Context Navigation

lscp_mod.F90 @ 4390

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

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

Download in other formats:

Original Format