source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 @ 5390

Last change on this file since 5390 was 5384, checked in by evignon, 7 weeks ago

ajout de diagnostiques des tendances de precip, Audran Borella
File size: 71.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,ptconv,ratqs,sigma_qtherm, &
10     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
11     pfraclr, pfracld,                                  &
12     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
13     radocond, radicefrac, rain, snow,                  &
14     frac_impa, frac_nucl, beta,                        &
15     prfl, psfl, rhcl, qta, fraca,                      &
16     tv, pspsk, tla, thl, iflag_cld_th,                 &
17     iflag_ice_thermo, distcltop, temp_cltop,           &
18     tke, tke_dissip,                                   &
19     cell_area,                                         &
20     cf_seri, rvc_seri, u_seri, v_seri,                 &
21     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
22     ratio_qi_qtot, dcf_sub, dcf_con, dcf_mix,          &
23     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
24     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
25     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
26     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
27     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
28     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
29     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
30     dqsmelt, dqsfreez)
31
32!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
34!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
35!--------------------------------------------------------------------------------
36! Date: 01/2021
37!--------------------------------------------------------------------------------
38! Aim: Large Scale Clouds and Precipitation (LSCP)
39!
40! This code is a new version of the fisrtilp.F90 routine, which is itself a
41! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
42! and 'ilp' (il pleut, L. Li)
43!
44! Compared to the original fisrtilp code, lscp
45! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
46! -> consider always precipitation thermalisation (fl_cor_ebil>0)
47! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
48! -> option "oldbug" by JYG has been removed
49! -> iflag_t_glace >0 always
50! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
51! -> rectangular distribution from L. Li no longer available
52! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
53!--------------------------------------------------------------------------------
54! References:
55!
56! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
57! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
58! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
59! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
60! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
61! - Touzze-Peifert Ludo, PhD thesis, p117-124
62! -------------------------------------------------------------------------------
63! Code structure:
64!
65! P0>     Thermalization of the precipitation coming from the overlying layer
66! P1>     Evaporation of the precipitation (falling from the k+1 level)
67! P2>     Cloud formation (at the k level)
68! P2.A.1> With the PDFs, calculation of cloud properties using the inital
69!         values of T and Q
70! P2.A.2> Coupling between condensed water and temperature
71! P2.A.3> Calculation of final quantities associated with cloud formation
72! P2.B>   Release of Latent heat after cloud formation
73! P3>     Autoconversion to precipitation (k-level)
74! P4>     Wet scavenging
75!------------------------------------------------------------------------------
76! Some preliminary comments (JBM) :
77!
78! The cloud water that the radiation scheme sees is not the same that the cloud
79! water used in the physics and the dynamics
80!
81! During the autoconversion to precipitation (P3 step), radocond (cloud water used
82! by the radiation scheme) is calculated as an average of the water that remains
83! in the cloud during the precipitation and not the water remaining at the end
84! of the time step. The latter is used in the rest of the physics and advected
85! by the dynamics.
86!
87! In summary:
88!
89! Radiation:
90! xflwc(newmicro)+xfiwc(newmicro) =
91!   radocond=lwcon(nc)+iwcon(nc)
92!
93! Notetheless, be aware of:
94!
95! radocond .NE. ocond(nc)
96! i.e.:
97! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
98! but oliq+(ocond-oliq) .EQ. ocond
99! (which is not trivial)
100!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101!
102
103! USE de modules contenant des fonctions.
104USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
105USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
106USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
107USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
108USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
109USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
110
111! Use du module lmdz_lscp_ini contenant les constantes
112USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
113USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
114USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
115USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
116USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
117USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
118USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
119USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
120USE lmdz_lscp_ini, ONLY : ok_poprecip
121USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
122
123IMPLICIT NONE
124
125!===============================================================================
126! VARIABLES DECLARATION
127!===============================================================================
128
129! INPUT VARIABLES:
130!-----------------
131
132  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
133  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
134  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
135
136  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
137  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
139  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
140  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
141  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
142  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
143                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
144  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
145
146  !Inputs associated with thermal plumes
147
148  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
153  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
154  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
155
156  ! INPUT/OUTPUT variables
157  !------------------------
158 
159  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
160  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
161
162
163  ! INPUT/OUTPUT condensation and ice supersaturation
164  !--------------------------------------------------
165  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
166  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratio_qi_qtot    ! solid specific water to total specific water ratio [-]
167  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
170  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
171
172  ! INPUT/OUTPUT aviation
173  !--------------------------------------------------
174  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
175  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
176 
177  ! OUTPUT variables
178  !-----------------
179
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
186  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
187  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
188  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
191  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
192  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
193  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
194  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
195  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
196  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
197  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
201
202  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
203 
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
206 
207  ! for condensation and ice supersaturation
208
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
228
229  ! for contrails and aviation
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
239
240
241  ! for POPRECIP
242
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
256
257  ! for thermals
258
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
263
264
265  ! LOCAL VARIABLES:
266  !----------------
267  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
268  REAL,DIMENSION(klon) :: qice_ini
269  REAL :: zct, zcl,zexpo
270  REAL, DIMENSION(klon,klev) :: ctot
271  REAL, DIMENSION(klon,klev) :: ctot_vol
272  REAL, DIMENSION(klon) :: zqs, zdqs
273  REAL :: zdelta, zcor, zcvm5
274  REAL, DIMENSION(klon) :: zdqsdT_raw
275  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
276  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
277  REAL :: num,denom
278  REAL :: cste
279  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
280  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
281  REAL :: erf
282  REAL, DIMENSION(klon) :: zfice_th
283  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
284  REAL, DIMENSION(klon) :: zrfl, zrfln
285  REAL :: zqev, zqevt
286  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
287  REAL :: zqev0,zqevi, zqevti
288  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
289  REAL, DIMENSION(klon) :: zoliql, zoliqi
290  REAL, DIMENSION(klon) :: zt
291  REAL, DIMENSION(klon,klev) :: zrho
292  REAL, DIMENSION(klon) :: zdz,iwc
293  REAL :: zchau,zfroi
294  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
295  REAL :: zmelt,zrain,zsnow,zprecip
296  REAL, DIMENSION(klon) :: dzfice
297  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
298  REAL :: zsolid
299  REAL, DIMENSION(klon) :: qtot, qzero
300  REAL, DIMENSION(klon) :: dqsl,dqsi
301  REAL :: smallestreal
302  !  Variables for Bergeron process
303  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
304  REAL, DIMENSION(klon) :: zqpreci, zqprecl
305  ! Variables precipitation energy conservation
306  REAL, DIMENSION(klon) :: zmqc
307  REAL :: zalpha_tr
308  REAL :: zfrac_lessi
309  REAL, DIMENSION(klon) :: zprec_cond
310  REAL :: zmair
311  REAL :: zcpair, zcpeau
312  REAL, DIMENSION(klon) :: zlh_solid
313  REAL, DIMENSION(klon) :: ztupnew
314  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
315  REAL :: zm_solid         ! for liquid -> solid conversion
316  REAL, DIMENSION(klon) :: zrflclr, zrflcld
317  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
318  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
319  REAL, DIMENSION(klon) :: ziflclr, ziflcld
320  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
321  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
322  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
323  REAL, DIMENSION(klon,klev) :: velo
324  REAL :: vr, ffallv
325  REAL :: qlmpc, qimpc, rnebmpc
326  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
327  REAL :: effective_zneb
328  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
329  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
330 
331  ! for condensation and ice supersaturation
332  REAL, DIMENSION(klon) :: qvc, shear
333  REAL :: delta_z
334  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
335  ! Constants used for calculating ratios that are advected (using a parent-child
336  ! formalism). This is not done in the dynamical core because at this moment,
337  ! only isotopes can use this parent-child formalism. Note that the two constants
338  ! are the same as the one use in the dynamical core, being also defined in
339  ! dyn3d_common/infotrac.F90
340  REAL :: min_qParent, min_ratio
341
342  INTEGER i, k, n, kk, iter
343  INTEGER, DIMENSION(klon) :: n_i
344  INTEGER ncoreczq
345  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
346  LOGICAL iftop
347
348  LOGICAL, DIMENSION(klon) :: lognormale
349  LOGICAL, DIMENSION(klon) :: keepgoing
350
351  CHARACTER (len = 20) :: modname = 'lscp'
352  CHARACTER (len = 80) :: abort_message
353 
354
355!===============================================================================
356! INITIALISATION
357!===============================================================================
358
359! Few initial checks
360
361
362IF (iflag_fisrtilp_qsat .LT. 0) THEN
363    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
364    CALL abort_physic(modname,abort_message,1)
365ENDIF
366
367! Few initialisations
368
369znebprecip(:)=0.0
370ctot_vol(1:klon,1:klev)=0.0
371rneblsvol(1:klon,1:klev)=0.0
372smallestreal=1.e-9
373znebprecipclr(:)=0.0
374znebprecipcld(:)=0.0
375mpc_bl_points(:,:)=0
376
377IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
378
379! AA for 'safety' reasons
380zalpha_tr   = 0.
381zfrac_lessi = 0.
382beta(:,:)= 0.
383
384! Initialisation of variables:
385
386prfl(:,:) = 0.0
387psfl(:,:) = 0.0
388d_t(:,:) = 0.0
389d_q(:,:) = 0.0
390d_ql(:,:) = 0.0
391d_qi(:,:) = 0.0
392rneb(:,:) = 0.0
393pfraclr(:,:)=0.0
394pfracld(:,:)=0.0
395cldfraliq(:,:)=0.
396sigma2_icefracturb(:,:)=0.
397mean_icefracturb(:,:)=0.
398radocond(:,:) = 0.0
399radicefrac(:,:) = 0.0
400frac_nucl(:,:) = 1.0
401frac_impa(:,:) = 1.0
402rain(:) = 0.0
403snow(:) = 0.0
404zoliq(:)=0.0
405zfice(:)=0.0
406dzfice(:)=0.0
407zfice_turb(:)=0.0
408dzfice_turb(:)=0.0
409zqprecl(:)=0.0
410zqpreci(:)=0.0
411zrfl(:) = 0.0
412zifl(:) = 0.0
413ziflprev(:)=0.0
414zneb(:) = seuil_neb
415zrflclr(:) = 0.0
416ziflclr(:) = 0.0
417zrflcld(:) = 0.0
418ziflcld(:) = 0.0
419tot_zneb(:) = 0.0
420tot_znebn(:) = 0.0
421d_tot_zneb(:) = 0.0
422qzero(:) = 0.0
423zdistcltop(:)=0.0
424ztemp_cltop(:) = 0.0
425ztupnew(:)=0.0
426
427distcltop(:,:)=0.
428temp_cltop(:,:)=0.
429
430!--Ice supersaturation
431gamma_cond(:,:) = 1.
432qissr(:,:)      = 0.
433issrfra(:,:)    = 0.
434dcf_sub(:,:)    = 0.
435dcf_con(:,:)    = 0.
436dcf_mix(:,:)    = 0.
437dqi_adj(:,:)    = 0.
438dqi_sub(:,:)    = 0.
439dqi_con(:,:)    = 0.
440dqi_mix(:,:)    = 0.
441dqvc_adj(:,:)   = 0.
442dqvc_sub(:,:)   = 0.
443dqvc_con(:,:)   = 0.
444dqvc_mix(:,:)   = 0.
445fcontrN(:,:)    = 0.
446fcontrP(:,:)    = 0.
447Tcontr(:,:)     = missing_val
448qcontr(:,:)     = missing_val
449qcontr2(:,:)    = missing_val
450dcf_avi(:,:)    = 0.
451dqi_avi(:,:)    = 0.
452dqvc_avi(:,:)   = 0.
453qvc(:)          = 0.
454shear(:)        = 0.
455min_qParent     = 1.e-30
456min_ratio       = 1.e-16
457
458!-- poprecip
459qraindiag(:,:)= 0.
460qsnowdiag(:,:)= 0.
461dqreva(:,:)   = 0.
462dqrauto(:,:)  = 0.
463dqrmelt(:,:)  = 0.
464dqrfreez(:,:) = 0.
465dqrcol(:,:)   = 0.
466dqssub(:,:)   = 0.
467dqsauto(:,:)  = 0.
468dqsrim(:,:)   = 0.
469dqsagg(:,:)   = 0.
470dqsfreez(:,:) = 0.
471dqsmelt(:,:)  = 0.
472zqupnew(:)    = 0.
473zqvapclr(:)   = 0.
474
475
476
477!c_iso: variable initialisation for iso
478
479
480!===============================================================================
481!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
482!===============================================================================
483
484ncoreczq=0
485
486DO k = klev, 1, -1
487
488    IF (k.LE.klev-1) THEN
489       iftop=.false.
490    ELSE
491       iftop=.true.
492    ENDIF
493
494    ! Initialisation temperature and specific humidity
495    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
496    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
497    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
498    ! d_t = temperature tendency due to lscp
499    ! The temperature of the overlying layer is updated here because needed for thermalization
500    DO i = 1, klon
501        zt(i)=temp(i,k)
502        zq(i)=qt(i,k)
503        IF (.not. iftop) THEN
504           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
505           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
506           !--zqs(i) is the saturation specific humidity in the layer above
507           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
508        ENDIF
509        !c_iso init of iso
510    ENDDO
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                              zrfl, zrflclr, zrflcld, &
520                              zifl, ziflclr, ziflcld, &
521                              dqreva(:,k),dqssub(:,k) &
522                             )
523
524    !================================================================
525    ELSE
526
527    ! --------------------------------------------------------------------
528    ! P1> Thermalization of precipitation falling from the overlying layer
529    ! --------------------------------------------------------------------
530    ! Computes air temperature variation due to enthalpy transported by
531    ! precipitation. Precipitation is then thermalized with the air in the
532    ! layer.
533    ! The precipitation should remain thermalized throughout the different
534    ! thermodynamical transformations.
535    ! The corresponding water mass should
536    ! be added when calculating the layer's enthalpy change with
537    ! temperature
538    ! See lmdzpedia page todoan
539    ! todoan: check consistency with ice phase
540    ! todoan: understand why several steps
541    ! ---------------------------------------------------------------------
542   
543    IF (iftop) THEN
544
545        DO i = 1, klon
546            zmqc(i) = 0.
547        ENDDO
548
549    ELSE
550
551        DO i = 1, klon
552 
553            zmair=(paprs(i,k)-paprs(i,k+1))/RG
554            ! no condensed water so cp=cp(vapor+dry air)
555            ! RVTMP2=rcpv/rcpd-1
556            zcpair=RCPD*(1.0+RVTMP2*zq(i))
557            zcpeau=RCPD*RVTMP2
558
559            ! zmqc: precipitation mass that has to be thermalized with
560            ! layer's air so that precipitation at the ground has the
561            ! same temperature as the lowermost layer
562            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
563            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
564            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
565             / (zcpair + zmqc(i)*zcpeau)
566
567        ENDDO
568 
569    ENDIF
570
571    qice_ini = zq(:)*ratio_qi_qtot(:,k)
572    ! --------------------------------------------------------------------
573    ! P2> Precipitation evaporation/sublimation/melting
574    ! --------------------------------------------------------------------
575    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
576    ! --------------------------------------------------------------------
577   
578    IF (iflag_evap_prec.GE.1) THEN
579
580        ! Calculation of saturation specific humidity
581        ! depending on temperature:
582        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,0,.false.,zqs,zdqs)
583        ! wrt liquid water
584        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsl,dqsl)
585        ! wrt ice
586        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsi,dqsi)
587
588        DO i = 1, klon
589
590            ! if precipitation
591            IF (zrfl(i)+zifl(i).GT.0.) THEN
592
593            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
594            ! c_iso: likely important to distinguish cs from neb isotope precipitation
595
596                IF (iflag_evap_prec.GE.4) THEN
597                    zrfl(i) = zrflclr(i)
598                    zifl(i) = ziflclr(i)
599                ENDIF
600
601                IF (iflag_evap_prec.EQ.1) THEN
602                    znebprecip(i)=zneb(i)
603                ELSE
604                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
605                ENDIF
606
607                IF (iflag_evap_prec.GT.4) THEN
608                ! Max evaporation not to saturate the clear sky precip fraction
609                ! i.e. the fraction where evaporation occurs
610                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
611                ELSEIF (iflag_evap_prec .EQ. 4) THEN
612                ! Max evaporation not to saturate the whole mesh
613                ! Pay attention -> lead to unrealistic and excessive evaporation
614                    zqev0 = MAX(0.0, zqs(i)-zq(i))
615                ELSE
616                ! Max evap not to saturate the fraction below the cloud
617                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
618                ENDIF
619
620                ! Evaporation of liquid precipitation coming from above
621                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
622                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
623                ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
624
625                IF (iflag_evap_prec.EQ.3) THEN
626                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
627                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
628                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
629                ELSE IF (iflag_evap_prec.GE.4) THEN
630                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
631                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
632                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
633                ELSE
634                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
635                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
636                ENDIF
637
638                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
639                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
640
641                ! sublimation of the solid precipitation coming from above
642                IF (iflag_evap_prec.EQ.3) THEN
643                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
644                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
645                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
646                ELSE IF (iflag_evap_prec.GE.4) THEN
647                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
648                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
649                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
650                ELSE
651                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
652                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
653                ENDIF
654
655                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
656                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
657
658                ! A. JAM
659                ! Evaporation limit: we ensure that the layer's fraction below
660                ! the cloud or the whole mesh (depending on iflag_evap_prec)
661                ! does not reach saturation. In this case, we
662                ! redistribute zqev0 conserving the ratio liquid/ice
663
664                IF (zqevt+zqevti.GT.zqev0) THEN
665                    zqev=zqev0*zqevt/(zqevt+zqevti)
666                    zqevi=zqev0*zqevti/(zqevt+zqevti)
667                ELSE
668                    zqev=zqevt
669                    zqevi=zqevti
670                ENDIF
671
672                !--Diagnostics
673                dqreva(i,k) = - zqev / dtime
674                dqssub(i,k) = - zqevti / dtime
675
676                ! New solid and liquid precipitation fluxes after evap and sublimation
677                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
678                         /RG/dtime)
679                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
680                         /RG/dtime)
681               
682
683                ! vapor, temperature, precip fluxes update
684                ! vapor is updated after evaporation/sublimation (it is increased)
685                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
686                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
687                ! zmqc is the total condensed water in the precip flux (it is decreased)
688                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
689                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
690                ! air and precip temperature (i.e., gridbox temperature)
691                ! is updated due to latent heat cooling
692                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
693                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
694                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
695                + (zifln(i)-zifl(i))                      &
696                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
697                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
698
699                ! New values of liquid and solid precipitation
700                zrfl(i) = zrfln(i)
701                zifl(i) = zifln(i)
702
703                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
704                !                         due to evap + sublim                                     
705                                           
706
707                IF (iflag_evap_prec.GE.4) THEN
708                    zrflclr(i) = zrfl(i)
709                    ziflclr(i) = zifl(i)
710                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
711                        znebprecipclr(i) = 0.0
712                    ENDIF
713                    zrfl(i) = zrflclr(i) + zrflcld(i)
714                    zifl(i) = ziflclr(i) + ziflcld(i)
715                ENDIF
716
717                ! c_iso duplicate for isotopes or loop on isotopes
718
719                ! Melting:
720                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
721                ! precip fraction that is melted
722                zmelt = MIN(MAX(zmelt,0.),1.)
723
724                ! update of rainfall and snowfall due to melting
725                IF (iflag_evap_prec.GE.4) THEN
726                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
727                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
728                    zrfl(i)=zrflclr(i)+zrflcld(i)
729                ELSE
730                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
731                ENDIF
732
733
734                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
735
736                ! Latent heat of melting because of precipitation melting
737                ! NB: the air + precip temperature is simultaneously updated
738                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
739                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
740               
741                IF (iflag_evap_prec.GE.4) THEN
742                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
743                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
744                    zifl(i)=ziflclr(i)+ziflcld(i)
745                ELSE
746                    zifl(i)=zifl(i)*(1.-zmelt)
747                ENDIF
748
749            ELSE
750                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
751                znebprecip(i)=0.
752
753            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
754
755        ENDDO ! loop on klon
756
757    ENDIF ! (iflag_evap_prec>=1)
758
759    ENDIF ! (ok_poprecip)
760
761    ! --------------------------------------------------------------------
762    ! End precip evaporation
763    ! --------------------------------------------------------------------
764   
765    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
766    !-------------------------------------------------------
767
768     qtot(:)=zq(:)+zmqc(:)
769     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
770     DO i = 1, klon
771            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
772            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
773            IF (zq(i) .LT. 1.e-15) THEN
774                ncoreczq=ncoreczq+1
775                zq(i)=1.e-15
776            ENDIF
777            ! c_iso: do something similar for isotopes
778
779     ENDDO
780   
781    ! --------------------------------------------------------------------
782    ! P2> Cloud formation
783    !---------------------------------------------------------------------
784    !
785    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
786    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
787    ! is prescribed and depends on large scale variables and boundary layer
788    ! properties)
789    ! The decrease in condensed part due tu latent heating is taken into
790    ! account
791    ! -------------------------------------------------------------------
792
793        ! P2.1> With the PDFs (log-normal, bigaussian)
794        ! cloud properties calculation with the initial values of t and q
795        ! ----------------------------------------------------------------
796
797        ! initialise gammasat and qincloud_mpc
798        gammasat(:)=1.
799        qincloud_mpc(:)=0.
800
801        IF (iflag_cld_th.GE.5) THEN
802               ! Cloud cover and content in meshes affected by shallow convection,
803               ! are retrieved from a bi-gaussian distribution of the saturation deficit
804               ! following Jam et al. 2013
805
806               IF (iflag_cloudth_vert.LE.2) THEN
807                  ! Old version of Arnaud Jam
808
809                    CALL cloudth(klon,klev,k,tv,                  &
810                         zq,qta,fraca,                            &
811                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
812                         ratqs,zqs,temp,                              &
813                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
814
815
816                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
817                   ! Default version of Arnaud Jam
818                       
819                    CALL cloudth_v3(klon,klev,k,tv,                        &
820                         zq,qta,fraca,                                     &
821                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
822                         ratqs,sigma_qtherm,zqs,temp, &
823                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
824
825
826                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
827                   ! Jean Jouhaud's version, with specific separation between surface and volume
828                   ! cloud fraction Decembre 2018
829
830                    CALL cloudth_v6(klon,klev,k,tv,                        &
831                         zq,qta,fraca,                                     &
832                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
833                         ratqs,zqs,temp, &
834                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
835
836                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
837                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
838                   ! for boundary-layer mixed phase clouds
839                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
840                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
841                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
842                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
843
844                ENDIF
845
846
847                DO i=1,klon
848                    rneb(i,k)=ctot(i,k)
849                    rneblsvol(i,k)=ctot_vol(i,k)
850                    zqn(i)=qcloud(i)
851                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
852                    qvc(i) = rneb(i,k) * zqs(i)
853                ENDDO
854
855        ENDIF
856
857        IF (iflag_cld_th .LE. 4) THEN
858           
859                ! lognormal
860            lognormale(:) = .TRUE.
861
862        ELSEIF (iflag_cld_th .GE. 6) THEN
863           
864                ! lognormal distribution when no thermals
865            lognormale(:) = fraca(:,k) < min_frac_th_cld
866
867        ELSE
868                ! When iflag_cld_th=5, we always assume
869                ! bi-gaussian distribution
870            lognormale(:) = .FALSE.
871       
872        ENDIF
873
874        DT(:) = 0.
875        n_i(:)=0
876        Tbef(:)=zt(:)
877        qlbef(:)=0.
878
879        ! Treatment of non-boundary layer clouds (lognormale)
880        ! We iterate here to take into account the change in qsat(T) and ice fraction
881        ! during the condensation process
882        ! the increment in temperature is calculated using the first principle of
883        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
884        ! and a clear sky fraction)
885        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
886
887        DO iter=1,iflag_fisrtilp_qsat+1
888
889                keepgoing(:) = .FALSE.
890
891                DO i=1,klon
892
893                    ! keepgoing = .true. while convergence is not satisfied
894
895                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
896
897                        ! if not convergence:
898                        ! we calculate a new iteration
899                        keepgoing(i) = .TRUE.
900
901                        ! P2.2.1> cloud fraction and condensed water mass calculation
902                        ! Calculated variables:
903                        ! rneb : cloud fraction
904                        ! zqn : total water within the cloud
905                        ! zcond: mean condensed water within the mesh
906                        ! rhcl: clear-sky relative humidity
907                        !---------------------------------------------------------------
908
909                        ! new temperature that only serves in the iteration process:
910                        Tbef(i)=Tbef(i)+DT(i)
911
912                        ! Rneb, qzn and zcond for lognormal PDFs
913                        qtot(i)=zq(i)+zmqc(i)
914
915                      ENDIF
916
917                  ENDDO
918
919                  ! Calculation of saturation specific humidity and ice fraction
920                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
921                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
922                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
923                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
924                  ! cloud phase determination
925                  IF (iflag_t_glace.GE.4) THEN
926                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
927                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
928                  ENDIF
929
930                  CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
931
932                  !--AB Activates a condensation scheme that allows for
933                  !--ice supersaturation and contrails evolution from aviation
934                  IF (ok_ice_supersat) THEN
935
936                    !--Calculate the shear value (input for condensation and ice supersat)
937                    DO i = 1, klon
938                      !--Cell thickness [m]
939                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
940                      IF ( iftop ) THEN
941                        ! top
942                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
943                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
944                      ELSEIF ( k .EQ. 1 ) THEN
945                        ! surface
946                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
947                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
948                      ELSE
949                        ! other layers
950                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
951                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
952                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
953                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
954                      ENDIF
955                    ENDDO
956
957                    !---------------------------------------------
958                    !--   CONDENSATION AND ICE SUPERSATURATION  --
959                    !---------------------------------------------
960
961                    CALL condensation_ice_supersat( &
962                        klon, dtime, missing_val, &
963                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
964                        cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
965                        shear(:), tke_dissip(:,k), cell_area(:), &
966                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
967                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
968                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
969                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
970                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
971                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
972                        flight_dist(:,k), flight_h2o(:,k), &
973                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
974
975
976                  ELSE
977                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
978
979                   CALL condensation_lognormal( &
980                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
981                       keepgoing, rneb(:,k), zqn, qvc)
982
983
984                  ENDIF ! .NOT. ok_ice_supersat
985
986                  DO i=1,klon
987                      IF (keepgoing(i)) THEN
988
989                        ! If vertical heterogeneity, change fraction by volume as well
990                        IF (iflag_cloudth_vert.GE.3) THEN
991                            ctot_vol(i,k)=rneb(i,k)
992                            rneblsvol(i,k)=ctot_vol(i,k)
993                        ENDIF
994
995
996                       ! P2.2.2> Approximative calculation of temperature variation DT
997                       !           due to condensation.
998                       ! Calculated variables:
999                       ! dT : temperature change due to condensation
1000                       !---------------------------------------------------------------
1001
1002               
1003                        IF (zfice(i).LT.1) THEN
1004                            cste=RLVTT
1005                        ELSE
1006                            cste=RLSTT
1007                        ENDIF
1008                       
1009                        ! LEA_R : check formule
1010                        IF ( ok_unadjusted_clouds ) THEN
1011                          !--AB We relax the saturation adjustment assumption
1012                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1013                          IF ( rneb(i,k) .GT. eps ) THEN
1014                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
1015                          ELSE
1016                            qlbef(i) = 0.
1017                          ENDIF
1018                        ELSE
1019                          qlbef(i)=max(0.,zqn(i)-zqs(i))
1020                        ENDIF
1021
1022                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1023                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
1024                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1025                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
1026                              *qlbef(i)*dzfice(i)
1027                        ! here we update a provisory temperature variable that only serves in the iteration
1028                        ! process
1029                        DT(i)=num/denom
1030                        n_i(i)=n_i(i)+1
1031
1032                    ENDIF ! end keepgoing
1033
1034                ENDDO     ! end loop on i
1035
1036        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
1037
1038        ! P2.3> Final quantities calculation
1039        ! Calculated variables:
1040        !   rneb : cloud fraction
1041        !   zcond: mean condensed water in the mesh
1042        !   zqn  : mean water vapor in the mesh
1043        !   zfice: ice fraction in clouds
1044        !   zt   : temperature
1045        !   rhcl : clear-sky relative humidity
1046        ! ----------------------------------------------------------------
1047
1048
1049        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
1050        IF (iflag_t_glace.GE.4) THEN
1051           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
1052           distcltop(:,k)=zdistcltop(:)
1053           temp_cltop(:,k)=ztemp_cltop(:)
1054        ENDIF
1055
1056        ! Partition function depending on temperature
1057        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
1058
1059        ! Partition function depending on tke for non shallow-convective clouds
1060        IF (iflag_icefrac .GE. 1) THEN
1061           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
1062           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
1063        ENDIF
1064
1065        ! Water vapor update, Phase determination and subsequent latent heat exchange
1066        DO i=1, klon
1067            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
1068            ! iflag_cloudth_vert=7 and specific param is activated
1069            IF (mpc_bl_points(i,k) .GT. 0) THEN
1070                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
1071                ! following line is very strange and probably wrong
1072                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
1073                ! water vapor update and partition function if thermals
1074                zq(i) = zq(i) - zcond(i)       
1075                zfice(i)=zfice_th(i)
1076            ELSE
1077                ! Checks on rneb, rhcl and zqn
1078                IF (rneb(i,k) .LE. 0.0) THEN
1079                    zqn(i) = 0.0
1080                    rneb(i,k) = 0.0
1081                    zcond(i) = 0.0
1082                    rhcl(i,k)=zq(i)/zqs(i)
1083                ELSE IF (rneb(i,k) .GE. 1.0) THEN
1084                    zqn(i) = zq(i)
1085                    rneb(i,k) = 1.0
1086                    IF ( ok_unadjusted_clouds ) THEN
1087                      !--AB We relax the saturation adjustment assumption
1088                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1089                      zcond(i) = MAX(0., zqn(i) - qvc(i))
1090                    ELSE
1091                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1092                    ENDIF
1093                    rhcl(i,k)=1.0
1094                ELSE
1095                    IF ( ok_unadjusted_clouds ) THEN
1096                      !--AB We relax the saturation adjustment assumption
1097                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1098                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
1099                    ELSE
1100                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1101                    ENDIF
1102                    ! following line is very strange and probably wrong:
1103                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
1104                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
1105                    IF (iflag_icefrac .GE. 1) THEN
1106                        IF (lognormale(i)) THEN 
1107                           zcond(i)  = zqliq(i) + zqice(i)
1108                           zfice(i)=zfice_turb(i)
1109                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
1110                        ENDIF
1111                    ENDIF
1112                ENDIF
1113
1114                ! water vapor update
1115                zq(i) = zq(i) - zcond(i)
1116
1117            ENDIF
1118
1119            ! c_iso : routine that computes in-cloud supersaturation
1120            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
1121
1122            ! temperature update due to phase change
1123            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1124            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1125                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1126        ENDDO
1127
1128        ! If vertical heterogeneity, change volume fraction
1129        IF (iflag_cloudth_vert .GE. 3) THEN
1130          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1131          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1132        ENDIF
1133   
1134    !--AB Write diagnostics and tracers for ice supersaturation
1135    IF ( ok_ice_supersat ) THEN
1136      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1137      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1138
1139      DO i = 1, klon
1140
1141        cf_seri(i,k) = rneb(i,k)
1142
1143        IF ( .NOT. ok_unadjusted_clouds ) THEN
1144          qvc(i) = zqs(i) * rneb(i,k)
1145        ENDIF
1146        IF ( zq(i) .GT. min_qParent ) THEN
1147          rvc_seri(i,k) = qvc(i) / zq(i)
1148        ELSE
1149          rvc_seri(i,k) = min_ratio
1150        ENDIF
1151        !--The MIN barrier is NEEDED because of:
1152        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1153        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1154        !--The MAX barrier is a safeguard that should not be activated
1155        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1156
1157        !--Diagnostics
1158        gamma_cond(i,k) = gammasat(i)
1159        IF ( issrfra(i,k) .LT. eps ) THEN
1160          issrfra(i,k) = 0.
1161          qissr(i,k) = 0.
1162        ENDIF
1163        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1164        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1165        IF ( subfra(i,k) .LT. eps ) THEN
1166          subfra(i,k) = 0.
1167          qsub(i,k) = 0.
1168        ENDIF
1169        qcld(i,k) = qvc(i) + zcond(i)
1170        IF ( cf_seri(i,k) .LT. eps ) THEN
1171          qcld(i,k) = 0.
1172        ENDIF
1173      ENDDO
1174    ENDIF
1175
1176
1177    ! ----------------------------------------------------------------
1178    ! P3> Precipitation formation
1179    ! ----------------------------------------------------------------
1180
1181    !================================================================
1182    IF (ok_poprecip) THEN
1183
1184      DO i = 1, klon
1185        zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1186        zoliqi(i) = zcond(i) * zfice(i)
1187      ENDDO
1188
1189      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1190                            ctot_vol(:,k), ptconv(:,k), &
1191                            zt, zq, zoliql, zoliqi, zfice, &
1192                            rneb(:,k), znebprecipclr, znebprecipcld, &
1193                            zrfl, zrflclr, zrflcld, &
1194                            zifl, ziflclr, ziflcld, &
1195                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1196                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1197                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1198                            dqsmelt(:,k), dqsfreez(:,k) &
1199                            )
1200
1201      DO i = 1, klon
1202        zoliq(i) = zoliql(i) + zoliqi(i)
1203        IF ( zoliq(i) .GT. 0. ) THEN
1204                zfice(i) = zoliqi(i)/zoliq(i)
1205        ELSE 
1206                zfice(i) = 0.0
1207        ENDIF
1208
1209        ! calculation of specific content of condensates seen by radiative scheme
1210        IF (ok_radocond_snow) THEN
1211           radocond(i,k) = zoliq(i)
1212           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1213           radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
1214        ELSE
1215           radocond(i,k) = zoliq(i)
1216           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
1217           radocondi(i,k)= radocond(i,k)*zfice(i)
1218        ENDIF
1219      ENDDO
1220
1221    !================================================================
1222    ELSE
1223
1224    ! LTP:
1225    IF (iflag_evap_prec .GE. 4) THEN
1226
1227        !Partitionning between precipitation coming from clouds and that coming from CS
1228
1229        !0) Calculate tot_zneb, total cloud fraction above the cloud
1230        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
1231       
1232        DO i=1, klon
1233                tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
1234                        /(1.-min(zneb(i),1.-smallestreal))
1235                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
1236                tot_zneb(i) = tot_znebn(i)
1237
1238
1239                !1) Cloudy to clear air
1240                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
1241                IF (znebprecipcld(i) .GT. 0.) THEN
1242                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
1243                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
1244                ELSE
1245                        d_zrfl_cld_clr(i) = 0.
1246                        d_zifl_cld_clr(i) = 0.
1247                ENDIF
1248
1249                !2) Clear to cloudy air
1250                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
1251                IF (znebprecipclr(i) .GT. 0) THEN
1252                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
1253                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
1254                ELSE
1255                        d_zrfl_clr_cld(i) = 0.
1256                        d_zifl_clr_cld(i) = 0.
1257                ENDIF
1258
1259                !Update variables
1260                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
1261                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
1262                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
1263                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
1264                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
1265                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
1266
1267                ! c_iso  do the same thing for isotopes precip
1268        ENDDO
1269    ENDIF
1270
1271
1272    ! Autoconversion
1273    ! -------------------------------------------------------------------------------
1274
1275
1276    ! Initialisation of zoliq and radocond variables
1277
1278    DO i = 1, klon
1279            zoliq(i) = zcond(i)
1280            zoliqi(i)= zoliq(i)*zfice(i)
1281            zoliql(i)= zoliq(i)*(1.-zfice(i))
1282            ! c_iso : initialisation of zoliq* also for isotopes
1283            zrho(i,k)  = pplay(i,k) / zt(i) / RD
1284            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
1285            iwc(i)   = 0.
1286            zneb(i)  = MAX(rneb(i,k),seuil_neb)
1287            radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
1288            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
1289            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
1290    ENDDO
1291
1292       
1293    DO n = 1, niter_lscp
1294
1295        DO i=1,klon
1296            IF (rneb(i,k).GT.0.0) THEN
1297                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
1298            ENDIF
1299        ENDDO
1300
1301        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
1302
1303        DO i = 1, klon
1304
1305            IF (rneb(i,k).GT.0.0) THEN
1306
1307                ! Initialization of zrain, zsnow and zprecip:
1308                zrain=0.
1309                zsnow=0.
1310                zprecip=0.
1311                ! c_iso same init for isotopes. Externalisation?
1312
1313                IF (zneb(i).EQ.seuil_neb) THEN
1314                    zprecip = 0.0
1315                    zsnow = 0.0
1316                    zrain= 0.0
1317                ELSE
1318
1319                    IF (ptconv(i,k)) THEN ! if convective point
1320                        zcl=cld_lc_con
1321                        zct=1./cld_tau_con
1322                        zexpo=cld_expo_con
1323                        ffallv=ffallv_con
1324                    ELSE
1325                        zcl=cld_lc_lsc
1326                        zct=1./cld_tau_lsc
1327                        zexpo=cld_expo_lsc
1328                        ffallv=ffallv_lsc
1329                    ENDIF
1330
1331
1332                    ! if vertical heterogeneity is taken into account, we use
1333                    ! the "true" volume fraction instead of a modified
1334                    ! surface fraction (which is larger and artificially
1335                    ! reduces the in-cloud water).
1336
1337                    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
1338                    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
1339                    !.........................................................
1340                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
1341
1342                    ! if vertical heterogeneity is taken into account, we use
1343                    ! the "true" volume fraction instead of a modified
1344                    ! surface fraction (which is larger and artificially
1345                    ! reduces the in-cloud water).
1346                        effective_zneb=ctot_vol(i,k)
1347                    ELSE
1348                        effective_zneb=zneb(i)
1349                    ENDIF
1350
1351
1352                    IF (iflag_autoconversion .EQ. 2) THEN
1353                    ! two-steps resolution with niter_lscp=1 sufficient
1354                    ! we first treat the second term (with exponential) in an explicit way
1355                    ! and then treat the first term (-q/tau) in an exact way
1356                        zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
1357                        *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
1358                    ELSE
1359                    ! old explicit resolution with subtimesteps
1360                        zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
1361                        *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
1362                    ENDIF
1363
1364
1365                    ! Ice water quantity to remove (Zender & Kiehl, 1997)
1366                    ! dqice/dt=1/rho*d(rho*wice*qice)/dz
1367                    !....................................
1368                    IF (iflag_autoconversion .EQ. 2) THEN
1369                    ! exact resolution, niter_lscp=1 is sufficient but works only
1370                    ! with iflag_vice=0
1371                       IF (zoliqi(i) .GT. 0.) THEN
1372                          zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
1373                          +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
1374                       ELSE
1375                          zfroi=0.
1376                       ENDIF
1377                    ELSE
1378                    ! old explicit resolution with subtimesteps
1379                       zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k)
1380                    ENDIF
1381
1382                    zrain   = MIN(MAX(zchau,0.0),zoliql(i))
1383                    zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
1384                    zprecip = MAX(zrain + zsnow,0.0)
1385
1386                ENDIF
1387
1388
1389                IF (iflag_autoconversion .GE. 1) THEN
1390                   ! debugged version with phase conservation through the autoconversion process
1391                   zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
1392                   zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
1393                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1394                ELSE
1395                   ! bugged version with phase resetting
1396                   zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
1397                   zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
1398                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1399                ENDIF
1400
1401                !--Diagnostics
1402                dqrauto(i,k) = dqrauto(i,k) + zrain / dtime
1403                dqsauto(i,k) = dqsauto(i,k) + zsnow / dtime
1404
1405
1406                ! c_iso: call isotope_conversion (for readibility) or duplicate
1407
1408                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
1409                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
1410                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
1411
1412            ENDIF ! rneb >0
1413
1414        ENDDO  ! i = 1,klon
1415
1416    ENDDO      ! n = 1,niter
1417
1418    ! Precipitation flux calculation
1419
1420    DO i = 1, klon
1421           
1422            IF (iflag_evap_prec.GE.4) THEN
1423                ziflprev(i)=ziflcld(i)
1424            ELSE
1425                ziflprev(i)=zifl(i)*zneb(i)
1426            ENDIF
1427
1428            IF (rneb(i,k) .GT. 0.0) THEN
1429
1430               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1431               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1432               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1433               ! taken into account through a linearization of the equations and by approximating
1434               ! the liquid precipitation process with a "threshold" process. We assume that
1435               ! condensates are not modified during this operation. Liquid precipitation is
1436               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1437               ! Water vapor increases as well
1438               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1439
1440                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1441                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1442                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1443                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1444                    ! Computation of DT if all the liquid precip freezes
1445                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1446                    ! T should not exceed the freezing point
1447                    ! that is Delta > RTT-zt(i)
1448                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1449                    zt(i)  = zt(i) + DeltaT
1450                    ! water vaporization due to temp. increase
1451                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1452                    ! we add this vaporized water to the total vapor and we remove it from the precip
1453                    zq(i) = zq(i) +  Deltaq
1454                    ! The three "max" lines herebelow protect from rounding errors
1455                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1456                    ! liquid precipitation converted to ice precip
1457                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1458                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1459                    ! iced water budget
1460                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1461
1462                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1463
1464                IF (iflag_evap_prec.GE.4) THEN
1465                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1466                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1467                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1468                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1469                    znebprecipcld(i) = rneb(i,k)
1470                    zrfl(i) = zrflcld(i) + zrflclr(i)
1471                    zifl(i) = ziflcld(i) + ziflclr(i)
1472                ELSE
1473                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1474                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1475                    zifl(i) = zifl(i)+ zqpreci(i)        &
1476                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1477                ENDIF
1478                ! c_iso : same for isotopes
1479 
1480           ENDIF ! rneb>0
1481
1482   ENDDO
1483
1484    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1485    ! if iflag_evap_prec>=4
1486    IF (iflag_evap_prec.GE.4) THEN
1487
1488        DO i=1,klon
1489
1490            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1491                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1492                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1493            ELSE
1494                znebprecipclr(i)=0.0
1495            ENDIF
1496
1497            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1498                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1499                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1500            ELSE
1501                znebprecipcld(i)=0.0
1502            ENDIF
1503        ENDDO
1504
1505    ENDIF
1506
1507
1508    ENDIF ! ok_poprecip
1509
1510    ! End of precipitation formation
1511    ! --------------------------------
1512
1513
1514    ! Calculation of cloud condensates amount seen by radiative scheme
1515    !-----------------------------------------------------------------
1516
1517    ! Calculation of the concentration of condensates seen by the radiation scheme
1518    ! for liquid phase, we take radocondl
1519    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
1520    ! we recalculate radocondi to account for contributions from the precipitation flux
1521    ! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
1522
1523    IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
1524        ! for the solid phase (crystals + snowflakes)
1525        ! we recalculate radocondi by summing
1526        ! the ice content calculated in the mesh
1527        ! + the contribution of the non-evaporated snowfall
1528        ! from the overlying layer
1529        DO i=1,klon
1530           IF (ziflprev(i) .NE. 0.0) THEN
1531              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1532           ELSE
1533              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1534           ENDIF
1535           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
1536        ENDDO
1537    ENDIF
1538
1539    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
1540    DO i=1,klon
1541        IF (radocond(i,k) .GT. 0.) THEN
1542            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
1543        ENDIF
1544    ENDDO
1545
1546    ! ----------------------------------------------------------------
1547    ! P4> Wet scavenging
1548    ! ----------------------------------------------------------------
1549
1550    !Scavenging through nucleation in the layer
1551
1552    DO i = 1,klon
1553       
1554        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1555            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1556        ELSE
1557            beta(i,k) = 0.
1558        ENDIF
1559
1560        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1561
1562        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1563
1564            IF (temp(i,k) .GE. t_glace_min) THEN
1565                zalpha_tr = a_tr_sca(3)
1566            ELSE
1567                zalpha_tr = a_tr_sca(4)
1568            ENDIF
1569
1570            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1571            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1572
1573            ! Nucleation with a factor of -1 instead of -0.5
1574            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1575
1576        ENDIF
1577
1578    ENDDO
1579
1580    ! Scavenging through impaction in the underlying layer
1581
1582    DO kk = k-1, 1, -1
1583
1584        DO i = 1, klon
1585
1586            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1587
1588                IF (temp(i,kk) .GE. t_glace_min) THEN
1589                    zalpha_tr = a_tr_sca(1)
1590                ELSE
1591                    zalpha_tr = a_tr_sca(2)
1592                ENDIF
1593
1594              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1595              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1596
1597            ENDIF
1598
1599        ENDDO
1600
1601    ENDDO
1602
1603    ! Outputs:
1604    !-------------------------------
1605    ! Precipitation fluxes at layer interfaces
1606    ! + precipitation fractions +
1607    ! temperature and water species tendencies
1608    DO i = 1, klon
1609        psfl(i,k)=zifl(i)
1610        prfl(i,k)=zrfl(i)
1611        pfraclr(i,k)=znebprecipclr(i)
1612        pfracld(i,k)=znebprecipcld(i)
1613        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1614        d_qi(i,k) = zfice(i)*zoliq(i)
1615        d_q(i,k) = zq(i) - qt(i,k)
1616        ! c_iso: same for isotopes
1617        d_t(i,k) = zt(i) - temp(i,k)
1618    ENDDO
1619
1620
1621ENDDO
1622
1623
1624  ! Rain or snow at the surface (depending on the first layer temperature)
1625  DO i = 1, klon
1626      snow(i) = zifl(i)
1627      rain(i) = zrfl(i)
1628      ! c_iso final output
1629  ENDDO
1630
1631  IF (ncoreczq>0) THEN
1632      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1633  ENDIF
1634
1635
1636RETURN
1637
1638END SUBROUTINE lscp
1639!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1640
1641END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.