Context Navigation

lmdz_lscp.F90 @ 5136

Last change on this file since 5136 was 5117, checked in by abarral, 4 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 63.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format