Context Navigation

cloudth_mod.F90 @ 4122

Last change on this file since 4122 was 2960, checked in by fhourdin, 7 years ago
Un lievre leve, sans doute. Un tout petit seuil sur les largeurs de distributions sous maille du schéma de nuage, qui produit un gros effet sur le contraste strato-cumulus / cumulus ...
File size: 42.1 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, &
	9	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	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 "thermcell.h"
	27	#include "nuage.h"
	28
	29	INTEGER itap,ind1,ind2
	30	INTEGER ngrid,klev,klon,l,ig
	31
	32	REAL ztv(ngrid,klev)
	33	REAL po(ngrid)
	34	REAL zqenv(ngrid)
	35	REAL zqta(ngrid,klev)
	36
	37	REAL fraca(ngrid,klev+1)
	38	REAL zpspsk(ngrid,klev)
	39	REAL paprs(ngrid,klev+1)
	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, &
	80	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	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, &
	253	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	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 "thermcell.h"
	271	#include "nuage.h"
	272
	273	INTEGER itap,ind1,ind2
	274	INTEGER ngrid,klev,klon,l,ig
	275
	276	REAL ztv(ngrid,klev)
	277	REAL po(ngrid)
	278	REAL zqenv(ngrid)
	279	REAL zqta(ngrid,klev)
	280
	281	REAL fraca(ngrid,klev+1)
	282	REAL zpspsk(ngrid,klev)
	283	REAL paprs(ngrid,klev+1)
	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
	587	SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
	588	& ztv,po,zqta,fraca, &
[2945]	589	& qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
[2686]	590	& ratqs,zqs,t)
	591
	592
	593	IMPLICIT NONE
	594
	595
	596	!===========================================================================
	597	! Author : Arnaud Octavio Jam (LMD/CNRS)
	598	! Date : 25 Mai 2010
	599	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	600	!===========================================================================
	601
	602
	603	#include "YOMCST.h"
	604	#include "YOETHF.h"
	605	#include "FCTTRE.h"
	606	#include "thermcell.h"
	607	#include "nuage.h"
	608
	609	INTEGER itap,ind1,ind2
	610	INTEGER ngrid,klev,klon,l,ig
	611
	612	REAL ztv(ngrid,klev)
	613	REAL po(ngrid)
	614	REAL zqenv(ngrid)
	615	REAL zqta(ngrid,klev)
	616
	617	REAL fraca(ngrid,klev+1)
	618	REAL zpspsk(ngrid,klev)
	619	REAL paprs(ngrid,klev+1)
	620	REAL ztla(ngrid,klev)
	621	REAL zthl(ngrid,klev)
	622
	623	REAL zqsatth(ngrid,klev)
	624	REAL zqsatenv(ngrid,klev)
	625
[2945]	626	REAL sigma1(ngrid,klev)
[2686]	627	REAL sigma2(ngrid,klev)
	628	REAL qlth(ngrid,klev)
	629	REAL qlenv(ngrid,klev)
	630	REAL qltot(ngrid,klev)
[2945]	631	REAL cth(ngrid,klev)
[2686]	632	REAL cenv(ngrid,klev)
	633	REAL ctot(ngrid,klev)
[2945]	634	REAL cth_vol(ngrid,klev)
	635	REAL cenv_vol(ngrid,klev)
	636	REAL ctot_vol(ngrid,klev)
	637	REAL rneb(ngrid,klev)
[2686]	638	REAL t(ngrid,klev)
	639	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
	640	REAL rdd,cppd,Lv
	641	REAL alth,alenv,ath,aenv
	642	REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
	643	REAL Tbef,zdelta,qsatbef,zcor
	644	REAL qlbef
	645	REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
	646	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	647	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	648	REAL zqs(ngrid), qcloud(ngrid)
	649	REAL erf
	650
	651
	652
	653	IF (iflag_cloudth_vert.GE.1) THEN
	654	CALL cloudth_vert_v3(ngrid,klev,ind2, &
	655	& ztv,po,zqta,fraca, &
[2945]	656	& qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
[2686]	657	& ratqs,zqs,t)
	658	RETURN
	659	ENDIF
	660	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	661
	662
	663	!-------------------------------------------------------------------------------
	664	! Initialisation des variables r?elles
	665	!-------------------------------------------------------------------------------
	666	sigma1(:,:)=0.
	667	sigma2(:,:)=0.
	668	qlth(:,:)=0.
	669	qlenv(:,:)=0.
	670	qltot(:,:)=0.
	671	rneb(:,:)=0.
	672	qcloud(:)=0.
	673	cth(:,:)=0.
	674	cenv(:,:)=0.
	675	ctot(:,:)=0.
