Context Navigation

lmdz_cloudth.F90 @ 5136

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

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

Download in other formats:

Original Format