Context Navigation

lscp_mod.F90 @ 4076

Last change on this file since 4076 was 4072, checked in by evignon, 4 years ago

Mise à jour de la routine de nuages lscp
les principaux changements consistent en:

des corrections de bug (déclaration et division

par zero) pour que la routine compile avec les modifications d'OB

des restructurations pour que les fonctions de lscp_tools_mod

travaillent sur des tableaux et non des scalaires

une séparation des coefficients d'évaporation de la neige et de la pluie

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

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

Download in other formats:

Original Format