Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_cloudth.F90 @ 5451

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

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

Download in other formats:

Original Format