Context Navigation

lmdz_cloudth.f90 @ 5479

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