Context Navigation

cloudth.F90 @ 2602

Last change on this file since 2602 was 2586, checked in by jbmadeleine, 8 years ago
Added a parameter called cloudth_vert_alpha defined in physiq.def; it can be used to change the width of the PDF used for the cloud vertical subgrid scale heterogeneity.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 19.5 KB

Rev	Line
[1399]	1	SUBROUTINE cloudth(ngrid,klev,ind2, &
	2	& ztv,po,zqta,fraca, &
	3	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	4	& ratqs,zqs,t)
	5
	6
	7	IMPLICIT NONE
	8
	9
	10	!===========================================================================
	11	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	12	! Date : 25 Mai 2010
	13	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	14	!===========================================================================
	15
	16
	17	#include "YOMCST.h"
	18	#include "YOETHF.h"
	19	#include "FCTTRE.h"
	20	#include "thermcell.h"
[2547]	21	#include "nuage.h"
[1399]	22
	23	INTEGER itap,ind1,ind2
	24	INTEGER ngrid,klev,klon,l,ig
	25
	26	REAL ztv(ngrid,klev)
	27	REAL po(ngrid)
	28	REAL zqenv(ngrid)
	29	REAL zqta(ngrid,klev)
	30
	31	REAL fraca(ngrid,klev+1)
	32	REAL zpspsk(ngrid,klev)
	33	REAL paprs(ngrid,klev+1)
	34	REAL ztla(ngrid,klev)
	35	REAL zthl(ngrid,klev)
	36
	37	REAL zqsatth(ngrid,klev)
	38	REAL zqsatenv(ngrid,klev)
	39
	40
[2267]	41	REAL sigma1(ngrid,klev)
[1399]	42	REAL sigma2(ngrid,klev)
	43	REAL qlth(ngrid,klev)
	44	REAL qlenv(ngrid,klev)
	45	REAL qltot(ngrid,klev)
	46	REAL cth(ngrid,klev)
	47	REAL cenv(ngrid,klev)
	48	REAL ctot(ngrid,klev)
	49	REAL rneb(ngrid,klev)
[2267]	50	REAL t(ngrid,klev)
[1399]	51	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
	52	REAL rdd,cppd,Lv
	53	REAL alth,alenv,ath,aenv
	54	REAL sth,senv,sigma1s,sigma2s,xth,xenv
	55	REAL Tbef,zdelta,qsatbef,zcor
[2586]	56	REAL qlbef
[1399]	57	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
	58
	59	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	60	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	61	REAL zqs(ngrid), qcloud(ngrid)
	62	REAL erf
	63
	64
	65
	66
[2267]	67	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2547]	68	! Gestion de deux versions de cloudth
[2267]	69	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2547]	70
	71	IF (iflag_cloudth_vert.GE.1) THEN
	72	CALL cloudth_vert(ngrid,klev,ind2, &
[2267]	73	& ztv,po,zqta,fraca, &
	74	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	75	& ratqs,zqs,t)
[2547]	76	RETURN
	77	ENDIF
[2267]	78	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	79
	80
	81	!-------------------------------------------------------------------------------
[2140]	82	! Initialisation des variables r?elles
[2267]	83	!-------------------------------------------------------------------------------
[1399]	84	sigma1(:,:)=0.
	85	sigma2(:,:)=0.
	86	qlth(:,:)=0.
	87	qlenv(:,:)=0.
	88	qltot(:,:)=0.
	89	rneb(:,:)=0.
	90	qcloud(:)=0.
	91	cth(:,:)=0.
	92	cenv(:,:)=0.
	93	ctot(:,:)=0.
	94	qsatmmussig1=0.
	95	qsatmmussig2=0.
	96	rdd=287.04
	97	cppd=1005.7
	98	pi=3.14159
	99	Lv=2.5e6
	100	sqrt2pi=sqrt(2.*pi)
	101
	102
	103
[2267]	104	!-------------------------------------------------------------------------------
[2140]	105	! Calcul de la fraction du thermique et des ?cart-types des distributions
[2267]	106	!-------------------------------------------------------------------------------
[1399]	107	do ind1=1,ngrid
	108
	109	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
	110
	111	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	112
	113
	114	! zqenv(ind1)=po(ind1)
	115	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	116	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	117	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	118	qsatbef=MIN(0.5,qsatbef)
	119	zcor=1./(1.-retv*qsatbef)
	120	qsatbef=qsatbef*zcor
	121	zqsatenv(ind1,ind2)=qsatbef
	122
	123
	124
	125
	126	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	127	aenv=1./(1.+(alenv*Lv/cppd))
	128	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	129
	130
	131
	132
	133	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	134	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	135	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	136	qsatbef=MIN(0.5,qsatbef)
	137	zcor=1./(1.-retv*qsatbef)
	138	qsatbef=qsatbef*zcor
	139	zqsatth(ind1,ind2)=qsatbef
	140
	141	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
	142	ath=1./(1.+(alth*Lv/cppd))
	143	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
	144
	145
	146
