source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_main.f90 @ 5728

Last change on this file since 5728 was 5728, checked in by aborella, 5 months ago

Correction to rare floating point errors + added support for the sedimentation of contrails
File size: 74.6 KB
Line 
1MODULE lmdz_lscp_main
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon, klev, dtime, missing_val,         &
9     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
10     ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, &
11     qccon_old, cfcon, qvcon, qccon,                    &
12     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
13     pfraclr, pfracld,                                  &
14     cldfraliq, cldfraliqth,                            &
15     sigma2_icefracturb,sigma2_icefracturbth,           &
16     mean_icefracturb,mean_icefracturbth,               &
17     radocond, radicefrac, rain, snow,                  &
18     frac_impa, frac_nucl, beta,                        &
19     prfl, psfl, rhcl, qta, fraca,                      &
20     tv, pspsk, tla, thl, wth, iflag_cld_th,            &
21     iflag_ice_thermo, distcltop, temp_cltop,           &
22     tke, tke_dissip,                                   &
23     entr_therm, detr_therm,                            &
24     cell_area, stratomask,                             &
25     cf_seri, qvc_seri, u_seri, v_seri,                 &
26     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
27     dcf_sub, dcf_con, dcf_mix, dqised, dcfsed, dqvcsed,&
28     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
29     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
30     cfl_seri, cfc_seri, qtl_seri, qtc_seri,            &
31     qice_lincont, qice_circont, flight_dist,           &
32     flight_h2o, qradice_lincont, qradice_circont,      &
33     Tcritcont, qcritcont, potcontfraP, potcontfraNP,   &
34     cloudth_sth,                                       &
35     cloudth_senv, cloudth_sigmath, cloudth_sigmaenv,   &
36     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
37     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
38     dqsmelt, dqsfreez)
39
40!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
42!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
43!--------------------------------------------------------------------------------
44! Date: 01/2021
45!--------------------------------------------------------------------------------
46! Aim: Large Scale Clouds and Precipitation (LSCP)
47!
48! This code is a new version of the fisrtilp.F90 routine, which is itself a
49! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
50! and 'ilp' (il pleut, L. Li)
51!
52! Compared to the original fisrtilp code, lscp
53! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
54! -> consider always precipitation thermalisation (fl_cor_ebil>0)
55! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
56! -> option "oldbug" by JYG has been removed
57! -> iflag_t_glace >0 always
58! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
59! -> rectangular distribution from L. Li no longer available
60! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
61!--------------------------------------------------------------------------------
62! References:
63!
64! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
65! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
66! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
67! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
68! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
69! - Touzze-Peifert Ludo, PhD thesis, p117-124
70! -------------------------------------------------------------------------------
71! Code structure:
72!
73! P0>     Thermalization of the precipitation coming from the overlying layer
74! P1>     Evaporation of the precipitation (falling from the k+1 level)
75! P2>     Cloud formation (at the k level)
76! P2.A.1> With the PDFs, calculation of cloud properties using the inital
77!         values of T and Q
78! P2.A.2> Coupling between condensed water and temperature
79! P2.A.3> Calculation of final quantities associated with cloud formation
80! P2.B>   Release of Latent heat after cloud formation
81! P3>     Autoconversion to precipitation (k-level)
82! P4>     Wet scavenging
83!------------------------------------------------------------------------------
84! Some preliminary comments (JBM) :
85!
86! The cloud water that the radiation scheme sees is not the same that the cloud
87! water used in the physics and the dynamics
88!
89! During the autoconversion to precipitation (P3 step), radocond (cloud water used
90! by the radiation scheme) is calculated as an average of the water that remains
91! in the cloud during the precipitation and not the water remaining at the end
92! of the time step. The latter is used in the rest of the physics and advected
93! by the dynamics.
94!
95! In summary:
96!
97! Radiation:
98! xflwc(newmicro)+xfiwc(newmicro) =
99!   radocond=lwcon(nc)+iwcon(nc)
100!
101! Notetheless, be aware of:
102!
103! radocond .NE. ocond(nc)
104! i.e.:
105! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
106! but oliq+(ocond-oliq) .EQ. ocond
107! (which is not trivial)
108!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
109!
110
111! USE de modules contenant des fonctions.
112USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
113USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
114USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
115USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth
116USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
117USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
118
119! Use du module lmdz_lscp_ini contenant les constantes
120USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
121USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
122USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
123USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
124USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater
125USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
126USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
127USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
128USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato
129USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_precip_contrails, ok_ice_sedim
130USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_nodeep_lscp_rad
131USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth
132
133! Temporary call for Lamquin et al (2012) diagnostics
134USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
135USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
136USE phys_local_var_mod, ONLY : dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub
137USE phys_local_var_mod, ONLY : dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix
138USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix
139USE phys_local_var_mod, ONLY : dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed
140USE geometry_mod, ONLY: longitude_deg, latitude_deg
141
142IMPLICIT NONE
143
144!===============================================================================
145! VARIABLES DECLARATION
146!===============================================================================
147
148! INPUT VARIABLES:
149!-----------------
150
151  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
152  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
153 
154
155  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
156  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
157
158  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
159  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
160  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
161  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
162  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
163  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
164  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
165  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke             ! turbulent kinetic energy [m2/s2]
166  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip      ! TKE dissipation [m2/s3]
167  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: entr_therm      ! thermal plume entrainment rate * dz [kg/s/m2]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: detr_therm      ! thermal plume detrainment rate * dz [kg/s/m2]
169
170
171 
172  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
173  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
174  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
175  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
176  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
177  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
178  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon           ! in-cloud condensed specific humidity from deep convection [kg/kg]
179
180  !Inputs associated with thermal plumes
181
182  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
183  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
184  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
185  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
186  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
187  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid potential temperature within thermals [K]
188  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: wth                 ! vertical velocity within thermals [m/s]
189  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: sigma_qtherm        ! controls the saturation deficit distribution width in thermals [-]
190
191
192
193  ! INPUT/OUTPUT variables
194  !------------------------
195 
196  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
197  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function that sets the large-scale water distribution width
198
199
200  ! INPUT/OUTPUT condensation and ice supersaturation
201  !--------------------------------------------------
202  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
203  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvc_seri         ! cloudy water vapor [kg/kg]
204  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
205  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
206  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
207  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: stratomask       ! fraction of stratosphere (0 or 1)
208
209  ! INPUT/OUTPUT aviation
210  !--------------------------------------------------
211  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfl_seri         ! linear contrails fraction [-]
212  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfc_seri         ! contrail cirrus fraction [-]
213  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtl_seri         ! linear contrails total specific humidity [kg/kg]
214  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtc_seri         ! contrail cirrus total specific humidity [kg/kg]
215  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
216  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
217
218  ! OUTPUT variables
219  !-----------------
220
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
228  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
229  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq        ! liquid fraction of cloud fraction [-]
230  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliqth      ! liquid fraction of cloud fraction in thermals [-]
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturbth   ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
238  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
239  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
240  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
241  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
242  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
245
246  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
247 
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
250 
251  ! for condensation and ice supersaturation
252
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
256  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
257  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
258  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
263  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
269  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
270  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqised         !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
271  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
272  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
273  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
274  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
275
276  ! for contrails and aviation
277
278  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_lincont   !--condensed water in linear contrails [kg/kg]
279  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_circont   !--condensed water in contrail cirrus [kg/kg]
280  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_lincont!--condensed water in linear contrails used in the radiation scheme [kg/kg]
281  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_circont!--condensed water in contrail cirrus used in the radiation scheme [kg/kg]
282  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
283  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
284  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraP    !--potential persistent contrail fraction [-]
285  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraNP   !--potential non-persistent contrail fraction [-]
286
287
288  ! for POPRECIP
289
290  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
291  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
292  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
293  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
294  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
295  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
296  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
297  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
298  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
299  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
300  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
301  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
302  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
303
304  ! for thermals
305
306  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
307  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
308  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
309  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
310
311
312  ! LOCAL VARIABLES:
313  !----------------
314  REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in
315  REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol
316  REAL, DIMENSION(klon,klev) :: wls                                 !-- large scalce vertical velocity [m/s]
317  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
318  REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith
319  REAL :: zdelta
320  REAL, DIMENSION(klon) :: zdqsdT_raw
321  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
322  REAL, DIMENSION(klon) :: Tbef,Tbefth,Tbefthm1,qlibef,DT                  ! temperature, humidity and temp. variation during condensation iteration
323  REAL :: num,denom
324  REAL :: cste
325  REAL, DIMENSION(klon) :: qincloud
326  REAL, DIMENSION(klon) :: zrfl, zifl
327  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl
328  REAL, DIMENSION(klon) :: zoliql, zoliqi
329  REAL, DIMENSION(klon) :: zt, zp
330  REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zsnow
331  REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv
332  REAL, DIMENSION(klon) :: qtot, zeroklon
333  ! Variables precipitation energy conservation
334  REAL, DIMENSION(klon) :: zmqc
335  REAL :: zalpha_tr
336  REAL :: zfrac_lessi
337  REAL, DIMENSION(klon) :: zprec_cond
338  REAL, DIMENSION(klon) :: zlh_solid
339  REAL, DIMENSION(klon) :: ztupnew
340  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
341  REAL, DIMENSION(klon) :: zrflclr, zrflcld
342  REAL, DIMENSION(klon) :: ziflclr, ziflcld
343  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
344  REAL, DIMENSION(klon) :: tot_zneb
345  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz
346  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb
347  ! for ice sedimentation
348  REAL, DIMENSION(klon) :: dzsed, flsed, cfsed
349  REAL, DIMENSION(klon) :: dzsed_abv, flsed_abv, cfsed_abv
350  REAL :: qice_sedim
351
352  ! for quantity of condensates seen by radiation
353  REAL, DIMENSION(klon) :: zradocond, zradoice
354  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
355 
356  ! for condensation and ice supersaturation
357  REAL, DIMENSION(klon) :: qvc, qvcl, shear
358  REAL :: delta_z, deepconv_coef
359  ! for contrails
360  REAL, DIMENSION(klon) :: lincontfra, circontfra, qlincont, qcircont
361  REAL, DIMENSION(klon) :: totfra_in, qtot_in
362  LOGICAL, DIMENSION(klon) :: pt_pron_clds
363  REAL, DIMENSION(klon) :: dzsed_lincont, flsed_lincont, cfsed_lincont
364  REAL, DIMENSION(klon) :: dzsed_circont, flsed_circont, cfsed_circont
365  REAL, DIMENSION(klon) :: dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv
366  REAL, DIMENSION(klon) :: dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv
367  REAL :: qice_cont
368  !--for Lamquin et al 2012 diagnostics
369  REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
370  REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP
371
372  INTEGER i, k, kk, iter
373  INTEGER, DIMENSION(klon) :: n_i
374  INTEGER ncoreczq
375  LOGICAL iftop
376
377  LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb
378  LOGICAL, DIMENSION(klon) :: keepgoing
379
380  CHARACTER (len = 20) :: modname = 'lscp'
381  CHARACTER (len = 80) :: abort_message
382 
383
384!===============================================================================
385! INITIALISATION
386!===============================================================================
387
388! Few initial checks
389
390
391IF (iflag_fisrtilp_qsat .LT. 0) THEN
392    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
393    CALL abort_physic(modname,abort_message,1)
394ENDIF
395
396! AA for 'safety' reasons
397zalpha_tr   = 0.
398zfrac_lessi = 0.
399beta(:,:)= 0.
400
401! Initialisation of variables:
402
403prfl(:,:) = 0.0
404psfl(:,:) = 0.0
405d_t(:,:) = 0.0
406d_q(:,:) = 0.0
407d_ql(:,:) = 0.0
408d_qi(:,:) = 0.0
409rneb(:,:) = 0.0
410rnebth(:,:)=0.0
411pfraclr(:,:)=0.0
412pfracld(:,:)=0.0
413cldfraliq(:,:)=0.
414sigma2_icefracturb(:,:)=0.
415mean_icefracturb(:,:)=0.
416cldfraliqth(:,:)=0.
417sigma2_icefracturbth(:,:)=0.
418mean_icefracturbth(:,:)=0.
419radocond(:,:) = 0.0
420radicefrac(:,:) = 0.0
421frac_nucl(:,:) = 1.0
422frac_impa(:,:) = 1.0
423rain(:) = 0.0
424snow(:) = 0.0
425zrfl(:) = 0.0
426zifl(:) = 0.0
427zneb(:) = seuil_neb
428zrflclr(:) = 0.0
429ziflclr(:) = 0.0
430zrflcld(:) = 0.0
431ziflcld(:) = 0.0
432tot_zneb(:) = 0.0
433zeroklon(:) = 0.0
434zdistcltop(:)=0.0
435ztemp_cltop(:) = 0.0
436ztupnew(:)=0.0
437ctot_vol(:,:)=0.0
438rneblsvol(:,:)=0.0
439znebprecip(:)=0.0
440znebprecipclr(:)=0.0
441znebprecipcld(:)=0.0
442distcltop(:,:)=0.
443temp_cltop(:,:)=0.
444
445
446!--Ice supersaturation
447gamma_cond(:,:) = 1.
448qissr(:,:)      = 0.
449issrfra(:,:)    = 0.
450dcf_sub(:,:)    = 0.
451dcf_con(:,:)    = 0.
452dcf_mix(:,:)    = 0.
453dcfsed(:,:)     = 0.
454dqi_adj(:,:)    = 0.
455dqi_sub(:,:)    = 0.
456dqi_con(:,:)    = 0.
457dqi_mix(:,:)    = 0.
458dqised(:,:)     = 0.
459dqvc_adj(:,:)   = 0.
460dqvc_sub(:,:)   = 0.
461dqvc_con(:,:)   = 0.
462dqvc_mix(:,:)   = 0.
463dqvcsed(:,:)    = 0.
464qvc(:)          = 0.
465shear(:)        = 0.
466flsed(:)        = 0.
467flsed_lincont(:)= 0.
468flsed_circont(:)= 0.
469pt_pron_clds(:) = .FALSE.
470
471!--for Lamquin et al (2012) diagnostics
472issrfra100to150(:)   = 0.
473issrfra100to150UP(:) = 0.
474issrfra150to200(:)   = 0.
475issrfra150to200UP(:) = 0.
476issrfra200to250(:)   = 0.
477issrfra200to250UP(:) = 0.
478issrfra250to300(:)   = 0.
479issrfra250to300UP(:) = 0.
480issrfra300to400(:)   = 0.
481issrfra300to400UP(:) = 0.
482issrfra400to500(:)   = 0.
483issrfra400to500UP(:) = 0.
484
485!-- poprecip
486qraindiag(:,:)= 0.
487qsnowdiag(:,:)= 0.
488dqreva(:,:)   = 0.
489dqrauto(:,:)  = 0.
490dqrmelt(:,:)  = 0.
491dqrfreez(:,:) = 0.
492dqrcol(:,:)   = 0.
493dqssub(:,:)   = 0.
494dqsauto(:,:)  = 0.
495dqsrim(:,:)   = 0.
496dqsagg(:,:)   = 0.
497dqsfreez(:,:) = 0.
498dqsmelt(:,:)  = 0.
499zqupnew(:)    = 0.
500zqvapclr(:)   = 0.
501
502!-- cloud phase useful variables
503wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:))
504zqliq(:)=0.
505zqice(:)=0.
506zqvapcl(:)=0.
507zqliqth(:)=0.
508zqiceth(:)=0.
509zqvapclth(:)=0.
510sursat_e(:)=0.
511invtau_e(:)=0.
512pticefracturb(:)=.FALSE.
513
514
515!===============================================================================
516!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
517!===============================================================================
518
519ncoreczq=0
520
521DO k = klev, 1, -1
522
523    IF (k.LE.klev-1) THEN
524       iftop=.false.
525    ELSE
526       iftop=.true.
527    ENDIF
528
529    ! Initialisation temperature and specific humidity
530    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
531    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
532    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
533    ! d_t = temperature tendency due to lscp
534    ! The temperature of the overlying layer is updated here because needed for thermalization
535    DO i = 1, klon
536        zt(i)=temp(i,k)
537        zq(i)=qt(i,k)
538        zp(i)=pplay(i,k)
539        qliq_in(i) = ql_seri(i,k)
540        qice_in(i) = qi_seri(i,k)
541        zcf(i) = 0.
542        zfice(i) = 1.0   ! initialized at 1 as by default we assume mpc to be at ice saturation
543        dzfice(i) = 0.0
544        zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
545        dzficeth(i) = 0.0
546        zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
547        dzficeenv(i) = 0.0
548       
549
550        IF (.not. iftop) THEN
551           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
552           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
553           !--zqs(i) is the saturation specific humidity in the layer above
554           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
555        ENDIF
556        !c_iso init of iso
557    ENDDO
558    IF ( ok_ice_supersat ) THEN
559      cldfra_in(:) = cf_seri(:,k)
560      qvc_in(:) = qvc_seri(:,k)
561    ENDIF
562
563    ! --------------------------------------------------------------------
564    ! P1> Precipitation processes, before cloud formation:
565    !     Thermalization of precipitation falling from the overlying layer AND
566    !     Precipitation evaporation/sublimation/melting
567    !---------------------------------------------------------------------
568   
569    !================================================================
570    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
571    IF ( ok_poprecip ) THEN
572
573      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
574                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
575                        zqvapclr, zqupnew, flsed, flsed_lincont, flsed_circont, &
576                        cldfra_in, qvc_in, qliq_in, qice_in, &
577                        zrfl, zrflclr, zrflcld, &
578                        zifl, ziflclr, ziflcld, &
579                        dqreva(:,k), dqssub(:,k) &
580                        )
581
582    !================================================================
583    ELSE
584
585      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
586                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
587                        flsed, flsed_lincont, flsed_circont, &
588                        zrfl, zrflclr, zrflcld, &
589                        zifl, ziflclr, ziflcld, &
590                        dqreva(:,k), dqssub(:,k) &
591                        )
592
593    ENDIF ! (ok_poprecip)
594   
595    ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter
596    !-------------------------------------------------------
597
598     qtot(:)=zq(:)+zmqc(:)
599     CALL calc_qsat_ecmwf(klon,zt,qtot,zp,RTT,0,.false.,zqs,zdqs)
600     DO i = 1, klon
601            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
602            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
603            IF (zq(i) .LT. 1.e-15) THEN
604                ncoreczq=ncoreczq+1
605                zq(i)=1.e-15
606            ENDIF
607            ! c_iso: do something similar for isotopes
608
609     ENDDO
610   
611    ! -------------------------------------------------------------------------
612    ! P2> Cloud formation including condensation and cloud phase determination
613    !--------------------------------------------------------------------------
614    !
615    ! We always assume a 'fractional cloud' approach
616    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
617    ! is prescribed and depends on large scale variables and boundary layer
618    ! properties)
619    ! The decrease in condensed part due tu latent heating is taken into
620    ! account
621    ! -------------------------------------------------------------------
622
623        ! P2.1> With the PDFs (log-normal, bigaussian)
624        ! cloud properties calculation with the initial values of t and q
625        ! ----------------------------------------------------------------
626
627        ! initialise gammasat and stratiform_or_all_distrib
628        stratiform_or_all_distrib(:)=.TRUE.
629        gammasat(:)=1.
630
631        ! part of the code that is supposed to become obsolete soon
632        !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
633        IF (.NOT. ok_lscp_mergecond) THEN
634          IF (iflag_cld_th.GE.5) THEN
635               ! Cloud cover and content in meshes affected by shallow convection,
636               ! are retrieved from a bi-gaussian distribution of the saturation deficit
637               ! following Jam et al. 2013
638
639               IF (iflag_cloudth_vert.LE.2) THEN
640                  ! Old version of Arnaud Jam
641
642                    CALL cloudth(klon,klev,k,tv,                  &
643                         zq,qta,fraca,                            &
644                         qincloud,ctot,pspsk,paprs,pplay,tla,thl, &
645                         ratqs,zqs,temp,                              &
646                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
647
648
649                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
650                   ! Default version of Arnaud Jam
651                       
652                    CALL cloudth_v3(klon,klev,k,tv,                        &
653                         zq,qta,fraca,                                     &
654                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
655                         ratqs,sigma_qtherm,zqs,temp, &
656                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
657
658
659                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
660                   ! Jean Jouhaud's version, with specific separation between surface and volume
661                   ! cloud fraction Decembre 2018
662
663                    CALL cloudth_v6(klon,klev,k,tv,                        &
664                         zq,qta,fraca,                                     &
665                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
666                         ratqs,zqs,temp, &
667                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
668
669                ENDIF
670
671
672                DO i=1,klon
673                    rneb(i,k)=ctot(i,k)
674                    ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
675                    rneblsvol(i,k)=ctot_vol(i,k)
676                    zqn(i)=qincloud(i)
677                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
678                    qvc(i) = rneb(i,k) * zqs(i)
679                ENDDO
680
681               ! Cloud phase final determination for clouds after "old" cloudth calls
682               ! for those clouds, only temperature dependent phase partitioning (eventually modulated by
683               ! distance to cloud top) is available
684               ! compute distance to cloud top when cloud phase is computed depending on temperature
685               ! and distance to cloud top
686               IF (iflag_t_glace .GE. 4) THEN
687                    CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
688               ENDIF
689               CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
690
691          ENDIF
692
693          IF (iflag_cld_th .EQ. 5) THEN
694               ! When iflag_cld_th=5, we always assume
695                ! bi-gaussian distribution
696            stratiform_or_all_distrib(:) = .FALSE.
697
698          ELSEIF (iflag_cld_th .GE. 6) THEN
699                ! stratiform distribution only when no thermals
700            stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld
701
702          ENDIF
703
704        ENDIF ! .not. ok_lscp_mergecond
705        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
706
707
708        IF ( ok_ice_supersat ) THEN
709
710          !--Initialisation
711          IF ( ok_plane_contrail ) THEN
712            IF ( iftop ) THEN
713              dzsed_lincont_abv(:) = 0.
714              flsed_lincont_abv(:) = 0.
715              cfsed_lincont_abv(:) = 0.
716              dzsed_circont_abv(:) = 0.
717              flsed_circont_abv(:) = 0.
718              cfsed_circont_abv(:) = 0.
719            ELSE
720              dzsed_lincont_abv(:) = dzsed_lincont(:)
721              flsed_lincont_abv(:) = flsed_lincont(:)
722              cfsed_lincont_abv(:) = cfsed_lincont(:)
723              dzsed_circont_abv(:) = dzsed_circont(:)
724              flsed_circont_abv(:) = flsed_circont(:)
725              cfsed_circont_abv(:) = cfsed_circont(:)
726            ENDIF
727            dzsed_lincont(:) = 0.
728            flsed_lincont(:) = 0.
729            cfsed_lincont(:) = 0.
730            dzsed_circont(:) = 0.
731            flsed_circont(:) = 0.
732            cfsed_circont(:) = 0.
733            lincontfra(:)    = 0.
734            circontfra(:)    = 0.
735            qlincont(:)      = 0.
736            qcircont(:)      = 0.
737          ENDIF
738
739          IF ( iftop ) THEN
740            dzsed_abv(:) = 0.
741            flsed_abv(:) = 0.
742            cfsed_abv(:) = 0.
743          ELSE
744            dzsed_abv(:) = dzsed(:)
745            flsed_abv(:) = flsed(:)
746            cfsed_abv(:) = cfsed(:)
747          ENDIF
748          dzsed(:) = 0.
749          flsed(:) = 0.
750          cfsed(:) = 0.
751
752          DO i = 1, klon
753            pt_pron_clds(i) = ( cfcon(i,k) .LT. ( 1. - eps ) )
754          ENDDO
755          IF ( .NOT. ok_weibull_warm_clouds ) THEN
756            DO i = 1, klon
757              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( zt(i) .LE. temp_nowater )
758            ENDDO
759          ENDIF
760          IF ( ok_no_issr_strato ) THEN
761            DO i = 1, klon
762              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( stratomask(i,k) .EQ. 0. )
763            ENDDO
764          ENDIF
765
766          totfra_in(:) = 1.
767          qtot_in(:) = zq(:)
768
769          IF ( ok_nodeep_lscp ) THEN
770            DO i = 1, klon
771              !--If deep convection is activated, the condensation scheme activates
772              !--only in the environment. NB. the clear sky fraction will the be
773              !--maximised by 1. - cfcon(i,k)
774              IF ( pt_pron_clds(i) .AND. ptconv(i,k) ) THEN
775                totfra_in(i) = 1. - cfcon(i,k)
776                qtot_in(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
777              ENDIF
778            ENDDO
779          ENDIF
780
781          DO i = 1, klon
782            IF ( pt_pron_clds(i) ) THEN
783              IF ( cfcon(i,k) .LT. cfcon_old(i,k) ) THEN
784                !--If deep convection is weakening, we add the clouds that are not anymore
785                !--'in' deep convection to the advected clouds
786                cldfra_in(i) = cldfra_in(i) + ( cfcon_old(i,k) - cfcon(i,k) )
787                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
788                qice_in(i) = qice_in(i) + qccon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
789              ELSE
790                !--Else if deep convection is strengthening, it consumes the existing cloud
791                !--fraction (which does not at this moment represent deep convection)
792                deepconv_coef = 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) )
793                cldfra_in(i) = cldfra_in(i) * deepconv_coef
794                qvc_in(i)    = qvc_in(i)    * deepconv_coef
795                qice_in(i)   = qice_in(i)   * deepconv_coef
796                IF ( ok_plane_contrail ) THEN
797                  !--If contrails are activated, their fraction is also reduced when deep
798                  !--convection is active
799                  cfl_seri(i,k) = cfl_seri(i,k) * deepconv_coef
800                  qtl_seri(i,k) = qtl_seri(i,k) * deepconv_coef
801                  cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef
802                  qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef
803                ENDIF
804              ENDIF
805
806              !--Calculate the shear value (input for condensation and ice supersat)
807              !--Cell thickness [m]
808              delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
809              IF ( iftop ) THEN
810                ! top
811                shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
812                               + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
813              ELSEIF ( k .EQ. 1 ) THEN
814                ! surface
815                shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
816                               + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
817              ELSE
818                ! other layers
819                shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
820                                   - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
821                               + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
822                                   - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
823              ENDIF
824            ENDIF
825          ENDDO
826        ENDIF
827
828
829        DT(:) = 0.
830        n_i(:)=0
831        Tbef(:)=zt(:)
832        qlibef(:)=0.
833        Tbefth(:)=tla(:,k)*pspsk(:,k)
834        IF (k .GT. 1) THEN
835         Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1)
836        ELSE
837         Tbefthm1(:)=Tbefth(:)       
838        ENDIF
839        zqth(:)=qta(:,k)
840        zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:)
841
842        ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth
843        ! in case of ok_lscp_mergecond)
844        ! We iterate here to take into account the change in qsat(T) and ice fraction
845        ! during the condensation process
846        ! the increment in temperature is calculated using the first principle of
847        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
848        ! and a clear sky fraction)
849        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
850
851        DO iter=1,iflag_fisrtilp_qsat+1
852
853                keepgoing(:) = .FALSE.
854
855                DO i=1,klon
856
857                    ! keepgoing = .true. while convergence is not satisfied
858
859                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN
860
861                        ! if not convergence:
862                        ! we calculate a new iteration
863                        keepgoing(i) = .TRUE.
864
865                        ! P2.2.1> cloud fraction and condensed water mass calculation
866                        ! Calculated variables:
867                        ! rneb : cloud fraction
868                        ! zqn : total water within the cloud
869                        ! zcond: mean condensed water within the mesh
870                        ! rhcl: clear-sky relative humidity
871                        !---------------------------------------------------------------
872
873                        ! new temperature that only serves in the iteration process:
874                        Tbef(i)=Tbef(i)+DT(i)
875
876                        ! total water including mass of precip
877                        qtot(i)=zq(i)+zmqc(i)
878
879                      ENDIF
880
881                  ENDDO
882
883                  ! Calculation of saturation specific humidity
884                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs)
885                  ! also in thermals
886                  CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth)
887
888                  IF (iflag_icefrac .GE. 1) THEN
889                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
890                  ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD
891                  ! we make this option works only for the physically-based and tke-dependent param from Raillard et al.
892                  ! (iflag_icefrac>=1)
893                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl)
894                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi)
895                      DO i=1,klon
896                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
897                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
898                      ENDDO
899                  ENDIF
900                  IF (iflag_icefrac .GE. 2) THEN
901                  ! same adjustment for thermals
902                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth)
903                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith)
904                      DO i=1,klon
905                           zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i)
906                           zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i)
907                      ENDDO
908                  ENDIF
909
910                  CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt)
911                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
912                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
913
914                  ! Cloud condensation based on subgrid pdf
915                  !--AB Activates a condensation scheme that allows for
916                  !--ice supersaturation and contrails evolution from aviation
917                  IF (ok_ice_supersat) THEN
918
919                    !---------------------------------------------
920                    !--   CONDENSATION AND ICE SUPERSATURATION  --
921                    !---------------------------------------------
922
923                    CALL condensation_ice_supersat( &
924                        klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
925                        totfra_in, cldfra_in, qvc_in, qliq_in, qice_in, &
926                        shear, tke_dissip(:,k), cell_area, Tbef, qtot_in, zqs, &
927                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
928                        dzsed_abv, flsed_abv, cfsed_abv, &
929                        dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, &
930                        dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, &
931                        dzsed, flsed, cfsed, dzsed_lincont, flsed_lincont, cfsed_lincont, &
932                        dzsed_circont, flsed_circont, cfsed_circont, &
933                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
934                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
935                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), dqised(:,k), &
936                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), dqvcsed(:,k), &
937                        cfl_seri(:,k), cfc_seri(:,k), qtl_seri(:,k), qtc_seri(:,k), &
938                        flight_dist(:,k), flight_h2o(:,k), &
939                        lincontfra, circontfra, qlincont, qcircont, &
940                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
941                        dcfl_ini(:,k), dqil_ini(:,k), dqtl_ini(:,k), &
942                        dcfl_sub(:,k), dqil_sub(:,k), dqtl_sub(:,k), &
943                        dcfl_cir(:,k), dqtl_cir(:,k), &
944                        dcfl_mix(:,k), dqil_mix(:,k), dqtl_mix(:,k), &
945                        dcfl_sed(:,k), dqil_sed(:,k), dqtl_sed(:,k), &
946                        dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), &
947                        dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), &
948                        dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k))
949
950                    IF ( ok_nodeep_lscp ) THEN
951                      DO i = 1, klon
952                        !--If prognostic clouds are activated, deep convection vapor is
953                        !--re-added to the total water vapor
954                        IF ( keepgoing(i) .AND. ptconv(i,k) .AND. pt_pron_clds(i) ) THEN
955                          IF ( ( rneb(i,k) + cfcon(i,k) ) .GT. eps ) THEN
956                            zqn(i) = ( zqn(i) * rneb(i,k) &
957                                + ( qccon(i,k) + qvcon(i,k) ) * cfcon(i,k) ) &
958                                / ( rneb(i,k) + cfcon(i,k) )
959                          ELSE
960                            zqn(i) = 0.
961                          ENDIF
962                          rneb(i,k) = rneb(i,k) + cfcon(i,k)
963                          qvc(i) = qvc(i) + qvcon(i,k) * cfcon(i,k)
964                        ENDIF
965                      ENDDO
966                    ENDIF
967
968                  ELSE
969                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
970
971                  CALL condensation_lognormal( &
972                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
973                       keepgoing, rneb(:,k), zqn, qvc)
974
975
976                  ENDIF ! .NOT. ok_ice_supersat
977
978                  ! Volume cloud fraction
979                  rneblsvol(:,k)=rneb(:,k)
980
981
982                  IF  (ok_lscp_mergecond) THEN
983                      ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only
984                      CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), &
985                                     pspsk(:,k),zp,tla(:,k), &
986                                     ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), &
987                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
988                  ENDIF
989
990
991
992                  ! Cloud phase determination
993
994
995                  IF (iflag_icefrac .LE. 1) THEN
996                     ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere
997                     IF (iflag_t_glace .GE. 4) THEN
998                        CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
999                     ENDIF
1000                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
1001                  ENDIF
1002
1003                  IF (iflag_icefrac .EQ. 1) THEN
1004                     ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics
1005                     ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025)
1006                     ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds
1007                     pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
1008                     DO i=1,klon
1009                        qincloud(i)=zqn(i)
1010                        zcf(i)=MIN(MAX(rneb(i,k), 0.),1.)
1011                        sursat_e(i) = 0.
1012                        invtau_e(i) = 0.
1013                     ENDDO
1014                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in,    &
1015                     ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
1016                     zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k))                     
1017                     DO i=1,klon
1018                        IF (pticefracturb(i)) THEN
1019                          zfice(i)=zficeenv(i)
1020                          dzfice(i)=dzficeenv(i)
1021                        ENDIF
1022                     ENDDO
1023                   
1024
1025                  ELSEIF (iflag_icefrac .GE. 2) THEN
1026                     ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption)
1027                     ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment
1028                     ! hence we assume that no snow falls in thermals.
1029                     ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T
1030                     pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
1031                     !Thermals
1032                     DO i=1,klon
1033                        IF (fraca(i,k) .GT. min_frac_th_cld) THEN
1034                           zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)/fraca(i,k)
1035                           qincloud(i)=zqn(i)*fraca(i,k)
1036                        ELSE
1037                           zcf(i) = 0.
1038                           qincloud(i) = 0.
1039                        ENDIF
1040                        sursat_e(i)=cloudth_senv(i,k)/zqsi(i)
1041                        invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
1042                     ENDDO
1043                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), qice_in,  &
1044                     zeroklon, znebprecipcld, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, &
1045                     zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k))
1046                     !Environment
1047                     DO i=1,klon
1048                        qincloud(i)=zqn(i)*(1.-fraca(i,k))
1049                        zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)/(1.-fraca(i,k))
1050                        IF (k .GT. 1) THEN
1051                           ! evaluate the mixing sursaturation using saturation deficit at level below
1052                           ! as air pacels detraining into clouds have not (less) seen yet entrainement from above
1053                           sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i)))
1054                           ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho
1055                           invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
1056                        ELSE
1057                           sursat_e(i)=0.
1058                           invtau_e(i)=0.
1059                        ENDIF
1060                     ENDDO
1061                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in,    &
1062                     ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
1063                     zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
1064 
1065                    ! adjust zfice to account for condensates in thermals'fraction
1066                     DO i=1,klon
1067                        IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN
1068                          zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i)))
1069                          dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb
1070                        ENDIF
1071                     ENDDO
1072
1073                  ENDIF ! iflag_icefrac
1074
1075
1076
1077
1078                  DO i=1,klon
1079                      IF (keepgoing(i)) THEN
1080
1081                       ! P2.2.2> Approximative calculation of temperature variation DT
1082                       !           due to condensation.
1083                       ! Calculated variables:
1084                       ! dT : temperature change due to condensation
1085                       !---------------------------------------------------------------
1086
1087               
1088                        IF (zfice(i).LT.1) THEN
1089                            cste=RLVTT
1090                        ELSE
1091                            cste=RLSTT
1092                        ENDIF
1093                       
1094                        IF ( ok_unadjusted_clouds ) THEN
1095                          !--AB We relax the saturation adjustment assumption
1096                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1097                          IF ( rneb(i,k) .GT. eps ) THEN
1098                            qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
1099                          ELSE
1100                            qlibef(i) = 0.
1101                          ENDIF
1102                        ELSE
1103                          qlibef(i)=max(0.,zqn(i)-zqs(i))
1104                        ENDIF
1105
1106                        IF ( ok_ice_sedim ) THEN
1107                          qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
1108                              / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
1109                          ! Add the ice that was sedimented, as it is not included in zqn
1110                          qlibef(i) = qlibef(i) + qice_sedim
1111                        ENDIF
1112
1113                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1114                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(i)
1115                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1116                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
1117                              *qlibef(i)*dzfice(i)
1118                        ! here we update a provisory temperature variable that only serves in the iteration
1119                        ! process
1120                        DT(i)=num/denom
1121                        n_i(i)=n_i(i)+1
1122
1123                    ENDIF ! end keepgoing
1124
1125                ENDDO     ! end loop on i
1126
1127        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
1128
1129        ! P2.2> Final quantities calculation
1130        ! Calculated variables:
1131        !   rneb : cloud fraction
1132        !   zcond: mean condensed water in the mesh
1133        !   zqn  : mean water vapor in the mesh
1134        !   zfice: ice fraction in clouds
1135        !   zt   : temperature
1136        !   rhcl : clear-sky relative humidity
1137        ! ----------------------------------------------------------------
1138
1139
1140        ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type
1141        !---------------------------------------------------------------------------------------------
1142        DO i=1, klon
1143                ! Checks on rneb, rhcl and zqn
1144                IF (rneb(i,k) .LE. 0.0) THEN
1145                    zqn(i) = 0.0
1146                    rneb(i,k) = 0.0
1147                    zcond(i) = 0.0
1148                    rhcl(i,k)=zq(i)/zqs(i)
1149                ELSE IF (rneb(i,k) .GE. 1.0) THEN
1150                    zqn(i) = zq(i)
1151                    rneb(i,k) = 1.0
1152                    IF ( ok_unadjusted_clouds ) THEN
1153                      !--AB We relax the saturation adjustment assumption
1154                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1155                      zcond(i) = MAX(0., zqn(i) - qvc(i))
1156                    ELSE
1157                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1158                    ENDIF
1159                    rhcl(i,k)=1.0
1160                ELSE
1161                    IF ( ok_unadjusted_clouds ) THEN
1162                      !--AB We relax the saturation adjustment assumption
1163     
1164                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1165                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
1166                    ELSE
1167                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1168                    ENDIF
1169
1170                    ! following line is very strange and probably wrong
1171                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
1172                    ! Correct calculation of clear-sky relative humidity. To activate once
1173                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
1174                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
1175                    !--This is slighly approximated, the actual formula is
1176                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
1177                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
1178                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
1179                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
1180                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
1181
1182                ENDIF
1183
1184                ! water vapor update
1185                zq(i) = zq(i) - zcond(i)
1186
1187                IF ( ok_ice_sedim ) THEN
1188                  qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
1189                      / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
1190                  ! Remove the ice that was sedimented. As it is not included in zqn,
1191                  ! we only remove it from the total water
1192                  zq(i) = zq(i) - qice_sedim
1193                  ! Temperature update due to phase change (sedimented ice was condensed)
1194                  zt(i) = zt(i) + qice_sedim &
1195                        * RLSTT / RCPD / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) + zcond(i) ) )
1196                ENDIF
1197                       
1198                ! temperature update due to phase change
1199                zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1200                &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1201                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1202        ENDDO
1203
1204    ! ----------------------------------------------------------------
1205    ! P3> Precipitation processes, after cloud formation
1206    !     - precipitation formation, melting/freezing
1207    ! ----------------------------------------------------------------
1208
1209    ! Initiate the quantity of liquid and solid condensates
1210    ! Note that in the following, zcond is the total amount of condensates
1211    ! including newly formed precipitations (i.e., condensates formed by the
1212    ! condensation process), while zoliq is the total amount of condensates in
1213    ! the cloud (i.e., on which precipitation processes have applied)
1214    DO i = 1, klon
1215      zoliq(i)  = zcond(i)
1216      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1217      zoliqi(i) = zcond(i) * zfice(i)
1218    ENDDO
1219
1220    IF (ok_plane_contrail) THEN
1221
1222      !--Ice water content of contrails
1223      qice_lincont(:,k) = qlincont(:) - zqs(:) * lincontfra(:)
1224      qice_circont(:,k) = qcircont(:) - zqs(:) * circontfra(:)
1225
1226      !--Contrails precipitate as natural clouds. We save the partition of ice
1227      !--between natural clouds and contrails
1228      !--NB. we use qlincont / qcircont as a temporary variable to save this partition
1229      IF ( ok_precip_contrails ) THEN
1230        DO i = 1, klon
1231          IF ( zoliqi(i) .GT. 0. ) THEN
1232            qlincont(i) = qice_lincont(i,k) / zoliqi(i)
1233            qcircont(i) = qice_circont(i,k) / zoliqi(i)
1234          ELSE
1235            qlincont(i) = 0.
1236            qcircont(i) = 0.
1237          ENDIF
1238        ENDDO
1239      ELSE
1240        !--If linear contrails do not precipitate, they are removed temporarily from
1241        !--the cloud variables
1242        DO i = 1, klon
1243          qice_cont = qice_lincont(i,k) + qice_circont(i,k)
1244          rneb(i,k) = rneb(i,k) - ( lincontfra(i) + circontfra(i) )
1245          zoliq(i) = zoliq(i) - qice_cont
1246          zoliqi(i) = zoliqi(i) - qice_cont
1247        ENDDO
1248      ENDIF
1249    ENDIF
1250
1251    !================================================================
1252    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1253    IF (ok_poprecip) THEN
1254
1255      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, &
1256                            ctot_vol(:,k), ptconv(:,k), &
1257                            zt, zq, zoliql, zoliqi, zfice, &
1258                            rneb(:,k), znebprecipclr, znebprecipcld, &
1259                            zrfl, zrflclr, zrflcld, &
1260                            zifl, ziflclr, ziflcld, &
1261                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1262                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1263                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1264                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
1265                            )
1266      DO i = 1, klon
1267          zoliq(i) = zoliql(i) + zoliqi(i)
1268      ENDDO
1269
1270    !================================================================
1271    ELSE
1272
1273      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
1274                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
1275                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
1276                            rneb(:,k), znebprecipclr, znebprecipcld, &
1277                            zneb, tot_zneb, zrho_up, zvelo_up, &
1278                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1279                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
1280                            )
1281
1282    ENDIF ! ok_poprecip
1283
1284    IF ( ok_plane_contrail ) THEN
1285      !--Contrails fraction is left unchanged, but contrails water has changed
1286      !--We alse compute the ice content that will be seen by radiation
1287      !--(qradice_lincont/circont)
1288      IF ( ok_precip_contrails ) THEN
1289        DO i = 1, klon
1290          IF ( zoliqi(i) .GT. 0. ) THEN
1291            qradice_lincont(i,k) = zradocond(i) * qlincont(i)
1292            qlincont(i) = zqs(i) * lincontfra(i) + zoliqi(i) * qlincont(i)
1293            qradice_circont(i,k) = zradocond(i) * qcircont(i)
1294            qcircont(i) = zqs(i) * circontfra(i) + zoliqi(i) * qcircont(i)
1295          ELSE
1296            qradice_lincont(i,k) = 0.
1297            lincontfra(i) = 0.
1298            qlincont(i) = 0.
1299            qradice_circont(i,k) = 0.
1300            circontfra(i) = 0.
1301            qcircont(i) = 0.
1302          ENDIF
1303        ENDDO
1304      ELSE
1305        !--If contrails do not precipitate, they are put back into
1306        !--the cloud variables
1307        DO i = 1, klon
1308          rneb(i,k) = rneb(i,k) + ( lincontfra(i) + circontfra(i) )
1309          qice_cont = qice_lincont(i,k) + qice_circont(i,k)
1310          zoliq(i) = zoliq(i) + qice_cont
1311          zoliqi(i) = zoliqi(i) + qice_cont
1312          zradocond(i) = zradocond(i) + qice_cont
1313          zradoice(i) = zradoice(i) + qice_cont
1314          qradice_lincont(i,k) = qice_lincont(i,k)
1315          qradice_circont(i,k) = qice_circont(i,k)
1316        ENDDO
1317      ENDIF
1318    ENDIF
1319
1320    ! End of precipitation processes after cloud formation
1321    ! ----------------------------------------------------
1322
1323    !----------------------------------------------------------------------
1324    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1325    !----------------------------------------------------------------------
1326
1327    DO i=1,klon
1328
1329      IF (ok_poprecip) THEN
1330        IF (ok_radocond_snow) THEN
1331          radocond(i,k) = zoliq(i)
1332          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1333        ELSE
1334          radocond(i,k) = zoliq(i)
1335          zradoice(i) = zoliqi(i)
1336        ENDIF
1337      ELSE
1338        radocond(i,k) = zradocond(i)
1339      ENDIF
1340
1341      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1342      ! by the radiation scheme
1343      IF (radocond(i,k) .GT. 0.) THEN
1344        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1345      ENDIF
1346    ENDDO
1347
1348    ! ----------------------------------------------------------------
1349    ! P5> Wet scavenging
1350    ! ----------------------------------------------------------------
1351
1352    !Scavenging through nucleation in the layer
1353
1354    DO i = 1,klon
1355       
1356        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1357            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1358        ELSE
1359            beta(i,k) = 0.
1360        ENDIF
1361
1362        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1363
1364        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1365
1366            IF (temp(i,k) .GE. temp_nowater) THEN
1367                zalpha_tr = a_tr_sca(3)
1368            ELSE
1369                zalpha_tr = a_tr_sca(4)
1370            ENDIF
1371
1372            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1373            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1374
1375            ! Nucleation with a factor of -1 instead of -0.5
1376            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1377
1378        ENDIF
1379
1380    ENDDO
1381
1382    ! Scavenging through impaction in the underlying layer
1383
1384    DO kk = k-1, 1, -1
1385
1386        DO i = 1, klon
1387
1388            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1389
1390                IF (temp(i,kk) .GE. temp_nowater) THEN
1391                    zalpha_tr = a_tr_sca(1)
1392                ELSE
1393                    zalpha_tr = a_tr_sca(2)
1394                ENDIF
1395
1396              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1397              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1398
1399            ENDIF
1400
1401        ENDDO
1402
1403    ENDDO
1404   
1405    !------------------------------------------------------------
1406    ! P6 > write diagnostics and outputs
1407    !------------------------------------------------------------
1408
1409    CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs)
1410    CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs)
1411   
1412    !--AB Write diagnostics and tracers for ice supersaturation
1413    IF ( ok_plane_contrail ) THEN
1414      DO i = 1, klon
1415        IF ( zoliq(i) .LE. 0. ) THEN
1416          lincontfra(i) = 0.
1417          circontfra(i) = 0.
1418          qlincont(i) = 0.
1419          qcircont(i) = 0.
1420        ENDIF
1421      ENDDO
1422      cfl_seri(:,k) = lincontfra(:)
1423      cfc_seri(:,k) = circontfra(:)
1424      qtl_seri(:,k) = qlincont(:)
1425      qtc_seri(:,k) = qcircont(:)
1426    ENDIF
1427
1428    IF ( ok_ice_supersat ) THEN
1429
1430      DO i = 1, klon
1431
1432        !--We save the cloud properties that will be advected
1433        cf_seri(i,k) = rneb(i,k)
1434        qvc_seri(i,k) = qvc(i)
1435
1436        !--We keep convective clouds properties in memory, and account for
1437        !--the sink of condensed water from precipitation
1438        IF ( ptconv(i,k) ) THEN
1439          IF ( zcond(i) .GT. 0. ) THEN
1440            qvcon_old(i,k) = qvcon(i,k)
1441            qccon_old(i,k) = qccon(i,k) * zoliq(i) / zcond(i)
1442          ELSE
1443            qvcon_old(i,k) = 0.
1444            qccon_old(i,k) = 0.
1445          ENDIF
1446        ELSE
1447          qvcon_old(i,k) = 0.
1448          qccon_old(i,k) = 0.
1449        ENDIF
1450
1451        !--Deep convection clouds properties are not advected
1452        IF ( ptconv(i,k) .AND. pt_pron_clds(i) .AND. ok_nodeep_lscp ) THEN
1453          cf_seri(i,k) = MAX(0., cf_seri(i,k) - cfcon(i,k))
1454          qvc_seri(i,k) = MAX(0., qvc_seri(i,k) - qvcon_old(i,k) * cfcon(i,k))
1455          zoliq(i) = MAX(0., zoliq(i) - qccon_old(i,k) * cfcon(i,k))
1456          zoliqi(i) = MAX(0., zoliqi(i) - qccon_old(i,k) * cfcon(i,k))
1457        ENDIF
1458        !--Deep convection clouds properties are removed from radiative properties
1459        !--outputed from lscp (NB. rneb and radocond are only used for the radiative
1460        !--properties and are NOT prognostics)
1461        !--We must have iflag_coupl == 5 for this coupling to work
1462        IF ( ptconv(i,k) .AND. pt_pron_clds(i) .AND. ok_nodeep_lscp_rad ) THEN
1463          rneb(i,k) = MAX(0., rneb(i,k) - cfcon(i,k))
1464          radocond(i,k) = MAX(0., radocond(i,k) - qccon_old(i,k) * cfcon(i,k))
1465        ENDIF
1466
1467        !--If everything was precipitated, the remaining empty cloud is dissipated
1468        !--and everything is transfered to the subsaturated clear sky region
1469        !--NB. we do not change rneb, as it is a diagnostic only
1470        IF ( zoliq(i) .LE. 0. ) THEN
1471          cf_seri(i,k) = 0.
1472          qvc_seri(i,k) = 0.
1473          qvc(i) = 0.
1474        ENDIF
1475
1476        !--Diagnostics
1477        gamma_cond(i,k) = gammasat(i)
1478        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1479        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1480        qcld(i,k) = qvc(i) + zoliq(i)
1481
1482        !--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
1483        !--methodology: in each layer, we make a maximum random overlap assumption for
1484        !--ice supersaturation
1485        IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
1486                IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
1487                        issrfra100to150(i) = 1.
1488                ELSE
1489                        issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
1490                                ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
1491                              / ( 1. - issrfra100to150UP(i) )
1492                        issrfra100to150UP(i) = issrfra(i,k)
1493                ENDIF
1494        ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
1495                IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
1496                        issrfra150to200(i) = 1.
1497                ELSE
1498                        issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
1499                                ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
1500                              / ( 1. - issrfra150to200UP(i) )
1501                        issrfra150to200UP(i) = issrfra(i,k)
1502                ENDIF
1503        ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
1504                IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
1505                        issrfra200to250(i) = 1.
1506                ELSE
1507                        issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
1508                                ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
1509                              / ( 1. - issrfra200to250UP(i) )
1510                        issrfra200to250UP(i) = issrfra(i,k)
1511                ENDIF
1512        ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
1513                IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
1514                        issrfra250to300(i) = 1.
1515                ELSE
1516                        issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
1517                                ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
1518                              / ( 1. - issrfra250to300UP(i) )
1519                        issrfra250to300UP(i) = issrfra(i,k)
1520                ENDIF
1521        ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
1522                IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
1523                        issrfra300to400(i) = 1.
1524                ELSE
1525                        issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
1526                                ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
1527                              / ( 1. - issrfra300to400UP(i) )
1528                        issrfra300to400UP(i) = issrfra(i,k)
1529                ENDIF
1530        ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
1531                IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
1532                        issrfra400to500(i) = 1.
1533                ELSE
1534                        issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
1535                                ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
1536                              / ( 1. - issrfra400to500UP(i) )
1537                        issrfra400to500UP(i) = issrfra(i,k)
1538                ENDIF
1539        ENDIF
1540
1541      ENDDO
1542    ENDIF
1543
1544    ! Outputs:
1545    !-------------------------------
1546    ! Precipitation fluxes at layer interfaces
1547    ! + precipitation fractions +
1548    ! temperature and water species tendencies
1549    DO i = 1, klon
1550        psfl(i,k)=zifl(i)
1551        prfl(i,k)=zrfl(i)
1552        pfraclr(i,k)=znebprecipclr(i)
1553        pfracld(i,k)=znebprecipcld(i)
1554        distcltop(i,k)=zdistcltop(i)
1555        temp_cltop(i,k)=ztemp_cltop(i)
1556        d_q(i,k) = zq(i) - qt(i,k)
1557        d_t(i,k) = zt(i) - temp(i,k)
1558
1559        IF (ok_bug_phase_lscp) THEN
1560           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1561           d_qi(i,k) = zfice(i)*zoliq(i)
1562        ELSE
1563           d_ql(i,k) = zoliql(i)
1564           d_qi(i,k) = zoliqi(i)   
1565        ENDIF
1566
1567    ENDDO
1568
1569
1570ENDDO ! loop on k from top to bottom
1571
1572
1573  ! Rain or snow at the surface (depending on the first layer temperature)
1574  DO i = 1, klon
1575      snow(i) = zifl(i)
1576      rain(i) = zrfl(i)
1577      ! c_iso final output
1578  ENDDO
1579
1580  IF ( ok_ice_sedim ) THEN
1581    DO i = 1, klon
1582      snow(i) = snow(i) + flsed(i) + flsed_lincont(i) + flsed_circont(i)
1583    ENDDO
1584  ENDIF
1585
1586  IF (ncoreczq>0) THEN
1587      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1588  ENDIF
1589
1590
1591RETURN
1592
1593END SUBROUTINE lscp
1594!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1595
1596END MODULE lmdz_lscp_main
Note: See TracBrowser for help on using the repository browser.