Context Navigation

cloudth_mod.F90 @ 2912

Last change on this file since 2912 was 2911, checked in by jbmadeleine, 7 years ago

Changes in the vertical PDF used for shallow convective clouds.
Added vert_alpha_th which is the variance of the vertical PDF used for thermals.
If nothing is specified in the def files, it is set equal to vert_alpha.
Otherwise, it is now possible to have a different vertical PDF in the thermals
and in the environment.
Also added a fourth option iflag_cloudth_vert = 4 where
deltas=vert_alpha*sigmas instead of deltas=a*vert_alpha*sigmas.

File size: 38.9 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, &
	589	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	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
	626
	627	REAL sigma1(ngrid,klev)
	628	REAL sigma2(ngrid,klev)
	629	REAL qlth(ngrid,klev)
	630	REAL qlenv(ngrid,klev)
	631	REAL qltot(ngrid,klev)
	632	REAL cth(ngrid,klev)
	633	REAL cenv(ngrid,klev)
	634	REAL ctot(ngrid,klev)
	635	REAL rneb(ngrid,klev)
	636	REAL t(ngrid,klev)
	637	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
	638	REAL rdd,cppd,Lv
	639	REAL alth,alenv,ath,aenv
	640	REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
	641	REAL Tbef,zdelta,qsatbef,zcor
	642	REAL qlbef
	643	REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
	644	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	645	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	646	REAL zqs(ngrid), qcloud(ngrid)
	647	REAL erf
	648
	649
	650
	651	IF (iflag_cloudth_vert.GE.1) THEN
	652	CALL cloudth_vert_v3(ngrid,klev,ind2, &
	653	& ztv,po,zqta,fraca, &
	654	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	655	& ratqs,zqs,t)
	656	RETURN
	657	ENDIF
	658	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	659
	660
	661	!-------------------------------------------------------------------------------
	662	! Initialisation des variables r?elles
	663	!-------------------------------------------------------------------------------
	664	sigma1(:,:)=0.
	665	sigma2(:,:)=0.
	666	qlth(:,:)=0.
	667	qlenv(:,:)=0.
	668	qltot(:,:)=0.
	669	rneb(:,:)=0.
	670	qcloud(:)=0.
	671	cth(:,:)=0.
	672	cenv(:,:)=0.
	673	ctot(:,:)=0.
	674	qsatmmussig1=0.
	675	qsatmmussig2=0.
	676	rdd=287.04
	677	cppd=1005.7
	678	pi=3.14159
	679	Lv=2.5e6
	680	sqrt2pi=sqrt(2.*pi)
	681	sqrt2=sqrt(2.)
	682	sqrtpi=sqrt(pi)
	683
	684
	685	!-------------------------------------------------------------------------------
	686	! Cloud fraction in the thermals and standard deviation of the PDFs
	687	!-------------------------------------------------------------------------------
	688	do ind1=1,ngrid
	689
	690	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
	691
	692	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	693
	694	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	695	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	696	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	697	qsatbef=MIN(0.5,qsatbef)
	698	zcor=1./(1.-retv*qsatbef)
	699	qsatbef=qsatbef*zcor
	700	zqsatenv(ind1,ind2)=qsatbef
	701
	702
	703	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	704	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	705	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	706
	707	!po = qt de l'environnement ET des thermique
	708	!zqenv = qt environnement
	709	!zqsatenv = qsat environnement
	710	!zthl = Tl environnement
	711
	712
	713	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	714	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	715	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	716	qsatbef=MIN(0.5,qsatbef)
	717	zcor=1./(1.-retv*qsatbef)
	718	qsatbef=qsatbef*zcor
	719	zqsatth(ind1,ind2)=qsatbef
	720
	721	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	722	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	723	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	724
	725	!zqta = qt thermals
	726	!zqsatth = qsat thermals
	727	!ztla = Tl thermals
	728
	729	!------------------------------------------------------------------------------
	730	! s standard deviations
	731	!------------------------------------------------------------------------------
	732
	733	! tests
	734	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	735	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
	736	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	737	! final option
	738	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	739	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	740
	741	!------------------------------------------------------------------------------
	742	! Condensed water and cloud cover
	743	!------------------------------------------------------------------------------
	744	xth=sth/(sqrt2*sigma2s)
	745	xenv=senv/(sqrt2*sigma1s)
	746	cth(ind1,ind2)=0.5(1.+1.erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	747	cenv(ind1,ind2)=0.5(1.+1.erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	748	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	749
	750	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt2*cth(ind1,ind2))
	751	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	752	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	753
	754	if (ctot(ind1,ind2).lt.1.e-10) then
	755	ctot(ind1,ind2)=0.
	756	qcloud(ind1)=zqsatenv(ind1,ind2)
	757	else
	758	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	759	endif
	760
	761	else ! Environnement only, follow the if l.110
	762
	763	zqenv(ind1)=po(ind1)
	764	Tbef=t(ind1,ind2)
	765	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	766	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	767	qsatbef=MIN(0.5,qsatbef)
	768	zcor=1./(1.-retv*qsatbef)
	769	qsatbef=qsatbef*zcor
	770	zqsatenv(ind1,ind2)=qsatbef
	771
	772	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	773	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	774	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	775	aenv=1./(1.+(alenv*Lv/cppd))
	776	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	777
	778	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	779
	780	xenv=senv/(sqrt2*sigma1s)
	781	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	782	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	783
	784	if (ctot(ind1,ind2).lt.1.e-3) then
	785	ctot(ind1,ind2)=0.
	786	qcloud(ind1)=zqsatenv(ind1,ind2)
	787	else
	788	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	789	endif
	790
	791
	792	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
	793	enddo ! from the loop on ngrid l.108
	794	return
	795	! end
	796	END SUBROUTINE cloudth_v3
	797
	798
	799
	800	!===========================================================================
	801	SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
	802	& ztv,po,zqta,fraca, &
	803	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	804	& ratqs,zqs,t)
	805
	806	!===========================================================================
	807	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	808	! Date : 25 Mai 2010
	809	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	810	!===========================================================================
	811
	812
	813	USE ioipsl_getin_p_mod, ONLY : getin_p
	814
	815	IMPLICIT NONE
	816
	817	#include "YOMCST.h"
	818	#include "YOETHF.h"
	819	#include "FCTTRE.h"
	820	#include "thermcell.h"
	821	#include "nuage.h"
	822
	823	INTEGER itap,ind1,ind2
	824	INTEGER ngrid,klev,klon,l,ig
	825
	826	REAL ztv(ngrid,klev)
	827	REAL po(ngrid)
	828	REAL zqenv(ngrid)
	829	REAL zqta(ngrid,klev)
	830
	831	REAL fraca(ngrid,klev+1)
	832	REAL zpspsk(ngrid,klev)
	833	REAL paprs(ngrid,klev+1)
	834	REAL ztla(ngrid,klev)
	835	REAL zthl(ngrid,klev)
	836
	837	REAL zqsatth(ngrid,klev)
	838	REAL zqsatenv(ngrid,klev)
	839
	840
	841	REAL sigma1(ngrid,klev)
	842	REAL sigma2(ngrid,klev)
	843	REAL qlth(ngrid,klev)
	844	REAL qlenv(ngrid,klev)
	845	REAL qltot(ngrid,klev)
	846	REAL cth(ngrid,klev)
	847	REAL cenv(ngrid,klev)
	848	REAL ctot(ngrid,klev)
	849	REAL rneb(ngrid,klev)
	850	REAL t(ngrid,klev)
	851	REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
	852	REAL rdd,cppd,Lv
	853	REAL alth,alenv,ath,aenv
	854	REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
	855	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
	856	REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
	857	REAL Tbef,zdelta,qsatbef,zcor
	858	REAL qlbef
	859	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	860	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	861	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2911]	862	REAL,SAVE :: vert_alpha, vert_alpha_th
	863	!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
