Context Navigation

cloudth.F90 @ 2212

Last change on this file since 2212 was 2140, checked in by fhourdin, 10 years ago

Nouvelle version des thermiques permettant de retrouver a 79 niveaux

des comportements proches de ceux obtenus dans la thèse d'Arnaud
Jam sur les cas de strato-cumulus.
thermcell_plume et cloudth modifiés.

New versions of thermal plume model for stratocumulus cases.

Arnaud Jam

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: 7.8 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	!-----------------------------------------------------------------------------------------------------------------
[2140]	76	! Initialisation des variables r?elles
[1399]	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	!------------------------------------------------------------------------------------------------------------------
[2140]	99	! Calcul de la fraction du thermique et des ?cart-types des distributions
[1399]	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	!-----------------------------------------------------------------------------------------------------------------
[2140]	142	! Calcul des ?cart-types pour s
[1399]	143	!-----------------------------------------------------------------------------------------------------------------
	144
[2140]	145	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	146	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
	147	! if (paprs(ind1,ind2).gt.90000) then
	148	! ratqs(ind1,ind2)=0.002
	149	! else
	150	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	151	! endif
	152	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	153	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	154	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	155	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[1399]	156
	157	!-----------------------------------------------------------------------------------------------------------------
[2140]	158	! Calcul de l'eau condens?e et de la couverture nuageuse
[1399]	159	!-----------------------------------------------------------------------------------------------------------------
	160	sqrt2pi=sqrt(2.*pi)
	161	xth=sth/(sqrt(2.)*sigma2s)
	162	xenv=senv/(sqrt(2.)*sigma1s)
	163	cth(ind1,ind2)=0.5(1.+1.erf(xth))
	164	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
	165	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	166	! ctot(ind1,ind2)=alphacth(ind1,ind2)+(1.-1.alpha)*cenv(ind1,ind2)
	167
	168
	169
	170	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
	171	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	172	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
	173	! qltot(ind1,ind2)=alphaqlth(ind1,ind2)+(1.-1.alpha)*qlenv(ind1,ind2)
	174
	175
	176	! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
	177
	178	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	179	if (ctot(ind1,ind2).lt.1.e-10) then
	180	ctot(ind1,ind2)=0.
	181	qcloud(ind1)=zqsatenv(ind1,ind2)
	182
	183	else
	184
	185	ctot(ind1,ind2)=ctot(ind1,ind2)
	186	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
	187
	188	endif
	189
	190
	191	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	192
	193
	194	else ! gaussienne environnement seule
	195
	196	zqenv(ind1)=po(ind1)
	197	Tbef=t(ind1,ind2)
	198	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	199	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
	200	qsatbef=MIN(0.5,qsatbef)
	201	zcor=1./(1.-retv*qsatbef)
	202	qsatbef=qsatbef*zcor
	203	zqsatenv(ind1,ind2)=qsatbef
	204
	205
	206	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	207	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
	208	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
	209	aenv=1./(1.+(alenv*Lv/cppd))
	210	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
	211
	212
[1411]	213	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
[1399]	214
	215	sqrt2pi=sqrt(2.*pi)
	216	xenv=senv/(sqrt(2.)*sigma1s)
	217	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
	218	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
	219
	220	if (ctot(ind1,ind2).lt.1.e-3) then
	221	ctot(ind1,ind2)=0.
	222	qcloud(ind1)=zqsatenv(ind1,ind2)
	223
	224	else
	225
	226	ctot(ind1,ind2)=ctot(ind1,ind2)
	227	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
	228
	229	endif
	230
	231
	232
	233
	234
	235
	236	endif
	237	enddo
	238
	239	return
	240	end
	241
	242
	243
	244
	245
	246
	247
	248
	249
	250

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

Download in other formats:

Original Format