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

Last change on this file since 5396 was 5396, checked in by evignon, 5 weeks ago

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