Context Navigation

lmdz_lscp.F90 @ 5441

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

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

Download in other formats:

Original Format