source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 @ 5209

Last change on this file since 5209 was 5208, checked in by Laurent Fairhead, 7 hours ago

Louis d'Alencon modifications:
Il s'agit des modifs permettant de choisir de calculer la variance dans les thermiques par le nouveau modèle (iflag_ratqs = 10) ou de
continuer à calculer cette variance par la param d'Arnaud ce qui fait une version hybride avec variance pronostique dans l'environnement mais
variance diagnostique dans les thermiques (iflag_ratqs = 11). Le iflag_ratqs = 4 version standard continue de faire toutes les variances de façon
diagnostique.
File size: 73.0 KB
Line 
1MODULE lmdz_lscp
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs,pplay,temp,qt,qice_save,ptconv,ratqs,sigma_qtherm, &
10     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
11     pfraclr, pfracld,                                  &
12     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
13     radocond, radicefrac, rain, snow,                  &
14     frac_impa, frac_nucl, beta,                        &
15     prfl, psfl, rhcl, qta, fraca,                      &
16     tv, pspsk, tla, thl, iflag_cld_th,                 &
17     iflag_ice_thermo, distcltop, temp_cltop,           &
18     tke, tke_dissip,                                   &
19     cell_area,                                         &
20     cf_seri, rvc_seri, u_seri, v_seri,                 &
21     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
22     ratio_qi_qtot, dcf_sub, dcf_con, dcf_mix,          &
23     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
24     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
25     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
26     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
27     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
28     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
29     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
30     dqsmelt, dqsfreez)
31
32!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
34!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
35!--------------------------------------------------------------------------------
36! Date: 01/2021
37!--------------------------------------------------------------------------------
38! Aim: Large Scale Clouds and Precipitation (LSCP)
39!
40! This code is a new version of the fisrtilp.F90 routine, which is itself a
41! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
42! and 'ilp' (il pleut, L. Li)
43!
44! Compared to the original fisrtilp code, lscp
45! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
46! -> consider always precipitation thermalisation (fl_cor_ebil>0)
47! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
48! -> option "oldbug" by JYG has been removed
49! -> iflag_t_glace >0 always
50! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
51! -> rectangular distribution from L. Li no longer available
52! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
53!--------------------------------------------------------------------------------
54! References:
55!
56! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
57! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
58! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
59! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
60! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
61! - Touzze-Peifert Ludo, PhD thesis, p117-124
62! -------------------------------------------------------------------------------
63! Code structure:
64!
65! P0>     Thermalization of the precipitation coming from the overlying layer
66! P1>     Evaporation of the precipitation (falling from the k+1 level)
67! P2>     Cloud formation (at the k level)
68! P2.A.1> With the PDFs, calculation of cloud properties using the inital
69!         values of T and Q
70! P2.A.2> Coupling between condensed water and temperature
71! P2.A.3> Calculation of final quantities associated with cloud formation
72! P2.B>   Release of Latent heat after cloud formation
73! P3>     Autoconversion to precipitation (k-level)
74! P4>     Wet scavenging
75!------------------------------------------------------------------------------
76! Some preliminary comments (JBM) :
77!
78! The cloud water that the radiation scheme sees is not the same that the cloud
79! water used in the physics and the dynamics
80!
81! During the autoconversion to precipitation (P3 step), radocond (cloud water used
82! by the radiation scheme) is calculated as an average of the water that remains
83! in the cloud during the precipitation and not the water remaining at the end
84! of the time step. The latter is used in the rest of the physics and advected
85! by the dynamics.
86!
87! In summary:
88!
89! Radiation:
90! xflwc(newmicro)+xfiwc(newmicro) =
91!   radocond=lwcon(nc)+iwcon(nc)
92!
93! Notetheless, be aware of:
94!
95! radocond .NE. ocond(nc)
96! i.e.:
97! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
98! but oliq+(ocond-oliq) .EQ. ocond
99! (which is not trivial)
100!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101!
102
103! USE de modules contenant des fonctions.
104USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
105USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
106USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
107USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
108USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
109USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
110
111! Use du module lmdz_lscp_ini contenant les constantes
112USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
113USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
114USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
115USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
116USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
117USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
118USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
119USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
120USE lmdz_lscp_ini, ONLY : ok_poprecip
121USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
122
123IMPLICIT NONE
124
125!===============================================================================
126! VARIABLES DECLARATION
127!===============================================================================
128
129! INPUT VARIABLES:
130!-----------------
131
132  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
133  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
134  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
135
136  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
137  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
139  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
140  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [kg/kg]
141  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
142  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
143                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
144  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
145
146  !Inputs associated with thermal plumes
147
148  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
153  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
154  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
155
156  ! INPUT/OUTPUT variables
157  !------------------------
158 
159  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
160  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
161
162
163  ! INPUT/OUTPUT condensation and ice supersaturation
164  !--------------------------------------------------
165  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
166  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratio_qi_qtot    ! solid specific water to total specific water ratio [-]
167  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
170  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
171
172  ! INPUT/OUTPUT aviation
173  !--------------------------------------------------
174  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
175  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
176 
177  ! OUTPUT variables
178  !-----------------
179
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
186  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
187  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
188  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
191  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
192  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
193  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
194  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
195  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
196  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
197  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
201
202  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
203 
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
206 
207  ! for condensation and ice supersaturation
208
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
228
229  ! for contrails and aviation
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
239
240
241  ! for POPRECIP
242
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
256
257  ! for thermals
258
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
263
264
265  ! LOCAL VARIABLES:
266  !----------------
267  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
268  REAL :: zct, zcl,zexpo
269  REAL, DIMENSION(klon,klev) :: ctot
270  REAL, DIMENSION(klon,klev) :: ctot_vol
271  REAL, DIMENSION(klon) :: zqs, zdqs
272  REAL :: zdelta, zcor, zcvm5
273  REAL, DIMENSION(klon) :: zdqsdT_raw
274  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
275  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
276  REAL :: num,denom
277  REAL :: cste
278  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
279  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
280  REAL :: erf
281  REAL, DIMENSION(klon) :: zfice_th
282  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
283  REAL, DIMENSION(klon) :: zrfl, zrfln
284  REAL :: zqev, zqevt
285  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
286  REAL :: zqev0,zqevi, zqevti
287  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
288  REAL, DIMENSION(klon) :: zoliql, zoliqi
289  REAL, DIMENSION(klon) :: zt
290  REAL, DIMENSION(klon,klev) :: zrho
291  REAL, DIMENSION(klon) :: zdz,iwc
292  REAL :: zchau,zfroi
293  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
294  REAL :: zmelt,zrain,zsnow,zprecip
295  REAL, DIMENSION(klon) :: dzfice
296  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
297  REAL :: zsolid
298  REAL, DIMENSION(klon) :: qtot, qzero
299  REAL, DIMENSION(klon) :: dqsl,dqsi
300  REAL :: smallestreal
301  !  Variables for Bergeron process
302  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
303  REAL, DIMENSION(klon) :: zqpreci, zqprecl
304  ! Variables precipitation energy conservation
305  REAL, DIMENSION(klon) :: zmqc
306  REAL :: zalpha_tr
307  REAL :: zfrac_lessi
308  REAL, DIMENSION(klon) :: zprec_cond
309  REAL :: zmair
310  REAL :: zcpair, zcpeau
311  REAL, DIMENSION(klon) :: zlh_solid
312  REAL, DIMENSION(klon) :: ztupnew
313  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
314  REAL :: zm_solid         ! for liquid -> solid conversion
315  REAL, DIMENSION(klon) :: zrflclr, zrflcld
316  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
317  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
318  REAL, DIMENSION(klon) :: ziflclr, ziflcld
319  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
320  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
321  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
322  REAL, DIMENSION(klon,klev) :: velo
323  REAL :: vr, ffallv
324  REAL :: qlmpc, qimpc, rnebmpc
325  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
326  REAL :: effective_zneb
327  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
328  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
329 
330  ! for condensation and ice supersaturation
331  REAL, DIMENSION(klon) :: qvc, shear
332  REAL :: delta_z
333  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
334  ! Constants used for calculating ratios that are advected (using a parent-child
335  ! formalism). This is not done in the dynamical core because at this moment,
336  ! only isotopes can use this parent-child formalism. Note that the two constants
337  ! are the same as the one use in the dynamical core, being also defined in
338  ! dyn3d_common/infotrac.F90
339  REAL :: min_qParent, min_ratio
340
341  INTEGER i, k, n, kk, iter
342  INTEGER, DIMENSION(klon) :: n_i
343  INTEGER ncoreczq
344  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
345  LOGICAL iftop
346
347  LOGICAL, DIMENSION(klon) :: lognormale
348  LOGICAL, DIMENSION(klon) :: keepgoing
349
350  CHARACTER (len = 20) :: modname = 'lscp'
351  CHARACTER (len = 80) :: abort_message
352 
353
354!===============================================================================
355! INITIALISATION
356!===============================================================================
357
358! Few initial checks
359
360
361IF (iflag_fisrtilp_qsat .LT. 0) THEN
362    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
363    CALL abort_physic(modname,abort_message,1)
364ENDIF
365
366! Few initialisations
367
368znebprecip(:)=0.0
369ctot_vol(1:klon,1:klev)=0.0
370rneblsvol(1:klon,1:klev)=0.0
371smallestreal=1.e-9
372znebprecipclr(:)=0.0
373znebprecipcld(:)=0.0
374mpc_bl_points(:,:)=0
375
376IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
377
378! AA for 'safety' reasons
379zalpha_tr   = 0.
380zfrac_lessi = 0.
381beta(:,:)= 0.
382
383! Initialisation of variables:
384
385prfl(:,:) = 0.0
386psfl(:,:) = 0.0
387d_t(:,:) = 0.0
388d_q(:,:) = 0.0
389d_ql(:,:) = 0.0
390d_qi(:,:) = 0.0
391rneb(:,:) = 0.0
392pfraclr(:,:)=0.0
393pfracld(:,:)=0.0
394cldfraliq(:,:)=0.
395sigma2_icefracturb(:,:)=0.
396mean_icefracturb(:,:)=0.
397radocond(:,:) = 0.0
398radicefrac(:,:) = 0.0
399frac_nucl(:,:) = 1.0
400frac_impa(:,:) = 1.0
401rain(:) = 0.0
402snow(:) = 0.0
403zoliq(:)=0.0
404zfice(:)=0.0
405dzfice(:)=0.0
406zfice_turb(:)=0.0
407dzfice_turb(:)=0.0
408zqprecl(:)=0.0
409zqpreci(:)=0.0
410zrfl(:) = 0.0
411zifl(:) = 0.0
412ziflprev(:)=0.0
413zneb(:) = seuil_neb
414zrflclr(:) = 0.0
415ziflclr(:) = 0.0
416zrflcld(:) = 0.0
417ziflcld(:) = 0.0
418tot_zneb(:) = 0.0
419tot_znebn(:) = 0.0
420d_tot_zneb(:) = 0.0
421qzero(:) = 0.0
422zdistcltop(:)=0.0
423ztemp_cltop(:) = 0.0
424ztupnew(:)=0.0
425
426distcltop(:,:)=0.
427temp_cltop(:,:)=0.
428
429!--Ice supersaturation
430gamma_cond(:,:) = 1.
431qissr(:,:)      = 0.
432issrfra(:,:)    = 0.
433dcf_sub(:,:)    = 0.
434dcf_con(:,:)    = 0.
435dcf_mix(:,:)    = 0.
436dqi_adj(:,:)    = 0.
437dqi_sub(:,:)    = 0.
438dqi_con(:,:)    = 0.
439dqi_mix(:,:)    = 0.
440dqvc_adj(:,:)   = 0.
441dqvc_sub(:,:)   = 0.
442dqvc_con(:,:)   = 0.
443dqvc_mix(:,:)   = 0.
444fcontrN(:,:)    = 0.
445fcontrP(:,:)    = 0.
446Tcontr(:,:)     = missing_val
447qcontr(:,:)     = missing_val
448qcontr2(:,:)    = missing_val
449dcf_avi(:,:)    = 0.
450dqi_avi(:,:)    = 0.
451dqvc_avi(:,:)   = 0.
452qvc(:)          = 0.
453shear(:)        = 0.
454min_qParent     = 1.e-30
455min_ratio       = 1.e-16
456
457!-- poprecip
458qraindiag(:,:)= 0.
459qsnowdiag(:,:)= 0.
460dqreva(:,:)   = 0.
461dqrauto(:,:)  = 0.
462dqrmelt(:,:)  = 0.
463dqrfreez(:,:) = 0.
464dqrcol(:,:)   = 0.
465dqssub(:,:)   = 0.
466dqsauto(:,:)  = 0.
467dqsrim(:,:)   = 0.
468dqsagg(:,:)   = 0.
469dqsfreez(:,:) = 0.
470dqsmelt(:,:)  = 0.
471zqupnew(:)    = 0.
472zqvapclr(:)   = 0.
473
474
475
476!c_iso: variable initialisation for iso
477
478
479!===============================================================================
480!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
481!===============================================================================
482
483ncoreczq=0
484
485DO k = klev, 1, -1
486
487    IF (k.LE.klev-1) THEN
488       iftop=.false.
489    ELSE
490       iftop=.true.
491    ENDIF
492
493    ! Initialisation temperature and specific humidity
494    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
495    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
496    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
497    ! d_t = temperature tendency due to lscp
498    ! The temperature of the overlying layer is updated here because needed for thermalization
499    DO i = 1, klon
500        zt(i)=temp(i,k)
501        zq(i)=qt(i,k)
502        IF (.not. iftop) THEN
503           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
504           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
505           !--zqs(i) is the saturation specific humidity in the layer above
506           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
507        ENDIF
508        !c_iso init of iso
509    ENDDO
510   
511    !================================================================
512    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
513    IF (ok_poprecip) THEN
514
515            CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
516                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
517                              zqvapclr, zqupnew, &
518                              zrfl, zrflclr, zrflcld, &
519                              zifl, ziflclr, ziflcld, &
520                              dqreva(:,k),dqssub(:,k) &
521                             )
522
523    !================================================================
524    ELSE
525
526    ! --------------------------------------------------------------------
527    ! P1> Thermalization of precipitation falling from the overlying layer
528    ! --------------------------------------------------------------------
529    ! Computes air temperature variation due to enthalpy transported by
530    ! precipitation. Precipitation is then thermalized with the air in the
531    ! layer.
532    ! The precipitation should remain thermalized throughout the different
533    ! thermodynamical transformations.
534    ! The corresponding water mass should
535    ! be added when calculating the layer's enthalpy change with
536    ! temperature
537    ! See lmdzpedia page todoan
538    ! todoan: check consistency with ice phase
539    ! todoan: understand why several steps
540    ! ---------------------------------------------------------------------
541   
542    IF (iftop) THEN
543
544        DO i = 1, klon
545            zmqc(i) = 0.
546        ENDDO
547
548    ELSE
549
550        DO i = 1, klon
551 
552            zmair=(paprs(i,k)-paprs(i,k+1))/RG
553            ! no condensed water so cp=cp(vapor+dry air)
554            ! RVTMP2=rcpv/rcpd-1
555            zcpair=RCPD*(1.0+RVTMP2*zq(i))
556            zcpeau=RCPD*RVTMP2
557
558            ! zmqc: precipitation mass that has to be thermalized with
559            ! layer's air so that precipitation at the ground has the
560            ! same temperature as the lowermost layer
561            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
562            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
563            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
564             / (zcpair + zmqc(i)*zcpeau)
565
566        ENDDO
567 
568    ENDIF
569
570    ! --------------------------------------------------------------------
571    ! P2> Precipitation evaporation/sublimation/melting
572    ! --------------------------------------------------------------------
573    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
574    ! --------------------------------------------------------------------
575   
576    IF (iflag_evap_prec.GE.1) THEN
577
578        ! Calculation of saturation specific humidity
579        ! depending on temperature:
580        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
581        ! wrt liquid water
582        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
583        ! wrt ice
584        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
585
586        DO i = 1, klon
587
588            ! if precipitation
589            IF (zrfl(i)+zifl(i).GT.0.) THEN
590
591            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
592            ! c_iso: likely important to distinguish cs from neb isotope precipitation
593
594                IF (iflag_evap_prec.GE.4) THEN
595                    zrfl(i) = zrflclr(i)
596                    zifl(i) = ziflclr(i)
597                ENDIF
598
599                IF (iflag_evap_prec.EQ.1) THEN
600                    znebprecip(i)=zneb(i)
601                ELSE
602                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
603                ENDIF
604
605                IF (iflag_evap_prec.GT.4) THEN
606                ! Max evaporation not to saturate the clear sky precip fraction
607                ! i.e. the fraction where evaporation occurs
608                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
609                ELSEIF (iflag_evap_prec .EQ. 4) THEN
610                ! Max evaporation not to saturate the whole mesh
611                ! Pay attention -> lead to unrealistic and excessive evaporation
612                    zqev0 = MAX(0.0, zqs(i)-zq(i))
613                ELSE
614                ! Max evap not to saturate the fraction below the cloud
615                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
616                ENDIF
617
618                ! Evaporation of liquid precipitation coming from above
619                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
620                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
621                ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
622
623                IF (iflag_evap_prec.EQ.3) THEN
624                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
625                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
626                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
627                ELSE IF (iflag_evap_prec.GE.4) THEN
628                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
629                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
630                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
631                ELSE
632                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
633                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
634                ENDIF
635
636                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
637                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
638
639                ! sublimation of the solid precipitation coming from above
640                IF (iflag_evap_prec.EQ.3) THEN
641                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
642                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
643                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
644                ELSE IF (iflag_evap_prec.GE.4) THEN
645                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
646                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
647                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
648                ELSE
649                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
650                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
651                ENDIF
652
653                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
654                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
655
656                ! A. JAM
657                ! Evaporation limit: we ensure that the layer's fraction below
658                ! the cloud or the whole mesh (depending on iflag_evap_prec)
659                ! does not reach saturation. In this case, we
660                ! redistribute zqev0 conserving the ratio liquid/ice
661
662                IF (zqevt+zqevti.GT.zqev0) THEN
663                    zqev=zqev0*zqevt/(zqevt+zqevti)
664                    zqevi=zqev0*zqevti/(zqevt+zqevti)
665                ELSE
666                    zqev=zqevt
667                    zqevi=zqevti
668                ENDIF
669
670
671                ! New solid and liquid precipitation fluxes after evap and sublimation
672                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
673                         /RG/dtime)
674                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
675                         /RG/dtime)
676
677
678                ! vapor, temperature, precip fluxes update
679                ! vapor is updated after evaporation/sublimation (it is increased)
680                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
681                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
682                ! zmqc is the total condensed water in the precip flux (it is decreased)
683                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
684                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
685                ! air and precip temperature (i.e., gridbox temperature)
686                ! is updated due to latent heat cooling
687                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
688                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
689                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
690                + (zifln(i)-zifl(i))                      &
691                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
692                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
693
694                ! New values of liquid and solid precipitation
695                zrfl(i) = zrfln(i)
696                zifl(i) = zifln(i)
697
698                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
699                !                         due to evap + sublim                                     
700                                           
701
702                IF (iflag_evap_prec.GE.4) THEN
703                    zrflclr(i) = zrfl(i)
704                    ziflclr(i) = zifl(i)
705                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
706                        znebprecipclr(i) = 0.0
707                    ENDIF
708                    zrfl(i) = zrflclr(i) + zrflcld(i)
709                    zifl(i) = ziflclr(i) + ziflcld(i)
710                ENDIF
711
712                ! c_iso duplicate for isotopes or loop on isotopes
713
714                ! Melting:
715                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
716                ! precip fraction that is melted
717                zmelt = MIN(MAX(zmelt,0.),1.)
718
719                ! update of rainfall and snowfall due to melting
720                IF (iflag_evap_prec.GE.4) THEN
721                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
722                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
723                    zrfl(i)=zrflclr(i)+zrflcld(i)
724                ELSE
725                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
726                ENDIF
727
728
729                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
730
731                ! Latent heat of melting because of precipitation melting
732                ! NB: the air + precip temperature is simultaneously updated
733                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
734                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
735               
736                IF (iflag_evap_prec.GE.4) THEN
737                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
738                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
739                    zifl(i)=ziflclr(i)+ziflcld(i)
740                ELSE
741                    zifl(i)=zifl(i)*(1.-zmelt)
742                ENDIF
743
744            ELSE
745                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
746                znebprecip(i)=0.
747
748            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
749
750        ENDDO ! loop on klon
751
752    ENDIF ! (iflag_evap_prec>=1)
753
754    ENDIF ! (ok_poprecip)
755
756    ! --------------------------------------------------------------------
757    ! End precip evaporation
758    ! --------------------------------------------------------------------
759   
760    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
761    !-------------------------------------------------------
762
763     qtot(:)=zq(:)+zmqc(:)
764     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
765     DO i = 1, klon
766            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
767            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
768            IF (zq(i) .LT. 1.e-15) THEN
769                ncoreczq=ncoreczq+1
770                zq(i)=1.e-15
771            ENDIF
772            ! c_iso: do something similar for isotopes
773
774     ENDDO
775   
776    ! --------------------------------------------------------------------
777    ! P2> Cloud formation
778    !---------------------------------------------------------------------
779    !
780    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
781    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
782    ! is prescribed and depends on large scale variables and boundary layer
783    ! properties)
784    ! The decrease in condensed part due tu latent heating is taken into
785    ! account
786    ! -------------------------------------------------------------------
787
788        ! P2.1> With the PDFs (log-normal, bigaussian)
789        ! cloud properties calculation with the initial values of t and q
790        ! ----------------------------------------------------------------
791
792        ! initialise gammasat and qincloud_mpc
793        gammasat(:)=1.
794        qincloud_mpc(:)=0.
795
796        IF (iflag_cld_th.GE.5) THEN
797               ! Cloud cover and content in meshes affected by shallow convection,
798               ! are retrieved from a bi-gaussian distribution of the saturation deficit
799               ! following Jam et al. 2013
800
801               IF (iflag_cloudth_vert.LE.2) THEN
802                  ! Old version of Arnaud Jam
803
804                    CALL cloudth(klon,klev,k,tv,                  &
805                         zq,qta,fraca,                            &
806                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
807                         ratqs,zqs,temp,                              &
808                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
809
810
811                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
812                   ! Default version of Arnaud Jam
813                       
814                    CALL cloudth_v3(klon,klev,k,tv,                        &
815                         zq,qta,fraca,                                     &
816                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
817                         ratqs,sigma_qtherm,zqs,temp, &
818                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
819
820
821                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
822                   ! Jean Jouhaud's version, with specific separation between surface and volume
823                   ! cloud fraction Decembre 2018
824
825                    CALL cloudth_v6(klon,klev,k,tv,                        &
826                         zq,qta,fraca,                                     &
827                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
828                         ratqs,zqs,temp, &
829                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
830
831                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
832                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
833                   ! for boundary-layer mixed phase clouds
834                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
835                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
836                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
837                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
838
839                ENDIF
840
841
842                DO i=1,klon
843                    rneb(i,k)=ctot(i,k)
844                    rneblsvol(i,k)=ctot_vol(i,k)
845                    zqn(i)=qcloud(i)
846                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
847                    qvc(i) = rneb(i,k) * zqs(i)
848                ENDDO
849
850        ENDIF
851
852        IF (iflag_cld_th .LE. 4) THEN
853           
854                ! lognormal
855            lognormale(:) = .TRUE.
856
857        ELSEIF (iflag_cld_th .GE. 6) THEN
858           
859                ! lognormal distribution when no thermals
860            lognormale(:) = fraca(:,k) < min_frac_th_cld
861
862        ELSE
863                ! When iflag_cld_th=5, we always assume
864                ! bi-gaussian distribution
865            lognormale(:) = .FALSE.
866       
867        ENDIF
868
869        DT(:) = 0.
870        n_i(:)=0
871        Tbef(:)=zt(:)
872        qlbef(:)=0.
873
874        ! Treatment of non-boundary layer clouds (lognormale)
875        ! condensation with qsat(T) variation (adaptation)
876        ! Iterative resolution to converge towards qsat
877        ! with update of temperature, ice fraction and qsat at
878        ! each iteration
879
880        ! todoan -> sensitivity to iflag_fisrtilp_qsat
881        DO iter=1,iflag_fisrtilp_qsat+1
882
883                keepgoing(:) = .FALSE.
884
885                DO i=1,klon
886
887                    ! keepgoing = .true. while convergence is not satisfied
888
889                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
890
891                        ! if not convergence:
892                        ! we calculate a new iteration
893                        keepgoing(i) = .TRUE.
894
895                        ! P2.2.1> cloud fraction and condensed water mass calculation
896                        ! Calculated variables:
897                        ! rneb : cloud fraction
898                        ! zqn : total water within the cloud
899                        ! zcond: mean condensed water within the mesh
900                        ! rhcl: clear-sky relative humidity
901                        !---------------------------------------------------------------
902
903                        ! new temperature that only serves in the iteration process:
904                        Tbef(i)=Tbef(i)+DT(i)
905
906                        ! Rneb, qzn and zcond for lognormal PDFs
907                        qtot(i)=zq(i)+zmqc(i)
908
909                      ENDIF
910
911                  ENDDO
912
913                  ! Calculation of saturation specific humidity and ice fraction
914                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
915                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
916                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
917                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
918                  ! cloud phase determination
919                  IF (iflag_t_glace.GE.4) THEN
920                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
921                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
922                  ENDIF
923
924                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
925
926                  !--AB Activates a condensation scheme that allows for
927                  !--ice supersaturation and contrails evolution from aviation
928                  IF (ok_ice_supersat) THEN
929
930                    !--Calculate the shear value (input for condensation and ice supersat)
931                    DO i = 1, klon
932                      !--Cell thickness [m]
933                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
934                      IF ( iftop ) THEN
935                        ! top
936                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
937                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
938                      ELSEIF ( k .EQ. 1 ) THEN
939                        ! surface
940                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
941                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
942                      ELSE
943                        ! other layers
944                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
945                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
946                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
947                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
948                      ENDIF
949                    ENDDO
950
951                    !---------------------------------------------
952                    !--   CONDENSATION AND ICE SUPERSATURATION  --
953                    !---------------------------------------------
954
955                    CALL condensation_ice_supersat( &
956                        klon, dtime, missing_val, &
957                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
958                        cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
959                        shear(:), tke_dissip(:,k), cell_area(:), &
960                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
961                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
962                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
963                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
964                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
965                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
966                        flight_dist(:,k), flight_h2o(:,k), &
967                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
968
969
970                  !--If .TRUE., calls an externalised version of the generalised
971                  !--lognormal condensation scheme (Bony and Emanuel 2001)
972                  !--Numerically, convergence is conserved with this option
973                  !--The objective is to simplify LSCP
974                  ELSEIF ( ok_external_lognormal ) THEN
975                         
976                   CALL condensation_lognormal( &
977                       klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &
978                       keepgoing(:), rneb(:,k), zqn(:), qvc(:))
979
980
981                 ELSE !--Generalised lognormal (Bony and Emanuel 2001)
982
983                  DO i=1,klon !todoan : check if loop in i is needed
984
985                      IF (keepgoing(i)) THEN
986
987                        zpdf_sig(i)=ratqs(i,k)*zq(i)
988                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
989                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
990                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
991                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
992                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
993                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
994                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
995                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
996                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
997                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
998
999                          IF (zpdf_e1(i).LT.1.e-10) THEN
1000                              rneb(i,k)=0.
1001                              zqn(i)=zqs(i)
1002                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
1003                              qvc(i) = 0.
1004                          ELSE
1005                              rneb(i,k)=0.5*zpdf_e1(i)
1006                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
1007                              !--AB grid-mean vapor in the cloud - we assume saturation adjustment
1008                              qvc(i) = rneb(i,k) * zqs(i)
1009                          ENDIF
1010
1011                      ENDIF ! keepgoing
1012                  ENDDO ! loop on klon
1013
1014                  ENDIF ! .NOT. ok_ice_supersat
1015
1016                  DO i=1,klon
1017                      IF (keepgoing(i)) THEN
1018
1019                        ! If vertical heterogeneity, change fraction by volume as well
1020                        IF (iflag_cloudth_vert.GE.3) THEN
1021                            ctot_vol(i,k)=rneb(i,k)
1022                            rneblsvol(i,k)=ctot_vol(i,k)
1023                        ENDIF
1024
1025
1026                       ! P2.2.2> Approximative calculation of temperature variation DT
1027                       !           due to condensation.
1028                       ! Calculated variables:
1029                       ! dT : temperature change due to condensation
1030                       !---------------------------------------------------------------
1031
1032               
1033                        IF (zfice(i).LT.1) THEN
1034                            cste=RLVTT
1035                        ELSE
1036                            cste=RLSTT
1037                        ENDIF
1038                       
1039                        ! LEA_R : check formule
1040                        IF ( ok_unadjusted_clouds ) THEN
1041                          !--AB We relax the saturation adjustment assumption
1042                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1043                          IF ( rneb(i,k) .GT. eps ) THEN
1044                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
1045                          ELSE
1046                            qlbef(i) = 0.
1047                          ENDIF
1048                        ELSE
1049                          qlbef(i)=max(0.,zqn(i)-zqs(i))
1050                        ENDIF
1051
1052                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1053                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
1054                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1055                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
1056                              *qlbef(i)*dzfice(i)
1057                        ! here we update a provisory temperature variable that only serves in the iteration
1058                        ! process
1059                        DT(i)=num/denom
1060                        n_i(i)=n_i(i)+1
1061
1062                    ENDIF ! end keepgoing
1063
1064                ENDDO     ! end loop on i
1065
1066        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
1067
1068        ! P2.3> Final quantities calculation
1069        ! Calculated variables:
1070        !   rneb : cloud fraction
1071        !   zcond: mean condensed water in the mesh
1072        !   zqn  : mean water vapor in the mesh
1073        !   zfice: ice fraction in clouds
1074        !   zt   : temperature
1075        !   rhcl : clear-sky relative humidity
1076        ! ----------------------------------------------------------------
1077
1078
1079        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
1080        IF (iflag_t_glace.GE.4) THEN
1081           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
1082           distcltop(:,k)=zdistcltop(:)
1083           temp_cltop(:,k)=ztemp_cltop(:)
1084        ENDIF
1085
1086        ! Partition function depending on temperature
1087        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
1088
1089        ! Partition function depending on tke for non shallow-convective clouds
1090        IF (iflag_icefrac .GE. 1) THEN
1091
1092           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
1093           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
1094
1095        ENDIF
1096
1097        ! Water vapor update, Phase determination and subsequent latent heat exchange
1098        DO i=1, klon
1099            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
1100            ! iflag_cloudth_vert=7 and specific param is activated
1101            IF (mpc_bl_points(i,k) .GT. 0) THEN
1102                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
1103                ! following line is very strange and probably wrong
1104                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
1105                ! water vapor update and partition function if thermals
1106                zq(i) = zq(i) - zcond(i)       
1107                zfice(i)=zfice_th(i)
1108            ELSE
1109                ! Checks on rneb, rhcl and zqn
1110                IF (rneb(i,k) .LE. 0.0) THEN
1111                    zqn(i) = 0.0
1112                    rneb(i,k) = 0.0
1113                    zcond(i) = 0.0
1114                    rhcl(i,k)=zq(i)/zqs(i)
1115                ELSE IF (rneb(i,k) .GE. 1.0) THEN
1116                    zqn(i) = zq(i)
1117                    rneb(i,k) = 1.0
1118                    IF ( ok_unadjusted_clouds ) THEN
1119                      !--AB We relax the saturation adjustment assumption
1120                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1121                      zcond(i) = MAX(0., zqn(i) - qvc(i))
1122                    ELSE
1123                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1124                    ENDIF
1125                    rhcl(i,k)=1.0
1126                ELSE
1127                    IF ( ok_unadjusted_clouds ) THEN
1128                      !--AB We relax the saturation adjustment assumption
1129                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1130                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
1131                    ELSE
1132                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1133                    ENDIF
1134                    ! following line is very strange and probably wrong:
1135                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
1136                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
1137                    IF (iflag_icefrac .GE. 1) THEN
1138                        IF (lognormale(i)) THEN 
1139                           zcond(i)  = zqliq(i) + zqice(i)
1140                           zfice(i)=zfice_turb(i)
1141                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
1142                        ENDIF
1143                    ENDIF
1144                ENDIF
1145
1146                ! water vapor update
1147                zq(i) = zq(i) - zcond(i)
1148
1149            ENDIF
1150
1151            ! c_iso : routine that computes in-cloud supersaturation
1152            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
1153
1154            ! temperature update due to phase change
1155            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1156            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1157                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1158        ENDDO
1159
1160        ! If vertical heterogeneity, change volume fraction
1161        IF (iflag_cloudth_vert .GE. 3) THEN
1162          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1163          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1164        ENDIF
1165   
1166    !--AB Write diagnostics and tracers for ice supersaturation
1167    IF ( ok_ice_supersat ) THEN
1168      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1169      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1170
1171      DO i = 1, klon
1172
1173        cf_seri(i,k) = rneb(i,k)
1174
1175        IF ( .NOT. ok_unadjusted_clouds ) THEN
1176          qvc(i) = zqs(i) * rneb(i,k)
1177        ENDIF
1178        IF ( zq(i) .GT. min_qParent ) THEN
1179          rvc_seri(i,k) = qvc(i) / zq(i)
1180        ELSE
1181          rvc_seri(i,k) = min_ratio
1182        ENDIF
1183        !--The MIN barrier is NEEDED because of:
1184        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1185        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1186        !--The MAX barrier is a safeguard that should not be activated
1187        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1188
1189        !--Diagnostics
1190        gamma_cond(i,k) = gammasat(i)
1191        IF ( issrfra(i,k) .LT. eps ) THEN
1192          issrfra(i,k) = 0.
1193          qissr(i,k) = 0.
1194        ENDIF
1195        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1196        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1197        IF ( subfra(i,k) .LT. eps ) THEN
1198          subfra(i,k) = 0.
1199          qsub(i,k) = 0.
1200        ENDIF
1201        qcld(i,k) = qvc(i) + zcond(i)
1202        IF ( cf_seri(i,k) .LT. eps ) THEN
1203          qcld(i,k) = 0.
1204        ENDIF
1205      ENDDO
1206    ENDIF
1207
1208
1209    ! ----------------------------------------------------------------
1210    ! P3> Precipitation formation
1211    ! ----------------------------------------------------------------
1212
1213    !================================================================
1214    IF (ok_poprecip) THEN
1215
1216      DO i = 1, klon
1217        zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1218        zoliqi(i) = zcond(i) * zfice(i)
1219      ENDDO
1220
1221      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1222                            ctot_vol(:,k), ptconv(:,k), &
1223                            zt, zq, zoliql, zoliqi, zfice, &
1224                            rneb(:,k), znebprecipclr, znebprecipcld, &
1225                            zrfl, zrflclr, zrflcld, &
1226                            zifl, ziflclr, ziflcld, &
1227                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1228                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1229                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1230                            dqsmelt(:,k), dqsfreez(:,k) &
1231                            )
1232
1233      DO i = 1, klon
1234        zoliq(i) = zoliql(i) + zoliqi(i)
1235        IF ( zoliq(i) .GT. 0. ) THEN
1236                zfice(i) = zoliqi(i)/zoliq(i)
1237        ELSE 
1238                zfice(i) = 0.0
1239        ENDIF
1240
1241        ! calculation of specific content of condensates seen by radiative scheme
1242        IF (ok_radocond_snow) THEN
1243           radocond(i,k) = zoliq(i)
1244           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1245           radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
1246        ELSE
1247           radocond(i,k) = zoliq(i)
1248           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1249           radocondi(i,k)= radocond(i,k)*zfice(i)
1250        ENDIF
1251      ENDDO
1252
1253    !================================================================
1254    ELSE
1255
1256    ! LTP:
1257    IF (iflag_evap_prec .GE. 4) THEN
1258
1259        !Partitionning between precipitation coming from clouds and that coming from CS
1260
1261        !0) Calculate tot_zneb, total cloud fraction above the cloud
1262        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
1263       
1264        DO i=1, klon
1265                tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
1266                        /(1.-min(zneb(i),1.-smallestreal))
1267                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
1268                tot_zneb(i) = tot_znebn(i)
1269
1270
1271                !1) Cloudy to clear air
1272                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
1273                IF (znebprecipcld(i) .GT. 0.) THEN
1274                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
1275                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
1276                ELSE
1277                        d_zrfl_cld_clr(i) = 0.
1278                        d_zifl_cld_clr(i) = 0.
1279                ENDIF
1280
1281                !2) Clear to cloudy air
1282                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
1283                IF (znebprecipclr(i) .GT. 0) THEN
1284                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
1285                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
1286                ELSE
1287                        d_zrfl_clr_cld(i) = 0.
1288                        d_zifl_clr_cld(i) = 0.
1289                ENDIF
1290
1291                !Update variables
1292                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
1293                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
1294                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
1295                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
1296                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
1297                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
1298
1299                ! c_iso  do the same thing for isotopes precip
1300        ENDDO
1301    ENDIF
1302
1303
1304    ! Autoconversion
1305    ! -------------------------------------------------------------------------------
1306
1307
1308    ! Initialisation of zoliq and radocond variables
1309
1310    DO i = 1, klon
1311            zoliq(i) = zcond(i)
1312            zoliqi(i)= zoliq(i)*zfice(i)
1313            zoliql(i)= zoliq(i)*(1.-zfice(i))
1314            ! c_iso : initialisation of zoliq* also for isotopes
1315            zrho(i,k)  = pplay(i,k) / zt(i) / RD
1316            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
1317            iwc(i)   = 0.
1318            zneb(i)  = MAX(rneb(i,k),seuil_neb)
1319            radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
1320            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
1321            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
1322    ENDDO
1323
1324       
1325    DO n = 1, niter_lscp
1326
1327        DO i=1,klon
1328            IF (rneb(i,k).GT.0.0) THEN
1329                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
1330            ENDIF
1331        ENDDO
1332
1333        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
1334
1335        DO i = 1, klon
1336
1337            IF (rneb(i,k).GT.0.0) THEN
1338
1339                ! Initialization of zrain, zsnow and zprecip:
1340                zrain=0.
1341                zsnow=0.
1342                zprecip=0.
1343                ! c_iso same init for isotopes. Externalisation?
1344
1345                IF (zneb(i).EQ.seuil_neb) THEN
1346                    zprecip = 0.0
1347                    zsnow = 0.0
1348                    zrain= 0.0
1349                ELSE
1350
1351                    IF (ptconv(i,k)) THEN ! if convective point
1352                        zcl=cld_lc_con
1353                        zct=1./cld_tau_con
1354                        zexpo=cld_expo_con
1355                        ffallv=ffallv_con
1356                    ELSE
1357                        zcl=cld_lc_lsc
1358                        zct=1./cld_tau_lsc
1359                        zexpo=cld_expo_lsc
1360                        ffallv=ffallv_lsc
1361                    ENDIF
1362
1363
1364                    ! if vertical heterogeneity is taken into account, we use
1365                    ! the "true" volume fraction instead of a modified
1366                    ! surface fraction (which is larger and artificially
1367                    ! reduces the in-cloud water).
1368
1369                    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
1370                    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
1371                    !.........................................................
1372                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
1373
1374                    ! if vertical heterogeneity is taken into account, we use
1375                    ! the "true" volume fraction instead of a modified
1376                    ! surface fraction (which is larger and artificially
1377                    ! reduces the in-cloud water).
1378                        effective_zneb=ctot_vol(i,k)
1379                    ELSE
1380                        effective_zneb=zneb(i)
1381                    ENDIF
1382
1383
1384                    IF (iflag_autoconversion .EQ. 2) THEN
1385                    ! two-steps resolution with niter_lscp=1 sufficient
1386                    ! we first treat the second term (with exponential) in an explicit way
1387                    ! and then treat the first term (-q/tau) in an exact way
1388                        zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
1389                        *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
1390                    ELSE
1391                    ! old explicit resolution with subtimesteps
1392                        zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
1393                        *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
1394                    ENDIF
1395
1396
1397                    ! Ice water quantity to remove (Zender & Kiehl, 1997)
1398                    ! dqice/dt=1/rho*d(rho*wice*qice)/dz
1399                    !....................................
1400                    IF (iflag_autoconversion .EQ. 2) THEN
1401                    ! exact resolution, niter_lscp=1 is sufficient but works only
1402                    ! with iflag_vice=0
1403                       IF (zoliqi(i) .GT. 0.) THEN
1404                          zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
1405                          +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
1406                       ELSE
1407                          zfroi=0.
1408                       ENDIF
1409                    ELSE
1410                    ! old explicit resolution with subtimesteps
1411                       zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k)
1412                    ENDIF
1413
1414                    zrain   = MIN(MAX(zchau,0.0),zoliql(i))
1415                    zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
1416                    zprecip = MAX(zrain + zsnow,0.0)
1417
1418                ENDIF
1419
1420                IF (iflag_autoconversion .GE. 1) THEN
1421                   ! debugged version with phase conservation through the autoconversion process
1422                   zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
1423                   zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
1424                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1425                ELSE
1426                   ! bugged version with phase resetting
1427                   zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
1428                   zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
1429                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1430                ENDIF
1431
1432                ! c_iso: call isotope_conversion (for readibility) or duplicate
1433
1434                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
1435                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
1436                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
1437
1438            ENDIF ! rneb >0
1439
1440        ENDDO  ! i = 1,klon
1441
1442    ENDDO      ! n = 1,niter
1443
1444    ! Precipitation flux calculation
1445
1446    DO i = 1, klon
1447           
1448            IF (iflag_evap_prec.GE.4) THEN
1449                ziflprev(i)=ziflcld(i)
1450            ELSE
1451                ziflprev(i)=zifl(i)*zneb(i)
1452            ENDIF
1453
1454            IF (rneb(i,k) .GT. 0.0) THEN
1455
1456               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1457               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1458               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1459               ! taken into account through a linearization of the equations and by approximating
1460               ! the liquid precipitation process with a "threshold" process. We assume that
1461               ! condensates are not modified during this operation. Liquid precipitation is
1462               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1463               ! Water vapor increases as well
1464               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1465
1466                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1467                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1468                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1469                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1470                    ! Computation of DT if all the liquid precip freezes
1471                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1472                    ! T should not exceed the freezing point
1473                    ! that is Delta > RTT-zt(i)
1474                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1475                    zt(i)  = zt(i) + DeltaT
1476                    ! water vaporization due to temp. increase
1477                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1478                    ! we add this vaporized water to the total vapor and we remove it from the precip
1479                    zq(i) = zq(i) +  Deltaq
1480                    ! The three "max" lines herebelow protect from rounding errors
1481                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1482                    ! liquid precipitation converted to ice precip
1483                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1484                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1485                    ! iced water budget
1486                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1487
1488                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1489
1490                IF (iflag_evap_prec.GE.4) THEN
1491                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1492                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1493                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1494                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1495                    znebprecipcld(i) = rneb(i,k)
1496                    zrfl(i) = zrflcld(i) + zrflclr(i)
1497                    zifl(i) = ziflcld(i) + ziflclr(i)
1498                ELSE
1499                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1500                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1501                    zifl(i) = zifl(i)+ zqpreci(i)        &
1502                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1503                ENDIF
1504                ! c_iso : same for isotopes
1505 
1506           ENDIF ! rneb>0
1507
1508   ENDDO
1509
1510    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1511    ! if iflag_evap_prec>=4
1512    IF (iflag_evap_prec.GE.4) THEN
1513
1514        DO i=1,klon
1515
1516            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1517                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1518                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1519            ELSE
1520                znebprecipclr(i)=0.0
1521            ENDIF
1522
1523            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1524                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1525                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1526            ELSE
1527                znebprecipcld(i)=0.0
1528            ENDIF
1529        ENDDO
1530
1531    ENDIF
1532
1533
1534    ENDIF ! ok_poprecip
1535
1536    ! End of precipitation formation
1537    ! --------------------------------
1538
1539
1540    ! Calculation of cloud condensates amount seen by radiative scheme
1541    !-----------------------------------------------------------------
1542
1543    ! Calculation of the concentration of condensates seen by the radiation scheme
1544    ! for liquid phase, we take radocondl
1545    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
1546    ! we recalculate radocondi to account for contributions from the precipitation flux
1547    ! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
1548
1549    IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
1550        ! for the solid phase (crystals + snowflakes)
1551        ! we recalculate radocondi by summing
1552        ! the ice content calculated in the mesh
1553        ! + the contribution of the non-evaporated snowfall
1554        ! from the overlying layer
1555        DO i=1,klon
1556           IF (ziflprev(i) .NE. 0.0) THEN
1557              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1558           ELSE
1559              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1560           ENDIF
1561           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
1562        ENDDO
1563    ENDIF
1564
1565    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
1566    DO i=1,klon
1567        IF (radocond(i,k) .GT. 0.) THEN
1568            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
1569        ENDIF
1570    ENDDO
1571
1572    ! ----------------------------------------------------------------
1573    ! P4> Wet scavenging
1574    ! ----------------------------------------------------------------
1575
1576    !Scavenging through nucleation in the layer
1577
1578    DO i = 1,klon
1579       
1580        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1581            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1582        ELSE
1583            beta(i,k) = 0.
1584        ENDIF
1585
1586        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1587
1588        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1589
1590            IF (temp(i,k) .GE. t_glace_min) THEN
1591                zalpha_tr = a_tr_sca(3)
1592            ELSE
1593                zalpha_tr = a_tr_sca(4)
1594            ENDIF
1595
1596            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1597            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1598
1599            ! Nucleation with a factor of -1 instead of -0.5
1600            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1601
1602        ENDIF
1603
1604    ENDDO
1605
1606    ! Scavenging through impaction in the underlying layer
1607
1608    DO kk = k-1, 1, -1
1609
1610        DO i = 1, klon
1611
1612            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1613
1614                IF (temp(i,kk) .GE. t_glace_min) THEN
1615                    zalpha_tr = a_tr_sca(1)
1616                ELSE
1617                    zalpha_tr = a_tr_sca(2)
1618                ENDIF
1619
1620              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1621              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1622
1623            ENDIF
1624
1625        ENDDO
1626
1627    ENDDO
1628
1629    ! Outputs:
1630    !-------------------------------
1631    ! Precipitation fluxes at layer interfaces
1632    ! + precipitation fractions +
1633    ! temperature and water species tendencies
1634    DO i = 1, klon
1635        psfl(i,k)=zifl(i)
1636        prfl(i,k)=zrfl(i)
1637        pfraclr(i,k)=znebprecipclr(i)
1638        pfracld(i,k)=znebprecipcld(i)
1639        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1640        d_qi(i,k) = zfice(i)*zoliq(i)
1641        d_q(i,k) = zq(i) - qt(i,k)
1642        ! c_iso: same for isotopes
1643        d_t(i,k) = zt(i) - temp(i,k)
1644    ENDDO
1645
1646
1647ENDDO
1648
1649
1650  ! Rain or snow at the surface (depending on the first layer temperature)
1651  DO i = 1, klon
1652      snow(i) = zifl(i)
1653      rain(i) = zrfl(i)
1654      ! c_iso final output
1655  ENDDO
1656
1657  IF (ncoreczq>0) THEN
1658      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1659  ENDIF
1660
1661
1662RETURN
1663
1664END SUBROUTINE lscp
1665!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1666
1667END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.