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