[2945]	676	cth_vol(:,:)=0.
	677	cenv_vol(:,:)=0.
	678	ctot_vol(:,:)=0.
[2686]	679	qsatmmussig1=0.
	680	qsatmmussig2=0.
	681	rdd=287.04
	682	cppd=1005.7
	683	pi=3.14159
	684	Lv=2.5e6
	685	sqrt2pi=sqrt(2.*pi)
	686	sqrt2=sqrt(2.)
	687	sqrtpi=sqrt(pi)
	688
	689
	690	!-------------------------------------------------------------------------------
	691	! Cloud fraction in the thermals and standard deviation of the PDFs
	692	!-------------------------------------------------------------------------------
	693	do ind1=1,ngrid
	694
	695	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
	696
	697	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	698
	699	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	700	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	701	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	702	qsatbef=MIN(0.5,qsatbef)
	703	zcor=1./(1.-retv*qsatbef)
	704	qsatbef=qsatbef*zcor
	705	zqsatenv(ind1,ind2)=qsatbef
	706
	707
	708	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	709	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	710	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	711
	712	!po = qt de l'environnement ET des thermique
	713	!zqenv = qt environnement
	714	!zqsatenv = qsat environnement
	715	!zthl = Tl environnement
	716
	717
	718	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	719	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	720	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	721	qsatbef=MIN(0.5,qsatbef)
	722	zcor=1./(1.-retv*qsatbef)
	723	qsatbef=qsatbef*zcor
	724	zqsatth(ind1,ind2)=qsatbef
	725
	726	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	727	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	728	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	729
	730	!zqta = qt thermals
	731	!zqsatth = qsat thermals
	732	!ztla = Tl thermals
	733
	734	!------------------------------------------------------------------------------
	735	! s standard deviations
	736	!------------------------------------------------------------------------------
	737
	738	! tests
	739	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	740	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
	741	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	742	! final option
	743	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	744	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	745
	746	!------------------------------------------------------------------------------
	747	! Condensed water and cloud cover
	748	!------------------------------------------------------------------------------
	749	xth=sth/(sqrt2*sigma2s)
	750	xenv=senv/(sqrt2*sigma1s)
	751	cth(ind1,ind2)=0.5(1.+1.erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	752	cenv(ind1,ind2)=0.5(1.+1.erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	753	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[2945]	754	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	755
	756	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth(ind1,ind2))
	757	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	758	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	759
	760	if (ctot(ind1,ind2).lt.1.e-10) then
	761	ctot(ind1,ind2)=0.
	762	qcloud(ind1)=zqsatenv(ind1,ind2)
	763	else
	764	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	765	endif
	766
	767	else ! Environnement only, follow the if l.110
	768
	769	zqenv(ind1)=po(ind1)
	770	Tbef=t(ind1,ind2)
	771	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	772	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	773	qsatbef=MIN(0.5,qsatbef)
	774	zcor=1./(1.-retv*qsatbef)
	775	qsatbef=qsatbef*zcor
	776	zqsatenv(ind1,ind2)=qsatbef
	777
	778	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	779	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	780	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	781	aenv=1./(1.+(alenv*Lv/cppd))
	782	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	783
	784	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	785
	786	xenv=senv/(sqrt2*sigma1s)
	787	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
[2945]	788	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	789	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	790
	791	if (ctot(ind1,ind2).lt.1.e-3) then
	792	ctot(ind1,ind2)=0.
	793	qcloud(ind1)=zqsatenv(ind1,ind2)
	794	else
	795	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	796	endif
	797
	798
	799	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
	800	enddo ! from the loop on ngrid l.108
	801	return
	802	! end
	803	END SUBROUTINE cloudth_v3
	804
	805
	806
	807	!===========================================================================
	808	SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
	809	& ztv,po,zqta,fraca, &
[2945]	810	& qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
[2686]	811	& ratqs,zqs,t)
	812
	813	!===========================================================================
	814	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	815	! Date : 25 Mai 2010
	816	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	817	!===========================================================================
	818
	819
	820	USE ioipsl_getin_p_mod, ONLY : getin_p
