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

Last change on this file was 5765, checked in by rkazeroni, 5 weeks ago

For GPU porting:

  • Add one Fortran comment for each ported param to specify possible names for horizontal variables used in the param
  • Add Fortran comment to exclude a file and its routines from GPU porting

Note:

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