Context Navigation

lscp_mod.F90 @ 4086

Last change on this file since 4086 was 4013, checked in by lguez, 3 years ago
Sync latest trunk changes to Ocean_skin
File size: 51.5 KB

Line
1	MODULE LSCP_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8	SUBROUTINE LSCP(dtime,paprs,pplay,t,q,ptconv,ratqs, &
9	d_t, d_q, d_ql, d_qi, rneb, radliq, radicefrac, &
10	rain, snow, &
11	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
12	frac_impa, frac_nucl, beta, &
13	prfl, psfl, rhcl, zqta, fraca, &
14	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
15	iflag_ice_thermo)
16
17	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
19	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
20	!--------------------------------------------------------------------------------
21	! Date: 01/2021
22	!--------------------------------------------------------------------------------
23	! Aim: Large Scale Clouds and Precipitation (LSCP)
24	!
25	! This code is a new version of the fisrtilp.F90 routine, which is itself a
26	! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
27	! and 'ilp' (il pleut, L. Li)
28	!
29	! Compared to fisrtilp, LSCP
30	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
31	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
32	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
33	! -> option "oldbug" by JYG has been removed
34	! -> iflag_t_glace >0 always
35	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
36	! -> rectangular distribution from L. Li no longer available
37	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
38	!--------------------------------------------------------------------------------
39	! References:
40	!
41	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
42	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
43	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
44	! -------------------------------------------------------------------------------
45	! Code structure:
46	!
47	! P0> Thermalization of the precipitation coming from the overlying layer
48	! P1> Evaporation of the precipitation (falling from the k+1 level)
49	! P2> Cloud formation (at the k level)
50	! P2.A.0> Cloud properties calculation from a rectangular pdf
51	! P2.A.1> With the new PDFs, calculation of cloud properties using the inital
52	! values of T and Q
53	! P2.A.2> Coupling between condensed water and temperature
54	! P2.A.3> Calculation of final quantities associated with cloud formation
55	! P2.B> 'All or nothing' cloud
56	! P2.C> Release of Latent heat after cloud formation
57	! P3> Autoconversion to precipitation (k-level)
58	! P4> Wet scavenging
59	!------------------------------------------------------------------------------
60	! Some preliminary comments (JBM) :
61	!
62	! The cloud water that the radiation scheme sees is not the same that the cloud
63	! water used in the physics and the dynamics
64	!
65	! During the autoconversion to precipitation (P3 step), radliq (cloud water used
66	! by the radiation scheme) is calculated as an average of the water that remains
67	! in the cloud during the precipitation and not the water remaining at the end
68	! of the time step. The latter is used in the rest of the physics and advected
69	! by the dynamics.
70	!
71	! In summary:
72	!
73	! Radiation:
74	! xflwc(newmicro)+xfiwc(newmicro) =
75	! cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc)
76	!
77	! Physics/Dynamics:
78	! ql_seri(physiq)+qs_seri(physiq)=ocond(nc)
79	!
80	! Notetheless, be aware of:
81	!
82	! radliq(fisrt) .NE. ocond(nc)
83	! i.e.:
84	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
85	! (which is not trivial)
86	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87
88	USE dimphy
89	USE print_control_mod, ONLY: prt_level, lunout
90	USE cloudth_mod
91	USE ioipsl_getin_p_mod, ONLY : getin_p
92	USE phys_local_var_mod, ONLY: ql_seri,qs_seri
93	USE phys_local_var_mod, ONLY: rneblsvol
94	USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY
95
96
97	IMPLICIT NONE
98
99
100
101	!===============================================================================
102	! VARIABLE DECLARATION
103	!===============================================================================
104
105	include "YOMCST.h"
106	include "YOETHF.h"
107	include "FCTTRE.h"
108	include "fisrtilp.h"
109	include "nuage.h"
110
111
112	! INPUT VARIABLES:
113	!-----------------
114
115	REAL, INTENT(IN) :: dtime ! time step [s]
116	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
117	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
118	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
119	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg]
120	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
121	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
122	! CR: if iflag_ice_thermo=2, only convection is active
123	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
124
125
126	!Inputs associated with thermal plumes
127
128	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K]
129	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg]
130	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
131	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp)
132	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K]
133	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K]
134
135	! INPUT/OUTPUT variables
136	!------------------------
137
138	REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale
139	! cloud PDF (sigma=ratqs*qt)
140
141	! OUTPUT variables
142	!-----------------
143
144	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
145	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
146	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
147	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
148	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
149	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! condensed water used in the radiation scheme [kg/kg]
150	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
151	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
152	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! large-scale rainfall [kg/s/m2]
153	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! large-scale snowfall [kg/s/m2]
154	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
155	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
156
157
158	! cofficients of scavenging fraction (for off-line)
159
160	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl ! scavenging fraction due tu nucleation [-]
161	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl ! scavenging fraction due tu nucleation with a -1 factor [-]
162	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa ! scavening fraction due to impaction [-]
163
164	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
165
166	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
167	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
168
169
170
171
172	! PROGRAM PARAMETERS:
173	!--------------------
174
175
176	REAL, SAVE :: seuil_neb=0.001 ! cloud fraction threshold: a cloud really exists when exceeded
177	!$OMP THREADPRIVATE(seuil_neb)
178
179	INTEGER, SAVE :: ninter=5 ! number of iterations to calculate autoconversion to precipitation
180	!$OMP THREADPRIVATE(ninter)
181
182	INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evap. flag. 0: nothing, 1: "old way", 2: Max cloud fraction above to calculate the max of reevaporation
183	! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
184	!$OMP THREADPRIVATE(iflag_evap_prec)
185
186
187	REAL t_coup ! temperature threshold which determines the phase
188	PARAMETER (t_coup=234.0) ! for which the saturation vapor pressure is calculated
189
190	REAL DDT0 ! iteration parameter
191	PARAMETER (DDT0=.01)
192
193	REAL ztfondue ! parameter to calculate melting fraction of precipitation
194	PARAMETER (ztfondue=278.15)
195
196	REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associates precipitation fraction
197	!$OMP THREADPRIVATE(rain_int_min)
198
199
200
201	! LOCAL VARIABLES:
202	!----------------
203
204
205	REAL qsl, qsi
206	REAL zct, zcl
207	REAL ctot(klon,klev)
208	REAL ctot_vol(klon,klev)
209	INTEGER mpc_bl_points(klon,klev)
210	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
211	REAL zdqsdT_raw(klon)
212	REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T
213	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
214	REAL cste
215	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
216	REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
217	REAL erf
218	REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
219	REAL zrfl(klon), zrfln(klon), zqev, zqevt
220	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
221	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
222	REAL zoliqp(klon), zoliqi(klon)
223	REAL zt(klon)
224	REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
225	REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
226	REAL zmelt,zpluie,zice
227	REAL dzfice(klon)
228	REAL zsolid,qtot,dqsl,dqsi
229	REAL smallestreal
230	! Variables for Bergeron process
231	REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
232	REAL zqpreci(klon), zqprecl(klon)
233	! Variables precipitation energy conservation
234	REAL zmqc(klon)
235	! Variables for tracers
236	! temporary: alpha parameter for scavenging
237	! 4 possible scavenging processes
238	REAL a_tr_sca(4)
239	save a_tr_sca
240	!$OMP THREADPRIVATE(a_tr_sca)
241	REAL zalpha_tr
242	REAL zfrac_lessi
243	REAL zprec_cond(klon)
244	REAL beta(klon,klev) ! conversion rate of condensed water
245	REAL zmair(klon), zcpair, zcpeau
246	REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion
247	REAL zrflclr(klon), zrflcld(klon)
248	REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
249	REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
250	REAL ziflclr(klon), ziflcld(klon)
251	REAL znebprecipclr(klon), znebprecipcld(klon)
252	REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
253	REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
254	REAL velo(klon)
255	REAL qlmpc, qimpc, rnebmpc
256	REAL radliqi(klon,klev)
257
258	INTEGER i, k, n, kk
259	INTEGER n_i(klon), iter
260	INTEGER ncoreczq
261	INTEGER,SAVE :: itap=0
262	!$OMP THREADPRIVATE(itap)
263
264
265	LOGICAL lognormale(klon)
266	LOGICAL convergence(klon)
267	LOGICAL appel1er
268	SAVE appel1er
269	!$OMP THREADPRIVATE(appel1er)
270
271
272	CHARACTER (len = 20) :: modname = 'lscp'
273	CHARACTER (len = 80) :: abort_message
274
275
276	DATA appel1er /.TRUE./
277
278
279
280
281	!===============================================================================
282	! INITIALISATION
283	!===============================================================================
284
285	! Few initial checksS
286
287	IF (iflag_t_glace.EQ.0) THEN
288	abort_message = 'lscp cannot be used if iflag_t_glace=0'
289	CALL abort_physic(modname,abort_message,1)
290	ENDIF
291
292	IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
293	abort_message = 'lscp cannot be used without ice thermodynamics'
294	CALL abort_physic(modname,abort_message,1)
295	ENDIF
296
297	IF (.NOT.thermcep) THEN
298	abort_message = 'lscp cannot be used when thermcep=false'
299	CALL abort_physic(modname,abort_message,1)
300	ENDIF
301
302	IF (iflag_fisrtilp_qsat .LT. 0) THEN
303	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
304	CALL abort_physic(modname,abort_message,1)
305	ENDIF
306
307
308	! Few initialisations
309
310	itap=itap+1
311	znebprecip(:)=0.
312	ctot_vol(1:klon,1:klev)=0.0
313	rneblsvol(1:klon,1:klev)=0.0
314	smallestreal=1.e-9
315	znebprecipclr(:)=0.
316	znebprecipcld(:)=0.
317	mpc_bl_points(:,:)=0
318
319	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
320
321	IF (appel1er) THEN
322	CALL getin_p('ninter',ninter)
323	CALL getin_p('iflag_evap_prec',iflag_evap_prec)
324	CALL getin_p('seuil_neb',seuil_neb)
325	CALL getin_p('rain_int_min',rain_int_min)
326	WRITE(lunout,*) 'lscp, ninter:', ninter
327	WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
328	WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
329	WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min
330
331	! check for precipitation sub-time steps
332	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
333	WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime
334	WRITE(lunout,*) 'I would prefer a 6 min sub-timestep'
335	ENDIF
336
337	!AA Temporary initialisation
338	a_tr_sca(1) = -0.5
339	a_tr_sca(2) = -0.5
340	a_tr_sca(3) = -0.5
341	a_tr_sca(4) = -0.5
342
343	!AA Set scavenged fractions to 1
344	DO k = 1, klev
345	DO i = 1, klon
346	pfrac_nucl(i,k)=1.
347	pfrac_1nucl(i,k)=1.
348	pfrac_impa(i,k)=1.
349	beta(i,k)=0.
350	ENDDO
351	ENDDO
352
353	appel1er = .FALSE.
354
355	ENDIF !(appel1er)
356
357	! AA for 'safety' reasons
358	zalpha_tr = 0.
359	zfrac_lessi = 0.
360
361
362	! Initialisation of output variables (JYG):
363	prfl(:,:) = 0.0
364	psfl(:,:) = 0.0
365
366	d_t(:,:) = 0.0
367	d_q(:,:) = 0.0
368	d_ql(:,:) = 0.0
369	d_qi(:,:) = 0.0
370	rneb(:,:) = 0.0
371	radliq(:,:) = 0.0
372	radicefrac(:,:) = 0.0
373	frac_nucl(:,:) = 1.
374	frac_impa(:,:) = 1.
375
376	rain(:) = 0.0
377	snow(:) = 0.0
378	zoliq(:)=0.
379	zrfl(:) = 0.0
380	zifl(:) = 0.0
381	zneb(:) = seuil_neb
382	zrflclr(:) = 0.0
383	ziflclr(:) = 0.0
384	zrflcld(:) = 0.0
385	ziflcld(:) = 0.0
386	tot_zneb(:) = 0.0
387	tot_znebn(:) = 0.0
388	d_tot_zneb(:) = 0.0
389
390
391	!===============================================================================
392	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
393	!===============================================================================
394
395	ncoreczq=0
396
397	DO k = klev, 1, -1
398
399
400	! Initialisation temperature and specific humidity
401	DO i = 1, klon
402	zt(i)=t(i,k)
403	zq(i)=q(i,k)
404	ENDDO
405
406	! --------------------------------------------------------------------
407	! P0> Thermalization of precipitation falling from the overlying layer
408	! --------------------------------------------------------------------
409	! Computes air temperature variation due to latent heat transported by
410	! precipitation. Precipitation is then thermalized with the air in the
411	! layer.
412	! The precipitation should remain thermalized throughout the different
413	! thermodynamical transformations. The corresponding water mass should
414	! be added when calculating the layer's enthalpy change with
415	! temperature
416	! ---------------------------------------------------------------------
417
418	IF (k.LE.klevm1) THEN
419
420	DO i = 1, klon
421
422	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
423	! no condensed water so cp=cp(vapor+dry air)
424	! RVTMP2=rcpv/rcpd-1
425	zcpair=RCPD(1.0+RVTMP2zq(i))
426	zcpeau=RCPD*RVTMP2
427
428	! zmqc: precipitation mass that has to be thermalized with
429	! layer's air so that precipitation at the ground has the
430	! same temperature as the lowermost layer
431	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
432	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
433	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zmqc(i)zcpeau + zcpair*zt(i) ) &
434	/ (zcpair + zmqc(i)*zcpeau)
435
436	ENDDO
437
438	ELSE
439
440	DO i = 1, klon
441	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
442	zmqc(i) = 0.
443	ENDDO
444
445	ENDIF
446
447	! --------------------------------------------------------------------
448	! P1> Precipitation evaporation/sublimation
449	! --------------------------------------------------------------------
450	! A part of the precipitation coming from above is evaporated/sublimated.
451	!
452	! --------------------------------------------------------------------
453
454
455	IF (iflag_evap_prec.GE.1) THEN
456
457
458	DO i = 1, klon
459
460	! if precipitation
461	IF (zrfl(i)+zifl(i).GT.0.) THEN
462
463	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
464
465	! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
466	IF (iflag_evap_prec.EQ.4) THEN
467	zrfl(i) = zrflclr(i)
468	zifl(i) = ziflclr(i)
469	ENDIF
470
471	IF (iflag_evap_prec.EQ.1) THEN
472	znebprecip(i)=zneb(i)
473	ELSE
474	znebprecip(i)=MAX(zneb(i),znebprecip(i))
475	ENDIF
476
477
478	IF (iflag_evap_prec.EQ.4) THEN
479	! Max evaporation not to saturate the whole mesh
480	zqev0 = MAX(0.0, zqs(i)-zq(i))
481	ELSE
482	! Max evap not to saturate the fraction below the cloud
483	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
484	ENDIF
485
486	! A. JAM
487	! We consider separately qsatice and qstal
488	! qsat wrt liquid phase according to thermcep
489	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl)
490
491
492	! Evaporation of liquid precipitation coming from above
493	! dP/dz=beta(1-q/qsat)sqrt(P)
494	! formula from Sundquist 1988, Klemp & Wilhemson 1978
495	! LTP: evaporation only in the clear sky part (iflag_evap_prec=4)
496
497	IF (iflag_evap_prec.EQ.3) THEN
498	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl) &
499	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
500	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
501	ELSE IF (iflag_evap_prec.EQ.4) THEN
502	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl) &
503	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
504	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
505	ELSE
506	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
507	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
508	ENDIF
509
510
511	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
512	RGdtime/(paprs(i,k)-paprs(i,k+1))
513
514	! qsat wrt ice according to thermcep
515	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)
516
517	! sublimation of the solid precipitation coming from above
518	IF (iflag_evap_prec.EQ.3) THEN
519	zqevti = znebprecip(i)coef_eva(1.0-zq(i)/qsi) &
520	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
521	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
522	ELSE IF (iflag_evap_prec.EQ.4) THEN
523	zqevti = znebprecipclr(i)coef_eva(1.0-zq(i)/qsi) &
524	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
525	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
526	ELSE
527	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
528	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
529	ENDIF
530
531	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
532	RGdtime/(paprs(i,k)-paprs(i,k+1))
533
534	! A. JAM
535	! Evaporation limit: we ensure that the layer's fraction below
536	! the cloud or the whole mesh (depending on iflag_evap_prec)
537	! does not reach saturation. In this case, we
538	! redistribute zqev0 conserving the ratio liquid/ice
539
540	IF (zqevt+zqevti.GT.zqev0) THEN
541	zqev=zqev0*zqevt/(zqevt+zqevti)
542	zqevi=zqev0*zqevti/(zqevt+zqevti)
543	ELSE
544	zqev=zqevt
545	zqevi=zqevti
546	ENDIF
547
548
549	! EV: agrees with JL's comments below so correct and comment
550	! JLD: I do not understand the lines below. We distribute the liquid
551	! and solid parts of the precipitation even if the layer does not get
552	! saturated. To my opinion, we should all replace with:
553	! zqev=zqevt
554	! zqevi=zqevti
555	! IF (zqevt+zqevti.GT.0.) THEN
556	! zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
557	! zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
558	! ELSE
559	! zqev=0.
560	! zqevi=0.
561	! ENDIF
562	! ENDIF
563
564	! New solid and liquid precipitation fluxes
565	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
566	/RG/dtime)
567	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
568	/RG/dtime)
569
570	! vapor, temperature, precip fluxes update
571	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
572	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
573	! precip thermalization
574	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
575	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
576	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
577	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
578	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
579	+ (zifln(i)-zifl(i)) &
580	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
581	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
582
583	! New values of liquid and solid precipitation
584	zrfl(i) = zrfln(i)
585	zifl(i) = zifln(i)
586
587	IF (iflag_evap_prec.EQ.4) THEN
588	zrflclr(i) = zrfl(i)
589	ziflclr(i) = zifl(i)
590	IF(zrflclr(i) + ziflclr(i).LE.0) THEN
591	znebprecipclr(i) = 0.
592	ENDIF
593	zrfl(i) = zrflclr(i) + zrflcld(i)
594	zifl(i) = ziflclr(i) + ziflcld(i)
595	ENDIF
596
597
598	! CR Be careful: ice melting has been moved
599	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
600	! precip fraction that is melted
601	zmelt = MIN(MAX(zmelt,0.),1.)
602	! Ice melting
603	IF (iflag_evap_prec.EQ.4) THEN
604	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
605	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
606	zrfl(i)=zrflclr(i)+zrflcld(i)
607	ELSE
608	zrfl(i)=zrfl(i)+zmelt*zifl(i)
609	ENDIF
610
611	! Latent heat of melting with precipitation thermalization
612	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
613	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
614
615	IF (iflag_evap_prec.EQ.4) THEN
616	ziflclr(i)=ziflclr(i)*(1.-zmelt)
617	ziflcld(i)=ziflcld(i)*(1.-zmelt)
618	zifl(i)=ziflclr(i)+ziflcld(i)
619	ELSE
620	zifl(i)=zifl(i)*(1.-zmelt)
621	ENDIF
622
623	ELSE
624	! if no precip, we reinitialize the cloud fraction used for the precip to 0
625	znebprecip(i)=0.
626
627	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
628
629	ENDDO
630
631	ENDIF ! (iflag_evap_prec>=1)
632
633	! --------------------------------------------------------------------
634	! End precip evaporation
635	! --------------------------------------------------------------------
636
637
638
639
640	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
641	!-------------------------------------------------------
642
643	DO i = 1, klon
644	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
645	qtot=zq(i)+zmqc(i)
646	CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
647	zdqsdT_raw(i) = zdqs(i)* &
648	& RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTTzdelta)
649
650	IF (zq(i) .LT. 1.e-15) THEN
651	ncoreczq=ncoreczq+1
652	zq(i)=1.e-15
653	ENDIF
654
655	ENDDO
656
657
658	! --------------------------------------------------------------------
659	! P2> Cloud formation
660	!---------------------------------------------------------------------
661	!
662	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
663	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
664	! is prescribed and depends on large scale variables and boundary layer
665	! properties)
666	! The decrease in condensed part due tu latent heating is taken into
667	! account
668	! -------------------------------------------------------------------
669
670
671	! P2.1> With the PDFs (log-normal, bigaussian)
672	! cloud propertues calculation with the initial values of t and q
673	! ----------------------------------------------------------------
674
675	! initialise gammasat and qincloud_mpc
676	gammasat(:)=1.
677	qincloud_mpc(:)=0.
678
679
680
681	IF (iflag_cld_th.GE.5) THEN
682
683	IF (iflag_mpc_bl .LT. 1) THEN
684
685	IF (iflag_cloudth_vert.LE.2) THEN
686
687	CALL cloudth(klon,klev,k,ztv, &
688	zq,zqta,fraca, &
689	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
690	ratqs,zqs,t)
691
692	ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
693
694	CALL cloudth_v3(klon,klev,k,ztv, &
695	zq,zqta,fraca, &
696	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
697	ratqs,zqs,t)
698
699
700	!Jean Jouhaud's version, Decembre 2018
701	!-------------------------------------
702
703	ELSEIF (iflag_cloudth_vert.EQ.6) THEN
704
705	CALL cloudth_v6(klon,klev,k,ztv, &
706	zq,zqta,fraca, &
707	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
708	ratqs,zqs,t)
709
710	ENDIF
711
712	ELSE
713	! cloudth_v3 + cold microphysical considerations
714	! + stationary mixed-phase cloud model
715
716	CALL cloudth_mpc(klon,klev,k,mpc_bl_points, &
717	t,ztv,zq,zqta,fraca, zpspsk, &
718	paprs,pplay,ztla,zthl,ratqs,zqs,psfl, &
719	qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
720
721	ENDIF ! iflag_mpc_bl
722
723	DO i=1,klon
724	rneb(i,k)=ctot(i,k)
725	rneblsvol(i,k)=ctot_vol(i,k)
726	zqn(i)=qcloud(i)
727	ENDDO
728
729	ENDIF
730
731	IF (iflag_cld_th .LE. 4) THEN
732
733	! lognormal
734	lognormale = .TRUE.
735
736	ELSEIF (iflag_cld_th .GE. 6) THEN
737
738	! lognormal distribution when no thermals
739	lognormale = fraca(:,k) < 1e-10
740
741	ELSE
742	! When iflag_cld_th=5, we always assume
743	! bi-gaussian distribution
744	lognormale = .FALSE.
745
746	ENDIF
747
748
749	DT(:) = 0.
750	n_i(:)=0
751	Tbef(:)=zt(:)
752	qlbef(:)=0.
753
754	! Treatment of non-boundary layer clouds (lognormale)
755	! condensation with qsat(T) variation (adaptation)
756	! Iterative Loop to converge towards qsat
757
758
759	DO iter=1,iflag_fisrtilp_qsat+1
760
761	DO i=1,klon
762
763	! convergence = .true. since convergence is not satisfied
764	convergence(i)=ABS(DT(i)).GT.DDT0
765
766	IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
767
768	! if not convergence:
769
770	! P2.2.1> cloud fraction and condensed water mass calculation
771	! Calculated variables:
772	! rneb : cloud fraction
773	! zqn : total water within the cloud
774	! zcond: mean condensed water within the mesh
775	! rhcl: clear-sky relative humidity
776	!---------------------------------------------------------------
777
778	! new temperature:
779	Tbef(i)=Tbef(i)+DT(i)
780
781
782	! Rneb, qzn and zcond for lognormal PDFs
783	qtot=zq(i)+zmqc(i)
784	CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
785	CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i))
786
787	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
788	zdqs(i) = gammasat(i)zdqs(i)+zqs(i)dgammasatdt(i)
789
790	zpdf_sig(i)=ratqs(i,k)*zq(i)
791	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
792	zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
793	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
794	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
795	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
796	zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
797	zpdf_e1(i)=1.-erf(zpdf_e1(i))
798	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
799	zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
800	zpdf_e2(i)=1.-erf(zpdf_e2(i))
801
802	IF (zpdf_e1(i).LT.1.e-10) THEN
803	rneb(i,k)=0.
804	zqn(i)=gammasat(i)*zqs(i)
805	ELSE
806	rneb(i,k)=0.5*zpdf_e1(i)
807	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
808	ENDIF
809
810	! If vertical heterogeneity, change fraction by volume as well
811	IF (iflag_cloudth_vert.GE.3) THEN
812	ctot_vol(i,k)=rneb(i,k)
813	rneblsvol(i,k)=ctot_vol(i,k)
814	ENDIF
815
816
817	! P2.2.2> Approximative calculation of temperature variation DT
818	! due to condensation.
819	! Calculated variables:
820	! dT : temperature change due to condensation
821	! ---------------------------------------------------------------
822
823	! EV: calculation of icefrac in one sole function
824	CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
825
826
827	IF (zfice(i).LT.1) THEN
828	cste=RLVTT
829	ELSE
830	cste=RLSTT
831	ENDIF
832
833	qlbef(i)=max(0.,zqn(i)-zqs(i))
834	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
835	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
836	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
837	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
838	& qlbef(i)dzfice(i)
839	DT(i)=num/denom
840	n_i(i)=n_i(i)+1
841
842	ENDIF ! end convergence
843
844	ENDDO ! end loop on i
845
846	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
847
848
849
850	! P2.3> Final quantities calculation
851	! Calculated variables:
852	! rneb : cloud fraction
853	! zcond: mean condensed water in the mesh
854	! zqn : mean water vapor in the mesh
855	! zt : temperature
856	! rhcl : clear-sky relative humidity
857	! ----------------------------------------------------------------
858
859
860
861	! Water vapor update, Phase determination and subsequent latent heat exchange
862
863	! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
864	CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
865
866	DO i=1, klon
867
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
897	! water vapor update
898	zq(i) = zq(i) - zcond(i)
899
900	ENDIF
901
902	! temperature update due to phase change
903	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
904	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
905	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
906	ENDDO
907
908
909	! If vertical heterogeneity, change volume fraction
910	IF (iflag_cloudth_vert .GE. 3) THEN
911	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
912	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
913	ENDIF
914
915	! ----------------------------------------------------------------
916	! P3> Precipitation formation and calculation of condensed
917	! water seen by the radiation scheme
918	! ----------------------------------------------------------------
919
920	! LTP:
921	IF (iflag_evap_prec==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) &
947	- d_tot_zneb(i) - zneb(i)))
948	IF (znebprecipclr(i) .GT. 0) THEN
949	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
950	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
951	ELSE
952	d_zrfl_clr_cld(i) = 0.
953	d_zifl_clr_cld(i) = 0.
954	ENDIF
955
956	!Update variables
957	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
958	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
959	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
960	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
961	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
962	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
963
964	ENDDO
965	ENDIF
966
967	! Initialisation of zoliq and radliq variables
968
969	DO i = 1, klon
970	IF (rneb(i,k).GT.0.0) THEN
971	zoliq(i) = zcond(i)
972	zrho(i) = pplay(i,k) / zt(i) / RD
973	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
974
975	zneb(i) = MAX(rneb(i,k), seuil_neb)
976	radliq(i,k) = zoliq(i)/REAL(ninter+1)
977	radicefrac(i,k) = zfice(i)
978	radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
979	ENDIF
980	ENDDO
981
982	! -------------------------------------------------------------------------------
983	! Iterations to calculate a "mean" cloud water content during the precipitation
984	! Comment: it is not the remaining water after precipitation but a mean
985	! remaining water in the cloud during the time step that is seen by the radiation
986	! -------------------------------------------------------------------------------
987
988
989	DO n = 1, ninter
990
991	DO i=1,klon
992	IF (rneb(i,k).GT.0.0) THEN
993	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
994	ENDIF
995	ENDDO
996
997	CALL FALLICE_VELOCITY(klon,zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))
998
999	DO i = 1, klon
1000
1001	IF (rneb(i,k).GT.0.0) THEN
1002
1003	! Initialization of zpluie, zice and ztot:
1004	zpluie=0.
1005	zice=0.
1006	ztot=0.
1007
1008	IF (zneb(i).EQ.seuil_neb) THEN
1009	ztot = 0.0
1010	zice = 0.0
1011	zpluie= 0.0
1012	ELSE
1013	! water quantity to remove: zchau (Sundqvist, 1978)
1014	! same thing for the ice: zfroi (Zender & Kiehl, 1997)
1015	IF (ptconv(i,k)) THEN ! if convective point
1016	zcl=cld_lc_con
1017	zct=1./cld_tau_con
1018	zfroi = dtime/REAL(ninter)/zdz(i)zoliq(i)velo(i)*zfice(i)
1019	ELSE
1020	zcl=cld_lc_lsc
1021	zct=1./cld_tau_lsc
1022	zfroi = dtime/REAL(ninter)/zdz(i)zoliq(i) & ! dqice/dt=1/rhod(rhowiceqice)/dz
1023	velo(i) zfice(i)
1024	ENDIF
1025
1026	! if vertical heterogeneity is taken into account, we use
1027	! the "true" volume fraction instead of a modified
1028	! surface fraction (which is larger and artificially
1029	! reduces the in-cloud water).
1030	IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
1031	zchau = zct dtime/REAL(ninter) zoliq(i) &
1032	(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)2)) (1.-zfice(i))
1033	ELSE
1034	zchau = zct dtime/REAL(ninter) zoliq(i) &
1035	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)2)) (1.-zfice(i)) ! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
1036	ENDIF
1037
1038	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1039	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1040	ztot = MAX(zpluie + zice,0.0)
1041
1042	ENDIF
1043
1044	ztot = MIN(ztot,zoliq(i))
1045	zice = MIN(zice,ztot)
1046	zpluie = MIN(zpluie,ztot)
1047
1048	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
1049	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
1050	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
1051
1052	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1053	radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
1054
1055	ENDIF
1056
1057	ENDDO ! i = 1,klon
1058
1059	ENDDO ! n = 1,niter
1060
1061
1062	! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
1063	DO i=1,klon
1064	IF (radliq(i,k) .GT. 0) THEN
1065	radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
1066	ENDIF
1067	ENDDO
1068
1069
1070	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1071	! If T<0C, liquid precip are converted into ice, which leads to an increase in
1072	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1073	! taken into account through a linearization of the equations and by approximating
1074	! the liquid precipitation process with a "threshold" process. We assume that
1075	! condensates are not modified during this operation. Liquid precipitation is
1076	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1077	! Water vapor increases as well
1078	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1079
1080
1081
1082	DO i = 1, klon
1083
1084	IF (rneb(i,k) .GT. 0.0) THEN
1085
1086	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1087	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1088	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
1089	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
1090	! Computation of DT if all the liquid precip freezes
1091	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
1092	! T should not exceed the freezing point
1093	! that is Delta > RTT-zt(i)
1094	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1095	zt(i) = zt(i) + DeltaT
1096	! water vaporization due to temp. increase
1097	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
1098	! we add this vaporized water to the total vapor and we remove it from the precip
1099	zq(i) = zq(i) + Deltaq
1100	! The three "max" lines herebelow protect from rounding errors
1101	zcond(i) = max( zcond(i)- Deltaq, 0. )
1102	! liquid precipitation freezes oe evaporates
1103	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
1104	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1105	! iced water budget
1106	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1107
1108	d_ql(i,k) = (1-zfice(i))*zoliq(i)
1109	d_qi(i,k) = zfice(i)*zoliq(i)
1110
1111	IF (iflag_evap_prec.EQ.4) THEN
1112	zrflcld(i) = zrflcld(i)+zqprecl(i) &
1113	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1114	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1115	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1116	znebprecipcld(i) = rneb(i,k)
1117	zrfl(i) = zrflcld(i) + zrflclr(i)
1118	zifl(i) = ziflcld(i) + ziflclr(i)
1119	ELSE
1120	zrfl(i) = zrfl(i)+ zqprecl(i) &
1121	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1122	zifl(i) = zifl(i)+ zqpreci(i) &
1123	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
1124	ENDIF
1125
1126	ENDIF ! rneb>0
1127
1128	ENDDO
1129
1130
1131
1132
1133	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1134	! if iflag_evap_pre=4
1135	IF (iflag_evap_prec.EQ.4) THEN
1136
1137	DO i=1, klon
1138
1139
1140	IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1141	znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1142	(MAX(znebprecipclr(i),seuil_neb)rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)rain_int_min)))
1143	ELSE
1144	znebprecipclr(i)=0.
1145	ENDIF
1146
1147
1148	IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1149	znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1150	(MAX(znebprecipcld(i),seuil_neb)rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)rain_int_min)))
1151	ELSE
1152	znebprecipcld(i)=0.
1153	ENDIF
1154
1155
1156	ENDDO
1157
1158	ENDIF
1159
1160	! --------------------------------
1161	! End of precipitation formation
1162	! --------------------------------
1163
1164
1165	! Outputs:
1166	! Precipitation fluxes at layer interfaces
1167	! and temperature and vapor tendencies
1168	DO i = 1, klon
1169	psfl(i,k)=zifl(i)
1170	prfl(i,k)=zrfl(i)
1171	d_q(i,k) = zq(i) - q(i,k)
1172	d_t(i,k) = zt(i) - t(i,k)
1173	ENDDO
1174
1175
1176
1177	! ----------------------------------------------------------------
1178	! P4> Wet scavenging
1179	! ----------------------------------------------------------------
1180
1181
1182	!Scavenging through nucleation in the layer
1183
1184	DO i = 1,klon
1185
1186	IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1187	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1188	ELSE
1189	beta(i,k) = 0.
1190	ENDIF
1191
1192	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1193	* (paprs(i,k)-paprs(i,k+1))/RG
1194
1195
1196	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1197
1198
1199	IF (t(i,k) .GE. t_glace_min) THEN
1200	zalpha_tr = a_tr_sca(3)
1201	ELSE
1202	zalpha_tr = a_tr_sca(4)
1203	ENDIF
1204
1205
1206	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1207	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
1208	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1209
1210	! Nucleation with a factor of -1 instead of -0.5
1211	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1212	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
1213
1214	ENDIF
1215
1216	ENDDO
1217
1218
1219	! Scavenging through impaction in the underlying layer
1220
1221	DO kk = k-1, 1, -1
1222
1223	DO i = 1, klon
1224
1225	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1226
1227	IF (t(i,kk) .GE. t_glace_min) THEN
1228	zalpha_tr = a_tr_sca(1)
1229	ELSE
1230	zalpha_tr = a_tr_sca(2)
1231	ENDIF
1232
1233	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1234	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
1235	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1236
1237	ENDIF
1238
1239	ENDDO
1240
1241	ENDDO
1242
1243
1244	END DO
1245
1246	!======================================================================
1247	! END OF VERTICAL LOOP
1248	!======================================================================
1249
1250
1251	! Rain or snow at the surface (depending on the first layer temperature)
1252	DO i = 1, klon
1253	snow(i) = zifl(i)
1254	rain(i) = zrfl(i)
1255	ENDDO
1256
1257
1258
1259	IF (ncoreczq>0) THEN
1260	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1261	ENDIF
1262
1263
1264
1265
1266	END SUBROUTINE LSCP
1267	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1268
1269	END MODULE LSCP_MOD

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

Download in other formats:

Original Format