source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90 @ 5647

Last change on this file since 5647 was 5647, checked in by evignon, 2 months ago

petite correction pour la phase mixte dans les thermiques, nouvelle param Lea Etienne
File size: 56.9 KB
RevLine 
[5614]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     ptconv, ratqs, sigma_qtherm,                       &
11     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
12     pfraclr, pfracld,                                  &
13     cldfraliq, cldfraliqth,                            &
14     sigma2_icefracturb,sigma2_icefracturbth,           &
15     mean_icefracturb,mean_icefracturbth,               &
16     radocond, radicefrac, rain, snow,                  &
17     frac_impa, frac_nucl, beta,                        &
18     prfl, psfl, rhcl, qta, fraca,                      &
19     tv, pspsk, tla, thl, wth, iflag_cld_th,            &
20     iflag_ice_thermo, distcltop, temp_cltop,           &
21     tke, tke_dissip,                                   &
22     entr_therm, detr_therm,                            &
23     cell_area,                                         &
24     cf_seri, rvc_seri, u_seri, v_seri,                 &
25     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
26     dcf_sub, dcf_con, dcf_mix,          &
27     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
28     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
29     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
30     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
31     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
32     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
33     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
34     dqsmelt, dqsfreez)
35
36!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
38!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
39!--------------------------------------------------------------------------------
40! Date: 01/2021
41!--------------------------------------------------------------------------------
42! Aim: Large Scale Clouds and Precipitation (LSCP)
43!
44! This code is a new version of the fisrtilp.F90 routine, which is itself a
45! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
46! and 'ilp' (il pleut, L. Li)
47!
48! Compared to the original fisrtilp code, lscp
49! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
50! -> consider always precipitation thermalisation (fl_cor_ebil>0)
51! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
52! -> option "oldbug" by JYG has been removed
53! -> iflag_t_glace >0 always
54! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
55! -> rectangular distribution from L. Li no longer available
56! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
57!--------------------------------------------------------------------------------
58! References:
59!
60! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
61! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
62! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
63! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
64! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
65! - Touzze-Peifert Ludo, PhD thesis, p117-124
66! -------------------------------------------------------------------------------
67! Code structure:
68!
69! P0>     Thermalization of the precipitation coming from the overlying layer
70! P1>     Evaporation of the precipitation (falling from the k+1 level)
71! P2>     Cloud formation (at the k level)
72! P2.A.1> With the PDFs, calculation of cloud properties using the inital
73!         values of T and Q
74! P2.A.2> Coupling between condensed water and temperature
75! P2.A.3> Calculation of final quantities associated with cloud formation
76! P2.B>   Release of Latent heat after cloud formation
77! P3>     Autoconversion to precipitation (k-level)
78! P4>     Wet scavenging
79!------------------------------------------------------------------------------
80! Some preliminary comments (JBM) :
81!
82! The cloud water that the radiation scheme sees is not the same that the cloud
83! water used in the physics and the dynamics
84!
85! During the autoconversion to precipitation (P3 step), radocond (cloud water used
86! by the radiation scheme) is calculated as an average of the water that remains
87! in the cloud during the precipitation and not the water remaining at the end
88! of the time step. The latter is used in the rest of the physics and advected
89! by the dynamics.
90!
91! In summary:
92!
93! Radiation:
94! xflwc(newmicro)+xfiwc(newmicro) =
95!   radocond=lwcon(nc)+iwcon(nc)
96!
97! Notetheless, be aware of:
98!
99! radocond .NE. ocond(nc)
100! i.e.:
101! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
102! but oliq+(ocond-oliq) .EQ. ocond
103! (which is not trivial)
104!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105!
106
107! USE de modules contenant des fonctions.
108USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
109USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
110USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
111USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth
112USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
113USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
114
115! Use du module lmdz_lscp_ini contenant les constantes
116USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
117USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
118USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
119USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
120USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater
121USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
122USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
123USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
124USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth
125
126IMPLICIT NONE
127
128!===============================================================================
129! VARIABLES DECLARATION
130!===============================================================================
131
132! INPUT VARIABLES:
133!-----------------
134
135  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
136  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
137 
138
139  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
140  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
141
142  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
143  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
144  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
145  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
146  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
147  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
148  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
149  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke             ! turbulent kinetic energy [m2/s2]
150  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip      ! TKE dissipation [m2/s3]
[5647]151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: entr_therm      ! thermal plume entrainment rate * dz [kg/s/m2]
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: detr_therm      ! thermal plume detrainment rate * dz [kg/s/m2]
[5614]153
154
155 
156  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
157
158  !Inputs associated with thermal plumes
159
160  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
161  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
162  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
163  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
164  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
165  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid potential temperature within thermals [K]
166  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: wth                 ! vertical velocity within thermals [m/s]
167  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: sigma_qtherm        ! controls the saturation deficit distribution width in thermals [-]
168
169
170
171  ! INPUT/OUTPUT variables
172  !------------------------
173 
174  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
175  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function that sets the large-scale water distribution width
176
177
178  ! INPUT/OUTPUT condensation and ice supersaturation
179  !--------------------------------------------------
180  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
181  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
182  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
183  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
184  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
185
186  ! INPUT/OUTPUT aviation
187  !--------------------------------------------------
188  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
189  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
190 
191  ! OUTPUT variables
192  !-----------------
193
194  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
195  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
196  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
197  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
201  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
202  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq        ! liquid fraction of cloud fraction [-]
203  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliqth      ! liquid fraction of cloud fraction in thermals [-]
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
206  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
207  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturbth   ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
208  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
211  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
212  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
213  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
214  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
218
219  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
220 
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
223 
224  ! for condensation and ice supersaturation
225
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
228  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
229  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
230  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
239  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
240  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
241  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
242  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
245
246  ! for contrails and aviation
247
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
256
257
258  ! for POPRECIP
259
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
263  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
269  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
270  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
271  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
272  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
273
274  ! for thermals
275
276  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
277  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
278  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
279  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
280
281
282  ! LOCAL VARIABLES:
283  !----------------
284  REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol
285  REAL, DIMENSION(klon,klev) :: wls                                 !-- large scalce vertical velocity [m/s]
286  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
287  REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith
288  REAL :: zdelta
289  REAL, DIMENSION(klon) :: zdqsdT_raw
290  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
[5633]291  REAL, DIMENSION(klon) :: Tbef,Tbefth,Tbefthm1,qlibef,DT                  ! temperature, humidity and temp. variation during condensation iteration
[5614]292  REAL :: num,denom
293  REAL :: cste
294  REAL, DIMENSION(klon) :: qincloud
295  REAL, DIMENSION(klon) :: zrfl, zifl
296  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl
297  REAL, DIMENSION(klon) :: zoliql, zoliqi
298  REAL, DIMENSION(klon) :: zt, zp
299  REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zqi_ini, zsnow
300  REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv
301  REAL, DIMENSION(klon) :: qtot, zeroklon
302  ! Variables precipitation energy conservation
303  REAL, DIMENSION(klon) :: zmqc
304  REAL :: zalpha_tr
305  REAL :: zfrac_lessi
306  REAL, DIMENSION(klon) :: zprec_cond
307  REAL, DIMENSION(klon) :: zlh_solid
308  REAL, DIMENSION(klon) :: ztupnew
309  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
310  REAL, DIMENSION(klon) :: zrflclr, zrflcld
311  REAL, DIMENSION(klon) :: ziflclr, ziflcld
312  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
313  REAL, DIMENSION(klon) :: tot_zneb
[5647]314  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz
[5614]315  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb
316
317  ! for quantity of condensates seen by radiation
318  REAL, DIMENSION(klon) :: zradocond, zradoice
319  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
320 
321  ! for condensation and ice supersaturation
322  REAL, DIMENSION(klon) :: qvc, qvcl, shear
323  REAL :: delta_z
324  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
325  ! Constants used for calculating ratios that are advected (using a parent-child
326  ! formalism). This is not done in the dynamical core because at this moment,
327  ! only isotopes can use this parent-child formalism. Note that the two constants
328  ! are the same as the one use in the dynamical core, being also defined in
329  ! dyn3d_common/infotrac.F90
330  REAL :: min_qParent, min_ratio
331
332  INTEGER i, k, kk, iter
333  INTEGER, DIMENSION(klon) :: n_i
334  INTEGER ncoreczq
335  LOGICAL iftop
336
337  LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb
338  LOGICAL, DIMENSION(klon) :: keepgoing
339
340  CHARACTER (len = 20) :: modname = 'lscp'
341  CHARACTER (len = 80) :: abort_message
342 
343
344!===============================================================================
345! INITIALISATION
346!===============================================================================
347
348! Few initial checks
349
350
351IF (iflag_fisrtilp_qsat .LT. 0) THEN
352    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
353    CALL abort_physic(modname,abort_message,1)
354ENDIF
355
356! AA for 'safety' reasons
357zalpha_tr   = 0.
358zfrac_lessi = 0.
359beta(:,:)= 0.
360
361! Initialisation of variables:
362
363prfl(:,:) = 0.0
364psfl(:,:) = 0.0
365d_t(:,:) = 0.0
366d_q(:,:) = 0.0
367d_ql(:,:) = 0.0
368d_qi(:,:) = 0.0
369rneb(:,:) = 0.0
370rnebth(:,:)=0.0
371pfraclr(:,:)=0.0
372pfracld(:,:)=0.0
373cldfraliq(:,:)=0.
374sigma2_icefracturb(:,:)=0.
375mean_icefracturb(:,:)=0.
376cldfraliqth(:,:)=0.
377sigma2_icefracturbth(:,:)=0.
378mean_icefracturbth(:,:)=0.
379radocond(:,:) = 0.0
380radicefrac(:,:) = 0.0
381frac_nucl(:,:) = 1.0
382frac_impa(:,:) = 1.0
383rain(:) = 0.0
384snow(:) = 0.0
385zrfl(:) = 0.0
386zifl(:) = 0.0
387zneb(:) = seuil_neb
388zrflclr(:) = 0.0
389ziflclr(:) = 0.0
390zrflcld(:) = 0.0
391ziflcld(:) = 0.0
392tot_zneb(:) = 0.0
393zeroklon(:) = 0.0
394zdistcltop(:)=0.0
395ztemp_cltop(:) = 0.0
396ztupnew(:)=0.0
397ctot_vol(:,:)=0.0
398rneblsvol(:,:)=0.0
399znebprecip(:)=0.0
400znebprecipclr(:)=0.0
401znebprecipcld(:)=0.0
402distcltop(:,:)=0.
403temp_cltop(:,:)=0.
404
405
406!--Ice supersaturation
407gamma_cond(:,:) = 1.
408qissr(:,:)      = 0.
409issrfra(:,:)    = 0.
410dcf_sub(:,:)    = 0.
411dcf_con(:,:)    = 0.
412dcf_mix(:,:)    = 0.
413dqi_adj(:,:)    = 0.
414dqi_sub(:,:)    = 0.
415dqi_con(:,:)    = 0.
416dqi_mix(:,:)    = 0.
417dqvc_adj(:,:)   = 0.
418dqvc_sub(:,:)   = 0.
419dqvc_con(:,:)   = 0.
420dqvc_mix(:,:)   = 0.
421fcontrN(:,:)    = 0.
422fcontrP(:,:)    = 0.
423Tcontr(:,:)     = missing_val
424qcontr(:,:)     = missing_val
425qcontr2(:,:)    = missing_val
426dcf_avi(:,:)    = 0.
427dqi_avi(:,:)    = 0.
428dqvc_avi(:,:)   = 0.
429qvc(:)          = 0.
430shear(:)        = 0.
431min_qParent     = 1.e-30
432min_ratio       = 1.e-16
433
434!-- poprecip
435qraindiag(:,:)= 0.
436qsnowdiag(:,:)= 0.
437dqreva(:,:)   = 0.
438dqrauto(:,:)  = 0.
439dqrmelt(:,:)  = 0.
440dqrfreez(:,:) = 0.
441dqrcol(:,:)   = 0.
442dqssub(:,:)   = 0.
443dqsauto(:,:)  = 0.
444dqsrim(:,:)   = 0.
445dqsagg(:,:)   = 0.
446dqsfreez(:,:) = 0.
447dqsmelt(:,:)  = 0.
448zqupnew(:)    = 0.
449zqvapclr(:)   = 0.
450
451!-- cloud phase useful variables
452wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:))
453zqliq(:)=0.
454zqice(:)=0.
455zqvapcl(:)=0.
456zqliqth(:)=0.
457zqiceth(:)=0.
458zqvapclth(:)=0.
459sursat_e(:)=0.
460invtau_e(:)=0.
461pticefracturb(:)=.FALSE.
462
463
464!===============================================================================
465!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
466!===============================================================================
467
468ncoreczq=0
469
470DO k = klev, 1, -1
471
472    IF (k.LE.klev-1) THEN
473       iftop=.false.
474    ELSE
475       iftop=.true.
476    ENDIF
477
478    ! Initialisation temperature and specific humidity
479    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
480    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
481    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
482    ! d_t = temperature tendency due to lscp
483    ! The temperature of the overlying layer is updated here because needed for thermalization
484    DO i = 1, klon
485        zt(i)=temp(i,k)
486        zq(i)=qt(i,k)
487        zp(i)=pplay(i,k)
488        zqi_ini(i)=qi_seri(i,k)
489        zcf(i) = 0.
490        zfice(i) = 1.0   ! initialized at 1 as by default we assume mpc to be at ice saturation
491        dzfice(i) = 0.0
492        zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
493        dzficeth(i) = 0.0
494        zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
495        dzficeenv(i) = 0.0
496       
497
498        IF (.not. iftop) THEN
499           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
500           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
501           !--zqs(i) is the saturation specific humidity in the layer above
502           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
503        ENDIF
504        !c_iso init of iso
505    ENDDO
506
507    ! --------------------------------------------------------------------
508    ! P1> Precipitation processes, before cloud formation:
509    !     Thermalization of precipitation falling from the overlying layer AND
510    !     Precipitation evaporation/sublimation/melting
511    !---------------------------------------------------------------------
512   
513    !================================================================
514    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
515    IF ( ok_poprecip ) THEN
516
517      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
518                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
519                        zqvapclr, zqupnew, &
520                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
521                        zrfl, zrflclr, zrflcld, &
522                        zifl, ziflclr, ziflcld, &
523                        dqreva(:,k), dqssub(:,k) &
524                        )
525
526    !================================================================
527    ELSE
528
529      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
530                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
531                        zrfl, zrflclr, zrflcld, &
532                        zifl, ziflclr, ziflcld, &
533                        dqreva(:,k), dqssub(:,k) &
534                        )
535
536    ENDIF ! (ok_poprecip)
537   
[5633]538    ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter
[5614]539    !-------------------------------------------------------
540
541     qtot(:)=zq(:)+zmqc(:)
542     CALL calc_qsat_ecmwf(klon,zt,qtot,zp,RTT,0,.false.,zqs,zdqs)
543     DO i = 1, klon
544            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
545            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
546            IF (zq(i) .LT. 1.e-15) THEN
547                ncoreczq=ncoreczq+1
548                zq(i)=1.e-15
549            ENDIF
550            ! c_iso: do something similar for isotopes
551
552     ENDDO
553   
554    ! -------------------------------------------------------------------------
555    ! P2> Cloud formation including condensation and cloud phase determination
556    !--------------------------------------------------------------------------
557    !
558    ! We always assume a 'fractional cloud' approach
559    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
560    ! is prescribed and depends on large scale variables and boundary layer
561    ! properties)
562    ! The decrease in condensed part due tu latent heating is taken into
563    ! account
564    ! -------------------------------------------------------------------
565
566        ! P2.1> With the PDFs (log-normal, bigaussian)
567        ! cloud properties calculation with the initial values of t and q
568        ! ----------------------------------------------------------------
569
570        ! initialise gammasat and stratiform_or_all_distrib
571        stratiform_or_all_distrib(:)=.TRUE.
572        gammasat(:)=1.
573
574        ! part of the code that is supposed to become obsolete soon
575        !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
576        IF (.NOT. ok_lscp_mergecond) THEN
577          IF (iflag_cld_th.GE.5) THEN
578               ! Cloud cover and content in meshes affected by shallow convection,
579               ! are retrieved from a bi-gaussian distribution of the saturation deficit
580               ! following Jam et al. 2013
581
582               IF (iflag_cloudth_vert.LE.2) THEN
583                  ! Old version of Arnaud Jam
584
585                    CALL cloudth(klon,klev,k,tv,                  &
586                         zq,qta,fraca,                            &
587                         qincloud,ctot,pspsk,paprs,pplay,tla,thl, &
588                         ratqs,zqs,temp,                              &
589                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
590
591
592                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
593                   ! Default version of Arnaud Jam
594                       
595                    CALL cloudth_v3(klon,klev,k,tv,                        &
596                         zq,qta,fraca,                                     &
597                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
598                         ratqs,sigma_qtherm,zqs,temp, &
599                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
600
601
602                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
603                   ! Jean Jouhaud's version, with specific separation between surface and volume
604                   ! cloud fraction Decembre 2018
605
606                    CALL cloudth_v6(klon,klev,k,tv,                        &
607                         zq,qta,fraca,                                     &
608                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
609                         ratqs,zqs,temp, &
610                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
611
612                ENDIF
613
614
615                DO i=1,klon
616                    rneb(i,k)=ctot(i,k)
617                    ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
618                    rneblsvol(i,k)=ctot_vol(i,k)
619                    zqn(i)=qincloud(i)
620                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
621                    qvc(i) = rneb(i,k) * zqs(i)
622                ENDDO
623
624               ! Cloud phase final determination for clouds after "old" cloudth calls
625               ! for those clouds, only temperature dependent phase partitioning (eventually modulated by
626               ! distance to cloud top) is available
627               ! compute distance to cloud top when cloud phase is computed depending on temperature
628               ! and distance to cloud top
629               IF (iflag_t_glace .GE. 4) THEN
630                    CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
631               ENDIF
632               CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
633
634          ENDIF
635
636          IF (iflag_cld_th .EQ. 5) THEN
637               ! When iflag_cld_th=5, we always assume
638                ! bi-gaussian distribution
639            stratiform_or_all_distrib(:) = .FALSE.
640
641          ELSEIF (iflag_cld_th .GE. 6) THEN
642                ! stratiform distribution only when no thermals
643            stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld
644
645          ENDIF
646
647        ENDIF ! .not. ok_lscp_mergecond
648        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
649       
650        DT(:) = 0.
651        n_i(:)=0
652        Tbef(:)=zt(:)
653        qlibef(:)=0.
654        Tbefth(:)=tla(:,k)*pspsk(:,k)
[5633]655        IF (k .GT. 1) THEN
656         Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1)
657        ELSE
658         Tbefthm1(:)=Tbefth(:)       
659        ENDIF
[5647]660        zqth(:)=qta(:,k)
661        zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:)
[5614]662
663        ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth
664        ! in case of ok_lscp_mergecond)
665        ! We iterate here to take into account the change in qsat(T) and ice fraction
666        ! during the condensation process
667        ! the increment in temperature is calculated using the first principle of
668        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
669        ! and a clear sky fraction)
670        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
671
672        DO iter=1,iflag_fisrtilp_qsat+1
673
674                keepgoing(:) = .FALSE.
675
676                DO i=1,klon
677
678                    ! keepgoing = .true. while convergence is not satisfied
679
680                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN
681
682                        ! if not convergence:
683                        ! we calculate a new iteration
684                        keepgoing(i) = .TRUE.
685
686                        ! P2.2.1> cloud fraction and condensed water mass calculation
687                        ! Calculated variables:
688                        ! rneb : cloud fraction
689                        ! zqn : total water within the cloud
690                        ! zcond: mean condensed water within the mesh
691                        ! rhcl: clear-sky relative humidity
692                        !---------------------------------------------------------------
693
694                        ! new temperature that only serves in the iteration process:
695                        Tbef(i)=Tbef(i)+DT(i)
696
697                        ! total water including mass of precip
698                        qtot(i)=zq(i)+zmqc(i)
699
700                      ENDIF
701
702                  ENDDO
703
704                  ! Calculation of saturation specific humidity
705                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs)
706                  ! also in thermals
707                  CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth)
708
709                  IF (iflag_icefrac .GE. 1) THEN
710                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
711                  ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD
712                  ! we make this option works only for the physically-based and tke-dependent param from Raillard et al.
713                  ! (iflag_icefrac>=1)
714                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl)
715                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi)
716                      DO i=1,klon
717                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
718                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
719                      ENDDO
720                  ENDIF
721                  IF (iflag_icefrac .GE. 2) THEN
722                  ! same adjustment for thermals
723                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth)
724                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith)
725                      DO i=1,klon
726                           zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i)
727                           zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i)
728                      ENDDO
729                  ENDIF
730
731                  CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt)
732                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
733                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
734
735                  ! Cloud condensation based on subgrid pdf
736                  !--AB Activates a condensation scheme that allows for
737                  !--ice supersaturation and contrails evolution from aviation
738                  IF (ok_ice_supersat) THEN
739
740                    !--Calculate the shear value (input for condensation and ice supersat)
741                    DO i = 1, klon
742                      !--Cell thickness [m]
743                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
744                      IF ( iftop ) THEN
745                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
746                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
747                      ELSEIF ( k .EQ. 1 ) THEN
748                        ! surface
749                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
750                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
751                      ELSE
752                        ! other layers
753                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
754                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
755                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
756                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
757                      ENDIF
758                    ENDDO
759
760                    !---------------------------------------------
761                    !--   CONDENSATION AND ICE SUPERSATURATION  --
762                    !---------------------------------------------
763
764                    CALL condensation_ice_supersat( &
765                        klon, dtime, missing_val, &
766                        zp, paprs(:,k), paprs(:,k+1), &
767                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
768                        shear, tke_dissip(:,k), cell_area, &
769                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
770                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
771                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
772                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
773                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
774                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
775                        flight_dist(:,k), flight_h2o(:,k), &
776                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
777
778
779                  ELSE
780                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
781
782                  CALL condensation_lognormal( &
783                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
784                       keepgoing, rneb(:,k), zqn, qvc)
785
786
787                  ENDIF ! .NOT. ok_ice_supersat
788
789                  ! Volume cloud fraction
790                  rneblsvol(:,k)=rneb(:,k)
791
792
793                  IF  (ok_lscp_mergecond) THEN
794                      ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only
795                      CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), &
796                                     pspsk(:,k),zp,tla(:,k), &
797                                     ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), &
798                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
799                  ENDIF
800
801
802
803                  ! Cloud phase determination
804
805
806                  IF (iflag_icefrac .LE. 1) THEN
807                     ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere
808                     IF (iflag_t_glace .GE. 4) THEN
809                        CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
810                     ENDIF
811                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
812                  ENDIF
813
814                  IF (iflag_icefrac .EQ. 1) THEN
815                     ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics
816                     ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025)
817                     ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds
818                     pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
819                     DO i=1,klon
820                        qincloud(i)=zqn(i)
821                        zcf(i)=MIN(MAX(rneb(i,k), 0.),1.)
822                        sursat_e(i) = 0.
823                        invtau_e(i) = 0.
824                     ENDDO
[5633]825                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini,    &
826                     ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
827                     zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k))                     
[5614]828                     DO i=1,klon
829                        IF (pticefracturb(i)) THEN
830                          zfice(i)=zficeenv(i)
831                          dzfice(i)=dzficeenv(i)
832                        ENDIF
833                     ENDDO
834                   
835
836                  ELSEIF (iflag_icefrac .GE. 2) THEN
837                     ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption)
838                     ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment
839                     ! hence we assume that no snow falls in thermals.
840                     ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T
841                     pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
842                     !Thermals
843                     DO i=1,klon
844                        IF (fraca(i,k) .GT. min_frac_th_cld) THEN
845                           zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)*fraca(i,k)
846                           qincloud(i)=zqn(i)/fraca(i,k)
847                        ELSE
848                           zcf(i) = 0.
849                           qincloud(i) = 0.
850                        ENDIF
851                        sursat_e(i)=cloudth_senv(i,k)/zqsi(i)
[5647]852                        invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
[5614]853                     ENDDO
[5633]854                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), zqi_ini,  &
855                     zeroklon, znebprecipcld, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, &
856                     zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k))
[5614]857                     !Environment
858                     DO i=1,klon
859                        qincloud(i)=zqn(i)/(1.-fraca(i,k))
860                        zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)*(1.-fraca(i,k))
[5633]861                        IF (k .GT. 1) THEN
862                           ! evaluate the mixing sursaturation using saturation deficit at level below
[5647]863                           ! as air pacels detraining into clouds have not (less) seen yet entrainement from above
[5633]864                           sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i)))
[5647]865                           ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho
866                           invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
[5633]867                        ELSE
868                           sursat_e(i)=0.
869                           invtau_e(i)=0.
870                        ENDIF
[5614]871                     ENDDO
[5633]872                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini,    &
873                     ziflcld, znebprecipcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
874                     zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
[5614]875 
876                    ! adjust zfice to account for condensates in thermals'fraction
877                     DO i=1,klon
878                        IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN
879                          zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i)))
880                          dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb
881                        ENDIF
882                     ENDDO
883
884                  ENDIF ! iflag_icefrac
885
886
887
888
889                  DO i=1,klon
890                      IF (keepgoing(i)) THEN
891
892                       ! P2.2.2> Approximative calculation of temperature variation DT
893                       !           due to condensation.
894                       ! Calculated variables:
895                       ! dT : temperature change due to condensation
896                       !---------------------------------------------------------------
897
898               
899                        IF (zfice(i).LT.1) THEN
900                            cste=RLVTT
901                        ELSE
902                            cste=RLSTT
903                        ENDIF
904                       
905                        IF ( ok_unadjusted_clouds ) THEN
906                          !--AB We relax the saturation adjustment assumption
907                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
908                          IF ( rneb(i,k) .GT. eps ) THEN
909                            qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
910                          ELSE
911                            qlibef(i) = 0.
912                          ENDIF
913                        ELSE
914                          qlibef(i)=max(0.,zqn(i)-zqs(i))
915                        ENDIF
916
917                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
918                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(i)
919                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
920                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
921                              *qlibef(i)*dzfice(i)
922                        ! here we update a provisory temperature variable that only serves in the iteration
923                        ! process
924                        DT(i)=num/denom
925                        n_i(i)=n_i(i)+1
926
927                    ENDIF ! end keepgoing
928
929                ENDDO     ! end loop on i
930
931        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
932
933        ! P2.2> Final quantities calculation
934        ! Calculated variables:
935        !   rneb : cloud fraction
936        !   zcond: mean condensed water in the mesh
937        !   zqn  : mean water vapor in the mesh
938        !   zfice: ice fraction in clouds
939        !   zt   : temperature
940        !   rhcl : clear-sky relative humidity
941        ! ----------------------------------------------------------------
942
943
944        ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type
945        !---------------------------------------------------------------------------------------------
946        DO i=1, klon
947                ! Checks on rneb, rhcl and zqn
948                IF (rneb(i,k) .LE. 0.0) THEN
949                    zqn(i) = 0.0
950                    rneb(i,k) = 0.0
951                    zcond(i) = 0.0
952                    rhcl(i,k)=zq(i)/zqs(i)
953                ELSE IF (rneb(i,k) .GE. 1.0) THEN
954                    zqn(i) = zq(i)
955                    rneb(i,k) = 1.0
956                    IF ( ok_unadjusted_clouds ) THEN
957                      !--AB We relax the saturation adjustment assumption
958                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
959                      zcond(i) = MAX(0., zqn(i) - qvc(i))
960                    ELSE
961                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
962                    ENDIF
963                    rhcl(i,k)=1.0
964                ELSE
965                    IF ( ok_unadjusted_clouds ) THEN
966                      !--AB We relax the saturation adjustment assumption
967     
968                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
969                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
970                    ELSE
971                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
972                    ENDIF
973
974                    ! following line is very strange and probably wrong
975                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
976                    ! Correct calculation of clear-sky relative humidity. To activate once
977                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
978                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
979                    !--This is slighly approximated, the actual formula is
980                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
981                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
982                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
983                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
984                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
985
986                ENDIF
987
988                ! water vapor update
989                zq(i) = zq(i) - zcond(i)
990
991                       
992                ! temperature update due to phase change
993                zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
994                &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
995                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
996        ENDDO
997
998    ! ----------------------------------------------------------------
999    ! P3> Precipitation processes, after cloud formation
1000    !     - precipitation formation, melting/freezing
1001    ! ----------------------------------------------------------------
1002
1003    ! Initiate the quantity of liquid and solid condensates
1004    ! Note that in the following, zcond is the total amount of condensates
1005    ! including newly formed precipitations (i.e., condensates formed by the
1006    ! condensation process), while zoliq is the total amount of condensates in
1007    ! the cloud (i.e., on which precipitation processes have applied)
1008    DO i = 1, klon
1009      zoliq(i)  = zcond(i)
1010      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1011      zoliqi(i) = zcond(i) * zfice(i)
1012    ENDDO
1013
1014    !================================================================
1015    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1016    IF (ok_poprecip) THEN
1017
1018      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, &
1019                            ctot_vol(:,k), ptconv(:,k), &
1020                            zt, zq, zoliql, zoliqi, zfice, &
1021                            rneb(:,k), znebprecipclr, znebprecipcld, &
1022                            zrfl, zrflclr, zrflcld, &
1023                            zifl, ziflclr, ziflcld, &
1024                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1025                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1026                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1027                            dqsmelt(:,k), dqsfreez(:,k) &
1028                            )
1029      DO i = 1, klon
1030          zoliq(i) = zoliql(i) + zoliqi(i)
1031      ENDDO
1032
1033    !================================================================
1034    ELSE
1035
1036      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
1037                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
1038                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
1039                            rneb(:,k), znebprecipclr, znebprecipcld, &
1040                            zneb, tot_zneb, zrho_up, zvelo_up, &
1041                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1042                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
1043                            )
1044
1045    ENDIF ! ok_poprecip
1046
1047    ! End of precipitation processes after cloud formation
1048    ! ----------------------------------------------------
1049
1050    !----------------------------------------------------------------------
1051    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1052    !----------------------------------------------------------------------
1053
1054    DO i=1,klon
1055
1056      IF (ok_poprecip) THEN
1057        IF (ok_radocond_snow) THEN
1058          radocond(i,k) = zoliq(i)
1059          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1060        ELSE
1061          radocond(i,k) = zoliq(i)
1062          zradoice(i) = zoliqi(i)
1063        ENDIF
1064      ELSE
1065        radocond(i,k) = zradocond(i)
1066      ENDIF
1067
1068      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1069      ! by the radiation scheme
1070      IF (radocond(i,k) .GT. 0.) THEN
1071        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1072      ENDIF
1073    ENDDO
1074
1075    ! ----------------------------------------------------------------
1076    ! P5> Wet scavenging
1077    ! ----------------------------------------------------------------
1078
1079    !Scavenging through nucleation in the layer
1080
1081    DO i = 1,klon
1082       
1083        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1084            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1085        ELSE
1086            beta(i,k) = 0.
1087        ENDIF
1088
1089        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1090
1091        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1092
1093            IF (temp(i,k) .GE. temp_nowater) THEN
1094                zalpha_tr = a_tr_sca(3)
1095            ELSE
1096                zalpha_tr = a_tr_sca(4)
1097            ENDIF
1098
1099            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1100            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1101
1102            ! Nucleation with a factor of -1 instead of -0.5
1103            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1104
1105        ENDIF
1106
1107    ENDDO
1108
1109    ! Scavenging through impaction in the underlying layer
1110
1111    DO kk = k-1, 1, -1
1112
1113        DO i = 1, klon
1114
1115            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1116
1117                IF (temp(i,kk) .GE. temp_nowater) THEN
1118                    zalpha_tr = a_tr_sca(1)
1119                ELSE
1120                    zalpha_tr = a_tr_sca(2)
1121                ENDIF
1122
1123              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1124              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1125
1126            ENDIF
1127
1128        ENDDO
1129
1130    ENDDO
1131   
1132    !------------------------------------------------------------
1133    ! P6 > write diagnostics and outputs
1134    !------------------------------------------------------------
1135   
1136    !--AB Write diagnostics and tracers for ice supersaturation
1137    IF ( ok_ice_supersat ) THEN
1138      CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs)
1139      CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs)
1140
1141      DO i = 1, klon
1142
1143        IF ( zoliq(i) .LE. 0. ) THEN
1144          !--If everything was precipitated, the remaining empty cloud is dissipated
1145          !--and everything is transfered to the subsaturated clear sky region
1146          rneb(i,k) = 0.
1147          qvc(i) = 0.
1148        ENDIF
1149
1150        cf_seri(i,k) = rneb(i,k)
1151
1152        IF ( .NOT. ok_unadjusted_clouds ) THEN
1153          qvc(i) = zqs(i) * rneb(i,k)
1154        ENDIF
1155        IF ( zq(i) .GT. min_qParent ) THEN
1156          rvc_seri(i,k) = qvc(i) / zq(i)
1157        ELSE
1158          rvc_seri(i,k) = min_ratio
1159        ENDIF
1160        !--The MIN barrier is NEEDED because of:
1161        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1162        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1163        !--The MAX barrier is a safeguard that should not be activated
1164        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1165
1166        !--Diagnostics
1167        gamma_cond(i,k) = gammasat(i)
1168        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1169        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1170        qcld(i,k) = qvc(i) + zoliq(i)
1171      ENDDO
1172    ENDIF
1173
1174    ! Outputs:
1175    !-------------------------------
1176    ! Precipitation fluxes at layer interfaces
1177    ! + precipitation fractions +
1178    ! temperature and water species tendencies
1179    DO i = 1, klon
1180        psfl(i,k)=zifl(i)
1181        prfl(i,k)=zrfl(i)
1182        pfraclr(i,k)=znebprecipclr(i)
1183        pfracld(i,k)=znebprecipcld(i)
1184        distcltop(i,k)=zdistcltop(i)
1185        temp_cltop(i,k)=ztemp_cltop(i)
1186        d_q(i,k) = zq(i) - qt(i,k)
1187        d_t(i,k) = zt(i) - temp(i,k)
1188
1189        IF (ok_bug_phase_lscp) THEN
1190           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1191           d_qi(i,k) = zfice(i)*zoliq(i)
1192        ELSE
1193           d_ql(i,k) = zoliql(i)
1194           d_qi(i,k) = zoliqi(i)   
1195        ENDIF
1196
1197    ENDDO
1198
1199
1200ENDDO ! loop on k from top to bottom
1201
1202
1203  ! Rain or snow at the surface (depending on the first layer temperature)
1204  DO i = 1, klon
1205      snow(i) = zifl(i)
1206      rain(i) = zrfl(i)
1207      ! c_iso final output
1208  ENDDO
1209
1210  IF (ncoreczq>0) THEN
1211      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1212  ENDIF
1213
1214
1215RETURN
1216
1217END SUBROUTINE lscp
1218!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1219
1220END MODULE lmdz_lscp_main
Note: See TracBrowser for help on using the repository browser.