Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_lscp.F90 @ 4743

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

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

Download in other formats:

Original Format