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

Last change on this file since 5573 was 5573, checked in by aborella, 4 months ago

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