Context Navigation

lmdz_cloudth.F90 @ 5097

Last change on this file since 5097 was 5087, checked in by abarral, 4 months ago
remove fixed-form ^{\s+& remaining in .f90,.F90}
File size: 86.8 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, &
[5087]	9	ztv,po,zqta,fraca, &
	10	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
	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
[2686]	17	IMPLICIT NONE
	18
	19
	20	!===========================================================================
	21	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	22	! Date : 25 Mai 2010
	23	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	24	!===========================================================================
	25
	26
[4593]	27	INCLUDE "YOMCST.h"
	28	INCLUDE "YOETHF.h"
	29	INCLUDE "FCTTRE.h"
[2686]	30
	31	INTEGER itap,ind1,ind2
	32	INTEGER ngrid,klev,klon,l,ig
[4651]	33	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
[2686]	34
	35	REAL ztv(ngrid,klev)
	36	REAL po(ngrid)
	37	REAL zqenv(ngrid)
	38	REAL zqta(ngrid,klev)
	39
	40	REAL fraca(ngrid,klev+1)
	41	REAL zpspsk(ngrid,klev)
	42	REAL paprs(ngrid,klev+1)
[3493]	43	REAL pplay(ngrid,klev)
[2686]	44	REAL ztla(ngrid,klev)
	45	REAL zthl(ngrid,klev)
	46
	47	REAL zqsatth(ngrid,klev)
	48	REAL zqsatenv(ngrid,klev)
	49
	50
	51	REAL sigma1(ngrid,klev)
	52	REAL sigma2(ngrid,klev)
	53	REAL qlth(ngrid,klev)
	54	REAL qlenv(ngrid,klev)
	55	REAL qltot(ngrid,klev)
	56	REAL cth(ngrid,klev)
	57	REAL cenv(ngrid,klev)
	58	REAL ctot(ngrid,klev)
	59	REAL rneb(ngrid,klev)
	60	REAL t(ngrid,klev)
	61	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
	62	REAL rdd,cppd,Lv
	63	REAL alth,alenv,ath,aenv
	64	REAL sth,senv,sigma1s,sigma2s,xth,xenv
	65	REAL Tbef,zdelta,qsatbef,zcor
	66	REAL qlbef
	67	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	68
	69	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	70	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	71	REAL zqs(ngrid), qcloud(ngrid)
	72	REAL erf
	73
	74
	75
	76
	77	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	78	! Gestion de deux versions de cloudth
	79	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	80
[5082]	81	IF (iflag_cloudth_vert>=1) THEN
[2686]	82	CALL cloudth_vert(ngrid,klev,ind2, &
[5087]	83	ztv,po,zqta,fraca, &
	84	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
	85	ratqs,zqs,t)
[2686]	86	RETURN
	87	ENDIF
	88	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	89
	90
	91	!-------------------------------------------------------------------------------
	92	! Initialisation des variables r?elles
	93	!-------------------------------------------------------------------------------
	94	sigma1(:,:)=0.
	95	sigma2(:,:)=0.
	96	qlth(:,:)=0.
	97	qlenv(:,:)=0.
	98	qltot(:,:)=0.
	99	rneb(:,:)=0.
	100	qcloud(:)=0.
	101	cth(:,:)=0.
	102	cenv(:,:)=0.
	103	ctot(:,:)=0.
	104	qsatmmussig1=0.
	105	qsatmmussig2=0.
	106	rdd=287.04
	107	cppd=1005.7
	108	pi=3.14159
	109	Lv=2.5e6
	110	sqrt2pi=sqrt(2.*pi)
	111
	112
	113
	114	!-------------------------------------------------------------------------------
	115	! Calcul de la fraction du thermique et des ?cart-types des distributions
	116	!-------------------------------------------------------------------------------
	117	do ind1=1,ngrid
	118
[5082]	119	if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then
[2686]	120
	121	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	122
	123
	124	! zqenv(ind1)=po(ind1)
	125	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	126	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	127	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	128	qsatbef=MIN(0.5,qsatbef)
	129	zcor=1./(1.-retv*qsatbef)
	130	qsatbef=qsatbef*zcor
	131	zqsatenv(ind1,ind2)=qsatbef
	132
	133
	134
	135
	136	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	137	aenv=1./(1.+(alenv*Lv/cppd))
	138	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	139
	140
	141
	142
	143	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	144	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	145	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	146	qsatbef=MIN(0.5,qsatbef)
	147	zcor=1./(1.-retv*qsatbef)
	148	qsatbef=qsatbef*zcor
	149	zqsatth(ind1,ind2)=qsatbef
	150
	151	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
	152	ath=1./(1.+(alth*Lv/cppd))
	153	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
	154
	155
	156
	157	!------------------------------------------------------------------------------
	158	! Calcul des ?cart-types pour s
	159	!------------------------------------------------------------------------------
	160
	161	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	162	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
	163	! if (paprs(ind1,ind2).gt.90000) then
	164	! ratqs(ind1,ind2)=0.002
	165	! else
	166	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	167	! endif
	168	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	169	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	170	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	171	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	172
	173	!------------------------------------------------------------------------------
	174	! Calcul de l'eau condens?e et de la couverture nuageuse
	175	!------------------------------------------------------------------------------
	176	sqrt2pi=sqrt(2.*pi)
	177	xth=sth/(sqrt(2.)*sigma2s)
	178	xenv=senv/(sqrt(2.)*sigma1s)
	179	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	180	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	181	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	182
	183	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	184	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	185	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	186
	187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[5082]	188	if (ctot(ind1,ind2)<1.e-10) then
[2686]	189	ctot(ind1,ind2)=0.
	190	qcloud(ind1)=zqsatenv(ind1,ind2)
	191
	192	else
	193
	194	ctot(ind1,ind2)=ctot(ind1,ind2)
	195	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	196
	197	endif
	198
	199
	200	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	201
	202
	203	else ! gaussienne environnement seule
	204
	205	zqenv(ind1)=po(ind1)
	206	Tbef=t(ind1,ind2)
	207	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	208	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	209	qsatbef=MIN(0.5,qsatbef)
	210	zcor=1./(1.-retv*qsatbef)
	211	qsatbef=qsatbef*zcor
	212	zqsatenv(ind1,ind2)=qsatbef
	213
	214
	215	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	216	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	217	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	218	aenv=1./(1.+(alenv*Lv/cppd))
	219	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	220
	221
	222	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	223
	224	sqrt2pi=sqrt(2.*pi)
	225	xenv=senv/(sqrt(2.)*sigma1s)
	226	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	227	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	228
[5082]	229	if (ctot(ind1,ind2)<1.e-3) then
[2686]	230	ctot(ind1,ind2)=0.
	231	qcloud(ind1)=zqsatenv(ind1,ind2)
	232
	233	else
	234
	235	ctot(ind1,ind2)=ctot(ind1,ind2)
	236	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	237
	238	endif
	239
	240
	241
	242
	243
	244
	245	endif
	246	enddo
	247
	248	return
	249	! end
	250	END SUBROUTINE cloudth
	251
	252
	253
	254	!===========================================================================
	255	SUBROUTINE cloudth_vert(ngrid,klev,ind2, &
[5087]	256	ztv,po,zqta,fraca, &
	257	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
	258	ratqs,zqs,t)
[2686]	259
	260	!===========================================================================
	261	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	262	! Date : 25 Mai 2010
	263	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	264	!===========================================================================
	265
	266
