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

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

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