Context Navigation

lmdz_cloudth.F90 @ 5133

Last change on this file since 5133 was 5117, checked in by abarral, 5 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 86.7 KB

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

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

Download in other formats:

Original Format