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