source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90 @ 5575

Last change on this file since 5575 was 5575, checked in by aborella, 3 months ago

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