Context Navigation

cloudth.F90 @ 2117

Last change on this file since 2117 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format