[2267]	147	!------------------------------------------------------------------------------
[2140]	148	! Calcul des ?cart-types pour s
[2267]	149	!------------------------------------------------------------------------------
[1399]	150
[2140]	151	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	152	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
	153	! if (paprs(ind1,ind2).gt.90000) then
	154	! ratqs(ind1,ind2)=0.002
	155	! else
	156	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	157	! endif
	158	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	159	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	160	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	161	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[1399]	162
[2267]	163	!------------------------------------------------------------------------------
[2140]	164	! Calcul de l'eau condens?e et de la couverture nuageuse
[2267]	165	!------------------------------------------------------------------------------
[1399]	166	sqrt2pi=sqrt(2.*pi)
	167	xth=sth/(sqrt(2.)*sigma2s)
	168	xenv=senv/(sqrt(2.)*sigma1s)
	169	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	170	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	171	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	172
	173	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	174	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	175	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	176
[2267]	177	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1399]	178	if (ctot(ind1,ind2).lt.1.e-10) then
	179	ctot(ind1,ind2)=0.
	180	qcloud(ind1)=zqsatenv(ind1,ind2)
	181
	182	else
	183
	184	ctot(ind1,ind2)=ctot(ind1,ind2)
	185	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	186
	187	endif
	188
	189
	190	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	191
	192
	193	else ! gaussienne environnement seule
	194
	195	zqenv(ind1)=po(ind1)
	196	Tbef=t(ind1,ind2)
	197	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	198	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	199	qsatbef=MIN(0.5,qsatbef)
	200	zcor=1./(1.-retv*qsatbef)
	201	qsatbef=qsatbef*zcor
	202	zqsatenv(ind1,ind2)=qsatbef
	203
	204
	205	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	206	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	207	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	208	aenv=1./(1.+(alenv*Lv/cppd))
	209	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	210
	211
[1411]	212	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[1399]	213
	214	sqrt2pi=sqrt(2.*pi)
	215	xenv=senv/(sqrt(2.)*sigma1s)
	216	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	217	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	218
	219	if (ctot(ind1,ind2).lt.1.e-3) then
	220	ctot(ind1,ind2)=0.
	221	qcloud(ind1)=zqsatenv(ind1,ind2)
	222
	223	else
	224
	225	ctot(ind1,ind2)=ctot(ind1,ind2)
	226	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	227
	228	endif
	229
	230
	231
	232
	233
	234
	235	endif
	236	enddo
	237
	238	return
	239	end
	240
	241
	242
[2267]	243	!===========================================================================
	244	SUBROUTINE cloudth_vert(ngrid,klev,ind2, &
	245	& ztv,po,zqta,fraca, &
	246	& qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	247	& ratqs,zqs,t)
[1399]	248
[2267]	249	!===========================================================================
	250	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	251	! Date : 25 Mai 2010
	252	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	253	!===========================================================================
[1399]	254
	255
[2586]	256	USE ioipsl_getin_p_mod, ONLY : getin_p
	257
	258	IMPLICIT NONE
	259
[2267]	260	#include "YOMCST.h"
	261	#include "YOETHF.h"
	262	#include "FCTTRE.h"
	263	#include "thermcell.h"
[2547]	264	#include "nuage.h"
	265
[2267]	266	INTEGER itap,ind1,ind2
	267	INTEGER ngrid,klev,klon,l,ig
	268
	269	REAL ztv(ngrid,klev)
	270	REAL po(ngrid)
	271	REAL zqenv(ngrid)
	272	REAL zqta(ngrid,klev)
	273
	274	REAL fraca(ngrid,klev+1)
	275	REAL zpspsk(ngrid,klev)
	276	REAL paprs(ngrid,klev+1)
	277	REAL ztla(ngrid,klev)
	278	REAL zthl(ngrid,klev)
	279
	280	REAL zqsatth(ngrid,klev)
	281	REAL zqsatenv(ngrid,klev)
	282
	283
	284	REAL sigma1(ngrid,klev)
	285	REAL sigma2(ngrid,klev)
	286	REAL qlth(ngrid,klev)
	287	REAL qlenv(ngrid,klev)
	288	REAL qltot(ngrid,klev)
	289	REAL cth(ngrid,klev)
	290	REAL cenv(ngrid,klev)
	291	REAL ctot(ngrid,klev)
	292	REAL rneb(ngrid,klev)
	293	REAL t(ngrid,klev)
	294	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
	295	REAL rdd,cppd,Lv,sqrt2,sqrtpi
	296	REAL alth,alenv,ath,aenv
	297	REAL sth,senv,sigma1s,sigma2s,xth,xenv
	298	REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
	299	REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
	300	REAL Tbef,zdelta,qsatbef,zcor
