Context Navigation

cloudth.F90 @ 2558

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

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

Download in other formats:

Original Format