Context Navigation

cloudth.F90 @ 2273

Last change on this file since 2273 was 2267, checked in by fhourdin, 10 years ago

Modification du modèle du thermique pour diminuer entre autres la
sensibilité au pas de temps. Convergence numérique si iflag_thermals_ed=7.
Ne concerne que la version iflag_thermals_ed=8

Premier essaie de prise en compte d'une variabilité verticale sous maille
de l'eau pour accroître la converture nuageuse.
Sous la clé iflag_coudth_vert=1 (1 par défaut)

Modification of the thermal plume model in case iflag_thermals_ed=8
Modification of the bigaussian scheme for clouds associated with

thermals. Activated if iflag_cloudth_vert=1

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

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

Download in other formats:

Original Format