Context Navigation

lscp_mod.F90 @ 4059

Last change on this file since 4059 was 4059, checked in by oboucher, 2 years ago
Audran Borella's parametrisation for ice supersaturation activated with flag_ice_sursat (FALSE by default)
File size: 53.9 KB

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

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

Download in other formats:

Original Format