source: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90 @ 5041

Last change on this file since 5041 was 4951, checked in by aborella, 4 months ago

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