[2686]	864	LOGICAL, SAVE :: firstcall = .TRUE.
[2910]	865	!$OMP THREADPRIVATE(firstcall)
[2686]	866
	867	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	868	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	869	REAL zqs(ngrid), qcloud(ngrid)
	870	REAL erf
	871
	872	!------------------------------------------------------------------------------
	873	! Initialize
	874	!------------------------------------------------------------------------------
	875	sigma1(:,:)=0.
	876	sigma2(:,:)=0.
	877	qlth(:,:)=0.
	878	qlenv(:,:)=0.
	879	qltot(:,:)=0.
	880	rneb(:,:)=0.
	881	qcloud(:)=0.
	882	cth(:,:)=0.
	883	cenv(:,:)=0.
	884	ctot(:,:)=0.
	885	qsatmmussig1=0.
	886	qsatmmussig2=0.
	887	rdd=287.04
	888	cppd=1005.7
	889	pi=3.14159
	890	Lv=2.5e6
	891	sqrt2pi=sqrt(2.*pi)
	892	sqrt2=sqrt(2.)
	893	sqrtpi=sqrt(pi)
	894
	895	IF (firstcall) THEN
	896	vert_alpha=0.5
	897	CALL getin_p('cloudth_vert_alpha',vert_alpha)
	898	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
[2911]	899	! The factor used for the thermal is equal to that of the environment
	900	! if nothing is explicitly specified in the def file
	901	vert_alpha_th=vert_alpha
	902	CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
	903	WRITE(,) 'cloudth_vert_alpha_th = ', vert_alpha_th
