source: LMDZ6/branches/ICOLMDZISO/libf/phylmd/lmdz_lscp.F90 @ 5875

Last change on this file since 5875 was 5211, checked in by evignon, 14 months ago

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