[4651]	267	use lmdz_cloudth_ini, only: iflag_cloudth_vert, vert_alpha
[2686]	268
	269	IMPLICIT NONE
	270
[4593]	271	INCLUDE "YOMCST.h"
	272	INCLUDE "YOETHF.h"
	273	INCLUDE "FCTTRE.h"
[2686]	274
	275	INTEGER itap,ind1,ind2
	276	INTEGER ngrid,klev,klon,l,ig
	277
	278	REAL ztv(ngrid,klev)
	279	REAL po(ngrid)
	280	REAL zqenv(ngrid)
	281	REAL zqta(ngrid,klev)
	282
	283	REAL fraca(ngrid,klev+1)
	284	REAL zpspsk(ngrid,klev)
	285	REAL paprs(ngrid,klev+1)
[3493]	286	REAL pplay(ngrid,klev)
[2686]	287	REAL ztla(ngrid,klev)
	288	REAL zthl(ngrid,klev)
	289
	290	REAL zqsatth(ngrid,klev)
	291	REAL zqsatenv(ngrid,klev)
	292
	293
	294	REAL sigma1(ngrid,klev)
	295	REAL sigma2(ngrid,klev)
	296	REAL qlth(ngrid,klev)
	297	REAL qlenv(ngrid,klev)
	298	REAL qltot(ngrid,klev)
	299	REAL cth(ngrid,klev)
	300	REAL cenv(ngrid,klev)
	301	REAL ctot(ngrid,klev)
	302	REAL rneb(ngrid,klev)
	303	REAL t(ngrid,klev)
	304	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
	305	REAL rdd,cppd,Lv,sqrt2,sqrtpi
	306	REAL alth,alenv,ath,aenv
	307	REAL sth,senv,sigma1s,sigma2s,xth,xenv
	308	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
	309	REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
	310	REAL Tbef,zdelta,qsatbef,zcor
	311	REAL qlbef
	312	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	313	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	314	! (J Jouhaud, JL Dufresne, JB Madeleine)
	315
	316	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	317	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	318	REAL zqs(ngrid), qcloud(ngrid)
	319	REAL erf
	320
	321	!------------------------------------------------------------------------------
	322	! Initialisation des variables r?elles
	323	!------------------------------------------------------------------------------
	324	sigma1(:,:)=0.
	325	sigma2(:,:)=0.
	326	qlth(:,:)=0.
	327	qlenv(:,:)=0.
	328	qltot(:,:)=0.
	329	rneb(:,:)=0.
	330	qcloud(:)=0.
	331	cth(:,:)=0.
	332	cenv(:,:)=0.
	333	ctot(:,:)=0.
	334	qsatmmussig1=0.
	335	qsatmmussig2=0.
	336	rdd=287.04
	337	cppd=1005.7
	338	pi=3.14159
	339	Lv=2.5e6
	340	sqrt2pi=sqrt(2.*pi)
	341	sqrt2=sqrt(2.)
	342	sqrtpi=sqrt(pi)
	343
	344	!-------------------------------------------------------------------------------
	345	! Calcul de la fraction du thermique et des ?cart-types des distributions
	346	!-------------------------------------------------------------------------------
	347	do ind1=1,ngrid
	348
[5082]	349	if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then
[2686]	350
	351	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	352
	353
	354	! zqenv(ind1)=po(ind1)
	355	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	356	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	357	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	358	qsatbef=MIN(0.5,qsatbef)
	359	zcor=1./(1.-retv*qsatbef)
	360	qsatbef=qsatbef*zcor
	361	zqsatenv(ind1,ind2)=qsatbef
	362
	363
	364
	365
	366	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	367	aenv=1./(1.+(alenv*Lv/cppd))
	368	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	369
	370
	371
	372
	373	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	374	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	375	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	376	qsatbef=MIN(0.5,qsatbef)
	377	zcor=1./(1.-retv*qsatbef)
	378	qsatbef=qsatbef*zcor
	379	zqsatth(ind1,ind2)=qsatbef
	380
	381	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
	382	ath=1./(1.+(alth*Lv/cppd))
	383	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
	384
	385
	386
	387	!------------------------------------------------------------------------------
	388	! Calcul des ?cart-types pour s
	389	!------------------------------------------------------------------------------
	390
	391	sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	392	sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	393	! if (paprs(ind1,ind2).gt.90000) then
	394	! ratqs(ind1,ind2)=0.002
	395	! else
	396	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	397	! endif
	398	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	399	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	400	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	401	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	402
	403	!------------------------------------------------------------------------------
	404	! Calcul de l'eau condens?e et de la couverture nuageuse
	405	!------------------------------------------------------------------------------
	406	sqrt2pi=sqrt(2.*pi)
	407	xth=sth/(sqrt(2.)*sigma2s)
	408	xenv=senv/(sqrt(2.)*sigma1s)
	409	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	410	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	411	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	412
	413	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	414	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	415	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	416
	417	IF (iflag_cloudth_vert == 1) THEN
	418	!-------------------------------------------------------------------------------
	419	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
	420	!-------------------------------------------------------------------------------
	421	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	422	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	423	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	424	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	425	! deltasenv=aenv0.01po(ind1)
	426	! deltasth=ath0.01zqta(ind1,ind2)
	427	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
	428	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
	429	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
	430	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
	431	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	432	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	433
	434	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
	435	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
	436	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	437
	438	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
	439	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	440	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
	441	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
	442
	443	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	444	! qlenv(ind1,ind2)=IntJ
	445	! print*, qlenv(ind1,ind2),'VERIF EAU'
	446
	447
	448	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
	449	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
	450	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
	451	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	452	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
	453	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
	454	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	455	! qlth(ind1,ind2)=IntJ
	456	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	457	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	458
	459	ELSE IF (iflag_cloudth_vert == 2) THEN
	460
	461	!-------------------------------------------------------------------------------
	462	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
	463	!-------------------------------------------------------------------------------
	464	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	465	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	466	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	467	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	468	deltasenv=aenvvert_alphasigma1s
	469	deltasth=athvert_alphasigma2s
	470
	471	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	472	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	473	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	474	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	475	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	476	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	477
	478	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	479	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
	480	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	481
	482	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)exp(-1.xenv2**2)
	483	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	484	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	485	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))(exp(-1.xenv1*2)-exp(-1.xenv2**2))
	486
	487	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
	488	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
	489	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
	490
	491	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	492	! qlenv(ind1,ind2)=IntJ
	493	! print*, qlenv(ind1,ind2),'VERIF EAU'
	494
	495	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)exp(-1.xth2**2)
	496	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	497	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	498	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))(exp(-1.xth1*2)-exp(-1.xth2**2))
	499
	500	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	501	! qlth(ind1,ind2)=IntJ
	502	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	503	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	504
	505
	506
	507
	508	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
	509
	510	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	511
[5082]	512	if (cenv(ind1,ind2)<1.e-10.or.cth(ind1,ind2)<1.e-10) then
[2686]	513	ctot(ind1,ind2)=0.
	514	qcloud(ind1)=zqsatenv(ind1,ind2)
	515
	516	else
	517
	518	ctot(ind1,ind2)=ctot(ind1,ind2)
	519	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	520	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	521	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	522
	523	endif
	524
	525
	526
	527	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	528
	529
	530	else ! gaussienne environnement seule
	531
	532	zqenv(ind1)=po(ind1)
	533	Tbef=t(ind1,ind2)
	534	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	535	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	536	qsatbef=MIN(0.5,qsatbef)
	537	zcor=1./(1.-retv*qsatbef)
	538	qsatbef=qsatbef*zcor
	539	zqsatenv(ind1,ind2)=qsatbef
	540
	541
	542	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	543	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	544	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	545	aenv=1./(1.+(alenv*Lv/cppd))
	546	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	547
	548
	549	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	550
	551	sqrt2pi=sqrt(2.*pi)
	552	xenv=senv/(sqrt(2.)*sigma1s)
	553	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	554	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	555
