Context Navigation

lscp_mod.F90 @ 4099

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

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

Download in other formats:

Original Format