Context Navigation

cloudth_mod.F90 @ 3847

Last change on this file since 3847 was 3570, checked in by fhourdin, 5 years ago
Sortie de nouveaux parametres pour controler la distribution des nuages.
File size: 57.1 KB

Rev	Line
[2686]	1	MODULE cloudth_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	SUBROUTINE cloudth(ngrid,klev,ind2, &
	8	& ztv,po,zqta,fraca, &
[3493]	9	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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)
[3493]	40	REAL pplay(ngrid,klev)
[2686]	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, &
[3493]	81	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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, &
[3493]	254	& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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)
[3493]	285	REAL pplay(ngrid,klev)
[2686]	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
[2909]	315	!$OMP THREADPRIVATE(vert_alpha)
[2686]	316	LOGICAL, SAVE :: firstcall = .TRUE.
[2909]	317	!$OMP THREADPRIVATE(firstcall)
[2686]	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
[3493]	589
	590
	591
[2686]	592	SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
	593	& ztv,po,zqta,fraca, &
[3493]	594	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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)
[3493]	625	REAL pplay(ngrid,klev)
[2686]	626	REAL ztla(ngrid,klev)
	627	REAL zthl(ngrid,klev)
	628
	629	REAL zqsatth(ngrid,klev)
	630	REAL zqsatenv(ngrid,klev)
	631
[2945]	632	REAL sigma1(ngrid,klev)
[2686]	633	REAL sigma2(ngrid,klev)
	634	REAL qlth(ngrid,klev)
	635	REAL qlenv(ngrid,klev)
	636	REAL qltot(ngrid,klev)
[2945]	637	REAL cth(ngrid,klev)
[2686]	638	REAL cenv(ngrid,klev)
	639	REAL ctot(ngrid,klev)
[2945]	640	REAL cth_vol(ngrid,klev)
	641	REAL cenv_vol(ngrid,klev)
	642	REAL ctot_vol(ngrid,klev)
	643	REAL rneb(ngrid,klev)
[2686]	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
[3493]	649	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2686]	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, &
[3493]	663	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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.
[2945]	683	cth_vol(:,:)=0.
	684	cenv_vol(:,:)=0.
	685	ctot_vol(:,:)=0.
[2686]	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)
[2945]	761	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	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))
[2945]	795	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	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, &
[3493]	817	& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	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
[2958]	828	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
	829	& cloudth_sigmath,cloudth_sigmaenv
[2686]	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)
[3493]	850	REAL pplay(ngrid,klev)
[2686]	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)
[2945]	862	REAL cth(ngrid,klev)
[2686]	863	REAL cenv(ngrid,klev)
	864	REAL ctot(ngrid,klev)
[2945]	865	REAL cth_vol(ngrid,klev)
	866	REAL cenv_vol(ngrid,klev)
	867	REAL ctot_vol(ngrid,klev)
[2686]	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
[2957]	873	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
[3493]	874	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2957]	875	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
[2686]	876	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
[2945]	877	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
[2686]	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)
[2911]	883	REAL,SAVE :: vert_alpha, vert_alpha_th
	884	!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
[2957]	885	REAL,SAVE :: sigma1s_factor=1.1
	886	REAL,SAVE :: sigma1s_power=0.6
[3570]	887	REAL,SAVE :: sigma2s_factor=0.09
	888	REAL,SAVE :: sigma2s_power=0.5
[2960]	889	REAL,SAVE :: cloudth_ratqsmin=-1.
[3570]	890	!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
[2957]	891	INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
	892	!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
	893
[2686]	894	LOGICAL, SAVE :: firstcall = .TRUE.
[2910]	895	!$OMP THREADPRIVATE(firstcall)
[2686]	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
[3493]	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
[2686]	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.
[2945]	927	cth_vol(:,:)=0.
	928	cenv_vol(:,:)=0.
	929	ctot_vol(:,:)=0.
