| blob | ed8b608466507285e198c2fa2d20a3089333538e |
1 #if ( RWORDSIZE == 4 )
2 # define VREC vsrec
3 # define VSQRT vssqrt
4 #else
5 # define VREC vrec
6 # define VSQRT vsqrt
7 #endif
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
129 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
130 INTENT(INOUT) :: snow, &
131 snowncv
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
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 INTEGER, DIMENSION( its:ite ) :: kwork1, kwork2
232 REAL, DIMENSION( its:ite , kts:kte ) :: &
233 falkc, work1c, work2c, fallc
234 ! variables for optimization
235 REAL, DIMENSION( its:ite ) :: tvec1
236 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
237 LOGICAL, DIMENSION( its:ite ) :: flgcld
238 REAL :: pi, &
239 cpmcal, xlcal, lamdar, lamdas, diffus, &
240 viscos, xka, venfac, conden, diffac, &
241 x, y, z, a, b, c, d, e, &
242 fallsum, fallsum_qsi, vt2i,vt2s,acrfac, &
243 qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
244 snomlt, hold, holdrs, facqci, supcol, coeres, &
245 supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
246 qimax, diameter, xni0, roqi0, supice
247 REAL :: holdc, holdci
248 INTEGER :: i, j, k, mstepmax, &
249 iprt, latd, lond, loop, loops, ifsat, kk, n
250 ! Temporaries used for inlining fpvs function
251 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
252 !
253 !=================================================================
254 ! compute internal functions
255 !
256 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
257 xlcal(x) = xlv0-xlv1*(x-t0c)
258 !----------------------------------------------------------------
259 ! size distributions: (x=mixing ratio, y=air density):
260 ! valid for mixing ratio > 1.e-9 kg/kg.
261 !
262 ! Optimizatin : A**B => exp(log(A)*(B))
263 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
264 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
265 !
266 !----------------------------------------------------------------
267 ! diffus: diffusion coefficient of the water vapor
268 ! viscos: kinematic viscosity(m2s-1)
269 !
270 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
271 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
272 xka(x,y) = 1.414e3*viscos(x,y)*y
273 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
274 ! venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
275 ! /viscos(b,c)**(.5)*(den0/c)**0.25
276 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
277 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
278 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
279 !
280 pi = 4. * atan(1.)
281 !
282 !----------------------------------------------------------------
283 ! paddint 0 for negative values generated by dynamics
284 !
285 do k = kts, kte
286 do i = its, ite
287 qci(i,k) = max(qci(i,k),0.0)
288 qrs(i,k) = max(qrs(i,k),0.0)
289 enddo
290 enddo
291 !
292 !----------------------------------------------------------------
293 ! latent heat for phase changes and heat capacity. neglect the
294 ! changes during microphysical process calculation
295 ! emanuel(1994)
296 !
297 do k = kts, kte
298 do i = its, ite
299 cpm(i,k) = cpmcal(q(i,k))
300 xl(i,k) = xlcal(t(i,k))
301 enddo
302 enddo
303 !
304 !----------------------------------------------------------------
305 ! compute the minor time steps.
306 !
307 loops = max(nint(delt/dtcldcr),1)
308 dtcld = delt/loops
309 if(delt.le.dtcldcr) dtcld = delt
310 !
311 do loop = 1,loops
312 !
313 !----------------------------------------------------------------
314 ! initialize the large scale variables
315 !
316 do i = its, ite
317 mstep(i) = 1
318 flgcld(i) = .true.
319 enddo
320 !
321 ! do k = kts, kte
322 ! do i = its, ite
323 ! denfac(i,k) = sqrt(den0/den(i,k))
324 ! enddo
325 ! enddo
326 do k = kts, kte
327 CALL VREC( tvec1(its), den(its,k), ite-its+1)
328 do i = its, ite
329 tvec1(i) = tvec1(i)*den0
330 enddo
331 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
332 enddo
333 !
334 ! Inline expansion for fpvs
335 ! qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
336 ! qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
337 cvap = cpv
338 hvap=xlv0
339 hsub=xls
340 ttp=t0c+0.01
341 dldt=cvap-cliq
342 xa=-dldt/rv
343 xb=xa+hvap/(rv*ttp)
344 dldti=cvap-cice
345 xai=-dldti/rv
346 xbi=xai+hsub/(rv*ttp)
347 do k = kts, kte
348 do i = its, ite
349 ! tr=ttp/t(i,k)
350 ! if(t(i,k).lt.ttp) then
351 ! qs(i,k) =psat*(tr**xai)*exp(xbi*(1.-tr))
352 ! else
353 ! qs(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
354 ! endif
355 ! qs0(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
356 tr=ttp/t(i,k)
357 if(t(i,k).lt.ttp) then
358 qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
359 else
360 qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
361 endif
362 qs0(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
363 qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
364 qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
365 qs(i,k) = max(qs(i,k),qmin)
366 rh(i,k) = max(q(i,k) / qs(i,k),qmin)
367 enddo
368 enddo
369 !
370 !----------------------------------------------------------------
371 ! initialize the variables for microphysical physics
372 !
373 !
374 do k = kts, kte
375 do i = its, ite
376 pres(i,k) = 0.
377 paut(i,k) = 0.
378 pacr(i,k) = 0.
379 pgen(i,k) = 0.
380 pisd(i,k) = 0.
381 pcon(i,k) = 0.
382 fall(i,k) = 0.
383 falk(i,k) = 0.
384 fallc(i,k) = 0.
385 falkc(i,k) = 0.
386 xni(i,k) = 1.e3
387 enddo
388 enddo
389 !
390 !----------------------------------------------------------------
391 ! compute the fallout term:
392 ! first, vertical terminal velosity for minor loops
393 !---------------------------------------------------------------
394 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
395 !---------------------------------------------------------------
396 do k = kts, kte
397 do i = its, ite
398 supcol = t0c-t(i,k)
399 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
400 if(t(i,k).ge.t0c) then
401 if(qrs(i,k).le.qcrmin)then
402 rslope(i,k) = rslopermax
403 rslopeb(i,k) = rsloperbmax
404 rslope2(i,k) = rsloper2max
405 rslope3(i,k) = rsloper3max
406 else
407 rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
408 ! rslopeb(i,k) = rslope(i,k)**bvtr
409 rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
410 rslope2(i,k) = rslope(i,k)*rslope(i,k)
411 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
412 endif
413 else
414 if(qrs(i,k).le.qcrmin)then
415 rslope(i,k) = rslopesmax
416 rslopeb(i,k) = rslopesbmax
417 rslope2(i,k) = rslopes2max
418 rslope3(i,k) = rslopes3max
419 else
420 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
421 ! rslopeb(i,k) = rslope(i,k)**bvts
422 rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
423 rslope2(i,k) = rslope(i,k)*rslope(i,k)
424 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
425 endif
426 endif
427 !-------------------------------------------------------------
428 ! Ni: ice crystal number concentraiton [HDC 5c]
429 !-------------------------------------------------------------
430 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
431 ! *max(qci(i,k),qmin))**0.75,1.e3),1.e6)
432 xni(i,k) = min(max(5.38e7*exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
433 enddo
434 enddo
435 !
436 mstepmax = 1
437 numdt = 1
438 do k = kte, kts, -1
439 do i = its, ite
440 if(t(i,k).lt.t0c) then
441 pvt = pvts
442 else
443 pvt = pvtr
444 endif
445 work1(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
446 work2(i,k) = work1(i,k)/delz(i,k)
447 numdt(i) = max(nint(work2(i,k)*dtcld+.5),1)
448 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
449 enddo
450 enddo
451 do i = its, ite
452 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
453 enddo
454 !
455 do n = 1, mstepmax
456 k = kte
457 do i = its, ite
458 if(n.le.mstep(i)) then
459 falk(i,k) = den(i,k)*qrs(i,k)*work2(i,k)/mstep(i)
460 hold = falk(i,k)
461 fall(i,k) = fall(i,k)+falk(i,k)
462 holdrs = qrs(i,k)
463 qrs(i,k) = max(qrs(i,k)-falk(i,k)*dtcld/den(i,k),0.)
464 endif
465 enddo
466 do k = kte-1, kts, -1
467 do i = its, ite
468 if(n.le.mstep(i)) then
469 falk(i,k) = den(i,k)*qrs(i,k)*work2(i,k)/mstep(i)
470 hold = falk(i,k)
471 fall(i,k) = fall(i,k)+falk(i,k)
472 holdrs = qrs(i,k)
473 qrs(i,k) = max(qrs(i,k)-(falk(i,k) &
474 -falk(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
475 endif
476 enddo
477 enddo
478 enddo
479 !---------------------------------------------------------------
480 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
481 !---------------------------------------------------------------
482 mstepmax = 1
483 mstep = 1
484 numdt = 1
485 do k = kte, kts, -1
486 do i = its, ite
487 if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
488 xmi = den(i,k)*qci(i,k)/xni(i,k)
489 ! diameter = dicon * sqrt(xmi)
490 ! work1c(i,k) = 1.49e4*diameter**1.31
491 diameter = max(dicon * sqrt(xmi), 1.e-25)
492 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
493 else
494 work1c(i,k) = 0.
495 endif
496 if(qci(i,k).le.0.) then
497 work2c(i,k) = 0.
498 else
499 work2c(i,k) = work1c(i,k)/delz(i,k)
500 endif
501 numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
502 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
503 enddo
504 enddo
505 do i = its, ite
506 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
507 enddo
508 !
509 do n = 1, mstepmax
510 k = kte
511 do i = its, ite
512 if (n.le.mstep(i)) then
513 falkc(i,k) = den(i,k)*qci(i,k)*work2c(i,k)/mstep(i)
514 holdc = falkc(i,k)
515 fallc(i,k) = fallc(i,k)+falkc(i,k)
516 holdci = qci(i,k)
517 qci(i,k) = max(qci(i,k)-falkc(i,k)*dtcld/den(i,k),0.)
518 endif
519 enddo
520 do k = kte-1, kts, -1
521 do i = its, ite
522 if (n.le.mstep(i)) then
523 falkc(i,k) = den(i,k)*qci(i,k)*work2c(i,k)/mstep(i)
524 holdc = falkc(i,k)
525 fallc(i,k) = fallc(i,k)+falkc(i,k)
526 holdci = qci(i,k)
527 qci(i,k) = max(qci(i,k)-(falkc(i,k) &
528 -falkc(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
529 endif
530 enddo
531 enddo
532 enddo
533 !
534 !----------------------------------------------------------------
535 ! compute the freezing/melting term. [D89 B16-B17]
536 ! freezing occurs one layer above the melting level
537 !
538 do i = its, ite
539 mstep(i) = 0
540 enddo
541 do k = kts, kte
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 kwork2(i) = mstep(i)
551 kwork1(i) = mstep(i)
552 if(mstep(i).ne.0) then
553 if (w(i,mstep(i)).gt.0.) then
554 kwork1(i) = mstep(i) + 1
555 endif
556 endif
557 enddo
558 !
559 do i = its, ite
560 k = kwork1(i)
561 kk = kwork2(i)
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