[2958]	821	USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
	822	& cloudth_sigmath,cloudth_sigmaenv
[2686]	823
	824	IMPLICIT NONE
	825
	826	#include "YOMCST.h"
	827	#include "YOETHF.h"
	828	#include "FCTTRE.h"
	829	#include "thermcell.h"
	830	#include "nuage.h"
	831
	832	INTEGER itap,ind1,ind2
	833	INTEGER ngrid,klev,klon,l,ig
	834
	835	REAL ztv(ngrid,klev)
	836	REAL po(ngrid)
	837	REAL zqenv(ngrid)
	838	REAL zqta(ngrid,klev)
	839
	840	REAL fraca(ngrid,klev+1)
	841	REAL zpspsk(ngrid,klev)
	842	REAL paprs(ngrid,klev+1)
	843	REAL ztla(ngrid,klev)
	844	REAL zthl(ngrid,klev)
	845
	846	REAL zqsatth(ngrid,klev)
	847	REAL zqsatenv(ngrid,klev)
	848
	849	REAL sigma1(ngrid,klev)
	850	REAL sigma2(ngrid,klev)
	851	REAL qlth(ngrid,klev)
	852	REAL qlenv(ngrid,klev)
	853	REAL qltot(ngrid,klev)
[2945]	854	REAL cth(ngrid,klev)
[2686]	855	REAL cenv(ngrid,klev)
	856	REAL ctot(ngrid,klev)
[2945]	857	REAL cth_vol(ngrid,klev)
	858	REAL cenv_vol(ngrid,klev)
	859	REAL ctot_vol(ngrid,klev)
[2686]	860	REAL rneb(ngrid,klev)
	861	REAL t(ngrid,klev)
	862	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
	863	REAL rdd,cppd,Lv
	864	REAL alth,alenv,ath,aenv
[2957]	865	REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
	866	REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
[2686]	867	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
[2945]	868	REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
[2686]	869	REAL Tbef,zdelta,qsatbef,zcor
	870	REAL qlbef
	871	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	872	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	873	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2911]	874	REAL,SAVE :: vert_alpha, vert_alpha_th
	875	!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
[2957]	876	REAL,SAVE :: sigma1s_factor=1.1
	877	REAL,SAVE :: sigma1s_power=0.6
[2960]	878	REAL,SAVE :: cloudth_ratqsmin=-1.
	879	!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,cloudth_ratqsmin)
[2957]	880	INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
	881	!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
	882
[2686]	883	LOGICAL, SAVE :: firstcall = .TRUE.
[2910]	884	!$OMP THREADPRIVATE(firstcall)
[2686]	885
	886	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	887	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	888	REAL zqs(ngrid), qcloud(ngrid)
	889	REAL erf
	890
	891	!------------------------------------------------------------------------------
	892	! Initialize
	893	!------------------------------------------------------------------------------
	894	sigma1(:,:)=0.
	895	sigma2(:,:)=0.
	896	qlth(:,:)=0.
	897	qlenv(:,:)=0.
	898	qltot(:,:)=0.
	899	rneb(:,:)=0.
	900	qcloud(:)=0.
	901	cth(:,:)=0.
	902	cenv(:,:)=0.
	903	ctot(:,:)=0.
[2945]	904	cth_vol(:,:)=0.
	905	cenv_vol(:,:)=0.
	906	ctot_vol(:,:)=0.
[2686]	907	qsatmmussig1=0.
	908	qsatmmussig2=0.
	909	rdd=287.04
	910	cppd=1005.7
	911	pi=3.14159
	912	Lv=2.5e6
	913	sqrt2pi=sqrt(2.*pi)
	914	sqrt2=sqrt(2.)
	915	sqrtpi=sqrt(pi)
	916
	917	IF (firstcall) THEN
	918	vert_alpha=0.5
	919	CALL getin_p('cloudth_vert_alpha',vert_alpha)
	920	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
