Context Navigation

lmdz_cloudth.F90 @ 5103

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

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

Download in other formats:

Original Format