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