Context Navigation

cloudth.F90 @ 3069

Last change on this file since 3069 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
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

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 "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
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
156
157	!-----------------------------------------------------------------------------------------------------------------
158	! Calcul de l'eau condens?e et de la couverture nuageuse
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
213	sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
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