Context Navigation

cloudth_mod.F90 @ 4380

Last change on this file since 4380 was 4380, checked in by evignon, 17 months ago
premier commit d'un travail en cours sur l'externalisation de la routine lscp pour l'utilisation du replay + nettoyage
File size: 92.6 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) &
1149	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1150	+ratqs(ind1,ind2)*po(ind1) !Environment
1151	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
1152	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1153	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1154	xth=sth/(sqrt(2.)*sigma2s)
1155	xenv=senv/(sqrt(2.)*sigma1s)
1156
1157	!Volumique
1158	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1159	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1160	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1161	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1162
1163	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1164	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1165	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1166
1167	!Surfacique
1168	!Neggers
1169	!beta=0.0044
1170	!inverse_rho=1.+beta*dz(ind1,ind2)
1171	!print *,'jeanjean : beta=',beta
1172	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1173	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1174	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1175
1176	!Brooks
1177	a_Brooks=0.6694
1178	b_Brooks=0.1882
1179	A_Maj_Brooks=0.1635 !-- sans shear
1180	!A_Maj_Brooks=0.17 !-- ARM LES
1181	!A_Maj_Brooks=0.18 !-- RICO LES
1182	!A_Maj_Brooks=0.19 !-- BOMEX LES
1183	Dx_Brooks=200000.
1184	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1185	!print *,'jeanjean_f=',f_Brooks
1186
1187	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1188	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1189	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1190	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1191
1192
1193
1194
1195
1196	ENDIF ! of if (iflag_cloudth_vert<5)
1197	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1198
1199	! if (ctot(ind1,ind2).lt.1.e-10) then
1200	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1201	ctot(ind1,ind2)=0.
1202	ctot_vol(ind1,ind2)=0.
1203	qcloud(ind1)=zqsatenv(ind1,ind2)
1204
1205	else
1206
1207	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1208	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1209	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1210
1211	endif
1212
1213	else ! gaussienne environnement seule
1214
1215
1216	zqenv(ind1)=po(ind1)
1217	Tbef=t(ind1,ind2)
1218	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1219	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1220	qsatbef=MIN(0.5,qsatbef)
1221	zcor=1./(1.-retv*qsatbef)
1222	qsatbef=qsatbef*zcor
1223	zqsatenv(ind1,ind2)=qsatbef
1224
1225
1226	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1227	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1228	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1229	aenv=1./(1.+(alenv*Lv/cppd))
1230	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1231	sth=0.
1232
1233
1234	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1235	sigma2s=0.
1236
1237	sqrt2pi=sqrt(2.*pi)
1238	xenv=senv/(sqrt(2.)*sigma1s)
1239	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1240	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1241	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1242
1243	if (ctot(ind1,ind2).lt.1.e-3) then
1244	ctot(ind1,ind2)=0.
1245	qcloud(ind1)=zqsatenv(ind1,ind2)
1246
1247	else
1248
1249	! ctot(ind1,ind2)=ctot(ind1,ind2)
1250	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1251
1252	endif
1253
1254
1255
1256
1257	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1258	! Outputs used to check the PDFs
1259	cloudth_senv(ind1,ind2) = senv
1260	cloudth_sth(ind1,ind2) = sth
1261	cloudth_sigmaenv(ind1,ind2) = sigma1s
1262	cloudth_sigmath(ind1,ind2) = sigma2s
1263
1264	enddo ! from the loop on ngrid l.333
1265	return
1266	! end
1267	END SUBROUTINE cloudth_vert_v3
1268	!
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1281	& ztv,po,zqta,fraca, &
1282	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1283	& ratqs,zqs,T)
1284
1285
1286	USE ioipsl_getin_p_mod, ONLY : getin_p
1287	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1288	& cloudth_sigmath,cloudth_sigmaenv
1289
1290	IMPLICIT NONE
1291
1292	#include "YOMCST.h"
1293	#include "YOETHF.h"
1294	#include "FCTTRE.h"
1295	#include "nuage.h"
1296
1297
1298	!Domain variables
1299	INTEGER ngrid !indice Max lat-lon
1300	INTEGER klev !indice Max alt
1301	INTEGER ind1 !indice in [1:ngrid]
1302	INTEGER ind2 !indice in [1:klev]
1303	!thermal plume fraction
1304	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1305	!temperatures
1306	REAL T(ngrid,klev) !temperature
1307	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1308	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1309	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1310	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1311	!pressure
1312	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1313	REAL pplay(ngrid,klev) !pressure at the middle of the level
1314	!humidity
1315	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1316	REAL po(ngrid) !total water (qt)
1317	REAL zqenv(ngrid) !total water in the environment (qt_env)
1318	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1319	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1320	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1321	REAL qlth(ngrid,klev) !condensed water in the thermals
1322	REAL qlenv(ngrid,klev) !condensed water in the environment
1323	REAL qltot(ngrid,klev) !condensed water in the gridbox
1324	!cloud fractions
1325	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1326	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1327	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1328	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1329	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1330	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1331	!PDF of saturation deficit variables
1332	REAL rdd,cppd,Lv
1333	REAL Tbef,zdelta,qsatbef,zcor
1334	REAL alth,alenv,ath,aenv
1335	REAL sth,senv !saturation deficits in the thermals and environment
1336	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1337	!cloud fraction variables
1338	REAL xth,xenv
1339	REAL inverse_rho,beta !Neggers et al. (2011) method
1340	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1341	!Incloud total water variables
1342	REAL zqs(ngrid) !q_sat
1343	REAL qcloud(ngrid) !eau totale dans le nuage
1344	!Some arithmetic variables
1345	REAL erf,pi,sqrt2,sqrt2pi
1346	!Depth of the layer
1347	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1348	REAL rhodz(ngrid,klev)
1349	REAL zrho(ngrid,klev)
1350	DO ind1 = 1, ngrid
1351	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1352	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1353	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1354	END DO
1355
1356	!------------------------------------------------------------------------------
1357	! Initialization
1358	!------------------------------------------------------------------------------
1359	qlth(:,:)=0.
1360	qlenv(:,:)=0.
1361	qltot(:,:)=0.
1362	cth_vol(:,:)=0.
1363	cenv_vol(:,:)=0.
1364	ctot_vol(:,:)=0.
1365	cth_surf(:,:)=0.
1366	cenv_surf(:,:)=0.
1367	ctot_surf(:,:)=0.
1368	qcloud(:)=0.
1369	rdd=287.04
1370	cppd=1005.7
1371	pi=3.14159
1372	Lv=2.5e6
1373	sqrt2=sqrt(2.)
1374	sqrt2pi=sqrt(2.*pi)
1375
1376
1377	DO ind1=1,ngrid
1378	!-------------------------------------------------------------------------------
1379	!Both thermal and environment in the gridbox
1380	!-------------------------------------------------------------------------------
1381	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1382	!--------------------------------------------
1383	!calcul de qsat_env
1384	!--------------------------------------------
1385	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1386	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1387	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1388	qsatbef=MIN(0.5,qsatbef)
1389	zcor=1./(1.-retv*qsatbef)
1390	qsatbef=qsatbef*zcor
1391	zqsatenv(ind1,ind2)=qsatbef
1392	!--------------------------------------------
1393	!calcul de s_env
1394	!--------------------------------------------
1395	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1396	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1397	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1398	!--------------------------------------------
1399	!calcul de qsat_th
1400	!--------------------------------------------
1401	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1402	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1403	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1404	qsatbef=MIN(0.5,qsatbef)
1405	zcor=1./(1.-retv*qsatbef)
1406	qsatbef=qsatbef*zcor
1407	zqsatth(ind1,ind2)=qsatbef
1408	!--------------------------------------------
1409	!calcul de s_th
1410	!--------------------------------------------
1411	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1412	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1413	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1414	!--------------------------------------------
1415	!calcul standard deviations bi-Gaussian PDF
1416	!--------------------------------------------
1417	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)
1418	sigma_env=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
1419	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
1420	+ratqs(ind1,ind2)*po(ind1)
1421	xth=sth/(sqrt2*sigma_th)
1422	xenv=senv/(sqrt2*sigma_env)
1423	!--------------------------------------------
1424	!Cloud fraction by volume CF_vol
1425	!--------------------------------------------
1426	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1427	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1428	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1429	!--------------------------------------------
1430	!Condensed water qc
1431	!--------------------------------------------
1432	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1433	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1434	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1435	!--------------------------------------------
1436	!Cloud fraction by surface CF_surf
1437	!--------------------------------------------
1438	!Method Neggers et al. (2011) : ok for cumulus clouds only
1439	!beta=0.0044 (Jouhaud et al.2018)
1440	!inverse_rho=1.+beta*dz(ind1,ind2)
1441	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1442	!Method Brooks et al. (2005) : ok for all types of clouds
1443	a_Brooks=0.6694
1444	b_Brooks=0.1882
1445	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1446	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1447	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1448	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1449	!--------------------------------------------
1450	!Incloud Condensed water qcloud
1451	!--------------------------------------------
1452	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1453	ctot_vol(ind1,ind2)=0.
1454	ctot_surf(ind1,ind2)=0.
1455	qcloud(ind1)=zqsatenv(ind1,ind2)
1456	else
1457	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1458	endif
1459
1460
1461
1462	!-------------------------------------------------------------------------------
1463	!Environment only in the gridbox
1464	!-------------------------------------------------------------------------------
1465	ELSE
1466	!--------------------------------------------
1467	!calcul de qsat_env
1468	!--------------------------------------------
1469	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1470	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1471	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1472	qsatbef=MIN(0.5,qsatbef)
1473	zcor=1./(1.-retv*qsatbef)
1474	qsatbef=qsatbef*zcor
1475	zqsatenv(ind1,ind2)=qsatbef
1476	!--------------------------------------------
1477	!calcul de s_env
1478	!--------------------------------------------
1479	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1480	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1481	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1482	!--------------------------------------------
1483	!calcul standard deviations Gaussian PDF
1484	!--------------------------------------------
1485	zqenv(ind1)=po(ind1)
1486	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1487	xenv=senv/(sqrt2*sigma_env)
1488	!--------------------------------------------
1489	!Cloud fraction by volume CF_vol
1490	!--------------------------------------------
1491	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1492	!--------------------------------------------
1493	!Condensed water qc
1494	!--------------------------------------------
1495	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1496	!--------------------------------------------
1497	!Cloud fraction by surface CF_surf
1498	!--------------------------------------------
1499	!Method Neggers et al. (2011) : ok for cumulus clouds only
1500	!beta=0.0044 (Jouhaud et al.2018)
1501	!inverse_rho=1.+beta*dz(ind1,ind2)
1502	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1503	!Method Brooks et al. (2005) : ok for all types of clouds
1504	a_Brooks=0.6694
1505	b_Brooks=0.1882
1506	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1507	Dx_Brooks=200000.
1508	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1509	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1510	!--------------------------------------------
1511	!Incloud Condensed water qcloud
1512	!--------------------------------------------
1513	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1514	ctot_vol(ind1,ind2)=0.
1515	ctot_surf(ind1,ind2)=0.
1516	qcloud(ind1)=zqsatenv(ind1,ind2)
1517	else
1518	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1519	endif
1520
1521
1522	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1523
1524	! Outputs used to check the PDFs
1525	cloudth_senv(ind1,ind2) = senv
1526	cloudth_sth(ind1,ind2) = sth
1527	cloudth_sigmaenv(ind1,ind2) = sigma_env
1528	cloudth_sigmath(ind1,ind2) = sigma_th
1529
1530	END DO ! From the loop on ngrid
1531	return
1532
1533	END SUBROUTINE cloudth_v6
1534
1535
1536
1537
1538
1539	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1540	SUBROUTINE cloudth_mpc(klon,klev,ind2,iflag_mpc_bl,mpc_bl_points, &
1541	& temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, &
1542	& ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol)
1543	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1544	! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS)
1545	! Date: Adapted from cloudth_vert_v3 in 2021
1546	! Aim : computes qc and rneb in thermals with cold microphysical considerations
1547	! + for mixed phase boundary layer clouds, calculate ql and qi from
1548	! a stationary MPC model
1549	! IMPORTANT NOTE: we assume iflag_clouth_vert=3
1550	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1551
1552
1553	USE ioipsl_getin_p_mod, ONLY : getin_p
1554	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
1555	USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY
1556	USE phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth
1557
1558	IMPLICIT NONE
1559
1560	#include "YOMCST.h"
1561	#include "YOETHF.h"
1562	#include "FCTTRE.h"
1563	#include "nuage.h"
1564
1565
1566	!------------------------------------------------------------------------------
1567	! Declaration
1568	!------------------------------------------------------------------------------
1569
1570	! INPUT/OUTPUT
1571
1572	INTEGER, INTENT(IN) :: klon,klev,ind2
1573
1574
1575	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K]
1576	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! Virtual potential temp [K]
1577	REAL, DIMENSION(klon), INTENT(IN) :: po ! specific humidity [kg/kg]
1578	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg]
1579	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: fraca ! Fraction of the mesh covered by thermals [0-1]
1580	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk
1581	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! Pressure at layer interfaces [Pa]
1582	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! Pressure at the center of layers [Pa]
1583	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! Liquid temp [K]
1584	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! Liquid potential temp [K]
1585	REAL, DIMENSION(klon,klev), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
1586	REAL, DIMENSION(klon), INTENT(IN) :: zqs ! Saturation specific humidity in the mesh [kg/kg]
1587	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s]
1588	INTEGER, INTENT(IN) :: iflag_mpc_bl ! option flag for mpc boundary layer clouds param.
1589
1590
1591	INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
1592
1593	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
1594	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
1595	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
1596	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
1597	REAL, DIMENSION(klon,klev), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
1598
1599
1600	! LOCAL VARIABLES
1601
1602	INTEGER itap,ind1,l,ig,iter,k
1603	INTEGER iflag_topthermals, niter
1604	LOGICAL falseklon(klon)
1605
1606	REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
1607	REAL sigma1(klon,klev)
1608	REAL sigma2(klon,klev)
1609	REAL qcth(klon,klev)
1610	REAL qcenv(klon,klev)
1611	REAL qctot(klon,klev)
1612	REAL cth(klon,klev)
1613	REAL cenv(klon,klev)
1614	REAL cth_vol(klon,klev)
1615	REAL cenv_vol(klon,klev)
1616	REAL rneb(klon,klev)
1617	REAL zqenv(klon), zqth(klon), zqenvonly(klon)
1618	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
1619	REAL rdd,cppd,Lv
1620	REAL alth,alenv,ath,aenv
1621	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1622	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1623	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1624	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1625	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1626	REAL zdelta,qsatbef,zcor
1627	REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon)
1628	REAL qlbef
1629	REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon)
1630	REAL erf
1631	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
1632	REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
1633	REAL rhodz(klon,klev)
1634	REAL zrho(klon,klev)
1635	REAL dz(klon,klev)
1636	REAL qslenv(klon), qslth(klon)
1637	REAL alenvl, aenvl
1638	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
1639	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
1640	REAL qmax, ttarget, stmp, cout, coutref, fraci
1641	REAL maxi, mini, pas, temp_lim
1642	REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
1643
1644	! Modifty the saturation deficit PDF in thermals
1645	! in the presence of ice crystals
1646	REAL,SAVE :: C_mpc
1647	!$OMP THREADPRIVATE(C_mpc)
1648	REAL, SAVE :: Ni,C_cap,Ei,d_top
1649	!$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top)
1650	! Change the width of the PDF used for vertical subgrid scale heterogeneity
1651	! (J Jouhaud, JL Dufresne, JB Madeleine)
1652	REAL,SAVE :: vert_alpha, vert_alpha_th
1653	!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
1654	REAL,SAVE :: sigma1s_factor=1.1
1655	REAL,SAVE :: sigma1s_power=0.6
1656	REAL,SAVE :: sigma2s_factor=0.09
1657	REAL,SAVE :: sigma2s_power=0.5
1658	REAL,SAVE :: cloudth_ratqsmin=-1.
1659	!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
1660	INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
1661	!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
1662	LOGICAL, SAVE :: firstcall = .TRUE.
1663	!$OMP THREADPRIVATE(firstcall)
1664
1665	CHARACTER (len = 80) :: abort_message
1666	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
1667
1668
1669	!------------------------------------------------------------------------------
1670	! Initialisation
1671	!------------------------------------------------------------------------------
1672
1673
1674	! Few initial checksS
1675
1676	IF (iflag_cloudth_vert.NE.3) THEN
1677	abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3'
1678	CALL abort_physic(routname,abort_message,1)
1679	ENDIF
1680
1681	DO k = 1,klev
1682	DO ind1 = 1, klon
1683	rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2
1684	zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3
1685	dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre
1686	END DO
1687	END DO
1688
1689
1690	sigma1(:,:)=0.
1691	sigma2(:,:)=0.
1692	qcth(:,:)=0.
1693	qcenv(:,:)=0.
1694	qctot(:,:)=0.
1695	qlth(:,ind2)=0.
1696	qith(:,ind2)=0.
1697	wiceth(:,ind2)=0.
1698	rneb(:,:)=0.
1699	qcloud(:)=0.
1700	cth(:,:)=0.
1701	cenv(:,:)=0.
1702	ctot(:,:)=0.
1703	cth_vol(:,:)=0.
1704	cenv_vol(:,:)=0.
1705	ctot_vol(:,:)=0.
1706	falseklon(:)=.false.
1707	qsatmmussig1=0.
1708	qsatmmussig2=0.
1709	rdd=287.04
1710	cppd=1005.7
1711	pi=3.14159
1712	sqrt2pi=sqrt(2.*pi)
1713	sqrt2=sqrt(2.)
1714	sqrtpi=sqrt(pi)
1715	icefrac(:,ind2)=0.
1716	althl=0.
1717	althi=0.
1718	athl=0.
1719	athi=0.
1720	senvl=0.
1721	senvi=0.
1722	sthi=0.
1723	sthl=0.
1724	sthil=0.
1725
1726
1727
1728	IF (firstcall) THEN
1729
1730	vert_alpha=0.5
1731	CALL getin_p('cloudth_vert_alpha',vert_alpha)
1732	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
1733	! The factor used for the thermal is equal to that of the environment
1734	! if nothing is explicitly specified in the def file
1735	vert_alpha_th=vert_alpha
1736	CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
1737	WRITE(,) 'cloudth_vert_alpha_th = ', vert_alpha_th
1738	! Factor used in the calculation of sigma1s
1739	CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
1740	WRITE(,) 'cloudth_sigma1s_factor = ', sigma1s_factor
1741	! Power used in the calculation of sigma1s
1742	CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
1743	WRITE(,) 'cloudth_sigma1s_power = ', sigma1s_power
1744	! Factor used in the calculation of sigma2s
1745	CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
1746	WRITE(,) 'cloudth_sigma2s_factor = ', sigma2s_factor
1747	! Power used in the calculation of sigma2s
1748	CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
1749	WRITE(,) 'cloudth_sigma2s_power = ', sigma2s_power
1750	! Minimum value for the environmental air subgrid water distrib
1751	CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
1752	WRITE(,) 'cloudth_ratqsmin = ', cloudth_ratqsmin
1753	! Remove the dependency to ratqs from the variance of the vertical PDF
1754	CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
1755	WRITE(,) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
1756	! Modifies the PDF in thermals when ice crystals are present
1757	C_mpc=1.e4
1758	CALL getin_p('C_mpc',C_mpc)
1759	WRITE(,) 'C_mpc = ', C_mpc
1760	Ni=2.0e3
1761	CALL getin_p('Ni', Ni)
1762	WRITE(,) 'Ni = ', Ni
1763	Ei=0.5
1764	CALL getin_p('Ei', Ei)
1765	WRITE(,) 'Ei = ', Ei
1766	C_cap=0.5
1767	CALL getin_p('C_cap', C_cap)
1768	WRITE(,) 'C_cap = ', C_cap
1769	d_top=1.2
1770	CALL getin_p('d_top', d_top)
1771	WRITE(,) 'd_top = ', d_top
1772
1773	firstcall=.FALSE.
1774
1775	ENDIF
1776
1777
1778
1779	!-------------------------------------------------------------------------------
1780	! Identify grid points with potential mixed-phase conditions
1781	!-------------------------------------------------------------------------------
1782
1783	temp_lim=RTT-40.0
1784
1785	DO ind1=1,klon
1786	IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) &
1787	.AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2) &
1788	.AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
1789	mpc_bl_points(ind1,ind2)=1
1790	ELSE
1791	mpc_bl_points(ind1,ind2)=0
1792	ENDIF
1793	ENDDO
1794
1795
1796	!-------------------------------------------------------------------------------
1797	! Thermal fraction calculation and standard deviation of the distribution
1798	!-------------------------------------------------------------------------------
1799
1800	! calculation of temperature, humidity and saturation specific humidity
1801
1802	Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
1803	Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
1804	Tbefenvonly(:)=temp(:,ind2)
1805	rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD
1806
1807	zqenv(:)=(po(:)-fraca(:,ind2)zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = aqtth + (1-a)*qtenv
1808	zqth(:)=zqta(:,ind2)
1809	zqenvonly(:)=po(:)
1810
1811
1812	CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
1813
1814	CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
1815	CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
1816
1817	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
1818	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
1819	CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
1820
1821
1822	DO ind1=1,klon
1823
1824
1825	IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
1826
1827
1828	! Environment:
1829
1830
1831	IF (Tbefenv(ind1) .GE. RTT) THEN
1832	Lv=RLVTT
1833	ELSE
1834	Lv=RLSTT
1835	ENDIF
1836
1837
1838	alenv=(0.622Lvzqsatenv(ind1))/(rddzthl(ind1,ind2)*2) !qsl, p84
1839	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
1840	senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84
1841
1842	! For MPCs:
1843	IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
1844	alenvl=(0.622RLVTTqslenv(ind1))/(rddzthl(ind1,ind2)*2)
1845	aenvl=1./(1.+(alenv*Lv/cppd))
1846	senvl=aenvl*(po(ind1)-qslenv(ind1))
1847	senvi=senv
1848	ENDIF
1849
1850
1851	! Thermals:
1852
1853
1854	IF (Tbefth(ind1) .GE. RTT) THEN
1855	Lv=RLVTT
1856	ELSE
1857	Lv=RLSTT
1858	ENDIF
1859
1860
1861	alth=(0.622Lvzqsatth(ind1))/(rddztla(ind1,ind2)*2)
1862	ath=1./(1.+(alth*Lv/cppd))
1863	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))
1864
1865	! For MPCs:
1866	IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
1867	althi=alth
1868	althl=(0.622RLVTTqslth(ind1))/(rddztla(ind1,ind2)*2)
1869	athl=1./(1.+(alth*RLVTT/cppd))
1870	athi=alth
1871	sthl=athl*(zqta(ind1,ind2)-qslth(ind1))
1872	sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2))
1873	sthil=athi*(zqta(ind1,ind2)-qslth(ind1))
1874	ENDIF
1875
1876
1877
1878	!-------------------------------------------------------------------------------
1879	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1880	! Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3
1881	!-------------------------------------------------------------------------------
1882
1883	IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
1884
1885	! Standard deviation of the distributions
1886
1887	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
1888	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
1889
1890	IF (cloudth_ratqsmin>0.) THEN
1891	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1892	ELSE
1893	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1894	ENDIF
1895
1896	sigma1s = sigma1s_fraca + sigma1s_ratqs
1897	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1898
1899
1900	deltasenv=aenvvert_alphasigma1s
1901	deltasth=athvert_alpha_thsigma2s
1902
1903	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1904	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1905	exp_xenv1 = exp(-1.xenv1*2)
1906	exp_xenv2 = exp(-1.xenv2*2)
1907	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1908	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1909	exp_xth1 = exp(-1.xth1*2)
1910	exp_xth2 = exp(-1.xth2*2)
1911
1912	!surface CF
1913
1914	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1915	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1916	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1917
1918
1919	!volume CF and condensed water
1920
1921	!environnement
1922
1923	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1924	IntJ_CF=0.5(1.-1.erf(xenv2))
1925
1926	IF (deltasenv .LT. 1.e-10) THEN
1927	qcenv(ind1,ind2)=IntJ
1928	cenv_vol(ind1,ind2)=IntJ_CF
1929	ELSE
1930	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1931	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1932	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1933	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1934	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1935	qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1936	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1937	ENDIF
1938
1939
1940
1941	!thermals
1942
1943	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1944	IntJ_CF=0.5(1.-1.erf(xth2))
1945
1946	IF (deltasth .LT. 1.e-10) THEN
1947	qcth(ind1,ind2)=IntJ
1948	cth_vol(ind1,ind2)=IntJ_CF
1949	ELSE
1950	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1951	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1952	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1953	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1954	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1955	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1956	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1957	ENDIF
1958
1959	qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
1960	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1961
1962
1963	IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
1964	ctot(ind1,ind2)=0.
1965	ctot_vol(ind1,ind2)=0.
1966	qcloud(ind1)=zqsatenv(ind1)
1967	qincloud(ind1)=0.
1968	ELSE
1969	qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1970	qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
1971	ENDIF
1972
1973
1974	ELSE ! mpc_bl_points>0
1975
1976	! Treat boundary layer mixed phase clouds
1977
1978	! thermals
1979	!=========
1980
1981	! ice phase
1982	!...........
1983
1984	qiceth=0.
1985	deltazlev_mpc=dz(ind1,:)
1986	temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
1987	pres_mpc=pplay(ind1,:)
1988	fraca_mpc=fraca(ind1,:)
1989	snowf_mpc=snowflux(ind1,:)
1990	iflag_topthermals=0
1991	IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
1992	iflag_topthermals = 1
1993	ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
1994	.AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
1995	iflag_topthermals = 2
1996	ELSE
1997	iflag_topthermals = 0
1998	ENDIF
1999
2000	CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
2001	qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
2002
2003	! qmax calculation
2004	sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
2005	deltasth=athlvert_alpha_thsigma2s
2006	xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
2007	xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
2008	exp_xth1 = exp(-1.xth1*2)
2009	exp_xth2 = exp(-1.xth2*2)
2010	IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
2011	IntJ_CF=0.5(1.-1.erf(xth2))
2012	IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
2013	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
2014	IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
2015	IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
2016	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
2017	qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
2018
2019
2020	! Liquid phase
2021	!................
2022	! We account for the effect of ice crystals in thermals on sthl
2023	! and on the width of the distribution
2024
2025	sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
2026	+ (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
2027
2028	sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
2029	/((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
2030	+0.002*zqta(ind1,ind2)
2031	deltasthc=athlvert_alpha_thsigma2sc
2032
2033
2034	xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
2035	xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
2036	exp_xth1 = exp(-1.xth1*2)
2037	exp_xth2 = exp(-1.xth2*2)
2038	IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
2039	IntJ_CF=0.5(1.-1.erf(xth2))
2040	IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
2041	IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
2042	IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
2043	IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
2044	IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
2045	qliqth=IntJ+IntI1+IntI2+IntI3
2046
2047	qlth(ind1,ind2)=MAX(0.,qliqth)
2048
2049	! Condensed water
2050
2051	qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
2052
2053
2054	! consistency with subgrid distribution
2055
2056	IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
2057	fraci=qith(ind1,ind2)/qcth(ind1,ind2)
2058	qcth(ind1,ind2)=qmax
2059	qith(ind1,ind2)=fraci*qmax
2060	qlth(ind1,ind2)=(1.-fraci)*qmax
2061	ENDIF
2062
2063	! Cloud Fraction
2064	!...............
2065	! calculation of qbase which is the value of the water vapor within mixed phase clouds
2066	! such that the total water in cloud = qbase+qliqth+qiceth
2067	! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
2068	! sbase and qbase calculation (note that sbase is wrt liq so negative)
2069	! look for an approximate solution with iteration
2070
2071	ttarget=qcth(ind1,ind2)
2072	mini= athl*(qsith(ind1,ind2)-qslth(ind1))
2073	maxi= 0. !athl(3.sqrt(sigma2s))
2074	niter=20
2075	pas=(maxi-mini)/niter
2076	stmp=mini
2077	sbase=stmp
2078	coutref=1.E6
2079	DO iter=1,niter
2080	cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
2081	+ stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
2082	IF (cout .LT. coutref) THEN
2083	sbase=stmp
2084	coutref=cout
2085	ELSE
2086	stmp=stmp+pas
2087	ENDIF
2088	ENDDO
2089	qbase=MAX(0., sbase/athl+qslth(ind1))
2090
2091	! surface cloud fraction in thermals
2092	cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
2093	cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
2094
2095
2096	!volume cloud fraction in thermals
2097	!to be checked
2098	xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
2099	xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
2100	exp_xth1 = exp(-1.xth1*2)
2101	exp_xth2 = exp(-1.xth2*2)
2102
2103	IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
2104	IntJ_CF=0.5(1.-1.erf(xth2))
2105
2106	IF (deltasth .LT. 1.e-10) THEN
2107	cth_vol(ind1,ind2)=IntJ_CF
2108	ELSE
2109	IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
2110	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
2111	IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
2112	IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
2113	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
2114	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2115	ENDIF
2116	cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
2117
2118
2119
2120	! Environment
2121	!=============
2122	! In the environment/downdrafts, only liquid clouds
2123	! See Shupe et al. 2008, JAS
2124
2125	! standard deviation of the distribution in the environment
2126	sth=sthl
2127	senv=senvl
2128	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
2129	& (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
2130	! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
2131	! in the environement
2132
2133	sigma1s_ratqs=1E-10
2134	IF (cloudth_ratqsmin>0.) THEN
2135	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
2136	ENDIF
2137
2138	sigma1s = sigma1s_fraca + sigma1s_ratqs
2139	deltasenv=aenvlvert_alphasigma1s
2140	xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
2141	xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
2142	exp_xenv1 = exp(-1.xenv1*2)
2143	exp_xenv2 = exp(-1.xenv2*2)
2144
2145	!surface CF
2146	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
2147
2148	!volume CF and condensed water
2149	IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
2150	IntJ_CF=0.5(1.-1.erf(xenv2))
2151
2152	IF (deltasenv .LT. 1.e-10) THEN
2153	qcenv(ind1,ind2)=IntJ
2154	cenv_vol(ind1,ind2)=IntJ_CF
2155	ELSE
2156	IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
2157	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
2158	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
2159	IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
2160	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
2161	qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
2162	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
2163	ENDIF
2164
2165	qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
2166	cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
2167
2168
2169
2170	! Thermals + environment
2171	! =======================
2172	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
2173	qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
2174	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2175	IF (qcth(ind1,ind2) .GT. 0) THEN
2176	icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
2177	/(fraca(ind1,ind2)*qcth(ind1,ind2) &
2178	+(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
2179	icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
2180	ELSE
2181	icefrac(ind1,ind2)=0.
2182	ENDIF
2183
2184	IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
2185	ctot(ind1,ind2)=0.
2186	ctot_vol(ind1,ind2)=0.
2187	qincloud(ind1)=0.
2188	qcloud(ind1)=zqsatenv(ind1)
2189	ELSE
2190	qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
2191	+(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
2192	qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
2193	+(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
2194	ENDIF
2195
2196	ENDIF ! mpc_bl_points
2197
2198
2199	ELSE ! gaussian for environment only
2200
2201
2202
2203
2204	IF (Tbefenvonly(ind1) .GE. RTT) THEN
2205	Lv=RLVTT
2206	ELSE
2207	Lv=RLSTT
2208	ENDIF
2209
2210
2211	zthl(ind1,ind2)=temp(ind1,ind2)(101325./paprs(ind1,ind2))*(rdd/cppd)
2212	alenv=(0.622Lvzqsatenvonly(ind1))/(rddzthl(ind1,ind2)*2)
2213	aenv=1./(1.+(alenv*Lv/cppd))
2214	senv=aenv*(po(ind1)-zqsatenvonly(ind1))
2215	sth=0.
2216
2217	sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
2218	sigma2s=0.
2219
2220	sqrt2pi=sqrt(2.*pi)
2221	xenv=senv/(sqrt(2.)*sigma1s)
2222	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
2223	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
2224	qctot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
2225
2226	IF (ctot(ind1,ind2).LT.1.e-3) THEN
2227	ctot(ind1,ind2)=0.
2228	qcloud(ind1)=zqsatenvonly(ind1)
2229	qincloud(ind1)=0.
2230	ELSE
2231	qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
2232	qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
2233	ENDIF
2234
2235
2236	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
2237
2238	! Outputs used to check the PDFs
2239	cloudth_senv(ind1,ind2) = senv
2240	cloudth_sth(ind1,ind2) = sth
2241	cloudth_sigmaenv(ind1,ind2) = sigma1s
2242	cloudth_sigmath(ind1,ind2) = sigma2s
2243
2244
2245	ENDDO !loop on klon
2246
2247	! Calcule ice fall velocity in thermals
2248
2249	CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2))
2250
2251	RETURN
2252
2253
2254	END SUBROUTINE cloudth_mpc
2255
2256	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2257	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)
2258
2259	! parameterization of ice for boundary
2260	! layer mixed-phase clouds assuming a stationary system
2261	!
2262	! Note that vapor deposition on ice crystals and riming of liquid droplets
2263	! depend on the ice number concentration Ni
2264	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
2265	! and use values from Hong et al. 2004, MWR for instance
2266	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
2267	! One could also think of a more complex expression of Ni;
2268	! function of qi, T, the concentration in aerosols or INP ..
2269	! Here we prefer fixing Ni to a tuning parameter
2270	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
2271	! in Mioche et al. 2017
2272	!
2273	!
2274	! References:
2275	!------------
2276	! This parameterization is thoroughly described in Vignon et al.
2277	!
2278	! More specifically
2279	! for the Water vapor deposition process:
2280	!
2281	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
2282	! ment of stratiform cloudfs and precipitation in large-scale
2283	! models. I: Description and evaluation of the microphysical
2284	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
2285	!
2286	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
2287	! stratiform cloud microphysics scheme in the NCAR Com-
2288	! munity Atmosphere Model (CAM3). Part I: Description and
2289	! numerical tests. J. Climate, 21, 3642–3659
2290	!
2291	! for the Riming process:
2292	!
2293	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
2294	! scale structure and organization of clouds and precipitation in
2295	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
2296	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
2297	!
2298	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
2299	! forecasts of winter precipitation using an improved bulk
2300	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
2301	! forecasts of winter precipitation using an improved bulk
2302	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
2303	!
2304	! For the formation of clouds by thermals:
2305	!
2306	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
2307	! the Atmospheric Sciences, 65, 407–425.
2308	!
2309	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
2310	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
2311	!
2312	!
2313	!
2314	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
2315	!=============================================================================
2316
2317	USE ioipsl_getin_p_mod, ONLY : getin_p
2318	USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
2319
2320	IMPLICIT none
2321
2322	INCLUDE "YOMCST.h"
2323	INCLUDE "nuage.h"
2324
2325	INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
2326	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
2327	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
2328	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
2329	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
2330	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
2331	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
2332	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
2333	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
2334	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
2335	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
2336	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
2337	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
2338	REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
2339	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
2340
2341
2342	INTEGER ind2p1,ind2p2
2343	REAL rho(klev)
2344	REAL unsurtaudet, unsurtaustardep, unsurtaurim
2345	REAL qi, AA, BB, Ka, Dv, rhoi
2346	REAL p0, t0, fp1, fp2
2347	REAL alpha, flux_term
2348	REAL det_term, precip_term, rim_term, dep_term
2349
2350
2351	ind2p1=ind2+1
2352	ind2p2=ind2+2
2353
2354	rho=pres/temp/RD ! air density kg/m3
2355
2356	Ka=2.4e-2 ! thermal conductivity of the air, SI
2357	p0=101325.0 ! ref pressure
2358	T0=273.15 ! ref temp
2359	rhoi=500.0 ! cloud ice density following Reisner et al. 1998
2360	alpha=700. ! fallvelocity param
2361
2362
2363	IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
2364
2365	Dv=0.00010.211(p0/pres(ind2))((temp(ind2)/T0)*1.94) ! water vapor diffusivity in air, SI
2366
2367	! Detrainment term:
2368
2369	unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
2370
2371
2372	! Deposition term
2373	AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
2374	BB=1./(rho(ind2)Dvqsith(ind2))
2375	unsurtaustardep=C_cap(Ni0.66)(qth(ind2)-qsith(ind2))/qsith(ind2) &
2376	4.RPI/(AA+BB)(6.rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
2377
2378	! Riming term neglected at this level
2379	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
2380
2381	qi=fraca(ind2)unsurtaustardep/MAX((d_topunsurtaudet),1E-12)
2382	qi=MAX(qi,0.)**(3./2.)
2383
2384	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
2385
2386	Dv=0.00010.211(p0/pres(ind2p1))((temp(ind2p1)/T0)*1.94) ! water vapor diffusivity in air, SI
2387
2388	! Detrainment term:
2389
2390	unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
2391	det_term=-unsurtaudetqith(ind2p1)rho(ind2p1)
2392
2393
2394	! Deposition term
2395
2396	AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
2397	BB=1./(rho(ind2p1)Dvqsith(ind2p1))
2398	unsurtaustardep=C_cap(Ni0.66)(qth(ind2p1)-qsith(ind2p1)) &
2399	/qsith(ind2p1)4.RPI/(AA+BB) &
2400	(6.rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
2401	dep_term=rho(ind2p1)fraca(ind2p1)(qith(ind2p1)*0.33)unsurtaustardep
2402
2403	! Riming term
2404
2405	unsurtaurim=rho(ind2p1)alpha3./rhoi/2.Eiqlth(ind2p1)((p0/pres(ind2p1))*0.4)
2406	rim_term=rho(ind2p1)fraca(ind2p1)qith(ind2p1)*unsurtaurim
2407
2408	! Precip term
2409
2410	! We assume that there is no solid precipitation outside thermals
2411	! so the precipitation flux within thermals is equal to the precipitation flux
2412	! at mesh-scale divided by thermals fraction
2413
2414	fp2=0.
2415	fp1=0.
2416	IF (fraca(ind2p1) .GT. 0.) THEN
2417	fp2=-qith(ind2p2)rho(ind2p2)vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
2418	fp1=-qith(ind2p1)rho(ind2p1)vith(ind2p1)*fraca(ind2p1)
2419	ENDIF
2420
2421	precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)
2422
2423	! Calculation in a top-to-bottom loop
2424
2425	IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
2426	qi= 1./fm_therm(ind1,ind2p1)* &
2427	(deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
2428	fm_therm(ind1,ind2p2)*(qith(ind2p1)))
2429	END IF
2430
2431	ENDIF ! top thermals
2432
2433	qith(ind2)=MAX(0.,qi)
2434
2435	RETURN
2436
2437	END SUBROUTINE ICE_MPC_BL_CLOUDS
2438
2439
2440
2441
2442	END MODULE cloudth_mod
2443
2444
2445
2446

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

Download in other formats:

Original Format