Context Navigation

lmdz_cloudth.F90 @ 5157

Last change on this file since 5157 was 5153, checked in by abarral, 3 months ago
Revert FCTTRE to INCLUDE to assess impact of inlining
File size: 92.5 KB

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

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

Download in other formats:

Original Format