[5082]	556	if (ctot(ind1,ind2)<1.e-3) then
[2686]	557	ctot(ind1,ind2)=0.
	558	qcloud(ind1)=zqsatenv(ind1,ind2)
	559
	560	else
	561
	562	ctot(ind1,ind2)=ctot(ind1,ind2)
	563	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	564
	565	endif
	566
	567
	568
	569
	570
	571
	572	endif
	573	enddo
	574
	575	return
	576	! end
	577	END SUBROUTINE cloudth_vert
	578
[3493]	579
	580
	581
[2686]	582	SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
[5087]	583	ztv,po,zqta,fraca, &
	584	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	585	ratqs,zqs,t, &
	586	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[2686]	587
[4651]	588	use lmdz_cloudth_ini, only: iflag_cloudth_vert
[2686]	589
	590	IMPLICIT NONE
	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	!===========================================================================
	598
	599
[4593]	600	INCLUDE "YOMCST.h"
	601	INCLUDE "YOETHF.h"
	602	INCLUDE "FCTTRE.h"
[2686]	603
[4674]	604	integer, intent(in) :: ind2
	605	integer, intent(in) :: ngrid,klev
	606
	607	real, dimension(ngrid,klev), intent(in) :: ztv
	608	real, dimension(ngrid), intent(in) :: po
	609	real, dimension(ngrid,klev), intent(in) :: zqta
	610	real, dimension(ngrid,klev+1), intent(in) :: fraca
	611	real, dimension(ngrid), intent(out) :: qcloud
	612	real, dimension(ngrid,klev), intent(out) :: ctot
	613	real, dimension(ngrid,klev), intent(out) :: ctot_vol
	614	real, dimension(ngrid,klev), intent(in) :: zpspsk
	615	real, dimension(ngrid,klev+1), intent(in) :: paprs
	616	real, dimension(ngrid,klev), intent(in) :: pplay
	617	real, dimension(ngrid,klev), intent(in) :: ztla
	618	real, dimension(ngrid,klev), intent(inout) :: zthl
	619	real, dimension(ngrid,klev), intent(in) :: ratqs
	620	real, dimension(ngrid), intent(in) :: zqs
	621	real, dimension(ngrid,klev), intent(in) :: t
[4651]	622	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
[4674]	623
	624
[2686]	625	REAL zqenv(ngrid)
	626	REAL zqsatth(ngrid,klev)
	627	REAL zqsatenv(ngrid,klev)
	628
[2945]	629	REAL sigma1(ngrid,klev)
[2686]	630	REAL sigma2(ngrid,klev)
	631	REAL qlth(ngrid,klev)
	632	REAL qlenv(ngrid,klev)
	633	REAL qltot(ngrid,klev)
[2945]	634	REAL cth(ngrid,klev)
[2686]	635	REAL cenv(ngrid,klev)
[2945]	636	REAL cth_vol(ngrid,klev)
	637	REAL cenv_vol(ngrid,klev)
	638	REAL rneb(ngrid,klev)
[2686]	639	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
	640	REAL rdd,cppd,Lv
	641	REAL alth,alenv,ath,aenv
	642	REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
[3493]	643	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2686]	644	REAL Tbef,zdelta,qsatbef,zcor
	645	REAL qlbef
	646	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	647	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	648	REAL erf
	649
	650
[4674]	651	INTEGER :: ind1,l, ig
	652
[5082]	653	IF (iflag_cloudth_vert>=1) THEN
[2686]	654	CALL cloudth_vert_v3(ngrid,klev,ind2, &
[5087]	655	ztv,po,zqta,fraca, &
	656	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	657	ratqs,zqs,t, &
	658	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[2686]	659	RETURN
	660	ENDIF
	661	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	662
	663
	664	!-------------------------------------------------------------------------------
	665	! Initialisation des variables r?elles
	666	!-------------------------------------------------------------------------------
	667	sigma1(:,:)=0.
	668	sigma2(:,:)=0.
	669	qlth(:,:)=0.
	670	qlenv(:,:)=0.
	671	qltot(:,:)=0.
	672	rneb(:,:)=0.
	673	qcloud(:)=0.
	674	cth(:,:)=0.
	675	cenv(:,:)=0.
	676	ctot(:,:)=0.
[2945]	677	cth_vol(:,:)=0.
	678	cenv_vol(:,:)=0.
	679	ctot_vol(:,:)=0.
[2686]	680	qsatmmussig1=0.
	681	qsatmmussig2=0.
	682	rdd=287.04
	683	cppd=1005.7
	684	pi=3.14159
	685	Lv=2.5e6
	686	sqrt2pi=sqrt(2.*pi)
	687	sqrt2=sqrt(2.)
	688	sqrtpi=sqrt(pi)
	689
	690
	691	!-------------------------------------------------------------------------------
	692	! Cloud fraction in the thermals and standard deviation of the PDFs
	693	!-------------------------------------------------------------------------------
	694	do ind1=1,ngrid
	695
[5082]	696	if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then
[2686]	697
	698	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	699
	700	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	701	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	702	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	703	qsatbef=MIN(0.5,qsatbef)
	704	zcor=1./(1.-retv*qsatbef)
	705	qsatbef=qsatbef*zcor
	706	zqsatenv(ind1,ind2)=qsatbef
	707
	708
	709	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	710	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	711	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	712
	713	!po = qt de l'environnement ET des thermique
	714	!zqenv = qt environnement
	715	!zqsatenv = qsat environnement
	716	!zthl = Tl environnement
	717
	718
	719	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	720	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	721	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	722	qsatbef=MIN(0.5,qsatbef)
	723	zcor=1./(1.-retv*qsatbef)
	724	qsatbef=qsatbef*zcor
	725	zqsatth(ind1,ind2)=qsatbef
	726
	727	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	728	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	729	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	730
	731	!zqta = qt thermals
	732	!zqsatth = qsat thermals
	733	!ztla = Tl thermals
	734
	735	!------------------------------------------------------------------------------
	736	! s standard deviations
	737	!------------------------------------------------------------------------------
	738
	739	! tests
	740	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	741	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
	742	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	743	! final option
	744	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	745	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	746
	747	!------------------------------------------------------------------------------
	748	! Condensed water and cloud cover
	749	!------------------------------------------------------------------------------
	750	xth=sth/(sqrt2*sigma2s)
	751	xenv=senv/(sqrt2*sigma1s)
	752	cth(ind1,ind2)=0.5(1.+1.erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	753	cenv(ind1,ind2)=0.5(1.+1.erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	754	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[2945]	755	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	756
	757	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth(ind1,ind2))
	758	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	759	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	760
[5082]	761	if (ctot(ind1,ind2)<1.e-10) then
[2686]	762	ctot(ind1,ind2)=0.
	763	qcloud(ind1)=zqsatenv(ind1,ind2)
	764	else
	765	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	766	endif
	767
	768	else ! Environnement only, follow the if l.110
	769
	770	zqenv(ind1)=po(ind1)
	771	Tbef=t(ind1,ind2)
	772	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	773	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	774	qsatbef=MIN(0.5,qsatbef)
	775	zcor=1./(1.-retv*qsatbef)
	776	qsatbef=qsatbef*zcor
	777	zqsatenv(ind1,ind2)=qsatbef
	778
	779	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	780	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	781	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	782	aenv=1./(1.+(alenv*Lv/cppd))
	783	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	784
	785	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	786
	787	xenv=senv/(sqrt2*sigma1s)
	788	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
[2945]	789	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	790	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	791
[5082]	792	if (ctot(ind1,ind2)<1.e-3) then
[2686]	793	ctot(ind1,ind2)=0.
	794	qcloud(ind1)=zqsatenv(ind1,ind2)
	795	else
	796	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	797	endif
	798
	799
	800	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
	801	enddo ! from the loop on ngrid l.108
	802	return
	803	! end
	804	END SUBROUTINE cloudth_v3
	805
	806
	807
	808	!===========================================================================
	809	SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
[5087]	810	ztv,po,zqta,fraca, &
	811	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	812	ratqs,zqs,t, &
	813	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[2686]	814
	815	!===========================================================================
	816	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	817	! Date : 25 Mai 2010
	818	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	819	!===========================================================================
	820
[4651]	821	use lmdz_cloudth_ini, only : iflag_cloudth_vert,iflag_ratqs
	822	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]	823
	824	IMPLICIT NONE
	825
[4651]	826
	827
[4593]	828	INCLUDE "YOMCST.h"
	829	INCLUDE "YOETHF.h"
	830	INCLUDE "FCTTRE.h"
[2686]	831
	832	INTEGER itap,ind1,ind2
	833	INTEGER ngrid,klev,klon,l,ig
[4651]	834	real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
[2686]	835
	836	REAL ztv(ngrid,klev)
	837	REAL po(ngrid)
	838	REAL zqenv(ngrid)
	839	REAL zqta(ngrid,klev)
	840
	841	REAL fraca(ngrid,klev+1)
	842	REAL zpspsk(ngrid,klev)
	843	REAL paprs(ngrid,klev+1)
[3493]	844	REAL pplay(ngrid,klev)
[2686]	845	REAL ztla(ngrid,klev)
	846	REAL zthl(ngrid,klev)
	847
	848	REAL zqsatth(ngrid,klev)
	849	REAL zqsatenv(ngrid,klev)
	850
	851	REAL sigma1(ngrid,klev)
	852	REAL sigma2(ngrid,klev)
	853	REAL qlth(ngrid,klev)
	854	REAL qlenv(ngrid,klev)
	855	REAL qltot(ngrid,klev)
[2945]	856	REAL cth(ngrid,klev)
[2686]	857	REAL cenv(ngrid,klev)
	858	REAL ctot(ngrid,klev)
[2945]	859	REAL cth_vol(ngrid,klev)
	860	REAL cenv_vol(ngrid,klev)
	861	REAL ctot_vol(ngrid,klev)
[2686]	862	REAL rneb(ngrid,klev)
	863	REAL t(ngrid,klev)
	864	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
	865	REAL rdd,cppd,Lv
	866	REAL alth,alenv,ath,aenv
[2957]	867	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
[3493]	868	REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
[2957]	869	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
[2686]	870	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
[2945]	871	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
[2686]	872	REAL Tbef,zdelta,qsatbef,zcor
	873	REAL qlbef
	874	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	875	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	876	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2957]	877
