Context Navigation

lmdz_cloudth.F90 @ 5143

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

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

Download in other formats:

Original Format