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