[2911]	921	! The factor used for the thermal is equal to that of the environment
	922	! if nothing is explicitly specified in the def file
	923	vert_alpha_th=vert_alpha
	924	CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
	925	WRITE(,) 'cloudth_vert_alpha_th = ', vert_alpha_th
[2957]	926	! Factor used in the calculation of sigma1s
	927	CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
	928	WRITE(,) 'cloudth_sigma1s_factor = ', sigma1s_factor
	929	! Power used in the calculation of sigma1s
	930	CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
	931	WRITE(,) 'cloudth_sigma1s_power = ', sigma1s_power
[2960]	932	! Minimum value for the environmental air subgrid water distrib
	933	CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
	934	WRITE(,) 'cloudth_ratqsmin = ', cloudth_ratqsmin
[2957]	935	! Remove the dependency to ratqs from the variance of the vertical PDF
	936	CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
	937	WRITE(,) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
	938
[2686]	939	firstcall=.FALSE.
	940	ENDIF
	941
	942	!-------------------------------------------------------------------------------
	943	! Calcul de la fraction du thermique et des ecart-types des distributions
	944	!-------------------------------------------------------------------------------
	945	do ind1=1,ngrid
	946
	947	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
	948
	949	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
	950
	951
	952	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	953	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	954	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	955	qsatbef=MIN(0.5,qsatbef)
	956	zcor=1./(1.-retv*qsatbef)
	957	qsatbef=qsatbef*zcor
	958	zqsatenv(ind1,ind2)=qsatbef
	959
	960
	961	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	962	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	963	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	964
	965	!zqenv = qt environnement
	966	!zqsatenv = qsat environnement
	967	!zthl = Tl environnement
	968
	969
	970	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	971	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	972	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	973	qsatbef=MIN(0.5,qsatbef)
	974	zcor=1./(1.-retv*qsatbef)
	975	qsatbef=qsatbef*zcor
	976	zqsatth(ind1,ind2)=qsatbef
	977
	978	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	979	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	980	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	981
	982
	983	!zqta = qt thermals
	984	!zqsatth = qsat thermals
	985	!ztla = Tl thermals
	986
	987	!------------------------------------------------------------------------------
	988	! s standard deviation
	989	!------------------------------------------------------------------------------
	990
[2957]	991	sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
	992	& (1-fraca(ind1,ind2))((sth-senv)2)*0.5
	993	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
[2960]	994	IF (cloudth_ratqsmin>0.) THEN
	995	sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
	996	ELSE
	997	sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
	998	ENDIF
[2957]	999	sigma1s = sigma1s_fraca + sigma1s_ratqs
[2686]	1000	sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	1001	! tests
	1002	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	1003	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
	1004	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	1005	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
	1006	! if (paprs(ind1,ind2).gt.90000) then
	1007	! ratqs(ind1,ind2)=0.002
	1008	! else
	1009	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	1010	! endif
	1011	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	1012	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	1013	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	1014	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	1015
	1016	IF (iflag_cloudth_vert == 1) THEN
	1017	!-------------------------------------------------------------------------------
	1018	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
	1019	!-------------------------------------------------------------------------------
	1020
	1021	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	1022	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	1023
	1024	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1025	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1026	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
	1027	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
	1028	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	1029	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	1030
	1031	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
	1032	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
	1033	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1034
	1035	! Environment
	1036	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
	1037	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	1038	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
	1039	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
	1040
	1041	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1042
	1043	! Thermal
	1044	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
	1045	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	1046	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
	1047	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
	1048	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1049	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	1050
[2911]	1051	ELSE IF (iflag_cloudth_vert >= 3) THEN
[2686]	1052
	1053	!-------------------------------------------------------------------------------
	1054	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
	1055	!-------------------------------------------------------------------------------
	1056	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	1057	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	1058	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	1059	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
[2911]	1060	IF (iflag_cloudth_vert == 3) THEN
	1061	deltasenv=aenvvert_alphasigma1s
	1062	deltasth=athvert_alpha_thsigma2s
	1063	ELSE IF (iflag_cloudth_vert == 4) THEN
