Context Navigation

cloudth.F90 @ 3882

Last change on this file since 3882 was 3817, checked in by millour, 10 years ago
Further cleanup and removal of references to iniprint.h. Also added bench testcase 48x36x19. EM
File size: 7.8 KB

Line
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 "thermcell.h"
21
22	INTEGER itap,ind1,ind2
23	INTEGER ngrid,klev,klon,l,ig
24
25	REAL ztv(ngrid,klev)
26	REAL po(ngrid)
27	REAL zqenv(ngrid)
28	REAL zqta(ngrid,klev)
29
30	REAL fraca(ngrid,klev+1)
31	REAL zpspsk(ngrid,klev)
32	REAL paprs(ngrid,klev+1)
33	REAL ztla(ngrid,klev)
34	REAL zthl(ngrid,klev)
35
36	REAL zqsatth(ngrid,klev)
37	REAL zqsatenv(ngrid,klev)
38
39
40	REAL sigma1(ngrid,klev)
41	REAL sigma2(ngrid,klev)
42	REAL qlth(ngrid,klev)
43	REAL qlenv(ngrid,klev)
44	REAL qltot(ngrid,klev)
45	REAL cth(ngrid,klev)
46	REAL cenv(ngrid,klev)
47	REAL ctot(ngrid,klev)
48	REAL rneb(ngrid,klev)
49	REAL t(ngrid,klev)
50	REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
51	REAL rdd,cppd,Lv
52	REAL alth,alenv,ath,aenv
53	REAL sth,senv,sigma1s,sigma2s,xth,xenv
54	REAL Tbef,zdelta,qsatbef,zcor
55	REAL alpha,qlbef
56	REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
57
58	REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
59	REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
60	REAL zqs(ngrid), qcloud(ngrid)
61	REAL erf
62
63
64
65
66
67	! print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), &
68	! & fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), &
69	! & 'verif'
70
71
72	! LOGICAL active(ngrid)
73
74	!-----------------------------------------------------------------------------------------------------------------
75	! Initialisation des variables r?elles
76	!-----------------------------------------------------------------------------------------------------------------
77	sigma1(:,:)=0.
78	sigma2(:,:)=0.
79	qlth(:,:)=0.
80	qlenv(:,:)=0.
81	qltot(:,:)=0.
82	rneb(:,:)=0.
83	qcloud(:)=0.
84	cth(:,:)=0.
85	cenv(:,:)=0.
86	ctot(:,:)=0.
87	qsatmmussig1=0.
88	qsatmmussig2=0.
89	rdd=287.04
90	cppd=1005.7
91	pi=3.14159
92	Lv=2.5e6
93	sqrt2pi=sqrt(2.*pi)
94
95
96
97	!------------------------------------------------------------------------------------------------------------------
98	! Calcul de la fraction du thermique et des ?cart-types des distributions
99	!------------------------------------------------------------------------------------------------------------------
100	do ind1=1,ngrid
101
102	if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
103
104	zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
105
106
107	! zqenv(ind1)=po(ind1)
108	Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
109	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
110	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
111	qsatbef=MIN(0.5,qsatbef)
112	zcor=1./(1.-retv*qsatbef)
113	qsatbef=qsatbef*zcor
114	zqsatenv(ind1,ind2)=qsatbef
115
116
117
118
119	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
120	aenv=1./(1.+(alenv*Lv/cppd))
121	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
122
123
124
125
126	Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
127	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
128	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
129	qsatbef=MIN(0.5,qsatbef)
130	zcor=1./(1.-retv*qsatbef)
131	qsatbef=qsatbef*zcor
132	zqsatth(ind1,ind2)=qsatbef
133
134	alth=(0.622Lvzqsatth(ind1,ind2))/(rddztla(ind1,ind2)*2)
135	ath=1./(1.+(alth*Lv/cppd))
136	sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
137
138
139
140	!-----------------------------------------------------------------------------------------------------------------
141	! Calcul des ?cart-types pour s
142	!-----------------------------------------------------------------------------------------------------------------
143
144	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
145	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
146	! if (paprs(ind1,ind2).gt.90000) then
147	! ratqs(ind1,ind2)=0.002
148	! else
149	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
150	! endif
151	sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
152	sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
153	! sigma1s=ratqs(ind1,ind2)*po(ind1)
154	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
155
156	!-----------------------------------------------------------------------------------------------------------------
157	! Calcul de l'eau condens?e et de la couverture nuageuse
158	!-----------------------------------------------------------------------------------------------------------------
159	sqrt2pi=sqrt(2.*pi)
160	xth=sth/(sqrt(2.)*sigma2s)
161	xenv=senv/(sqrt(2.)*sigma1s)
162	cth(ind1,ind2)=0.5(1.+1.erf(xth))
163	cenv(ind1,ind2)=0.5(1.+1.erf(xenv))
164	ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
165	! ctot(ind1,ind2)=alphacth(ind1,ind2)+(1.-1.alpha)*cenv(ind1,ind2)
166
167
168
169	qlth(ind1,ind2)=sigma2s((exp(-1.xth*2)/sqrt2pi)+xthsqrt(2.)*cth(ind1,ind2))
170	qlenv(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
171	qltot(ind1,ind2)=fraca(ind1,ind2)qlth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qlenv(ind1,ind2)
172	! qltot(ind1,ind2)=alphaqlth(ind1,ind2)+(1.-1.alpha)*qlenv(ind1,ind2)
173
174
175	! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
176
177	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178	if (ctot(ind1,ind2).lt.1.e-10) then
179	ctot(ind1,ind2)=0.
180	qcloud(ind1)=zqsatenv(ind1,ind2)
181
182	else
183
184	ctot(ind1,ind2)=ctot(ind1,ind2)
185	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
186
187	endif
188
189
190	! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
191
192
193	else ! gaussienne environnement seule
194
195	zqenv(ind1)=po(ind1)
196	Tbef=t(ind1,ind2)
197	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
198	qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
199	qsatbef=MIN(0.5,qsatbef)
200	zcor=1./(1.-retv*qsatbef)
201	qsatbef=qsatbef*zcor
202	zqsatenv(ind1,ind2)=qsatbef
203
204
205	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
206	zthl(ind1,ind2)=t(ind1,ind2)(101325/paprs(ind1,ind2))*(rdd/cppd)
207	alenv=(0.622Lvzqsatenv(ind1,ind2))/(rddzthl(ind1,ind2)*2)
208	aenv=1./(1.+(alenv*Lv/cppd))
209	senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
210
211
212	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
213
214	sqrt2pi=sqrt(2.*pi)
215	xenv=senv/(sqrt(2.)*sigma1s)
216	ctot(ind1,ind2)=0.5(1.+1.erf(xenv))
217	qltot(ind1,ind2)=sigma1s((exp(-1.xenv*2)/sqrt2pi)+xenvsqrt(2.)*cenv(ind1,ind2))
218
219	if (ctot(ind1,ind2).lt.1.e-3) then
220	ctot(ind1,ind2)=0.
221	qcloud(ind1)=zqsatenv(ind1,ind2)
222
223	else
224
225	ctot(ind1,ind2)=ctot(ind1,ind2)
226	qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
227
228	endif
229
230
231
232
233
234
235	endif
236	enddo
237
238	return
239	end
240
241
242
243
244
245
246
247
248
249

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

Download in other formats:

Original Format