Context Navigation

cloudth_mod.F90 @ 5454

Last change on this file since 5454 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
File size: 92.5 KB

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

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

Download in other formats:

Original Format