Context Navigation

lmdz_cloudth.F90 @ 4664

Last change on this file since 4664 was 4664, checked in by fhourdin, 9 months ago

standardisatio des noms pour lscp et fisrtilp

fisrtilp passe dans le module lmdz_lscp_old.F90
Prepartation de la replaysation de fisrtilp (deja fait pour lscp)

File size: 88.7 KB

Line
1	MODULE lmdz_cloudth
2
3
4	IMPLICIT NONE
5
6	CONTAINS
7
8	SUBROUTINE cloudth(ngrid,klev,ind2, &
9	& ztv,po,zqta,fraca, &
10	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
11	& ratqs,zqs,t, &
12	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
13
14
15	use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
16
17	IMPLICIT NONE
18
19
20	!===========================================================================
21	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
22	! Date : 25 Mai 2010
23	! Objet : calcule les valeurs de qc et rneb dans les thermiques
24	!===========================================================================
25
26
27	INCLUDE "YOMCST.h"
28	INCLUDE "YOETHF.h"
29	INCLUDE "FCTTRE.h"
30
31	INTEGER itap,ind1,ind2
32	INTEGER ngrid,klev,klon,l,ig
33	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
34
35	REAL ztv(ngrid,klev)
36	REAL po(ngrid)
37	REAL zqenv(ngrid)
38	REAL zqta(ngrid,klev)
39
40	REAL fraca(ngrid,klev+1)
41	REAL zpspsk(ngrid,klev)
42	REAL paprs(ngrid,klev+1)
43	REAL pplay(ngrid,klev)
44	REAL ztla(ngrid,klev)
45	REAL zthl(ngrid,klev)
46
47	REAL zqsatth(ngrid,klev)
48	REAL zqsatenv(ngrid,klev)
49
50
51	REAL sigma1(ngrid,klev)
52	REAL sigma2(ngrid,klev)
53	REAL qlth(ngrid,klev)
54	REAL qlenv(ngrid,klev)
55	REAL qltot(ngrid,klev)
56	REAL cth(ngrid,klev)
57	REAL cenv(ngrid,klev)
58	REAL ctot(ngrid,klev)
59	REAL rneb(ngrid,klev)
60	REAL t(ngrid,klev)
61	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
62	REAL rdd,cppd,Lv
63	REAL alth,alenv,ath,aenv
64	REAL sth,senv,sigma1s,sigma2s,xth,xenv
65	REAL Tbef,zdelta,qsatbef,zcor
66	REAL qlbef
67	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
68
69	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
70	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
71	REAL zqs(ngrid), qcloud(ngrid)
72	REAL erf
73
74
75
76
77	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78	! Gestion de deux versions de cloudth
79	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80
81	IF (iflag_cloudth_vert.GE.1) THEN
82	CALL cloudth_vert(ngrid,klev,ind2, &
83	& ztv,po,zqta,fraca, &
84	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
85	& ratqs,zqs,t)
86	RETURN
87	ENDIF
88	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
89
90
91	!-------------------------------------------------------------------------------
92	! Initialisation des variables r?elles
93	!-------------------------------------------------------------------------------
94	sigma1(:,:)=0.
95	sigma2(:,:)=0.
96	qlth(:,:)=0.
97	qlenv(:,:)=0.
98	qltot(:,:)=0.
99	rneb(:,:)=0.
100	qcloud(:)=0.
101	cth(:,:)=0.
102	cenv(:,:)=0.
103	ctot(:,:)=0.
104	qsatmmussig1=0.
105	qsatmmussig2=0.
106	rdd=287.04
107	cppd=1005.7
108	pi=3.14159
109	Lv=2.5e6
110	sqrt2pi=sqrt(2.*pi)
111
112
113
114	!-------------------------------------------------------------------------------
115	! Calcul de la fraction du thermique et des ?cart-types des distributions
116	!-------------------------------------------------------------------------------
117	do ind1=1,ngrid
118
119	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
120
121	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
122
123
124	! zqenv(ind1)=po(ind1)
125	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
126	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
127	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
128	qsatbef=MIN(0.5,qsatbef)
129	zcor=1./(1.-retv*qsatbef)
130	qsatbef=qsatbef*zcor
131	zqsatenv(ind1,ind2)=qsatbef
132
133
134
135
136	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
137	aenv=1./(1.+(alenv*Lv/cppd))
138	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
139
140
141
142
143	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
144	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
145	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
146	qsatbef=MIN(0.5,qsatbef)
147	zcor=1./(1.-retv*qsatbef)
148	qsatbef=qsatbef*zcor
149	zqsatth(ind1,ind2)=qsatbef
150
151	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
152	ath=1./(1.+(alth*Lv/cppd))
153	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
154
155
156
157	!------------------------------------------------------------------------------
158	! Calcul des ?cart-types pour s
159	!------------------------------------------------------------------------------
160
161	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
162	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
163	! if (paprs(ind1,ind2).gt.90000) then
164	! ratqs(ind1,ind2)=0.002
165	! else
166	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
167	! endif
168	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
169	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
170	! sigma1s=ratqs(ind1,ind2)*po(ind1)
171	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
172
173	!------------------------------------------------------------------------------
174	! Calcul de l'eau condens?e et de la couverture nuageuse
175	!------------------------------------------------------------------------------
176	sqrt2pi=sqrt(2.*pi)
177	xth=sth/(sqrt(2.)*sigma2s)
178	xenv=senv/(sqrt(2.)*sigma1s)
179	cth(ind1,ind2)=0.5(1.+1.erf(xth))
180	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
181	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
182
183	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
184	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
185	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
186
187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188	if (ctot(ind1,ind2).lt.1.e-10) then
189	ctot(ind1,ind2)=0.
190	qcloud(ind1)=zqsatenv(ind1,ind2)
191
192	else
193
194	ctot(ind1,ind2)=ctot(ind1,ind2)
195	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
196
197	endif
198
199
200	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
201
202
203	else ! gaussienne environnement seule
204
205	zqenv(ind1)=po(ind1)
206	Tbef=t(ind1,ind2)
207	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
208	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
209	qsatbef=MIN(0.5,qsatbef)
210	zcor=1./(1.-retv*qsatbef)
211	qsatbef=qsatbef*zcor
212	zqsatenv(ind1,ind2)=qsatbef
213
214
215	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
216	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
217	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
218	aenv=1./(1.+(alenv*Lv/cppd))
219	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
220
221
222	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
223
224	sqrt2pi=sqrt(2.*pi)
225	xenv=senv/(sqrt(2.)*sigma1s)
226	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
227	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
228
229	if (ctot(ind1,ind2).lt.1.e-3) then
230	ctot(ind1,ind2)=0.
231	qcloud(ind1)=zqsatenv(ind1,ind2)
232
233	else
234
235	ctot(ind1,ind2)=ctot(ind1,ind2)
236	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
237
238	endif
239
240
241
242
243
244
245	endif
246	enddo
247
248	return
249	! end
250	END SUBROUTINE cloudth
251
252
253
254	!===========================================================================
255	SUBROUTINE cloudth_vert(ngrid,klev,ind2, &
256	& ztv,po,zqta,fraca, &
257	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
258	& ratqs,zqs,t)
259
260	!===========================================================================
261	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
262	! Date : 25 Mai 2010
263	! Objet : calcule les valeurs de qc et rneb dans les thermiques
264	!===========================================================================
265
266
267	use lmdz_cloudth_ini, only: iflag_cloudth_vert, vert_alpha
268
269	IMPLICIT NONE
270
271	INCLUDE "YOMCST.h"
272	INCLUDE "YOETHF.h"
273	INCLUDE "FCTTRE.h"
274
275	INTEGER itap,ind1,ind2
276	INTEGER ngrid,klev,klon,l,ig
277
278	REAL ztv(ngrid,klev)
279	REAL po(ngrid)
280	REAL zqenv(ngrid)
281	REAL zqta(ngrid,klev)
282
283	REAL fraca(ngrid,klev+1)
284	REAL zpspsk(ngrid,klev)
285	REAL paprs(ngrid,klev+1)
286	REAL pplay(ngrid,klev)
287	REAL ztla(ngrid,klev)
288	REAL zthl(ngrid,klev)
289
290	REAL zqsatth(ngrid,klev)
291	REAL zqsatenv(ngrid,klev)
292
293
294	REAL sigma1(ngrid,klev)
295	REAL sigma2(ngrid,klev)
296	REAL qlth(ngrid,klev)
297	REAL qlenv(ngrid,klev)
298	REAL qltot(ngrid,klev)
299	REAL cth(ngrid,klev)
300	REAL cenv(ngrid,klev)
301	REAL ctot(ngrid,klev)
302	REAL rneb(ngrid,klev)
303	REAL t(ngrid,klev)
304	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
305	REAL rdd,cppd,Lv,sqrt2,sqrtpi
306	REAL alth,alenv,ath,aenv
307	REAL sth,senv,sigma1s,sigma2s,xth,xenv
308	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
309	REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
310	REAL Tbef,zdelta,qsatbef,zcor
311	REAL qlbef
312	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
313	! Change the width of the PDF used for vertical subgrid scale heterogeneity
314	! (J Jouhaud, JL Dufresne, JB Madeleine)
315
316	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
317	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
318	REAL zqs(ngrid), qcloud(ngrid)
319	REAL erf
320
321	!------------------------------------------------------------------------------
322	! Initialisation des variables r?elles
323	!------------------------------------------------------------------------------
324	sigma1(:,:)=0.
325	sigma2(:,:)=0.
326	qlth(:,:)=0.
327	qlenv(:,:)=0.
328	qltot(:,:)=0.
329	rneb(:,:)=0.
330	qcloud(:)=0.
331	cth(:,:)=0.
332	cenv(:,:)=0.
333	ctot(:,:)=0.
334	qsatmmussig1=0.
335	qsatmmussig2=0.
336	rdd=287.04
337	cppd=1005.7
338	pi=3.14159
339	Lv=2.5e6
340	sqrt2pi=sqrt(2.*pi)
341	sqrt2=sqrt(2.)
342	sqrtpi=sqrt(pi)
343
344	!-------------------------------------------------------------------------------
345	! Calcul de la fraction du thermique et des ?cart-types des distributions
346	!-------------------------------------------------------------------------------
347	do ind1=1,ngrid
348
349	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
350
351	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
352
353
354	! zqenv(ind1)=po(ind1)
355	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
356	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
357	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
358	qsatbef=MIN(0.5,qsatbef)
359	zcor=1./(1.-retv*qsatbef)
360	qsatbef=qsatbef*zcor
361	zqsatenv(ind1,ind2)=qsatbef
362
363
364
365
366	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
367	aenv=1./(1.+(alenv*Lv/cppd))
368	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
369
370
371
372
373	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
374	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
375	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
376	qsatbef=MIN(0.5,qsatbef)
377	zcor=1./(1.-retv*qsatbef)
378	qsatbef=qsatbef*zcor
379	zqsatth(ind1,ind2)=qsatbef
380
381	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
382	ath=1./(1.+(alth*Lv/cppd))
383	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
384
385
386
387	!------------------------------------------------------------------------------
388	! Calcul des ?cart-types pour s
389	!------------------------------------------------------------------------------
390
391	sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
392	sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
393	! if (paprs(ind1,ind2).gt.90000) then
394	! ratqs(ind1,ind2)=0.002
395	! else
396	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
397	! endif
398	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
399	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
400	! sigma1s=ratqs(ind1,ind2)*po(ind1)
401	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
402
403	!------------------------------------------------------------------------------
404	! Calcul de l'eau condens?e et de la couverture nuageuse
405	!------------------------------------------------------------------------------
406	sqrt2pi=sqrt(2.*pi)
407	xth=sth/(sqrt(2.)*sigma2s)
408	xenv=senv/(sqrt(2.)*sigma1s)
409	cth(ind1,ind2)=0.5(1.+1.erf(xth))
410	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
411	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
412
413	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
414	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
415	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
416
417	IF (iflag_cloudth_vert == 1) THEN
418	!-------------------------------------------------------------------------------
419	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
420	!-------------------------------------------------------------------------------
421	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
422	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
423	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
424	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
425	! deltasenv=aenv0.01po(ind1)
426	! deltasth=ath0.01zqta(ind1,ind2)
427	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
428	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
429	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
430	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
431	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
432	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
433
434	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
435	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
436	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
437
438	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
439	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
440	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
441	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
442
443	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
444	! qlenv(ind1,ind2)=IntJ
445	! print*, qlenv(ind1,ind2),'VERIF EAU'
446
447
448	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
449	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
450	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
451	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
452	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
453	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
454	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
455	! qlth(ind1,ind2)=IntJ
456	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
457	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
458
459	ELSE IF (iflag_cloudth_vert == 2) THEN
460
461	!-------------------------------------------------------------------------------
462	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
463	!-------------------------------------------------------------------------------
464	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
465	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
466	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
467	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
468	deltasenv=aenvvert_alphasigma1s
469	deltasth=athvert_alphasigma2s
470
471	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
472	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
473	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
474	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
475	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
476	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
477
478	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
479	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
480	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
481
482	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)exp(-1.xenv2**2)
483	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
484	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
485	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))(exp(-1.xenv1*2)-exp(-1.xenv2**2))
486
487	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
488	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
489	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
490
491	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
492	! qlenv(ind1,ind2)=IntJ
493	! print*, qlenv(ind1,ind2),'VERIF EAU'
494
495	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)exp(-1.xth2**2)
496	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
497	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
498	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))(exp(-1.xth1*2)-exp(-1.xth2**2))
499
500	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
501	! qlth(ind1,ind2)=IntJ
502	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
503	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
504
505
506
507
508	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
509
510	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
511
512	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
513	ctot(ind1,ind2)=0.
514	qcloud(ind1)=zqsatenv(ind1,ind2)
515
516	else
517
518	ctot(ind1,ind2)=ctot(ind1,ind2)
519	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
520	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
521	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
522
523	endif
524
525
526
527	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
528
529
530	else ! gaussienne environnement seule
531
532	zqenv(ind1)=po(ind1)
533	Tbef=t(ind1,ind2)
534	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
535	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
536	qsatbef=MIN(0.5,qsatbef)
537	zcor=1./(1.-retv*qsatbef)
538	qsatbef=qsatbef*zcor
539	zqsatenv(ind1,ind2)=qsatbef
540
541
542	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
543	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
544	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
545	aenv=1./(1.+(alenv*Lv/cppd))
546	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
547
548
549	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
550
551	sqrt2pi=sqrt(2.*pi)
552	xenv=senv/(sqrt(2.)*sigma1s)
553	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
554	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
555
556	if (ctot(ind1,ind2).lt.1.e-3) then
557	ctot(ind1,ind2)=0.
558	qcloud(ind1)=zqsatenv(ind1,ind2)
559
560	else
561
562	ctot(ind1,ind2)=ctot(ind1,ind2)
563	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
564
565	endif
566
567
568
569
570
571
572	endif
573	enddo
574
575	return
576	! end
577	END SUBROUTINE cloudth_vert
578
579
580
581
582	SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
583	& ztv,po,zqta,fraca, &
584	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
585	& ratqs,zqs,t, &
586	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
587
588	use lmdz_cloudth_ini, only: iflag_cloudth_vert
589
590	IMPLICIT NONE
591
592
593	!===========================================================================
594	! Author : Arnaud Octavio Jam (LMD/CNRS)
595	! Date : 25 Mai 2010
596	! Objet : calcule les valeurs de qc et rneb dans les thermiques
597	!===========================================================================
598
599
600	INCLUDE "YOMCST.h"
601	INCLUDE "YOETHF.h"
602	INCLUDE "FCTTRE.h"
603
604	INTEGER itap,ind1,ind2
605	INTEGER ngrid,klev,klon,l,ig
606	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
607
608	REAL ztv(ngrid,klev)
609	REAL po(ngrid)
610	REAL zqenv(ngrid)
611	REAL zqta(ngrid,klev)
612
613	REAL fraca(ngrid,klev+1)
614	REAL zpspsk(ngrid,klev)
615	REAL paprs(ngrid,klev+1)
616	REAL pplay(ngrid,klev)
617	REAL ztla(ngrid,klev)
618	REAL zthl(ngrid,klev)
619
620	REAL zqsatth(ngrid,klev)
621	REAL zqsatenv(ngrid,klev)
622
623	REAL sigma1(ngrid,klev)
624	REAL sigma2(ngrid,klev)
625	REAL qlth(ngrid,klev)
626	REAL qlenv(ngrid,klev)
627	REAL qltot(ngrid,klev)
628	REAL cth(ngrid,klev)
629	REAL cenv(ngrid,klev)
630	REAL ctot(ngrid,klev)
631	REAL cth_vol(ngrid,klev)
632	REAL cenv_vol(ngrid,klev)
633	REAL ctot_vol(ngrid,klev)
634	REAL rneb(ngrid,klev)
635	REAL t(ngrid,klev)
636	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
637	REAL rdd,cppd,Lv
638	REAL alth,alenv,ath,aenv
639	REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
640	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
641	REAL Tbef,zdelta,qsatbef,zcor
642	REAL qlbef
643	REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
644	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
645	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
646	REAL zqs(ngrid), qcloud(ngrid)
647	REAL erf
648
649
650	IF (iflag_cloudth_vert.GE.1) THEN
651	CALL cloudth_vert_v3(ngrid,klev,ind2, &
652	& ztv,po,zqta,fraca, &
653	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
654	& ratqs,zqs,t, &
655	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
656	RETURN
657	ENDIF
658	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
659
660
661	!-------------------------------------------------------------------------------
662	! Initialisation des variables r?elles
663	!-------------------------------------------------------------------------------
664	sigma1(:,:)=0.
665	sigma2(:,:)=0.
666	qlth(:,:)=0.
667	qlenv(:,:)=0.
668	qltot(:,:)=0.
669	rneb(:,:)=0.
670	qcloud(:)=0.
671	cth(:,:)=0.
672	cenv(:,:)=0.
673	ctot(:,:)=0.
674	cth_vol(:,:)=0.
675	cenv_vol(:,:)=0.
676	ctot_vol(:,:)=0.
677	qsatmmussig1=0.
678	qsatmmussig2=0.
679	rdd=287.04
680	cppd=1005.7
681	pi=3.14159
682	Lv=2.5e6
683	sqrt2pi=sqrt(2.*pi)
684	sqrt2=sqrt(2.)
685	sqrtpi=sqrt(pi)
686
687
688	!-------------------------------------------------------------------------------
689	! Cloud fraction in the thermals and standard deviation of the PDFs
690	!-------------------------------------------------------------------------------
691	do ind1=1,ngrid
692
693	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
694
695	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
696
697	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
698	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
699	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
700	qsatbef=MIN(0.5,qsatbef)
701	zcor=1./(1.-retv*qsatbef)
702	qsatbef=qsatbef*zcor
703	zqsatenv(ind1,ind2)=qsatbef
704
705
706	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
707	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
708	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
709
710	!po = qt de l'environnement ET des thermique
711	!zqenv = qt environnement
712	!zqsatenv = qsat environnement
713	!zthl = Tl environnement
714
715
716	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
717	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
718	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
719	qsatbef=MIN(0.5,qsatbef)
720	zcor=1./(1.-retv*qsatbef)
721	qsatbef=qsatbef*zcor
722	zqsatth(ind1,ind2)=qsatbef
723
724	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
725	ath=1./(1.+(alth*Lv/cppd)) !al, p84
726	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
727
728	!zqta = qt thermals
729	!zqsatth = qsat thermals
730	!ztla = Tl thermals
731
732	!------------------------------------------------------------------------------
733	! s standard deviations
734	!------------------------------------------------------------------------------
735
736	! tests
737	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
738	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
739	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
740	! final option
741	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
742	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
743
744	!------------------------------------------------------------------------------
745	! Condensed water and cloud cover
746	!------------------------------------------------------------------------------
747	xth=sth/(sqrt2*sigma2s)
748	xenv=senv/(sqrt2*sigma1s)
749	cth(ind1,ind2)=0.5(1.+1.erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
750	cenv(ind1,ind2)=0.5(1.+1.erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
751	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
752	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
753
754	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth(ind1,ind2))
755	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
756	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
757
758	if (ctot(ind1,ind2).lt.1.e-10) then
759	ctot(ind1,ind2)=0.
760	qcloud(ind1)=zqsatenv(ind1,ind2)
761	else
762	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
763	endif
764
765	else ! Environnement only, follow the if l.110
766
767	zqenv(ind1)=po(ind1)
768	Tbef=t(ind1,ind2)
769	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
770	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
771	qsatbef=MIN(0.5,qsatbef)
772	zcor=1./(1.-retv*qsatbef)
773	qsatbef=qsatbef*zcor
774	zqsatenv(ind1,ind2)=qsatbef
775
776	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
777	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
778	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
779	aenv=1./(1.+(alenv*Lv/cppd))
780	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
781
782	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
783
784	xenv=senv/(sqrt2*sigma1s)
785	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
786	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
787	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
788
789	if (ctot(ind1,ind2).lt.1.e-3) then
790	ctot(ind1,ind2)=0.
791	qcloud(ind1)=zqsatenv(ind1,ind2)
792	else
793	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
794	endif
795
796
797	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
798	enddo ! from the loop on ngrid l.108
799	return
800	! end
801	END SUBROUTINE cloudth_v3
802
803
804
805	!===========================================================================
806	SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
807	& ztv,po,zqta,fraca, &
808	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
809	& ratqs,zqs,t, &
810	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
811
812	!===========================================================================
813	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
814	! Date : 25 Mai 2010
815	! Objet : calcule les valeurs de qc et rneb dans les thermiques
816	!===========================================================================
817
818	use lmdz_cloudth_ini, only : iflag_cloudth_vert,iflag_ratqs
819	use lmdz_cloudth_ini, only : vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs
820
821	IMPLICIT NONE
822
823
824
825	INCLUDE "YOMCST.h"
826	INCLUDE "YOETHF.h"
827	INCLUDE "FCTTRE.h"
828
829	INTEGER itap,ind1,ind2
830	INTEGER ngrid,klev,klon,l,ig
831	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
832
833	REAL ztv(ngrid,klev)
834	REAL po(ngrid)
835	REAL zqenv(ngrid)
836	REAL zqta(ngrid,klev)
837
838	REAL fraca(ngrid,klev+1)
839	REAL zpspsk(ngrid,klev)
840	REAL paprs(ngrid,klev+1)
841	REAL pplay(ngrid,klev)
842	REAL ztla(ngrid,klev)
843	REAL zthl(ngrid,klev)
844
845	REAL zqsatth(ngrid,klev)
846	REAL zqsatenv(ngrid,klev)
847
848	REAL sigma1(ngrid,klev)
849	REAL sigma2(ngrid,klev)
850	REAL qlth(ngrid,klev)
851	REAL qlenv(ngrid,klev)
852	REAL qltot(ngrid,klev)
853	REAL cth(ngrid,klev)
854	REAL cenv(ngrid,klev)
855	REAL ctot(ngrid,klev)
856	REAL cth_vol(ngrid,klev)
857	REAL cenv_vol(ngrid,klev)
858	REAL ctot_vol(ngrid,klev)
859	REAL rneb(ngrid,klev)
860	REAL t(ngrid,klev)
861	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
862	REAL rdd,cppd,Lv
863	REAL alth,alenv,ath,aenv
864	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
865	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
866	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
867	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
868	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
869	REAL Tbef,zdelta,qsatbef,zcor
870	REAL qlbef
871	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
872	! Change the width of the PDF used for vertical subgrid scale heterogeneity
873	! (J Jouhaud, JL Dufresne, JB Madeleine)
874
875	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
876	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
877	REAL zqs(ngrid), qcloud(ngrid)
878	REAL erf
879
880	REAL rhodz(ngrid,klev)
881	REAL zrho(ngrid,klev)
882	REAL dz(ngrid,klev)
883
884	DO ind1 = 1, ngrid
885	!Layer calculation
886	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
887	zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
888	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
889	END DO
890
891	!------------------------------------------------------------------------------
892	! Initialize
893	!------------------------------------------------------------------------------
894
895	sigma1(:,:)=0.
896	sigma2(:,:)=0.
897	qlth(:,:)=0.
898	qlenv(:,:)=0.
899	qltot(:,:)=0.
900	rneb(:,:)=0.
901	qcloud(:)=0.
902	cth(:,:)=0.
903	cenv(:,:)=0.
904	ctot(:,:)=0.
905	cth_vol(:,:)=0.
906	cenv_vol(:,:)=0.
907	ctot_vol(:,:)=0.
908	qsatmmussig1=0.
909	qsatmmussig2=0.
910	rdd=287.04
911	cppd=1005.7
912	pi=3.14159
913	Lv=2.5e6
914	sqrt2pi=sqrt(2.*pi)
915	sqrt2=sqrt(2.)
916	sqrtpi=sqrt(pi)
917
918
919
920	!-------------------------------------------------------------------------------
921	! Calcul de la fraction du thermique et des ecart-types des distributions
922	!-------------------------------------------------------------------------------
923	do ind1=1,ngrid
924
925	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
926
927	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
928
929
930	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
931	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
932	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
933	qsatbef=MIN(0.5,qsatbef)
934	zcor=1./(1.-retv*qsatbef)
935	qsatbef=qsatbef*zcor
936	zqsatenv(ind1,ind2)=qsatbef
937
938
939	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
940	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
941	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
942
943	!zqenv = qt environnement
944	!zqsatenv = qsat environnement
945	!zthl = Tl environnement
946
947
948	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
949	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
950	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
951	qsatbef=MIN(0.5,qsatbef)
952	zcor=1./(1.-retv*qsatbef)
953	qsatbef=qsatbef*zcor
954	zqsatth(ind1,ind2)=qsatbef
955
956	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
957	ath=1./(1.+(alth*Lv/cppd)) !al, p84
958	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
959
960
961	!zqta = qt thermals
962	!zqsatth = qsat thermals
963	!ztla = Tl thermals
964	!------------------------------------------------------------------------------
965	! s standard deviation
966	!------------------------------------------------------------------------------
967
968	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
969	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
970	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
971	IF (cloudth_ratqsmin>0.) THEN
972	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
973	ELSE
974	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
975	ENDIF
976	sigma1s = sigma1s_fraca + sigma1s_ratqs
977	IF (iflag_ratqs.eq.11) then
978	sigma1s = ratqs(ind1,ind2)po(ind1)aenv
979	ENDIF
980	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
981	! tests
982	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
983	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
984	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
985	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
986	! if (paprs(ind1,ind2).gt.90000) then
987	! ratqs(ind1,ind2)=0.002
988	! else
989	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
990	! endif
991	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
992	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
993	! sigma1s=ratqs(ind1,ind2)*po(ind1)
994	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
995
996	IF (iflag_cloudth_vert == 1) THEN
997	!-------------------------------------------------------------------------------
998	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
999	!-------------------------------------------------------------------------------
1000
1001	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1002	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1003
1004	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1005	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1006	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1007	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1008	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
1009	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
1010
1011	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
1012	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
1013	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1014
1015	! Environment
1016	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
1017	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
1018	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
1019	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
1020
1021	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1022
1023	! Thermal
1024	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
1025	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
1026	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
1027	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
1028	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1029	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1030
1031	ELSE IF (iflag_cloudth_vert >= 3) THEN
1032	IF (iflag_cloudth_vert < 5) THEN
1033	!-------------------------------------------------------------------------------
1034	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1035	!-------------------------------------------------------------------------------
1036	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
1037	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
1038	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1039	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1040	IF (iflag_cloudth_vert == 3) THEN
1041	deltasenv=aenvvert_alphasigma1s
1042	deltasth=athvert_alpha_thsigma2s
1043	ELSE IF (iflag_cloudth_vert == 4) THEN
1044	IF (iflag_cloudth_vert_noratqs == 1) THEN
1045	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1046	deltasth=vert_alpha_th*sigma2s
1047	ELSE
1048	deltasenv=vert_alpha*sigma1s
1049	deltasth=vert_alpha_th*sigma2s
1050	ENDIF
1051	ENDIF
1052
1053	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1054	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1055	exp_xenv1 = exp(-1.xenv1*2)
1056	exp_xenv2 = exp(-1.xenv2*2)
1057	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1058	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1059	exp_xth1 = exp(-1.xth1*2)
1060	exp_xth2 = exp(-1.xth2*2)
1061
1062	!CF_surfacique
1063	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1064	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1065	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1066
1067
1068	!CF_volumique & eau condense
1069	!environnement
1070	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1071	IntJ_CF=0.5(1.-1.erf(xenv2))
1072	if (deltasenv .lt. 1.e-10) then
1073	qlenv(ind1,ind2)=IntJ
1074	cenv_vol(ind1,ind2)=IntJ_CF
1075	else
1076	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1077	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1078	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1079	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1080	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1081	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1082	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1083	endif
1084
1085	!thermique
1086	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1087	IntJ_CF=0.5(1.-1.erf(xth2))
1088	if (deltasth .lt. 1.e-10) then
1089	qlth(ind1,ind2)=IntJ
1090	cth_vol(ind1,ind2)=IntJ_CF
1091	else
1092	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1093	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1094	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1095	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1096	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1097	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1098	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1099	endif
1100
1101	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1102	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1103
1104	ELSE IF (iflag_cloudth_vert == 5) THEN
1105	sigma1s=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
1106	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1107	+ratqs(ind1,ind2)*po(ind1) !Environment
1108	sigma2s=(0.03218+0.000092655dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)0.5)(((sth-senv)2)0.5)+0.002*zqta(ind1,ind2) !Thermals
1109	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1110	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1111	xth=sth/(sqrt(2.)*sigma2s)
1112	xenv=senv/(sqrt(2.)*sigma1s)
1113
1114	!Volumique
1115	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1116	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1117	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1118	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1119
1120	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1121	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1122	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1123
1124	!Surfacique
1125	!Neggers
1126	!beta=0.0044
1127	!inverse_rho=1.+beta*dz(ind1,ind2)
1128	!print *,'jeanjean : beta=',beta
1129	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1130	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1131	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1132
1133	!Brooks
1134	a_Brooks=0.6694
1135	b_Brooks=0.1882
1136	A_Maj_Brooks=0.1635 !-- sans shear
1137	!A_Maj_Brooks=0.17 !-- ARM LES
1138	!A_Maj_Brooks=0.18 !-- RICO LES
1139	!A_Maj_Brooks=0.19 !-- BOMEX LES
1140	Dx_Brooks=200000.
1141	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1142	!print *,'jeanjean_f=',f_Brooks
1143
1144	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1145	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1146	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1147	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1148
1149
1150
1151
1152
1153	ENDIF ! of if (iflag_cloudth_vert<5)
1154	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1155
1156	! if (ctot(ind1,ind2).lt.1.e-10) then
1157	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1158	ctot(ind1,ind2)=0.
1159	ctot_vol(ind1,ind2)=0.
1160	qcloud(ind1)=zqsatenv(ind1,ind2)
1161
1162	else
1163
1164	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1165	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1166	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1167
1168	endif
1169
1170	else ! gaussienne environnement seule
1171
1172
1173	zqenv(ind1)=po(ind1)
1174	Tbef=t(ind1,ind2)
1175	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1176	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1177	qsatbef=MIN(0.5,qsatbef)
1178	zcor=1./(1.-retv*qsatbef)
1179	qsatbef=qsatbef*zcor
1180	zqsatenv(ind1,ind2)=qsatbef
1181
1182
1183	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1184	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1185	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1186	aenv=1./(1.+(alenv*Lv/cppd))
1187	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1188	sth=0.
1189
1190
1191	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1192	sigma2s=0.
1193
1194	sqrt2pi=sqrt(2.*pi)
1195	xenv=senv/(sqrt(2.)*sigma1s)
1196	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1197	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1198	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1199
1200	if (ctot(ind1,ind2).lt.1.e-3) then
1201	ctot(ind1,ind2)=0.
1202	qcloud(ind1)=zqsatenv(ind1,ind2)
1203
1204	else
1205
1206	! ctot(ind1,ind2)=ctot(ind1,ind2)
1207	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1208
1209	endif
1210
1211
1212
1213
1214	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1215	! Outputs used to check the PDFs
1216	cloudth_senv(ind1,ind2) = senv
1217	cloudth_sth(ind1,ind2) = sth
1218	cloudth_sigmaenv(ind1,ind2) = sigma1s
1219	cloudth_sigmath(ind1,ind2) = sigma2s
1220
1221	enddo ! from the loop on ngrid l.333
1222	return
1223	! end
1224	END SUBROUTINE cloudth_vert_v3
1225	!
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1238	& ztv,po,zqta,fraca, &
1239	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1240	& ratqs,zqs,T, &
1241	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1242
1243	use lmdz_cloudth_ini, only: iflag_cloudth_vert
1244
1245	IMPLICIT NONE
1246
1247
1248	INCLUDE "YOMCST.h"
1249	INCLUDE "YOETHF.h"
1250	INCLUDE "FCTTRE.h"
1251
1252
1253	!Domain variables
1254	INTEGER ngrid !indice Max lat-lon
1255	INTEGER klev !indice Max alt
1256	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1257	INTEGER ind1 !indice in [1:ngrid]
1258	INTEGER ind2 !indice in [1:klev]
1259	!thermal plume fraction
1260	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1261	!temperatures
1262	REAL T(ngrid,klev) !temperature
1263	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1264	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1265	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1266	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1267	!pressure
1268	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1269	REAL pplay(ngrid,klev) !pressure at the middle of the level
1270	!humidity
1271	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1272	REAL po(ngrid) !total water (qt)
1273	REAL zqenv(ngrid) !total water in the environment (qt_env)
1274	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1275	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1276	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1277	REAL qlth(ngrid,klev) !condensed water in the thermals
1278	REAL qlenv(ngrid,klev) !condensed water in the environment
1279	REAL qltot(ngrid,klev) !condensed water in the gridbox
1280	!cloud fractions
1281	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1282	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1283	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1284	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1285	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1286	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1287	!PDF of saturation deficit variables
1288	REAL rdd,cppd,Lv
1289	REAL Tbef,zdelta,qsatbef,zcor
1290	REAL alth,alenv,ath,aenv
1291	REAL sth,senv !saturation deficits in the thermals and environment
1292	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1293	!cloud fraction variables
1294	REAL xth,xenv
1295	REAL inverse_rho,beta !Neggers et al. (2011) method
1296	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1297	!Incloud total water variables
1298	REAL zqs(ngrid) !q_sat
1299	REAL qcloud(ngrid) !eau totale dans le nuage
1300	!Some arithmetic variables
1301	REAL erf,pi,sqrt2,sqrt2pi
1302	!Depth of the layer
1303	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1304	REAL rhodz(ngrid,klev)
1305	REAL zrho(ngrid,klev)
1306	DO ind1 = 1, ngrid
1307	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1308	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1309	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1310	END DO
1311
1312	!------------------------------------------------------------------------------
1313	! Initialization
1314	!------------------------------------------------------------------------------
1315	qlth(:,:)=0.
1316	qlenv(:,:)=0.
1317	qltot(:,:)=0.
1318	cth_vol(:,:)=0.
1319	cenv_vol(:,:)=0.
1320	ctot_vol(:,:)=0.
1321	cth_surf(:,:)=0.
1322	cenv_surf(:,:)=0.
1323	ctot_surf(:,:)=0.
1324	qcloud(:)=0.
1325	rdd=287.04
1326	cppd=1005.7
1327	pi=3.14159
1328	Lv=2.5e6
1329	sqrt2=sqrt(2.)
1330	sqrt2pi=sqrt(2.*pi)
1331
1332
1333	DO ind1=1,ngrid
1334	!-------------------------------------------------------------------------------
1335	!Both thermal and environment in the gridbox
1336	!-------------------------------------------------------------------------------
1337	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1338	!--------------------------------------------
1339	!calcul de qsat_env
1340	!--------------------------------------------
1341	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1342	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1343	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1344	qsatbef=MIN(0.5,qsatbef)
1345	zcor=1./(1.-retv*qsatbef)
1346	qsatbef=qsatbef*zcor
1347	zqsatenv(ind1,ind2)=qsatbef
1348	!--------------------------------------------
1349	!calcul de s_env
1350	!--------------------------------------------
1351	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1352	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1353	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1354	!--------------------------------------------
1355	!calcul de qsat_th
1356	!--------------------------------------------
1357	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1358	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1359	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1360	qsatbef=MIN(0.5,qsatbef)
1361	zcor=1./(1.-retv*qsatbef)
1362	qsatbef=qsatbef*zcor
1363	zqsatth(ind1,ind2)=qsatbef
1364	!--------------------------------------------
1365	!calcul de s_th
1366	!--------------------------------------------
1367	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1368	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1369	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1370	!--------------------------------------------
1371	!calcul standard deviations bi-Gaussian PDF
1372	!--------------------------------------------
1373	sigma_th=(0.03218+0.000092655dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)0.5)(((sth-senv)2)0.5)+0.002*zqta(ind1,ind2)
1374	sigma_env=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
1375	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1376	+ratqs(ind1,ind2)*po(ind1)
1377	xth=sth/(sqrt2*sigma_th)
1378	xenv=senv/(sqrt2*sigma_env)
1379	!--------------------------------------------
1380	!Cloud fraction by volume CF_vol
1381	!--------------------------------------------
1382	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1383	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1384	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1385	!--------------------------------------------
1386	!Condensed water qc
1387	!--------------------------------------------
1388	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1389	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1390	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1391	!--------------------------------------------
1392	!Cloud fraction by surface CF_surf
1393	!--------------------------------------------
1394	!Method Neggers et al. (2011) : ok for cumulus clouds only
1395	!beta=0.0044 (Jouhaud et al.2018)
1396	!inverse_rho=1.+beta*dz(ind1,ind2)
1397	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1398	!Method Brooks et al. (2005) : ok for all types of clouds
1399	a_Brooks=0.6694
1400	b_Brooks=0.1882
1401	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1402	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1403	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1404	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1405	!--------------------------------------------
1406	!Incloud Condensed water qcloud
1407	!--------------------------------------------
1408	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1409	ctot_vol(ind1,ind2)=0.
1410	ctot_surf(ind1,ind2)=0.
1411	qcloud(ind1)=zqsatenv(ind1,ind2)
1412	else
1413	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1414	endif
1415
1416
1417
1418	!-------------------------------------------------------------------------------
1419	!Environment only in the gridbox
1420	!-------------------------------------------------------------------------------
1421	ELSE
1422	!--------------------------------------------
1423	!calcul de qsat_env
1424	!--------------------------------------------
1425	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1426	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1427	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1428	qsatbef=MIN(0.5,qsatbef)
1429	zcor=1./(1.-retv*qsatbef)
1430	qsatbef=qsatbef*zcor
1431	zqsatenv(ind1,ind2)=qsatbef
1432	!--------------------------------------------
1433	!calcul de s_env
1434	!--------------------------------------------
1435	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1436	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1437	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1438	!--------------------------------------------
1439	!calcul standard deviations Gaussian PDF
1440	!--------------------------------------------
1441	zqenv(ind1)=po(ind1)
1442	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1443	xenv=senv/(sqrt2*sigma_env)
1444	!--------------------------------------------
1445	!Cloud fraction by volume CF_vol
1446	!--------------------------------------------
1447	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1448	!--------------------------------------------
1449	!Condensed water qc
1450	!--------------------------------------------
1451	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1452	!--------------------------------------------
1453	!Cloud fraction by surface CF_surf
1454	!--------------------------------------------
1455	!Method Neggers et al. (2011) : ok for cumulus clouds only
1456	!beta=0.0044 (Jouhaud et al.2018)
1457	!inverse_rho=1.+beta*dz(ind1,ind2)
1458	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1459	!Method Brooks et al. (2005) : ok for all types of clouds
1460	a_Brooks=0.6694
1461	b_Brooks=0.1882
1462	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1463	Dx_Brooks=200000.
1464	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1465	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1466	!--------------------------------------------
1467	!Incloud Condensed water qcloud
1468	!--------------------------------------------
1469	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1470	ctot_vol(ind1,ind2)=0.
1471	ctot_surf(ind1,ind2)=0.
1472	qcloud(ind1)=zqsatenv(ind1,ind2)
1473	else
1474	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1475	endif
1476
1477
1478	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1479
1480	! Outputs used to check the PDFs
1481	cloudth_senv(ind1,ind2) = senv
1482	cloudth_sth(ind1,ind2) = sth
1483	cloudth_sigmaenv(ind1,ind2) = sigma_env
1484	cloudth_sigmath(ind1,ind2) = sigma_th
1485
1486	END DO ! From the loop on ngrid
1487	return
1488
1489	END SUBROUTINE cloudth_v6
1490
1491
1492
1493
1494
1495	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1496	SUBROUTINE cloudth_mpc(klon,klev,ind2,iflag_mpc_bl,mpc_bl_points, &
1497	& temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, &
1498	& ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol, &
1499	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1500	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1501	! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS)
1502	! Date: Adapted from cloudth_vert_v3 in 2021
1503	! Aim : computes qc and rneb in thermals with cold microphysical considerations
1504	! + for mixed phase boundary layer clouds, calculate ql and qi from
1505	! a stationary MPC model
1506	! IMPORTANT NOTE: we assume iflag_clouth_vert=3
1507	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1508
1509	use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
1510	use lmdz_cloudth_ini, only: C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
1511	use lmdz_lscp_tools, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY
1512
1513	use phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth
1514
1515	IMPLICIT NONE
1516
1517	INCLUDE "YOMCST.h"
1518	INCLUDE "YOETHF.h"
1519	INCLUDE "FCTTRE.h"
1520
1521
1522	!------------------------------------------------------------------------------
1523	! Declaration
1524	!------------------------------------------------------------------------------
1525
1526	! INPUT/OUTPUT
1527
1528	INTEGER, INTENT(IN) :: klon,klev,ind2
1529
1530
1531	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K]
1532	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! Virtual potential temp [K]
1533	REAL, DIMENSION(klon), INTENT(IN) :: po ! specific humidity [kg/kg]
1534	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg]
1535	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: fraca ! Fraction of the mesh covered by thermals [0-1]
1536	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk
1537	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! Pressure at layer interfaces [Pa]
1538	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! Pressure at the center of layers [Pa]
1539	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! Liquid temp [K]
1540	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! Liquid potential temp [K]
1541	REAL, DIMENSION(klon,klev), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
1542	REAL, DIMENSION(klon), INTENT(IN) :: zqs ! Saturation specific humidity in the mesh [kg/kg]
1543	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s]
1544	INTEGER, INTENT(IN) :: iflag_mpc_bl ! option flag for mpc boundary layer clouds param.
1545
1546
1547	INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
1548
1549	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
1550	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
1551	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
1552	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
1553	REAL, DIMENSION(klon,klev), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
1554	real, dimension(klon,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1555
1556
1557	! LOCAL VARIABLES
1558
1559	INTEGER itap,ind1,l,ig,iter,k
1560	INTEGER iflag_topthermals, niter
1561	LOGICAL falseklon(klon)
1562
1563	REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
1564	REAL sigma1(klon,klev)
1565	REAL sigma2(klon,klev)
1566	REAL qcth(klon,klev)
1567	REAL qcenv(klon,klev)
1568	REAL qctot(klon,klev)
1569	REAL cth(klon,klev)
1570	REAL cenv(klon,klev)
1571	REAL cth_vol(klon,klev)
1572	REAL cenv_vol(klon,klev)
1573	REAL rneb(klon,klev)
1574	REAL zqenv(klon), zqth(klon), zqenvonly(klon)
1575	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
1576	REAL rdd,cppd,Lv
1577	REAL alth,alenv,ath,aenv
1578	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1579	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1580	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1581	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1582	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1583	REAL zdelta,qsatbef,zcor
1584	REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon)
1585	REAL qlbef
1586	REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon)
1587	REAL erf
1588	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1589	REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1590	REAL rhodz(klon,klev)
1591	REAL zrho(klon,klev)
1592	REAL dz(klon,klev)
1593	REAL qslenv(klon), qslth(klon)
1594	REAL alenvl, aenvl
1595	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1596	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
1597	REAL qmax, ttarget, stmp, cout, coutref, fraci
1598	REAL maxi, mini, pas, temp_lim
1599	REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
1600
1601	! Modifty the saturation deficit PDF in thermals
1602	! in the presence of ice crystals
1603	CHARACTER (len = 80) :: abort_message
1604	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
1605
1606
1607	!------------------------------------------------------------------------------
1608	! Initialisation
1609	!------------------------------------------------------------------------------
1610
1611
1612	! Few initial checksS
1613
1614	IF (iflag_cloudth_vert.NE.3) THEN
1615	abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3'
1616	CALL abort_physic(routname,abort_message,1)
1617	ENDIF
1618
1619	DO k = 1,klev
1620	DO ind1 = 1, klon
1621	rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2
1622	zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3
1623	dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre
1624	END DO
1625	END DO
1626
1627
1628	sigma1(:,:)=0.
1629	sigma2(:,:)=0.
1630	qcth(:,:)=0.
1631	qcenv(:,:)=0.
1632	qctot(:,:)=0.
1633	qlth(:,ind2)=0.
1634	qith(:,ind2)=0.
1635	wiceth(:,ind2)=0.
1636	rneb(:,:)=0.
1637	qcloud(:)=0.
1638	cth(:,:)=0.
1639	cenv(:,:)=0.
1640	ctot(:,:)=0.
1641	cth_vol(:,:)=0.
1642	cenv_vol(:,:)=0.
1643	ctot_vol(:,:)=0.
1644	falseklon(:)=.false.
1645	qsatmmussig1=0.
1646	qsatmmussig2=0.
1647	rdd=287.04
1648	cppd=1005.7
1649	pi=3.14159
1650	sqrt2pi=sqrt(2.*pi)
1651	sqrt2=sqrt(2.)
1652	sqrtpi=sqrt(pi)
1653	icefrac(:,ind2)=0.
1654	althl=0.
1655	althi=0.
1656	athl=0.
1657	athi=0.
1658	senvl=0.
1659	senvi=0.
1660	sthi=0.
1661	sthl=0.
1662	sthil=0.
1663
1664	!-------------------------------------------------------------------------------
1665	! Identify grid points with potential mixed-phase conditions
1666	!-------------------------------------------------------------------------------
1667
1668	temp_lim=RTT-40.0
1669
1670	DO ind1=1,klon
1671	IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) &
1672	.AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2) &
1673	.AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
1674	mpc_bl_points(ind1,ind2)=1
1675	ELSE
1676	mpc_bl_points(ind1,ind2)=0
1677	ENDIF
1678	ENDDO
1679
1680
1681	!-------------------------------------------------------------------------------
1682	! Thermal fraction calculation and standard deviation of the distribution
1683	!-------------------------------------------------------------------------------
1684
1685	! calculation of temperature, humidity and saturation specific humidity
1686
1687	Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
1688	Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
1689	Tbefenvonly(:)=temp(:,ind2)
1690	rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD
1691
1692	zqenv(:)=(po(:)-fraca(:,ind2)zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = aqtth + (1-a)*qtenv
1693	zqth(:)=zqta(:,ind2)
1694	zqenvonly(:)=po(:)
1695
1696
1697	CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
1698
1699	CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
1700	CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
1701
1702	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
1703	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
1704	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
1705
1706
1707	DO ind1=1,klon
1708
1709
1710	IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
1711
1712
1713	! Environment:
1714
1715
1716	IF (Tbefenv(ind1) .GE. RTT) THEN
1717	Lv=RLVTT
1718	ELSE
1719	Lv=RLSTT
1720	ENDIF
1721
1722
1723	alenv=(0.622Lvzqsatenv(ind1))/(rddzthl(ind1,ind2)*2) !qsl, p84
1724	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
1725	senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84
1726
1727	! For MPCs:
1728	IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
1729	alenvl=(0.622RLVTTqslenv(ind1))/(rddzthl(ind1,ind2)*2)
1730	aenvl=1./(1.+(alenv*Lv/cppd))
1731	senvl=aenvl*(po(ind1)-qslenv(ind1))
1732	senvi=senv
1733	ENDIF
1734
1735
1736	! Thermals:
1737
1738
1739	IF (Tbefth(ind1) .GE. RTT) THEN
1740	Lv=RLVTT
1741	ELSE
1742	Lv=RLSTT
1743	ENDIF
1744
1745
1746	alth=(0.622Lvzqsatth(ind1))/(rddztla(ind1,ind2)*2)
1747	ath=1./(1.+(alth*Lv/cppd))
1748	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))
1749
1750	! For MPCs:
1751	IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
1752	althi=alth
1753	althl=(0.622RLVTTqslth(ind1))/(rddztla(ind1,ind2)*2)
1754	athl=1./(1.+(alth*RLVTT/cppd))
1755	athi=alth
1756	sthl=athl*(zqta(ind1,ind2)-qslth(ind1))
1757	sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2))
1758	sthil=athi*(zqta(ind1,ind2)-qslth(ind1))
1759	ENDIF
1760
1761
1762
1763	!-------------------------------------------------------------------------------
1764	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1765	! Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3
1766	!-------------------------------------------------------------------------------
1767
1768	IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
1769
1770	! Standard deviation of the distributions
1771
1772	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
1773	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
1774
1775	IF (cloudth_ratqsmin>0.) THEN
1776	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1777	ELSE
1778	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1779	ENDIF
1780
1781	sigma1s = sigma1s_fraca + sigma1s_ratqs
1782	IF (iflag_ratqs.eq.11) then
1783	sigma1s = ratqs(ind1,ind2)po(ind1)aenv
1784	ENDIF
1785	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1786
1787
1788	deltasenv=aenvvert_alphasigma1s
1789	deltasth=athvert_alpha_thsigma2s
1790
1791	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1792	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1793	exp_xenv1 = exp(-1.xenv1*2)
1794	exp_xenv2 = exp(-1.xenv2*2)
1795	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1796	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1797	exp_xth1 = exp(-1.xth1*2)
1798	exp_xth2 = exp(-1.xth2*2)
1799
1800	!surface CF
1801
1802	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1803	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1804	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1805
1806
1807	!volume CF and condensed water
1808
1809	!environnement
1810
1811	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1812	IntJ_CF=0.5(1.-1.erf(xenv2))
1813
1814	IF (deltasenv .LT. 1.e-10) THEN
1815	qcenv(ind1,ind2)=IntJ
1816	cenv_vol(ind1,ind2)=IntJ_CF
1817	ELSE
1818	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1819	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1820	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1821	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1822	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1823	qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1824	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1825	ENDIF
1826
1827
1828
1829	!thermals
1830
1831	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1832	IntJ_CF=0.5(1.-1.erf(xth2))
1833
1834	IF (deltasth .LT. 1.e-10) THEN
1835	qcth(ind1,ind2)=IntJ
1836	cth_vol(ind1,ind2)=IntJ_CF
1837	ELSE
1838	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1839	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1840	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1841	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1842	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1843	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1844	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1845	ENDIF
1846
1847	qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
1848	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1849
1850
1851	IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
1852	ctot(ind1,ind2)=0.
1853	ctot_vol(ind1,ind2)=0.
1854	qcloud(ind1)=zqsatenv(ind1)
1855	qincloud(ind1)=0.
1856	ELSE
1857	qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1858	qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
1859	ENDIF
1860
1861
1862	ELSE ! mpc_bl_points>0
1863
1864	! Treat boundary layer mixed phase clouds
1865
1866	! thermals
1867	!=========
1868
1869	! ice phase
1870	!...........
1871
1872	qiceth=0.
1873	deltazlev_mpc=dz(ind1,:)
1874	temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
1875	pres_mpc=pplay(ind1,:)
1876	fraca_mpc=fraca(ind1,:)
1877	snowf_mpc=snowflux(ind1,:)
1878	iflag_topthermals=0
1879	IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
1880	iflag_topthermals = 1
1881	ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
1882	.AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
1883	iflag_topthermals = 2
1884	ELSE
1885	iflag_topthermals = 0
1886	ENDIF
1887
1888	CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
1889	qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
1890
1891	! qmax calculation
1892	sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1893	deltasth=athlvert_alpha_thsigma2s
1894	xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
1895	xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
1896	exp_xth1 = exp(-1.xth1*2)
1897	exp_xth2 = exp(-1.xth2*2)
1898	IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1899	IntJ_CF=0.5(1.-1.erf(xth2))
1900	IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1901	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1902	IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1903	IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1904	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1905	qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
1906
1907
1908	! Liquid phase
1909	!................
1910	! We account for the effect of ice crystals in thermals on sthl
1911	! and on the width of the distribution
1912
1913	sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
1914	+ (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
1915
1916	sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
1917	/((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
1918	+0.002*zqta(ind1,ind2)
1919	deltasthc=athlvert_alpha_thsigma2sc
1920
1921
1922	xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
1923	xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
1924	exp_xth1 = exp(-1.xth1*2)
1925	exp_xth2 = exp(-1.xth2*2)
1926	IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
1927	IntJ_CF=0.5(1.-1.erf(xth2))
1928	IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
1929	IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1930	IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
1931	IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
1932	IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
1933	qliqth=IntJ+IntI1+IntI2+IntI3
1934
1935	qlth(ind1,ind2)=MAX(0.,qliqth)
1936
1937	! Condensed water
1938
1939	qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
1940
1941
1942	! consistency with subgrid distribution
1943
1944	IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
1945	fraci=qith(ind1,ind2)/qcth(ind1,ind2)
1946	qcth(ind1,ind2)=qmax
1947	qith(ind1,ind2)=fraci*qmax
1948	qlth(ind1,ind2)=(1.-fraci)*qmax
1949	ENDIF
1950
1951	! Cloud Fraction
1952	!...............
1953	! calculation of qbase which is the value of the water vapor within mixed phase clouds
1954	! such that the total water in cloud = qbase+qliqth+qiceth
1955	! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
1956	! sbase and qbase calculation (note that sbase is wrt liq so negative)
1957	! look for an approximate solution with iteration
1958
1959	ttarget=qcth(ind1,ind2)
1960	mini= athl*(qsith(ind1,ind2)-qslth(ind1))
1961	maxi= 0. !athl(3.sqrt(sigma2s))
1962	niter=20
1963	pas=(maxi-mini)/niter
1964	stmp=mini
1965	sbase=stmp
1966	coutref=1.E6
1967	DO iter=1,niter
1968	cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
1969	+ stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
1970	IF (cout .LT. coutref) THEN
1971	sbase=stmp
1972	coutref=cout
1973	ELSE
1974	stmp=stmp+pas
1975	ENDIF
1976	ENDDO
1977	qbase=MAX(0., sbase/athl+qslth(ind1))
1978
1979	! surface cloud fraction in thermals
1980	cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
1981	cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
1982
1983
1984	!volume cloud fraction in thermals
1985	!to be checked
1986	xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
1987	xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
1988	exp_xth1 = exp(-1.xth1*2)
1989	exp_xth2 = exp(-1.xth2*2)
1990
1991	IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1992	IntJ_CF=0.5(1.-1.erf(xth2))
1993
1994	IF (deltasth .LT. 1.e-10) THEN
1995	cth_vol(ind1,ind2)=IntJ_CF
1996	ELSE
1997	IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1998	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1999	IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
2000	IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
2001	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
2002	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2003	ENDIF
2004	cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
2005
2006
2007
2008	! Environment
2009	!=============
2010	! In the environment/downdrafts, only liquid clouds
2011	! See Shupe et al. 2008, JAS
2012
2013	! standard deviation of the distribution in the environment
2014	sth=sthl
2015	senv=senvl
2016	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
2017	& (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
2018	! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
2019	! in the environement
2020
2021	sigma1s_ratqs=1E-10
2022	IF (cloudth_ratqsmin>0.) THEN
2023	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2024	ENDIF
2025
2026	sigma1s = sigma1s_fraca + sigma1s_ratqs
2027	IF (iflag_ratqs.eq.11) then
2028	sigma1s = ratqs(ind1,ind2)po(ind1)aenv
2029	ENDIF
2030	IF (iflag_ratqs.eq.11) then
2031	sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
2032	ENDIF
2033	deltasenv=aenvlvert_alphasigma1s
2034	xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
2035	xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
2036	exp_xenv1 = exp(-1.xenv1*2)
2037	exp_xenv2 = exp(-1.xenv2*2)
2038
2039	!surface CF
2040	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
2041
2042	!volume CF and condensed water
2043	IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2044	IntJ_CF=0.5(1.-1.erf(xenv2))
2045
2046	IF (deltasenv .LT. 1.e-10) THEN
2047	qcenv(ind1,ind2)=IntJ
2048	cenv_vol(ind1,ind2)=IntJ_CF
2049	ELSE
2050	IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
2051	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
2052	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
2053	IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
2054	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
2055	qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
2056	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2057	ENDIF
2058
2059	qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
2060	cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
2061
2062
2063
2064	! Thermals + environment
2065	! =======================
2066	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
2067	qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
2068	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2069	IF (qcth(ind1,ind2) .GT. 0) THEN
2070	icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
2071	/(fraca(ind1,ind2)*qcth(ind1,ind2) &
2072	+(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
2073	icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
2074	ELSE
2075	icefrac(ind1,ind2)=0.
2076	ENDIF
2077
2078	IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
2079	ctot(ind1,ind2)=0.
2080	ctot_vol(ind1,ind2)=0.
2081	qincloud(ind1)=0.
2082	qcloud(ind1)=zqsatenv(ind1)
2083	ELSE
2084	qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
2085	+(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
2086	qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
2087	+(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
2088	ENDIF
2089
2090	ENDIF ! mpc_bl_points
2091
2092
2093	ELSE ! gaussian for environment only
2094
2095
2096
2097
2098	IF (Tbefenvonly(ind1) .GE. RTT) THEN
2099	Lv=RLVTT
2100	ELSE
2101	Lv=RLSTT
2102	ENDIF
2103
2104
2105	zthl(ind1,ind2)=temp(ind1,ind2)(101325./paprs(ind1,ind2))*(rdd/cppd)
2106	alenv=(0.622Lvzqsatenvonly(ind1))/(rddzthl(ind1,ind2)*2)
2107	aenv=1./(1.+(alenv*Lv/cppd))
2108	senv=aenv*(po(ind1)-zqsatenvonly(ind1))
2109	sth=0.
2110
2111	sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
2112	sigma2s=0.
2113
2114	sqrt2pi=sqrt(2.*pi)
2115	xenv=senv/(sqrt(2.)*sigma1s)
2116	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
2117	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
2118	qctot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
2119
2120	IF (ctot(ind1,ind2).LT.1.e-3) THEN
2121	ctot(ind1,ind2)=0.
2122	qcloud(ind1)=zqsatenvonly(ind1)
2123	qincloud(ind1)=0.
2124	ELSE
2125	qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
2126	qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
2127	ENDIF
2128
2129
2130	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
2131
2132	! Outputs used to check the PDFs
2133	cloudth_senv(ind1,ind2) = senv
2134	cloudth_sth(ind1,ind2) = sth
2135	cloudth_sigmaenv(ind1,ind2) = sigma1s
2136	cloudth_sigmath(ind1,ind2) = sigma2s
2137
2138
2139	ENDDO !loop on klon
2140
2141	! Calcule ice fall velocity in thermals
2142
2143	CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2))
2144
2145	RETURN
2146
2147
2148	END SUBROUTINE cloudth_mpc
2149
2150	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2151	SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,vith,fraca,qith)
2152
2153	! parameterization of ice for boundary
2154	! layer mixed-phase clouds assuming a stationary system
2155	!
2156	! Note that vapor deposition on ice crystals and riming of liquid droplets
2157	! depend on the ice number concentration Ni
2158	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
2159	! and use values from Hong et al. 2004, MWR for instance
2160	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
2161	! One could also think of a more complex expression of Ni;
2162	! function of qi, T, the concentration in aerosols or INP ..
2163	! Here we prefer fixing Ni to a tuning parameter
2164	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
2165	! in Mioche et al. 2017
2166	!
2167	!
2168	! References:
2169	!------------
2170	! This parameterization is thoroughly described in Vignon et al.
2171	!
2172	! More specifically
2173	! for the Water vapor deposition process:
2174	!
2175	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
2176	! ment of stratiform cloudfs and precipitation in large-scale
2177	! models. I: Description and evaluation of the microphysical
2178	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
2179	!
2180	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
2181	! stratiform cloud microphysics scheme in the NCAR Com-
2182	! munity Atmosphere Model (CAM3). Part I: Description and
2183	! numerical tests. J. Climate, 21, 3642–3659
2184	!
2185	! for the Riming process:
2186	!
2187	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
2188	! scale structure and organization of clouds and precipitation in
2189	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
2190	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
2191	!
2192	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
2193	! forecasts of winter precipitation using an improved bulk
2194	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
2195	! forecasts of winter precipitation using an improved bulk
2196	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
2197	!
2198	! For the formation of clouds by thermals:
2199	!
2200	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
2201	! the Atmospheric Sciences, 65, 407–425.
2202	!
2203	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
2204	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
2205	!
2206	!
2207	!
2208	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
2209	!=============================================================================
2210
2211	use phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
2212
2213	IMPLICIT none
2214
2215	INCLUDE "YOMCST.h"
2216
2217	INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
2218	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
2219	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
2220	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
2221	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
2222	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
2223	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
2224	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
2225	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
2226	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
2227	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
2228	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
2229	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
2230	REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
2231	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
2232
2233
2234	INTEGER ind2p1,ind2p2
2235	REAL rho(klev)
2236	REAL unsurtaudet, unsurtaustardep, unsurtaurim
2237	REAL qi, AA, BB, Ka, Dv, rhoi
2238	REAL p0, t0, fp1, fp2
2239	REAL alpha, flux_term
2240	REAL det_term, precip_term, rim_term, dep_term
2241
2242
2243	ind2p1=ind2+1
2244	ind2p2=ind2+2
2245
2246	rho=pres/temp/RD ! air density kg/m3
2247
2248	Ka=2.4e-2 ! thermal conductivity of the air, SI
2249	p0=101325.0 ! ref pressure
2250	T0=273.15 ! ref temp
2251	rhoi=500.0 ! cloud ice density following Reisner et al. 1998
2252	alpha=700. ! fallvelocity param
2253
2254
2255	IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
2256
2257	Dv=0.00010.211(p0/pres(ind2))((temp(ind2)/T0)*1.94) ! water vapor diffusivity in air, SI
2258
2259	! Detrainment term:
2260
2261	unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
2262
2263
2264	! Deposition term
2265	AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
2266	BB=1./(rho(ind2)Dvqsith(ind2))
2267	unsurtaustardep=C_cap(Ni0.66)(qth(ind2)-qsith(ind2))/qsith(ind2) &
2268	4.RPI/(AA+BB)(6.rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
2269
2270	! Riming term neglected at this level
2271	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
2272
2273	qi=fraca(ind2)unsurtaustardep/MAX((d_topunsurtaudet),1E-12)
2274	qi=MAX(qi,0.)**(3./2.)
2275
2276	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
2277
2278	Dv=0.00010.211(p0/pres(ind2p1))((temp(ind2p1)/T0)*1.94) ! water vapor diffusivity in air, SI
2279
2280	! Detrainment term:
2281
2282	unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
2283	det_term=-unsurtaudetqith(ind2p1)rho(ind2p1)
2284
2285
2286	! Deposition term
2287
2288	AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
2289	BB=1./(rho(ind2p1)Dvqsith(ind2p1))
2290	unsurtaustardep=C_cap(Ni0.66)(qth(ind2p1)-qsith(ind2p1)) &
2291	/qsith(ind2p1)4.RPI/(AA+BB) &
2292	(6.rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
2293	dep_term=rho(ind2p1)fraca(ind2p1)(qith(ind2p1)*0.33)unsurtaustardep
2294
2295	! Riming term
2296
2297	unsurtaurim=rho(ind2p1)alpha3./rhoi/2.Eiqlth(ind2p1)((p0/pres(ind2p1))*0.4)
2298	rim_term=rho(ind2p1)fraca(ind2p1)qith(ind2p1)*unsurtaurim
2299
2300	! Precip term
2301
2302	! We assume that there is no solid precipitation outside thermals
2303	! so the precipitation flux within thermals is equal to the precipitation flux
2304	! at mesh-scale divided by thermals fraction
2305
2306	fp2=0.
2307	fp1=0.
2308	IF (fraca(ind2p1) .GT. 0.) THEN
2309	fp2=-qith(ind2p2)rho(ind2p2)vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
2310	fp1=-qith(ind2p1)rho(ind2p1)vith(ind2p1)*fraca(ind2p1)
2311	ENDIF
2312
2313	precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
2314
2315	! Calculation in a top-to-bottom loop
2316
2317	IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
2318	qi= 1./fm_therm(ind1,ind2p1)* &
2319	(deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
2320	fm_therm(ind1,ind2p2)*(qith(ind2p1)))
2321	END IF
2322
2323	ENDIF ! top thermals
2324
2325	qith(ind2)=MAX(0.,qi)
2326
2327	RETURN
2328
2329	END SUBROUTINE ICE_MPC_BL_CLOUDS
2330
2331	END MODULE lmdz_cloudth

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

Download in other formats:

Original Format