[2686]	878	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	879	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	880	REAL zqs(ngrid), qcloud(ngrid)
	881	REAL erf
	882
[3493]	883	REAL rhodz(ngrid,klev)
	884	REAL zrho(ngrid,klev)
	885	REAL dz(ngrid,klev)
	886
	887	DO ind1 = 1, ngrid
	888	!Layer calculation
	889	rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
	890	zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
	891	dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
	892	END DO
	893
[2686]	894	!------------------------------------------------------------------------------
	895	! Initialize
	896	!------------------------------------------------------------------------------
[3999]	897
[2686]	898	sigma1(:,:)=0.
	899	sigma2(:,:)=0.
	900	qlth(:,:)=0.
	901	qlenv(:,:)=0.
	902	qltot(:,:)=0.
	903	rneb(:,:)=0.
	904	qcloud(:)=0.
	905	cth(:,:)=0.
	906	cenv(:,:)=0.
	907	ctot(:,:)=0.
[2945]	908	cth_vol(:,:)=0.
	909	cenv_vol(:,:)=0.
	910	ctot_vol(:,:)=0.
[2686]	911	qsatmmussig1=0.
	912	qsatmmussig2=0.
	913	rdd=287.04
	914	cppd=1005.7
	915	pi=3.14159
	916	Lv=2.5e6
	917	sqrt2pi=sqrt(2.*pi)
	918	sqrt2=sqrt(2.)
	919	sqrtpi=sqrt(pi)
	920
[2957]	921
[2686]	922
	923	!-------------------------------------------------------------------------------
	924	! Calcul de la fraction du thermique et des ecart-types des distributions
	925	!-------------------------------------------------------------------------------
	926	do ind1=1,ngrid
	927
[5082]	928	if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then !Thermal and environnement
[2686]	929
	930	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
	931
	932
	933	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	934	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	935	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	936	qsatbef=MIN(0.5,qsatbef)
	937	zcor=1./(1.-retv*qsatbef)
	938	qsatbef=qsatbef*zcor
	939	zqsatenv(ind1,ind2)=qsatbef
	940
	941
	942	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	943	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	944	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	945
	946	!zqenv = qt environnement
	947	!zqsatenv = qsat environnement
	948	!zthl = Tl environnement
	949
	950
	951	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	952	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	953	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	954	qsatbef=MIN(0.5,qsatbef)
	955	zcor=1./(1.-retv*qsatbef)
	956	qsatbef=qsatbef*zcor
	957	zqsatth(ind1,ind2)=qsatbef
	958
	959	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	960	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	961	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	962
	963
	964	!zqta = qt thermals
	965	!zqsatth = qsat thermals
	966	!ztla = Tl thermals
	967	!------------------------------------------------------------------------------
	968	! s standard deviation
	969	!------------------------------------------------------------------------------
	970
[2957]	971	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
[5087]	972	(1-fraca(ind1,ind2))((sth-senv)2)*0.5
[2957]	973	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
[2960]	974	IF (cloudth_ratqsmin>0.) THEN
	975	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
	976	ELSE
	977	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
	978	ENDIF
[2957]	979	sigma1s = sigma1s_fraca + sigma1s_ratqs
[5082]	980	IF (iflag_ratqs==11) then
[4623]	981	sigma1s = ratqs(ind1,ind2)po(ind1)aenv
	982	ENDIF
[3570]	983	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
[2686]	984	! tests
	985	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	986	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
	987	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	988	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
	989	! if (paprs(ind1,ind2).gt.90000) then
	990	! ratqs(ind1,ind2)=0.002
	991	! else
	992	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	993	! endif
	994	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	995	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	996	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	997	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	998
	999	IF (iflag_cloudth_vert == 1) THEN
	1000	!-------------------------------------------------------------------------------
	1001	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
	1002	!-------------------------------------------------------------------------------
	1003
	1004	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	1005	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	1006
	1007	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1008	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1009	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
	1010	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
	1011	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	1012	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	1013
	1014	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
	1015	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
	1016	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1017
	1018	! Environment
	1019	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
	1020	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	1021	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
	1022	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
	1023
	1024	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1025
	1026	! Thermal
	1027	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
	1028	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	1029	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
	1030	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
	1031	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1032	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	1033
