Context Navigation

lmdz_cloudth.f90 @ 5285

Last change on this file since 5285 was 5285, checked in by abarral, 10 hours ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 87.3 KB

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

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

Download in other formats:

Original Format