Context Navigation

cloudth_mod.F90 @ 4114

Last change on this file since 4114 was 4114, checked in by evignon, 2 years ago
option pour inclure les effets radiatifs des chutes de neige plus corrections du schema des nuages de phase mixte
File size: 92.3 KB

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

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

Download in other formats:

Original Format