module_mp_wsm5.F
References to this file elsewhere.
1 #if ( RWORDSIZE == 4 )
2 # define VREC vsrec
3 # define VSQRT vssqrt
4 #else
5 # define VREC vrec
6 # define VSQRT vsqrt
7 #endif
8
9 !Including inline expansion statistical function
10 MODULE module_mp_wsm5
11 !
12 !
13 REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
14 REAL, PARAMETER, PRIVATE :: n0r = 8.e6
15 REAL, PARAMETER, PRIVATE :: avtr = 841.9
16 REAL, PARAMETER, PRIVATE :: bvtr = 0.8
17 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
18 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
19 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
20 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
21 REAL, PARAMETER, PRIVATE :: avts = 11.72
22 REAL, PARAMETER, PRIVATE :: bvts = .41
23 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
24 REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
25 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
26 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
27 REAL, PARAMETER, PRIVATE :: betai = .6
28 REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
29 REAL, PARAMETER, PRIVATE :: dicon = 11.9
30 REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
31 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
32 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
33 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
34 REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
35 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
36 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
37 REAL, PARAMETER, PRIVATE :: t40c = 233.16
38 REAL, PARAMETER, PRIVATE :: eacrc = 1.0
39 REAL, SAVE :: &
40 qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
41 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
42 precr1,precr2,xm0,xmmax,roqimax,bvts1, &
43 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
44 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
45 pidn0s,xlv1,pacrc, &
46 rslopermax,rslopesmax,rslopegmax, &
47 rsloperbmax,rslopesbmax,rslopegbmax, &
48 rsloper2max,rslopes2max,rslopeg2max, &
49 rsloper3max,rslopes3max,rslopeg3max
50 !
51 ! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
52 !
53 CONTAINS
54 !===================================================================
55 !
56 SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
57 ,den, pii, p, delz &
58 ,delt,g, cpd, cpv, rd, rv, t0c &
59 ,ep1, ep2, qmin &
60 ,XLS, XLV0, XLF0, den0, denr &
61 ,cliq,cice,psat &
62 ,rain, rainncv &
63 ,snow, snowncv &
64 ,sr &
65 ,ids,ide, jds,jde, kds,kde &
66 ,ims,ime, jms,jme, kms,kme &
67 ,its,ite, jts,jte, kts,kte &
68 )
69 !-------------------------------------------------------------------
70 IMPLICIT NONE
71 !-------------------------------------------------------------------
72 !
73 ! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
74 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
75 ! number concentration is a function of temperature, and seperate assumption
76 ! is developed, in which ice crystal number concentration is a function
77 ! of ice amount. A theoretical background of the ice-microphysics and related
78 ! processes in the WSMMPs are described in Hong et al. (2004).
79 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
80 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
81 !
82 ! WSM5 cloud scheme
83 !
84 ! Coded by Song-You Hong (Yonsei Univ.)
85 ! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
86 ! Summer 2002
87 !
88 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
89 ! Summer 2003
90 !
91 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
92 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
93 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
94 !
95 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
96 ims,ime, jms,jme, kms,kme , &
97 its,ite, jts,jte, kts,kte
98 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
99 INTENT(INOUT) :: &
100 th, &
101 q, &
102 qc, &
103 qi, &
104 qr, &
105 qs
106 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
107 INTENT(IN ) :: &
108 den, &
109 pii, &
110 p, &
111 delz
112 REAL, INTENT(IN ) :: delt, &
113 g, &
114 rd, &
115 rv, &
116 t0c, &
117 den0, &
118 cpd, &
119 cpv, &
120 ep1, &
121 ep2, &
122 qmin, &
123 XLS, &
124 XLV0, &
125 XLF0, &
126 cliq, &
127 cice, &
128 psat, &
129 denr
130 REAL, DIMENSION( ims:ime , jms:jme ), &
131 INTENT(INOUT) :: rain, &
132 rainncv, &
133 sr
134
135 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
136 INTENT(INOUT) :: snow, &
137 snowncv
138
139 ! LOCAL VAR
140 REAL, DIMENSION( its:ite , kts:kte ) :: t
141 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
142 INTEGER :: i,j,k
143 !-------------------------------------------------------------------
144 DO j=jts,jte
145 DO k=kts,kte
146 DO i=its,ite
147 t(i,k)=th(i,k,j)*pii(i,k,j)
148 qci(i,k,1) = qc(i,k,j)
149 qci(i,k,2) = qi(i,k,j)
150 qrs(i,k,1) = qr(i,k,j)
151 qrs(i,k,2) = qs(i,k,j)
152 ENDDO
153 ENDDO
154 CALL wsm52D(t, q(ims,kms,j), qci, qrs &
155 ,den(ims,kms,j) &
156 ,p(ims,kms,j), delz(ims,kms,j) &
157 ,delt,g, cpd, cpv, rd, rv, t0c &
158 ,ep1, ep2, qmin &
159 ,XLS, XLV0, XLF0, den0, denr &
160 ,cliq,cice,psat &
161 ,j &
162 ,rain(ims,j),rainncv(ims,j) &
163 ,sr(ims,j) &
164 ,ids,ide, jds,jde, kds,kde &
165 ,ims,ime, jms,jme, kms,kme &
166 ,its,ite, jts,jte, kts,kte &
167 ,snow(ims,j),snowncv(ims,j) &
168 )
169 DO K=kts,kte
170 DO I=its,ite
171 th(i,k,j)=t(i,k)/pii(i,k,j)
172 qc(i,k,j) = qci(i,k,1)
173 qi(i,k,j) = qci(i,k,2)
174 qr(i,k,j) = qrs(i,k,1)
175 qs(i,k,j) = qrs(i,k,2)
176 ENDDO
177 ENDDO
178 ENDDO
179 END SUBROUTINE wsm5
180 !===================================================================
181 !
182 SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz &
183 ,delt,g, cpd, cpv, rd, rv, t0c &
184 ,ep1, ep2, qmin &
185 ,XLS, XLV0, XLF0, den0, denr &
186 ,cliq,cice,psat &
187 ,lat &
188 ,rain,rainncv &
189 ,sr &
190 ,ids,ide, jds,jde, kds,kde &
191 ,ims,ime, jms,jme, kms,kme &
192 ,its,ite, jts,jte, kts,kte &
193 ,snow,snowncv &
194 )
195 !-------------------------------------------------------------------
196 IMPLICIT NONE
197 !-------------------------------------------------------------------
198 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
199 ims,ime, jms,jme, kms,kme , &
200 its,ite, jts,jte, kts,kte, &
201 lat
202 REAL, DIMENSION( its:ite , kts:kte ), &
203 INTENT(INOUT) :: &
204 t
205 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
206 INTENT(INOUT) :: &
207 qci, &
208 qrs
209 REAL, DIMENSION( ims:ime , kms:kme ), &
210 INTENT(INOUT) :: &
211 q
212 REAL, DIMENSION( ims:ime , kms:kme ), &
213 INTENT(IN ) :: &
214 den, &
215 p, &
216 delz
217 REAL, INTENT(IN ) :: delt, &
218 g, &
219 cpd, &
220 cpv, &
221 t0c, &
222 den0, &
223 rd, &
224 rv, &
225 ep1, &
226 ep2, &
227 qmin, &
228 XLS, &
229 XLV0, &
230 XLF0, &
231 cliq, &
232 cice, &
233 psat, &
234 denr
235 REAL, DIMENSION( ims:ime ), &
236 INTENT(INOUT) :: rain, &
237 rainncv, &
238 sr
239
240 REAL, DIMENSION( ims:ime ), OPTIONAL, &
241 INTENT(INOUT) :: snow, &
242 snowncv
243
244 ! LOCAL VAR
245 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
246 rh, qs, rslope, rslope2, rslope3, rslopeb, &
247 falk, fall, work1
248 REAL, DIMENSION( its:ite , kts:kte ) :: &
249 falkc, work1c, work2c, fallc
250 REAL, DIMENSION( its:ite , kts:kte ) :: &
251 praut, psaut, prevp, psdep, pracw, psaci, psacw, &
252 pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, xni,&
253 n0sfac
254 ! variables for optimization
255 REAL, DIMENSION( its:ite ) :: tvec1
256 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
257 REAL, DIMENSION(its:ite) :: rmstep
258 REAL dtcldden, rdelz, rdtcld
259 LOGICAL, DIMENSION( its:ite ) :: flgcld
260 REAL :: pi, &
261 cpmcal, xlcal, lamdar, lamdas, diffus, &
262 viscos, xka, venfac, conden, diffac, &
263 x, y, z, a, b, c, d, e, &
264 qdt, holdrr, holdrs, supcol, pvt, &
265 coeres, supsat, dtcld, xmi, eacrs, satdt, &
266 vt2i,vt2s,acrfac, &
267 qimax, diameter, xni0, roqi0, &
268 fallsum, fallsum_qsi, xlwork2, factor, source, &
269 value, xlf, pfrzdtc, pfrzdtr, supice
270 REAL :: temp
271 REAL :: holdc, holdci
272 INTEGER :: i, j, k, mstepmax, &
273 iprt, latd, lond, loop, loops, ifsat, n
274 ! Temporaries used for inlining fpvs function
275 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
276 !
277 !=================================================================
278 ! compute internal functions
279 !
280 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
281 xlcal(x) = xlv0-xlv1*(x-t0c)
282 !----------------------------------------------------------------
283 ! size distributions: (x=mixing ratio, y=air density):
284 ! valid for mixing ratio > 1.e-9 kg/kg.
285 !
286 ! Optimizatin : A**B => exp(log(A)*(B))
287 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
288 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
289 !
290 !----------------------------------------------------------------
291 ! diffus: diffusion coefficient of the water vapor
292 ! viscos: kinematic viscosity(m2s-1)
293 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
294 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
295 ! xka(x,y) = 1.414e3*viscos(x,y)*y
296 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
297 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
298 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
299 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
300 !
301 !
302 pi = 4. * atan(1.)
303 !
304 !----------------------------------------------------------------
305 ! paddint 0 for negative values generated by dynamics
306 !
307 do k = kts, kte
308 do i = its, ite
309 qci(i,k,1) = max(qci(i,k,1),0.0)
310 qrs(i,k,1) = max(qrs(i,k,1),0.0)
311 qci(i,k,2) = max(qci(i,k,2),0.0)
312 qrs(i,k,2) = max(qrs(i,k,2),0.0)
313 enddo
314 enddo
315 !
316 !----------------------------------------------------------------
317 ! latent heat for phase changes and heat capacity. neglect the
318 ! changes during microphysical process calculation
319 ! emanuel(1994)
320 !
321 do k = kts, kte
322 do i = its, ite
323 cpm(i,k) = cpmcal(q(i,k))
324 xl(i,k) = xlcal(t(i,k))
325 enddo
326 enddo
327 !
328 !----------------------------------------------------------------
329 ! compute the minor time steps.
330 !
331 loops = max(nint(delt/dtcldcr),1)
332 dtcld = delt/loops
333 if(delt.le.dtcldcr) dtcld = delt
334 !
335 do loop = 1,loops
336 !
337 !----------------------------------------------------------------
338 ! initialize the large scale variables
339 !
340 do i = its, ite
341 mstep(i) = 1
342 flgcld(i) = .true.
343 enddo
344 !
345 ! do k = kts, kte
346 ! do i = its, ite
347 ! denfac(i,k) = sqrt(den0/den(i,k))
348 ! enddo
349 ! enddo
350 do k = kts, kte
351 CALL VREC( tvec1(its), den(its,k), ite-its+1)
352 do i = its, ite
353 tvec1(i) = tvec1(i)*den0
354 enddo
355 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
356 enddo
357 !
358 ! Inline expansion for fpvs
359 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
360 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
361 hsub = xls
362 hvap = xlv0
363 cvap = cpv
364 ttp=t0c+0.01
365 dldt=cvap-cliq
366 xa=-dldt/rv
367 xb=xa+hvap/(rv*ttp)
368 dldti=cvap-cice
369 xai=-dldti/rv
370 xbi=xai+hsub/(rv*ttp)
371 do k = kts, kte
372 do i = its, ite
373 tr=ttp/t(i,k)
374 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
375 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
376 qs(i,k,1) = max(qs(i,k,1),qmin)
377 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
378 tr=ttp/t(i,k)
379 if(t(i,k).lt.ttp) then
380 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
381 else
382 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
383 endif
384 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
385 qs(i,k,2) = max(qs(i,k,2),qmin)
386 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
387 enddo
388 enddo
389 !
390 !----------------------------------------------------------------
391 ! initialize the variables for microphysical physics
392 !
393 !
394 do k = kts, kte
395 do i = its, ite
396 prevp(i,k) = 0.
397 psdep(i,k) = 0.
398 praut(i,k) = 0.
399 psaut(i,k) = 0.
400 pracw(i,k) = 0.
401 psaci(i,k) = 0.
402 psacw(i,k) = 0.
403 pigen(i,k) = 0.
404 pidep(i,k) = 0.
405 pcond(i,k) = 0.
406 psmlt(i,k) = 0.
407 psevp(i,k) = 0.
408 falk(i,k,1) = 0.
409 falk(i,k,2) = 0.
410 fall(i,k,1) = 0.
411 fall(i,k,2) = 0.
412 fallc(i,k) = 0.
413 falkc(i,k) = 0.
414 xni(i,k) = 1.e3
415 enddo
416 enddo
417 !
418 !----------------------------------------------------------------
419 ! compute the fallout term:
420 ! first, vertical terminal velosity for minor loops
421 !
422 do k = kts, kte
423 do i = its, ite
424 supcol = t0c-t(i,k)
425 !---------------------------------------------------------------
426 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
427 !---------------------------------------------------------------
428 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
429 if(qrs(i,k,1).le.qcrmin)then
430 rslope(i,k,1) = rslopermax
431 rslopeb(i,k,1) = rsloperbmax
432 rslope2(i,k,1) = rsloper2max
433 rslope3(i,k,1) = rsloper3max
434 else
435 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
436 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
437 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
438 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
439 endif
440 if(qrs(i,k,2).le.qcrmin)then
441 rslope(i,k,2) = rslopesmax
442 rslopeb(i,k,2) = rslopesbmax
443 rslope2(i,k,2) = rslopes2max
444 rslope3(i,k,2) = rslopes3max
445 else
446 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
447 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
448 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
449 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
450 endif
451 !-------------------------------------------------------------
452 ! Ni: ice crystal number concentraiton [HDC 5c]
453 !-------------------------------------------------------------
454 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
455 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
456 temp = (den(i,k)*max(qci(i,k,2),qmin))
457 temp = sqrt(sqrt(temp*temp*temp))
458 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
459 enddo
460 enddo
461 !
462 mstepmax = 1
463 numdt = 1
464 do k = kte, kts, -1
465 do i = its, ite
466 work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
467 work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
468 numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
469 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
470 enddo
471 enddo
472 do i = its, ite
473 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
474 rmstep(i) = 1./mstep(i)
475 enddo
476 !
477 do n = 1, mstepmax
478 k = kte
479 do i = its, ite
480 if(n.le.mstep(i)) then
481 ! falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
482 ! falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
483 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
484 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
485 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
486 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
487 ! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
488 ! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
489 dtcldden = dtcld/den(i,k)
490 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
491 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
492 endif
493 enddo
494 do k = kte-1, kts, -1
495 do i = its, ite
496 if(n.le.mstep(i)) then
497 ! falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
498 ! falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
499 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
500 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
501 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
502 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
503 ! qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
504 ! *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
505 ! qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
506 ! *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
507 dtcldden = dtcld/den(i,k)
508 rdelz = 1./delz(i,k)
509 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
510 *delz(i,k+1)*rdelz)*dtcldden,0.)
511 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
512 *delz(i,k+1)*rdelz)*dtcldden,0.)
513 endif
514 enddo
515 enddo
516 do k = kte, kts, -1
517 do i = its, ite
518 if(n.le.mstep(i)) then
519 if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
520 !----------------------------------------------------------------
521 ! psmlt: melting of snow [HL A33] [RH83 A25]
522 ! (T>T0: S->R)
523 !----------------------------------------------------------------
524 xlf = xlf0
525 ! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
526 work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
527 /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
528 *exp(log(t(i,k))*(1.81))/p(i,k)))) &
529 *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
530 *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
531 *sqrt(sqrt(den0/(den(i,k)))))
532 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
533 ! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
534 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
535 ! *work2(i,k)*coeres)
536 psmlt(i,k) = &
537 (1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k))) /((t(i,k))+120.)/(den(i,k)) )*(den(i,k)))&
538 /xlf*(t0c-t(i,k))*pi/2. &
539 *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
540 *work2(i,k)*coeres)
541 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
542 -qrs(i,k,2)/mstep(i)),0.)
543 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
544 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
545 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
546 endif
547 endif
548 enddo
549 enddo
550 enddo
551 !---------------------------------------------------------------
552 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
553 !---------------------------------------------------------------
554 mstepmax = 1
555 mstep = 1
556 numdt = 1
557 do k = kte, kts, -1
558 do i = its, ite
559 if(qci(i,k,2).le.0.) then
560 work2c(i,k) = 0.
561 else
562 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
563 ! diameter = min(dicon * sqrt(xmi),dimax)
564 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
565 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
566 work2c(i,k) = work1c(i,k)/delz(i,k)
567 endif
568 numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
569 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
570 enddo
571 enddo
572 do i = its, ite
573 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
574 enddo
575 !
576 do n = 1, mstepmax
577 k = kte
578 do i = its, ite
579 if(n.le.mstep(i)) then
580 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
581 holdc = falkc(i,k)
582 fallc(i,k) = fallc(i,k)+falkc(i,k)
583 holdci = qci(i,k,2)
584 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
585 endif
586 enddo
587 do k = kte-1, kts, -1
588 do i = its, ite
589 if(n.le.mstep(i)) then
590 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
591 holdc = falkc(i,k)
592 fallc(i,k) = fallc(i,k)+falkc(i,k)
593 holdci = qci(i,k,2)
594 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
595 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
596 endif
597 enddo
598 enddo
599 enddo
600 !
601 !
602 !----------------------------------------------------------------
603 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
604 !
605 do i = its, ite
606 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
607 fallsum_qsi = fall(i,1,2)+fallc(i,1)
608 rainncv(i) = 0.
609 if(fallsum.gt.0.) then
610 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
611 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
612 endif
613 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
614 snowncv(i) = 0.
615 if(fallsum_qsi.gt.0.) then
616 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
617 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
618 endif
619 ENDIF
620 sr(i) = 0.
621 if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000./(rainncv(i)+1.e-12)
622 enddo
623 !
624 !---------------------------------------------------------------
625 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
626 ! (T>T0: I->C)
627 !---------------------------------------------------------------
628 do k = kts, kte
629 do i = its, ite
630 supcol = t0c-t(i,k)
631 xlf = xls-xl(i,k)
632 if(supcol.lt.0.) xlf = xlf0
633 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
634 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
635 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
636 qci(i,k,2) = 0.
637 endif
638 !---------------------------------------------------------------
639 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
640 ! (T<-40C: C->I)
641 !---------------------------------------------------------------
642 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
643 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
644 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
645 qci(i,k,1) = 0.
646 endif
647 !---------------------------------------------------------------
648 ! pihtf: heterogeneous freezing of cloud water [HL A44]
649 ! (T0>T>-40C: C->I)
650 !---------------------------------------------------------------
651 if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
652 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
653 ! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
654 pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
655 *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
656 qci(i,k,2) = qci(i,k,2) + pfrzdtc
657 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
658 qci(i,k,1) = qci(i,k,1)-pfrzdtc
659 endif
660 !---------------------------------------------------------------
661 ! psfrz: freezing of rain water [HL A20] [LFO 45]
662 ! (T<T0, R->S)
663 !---------------------------------------------------------------
664 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
665 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
666 ! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
667 ! qrs(i,k,1))
668 temp = rslope(i,k,1)
669 temp = temp*temp*temp*temp*temp*temp*temp
670 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
671 *(exp(pfrz2*supcol)-1.)*temp*dtcld, &
672 qrs(i,k,1))
673 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
674 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
675 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
676 endif
677 enddo
678 enddo
679 !
680 !----------------------------------------------------------------
681 ! rsloper: reverse of the slope parameter of the rain(m)
682 ! xka: thermal conductivity of air(jm-1s-1k-1)
683 ! work1: the thermodynamic term in the denominator associated with
684 ! heat conduction and vapor diffusion
685 ! (ry88, y93, h85)
686 ! work2: parameter associated with the ventilation effects(y93)
687 !
688 do k = kts, kte
689 do i = its, ite
690 if(qrs(i,k,1).le.qcrmin)then
691 rslope(i,k,1) = rslopermax
692 rslopeb(i,k,1) = rsloperbmax
693 rslope2(i,k,1) = rsloper2max
694 rslope3(i,k,1) = rsloper3max
695 else
696 ! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
697 rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
698 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
699 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
700 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
701 endif
702 if(qrs(i,k,2).le.qcrmin)then
703 rslope(i,k,2) = rslopesmax
704 rslopeb(i,k,2) = rslopesbmax
705 rslope2(i,k,2) = rslopes2max
706 rslope3(i,k,2) = rslopes3max
707 else
708 ! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
709 rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2))*(den(i,k))))))
710 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
711 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
712 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
713 endif
714 enddo
715 enddo
716 !
717 do k = kts, kte
718 do i = its, ite
719 ! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
720 work1(i,k,1) = &
721 ((((den(i,k))*(xl(i,k))*(xl(i,k))) * ((t(i,k))+120.) * (den(i,k))) &
722 / &
723 ( 1.414e3 * (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) * (den(i,k)) * &
724 (rv*(t(i,k))*(t(i,k))))) &
725 + &
726 p(i,k) / ( (qs(i,k,1)) * ( 8.794e-5 * exp(log(t(i,k))*(1.81)) ) )
727 ! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
728 work1(i,k,2) = &
729 ( &
730 (((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k))) &
731 / &
732 ( &
733 1.414e3 * (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) * (den(i,k)) * &
734 (rv*(t(i,k))*(t(i,k))) &
735 ) &
736 + &
737 p(i,k) &
738 / &
739 ( qs(i,k,2) * (8.794e-5 * exp(log(t(i,k))*(1.81)))) &
740 )
741 ! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
742 work2(i,k) = &
743 ( &
744 exp(.3333333*log( &
745 ((1.496e-6 * ((t(i,k))*sqrt(t(i,k))))*p(i,k)) &
746 / &
747 (((t(i,k))+120.)*den(i,k)*(8.794e-5 * exp(log(t(i,k))*(1.81)))) &
748 )) &
749 * &
750 sqrt(sqrt(den0/(den(i,k)))) &
751 ) &
752 / &
753 sqrt( &
754 (1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
755 / &
756 ( &
757 ((t(i,k))+120.) * den(i,k) &
758 ) &
759 )
760 ENDDO
761 ENDDO
762 !
763 !===============================================================
764 !
765 ! warm rain processes
766 !
767 ! - follows the processes in RH83 and LFO except for autoconcersion
768 !
769 !===============================================================
770 !
771 do k = kts, kte
772 do i = its, ite
773 supsat = max(q(i,k),qmin)-qs(i,k,1)
774 satdt = supsat/dtcld
775 !---------------------------------------------------------------
776 ! praut: auto conversion rate from cloud to rain [HDC 16]
777 ! (C->R)
778 !---------------------------------------------------------------
779 if(qci(i,k,1).gt.qc0) then
780 praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
781 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
782 endif
783 !---------------------------------------------------------------
784 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
785 ! (C->R)
786 !---------------------------------------------------------------
787 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
788 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
789 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
790 endif
791 !---------------------------------------------------------------
792 ! prevp: evaporation/condensation rate of rain [HDC 14]
793 ! (V->R or R->V)
794 !---------------------------------------------------------------
795 if(qrs(i,k,1).gt.0.) then
796 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
797 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
798 +precr2*work2(i,k)*coeres)/work1(i,k,1)
799 if(prevp(i,k).lt.0.) then
800 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
801 prevp(i,k) = max(prevp(i,k),satdt/2)
802 else
803 prevp(i,k) = min(prevp(i,k),satdt/2)
804 endif
805 endif
806 enddo
807 enddo
808 !
809 !===============================================================
810 !
811 ! cold rain processes
812 !
813 ! - follows the revised ice microphysics processes in HDC
814 ! - the processes same as in RH83 and RH84 and LFO behave
815 ! following ice crystal hapits defined in HDC, inclduing
816 ! intercept parameter for snow (n0s), ice crystal number
817 ! concentration (ni), ice nuclei number concentration
818 ! (n0i), ice diameter (d)
819 !
820 !===============================================================
821 !
822 rdtcld = 1./dtcld
823 do k = kts, kte
824 do i = its, ite
825 supcol = t0c-t(i,k)
826 supsat = max(q(i,k),qmin)-qs(i,k,2)
827 satdt = supsat/dtcld
828 ifsat = 0
829 !-------------------------------------------------------------
830 ! Ni: ice crystal number concentraiton [HDC 5c]
831 !-------------------------------------------------------------
832 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
833 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
834 temp = (den(i,k)*max(qci(i,k,2),qmin))
835 temp = sqrt(sqrt(temp*temp*temp))
836 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
837 eacrs = exp(0.07*(-supcol))
838 !
839 if(supcol.gt.0) then
840 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
841 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
842 diameter = min(dicon * sqrt(xmi),dimax)
843 vt2i = 1.49e4*diameter**1.31
844 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
845 !-------------------------------------------------------------
846 ! psaci: Accretion of cloud ice by rain [HDC 10]
847 ! (T<T0: I->S)
848 !-------------------------------------------------------------
849 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
850 +diameter**2*rslope(i,k,2)
851 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
852 *abs(vt2s-vt2i)*acrfac/4.
853 endif
854 !-------------------------------------------------------------
855 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
856 ! (T<T0: C->S, and T>=T0: C->R)
857 !-------------------------------------------------------------
858 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
859 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
860 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
861 ! ,qci(i,k,1)/dtcld)
862 ,qci(i,k,1)*rdtcld)
863 endif
864 !-------------------------------------------------------------
865 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
866 ! (T<T0: V->I or I->V)
867 !-------------------------------------------------------------
868 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
869 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
870 diameter = dicon * sqrt(xmi)
871 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
872 supice = satdt-prevp(i,k)
873 if(pidep(i,k).lt.0.) then
874 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
875 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
876 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
877 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
878 else
879 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
880 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
881 endif
882 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
883 endif
884 endif
885 !-------------------------------------------------------------
886 ! psdep: deposition/sublimation rate of snow [HDC 14]
887 ! (V->S or S->V)
888 !-------------------------------------------------------------
889 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
890 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
891 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
892 *(precs1*rslope2(i,k,2)+precs2 &
893 *work2(i,k)*coeres)/work1(i,k,2)
894 supice = satdt-prevp(i,k)-pidep(i,k)
895 if(psdep(i,k).lt.0.) then
896 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
897 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
898 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
899 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
900 else
901 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
902 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
903 endif
904 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
905 ifsat = 1
906 endif
907 !-------------------------------------------------------------
908 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
909 ! (T<T0: V->I)
910 !-------------------------------------------------------------
911 if(supcol.gt.0) then
912 if(supsat.gt.0.and.ifsat.ne.1) then
913 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
914 xni0 = 1.e3*exp(0.1*supcol)
915 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
916 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
917 ! /dtcld)
918 *rdtcld)
919 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
920 endif
921 !
922 !-------------------------------------------------------------
923 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
924 ! (T<T0: I->S)
925 !-------------------------------------------------------------
926 if(qci(i,k,2).gt.0.) then
927 qimax = roqimax/den(i,k)
928 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
929 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
930 endif
931 endif
932 !-------------------------------------------------------------
933 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
934 ! (T>T0: S->V)
935 !-------------------------------------------------------------
936 if(supcol.lt.0.) then
937 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
938 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
939 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
940 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
941 endif
942 enddo
943 enddo
944 !
945 !
946 !----------------------------------------------------------------
947 ! check mass conservation of generation terms and feedback to the
948 ! large scale
949 !
950 do k = kts, kte
951 do i = its, ite
952 if(t(i,k).le.t0c) then
953 !
954 ! cloud water
955 !
956 value = max(qmin,qci(i,k,1))
957 source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
958 if (source.gt.value) then
959 factor = value/source
960 praut(i,k) = praut(i,k)*factor
961 pracw(i,k) = pracw(i,k)*factor
962 psacw(i,k) = psacw(i,k)*factor
963 endif
964 !
965 ! cloud ice
966 !
967 value = max(qmin,qci(i,k,2))
968 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
969 if (source.gt.value) then
970 factor = value/source
971 psaut(i,k) = psaut(i,k)*factor
972 psaci(i,k) = psaci(i,k)*factor
973 pigen(i,k) = pigen(i,k)*factor
974 pidep(i,k) = pidep(i,k)*factor
975 endif
976 !
977 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
978 ! update
979 q(i,k) = q(i,k)+work2(i,k)*dtcld
980 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
981 +psacw(i,k))*dtcld,0.)
982 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
983 +prevp(i,k))*dtcld,0.)
984 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
985 -pigen(i,k)-pidep(i,k))*dtcld,0.)
986 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
987 +psaci(i,k)+psacw(i,k))*dtcld,0.)
988 xlf = xls-xl(i,k)
989 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
990 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
991 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
992 else
993 !
994 ! cloud water
995 !
996 value = max(qmin,qci(i,k,1))
997 source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
998 if (source.gt.value) then
999 factor = value/source
1000 praut(i,k) = praut(i,k)*factor
1001 pracw(i,k) = pracw(i,k)*factor
1002 psacw(i,k) = psacw(i,k)*factor
1003 endif
1004 !
1005 ! snow
1006 !
1007 value = max(qcrmin,qrs(i,k,2))
1008 source=(-psevp(i,k))*dtcld
1009 if (source.gt.value) then
1010 factor = value/source
1011 psevp(i,k) = psevp(i,k)*factor
1012 endif
1013 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1014 ! update
1015 q(i,k) = q(i,k)+work2(i,k)*dtcld
1016 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1017 +psacw(i,k))*dtcld,0.)
1018 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1019 +prevp(i,k) +psacw(i,k))*dtcld,0.)
1020 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1021 xlf = xls-xl(i,k)
1022 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1023 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1024 endif
1025 enddo
1026 enddo
1027 !
1028 ! Inline expansion for fpvs
1029 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1030 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1031 hsub = xls
1032 hvap = xlv0
1033 cvap = cpv
1034 ttp=t0c+0.01
1035 dldt=cvap-cliq
1036 xa=-dldt/rv
1037 xb=xa+hvap/(rv*ttp)
1038 dldti=cvap-cice
1039 xai=-dldti/rv
1040 xbi=xai+hsub/(rv*ttp)
1041 do k = kts, kte
1042 do i = its, ite
1043 tr=ttp/t(i,k)
1044 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1045 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1046 qs(i,k,1) = max(qs(i,k,1),qmin)
1047 tr=ttp/t(i,k)
1048 if(t(i,k).lt.ttp) then
1049 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1050 else
1051 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1052 endif
1053 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1054 qs(i,k,2) = max(qs(i,k,2),qmin)
1055 enddo
1056 enddo
1057 !
1058 !----------------------------------------------------------------
1059 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1060 ! if there exists additional water vapor condensated/if
1061 ! evaporation of cloud water is not enough to remove subsaturation
1062 !
1063 do k = kts, kte
1064 do i = its, ite
1065 ! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1066 work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/ &
1067 (1.+(xl(i,k))*(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k))*(t(i,k)))))
1068 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1069 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1070 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1071 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1072 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1073 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1074 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1075 enddo
1076 enddo
1077 !
1078 !
1079 !----------------------------------------------------------------
1080 ! padding for small values
1081 !
1082 do k = kts, kte
1083 do i = its, ite
1084 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1085 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1086 enddo
1087 enddo
1088 enddo ! big loops
1089 END SUBROUTINE wsm52d
1090 ! ...................................................................
1091 REAL FUNCTION rgmma(x)
1092 !-------------------------------------------------------------------
1093 IMPLICIT NONE
1094 !-------------------------------------------------------------------
1095 ! rgmma function: use infinite product form
1096 REAL :: euler
1097 PARAMETER (euler=0.577215664901532)
1098 REAL :: x, y
1099 INTEGER :: i
1100 if(x.eq.1.)then
1101 rgmma=0.
1102 else
1103 rgmma=x*exp(euler*x)
1104 do i=1,10000
1105 y=float(i)
1106 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1107 enddo
1108 rgmma=1./rgmma
1109 endif
1110 END FUNCTION rgmma
1111 !
1112 !--------------------------------------------------------------------------
1113 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1114 !--------------------------------------------------------------------------
1115 IMPLICIT NONE
1116 !--------------------------------------------------------------------------
1117 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1118 xai,xbi,ttp,tr
1119 INTEGER ice
1120 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1121 ttp=t0c+0.01
1122 dldt=cvap-cliq
1123 xa=-dldt/rv
1124 xb=xa+hvap/(rv*ttp)
1125 dldti=cvap-cice
1126 xai=-dldti/rv
1127 xbi=xai+hsub/(rv*ttp)
1128 tr=ttp/t
1129 if(t.lt.ttp.and.ice.eq.1) then
1130 fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1131 else
1132 fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1133 endif
1134 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1135 END FUNCTION fpvs
1136 !-------------------------------------------------------------------
1137 SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
1138 !-------------------------------------------------------------------
1139 IMPLICIT NONE
1140 !-------------------------------------------------------------------
1141 !.... constants which may not be tunable
1142 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1143 LOGICAL, INTENT(IN) :: allowed_to_read
1144 REAL :: pi
1145 !
1146 pi = 4.*atan(1.)
1147 xlv1 = cl-cpv
1148 !
1149 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1150 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1151 !
1152 bvtr1 = 1.+bvtr
1153 bvtr2 = 2.5+.5*bvtr
1154 bvtr3 = 3.+bvtr
1155 bvtr4 = 4.+bvtr
1156 g1pbr = rgmma(bvtr1)
1157 g3pbr = rgmma(bvtr3)
1158 g4pbr = rgmma(bvtr4) ! 17.837825
1159 g5pbro2 = rgmma(bvtr2) ! 1.8273
1160 pvtr = avtr*g4pbr/6.
1161 eacrr = 1.0
1162 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1163 precr1 = 2.*pi*n0r*.78
1164 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
1165 xm0 = (di0/dicon)**2
1166 xmmax = (dimax/dicon)**2
1167 roqimax = 2.08e22*dimax**8
1168 !
1169 bvts1 = 1.+bvts
1170 bvts2 = 2.5+.5*bvts
1171 bvts3 = 3.+bvts
1172 bvts4 = 4.+bvts
1173 g1pbs = rgmma(bvts1) !.8875
1174 g3pbs = rgmma(bvts3)
1175 g4pbs = rgmma(bvts4) ! 12.0786
1176 g5pbso2 = rgmma(bvts2)
1177 pvts = avts*g4pbs/6.
1178 pacrs = pi*n0s*avts*g3pbs*.25
1179 precs1 = 4.*n0s*.65
1180 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1181 pidn0r = pi*denr*n0r
1182 pidn0s = pi*dens*n0s
1183 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1184 !
1185 rslopermax = 1./lamdarmax
1186 rslopesmax = 1./lamdasmax
1187 rsloperbmax = rslopermax ** bvtr
1188 rslopesbmax = rslopesmax ** bvts
1189 rsloper2max = rslopermax * rslopermax
1190 rslopes2max = rslopesmax * rslopesmax
1191 rsloper3max = rsloper2max * rslopermax
1192 rslopes3max = rslopes2max * rslopesmax
1193 !
1194 END SUBROUTINE wsm5init
1195 END MODULE module_mp_wsm5