Context Navigation

cloudth.F90 @ 2404

Last change on this file since 2404 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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

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

Download in other formats:

Original Format