Context Navigation

lmdz_cloudth.F90 @ 5097

Last change on this file since 5097 was 5087, checked in by abarral, 4 months ago
remove fixed-form ^{\s+& remaining in .f90,.F90}
File size: 86.8 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>=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)>ztv(ind1,2)).and.(fraca(ind1,ind2)>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)<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)<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)>ztv(ind1,2)).and.(fraca(ind1,ind2)>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)<1.e-10.or.cth(ind1,ind2)<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)<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>=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)>ztv(ind1,2)).and.(fraca(ind1,ind2)>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)<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)<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)>ztv(ind1,2)).and.(fraca(ind1,ind2)>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==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 < 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 < 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)<1.e-10.or.cth(ind1,ind2)<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)<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)>ztv(ind1,2)).and.(fraca(ind1,ind2)>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) < 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) < 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,mpc_bl_points, &
1500	& temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th, &
1501	& ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol, &
1502	& cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
1503	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1504	! Author : Etienne Vignon (LMDZ/CNRS)
1505	! Date: April 2024
1506	! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
1507	! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
1508	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1509
1510	use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
1511	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
1512	use lmdz_lscp_tools, only: CALC_QSAT_ECMWF
1513	use lmdz_lscp_ini, only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
1514
1515	IMPLICIT NONE
1516
1517
1518	!------------------------------------------------------------------------------
1519	! Declaration
1520	!------------------------------------------------------------------------------
1521
1522	! INPUT/OUTPUT
1523
1524	INTEGER, INTENT(IN) :: klon,klev,ind2
1525
1526
1527	REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
1528	REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
1529	REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
1530	REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
1531	REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1]
1532	REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential
1533	REAL, DIMENSION(klon), INTENT(IN) :: paprsup ! Pressure at top layer interface [Pa]
1534	REAL, DIMENSION(klon), INTENT(IN) :: paprsdn ! Pressure at bottom layer interface [Pa]
1535	REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa]
1536	REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
1537
1538
1539	INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
1540
1541	REAL, DIMENSION(klon), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
1542	REAL, DIMENSION(klon), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
1543	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
1544	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
1545	REAL, DIMENSION(klon), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
1546	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
1547	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
1548	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
1549	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
1550
1551
1552	! LOCAL VARIABLES
1553
1554	INTEGER itap,ind1,l,ig,iter,k
1555	INTEGER iflag_topthermals, niter
1556
1557	REAL qcth(klon)
1558	REAL qcenv(klon)
1559	REAL qctot(klon)
1560	REAL cth(klon)
1561	REAL cenv(klon)
1562	REAL cth_vol(klon)
1563	REAL cenv_vol(klon)
1564	REAL qt_env(klon), thetal_env(klon)
1565	REAL icefracenv(klon), icefracth(klon)
1566	REAL sqrtpi,sqrt2,sqrt2pi
1567	REAL alth,alenv,ath,aenv
1568	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1569	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1570	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1571	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1572	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1573	REAL zdelta,qsatbef,zcor
1574	REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
1575	REAL qlbef
1576	REAL dqsatenv(klon), dqsatth(klon)
1577	REAL erf
1578	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1579	REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1580	REAL rhodz(klon)
1581	REAL rho(klon)
1582	REAL dz(klon)
1583	REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
1584	REAL alenvl, aenvl
1585	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1586	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
1587	REAL qmax, ttarget, stmp, cout, coutref, fraci
1588	REAL maxi, mini, pas
1589
1590	! Modifty the saturation deficit PDF in thermals
1591	! in the presence of ice crystals
1592	CHARACTER (len = 80) :: abort_message
1593	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
1594
1595
1596	!------------------------------------------------------------------------------
1597	! Initialisation
1598	!------------------------------------------------------------------------------
1599
1600
1601	! Few initial checksS
1602
1603	DO k = 1,klev
1604	DO ind1 = 1, klon
1605	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
1606	zrho(ind1) = play(ind1)/temp(ind1)/rd !kg/m3
1607	dz(ind1) = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
1608	END DO
1609	END DO
1610
1611
1612	icefracth(:)=0.
1613	icefracenv(:)=0.
1614	sqrt2pi=sqrt(2.*rpi)
1615	sqrt2=sqrt(2.)
1616	sqrtpi=sqrt(rpi)
1617	icefrac(:)=0.
1618
1619	!-------------------------------------------------------------------------------
1620	! Identify grid points with potential mixed-phase conditions
1621	!-------------------------------------------------------------------------------
1622
1623
1624	DO ind1=1,klon
1625	IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) &
1626	.AND. (ind2<=klev-2) &
1627	.AND. (frac_th(ind1)>min_frac_th_cld)) THEN
1628	mpc_bl_points(ind1,ind2)=1
1629	ELSE
1630	mpc_bl_points(ind1,ind2)=0
1631	ENDIF
1632	ENDDO
1633
1634
1635	!-------------------------------------------------------------------------------
1636	! Thermal fraction calculation and standard deviation of the distribution
1637	!-------------------------------------------------------------------------------
1638
1639	! initialisations and calculation of temperature, humidity and saturation specific humidity
1640
1641	DO ind1=1,klon
1642
1643	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg ! kg/m2
1644	rho(ind1) = play(ind1)/temp(ind1)/rd ! kg/m3
1645	dz(ind1) = rhodz(ind1)/rho(ind1) ! m : epaisseur de la couche en metre
1646	Tbefenv(ind1) = temp(ind1)
1647	thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
1648	Tbefth(ind1) = thetal_th(ind1)*zpspsk(ind1)
1649	rhoth(ind1) = play(ind1)/Tbefth(ind1)/RD
1650	qt_env(ind1) = (qt(ind1)-frac_th(ind1)qt_th(ind1))/(1.-frac_th(ind1)) !qt = aqtth + (1-a)*qtenv
1651
1652	ENDDO
1653
1654	! Calculation of saturation specific humidity
1655	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
1656	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
1657	CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
1658
1659
1660	DO ind1=1,klon
1661
1662
1663	IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement
1664
1665	! unlike in the other cloudth routine,
1666	! We consider distributions of the saturation deficit WRT liquid
1667	! at positive AND negative celcius temperature
1668	! subsequently, cloud fraction corresponds to the part of the pdf corresponding
1669	! to superstauration with respect to liquid WHATEVER the temperature
1670
1671	! Environment:
1672
1673
1674	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
1675	aenv=1./(1.+(alenv*RLVTT/rcpd))
1676	senv=aenv*(qt_env(ind1)-qslenv(ind1))
1677
1678
1679	! Thermals:
1680
1681
1682	alth=(0.622RLVTTqslth(ind1))/(rdthetal_th(ind1)*2)
1683	ath=1./(1.+(alth*RLVTT/rcpd))
1684	sth=ath*(qt_th(ind1)-qslth(ind1))
1685
1686
1687	! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
1688
1689
1690	! Standard deviation of the distributions
1691
1692	sigma1s_fraca = (sigma1s_factor*0.5)(frac_th(ind1)**sigma1s_power) / &
1693	(1-frac_th(ind1))((sth-senv)2)*0.5
1694
1695	IF (cloudth_ratqsmin>0.) THEN
1696	sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
1697	ELSE
1698	sigma1s_ratqs = ratqs(ind1)*qt(ind1)
1699	ENDIF
1700
1701	sigma1s = sigma1s_fraca + sigma1s_ratqs
1702	IF (iflag_ratqs==11) then
1703	sigma1s = ratqs(ind1)qt(ind1)aenv
1704	ENDIF
1705	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((frac_th(ind1)+0.02)sigma2s_power))+0.002qt_th(ind1)
1706
1707
1708	deltasenv=aenvvert_alphasigma1s
1709	deltasth=athvert_alpha_thsigma2s
1710
1711	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1712	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1713	exp_xenv1 = exp(-1.xenv1*2)
1714	exp_xenv2 = exp(-1.xenv2*2)
1715	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1716	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1717	exp_xth1 = exp(-1.xth1*2)
1718	exp_xth2 = exp(-1.xth2*2)
1719
1720	!surface CF
1721
1722	cth(ind1)=0.5(1.-1.erf(xth1))
1723	cenv(ind1)=0.5(1.-1.erf(xenv1))
1724	ctot(ind1)=frac_th(ind1)cth(ind1)+(1.-1.frac_th(ind1))*cenv(ind1)
1725
1726
1727	!volume CF, condensed water and ice fraction
1728
1729	!environnement
1730
1731	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1732	IntJ_CF=0.5(1.-1.erf(xenv2))
1733
1734	IF (deltasenv < 1.e-10) THEN
1735	qcenv(ind1)=IntJ
1736	cenv_vol(ind1)=IntJ_CF
1737	ELSE
1738	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1739	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1740	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1741	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1742	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1743	qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
1744	cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1745	IF (Tbefenv(ind1) < temp_nowater) THEN
1746	! freeze all droplets in cirrus temperature regime
1747	icefracenv(ind1)=1.
1748	ENDIF
1749	ENDIF
1750
1751
1752
1753	!thermals
1754
1755	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1756	IntJ_CF=0.5(1.-1.erf(xth2))
1757
1758	IF (deltasth < 1.e-10) THEN
1759	qcth(ind1)=IntJ
1760	cth_vol(ind1)=IntJ_CF
1761	ELSE
1762	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1763	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1764	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1765	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1766	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1767	qcth(ind1)=IntJ+IntI1+IntI2+IntI3
1768	cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
1769	IF (Tbefth(ind1) < temp_nowater) THEN
1770	! freeze all droplets in cirrus temperature regime
1771	icefracth(ind1)=1.
1772	ENDIF
1773
1774	ENDIF
1775
1776	qctot(ind1)=frac_th(ind1)qcth(ind1)+(1.-1.frac_th(ind1))*qcenv(ind1)
1777	ctot_vol(ind1)=frac_th(ind1)cth_vol(ind1)+(1.-1.frac_th(ind1))*cenv_vol(ind1)
1778
1779	IF (cenv(ind1)<min_neb_th.and.cth(ind1)<min_neb_th) THEN
1780	ctot(ind1)=0.
1781	ctot_vol(ind1)=0.
1782	qcloud(ind1)=qslenv(ind1)
1783	qincloud(ind1)=0.
1784	icefrac(ind1)=0.
1785	ELSE
1786	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1787	qincloud(ind1)=qctot(ind1)/ctot(ind1)
1788	IF (qctot(ind1) > 0) THEN
1789	icefrac(ind1)=(frac_th(ind1)qcth(ind1)icefracth(ind1)+(1.-1.frac_th(ind1))qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
1790	icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
1791	ENDIF
1792	ENDIF
1793
1794
1795	! ELSE ! mpc_bl_points>0
1796
1797	ELSE ! gaussian for environment only
1798
1799
1800	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
1801	aenv=1./(1.+(alenv*RLVTT/rcpd))
1802	senv=aenv*(qt_env(ind1)-qslenv(ind1))
1803	sth=0.
1804	sigma1s=ratqs(ind1)*qt_env(ind1)
1805	sigma2s=0.
1806
1807	xenv=senv/(sqrt(2.)*sigma1s)
1808	cenv(ind1)=0.5(1.-1.erf(xenv))
1809	ctot(ind1)=0.5(1.+1.erf(xenv))
1810	ctot_vol(ind1)=ctot(ind1)
1811	qctot(ind1)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1))
1812
1813
1814	IF (ctot(ind1)<min_neb_th) THEN
1815	ctot(ind1)=0.
1816	qcloud(ind1)=qslenv(ind1)
1817	qincloud(ind1)=0.
1818	ELSE
1819	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
1820	qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
1821	ENDIF
1822
1823
1824	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
1825
1826	! Outputs used to check the PDFs
1827	cloudth_senv(ind1) = senv
1828	cloudth_sth(ind1) = sth
1829	cloudth_sigmaenv(ind1) = sigma1s
1830	cloudth_sigmath(ind1) = sigma2s
1831
1832
1833	ENDDO !loop on klon
1834
1835
1836	RETURN
1837
1838
1839	END SUBROUTINE cloudth_mpc
1840
1841
1842
1843	! ELSE ! mpc_bl_points>0
1844	!
1845	! ! Treat boundary layer mixed phase clouds
1846	!
1847	! ! thermals
1848	! !=========
1849	!
1850	! ! ice phase
1851	! !...........
1852	!
1853	! qiceth=0.
1854	! deltazlev_mpc=dz(ind1,:)
1855	! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
1856	! pres_mpc=pplay(ind1,:)
1857	! fraca_mpc=fraca(ind1,:)
1858	! snowf_mpc=snowflux(ind1,:)
1859	! iflag_topthermals=0
1860	! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
1861	! iflag_topthermals = 1
1862	! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
1863	! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
1864	! iflag_topthermals = 2
1865	! ELSE
1866	! iflag_topthermals = 0
1867	! ENDIF
1868	!
1869	! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
1870	! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
1871	!
1872	! ! qmax calculation
1873	! sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1874	! deltasth=athlvert_alpha_thsigma2s
1875	! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
1876	! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
1877	! exp_xth1 = exp(-1.xth1*2)
1878	! exp_xth2 = exp(-1.xth2*2)
1879	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1880	! IntJ_CF=0.5(1.-1.erf(xth2))
1881	! IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1882	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1883	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1884	! IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1885	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1886	! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
1887	!
1888	!
1889	! ! Liquid phase
1890	! !................
1891	! ! We account for the effect of ice crystals in thermals on sthl
1892	! ! and on the width of the distribution
1893	!
1894	! sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
1895	! + (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
1896	!
1897	! sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
1898	! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
1899	! +0.002*zqta(ind1,ind2)
1900	! deltasthc=athlvert_alpha_thsigma2sc
1901	!
1902	!
1903	! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
1904	! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
1905	! exp_xth1 = exp(-1.xth1*2)
1906	! exp_xth2 = exp(-1.xth2*2)
1907	! IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
1908	! IntJ_CF=0.5(1.-1.erf(xth2))
1909	! IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
1910	! IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1911	! IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
1912	! IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
1913	! IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
1914	! qliqth=IntJ+IntI1+IntI2+IntI3
1915	!
1916	! qlth(ind1,ind2)=MAX(0.,qliqth)
1917	!
1918	! ! Condensed water
1919	!
1920	! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
1921	!
1922	!
1923	! ! consistency with subgrid distribution
1924	!
1925	! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
1926	! fraci=qith(ind1,ind2)/qcth(ind1,ind2)
1927	! qcth(ind1,ind2)=qmax
1928	! qith(ind1,ind2)=fraci*qmax
1929	! qlth(ind1,ind2)=(1.-fraci)*qmax
1930	! ENDIF
1931	!
1932	! ! Cloud Fraction
1933	! !...............
1934	! ! calculation of qbase which is the value of the water vapor within mixed phase clouds
1935	! ! such that the total water in cloud = qbase+qliqth+qiceth
1936	! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
1937	! ! sbase and qbase calculation (note that sbase is wrt liq so negative)
1938	! ! look for an approximate solution with iteration
1939	!
1940	! ttarget=qcth(ind1,ind2)
1941	! mini= athl*(qsith(ind1,ind2)-qslth(ind1))
1942	! maxi= 0. !athl(3.sqrt(sigma2s))
1943	! niter=20
1944	! pas=(maxi-mini)/niter
1945	! stmp=mini
1946	! sbase=stmp
1947	! coutref=1.E6
1948	! DO iter=1,niter
1949	! cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
1950	! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
1951	! IF (cout .LT. coutref) THEN
1952	! sbase=stmp
1953	! coutref=cout
1954	! ELSE
1955	! stmp=stmp+pas
1956	! ENDIF
1957	! ENDDO
1958	! qbase=MAX(0., sbase/athl+qslth(ind1))
1959	!
1960	! ! surface cloud fraction in thermals
1961	! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
1962	! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
1963	!
1964	!
1965	! !volume cloud fraction in thermals
1966	! !to be checked
1967	! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
1968	! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
1969	! exp_xth1 = exp(-1.xth1*2)
1970	! exp_xth2 = exp(-1.xth2*2)
1971	!
1972	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1973	! IntJ_CF=0.5(1.-1.erf(xth2))
1974	!
1975	! IF (deltasth .LT. 1.e-10) THEN
1976	! cth_vol(ind1,ind2)=IntJ_CF
1977	! ELSE
1978	! IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1979	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1980	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1981	! IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1982	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1983	! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1984	! ENDIF
1985	! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
1986	!
1987	!
1988	!
1989	! ! Environment
1990	! !=============
1991	! ! In the environment/downdrafts, only liquid clouds
1992	! ! See Shupe et al. 2008, JAS
1993	!
1994	! ! standard deviation of the distribution in the environment
1995	! sth=sthl
1996	! senv=senvl
1997	! sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
1998	! & (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
1999	! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
2000	! ! in the environement
2001	!
2002	! sigma1s_ratqs=1E-10
2003	! IF (cloudth_ratqsmin>0.) THEN
2004	! sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2005	! ENDIF
2006	!
2007	! sigma1s = sigma1s_fraca + sigma1s_ratqs
2008	! IF (iflag_ratqs.eq.11) then
2009	! sigma1s = ratqs(ind1,ind2)po(ind1)aenv
2010	! ENDIF
2011	! IF (iflag_ratqs.eq.11) then
2012	! sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
2013	! ENDIF
2014	! deltasenv=aenvlvert_alphasigma1s
2015	! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
2016	! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
2017	! exp_xenv1 = exp(-1.xenv1*2)
2018	! exp_xenv2 = exp(-1.xenv2*2)
2019	!
2020	! !surface CF
2021	! cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
2022	!
2023	! !volume CF and condensed water
2024	! IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2025	! IntJ_CF=0.5(1.-1.erf(xenv2))
2026	!
2027	! IF (deltasenv .LT. 1.e-10) THEN
2028	! qcenv(ind1,ind2)=IntJ
2029	! cenv_vol(ind1,ind2)=IntJ_CF
2030	! ELSE
2031	! IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
2032	! IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
2033	! IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
2034	! IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
2035	! IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
2036	! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
2037	! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2038	! ENDIF
2039	!
2040	! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
2041	! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
2042	!
2043	!
2044	!
2045	! ! Thermals + environment
2046	! ! =======================
2047	! ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
2048	! qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
2049	! ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2050	! IF (qcth(ind1,ind2) .GT. 0) THEN
2051	! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
2052	! /(fraca(ind1,ind2)*qcth(ind1,ind2) &
2053	! +(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
2054	! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
2055	! ELSE
2056	! icefrac(ind1,ind2)=0.
2057	! ENDIF
2058	!
2059	! IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
2060	! ctot(ind1,ind2)=0.
2061	! ctot_vol(ind1,ind2)=0.
2062	! qincloud(ind1)=0.
2063	! qcloud(ind1)=zqsatenv(ind1)
2064	! ELSE
2065	! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
2066	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
2067	! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
2068	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
2069	! ENDIF
2070	!
2071	! ENDIF ! mpc_bl_points
2072
2073
2074
2075	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2076	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)
2077
2078	! parameterization of ice for boundary
2079	! layer mixed-phase clouds assuming a stationary system
2080	!
2081	! Note that vapor deposition on ice crystals and riming of liquid droplets
2082	! depend on the ice number concentration Ni
2083	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
2084	! and use values from Hong et al. 2004, MWR for instance
2085	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
2086	! One could also think of a more complex expression of Ni;
2087	! function of qi, T, the concentration in aerosols or INP ..
2088	! Here we prefer fixing Ni to a tuning parameter
2089	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
2090	! in Mioche et al. 2017
2091	!
2092	!
2093	! References:
2094	!------------
2095	! This parameterization is thoroughly described in Vignon et al.
2096	!
2097	! More specifically
2098	! for the Water vapor deposition process:
2099	!
2100	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
2101	! ment of stratiform cloudfs and precipitation in large-scale
2102	! models. I: Description and evaluation of the microphysical
2103	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
2104	!
2105	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
2106	! stratiform cloud microphysics scheme in the NCAR Com-
2107	! munity Atmosphere Model (CAM3). Part I: Description and
2108	! numerical tests. J. Climate, 21, 3642–3659
2109	!
2110	! for the Riming process:
2111	!
2112	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
2113	! scale structure and organization of clouds and precipitation in
2114	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
2115	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
2116	!
2117	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
2118	! forecasts of winter precipitation using an improved bulk
2119	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
2120	! forecasts of winter precipitation using an improved bulk
2121	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
2122	!
2123	! For the formation of clouds by thermals:
2124	!
2125	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
2126	! the Atmospheric Sciences, 65, 407–425.
2127	!
2128	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
2129	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
2130	!
2131	!
2132	!
2133	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
2134	!=============================================================================
2135
2136	use phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
2137
2138	IMPLICIT none
2139
2140	INCLUDE "YOMCST.h"
2141
2142	INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
2143	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
2144	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
2145	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
2146	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
2147	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
2148	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
2149	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
2150	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
2151	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
2152	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
2153	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
2154	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
2155	REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
2156	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
2157
2158
2159	INTEGER ind2p1,ind2p2
2160	REAL rho(klev)
2161	REAL unsurtaudet, unsurtaustardep, unsurtaurim
2162	REAL qi, AA, BB, Ka, Dv, rhoi
2163	REAL p0, t0, fp1, fp2
2164	REAL alpha, flux_term
2165	REAL det_term, precip_term, rim_term, dep_term
2166
2167
2168	ind2p1=ind2+1
2169	ind2p2=ind2+2
2170
2171	rho=pres/temp/RD ! air density kg/m3
2172
2173	Ka=2.4e-2 ! thermal conductivity of the air, SI
2174	p0=101325.0 ! ref pressure
2175	T0=273.15 ! ref temp
2176	rhoi=500.0 ! cloud ice density following Reisner et al. 1998
2177	alpha=700. ! fallvelocity param
2178
2179
2180	IF (iflag_topthermals > 0) THEN ! uppermost thermals levels
2181
2182	Dv=0.00010.211(p0/pres(ind2))((temp(ind2)/T0)*1.94) ! water vapor diffusivity in air, SI
2183
2184	! Detrainment term:
2185
2186	unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
2187
2188
2189	! Deposition term
2190	AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
2191	BB=1./(rho(ind2)Dvqsith(ind2))
2192	unsurtaustardep=C_cap(Ni0.66)(qth(ind2)-qsith(ind2))/qsith(ind2) &
2193	4.RPI/(AA+BB)(6.rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
2194
2195	! Riming term neglected at this level
2196	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
2197
2198	qi=fraca(ind2)unsurtaustardep/MAX((d_topunsurtaudet),1E-12)
2199	qi=MAX(qi,0.)**(3./2.)
2200
2201	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
2202
2203	Dv=0.00010.211(p0/pres(ind2p1))((temp(ind2p1)/T0)*1.94) ! water vapor diffusivity in air, SI
2204
2205	! Detrainment term:
2206
2207	unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
2208	det_term=-unsurtaudetqith(ind2p1)rho(ind2p1)
2209
2210
2211	! Deposition term
2212
2213	AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
2214	BB=1./(rho(ind2p1)Dvqsith(ind2p1))
2215	unsurtaustardep=C_cap(Ni0.66)(qth(ind2p1)-qsith(ind2p1)) &
2216	/qsith(ind2p1)4.RPI/(AA+BB) &
2217	(6.rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
2218	dep_term=rho(ind2p1)fraca(ind2p1)(qith(ind2p1)*0.33)unsurtaustardep
2219
2220	! Riming term
2221
2222	unsurtaurim=rho(ind2p1)alpha3./rhoi/2.Eiqlth(ind2p1)((p0/pres(ind2p1))*0.4)
2223	rim_term=rho(ind2p1)fraca(ind2p1)qith(ind2p1)*unsurtaurim
2224
2225	! Precip term
2226
2227	! We assume that there is no solid precipitation outside thermals
2228	! so the precipitation flux within thermals is equal to the precipitation flux
2229	! at mesh-scale divided by thermals fraction
2230
2231	fp2=0.
2232	fp1=0.
2233	IF (fraca(ind2p1) > 0.) THEN
2234	fp2=-qith(ind2p2)rho(ind2p2)vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
2235	fp1=-qith(ind2p1)rho(ind2p1)vith(ind2p1)*fraca(ind2p1)
2236	ENDIF
2237
2238	precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
2239
2240	! Calculation in a top-to-bottom loop
2241
2242	IF (fm_therm(ind1,ind2p1) > 0.) THEN
2243	qi= 1./fm_therm(ind1,ind2p1)* &
2244	(deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
2245	fm_therm(ind1,ind2p2)*(qith(ind2p1)))
2246	END IF
2247
2248	ENDIF ! top thermals
2249
2250	qith(ind2)=MAX(0.,qi)
2251
2252	RETURN
2253
2254	END SUBROUTINE ICE_MPC_BL_CLOUDS
2255
2256	END MODULE lmdz_cloudth

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

Download in other formats:

Original Format