1 |
|
|
MODULE cloudth_mod |
2 |
|
|
|
3 |
|
|
|
4 |
|
|
IMPLICIT NONE |
5 |
|
|
|
6 |
|
|
CONTAINS |
7 |
|
|
|
8 |
|
|
SUBROUTINE cloudth(ngrid,klev,ind2, & |
9 |
|
|
& ztv,po,zqta,fraca, & |
10 |
|
|
& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
11 |
|
|
& ratqs,zqs,t) |
12 |
|
|
|
13 |
|
|
|
14 |
|
|
use lscp_ini_mod, only: iflag_cloudth_vert |
15 |
|
|
|
16 |
|
|
IMPLICIT NONE |
17 |
|
|
|
18 |
|
|
|
19 |
|
|
!=========================================================================== |
20 |
|
|
! Auteur : Arnaud Octavio Jam (LMD/CNRS) |
21 |
|
|
! Date : 25 Mai 2010 |
22 |
|
|
! Objet : calcule les valeurs de qc et rneb dans les thermiques |
23 |
|
|
!=========================================================================== |
24 |
|
|
|
25 |
|
|
|
26 |
|
|
INCLUDE "YOMCST.h" |
27 |
|
|
INCLUDE "YOETHF.h" |
28 |
|
|
INCLUDE "FCTTRE.h" |
29 |
|
|
|
30 |
|
|
INTEGER itap,ind1,ind2 |
31 |
|
|
INTEGER ngrid,klev,klon,l,ig |
32 |
|
|
|
33 |
|
|
REAL ztv(ngrid,klev) |
34 |
|
|
REAL po(ngrid) |
35 |
|
|
REAL zqenv(ngrid) |
36 |
|
|
REAL zqta(ngrid,klev) |
37 |
|
|
|
38 |
|
|
REAL fraca(ngrid,klev+1) |
39 |
|
|
REAL zpspsk(ngrid,klev) |
40 |
|
|
REAL paprs(ngrid,klev+1) |
41 |
|
|
REAL pplay(ngrid,klev) |
42 |
|
|
REAL ztla(ngrid,klev) |
43 |
|
|
REAL zthl(ngrid,klev) |
44 |
|
|
|
45 |
|
|
REAL zqsatth(ngrid,klev) |
46 |
|
|
REAL zqsatenv(ngrid,klev) |
47 |
|
|
|
48 |
|
|
|
49 |
|
|
REAL sigma1(ngrid,klev) |
50 |
|
|
REAL sigma2(ngrid,klev) |
51 |
|
|
REAL qlth(ngrid,klev) |
52 |
|
|
REAL qlenv(ngrid,klev) |
53 |
|
|
REAL qltot(ngrid,klev) |
54 |
|
|
REAL cth(ngrid,klev) |
55 |
|
|
REAL cenv(ngrid,klev) |
56 |
|
|
REAL ctot(ngrid,klev) |
57 |
|
|
REAL rneb(ngrid,klev) |
58 |
|
|
REAL t(ngrid,klev) |
59 |
|
|
REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi |
60 |
|
|
REAL rdd,cppd,Lv |
61 |
|
|
REAL alth,alenv,ath,aenv |
62 |
|
|
REAL sth,senv,sigma1s,sigma2s,xth,xenv |
63 |
|
|
REAL Tbef,zdelta,qsatbef,zcor |
64 |
|
|
REAL qlbef |
65 |
|
|
REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur |
66 |
|
|
|
67 |
|
|
REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) |
68 |
|
|
REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) |
69 |
|
|
REAL zqs(ngrid), qcloud(ngrid) |
70 |
|
|
REAL erf |
71 |
|
|
|
72 |
|
|
|
73 |
|
|
|
74 |
|
|
|
75 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
76 |
|
|
! Gestion de deux versions de cloudth |
77 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
78 |
|
|
|
79 |
|
|
IF (iflag_cloudth_vert.GE.1) THEN |
80 |
|
|
CALL cloudth_vert(ngrid,klev,ind2, & |
81 |
|
|
& ztv,po,zqta,fraca, & |
82 |
|
|
& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
83 |
|
|
& ratqs,zqs,t) |
84 |
|
|
RETURN |
85 |
|
|
ENDIF |
86 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
87 |
|
|
|
88 |
|
|
|
89 |
|
|
!------------------------------------------------------------------------------- |
90 |
|
|
! Initialisation des variables r?elles |
91 |
|
|
!------------------------------------------------------------------------------- |
92 |
|
|
sigma1(:,:)=0. |
93 |
|
|
sigma2(:,:)=0. |
94 |
|
|
qlth(:,:)=0. |
95 |
|
|
qlenv(:,:)=0. |
96 |
|
|
qltot(:,:)=0. |
97 |
|
|
rneb(:,:)=0. |
98 |
|
|
qcloud(:)=0. |
99 |
|
|
cth(:,:)=0. |
100 |
|
|
cenv(:,:)=0. |
101 |
|
|
ctot(:,:)=0. |
102 |
|
|
qsatmmussig1=0. |
103 |
|
|
qsatmmussig2=0. |
104 |
|
|
rdd=287.04 |
105 |
|
|
cppd=1005.7 |
106 |
|
|
pi=3.14159 |
107 |
|
|
Lv=2.5e6 |
108 |
|
|
sqrt2pi=sqrt(2.*pi) |
109 |
|
|
|
110 |
|
|
|
111 |
|
|
|
112 |
|
|
!------------------------------------------------------------------------------- |
113 |
|
|
! Calcul de la fraction du thermique et des ?cart-types des distributions |
114 |
|
|
!------------------------------------------------------------------------------- |
115 |
|
|
do ind1=1,ngrid |
116 |
|
|
|
117 |
|
|
if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then |
118 |
|
|
|
119 |
|
|
zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) |
120 |
|
|
|
121 |
|
|
|
122 |
|
|
! zqenv(ind1)=po(ind1) |
123 |
|
|
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
124 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
125 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
126 |
|
|
qsatbef=MIN(0.5,qsatbef) |
127 |
|
|
zcor=1./(1.-retv*qsatbef) |
128 |
|
|
qsatbef=qsatbef*zcor |
129 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
130 |
|
|
|
131 |
|
|
|
132 |
|
|
|
133 |
|
|
|
134 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
135 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
136 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
137 |
|
|
|
138 |
|
|
|
139 |
|
|
|
140 |
|
|
|
141 |
|
|
Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
142 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
143 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
144 |
|
|
qsatbef=MIN(0.5,qsatbef) |
145 |
|
|
zcor=1./(1.-retv*qsatbef) |
146 |
|
|
qsatbef=qsatbef*zcor |
147 |
|
|
zqsatth(ind1,ind2)=qsatbef |
148 |
|
|
|
149 |
|
|
alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) |
150 |
|
|
ath=1./(1.+(alth*Lv/cppd)) |
151 |
|
|
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) |
152 |
|
|
|
153 |
|
|
|
154 |
|
|
|
155 |
|
|
!------------------------------------------------------------------------------ |
156 |
|
|
! Calcul des ?cart-types pour s |
157 |
|
|
!------------------------------------------------------------------------------ |
158 |
|
|
|
159 |
|
|
! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) |
160 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) |
161 |
|
|
! if (paprs(ind1,ind2).gt.90000) then |
162 |
|
|
! ratqs(ind1,ind2)=0.002 |
163 |
|
|
! else |
164 |
|
|
! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 |
165 |
|
|
! endif |
166 |
|
|
sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) |
167 |
|
|
sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) |
168 |
|
|
! sigma1s=ratqs(ind1,ind2)*po(ind1) |
169 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 |
170 |
|
|
|
171 |
|
|
!------------------------------------------------------------------------------ |
172 |
|
|
! Calcul de l'eau condens?e et de la couverture nuageuse |
173 |
|
|
!------------------------------------------------------------------------------ |
174 |
|
|
sqrt2pi=sqrt(2.*pi) |
175 |
|
|
xth=sth/(sqrt(2.)*sigma2s) |
176 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
177 |
|
|
cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) |
178 |
|
|
cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
179 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
180 |
|
|
|
181 |
|
|
qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) |
182 |
|
|
qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
183 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
184 |
|
|
|
185 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
186 |
|
|
if (ctot(ind1,ind2).lt.1.e-10) then |
187 |
|
|
ctot(ind1,ind2)=0. |
188 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
189 |
|
|
|
190 |
|
|
else |
191 |
|
|
|
192 |
|
|
ctot(ind1,ind2)=ctot(ind1,ind2) |
193 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
194 |
|
|
|
195 |
|
|
endif |
196 |
|
|
|
197 |
|
|
|
198 |
|
|
! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' |
199 |
|
|
|
200 |
|
|
|
201 |
|
|
else ! gaussienne environnement seule |
202 |
|
|
|
203 |
|
|
zqenv(ind1)=po(ind1) |
204 |
|
|
Tbef=t(ind1,ind2) |
205 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
206 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
207 |
|
|
qsatbef=MIN(0.5,qsatbef) |
208 |
|
|
zcor=1./(1.-retv*qsatbef) |
209 |
|
|
qsatbef=qsatbef*zcor |
210 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
211 |
|
|
|
212 |
|
|
|
213 |
|
|
! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) |
214 |
|
|
zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) |
215 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
216 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
217 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
218 |
|
|
|
219 |
|
|
|
220 |
|
|
sigma1s=ratqs(ind1,ind2)*zqenv(ind1) |
221 |
|
|
|
222 |
|
|
sqrt2pi=sqrt(2.*pi) |
223 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
224 |
|
|
ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
225 |
|
|
qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
226 |
|
|
|
227 |
|
|
if (ctot(ind1,ind2).lt.1.e-3) then |
228 |
|
|
ctot(ind1,ind2)=0. |
229 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
230 |
|
|
|
231 |
|
|
else |
232 |
|
|
|
233 |
|
|
ctot(ind1,ind2)=ctot(ind1,ind2) |
234 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) |
235 |
|
|
|
236 |
|
|
endif |
237 |
|
|
|
238 |
|
|
|
239 |
|
|
|
240 |
|
|
|
241 |
|
|
|
242 |
|
|
|
243 |
|
|
endif |
244 |
|
|
enddo |
245 |
|
|
|
246 |
|
|
return |
247 |
|
|
! end |
248 |
|
|
END SUBROUTINE cloudth |
249 |
|
|
|
250 |
|
|
|
251 |
|
|
|
252 |
|
|
!=========================================================================== |
253 |
|
|
SUBROUTINE cloudth_vert(ngrid,klev,ind2, & |
254 |
|
|
& ztv,po,zqta,fraca, & |
255 |
|
|
& qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
256 |
|
|
& ratqs,zqs,t) |
257 |
|
|
|
258 |
|
|
!=========================================================================== |
259 |
|
|
! Auteur : Arnaud Octavio Jam (LMD/CNRS) |
260 |
|
|
! Date : 25 Mai 2010 |
261 |
|
|
! Objet : calcule les valeurs de qc et rneb dans les thermiques |
262 |
|
|
!=========================================================================== |
263 |
|
|
|
264 |
|
|
|
265 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
266 |
|
|
use lscp_ini_mod, only: iflag_cloudth_vert |
267 |
|
|
|
268 |
|
|
IMPLICIT NONE |
269 |
|
|
|
270 |
|
|
INCLUDE "YOMCST.h" |
271 |
|
|
INCLUDE "YOETHF.h" |
272 |
|
|
INCLUDE "FCTTRE.h" |
273 |
|
|
|
274 |
|
|
INTEGER itap,ind1,ind2 |
275 |
|
|
INTEGER ngrid,klev,klon,l,ig |
276 |
|
|
|
277 |
|
|
REAL ztv(ngrid,klev) |
278 |
|
|
REAL po(ngrid) |
279 |
|
|
REAL zqenv(ngrid) |
280 |
|
|
REAL zqta(ngrid,klev) |
281 |
|
|
|
282 |
|
|
REAL fraca(ngrid,klev+1) |
283 |
|
|
REAL zpspsk(ngrid,klev) |
284 |
|
|
REAL paprs(ngrid,klev+1) |
285 |
|
|
REAL pplay(ngrid,klev) |
286 |
|
|
REAL ztla(ngrid,klev) |
287 |
|
|
REAL zthl(ngrid,klev) |
288 |
|
|
|
289 |
|
|
REAL zqsatth(ngrid,klev) |
290 |
|
|
REAL zqsatenv(ngrid,klev) |
291 |
|
|
|
292 |
|
|
|
293 |
|
|
REAL sigma1(ngrid,klev) |
294 |
|
|
REAL sigma2(ngrid,klev) |
295 |
|
|
REAL qlth(ngrid,klev) |
296 |
|
|
REAL qlenv(ngrid,klev) |
297 |
|
|
REAL qltot(ngrid,klev) |
298 |
|
|
REAL cth(ngrid,klev) |
299 |
|
|
REAL cenv(ngrid,klev) |
300 |
|
|
REAL ctot(ngrid,klev) |
301 |
|
|
REAL rneb(ngrid,klev) |
302 |
|
|
REAL t(ngrid,klev) |
303 |
|
|
REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi |
304 |
|
|
REAL rdd,cppd,Lv,sqrt2,sqrtpi |
305 |
|
|
REAL alth,alenv,ath,aenv |
306 |
|
|
REAL sth,senv,sigma1s,sigma2s,xth,xenv |
307 |
|
|
REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv |
308 |
|
|
REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth |
309 |
|
|
REAL Tbef,zdelta,qsatbef,zcor |
310 |
|
|
REAL qlbef |
311 |
|
|
REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur |
312 |
|
|
! Change the width of the PDF used for vertical subgrid scale heterogeneity |
313 |
|
|
! (J Jouhaud, JL Dufresne, JB Madeleine) |
314 |
|
|
REAL,SAVE :: vert_alpha |
315 |
|
|
!$OMP THREADPRIVATE(vert_alpha) |
316 |
|
|
LOGICAL, SAVE :: firstcall = .TRUE. |
317 |
|
|
!$OMP THREADPRIVATE(firstcall) |
318 |
|
|
|
319 |
|
|
REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) |
320 |
|
|
REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) |
321 |
|
|
REAL zqs(ngrid), qcloud(ngrid) |
322 |
|
|
REAL erf |
323 |
|
|
|
324 |
|
|
!------------------------------------------------------------------------------ |
325 |
|
|
! Initialisation des variables r?elles |
326 |
|
|
!------------------------------------------------------------------------------ |
327 |
|
|
sigma1(:,:)=0. |
328 |
|
|
sigma2(:,:)=0. |
329 |
|
|
qlth(:,:)=0. |
330 |
|
|
qlenv(:,:)=0. |
331 |
|
|
qltot(:,:)=0. |
332 |
|
|
rneb(:,:)=0. |
333 |
|
|
qcloud(:)=0. |
334 |
|
|
cth(:,:)=0. |
335 |
|
|
cenv(:,:)=0. |
336 |
|
|
ctot(:,:)=0. |
337 |
|
|
qsatmmussig1=0. |
338 |
|
|
qsatmmussig2=0. |
339 |
|
|
rdd=287.04 |
340 |
|
|
cppd=1005.7 |
341 |
|
|
pi=3.14159 |
342 |
|
|
Lv=2.5e6 |
343 |
|
|
sqrt2pi=sqrt(2.*pi) |
344 |
|
|
sqrt2=sqrt(2.) |
345 |
|
|
sqrtpi=sqrt(pi) |
346 |
|
|
|
347 |
|
|
IF (firstcall) THEN |
348 |
|
|
vert_alpha=0.5 |
349 |
|
|
CALL getin_p('cloudth_vert_alpha',vert_alpha) |
350 |
|
|
WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha |
351 |
|
|
firstcall=.FALSE. |
352 |
|
|
ENDIF |
353 |
|
|
|
354 |
|
|
!------------------------------------------------------------------------------- |
355 |
|
|
! Calcul de la fraction du thermique et des ?cart-types des distributions |
356 |
|
|
!------------------------------------------------------------------------------- |
357 |
|
|
do ind1=1,ngrid |
358 |
|
|
|
359 |
|
|
if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then |
360 |
|
|
|
361 |
|
|
zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) |
362 |
|
|
|
363 |
|
|
|
364 |
|
|
! zqenv(ind1)=po(ind1) |
365 |
|
|
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
366 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
367 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
368 |
|
|
qsatbef=MIN(0.5,qsatbef) |
369 |
|
|
zcor=1./(1.-retv*qsatbef) |
370 |
|
|
qsatbef=qsatbef*zcor |
371 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
372 |
|
|
|
373 |
|
|
|
374 |
|
|
|
375 |
|
|
|
376 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
377 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
378 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
379 |
|
|
|
380 |
|
|
|
381 |
|
|
|
382 |
|
|
|
383 |
|
|
Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
384 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
385 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
386 |
|
|
qsatbef=MIN(0.5,qsatbef) |
387 |
|
|
zcor=1./(1.-retv*qsatbef) |
388 |
|
|
qsatbef=qsatbef*zcor |
389 |
|
|
zqsatth(ind1,ind2)=qsatbef |
390 |
|
|
|
391 |
|
|
alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) |
392 |
|
|
ath=1./(1.+(alth*Lv/cppd)) |
393 |
|
|
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) |
394 |
|
|
|
395 |
|
|
|
396 |
|
|
|
397 |
|
|
!------------------------------------------------------------------------------ |
398 |
|
|
! Calcul des ?cart-types pour s |
399 |
|
|
!------------------------------------------------------------------------------ |
400 |
|
|
|
401 |
|
|
sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) |
402 |
|
|
sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) |
403 |
|
|
! if (paprs(ind1,ind2).gt.90000) then |
404 |
|
|
! ratqs(ind1,ind2)=0.002 |
405 |
|
|
! else |
406 |
|
|
! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 |
407 |
|
|
! endif |
408 |
|
|
! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) |
409 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) |
410 |
|
|
! sigma1s=ratqs(ind1,ind2)*po(ind1) |
411 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 |
412 |
|
|
|
413 |
|
|
!------------------------------------------------------------------------------ |
414 |
|
|
! Calcul de l'eau condens?e et de la couverture nuageuse |
415 |
|
|
!------------------------------------------------------------------------------ |
416 |
|
|
sqrt2pi=sqrt(2.*pi) |
417 |
|
|
xth=sth/(sqrt(2.)*sigma2s) |
418 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
419 |
|
|
cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) |
420 |
|
|
cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
421 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
422 |
|
|
|
423 |
|
|
qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) |
424 |
|
|
qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
425 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
426 |
|
|
|
427 |
|
|
IF (iflag_cloudth_vert == 1) THEN |
428 |
|
|
!------------------------------------------------------------------------------- |
429 |
|
|
! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs |
430 |
|
|
!------------------------------------------------------------------------------- |
431 |
|
|
! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) |
432 |
|
|
! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) |
433 |
|
|
deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) |
434 |
|
|
deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) |
435 |
|
|
! deltasenv=aenv*0.01*po(ind1) |
436 |
|
|
! deltasth=ath*0.01*zqta(ind1,ind2) |
437 |
|
|
xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) |
438 |
|
|
xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) |
439 |
|
|
xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) |
440 |
|
|
xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) |
441 |
|
|
coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) |
442 |
|
|
coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) |
443 |
|
|
|
444 |
|
|
cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) |
445 |
|
|
cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) |
446 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
447 |
|
|
|
448 |
|
|
IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) |
449 |
|
|
IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) |
450 |
|
|
IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) |
451 |
|
|
IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) |
452 |
|
|
|
453 |
|
|
qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
454 |
|
|
! qlenv(ind1,ind2)=IntJ |
455 |
|
|
! print*, qlenv(ind1,ind2),'VERIF EAU' |
456 |
|
|
|
457 |
|
|
|
458 |
|
|
IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) |
459 |
|
|
! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2)) |
460 |
|
|
! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1)) |
461 |
|
|
IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) |
462 |
|
|
IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) |
463 |
|
|
IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) |
464 |
|
|
qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
465 |
|
|
! qlth(ind1,ind2)=IntJ |
466 |
|
|
! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' |
467 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
468 |
|
|
|
469 |
|
|
ELSE IF (iflag_cloudth_vert == 2) THEN |
470 |
|
|
|
471 |
|
|
!------------------------------------------------------------------------------- |
472 |
|
|
! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s |
473 |
|
|
!------------------------------------------------------------------------------- |
474 |
|
|
! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) |
475 |
|
|
! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) |
476 |
|
|
! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) |
477 |
|
|
! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) |
478 |
|
|
deltasenv=aenv*vert_alpha*sigma1s |
479 |
|
|
deltasth=ath*vert_alpha*sigma2s |
480 |
|
|
|
481 |
|
|
xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) |
482 |
|
|
xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) |
483 |
|
|
xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) |
484 |
|
|
xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) |
485 |
|
|
! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) |
486 |
|
|
! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) |
487 |
|
|
|
488 |
|
|
cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) |
489 |
|
|
cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) |
490 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
491 |
|
|
|
492 |
|
|
IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2) |
493 |
|
|
IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) |
494 |
|
|
IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) |
495 |
|
|
IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2)) |
496 |
|
|
|
497 |
|
|
! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) |
498 |
|
|
! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) |
499 |
|
|
! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) |
500 |
|
|
|
501 |
|
|
qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
502 |
|
|
! qlenv(ind1,ind2)=IntJ |
503 |
|
|
! print*, qlenv(ind1,ind2),'VERIF EAU' |
504 |
|
|
|
505 |
|
|
IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2) |
506 |
|
|
IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) |
507 |
|
|
IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) |
508 |
|
|
IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2)) |
509 |
|
|
|
510 |
|
|
qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
511 |
|
|
! qlth(ind1,ind2)=IntJ |
512 |
|
|
! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' |
513 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
514 |
|
|
|
515 |
|
|
|
516 |
|
|
|
517 |
|
|
|
518 |
|
|
ENDIF ! of if (iflag_cloudth_vert==1 or 2) |
519 |
|
|
|
520 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
521 |
|
|
|
522 |
|
|
if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then |
523 |
|
|
ctot(ind1,ind2)=0. |
524 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
525 |
|
|
|
526 |
|
|
else |
527 |
|
|
|
528 |
|
|
ctot(ind1,ind2)=ctot(ind1,ind2) |
529 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
530 |
|
|
! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & |
531 |
|
|
! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) |
532 |
|
|
|
533 |
|
|
endif |
534 |
|
|
|
535 |
|
|
|
536 |
|
|
|
537 |
|
|
! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' |
538 |
|
|
|
539 |
|
|
|
540 |
|
|
else ! gaussienne environnement seule |
541 |
|
|
|
542 |
|
|
zqenv(ind1)=po(ind1) |
543 |
|
|
Tbef=t(ind1,ind2) |
544 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
545 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
546 |
|
|
qsatbef=MIN(0.5,qsatbef) |
547 |
|
|
zcor=1./(1.-retv*qsatbef) |
548 |
|
|
qsatbef=qsatbef*zcor |
549 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
550 |
|
|
|
551 |
|
|
|
552 |
|
|
! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) |
553 |
|
|
zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) |
554 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
555 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
556 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
557 |
|
|
|
558 |
|
|
|
559 |
|
|
sigma1s=ratqs(ind1,ind2)*zqenv(ind1) |
560 |
|
|
|
561 |
|
|
sqrt2pi=sqrt(2.*pi) |
562 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
563 |
|
|
ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
564 |
|
|
qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
565 |
|
|
|
566 |
|
|
if (ctot(ind1,ind2).lt.1.e-3) then |
567 |
|
|
ctot(ind1,ind2)=0. |
568 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
569 |
|
|
|
570 |
|
|
else |
571 |
|
|
|
572 |
|
|
ctot(ind1,ind2)=ctot(ind1,ind2) |
573 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) |
574 |
|
|
|
575 |
|
|
endif |
576 |
|
|
|
577 |
|
|
|
578 |
|
|
|
579 |
|
|
|
580 |
|
|
|
581 |
|
|
|
582 |
|
|
endif |
583 |
|
|
enddo |
584 |
|
|
|
585 |
|
|
return |
586 |
|
|
! end |
587 |
|
|
END SUBROUTINE cloudth_vert |
588 |
|
|
|
589 |
|
|
|
590 |
|
|
|
591 |
|
|
|
592 |
|
11232 |
SUBROUTINE cloudth_v3(ngrid,klev,ind2, & |
593 |
|
11232 |
& ztv,po,zqta,fraca, & |
594 |
|
|
& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
595 |
|
|
& ratqs,zqs,t) |
596 |
|
|
|
597 |
|
|
use lscp_ini_mod, only: iflag_cloudth_vert |
598 |
|
|
|
599 |
|
|
IMPLICIT NONE |
600 |
|
|
|
601 |
|
|
|
602 |
|
|
!=========================================================================== |
603 |
|
|
! Author : Arnaud Octavio Jam (LMD/CNRS) |
604 |
|
|
! Date : 25 Mai 2010 |
605 |
|
|
! Objet : calcule les valeurs de qc et rneb dans les thermiques |
606 |
|
|
!=========================================================================== |
607 |
|
|
|
608 |
|
|
|
609 |
|
|
INCLUDE "YOMCST.h" |
610 |
|
|
INCLUDE "YOETHF.h" |
611 |
|
|
INCLUDE "FCTTRE.h" |
612 |
|
|
|
613 |
|
|
INTEGER itap,ind1,ind2 |
614 |
|
|
INTEGER ngrid,klev,klon,l,ig |
615 |
|
|
|
616 |
|
|
REAL ztv(ngrid,klev) |
617 |
|
|
REAL po(ngrid) |
618 |
|
22464 |
REAL zqenv(ngrid) |
619 |
|
|
REAL zqta(ngrid,klev) |
620 |
|
|
|
621 |
|
|
REAL fraca(ngrid,klev+1) |
622 |
|
|
REAL zpspsk(ngrid,klev) |
623 |
|
|
REAL paprs(ngrid,klev+1) |
624 |
|
|
REAL pplay(ngrid,klev) |
625 |
|
|
REAL ztla(ngrid,klev) |
626 |
|
|
REAL zthl(ngrid,klev) |
627 |
|
|
|
628 |
|
22464 |
REAL zqsatth(ngrid,klev) |
629 |
|
22464 |
REAL zqsatenv(ngrid,klev) |
630 |
|
|
|
631 |
|
22464 |
REAL sigma1(ngrid,klev) |
632 |
|
22464 |
REAL sigma2(ngrid,klev) |
633 |
|
22464 |
REAL qlth(ngrid,klev) |
634 |
|
22464 |
REAL qlenv(ngrid,klev) |
635 |
|
22464 |
REAL qltot(ngrid,klev) |
636 |
|
22464 |
REAL cth(ngrid,klev) |
637 |
|
22464 |
REAL cenv(ngrid,klev) |
638 |
|
|
REAL ctot(ngrid,klev) |
639 |
|
22464 |
REAL cth_vol(ngrid,klev) |
640 |
|
22464 |
REAL cenv_vol(ngrid,klev) |
641 |
|
|
REAL ctot_vol(ngrid,klev) |
642 |
|
22464 |
REAL rneb(ngrid,klev) |
643 |
|
|
REAL t(ngrid,klev) |
644 |
|
|
REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi |
645 |
|
|
REAL rdd,cppd,Lv |
646 |
|
|
REAL alth,alenv,ath,aenv |
647 |
|
|
REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2 |
648 |
|
|
REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks |
649 |
|
|
REAL Tbef,zdelta,qsatbef,zcor |
650 |
|
|
REAL qlbef |
651 |
|
|
REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution |
652 |
|
|
REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) |
653 |
|
|
REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) |
654 |
|
|
REAL zqs(ngrid), qcloud(ngrid) |
655 |
|
|
REAL erf |
656 |
|
|
|
657 |
|
|
|
658 |
✓✗ |
11232 |
IF (iflag_cloudth_vert.GE.1) THEN |
659 |
|
|
CALL cloudth_vert_v3(ngrid,klev,ind2, & |
660 |
|
|
& ztv,po,zqta,fraca, & |
661 |
|
|
& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
662 |
|
11232 |
& ratqs,zqs,t) |
663 |
|
11232 |
RETURN |
664 |
|
|
ENDIF |
665 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
666 |
|
|
|
667 |
|
|
|
668 |
|
|
!------------------------------------------------------------------------------- |
669 |
|
|
! Initialisation des variables r?elles |
670 |
|
|
!------------------------------------------------------------------------------- |
671 |
|
|
sigma1(:,:)=0. |
672 |
|
|
sigma2(:,:)=0. |
673 |
|
|
qlth(:,:)=0. |
674 |
|
|
qlenv(:,:)=0. |
675 |
|
|
qltot(:,:)=0. |
676 |
|
|
rneb(:,:)=0. |
677 |
|
|
qcloud(:)=0. |
678 |
|
|
cth(:,:)=0. |
679 |
|
|
cenv(:,:)=0. |
680 |
|
|
ctot(:,:)=0. |
681 |
|
|
cth_vol(:,:)=0. |
682 |
|
|
cenv_vol(:,:)=0. |
683 |
|
|
ctot_vol(:,:)=0. |
684 |
|
|
qsatmmussig1=0. |
685 |
|
|
qsatmmussig2=0. |
686 |
|
|
rdd=287.04 |
687 |
|
|
cppd=1005.7 |
688 |
|
|
pi=3.14159 |
689 |
|
|
Lv=2.5e6 |
690 |
|
|
sqrt2pi=sqrt(2.*pi) |
691 |
|
|
sqrt2=sqrt(2.) |
692 |
|
|
sqrtpi=sqrt(pi) |
693 |
|
|
|
694 |
|
|
|
695 |
|
|
!------------------------------------------------------------------------------- |
696 |
|
|
! Cloud fraction in the thermals and standard deviation of the PDFs |
697 |
|
|
!------------------------------------------------------------------------------- |
698 |
|
|
do ind1=1,ngrid |
699 |
|
|
|
700 |
|
|
if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then |
701 |
|
|
|
702 |
|
|
zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) |
703 |
|
|
|
704 |
|
|
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
705 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
706 |
|
|
qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
707 |
|
|
qsatbef=MIN(0.5,qsatbef) |
708 |
|
|
zcor=1./(1.-retv*qsatbef) |
709 |
|
|
qsatbef=qsatbef*zcor |
710 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
711 |
|
|
|
712 |
|
|
|
713 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 |
714 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 |
715 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 |
716 |
|
|
|
717 |
|
|
!po = qt de l'environnement ET des thermique |
718 |
|
|
!zqenv = qt environnement |
719 |
|
|
!zqsatenv = qsat environnement |
720 |
|
|
!zthl = Tl environnement |
721 |
|
|
|
722 |
|
|
|
723 |
|
|
Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
724 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
725 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
726 |
|
|
qsatbef=MIN(0.5,qsatbef) |
727 |
|
|
zcor=1./(1.-retv*qsatbef) |
728 |
|
|
qsatbef=qsatbef*zcor |
729 |
|
|
zqsatth(ind1,ind2)=qsatbef |
730 |
|
|
|
731 |
|
|
alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 |
732 |
|
|
ath=1./(1.+(alth*Lv/cppd)) !al, p84 |
733 |
|
|
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 |
734 |
|
|
|
735 |
|
|
!zqta = qt thermals |
736 |
|
|
!zqsatth = qsat thermals |
737 |
|
|
!ztla = Tl thermals |
738 |
|
|
|
739 |
|
|
!------------------------------------------------------------------------------ |
740 |
|
|
! s standard deviations |
741 |
|
|
!------------------------------------------------------------------------------ |
742 |
|
|
|
743 |
|
|
! tests |
744 |
|
|
! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) |
745 |
|
|
! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) |
746 |
|
|
! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) |
747 |
|
|
! final option |
748 |
|
|
sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) |
749 |
|
|
sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) |
750 |
|
|
|
751 |
|
|
!------------------------------------------------------------------------------ |
752 |
|
|
! Condensed water and cloud cover |
753 |
|
|
!------------------------------------------------------------------------------ |
754 |
|
|
xth=sth/(sqrt2*sigma2s) |
755 |
|
|
xenv=senv/(sqrt2*sigma1s) |
756 |
|
|
cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam |
757 |
|
|
cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam |
758 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
759 |
|
|
ctot_vol(ind1,ind2)=ctot(ind1,ind2) |
760 |
|
|
|
761 |
|
|
qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) |
762 |
|
|
qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) |
763 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
764 |
|
|
|
765 |
|
|
if (ctot(ind1,ind2).lt.1.e-10) then |
766 |
|
|
ctot(ind1,ind2)=0. |
767 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
768 |
|
|
else |
769 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
770 |
|
|
endif |
771 |
|
|
|
772 |
|
|
else ! Environnement only, follow the if l.110 |
773 |
|
|
|
774 |
|
|
zqenv(ind1)=po(ind1) |
775 |
|
|
Tbef=t(ind1,ind2) |
776 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
777 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
778 |
|
|
qsatbef=MIN(0.5,qsatbef) |
779 |
|
|
zcor=1./(1.-retv*qsatbef) |
780 |
|
|
qsatbef=qsatbef*zcor |
781 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
782 |
|
|
|
783 |
|
|
! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) |
784 |
|
|
zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) |
785 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
786 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
787 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
788 |
|
|
|
789 |
|
|
sigma1s=ratqs(ind1,ind2)*zqenv(ind1) |
790 |
|
|
|
791 |
|
|
xenv=senv/(sqrt2*sigma1s) |
792 |
|
|
ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
793 |
|
|
ctot_vol(ind1,ind2)=ctot(ind1,ind2) |
794 |
|
|
qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) |
795 |
|
|
|
796 |
|
|
if (ctot(ind1,ind2).lt.1.e-3) then |
797 |
|
|
ctot(ind1,ind2)=0. |
798 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
799 |
|
|
else |
800 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) |
801 |
|
|
endif |
802 |
|
|
|
803 |
|
|
|
804 |
|
|
endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 |
805 |
|
|
enddo ! from the loop on ngrid l.108 |
806 |
|
|
return |
807 |
|
|
! end |
808 |
|
|
END SUBROUTINE cloudth_v3 |
809 |
|
|
|
810 |
|
|
|
811 |
|
|
|
812 |
|
|
!=========================================================================== |
813 |
|
11232 |
SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, & |
814 |
|
|
& ztv,po,zqta,fraca, & |
815 |
|
11232 |
& qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
816 |
|
|
& ratqs,zqs,t) |
817 |
|
|
|
818 |
|
|
!=========================================================================== |
819 |
|
|
! Auteur : Arnaud Octavio Jam (LMD/CNRS) |
820 |
|
|
! Date : 25 Mai 2010 |
821 |
|
|
! Objet : calcule les valeurs de qc et rneb dans les thermiques |
822 |
|
|
!=========================================================================== |
823 |
|
|
|
824 |
|
|
use lscp_ini_mod, only: iflag_cloudth_vert |
825 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
826 |
|
|
USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, & |
827 |
|
|
& cloudth_sigmath,cloudth_sigmaenv |
828 |
|
|
|
829 |
|
|
IMPLICIT NONE |
830 |
|
|
|
831 |
|
|
INCLUDE "YOMCST.h" |
832 |
|
|
INCLUDE "YOETHF.h" |
833 |
|
|
INCLUDE "FCTTRE.h" |
834 |
|
|
|
835 |
|
|
INTEGER itap,ind1,ind2 |
836 |
|
|
INTEGER ngrid,klev,klon,l,ig |
837 |
|
|
|
838 |
|
|
REAL ztv(ngrid,klev) |
839 |
|
|
REAL po(ngrid) |
840 |
|
22464 |
REAL zqenv(ngrid) |
841 |
|
|
REAL zqta(ngrid,klev) |
842 |
|
|
|
843 |
|
|
REAL fraca(ngrid,klev+1) |
844 |
|
|
REAL zpspsk(ngrid,klev) |
845 |
|
|
REAL paprs(ngrid,klev+1) |
846 |
|
|
REAL pplay(ngrid,klev) |
847 |
|
|
REAL ztla(ngrid,klev) |
848 |
|
|
REAL zthl(ngrid,klev) |
849 |
|
|
|
850 |
|
22464 |
REAL zqsatth(ngrid,klev) |
851 |
|
22464 |
REAL zqsatenv(ngrid,klev) |
852 |
|
|
|
853 |
|
22464 |
REAL sigma1(ngrid,klev) |
854 |
|
22464 |
REAL sigma2(ngrid,klev) |
855 |
|
22464 |
REAL qlth(ngrid,klev) |
856 |
|
22464 |
REAL qlenv(ngrid,klev) |
857 |
|
22464 |
REAL qltot(ngrid,klev) |
858 |
|
22464 |
REAL cth(ngrid,klev) |
859 |
|
22464 |
REAL cenv(ngrid,klev) |
860 |
|
|
REAL ctot(ngrid,klev) |
861 |
|
22464 |
REAL cth_vol(ngrid,klev) |
862 |
|
22464 |
REAL cenv_vol(ngrid,klev) |
863 |
|
|
REAL ctot_vol(ngrid,klev) |
864 |
|
22464 |
REAL rneb(ngrid,klev) |
865 |
|
|
REAL t(ngrid,klev) |
866 |
|
|
REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi |
867 |
|
|
REAL rdd,cppd,Lv |
868 |
|
|
REAL alth,alenv,ath,aenv |
869 |
|
|
REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs |
870 |
|
|
REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks |
871 |
|
|
REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 |
872 |
|
|
REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv |
873 |
|
|
REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth |
874 |
|
|
REAL Tbef,zdelta,qsatbef,zcor |
875 |
|
|
REAL qlbef |
876 |
|
|
REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur |
877 |
|
|
! Change the width of the PDF used for vertical subgrid scale heterogeneity |
878 |
|
|
! (J Jouhaud, JL Dufresne, JB Madeleine) |
879 |
|
|
REAL,SAVE :: vert_alpha, vert_alpha_th |
880 |
|
|
!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th) |
881 |
|
|
REAL,SAVE :: sigma1s_factor=1.1 |
882 |
|
|
REAL,SAVE :: sigma1s_power=0.6 |
883 |
|
|
REAL,SAVE :: sigma2s_factor=0.09 |
884 |
|
|
REAL,SAVE :: sigma2s_power=0.5 |
885 |
|
|
REAL,SAVE :: cloudth_ratqsmin=-1. |
886 |
|
|
!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin) |
887 |
|
|
INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0 |
888 |
|
|
!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs) |
889 |
|
|
|
890 |
|
|
LOGICAL, SAVE :: firstcall = .TRUE. |
891 |
|
|
!$OMP THREADPRIVATE(firstcall) |
892 |
|
|
|
893 |
|
|
REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) |
894 |
|
|
REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) |
895 |
|
|
REAL zqs(ngrid), qcloud(ngrid) |
896 |
|
|
REAL erf |
897 |
|
|
|
898 |
|
22464 |
REAL rhodz(ngrid,klev) |
899 |
|
22464 |
REAL zrho(ngrid,klev) |
900 |
|
11232 |
REAL dz(ngrid,klev) |
901 |
|
|
|
902 |
✓✓ |
11175840 |
DO ind1 = 1, ngrid |
903 |
|
|
!Layer calculation |
904 |
|
11164608 |
rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2 |
905 |
|
11164608 |
zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3 |
906 |
|
11232 |
dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre |
907 |
|
|
END DO |
908 |
|
|
|
909 |
|
|
!------------------------------------------------------------------------------ |
910 |
|
|
! Initialize |
911 |
|
|
!------------------------------------------------------------------------------ |
912 |
|
|
|
913 |
✓✓✓✓
|
435868992 |
sigma1(:,:)=0. |
914 |
✓✓✓✓
|
435868992 |
sigma2(:,:)=0. |
915 |
✓✓✓✓
|
435868992 |
qlth(:,:)=0. |
916 |
✓✓✓✓
|
435868992 |
qlenv(:,:)=0. |
917 |
✓✓✓✓
|
435868992 |
qltot(:,:)=0. |
918 |
✓✓✓✓
|
435868992 |
rneb(:,:)=0. |
919 |
✓✓ |
11175840 |
qcloud(:)=0. |
920 |
✓✓✓✓
|
435868992 |
cth(:,:)=0. |
921 |
✓✓✓✓
|
435868992 |
cenv(:,:)=0. |
922 |
✓✓✓✓
|
435868992 |
ctot(:,:)=0. |
923 |
✓✓✓✓
|
435868992 |
cth_vol(:,:)=0. |
924 |
✓✓✓✓
|
435868992 |
cenv_vol(:,:)=0. |
925 |
✓✓✓✓
|
435868992 |
ctot_vol(:,:)=0. |
926 |
|
|
qsatmmussig1=0. |
927 |
|
|
qsatmmussig2=0. |
928 |
|
|
rdd=287.04 |
929 |
|
|
cppd=1005.7 |
930 |
|
|
pi=3.14159 |
931 |
|
|
Lv=2.5e6 |
932 |
|
|
sqrt2pi=sqrt(2.*pi) |
933 |
|
|
sqrt2=sqrt(2.) |
934 |
|
|
sqrtpi=sqrt(pi) |
935 |
|
|
|
936 |
✓✓ |
11232 |
IF (firstcall) THEN |
937 |
|
1 |
vert_alpha=0.5 |
938 |
|
1 |
CALL getin_p('cloudth_vert_alpha',vert_alpha) |
939 |
|
1 |
WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha |
940 |
|
|
! The factor used for the thermal is equal to that of the environment |
941 |
|
|
! if nothing is explicitly specified in the def file |
942 |
|
1 |
vert_alpha_th=vert_alpha |
943 |
|
1 |
CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th) |
944 |
|
1 |
WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th |
945 |
|
|
! Factor used in the calculation of sigma1s |
946 |
|
1 |
CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor) |
947 |
|
1 |
WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor |
948 |
|
|
! Power used in the calculation of sigma1s |
949 |
|
1 |
CALL getin_p('cloudth_sigma1s_power',sigma1s_power) |
950 |
|
1 |
WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power |
951 |
|
|
! Factor used in the calculation of sigma2s |
952 |
|
1 |
CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor) |
953 |
|
1 |
WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor |
954 |
|
|
! Power used in the calculation of sigma2s |
955 |
|
1 |
CALL getin_p('cloudth_sigma2s_power',sigma2s_power) |
956 |
|
1 |
WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power |
957 |
|
|
! Minimum value for the environmental air subgrid water distrib |
958 |
|
1 |
CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin) |
959 |
|
1 |
WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin |
960 |
|
|
! Remove the dependency to ratqs from the variance of the vertical PDF |
961 |
|
1 |
CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs) |
962 |
|
1 |
WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs |
963 |
|
|
|
964 |
|
1 |
firstcall=.FALSE. |
965 |
|
|
ENDIF |
966 |
|
|
|
967 |
|
|
!------------------------------------------------------------------------------- |
968 |
|
|
! Calcul de la fraction du thermique et des ecart-types des distributions |
969 |
|
|
!------------------------------------------------------------------------------- |
970 |
✓✓ |
11175840 |
do ind1=1,ngrid |
971 |
|
|
|
972 |
✓✓✓✓
|
11164608 |
if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement |
973 |
|
|
|
974 |
|
741354 |
zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv |
975 |
|
|
|
976 |
|
|
|
977 |
|
741354 |
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
978 |
|
741354 |
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
979 |
|
741354 |
qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
980 |
|
741354 |
qsatbef=MIN(0.5,qsatbef) |
981 |
|
741354 |
zcor=1./(1.-retv*qsatbef) |
982 |
|
741354 |
qsatbef=qsatbef*zcor |
983 |
|
741354 |
zqsatenv(ind1,ind2)=qsatbef |
984 |
|
|
|
985 |
|
|
|
986 |
|
741354 |
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 |
987 |
|
741354 |
aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 |
988 |
|
741354 |
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 |
989 |
|
|
|
990 |
|
|
!zqenv = qt environnement |
991 |
|
|
!zqsatenv = qsat environnement |
992 |
|
|
!zthl = Tl environnement |
993 |
|
|
|
994 |
|
|
|
995 |
|
741354 |
Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
996 |
|
741354 |
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
997 |
|
741354 |
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
998 |
|
741354 |
qsatbef=MIN(0.5,qsatbef) |
999 |
|
741354 |
zcor=1./(1.-retv*qsatbef) |
1000 |
|
741354 |
qsatbef=qsatbef*zcor |
1001 |
|
741354 |
zqsatth(ind1,ind2)=qsatbef |
1002 |
|
|
|
1003 |
|
741354 |
alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 |
1004 |
|
741354 |
ath=1./(1.+(alth*Lv/cppd)) !al, p84 |
1005 |
|
741354 |
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 |
1006 |
|
|
|
1007 |
|
|
|
1008 |
|
|
!zqta = qt thermals |
1009 |
|
|
!zqsatth = qsat thermals |
1010 |
|
|
!ztla = Tl thermals |
1011 |
|
|
!------------------------------------------------------------------------------ |
1012 |
|
|
! s standard deviation |
1013 |
|
|
!------------------------------------------------------------------------------ |
1014 |
|
|
|
1015 |
|
|
sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & |
1016 |
|
741354 |
& (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 |
1017 |
|
|
! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 |
1018 |
|
741354 |
IF (cloudth_ratqsmin>0.) THEN |
1019 |
|
|
sigma1s_ratqs = cloudth_ratqsmin*po(ind1) |
1020 |
|
|
ELSE |
1021 |
|
741354 |
sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) |
1022 |
|
|
ENDIF |
1023 |
|
741354 |
sigma1s = sigma1s_fraca + sigma1s_ratqs |
1024 |
|
741354 |
sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) |
1025 |
|
|
! tests |
1026 |
|
|
! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) |
1027 |
|
|
! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) |
1028 |
|
|
! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) |
1029 |
|
|
! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) |
1030 |
|
|
! if (paprs(ind1,ind2).gt.90000) then |
1031 |
|
|
! ratqs(ind1,ind2)=0.002 |
1032 |
|
|
! else |
1033 |
|
|
! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 |
1034 |
|
|
! endif |
1035 |
|
|
! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) |
1036 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) |
1037 |
|
|
! sigma1s=ratqs(ind1,ind2)*po(ind1) |
1038 |
|
|
! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 |
1039 |
|
|
|
1040 |
✗✓ |
741354 |
IF (iflag_cloudth_vert == 1) THEN |
1041 |
|
|
!------------------------------------------------------------------------------- |
1042 |
|
|
! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs |
1043 |
|
|
!------------------------------------------------------------------------------- |
1044 |
|
|
|
1045 |
|
|
deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) |
1046 |
|
|
deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) |
1047 |
|
|
|
1048 |
|
|
xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) |
1049 |
|
|
xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) |
1050 |
|
|
xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) |
1051 |
|
|
xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) |
1052 |
|
|
coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) |
1053 |
|
|
coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) |
1054 |
|
|
|
1055 |
|
|
cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) |
1056 |
|
|
cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) |
1057 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
1058 |
|
|
|
1059 |
|
|
! Environment |
1060 |
|
|
IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) |
1061 |
|
|
IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) |
1062 |
|
|
IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) |
1063 |
|
|
IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) |
1064 |
|
|
|
1065 |
|
|
qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1066 |
|
|
|
1067 |
|
|
! Thermal |
1068 |
|
|
IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) |
1069 |
|
|
IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) |
1070 |
|
|
IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) |
1071 |
|
|
IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) |
1072 |
|
|
qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1073 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
1074 |
|
|
|
1075 |
✓✗ |
741354 |
ELSE IF (iflag_cloudth_vert >= 3) THEN |
1076 |
✓✗ |
741354 |
IF (iflag_cloudth_vert < 5) THEN |
1077 |
|
|
!------------------------------------------------------------------------------- |
1078 |
|
|
! Version 3: Changes by J. Jouhaud; condensation for q > -delta s |
1079 |
|
|
!------------------------------------------------------------------------------- |
1080 |
|
|
! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) |
1081 |
|
|
! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) |
1082 |
|
|
! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) |
1083 |
|
|
! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) |
1084 |
✓✗ |
741354 |
IF (iflag_cloudth_vert == 3) THEN |
1085 |
|
741354 |
deltasenv=aenv*vert_alpha*sigma1s |
1086 |
|
741354 |
deltasth=ath*vert_alpha_th*sigma2s |
1087 |
|
|
ELSE IF (iflag_cloudth_vert == 4) THEN |
1088 |
|
|
IF (iflag_cloudth_vert_noratqs == 1) THEN |
1089 |
|
|
deltasenv=vert_alpha*max(sigma1s_fraca,1e-10) |
1090 |
|
|
deltasth=vert_alpha_th*sigma2s |
1091 |
|
|
ELSE |
1092 |
|
|
deltasenv=vert_alpha*sigma1s |
1093 |
|
|
deltasth=vert_alpha_th*sigma2s |
1094 |
|
|
ENDIF |
1095 |
|
|
ENDIF |
1096 |
|
|
|
1097 |
|
741354 |
xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) |
1098 |
|
741354 |
xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) |
1099 |
|
741354 |
exp_xenv1 = exp(-1.*xenv1**2) |
1100 |
|
741354 |
exp_xenv2 = exp(-1.*xenv2**2) |
1101 |
|
741354 |
xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) |
1102 |
|
741354 |
xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) |
1103 |
|
741354 |
exp_xth1 = exp(-1.*xth1**2) |
1104 |
|
741354 |
exp_xth2 = exp(-1.*xth2**2) |
1105 |
|
|
|
1106 |
|
|
!CF_surfacique |
1107 |
|
741354 |
cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) |
1108 |
|
741354 |
cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) |
1109 |
|
741354 |
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
1110 |
|
|
|
1111 |
|
|
|
1112 |
|
|
!CF_volumique & eau condense |
1113 |
|
|
!environnement |
1114 |
|
741354 |
IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 |
1115 |
|
741354 |
IntJ_CF=0.5*(1.-1.*erf(xenv2)) |
1116 |
✗✓ |
741354 |
if (deltasenv .lt. 1.e-10) then |
1117 |
|
|
qlenv(ind1,ind2)=IntJ |
1118 |
|
|
cenv_vol(ind1,ind2)=IntJ_CF |
1119 |
|
|
else |
1120 |
|
741354 |
IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) |
1121 |
|
741354 |
IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) |
1122 |
|
741354 |
IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) |
1123 |
|
741354 |
IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) |
1124 |
|
741354 |
IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) |
1125 |
|
741354 |
qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1126 |
|
741354 |
cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
1127 |
|
|
endif |
1128 |
|
|
|
1129 |
|
|
!thermique |
1130 |
|
741354 |
IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 |
1131 |
|
741354 |
IntJ_CF=0.5*(1.-1.*erf(xth2)) |
1132 |
✗✓ |
741354 |
if (deltasth .lt. 1.e-10) then |
1133 |
|
|
qlth(ind1,ind2)=IntJ |
1134 |
|
|
cth_vol(ind1,ind2)=IntJ_CF |
1135 |
|
|
else |
1136 |
|
741354 |
IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) |
1137 |
|
741354 |
IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) |
1138 |
|
741354 |
IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) |
1139 |
|
741354 |
IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) |
1140 |
|
741354 |
IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) |
1141 |
|
741354 |
qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1142 |
|
741354 |
cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
1143 |
|
|
endif |
1144 |
|
|
|
1145 |
|
741354 |
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
1146 |
|
741354 |
ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) |
1147 |
|
|
|
1148 |
|
|
ELSE IF (iflag_cloudth_vert == 5) THEN |
1149 |
|
|
sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & |
1150 |
|
|
/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & |
1151 |
|
|
+ratqs(ind1,ind2)*po(ind1) !Environment |
1152 |
|
|
sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals |
1153 |
|
|
!sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) |
1154 |
|
|
!sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) |
1155 |
|
|
xth=sth/(sqrt(2.)*sigma2s) |
1156 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
1157 |
|
|
|
1158 |
|
|
!Volumique |
1159 |
|
|
cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) |
1160 |
|
|
cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
1161 |
|
|
ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) |
1162 |
|
|
!print *,'jeanjean_CV=',ctot_vol(ind1,ind2) |
1163 |
|
|
|
1164 |
|
|
qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2)) |
1165 |
|
|
qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) |
1166 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
1167 |
|
|
|
1168 |
|
|
!Surfacique |
1169 |
|
|
!Neggers |
1170 |
|
|
!beta=0.0044 |
1171 |
|
|
!inverse_rho=1.+beta*dz(ind1,ind2) |
1172 |
|
|
!print *,'jeanjean : beta=',beta |
1173 |
|
|
!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho |
1174 |
|
|
!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho |
1175 |
|
|
!ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
1176 |
|
|
|
1177 |
|
|
!Brooks |
1178 |
|
|
a_Brooks=0.6694 |
1179 |
|
|
b_Brooks=0.1882 |
1180 |
|
|
A_Maj_Brooks=0.1635 !-- sans shear |
1181 |
|
|
!A_Maj_Brooks=0.17 !-- ARM LES |
1182 |
|
|
!A_Maj_Brooks=0.18 !-- RICO LES |
1183 |
|
|
!A_Maj_Brooks=0.19 !-- BOMEX LES |
1184 |
|
|
Dx_Brooks=200000. |
1185 |
|
|
f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) |
1186 |
|
|
!print *,'jeanjean_f=',f_Brooks |
1187 |
|
|
|
1188 |
|
|
cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.)) |
1189 |
|
|
cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.)) |
1190 |
|
|
ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) |
1191 |
|
|
!print *,'JJ_ctot_1',ctot(ind1,ind2) |
1192 |
|
|
|
1193 |
|
|
|
1194 |
|
|
|
1195 |
|
|
|
1196 |
|
|
|
1197 |
|
|
ENDIF ! of if (iflag_cloudth_vert<5) |
1198 |
|
|
ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) |
1199 |
|
|
|
1200 |
|
|
! if (ctot(ind1,ind2).lt.1.e-10) then |
1201 |
✓✓✓✓
|
741354 |
if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then |
1202 |
|
455512 |
ctot(ind1,ind2)=0. |
1203 |
|
455512 |
ctot_vol(ind1,ind2)=0. |
1204 |
|
455512 |
qcloud(ind1)=zqsatenv(ind1,ind2) |
1205 |
|
|
|
1206 |
|
|
else |
1207 |
|
|
|
1208 |
|
285842 |
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
1209 |
|
|
! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & |
1210 |
|
|
! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) |
1211 |
|
|
|
1212 |
|
|
endif |
1213 |
|
|
|
1214 |
|
|
else ! gaussienne environnement seule |
1215 |
|
|
|
1216 |
|
|
|
1217 |
|
10423254 |
zqenv(ind1)=po(ind1) |
1218 |
|
10423254 |
Tbef=t(ind1,ind2) |
1219 |
|
10423254 |
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
1220 |
|
10423254 |
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
1221 |
|
10423254 |
qsatbef=MIN(0.5,qsatbef) |
1222 |
|
10423254 |
zcor=1./(1.-retv*qsatbef) |
1223 |
|
10423254 |
qsatbef=qsatbef*zcor |
1224 |
|
10423254 |
zqsatenv(ind1,ind2)=qsatbef |
1225 |
|
|
|
1226 |
|
|
|
1227 |
|
|
! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) |
1228 |
|
10423254 |
zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) |
1229 |
|
10423254 |
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
1230 |
|
10423254 |
aenv=1./(1.+(alenv*Lv/cppd)) |
1231 |
|
10423254 |
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
1232 |
|
|
sth=0. |
1233 |
|
|
|
1234 |
|
|
|
1235 |
|
10423254 |
sigma1s=ratqs(ind1,ind2)*zqenv(ind1) |
1236 |
|
|
sigma2s=0. |
1237 |
|
|
|
1238 |
|
|
sqrt2pi=sqrt(2.*pi) |
1239 |
|
10423254 |
xenv=senv/(sqrt(2.)*sigma1s) |
1240 |
|
10423254 |
ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
1241 |
|
10423254 |
ctot_vol(ind1,ind2)=ctot(ind1,ind2) |
1242 |
|
10423254 |
qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
1243 |
|
|
|
1244 |
|
10423254 |
if (ctot(ind1,ind2).lt.1.e-3) then |
1245 |
|
8258508 |
ctot(ind1,ind2)=0. |
1246 |
|
8258508 |
qcloud(ind1)=zqsatenv(ind1,ind2) |
1247 |
|
|
|
1248 |
|
|
else |
1249 |
|
|
|
1250 |
|
|
! ctot(ind1,ind2)=ctot(ind1,ind2) |
1251 |
|
2164746 |
qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) |
1252 |
|
|
|
1253 |
|
|
endif |
1254 |
|
|
|
1255 |
|
|
|
1256 |
|
|
|
1257 |
|
|
|
1258 |
|
|
endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 |
1259 |
|
|
! Outputs used to check the PDFs |
1260 |
|
11164608 |
cloudth_senv(ind1,ind2) = senv |
1261 |
|
11164608 |
cloudth_sth(ind1,ind2) = sth |
1262 |
|
11164608 |
cloudth_sigmaenv(ind1,ind2) = sigma1s |
1263 |
|
11175840 |
cloudth_sigmath(ind1,ind2) = sigma2s |
1264 |
|
|
|
1265 |
|
|
enddo ! from the loop on ngrid l.333 |
1266 |
|
11232 |
return |
1267 |
|
|
! end |
1268 |
|
|
END SUBROUTINE cloudth_vert_v3 |
1269 |
|
|
! |
1270 |
|
|
|
1271 |
|
|
|
1272 |
|
|
|
1273 |
|
|
|
1274 |
|
|
|
1275 |
|
|
|
1276 |
|
|
|
1277 |
|
|
|
1278 |
|
|
|
1279 |
|
|
|
1280 |
|
|
|
1281 |
|
|
SUBROUTINE cloudth_v6(ngrid,klev,ind2, & |
1282 |
|
|
& ztv,po,zqta,fraca, & |
1283 |
|
|
& qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
1284 |
|
|
& ratqs,zqs,T) |
1285 |
|
|
|
1286 |
✗✓✓✓
|
22329216 |
use lscp_ini_mod, only: iflag_cloudth_vert |
1287 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
1288 |
|
|
USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, & |
1289 |
|
|
& cloudth_sigmath,cloudth_sigmaenv |
1290 |
|
|
|
1291 |
|
|
IMPLICIT NONE |
1292 |
|
|
|
1293 |
|
|
INCLUDE "YOMCST.h" |
1294 |
|
|
INCLUDE "YOETHF.h" |
1295 |
|
|
INCLUDE "FCTTRE.h" |
1296 |
|
|
|
1297 |
|
|
|
1298 |
|
|
!Domain variables |
1299 |
|
|
INTEGER ngrid !indice Max lat-lon |
1300 |
|
|
INTEGER klev !indice Max alt |
1301 |
|
|
INTEGER ind1 !indice in [1:ngrid] |
1302 |
|
|
INTEGER ind2 !indice in [1:klev] |
1303 |
|
|
!thermal plume fraction |
1304 |
|
|
REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox |
1305 |
|
|
!temperatures |
1306 |
|
|
REAL T(ngrid,klev) !temperature |
1307 |
|
|
REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables) |
1308 |
|
|
REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90) |
1309 |
|
|
REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th) |
1310 |
|
|
REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env) |
1311 |
|
|
!pressure |
1312 |
|
|
REAL paprs(ngrid,klev+1) !pressure at the interface of levels |
1313 |
|
|
REAL pplay(ngrid,klev) !pressure at the middle of the level |
1314 |
|
|
!humidity |
1315 |
|
|
REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution |
1316 |
|
|
REAL po(ngrid) !total water (qt) |
1317 |
|
|
REAL zqenv(ngrid) !total water in the environment (qt_env) |
1318 |
|
|
REAL zqta(ngrid,klev) !total water in the thermals (qt_th) |
1319 |
|
|
REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th) |
1320 |
|
|
REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env) |
1321 |
|
|
REAL qlth(ngrid,klev) !condensed water in the thermals |
1322 |
|
|
REAL qlenv(ngrid,klev) !condensed water in the environment |
1323 |
|
|
REAL qltot(ngrid,klev) !condensed water in the gridbox |
1324 |
|
|
!cloud fractions |
1325 |
|
|
REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals |
1326 |
|
|
REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment |
1327 |
|
|
REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox |
1328 |
|
|
REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals |
1329 |
|
|
REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment |
1330 |
|
|
REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox |
1331 |
|
|
!PDF of saturation deficit variables |
1332 |
|
|
REAL rdd,cppd,Lv |
1333 |
|
|
REAL Tbef,zdelta,qsatbef,zcor |
1334 |
|
|
REAL alth,alenv,ath,aenv |
1335 |
|
|
REAL sth,senv !saturation deficits in the thermals and environment |
1336 |
|
|
REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF |
1337 |
|
|
!cloud fraction variables |
1338 |
|
|
REAL xth,xenv |
1339 |
|
|
REAL inverse_rho,beta !Neggers et al. (2011) method |
1340 |
|
|
REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method |
1341 |
|
|
!Incloud total water variables |
1342 |
|
|
REAL zqs(ngrid) !q_sat |
1343 |
|
|
REAL qcloud(ngrid) !eau totale dans le nuage |
1344 |
|
|
!Some arithmetic variables |
1345 |
|
|
REAL erf,pi,sqrt2,sqrt2pi |
1346 |
|
|
!Depth of the layer |
1347 |
|
|
REAL dz(ngrid,klev) !epaisseur de la couche en metre |
1348 |
|
|
REAL rhodz(ngrid,klev) |
1349 |
|
|
REAL zrho(ngrid,klev) |
1350 |
|
|
DO ind1 = 1, ngrid |
1351 |
|
|
rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2] |
1352 |
|
|
zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3] |
1353 |
|
|
dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m] |
1354 |
|
|
END DO |
1355 |
|
|
|
1356 |
|
|
!------------------------------------------------------------------------------ |
1357 |
|
|
! Initialization |
1358 |
|
|
!------------------------------------------------------------------------------ |
1359 |
|
|
qlth(:,:)=0. |
1360 |
|
|
qlenv(:,:)=0. |
1361 |
|
|
qltot(:,:)=0. |
1362 |
|
|
cth_vol(:,:)=0. |
1363 |
|
|
cenv_vol(:,:)=0. |
1364 |
|
|
ctot_vol(:,:)=0. |
1365 |
|
|
cth_surf(:,:)=0. |
1366 |
|
|
cenv_surf(:,:)=0. |
1367 |
|
|
ctot_surf(:,:)=0. |
1368 |
|
|
qcloud(:)=0. |
1369 |
|
|
rdd=287.04 |
1370 |
|
|
cppd=1005.7 |
1371 |
|
|
pi=3.14159 |
1372 |
|
|
Lv=2.5e6 |
1373 |
|
|
sqrt2=sqrt(2.) |
1374 |
|
|
sqrt2pi=sqrt(2.*pi) |
1375 |
|
|
|
1376 |
|
|
|
1377 |
|
|
DO ind1=1,ngrid |
1378 |
|
|
!------------------------------------------------------------------------------- |
1379 |
|
|
!Both thermal and environment in the gridbox |
1380 |
|
|
!------------------------------------------------------------------------------- |
1381 |
|
|
IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN |
1382 |
|
|
!-------------------------------------------- |
1383 |
|
|
!calcul de qsat_env |
1384 |
|
|
!-------------------------------------------- |
1385 |
|
|
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
1386 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
1387 |
|
|
qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
1388 |
|
|
qsatbef=MIN(0.5,qsatbef) |
1389 |
|
|
zcor=1./(1.-retv*qsatbef) |
1390 |
|
|
qsatbef=qsatbef*zcor |
1391 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
1392 |
|
|
!-------------------------------------------- |
1393 |
|
|
!calcul de s_env |
1394 |
|
|
!-------------------------------------------- |
1395 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam |
1396 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam |
1397 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam |
1398 |
|
|
!-------------------------------------------- |
1399 |
|
|
!calcul de qsat_th |
1400 |
|
|
!-------------------------------------------- |
1401 |
|
|
Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
1402 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
1403 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
1404 |
|
|
qsatbef=MIN(0.5,qsatbef) |
1405 |
|
|
zcor=1./(1.-retv*qsatbef) |
1406 |
|
|
qsatbef=qsatbef*zcor |
1407 |
|
|
zqsatth(ind1,ind2)=qsatbef |
1408 |
|
|
!-------------------------------------------- |
1409 |
|
|
!calcul de s_th |
1410 |
|
|
!-------------------------------------------- |
1411 |
|
|
alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 these Arnaud Jam |
1412 |
|
|
ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam |
1413 |
|
|
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam |
1414 |
|
|
!-------------------------------------------- |
1415 |
|
|
!calcul standard deviations bi-Gaussian PDF |
1416 |
|
|
!-------------------------------------------- |
1417 |
|
|
sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) |
1418 |
|
|
sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & |
1419 |
|
|
/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & |
1420 |
|
|
+ratqs(ind1,ind2)*po(ind1) |
1421 |
|
|
xth=sth/(sqrt2*sigma_th) |
1422 |
|
|
xenv=senv/(sqrt2*sigma_env) |
1423 |
|
|
!-------------------------------------------- |
1424 |
|
|
!Cloud fraction by volume CF_vol |
1425 |
|
|
!-------------------------------------------- |
1426 |
|
|
cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) |
1427 |
|
|
cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
1428 |
|
|
ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) |
1429 |
|
|
!-------------------------------------------- |
1430 |
|
|
!Condensed water qc |
1431 |
|
|
!-------------------------------------------- |
1432 |
|
|
qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2)) |
1433 |
|
|
qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) |
1434 |
|
|
qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
1435 |
|
|
!-------------------------------------------- |
1436 |
|
|
!Cloud fraction by surface CF_surf |
1437 |
|
|
!-------------------------------------------- |
1438 |
|
|
!Method Neggers et al. (2011) : ok for cumulus clouds only |
1439 |
|
|
!beta=0.0044 (Jouhaud et al.2018) |
1440 |
|
|
!inverse_rho=1.+beta*dz(ind1,ind2) |
1441 |
|
|
!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho |
1442 |
|
|
!Method Brooks et al. (2005) : ok for all types of clouds |
1443 |
|
|
a_Brooks=0.6694 |
1444 |
|
|
b_Brooks=0.1882 |
1445 |
|
|
A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent |
1446 |
|
|
Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote |
1447 |
|
|
f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) |
1448 |
|
|
ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) |
1449 |
|
|
!-------------------------------------------- |
1450 |
|
|
!Incloud Condensed water qcloud |
1451 |
|
|
!-------------------------------------------- |
1452 |
|
|
if (ctot_surf(ind1,ind2) .lt. 1.e-10) then |
1453 |
|
|
ctot_vol(ind1,ind2)=0. |
1454 |
|
|
ctot_surf(ind1,ind2)=0. |
1455 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
1456 |
|
|
else |
1457 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1) |
1458 |
|
|
endif |
1459 |
|
|
|
1460 |
|
|
|
1461 |
|
|
|
1462 |
|
|
!------------------------------------------------------------------------------- |
1463 |
|
|
!Environment only in the gridbox |
1464 |
|
|
!------------------------------------------------------------------------------- |
1465 |
|
|
ELSE |
1466 |
|
|
!-------------------------------------------- |
1467 |
|
|
!calcul de qsat_env |
1468 |
|
|
!-------------------------------------------- |
1469 |
|
|
Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
1470 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
1471 |
|
|
qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
1472 |
|
|
qsatbef=MIN(0.5,qsatbef) |
1473 |
|
|
zcor=1./(1.-retv*qsatbef) |
1474 |
|
|
qsatbef=qsatbef*zcor |
1475 |
|
|
zqsatenv(ind1,ind2)=qsatbef |
1476 |
|
|
!-------------------------------------------- |
1477 |
|
|
!calcul de s_env |
1478 |
|
|
!-------------------------------------------- |
1479 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam |
1480 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam |
1481 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam |
1482 |
|
|
!-------------------------------------------- |
1483 |
|
|
!calcul standard deviations Gaussian PDF |
1484 |
|
|
!-------------------------------------------- |
1485 |
|
|
zqenv(ind1)=po(ind1) |
1486 |
|
|
sigma_env=ratqs(ind1,ind2)*zqenv(ind1) |
1487 |
|
|
xenv=senv/(sqrt2*sigma_env) |
1488 |
|
|
!-------------------------------------------- |
1489 |
|
|
!Cloud fraction by volume CF_vol |
1490 |
|
|
!-------------------------------------------- |
1491 |
|
|
ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
1492 |
|
|
!-------------------------------------------- |
1493 |
|
|
!Condensed water qc |
1494 |
|
|
!-------------------------------------------- |
1495 |
|
|
qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2)) |
1496 |
|
|
!-------------------------------------------- |
1497 |
|
|
!Cloud fraction by surface CF_surf |
1498 |
|
|
!-------------------------------------------- |
1499 |
|
|
!Method Neggers et al. (2011) : ok for cumulus clouds only |
1500 |
|
|
!beta=0.0044 (Jouhaud et al.2018) |
1501 |
|
|
!inverse_rho=1.+beta*dz(ind1,ind2) |
1502 |
|
|
!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho |
1503 |
|
|
!Method Brooks et al. (2005) : ok for all types of clouds |
1504 |
|
|
a_Brooks=0.6694 |
1505 |
|
|
b_Brooks=0.1882 |
1506 |
|
|
A_Maj_Brooks=0.1635 !-- sans dependence au shear |
1507 |
|
|
Dx_Brooks=200000. |
1508 |
|
|
f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) |
1509 |
|
|
ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) |
1510 |
|
|
!-------------------------------------------- |
1511 |
|
|
!Incloud Condensed water qcloud |
1512 |
|
|
!-------------------------------------------- |
1513 |
|
|
if (ctot_surf(ind1,ind2) .lt. 1.e-8) then |
1514 |
|
|
ctot_vol(ind1,ind2)=0. |
1515 |
|
|
ctot_surf(ind1,ind2)=0. |
1516 |
|
|
qcloud(ind1)=zqsatenv(ind1,ind2) |
1517 |
|
|
else |
1518 |
|
|
qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2) |
1519 |
|
|
endif |
1520 |
|
|
|
1521 |
|
|
|
1522 |
|
|
END IF ! From the separation (thermal/envrionnement) et (environnement only) |
1523 |
|
|
|
1524 |
|
|
! Outputs used to check the PDFs |
1525 |
|
|
cloudth_senv(ind1,ind2) = senv |
1526 |
|
|
cloudth_sth(ind1,ind2) = sth |
1527 |
|
|
cloudth_sigmaenv(ind1,ind2) = sigma_env |
1528 |
|
|
cloudth_sigmath(ind1,ind2) = sigma_th |
1529 |
|
|
|
1530 |
|
|
END DO ! From the loop on ngrid |
1531 |
|
|
return |
1532 |
|
|
|
1533 |
|
|
END SUBROUTINE cloudth_v6 |
1534 |
|
|
|
1535 |
|
|
|
1536 |
|
|
|
1537 |
|
|
|
1538 |
|
|
|
1539 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1540 |
|
|
SUBROUTINE cloudth_mpc(klon,klev,ind2,iflag_mpc_bl,mpc_bl_points, & |
1541 |
|
|
& temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, & |
1542 |
|
|
& ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol) |
1543 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1544 |
|
|
! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS) |
1545 |
|
|
! Date: Adapted from cloudth_vert_v3 in 2021 |
1546 |
|
|
! Aim : computes qc and rneb in thermals with cold microphysical considerations |
1547 |
|
|
! + for mixed phase boundary layer clouds, calculate ql and qi from |
1548 |
|
|
! a stationary MPC model |
1549 |
|
|
! IMPORTANT NOTE: we assume iflag_clouth_vert=3 |
1550 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1551 |
|
|
|
1552 |
|
|
use lscp_ini_mod, only: iflag_cloudth_vert |
1553 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
1554 |
|
|
USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv |
1555 |
|
|
USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY |
1556 |
|
|
USE phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth |
1557 |
|
|
|
1558 |
|
|
IMPLICIT NONE |
1559 |
|
|
|
1560 |
|
|
INCLUDE "YOMCST.h" |
1561 |
|
|
INCLUDE "YOETHF.h" |
1562 |
|
|
INCLUDE "FCTTRE.h" |
1563 |
|
|
|
1564 |
|
|
|
1565 |
|
|
!------------------------------------------------------------------------------ |
1566 |
|
|
! Declaration |
1567 |
|
|
!------------------------------------------------------------------------------ |
1568 |
|
|
|
1569 |
|
|
! INPUT/OUTPUT |
1570 |
|
|
|
1571 |
|
|
INTEGER, INTENT(IN) :: klon,klev,ind2 |
1572 |
|
|
|
1573 |
|
|
|
1574 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K] |
1575 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! Virtual potential temp [K] |
1576 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: po ! specific humidity [kg/kg] |
1577 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] |
1578 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: fraca ! Fraction of the mesh covered by thermals [0-1] |
1579 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk |
1580 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! Pressure at layer interfaces [Pa] |
1581 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! Pressure at the center of layers [Pa] |
1582 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! Liquid temp [K] |
1583 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! Liquid potential temp [K] |
1584 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib. |
1585 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: zqs ! Saturation specific humidity in the mesh [kg/kg] |
1586 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s] |
1587 |
|
|
INTEGER, INTENT(IN) :: iflag_mpc_bl ! option flag for mpc boundary layer clouds param. |
1588 |
|
|
|
1589 |
|
|
|
1590 |
|
|
INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds |
1591 |
|
|
|
1592 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1] |
1593 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1] |
1594 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg] |
1595 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg] |
1596 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1] |
1597 |
|
|
|
1598 |
|
|
|
1599 |
|
|
! LOCAL VARIABLES |
1600 |
|
|
|
1601 |
|
|
INTEGER itap,ind1,l,ig,iter,k |
1602 |
|
|
INTEGER iflag_topthermals, niter |
1603 |
|
|
LOGICAL falseklon(klon) |
1604 |
|
|
|
1605 |
|
|
REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon) |
1606 |
|
|
REAL sigma1(klon,klev) |
1607 |
|
|
REAL sigma2(klon,klev) |
1608 |
|
|
REAL qcth(klon,klev) |
1609 |
|
|
REAL qcenv(klon,klev) |
1610 |
|
|
REAL qctot(klon,klev) |
1611 |
|
|
REAL cth(klon,klev) |
1612 |
|
|
REAL cenv(klon,klev) |
1613 |
|
|
REAL cth_vol(klon,klev) |
1614 |
|
|
REAL cenv_vol(klon,klev) |
1615 |
|
|
REAL rneb(klon,klev) |
1616 |
|
|
REAL zqenv(klon), zqth(klon), zqenvonly(klon) |
1617 |
|
|
REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi |
1618 |
|
|
REAL rdd,cppd,Lv |
1619 |
|
|
REAL alth,alenv,ath,aenv |
1620 |
|
|
REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs |
1621 |
|
|
REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks |
1622 |
|
|
REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 |
1623 |
|
|
REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv |
1624 |
|
|
REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth |
1625 |
|
|
REAL zdelta,qsatbef,zcor |
1626 |
|
|
REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon) |
1627 |
|
|
REAL qlbef |
1628 |
|
|
REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon) |
1629 |
|
|
REAL erf |
1630 |
|
|
REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) |
1631 |
|
|
REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) |
1632 |
|
|
REAL rhodz(klon,klev) |
1633 |
|
|
REAL zrho(klon,klev) |
1634 |
|
|
REAL dz(klon,klev) |
1635 |
|
|
REAL qslenv(klon), qslth(klon) |
1636 |
|
|
REAL alenvl, aenvl |
1637 |
|
|
REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc |
1638 |
|
|
REAL senvi, senvl, qbase, sbase, qliqth, qiceth |
1639 |
|
|
REAL qmax, ttarget, stmp, cout, coutref, fraci |
1640 |
|
|
REAL maxi, mini, pas, temp_lim |
1641 |
|
|
REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1) |
1642 |
|
|
|
1643 |
|
|
! Modifty the saturation deficit PDF in thermals |
1644 |
|
|
! in the presence of ice crystals |
1645 |
|
|
REAL,SAVE :: C_mpc |
1646 |
|
|
!$OMP THREADPRIVATE(C_mpc) |
1647 |
|
|
REAL, SAVE :: Ni,C_cap,Ei,d_top |
1648 |
|
|
!$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top) |
1649 |
|
|
! Change the width of the PDF used for vertical subgrid scale heterogeneity |
1650 |
|
|
! (J Jouhaud, JL Dufresne, JB Madeleine) |
1651 |
|
|
REAL,SAVE :: vert_alpha, vert_alpha_th |
1652 |
|
|
!$OMP THREADPRIVATE(vert_alpha, vert_alpha_th) |
1653 |
|
|
REAL,SAVE :: sigma1s_factor=1.1 |
1654 |
|
|
REAL,SAVE :: sigma1s_power=0.6 |
1655 |
|
|
REAL,SAVE :: sigma2s_factor=0.09 |
1656 |
|
|
REAL,SAVE :: sigma2s_power=0.5 |
1657 |
|
|
REAL,SAVE :: cloudth_ratqsmin=-1. |
1658 |
|
|
!$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin) |
1659 |
|
|
INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0 |
1660 |
|
|
!$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs) |
1661 |
|
|
LOGICAL, SAVE :: firstcall = .TRUE. |
1662 |
|
|
!$OMP THREADPRIVATE(firstcall) |
1663 |
|
|
|
1664 |
|
|
CHARACTER (len = 80) :: abort_message |
1665 |
|
|
CHARACTER (len = 20) :: routname = 'cloudth_mpc' |
1666 |
|
|
|
1667 |
|
|
|
1668 |
|
|
!------------------------------------------------------------------------------ |
1669 |
|
|
! Initialisation |
1670 |
|
|
!------------------------------------------------------------------------------ |
1671 |
|
|
|
1672 |
|
|
|
1673 |
|
|
! Few initial checksS |
1674 |
|
|
|
1675 |
|
|
IF (iflag_cloudth_vert.NE.3) THEN |
1676 |
|
|
abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3' |
1677 |
|
|
CALL abort_physic(routname,abort_message,1) |
1678 |
|
|
ENDIF |
1679 |
|
|
|
1680 |
|
|
DO k = 1,klev |
1681 |
|
|
DO ind1 = 1, klon |
1682 |
|
|
rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2 |
1683 |
|
|
zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3 |
1684 |
|
|
dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre |
1685 |
|
|
END DO |
1686 |
|
|
END DO |
1687 |
|
|
|
1688 |
|
|
|
1689 |
|
|
sigma1(:,:)=0. |
1690 |
|
|
sigma2(:,:)=0. |
1691 |
|
|
qcth(:,:)=0. |
1692 |
|
|
qcenv(:,:)=0. |
1693 |
|
|
qctot(:,:)=0. |
1694 |
|
|
qlth(:,ind2)=0. |
1695 |
|
|
qith(:,ind2)=0. |
1696 |
|
|
wiceth(:,ind2)=0. |
1697 |
|
|
rneb(:,:)=0. |
1698 |
|
|
qcloud(:)=0. |
1699 |
|
|
cth(:,:)=0. |
1700 |
|
|
cenv(:,:)=0. |
1701 |
|
|
ctot(:,:)=0. |
1702 |
|
|
cth_vol(:,:)=0. |
1703 |
|
|
cenv_vol(:,:)=0. |
1704 |
|
|
ctot_vol(:,:)=0. |
1705 |
|
|
falseklon(:)=.false. |
1706 |
|
|
qsatmmussig1=0. |
1707 |
|
|
qsatmmussig2=0. |
1708 |
|
|
rdd=287.04 |
1709 |
|
|
cppd=1005.7 |
1710 |
|
|
pi=3.14159 |
1711 |
|
|
sqrt2pi=sqrt(2.*pi) |
1712 |
|
|
sqrt2=sqrt(2.) |
1713 |
|
|
sqrtpi=sqrt(pi) |
1714 |
|
|
icefrac(:,ind2)=0. |
1715 |
|
|
althl=0. |
1716 |
|
|
althi=0. |
1717 |
|
|
athl=0. |
1718 |
|
|
athi=0. |
1719 |
|
|
senvl=0. |
1720 |
|
|
senvi=0. |
1721 |
|
|
sthi=0. |
1722 |
|
|
sthl=0. |
1723 |
|
|
sthil=0. |
1724 |
|
|
|
1725 |
|
|
|
1726 |
|
|
|
1727 |
|
|
IF (firstcall) THEN |
1728 |
|
|
|
1729 |
|
|
vert_alpha=0.5 |
1730 |
|
|
CALL getin_p('cloudth_vert_alpha',vert_alpha) |
1731 |
|
|
WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha |
1732 |
|
|
! The factor used for the thermal is equal to that of the environment |
1733 |
|
|
! if nothing is explicitly specified in the def file |
1734 |
|
|
vert_alpha_th=vert_alpha |
1735 |
|
|
CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th) |
1736 |
|
|
WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th |
1737 |
|
|
! Factor used in the calculation of sigma1s |
1738 |
|
|
CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor) |
1739 |
|
|
WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor |
1740 |
|
|
! Power used in the calculation of sigma1s |
1741 |
|
|
CALL getin_p('cloudth_sigma1s_power',sigma1s_power) |
1742 |
|
|
WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power |
1743 |
|
|
! Factor used in the calculation of sigma2s |
1744 |
|
|
CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor) |
1745 |
|
|
WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor |
1746 |
|
|
! Power used in the calculation of sigma2s |
1747 |
|
|
CALL getin_p('cloudth_sigma2s_power',sigma2s_power) |
1748 |
|
|
WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power |
1749 |
|
|
! Minimum value for the environmental air subgrid water distrib |
1750 |
|
|
CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin) |
1751 |
|
|
WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin |
1752 |
|
|
! Remove the dependency to ratqs from the variance of the vertical PDF |
1753 |
|
|
CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs) |
1754 |
|
|
WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs |
1755 |
|
|
! Modifies the PDF in thermals when ice crystals are present |
1756 |
|
|
C_mpc=1.e4 |
1757 |
|
|
CALL getin_p('C_mpc',C_mpc) |
1758 |
|
|
WRITE(*,*) 'C_mpc = ', C_mpc |
1759 |
|
|
Ni=2.0e3 |
1760 |
|
|
CALL getin_p('Ni', Ni) |
1761 |
|
|
WRITE(*,*) 'Ni = ', Ni |
1762 |
|
|
Ei=0.5 |
1763 |
|
|
CALL getin_p('Ei', Ei) |
1764 |
|
|
WRITE(*,*) 'Ei = ', Ei |
1765 |
|
|
C_cap=0.5 |
1766 |
|
|
CALL getin_p('C_cap', C_cap) |
1767 |
|
|
WRITE(*,*) 'C_cap = ', C_cap |
1768 |
|
|
d_top=1.2 |
1769 |
|
|
CALL getin_p('d_top', d_top) |
1770 |
|
|
WRITE(*,*) 'd_top = ', d_top |
1771 |
|
|
|
1772 |
|
|
firstcall=.FALSE. |
1773 |
|
|
|
1774 |
|
|
ENDIF |
1775 |
|
|
|
1776 |
|
|
|
1777 |
|
|
|
1778 |
|
|
!------------------------------------------------------------------------------- |
1779 |
|
|
! Identify grid points with potential mixed-phase conditions |
1780 |
|
|
!------------------------------------------------------------------------------- |
1781 |
|
|
|
1782 |
|
|
temp_lim=RTT-40.0 |
1783 |
|
|
|
1784 |
|
|
DO ind1=1,klon |
1785 |
|
|
IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) & |
1786 |
|
|
.AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2) & |
1787 |
|
|
.AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN |
1788 |
|
|
mpc_bl_points(ind1,ind2)=1 |
1789 |
|
|
ELSE |
1790 |
|
|
mpc_bl_points(ind1,ind2)=0 |
1791 |
|
|
ENDIF |
1792 |
|
|
ENDDO |
1793 |
|
|
|
1794 |
|
|
|
1795 |
|
|
!------------------------------------------------------------------------------- |
1796 |
|
|
! Thermal fraction calculation and standard deviation of the distribution |
1797 |
|
|
!------------------------------------------------------------------------------- |
1798 |
|
|
|
1799 |
|
|
! calculation of temperature, humidity and saturation specific humidity |
1800 |
|
|
|
1801 |
|
|
Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2) |
1802 |
|
|
Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2) |
1803 |
|
|
Tbefenvonly(:)=temp(:,ind2) |
1804 |
|
|
rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD |
1805 |
|
|
|
1806 |
|
|
zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv |
1807 |
|
|
zqth(:)=zqta(:,ind2) |
1808 |
|
|
zqenvonly(:)=po(:) |
1809 |
|
|
|
1810 |
|
|
|
1811 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv) |
1812 |
|
|
|
1813 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv) |
1814 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv) |
1815 |
|
|
|
1816 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth) |
1817 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth) |
1818 |
|
|
CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth) |
1819 |
|
|
|
1820 |
|
|
|
1821 |
|
|
DO ind1=1,klon |
1822 |
|
|
|
1823 |
|
|
|
1824 |
|
|
IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement |
1825 |
|
|
|
1826 |
|
|
|
1827 |
|
|
! Environment: |
1828 |
|
|
|
1829 |
|
|
|
1830 |
|
|
IF (Tbefenv(ind1) .GE. RTT) THEN |
1831 |
|
|
Lv=RLVTT |
1832 |
|
|
ELSE |
1833 |
|
|
Lv=RLSTT |
1834 |
|
|
ENDIF |
1835 |
|
|
|
1836 |
|
|
|
1837 |
|
|
alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 |
1838 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 |
1839 |
|
|
senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84 |
1840 |
|
|
|
1841 |
|
|
! For MPCs: |
1842 |
|
|
IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN |
1843 |
|
|
alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2) |
1844 |
|
|
aenvl=1./(1.+(alenv*Lv/cppd)) |
1845 |
|
|
senvl=aenvl*(po(ind1)-qslenv(ind1)) |
1846 |
|
|
senvi=senv |
1847 |
|
|
ENDIF |
1848 |
|
|
|
1849 |
|
|
|
1850 |
|
|
! Thermals: |
1851 |
|
|
|
1852 |
|
|
|
1853 |
|
|
IF (Tbefth(ind1) .GE. RTT) THEN |
1854 |
|
|
Lv=RLVTT |
1855 |
|
|
ELSE |
1856 |
|
|
Lv=RLSTT |
1857 |
|
|
ENDIF |
1858 |
|
|
|
1859 |
|
|
|
1860 |
|
|
alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2) |
1861 |
|
|
ath=1./(1.+(alth*Lv/cppd)) |
1862 |
|
|
sth=ath*(zqta(ind1,ind2)-zqsatth(ind1)) |
1863 |
|
|
|
1864 |
|
|
! For MPCs: |
1865 |
|
|
IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN |
1866 |
|
|
althi=alth |
1867 |
|
|
althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2) |
1868 |
|
|
athl=1./(1.+(alth*RLVTT/cppd)) |
1869 |
|
|
athi=alth |
1870 |
|
|
sthl=athl*(zqta(ind1,ind2)-qslth(ind1)) |
1871 |
|
|
sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2)) |
1872 |
|
|
sthil=athi*(zqta(ind1,ind2)-qslth(ind1)) |
1873 |
|
|
ENDIF |
1874 |
|
|
|
1875 |
|
|
|
1876 |
|
|
|
1877 |
|
|
!------------------------------------------------------------------------------- |
1878 |
|
|
! Version 3: Changes by J. Jouhaud; condensation for q > -delta s |
1879 |
|
|
! Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3 |
1880 |
|
|
!------------------------------------------------------------------------------- |
1881 |
|
|
|
1882 |
|
|
IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC |
1883 |
|
|
|
1884 |
|
|
! Standard deviation of the distributions |
1885 |
|
|
|
1886 |
|
|
sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & |
1887 |
|
|
& (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 |
1888 |
|
|
|
1889 |
|
|
IF (cloudth_ratqsmin>0.) THEN |
1890 |
|
|
sigma1s_ratqs = cloudth_ratqsmin*po(ind1) |
1891 |
|
|
ELSE |
1892 |
|
|
sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) |
1893 |
|
|
ENDIF |
1894 |
|
|
|
1895 |
|
|
sigma1s = sigma1s_fraca + sigma1s_ratqs |
1896 |
|
|
sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) |
1897 |
|
|
|
1898 |
|
|
|
1899 |
|
|
deltasenv=aenv*vert_alpha*sigma1s |
1900 |
|
|
deltasth=ath*vert_alpha_th*sigma2s |
1901 |
|
|
|
1902 |
|
|
xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) |
1903 |
|
|
xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) |
1904 |
|
|
exp_xenv1 = exp(-1.*xenv1**2) |
1905 |
|
|
exp_xenv2 = exp(-1.*xenv2**2) |
1906 |
|
|
xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) |
1907 |
|
|
xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) |
1908 |
|
|
exp_xth1 = exp(-1.*xth1**2) |
1909 |
|
|
exp_xth2 = exp(-1.*xth2**2) |
1910 |
|
|
|
1911 |
|
|
!surface CF |
1912 |
|
|
|
1913 |
|
|
cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) |
1914 |
|
|
cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) |
1915 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
1916 |
|
|
|
1917 |
|
|
|
1918 |
|
|
!volume CF and condensed water |
1919 |
|
|
|
1920 |
|
|
!environnement |
1921 |
|
|
|
1922 |
|
|
IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 |
1923 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xenv2)) |
1924 |
|
|
|
1925 |
|
|
IF (deltasenv .LT. 1.e-10) THEN |
1926 |
|
|
qcenv(ind1,ind2)=IntJ |
1927 |
|
|
cenv_vol(ind1,ind2)=IntJ_CF |
1928 |
|
|
ELSE |
1929 |
|
|
IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) |
1930 |
|
|
IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) |
1931 |
|
|
IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) |
1932 |
|
|
IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) |
1933 |
|
|
IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) |
1934 |
|
|
qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1935 |
|
|
cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
1936 |
|
|
ENDIF |
1937 |
|
|
|
1938 |
|
|
|
1939 |
|
|
|
1940 |
|
|
!thermals |
1941 |
|
|
|
1942 |
|
|
IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 |
1943 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xth2)) |
1944 |
|
|
|
1945 |
|
|
IF (deltasth .LT. 1.e-10) THEN |
1946 |
|
|
qcth(ind1,ind2)=IntJ |
1947 |
|
|
cth_vol(ind1,ind2)=IntJ_CF |
1948 |
|
|
ELSE |
1949 |
|
|
IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) |
1950 |
|
|
IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) |
1951 |
|
|
IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) |
1952 |
|
|
IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) |
1953 |
|
|
IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) |
1954 |
|
|
qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 |
1955 |
|
|
cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
1956 |
|
|
ENDIF |
1957 |
|
|
|
1958 |
|
|
qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) |
1959 |
|
|
ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) |
1960 |
|
|
|
1961 |
|
|
|
1962 |
|
|
IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN |
1963 |
|
|
ctot(ind1,ind2)=0. |
1964 |
|
|
ctot_vol(ind1,ind2)=0. |
1965 |
|
|
qcloud(ind1)=zqsatenv(ind1) |
1966 |
|
|
qincloud(ind1)=0. |
1967 |
|
|
ELSE |
1968 |
|
|
qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
1969 |
|
|
qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2) |
1970 |
|
|
ENDIF |
1971 |
|
|
|
1972 |
|
|
|
1973 |
|
|
ELSE ! mpc_bl_points>0 |
1974 |
|
|
|
1975 |
|
|
! Treat boundary layer mixed phase clouds |
1976 |
|
|
|
1977 |
|
|
! thermals |
1978 |
|
|
!========= |
1979 |
|
|
|
1980 |
|
|
! ice phase |
1981 |
|
|
!........... |
1982 |
|
|
|
1983 |
|
|
qiceth=0. |
1984 |
|
|
deltazlev_mpc=dz(ind1,:) |
1985 |
|
|
temp_mpc=ztla(ind1,:)*zpspsk(ind1,:) |
1986 |
|
|
pres_mpc=pplay(ind1,:) |
1987 |
|
|
fraca_mpc=fraca(ind1,:) |
1988 |
|
|
snowf_mpc=snowflux(ind1,:) |
1989 |
|
|
iflag_topthermals=0 |
1990 |
|
|
IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN |
1991 |
|
|
iflag_topthermals = 1 |
1992 |
|
|
ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & |
1993 |
|
|
.AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN |
1994 |
|
|
iflag_topthermals = 2 |
1995 |
|
|
ELSE |
1996 |
|
|
iflag_topthermals = 0 |
1997 |
|
|
ENDIF |
1998 |
|
|
|
1999 |
|
|
CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), & |
2000 |
|
|
qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:)) |
2001 |
|
|
|
2002 |
|
|
! qmax calculation |
2003 |
|
|
sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) |
2004 |
|
|
deltasth=athl*vert_alpha_th*sigma2s |
2005 |
|
|
xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s) |
2006 |
|
|
xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) |
2007 |
|
|
exp_xth1 = exp(-1.*xth1**2) |
2008 |
|
|
exp_xth2 = exp(-1.*xth2**2) |
2009 |
|
|
IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 |
2010 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xth2)) |
2011 |
|
|
IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) |
2012 |
|
|
IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) |
2013 |
|
|
IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) |
2014 |
|
|
IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) |
2015 |
|
|
IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) |
2016 |
|
|
qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.) |
2017 |
|
|
|
2018 |
|
|
|
2019 |
|
|
! Liquid phase |
2020 |
|
|
!................ |
2021 |
|
|
! We account for the effect of ice crystals in thermals on sthl |
2022 |
|
|
! and on the width of the distribution |
2023 |
|
|
|
2024 |
|
|
sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2)) & |
2025 |
|
|
+ (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) |
2026 |
|
|
|
2027 |
|
|
sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) & |
2028 |
|
|
/((fraca(ind1,ind2)+0.02)**sigma2s_power)) & |
2029 |
|
|
+0.002*zqta(ind1,ind2) |
2030 |
|
|
deltasthc=athl*vert_alpha_th*sigma2sc |
2031 |
|
|
|
2032 |
|
|
|
2033 |
|
|
xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc) |
2034 |
|
|
xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc) |
2035 |
|
|
exp_xth1 = exp(-1.*xth1**2) |
2036 |
|
|
exp_xth2 = exp(-1.*xth2**2) |
2037 |
|
|
IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2 |
2038 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xth2)) |
2039 |
|
|
IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1)) |
2040 |
|
|
IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) |
2041 |
|
|
IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2) |
2042 |
|
|
IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc) |
2043 |
|
|
IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc) |
2044 |
|
|
qliqth=IntJ+IntI1+IntI2+IntI3 |
2045 |
|
|
|
2046 |
|
|
qlth(ind1,ind2)=MAX(0.,qliqth) |
2047 |
|
|
|
2048 |
|
|
! Condensed water |
2049 |
|
|
|
2050 |
|
|
qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2) |
2051 |
|
|
|
2052 |
|
|
|
2053 |
|
|
! consistency with subgrid distribution |
2054 |
|
|
|
2055 |
|
|
IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN |
2056 |
|
|
fraci=qith(ind1,ind2)/qcth(ind1,ind2) |
2057 |
|
|
qcth(ind1,ind2)=qmax |
2058 |
|
|
qith(ind1,ind2)=fraci*qmax |
2059 |
|
|
qlth(ind1,ind2)=(1.-fraci)*qmax |
2060 |
|
|
ENDIF |
2061 |
|
|
|
2062 |
|
|
! Cloud Fraction |
2063 |
|
|
!............... |
2064 |
|
|
! calculation of qbase which is the value of the water vapor within mixed phase clouds |
2065 |
|
|
! such that the total water in cloud = qbase+qliqth+qiceth |
2066 |
|
|
! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction |
2067 |
|
|
! sbase and qbase calculation (note that sbase is wrt liq so negative) |
2068 |
|
|
! look for an approximate solution with iteration |
2069 |
|
|
|
2070 |
|
|
ttarget=qcth(ind1,ind2) |
2071 |
|
|
mini= athl*(qsith(ind1,ind2)-qslth(ind1)) |
2072 |
|
|
maxi= 0. !athl*(3.*sqrt(sigma2s)) |
2073 |
|
|
niter=20 |
2074 |
|
|
pas=(maxi-mini)/niter |
2075 |
|
|
stmp=mini |
2076 |
|
|
sbase=stmp |
2077 |
|
|
coutref=1.E6 |
2078 |
|
|
DO iter=1,niter |
2079 |
|
|
cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) & |
2080 |
|
|
+ stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget) |
2081 |
|
|
IF (cout .LT. coutref) THEN |
2082 |
|
|
sbase=stmp |
2083 |
|
|
coutref=cout |
2084 |
|
|
ELSE |
2085 |
|
|
stmp=stmp+pas |
2086 |
|
|
ENDIF |
2087 |
|
|
ENDDO |
2088 |
|
|
qbase=MAX(0., sbase/athl+qslth(ind1)) |
2089 |
|
|
|
2090 |
|
|
! surface cloud fraction in thermals |
2091 |
|
|
cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s)) |
2092 |
|
|
cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.) |
2093 |
|
|
|
2094 |
|
|
|
2095 |
|
|
!volume cloud fraction in thermals |
2096 |
|
|
!to be checked |
2097 |
|
|
xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s) |
2098 |
|
|
xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s) |
2099 |
|
|
exp_xth1 = exp(-1.*xth1**2) |
2100 |
|
|
exp_xth2 = exp(-1.*xth2**2) |
2101 |
|
|
|
2102 |
|
|
IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 |
2103 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xth2)) |
2104 |
|
|
|
2105 |
|
|
IF (deltasth .LT. 1.e-10) THEN |
2106 |
|
|
cth_vol(ind1,ind2)=IntJ_CF |
2107 |
|
|
ELSE |
2108 |
|
|
IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) |
2109 |
|
|
IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) |
2110 |
|
|
IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) |
2111 |
|
|
IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) |
2112 |
|
|
IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) |
2113 |
|
|
cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
2114 |
|
|
ENDIF |
2115 |
|
|
cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.) |
2116 |
|
|
|
2117 |
|
|
|
2118 |
|
|
|
2119 |
|
|
! Environment |
2120 |
|
|
!============= |
2121 |
|
|
! In the environment/downdrafts, only liquid clouds |
2122 |
|
|
! See Shupe et al. 2008, JAS |
2123 |
|
|
|
2124 |
|
|
! standard deviation of the distribution in the environment |
2125 |
|
|
sth=sthl |
2126 |
|
|
senv=senvl |
2127 |
|
|
sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & |
2128 |
|
|
& (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5 |
2129 |
|
|
! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution |
2130 |
|
|
! in the environement |
2131 |
|
|
|
2132 |
|
|
sigma1s_ratqs=1E-10 |
2133 |
|
|
IF (cloudth_ratqsmin>0.) THEN |
2134 |
|
|
sigma1s_ratqs = cloudth_ratqsmin*po(ind1) |
2135 |
|
|
ENDIF |
2136 |
|
|
|
2137 |
|
|
sigma1s = sigma1s_fraca + sigma1s_ratqs |
2138 |
|
|
deltasenv=aenvl*vert_alpha*sigma1s |
2139 |
|
|
xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s) |
2140 |
|
|
xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s) |
2141 |
|
|
exp_xenv1 = exp(-1.*xenv1**2) |
2142 |
|
|
exp_xenv2 = exp(-1.*xenv2**2) |
2143 |
|
|
|
2144 |
|
|
!surface CF |
2145 |
|
|
cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) |
2146 |
|
|
|
2147 |
|
|
!volume CF and condensed water |
2148 |
|
|
IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 |
2149 |
|
|
IntJ_CF=0.5*(1.-1.*erf(xenv2)) |
2150 |
|
|
|
2151 |
|
|
IF (deltasenv .LT. 1.e-10) THEN |
2152 |
|
|
qcenv(ind1,ind2)=IntJ |
2153 |
|
|
cenv_vol(ind1,ind2)=IntJ_CF |
2154 |
|
|
ELSE |
2155 |
|
|
IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) |
2156 |
|
|
IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) |
2157 |
|
|
IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) |
2158 |
|
|
IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) |
2159 |
|
|
IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) |
2160 |
|
|
qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment |
2161 |
|
|
cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF |
2162 |
|
|
ENDIF |
2163 |
|
|
|
2164 |
|
|
qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.) |
2165 |
|
|
cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.) |
2166 |
|
|
|
2167 |
|
|
|
2168 |
|
|
|
2169 |
|
|
! Thermals + environment |
2170 |
|
|
! ======================= |
2171 |
|
|
ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
2172 |
|
|
qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) |
2173 |
|
|
ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) |
2174 |
|
|
IF (qcth(ind1,ind2) .GT. 0) THEN |
2175 |
|
|
icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) & |
2176 |
|
|
/(fraca(ind1,ind2)*qcth(ind1,ind2) & |
2177 |
|
|
+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)) |
2178 |
|
|
icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.) |
2179 |
|
|
ELSE |
2180 |
|
|
icefrac(ind1,ind2)=0. |
2181 |
|
|
ENDIF |
2182 |
|
|
|
2183 |
|
|
IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN |
2184 |
|
|
ctot(ind1,ind2)=0. |
2185 |
|
|
ctot_vol(ind1,ind2)=0. |
2186 |
|
|
qincloud(ind1)=0. |
2187 |
|
|
qcloud(ind1)=zqsatenv(ind1) |
2188 |
|
|
ELSE |
2189 |
|
|
qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) & |
2190 |
|
|
+(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1)) |
2191 |
|
|
qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) & |
2192 |
|
|
+(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.) |
2193 |
|
|
ENDIF |
2194 |
|
|
|
2195 |
|
|
ENDIF ! mpc_bl_points |
2196 |
|
|
|
2197 |
|
|
|
2198 |
|
|
ELSE ! gaussian for environment only |
2199 |
|
|
|
2200 |
|
|
|
2201 |
|
|
|
2202 |
|
|
|
2203 |
|
|
IF (Tbefenvonly(ind1) .GE. RTT) THEN |
2204 |
|
|
Lv=RLVTT |
2205 |
|
|
ELSE |
2206 |
|
|
Lv=RLSTT |
2207 |
|
|
ENDIF |
2208 |
|
|
|
2209 |
|
|
|
2210 |
|
|
zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd) |
2211 |
|
|
alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2) |
2212 |
|
|
aenv=1./(1.+(alenv*Lv/cppd)) |
2213 |
|
|
senv=aenv*(po(ind1)-zqsatenvonly(ind1)) |
2214 |
|
|
sth=0. |
2215 |
|
|
|
2216 |
|
|
sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1) |
2217 |
|
|
sigma2s=0. |
2218 |
|
|
|
2219 |
|
|
sqrt2pi=sqrt(2.*pi) |
2220 |
|
|
xenv=senv/(sqrt(2.)*sigma1s) |
2221 |
|
|
ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
2222 |
|
|
ctot_vol(ind1,ind2)=ctot(ind1,ind2) |
2223 |
|
|
qctot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
2224 |
|
|
|
2225 |
|
|
IF (ctot(ind1,ind2).LT.1.e-3) THEN |
2226 |
|
|
ctot(ind1,ind2)=0. |
2227 |
|
|
qcloud(ind1)=zqsatenvonly(ind1) |
2228 |
|
|
qincloud(ind1)=0. |
2229 |
|
|
ELSE |
2230 |
|
|
qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1) |
2231 |
|
|
qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.) |
2232 |
|
|
ENDIF |
2233 |
|
|
|
2234 |
|
|
|
2235 |
|
|
ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492 |
2236 |
|
|
|
2237 |
|
|
! Outputs used to check the PDFs |
2238 |
|
|
cloudth_senv(ind1,ind2) = senv |
2239 |
|
|
cloudth_sth(ind1,ind2) = sth |
2240 |
|
|
cloudth_sigmaenv(ind1,ind2) = sigma1s |
2241 |
|
|
cloudth_sigmath(ind1,ind2) = sigma2s |
2242 |
|
|
|
2243 |
|
|
|
2244 |
|
|
ENDDO !loop on klon |
2245 |
|
|
|
2246 |
|
|
! Calcule ice fall velocity in thermals |
2247 |
|
|
|
2248 |
|
|
CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2)) |
2249 |
|
|
|
2250 |
|
|
RETURN |
2251 |
|
|
|
2252 |
|
|
|
2253 |
|
|
END SUBROUTINE cloudth_mpc |
2254 |
|
|
|
2255 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
2256 |
|
|
SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,vith,fraca,qith) |
2257 |
|
|
|
2258 |
|
|
! parameterization of ice for boundary |
2259 |
|
|
! layer mixed-phase clouds assuming a stationary system |
2260 |
|
|
! |
2261 |
|
|
! Note that vapor deposition on ice crystals and riming of liquid droplets |
2262 |
|
|
! depend on the ice number concentration Ni |
2263 |
|
|
! One could assume that Ni depends on qi, e.g., Ni=beta*(qi*rho)**xi |
2264 |
|
|
! and use values from Hong et al. 2004, MWR for instance |
2265 |
|
|
! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962 |
2266 |
|
|
! One could also think of a more complex expression of Ni; |
2267 |
|
|
! function of qi, T, the concentration in aerosols or INP .. |
2268 |
|
|
! Here we prefer fixing Ni to a tuning parameter |
2269 |
|
|
! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard |
2270 |
|
|
! in Mioche et al. 2017 |
2271 |
|
|
! |
2272 |
|
|
! |
2273 |
|
|
! References: |
2274 |
|
|
!------------ |
2275 |
|
|
! This parameterization is thoroughly described in Vignon et al. |
2276 |
|
|
! |
2277 |
|
|
! More specifically |
2278 |
|
|
! for the Water vapor deposition process: |
2279 |
|
|
! |
2280 |
|
|
! Rotstayn, L. D., 1997: A physically based scheme for the treat- |
2281 |
|
|
! ment of stratiform cloudfs and precipitation in large-scale |
2282 |
|
|
! models. I: Description and evaluation of the microphysical |
2283 |
|
|
! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282. |
2284 |
|
|
! |
2285 |
|
|
! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk |
2286 |
|
|
! stratiform cloud microphysics scheme in the NCAR Com- |
2287 |
|
|
! munity Atmosphere Model (CAM3). Part I: Description and |
2288 |
|
|
! numerical tests. J. Climate, 21, 3642–3659 |
2289 |
|
|
! |
2290 |
|
|
! for the Riming process: |
2291 |
|
|
! |
2292 |
|
|
! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro- |
2293 |
|
|
! scale structure and organization of clouds and precipitation in |
2294 |
|
|
! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’ |
2295 |
|
|
! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206 |
2296 |
|
|
! |
2297 |
|
|
! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit |
2298 |
|
|
! forecasts of winter precipitation using an improved bulk |
2299 |
|
|
! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit |
2300 |
|
|
! forecasts of winter precipitation using an improved bulk |
2301 |
|
|
! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542 |
2302 |
|
|
! |
2303 |
|
|
! For the formation of clouds by thermals: |
2304 |
|
|
! |
2305 |
|
|
! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of |
2306 |
|
|
! the Atmospheric Sciences, 65, 407–425. |
2307 |
|
|
! |
2308 |
|
|
! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a |
2309 |
|
|
! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3 |
2310 |
|
|
! |
2311 |
|
|
! |
2312 |
|
|
! |
2313 |
|
|
! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr |
2314 |
|
|
!============================================================================= |
2315 |
|
|
|
2316 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
2317 |
|
|
USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm |
2318 |
|
|
|
2319 |
|
|
IMPLICIT none |
2320 |
|
|
|
2321 |
|
|
INCLUDE "YOMCST.h" |
2322 |
|
|
|
2323 |
|
|
INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions |
2324 |
|
|
INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ? |
2325 |
|
|
REAL, INTENT(IN) :: Ni ! ice number concentration [m-3] |
2326 |
|
|
REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency |
2327 |
|
|
REAL, INTENT(IN) :: C_cap ! ice crystal capacitance |
2328 |
|
|
REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter |
2329 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals |
2330 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa] |
2331 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg] |
2332 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg] |
2333 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg] |
2334 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m] |
2335 |
|
|
REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s] |
2336 |
|
|
REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals |
2337 |
|
|
REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg] |
2338 |
|
|
|
2339 |
|
|
|
2340 |
|
|
INTEGER ind2p1,ind2p2 |
2341 |
|
|
REAL rho(klev) |
2342 |
|
|
REAL unsurtaudet, unsurtaustardep, unsurtaurim |
2343 |
|
|
REAL qi, AA, BB, Ka, Dv, rhoi |
2344 |
|
|
REAL p0, t0, fp1, fp2 |
2345 |
|
|
REAL alpha, flux_term |
2346 |
|
|
REAL det_term, precip_term, rim_term, dep_term |
2347 |
|
|
|
2348 |
|
|
|
2349 |
|
|
ind2p1=ind2+1 |
2350 |
|
|
ind2p2=ind2+2 |
2351 |
|
|
|
2352 |
|
|
rho=pres/temp/RD ! air density kg/m3 |
2353 |
|
|
|
2354 |
|
|
Ka=2.4e-2 ! thermal conductivity of the air, SI |
2355 |
|
|
p0=101325.0 ! ref pressure |
2356 |
|
|
T0=273.15 ! ref temp |
2357 |
|
|
rhoi=500.0 ! cloud ice density following Reisner et al. 1998 |
2358 |
|
|
alpha=700. ! fallvelocity param |
2359 |
|
|
|
2360 |
|
|
|
2361 |
|
|
IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels |
2362 |
|
|
|
2363 |
|
|
Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI |
2364 |
|
|
|
2365 |
|
|
! Detrainment term: |
2366 |
|
|
|
2367 |
|
|
unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2) |
2368 |
|
|
|
2369 |
|
|
|
2370 |
|
|
! Deposition term |
2371 |
|
|
AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.) |
2372 |
|
|
BB=1./(rho(ind2)*Dv*qsith(ind2)) |
2373 |
|
|
unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) & |
2374 |
|
|
*4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33) |
2375 |
|
|
|
2376 |
|
|
! Riming term neglected at this level |
2377 |
|
|
!unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4) |
2378 |
|
|
|
2379 |
|
|
qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12) |
2380 |
|
|
qi=MAX(qi,0.)**(3./2.) |
2381 |
|
|
|
2382 |
|
|
ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2 |
2383 |
|
|
|
2384 |
|
|
Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI |
2385 |
|
|
|
2386 |
|
|
! Detrainment term: |
2387 |
|
|
|
2388 |
|
|
unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1) |
2389 |
|
|
det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1) |
2390 |
|
|
|
2391 |
|
|
|
2392 |
|
|
! Deposition term |
2393 |
|
|
|
2394 |
|
|
AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.) |
2395 |
|
|
BB=1./(rho(ind2p1)*Dv*qsith(ind2p1)) |
2396 |
|
|
unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1)) & |
2397 |
|
|
/qsith(ind2p1)*4.*RPI/(AA+BB) & |
2398 |
|
|
*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33) |
2399 |
|
|
dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep |
2400 |
|
|
|
2401 |
|
|
! Riming term |
2402 |
|
|
|
2403 |
|
|
unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4) |
2404 |
|
|
rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim |
2405 |
|
|
|
2406 |
|
|
! Precip term |
2407 |
|
|
|
2408 |
|
|
! We assume that there is no solid precipitation outside thermals |
2409 |
|
|
! so the precipitation flux within thermals is equal to the precipitation flux |
2410 |
|
|
! at mesh-scale divided by thermals fraction |
2411 |
|
|
|
2412 |
|
|
fp2=0. |
2413 |
|
|
fp1=0. |
2414 |
|
|
IF (fraca(ind2p1) .GT. 0.) THEN |
2415 |
|
|
fp2=-qith(ind2p2)*rho(ind2p2)*vith(ind2p2)*fraca(ind2p2)! flux defined positive upward |
2416 |
|
|
fp1=-qith(ind2p1)*rho(ind2p1)*vith(ind2p1)*fraca(ind2p1) |
2417 |
|
|
ENDIF |
2418 |
|
|
|
2419 |
|
|
precip_term=-1./deltazlev(ind2p1)*(fp2-fp1) |
2420 |
|
|
|
2421 |
|
|
! Calculation in a top-to-bottom loop |
2422 |
|
|
|
2423 |
|
|
IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN |
2424 |
|
|
qi= 1./fm_therm(ind1,ind2p1)* & |
2425 |
|
|
(deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + & |
2426 |
|
|
fm_therm(ind1,ind2p2)*(qith(ind2p1))) |
2427 |
|
|
END IF |
2428 |
|
|
|
2429 |
|
|
ENDIF ! top thermals |
2430 |
|
|
|
2431 |
|
|
qith(ind2)=MAX(0.,qi) |
2432 |
|
|
|
2433 |
|
|
RETURN |
2434 |
|
|
|
2435 |
|
|
END SUBROUTINE ICE_MPC_BL_CLOUDS |
2436 |
|
|
|
2437 |
|
|
|
2438 |
|
|
|
2439 |
|
|
|
2440 |
|
|
END MODULE cloudth_mod |
2441 |
|
|
|
2442 |
|
|
|
2443 |
|
|
|
2444 |
|
|
|