[2957]	1064	IF (iflag_cloudth_vert_noratqs == 1) THEN
[2959]	1065	deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
[2957]	1066	deltasth=vert_alpha_th*sigma2s
	1067	ELSE
	1068	deltasenv=vert_alpha*sigma1s
	1069	deltasth=vert_alpha_th*sigma2s
	1070	ENDIF
[2911]	1071	ENDIF
[2686]	1072
	1073	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1074	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1075	exp_xenv1 = exp(-1.xenv1*2)
	1076	exp_xenv2 = exp(-1.xenv2*2)
	1077	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	1078	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	1079	exp_xth1 = exp(-1.xth1*2)
	1080	exp_xth2 = exp(-1.xth2*2)
	1081
[2945]	1082	!CF_surfacique
[2686]	1083	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	1084	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
[2945]	1085	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[2686]	1086
[2945]	1087
	1088	!CF_volumique & eau condense
[2686]	1089	!environnement
	1090	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
[2945]	1091	IntJ_CF=0.5(1.-1.erf(xenv2))
[2686]	1092	if (deltasenv .lt. 1.e-10) then
	1093	qlenv(ind1,ind2)=IntJ
[2945]	1094	cenv_vol(ind1,ind2)=IntJ_CF
[2686]	1095	else
	1096	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1097	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1098	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
[2945]	1099	IntI1_CF=((senv+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1100	IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
[2686]	1101	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1102	cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1103	endif
	1104
	1105	!thermique
	1106	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
[2945]	1107	IntJ_CF=0.5(1.-1.erf(xth2))
	1108	if (deltasth .lt. 1.e-10) then
[2686]	1109	qlth(ind1,ind2)=IntJ
[2945]	1110	cth_vol(ind1,ind2)=IntJ_CF
[2686]	1111	else
	1112	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1113	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1114	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
[2945]	1115	IntI1_CF=((sth+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1116	IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
[2686]	1117	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
[2945]	1118	cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
[2686]	1119	endif
	1120
[2945]	1121
[2686]	1122	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
[2945]	1123	ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
[2686]	1124
[2911]	1125	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1126
[2945]	1127
[2686]	1128	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
	1129	ctot(ind1,ind2)=0.
[2945]	1130	ctot_vol(ind1,ind2)=0.
[2686]	1131	qcloud(ind1)=zqsatenv(ind1,ind2)
	1132
	1133	else
	1134
	1135	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	1136	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	1137	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	1138
	1139	endif
	1140
	1141	else ! Environment only
	1142
	1143	zqenv(ind1)=po(ind1)
	1144	Tbef=t(ind1,ind2)
	1145	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	1146	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	1147	qsatbef=MIN(0.5,qsatbef)
	1148	zcor=1./(1.-retv*qsatbef)
	1149	qsatbef=qsatbef*zcor
	1150	zqsatenv(ind1,ind2)=qsatbef
	1151
	1152
	1153	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	1154	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	1155	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	1156	aenv=1./(1.+(alenv*Lv/cppd))
[2958]	1157	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	1158	sth=0.
[2686]	1159
	1160	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[2958]	1161	sigma2s=0.
[2686]	1162
	1163	xenv=senv/(sqrt2*sigma1s)
	1164	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
[2945]	1165	ctot_vol(ind1,ind2)=ctot(ind1,ind2)
[2686]	1166	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	1167
	1168	if (ctot(ind1,ind2).lt.1.e-3) then
	1169	ctot(ind1,ind2)=0.
	1170	qcloud(ind1)=zqsatenv(ind1,ind2)
	1171
	1172	else
	1173
	1174	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	1175
	1176	endif
	1177
	1178	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]	1179	! Outputs used to check the PDFs
	1180	cloudth_senv(ind1,ind2) = senv
	1181	cloudth_sth(ind1,ind2) = sth
	1182	cloudth_sigmaenv(ind1,ind2) = sigma1s
	1183	cloudth_sigmath(ind1,ind2) = sigma2s
	1184
[2686]	1185	enddo ! from the loop on ngrid l.333
	1186
	1187	return
	1188	! end
	1189	END SUBROUTINE cloudth_vert_v3
	1190	!
	1191	END MODULE cloudth_mod

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

Download in other formats:

Original Format