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

Last change on this file since 5456 was 5456, checked in by aborella, 37 hours ago

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