Context Navigation

lmdz_cloudth.f90 @ 5408

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

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

Download in other formats:

Original Format