[2686]	904	firstcall=.FALSE.
	905	ENDIF
	906
	907	!-------------------------------------------------------------------------------
	908	! Calcul de la fraction du thermique et des ecart-types des distributions
	909	!-------------------------------------------------------------------------------
	910	do ind1=1,ngrid
	911
	912	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
	913
	914	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = aqtth + (1-a)*qtenv
	915
	916
	917	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	918	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	919	qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	920	qsatbef=MIN(0.5,qsatbef)
	921	zcor=1./(1.-retv*qsatbef)
	922	qsatbef=qsatbef*zcor
	923	zqsatenv(ind1,ind2)=qsatbef
	924
	925
	926	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2) !qsl, p84
	927	aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
	928	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
	929
	930	!zqenv = qt environnement
	931	!zqsatenv = qsat environnement
	932	!zthl = Tl environnement
	933
	934
	935	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	936	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	937	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	938	qsatbef=MIN(0.5,qsatbef)
	939	zcor=1./(1.-retv*qsatbef)
	940	qsatbef=qsatbef*zcor
	941	zqsatth(ind1,ind2)=qsatbef
	942
	943	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2) !qsl, p84
	944	ath=1./(1.+(alth*Lv/cppd)) !al, p84
	945	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
	946
	947
	948	!zqta = qt thermals
	949	!zqsatth = qsat thermals
	950	!ztla = Tl thermals
	951
	952	!------------------------------------------------------------------------------
	953	! s standard deviation
	954	!------------------------------------------------------------------------------
	955
	956	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	957	sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	958	! tests
	959	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	960	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
	961	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	962	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
	963	! if (paprs(ind1,ind2).gt.90000) then
	964	! ratqs(ind1,ind2)=0.002
	965	! else
	966	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	967	! endif
	968	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	969	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	970	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	971	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	972
	973	IF (iflag_cloudth_vert == 1) THEN
	974	!-------------------------------------------------------------------------------
	975	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
	976	!-------------------------------------------------------------------------------
	977
	978	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	979	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	980
	981	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
	982	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
	983	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
	984	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
	985	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	986	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	987
	988	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
	989	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
	990	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	991
	992	! Environment
	993	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
	994	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	995	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
	996	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
	997
	998	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	999
	1000	! Thermal
	1001	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
	1002	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	1003	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
	1004	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
	1005	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1006	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	1007
[2911]	1008	ELSE IF (iflag_cloudth_vert >= 3) THEN
[2686]	1009
	1010	!-------------------------------------------------------------------------------
	1011	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
	1012	!-------------------------------------------------------------------------------
	1013	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	1014	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	1015	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	1016	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
[2911]	1017	IF (iflag_cloudth_vert == 3) THEN
	1018	deltasenv=aenvvert_alphasigma1s
	1019	deltasth=athvert_alpha_thsigma2s
	1020	ELSE IF (iflag_cloudth_vert == 4) THEN
	1021	deltasenv=vert_alpha*sigma1s
	1022	deltasth=vert_alpha_th*sigma2s
	1023	ENDIF
[2686]	1024
	1025	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	1026	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	1027	exp_xenv1 = exp(-1.xenv1*2)
	1028	exp_xenv2 = exp(-1.xenv2*2)
	1029	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	1030	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	1031	exp_xth1 = exp(-1.xth1*2)
	1032	exp_xth2 = exp(-1.xth2*2)
	1033
	1034	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	1035	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
	1036	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1037
	1038
	1039	!environnement
	1040	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
	1041	if (deltasenv .lt. 1.e-10) then
	1042	qlenv(ind1,ind2)=IntJ
	1043	else
	1044	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1045	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1046	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
	1047	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1048	endif
	1049
	1050
	1051	!thermique
	1052	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1053	if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
	1054	qlth(ind1,ind2)=IntJ
	1055	else
	1056	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1057	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1058	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1059	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	1060	endif
	1061
	1062	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	1063
	1064
[2911]	1065	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1066
	1067	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
	1068	ctot(ind1,ind2)=0.
	1069	qcloud(ind1)=zqsatenv(ind1,ind2)
	1070
	1071	else
	1072
	1073	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	1074	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	1075	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	1076
	1077	endif
	1078
	1079	else ! Environment only
	1080
	1081	zqenv(ind1)=po(ind1)
	1082	Tbef=t(ind1,ind2)
	1083	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	1084	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	1085	qsatbef=MIN(0.5,qsatbef)
	1086	zcor=1./(1.-retv*qsatbef)
	1087	qsatbef=qsatbef*zcor
	1088	zqsatenv(ind1,ind2)=qsatbef
	1089
	1090
	1091	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	1092	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	1093	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	1094	aenv=1./(1.+(alenv*Lv/cppd))
	1095	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	1096
	1097
	1098	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	1099
	1100	xenv=senv/(sqrt2*sigma1s)
	1101	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	1102	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt2*cenv(ind1,ind2))
	1103
	1104	if (ctot(ind1,ind2).lt.1.e-3) then
	1105	ctot(ind1,ind2)=0.
	1106	qcloud(ind1)=zqsatenv(ind1,ind2)
	1107
	1108	else
	1109
	1110	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	1111
	1112	endif
	1113
	1114	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
	1115	enddo ! from the loop on ngrid l.333
	1116
	1117	return
	1118	! end
	1119	END SUBROUTINE cloudth_vert_v3
	1120	!
	1121	END MODULE cloudth_mod

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

Download in other formats:

Original Format