Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cloudth.F90 @ 5360

Last change on this file since 5360 was 1399, checked in by idelkadi, 15 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

physiq.F : appel modifie a calltherm et fisrtilp
thermcell_main.F90 : changement dans les flag_thermals_ed
thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

cloudth.F90 : PDF d'eau bi-gaussiennes

File size: 8.0 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 "iniprint.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
	41	REAL sigma1(ngrid,klev)
	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)
	50	REAL t(ngrid,klev)
	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
	64
	65
	66
	67
	68	! print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), &
	69	! & fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), &
	70	! & 'verif'
	71
	72
	73	! LOGICAL active(ngrid)
	74
	75	!-----------------------------------------------------------------------------------------------------------------
	76	! Initialisation des variables réelles
	77	!-----------------------------------------------------------------------------------------------------------------
	78	sigma1(:,:)=0.
	79	sigma2(:,:)=0.
	80	qlth(:,:)=0.
	81	qlenv(:,:)=0.
	82	qltot(:,:)=0.
	83	rneb(:,:)=0.
	84	qcloud(:)=0.
	85	cth(:,:)=0.
	86	cenv(:,:)=0.
	87	ctot(:,:)=0.
	88	qsatmmussig1=0.
	89	qsatmmussig2=0.
	90	rdd=287.04
	91	cppd=1005.7
	92	pi=3.14159
	93	Lv=2.5e6
	94	sqrt2pi=sqrt(2.*pi)
	95
	96
	97
	98	!------------------------------------------------------------------------------------------------------------------
	99	! Calcul de la fraction du thermique et des écart-types des distributions
	100	!------------------------------------------------------------------------------------------------------------------
	101	do ind1=1,ngrid
	102
	103	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
	104
	105	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
	106
	107
	108	! zqenv(ind1)=po(ind1)
	109	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
	110	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	111	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	112	qsatbef=MIN(0.5,qsatbef)
	113	zcor=1./(1.-retv*qsatbef)
	114	qsatbef=qsatbef*zcor
	115	zqsatenv(ind1,ind2)=qsatbef
	116
	117
	118
	119
	120	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	121	aenv=1./(1.+(alenv*Lv/cppd))
	122	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	123
	124
	125
	126
	127	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
	128	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	129	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	130	qsatbef=MIN(0.5,qsatbef)
	131	zcor=1./(1.-retv*qsatbef)
	132	qsatbef=qsatbef*zcor
	133	zqsatth(ind1,ind2)=qsatbef
	134
	135	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
	136	ath=1./(1.+(alth*Lv/cppd))
	137	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
	138
	139
	140
	141	!-----------------------------------------------------------------------------------------------------------------
	142	! Calcul des écart-types pour s
	143	!-----------------------------------------------------------------------------------------------------------------
	144	! alpha=0.5*(fraca(ind1,ind2)+fraca(ind1,ind2-1))
	145	! sigma1s=(Max(1.2fraca(ind1,ind2)(sth-senv)2,(senv/100)2))**0.5
	146	! sigma2s=((0.021/(fraca(ind1,ind2)+0.0)*0.5)(sth-senv)2+(sth/100)2)**0.5
	147	! sigma1s=(1.5*0.5)(fraca(ind1,ind2)*0.65)(sth-senv)+0.000045
	148	! sigma2s=0.1265(sth-senv)/(fraca(ind1,ind2)+0.0)*0.35+0.000045
	149
	150	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.00003
	151	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
	152
	153	! sigma1s=(1.5*0.5)(alpha*0.65)(sth-senv)+0.000045
	154	! sigma2s=0.1265(sth-senv)/alpha*0.35+0.000045
	155
	156	! sigma1s=(2.8*0.5)(0.1*0.7)(sth-senv)+0.00002
	157	! sigma2s=((0.126/(0.1)*0.3)(sth-senv))+0.00002
	158
	159
	160
	161	!-----------------------------------------------------------------------------------------------------------------
	162	! Calcul de l'eau condensée et de la couverture nuageuse
	163	!-----------------------------------------------------------------------------------------------------------------
	164	sqrt2pi=sqrt(2.*pi)
	165	xth=sth/(sqrt(2.)*sigma2s)
	166	xenv=senv/(sqrt(2.)*sigma1s)
	167	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	168	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	169	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	170	! ctot(ind1,ind2)=alphacth(ind1,ind2)+(1.-1.alpha)*cenv(ind1,ind2)
	171
	172
	173
	174	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	175	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	176	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	177	! qltot(ind1,ind2)=alphaqlth(ind1,ind2)+(1.-1.alpha)*qlenv(ind1,ind2)
	178
	179
	180	! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
	181
	182	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	183	if (ctot(ind1,ind2).lt.1.e-10) then
	184	ctot(ind1,ind2)=0.
	185	qcloud(ind1)=zqsatenv(ind1,ind2)
	186
	187	else
	188
	189	ctot(ind1,ind2)=ctot(ind1,ind2)
	190	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	191
	192	endif
	193
	194
	195	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	196
	197
	198	else ! gaussienne environnement seule
	199
	200	zqenv(ind1)=po(ind1)
	201	Tbef=t(ind1,ind2)
	202	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	203	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	204	qsatbef=MIN(0.5,qsatbef)
	205	zcor=1./(1.-retv*qsatbef)
	206	qsatbef=qsatbef*zcor
	207	zqsatenv(ind1,ind2)=qsatbef
	208
	209
	210	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	211	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	212	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	213	aenv=1./(1.+(alenv*Lv/cppd))
	214	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	215
	216	! if (zqenv(ind1).gt.0.005) then
	217	sigma1s=0.005*zqenv(ind1)
	218	! else
	219	! sigma1s=0.0005
	220	! endif
	221
	222	! sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
	223
	224
	225
	226	! sigma1s=0.00003
	227	sqrt2pi=sqrt(2.*pi)
	228	xenv=senv/(sqrt(2.)*sigma1s)
	229	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	230	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	231
	232	if (ctot(ind1,ind2).lt.1.e-3) then
	233	ctot(ind1,ind2)=0.
	234	qcloud(ind1)=zqsatenv(ind1,ind2)
	235
	236	else
	237
	238	ctot(ind1,ind2)=ctot(ind1,ind2)
	239	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	240
	241	endif
	242
	243
	244
	245
	246
	247
	248	endif
	249	enddo
	250
	251	return
	252	end
	253
	254
	255
	256
	257
	258
	259
	260
	261
	262
	263

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

Download in other formats:

Original Format