Context Navigation

cloudth_mod.F90 @ 4639

Last change on this file since 4639 was 4623, checked in by fhourdin, 12 months ago

File size: 93.5 KB

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

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

Download in other formats:

Original Format