Context Navigation

cloudth_mod.F90 @ 3627

Last change on this file since 3627 was 3605, checked in by lguez, 4 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
File size: 57.1 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 :: sigma2s_factor=0.09
888	REAL,SAVE :: sigma2s_power=0.5
889	REAL,SAVE :: cloudth_ratqsmin=-1.
890	!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
891	INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
892	!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
893
894	LOGICAL, SAVE :: firstcall = .TRUE.
895	!$OMP THREADPRIVATE(firstcall)
896
897	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
898	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
899	REAL zqs(ngrid), qcloud(ngrid)
900	REAL erf
901
902	REAL rhodz(ngrid,klev)
903	REAL zrho(ngrid,klev)
904	REAL dz(ngrid,klev)
905
906	DO ind1 = 1, ngrid
907	!Layer calculation
908	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
909	zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
910	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
911	END DO
912
913
914	!------------------------------------------------------------------------------
915	! Initialize
916	!------------------------------------------------------------------------------
917	sigma1(:,:)=0.
918	sigma2(:,:)=0.
919	qlth(:,:)=0.
920	qlenv(:,:)=0.
921	qltot(:,:)=0.
922	rneb(:,:)=0.
923	qcloud(:)=0.
924	cth(:,:)=0.
925	cenv(:,:)=0.
926	ctot(:,:)=0.
927	cth_vol(:,:)=0.
928	cenv_vol(:,:)=0.
929	ctot_vol(:,:)=0.
930	qsatmmussig1=0.
931	qsatmmussig2=0.
932	rdd=287.04
933	cppd=1005.7
934	pi=3.14159
935	Lv=2.5e6
936	sqrt2pi=sqrt(2.*pi)
937	sqrt2=sqrt(2.)
938	sqrtpi=sqrt(pi)
939
940	IF (firstcall) THEN
941	vert_alpha=0.5
942	CALL getin_p('cloudth_vert_alpha',vert_alpha)
943	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
944	! The factor used for the thermal is equal to that of the environment
945	! if nothing is explicitly specified in the def file
946	vert_alpha_th=vert_alpha
947	CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
948	WRITE(,) 'cloudth_vert_alpha_th = ', vert_alpha_th
949	! Factor used in the calculation of sigma1s
950	CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
951	WRITE(,) 'cloudth_sigma1s_factor = ', sigma1s_factor
952	! Power used in the calculation of sigma1s
953	CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
954	WRITE(,) 'cloudth_sigma1s_power = ', sigma1s_power
955	! Factor used in the calculation of sigma2s
956	CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
957	WRITE(,) 'cloudth_sigma2s_factor = ', sigma2s_factor
958	! Power used in the calculation of sigma2s
959	CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
960	WRITE(,) 'cloudth_sigma2s_power = ', sigma2s_power
961	! Minimum value for the environmental air subgrid water distrib
962	CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
963	WRITE(,) 'cloudth_ratqsmin = ', cloudth_ratqsmin
964	! Remove the dependency to ratqs from the variance of the vertical PDF
965	CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
966	WRITE(,) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
967
968	firstcall=.FALSE.
969	ENDIF
970
971	!-------------------------------------------------------------------------------
972	! Calcul de la fraction du thermique et des ecart-types des distributions
973	!-------------------------------------------------------------------------------
974	do ind1=1,ngrid
975
976	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
977
978	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
979
980
981	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
982	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
983	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
984	qsatbef=MIN(0.5,qsatbef)
985	zcor=1./(1.-retv*qsatbef)
986	qsatbef=qsatbef*zcor
987	zqsatenv(ind1,ind2)=qsatbef
988
989
990	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
991	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
992	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
993
994	!zqenv = qt environnement
995	!zqsatenv = qsat environnement
996	!zthl = Tl environnement
997
998
999	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1000	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1001	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1002	qsatbef=MIN(0.5,qsatbef)
1003	zcor=1./(1.-retv*qsatbef)
1004	qsatbef=qsatbef*zcor
1005	zqsatth(ind1,ind2)=qsatbef
1006
1007	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
1008	ath=1./(1.+(alth*Lv/cppd)) !al, p84
1009	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
1010
1011
1012	!zqta = qt thermals
1013	!zqsatth = qsat thermals
1014	!ztla = Tl thermals
1015
1016	!------------------------------------------------------------------------------
1017	! s standard deviation
1018	!------------------------------------------------------------------------------
1019
1020	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
1021	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
1022	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
1023	IF (cloudth_ratqsmin>0.) THEN
1024	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1025	ELSE
1026	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1027	ENDIF
1028	sigma1s = sigma1s_fraca + sigma1s_ratqs
1029	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
1030	! tests
1031	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
1032	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
1033	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
1034	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
1035	! if (paprs(ind1,ind2).gt.90000) then
1036	! ratqs(ind1,ind2)=0.002
1037	! else
1038	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1039	! endif
1040	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1041	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1042	! sigma1s=ratqs(ind1,ind2)*po(ind1)
1043	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
1044
1045	IF (iflag_cloudth_vert == 1) THEN
1046	!-------------------------------------------------------------------------------
1047	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1048	!-------------------------------------------------------------------------------
1049
1050	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1051	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1052
1053	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1054	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1055	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1056	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1057	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
1058	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
1059
1060	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
1061	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
1062	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1063
1064	! Environment
1065	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
1066	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
1067	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
1068	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
1069
1070	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1071
1072	! Thermal
1073	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
1074	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
1075	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
1076	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
1077	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1078	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1079
1080	ELSE IF (iflag_cloudth_vert >= 3) THEN
1081	IF (iflag_cloudth_vert < 5) THEN
1082	!-------------------------------------------------------------------------------
1083	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1084	!-------------------------------------------------------------------------------
1085	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
1086	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
1087	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
1088	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
1089	IF (iflag_cloudth_vert == 3) THEN
1090	deltasenv=aenvvert_alphasigma1s
1091	deltasth=athvert_alpha_thsigma2s
1092	ELSE IF (iflag_cloudth_vert == 4) THEN
1093	IF (iflag_cloudth_vert_noratqs == 1) THEN
1094	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1095	deltasth=vert_alpha_th*sigma2s
1096	ELSE
1097	deltasenv=vert_alpha*sigma1s
1098	deltasth=vert_alpha_th*sigma2s
1099	ENDIF
1100	ENDIF
1101
1102	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1103	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1104	exp_xenv1 = exp(-1.xenv1*2)
1105	exp_xenv2 = exp(-1.xenv2*2)
1106	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1107	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1108	exp_xth1 = exp(-1.xth1*2)
1109	exp_xth2 = exp(-1.xth2*2)
1110
1111	!CF_surfacique
1112	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
1113	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
1114	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1115
1116
1117	!CF_volumique & eau condense
1118	!environnement
1119	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1120	IntJ_CF=0.5(1.-1.erf(xenv2))
1121	if (deltasenv .lt. 1.e-10) then
1122	qlenv(ind1,ind2)=IntJ
1123	cenv_vol(ind1,ind2)=IntJ_CF
1124	else
1125	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
1126	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
1127	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
1128	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
1129	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
1130	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1131	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1132	endif
1133
1134	!thermique
1135	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1136	IntJ_CF=0.5(1.-1.erf(xth2))
1137	if (deltasth .lt. 1.e-10) then
1138	qlth(ind1,ind2)=IntJ
1139	cth_vol(ind1,ind2)=IntJ_CF
1140	else
1141	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
1142	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
1143	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
1144	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
1145	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
1146	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1147	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1148	endif
1149
1150	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1151	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1152
1153	ELSE IF (iflag_cloudth_vert == 5) THEN
1154	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
1155	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
1156	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1157	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1158	xth=sth/(sqrt(2.)*sigma2s)
1159	xenv=senv/(sqrt(2.)*sigma1s)
1160
1161	!Volumique
1162	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1163	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1164	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1165	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1166
1167	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1168	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1169	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1170
1171	!Surfacique
1172	!Neggers
1173	!beta=0.0044
1174	!inverse_rho=1.+beta*dz(ind1,ind2)
1175	!print *,'jeanjean : beta=',beta
1176	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1177	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1178	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1179
1180	!Brooks
1181	a_Brooks=0.6694
1182	b_Brooks=0.1882
1183	A_Maj_Brooks=0.1635 !-- sans shear
1184	!A_Maj_Brooks=0.17 !-- ARM LES
1185	!A_Maj_Brooks=0.18 !-- RICO LES
1186	!A_Maj_Brooks=0.19 !-- BOMEX LES
1187	Dx_Brooks=200000.
1188	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1189	!print *,'jeanjean_f=',f_Brooks
1190
1191	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1192	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1193	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1194	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1195
1196
1197
1198
1199
1200	ENDIF ! of if (iflag_cloudth_vert<5)
1201	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1202
1203	! if (ctot(ind1,ind2).lt.1.e-10) then
1204	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1205	ctot(ind1,ind2)=0.
1206	ctot_vol(ind1,ind2)=0.
1207	qcloud(ind1)=zqsatenv(ind1,ind2)
1208
1209	else
1210
1211	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1212	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1213	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1214
1215	endif
1216
1217	else ! gaussienne environnement seule
1218
1219	zqenv(ind1)=po(ind1)
1220	Tbef=t(ind1,ind2)
1221	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1222	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1223	qsatbef=MIN(0.5,qsatbef)
1224	zcor=1./(1.-retv*qsatbef)
1225	qsatbef=qsatbef*zcor
1226	zqsatenv(ind1,ind2)=qsatbef
1227
1228
1229	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1230	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1231	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1232	aenv=1./(1.+(alenv*Lv/cppd))
1233	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1234	sth=0.
1235
1236
1237	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1238	sigma2s=0.
1239
1240	sqrt2pi=sqrt(2.*pi)
1241	xenv=senv/(sqrt(2.)*sigma1s)
1242	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1243	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1244	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1245
1246	if (ctot(ind1,ind2).lt.1.e-3) then
1247	ctot(ind1,ind2)=0.
1248	qcloud(ind1)=zqsatenv(ind1,ind2)
1249
1250	else
1251
1252	! ctot(ind1,ind2)=ctot(ind1,ind2)
1253	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1254
1255	endif
1256
1257
1258
1259
1260	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1261	! Outputs used to check the PDFs
1262	cloudth_senv(ind1,ind2) = senv
1263	cloudth_sth(ind1,ind2) = sth
1264	cloudth_sigmaenv(ind1,ind2) = sigma1s
1265	cloudth_sigmath(ind1,ind2) = sigma2s
1266
1267	enddo ! from the loop on ngrid l.333
1268	return
1269	! end
1270	END SUBROUTINE cloudth_vert_v3
1271	!
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1284	& ztv,po,zqta,fraca, &
1285	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1286	& ratqs,zqs,T)
1287
1288
1289	USE ioipsl_getin_p_mod, ONLY : getin_p
1290	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1291	& cloudth_sigmath,cloudth_sigmaenv
1292
1293	IMPLICIT NONE
1294
1295	#include "YOMCST.h"
1296	#include "YOETHF.h"
1297	#include "FCTTRE.h"
1298	#include "thermcell.h"
1299	#include "nuage.h"
1300
1301
1302	!Domain variables
1303	INTEGER ngrid !indice Max lat-lon
1304	INTEGER klev !indice Max alt
1305	INTEGER ind1 !indice in [1:ngrid]
1306	INTEGER ind2 !indice in [1:klev]
1307	!thermal plume fraction
1308	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1309	!temperatures
1310	REAL T(ngrid,klev) !temperature
1311	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1312	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1313	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1314	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1315	!pressure
1316	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1317	REAL pplay(ngrid,klev) !pressure at the middle of the level
1318	!humidity
1319	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1320	REAL po(ngrid) !total water (qt)
1321	REAL zqenv(ngrid) !total water in the environment (qt_env)
1322	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1323	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1324	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1325	REAL qlth(ngrid,klev) !condensed water in the thermals
1326	REAL qlenv(ngrid,klev) !condensed water in the environment
1327	REAL qltot(ngrid,klev) !condensed water in the gridbox
1328	!cloud fractions
1329	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1330	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1331	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1332	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1333	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1334	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1335	!PDF of saturation deficit variables
1336	REAL rdd,cppd,Lv
1337	REAL Tbef,zdelta,qsatbef,zcor
1338	REAL alth,alenv,ath,aenv
1339	REAL sth,senv !saturation deficits in the thermals and environment
1340	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1341	!cloud fraction variables
1342	REAL xth,xenv
1343	REAL inverse_rho,beta !Neggers et al. (2011) method
1344	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1345	!Incloud total water variables
1346	REAL zqs(ngrid) !q_sat
1347	REAL qcloud(ngrid) !eau totale dans le nuage
1348	!Some arithmetic variables
1349	REAL erf,pi,sqrt2,sqrt2pi
1350	!Depth of the layer
1351	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1352	REAL rhodz(ngrid,klev)
1353	REAL zrho(ngrid,klev)
1354	DO ind1 = 1, ngrid
1355	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1356	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1357	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1358	END DO
1359
1360	!------------------------------------------------------------------------------
1361	! Initialization
1362	!------------------------------------------------------------------------------
1363	qlth(:,:)=0.
1364	qlenv(:,:)=0.
1365	qltot(:,:)=0.
1366	cth_vol(:,:)=0.
1367	cenv_vol(:,:)=0.
1368	ctot_vol(:,:)=0.
1369	cth_surf(:,:)=0.
1370	cenv_surf(:,:)=0.
1371	ctot_surf(:,:)=0.
1372	qcloud(:)=0.
1373	rdd=287.04
1374	cppd=1005.7
1375	pi=3.14159
1376	Lv=2.5e6
1377	sqrt2=sqrt(2.)
1378	sqrt2pi=sqrt(2.*pi)
1379
1380
1381	DO ind1=1,ngrid
1382	!-------------------------------------------------------------------------------
1383	!Both thermal and environment in the gridbox
1384	!-------------------------------------------------------------------------------
1385	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1386	!--------------------------------------------
1387	!calcul de qsat_env
1388	!--------------------------------------------
1389	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1390	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1391	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1392	qsatbef=MIN(0.5,qsatbef)
1393	zcor=1./(1.-retv*qsatbef)
1394	qsatbef=qsatbef*zcor
1395	zqsatenv(ind1,ind2)=qsatbef
1396	!--------------------------------------------
1397	!calcul de s_env
1398	!--------------------------------------------
1399	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1400	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1401	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1402	!--------------------------------------------
1403	!calcul de qsat_th
1404	!--------------------------------------------
1405	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1406	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1407	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1408	qsatbef=MIN(0.5,qsatbef)
1409	zcor=1./(1.-retv*qsatbef)
1410	qsatbef=qsatbef*zcor
1411	zqsatth(ind1,ind2)=qsatbef
1412	!--------------------------------------------
1413	!calcul de s_th
1414	!--------------------------------------------
1415	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1416	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1417	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1418	!--------------------------------------------
1419	!calcul standard deviations bi-Gaussian PDF
1420	!--------------------------------------------
1421	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)
1422	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)
1423	xth=sth/(sqrt2*sigma_th)
1424	xenv=senv/(sqrt2*sigma_env)
1425	!--------------------------------------------
1426	!Cloud fraction by volume CF_vol
1427	!--------------------------------------------
1428	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1429	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1430	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1431	!--------------------------------------------
1432	!Condensed water qc
1433	!--------------------------------------------
1434	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1435	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1436	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1437	!--------------------------------------------
1438	!Cloud fraction by surface CF_surf
1439	!--------------------------------------------
1440	!Method Neggers et al. (2011) : ok for cumulus clouds only
1441	!beta=0.0044 (Jouhaud et al.2018)
1442	!inverse_rho=1.+beta*dz(ind1,ind2)
1443	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1444	!Method Brooks et al. (2005) : ok for all types of clouds
1445	a_Brooks=0.6694
1446	b_Brooks=0.1882
1447	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1448	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1449	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1450	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1451	!--------------------------------------------
1452	!Incloud Condensed water qcloud
1453	!--------------------------------------------
1454	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1455	ctot_vol(ind1,ind2)=0.
1456	ctot_surf(ind1,ind2)=0.
1457	qcloud(ind1)=zqsatenv(ind1,ind2)
1458	else
1459	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1460	endif
1461
1462
1463
1464	!-------------------------------------------------------------------------------
1465	!Environment only in the gridbox
1466	!-------------------------------------------------------------------------------
1467	ELSE
1468	!--------------------------------------------
1469	!calcul de qsat_env
1470	!--------------------------------------------
1471	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1472	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1473	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1474	qsatbef=MIN(0.5,qsatbef)
1475	zcor=1./(1.-retv*qsatbef)
1476	qsatbef=qsatbef*zcor
1477	zqsatenv(ind1,ind2)=qsatbef
1478	!--------------------------------------------
1479	!calcul de s_env
1480	!--------------------------------------------
1481	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1482	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1483	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1484	!--------------------------------------------
1485	!calcul standard deviations Gaussian PDF
1486	!--------------------------------------------
1487	zqenv(ind1)=po(ind1)
1488	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1489	xenv=senv/(sqrt2*sigma_env)
1490	!--------------------------------------------
1491	!Cloud fraction by volume CF_vol
1492	!--------------------------------------------
1493	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1494	!--------------------------------------------
1495	!Condensed water qc
1496	!--------------------------------------------
1497	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1498	!--------------------------------------------
1499	!Cloud fraction by surface CF_surf
1500	!--------------------------------------------
1501	!Method Neggers et al. (2011) : ok for cumulus clouds only
1502	!beta=0.0044 (Jouhaud et al.2018)
1503	!inverse_rho=1.+beta*dz(ind1,ind2)
1504	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1505	!Method Brooks et al. (2005) : ok for all types of clouds
1506	a_Brooks=0.6694
1507	b_Brooks=0.1882
1508	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1509	Dx_Brooks=200000.
1510	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1511	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1512	!--------------------------------------------
1513	!Incloud Condensed water qcloud
1514	!--------------------------------------------
1515	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1516	ctot_vol(ind1,ind2)=0.
1517	ctot_surf(ind1,ind2)=0.
1518	qcloud(ind1)=zqsatenv(ind1,ind2)
1519	else
1520	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1521	endif
1522
1523
1524	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1525
1526	! Outputs used to check the PDFs
1527	cloudth_senv(ind1,ind2) = senv
1528	cloudth_sth(ind1,ind2) = sth
1529	cloudth_sigmaenv(ind1,ind2) = sigma_env
1530	cloudth_sigmath(ind1,ind2) = sigma_th
1531
1532	END DO ! From the loop on ngrid
1533	return
1534
1535	END SUBROUTINE cloudth_v6
1536	END MODULE cloudth_mod
1537
1538
1539
1540

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

Download in other formats:

Original Format