Context Navigation

lmdz_cloudth.f90 @ 5480

Last change on this file since 5480 was 5479, checked in by yann meurdesoif, 2 weeks ago
Optimization : setting array to 0 only at computed vertical level. Very sensitive when running on gpu YM
File size: 87.4 KB

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

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

Download in other formats:

Original Format