[2586]	301	REAL qlbef
[2267]	302	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
[2586]	303	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	304	! (J Jouhaud, JL Dufresne, JB Madeleine)
	305	REAL,SAVE :: vert_alpha
	306	LOGICAL, SAVE :: firstcall = .TRUE.
[2267]	307
	308	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
	309	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
	310	REAL zqs(ngrid), qcloud(ngrid)
	311	REAL erf
	312
	313	!------------------------------------------------------------------------------
	314	! Initialisation des variables r?elles
	315	!------------------------------------------------------------------------------
	316	sigma1(:,:)=0.
	317	sigma2(:,:)=0.
	318	qlth(:,:)=0.
	319	qlenv(:,:)=0.
	320	qltot(:,:)=0.
	321	rneb(:,:)=0.
	322	qcloud(:)=0.
	323	cth(:,:)=0.
	324	cenv(:,:)=0.
	325	ctot(:,:)=0.
	326	qsatmmussig1=0.
	327	qsatmmussig2=0.
	328	rdd=287.04
	329	cppd=1005.7
	330	pi=3.14159
	331	Lv=2.5e6
	332	sqrt2pi=sqrt(2.*pi)
	333	sqrt2=sqrt(2.)
	334	sqrtpi=sqrt(pi)
	335
[2586]	336	IF (firstcall) THEN
	337	vert_alpha=0.5
	338	CALL getin_p('cloudth_vert_alpha',vert_alpha)
	339	WRITE(,) 'cloudth_vert_alpha = ', vert_alpha
	340	firstcall=.FALSE.
	341	ENDIF
[2267]	342
	343	!-------------------------------------------------------------------------------
	344	! Calcul de la fraction du thermique et des ?cart-types des distributions
	345	!-------------------------------------------------------------------------------
	346	do ind1=1,ngrid
	347
	348	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
	349
	350	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	351
	352
	353	! zqenv(ind1)=po(ind1)
	354	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	355	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	356	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	357	qsatbef=MIN(0.5,qsatbef)
	358	zcor=1./(1.-retv*qsatbef)
	359	qsatbef=qsatbef*zcor
	360	zqsatenv(ind1,ind2)=qsatbef
	361
	362
	363
	364
	365	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	366	aenv=1./(1.+(alenv*Lv/cppd))
	367	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	368
	369
	370
	371
	372	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	373	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	374	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	375	qsatbef=MIN(0.5,qsatbef)
	376	zcor=1./(1.-retv*qsatbef)
	377	qsatbef=qsatbef*zcor
	378	zqsatth(ind1,ind2)=qsatbef
	379
	380	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
	381	ath=1./(1.+(alth*Lv/cppd))
	382	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
	383
	384
	385
	386	!------------------------------------------------------------------------------
	387	! Calcul des ?cart-types pour s
	388	!------------------------------------------------------------------------------
	389
	390	sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	391	sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	392	! if (paprs(ind1,ind2).gt.90000) then
	393	! ratqs(ind1,ind2)=0.002
	394	! else
	395	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	396	! endif
	397	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	398	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	399	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	400	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	401
	402	!------------------------------------------------------------------------------
	403	! Calcul de l'eau condens?e et de la couverture nuageuse
	404	!------------------------------------------------------------------------------
	405	sqrt2pi=sqrt(2.*pi)
	406	xth=sth/(sqrt(2.)*sigma2s)
	407	xenv=senv/(sqrt(2.)*sigma1s)
	408	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	409	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	410	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	411
	412	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	413	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	414	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	415
