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

Last change on this file since 5614 was 5614, checked in by evignon, 8 weeks ago

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