Context Navigation

cloudth_mod.F90 @ 5304

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

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

Download in other formats:

Original Format