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

Last change on this file since 5602 was 5602, checked in by aborella, 2 months ago

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