[2911]	1034	ELSE IF (iflag_cloudth_vert >= 3) THEN
[3493]	1035	IF (iflag_cloudth_vert < 5) THEN
[2686]	1036	!-------------------------------------------------------------------------------
	1037	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
	1038	!-------------------------------------------------------------------------------
	1039	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	1040	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	1041	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	1042	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
[2911]	1043	IF (iflag_cloudth_vert == 3) THEN
	1044	deltasenv=aenvvert_alphasigma1s
	1045	deltasth=athvert_alpha_thsigma2s
	1046	ELSE IF (iflag_cloudth_vert == 4) THEN
[2957]	1047	IF (iflag_cloudth_vert_noratqs == 1) THEN
[2959]	1048	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
[2957]	1049	deltasth=vert_alpha_th*sigma2s
	1050	ELSE
	1051	deltasenv=vert_alpha*sigma1s
	1052	deltasth=vert_alpha_th*sigma2s
	1053	ENDIF
[2911]	1054	ENDIF
[2686]	1055
	1056	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1057	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1058	exp_xenv1 = exp(-1.xenv1*2)
	1059	exp_xenv2 = exp(-1.xenv2*2)
	1060	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	1061	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	1062	exp_xth1 = exp(-1.xth1*2)
	1063	exp_xth2 = exp(-1.xth2*2)
	1064
[2945]	1065	!CF_surfacique
[2686]	1066	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	1067	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
[2945]	1068	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[2686]	1069
[2945]	1070
	1071	!CF_volumique & eau condense
[2686]	1072	!environnement
	1073	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
[2945]	1074	IntJ_CF=0.5(1.-1.erf(xenv2))
[5082]	1075	if (deltasenv < 1.e-10) then
[2686]	1076	qlenv(ind1,ind2)=IntJ
[2945]	1077	cenv_vol(ind1,ind2)=IntJ_CF
[2686]	1078	else
	1079	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1080	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1081	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
[2945]	1082	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1083	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
[2686]	1084	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1085	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1086	endif
	1087
	1088	!thermique
	1089	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
[2945]	1090	IntJ_CF=0.5(1.-1.erf(xth2))
[5082]	1091	if (deltasth < 1.e-10) then
[2686]	1092	qlth(ind1,ind2)=IntJ
[2945]	1093	cth_vol(ind1,ind2)=IntJ_CF
[2686]	1094	else
	1095	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1096	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1097	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
[2945]	1098	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1099	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
[2686]	1100	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1101	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1102	endif
	1103
	1104	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
[2945]	1105	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
[2686]	1106
[3493]	1107	ELSE IF (iflag_cloudth_vert == 5) THEN
[4260]	1108	sigma1s=(0.71794+0.000498239dz(ind1,ind2))(fraca(ind1,ind2)**0.5) &
	1109	/(1-fraca(ind1,ind2))(((sth-senv)2)*0.5) &
	1110	+ratqs(ind1,ind2)*po(ind1) !Environment
[3493]	1111	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
	1112	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	1113	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	1114	xth=sth/(sqrt(2.)*sigma2s)
	1115	xenv=senv/(sqrt(2.)*sigma1s)
	1116
	1117	!Volumique
	1118	cth_vol(ind1,ind2)=0.5(1.+1.erf(xth))
	1119	cenv_vol(ind1,ind2)=0.5(1.+1.erf(xenv))
	1120	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
	1121	!print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
	1122
	1123	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth_vol(ind1,ind2))
	1124	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv_vol(ind1,ind2))
	1125	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	1126
	1127	!Surfacique
	1128	!Neggers
	1129	!beta=0.0044
	1130	!inverse_rho=1.+beta*dz(ind1,ind2)
	1131	!print *,'jeanjean : beta=',beta
	1132	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
	1133	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
	1134	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1135
	1136	!Brooks
	1137	a_Brooks=0.6694
	1138	b_Brooks=0.1882
	1139	A_Maj_Brooks=0.1635 !-- sans shear
	1140	!A_Maj_Brooks=0.17 !-- ARM LES
	1141	!A_Maj_Brooks=0.18 !-- RICO LES
	1142	!A_Maj_Brooks=0.19 !-- BOMEX LES
	1143	Dx_Brooks=200000.
	1144	f_Brooks=A_Maj_Brooks(dz(ind1,ind2)(a_Brooks))(Dx_Brooks**(-b_Brooks))
	1145	!print *,'jeanjean_f=',f_Brooks
	1146
	1147	cth(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
	1148	cenv(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
	1149	ctot(ind1,ind2)=1./(1.+exp(-1.f_Brooks)((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
	1150	!print *,'JJ_ctot_1',ctot(ind1,ind2)
	1151
	1152
	1153
	1154
	1155
	1156	ENDIF ! of if (iflag_cloudth_vert<5)
[2911]	1157	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1158
[3495]	1159	! if (ctot(ind1,ind2).lt.1.e-10) then
[5082]	1160	if (cenv(ind1,ind2)<1.e-10.or.cth(ind1,ind2)<1.e-10) then
[2686]	1161	ctot(ind1,ind2)=0.
[2945]	1162	ctot_vol(ind1,ind2)=0.
[3493]	1163	qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]	1164
[3493]	1165	else
[2686]	1166
	1167	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	1168	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	1169	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	1170
	1171	endif
	1172
[3493]	1173	else ! gaussienne environnement seule
[2686]	1174
[3999]	1175
[2686]	1176	zqenv(ind1)=po(ind1)
	1177	Tbef=t(ind1,ind2)
	1178	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	1179	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	1180	qsatbef=MIN(0.5,qsatbef)
	1181	zcor=1./(1.-retv*qsatbef)
	1182	qsatbef=qsatbef*zcor
	1183	zqsatenv(ind1,ind2)=qsatbef
	1184
	1185
[3493]	1186	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
[2686]	1187	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
[3495]	1188	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
[2686]	1189	aenv=1./(1.+(alenv*Lv/cppd))
[3495]	1190	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	1191	sth=0.
[3493]	1192
[2686]	1193
	1194	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[3495]	1195	sigma2s=0.
[2686]	1196
[3493]	1197	sqrt2pi=sqrt(2.*pi)
	1198	xenv=senv/(sqrt(2.)*sigma1s)
[2686]	1199	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
[3495]	1200	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[3493]	1201	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
[2686]	1202
[5082]	1203	if (ctot(ind1,ind2)<1.e-3) then
[2686]	1204	ctot(ind1,ind2)=0.
[3495]	1205	qcloud(ind1)=zqsatenv(ind1,ind2)
[2686]	1206
	1207	else
	1208
[3495]	1209	! ctot(ind1,ind2)=ctot(ind1,ind2)
[2686]	1210	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	1211
[3493]	1212	endif
[2686]	1213
[3493]	1214
	1215
	1216
[2686]	1217	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]	1218	! Outputs used to check the PDFs
	1219	cloudth_senv(ind1,ind2) = senv
	1220	cloudth_sth(ind1,ind2) = sth
	1221	cloudth_sigmaenv(ind1,ind2) = sigma1s
	1222	cloudth_sigmath(ind1,ind2) = sigma2s
	1223
[2686]	1224	enddo ! from the loop on ngrid l.333
	1225	return
	1226	! end
	1227	END SUBROUTINE cloudth_vert_v3
	1228	!
[3493]	1229
	1230
	1231
	1232
	1233
	1234
	1235
	1236
	1237
	1238
	1239
	1240	SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
[5087]	1241	ztv,po,zqta,fraca, &
	1242	qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	1243	ratqs,zqs,T, &
	1244	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3493]	1245
[4651]	1246	use lmdz_cloudth_ini, only: iflag_cloudth_vert
[3493]	1247
	1248	IMPLICIT NONE
	1249
[4651]	1250
[4593]	1251	INCLUDE "YOMCST.h"
	1252	INCLUDE "YOETHF.h"
	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
	1304	REAL erf,pi,sqrt2,sqrt2pi
	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	!------------------------------------------------------------------------------
	1318	qlth(:,:)=0.
	1319	qlenv(:,:)=0.
	1320	qltot(:,:)=0.
	1321	cth_vol(:,:)=0.
	1322	cenv_vol(:,:)=0.
	1323	ctot_vol(:,:)=0.
	1324	cth_surf(:,:)=0.
	1325	cenv_surf(:,:)=0.
	1326	ctot_surf(:,:)=0.
	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	!-------------------------------------------------------------------------------
[5082]	1340	IF ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) THEN
[3493]	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	!--------------------------------------------
[5082]	1411	if (ctot_surf(ind1,ind2) < 1.e-10) then
[3493]	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	!--------------------------------------------
[5082]	1472	if (ctot_surf(ind1,ind2) < 1.e-8) then
[3493]	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 erf
	1578	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
	1579	REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
[4910]	1580	REAL rhodz(klon)
	1581	REAL rho(klon)
	1582	REAL dz(klon)
	1583	REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
[4072]	1584	REAL alenvl, aenvl
	1585	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
[3999]	1586	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
[4072]	1587	REAL qmax, ttarget, stmp, cout, coutref, fraci
[4910]	1588	REAL maxi, mini, pas
[3999]	1589
	1590	! Modifty the saturation deficit PDF in thermals
	1591	! in the presence of ice crystals
	1592	CHARACTER (len = 80) :: abort_message
	1593	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
	1594
	1595
	1596	!------------------------------------------------------------------------------
	1597	! Initialisation
	1598	!------------------------------------------------------------------------------
	1599
	1600
	1601	! Few initial checksS
	1602
	1603	DO k = 1,klev
	1604	DO ind1 = 1, klon
[4910]	1605	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
	1606	zrho(ind1) = play(ind1)/temp(ind1)/rd !kg/m3
	1607	dz(ind1) = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
[3999]	1608	END DO
	1609	END DO
	1610
	1611
[4910]	1612	icefracth(:)=0.
	1613	icefracenv(:)=0.
	1614	sqrt2pi=sqrt(2.*rpi)
[3999]	1615	sqrt2=sqrt(2.)
[4910]	1616	sqrtpi=sqrt(rpi)
	1617	icefrac(:)=0.
[3999]	1618
	1619	!-------------------------------------------------------------------------------
	1620	! Identify grid points with potential mixed-phase conditions
	1621	!-------------------------------------------------------------------------------
	1622
	1623
	1624	DO ind1=1,klon
[5082]	1625	IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) &
[4910]	1626	.AND. (ind2<=klev-2) &
[5082]	1627	.AND. (frac_th(ind1)>min_frac_th_cld)) THEN
[3999]	1628	mpc_bl_points(ind1,ind2)=1
	1629	ELSE
	1630	mpc_bl_points(ind1,ind2)=0
	1631	ENDIF
	1632	ENDDO
	1633
	1634
	1635	!-------------------------------------------------------------------------------
	1636	! Thermal fraction calculation and standard deviation of the distribution
	1637	!-------------------------------------------------------------------------------
	1638
