Context Navigation

cloudth_mod.F90 @ 3493

Last change on this file since 3493 was 3493, checked in by idelkadi, 5 years ago
Amelioration de la representation des nuages bas (Jean Jouhaud) : Calcul des nouveaux ecarts types sigmas de la bi-Gaussienne en prenant en compte l'epaisseur des mailles, CF_vol avec ces sigmas, et CF_surf avec la methode de Brooks
File size: 56.5 KB

Line
1	MODULE cloudth_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE cloudth(ngrid,klev,ind2, &
8	& ztv,po,zqta,fraca, &
9	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
10	& ratqs,zqs,t)
11
12
13	IMPLICIT NONE
14
15
16	!===========================================================================
17	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
18	! Date : 25 Mai 2010
19	! Objet : calcule les valeurs de qc et rneb dans les thermiques
20	!===========================================================================
21
22
23	#include "YOMCST.h"
24	#include "YOETHF.h"
25	#include "FCTTRE.h"
26	#include "thermcell.h"
27	#include "nuage.h"
28
29	INTEGER itap,ind1,ind2
30	INTEGER ngrid,klev,klon,l,ig
31
32	REAL ztv(ngrid,klev)
33	REAL po(ngrid)
34	REAL zqenv(ngrid)
35	REAL zqta(ngrid,klev)
36
37	REAL fraca(ngrid,klev+1)
38	REAL zpspsk(ngrid,klev)
39	REAL paprs(ngrid,klev+1)
40	REAL pplay(ngrid,klev)
41	REAL ztla(ngrid,klev)
42	REAL zthl(ngrid,klev)
43
44	REAL zqsatth(ngrid,klev)
45	REAL zqsatenv(ngrid,klev)
46
47
48	REAL sigma1(ngrid,klev)
49	REAL sigma2(ngrid,klev)
50	REAL qlth(ngrid,klev)
51	REAL qlenv(ngrid,klev)
52	REAL qltot(ngrid,klev)
53	REAL cth(ngrid,klev)
54	REAL cenv(ngrid,klev)
55	REAL ctot(ngrid,klev)
56	REAL rneb(ngrid,klev)
57	REAL t(ngrid,klev)
58	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
59	REAL rdd,cppd,Lv
60	REAL alth,alenv,ath,aenv
61	REAL sth,senv,sigma1s,sigma2s,xth,xenv
62	REAL Tbef,zdelta,qsatbef,zcor
63	REAL qlbef
64	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
65
66	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
67	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
68	REAL zqs(ngrid), qcloud(ngrid)
69	REAL erf
70
71
72
73
74	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75	! Gestion de deux versions de cloudth
76	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77
78	IF (iflag_cloudth_vert.GE.1) THEN
79	CALL cloudth_vert(ngrid,klev,ind2, &
80	& ztv,po,zqta,fraca, &
81	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
82	& ratqs,zqs,t)
83	RETURN
84	ENDIF
85	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86
87
88	!-------------------------------------------------------------------------------
89	! Initialisation des variables r?elles
90	!-------------------------------------------------------------------------------
91	sigma1(:,:)=0.
92	sigma2(:,:)=0.
93	qlth(:,:)=0.
94	qlenv(:,:)=0.
95	qltot(:,:)=0.
96	rneb(:,:)=0.
97	qcloud(:)=0.
98	cth(:,:)=0.
99	cenv(:,:)=0.
100	ctot(:,:)=0.
101	qsatmmussig1=0.
102	qsatmmussig2=0.
103	rdd=287.04
104	cppd=1005.7
105	pi=3.14159
106	Lv=2.5e6
107	sqrt2pi=sqrt(2.*pi)
108
109
110
111	!-------------------------------------------------------------------------------
112	! Calcul de la fraction du thermique et des ?cart-types des distributions
113	!-------------------------------------------------------------------------------
114	do ind1=1,ngrid
115
116	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
117
118	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
119
120
121	! zqenv(ind1)=po(ind1)
122	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
123	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
124	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
125	qsatbef=MIN(0.5,qsatbef)
126	zcor=1./(1.-retv*qsatbef)
127	qsatbef=qsatbef*zcor
128	zqsatenv(ind1,ind2)=qsatbef
129
130
131
132
133	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
134	aenv=1./(1.+(alenv*Lv/cppd))
135	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
136
137
138
139
140	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
141	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
142	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
143	qsatbef=MIN(0.5,qsatbef)
144	zcor=1./(1.-retv*qsatbef)
145	qsatbef=qsatbef*zcor
146	zqsatth(ind1,ind2)=qsatbef
147
148	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
149	ath=1./(1.+(alth*Lv/cppd))
150	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
151
152
153
154	!------------------------------------------------------------------------------
155	! Calcul des ?cart-types pour s
156	!------------------------------------------------------------------------------
157
158	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
159	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
160	! if (paprs(ind1,ind2).gt.90000) then
161	! ratqs(ind1,ind2)=0.002
162	! else
163	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
164	! endif
165	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
166	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
167	! sigma1s=ratqs(ind1,ind2)*po(ind1)
168	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
169
170	!------------------------------------------------------------------------------
171	! Calcul de l'eau condens?e et de la couverture nuageuse
172	!------------------------------------------------------------------------------
173	sqrt2pi=sqrt(2.*pi)
174	xth=sth/(sqrt(2.)*sigma2s)
175	xenv=senv/(sqrt(2.)*sigma1s)
176	cth(ind1,ind2)=0.5(1.+1.erf(xth))
177	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
178	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
179
180	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
181	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
182	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
183
184	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185	if (ctot(ind1,ind2).lt.1.e-10) then
186	ctot(ind1,ind2)=0.
187	qcloud(ind1)=zqsatenv(ind1,ind2)
188
189	else
190
191	ctot(ind1,ind2)=ctot(ind1,ind2)
192	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
193
194	endif
195
196
197	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
198
199
200	else ! gaussienne environnement seule
201
202	zqenv(ind1)=po(ind1)
203	Tbef=t(ind1,ind2)
204	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
205	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
206	qsatbef=MIN(0.5,qsatbef)
207	zcor=1./(1.-retv*qsatbef)
208	qsatbef=qsatbef*zcor
209	zqsatenv(ind1,ind2)=qsatbef
210
211
212	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
213	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
214	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
215	aenv=1./(1.+(alenv*Lv/cppd))
216	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
217
218
219	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
220
221	sqrt2pi=sqrt(2.*pi)
222	xenv=senv/(sqrt(2.)*sigma1s)
223	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
224	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
225
226	if (ctot(ind1,ind2).lt.1.e-3) then
227	ctot(ind1,ind2)=0.
228	qcloud(ind1)=zqsatenv(ind1,ind2)
229
230	else
231
232	ctot(ind1,ind2)=ctot(ind1,ind2)
233	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
234
235	endif
236
237
238
239
240
241
242	endif
243	enddo
244
245	return
246	! end
247	END SUBROUTINE cloudth
248
249
250
251	!===========================================================================
252	SUBROUTINE cloudth_vert(ngrid,klev,ind2, &
253	& ztv,po,zqta,fraca, &
254	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
255	& ratqs,zqs,t)
256
257	!===========================================================================
258	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
259	! Date : 25 Mai 2010
260	! Objet : calcule les valeurs de qc et rneb dans les thermiques
261	!===========================================================================
262
263
264	USE ioipsl_getin_p_mod, ONLY : getin_p
265
266	IMPLICIT NONE
267
268	#include "YOMCST.h"
269	#include "YOETHF.h"
270	#include "FCTTRE.h"
271	#include "thermcell.h"
272	#include "nuage.h"
273
274	INTEGER itap,ind1,ind2
275	INTEGER ngrid,klev,klon,l,ig
276
277	REAL ztv(ngrid,klev)
278	REAL po(ngrid)
279	REAL zqenv(ngrid)
280	REAL zqta(ngrid,klev)
281
282	REAL fraca(ngrid,klev+1)
283	REAL zpspsk(ngrid,klev)
284	REAL paprs(ngrid,klev+1)
285	REAL pplay(ngrid,klev)
286	REAL ztla(ngrid,klev)
287	REAL zthl(ngrid,klev)
288
289	REAL zqsatth(ngrid,klev)
290	REAL zqsatenv(ngrid,klev)
291
292
293	REAL sigma1(ngrid,klev)
294	REAL sigma2(ngrid,klev)
295	REAL qlth(ngrid,klev)
296	REAL qlenv(ngrid,klev)
297	REAL qltot(ngrid,klev)
298	REAL cth(ngrid,klev)
299	REAL cenv(ngrid,klev)
300	REAL ctot(ngrid,klev)
301	REAL rneb(ngrid,klev)
302	REAL t(ngrid,klev)
303	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
304	REAL rdd,cppd,Lv,sqrt2,sqrtpi
305	REAL alth,alenv,ath,aenv
306	REAL sth,senv,sigma1s,sigma2s,xth,xenv
307	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
308	REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
309	REAL Tbef,zdelta,qsatbef,zcor
310	REAL qlbef
311	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
312	! Change the width of the PDF used for vertical subgrid scale heterogeneity
313	! (J Jouhaud, JL Dufresne, JB Madeleine)
314	REAL,SAVE :: vert_alpha
315	!$OMP THREADPRIVATE(vert_alpha)
316	LOGICAL, SAVE :: firstcall = .TRUE.
317	!$OMP THREADPRIVATE(firstcall)
318
319	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
320	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
321	REAL zqs(ngrid), qcloud(ngrid)
322	REAL erf
323
324	!------------------------------------------------------------------------------
325	! Initialisation des variables r?elles
326	!------------------------------------------------------------------------------
327	sigma1(:,:)=0.
328	sigma2(:,:)=0.
329	qlth(:,:)=0.
330	qlenv(:,:)=0.
331	qltot(:,:)=0.
332	rneb(:,:)=0.
333	qcloud(:)=0.
334	cth(:,:)=0.
335	cenv(:,:)=0.
336	ctot(:,:)=0.
337	qsatmmussig1=0.
338	qsatmmussig2=0.
339	rdd=287.04
340	cppd=1005.7
341	pi=3.14159
342	Lv=2.5e6
343	sqrt2pi=sqrt(2.*pi)
344	sqrt2=sqrt(2.)
345	sqrtpi=sqrt(pi)
346
347	IF (firstcall) THEN
348	vert_alpha=0.5
349	CALL getin_p('cloudth_vert_alpha',vert_alpha)
350	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
351	firstcall=.FALSE.
352	ENDIF
353
354	!-------------------------------------------------------------------------------
355	! Calcul de la fraction du thermique et des ?cart-types des distributions
356	!-------------------------------------------------------------------------------
357	do ind1=1,ngrid
358
359	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
360
361	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
362
363
364	! zqenv(ind1)=po(ind1)
365	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
366	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
367	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
368	qsatbef=MIN(0.5,qsatbef)
369	zcor=1./(1.-retv*qsatbef)
370	qsatbef=qsatbef*zcor
371	zqsatenv(ind1,ind2)=qsatbef
372
373
374
375
376	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
377	aenv=1./(1.+(alenv*Lv/cppd))
378	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
379
380
381
382
383	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
384	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
385	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
386	qsatbef=MIN(0.5,qsatbef)
387	zcor=1./(1.-retv*qsatbef)
388	qsatbef=qsatbef*zcor
389	zqsatth(ind1,ind2)=qsatbef
390
391	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
392	ath=1./(1.+(alth*Lv/cppd))
393	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
394
395
396
397	!------------------------------------------------------------------------------
398	! Calcul des ?cart-types pour s
399	!------------------------------------------------------------------------------
400
401	sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
402	sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
403	! if (paprs(ind1,ind2).gt.90000) then
404	! ratqs(ind1,ind2)=0.002
405	! else
406	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
407	! endif
408	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
409	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
410	! sigma1s=ratqs(ind1,ind2)*po(ind1)
411	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
412
413	!------------------------------------------------------------------------------
414	! Calcul de l'eau condens?e et de la couverture nuageuse
415	!------------------------------------------------------------------------------
416	sqrt2pi=sqrt(2.*pi)
417	xth=sth/(sqrt(2.)*sigma2s)
418	xenv=senv/(sqrt(2.)*sigma1s)
419	cth(ind1,ind2)=0.5(1.+1.erf(xth))
420	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
421	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
422
423	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
424	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
425	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
426
427	IF (iflag_cloudth_vert == 1) THEN
428	!-------------------------------------------------------------------------------
429	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
430	!-------------------------------------------------------------------------------
431	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
432	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
433	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
434	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
435	! deltasenv=aenv0.01po(ind1)
436	! deltasth=ath0.01zqta(ind1,ind2)
437	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
438	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
439	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
440	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
441	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
442	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
443
444	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
445	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
446	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
447
448	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
449	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
450	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
451	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
452
453	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
454	! qlenv(ind1,ind2)=IntJ
455	! print*, qlenv(ind1,ind2),'VERIF EAU'
456
457
458	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
459	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
460	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
461	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
462	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
463	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
464	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
465	! qlth(ind1,ind2)=IntJ
466	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
467	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
468
469	ELSE IF (iflag_cloudth_vert == 2) THEN
470
471	!-------------------------------------------------------------------------------
472	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
473	!-------------------------------------------------------------------------------
474	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
475	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
476	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
477	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
478	deltasenv=aenvvert_alphasigma1s
479	deltasth=athvert_alphasigma2s
480
481	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
482	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
483	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
484	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
485	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
486	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
487
488	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
489	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
490	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
491
492	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)exp(-1.xenv2**2)
493	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
494	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
495	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))(exp(-1.xenv1*2)-exp(-1.xenv2**2))
496
497	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
498	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
499	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
500
501	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
502	! qlenv(ind1,ind2)=IntJ
503	! print*, qlenv(ind1,ind2),'VERIF EAU'
504
505	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)exp(-1.xth2**2)
506	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
507	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
508	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))(exp(-1.xth1*2)-exp(-1.xth2**2))
509
510	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
511	! qlth(ind1,ind2)=IntJ
512	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
513	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
514
515
516
517
518	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
519
520	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
521
522	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
523	ctot(ind1,ind2)=0.
524	qcloud(ind1)=zqsatenv(ind1,ind2)
525
526	else
527
528	ctot(ind1,ind2)=ctot(ind1,ind2)
529	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
530	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
531	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
532
533	endif
534
535
536
537	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
538
539
540	else ! gaussienne environnement seule
541
542	zqenv(ind1)=po(ind1)
543	Tbef=t(ind1,ind2)
544	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
545	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
546	qsatbef=MIN(0.5,qsatbef)
547	zcor=1./(1.-retv*qsatbef)
548	qsatbef=qsatbef*zcor
549	zqsatenv(ind1,ind2)=qsatbef
550
551
552	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
553	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
554	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
555	aenv=1./(1.+(alenv*Lv/cppd))
556	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
557
558
559	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
560
561	sqrt2pi=sqrt(2.*pi)
562	xenv=senv/(sqrt(2.)*sigma1s)
563	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
564	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
565
566	if (ctot(ind1,ind2).lt.1.e-3) then
567	ctot(ind1,ind2)=0.
568	qcloud(ind1)=zqsatenv(ind1,ind2)
569
570	else
571
572	ctot(ind1,ind2)=ctot(ind1,ind2)
573	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
574
575	endif
576
577
578
579
580
581
582	endif
583	enddo
584
585	return
586	! end
587	END SUBROUTINE cloudth_vert
588
589
590
591
592	SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
593	& ztv,po,zqta,fraca, &
594	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
595	& ratqs,zqs,t)
596
597
598	IMPLICIT NONE
599
600
601	!===========================================================================
602	! Author : Arnaud Octavio Jam (LMD/CNRS)
603	! Date : 25 Mai 2010
604	! Objet : calcule les valeurs de qc et rneb dans les thermiques
605	!===========================================================================
606
607
608	#include "YOMCST.h"
609	#include "YOETHF.h"
610	#include "FCTTRE.h"
611	#include "thermcell.h"
612	#include "nuage.h"
613
614	INTEGER itap,ind1,ind2
615	INTEGER ngrid,klev,klon,l,ig
616
617	REAL ztv(ngrid,klev)
618	REAL po(ngrid)
619	REAL zqenv(ngrid)
620	REAL zqta(ngrid,klev)
621
622	REAL fraca(ngrid,klev+1)
623	REAL zpspsk(ngrid,klev)
624	REAL paprs(ngrid,klev+1)
625	REAL pplay(ngrid,klev)
626	REAL ztla(ngrid,klev)
627	REAL zthl(ngrid,klev)
628
629	REAL zqsatth(ngrid,klev)
630	REAL zqsatenv(ngrid,klev)
631
632	REAL sigma1(ngrid,klev)
633	REAL sigma2(ngrid,klev)
634	REAL qlth(ngrid,klev)
635	REAL qlenv(ngrid,klev)
636	REAL qltot(ngrid,klev)
637	REAL cth(ngrid,klev)
638	REAL cenv(ngrid,klev)
639	REAL ctot(ngrid,klev)
640	REAL cth_vol(ngrid,klev)
641	REAL cenv_vol(ngrid,klev)
642	REAL ctot_vol(ngrid,klev)
643	REAL rneb(ngrid,klev)
644	REAL t(ngrid,klev)
645	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
646	REAL rdd,cppd,Lv
647	REAL alth,alenv,ath,aenv
648	REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
649	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
650	REAL Tbef,zdelta,qsatbef,zcor
651	REAL qlbef
652	REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
653	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
654	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
655	REAL zqs(ngrid), qcloud(ngrid)
656	REAL erf
657
658
659
660	IF (iflag_cloudth_vert.GE.1) THEN
661	CALL cloudth_vert_v3(ngrid,klev,ind2, &
662	& ztv,po,zqta,fraca, &
663	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
664	& ratqs,zqs,t)
665	RETURN
666	ENDIF
667	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
668
669
670	!-------------------------------------------------------------------------------
671	! Initialisation des variables r?elles
672	!-------------------------------------------------------------------------------
673	sigma1(:,:)=0.
674	sigma2(:,:)=0.
675	qlth(:,:)=0.
676	qlenv(:,:)=0.
677	qltot(:,:)=0.
678	rneb(:,:)=0.
679	qcloud(:)=0.
680	cth(:,:)=0.
681	cenv(:,:)=0.
682	ctot(:,:)=0.
683	cth_vol(:,:)=0.
684	cenv_vol(:,:)=0.
685	ctot_vol(:,:)=0.
686	qsatmmussig1=0.
687	qsatmmussig2=0.
688	rdd=287.04
689	cppd=1005.7
690	pi=3.14159
691	Lv=2.5e6
692	sqrt2pi=sqrt(2.*pi)
693	sqrt2=sqrt(2.)
694	sqrtpi=sqrt(pi)
695
696
697	!-------------------------------------------------------------------------------
698	! Cloud fraction in the thermals and standard deviation of the PDFs
699	!-------------------------------------------------------------------------------
700	do ind1=1,ngrid
701
702	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
703
704	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
705
706	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
707	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
708	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
709	qsatbef=MIN(0.5,qsatbef)
710	zcor=1./(1.-retv*qsatbef)
711	qsatbef=qsatbef*zcor
712	zqsatenv(ind1,ind2)=qsatbef
713
714
715	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
716	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
717	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
718
719	!po = qt de l'environnement ET des thermique
720	!zqenv = qt environnement
721	!zqsatenv = qsat environnement
722	!zthl = Tl environnement
723
724
725	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
726	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
727	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
728	qsatbef=MIN(0.5,qsatbef)
729	zcor=1./(1.-retv*qsatbef)
730	qsatbef=qsatbef*zcor
731	zqsatth(ind1,ind2)=qsatbef
732
733	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
734	ath=1./(1.+(alth*Lv/cppd)) !al, p84
735	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
736
737	!zqta = qt thermals
738	!zqsatth = qsat thermals
739	!ztla = Tl thermals
740
741	!------------------------------------------------------------------------------
742	! s standard deviations
743	!------------------------------------------------------------------------------
744
745	! tests
746	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
747	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
748	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
749	! final option
750	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
751	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
752
753	!------------------------------------------------------------------------------
754	! Condensed water and cloud cover
755	!------------------------------------------------------------------------------
756	xth=sth/(sqrt2*sigma2s)
757	xenv=senv/(sqrt2*sigma1s)
758	cth(ind1,ind2)=0.5(1.+1.erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
759	cenv(ind1,ind2)=0.5(1.+1.erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
760	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
761	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
762
763	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth(ind1,ind2))
764	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
765	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
766
767	if (ctot(ind1,ind2).lt.1.e-10) then
768	ctot(ind1,ind2)=0.
769	qcloud(ind1)=zqsatenv(ind1,ind2)
770	else
771	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
772	endif
773
774	else ! Environnement only, follow the if l.110
775
776	zqenv(ind1)=po(ind1)
777	Tbef=t(ind1,ind2)
778	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
779	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
780	qsatbef=MIN(0.5,qsatbef)
781	zcor=1./(1.-retv*qsatbef)
782	qsatbef=qsatbef*zcor
783	zqsatenv(ind1,ind2)=qsatbef
784
785	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
786	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
787	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
788	aenv=1./(1.+(alenv*Lv/cppd))
789	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
790
791	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
792
793	xenv=senv/(sqrt2*sigma1s)
794	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
795	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
796	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
797
798	if (ctot(ind1,ind2).lt.1.e-3) then
799	ctot(ind1,ind2)=0.
800	qcloud(ind1)=zqsatenv(ind1,ind2)
801	else
802	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
803	endif
804
805
806	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
807	enddo ! from the loop on ngrid l.108
808	return
809	! end
810	END SUBROUTINE cloudth_v3
811
812
813
814	!===========================================================================
815	SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
816	& ztv,po,zqta,fraca, &
817	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
818	& ratqs,zqs,t)
819
820	!===========================================================================
821	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
822	! Date : 25 Mai 2010
823	! Objet : calcule les valeurs de qc et rneb dans les thermiques
824	!===========================================================================
825
826
827	USE ioipsl_getin_p_mod, ONLY : getin_p
828	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
829	& cloudth_sigmath,cloudth_sigmaenv
830
831	IMPLICIT NONE
832
833	#include "YOMCST.h"
834	#include "YOETHF.h"
835	#include "FCTTRE.h"
836	#include "thermcell.h"
837	#include "nuage.h"
838
839	INTEGER itap,ind1,ind2
840	INTEGER ngrid,klev,klon,l,ig
841
842	REAL ztv(ngrid,klev)
843	REAL po(ngrid)
844	REAL zqenv(ngrid)
845	REAL zqta(ngrid,klev)
846
847	REAL fraca(ngrid,klev+1)
848	REAL zpspsk(ngrid,klev)
849	REAL paprs(ngrid,klev+1)
850	REAL pplay(ngrid,klev)
851	REAL ztla(ngrid,klev)
852	REAL zthl(ngrid,klev)
853
854	REAL zqsatth(ngrid,klev)
855	REAL zqsatenv(ngrid,klev)
856
857	REAL sigma1(ngrid,klev)
858	REAL sigma2(ngrid,klev)
859	REAL qlth(ngrid,klev)
860	REAL qlenv(ngrid,klev)
861	REAL qltot(ngrid,klev)
862	REAL cth(ngrid,klev)
863	REAL cenv(ngrid,klev)
864	REAL ctot(ngrid,klev)
865	REAL cth_vol(ngrid,klev)
866	REAL cenv_vol(ngrid,klev)
867	REAL ctot_vol(ngrid,klev)
868	REAL rneb(ngrid,klev)
869	REAL t(ngrid,klev)
870	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
871	REAL rdd,cppd,Lv
872	REAL alth,alenv,ath,aenv
873	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
874	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
875	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
876	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
877	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
878	REAL Tbef,zdelta,qsatbef,zcor
879	REAL qlbef
880	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
881	! Change the width of the PDF used for vertical subgrid scale heterogeneity
882	! (J Jouhaud, JL Dufresne, JB Madeleine)
883	REAL,SAVE :: vert_alpha, vert_alpha_th
884	!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
885	REAL,SAVE :: sigma1s_factor=1.1
886	REAL,SAVE :: sigma1s_power=0.6
887	REAL,SAVE :: cloudth_ratqsmin=-1.
888	!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,cloudth_ratqsmin)
889	INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
890	!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
891
892	LOGICAL, SAVE :: firstcall = .TRUE.
893	!$OMP THREADPRIVATE(firstcall)
894
895	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
896	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
897	REAL zqs(ngrid), qcloud(ngrid)
898	REAL erf
899
900	REAL rhodz(ngrid,klev)
901	REAL zrho(ngrid,klev)
902	REAL dz(ngrid,klev)
903
904	DO ind1 = 1, ngrid
905	!Layer calculation
906	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
907	zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
908	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
909	END DO
910
911
912	!------------------------------------------------------------------------------
913	! Initialize
914	!------------------------------------------------------------------------------
915	sigma1(:,:)=0.
916	sigma2(:,:)=0.
917	qlth(:,:)=0.
918	qlenv(:,:)=0.
919	qltot(:,:)=0.
920	rneb(:,:)=0.
921	qcloud(:)=0.
922	cth(:,:)=0.
923	cenv(:,:)=0.
924	ctot(:,:)=0.
925	cth_vol(:,:)=0.
926	cenv_vol(:,:)=0.
927	ctot_vol(:,:)=0.
928	qsatmmussig1=0.
929	qsatmmussig2=0.
930	rdd=287.04
931	cppd=1005.7
932	pi=3.14159
933	Lv=2.5e6
934	sqrt2pi=sqrt(2.*pi)
935	sqrt2=sqrt(2.)
936	sqrtpi=sqrt(pi)
937
938	IF (firstcall) THEN
939	vert_alpha=0.5
940	CALL getin_p('cloudth_vert_alpha',vert_alpha)
941	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
942	! The factor used for the thermal is equal to that of the environment
943	! if nothing is explicitly specified in the def file
944	vert_alpha_th=vert_alpha
945	CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
946	WRITE(,) 'cloudth_vert_alpha_th = ', vert_alpha_th
947	! Factor used in the calculation of sigma1s
948	CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
949	WRITE(,) 'cloudth_sigma1s_factor = ', sigma1s_factor
950	! Power used in the calculation of sigma1s
951	CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
952	WRITE(,) 'cloudth_sigma1s_power = ', sigma1s_power
953	! Minimum value for the environmental air subgrid water distrib
954	CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
955	WRITE(,) 'cloudth_ratqsmin = ', cloudth_ratqsmin
956	! Remove the dependency to ratqs from the variance of the vertical PDF
957	CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
958	WRITE(,) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
959
960	firstcall=.FALSE.
961	ENDIF
962
963	!-------------------------------------------------------------------------------
964	! Calcul de la fraction du thermique et des ecart-types des distributions
965	!-------------------------------------------------------------------------------
966	do ind1=1,ngrid
967
968	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
969
970	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
971
972
973	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
974	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
975	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
976	qsatbef=MIN(0.5,qsatbef)
977	zcor=1./(1.-retv*qsatbef)
978	qsatbef=qsatbef*zcor
979	zqsatenv(ind1,ind2)=qsatbef
980
981
982	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
983	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
984	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
985
986	!zqenv = qt environnement
987	!zqsatenv = qsat environnement
988	!zthl = Tl environnement
989
990
991	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
992	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
993	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
994	qsatbef=MIN(0.5,qsatbef)
995	zcor=1./(1.-retv*qsatbef)
996	qsatbef=qsatbef*zcor
997	zqsatth(ind1,ind2)=qsatbef
998
999	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
1000	ath=1./(1.+(alth*Lv/cppd)) !al, p84
1001	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
1002
1003
1004	!zqta = qt thermals
1005	!zqsatth = qsat thermals
1006	!ztla = Tl thermals
1007
1008	!------------------------------------------------------------------------------
1009	! s standard deviation
1010	!------------------------------------------------------------------------------
1011
1012	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
1013	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
1014	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
1015	IF (cloudth_ratqsmin>0.) THEN
1016	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1017	ELSE
1018	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1019	ENDIF
1020	sigma1s = sigma1s_fraca + sigma1s_ratqs
1021	sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
1022	! tests
1023	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
1024	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
1025	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
1026	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
1027	! if (paprs(ind1,ind2).gt.90000) then
1028	! ratqs(ind1,ind2)=0.002
1029	! else
1030	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1031	! endif
1032	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1033	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1034	! sigma1s=ratqs(ind1,ind2)*po(ind1)
1035	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
1036
1037	IF (iflag_cloudth_vert == 1) THEN
1038	!-------------------------------------------------------------------------------
1039	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1040	!-------------------------------------------------------------------------------
1041
1042	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1043	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1044
1045	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1046	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1047	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1048	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1049	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
1050	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
1051
1052	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
1053	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
1054	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1055
1056	! Environment
1057	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
1058	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
1059	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
1060	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
1061
1062	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1063
1064	! Thermal
1065	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
1066	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
1067	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
1068	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
1069	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1070	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1071
1072	ELSE IF (iflag_cloudth_vert >= 3) THEN
1073	IF (iflag_cloudth_vert < 5) THEN
1074	!-------------------------------------------------------------------------------
1075	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1076	!-------------------------------------------------------------------------------
1077	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
1078	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
1079	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1080	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1081	IF (iflag_cloudth_vert == 3) THEN
1082	deltasenv=aenvvert_alphasigma1s
1083	deltasth=athvert_alpha_thsigma2s
1084	ELSE IF (iflag_cloudth_vert == 4) THEN
1085	IF (iflag_cloudth_vert_noratqs == 1) THEN
1086	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1087	deltasth=vert_alpha_th*sigma2s
1088	ELSE
1089	deltasenv=vert_alpha*sigma1s
1090	deltasth=vert_alpha_th*sigma2s
1091	ENDIF
1092	ENDIF
1093
1094	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1095	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1096	exp_xenv1 = exp(-1.xenv1*2)
1097	exp_xenv2 = exp(-1.xenv2*2)
1098	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1099	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1100	exp_xth1 = exp(-1.xth1*2)
1101	exp_xth2 = exp(-1.xth2*2)
1102
1103	!CF_surfacique
1104	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1105	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1106	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1107
1108
1109	!CF_volumique & eau condense
1110	!environnement
1111	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1112	IntJ_CF=0.5(1.-1.erf(xenv2))
1113	if (deltasenv .lt. 1.e-10) then
1114	qlenv(ind1,ind2)=IntJ
1115	cenv_vol(ind1,ind2)=IntJ_CF
1116	else
1117	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1118	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1119	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1120	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1121	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1122	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1123	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1124	endif
1125
1126	!thermique
1127	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1128	IntJ_CF=0.5(1.-1.erf(xth2))
1129	if (deltasth .lt. 1.e-10) then
1130	qlth(ind1,ind2)=IntJ
1131	cth_vol(ind1,ind2)=IntJ_CF
1132	else
1133	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1134	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1135	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1136	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1137	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1138	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1139	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1140	endif
1141
1142	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1143	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1144
1145
1146
1147
1148
1149
1150	ELSE IF (iflag_cloudth_vert == 5) THEN
1151	sigma1s=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5)+ratqs(ind1,ind2)*po(ind1) !Environment
1152	sigma2s=(0.03218+0.000092655dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)0.5)(((sth-senv)2)0.5)+0.002*zqta(ind1,ind2) !Thermals
1153	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1154	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1155	xth=sth/(sqrt(2.)*sigma2s)
1156	xenv=senv/(sqrt(2.)*sigma1s)
1157
1158	!Volumique
1159	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1160	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1161	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1162	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1163
1164	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1165	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1166	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1167
1168	!Surfacique
1169	!Neggers
1170	!beta=0.0044
1171	!inverse_rho=1.+beta*dz(ind1,ind2)
1172	!print *,'jeanjean : beta=',beta
1173	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1174	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1175	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1176
1177	!Brooks
1178	a_Brooks=0.6694
1179	b_Brooks=0.1882
1180	A_Maj_Brooks=0.1635 !-- sans shear
1181	!A_Maj_Brooks=0.17 !-- ARM LES
1182	!A_Maj_Brooks=0.18 !-- RICO LES
1183	!A_Maj_Brooks=0.19 !-- BOMEX LES
1184	Dx_Brooks=200000.
1185	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1186	!print *,'jeanjean_f=',f_Brooks
1187
1188	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1189	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1190	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1191	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1192
1193
1194
1195
1196
1197	ENDIF ! of if (iflag_cloudth_vert<5)
1198	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1199
1200	if (ctot(ind1,ind2).lt.1.e-10) then
1201	ctot(ind1,ind2)=0.
1202	ctot_vol(ind1,ind2)=0.
1203	qcloud(ind1)=zqsatenv(ind1,ind2)
1204
1205	else
1206
1207	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1208	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1209	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1210
1211	endif
1212
1213	else ! gaussienne environnement seule
1214
1215	zqenv(ind1)=po(ind1)
1216	Tbef=t(ind1,ind2)
1217	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1218	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1219	qsatbef=MIN(0.5,qsatbef)
1220	zcor=1./(1.-retv*qsatbef)
1221	qsatbef=qsatbef*zcor
1222	zqsatenv(ind1,ind2)=qsatbef
1223
1224
1225	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1226	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1227	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1228	aenv=1./(1.+(alenv*Lv/cppd))
1229	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1230
1231
1232	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1233
1234	sqrt2pi=sqrt(2.*pi)
1235	xenv=senv/(sqrt(2.)*sigma1s)
1236	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1237	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1238
1239	if (ctot(ind1,ind2).lt.1.e-3) then
1240	ctot(ind1,ind2)=0.
1241	qcloud(ind1)=zqsatenv(ind1,ind2)
1242
1243	else
1244
1245	ctot(ind1,ind2)=ctot(ind1,ind2)
1246	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1247
1248	endif
1249
1250
1251
1252
1253	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1254	! Outputs used to check the PDFs
1255	cloudth_senv(ind1,ind2) = senv
1256	cloudth_sth(ind1,ind2) = sth
1257	cloudth_sigmaenv(ind1,ind2) = sigma1s
1258	cloudth_sigmath(ind1,ind2) = sigma2s
1259
1260	enddo ! from the loop on ngrid l.333
1261	return
1262	! end
1263	END SUBROUTINE cloudth_vert_v3
1264	!
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1277	& ztv,po,zqta,fraca, &
1278	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1279	& ratqs,zqs,T)
1280
1281
1282	USE ioipsl_getin_p_mod, ONLY : getin_p
1283	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1284	& cloudth_sigmath,cloudth_sigmaenv
1285
1286	IMPLICIT NONE
1287
1288	#include "YOMCST.h"
1289	#include "YOETHF.h"
1290	#include "FCTTRE.h"
1291	#include "thermcell.h"
1292	#include "nuage.h"
1293
1294
1295	!Domain variables
1296	INTEGER ngrid !indice Max lat-lon
1297	INTEGER klev !indice Max alt
1298	INTEGER ind1 !indice in [1:ngrid]
1299	INTEGER ind2 !indice in [1:klev]
1300	!thermal plume fraction
1301	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1302	!temperatures
1303	REAL T(ngrid,klev) !temperature
1304	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1305	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1306	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1307	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1308	!pressure
1309	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1310	REAL pplay(ngrid,klev) !pressure at the middle of the level
1311	!humidity
1312	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1313	REAL po(ngrid) !total water (qt)
1314	REAL zqenv(ngrid) !total water in the environment (qt_env)
1315	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1316	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1317	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1318	REAL qlth(ngrid,klev) !condensed water in the thermals
1319	REAL qlenv(ngrid,klev) !condensed water in the environment
1320	REAL qltot(ngrid,klev) !condensed water in the gridbox
1321	!cloud fractions
1322	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1323	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1324	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1325	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1326	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1327	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1328	!PDF of saturation deficit variables
1329	REAL rdd,cppd,Lv
1330	REAL Tbef,zdelta,qsatbef,zcor
1331	REAL alth,alenv,ath,aenv
1332	REAL sth,senv !saturation deficits in the thermals and environment
1333	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1334	!cloud fraction variables
1335	REAL xth,xenv
1336	REAL inverse_rho,beta !Neggers et al. (2011) method
1337	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1338	!Incloud total water variables
1339	REAL zqs(ngrid) !q_sat
1340	REAL qcloud(ngrid) !eau totale dans le nuage
1341	!Some arithmetic variables
1342	REAL erf,pi,sqrt2,sqrt2pi
1343	!Depth of the layer
1344	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1345	REAL rhodz(ngrid,klev)
1346	REAL zrho(ngrid,klev)
1347	DO ind1 = 1, ngrid
1348	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1349	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1350	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1351	END DO
1352
1353	!------------------------------------------------------------------------------
1354	! Initialization
1355	!------------------------------------------------------------------------------
1356	qlth(:,:)=0.
1357	qlenv(:,:)=0.
1358	qltot(:,:)=0.
1359	cth_vol(:,:)=0.
1360	cenv_vol(:,:)=0.
1361	ctot_vol(:,:)=0.
1362	cth_surf(:,:)=0.
1363	cenv_surf(:,:)=0.
1364	ctot_surf(:,:)=0.
1365	qcloud(:)=0.
1366	rdd=287.04
1367	cppd=1005.7
1368	pi=3.14159
1369	Lv=2.5e6
1370	sqrt2=sqrt(2.)
1371	sqrt2pi=sqrt(2.*pi)
1372
1373
1374	DO ind1=1,ngrid
1375	!-------------------------------------------------------------------------------
1376	!Both thermal and environment in the gridbox
1377	!-------------------------------------------------------------------------------
1378	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1379	!--------------------------------------------
1380	!calcul de qsat_env
1381	!--------------------------------------------
1382	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1383	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1384	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1385	qsatbef=MIN(0.5,qsatbef)
1386	zcor=1./(1.-retv*qsatbef)
1387	qsatbef=qsatbef*zcor
1388	zqsatenv(ind1,ind2)=qsatbef
1389	!--------------------------------------------
1390	!calcul de s_env
1391	!--------------------------------------------
1392	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1393	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1394	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1395	!--------------------------------------------
1396	!calcul de qsat_th
1397	!--------------------------------------------
1398	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1399	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1400	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1401	qsatbef=MIN(0.5,qsatbef)
1402	zcor=1./(1.-retv*qsatbef)
1403	qsatbef=qsatbef*zcor
1404	zqsatth(ind1,ind2)=qsatbef
1405	!--------------------------------------------
1406	!calcul de s_th
1407	!--------------------------------------------
1408	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1409	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1410	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1411	!--------------------------------------------
1412	!calcul standard deviations bi-Gaussian PDF
1413	!--------------------------------------------
1414	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)
1415	sigma_env=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5)+ratqs(ind1,ind2)*po(ind1)
1416	xth=sth/(sqrt2*sigma_th)
1417	xenv=senv/(sqrt2*sigma_env)
1418	!--------------------------------------------
1419	!Cloud fraction by volume CF_vol
1420	!--------------------------------------------
1421	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1422	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1423	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1424	!--------------------------------------------
1425	!Condensed water qc
1426	!--------------------------------------------
1427	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1428	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1429	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1430	!--------------------------------------------
1431	!Cloud fraction by surface CF_surf
1432	!--------------------------------------------
1433	!Method Neggers et al. (2011) : ok for cumulus clouds only
1434	!beta=0.0044 (Jouhaud et al.2018)
1435	!inverse_rho=1.+beta*dz(ind1,ind2)
1436	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1437	!Method Brooks et al. (2005) : ok for all types of clouds
1438	a_Brooks=0.6694
1439	b_Brooks=0.1882
1440	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1441	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1442	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1443	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1444	!--------------------------------------------
1445	!Incloud Condensed water qcloud
1446	!--------------------------------------------
1447	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1448	ctot_vol(ind1,ind2)=0.
1449	ctot_surf(ind1,ind2)=0.
1450	qcloud(ind1)=zqsatenv(ind1,ind2)
1451	else
1452	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1453	endif
1454
1455
1456
1457	!-------------------------------------------------------------------------------
1458	!Environment only in the gridbox
1459	!-------------------------------------------------------------------------------
1460	ELSE
1461	!--------------------------------------------
1462	!calcul de qsat_env
1463	!--------------------------------------------
1464	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1465	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1466	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1467	qsatbef=MIN(0.5,qsatbef)
1468	zcor=1./(1.-retv*qsatbef)
1469	qsatbef=qsatbef*zcor
1470	zqsatenv(ind1,ind2)=qsatbef
1471	!--------------------------------------------
1472	!calcul de s_env
1473	!--------------------------------------------
1474	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1475	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1476	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1477	!--------------------------------------------
1478	!calcul standard deviations Gaussian PDF
1479	!--------------------------------------------
1480	zqenv(ind1)=po(ind1)
1481	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1482	xenv=senv/(sqrt2*sigma_env)
1483	!--------------------------------------------
1484	!Cloud fraction by volume CF_vol
1485	!--------------------------------------------
1486	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1487	!--------------------------------------------
1488	!Condensed water qc
1489	!--------------------------------------------
1490	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1491	!--------------------------------------------
1492	!Cloud fraction by surface CF_surf
1493	!--------------------------------------------
1494	!Method Neggers et al. (2011) : ok for cumulus clouds only
1495	!beta=0.0044 (Jouhaud et al.2018)
1496	!inverse_rho=1.+beta*dz(ind1,ind2)
1497	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1498	!Method Brooks et al. (2005) : ok for all types of clouds
1499	a_Brooks=0.6694
1500	b_Brooks=0.1882
1501	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1502	Dx_Brooks=200000.
1503	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1504	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1505	!--------------------------------------------
1506	!Incloud Condensed water qcloud
1507	!--------------------------------------------
1508	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1509	ctot_vol(ind1,ind2)=0.
1510	ctot_surf(ind1,ind2)=0.
1511	qcloud(ind1)=zqsatenv(ind1,ind2)
1512	else
1513	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1514	endif
1515
1516
1517	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1518
1519	! Outputs used to check the PDFs
1520	cloudth_senv(ind1,ind2) = senv
1521	cloudth_sth(ind1,ind2) = sth
1522	cloudth_sigmaenv(ind1,ind2) = sigma_env
1523	cloudth_sigmath(ind1,ind2) = sigma_th
1524
1525	END DO ! From the loop on ngrid
1526	return
1527
1528	END SUBROUTINE cloudth_v6
1529	END MODULE cloudth_mod
1530
1531
1532
1533

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

Download in other formats:

Original Format