[2547]	416	IF (iflag_cloudth_vert == 1) THEN
[2267]	417	!-------------------------------------------------------------------------------
	418	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
	419	!-------------------------------------------------------------------------------
	420	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	421	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	422	deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	423	deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	424	! deltasenv=aenv0.01po(ind1)
	425	! deltasth=ath0.01zqta(ind1,ind2)
	426	xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
	427	xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
	428	xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
	429	xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
	430	coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	431	coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	432
	433	cth(ind1,ind2)=0.5(1.+1.erf(xth2))
	434	cenv(ind1,ind2)=0.5(1.+1.erf(xenv2))
	435	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	436
	437	IntJ=sigma1s(exp(-1.xenv1*2)/sqrt2pi)+0.5senv*(1+erf(xenv1))
	438	IntI1=coeffqlenv0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	439	IntI2=coeffqlenvxenv2(exp(-1.xenv22)-exp(-1.xenv1**2))
	440	IntI3=coeffqlenv0.5sqrtpixenv22(erf(xenv2)-erf(xenv1))
	441
	442	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	443	! qlenv(ind1,ind2)=IntJ
	444	! print*, qlenv(ind1,ind2),'VERIF EAU'
	445
	446
	447	IntJ=sigma2s(exp(-1.xth1*2)/sqrt2pi)+0.5sth*(1+erf(xth1))
	448	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
	449	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
	450	IntI1=coeffqlth0.5(0.5sqrtpi(erf(xth2)-erf(xth1))+xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	451	IntI2=coeffqlthxth2(exp(-1.xth22)-exp(-1.xth1**2))
	452	IntI3=coeffqlth0.5sqrtpixth22(erf(xth2)-erf(xth1))
	453	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	454	! qlth(ind1,ind2)=IntJ
	455	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	456	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	457
[2547]	458	ELSE IF (iflag_cloudth_vert == 2) THEN
	459
	460	!-------------------------------------------------------------------------------
	461	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
	462	!-------------------------------------------------------------------------------
	463	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	464	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	465	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	466	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
[2586]	467	deltasenv=aenvvert_alphasigma1s
	468	deltasth=athvert_alphasigma2s
[2547]	469
	470	xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
	471	xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
	472	xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
	473	xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
	474	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	475	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
	476
	477	cth(ind1,ind2)=0.5(1.-1.erf(xth1))
	478	cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
	479	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	480
	481	IntJ=0.5senv(1-erf(xenv2))+(sigma1s/sqrt2pi)exp(-1.xenv2**2)
	482	IntI1=(((senv+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	483	IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp(-1.xenv1*2)-xenv2exp(-1.xenv2*2))
	484	IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))(exp(-1.xenv1*2)-exp(-1.xenv2**2))
	485
	486	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
	487	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
	488	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
	489
	490	qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	491	! qlenv(ind1,ind2)=IntJ
	492	! print*, qlenv(ind1,ind2),'VERIF EAU'
	493
	494	IntJ=0.5sth(1-erf(xth2))+(sigma2s/sqrt2pi)exp(-1.xth2**2)
	495	IntI1=(((sth+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	496	IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp(-1.xth1*2)-xth2exp(-1.xth2*2))
	497	IntI3=((sqrt2sigma2s(sth+deltasth))/(4sqrtpideltasth))(exp(-1.xth1*2)-exp(-1.xth2**2))
	498
	499	qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
	500	! qlth(ind1,ind2)=IntJ
	501	! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	502	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	503
	504
	505
	506
	507	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
	508
[2267]	509	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2547]	510
[2267]	511	if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
	512	ctot(ind1,ind2)=0.
	513	qcloud(ind1)=zqsatenv(ind1,ind2)
	514
	515	else
	516
	517	ctot(ind1,ind2)=ctot(ind1,ind2)
	518	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	519	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	520	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	521
	522	endif
	523
	524
	525
	526	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	527
	528
	529	else ! gaussienne environnement seule
	530
	531	zqenv(ind1)=po(ind1)
	532	Tbef=t(ind1,ind2)
	533	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	534	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	535	qsatbef=MIN(0.5,qsatbef)
	536	zcor=1./(1.-retv*qsatbef)
	537	qsatbef=qsatbef*zcor
	538	zqsatenv(ind1,ind2)=qsatbef
	539
	540
	541	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	542	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	543	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	544	aenv=1./(1.+(alenv*Lv/cppd))
	545	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	546
	547
	548	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	549
	550	sqrt2pi=sqrt(2.*pi)
	551	xenv=senv/(sqrt(2.)*sigma1s)
	552	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	553	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	554
	555	if (ctot(ind1,ind2).lt.1.e-3) then
	556	ctot(ind1,ind2)=0.
	557	qcloud(ind1)=zqsatenv(ind1,ind2)
	558
	559	else
	560
	561	ctot(ind1,ind2)=ctot(ind1,ind2)
	562	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	563
	564	endif
	565
	566
	567
	568
	569
	570
	571	endif
	572	enddo
	573
	574	return
	575	end

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

Download in other formats:

Original Format