source: LMDZ6/branches/Amaury_posttrusting/libf/phylmd/lmdz_lscp.f90 @ 5356

Last change on this file since 5356 was 5268, checked in by abarral, 13 months ago

.f90 <-> .F90 depending on cpp key use
File size: 71.4 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
108!USE lmdz_lscp_condensation, ONLY : condensation_lognormal_test, condensation_ice_supersat
109USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
110USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
111
112! Use du module lmdz_lscp_ini contenant les constantes
113USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
114USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
115USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
116USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
117USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
118USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
119USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
120USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
121USE lmdz_lscp_ini, ONLY : ok_poprecip
122USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
123
124IMPLICIT NONE
125
126!===============================================================================
127! VARIABLES DECLARATION
128!===============================================================================
129
130! INPUT VARIABLES:
131!-----------------
132
133  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
134  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
135  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
136
137  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
139  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
140  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
141  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [kg/kg]
142  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
143  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
144                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
145  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
146
147  !Inputs associated with thermal plumes
148
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
153  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
154  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
155  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
156
157  ! INPUT/OUTPUT variables
158  !------------------------
159 
160  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
161  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
162
163
164  ! INPUT/OUTPUT condensation and ice supersaturation
165  !--------------------------------------------------
166  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
167  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratio_qi_qtot    ! solid specific water to total specific water ratio [-]
168  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
170  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
171  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
172
173  ! INPUT/OUTPUT aviation
174  !--------------------------------------------------
175  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
176  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
177 
178  ! OUTPUT variables
179  !-----------------
180
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
186  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
187  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
188  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
191  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
192  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
193  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
194  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
195  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
196  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
197  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
198  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
201  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
202
203  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
204 
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
206  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
207 
208  ! for condensation and ice supersaturation
209
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
228  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
229
230  ! for contrails and aviation
231
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
239  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
240
241
242  ! for POPRECIP
243
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
256  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
257
258  ! for thermals
259
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
263  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
264
265
266  ! LOCAL VARIABLES:
267  !----------------
268  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
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    ! --------------------------------------------------------------------
572    ! P2> Precipitation evaporation/sublimation/melting
573    ! --------------------------------------------------------------------
574    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
575    ! --------------------------------------------------------------------
576   
577    IF (iflag_evap_prec.GE.1) THEN
578
579        ! Calculation of saturation specific humidity
580        ! depending on temperature:
581        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
582        ! wrt liquid water
583        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
584        ! wrt ice
585        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
586
587        DO i = 1, klon
588
589            ! if precipitation
590            IF (zrfl(i)+zifl(i).GT.0.) THEN
591
592            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
593            ! c_iso: likely important to distinguish cs from neb isotope precipitation
594
595                IF (iflag_evap_prec.GE.4) THEN
596                    zrfl(i) = zrflclr(i)
597                    zifl(i) = ziflclr(i)
598                ENDIF
599
600                IF (iflag_evap_prec.EQ.1) THEN
601                    znebprecip(i)=zneb(i)
602                ELSE
603                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
604                ENDIF
605
606                IF (iflag_evap_prec.GT.4) THEN
607                ! Max evaporation not to saturate the clear sky precip fraction
608                ! i.e. the fraction where evaporation occurs
609                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
610                ELSEIF (iflag_evap_prec .EQ. 4) THEN
611                ! Max evaporation not to saturate the whole mesh
612                ! Pay attention -> lead to unrealistic and excessive evaporation
613                    zqev0 = MAX(0.0, zqs(i)-zq(i))
614                ELSE
615                ! Max evap not to saturate the fraction below the cloud
616                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
617                ENDIF
618
619                ! Evaporation of liquid precipitation coming from above
620                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
621                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
622                ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
623
624                IF (iflag_evap_prec.EQ.3) THEN
625                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
626                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
627                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
628                ELSE IF (iflag_evap_prec.GE.4) THEN
629                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
630                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
631                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
632                ELSE
633                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
634                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
635                ENDIF
636
637                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
638                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
639
640                ! sublimation of the solid precipitation coming from above
641                IF (iflag_evap_prec.EQ.3) THEN
642                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
643                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
644                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
645                ELSE IF (iflag_evap_prec.GE.4) THEN
646                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
647                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
648                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
649                ELSE
650                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
651                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
652                ENDIF
653
654                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
655                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
656
657                ! A. JAM
658                ! Evaporation limit: we ensure that the layer's fraction below
659                ! the cloud or the whole mesh (depending on iflag_evap_prec)
660                ! does not reach saturation. In this case, we
661                ! redistribute zqev0 conserving the ratio liquid/ice
662
663                IF (zqevt+zqevti.GT.zqev0) THEN
664                    zqev=zqev0*zqevt/(zqevt+zqevti)
665                    zqevi=zqev0*zqevti/(zqevt+zqevti)
666                ELSE
667                    zqev=zqevt
668                    zqevi=zqevti
669                ENDIF
670
671
672                ! New solid and liquid precipitation fluxes after evap and sublimation
673                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
674                         /RG/dtime)
675                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
676                         /RG/dtime)
677
678
679                ! vapor, temperature, precip fluxes update
680                ! vapor is updated after evaporation/sublimation (it is increased)
681                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
682                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
683                ! zmqc is the total condensed water in the precip flux (it is decreased)
684                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
685                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
686                ! air and precip temperature (i.e., gridbox temperature)
687                ! is updated due to latent heat cooling
688                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
689                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
690                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
691                + (zifln(i)-zifl(i))                      &
692                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
693                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
694
695                ! New values of liquid and solid precipitation
696                zrfl(i) = zrfln(i)
697                zifl(i) = zifln(i)
698
699                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
700                !                         due to evap + sublim                                     
701                                           
702
703                IF (iflag_evap_prec.GE.4) THEN
704                    zrflclr(i) = zrfl(i)
705                    ziflclr(i) = zifl(i)
706                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
707                        znebprecipclr(i) = 0.0
708                    ENDIF
709                    zrfl(i) = zrflclr(i) + zrflcld(i)
710                    zifl(i) = ziflclr(i) + ziflcld(i)
711                ENDIF
712
713                ! c_iso duplicate for isotopes or loop on isotopes
714
715                ! Melting:
716                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
717                ! precip fraction that is melted
718                zmelt = MIN(MAX(zmelt,0.),1.)
719
720                ! update of rainfall and snowfall due to melting
721                IF (iflag_evap_prec.GE.4) THEN
722                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
723                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
724                    zrfl(i)=zrflclr(i)+zrflcld(i)
725                ELSE
726                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
727                ENDIF
728
729
730                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
731
732                ! Latent heat of melting because of precipitation melting
733                ! NB: the air + precip temperature is simultaneously updated
734                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
735                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
736               
737                IF (iflag_evap_prec.GE.4) THEN
738                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
739                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
740                    zifl(i)=ziflclr(i)+ziflcld(i)
741                ELSE
742                    zifl(i)=zifl(i)*(1.-zmelt)
743                ENDIF
744
745            ELSE
746                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
747                znebprecip(i)=0.
748
749            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
750
751        ENDDO ! loop on klon
752
753    ENDIF ! (iflag_evap_prec>=1)
754
755    ENDIF ! (ok_poprecip)
756
757    ! --------------------------------------------------------------------
758    ! End precip evaporation
759    ! --------------------------------------------------------------------
760   
761    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
762    !-------------------------------------------------------
763
764     qtot(:)=zq(:)+zmqc(:)
765     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
766     DO i = 1, klon
767            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
768            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
769            IF (zq(i) .LT. 1.e-15) THEN
770                ncoreczq=ncoreczq+1
771                zq(i)=1.e-15
772            ENDIF
773            ! c_iso: do something similar for isotopes
774
775     ENDDO
776   
777    ! --------------------------------------------------------------------
778    ! P2> Cloud formation
779    !---------------------------------------------------------------------
780    !
781    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
782    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
783    ! is prescribed and depends on large scale variables and boundary layer
784    ! properties)
785    ! The decrease in condensed part due tu latent heating is taken into
786    ! account
787    ! -------------------------------------------------------------------
788
789        ! P2.1> With the PDFs (log-normal, bigaussian)
790        ! cloud properties calculation with the initial values of t and q
791        ! ----------------------------------------------------------------
792
793        ! initialise gammasat and qincloud_mpc
794        gammasat(:)=1.
795        qincloud_mpc(:)=0.
796
797        IF (iflag_cld_th.GE.5) THEN
798               ! Cloud cover and content in meshes affected by shallow convection,
799               ! are retrieved from a bi-gaussian distribution of the saturation deficit
800               ! following Jam et al. 2013
801
802               IF (iflag_cloudth_vert.LE.2) THEN
803                  ! Old version of Arnaud Jam
804
805                    CALL cloudth(klon,klev,k,tv,                  &
806                         zq,qta,fraca,                            &
807                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
808                         ratqs,zqs,temp,                              &
809                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
810
811
812                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
813                   ! Default version of Arnaud Jam
814                       
815                    CALL cloudth_v3(klon,klev,k,tv,                        &
816                         zq,qta,fraca,                                     &
817                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
818                         ratqs,sigma_qtherm,zqs,temp, &
819                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
820
821
822                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
823                   ! Jean Jouhaud's version, with specific separation between surface and volume
824                   ! cloud fraction Decembre 2018
825
826                    CALL cloudth_v6(klon,klev,k,tv,                        &
827                         zq,qta,fraca,                                     &
828                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
829                         ratqs,zqs,temp, &
830                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
831
832                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
833                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
834                   ! for boundary-layer mixed phase clouds
835                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
836                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
837                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
838                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
839
840                ENDIF
841
842
843                DO i=1,klon
844                    rneb(i,k)=ctot(i,k)
845                    rneblsvol(i,k)=ctot_vol(i,k)
846                    zqn(i)=qcloud(i)
847                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
848                    qvc(i) = rneb(i,k) * zqs(i)
849                ENDDO
850
851        ENDIF
852
853        IF (iflag_cld_th .LE. 4) THEN
854           
855                ! lognormal
856            lognormale(:) = .TRUE.
857
858        ELSEIF (iflag_cld_th .GE. 6) THEN
859           
860                ! lognormal distribution when no thermals
861            lognormale(:) = fraca(:,k) < min_frac_th_cld
862
863        ELSE
864                ! When iflag_cld_th=5, we always assume
865                ! bi-gaussian distribution
866            lognormale(:) = .FALSE.
867       
868        ENDIF
869
870        DT(:) = 0.
871        n_i(:)=0
872        Tbef(:)=zt(:)
873        qlbef(:)=0.
874
875        ! Treatment of non-boundary layer clouds (lognormale)
876        ! condensation with qsat(T) variation (adaptation)
877        ! Iterative resolution to converge towards qsat
878        ! with update of temperature, ice fraction and qsat at
879        ! each iteration
880
881        ! todoan -> sensitivity to iflag_fisrtilp_qsat
882        DO iter=1,iflag_fisrtilp_qsat+1
883
884                keepgoing(:) = .FALSE.
885
886                DO i=1,klon
887
888                    ! keepgoing = .true. while convergence is not satisfied
889
890                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
891
892                        ! if not convergence:
893                        ! we calculate a new iteration
894                        keepgoing(i) = .TRUE.
895
896                        ! P2.2.1> cloud fraction and condensed water mass calculation
897                        ! Calculated variables:
898                        ! rneb : cloud fraction
899                        ! zqn : total water within the cloud
900                        ! zcond: mean condensed water within the mesh
901                        ! rhcl: clear-sky relative humidity
902                        !---------------------------------------------------------------
903
904                        ! new temperature that only serves in the iteration process:
905                        Tbef(i)=Tbef(i)+DT(i)
906
907                        ! Rneb, qzn and zcond for lognormal PDFs
908                        qtot(i)=zq(i)+zmqc(i)
909
910                      ENDIF
911
912                  ENDDO
913
914                  ! Calculation of saturation specific humidity and ice fraction
915                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
916                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
917                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
918                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
919                  ! cloud phase determination
920                  IF (iflag_t_glace.GE.4) THEN
921                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
922                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
923                  ENDIF
924
925                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
926
927                  !--AB Activates a condensation scheme that allows for
928                  !--ice supersaturation and contrails evolution from aviation
929                  IF (ok_ice_supersat) THEN
930
931                    !--Calculate the shear value (input for condensation and ice supersat)
932                    DO i = 1, klon
933                      !--Cell thickness [m]
934                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
935                      IF ( iftop ) THEN
936                        ! top
937                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
938                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
939                      ELSEIF ( k .EQ. 1 ) THEN
940                        ! surface
941                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
942                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
943                      ELSE
944                        ! other layers
945                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
946                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
947                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
948                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
949                      ENDIF
950                    ENDDO
951
952                    !---------------------------------------------
953                    !--   CONDENSATION AND ICE SUPERSATURATION  --
954                    !---------------------------------------------
955
956                    CALL condensation_ice_supersat( &
957                        klon, dtime, missing_val, &
958                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
959                        cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
960                        shear(:), tke_dissip(:,k), cell_area(:), &
961                        Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
962                        rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
963                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
964                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
965                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
966                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
967                        flight_dist(:,k), flight_h2o(:,k), &
968                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
969
970
971                  ELSE
972                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
973
974                   CALL condensation_lognormal( &
975                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
976                       keepgoing, rneb(:,k), zqn, qvc)
977!                   CALL condensation_lognormal_test( &
978!                       klon, klev, k, keepgoing, ratqs, gammasat, zq, zqs, &
979!                       zqn, qvc, rneb)
980
981
982                  ENDIF ! .NOT. ok_ice_supersat
983
984                  DO i=1,klon
985                      IF (keepgoing(i)) THEN
986
987                        ! If vertical heterogeneity, change fraction by volume as well
988                        IF (iflag_cloudth_vert.GE.3) THEN
989                            ctot_vol(i,k)=rneb(i,k)
990                            rneblsvol(i,k)=ctot_vol(i,k)
991                        ENDIF
992
993
994                       ! P2.2.2> Approximative calculation of temperature variation DT
995                       !           due to condensation.
996                       ! Calculated variables:
997                       ! dT : temperature change due to condensation
998                       !---------------------------------------------------------------
999
1000               
1001                        IF (zfice(i).LT.1) THEN
1002                            cste=RLVTT
1003                        ELSE
1004                            cste=RLSTT
1005                        ENDIF
1006                       
1007                        ! LEA_R : check formule
1008                        IF ( ok_unadjusted_clouds ) THEN
1009                          !--AB We relax the saturation adjustment assumption
1010                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1011                          IF ( rneb(i,k) .GT. eps ) THEN
1012                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
1013                          ELSE
1014                            qlbef(i) = 0.
1015                          ENDIF
1016                        ELSE
1017                          qlbef(i)=max(0.,zqn(i)-zqs(i))
1018                        ENDIF
1019
1020                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1021                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
1022                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1023                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
1024                              *qlbef(i)*dzfice(i)
1025                        ! here we update a provisory temperature variable that only serves in the iteration
1026                        ! process
1027                        DT(i)=num/denom
1028                        n_i(i)=n_i(i)+1
1029
1030                    ENDIF ! end keepgoing
1031
1032                ENDDO     ! end loop on i
1033
1034        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
1035
1036        ! P2.3> Final quantities calculation
1037        ! Calculated variables:
1038        !   rneb : cloud fraction
1039        !   zcond: mean condensed water in the mesh
1040        !   zqn  : mean water vapor in the mesh
1041        !   zfice: ice fraction in clouds
1042        !   zt   : temperature
1043        !   rhcl : clear-sky relative humidity
1044        ! ----------------------------------------------------------------
1045
1046
1047        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
1048        IF (iflag_t_glace.GE.4) THEN
1049           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
1050           distcltop(:,k)=zdistcltop(:)
1051           temp_cltop(:,k)=ztemp_cltop(:)
1052        ENDIF
1053
1054        ! Partition function depending on temperature
1055        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
1056
1057        ! Partition function depending on tke for non shallow-convective clouds
1058        IF (iflag_icefrac .GE. 1) THEN
1059
1060           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
1061           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
1062
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                IF (iflag_autoconversion .GE. 1) THEN
1389                   ! debugged version with phase conservation through the autoconversion process
1390                   zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
1391                   zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
1392                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1393                ELSE
1394                   ! bugged version with phase resetting
1395                   zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
1396                   zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
1397                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1398                ENDIF
1399
1400                ! c_iso: call isotope_conversion (for readibility) or duplicate
1401
1402                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
1403                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
1404                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
1405
1406            ENDIF ! rneb >0
1407
1408        ENDDO  ! i = 1,klon
1409
1410    ENDDO      ! n = 1,niter
1411
1412    ! Precipitation flux calculation
1413
1414    DO i = 1, klon
1415           
1416            IF (iflag_evap_prec.GE.4) THEN
1417                ziflprev(i)=ziflcld(i)
1418            ELSE
1419                ziflprev(i)=zifl(i)*zneb(i)
1420            ENDIF
1421
1422            IF (rneb(i,k) .GT. 0.0) THEN
1423
1424               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1425               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1426               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1427               ! taken into account through a linearization of the equations and by approximating
1428               ! the liquid precipitation process with a "threshold" process. We assume that
1429               ! condensates are not modified during this operation. Liquid precipitation is
1430               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1431               ! Water vapor increases as well
1432               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1433
1434                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1435                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1436                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1437                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1438                    ! Computation of DT if all the liquid precip freezes
1439                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1440                    ! T should not exceed the freezing point
1441                    ! that is Delta > RTT-zt(i)
1442                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1443                    zt(i)  = zt(i) + DeltaT
1444                    ! water vaporization due to temp. increase
1445                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1446                    ! we add this vaporized water to the total vapor and we remove it from the precip
1447                    zq(i) = zq(i) +  Deltaq
1448                    ! The three "max" lines herebelow protect from rounding errors
1449                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1450                    ! liquid precipitation converted to ice precip
1451                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1452                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1453                    ! iced water budget
1454                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1455
1456                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1457
1458                IF (iflag_evap_prec.GE.4) THEN
1459                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1460                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1461                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1462                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1463                    znebprecipcld(i) = rneb(i,k)
1464                    zrfl(i) = zrflcld(i) + zrflclr(i)
1465                    zifl(i) = ziflcld(i) + ziflclr(i)
1466                ELSE
1467                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1468                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1469                    zifl(i) = zifl(i)+ zqpreci(i)        &
1470                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1471                ENDIF
1472                ! c_iso : same for isotopes
1473 
1474           ENDIF ! rneb>0
1475
1476   ENDDO
1477
1478    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1479    ! if iflag_evap_prec>=4
1480    IF (iflag_evap_prec.GE.4) THEN
1481
1482        DO i=1,klon
1483
1484            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1485                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1486                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1487            ELSE
1488                znebprecipclr(i)=0.0
1489            ENDIF
1490
1491            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1492                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1493                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1494            ELSE
1495                znebprecipcld(i)=0.0
1496            ENDIF
1497        ENDDO
1498
1499    ENDIF
1500
1501
1502    ENDIF ! ok_poprecip
1503
1504    ! End of precipitation formation
1505    ! --------------------------------
1506
1507
1508    ! Calculation of cloud condensates amount seen by radiative scheme
1509    !-----------------------------------------------------------------
1510
1511    ! Calculation of the concentration of condensates seen by the radiation scheme
1512    ! for liquid phase, we take radocondl
1513    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
1514    ! we recalculate radocondi to account for contributions from the precipitation flux
1515    ! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
1516
1517    IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
1518        ! for the solid phase (crystals + snowflakes)
1519        ! we recalculate radocondi by summing
1520        ! the ice content calculated in the mesh
1521        ! + the contribution of the non-evaporated snowfall
1522        ! from the overlying layer
1523        DO i=1,klon
1524           IF (ziflprev(i) .NE. 0.0) THEN
1525              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1526           ELSE
1527              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1528           ENDIF
1529           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
1530        ENDDO
1531    ENDIF
1532
1533    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
1534    DO i=1,klon
1535        IF (radocond(i,k) .GT. 0.) THEN
1536            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
1537        ENDIF
1538    ENDDO
1539
1540    ! ----------------------------------------------------------------
1541    ! P4> Wet scavenging
1542    ! ----------------------------------------------------------------
1543
1544    !Scavenging through nucleation in the layer
1545
1546    DO i = 1,klon
1547       
1548        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1549            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1550        ELSE
1551            beta(i,k) = 0.
1552        ENDIF
1553
1554        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1555
1556        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1557
1558            IF (temp(i,k) .GE. t_glace_min) THEN
1559                zalpha_tr = a_tr_sca(3)
1560            ELSE
1561                zalpha_tr = a_tr_sca(4)
1562            ENDIF
1563
1564            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1565            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1566
1567            ! Nucleation with a factor of -1 instead of -0.5
1568            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1569
1570        ENDIF
1571
1572    ENDDO
1573
1574    ! Scavenging through impaction in the underlying layer
1575
1576    DO kk = k-1, 1, -1
1577
1578        DO i = 1, klon
1579
1580            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1581
1582                IF (temp(i,kk) .GE. t_glace_min) THEN
1583                    zalpha_tr = a_tr_sca(1)
1584                ELSE
1585                    zalpha_tr = a_tr_sca(2)
1586                ENDIF
1587
1588              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1589              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1590
1591            ENDIF
1592
1593        ENDDO
1594
1595    ENDDO
1596
1597    ! Outputs:
1598    !-------------------------------
1599    ! Precipitation fluxes at layer interfaces
1600    ! + precipitation fractions +
1601    ! temperature and water species tendencies
1602    DO i = 1, klon
1603        psfl(i,k)=zifl(i)
1604        prfl(i,k)=zrfl(i)
1605        pfraclr(i,k)=znebprecipclr(i)
1606        pfracld(i,k)=znebprecipcld(i)
1607        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1608        d_qi(i,k) = zfice(i)*zoliq(i)
1609        d_q(i,k) = zq(i) - qt(i,k)
1610        ! c_iso: same for isotopes
1611        d_t(i,k) = zt(i) - temp(i,k)
1612    ENDDO
1613
1614
1615ENDDO
1616
1617
1618  ! Rain or snow at the surface (depending on the first layer temperature)
1619  DO i = 1, klon
1620      snow(i) = zifl(i)
1621      rain(i) = zrfl(i)
1622      ! c_iso final output
1623  ENDDO
1624
1625  IF (ncoreczq>0) THEN
1626      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1627  ENDIF
1628
1629
1630RETURN
1631
1632END SUBROUTINE lscp
1633!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1634
1635END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.