Context Navigation

lmdz_cloudth.F90 @ 5222

Last change on this file since 5222 was 5208, checked in by Laurent Fairhead, 8 weeks ago

Louis d'Alencon modifications:
Il s'agit des modifs permettant de choisir de calculer la variance dans les thermiques par le nouveau modèle (iflag_ratqs = 10) ou de
continuer à calculer cette variance par la param d'Arnaud ce qui fait une version hybride avec variance pronostique dans l'environnement mais
variance diagnostique dans les thermiques (iflag_ratqs = 11). Le iflag_ratqs = 4 version standard continue de faire toutes les variances de façon
diagnostique.

File size: 87.3 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,sigma_qtherm,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,sigma_qtherm
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,sigma_qtherm,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,sigma_qtherm,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),sigma_qtherm(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	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
981	IF (iflag_ratqs.eq.10.or.iflag_ratqs.eq.11) then
982	sigma1s = ratqs(ind1,ind2)po(ind1)aenv
983	IF (iflag_ratqs.eq.10.and.sigma_qtherm(ind1,ind2).ne.0) then
984	sigma2s = sigma_qtherm(ind1,ind2)*ath
985	ENDIF
986	ENDIF
987
988	! tests
989	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
990	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
991	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
992	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
993	! if (paprs(ind1,ind2).gt.90000) then
994	! ratqs(ind1,ind2)=0.002
995	! else
996	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
997	! endif
998	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
999	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1000	! sigma1s=ratqs(ind1,ind2)*po(ind1)
1001	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
1002
1003	IF (iflag_cloudth_vert == 1) THEN
1004	!-------------------------------------------------------------------------------
1005	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1006	!-------------------------------------------------------------------------------
1007
1008	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1009	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1010
1011	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1012	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1013	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1014	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1015	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
1016	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
1017
1018	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
1019	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
1020	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1021
1022	! Environment
1023	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
1024	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
1025	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
1026	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
1027
1028	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1029
1030	! Thermal
1031	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
1032	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
1033	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
1034	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
1035	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1036	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1037
1038	ELSE IF (iflag_cloudth_vert >= 3) THEN
1039	IF (iflag_cloudth_vert < 5) THEN
1040	!-------------------------------------------------------------------------------
1041	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1042	!-------------------------------------------------------------------------------
1043	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
1044	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
1045	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1046	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1047	IF (iflag_cloudth_vert == 3) THEN
1048	deltasenv=aenvvert_alphasigma1s
1049	deltasth=athvert_alpha_thsigma2s
1050	ELSE IF (iflag_cloudth_vert == 4) THEN
1051	IF (iflag_cloudth_vert_noratqs == 1) THEN
1052	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1053	deltasth=vert_alpha_th*sigma2s
1054	ELSE
1055	deltasenv=vert_alpha*sigma1s
1056	deltasth=vert_alpha_th*sigma2s
1057	ENDIF
1058	ENDIF
1059
1060	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1061	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1062	exp_xenv1 = exp(-1.xenv1*2)
1063	exp_xenv2 = exp(-1.xenv2*2)
1064	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1065	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1066	exp_xth1 = exp(-1.xth1*2)
1067	exp_xth2 = exp(-1.xth2*2)
1068
1069	!CF_surfacique
1070	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1071	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1072	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1073
1074
1075	!CF_volumique & eau condense
1076	!environnement
1077	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1078	IntJ_CF=0.5(1.-1.erf(xenv2))
1079	if (deltasenv .lt. 1.e-10) then
1080	qlenv(ind1,ind2)=IntJ
1081	cenv_vol(ind1,ind2)=IntJ_CF
1082	else
1083	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1084	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1085	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1086	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1087	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1088	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1089	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1090	endif
1091
1092	!thermique
1093	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1094	IntJ_CF=0.5(1.-1.erf(xth2))
1095	if (deltasth .lt. 1.e-10) then
1096	qlth(ind1,ind2)=IntJ
1097	cth_vol(ind1,ind2)=IntJ_CF
1098	else
1099	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1100	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1101	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1102	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1103	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1104	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1105	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1106	endif
1107
1108	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1109	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1110
1111	ELSE IF (iflag_cloudth_vert == 5) THEN
1112	sigma1s=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
1113	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1114	+ratqs(ind1,ind2)*po(ind1) !Environment
1115	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
1116	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1117	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1118	xth=sth/(sqrt(2.)*sigma2s)
1119	xenv=senv/(sqrt(2.)*sigma1s)
1120
1121	!Volumique
1122	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1123	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1124	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1125	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1126
1127	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1128	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1129	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1130
1131	!Surfacique
1132	!Neggers
1133	!beta=0.0044
1134	!inverse_rho=1.+beta*dz(ind1,ind2)
1135	!print *,'jeanjean : beta=',beta
1136	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1137	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1138	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1139
1140	!Brooks
1141	a_Brooks=0.6694
1142	b_Brooks=0.1882
1143	A_Maj_Brooks=0.1635 !-- sans shear
1144	!A_Maj_Brooks=0.17 !-- ARM LES
1145	!A_Maj_Brooks=0.18 !-- RICO LES
1146	!A_Maj_Brooks=0.19 !-- BOMEX LES
1147	Dx_Brooks=200000.
1148	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1149	!print *,'jeanjean_f=',f_Brooks
1150
1151	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1152	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1153	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1154	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1155
1156
1157
1158
1159
1160	ENDIF ! of if (iflag_cloudth_vert<5)
1161	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1162
1163	! if (ctot(ind1,ind2).lt.1.e-10) then
1164	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1165	ctot(ind1,ind2)=0.
1166	ctot_vol(ind1,ind2)=0.
1167	qcloud(ind1)=zqsatenv(ind1,ind2)
1168
1169	else
1170
1171	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1172	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1173	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1174
1175	endif
1176
1177	else ! gaussienne environnement seule
1178
1179
1180	zqenv(ind1)=po(ind1)
1181	Tbef=t(ind1,ind2)
1182	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1183	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1184	qsatbef=MIN(0.5,qsatbef)
1185	zcor=1./(1.-retv*qsatbef)
1186	qsatbef=qsatbef*zcor
1187	zqsatenv(ind1,ind2)=qsatbef
1188
1189
1190	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1191	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1192	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1193	aenv=1./(1.+(alenv*Lv/cppd))
1194	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1195	sth=0.
1196
1197
1198	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1199	sigma2s=0.
1200
1201	sqrt2pi=sqrt(2.*pi)
1202	xenv=senv/(sqrt(2.)*sigma1s)
1203	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1204	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1205	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1206
1207	if (ctot(ind1,ind2).lt.1.e-3) then
1208	ctot(ind1,ind2)=0.
1209	qcloud(ind1)=zqsatenv(ind1,ind2)
1210
1211	else
1212
1213	! ctot(ind1,ind2)=ctot(ind1,ind2)
1214	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1215
1216	endif
1217
1218
1219
1220
1221	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1222	! Outputs used to check the PDFs
1223	cloudth_senv(ind1,ind2) = senv
1224	cloudth_sth(ind1,ind2) = sth
1225	cloudth_sigmaenv(ind1,ind2) = sigma1s
1226	cloudth_sigmath(ind1,ind2) = sigma2s
1227
1228	enddo ! from the loop on ngrid l.333
1229	return
1230	! end
1231	END SUBROUTINE cloudth_vert_v3
1232	!
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1245	& ztv,po,zqta,fraca, &
1246	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1247	& ratqs,zqs,T, &
1248	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1249
1250	use lmdz_cloudth_ini, only: iflag_cloudth_vert
1251
1252	IMPLICIT NONE
1253
1254
1255	INCLUDE "YOMCST.h"
1256	INCLUDE "YOETHF.h"
1257	INCLUDE "FCTTRE.h"
1258
1259
1260	!Domain variables
1261	INTEGER ngrid !indice Max lat-lon
1262	INTEGER klev !indice Max alt
1263	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1264	INTEGER ind1 !indice in [1:ngrid]
1265	INTEGER ind2 !indice in [1:klev]
1266	!thermal plume fraction
1267	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1268	!temperatures
1269	REAL T(ngrid,klev) !temperature
1270	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1271	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1272	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1273	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1274	!pressure
1275	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1276	REAL pplay(ngrid,klev) !pressure at the middle of the level
1277	!humidity
1278	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1279	REAL po(ngrid) !total water (qt)
1280	REAL zqenv(ngrid) !total water in the environment (qt_env)
1281	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1282	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1283	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1284	REAL qlth(ngrid,klev) !condensed water in the thermals
1285	REAL qlenv(ngrid,klev) !condensed water in the environment
1286	REAL qltot(ngrid,klev) !condensed water in the gridbox
1287	!cloud fractions
1288	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1289	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1290	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1291	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1292	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1293	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1294	!PDF of saturation deficit variables
1295	REAL rdd,cppd,Lv
1296	REAL Tbef,zdelta,qsatbef,zcor
1297	REAL alth,alenv,ath,aenv
1298	REAL sth,senv !saturation deficits in the thermals and environment
1299	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1300	!cloud fraction variables
1301	REAL xth,xenv
1302	REAL inverse_rho,beta !Neggers et al. (2011) method
1303	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1304	!Incloud total water variables
1305	REAL zqs(ngrid) !q_sat
1306	REAL qcloud(ngrid) !eau totale dans le nuage
1307	!Some arithmetic variables
1308	REAL erf,pi,sqrt2,sqrt2pi
1309	!Depth of the layer
1310	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1311	REAL rhodz(ngrid,klev)
1312	REAL zrho(ngrid,klev)
1313	DO ind1 = 1, ngrid
1314	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1315	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1316	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1317	END DO
1318
1319	!------------------------------------------------------------------------------
1320	! Initialization
1321	!------------------------------------------------------------------------------
1322	qlth(:,:)=0.
1323	qlenv(:,:)=0.
1324	qltot(:,:)=0.
1325	cth_vol(:,:)=0.
1326	cenv_vol(:,:)=0.
1327	ctot_vol(:,:)=0.
1328	cth_surf(:,:)=0.
1329	cenv_surf(:,:)=0.
1330	ctot_surf(:,:)=0.
1331	qcloud(:)=0.
1332	rdd=287.04
1333	cppd=1005.7
1334	pi=3.14159
1335	Lv=2.5e6
1336	sqrt2=sqrt(2.)
1337	sqrt2pi=sqrt(2.*pi)
1338
1339
1340	DO ind1=1,ngrid
1341	!-------------------------------------------------------------------------------
1342	!Both thermal and environment in the gridbox
1343	!-------------------------------------------------------------------------------
1344	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1345	!--------------------------------------------
1346	!calcul de qsat_env
1347	!--------------------------------------------
1348	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1349	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1350	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1351	qsatbef=MIN(0.5,qsatbef)
1352	zcor=1./(1.-retv*qsatbef)
1353	qsatbef=qsatbef*zcor
1354	zqsatenv(ind1,ind2)=qsatbef
1355	!--------------------------------------------
1356	!calcul de s_env
1357	!--------------------------------------------
1358	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1359	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1360	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1361	!--------------------------------------------
1362	!calcul de qsat_th
1363	!--------------------------------------------
1364	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1365	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1366	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1367	qsatbef=MIN(0.5,qsatbef)
1368	zcor=1./(1.-retv*qsatbef)
1369	qsatbef=qsatbef*zcor
1370	zqsatth(ind1,ind2)=qsatbef
1371	!--------------------------------------------
1372	!calcul de s_th
1373	!--------------------------------------------
1374	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1375	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1376	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1377	!--------------------------------------------
1378	!calcul standard deviations bi-Gaussian PDF
1379	!--------------------------------------------
1380	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)
1381	sigma_env=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
1382	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1383	+ratqs(ind1,ind2)*po(ind1)
1384	xth=sth/(sqrt2*sigma_th)
1385	xenv=senv/(sqrt2*sigma_env)
1386	!--------------------------------------------
1387	!Cloud fraction by volume CF_vol
1388	!--------------------------------------------
1389	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1390	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1391	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1392	!--------------------------------------------
1393	!Condensed water qc
1394	!--------------------------------------------
1395	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1396	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1397	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1398	!--------------------------------------------
1399	!Cloud fraction by surface CF_surf
1400	!--------------------------------------------
1401	!Method Neggers et al. (2011) : ok for cumulus clouds only
1402	!beta=0.0044 (Jouhaud et al.2018)
1403	!inverse_rho=1.+beta*dz(ind1,ind2)
1404	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1405	!Method Brooks et al. (2005) : ok for all types of clouds
1406	a_Brooks=0.6694
1407	b_Brooks=0.1882
1408	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1409	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1410	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1411	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1412	!--------------------------------------------
1413	!Incloud Condensed water qcloud
1414	!--------------------------------------------
1415	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1416	ctot_vol(ind1,ind2)=0.
1417	ctot_surf(ind1,ind2)=0.
1418	qcloud(ind1)=zqsatenv(ind1,ind2)
1419	else
1420	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1421	endif
1422
1423
1424
1425	!-------------------------------------------------------------------------------
1426	!Environment only in the gridbox
1427	!-------------------------------------------------------------------------------
1428	ELSE
1429	!--------------------------------------------
1430	!calcul de qsat_env
1431	!--------------------------------------------
1432	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1433	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1434	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1435	qsatbef=MIN(0.5,qsatbef)
1436	zcor=1./(1.-retv*qsatbef)
1437	qsatbef=qsatbef*zcor
1438	zqsatenv(ind1,ind2)=qsatbef
1439	!--------------------------------------------
1440	!calcul de s_env
1441	!--------------------------------------------
1442	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1443	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1444	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1445	!--------------------------------------------
1446	!calcul standard deviations Gaussian PDF
1447	!--------------------------------------------
1448	zqenv(ind1)=po(ind1)
1449	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1450	xenv=senv/(sqrt2*sigma_env)
1451	!--------------------------------------------
1452	!Cloud fraction by volume CF_vol
1453	!--------------------------------------------
1454	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1455	!--------------------------------------------
1456	!Condensed water qc
1457	!--------------------------------------------
1458	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1459	!--------------------------------------------
1460	!Cloud fraction by surface CF_surf
1461	!--------------------------------------------
1462	!Method Neggers et al. (2011) : ok for cumulus clouds only
1463	!beta=0.0044 (Jouhaud et al.2018)
1464	!inverse_rho=1.+beta*dz(ind1,ind2)
1465	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1466	!Method Brooks et al. (2005) : ok for all types of clouds
1467	a_Brooks=0.6694
1468	b_Brooks=0.1882
1469	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1470	Dx_Brooks=200000.
1471	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1472	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1473	!--------------------------------------------
1474	!Incloud Condensed water qcloud
1475	!--------------------------------------------
1476	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1477	ctot_vol(ind1,ind2)=0.
1478	ctot_surf(ind1,ind2)=0.
1479	qcloud(ind1)=zqsatenv(ind1,ind2)
1480	else
1481	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1482	endif
1483
1484
1485	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1486
1487	! Outputs used to check the PDFs
1488	cloudth_senv(ind1,ind2) = senv
1489	cloudth_sth(ind1,ind2) = sth
1490	cloudth_sigmaenv(ind1,ind2) = sigma_env
1491	cloudth_sigmath(ind1,ind2) = sigma_th
1492
1493	END DO ! From the loop on ngrid
1494	return
1495
1496	END SUBROUTINE cloudth_v6
1497
1498
1499
1500
1501
1502	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1503	SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, &
1504	& temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th, &
1505	& ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol, &
1506	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1507	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1508	! Author : Etienne Vignon (LMDZ/CNRS)
1509	! Date: April 2024
1510	! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
1511	! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
1512	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1513
1514	use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
1515	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
1516	use lmdz_lscp_tools, only: CALC_QSAT_ECMWF
1517	use lmdz_lscp_ini, only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
1518
1519	IMPLICIT NONE
1520
1521
1522	!------------------------------------------------------------------------------
1523	! Declaration
1524	!------------------------------------------------------------------------------
1525
1526	! INPUT/OUTPUT
1527
1528	INTEGER, INTENT(IN) :: klon,klev,ind2
1529
1530
1531	REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
1532	REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
1533	REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
1534	REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
1535	REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1]
1536	REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential
1537	REAL, DIMENSION(klon), INTENT(IN) :: paprsup ! Pressure at top layer interface [Pa]
1538	REAL, DIMENSION(klon), INTENT(IN) :: paprsdn ! Pressure at bottom layer interface [Pa]
1539	REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa]
1540	REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
1541
1542
1543	INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
1544
1545	REAL, DIMENSION(klon), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
1546	REAL, DIMENSION(klon), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
1547	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
1548	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
1549	REAL, DIMENSION(klon), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
1550	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
1551	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
1552	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
1553	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
1554
1555
1556	! LOCAL VARIABLES
1557
1558	INTEGER itap,ind1,l,ig,iter,k
1559	INTEGER iflag_topthermals, niter
1560
1561	REAL qcth(klon)
1562	REAL qcenv(klon)
1563	REAL qctot(klon)
1564	REAL cth(klon)
1565	REAL cenv(klon)
1566	REAL cth_vol(klon)
1567	REAL cenv_vol(klon)
1568	REAL qt_env(klon), thetal_env(klon)
1569	REAL icefracenv(klon), icefracth(klon)
1570	REAL sqrtpi,sqrt2,sqrt2pi
1571	REAL alth,alenv,ath,aenv
1572	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1573	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1574	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1575	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1576	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1577	REAL zdelta,qsatbef,zcor
1578	REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
1579	REAL qlbef
1580	REAL dqsatenv(klon), dqsatth(klon)
1581	REAL erf
1582	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1583	REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1584	REAL rhodz(klon)
1585	REAL rho(klon)
1586	REAL dz(klon)
1587	REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
1588	REAL alenvl, aenvl
1589	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1590	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
1591	REAL qmax, ttarget, stmp, cout, coutref, fraci
1592	REAL maxi, mini, pas
1593
1594	! Modifty the saturation deficit PDF in thermals
1595	! in the presence of ice crystals
1596	CHARACTER (len = 80) :: abort_message
1597	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
1598
1599
1600	!------------------------------------------------------------------------------
1601	! Initialisation
1602	!------------------------------------------------------------------------------
1603
1604
1605	! Few initial checksS
1606
1607	DO k = 1,klev
1608	DO ind1 = 1, klon
1609	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
1610	zrho(ind1) = play(ind1)/temp(ind1)/rd !kg/m3
1611	dz(ind1) = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
1612	END DO
1613	END DO
1614
1615
1616	icefracth(:)=0.
1617	icefracenv(:)=0.
1618	sqrt2pi=sqrt(2.*rpi)
1619	sqrt2=sqrt(2.)
1620	sqrtpi=sqrt(rpi)
1621	icefrac(:)=0.
1622
1623	!-------------------------------------------------------------------------------
1624	! Identify grid points with potential mixed-phase conditions
1625	!-------------------------------------------------------------------------------
1626
1627
1628	DO ind1=1,klon
1629	IF ((temp(ind1) .LT. RTT) .AND. (temp(ind1) .GT. temp_nowater) &
1630	.AND. (ind2<=klev-2) &
1631	.AND. (frac_th(ind1).GT.min_frac_th_cld)) THEN
1632	mpc_bl_points(ind1,ind2)=1
1633	ELSE
1634	mpc_bl_points(ind1,ind2)=0
1635	ENDIF
1636	ENDDO
1637
1638
1639	!-------------------------------------------------------------------------------
1640	! Thermal fraction calculation and standard deviation of the distribution
1641	!-------------------------------------------------------------------------------
1642
1643	! initialisations and calculation of temperature, humidity and saturation specific humidity
1644
1645	DO ind1=1,klon
1646
1647	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg ! kg/m2
1648	rho(ind1) = play(ind1)/temp(ind1)/rd ! kg/m3
1649	dz(ind1) = rhodz(ind1)/rho(ind1) ! m : epaisseur de la couche en metre
1650	Tbefenv(ind1) = temp(ind1)
1651	thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
1652	Tbefth(ind1) = thetal_th(ind1)*zpspsk(ind1)
1653	rhoth(ind1) = play(ind1)/Tbefth(ind1)/RD
1654	qt_env(ind1) = (qt(ind1)-frac_th(ind1)qt_th(ind1))/(1.-frac_th(ind1)) !qt = aqtth + (1-a)*qtenv
1655
1656	ENDDO
1657
1658	! Calculation of saturation specific humidity
1659	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
1660	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
1661	CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
1662
1663
1664	DO ind1=1,klon
1665
1666
1667	IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
1668
1669	! unlike in the other cloudth routine,
1670	! We consider distributions of the saturation deficit WRT liquid
1671	! at positive AND negative celcius temperature
1672	! subsequently, cloud fraction corresponds to the part of the pdf corresponding
1673	! to superstauration with respect to liquid WHATEVER the temperature
1674
1675	! Environment:
1676
1677
1678	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
1679	aenv=1./(1.+(alenv*RLVTT/rcpd))
1680	senv=aenv*(qt_env(ind1)-qslenv(ind1))
1681
1682
1683	! Thermals:
1684
1685
1686	alth=(0.622RLVTTqslth(ind1))/(rdthetal_th(ind1)*2)
1687	ath=1./(1.+(alth*RLVTT/rcpd))
1688	sth=ath*(qt_th(ind1)-qslth(ind1))
1689
1690
1691	! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
1692
1693
1694	! Standard deviation of the distributions
1695
1696	sigma1s_fraca = (sigma1s_factor*0.5)(frac_th(ind1)**sigma1s_power) / &
1697	& (1-frac_th(ind1))((sth-senv)2)*0.5
1698
1699	IF (cloudth_ratqsmin>0.) THEN
1700	sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
1701	ELSE
1702	sigma1s_ratqs = ratqs(ind1)*qt(ind1)
1703	ENDIF
1704
1705	sigma1s = sigma1s_fraca + sigma1s_ratqs
1706	IF (iflag_ratqs.eq.11) then
1707	sigma1s = ratqs(ind1)qt(ind1)aenv
1708	ENDIF
1709	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((frac_th(ind1)+0.02)sigma2s_power))+0.002qt_th(ind1)
1710
1711
1712	deltasenv=aenvvert_alphasigma1s
1713	deltasth=athvert_alpha_thsigma2s
1714
1715	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1716	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1717	exp_xenv1 = exp(-1.xenv1*2)
1718	exp_xenv2 = exp(-1.xenv2*2)
1719	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1720	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1721	exp_xth1 = exp(-1.xth1*2)
1722	exp_xth2 = exp(-1.xth2*2)
1723
1724	!surface CF
1725
1726	cth(ind1)=0.5(1.-1.erf(xth1))
1727	cenv(ind1)=0.5(1.-1.erf(xenv1))
1728	ctot(ind1)=frac_th(ind1)cth(ind1)+(1.-1.frac_th(ind1))*cenv(ind1)
1729
1730
1731	!volume CF, condensed water and ice fraction
1732
1733	!environnement
1734
1735	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1736	IntJ_CF=0.5(1.-1.erf(xenv2))
1737
1738	IF (deltasenv .LT. 1.e-10) THEN
1739	qcenv(ind1)=IntJ
1740	cenv_vol(ind1)=IntJ_CF
1741	ELSE
1742	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1743	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1744	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1745	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1746	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1747	qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
1748	cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1749	IF (Tbefenv(ind1) .LT. temp_nowater) THEN
1750	! freeze all droplets in cirrus temperature regime
1751	icefracenv(ind1)=1.
1752	ENDIF
1753	ENDIF
1754
1755
1756
1757	!thermals
1758
1759	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1760	IntJ_CF=0.5(1.-1.erf(xth2))
1761
1762	IF (deltasth .LT. 1.e-10) THEN
1763	qcth(ind1)=IntJ
1764	cth_vol(ind1)=IntJ_CF
1765	ELSE
1766	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1767	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1768	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1769	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1770	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1771	qcth(ind1)=IntJ+IntI1+IntI2+IntI3
1772	cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1773	IF (Tbefth(ind1) .LT. temp_nowater) THEN
1774	! freeze all droplets in cirrus temperature regime
1775	icefracth(ind1)=1.
1776	ENDIF
1777
1778	ENDIF
1779
1780	qctot(ind1)=frac_th(ind1)qcth(ind1)+(1.-1.frac_th(ind1))*qcenv(ind1)
1781	ctot_vol(ind1)=frac_th(ind1)cth_vol(ind1)+(1.-1.frac_th(ind1))*cenv_vol(ind1)
1782
1783	IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
1784	ctot(ind1)=0.
1785	ctot_vol(ind1)=0.
1786	qcloud(ind1)=qslenv(ind1)
1787	qincloud(ind1)=0.
1788	icefrac(ind1)=0.
1789	ELSE
1790	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1791	qincloud(ind1)=qctot(ind1)/ctot(ind1)
1792	IF (qctot(ind1) .GT. 0) THEN
1793	icefrac(ind1)=(frac_th(ind1)qcth(ind1)icefracth(ind1)+(1.-1.frac_th(ind1))qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
1794	icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
1795	ENDIF
1796	ENDIF
1797
1798
1799	! ELSE ! mpc_bl_points>0
1800
1801	ELSE ! gaussian for environment only
1802
1803
1804	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
1805	aenv=1./(1.+(alenv*RLVTT/rcpd))
1806	senv=aenv*(qt_env(ind1)-qslenv(ind1))
1807	sth=0.
1808	sigma1s=ratqs(ind1)*qt_env(ind1)
1809	sigma2s=0.
1810
1811	xenv=senv/(sqrt(2.)*sigma1s)
1812	cenv(ind1)=0.5(1.-1.erf(xenv))
1813	ctot(ind1)=0.5(1.+1.erf(xenv))
1814	ctot_vol(ind1)=ctot(ind1)
1815	qctot(ind1)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1))
1816
1817
1818	IF (ctot(ind1).LT.min_neb_th) THEN
1819	ctot(ind1)=0.
1820	qcloud(ind1)=qslenv(ind1)
1821	qincloud(ind1)=0.
1822	ELSE
1823	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1824	qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
1825	ENDIF
1826
1827
1828	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
1829
1830	! Outputs used to check the PDFs
1831	cloudth_senv(ind1) = senv
1832	cloudth_sth(ind1) = sth
1833	cloudth_sigmaenv(ind1) = sigma1s
1834	cloudth_sigmath(ind1) = sigma2s
1835
1836
1837	ENDDO !loop on klon
1838
1839
1840	RETURN
1841
1842
1843	END SUBROUTINE cloudth_mpc
1844
1845
1846
1847	! ELSE ! mpc_bl_points>0
1848	!
1849	! ! Treat boundary layer mixed phase clouds
1850	!
1851	! ! thermals
1852	! !=========
1853	!
1854	! ! ice phase
1855	! !...........
1856	!
1857	! qiceth=0.
1858	! deltazlev_mpc=dz(ind1,:)
1859	! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
1860	! pres_mpc=pplay(ind1,:)
1861	! fraca_mpc=fraca(ind1,:)
1862	! snowf_mpc=snowflux(ind1,:)
1863	! iflag_topthermals=0
1864	! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
1865	! iflag_topthermals = 1
1866	! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
1867	! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
1868	! iflag_topthermals = 2
1869	! ELSE
1870	! iflag_topthermals = 0
1871	! ENDIF
1872	!
1873	! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
1874	! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
1875	!
1876	! ! qmax calculation
1877	! sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1878	! deltasth=athlvert_alpha_thsigma2s
1879	! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
1880	! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
1881	! exp_xth1 = exp(-1.xth1*2)
1882	! exp_xth2 = exp(-1.xth2*2)
1883	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1884	! IntJ_CF=0.5(1.-1.erf(xth2))
1885	! IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1886	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1887	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1888	! IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1889	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1890	! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
1891	!
1892	!
1893	! ! Liquid phase
1894	! !................
1895	! ! We account for the effect of ice crystals in thermals on sthl
1896	! ! and on the width of the distribution
1897	!
1898	! sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
1899	! + (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
1900	!
1901	! sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
1902	! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
1903	! +0.002*zqta(ind1,ind2)
1904	! deltasthc=athlvert_alpha_thsigma2sc
1905	!
1906	!
1907	! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
1908	! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
1909	! exp_xth1 = exp(-1.xth1*2)
1910	! exp_xth2 = exp(-1.xth2*2)
1911	! IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
1912	! IntJ_CF=0.5(1.-1.erf(xth2))
1913	! IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
1914	! IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1915	! IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
1916	! IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
1917	! IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
1918	! qliqth=IntJ+IntI1+IntI2+IntI3
1919	!
1920	! qlth(ind1,ind2)=MAX(0.,qliqth)
1921	!
1922	! ! Condensed water
1923	!
1924	! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
1925	!
1926	!
1927	! ! consistency with subgrid distribution
1928	!
1929	! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
1930	! fraci=qith(ind1,ind2)/qcth(ind1,ind2)
1931	! qcth(ind1,ind2)=qmax
1932	! qith(ind1,ind2)=fraci*qmax
1933	! qlth(ind1,ind2)=(1.-fraci)*qmax
1934	! ENDIF
1935	!
1936	! ! Cloud Fraction
1937	! !...............
1938	! ! calculation of qbase which is the value of the water vapor within mixed phase clouds
1939	! ! such that the total water in cloud = qbase+qliqth+qiceth
1940	! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
1941	! ! sbase and qbase calculation (note that sbase is wrt liq so negative)
1942	! ! look for an approximate solution with iteration
1943	!
1944	! ttarget=qcth(ind1,ind2)
1945	! mini= athl*(qsith(ind1,ind2)-qslth(ind1))
1946	! maxi= 0. !athl(3.sqrt(sigma2s))
1947	! niter=20
1948	! pas=(maxi-mini)/niter
1949	! stmp=mini
1950	! sbase=stmp
1951	! coutref=1.E6
1952	! DO iter=1,niter
1953	! cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
1954	! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
1955	! IF (cout .LT. coutref) THEN
1956	! sbase=stmp
1957	! coutref=cout
1958	! ELSE
1959	! stmp=stmp+pas
1960	! ENDIF
1961	! ENDDO
1962	! qbase=MAX(0., sbase/athl+qslth(ind1))
1963	!
1964	! ! surface cloud fraction in thermals
1965	! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
1966	! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
1967	!
1968	!
1969	! !volume cloud fraction in thermals
1970	! !to be checked
1971	! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
1972	! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
1973	! exp_xth1 = exp(-1.xth1*2)
1974	! exp_xth2 = exp(-1.xth2*2)
1975	!
1976	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1977	! IntJ_CF=0.5(1.-1.erf(xth2))
1978	!
1979	! IF (deltasth .LT. 1.e-10) THEN
1980	! cth_vol(ind1,ind2)=IntJ_CF
1981	! ELSE
1982	! IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1983	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1984	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1985	! IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1986	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1987	! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1988	! ENDIF
1989	! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
1990	!
1991	!
1992	!
1993	! ! Environment
1994	! !=============
1995	! ! In the environment/downdrafts, only liquid clouds
1996	! ! See Shupe et al. 2008, JAS
1997	!
1998	! ! standard deviation of the distribution in the environment
1999	! sth=sthl
2000	! senv=senvl
2001	! sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
2002	! & (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
2003	! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
2004	! ! in the environement
2005	!
2006	! sigma1s_ratqs=1E-10
2007	! IF (cloudth_ratqsmin>0.) THEN
2008	! sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2009	! ENDIF
2010	!
2011	! sigma1s = sigma1s_fraca + sigma1s_ratqs
2012	! IF (iflag_ratqs.eq.11) then
2013	! sigma1s = ratqs(ind1,ind2)po(ind1)aenv
2014	! ENDIF
2015	! IF (iflag_ratqs.eq.11) then
2016	! sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
2017	! ENDIF
2018	! deltasenv=aenvlvert_alphasigma1s
2019	! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
2020	! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
2021	! exp_xenv1 = exp(-1.xenv1*2)
2022	! exp_xenv2 = exp(-1.xenv2*2)
2023	!
2024	! !surface CF
2025	! cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
2026	!
2027	! !volume CF and condensed water
2028	! IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2029	! IntJ_CF=0.5(1.-1.erf(xenv2))
2030	!
2031	! IF (deltasenv .LT. 1.e-10) THEN
2032	! qcenv(ind1,ind2)=IntJ
2033	! cenv_vol(ind1,ind2)=IntJ_CF
2034	! ELSE
2035	! IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
2036	! IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
2037	! IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
2038	! IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
2039	! IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
2040	! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
2041	! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2042	! ENDIF
2043	!
2044	! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
2045	! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
2046	!
2047	!
2048	!
2049	! ! Thermals + environment
2050	! ! =======================
2051	! ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
2052	! qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
2053	! ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2054	! IF (qcth(ind1,ind2) .GT. 0) THEN
2055	! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
2056	! /(fraca(ind1,ind2)*qcth(ind1,ind2) &
2057	! +(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
2058	! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
2059	! ELSE
2060	! icefrac(ind1,ind2)=0.
2061	! ENDIF
2062	!
2063	! IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
2064	! ctot(ind1,ind2)=0.
2065	! ctot_vol(ind1,ind2)=0.
2066	! qincloud(ind1)=0.
2067	! qcloud(ind1)=zqsatenv(ind1)
2068	! ELSE
2069	! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
2070	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
2071	! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
2072	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
2073	! ENDIF
2074	!
2075	! ENDIF ! mpc_bl_points
2076
2077
2078
2079	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2080	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)
2081
2082	! parameterization of ice for boundary
2083	! layer mixed-phase clouds assuming a stationary system
2084	!
2085	! Note that vapor deposition on ice crystals and riming of liquid droplets
2086	! depend on the ice number concentration Ni
2087	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
2088	! and use values from Hong et al. 2004, MWR for instance
2089	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
2090	! One could also think of a more complex expression of Ni;
2091	! function of qi, T, the concentration in aerosols or INP ..
2092	! Here we prefer fixing Ni to a tuning parameter
2093	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
2094	! in Mioche et al. 2017
2095	!
2096	!
2097	! References:
2098	!------------
2099	! This parameterization is thoroughly described in Vignon et al.
2100	!
2101	! More specifically
2102	! for the Water vapor deposition process:
2103	!
2104	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
2105	! ment of stratiform cloudfs and precipitation in large-scale
2106	! models. I: Description and evaluation of the microphysical
2107	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
2108	!
2109	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
2110	! stratiform cloud microphysics scheme in the NCAR Com-
2111	! munity Atmosphere Model (CAM3). Part I: Description and
2112	! numerical tests. J. Climate, 21, 3642–3659
2113	!
2114	! for the Riming process:
2115	!
2116	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
2117	! scale structure and organization of clouds and precipitation in
2118	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
2119	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
2120	!
2121	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
2122	! forecasts of winter precipitation using an improved bulk
2123	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
2124	! forecasts of winter precipitation using an improved bulk
2125	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
2126	!
2127	! For the formation of clouds by thermals:
2128	!
2129	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
2130	! the Atmospheric Sciences, 65, 407–425.
2131	!
2132	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
2133	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
2134	!
2135	!
2136	!
2137	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
2138	!=============================================================================
2139
2140	use phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
2141
2142	IMPLICIT none
2143
2144	INCLUDE "YOMCST.h"
2145
2146	INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
2147	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
2148	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
2149	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
2150	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
2151	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
2152	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
2153	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
2154	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
2155	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
2156	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
2157	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
2158	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
2159	REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
2160	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
2161
2162
2163	INTEGER ind2p1,ind2p2
2164	REAL rho(klev)
2165	REAL unsurtaudet, unsurtaustardep, unsurtaurim
2166	REAL qi, AA, BB, Ka, Dv, rhoi
2167	REAL p0, t0, fp1, fp2
2168	REAL alpha, flux_term
2169	REAL det_term, precip_term, rim_term, dep_term
2170
2171
2172	ind2p1=ind2+1
2173	ind2p2=ind2+2
2174
2175	rho=pres/temp/RD ! air density kg/m3
2176
2177	Ka=2.4e-2 ! thermal conductivity of the air, SI
2178	p0=101325.0 ! ref pressure
2179	T0=273.15 ! ref temp
2180	rhoi=500.0 ! cloud ice density following Reisner et al. 1998
2181	alpha=700. ! fallvelocity param
2182
2183
2184	IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
2185
2186	Dv=0.00010.211(p0/pres(ind2))((temp(ind2)/T0)*1.94) ! water vapor diffusivity in air, SI
2187
2188	! Detrainment term:
2189
2190	unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
2191
2192
2193	! Deposition term
2194	AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
2195	BB=1./(rho(ind2)Dvqsith(ind2))
2196	unsurtaustardep=C_cap(Ni0.66)(qth(ind2)-qsith(ind2))/qsith(ind2) &
2197	4.RPI/(AA+BB)(6.rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
2198
2199	! Riming term neglected at this level
2200	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
2201
2202	qi=fraca(ind2)unsurtaustardep/MAX((d_topunsurtaudet),1E-12)
2203	qi=MAX(qi,0.)**(3./2.)
2204
2205	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
2206
2207	Dv=0.00010.211(p0/pres(ind2p1))((temp(ind2p1)/T0)*1.94) ! water vapor diffusivity in air, SI
2208
2209	! Detrainment term:
2210
2211	unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
2212	det_term=-unsurtaudetqith(ind2p1)rho(ind2p1)
2213
2214
2215	! Deposition term
2216
2217	AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
2218	BB=1./(rho(ind2p1)Dvqsith(ind2p1))
2219	unsurtaustardep=C_cap(Ni0.66)(qth(ind2p1)-qsith(ind2p1)) &
2220	/qsith(ind2p1)4.RPI/(AA+BB) &
2221	(6.rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
2222	dep_term=rho(ind2p1)fraca(ind2p1)(qith(ind2p1)*0.33)unsurtaustardep
2223
2224	! Riming term
2225
2226	unsurtaurim=rho(ind2p1)alpha3./rhoi/2.Eiqlth(ind2p1)((p0/pres(ind2p1))*0.4)
2227	rim_term=rho(ind2p1)fraca(ind2p1)qith(ind2p1)*unsurtaurim
2228
2229	! Precip term
2230
2231	! We assume that there is no solid precipitation outside thermals
2232	! so the precipitation flux within thermals is equal to the precipitation flux
2233	! at mesh-scale divided by thermals fraction
2234
2235	fp2=0.
2236	fp1=0.
2237	IF (fraca(ind2p1) .GT. 0.) THEN
2238	fp2=-qith(ind2p2)rho(ind2p2)vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
2239	fp1=-qith(ind2p1)rho(ind2p1)vith(ind2p1)*fraca(ind2p1)
2240	ENDIF
2241
2242	precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
2243
2244	! Calculation in a top-to-bottom loop
2245
2246	IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
2247	qi= 1./fm_therm(ind1,ind2p1)* &
2248	(deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
2249	fm_therm(ind1,ind2p2)*(qith(ind2p1)))
2250	END IF
2251
2252	ENDIF ! top thermals
2253
2254	qith(ind2)=MAX(0.,qi)
2255
2256	RETURN
2257
2258	END SUBROUTINE ICE_MPC_BL_CLOUDS
2259
2260	END MODULE lmdz_cloudth

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

Download in other formats:

Original Format