[4910]	1639	! initialisations and calculation of temperature, humidity and saturation specific humidity
[3999]	1640
[4910]	1641	DO ind1=1,klon
	1642
	1643	rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg ! kg/m2
	1644	rho(ind1) = play(ind1)/temp(ind1)/rd ! kg/m3
	1645	dz(ind1) = rhodz(ind1)/rho(ind1) ! m : epaisseur de la couche en metre
	1646	Tbefenv(ind1) = temp(ind1)
	1647	thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
	1648	Tbefth(ind1) = thetal_th(ind1)*zpspsk(ind1)
	1649	rhoth(ind1) = play(ind1)/Tbefth(ind1)/RD
	1650	qt_env(ind1) = (qt(ind1)-frac_th(ind1)qt_th(ind1))/(1.-frac_th(ind1)) !qt = aqtth + (1-a)*qtenv
[3999]	1651
[4910]	1652	ENDDO
[4072]	1653
[4910]	1654	! Calculation of saturation specific humidity
	1655	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
	1656	CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
	1657	CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
[4072]	1658
	1659
[4910]	1660	DO ind1=1,klon
[4072]	1661
	1662
[5082]	1663	IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement
[4072]	1664
[4910]	1665	! unlike in the other cloudth routine,
	1666	! We consider distributions of the saturation deficit WRT liquid
	1667	! at positive AND negative celcius temperature
	1668	! subsequently, cloud fraction corresponds to the part of the pdf corresponding
	1669	! to superstauration with respect to liquid WHATEVER the temperature
[4072]	1670
[3999]	1671	! Environment:
	1672
	1673
[4910]	1674	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
	1675	aenv=1./(1.+(alenv*RLVTT/rcpd))
	1676	senv=aenv*(qt_env(ind1)-qslenv(ind1))
[3999]	1677
	1678
	1679	! Thermals:
	1680
	1681
[4910]	1682	alth=(0.622RLVTTqslth(ind1))/(rdthetal_th(ind1)*2)
	1683	ath=1./(1.+(alth*RLVTT/rcpd))
	1684	sth=ath*(qt_th(ind1)-qslth(ind1))
[3999]	1685
	1686
[4910]	1687	! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
[3999]	1688
	1689
[4910]	1690	! Standard deviation of the distributions
[3999]	1691
[4910]	1692	sigma1s_fraca = (sigma1s_factor*0.5)(frac_th(ind1)**sigma1s_power) / &
[5087]	1693	(1-frac_th(ind1))((sth-senv)2)*0.5
[3999]	1694
	1695	IF (cloudth_ratqsmin>0.) THEN
[4910]	1696	sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
[3999]	1697	ELSE
[4910]	1698	sigma1s_ratqs = ratqs(ind1)*qt(ind1)
[3999]	1699	ENDIF
	1700
	1701	sigma1s = sigma1s_fraca + sigma1s_ratqs
[5082]	1702	IF (iflag_ratqs==11) then
[4910]	1703	sigma1s = ratqs(ind1)qt(ind1)aenv
[4623]	1704	ENDIF
[4910]	1705	sigma2s=(sigma2s_factor(((sth-senv)2)0.5)/((frac_th(ind1)+0.02)sigma2s_power))+0.002qt_th(ind1)
[3999]	1706
	1707
	1708	deltasenv=aenvvert_alphasigma1s
	1709	deltasth=athvert_alpha_thsigma2s
	1710
	1711	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1712	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1713	exp_xenv1 = exp(-1.xenv1*2)
	1714	exp_xenv2 = exp(-1.xenv2*2)
	1715	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	1716	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	1717	exp_xth1 = exp(-1.xth1*2)
	1718	exp_xth2 = exp(-1.xth2*2)
	1719
