Context Navigation

cloudth_mod.F90 @ 3495

Last change on this file since 3495 was 3495, checked in by idelkadi, 6 years ago
Correction de l'erreur introduite dans la version svn3493
File size: 56.6 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	ELSE IF (iflag_cloudth_vert == 5) THEN
1146	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
1147	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
1148	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
1149	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
1150	xth=sth/(sqrt(2.)*sigma2s)
1151	xenv=senv/(sqrt(2.)*sigma1s)
1152
1153	!Volumique
1154	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1155	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1156	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1157	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1158
1159	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
1160	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
1161	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1162
1163	!Surfacique
1164	!Neggers
1165	!beta=0.0044
1166	!inverse_rho=1.+beta*dz(ind1,ind2)
1167	!print *,'jeanjean : beta=',beta
1168	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1169	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1170	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
1171
1172	!Brooks
1173	a_Brooks=0.6694
1174	b_Brooks=0.1882
1175	A_Maj_Brooks=0.1635 !-- sans shear
1176	!A_Maj_Brooks=0.17 !-- ARM LES
1177	!A_Maj_Brooks=0.18 !-- RICO LES
1178	!A_Maj_Brooks=0.19 !-- BOMEX LES
1179	Dx_Brooks=200000.
1180	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1181	!print *,'jeanjean_f=',f_Brooks
1182
1183	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1184	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1185	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1186	!print *,'JJ_ctot_1',ctot(ind1,ind2)
1187
1188
1189
1190
1191
1192	ENDIF ! of if (iflag_cloudth_vert<5)
1193	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1194
1195	! if (ctot(ind1,ind2).lt.1.e-10) then
1196	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1197	ctot(ind1,ind2)=0.
1198	ctot_vol(ind1,ind2)=0.
1199	qcloud(ind1)=zqsatenv(ind1,ind2)
1200
1201	else
1202
1203	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1204	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1205	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1206
1207	endif
1208
1209	else ! gaussienne environnement seule
1210
1211	zqenv(ind1)=po(ind1)
1212	Tbef=t(ind1,ind2)
1213	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1214	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1215	qsatbef=MIN(0.5,qsatbef)
1216	zcor=1./(1.-retv*qsatbef)
1217	qsatbef=qsatbef*zcor
1218	zqsatenv(ind1,ind2)=qsatbef
1219
1220
1221	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1222	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
1223	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
1224	aenv=1./(1.+(alenv*Lv/cppd))
1225	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1226	sth=0.
1227
1228
1229	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1230	sigma2s=0.
1231
1232	sqrt2pi=sqrt(2.*pi)
1233	xenv=senv/(sqrt(2.)*sigma1s)
1234	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
1235	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1236	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
1237
1238	if (ctot(ind1,ind2).lt.1.e-3) then
1239	ctot(ind1,ind2)=0.
1240	qcloud(ind1)=zqsatenv(ind1,ind2)
1241
1242	else
1243
1244	! ctot(ind1,ind2)=ctot(ind1,ind2)
1245	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1246
1247	endif
1248
1249
1250
1251
1252	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1253	! Outputs used to check the PDFs
1254	cloudth_senv(ind1,ind2) = senv
1255	cloudth_sth(ind1,ind2) = sth
1256	cloudth_sigmaenv(ind1,ind2) = sigma1s
1257	cloudth_sigmath(ind1,ind2) = sigma2s
1258
1259	enddo ! from the loop on ngrid l.333
1260	return
1261	! end
1262	END SUBROUTINE cloudth_vert_v3
1263	!
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1276	& ztv,po,zqta,fraca, &
1277	& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1278	& ratqs,zqs,T)
1279
1280
1281	USE ioipsl_getin_p_mod, ONLY : getin_p
1282	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1283	& cloudth_sigmath,cloudth_sigmaenv
1284
1285	IMPLICIT NONE
1286
1287	#include "YOMCST.h"
1288	#include "YOETHF.h"
1289	#include "FCTTRE.h"
1290	#include "thermcell.h"
1291	#include "nuage.h"
1292
1293
1294	!Domain variables
1295	INTEGER ngrid !indice Max lat-lon
1296	INTEGER klev !indice Max alt
1297	INTEGER ind1 !indice in [1:ngrid]
1298	INTEGER ind2 !indice in [1:klev]
1299	!thermal plume fraction
1300	REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
1301	!temperatures
1302	REAL T(ngrid,klev) !temperature
1303	REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
1304	REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
1305	REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
1306	REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
1307	!pressure
1308	REAL paprs(ngrid,klev+1) !pressure at the interface of levels
1309	REAL pplay(ngrid,klev) !pressure at the middle of the level
1310	!humidity
1311	REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
1312	REAL po(ngrid) !total water (qt)
1313	REAL zqenv(ngrid) !total water in the environment (qt_env)
1314	REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
1315	REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
1316	REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
1317	REAL qlth(ngrid,klev) !condensed water in the thermals
1318	REAL qlenv(ngrid,klev) !condensed water in the environment
1319	REAL qltot(ngrid,klev) !condensed water in the gridbox
1320	!cloud fractions
1321	REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
1322	REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
1323	REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
1324	REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
1325	REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
1326	REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
1327	!PDF of saturation deficit variables
1328	REAL rdd,cppd,Lv
1329	REAL Tbef,zdelta,qsatbef,zcor
1330	REAL alth,alenv,ath,aenv
1331	REAL sth,senv !saturation deficits in the thermals and environment
1332	REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
1333	!cloud fraction variables
1334	REAL xth,xenv
1335	REAL inverse_rho,beta !Neggers et al. (2011) method
1336	REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
1337	!Incloud total water variables
1338	REAL zqs(ngrid) !q_sat
1339	REAL qcloud(ngrid) !eau totale dans le nuage
1340	!Some arithmetic variables
1341	REAL erf,pi,sqrt2,sqrt2pi
1342	!Depth of the layer
1343	REAL dz(ngrid,klev) !epaisseur de la couche en metre
1344	REAL rhodz(ngrid,klev)
1345	REAL zrho(ngrid,klev)
1346	DO ind1 = 1, ngrid
1347	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
1348	zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
1349	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
1350	END DO
1351
1352	!------------------------------------------------------------------------------
1353	! Initialization
1354	!------------------------------------------------------------------------------
1355	qlth(:,:)=0.
1356	qlenv(:,:)=0.
1357	qltot(:,:)=0.
1358	cth_vol(:,:)=0.
1359	cenv_vol(:,:)=0.
1360	ctot_vol(:,:)=0.
1361	cth_surf(:,:)=0.
1362	cenv_surf(:,:)=0.
1363	ctot_surf(:,:)=0.
1364	qcloud(:)=0.
1365	rdd=287.04
1366	cppd=1005.7
1367	pi=3.14159
1368	Lv=2.5e6
1369	sqrt2=sqrt(2.)
1370	sqrt2pi=sqrt(2.*pi)
1371
1372
1373	DO ind1=1,ngrid
1374	!-------------------------------------------------------------------------------
1375	!Both thermal and environment in the gridbox
1376	!-------------------------------------------------------------------------------
1377	IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
1378	!--------------------------------------------
1379	!calcul de qsat_env
1380	!--------------------------------------------
1381	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1382	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1383	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1384	qsatbef=MIN(0.5,qsatbef)
1385	zcor=1./(1.-retv*qsatbef)
1386	qsatbef=qsatbef*zcor
1387	zqsatenv(ind1,ind2)=qsatbef
1388	!--------------------------------------------
1389	!calcul de s_env
1390	!--------------------------------------------
1391	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1392	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1393	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1394	!--------------------------------------------
1395	!calcul de qsat_th
1396	!--------------------------------------------
1397	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1398	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1399	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1400	qsatbef=MIN(0.5,qsatbef)
1401	zcor=1./(1.-retv*qsatbef)
1402	qsatbef=qsatbef*zcor
1403	zqsatth(ind1,ind2)=qsatbef
1404	!--------------------------------------------
1405	!calcul de s_th
1406	!--------------------------------------------
1407	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1408	ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
1409	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
1410	!--------------------------------------------
1411	!calcul standard deviations bi-Gaussian PDF
1412	!--------------------------------------------
1413	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)
1414	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)
1415	xth=sth/(sqrt2*sigma_th)
1416	xenv=senv/(sqrt2*sigma_env)
1417	!--------------------------------------------
1418	!Cloud fraction by volume CF_vol
1419	!--------------------------------------------
1420	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
1421	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1422	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1423	!--------------------------------------------
1424	!Condensed water qc
1425	!--------------------------------------------
1426	qlth(ind1,ind2)=sigma_th((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth_vol(ind1,ind2))
1427	qlenv(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv_vol(ind1,ind2))
1428	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
1429	!--------------------------------------------
1430	!Cloud fraction by surface CF_surf
1431	!--------------------------------------------
1432	!Method Neggers et al. (2011) : ok for cumulus clouds only
1433	!beta=0.0044 (Jouhaud et al.2018)
1434	!inverse_rho=1.+beta*dz(ind1,ind2)
1435	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1436	!Method Brooks et al. (2005) : ok for all types of clouds
1437	a_Brooks=0.6694
1438	b_Brooks=0.1882
1439	A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
1440	Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
1441	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1442	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1443	!--------------------------------------------
1444	!Incloud Condensed water qcloud
1445	!--------------------------------------------
1446	if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
1447	ctot_vol(ind1,ind2)=0.
1448	ctot_surf(ind1,ind2)=0.
1449	qcloud(ind1)=zqsatenv(ind1,ind2)
1450	else
1451	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
1452	endif
1453
1454
1455
1456	!-------------------------------------------------------------------------------
1457	!Environment only in the gridbox
1458	!-------------------------------------------------------------------------------
1459	ELSE
1460	!--------------------------------------------
1461	!calcul de qsat_env
1462	!--------------------------------------------
1463	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1464	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1465	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1466	qsatbef=MIN(0.5,qsatbef)
1467	zcor=1./(1.-retv*qsatbef)
1468	qsatbef=qsatbef*zcor
1469	zqsatenv(ind1,ind2)=qsatbef
1470	!--------------------------------------------
1471	!calcul de s_env
1472	!--------------------------------------------
1473	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84 these Arnaud Jam
1474	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
1475	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
1476	!--------------------------------------------
1477	!calcul standard deviations Gaussian PDF
1478	!--------------------------------------------
1479	zqenv(ind1)=po(ind1)
1480	sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
1481	xenv=senv/(sqrt2*sigma_env)
1482	!--------------------------------------------
1483	!Cloud fraction by volume CF_vol
1484	!--------------------------------------------
1485	ctot_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
1486	!--------------------------------------------
1487	!Condensed water qc
1488	!--------------------------------------------
1489	qltot(ind1,ind2)=sigma_env((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*ctot_vol(ind1,ind2))
1490	!--------------------------------------------
1491	!Cloud fraction by surface CF_surf
1492	!--------------------------------------------
1493	!Method Neggers et al. (2011) : ok for cumulus clouds only
1494	!beta=0.0044 (Jouhaud et al.2018)
1495	!inverse_rho=1.+beta*dz(ind1,ind2)
1496	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
1497	!Method Brooks et al. (2005) : ok for all types of clouds
1498	a_Brooks=0.6694
1499	b_Brooks=0.1882
1500	A_Maj_Brooks=0.1635 !-- sans dependence au shear
1501	Dx_Brooks=200000.
1502	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
1503	ctot_surf(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1504	!--------------------------------------------
1505	!Incloud Condensed water qcloud
1506	!--------------------------------------------
1507	if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
1508	ctot_vol(ind1,ind2)=0.
1509	ctot_surf(ind1,ind2)=0.
1510	qcloud(ind1)=zqsatenv(ind1,ind2)
1511	else
1512	qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
1513	endif
1514
1515
1516	END IF ! From the separation (thermal/envrionnement) et (environnement only)
1517
1518	! Outputs used to check the PDFs
1519	cloudth_senv(ind1,ind2) = senv
1520	cloudth_sth(ind1,ind2) = sth
1521	cloudth_sigmaenv(ind1,ind2) = sigma_env
1522	cloudth_sigmath(ind1,ind2) = sigma_th
1523
1524	END DO ! From the loop on ngrid
1525	return
1526
1527	END SUBROUTINE cloudth_v6
1528	END MODULE cloudth_mod
1529
1530
1531
1532

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

Download in other formats:

Original Format