module_mp_wsm6.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_wsm6
10 !
11 !
12 REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
13 REAL, PARAMETER, PRIVATE :: n0r = 8.e6
14 REAL, PARAMETER, PRIVATE :: n0g = 4.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 :: avtg = 330.
24 REAL, PARAMETER, PRIVATE :: bvtg = 0.8
25 REAL, PARAMETER, PRIVATE :: deng = 500.
26 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
27 REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
28 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
29 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
30 REAL, PARAMETER, PRIVATE :: betai = .6
31 REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
32 REAL, PARAMETER, PRIVATE :: dicon = 11.9
33 REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
34 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
35 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
36 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
37 REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
38 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
39 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
40 REAL, PARAMETER, PRIVATE :: t40c = 233.16
41 REAL, PARAMETER, PRIVATE :: eacrc = 1.0
42 REAL, PARAMETER, PRIVATE :: dens = 100.0
43 REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! pgaut
44 REAL, SAVE :: &
45 qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
46 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
47 bvtr6,g6pbr, &
48 precr1,precr2,xm0,xmmax,roqimax,bvts1, &
49 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
50 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
51 pidn0s,xlv1,pacrc, &
52 bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
53 g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
54 precg1,precg2,pidn0g, &
55 rslopermax,rslopesmax,rslopegmax, &
56 rsloperbmax,rslopesbmax,rslopegbmax, &
57 rsloper2max,rslopes2max,rslopeg2max, &
58 rsloper3max,rslopes3max,rslopeg3max
59 CONTAINS
60 !===================================================================
61 !
62 SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
63 ,den, pii, p, delz &
64 ,delt,g, cpd, cpv, rd, rv, t0c &
65 ,ep1, ep2, qmin &
66 ,XLS, XLV0, XLF0, den0, denr &
67 ,cliq,cice,psat &
68 ,rain, rainncv &
69 ,snow, snowncv &
70 ,graupel, graupelncv &
71 ,sr &
72 ,ids,ide, jds,jde, kds,kde &
73 ,ims,ime, jms,jme, kms,kme &
74 ,its,ite, jts,jte, kts,kte &
75 )
76 !-------------------------------------------------------------------
77 IMPLICIT NONE
78 !-------------------------------------------------------------------
79 !
80 ! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the WRF
81 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
82 ! number concentration is a function of temperature, and seperate assumption
83 ! is developed, in which ice crystal number concentration is a function
84 ! of ice amount. A theoretical background of the ice-microphysics and related
85 ! processes in the WSMMPs are described in Hong et al. (2004).
86 ! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
87 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
88 !
89 ! WSM6 cloud scheme
90 !
91 ! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
92 ! Summer 2003
93 !
94 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
95 ! Summer 2004
96 !
97 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
98 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
99 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
100 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
101 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
102 !
103 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
104 ims,ime, jms,jme, kms,kme , &
105 its,ite, jts,jte, kts,kte
106 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
107 INTENT(INOUT) :: &
108 th, &
109 q, &
110 qc, &
111 qi, &
112 qr, &
113 qs, &
114 qg
115 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
116 INTENT(IN ) :: &
117 den, &
118 pii, &
119 p, &
120 delz
121 REAL, INTENT(IN ) :: delt, &
122 g, &
123 rd, &
124 rv, &
125 t0c, &
126 den0, &
127 cpd, &
128 cpv, &
129 ep1, &
130 ep2, &
131 qmin, &
132 XLS, &
133 XLV0, &
134 XLF0, &
135 cliq, &
136 cice, &
137 psat, &
138 denr
139 REAL, DIMENSION( ims:ime , jms:jme ), &
140 INTENT(INOUT) :: rain, &
141 rainncv, &
142 sr
143
144 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
145 INTENT(INOUT) :: snow, &
146 snowncv
147
148 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
149 INTENT(INOUT) :: graupel, &
150 graupelncv
151 ! LOCAL VAR
152 REAL, DIMENSION( its:ite , kts:kte ) :: t
153 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
154 REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
155 INTEGER :: i,j,k
156 !-------------------------------------------------------------------
157 DO j=jts,jte
158 DO k=kts,kte
159 DO i=its,ite
160 t(i,k)=th(i,k,j)*pii(i,k,j)
161 qci(i,k,1) = qc(i,k,j)
162 qci(i,k,2) = qi(i,k,j)
163 qrs(i,k,1) = qr(i,k,j)
164 qrs(i,k,2) = qs(i,k,j)
165 qrs(i,k,3) = qg(i,k,j)
166 ENDDO
167 ENDDO
168 CALL wsm62D(t, q(ims,kms,j), qci, qrs &
169 ,den(ims,kms,j) &
170 ,p(ims,kms,j), delz(ims,kms,j) &
171 ,delt,g, cpd, cpv, rd, rv, t0c &
172 ,ep1, ep2, qmin &
173 ,XLS, XLV0, XLF0, den0, denr &
174 ,cliq,cice,psat &
175 ,j &
176 ,rain(ims,j),rainncv(ims,j) &
177 ,sr(ims,j) &
178 ,ids,ide, jds,jde, kds,kde &
179 ,ims,ime, jms,jme, kms,kme &
180 ,its,ite, jts,jte, kts,kte &
181 ,snow(ims,j),snowncv(ims,j) &
182 ,graupel(ims,j),graupelncv(ims,j) &
183 )
184 DO K=kts,kte
185 DO I=its,ite
186 th(i,k,j)=t(i,k)/pii(i,k,j)
187 qc(i,k,j) = qci(i,k,1)
188 qi(i,k,j) = qci(i,k,2)
189 qr(i,k,j) = qrs(i,k,1)
190 qs(i,k,j) = qrs(i,k,2)
191 qg(i,k,j) = qrs(i,k,3)
192 ENDDO
193 ENDDO
194 ENDDO
195 END SUBROUTINE wsm6
196 !===================================================================
197 !
198 SUBROUTINE wsm62D(t, q, qci, qrs, den, p, delz &
199 ,delt,g, cpd, cpv, rd, rv, t0c &
200 ,ep1, ep2, qmin &
201 ,XLS, XLV0, XLF0, den0, denr &
202 ,cliq,cice,psat &
203 ,lat &
204 ,rain,rainncv &
205 ,sr &
206 ,ids,ide, jds,jde, kds,kde &
207 ,ims,ime, jms,jme, kms,kme &
208 ,its,ite, jts,jte, kts,kte &
209 ,snow,snowncv &
210 ,graupel,graupelncv &
211 )
212 !-------------------------------------------------------------------
213 IMPLICIT NONE
214 !-------------------------------------------------------------------
215 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
216 ims,ime, jms,jme, kms,kme , &
217 its,ite, jts,jte, kts,kte, &
218 lat
219 REAL, DIMENSION( its:ite , kts:kte ), &
220 INTENT(INOUT) :: &
221 t
222 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
223 INTENT(INOUT) :: &
224 qci
225 REAL, DIMENSION( its:ite , kts:kte, 3 ), &
226 INTENT(INOUT) :: &
227 qrs
228 REAL, DIMENSION( ims:ime , kms:kme ), &
229 INTENT(INOUT) :: &
230 q
231 REAL, DIMENSION( ims:ime , kms:kme ), &
232 INTENT(IN ) :: &
233 den, &
234 p, &
235 delz
236 REAL, INTENT(IN ) :: delt, &
237 g, &
238 cpd, &
239 cpv, &
240 t0c, &
241 den0, &
242 rd, &
243 rv, &
244 ep1, &
245 ep2, &
246 qmin, &
247 XLS, &
248 XLV0, &
249 XLF0, &
250 cliq, &
251 cice, &
252 psat, &
253 denr
254 REAL, DIMENSION( ims:ime ), &
255 INTENT(INOUT) :: rain, &
256 rainncv, &
257 sr
258 REAL, DIMENSION( ims:ime ), OPTIONAL, &
259 INTENT(INOUT) :: snow, &
260 snowncv
261
262 REAL, DIMENSION( ims:ime ), OPTIONAL, &
263 INTENT(INOUT) :: graupel, &
264 graupelncv
265 ! LOCAL VAR
266 REAL, DIMENSION( its:ite , kts:kte , 3) :: &
267 rh, qs, rslope, rslope2, rslope3, rslopeb, &
268 falk, fall, work1
269 REAL, DIMENSION( its:ite , kts:kte , 1) :: &
270 worka
271 REAL, DIMENSION( its:ite , kts:kte ) :: &
272 falkc, work1c, work2c, fallc
273 REAL, DIMENSION( its:ite , kts:kte ) :: &
274 prevp, psdep, pgdep, praut, psaut, pgaut, &
275 pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci, &
276 piacr, pracs, psacr, pgaci, pseml, pgeml
277 REAL, DIMENSION( its:ite , kts:kte) :: qsum
278 REAL, DIMENSION( its:ite , kts:kte ) :: paacw
279 REAL, DIMENSION( its:ite , kts:kte ) :: &
280 pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, &
281 xni, pgevp,n0sfac
282 ! variables for optimization
283 REAL, DIMENSION( its:ite ) :: tvec1
284 REAL :: temp
285 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
286 LOGICAL, DIMENSION( its:ite ) :: flgcld
287 REAL :: pi, &
288 cpmcal, xlcal, lamdar, lamdas, lamdag, diffus, &
289 viscos, xka, venfac, conden, diffac, &
290 x, y, z, a, b, c, d, e, &
291 qdt, holdrr, holdrs, holdrg, supcol, pvt, &
292 coeres, supsat, dtcld, xmi, eacrs, satdt, &
293 qimax, diameter, xni0, roqi0, &
294 fallsum, fallsum_qsi, fallsum_qg, &
295 vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
296 xlwork2, factor, source, value, &
297 xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
298 REAL :: vt2ave
299 REAL :: holdc, holdci
300 INTEGER :: i, j, k, mstepmax, &
301 iprt, latd, lond, loop, loops, ifsat, n
302 ! Temporaries used for inlining fpvs function
303 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
304 !
305 !=================================================================
306 ! compute internal functions
307 !
308 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
309 xlcal(x) = xlv0-xlv1*(x-t0c)
310 !----------------------------------------------------------------
311 ! size distributions: (x=mixing ratio, y=air density):
312 ! valid for mixing ratio > 1.e-9 kg/kg.
313 !
314 ! Optimizatin : A**B => exp(log(A)*(B))
315 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
316 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
317 lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
318 !
319 !----------------------------------------------------------------
320 ! diffus: diffusion coefficient of the water vapor
321 ! viscos: kinematic viscosity(m2s-1)
322 !
323 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
324 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
325 xka(x,y) = 1.414e3*viscos(x,y)*y
326 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
327 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
328 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
329 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
330 !
331 pi = 4. * atan(1.)
332 !
333 !
334 !----------------------------------------------------------------
335 ! paddint 0 for negative values generated by dynamics
336 !
337 do k = kts, kte
338 do i = its, ite
339 qci(i,k,1) = max(qci(i,k,1),0.0)
340 qrs(i,k,1) = max(qrs(i,k,1),0.0)
341 qci(i,k,2) = max(qci(i,k,2),0.0)
342 qrs(i,k,2) = max(qrs(i,k,2),0.0)
343 qrs(i,k,3) = max(qrs(i,k,3),0.0)
344 enddo
345 enddo
346 !
347 !----------------------------------------------------------------
348 ! latent heat for phase changes and heat capacity. neglect the
349 ! changes during microphysical process calculation
350 ! emanuel(1994)
351 !
352 do k = kts, kte
353 do i = its, ite
354 cpm(i,k) = cpmcal(q(i,k))
355 xl(i,k) = xlcal(t(i,k))
356 enddo
357 enddo
358 !
359 !----------------------------------------------------------------
360 ! compute the minor time steps.
361 !
362 loops = max(nint(delt/dtcldcr),1)
363 dtcld = delt/loops
364 if(delt.le.dtcldcr) dtcld = delt
365 !
366 do loop = 1,loops
367 !
368 !----------------------------------------------------------------
369 ! initialize the large scale variables
370 !
371 do i = its, ite
372 mstep(i) = 1
373 flgcld(i) = .true.
374 enddo
375 !
376 ! do k = kts, kte
377 ! do i = its, ite
378 ! denfac(i,k) = sqrt(den0/den(i,k))
379 ! enddo
380 ! enddo
381 do k = kts, kte
382 CALL VREC( tvec1(its), den(its,k), ite-its+1)
383 do i = its, ite
384 tvec1(i) = tvec1(i)*den0
385 enddo
386 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
387 enddo
388 !
389 ! Inline expansion for fpvs
390 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
391 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
392 hsub = xls
393 hvap = xlv0
394 cvap = cpv
395 ttp=t0c+0.01
396 dldt=cvap-cliq
397 xa=-dldt/rv
398 xb=xa+hvap/(rv*ttp)
399 dldti=cvap-cice
400 xai=-dldti/rv
401 xbi=xai+hsub/(rv*ttp)
402 do k = kts, kte
403 do i = its, ite
404 tr=ttp/t(i,k)
405 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
406 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
407 qs(i,k,1) = max(qs(i,k,1),qmin)
408 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
409 tr=ttp/t(i,k)
410 if(t(i,k).lt.ttp) then
411 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
412 else
413 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
414 endif
415 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
416 qs(i,k,2) = max(qs(i,k,2),qmin)
417 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
418 enddo
419 enddo
420 !
421 !----------------------------------------------------------------
422 ! initialize the variables for microphysical physics
423 !
424 !
425 do k = kts, kte
426 do i = its, ite
427 prevp(i,k) = 0.
428 psdep(i,k) = 0.
429 pgdep(i,k) = 0.
430 praut(i,k) = 0.
431 psaut(i,k) = 0.
432 pgaut(i,k) = 0.
433 pracw(i,k) = 0.
434 praci(i,k) = 0.
435 piacr(i,k) = 0.
436 psaci(i,k) = 0.
437 psacw(i,k) = 0.
438 pracs(i,k) = 0.
439 psacr(i,k) = 0.
440 pgacw(i,k) = 0.
441 paacw(i,k) = 0.
442 pgaci(i,k) = 0.
443 pgacr(i,k) = 0.
444 pgacs(i,k) = 0.
445 pigen(i,k) = 0.
446 pidep(i,k) = 0.
447 pcond(i,k) = 0.
448 psmlt(i,k) = 0.
449 pgmlt(i,k) = 0.
450 pseml(i,k) = 0.
451 pgeml(i,k) = 0.
452 psevp(i,k) = 0.
453 pgevp(i,k) = 0.
454 falk(i,k,1) = 0.
455 falk(i,k,2) = 0.
456 falk(i,k,3) = 0.
457 fall(i,k,1) = 0.
458 fall(i,k,2) = 0.
459 fall(i,k,3) = 0.
460 fallc(i,k) = 0.
461 falkc(i,k) = 0.
462 xni(i,k) = 1.e3
463 enddo
464 enddo
465 !
466 !----------------------------------------------------------------
467 ! compute the fallout term:
468 ! first, vertical terminal velosity for minor loops
469 !
470 do k = kts, kte
471 do i = its, ite
472 supcol = t0c-t(i,k)
473 !---------------------------------------------------------------
474 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
475 !---------------------------------------------------------------
476 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
477 if(qrs(i,k,1).le.qcrmin)then
478 rslope(i,k,1) = rslopermax
479 rslopeb(i,k,1) = rsloperbmax
480 rslope2(i,k,1) = rsloper2max
481 rslope3(i,k,1) = rsloper3max
482 else
483 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
484 rslopeb(i,k,1) = rslope(i,k,1)**bvtr
485 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
486 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
487 endif
488 if(qrs(i,k,2).le.qcrmin)then
489 rslope(i,k,2) = rslopesmax
490 rslopeb(i,k,2) = rslopesbmax
491 rslope2(i,k,2) = rslopes2max
492 rslope3(i,k,2) = rslopes3max
493 else
494 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
495 rslopeb(i,k,2) = rslope(i,k,2)**bvts
496 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
497 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
498 endif
499 if(qrs(i,k,3).le.qcrmin)then
500 rslope(i,k,3) = rslopegmax
501 rslopeb(i,k,3) = rslopegbmax
502 rslope2(i,k,3) = rslopeg2max
503 rslope3(i,k,3) = rslopeg3max
504 else
505 rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
506 rslopeb(i,k,3) = rslope(i,k,3)**bvtg
507 rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
508 rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
509 endif
510 !-------------------------------------------------------------
511 ! Ni: ice crystal number concentraiton [HDC 5c]
512 !-------------------------------------------------------------
513 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
514 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
515 temp = (den(i,k)*max(qci(i,k,2),qmin))
516 temp = sqrt(sqrt(temp*temp*temp))
517 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
518 enddo
519 enddo
520 !
521 mstepmax = 1
522 numdt = 1
523 do k = kte, kts, -1
524 do i = its, ite
525 work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
526 work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
527 work1(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)/delz(i,k)
528 qsum(i,k) = qrs(i,k,2)+qrs(i,k,3)
529 if(qsum(i,k) .le. 0) then
530 qsum(i,k) = 1.E-15
531 endif
532 worka(i,k,1) = (pvts*rslopeb(i,k,2)*denfac(i,k)*qrs(i,k,2) &
533 +pvtg*rslopeb(i,k,3)*denfac(i,k)*qrs(i,k,3))/qsum(i,k)/delz(i,k)
534 numdt(i) = max(nint(max(work1(i,k,1),worka(i,k,1)) &
535 *dtcld+.5),1)
536 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
537 enddo
538 enddo
539 do i = its, ite
540 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
541 enddo
542 !
543 do n = 1, mstepmax
544 k = kte
545 do i = its, ite
546 if(n.le.mstep(i)) then
547 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
548 falk(i,k,2) = den(i,k)*qrs(i,k,2)*worka(i,k,1)/mstep(i)
549 falk(i,k,3) = den(i,k)*qrs(i,k,3)*worka(i,k,1)/mstep(i)
550 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
551 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
552 fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
553 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
554 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
555 qrs(i,k,3) = max(qrs(i,k,3)-falk(i,k,3)*dtcld/den(i,k),0.)
556 endif
557 enddo
558 do k = kte-1, kts, -1
559 do i = its, ite
560 if(n.le.mstep(i)) then
561 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
562 falk(i,k,2) = den(i,k)*qrs(i,k,2)*worka(i,k,1)/mstep(i)
563 falk(i,k,3) = den(i,k)*qrs(i,k,3)*worka(i,k,1)/mstep(i)
564 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
565 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
566 fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
567 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
568 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
569 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
570 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
571 qrs(i,k,3) = max(qrs(i,k,3)-(falk(i,k,3)-falk(i,k+1,3) &
572 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
573 endif
574 enddo
575 enddo
576 do k = kte, kts, -1
577 do i = its, ite
578 if(n.le.mstep(i).and.t(i,k).gt.t0c) then
579 !---------------------------------------------------------------
580 ! psmlt: melting of snow [HL A33] [RH83 A25]
581 ! (T>T0: S->R)
582 !---------------------------------------------------------------
583 xlf = xlf0
584 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
585 if(qrs(i,k,2).gt.0.) then
586 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
587 psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
588 *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
589 +precs2*work2(i,k)*coeres)
590 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
591 -qrs(i,k,2)/mstep(i)),0.)
592 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
593 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
594 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
595 endif
596 !---------------------------------------------------------------
597 ! pgmlt: melting of graupel [HL A23] [LFO 47]
598 ! (T>T0: G->R)
599 !---------------------------------------------------------------
600 if(qrs(i,k,3).gt.0.) then
601 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
602 pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
603 *(t0c-t(i,k))*(precg1*rslope2(i,k,3) &
604 +precg2*work2(i,k)*coeres)
605 pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
606 -qrs(i,k,3)/mstep(i)),0.)
607 qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
608 qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
609 t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
610 endif
611 endif
612 enddo
613 enddo
614 enddo
615 !---------------------------------------------------------------
616 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
617 !---------------------------------------------------------------
618 mstepmax = 1
619 mstep = 1
620 numdt = 1
621 do k = kte, kts, -1
622 do i = its, ite
623 if(qci(i,k,2).le.0.) then
624 work2c(i,k) = 0.
625 else
626 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
627 ! diameter = min(dicon * sqrt(xmi),dimax)
628 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
629 work1c(i,k) = 1.49e4*diameter**1.31
630 work2c(i,k) = work1c(i,k)/delz(i,k)
631 endif
632 numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
633 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
634 enddo
635 enddo
636 do i = its, ite
637 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
638 enddo
639 !
640 do n = 1, mstepmax
641 k = kte
642 do i = its, ite
643 if(n.le.mstep(i)) then
644 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
645 holdc = falkc(i,k)
646 fallc(i,k) = fallc(i,k)+falkc(i,k)
647 holdci = qci(i,k,2)
648 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
649 endif
650 enddo
651 do k = kte-1, kts, -1
652 do i = its, ite
653 if(n.le.mstep(i)) then
654 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
655 holdc = falkc(i,k)
656 fallc(i,k) = fallc(i,k)+falkc(i,k)
657 holdci = qci(i,k,2)
658 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
659 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
660 endif
661 enddo
662 enddo
663 enddo
664 !
665 !----------------------------------------------------------------
666 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
667 !
668 do i = its, ite
669 fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
670 fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
671 fallsum_qg = fall(i,kts,3)
672 rainncv(i) = 0.
673 if(fallsum.gt.0.) then
674 rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000.
675 rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
676 endif
677 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
678 snowncv(i) = 0.
679 if(fallsum_qsi.gt.0.) then
680 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
681 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
682 endif
683 ENDIF
684 IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
685 graupelncv(i) = 0.
686 if(fallsum_qg.gt.0.) then
687 graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.
688 graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
689 endif
690 ENDIF
691 sr(i) = 0.
692 if(fallsum.gt.0.)sr(i)=(fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + &
693 fallsum_qg*delz(i,kts)/denr*dtcld*1000.)/(rainncv(i)+1.e-12)
694 enddo
695 !
696 !---------------------------------------------------------------
697 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
698 ! (T>T0: I->C)
699 !---------------------------------------------------------------
700 do k = kts, kte
701 do i = its, ite
702 supcol = t0c-t(i,k)
703 xlf = xls-xl(i,k)
704 if(supcol.lt.0.) xlf = xlf0
705 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
706 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
707 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
708 qci(i,k,2) = 0.
709 endif
710 !---------------------------------------------------------------
711 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
712 ! (T<-40C: C->I)
713 !---------------------------------------------------------------
714 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
715 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
716 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
717 qci(i,k,1) = 0.
718 endif
719 !---------------------------------------------------------------
720 ! pihtf: heterogeneous freezing of cloud water [HL A44]
721 ! (T0>T>-40C: C->I)
722 !---------------------------------------------------------------
723 if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
724 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
725 ! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
726 pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
727 *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
728 qci(i,k,2) = qci(i,k,2) + pfrzdtc
729 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
730 qci(i,k,1) = qci(i,k,1)-pfrzdtc
731 endif
732 !---------------------------------------------------------------
733 ! pgfrz: freezing of rain water [HL A20] [LFO 45]
734 ! (T<T0, R->G)
735 !---------------------------------------------------------------
736 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
737 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
738 ! *(exp(pfrz2*supcol)-1.)*rslope3(i,k,1)**2 &
739 ! *rslope(i,k,1)*dtcld,qrs(i,k,1))
740 temp = rslope3(i,k,1)
741 temp = temp*temp*rslope(i,k,1)
742 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
743 *(exp(pfrz2*supcol)-1.)*temp*dtcld, &
744 qrs(i,k,1))
745 qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
746 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
747 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
748 endif
749 enddo
750 enddo
751 !
752 !
753 !----------------------------------------------------------------
754 ! rsloper: reverse of the slope parameter of the rain(m)
755 ! xka: thermal conductivity of air(jm-1s-1k-1)
756 ! work1: the thermodynamic term in the denominator associated with
757 ! heat conduction and vapor diffusion
758 ! (ry88, y93, h85)
759 ! work2: parameter associated with the ventilation effects(y93)
760 !
761 do k = kts, kte
762 do i = its, ite
763 if(qrs(i,k,1).le.qcrmin)then
764 rslope(i,k,1) = rslopermax
765 rslopeb(i,k,1) = rsloperbmax
766 rslope2(i,k,1) = rsloper2max
767 rslope3(i,k,1) = rsloper3max
768 else
769 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
770 rslopeb(i,k,1) = rslope(i,k,1)**bvtr
771 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
772 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
773 endif
774 if(qrs(i,k,2).le.qcrmin)then
775 rslope(i,k,2) = rslopesmax
776 rslopeb(i,k,2) = rslopesbmax
777 rslope2(i,k,2) = rslopes2max
778 rslope3(i,k,2) = rslopes3max
779 else
780 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
781 rslopeb(i,k,2) = rslope(i,k,2)**bvts
782 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
783 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
784 endif
785 if(qrs(i,k,3).le.qcrmin)then
786 rslope(i,k,3) = rslopegmax
787 rslopeb(i,k,3) = rslopegbmax
788 rslope2(i,k,3) = rslopeg2max
789 rslope3(i,k,3) = rslopeg3max
790 else
791 rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
792 rslopeb(i,k,3) = rslope(i,k,3)**bvtg
793 rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
794 rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
795 endif
796 enddo
797 enddo
798 !
799 do k = kts, kte
800 do i = its, ite
801 work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
802 work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
803 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
804 enddo
805 enddo
806 !
807 !===============================================================
808 !
809 ! warm rain processes
810 !
811 ! - follows the processes in RH83 and LFO except for autoconcersion
812 !
813 !===============================================================
814 !
815 do k = kts, kte
816 do i = its, ite
817 supsat = max(q(i,k),qmin)-qs(i,k,1)
818 satdt = supsat/dtcld
819 !---------------------------------------------------------------
820 ! praut: auto conversion rate from cloud to rain [HDC 16]
821 ! (C->R)
822 !---------------------------------------------------------------
823 if(qci(i,k,1).gt.qc0) then
824 praut(i,k) = qck1*qci(i,k,1)**(7./3.)
825 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
826 endif
827 !---------------------------------------------------------------
828 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
829 ! (C->R)
830 !---------------------------------------------------------------
831 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
832 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
833 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
834 endif
835 !---------------------------------------------------------------
836 ! prevp: evaporation/condensation rate of rain [HDC 14]
837 ! (V->R or R->V)
838 !---------------------------------------------------------------
839 if(qrs(i,k,1).gt.0.) then
840 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
841 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
842 +precr2*work2(i,k)*coeres)/work1(i,k,1)
843 if(prevp(i,k).lt.0.) then
844 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
845 prevp(i,k) = max(prevp(i,k),satdt/2)
846 else
847 prevp(i,k) = min(prevp(i,k),satdt/2)
848 endif
849 endif
850 enddo
851 enddo
852 !
853 !===============================================================
854 !
855 ! cold rain processes
856 !
857 ! - follows the revised ice microphysics processes in HDC
858 ! - the processes same as in RH83 and RH84 and LFO behave
859 ! following ice crystal hapits defined in HDC, inclduing
860 ! intercept parameter for snow (n0s), ice crystal number
861 ! concentration (ni), ice nuclei number concentration
862 ! (n0i), ice diameter (d)
863 !
864 !===============================================================
865 !
866 do k = kts, kte
867 do i = its, ite
868 supcol = t0c-t(i,k)
869 supsat = max(q(i,k),qmin)-qs(i,k,2)
870 satdt = supsat/dtcld
871 ifsat = 0
872 !-------------------------------------------------------------
873 ! Ni: ice crystal number concentraiton [HDC 5c]
874 !-------------------------------------------------------------
875 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
876 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
877 temp = (den(i,k)*max(qci(i,k,2),qmin))
878 temp = sqrt(sqrt(temp*temp*temp))
879 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
880 eacrs = exp(0.07*(-supcol))
881 !
882 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
883 diameter = min(dicon * sqrt(xmi),dimax)
884 vt2i = 1.49e4*diameter**1.31
885 vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
886 vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
887 vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
888 qsum(i,k) = qrs(i,k,2)+qrs(i,k,3)
889 if(qsum(i,k) .le. 0) then
890 qsum(i,k) = 1.E-15
891 endif
892 vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
893 if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
894 if(qrs(i,k,1).gt.qcrmin) then
895 !-------------------------------------------------------------
896 ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
897 ! (T<T0: I->R)
898 !-------------------------------------------------------------
899 acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1) &
900 +diameter**2*rslope(i,k,1)
901 praci(i,k) = pi*qci(i,k,2)*n0r*abs(vt2r-vt2i)*acrfac/4.
902 praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
903 !-------------------------------------------------------------
904 ! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
905 ! (T<T0: R->S or R->G)
906 !-------------------------------------------------------------
907 piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k) &
908 *g6pbr*rslope3(i,k,1)*rslope3(i,k,1) &
909 *rslopeb(i,k,1)/24./den(i,k)
910 piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
911 endif
912 !-------------------------------------------------------------
913 ! psaci: Accretion of cloud ice by snow [HDC 10]
914 ! (T<T0: I->S)
915 !-------------------------------------------------------------
916 if(qrs(i,k,2).gt.qcrmin) then
917 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
918 +diameter**2*rslope(i,k,2)
919 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
920 *abs(vt2ave-vt2i)*acrfac/4.
921 psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
922 endif
923 !-------------------------------------------------------------
924 ! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
925 ! (T<T0: I->G)
926 !-------------------------------------------------------------
927 if(qrs(i,k,3).gt.qcrmin) then
928 egi = exp(0.07*(-supcol))
929 acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
930 +diameter**2*rslope(i,k,3)
931 pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
932 pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
933 endif
934 endif
935 !-------------------------------------------------------------
936 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
937 ! (T<T0: C->S, and T>=T0: C->R)
938 !-------------------------------------------------------------
939 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
940 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
941 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
942 ,qci(i,k,1)/dtcld)
943 endif
944 !-------------------------------------------------------------
945 ! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
946 ! (T<T0: C->G, and T>=T0: C->R)
947 !-------------------------------------------------------------
948 if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
949 pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3) &
950 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
951 endif
952 !-------------------------------------------------------------
953 ! paacw: Accretion of cloud water by averaged snow/graupel
954 ! (T<T0: C->G or S, and T>=T0: C->R)
955 !-------------------------------------------------------------
956 if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
957 paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k))
958 endif
959 !-------------------------------------------------------------
960 ! pracs: Accretion of snow by rain [HL A11] [LFO 27]
961 ! (T<T0: S->G)
962 !-------------------------------------------------------------
963 if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
964 if(supcol.gt.0) then
965 acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1) &
966 +2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) &
967 +.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
968 pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
969 *(dens/den(i,k))*acrfac
970 pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
971 endif
972 !-------------------------------------------------------------
973 ! psacr: Accretion of rain by snow [HL A10] [LFO 28]
974 ! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
975 !-------------------------------------------------------------
976 acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2) &
977 +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
978 +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
979 psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
980 *(denr/den(i,k))*acrfac
981 psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
982 endif
983 !-------------------------------------------------------------
984 ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
985 ! (T<T0: R->G) (T>=T0: enhance melting of graupel)
986 !-------------------------------------------------------------
987 if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
988 acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3) &
989 +2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
990 +.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
991 pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
992 *acrfac
993 pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
994 endif
995 !
996 !-------------------------------------------------------------
997 ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
998 ! (S->G)
999 !-------------------------------------------------------------
1000 if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
1001 acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,3) &
1002 +2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,3) &
1003 +.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,3)
1004 if(supcol.gt.0) then
1005 egs = exp(-0.09*supcol)
1006 else
1007 egs = 1.
1008 endif
1009 pgacs(i,k) = pi**2*egs*n0s*n0sfac(i,k)*n0g*abs(vt2ave-vt2ave) &
1010 *(dens/den(i,k))*acrfac
1011 pgacs(i,k) = min(pgacs(i,k),qrs(i,k,2)/dtcld)
1012 endif
1013 if(supcol.le.0) then
1014 xlf = xlf0
1015 !-------------------------------------------------------------
1016 ! pseml: Enhanced melting of snow by accretion of water [HL A34]
1017 ! (T>=T0: S->R)
1018 !-------------------------------------------------------------
1019 if(qrs(i,k,2).gt.0.) &
1020 pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
1021 /xlf,-qrs(i,k,2)/dtcld),0.)
1022 !-------------------------------------------------------------
1023 ! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1024 ! (T>=T0: G->R)
1025 !-------------------------------------------------------------
1026 if(qrs(i,k,3).gt.0.) &
1027 pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k)) &
1028 /xlf,-qrs(i,k,3)/dtcld),0.)
1029 endif
1030 if(supcol.gt.0) then
1031 !-------------------------------------------------------------
1032 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1033 ! (T<T0: V->I or I->V)
1034 !-------------------------------------------------------------
1035 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1036 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1037 supice = satdt-prevp(i,k)
1038 if(pidep(i,k).lt.0.) then
1039 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1040 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1041 else
1042 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1043 endif
1044 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1045 endif
1046 !-------------------------------------------------------------
1047 ! psdep: deposition/sublimation rate of snow [HDC 14]
1048 ! (T<T0: V->S or S->V)
1049 !-------------------------------------------------------------
1050 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1051 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1052 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1 &
1053 *rslope2(i,k,2)+precs2*work2(i,k) &
1054 *coeres)/work1(i,k,2)
1055 supice = satdt-prevp(i,k)-pidep(i,k)
1056 if(psdep(i,k).lt.0.) then
1057 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1058 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1059 else
1060 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1061 endif
1062 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1063 ifsat = 1
1064 endif
1065 !-------------------------------------------------------------
1066 ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1067 ! (T<T0: V->G or G->V)
1068 !-------------------------------------------------------------
1069 if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
1070 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1071 pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
1072 +precg2*work2(i,k)*coeres)/work1(i,k,2)
1073 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1074 if(pgdep(i,k).lt.0.) then
1075 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1076 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1077 else
1078 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1079 endif
1080 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1081 abs(satdt)) ifsat = 1
1082 endif
1083 !-------------------------------------------------------------
1084 ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1085 ! (T<T0: V->I)
1086 !-------------------------------------------------------------
1087 if(supsat.gt.0.and.ifsat.ne.1) then
1088 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1089 xni0 = 1.e3*exp(0.1*supcol)
1090 roqi0 = 4.92e-11*xni0**1.33
1091 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
1092 /dtcld)
1093 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1094 endif
1095 !
1096 !-------------------------------------------------------------
1097 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1098 ! (T<T0: I->S)
1099 !-------------------------------------------------------------
1100 if(qci(i,k,2).gt.0.) then
1101 qimax = roqimax/den(i,k)
1102 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1103 endif
1104 !
1105 !-------------------------------------------------------------
1106 ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1107 ! (T<T0: S->G)
1108 !-------------------------------------------------------------
1109 if(qrs(i,k,2).gt.0.) then
1110 alpha2 = 1.e-3*exp(0.09*(-supcol))
1111 pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)) &
1112 ,qrs(i,k,2)/dtcld)
1113 endif
1114 endif
1115 !
1116 !-------------------------------------------------------------
1117 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1118 ! (T>=T0: S->V)
1119 !-------------------------------------------------------------
1120 if(supcol.lt.0.) then
1121 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
1122 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1123 psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1 &
1124 *rslope2(i,k,2)+precs2*work2(i,k) &
1125 *coeres)/work1(i,k,1)
1126 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1127 endif
1128 !-------------------------------------------------------------
1129 ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
1130 ! (T>=T0: G->V)
1131 !-------------------------------------------------------------
1132 if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
1133 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1134 pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
1135 +precg2*work2(i,k)*coeres)/work1(i,k,1)
1136 pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1137 endif
1138 endif
1139 enddo
1140 enddo
1141 !
1142 !
1143 !----------------------------------------------------------------
1144 ! check mass conservation of generation terms and feedback to the
1145 ! large scale
1146 !
1147 do k = kts, kte
1148 do i = its, ite
1149 !
1150 delta2=0.
1151 delta3=0.
1152 if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
1153 if(qrs(i,k,1).lt.1.e-4) delta3=1.
1154 if(t(i,k).le.t0c) then
1155 !
1156 ! cloud water
1157 !
1158 value = max(qmin,qci(i,k,1))
1159 source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1160 if (source.gt.value) then
1161 factor = value/source
1162 praut(i,k) = praut(i,k)*factor
1163 pracw(i,k) = pracw(i,k)*factor
1164 paacw(i,k) = paacw(i,k)*factor
1165 endif
1166 !
1167 ! cloud ice
1168 !
1169 value = max(qmin,qci(i,k,2))
1170 source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k) &
1171 +psaci(i,k)+pgaci(i,k))*dtcld
1172 if (source.gt.value) then
1173 factor = value/source
1174 psaut(i,k) = psaut(i,k)*factor
1175 pigen(i,k) = pigen(i,k)*factor
1176 pidep(i,k) = pidep(i,k)*factor
1177 praci(i,k) = praci(i,k)*factor
1178 psaci(i,k) = psaci(i,k)*factor
1179 pgaci(i,k) = pgaci(i,k)*factor
1180 endif
1181 !
1182 ! rain
1183 !
1184 value = max(qmin,qrs(i,k,1))
1185 source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
1186 +psacr(i,k)+pgacr(i,k))*dtcld
1187 if (source.gt.value) then
1188 factor = value/source
1189 praut(i,k) = praut(i,k)*factor
1190 prevp(i,k) = prevp(i,k)*factor
1191 pracw(i,k) = pracw(i,k)*factor
1192 piacr(i,k) = piacr(i,k)*factor
1193 psacr(i,k) = psacr(i,k)*factor
1194 pgacr(i,k) = pgacr(i,k)*factor
1195 endif
1196 !
1197 ! snow
1198 !
1199 value = max(qmin,qrs(i,k,2))
1200 source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
1201 +piacr(i,k)*delta3+praci(i,k)*delta3 &
1202 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2 &
1203 +psaci(i,k)-pgacs(i,k) )*dtcld
1204 if (source.gt.value) then
1205 factor = value/source
1206 psdep(i,k) = psdep(i,k)*factor
1207 psaut(i,k) = psaut(i,k)*factor
1208 pgaut(i,k) = pgaut(i,k)*factor
1209 paacw(i,k) = paacw(i,k)*factor
1210 piacr(i,k) = piacr(i,k)*factor
1211 praci(i,k) = praci(i,k)*factor
1212 psaci(i,k) = psaci(i,k)*factor
1213 pracs(i,k) = pracs(i,k)*factor
1214 psacr(i,k) = psacr(i,k)*factor
1215 pgacs(i,k) = pgacs(i,k)*factor
1216 endif
1217 !
1218 ! graupel
1219 !
1220 value = max(qmin,qrs(i,k,3))
1221 source = -(pgdep(i,k)+pgaut(i,k) &
1222 +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
1223 +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
1224 +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1225 if (source.gt.value) then
1226 factor = value/source
1227 pgdep(i,k) = pgdep(i,k)*factor
1228 pgaut(i,k) = pgaut(i,k)*factor
1229 piacr(i,k) = piacr(i,k)*factor
1230 praci(i,k) = praci(i,k)*factor
1231 psacr(i,k) = psacr(i,k)*factor
1232 pracs(i,k) = pracs(i,k)*factor
1233 paacw(i,k) = paacw(i,k)*factor
1234 pgaci(i,k) = pgaci(i,k)*factor
1235 pgacr(i,k) = pgacr(i,k)*factor
1236 pgacs(i,k) = pgacs(i,k)*factor
1237 endif
1238 !
1239 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k) &
1240 +pidep(i,k))
1241 ! update
1242 q(i,k) = q(i,k)+work2(i,k)*dtcld
1243 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1244 +paacw(i,k)+paacw(i,k))*dtcld,0.)
1245 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1246 +prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1247 -psacr(i,k))*dtcld,0.)
1248 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
1249 +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1250 *dtcld,0.)
1251 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1252 -pgaut(i,k)+piacr(i,k)*delta3 &
1253 +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1254 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
1255 *dtcld,0.)
1256 qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
1257 +piacr(i,k)*(1.-delta3) &
1258 +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2)&
1259 +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1260 +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1261 xlf = xls-xl(i,k)
1262 xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1263 -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
1264 +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1265 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1266 else
1267 !
1268 ! cloud water
1269 !
1270 value = max(qmin,qci(i,k,1))
1271 source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1272 if (source.gt.value) then
1273 factor = value/source
1274 praut(i,k) = praut(i,k)*factor
1275 pracw(i,k) = pracw(i,k)*factor
1276 paacw(i,k) = paacw(i,k)*factor
1277 endif
1278 !
1279 ! rain
1280 !
1281 value = max(qmin,qrs(i,k,1))
1282 source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
1283 -pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
1284 if (source.gt.value) then
1285 factor = value/source
1286 praut(i,k) = praut(i,k)*factor
1287 prevp(i,k) = prevp(i,k)*factor
1288 pracw(i,k) = pracw(i,k)*factor
1289 paacw(i,k) = paacw(i,k)*factor
1290 pseml(i,k) = pseml(i,k)*factor
1291 pgeml(i,k) = pgeml(i,k)*factor
1292 endif
1293 !
1294 ! snow
1295 !
1296 value = max(qcrmin,qrs(i,k,2))
1297 source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1298 if (source.gt.value) then
1299 factor = value/source
1300 pgacs(i,k) = pgacs(i,k)*factor
1301 psevp(i,k) = psevp(i,k)*factor
1302 pseml(i,k) = pseml(i,k)*factor
1303 endif
1304 !
1305 ! graupel
1306 !
1307 value = max(qcrmin,qrs(i,k,3))
1308 source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1309 if (source.gt.value) then
1310 factor = value/source
1311 pgacs(i,k) = pgacs(i,k)*factor
1312 pgevp(i,k) = pgevp(i,k)*factor
1313 pgeml(i,k) = pgeml(i,k)*factor
1314 endif
1315 work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1316 ! update
1317 q(i,k) = q(i,k)+work2(i,k)*dtcld
1318 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1319 +paacw(i,k)+paacw(i,k))*dtcld,0.)
1320 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1321 +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1322 -pgeml(i,k))*dtcld,0.)
1323 qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
1324 +pseml(i,k))*dtcld,0.)
1325 qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
1326 +pgeml(i,k))*dtcld,0.)
1327 xlf = xls-xl(i,k)
1328 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1329 -xlf*(pseml(i,k)+pgeml(i,k))
1330 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1331 endif
1332 enddo
1333 enddo
1334 !
1335 ! Inline expansion for fpvs
1336 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1337 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1338 hsub = xls
1339 hvap = xlv0
1340 cvap = cpv
1341 ttp=t0c+0.01
1342 dldt=cvap-cliq
1343 xa=-dldt/rv
1344 xb=xa+hvap/(rv*ttp)
1345 dldti=cvap-cice
1346 xai=-dldti/rv
1347 xbi=xai+hsub/(rv*ttp)
1348 do k = kts, kte
1349 do i = its, ite
1350 tr=ttp/t(i,k)
1351 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1352 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1353 qs(i,k,1) = max(qs(i,k,1),qmin)
1354 tr=ttp/t(i,k)
1355 if(t(i,k).lt.ttp) then
1356 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1357 else
1358 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1359 endif
1360 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1361 qs(i,k,2) = max(qs(i,k,2),qmin)
1362 enddo
1363 enddo
1364 !
1365 !----------------------------------------------------------------
1366 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1367 ! if there exists additional water vapor condensated/if
1368 ! evaporation of cloud water is not enough to remove subsaturation
1369 !
1370 do k = kts, kte
1371 do i = its, ite
1372 work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1373 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1374 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1375 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1376 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1377 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1378 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1379 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1380 enddo
1381 enddo
1382 !
1383 !
1384 !----------------------------------------------------------------
1385 ! padding for small values
1386 !
1387 do k = kts, kte
1388 do i = its, ite
1389 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1390 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1391 enddo
1392 enddo
1393 enddo ! big loops
1394 END SUBROUTINE wsm62d
1395 ! ...................................................................
1396 REAL FUNCTION rgmma(x)
1397 !-------------------------------------------------------------------
1398 IMPLICIT NONE
1399 !-------------------------------------------------------------------
1400 ! rgmma function: use infinite product form
1401 REAL :: euler
1402 PARAMETER (euler=0.577215664901532)
1403 REAL :: x, y
1404 INTEGER :: i
1405 if(x.eq.1.)then
1406 rgmma=0.
1407 else
1408 rgmma=x*exp(euler*x)
1409 do i=1,10000
1410 y=float(i)
1411 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1412 enddo
1413 rgmma=1./rgmma
1414 endif
1415 END FUNCTION rgmma
1416 !
1417 !--------------------------------------------------------------------------
1418 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1419 !--------------------------------------------------------------------------
1420 IMPLICIT NONE
1421 !--------------------------------------------------------------------------
1422 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1423 xai,xbi,ttp,tr
1424 INTEGER ice
1425 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1426 ttp=t0c+0.01
1427 dldt=cvap-cliq
1428 xa=-dldt/rv
1429 xb=xa+hvap/(rv*ttp)
1430 dldti=cvap-cice
1431 xai=-dldti/rv
1432 xbi=xai+hsub/(rv*ttp)
1433 tr=ttp/t
1434 if(t.lt.ttp.and.ice.eq.1) then
1435 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1436 else
1437 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1438 endif
1439 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1440 END FUNCTION fpvs
1441 !-------------------------------------------------------------------
1442 SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
1443 !-------------------------------------------------------------------
1444 IMPLICIT NONE
1445 !-------------------------------------------------------------------
1446 !.... constants which may not be tunable
1447 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1448 LOGICAL, INTENT(IN) :: allowed_to_read
1449 REAL :: pi
1450 !
1451 pi = 4.*atan(1.)
1452 xlv1 = cl-cpv
1453 !
1454 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1455 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1456 !
1457 bvtr1 = 1.+bvtr
1458 bvtr2 = 2.5+.5*bvtr
1459 bvtr3 = 3.+bvtr
1460 bvtr4 = 4.+bvtr
1461 bvtr6 = 6.+bvtr
1462 g1pbr = rgmma(bvtr1)
1463 g3pbr = rgmma(bvtr3)
1464 g4pbr = rgmma(bvtr4) ! 17.837825
1465 g6pbr = rgmma(bvtr6)
1466 g5pbro2 = rgmma(bvtr2) ! 1.8273
1467 pvtr = avtr*g4pbr/6.
1468 eacrr = 1.0
1469 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1470 precr1 = 2.*pi*n0r*.78
1471 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
1472 xm0 = (di0/dicon)**2
1473 xmmax = (dimax/dicon)**2
1474 roqimax = 2.08e22*dimax**8
1475 !
1476 bvts1 = 1.+bvts
1477 bvts2 = 2.5+.5*bvts
1478 bvts3 = 3.+bvts
1479 bvts4 = 4.+bvts
1480 g1pbs = rgmma(bvts1) !.8875
1481 g3pbs = rgmma(bvts3)
1482 g4pbs = rgmma(bvts4) ! 12.0786
1483 g5pbso2 = rgmma(bvts2)
1484 pvts = avts*g4pbs/6.
1485 pacrs = pi*n0s*avts*g3pbs*.25
1486 precs1 = 4.*n0s*.65
1487 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1488 pidn0r = pi*denr*n0r
1489 pidn0s = pi*dens*n0s
1490 !
1491 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1492 !
1493 bvtg1 = 1.+bvtg
1494 bvtg2 = 2.5+.5*bvtg
1495 bvtg3 = 3.+bvtg
1496 bvtg4 = 4.+bvtg
1497 g1pbg = rgmma(bvtg1)
1498 g3pbg = rgmma(bvtg3)
1499 g4pbg = rgmma(bvtg4)
1500 pacrg = pi*n0g*avtg*g3pbg*.25
1501 g5pbgo2 = rgmma(bvtg2)
1502 pvtg = avtg*g4pbg/6.
1503 precg1 = 2.*pi*n0g*.78
1504 precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
1505 pidn0g = pi*deng*n0g
1506 !
1507 rslopermax = 1./lamdarmax
1508 rslopesmax = 1./lamdasmax
1509 rslopegmax = 1./lamdagmax
1510 rsloperbmax = rslopermax ** bvtr
1511 rslopesbmax = rslopesmax ** bvts
1512 rslopegbmax = rslopegmax ** bvtg
1513 rsloper2max = rslopermax * rslopermax
1514 rslopes2max = rslopesmax * rslopesmax
1515 rslopeg2max = rslopegmax * rslopegmax
1516 rsloper3max = rsloper2max * rslopermax
1517 rslopes3max = rslopes2max * rslopesmax
1518 rslopeg3max = rslopeg2max * rslopegmax
1519 !
1520 END SUBROUTINE wsm6init
1521 END MODULE module_mp_wsm6