[4910]	1720	!surface CF
[3999]	1721
[4910]	1722	cth(ind1)=0.5(1.-1.erf(xth1))
	1723	cenv(ind1)=0.5(1.-1.erf(xenv1))
	1724	ctot(ind1)=frac_th(ind1)cth(ind1)+(1.-1.frac_th(ind1))*cenv(ind1)
[3999]	1725
	1726
[4910]	1727	!volume CF, condensed water and ice fraction
[3999]	1728
	1729	!environnement
	1730
	1731	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
	1732	IntJ_CF=0.5(1.-1.erf(xenv2))
	1733
[5082]	1734	IF (deltasenv < 1.e-10) THEN
[4910]	1735	qcenv(ind1)=IntJ
	1736	cenv_vol(ind1)=IntJ_CF
[3999]	1737	ELSE
	1738	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1739	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1740	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
	1741	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1742	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
[4910]	1743	qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
	1744	cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
[5082]	1745	IF (Tbefenv(ind1) < temp_nowater) THEN
[4910]	1746	! freeze all droplets in cirrus temperature regime
	1747	icefracenv(ind1)=1.
	1748	ENDIF
[3999]	1749	ENDIF
	1750
	1751
	1752
	1753	!thermals
	1754
	1755	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1756	IntJ_CF=0.5(1.-1.erf(xth2))
	1757
[5082]	1758	IF (deltasth < 1.e-10) THEN
[4910]	1759	qcth(ind1)=IntJ
	1760	cth_vol(ind1)=IntJ_CF
[3999]	1761	ELSE
	1762	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1763	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1764	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1765	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1766	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
[4910]	1767	qcth(ind1)=IntJ+IntI1+IntI2+IntI3
	1768	cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
[5082]	1769	IF (Tbefth(ind1) < temp_nowater) THEN
[4910]	1770	! freeze all droplets in cirrus temperature regime
	1771	icefracth(ind1)=1.
	1772	ENDIF
	1773
[3999]	1774	ENDIF
	1775
[4910]	1776	qctot(ind1)=frac_th(ind1)qcth(ind1)+(1.-1.frac_th(ind1))*qcenv(ind1)
	1777	ctot_vol(ind1)=frac_th(ind1)cth_vol(ind1)+(1.-1.frac_th(ind1))*cenv_vol(ind1)
[3999]	1778
[5082]	1779	IF (cenv(ind1)<min_neb_th.and.cth(ind1)<min_neb_th) THEN
[4910]	1780	ctot(ind1)=0.
	1781	ctot_vol(ind1)=0.
	1782	qcloud(ind1)=qslenv(ind1)
[3999]	1783	qincloud(ind1)=0.
[4910]	1784	icefrac(ind1)=0.
[3999]	1785	ELSE
[4910]	1786	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
	1787	qincloud(ind1)=qctot(ind1)/ctot(ind1)
[5082]	1788	IF (qctot(ind1) > 0) THEN
[4910]	1789	icefrac(ind1)=(frac_th(ind1)qcth(ind1)icefracth(ind1)+(1.-1.frac_th(ind1))qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
	1790	icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
[3999]	1791	ENDIF
	1792	ENDIF
	1793
	1794
[4910]	1795	! ELSE ! mpc_bl_points>0
[3999]	1796
	1797	ELSE ! gaussian for environment only
	1798
	1799
[4910]	1800	alenv=(0.622RLVTTqslenv(ind1))/(rdthetal_env(ind1)*2)
	1801	aenv=1./(1.+(alenv*RLVTT/rcpd))
	1802	senv=aenv*(qt_env(ind1)-qslenv(ind1))
[3999]	1803	sth=0.
[4910]	1804	sigma1s=ratqs(ind1)*qt_env(ind1)
[3999]	1805	sigma2s=0.
	1806
	1807	xenv=senv/(sqrt(2.)*sigma1s)
[4910]	1808	cenv(ind1)=0.5(1.-1.erf(xenv))
	1809	ctot(ind1)=0.5(1.+1.erf(xenv))
	1810	ctot_vol(ind1)=ctot(ind1)
	1811	qctot(ind1)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1))
	1812
	1813
[5082]	1814	IF (ctot(ind1)<min_neb_th) THEN
[4910]	1815	ctot(ind1)=0.
	1816	qcloud(ind1)=qslenv(ind1)
[3999]	1817	qincloud(ind1)=0.
	1818	ELSE
[4910]	1819	qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
	1820	qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
[3999]	1821	ENDIF
	1822
	1823
	1824	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
	1825
	1826	! Outputs used to check the PDFs
[4910]	1827	cloudth_senv(ind1) = senv
	1828	cloudth_sth(ind1) = sth
	1829	cloudth_sigmaenv(ind1) = sigma1s
	1830	cloudth_sigmath(ind1) = sigma2s
[3999]	1831
	1832
[4114]	1833	ENDDO !loop on klon
[3999]	1834
[4114]	1835
[3999]	1836	RETURN
	1837
	1838
	1839	END SUBROUTINE cloudth_mpc
	1840