[2686]	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
[2911]	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
[2957]	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
[3570]	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
[2960]	961	! Minimum value for the environmental air subgrid water distrib
	962	CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
	963	WRITE(,) 'cloudth_ratqsmin = ', cloudth_ratqsmin
[2957]	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
[2686]	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
[2957]	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
[2960]	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
[2957]	1028	sigma1s = sigma1s_fraca + sigma1s_ratqs
[3570]	1029	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
[2686]	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
[2911]	1080	ELSE IF (iflag_cloudth_vert >= 3) THEN
[3493]	1081	IF (iflag_cloudth_vert < 5) THEN
[2686]	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)
[2911]	1089	IF (iflag_cloudth_vert == 3) THEN
	1090	deltasenv=aenvvert_alphasigma1s
	1091	deltasth=athvert_alpha_thsigma2s
	1092	ELSE IF (iflag_cloudth_vert == 4) THEN
[2957]	1093	IF (iflag_cloudth_vert_noratqs == 1) THEN
[2959]	1094	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
[2957]	1095	deltasth=vert_alpha_th*sigma2s
	1096	ELSE
	1097	deltasenv=vert_alpha*sigma1s
	1098	deltasth=vert_alpha_th*sigma2s
	1099	ENDIF
[2911]	1100	ENDIF
[2686]	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
[2945]	1111	!CF_surfacique
[2686]	1112	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	1113	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
[2945]	1114	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[2686]	1115
[2945]	1116
	1117	!CF_volumique & eau condense
[2686]	1118	!environnement
	1119	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
[2945]	1120	IntJ_CF=0.5(1.-1.erf(xenv2))
[2686]	1121	if (deltasenv .lt. 1.e-10) then
	1122	qlenv(ind1,ind2)=IntJ
[2945]	1123	cenv_vol(ind1,ind2)=IntJ_CF
[2686]	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)
[2945]	1128	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1129	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
[2686]	1130	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1131	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1132	endif
	1133
	1134	!thermique
	1135	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
[2945]	1136	IntJ_CF=0.5(1.-1.erf(xth2))
	1137	if (deltasth .lt. 1.e-10) then
[2686]	1138	qlth(ind1,ind2)=IntJ
[2945]	1139	cth_vol(ind1,ind2)=IntJ_CF
[2686]	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)
[2945]	1144	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1145	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
[2686]	1146	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1147	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1148	endif
	1149
	1150	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
[2945]	1151	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
[2686]	1152
[3493]	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)
[2911]	1201	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1202
[3495]	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
[2686]	1205	ctot(ind1,ind2)=0.
[2945]	1206	ctot_vol(ind1,ind2)=0.
[3493]	1207	qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]	1208
[3493]	1209	else
[2686]	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
[3493]	1217	else ! gaussienne environnement seule
[2686]	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
[3493]	1229	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
[2686]	1230	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
[3495]	1231	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
[2686]	1232	aenv=1./(1.+(alenv*Lv/cppd))
[3495]	1233	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	1234	sth=0.
[3493]	1235
[2686]	1236
	1237	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[3495]	1238	sigma2s=0.
[2686]	1239
[3493]	1240	sqrt2pi=sqrt(2.*pi)
	1241	xenv=senv/(sqrt(2.)*sigma1s)
[2686]	1242	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
[3495]	1243	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[3493]	1244	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
[2686]	1245
	1246	if (ctot(ind1,ind2).lt.1.e-3) then
	1247	ctot(ind1,ind2)=0.
[3495]	1248	qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]	1249
	1250	else
	1251
[3495]	1252	! ctot(ind1,ind2)=ctot(ind1,ind2)
[2686]	1253	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	1254
[3493]	1255	endif
[2686]	1256
[3493]	1257
	1258
	1259
[2686]	1260	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]	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
[2686]	1267	enddo ! from the loop on ngrid l.333
	1268	return
	1269	! end
	1270	END SUBROUTINE cloudth_vert_v3
	1271	!
[3493]	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
[2686]	1536	END MODULE cloudth_mod
[3493]	1537
	1538
	1539
	1540

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

Download in other formats:

Original Format