Context Navigation

cloudth_mod.F90 @ 4190

Last change on this file since 4190 was 4013, checked in by lguez, 3 years ago
Sync latest trunk changes to Ocean_skin
File size: 91.2 KB

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

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

Download in other formats:

Original Format