[4910]	1841
	1842
	1843	! ELSE ! mpc_bl_points>0
	1844	!
	1845	! ! Treat boundary layer mixed phase clouds
	1846	!
	1847	! ! thermals
	1848	! !=========
	1849	!
	1850	! ! ice phase
	1851	! !...........
	1852	!
	1853	! qiceth=0.
	1854	! deltazlev_mpc=dz(ind1,:)
	1855	! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
	1856	! pres_mpc=pplay(ind1,:)
	1857	! fraca_mpc=fraca(ind1,:)
	1858	! snowf_mpc=snowflux(ind1,:)
	1859	! iflag_topthermals=0
	1860	! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
	1861	! iflag_topthermals = 1
	1862	! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
	1863	! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
	1864	! iflag_topthermals = 2
	1865	! ELSE
	1866	! iflag_topthermals = 0
	1867	! ENDIF
	1868	!
	1869	! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
	1870	! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
	1871	!
	1872	! ! qmax calculation
	1873	! sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
	1874	! deltasth=athlvert_alpha_thsigma2s
	1875	! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
	1876	! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
	1877	! exp_xth1 = exp(-1.xth1*2)
	1878	! exp_xth2 = exp(-1.xth2*2)
	1879	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1880	! IntJ_CF=0.5(1.-1.erf(xth2))
	1881	! IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1882	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1883	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1884	! IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1885	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1886	! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
	1887	!
	1888	!
	1889	! ! Liquid phase
	1890	! !................
	1891	! ! We account for the effect of ice crystals in thermals on sthl
	1892	! ! and on the width of the distribution
	1893	!
	1894	! sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
	1895	! + (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
	1896	!
	1897	! sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
	1898	! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
	1899	! +0.002*zqta(ind1,ind2)
	1900	! deltasthc=athlvert_alpha_thsigma2sc
	1901	!
	1902	!
	1903	! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
	1904	! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
	1905	! exp_xth1 = exp(-1.xth1*2)
	1906	! exp_xth2 = exp(-1.xth2*2)
	1907	! IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
	1908	! IntJ_CF=0.5(1.-1.erf(xth2))
	1909	! IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
	1910	! IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1911	! IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
	1912	! IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
	1913	! IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
	1914	! qliqth=IntJ+IntI1+IntI2+IntI3
	1915	!
	1916	! qlth(ind1,ind2)=MAX(0.,qliqth)
	1917	!
	1918	! ! Condensed water
	1919	!
	1920	! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
	1921	!
	1922	!
	1923	! ! consistency with subgrid distribution
	1924	!
	1925	! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
	1926	! fraci=qith(ind1,ind2)/qcth(ind1,ind2)
	1927	! qcth(ind1,ind2)=qmax
	1928	! qith(ind1,ind2)=fraci*qmax
	1929	! qlth(ind1,ind2)=(1.-fraci)*qmax
	1930	! ENDIF
	1931	!
	1932	! ! Cloud Fraction
	1933	! !...............
	1934	! ! calculation of qbase which is the value of the water vapor within mixed phase clouds
	1935	! ! such that the total water in cloud = qbase+qliqth+qiceth
	1936	! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
	1937	! ! sbase and qbase calculation (note that sbase is wrt liq so negative)
	1938	! ! look for an approximate solution with iteration
	1939	!
	1940	! ttarget=qcth(ind1,ind2)
	1941	! mini= athl*(qsith(ind1,ind2)-qslth(ind1))
	1942	! maxi= 0. !athl(3.sqrt(sigma2s))
	1943	! niter=20
	1944	! pas=(maxi-mini)/niter
	1945	! stmp=mini
	1946	! sbase=stmp
	1947	! coutref=1.E6
	1948	! DO iter=1,niter
	1949	! cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
	1950	! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
	1951	! IF (cout .LT. coutref) THEN
	1952	! sbase=stmp
	1953	! coutref=cout
	1954	! ELSE
	1955	! stmp=stmp+pas
	1956	! ENDIF
	1957	! ENDDO
	1958	! qbase=MAX(0., sbase/athl+qslth(ind1))
	1959	!
	1960	! ! surface cloud fraction in thermals
	1961	! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
	1962	! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
	1963	!
	1964	!
	1965	! !volume cloud fraction in thermals
	1966	! !to be checked
	1967	! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
	1968	! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
	1969	! exp_xth1 = exp(-1.xth1*2)
	1970	! exp_xth2 = exp(-1.xth2*2)
	1971	!
	1972	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1973	! IntJ_CF=0.5(1.-1.erf(xth2))
	1974	!
	1975	! IF (deltasth .LT. 1.e-10) THEN
	1976	! cth_vol(ind1,ind2)=IntJ_CF
	1977	! ELSE
	1978	! IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1979	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1980	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1981	! IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1982	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1983	! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	1984	! ENDIF
	1985	! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
	1986	!
	1987	!
	1988	!
	1989	! ! Environment
	1990	! !=============
	1991	! ! In the environment/downdrafts, only liquid clouds
	1992	! ! See Shupe et al. 2008, JAS
	1993	!
	1994	! ! standard deviation of the distribution in the environment
	1995	! sth=sthl
	1996	! senv=senvl
	1997	! sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
	1998	! & (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
	1999	! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
	2000	! ! in the environement
	2001	!
	2002	! sigma1s_ratqs=1E-10
	2003	! IF (cloudth_ratqsmin>0.) THEN
	2004	! sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
	2005	! ENDIF
	2006	!
	2007	! sigma1s = sigma1s_fraca + sigma1s_ratqs
	2008	! IF (iflag_ratqs.eq.11) then
	2009	! sigma1s = ratqs(ind1,ind2)po(ind1)aenv
	2010	! ENDIF
	2011	! IF (iflag_ratqs.eq.11) then
	2012	! sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
	2013	! ENDIF
	2014	! deltasenv=aenvlvert_alphasigma1s
	2015	! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
	2016	! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
	2017	! exp_xenv1 = exp(-1.xenv1*2)
	2018	! exp_xenv2 = exp(-1.xenv2*2)
	2019	!
	2020	! !surface CF
	2021	! cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
	2022	!
	2023	! !volume CF and condensed water
	2024	! IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
	2025	! IntJ_CF=0.5(1.-1.erf(xenv2))
	2026	!
	2027	! IF (deltasenv .LT. 1.e-10) THEN
	2028	! qcenv(ind1,ind2)=IntJ
	2029	! cenv_vol(ind1,ind2)=IntJ_CF
	2030	! ELSE
	2031	! IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	2032	! IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	2033	! IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
	2034	! IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	2035	! IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
	2036	! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
	2037	! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	2038	! ENDIF
	2039	!
	2040	! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
	2041	! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
	2042	!
	2043	!
	2044	!
	2045	! ! Thermals + environment
	2046	! ! =======================
	2047	! ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	2048	! qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
	2049	! ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
	2050	! IF (qcth(ind1,ind2) .GT. 0) THEN
	2051	! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
	2052	! /(fraca(ind1,ind2)*qcth(ind1,ind2) &
	2053	! +(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
	2054	! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
	2055	! ELSE
	2056	! icefrac(ind1,ind2)=0.
	2057	! ENDIF
	2058	!
	2059	! IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
	2060	! ctot(ind1,ind2)=0.
	2061	! ctot_vol(ind1,ind2)=0.
	2062	! qincloud(ind1)=0.
	2063	! qcloud(ind1)=zqsatenv(ind1)
	2064	! ELSE
	2065	! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
	2066	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
	2067	! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
	2068	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
	2069	! ENDIF
	2070	!
	2071	! ENDIF ! mpc_bl_points
	2072
	2073
	2074
[3999]	2075	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4118]	2076	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]	2077
	2078	! parameterization of ice for boundary
	2079	! layer mixed-phase clouds assuming a stationary system
	2080	!
	2081	! Note that vapor deposition on ice crystals and riming of liquid droplets
	2082	! depend on the ice number concentration Ni
	2083	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
	2084	! and use values from Hong et al. 2004, MWR for instance
	2085	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
	2086	! One could also think of a more complex expression of Ni;
	2087	! function of qi, T, the concentration in aerosols or INP ..
	2088	! Here we prefer fixing Ni to a tuning parameter
	2089	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
	2090	! in Mioche et al. 2017
	2091	!
	2092	!
	2093	! References:
	2094	!------------
	2095	! This parameterization is thoroughly described in Vignon et al.
	2096	!
	2097	! More specifically
	2098	! for the Water vapor deposition process:
	2099	!
	2100	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
	2101	! ment of stratiform cloudfs and precipitation in large-scale
	2102	! models. I: Description and evaluation of the microphysical
	2103	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
	2104	!
	2105	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
	2106	! stratiform cloud microphysics scheme in the NCAR Com-
	2107	! munity Atmosphere Model (CAM3). Part I: Description and
	2108	! numerical tests. J. Climate, 21, 3642–3659
	2109	!
	2110	! for the Riming process:
	2111	!
	2112	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
	2113	! scale structure and organization of clouds and precipitation in
	2114	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
	2115	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
	2116	!
	2117	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
	2118	! forecasts of winter precipitation using an improved bulk
	2119	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
	2120	! forecasts of winter precipitation using an improved bulk
	2121	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
	2122	!
	2123	! For the formation of clouds by thermals:
	2124	!
	2125	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
	2126	! the Atmospheric Sciences, 65, 407–425.
	2127	!
	2128	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
	2129	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
	2130	!
	2131	!
	2132	!
	2133	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
	2134	!=============================================================================
	2135
[4651]	2136	use phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
[3999]	2137
	2138	IMPLICIT none
	2139
	2140	INCLUDE "YOMCST.h"
	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
[5082]	2180	IF (iflag_topthermals > 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.
[5082]	2233	IF (fraca(ind2p1) > 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
[5082]	2242	IF (fm_therm(ind1,ind2p1) > 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