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
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, &
38 rslopermax,rslopesmax,rslopegmax, &
39 rsloperbmax,rslopesbmax,rslopegbmax, &
40 rsloper2max,rslopes2max,rslopeg2max, &
41 rsloper3max,rslopes3max,rslopeg3max
43 ! Specifies code-inlining of fpvs function in WSM32D below. JM 20040507
46 !===================================================================
48 SUBROUTINE wsm3(th, q, qci, qrs &
49 , w, den, pii, p, delz &
50 , delt,g, cpd, cpv, rd, rv, t0c &
52 , XLS, XLV0, XLF0, den0, denr &
57 , ids,ide, jds,jde, kds,kde &
58 , ims,ime, jms,jme, kms,kme &
59 , its,ite, jts,jte, kts,kte &
61 !-------------------------------------------------------------------
63 !-------------------------------------------------------------------
65 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
66 ims,ime, jms,jme, kms,kme , &
67 its,ite, jts,jte, kts,kte
68 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
74 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
80 REAL, INTENT(IN ) :: delt, &
98 REAL, DIMENSION( ims:ime , jms:jme ), &
99 INTENT(INOUT) :: rain, &
101 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
102 INTENT(INOUT) :: snow, &
106 REAL, DIMENSION( its:ite , kts:kte ) :: t
108 !-------------------------------------------------------------------
112 t(i,k)=th(i,k,j)*pii(i,k,j)
115 CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j) &
116 ,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j) &
117 ,p(ims,kms,j), delz(ims,kms,j) &
118 ,delt,g, cpd, cpv, rd, rv, t0c &
120 ,XLS, XLV0, XLF0, den0, denr &
123 ,rain(ims,j), rainncv(ims,j) &
124 ,snow(ims,j),snowncv(ims,j) &
126 ,ids,ide, jds,jde, kds,kde &
127 ,ims,ime, jms,jme, kms,kme &
128 ,its,ite, jts,jte, kts,kte &
132 th(i,k,j)=t(i,k)/pii(i,k,j)
137 !===================================================================
139 SUBROUTINE wsm32D(t, q &
140 ,qci, qrs,w, den, p, delz &
141 ,delt,g, cpd, cpv, rd, rv, t0c &
143 ,XLS, XLV0, XLF0, den0, denr &
149 ,ids,ide, jds,jde, kds,kde &
150 ,ims,ime, jms,jme, kms,kme &
151 ,its,ite, jts,jte, kts,kte &
153 !-------------------------------------------------------------------
155 !-------------------------------------------------------------------
157 ! This code is a 3-class simple ice microphyiscs scheme (WSM3) of the
158 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
159 ! number concentration is a function of temperature, and seperate assumption
160 ! is developed, in which ice crystal number concentration is a function
161 ! of ice amount. A theoretical background of the ice-microphysics and related
162 ! processes in the WSMMPs are described in Hong et al. (2004).
163 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
164 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
168 ! Developed by Song-You Hong (Yonsei Univ.), Jimy Dudhia (NCAR)
169 ! and Shu-Hua Chen (UC Davis)
172 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
175 ! History : semi-lagrangian scheme sedimentation(JH), and clean up
178 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
179 ! Dudhia (D89, 1989) J. Atmos. Sci.
180 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
181 ! Juang and Hong (JH, 2010) Mon. Wea. Rev.
183 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
184 ims,ime, jms,jme, kms,kme, &
185 its,ite, jts,jte, kts,kte, &
187 REAL, DIMENSION( its:ite , kts:kte ), &
190 REAL, DIMENSION( ims:ime , kms:kme ), &
195 REAL, DIMENSION( ims:ime , kms:kme ), &
200 REAL, INTENT(IN ) :: delt, &
218 REAL, DIMENSION( ims:ime ), &
219 INTENT(INOUT) :: rain, &
221 REAL, DIMENSION( ims:ime ), OPTIONAL, &
222 INTENT(INOUT) :: snow, &
226 REAL, DIMENSION( its:ite , kts:kte ) :: &
237 REAL, DIMENSION( its:ite , kts:kte ) :: &
244 REAL, DIMENSION( its:ite , kts:kte ) :: &
260 REAL, DIMENSION( its:ite ) :: delqrs,&
262 INTEGER, DIMENSION( its:ite ) :: kwork1,&
264 INTEGER, DIMENSION( its:ite ) :: mstep, &
266 LOGICAL, DIMENSION( its:ite ) :: flgcld
268 cpmcal, xlcal, diffus, &
269 viscos, xka, venfac, conden, diffac, &
270 x, y, z, a, b, c, d, e, &
271 fallsum, fallsum_qsi, vt2i,vt2s,acrfac, &
272 qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
273 snomlt, hold, holdrs, facqci, supcol, coeres, &
274 supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
275 qimax, diameter, xni0, roqi0, supice,holdc, holdci
276 INTEGER :: i, j, k, mstepmax, &
277 iprt, latd, lond, loop, loops, ifsat, kk, n, idim, kdim
278 ! Temporaries used for inlining fpvs function
279 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
280 ! variables for optimization
281 REAL, DIMENSION( its:ite ) :: tvec1
283 !=================================================================
284 ! compute internal functions
286 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
287 xlcal(x) = xlv0-xlv1*(x-t0c)
288 !----------------------------------------------------------------
289 ! diffus: diffusion coefficient of the water vapor
290 ! viscos: kinematic viscosity(m2s-1)
291 ! Optimizatin : A**B => exp(log(A)*(B))
293 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
294 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
295 xka(x,y) = 1.414e3*viscos(x,y)*y
296 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
297 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
298 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
299 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
304 !----------------------------------------------------------------
305 ! paddint 0 for negative values generated by dynamics
309 qci(i,k) = max(qci(i,k),0.0)
310 qrs(i,k) = max(qrs(i,k),0.0)
314 !----------------------------------------------------------------
315 ! latent heat for phase changes and heat capacity. neglect the
316 ! changes during microphysical process calculation
321 cpm(i,k) = cpmcal(q(i,k))
322 xl(i,k) = xlcal(t(i,k))
327 delz_tmp(i,k) = delz(i,k)
328 den_tmp(i,k) = den(i,k)
332 !----------------------------------------------------------------
333 ! compute the minor time steps.
335 loops = max(nint(delt/dtcldcr),1)
337 if(delt.le.dtcldcr) dtcld = delt
341 !----------------------------------------------------------------
342 ! initialize the large scale variables
349 CALL VREC( tvec1(its), den(its,k), ite-its+1)
351 tvec1(i) = tvec1(i)*den0
353 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
356 ! Inline expansion for fpvs
357 ! qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
358 ! qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
368 xbi=xai+hsub/(rv*ttp)
372 if(t(i,k).lt.ttp) then
373 qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
375 qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
377 qs0(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
378 qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
379 qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
380 qs(i,k) = max(qs(i,k),qmin)
381 rh(i,k) = max(q(i,k) / qs(i,k),qmin)
385 !----------------------------------------------------------------
386 ! initialize the variables for microphysical physics
404 !-------------------------------------------------------------
405 ! Ni: ice crystal number concentraiton [HDC 5c]
406 !-------------------------------------------------------------
409 xni(i,k) = min(max(5.38e7 &
410 *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
414 !----------------------------------------------------------------
415 ! compute the fallout term:
416 ! first, vertical terminal velosity for minor loops
417 !---------------------------------------------------------------
420 qrs_tmp(i,k) = qrs(i,k)
423 call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
424 work1,its,ite,kts,kte)
427 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
431 denqrs(i,k) = den(i,k)*qrs(i,k)
434 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1,denqrs, &
438 qrs(i,k) = max(denqrs(i,k)/den(i,k),0.)
439 fall(i,k) = denqrs(i,k)*work1(i,k)/delz(i,k)
443 fall(i,1) = delqrs(i)/delz(i,1)/dtcld
445 !---------------------------------------------------------------
446 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
447 !---------------------------------------------------------------
450 if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
451 xmi = den(i,k)*qci(i,k)/xni(i,k)
452 diameter = max(dicon * sqrt(xmi), 1.e-25)
453 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
460 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
464 denqci(i,k) = den(i,k)*qci(i,k)
467 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
471 qci(i,k) = max(denqci(i,k)/den(i,k),0.)
475 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
478 !----------------------------------------------------------------
479 ! compute the freezing/melting term. [D89 B16-B17]
480 ! freezing occurs one layer above the melting level
488 if(t(i,k).ge.t0c) then
497 if(mstep(i).ne.0) then
498 if (w(i,mstep(i)).gt.0.) then
499 kwork1(i) = mstep(i) + 1
508 qrsci = qrs(i,k) + qci(i,k)
509 if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
510 frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld), &
512 snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld), &
515 t(i,k) = t(i,k) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld
517 t(i,k) = t(i,k) - xlf0/cpm(i,k)*frzmlt*dtcld
518 t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)*snomlt*dtcld
524 !----------------------------------------------------------------
525 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
530 if((t0c-t(i,1)).gt.0) then
531 fallsum = fallsum+fallc(i,1)
532 fallsum_qsi = fall(i,1)+fallc(i,1)
535 if(fallsum.gt.0.) then
536 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
537 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
539 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
541 if(fallsum_qsi.gt.0.) then
542 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
543 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
547 if(fallsum.gt.0.) sr(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
551 !----------------------------------------------------------------
552 ! update the slope parameters for microphysics computation
556 qrs_tmp(i,k) = qrs(i,k)
559 call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
560 work1,its,ite,kts,kte)
562 ! work1: the thermodynamic term in the denominator associated with
563 ! heat conduction and vapor diffusion
564 ! work2: parameter associated with the ventilation effects(y93)
568 if(t(i,k).ge.t0c) then
569 work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
571 work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
573 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
579 supsat = max(q(i,k),qmin)-qs(i,k)
581 if(t(i,k).ge.t0c) then
583 !===============================================================
585 ! warm rain processes
587 ! - follows the processes in RH83 and LFO except for autoconcersion
589 !===============================================================
590 !---------------------------------------------------------------
591 ! praut: auto conversion rate from cloud to rain [HDC 16]
593 !---------------------------------------------------------------
594 if(qci(i,k).gt.qc0) then
595 ! paut(i,k) = qck1*qci(i,k)**(7./3.)
596 paut(i,k) = qck1*exp(log(qci(i,k))*((7./3.)))
597 paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
599 !---------------------------------------------------------------
600 ! pracw: accretion of cloud water by rain [HL A40] [D89 B15]
602 !---------------------------------------------------------------
603 if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
604 pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k) &
605 *qci(i,k)*denfac(i,k),qci(i,k)/dtcld)
607 !---------------------------------------------------------------
608 ! prevp: evaporation/condensation rate of rain [HDC 14]
610 !---------------------------------------------------------------
611 if(qrs(i,k).gt.0.) then
612 coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
613 pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k) &
614 +precr2*work2(i,k)*coeres)/work1(i,k)
615 if(pres(i,k).lt.0.) then
616 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
617 pres(i,k) = max(pres(i,k),satdt/2)
619 pres(i,k) = min(pres(i,k),satdt/2)
624 !===============================================================
626 ! cold rain processes
628 ! - follows the revised ice microphysics processes in HDC
629 ! - the processes same as in RH83 and LFO behave
630 ! following ice crystal hapits defined in HDC, inclduing
631 ! intercept parameter for snow (n0s), ice crystal number
632 ! concentration (ni), ice nuclei number concentration
633 ! (n0i), ice diameter (d)
635 !===============================================================
638 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
640 !-------------------------------------------------------------
641 ! Ni: ice crystal number concentraiton [HDC 5c]
642 !-------------------------------------------------------------
643 xni(i,k) = min(max(5.38e7 &
644 *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
645 eacrs = exp(0.07*(-supcol))
646 if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
647 xmi = den(i,k)*qci(i,k)/xni(i,k)
648 diameter = min(dicon * sqrt(xmi),dimax)
649 vt2i = 1.49e4*diameter**1.31
650 vt2s = pvts*rslopeb(i,k)*denfac(i,k)
651 !-------------------------------------------------------------
652 ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
654 !-------------------------------------------------------------
655 acrfac = 2.*rslope3(i,k)+2.*diameter*rslope2(i,k) &
656 +diameter**2*rslope(i,k)
657 pacr(i,k) = min(pi*qci(i,k)*eacrs*n0s*n0sfac(i,k) &
658 *abs(vt2s-vt2i)*acrfac/4.,qci(i,k)/dtcld)
660 !-------------------------------------------------------------
661 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
662 ! (T<T0: V->I or I->V)
663 !-------------------------------------------------------------
664 if(qci(i,k).gt.0.) then
665 xmi = den(i,k)*qci(i,k)/xni(i,k)
666 diameter = dicon * sqrt(xmi)
667 pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.)/work1(i,k)
668 if(pisd(i,k).lt.0.) then
669 pisd(i,k) = max(pisd(i,k),satdt/2)
670 pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
672 pisd(i,k) = min(pisd(i,k),satdt/2)
674 if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
676 !-------------------------------------------------------------
677 ! psdep: deposition/sublimation rate of snow [HDC 14]
679 !-------------------------------------------------------------
680 if(qrs(i,k).gt.0..and.ifsat.ne.1) then
681 coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
682 pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k) &
683 +precs2*work2(i,k)*coeres)/work1(i,k)
684 supice = satdt-pisd(i,k)
685 if(pres(i,k).lt.0.) then
686 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
687 pres(i,k) = max(max(pres(i,k),satdt/2),supice)
689 pres(i,k) = min(min(pres(i,k),satdt/2),supice)
691 if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
693 !-------------------------------------------------------------
694 ! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
696 !-------------------------------------------------------------
697 if(supsat.gt.0.and.ifsat.ne.1) then
698 supice = satdt-pisd(i,k)-pres(i,k)
699 xni0 = 1.e3*exp(0.1*supcol)
700 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
701 pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld)
702 pgen(i,k) = min(min(pgen(i,k),satdt),supice)
704 !-------------------------------------------------------------
705 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
707 !-------------------------------------------------------------
708 if(qci(i,k).gt.0.) then
709 qimax = roqimax/den(i,k)
710 paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld)
716 !----------------------------------------------------------------
717 ! check mass conservation of generation terms and feedback to the
722 qciik = max(qmin,qci(i,k))
723 delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
724 if(delqci.ge.qciik) then
725 facqci = qciik/delqci
726 paut(i,k) = paut(i,k)*facqci
727 pacr(i,k) = pacr(i,k)*facqci
728 pgen(i,k) = pgen(i,k)*facqci
729 pisd(i,k) = pisd(i,k)*facqci
731 qik = max(qmin,q(i,k))
732 delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
735 pres(i,k) = pres(i,k)*facq
736 pgen(i,k) = pgen(i,k)*facq
737 pisd(i,k) = pisd(i,k)*facq
739 work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
740 q(i,k) = q(i,k)+work2(i,k)*dtcld
741 qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k)) &
743 qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k)+pres(i,k))*dtcld,0.)
744 if(t(i,k).lt.t0c) then
745 t(i,k) = t(i,k)-xls*work2(i,k)/cpm(i,k)*dtcld
747 t(i,k) = t(i,k)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld
761 xbi=xai+hsub/(rv*ttp)
765 qs(i,k)=psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
766 qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
767 qs(i,k) = max(qs(i,k),qmin)
768 denfac(i,k) = sqrt(den0/den(i,k))
772 !----------------------------------------------------------------
773 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
774 ! if there exists additional water vapor condensated/if
775 ! evaporation of cloud water is not enough to remove subsaturation
779 work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
780 work2(i,k) = qci(i,k)+work1(i,k)
781 pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
782 if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c) &
783 pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
784 q(i,k) = q(i,k)-pcon(i,k)*dtcld
785 qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
786 t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
790 !----------------------------------------------------------------
791 ! padding for small values
795 if(qci(i,k).le.qmin) qci(i,k) = 0.0
796 if(qrs(i,k).le.qcrmin) qrs(i,k) = 0.0
801 END SUBROUTINE wsm32D
802 ! ...................................................................
803 REAL FUNCTION rgmma(x)
804 !-------------------------------------------------------------------
806 !-------------------------------------------------------------------
807 ! rgmma function: use infinite product form
809 PARAMETER (euler=0.577215664901532)
818 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
824 !--------------------------------------------------------------------------
825 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
826 !--------------------------------------------------------------------------
828 !--------------------------------------------------------------------------
829 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
832 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
839 xbi=xai+hsub/(rv*ttp)
841 if(t.lt.ttp.and.ice.eq.1) then
842 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
844 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
846 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
848 !-------------------------------------------------------------------
849 SUBROUTINE wsm3init(den0,denr,dens,cl,cpv,allowed_to_read)
850 !-------------------------------------------------------------------
852 !-------------------------------------------------------------------
853 !.... constants which may not be tunable
854 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
855 LOGICAL, INTENT(IN) :: allowed_to_read
860 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
861 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
869 g4pbr = rgmma(bvtr4) ! 17.837825
870 g5pbro2 = rgmma(bvtr2) ! 1.8273
873 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
874 precr1 = 2.*pi*n0r*.78
875 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
876 xmmax = (dimax/dicon)**2
877 roqimax = 2.08e22*dimax**8
883 g1pbs = rgmma(bvts1) !.8875
885 g4pbs = rgmma(bvts4) ! 12.0786
886 g5pbso2 = rgmma(bvts2)
888 pacrs = pi*n0s*avts*g3pbs*.25
890 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
894 rslopermax = 1./lamdarmax
895 rslopesmax = 1./lamdasmax
896 rsloperbmax = rslopermax ** bvtr
897 rslopesbmax = rslopesmax ** bvts
898 rsloper2max = rslopermax * rslopermax
899 rslopes2max = rslopesmax * rslopesmax
900 rsloper3max = rsloper2max * rslopermax
901 rslopes3max = rslopes2max * rslopesmax
903 END SUBROUTINE wsm3init
905 subroutine slope_wsm3(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,vt,its,ite,kts,kte)
907 INTEGER :: its,ite, jts,jte, kts,kte
908 REAL, DIMENSION( its:ite , kts:kte ) :: &
918 REAL, PARAMETER :: t0c = 273.15
919 REAL, DIMENSION( its:ite , kts:kte ) :: &
921 REAL :: lamdar,lamdas,x, y, z, supcol, pvt
923 !----------------------------------------------------------------
924 ! size distributions: (x=mixing ratio, y=air density):
925 ! valid for mixing ratio > 1.e-9 kg/kg.
927 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
928 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
932 if(t(i,k).ge.t0c) then
934 if(qrs(i,k).le.qcrmin)then
935 rslope(i,k) = rslopermax
936 rslopeb(i,k) = rsloperbmax
937 rslope2(i,k) = rsloper2max
938 rslope3(i,k) = rsloper3max
940 rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
941 rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
942 rslope2(i,k) = rslope(i,k)*rslope(i,k)
943 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
947 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
949 if(qrs(i,k).le.qcrmin)then
950 rslope(i,k) = rslopesmax
951 rslopeb(i,k) = rslopesbmax
952 rslope2(i,k) = rslopes2max
953 rslope3(i,k) = rslopes3max
955 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
956 rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
957 rslope2(i,k) = rslope(i,k)*rslope(i,k)
958 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
961 vt(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
962 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
965 END subroutine slope_wsm3
966 !-------------------------------------------------------------------
967 SUBROUTINE nislfv_rain_pcm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
968 !-------------------------------------------------------------------
970 ! for non-iteration semi-Lagrangain forward advection for cloud
971 ! with mass conservation and positive definite advection
972 ! 2nd order interpolation with monotonic piecewise linear method
973 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
975 ! dzl depth of model layer in meter
976 ! wwl terminal velocity at model layer m/s
977 ! rql cloud density*mixing ration
978 ! precip precipitation
980 ! id kind of precip: 0 test case; 1 raindrop
981 ! iter how many time to guess mean terminal velocity: 0 pure forward.
982 ! 0 : use departure wind for advection
983 ! 1 : use mean wind for advection
984 ! > 1 : use mean wind after iter-1 iterations
986 ! author: hann-ming henry juang <henry.juang@noaa.gov>
987 ! implemented by song-you hong
992 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
993 real denl(im,km),denfacl(im,km),tkl(im,km)
995 integer i,k,n,m,kk,kb,kt,iter
996 real tl,tl2,qql,dql,qqd
998 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
999 real zsumt,qsumt,zsumb,qsumb
1000 real allold, allnew, zz, dzamin, cflmax, decfl
1001 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1002 real den(km), denfac(km), tk(km)
1003 real wi(km+1), zi(km+1), za(km+1)
1004 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1005 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1010 ! -----------------------------------
1015 denfac(:) = denfacl(i,:)
1017 ! skip for no precipitation for all layers
1020 allold = allold + qq(k)
1022 if(allold.le.0.0) then
1026 ! compute interface values
1029 zi(k+1) = zi(k)+dz(k)
1032 ! save departure wind
1036 ! pcm is 1st order, we should use 2nd order wi
1037 ! 2nd order interpolation to get wi
1040 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1044 ! terminate of top of raingroup
1046 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1052 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1053 if( decfl .gt. con1 ) then
1054 wi(k) = wi(k+1) - con1*dz(k)/dt
1057 ! compute arrival point
1059 za(k) = zi(k) - wi(k)*dt
1063 dza(k) = za(k+1)-za(k)
1065 dza(km+1) = zi(km+1) - za(km+1)
1067 ! computer deformation at arrival point
1069 qa(k) = qq(k)*dz(k)/dza(k)
1070 qr(k) = qa(k)/den(k)
1073 ! call maxmin(km,1,qa,' arrival points ')
1075 ! compute arrival terminal velocity, and estimate mean terminal velocity
1076 ! then back to use mean terminal velocity
1077 if( n.le.iter ) then
1078 call slope_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1079 if( n.eq.2 ) wa(1:km) = 0.5*(wa(1:km)+was(1:km))
1082 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1084 ! mean wind is average of departure and new arrival winds
1085 ww(k) = 0.5* ( wd(k)+wa(k) )
1093 ! interpolation to regular point
1101 if( zi(k).ge.za(km+1) ) then
1104 find_kb : do kk=kb,km
1105 if( zi(k).le.za(kk+1) ) then
1112 find_kt : do kk=kt,km
1113 if( zi(k+1).le.za(kk) ) then
1120 ! compute q with piecewise constant method
1121 if( kt-kb.eq.1 ) then
1123 else if( kt-kb.ge.2 ) then
1124 zsumb = za(kb+1)-zi(k)
1125 qsumb = qa(kb) * zsumb
1126 zsumt = zi(k+1)-za(kt-1)
1127 qsumt = qa(kt-1) * zsumt
1130 if( kt-kb.ge.3 ) then
1132 qsum = qsum + qa(m) * dza(m)
1133 zsum = zsum + dza(m)
1136 qn(k) = (qsumb+qsum+qsumt)/(zsumb+zsum+zsumt)
1144 sum_precip: do k=1,km
1145 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1146 precip(i) = precip(i) + qa(k)*dza(k)
1148 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1149 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1155 ! replace the new values
1158 ! ----------------------------------
1161 END SUBROUTINE nislfv_rain_pcm
1162 !-------------------------------------------------------------------
1163 SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1164 !-------------------------------------------------------------------
1166 ! for non-iteration semi-Lagrangain forward advection for cloud
1167 ! with mass conservation and positive definite advection
1168 ! 2nd order interpolation with monotonic piecewise linear method
1169 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1171 ! dzl depth of model layer in meter
1172 ! wwl terminal velocity at model layer m/s
1173 ! rql cloud density*mixing ration
1174 ! precip precipitation
1176 ! id kind of precip: 0 test case; 1 raindrop
1177 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1178 ! 0 : use departure wind for advection
1179 ! 1 : use mean wind for advection
1180 ! > 1 : use mean wind after iter-1 iterations
1182 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1183 ! implemented by song-you hong
1188 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1189 real denl(im,km),denfacl(im,km),tkl(im,km)
1191 integer i,k,n,m,kk,kb,kt,iter
1192 real tl,tl2,qql,dql,qqd
1194 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1195 real allold, allnew, zz, dzamin, cflmax, decfl
1196 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1197 real den(km), denfac(km), tk(km)
1198 real wi(km+1), zi(km+1), za(km+1)
1199 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1200 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1205 ! -----------------------------------
1210 denfac(:) = denfacl(i,:)
1212 ! skip for no precipitation for all layers
1215 allold = allold + qq(k)
1217 if(allold.le.0.0) then
1221 ! compute interface values
1224 zi(k+1) = zi(k)+dz(k)
1227 ! save departure wind
1231 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1232 ! 2nd order interpolation to get wi
1236 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1238 ! 3rd order interpolation to get wi
1242 wi(2) = 0.5*(ww(2)+ww(1))
1244 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1246 wi(km) = 0.5*(ww(km)+ww(km-1))
1249 ! terminate of top of raingroup
1251 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1257 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1258 if( decfl .gt. con1 ) then
1259 wi(k) = wi(k+1) - con1*dz(k)/dt
1262 ! compute arrival point
1264 za(k) = zi(k) - wi(k)*dt
1268 dza(k) = za(k+1)-za(k)
1270 dza(km+1) = zi(km+1) - za(km+1)
1272 ! computer deformation at arrival point
1274 qa(k) = qq(k)*dz(k)/dza(k)
1275 qr(k) = qa(k)/den(k)
1278 ! call maxmin(km,1,qa,' arrival points ')
1280 ! compute arrival terminal velocity, and estimate mean terminal velocity
1281 ! then back to use mean terminal velocity
1282 if( n.le.iter ) then
1283 call slope_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1284 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1287 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1289 ! mean wind is average of departure and new arrival winds
1290 ww(k) = 0.5* ( wd(k)+wa(k) )
1297 ! estimate values at arrival cell interface with monotone
1299 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1300 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1301 if( dip*dim.le.0.0 ) then
1305 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1306 qmi(k)=2.0*qa(k)-qpi(k)
1307 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1318 ! interpolation to regular point
1326 if( zi(k).ge.za(km+1) ) then
1329 find_kb : do kk=kb,km
1330 if( zi(k).le.za(kk+1) ) then
1337 find_kt : do kk=kt,km
1338 if( zi(k+1).le.za(kk) ) then
1346 ! compute q with piecewise constant method
1348 tl=(zi(k)-za(kb))/dza(kb)
1349 th=(zi(k+1)-za(kb))/dza(kb)
1352 qqd=0.5*(qpi(kb)-qmi(kb))
1353 qqh=qqd*th2+qmi(kb)*th
1354 qql=qqd*tl2+qmi(kb)*tl
1355 qn(k) = (qqh-qql)/(th-tl)
1356 else if( kt.gt.kb ) then
1357 tl=(zi(k)-za(kb))/dza(kb)
1359 qqd=0.5*(qpi(kb)-qmi(kb))
1360 qql=qqd*tl2+qmi(kb)*tl
1362 zsum = (1.-tl)*dza(kb)
1364 if( kt-kb.gt.1 ) then
1366 zsum = zsum + dza(m)
1367 qsum = qsum + qa(m) * dza(m)
1370 th=(zi(k+1)-za(kt))/dza(kt)
1372 qqd=0.5*(qpi(kt)-qmi(kt))
1373 dqh=qqd*th2+qmi(kt)*th
1374 zsum = zsum + th*dza(kt)
1375 qsum = qsum + dqh*dza(kt)
1384 sum_precip: do k=1,km
1385 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1386 precip(i) = precip(i) + qa(k)*dza(k)
1388 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1389 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1395 ! replace the new values
1398 ! ----------------------------------
1401 END SUBROUTINE nislfv_rain_plm
1403 END MODULE module_mp_wsm3