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

Last change on this file since 5691 was 5691, checked in by aborella, 3 weeks ago

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