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

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

Multiple changes:

  • tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

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