9 !Including inline expansion statistical function
13 REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
14 REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
15 REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
16 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
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 ! a constant for terminal velocity of snow
22 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
23 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
24 REAL, PARAMETER, PRIVATE :: lamdacmax = 1.e10 ! limited maximum value for slope parameter of cloud water
25 REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e8 ! limited maximum value for slope parameter of rain
26 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
27 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
28 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
29 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
30 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
31 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
32 REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
33 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
34 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
35 REAL, PARAMETER, PRIVATE :: ncmin = 1.e1 ! minimum value for Nc
36 REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2 ! minimum value for Nr
37 REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
39 REAL, PARAMETER, PRIVATE :: satmax = 1.0048 ! maximum saturation value for CCN activation
40 ! 1.008 for maritime air mass /1.0048 for conti
41 REAL, PARAMETER, PRIVATE :: actk = 0.6 ! parameter for the CCN activation
42 REAL, PARAMETER, PRIVATE :: actr = 1.5 ! radius of activated CCN drops
43 REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3 ! Long's collection kernel coefficient
44 REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15 ! Long's collection kernel coefficient
45 REAL, PARAMETER, PRIVATE :: di100 = 1.e-4 ! parameter related with accretion and collection of cloud drops
46 REAL, PARAMETER, PRIVATE :: di600 = 6.e-4 ! parameter related with accretion and collection of cloud drops
47 REAL, PARAMETER, PRIVATE :: di2000 = 20.e-4 ! parameter related with accretion and collection of cloud drops
48 REAL, PARAMETER, PRIVATE :: di82 = 82.e-6 ! dimater related with raindrops evaporation
49 REAL, PARAMETER, PRIVATE :: di15 = 15.e-6 ! auto conversion takes place beyond this diameter
51 qc0, qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4, &
52 bvtr5,bvtr7,bvtr2o5,bvtr3o5,g1pbr,g2pbr, &
53 g3pbr,g4pbr,g5pbr,g7pbr,g5pbro2,g7pbro2, &
54 pvtr,pvtrn,eacrr,pacrr, pi, &
55 precr1,precr2,xmmax,roqimax,bvts1, &
56 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
57 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
58 pidn0s,pidnr,xlv1,pacrc, &
59 rslopecmax,rslopec2max,rslopec3max, &
60 rslopermax,rslopesmax,rslopegmax, &
61 rsloperbmax,rslopesbmax,rslopegbmax, &
62 rsloper2max,rslopes2max,rslopeg2max, &
63 rsloper3max,rslopes3max,rslopeg3max
65 ! Specifies code-inlining of fpvs function in WDM52D below. JM 20040507
68 !===================================================================
70 SUBROUTINE wdm5(th, q, qc, qr, qi, qs &
73 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
75 ,XLS, XLV0, XLF0, den0, denr &
81 ,ids,ide, jds,jde, kds,kde &
82 ,ims,ime, jms,jme, kms,kme &
83 ,its,ite, jts,jte, kts,kte &
85 !-------------------------------------------------------------------
87 !-------------------------------------------------------------------
89 ! This code is a WRF double-moment 5-class mixed ice
90 ! microphyiscs scheme (WDM5). The WDM microphysics scheme predicts
91 ! number concentrations for warm rain species including clouds and
92 ! rain. cloud condensation nuclei (CCN) is also predicted.
93 ! The cold rain species including ice, snow, graupel follow the
94 ! WRF single-moment 5-class microphysics (WSM5)
95 ! in which theoretical background for WSM ice phase microphysics is
96 ! based on Hong et al. (2004).
97 ! The WDM scheme is described in Lim and Hong (2009).
98 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
102 ! Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
104 ! Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
106 ! Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
107 ! Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
108 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
109 ! Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
110 ! Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
111 ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
113 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
114 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
115 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
117 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
118 ims,ime, jms,jme, kms,kme , &
119 its,ite, jts,jte, kts,kte
120 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
131 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
137 REAL, INTENT(IN ) :: delt, &
156 INTEGER, INTENT(IN ) :: itimestep
157 REAL, DIMENSION( ims:ime , jms:jme ), &
158 INTENT(INOUT) :: rain, &
161 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
162 INTENT(INOUT) :: snow, &
165 REAL, DIMENSION( its:ite , kts:kte ) :: t
166 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
167 REAL, DIMENSION( its:ite , kts:kte, 3 ) :: ncr
168 CHARACTER*256 :: emess
171 !-------------------------------------------------------------------
173 IF (itimestep .eq. 1) THEN
186 t(i,k)=th(i,k,j)*pii(i,k,j)
187 qci(i,k,1) = qc(i,k,j)
188 qci(i,k,2) = qi(i,k,j)
189 qrs(i,k,1) = qr(i,k,j)
190 qrs(i,k,2) = qs(i,k,j)
191 ncr(i,k,1) = nn(i,k,j)
192 ncr(i,k,2) = nc(i,k,j)
193 ncr(i,k,3) = nr(i,k,j)
196 ! Sending array starting locations of optional variables may cause
197 ! troubles, so we explicitly change the call.
198 CALL wdm52D(t, q(ims,kms,j), qci, qrs, ncr &
200 ,p(ims,kms,j), delz(ims,kms,j) &
201 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
203 ,XLS, XLV0, XLF0, den0, denr &
206 ,rain(ims,j),rainncv(ims,j) &
208 ,ids,ide, jds,jde, kds,kde &
209 ,ims,ime, jms,jme, kms,kme &
210 ,its,ite, jts,jte, kts,kte &
211 ,snow(ims,j),snowncv(ims,j) &
215 th(i,k,j)=t(i,k)/pii(i,k,j)
216 qc(i,k,j) = qci(i,k,1)
217 qi(i,k,j) = qci(i,k,2)
218 qr(i,k,j) = qrs(i,k,1)
219 qs(i,k,j) = qrs(i,k,2)
220 nn(i,k,j) = ncr(i,k,1)
221 nc(i,k,j) = ncr(i,k,2)
222 nr(i,k,j) = ncr(i,k,3)
227 CALL get_wsm5_gpu_levels ( mkx_test )
228 IF ( mkx_test .LT. kte ) THEN
229 WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ', &
231 CALL wrf_error_fatal(emess)
234 th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte) &
235 ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte) &
236 ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte) &
237 ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte) &
238 ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte) &
240 ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte) &
241 ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte) &
242 ,sr(its:ite,jts:jte) &
243 ,its, ite, jts, jte, kts, kte &
244 ,its, ite, jts, jte, kts, kte &
245 ,its, ite, jts, jte, kts, kte &
249 !===================================================================
251 SUBROUTINE wdm52D(t, q, qci, qrs, ncr, den, p, delz &
252 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
254 ,XLS, XLV0, XLF0, den0, denr &
259 ,ids,ide, jds,jde, kds,kde &
260 ,ims,ime, jms,jme, kms,kme &
261 ,its,ite, jts,jte, kts,kte &
264 !-------------------------------------------------------------------
266 !-------------------------------------------------------------------
267 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
268 ims,ime, jms,jme, kms,kme , &
269 its,ite, jts,jte, kts,kte, &
271 REAL, DIMENSION( its:ite , kts:kte ), &
274 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
278 REAL, DIMENSION( its:ite , kts:kte, 3 ), &
281 REAL, DIMENSION( ims:ime , kms:kme ), &
284 REAL, DIMENSION( ims:ime , kms:kme ), &
289 REAL, INTENT(IN ) :: delt, &
308 REAL, DIMENSION( ims:ime ), &
309 INTENT(INOUT) :: rain, &
312 REAL, DIMENSION( ims:ime ), OPTIONAL, &
313 INTENT(INOUT) :: snow, &
316 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
317 rh, qs, rslope, rslope2, rslope3, rslopeb, &
318 falk, fall, work1, qrs_tmp
319 REAL, DIMENSION( its:ite , kts:kte ) :: &
320 rslopec, rslopec2,rslopec3
321 REAL, DIMENSION( its:ite , kts:kte, 2) :: &
323 REAL, DIMENSION( its:ite , kts:kte ) :: &
325 REAL, DIMENSION( its:ite , kts:kte ) :: &
327 REAL, DIMENSION( its:ite , kts:kte ) :: &
328 den_tmp, delz_tmp, ncr_tmp
329 REAL, DIMENSION( its:ite , kts:kte ) :: &
330 falkc, work1c, work2c, fallc
331 REAL, DIMENSION( its:ite , kts:kte ) :: &
332 pcact, praut, psaut, prevp, psdep, pracw, psaci, psacw, &
333 pigen, pidep, pcond, prevp_s, &
334 xl, cpm, work2, psmlt, psevp, denfac, xni, &
335 n0sfac, denqrs2, denqci
336 REAL, DIMENSION( its:ite ) :: &
338 REAL, DIMENSION( its:ite , kts:kte ) :: &
339 nraut, nracw, nrevp, ncevp, nccol, nrcol, &
343 #define WSM_NO_CONDITIONAL_IN_VECTOR
344 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
345 REAL, DIMENSION(its:ite) :: xal, xbl
347 ! variables for optimization
348 REAL, DIMENSION( its:ite ) :: tvec1
349 INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
350 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
351 REAL, DIMENSION(its:ite) :: rmstep
352 REAL dtcldden, rdelz, rdtcld
353 LOGICAL, DIMENSION( its:ite ) :: flgcld
355 cpmcal, xlcal, lamdac, diffus, &
356 viscos, xka, venfac, conden, diffac, &
357 x, y, z, a, b, c, d, e, &
358 ndt, qdt, holdrr, holdrs, supcol, supcolt, pvt, &
359 coeres, supsat, dtcld, xmi, eacrs, satdt, &
360 vt2i,vt2s,acrfac, coecol, &
362 taucon, lencon, lenconcr, &
363 qimax, diameter, xni0, roqi0, &
364 fallsum, fallsum_qsi, xlwork2, factor, source, &
365 value, xlf, pfrzdtc, pfrzdtr, supice
367 REAL :: holdc, holdci
368 INTEGER :: i, j, k, mstepmax, &
369 iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
370 ! Temporaries used for inlining fpvs function
371 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
374 !=================================================================
375 ! compute internal functions
377 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
378 xlcal(x) = xlv0-xlv1*(x-t0c)
379 !----------------------------------------------------------------
380 ! size distributions: (x=mixing ratio, y=air density):
381 ! valid for mixing ratio > 1.e-9 kg/kg.
383 ! Optimizatin : A**B => exp(log(A)*(B))
384 lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
386 !----------------------------------------------------------------
387 ! diffus: diffusion coefficient of the water vapor
388 ! viscos: kinematic viscosity(m2s-1)
389 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y
390 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y
391 ! xka(x,y) = 1.414e3*viscos(x,y)*y
392 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
393 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
394 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
395 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
401 !----------------------------------------------------------------
402 ! paddint 0 for negative values generated by dynamics
406 qci(i,k,1) = max(qci(i,k,1),0.0)
407 qrs(i,k,1) = max(qrs(i,k,1),0.0)
408 qci(i,k,2) = max(qci(i,k,2),0.0)
409 qrs(i,k,2) = max(qrs(i,k,2),0.0)
410 ncr(i,k,1) = max(ncr(i,k,1),0.)
411 ncr(i,k,2) = max(ncr(i,k,2),0.)
412 ncr(i,k,3) = max(ncr(i,k,3),0.)
416 ! latent heat for phase changes and heat capacity. neglect the
417 ! changes during microphysical process calculation
422 cpm(i,k) = cpmcal(q(i,k))
423 xl(i,k) = xlcal(t(i,k))
424 delz_tmp(i,k) = delz(i,k)
425 den_tmp(i,k) = den(i,k)
429 !----------------------------------------------------------------
430 ! compute the minor time steps.
432 loops = max(nint(delt/dtcldcr),1)
434 if(delt.le.dtcldcr) dtcld = delt
438 !----------------------------------------------------------------
439 ! initialize the large scale variables
449 ! denfac(i,k) = sqrt(den0/den(i,k))
453 CALL VREC( tvec1(its), den(its,k), ite-its+1)
455 tvec1(i) = tvec1(i)*den0
457 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
460 ! Inline expansion for fpvs
461 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
462 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
472 xbi=xai+hsub/(rv*ttp)
473 ! this is for compilers where the conditional inhibits vectorization
474 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
477 if(t(i,k).lt.ttp) then
488 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
489 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
490 qs(i,k,1) = max(qs(i,k,1),qmin)
491 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
492 qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
493 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
494 qs(i,k,2) = max(qs(i,k,2),qmin)
495 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
503 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
504 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
505 qs(i,k,1) = max(qs(i,k,1),qmin)
506 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
507 if(t(i,k).lt.ttp) then
508 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
510 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
512 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
513 qs(i,k,2) = max(qs(i,k,2),qmin)
514 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
519 !----------------------------------------------------------------
520 ! initialize the variables for microphysical physics
560 !----------------------------------------------------------------
561 ! compute the fallout term:
562 ! first, vertical terminal velosity for minor loops
566 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin)then
567 rslopec(i,k) = rslopecmax
568 rslopec2(i,k) = rslopec2max
569 rslopec3(i,k) = rslopec3max
571 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
572 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
573 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
575 !-------------------------------------------------------------
576 ! Ni: ice crystal number concentraiton [HDC 5c]
577 !-------------------------------------------------------------
578 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
579 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
580 temp = (den(i,k)*max(qci(i,k,2),qmin))
581 temp = sqrt(sqrt(temp*temp*temp))
582 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
587 qrs_tmp(i,k,1) = qrs(i,k,1)
588 qrs_tmp(i,k,2) = qrs(i,k,2)
589 ncr_tmp(i,k) = ncr(i,k,3)
592 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
593 rslope3,work1,workn,its,ite,kts,kte)
594 !----------------------------------------------------------------
595 ! compute the fallout term:
596 ! first, vertical terminal velosity for minor loops
597 !----------------------------------------------------------------
599 ! vt update for qr and nr
604 work1(i,k,1) = work1(i,k,1)/delz(i,k)
605 workn(i,k) = workn(i,k)/delz(i,k)
606 numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
607 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
611 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
617 if(n.le.mstep(i)) then
618 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
619 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
620 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
621 falln(i,k) = falln(i,k)+falkn(i,k)
622 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
623 ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
626 do k = kte-1, kts, -1
628 if(n.le.mstep(i)) then
629 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
630 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
631 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
632 falln(i,k) = falln(i,k)+falkn(i,k)
633 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
634 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
635 ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
636 /delz(i,k))*dtcld,0.)
642 qrs_tmp(i,k,1) = qrs(i,k,1)
643 ncr_tmp(i,k) = ncr(i,k,3)
646 call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
647 rslope3,work1,workn,its,ite,kts,kte)
650 work1(i,k,1) = work1(i,k,1)/delz(i,k)
651 workn(i,k) = workn(i,k)/delz(i,k)
658 works(i,k) = work1(i,k,2)
659 denqrs2(i,k) = den(i,k)*qrs(i,k,2)
660 if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
663 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2, &
667 qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
668 fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
672 fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
676 qrs_tmp(i,k,1) = qrs(i,k,1)
677 qrs_tmp(i,k,2) = qrs(i,k,2)
678 ncr_tmp(i,k) = ncr(i,k,3)
681 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
682 rslope3,work1,workn,its,ite,kts,kte)
687 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
688 if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
689 !----------------------------------------------------------------
690 ! psmlt: melting of snow [HL A33] [RH83 A25]
692 !----------------------------------------------------------------
694 ! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
695 work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
696 /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
697 *exp(log(t(i,k))*(1.81))/p(i,k)))) &
698 *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
699 *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
700 *sqrt(sqrt(den0/(den(i,k)))))
701 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
702 ! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
703 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
704 ! *work2(i,k)*coeres)
705 psmlt(i,k) = (1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k))) &
706 /((t(i,k))+120.)/(den(i,k)))*(den(i,k)))/xlf &
707 *(t0c-t(i,k))*pi/2.*n0sfac(i,k) &
708 *(precs1*rslope2(i,k,2)+precs2*work2(i,k)*coeres)
709 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2) &
711 !-------------------------------------------------------------------
712 ! nsmlt: melgin of snow [LH A27]
714 !-------------------------------------------------------------------
715 if(qrs(i,k,2).gt.qcrmin) then
716 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)*mstep(i)/qrs(i,k,2)
717 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
719 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
720 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
721 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
725 !---------------------------------------------------------------
726 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
727 !---------------------------------------------------------------
730 if(qci(i,k,2).le.0.) then
733 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
734 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
735 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
740 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
744 denqci(i,k) = den(i,k)*qci(i,k,2)
747 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
751 qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
755 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
758 !----------------------------------------------------------------
759 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
762 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
763 fallsum_qsi = fall(i,1,2)+fallc(i,1)
765 if(fallsum.gt.0.) then
766 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
767 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000.+rain(i)
769 if (PRESENT (snowncv) .and. PRESENT (snow)) then
771 if(fallsum_qsi.gt.0.) then
772 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
773 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.+snow(i)
777 if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
781 !---------------------------------------------------------------
782 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
784 !---------------------------------------------------------------
789 if(supcol.lt.0.) xlf = xlf0
790 if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
791 qci(i,k,1) = qci(i,k,1)+qci(i,k,2)
792 !---------------------------------------------------------------
793 ! nimlt: instantaneous melting of cloud ice [LH A18]
795 !--------------------------------------------------------------
796 ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
797 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
800 !---------------------------------------------------------------
801 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
803 !---------------------------------------------------------------
804 if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
805 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
806 !---------------------------------------------------------------
807 ! nihmf: homogeneous of cloud water below -40c [LH A17]
809 !---------------------------------------------------------------
810 if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
811 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
814 !---------------------------------------------------------------
815 ! pihtf: heterogeneous freezing of cloud water [HL A44]
816 ! (T0>T>-40C: QC->QI)
817 !---------------------------------------------------------------
818 if(supcol.gt.0. .and. qci(i,k,1).gt.0.) then
819 supcolt=min(supcol,70.)
820 pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k) &
821 *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld,qci(i,k,1))
822 !---------------------------------------------------------------
823 ! nihtf: heterogeneous of cloud water [LH A16]
825 !---------------------------------------------------------------
826 if(ncr(i,k,2).gt.ncmin) then
827 nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2) &
828 *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
829 ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
831 qci(i,k,2) = qci(i,k,2) + pfrzdtc
832 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
833 qci(i,k,1) = qci(i,k,1)-pfrzdtc
835 !---------------------------------------------------------------
836 ! psfrz: freezing of rain water [HL A20] [LFO 45]
838 !---------------------------------------------------------------
839 if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
840 supcolt=min(supcol,70.)
841 pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k) &
842 *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1) &
844 !---------------------------------------------------------------
845 ! nsfrz: freezing of rain water [LH A26]
847 !---------------------------------------------------------------
848 if(ncr(i,k,3).gt.nrmin) then
849 nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.) &
850 *rslope3(i,k,1)*dtcld,ncr(i,k,3))
851 ncr(i,k,3) = ncr(i,k,3)-nfrzdtr
853 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
854 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
855 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
862 ncr(i,k,2) = max(ncr(i,k,2),0.0)
863 ncr(i,k,3) = max(ncr(i,k,3),0.0)
866 !----------------------------------------------------------------
867 ! update the slope parameters for microphysics computation
871 qrs_tmp(i,k,1) = qrs(i,k,1)
872 qrs_tmp(i,k,2) = qrs(i,k,2)
873 ncr_tmp(i,k) = ncr(i,k,3)
876 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
877 rslope3,work1,workn,its,ite,kts,kte)
880 !-----------------------------------------------------------------
881 ! compute the mean-volume drop diameter [LH A10]
882 ! for raindrop distribution
883 !-----------------------------------------------------------------
884 avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
885 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
886 rslopec(i,k) = rslopecmax
887 rslopec2(i,k) = rslopec2max
888 rslopec3(i,k) = rslopec3max
890 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
891 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
892 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
894 !--------------------------------------------------------------------
895 ! compute the mean-volume drop diameter [LH A7]
896 ! for cloud-droplet distribution
897 !--------------------------------------------------------------------
898 avedia(i,k,1) = rslopec(i,k)
901 !----------------------------------------------------------------
902 ! work1: the thermodynamic term in the denominator associated with
903 ! heat conduction and vapor diffusion
905 ! work2: parameter associated with the ventilation effects(y93)
909 ! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
910 work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.) &
911 *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
912 *(den(i,k))*(rv*(t(i,k))*(t(i,k))))) &
913 + p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
914 ! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
915 work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
916 /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
917 *(rv*(t(i,k))*(t(i,k)))) &
918 + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
919 ! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
920 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
921 *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5 &
922 *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
923 /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k)))) &
924 /(((t(i,k))+120.)*den(i,k)))
928 !===============================================================
930 ! warm rain processes
932 ! - follows the processes in RH83 and LFO except for autoconcersion
934 !===============================================================
938 supsat = max(q(i,k),qmin)-qs(i,k,1)
940 !---------------------------------------------------------------
941 ! praut: auto conversion rate from cloud to rain [LH 9] [CP 17]
943 !---------------------------------------------------------------
944 lencon = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k) &
946 lenconcr = max(1.2*lencon,qcrmin)
947 if(avedia(i,k,1).gt.di15) then
948 taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
949 praut(i,k) = lencon/taucon
950 praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
951 !---------------------------------------------------------------
952 ! nraut: auto conversion rate from cloud to rain [LH A6][CP 18 & 19]
954 !---------------------------------------------------------------
955 nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
956 if(qrs(i,k,1).gt.lenconcr) &
957 nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
958 nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
960 !---------------------------------------------------------------
961 ! pracw: accretion of cloud water by rain [LH 10][CP 22 & 23]
963 ! nracw: accretion of cloud water by rain [LH A9]
965 !---------------------------------------------------------------
966 if(qrs(i,k,1).ge.lenconcr) then
967 if(avedia(i,k,2).ge.di100) then
968 nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k) &
969 + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
970 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2) &
971 *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k) &
972 + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
974 nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k) &
975 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
976 *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
977 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2) &
978 *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k) &
979 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
980 *rslope3(i,k,1)),qci(i,k,1)/dtcld)
983 !----------------------------------------------------------------
984 ! nccol: self collection of cloud water [LH A8][CP 24 & 25]
986 !----------------------------------------------------------------
987 if(avedia(i,k,1).ge.di100) then
988 nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
990 nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) &
993 !----------------------------------------------------------------
994 ! nrcol: self collection of rain-drops and break-up [LH A21][CP 24 & 25]
996 !----------------------------------------------------------------
997 if(qrs(i,k,1).ge.lenconcr) then
998 if(avedia(i,k,2).lt.di100) then
999 nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) &
1001 elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
1002 nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1003 elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
1004 coecol = -2.5e3*(avedia(i,k,2)-di600)
1005 nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3) &
1011 !---------------------------------------------------------------
1012 ! prevp: evaporation/condensation rate of rain [HL A41]
1013 ! (QV->QR or QR->QV)
1014 !---------------------------------------------------------------
1015 if(qrs(i,k,1).gt.0.) then
1016 coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1017 prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1) &
1018 +precr2*work2(i,k)*coeres)/work1(i,k,1)
1019 if(prevp(i,k).lt.0.) then
1020 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1021 prevp(i,k) = max(prevp(i,k),satdt/2)
1022 !----------------------------------------------------------------
1023 ! Nrevp: evaporation/condensation rate of rain [LH A14]
1025 !----------------------------------------------------------------
1026 if(avedia(i,k,2).le.di82) then
1027 nrevp(i,k) = ncr(i,k,3)/dtcld
1028 !----------------------------------------------------------------
1029 ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
1031 !----------------------------------------------------------------
1032 prevp_s(i,k) = qrs(i,k,1)/dtcld
1036 prevp(i,k) = min(prevp(i,k),satdt/2)
1042 !===============================================================
1044 ! cold rain processes
1046 ! - follows the revised ice microphysics processes in HDC
1047 ! - the processes same as in RH83 and RH84 and LFO behave
1048 ! following ice crystal hapits defined in HDC, inclduing
1049 ! intercept parameter for snow (n0s), ice crystal number
1050 ! concentration (ni), ice nuclei number concentration
1051 ! (n0i), ice diameter (d)
1053 !===============================================================
1059 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1060 supsat = max(q(i,k),qmin)-qs(i,k,2)
1061 satdt = supsat/dtcld
1063 !-------------------------------------------------------------
1064 ! Ni: ice crystal number concentraiton [HDC 5c]
1065 !-------------------------------------------------------------
1066 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
1067 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1068 temp = (den(i,k)*max(qci(i,k,2),qmin))
1069 temp = sqrt(sqrt(temp*temp*temp))
1070 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1071 eacrs = exp(0.07*(-supcol))
1073 if(supcol.gt.0) then
1074 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,2).gt.qmin) then
1075 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1076 diameter = min(dicon * sqrt(xmi),dimax)
1077 vt2i = 1.49e4*diameter**1.31
1078 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
1079 !-------------------------------------------------------------
1080 ! psaci: Accretion of cloud ice by rain [HDC 10]
1082 !-------------------------------------------------------------
1083 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1084 + diameter**2*rslope(i,k,2)
1085 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)*abs(vt2s-vt2i) &
1089 !-------------------------------------------------------------
1090 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
1091 ! (T<T0: QC->QS, and T>=T0: QC->QR)
1092 !-------------------------------------------------------------
1093 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1094 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1095 *qci(i,k,1)*denfac(i,k),qci(i,k,1)*rdtcld)
1097 !-------------------------------------------------------------
1098 ! nsacw: Accretion of cloud water by snow [LH A12]
1100 !-------------------------------------------------------------
1101 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1102 nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1103 *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1105 if(supcol.le.0) then
1107 !--------------------------------------------------------------
1108 ! nseml: Enhanced melting of snow by accretion of water [LH A29]
1110 !--------------------------------------------------------------
1111 if (qrs(i,k,2).gt.qcrmin) then
1112 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1113 nseml(i,k) = -sfac*min(max(cliq*supcol*(psacw(i,k))/xlf &
1114 ,-qrs(i,k,2)/dtcld),0.)
1117 if(supcol.gt.0) then
1118 !-------------------------------------------------------------
1119 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1120 ! (T<T0: QV->QI or QI->QV)
1121 !-------------------------------------------------------------
1122 if(qci(i,k,2).gt.0 .and. ifsat.ne.1) then
1123 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1124 diameter = dicon * sqrt(xmi)
1125 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1126 supice = satdt-prevp(i,k)
1127 if(pidep(i,k).lt.0.) then
1128 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1129 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1130 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1131 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1133 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1134 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1136 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1138 !-------------------------------------------------------------
1139 ! psdep: deposition/sublimation rate of snow [HDC 14]
1140 ! (QV->QS or QS->QV)
1141 !-------------------------------------------------------------
1142 if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
1143 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1144 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1145 + precs2*work2(i,k)*coeres)/work1(i,k,2)
1146 supice = satdt-prevp(i,k)-pidep(i,k)
1147 if(psdep(i,k).lt.0.) then
1148 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1149 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1150 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1151 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1153 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1154 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1156 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1158 !-------------------------------------------------------------
1159 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1161 !-------------------------------------------------------------
1162 if(supsat.gt.0 .and. ifsat.ne.1) then
1163 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1164 xni0 = 1.e3*exp(0.1*supcol)
1165 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1166 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))*rdtcld)
1167 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1170 !-------------------------------------------------------------
1171 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1173 !-------------------------------------------------------------
1174 if(qci(i,k,2).gt.0.) then
1175 qimax = roqimax/den(i,k)
1176 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1177 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1180 !-------------------------------------------------------------
1181 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1183 !-------------------------------------------------------------
1184 if(supcol.lt.0.) then
1185 if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) &
1186 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1187 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1188 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1194 !----------------------------------------------------------------
1195 ! check mass conservation of generation terms and feedback to the
1200 if(t(i,k).le.t0c) then
1204 value = max(qmin,qci(i,k,1))
1205 source = (praut(i,k)+pracw(i,k)+psacw(i,k)-prevp_s(i,k))*dtcld
1206 if (source.gt.value) then
1207 factor = value/source
1208 praut(i,k) = praut(i,k)*factor
1209 pracw(i,k) = pracw(i,k)*factor
1210 psacw(i,k) = psacw(i,k)*factor
1211 prevp_s(i,k) = prevp_s(i,k)*factor
1216 value = max(qmin,qci(i,k,2))
1217 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1218 if (source.gt.value) then
1219 factor = value/source
1220 psaut(i,k) = psaut(i,k)*factor
1221 psaci(i,k) = psaci(i,k)*factor
1222 pigen(i,k) = pigen(i,k)*factor
1223 pidep(i,k) = pidep(i,k)*factor
1229 value = max(qmin,qrs(i,k,1))
1230 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)+prevp_s(i,k))*dtcld
1231 if (source.gt.value) then
1232 factor = value/source
1233 praut(i,k) = praut(i,k)*factor
1234 pracw(i,k) = pracw(i,k)*factor
1235 prevp(i,k) = prevp(i,k)*factor
1236 prevp_s(i,k) = prevp_s(i,k)*factor
1241 value = max(qmin,qrs(i,k,2))
1242 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1243 if (source.gt.value) then
1244 factor = value/source
1245 psdep(i,k) = psdep(i,k)*factor
1246 psaut(i,k) = psaut(i,k)*factor
1247 psaci(i,k) = psaci(i,k)*factor
1248 psacw(i,k) = psacw(i,k)*factor
1253 value = max(ncmin,ncr(i,k,2))
1254 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k) &
1256 if (source.gt.value) then
1257 factor = value/source
1258 nraut(i,k) = nraut(i,k)*factor
1259 nccol(i,k) = nccol(i,k)*factor
1260 nracw(i,k) = nracw(i,k)*factor
1261 nsacw(i,k) = nsacw(i,k)*factor
1262 nrevp(i,k) = nrevp(i,k)*factor
1267 value = max(nrmin,ncr(i,k,3))
1268 source = (-nraut(i,k)+nrcol(i,k)+nrevp(i,k))*dtcld
1269 if (source.gt.value) then
1270 factor = value/source
1271 nraut(i,k) = nraut(i,k)*factor
1272 nrcol(i,k) = nrcol(i,k)*factor
1273 nrevp(i,k) = nrevp(i,k)*factor
1276 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1278 q(i,k) = q(i,k)+work2(i,k)*dtcld
1279 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k) &
1280 +prevp_s(i,k))*dtcld,0.)
1281 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k) &
1282 -prevp_s(i,k))*dtcld,0.)
1283 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)-pigen(i,k) &
1284 -pidep(i,k))*dtcld,0.)
1285 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+psaci(i,k) &
1286 +psacw(i,k))*dtcld,0.)
1287 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1288 -nsacw(i,k)+nrevp(i,k))*dtcld,0.)
1289 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-nrevp(i,k)) &
1292 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1293 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1294 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1299 value = max(qmin,qci(i,k,1))
1300 source=(praut(i,k)+pracw(i,k)+psacw(i,k)-prevp_s(i,k))*dtcld
1301 if (source.gt.value) then
1302 factor = value/source
1303 praut(i,k) = praut(i,k)*factor
1304 pracw(i,k) = pracw(i,k)*factor
1305 psacw(i,k) = psacw(i,k)*factor
1306 prevp_s(i,k) = prevp_s(i,k)*factor
1311 value = max(qmin,qrs(i,k,1))
1312 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)+prevp_s(i,k) &
1314 if (source.gt.value) then
1315 factor = value/source
1316 praut(i,k) = praut(i,k)*factor
1317 pracw(i,k) = pracw(i,k)*factor
1318 prevp(i,k) = prevp(i,k)*factor
1319 psacw(i,k) = psacw(i,k)*factor
1320 prevp_s(i,k) = prevp_s(i,k)*factor
1325 value = max(qcrmin,qrs(i,k,2))
1326 source=(-psevp(i,k))*dtcld
1327 if (source.gt.value) then
1328 factor = value/source
1329 psevp(i,k) = psevp(i,k)*factor
1334 value = max(ncmin,ncr(i,k,2))
1335 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k) &
1337 if (source.gt.value) then
1338 factor = value/source
1339 nraut(i,k) = nraut(i,k)*factor
1340 nccol(i,k) = nccol(i,k)*factor
1341 nracw(i,k) = nracw(i,k)*factor
1342 nsacw(i,k) = nsacw(i,k)*factor
1343 nrevp(i,k) = nrevp(i,k)*factor
1348 value = max(nrmin,ncr(i,k,3))
1349 source = (-nraut(i,k)-nseml(i,k)+nrcol(i,k)+nrevp(i,k))*dtcld
1350 if (source.gt.value) then
1351 factor = value/source
1352 nraut(i,k) = nraut(i,k)*factor
1353 nseml(i,k) = nseml(i,k)*factor
1354 nrcol(i,k) = nrcol(i,k)*factor
1355 nrevp(i,k) = nrevp(i,k)*factor
1357 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1359 q(i,k) = q(i,k)+work2(i,k)*dtcld
1360 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k) &
1361 +prevp_s(i,k))*dtcld,0.)
1362 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k) &
1363 +psacw(i,k)-prevp_s(i,k))*dtcld,0.)
1364 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1365 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1366 -nsacw(i,k)+nrevp(i,k))*dtcld,0.)
1367 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)+nseml(i,k)-nrcol(i,k) &
1368 -nrevp(i,k))*dtcld,0.)
1370 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1371 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1376 ! Inline expansion for fpvs
1377 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1378 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1388 xbi=xai+hsub/(rv*ttp)
1393 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1394 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1395 qs(i,k,1) = max(qs(i,k,1),qmin)
1396 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
1402 !-------------------------------------------------------------------
1403 ! rate of change of cloud drop concentration due to CCN activation
1404 ! pcact: QV -> QC [LH 8] [KK 14]
1405 ! ncact: NCCN -> NC [LH A2] [KK 12]
1406 !-------------------------------------------------------------------
1407 if(rh(i,k,1).gt.1.) then
1408 ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) &
1409 *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
1410 ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
1411 pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/ &
1412 (3.*den(i,k)),max(q(i,k),0.)/dtcld)
1413 q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
1414 qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
1415 ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
1416 ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
1417 t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
1419 !---------------------------------------------------------------------
1420 ! pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1421 ! if there exists additional water vapor condensated/if
1422 ! evaporation of cloud water is not enough to remove subsaturation
1423 ! (QV->QC or QC->QV)
1424 !---------------------------------------------------------------------
1426 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1427 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1428 qs(i,k,1) = max(qs(i,k,1),qmin)
1429 work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1430 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k)) &
1432 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1433 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1434 if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) &
1435 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1436 !----------------------------------------------------------------------
1437 ! ncevp: evpration of Cloud number concentration [LH A3]
1439 !----------------------------------------------------------------------
1440 if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
1442 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
1445 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1446 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1447 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1451 !----------------------------------------------------------------
1452 ! padding for small values
1456 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1457 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1461 END SUBROUTINE wdm52d
1462 ! ...................................................................
1463 REAL FUNCTION rgmma(x)
1464 !-------------------------------------------------------------------
1466 !-------------------------------------------------------------------
1467 ! rgmma function: use infinite product form
1469 PARAMETER (euler=0.577215664901532)
1475 rgmma=x*exp(euler*x)
1478 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1484 !--------------------------------------------------------------------------
1485 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1486 !--------------------------------------------------------------------------
1488 !--------------------------------------------------------------------------
1489 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1492 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1499 xbi=xai+hsub/(rv*ttp)
1501 if(t.lt.ttp .and. ice.eq.1) then
1502 fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1504 fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1506 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1508 !-------------------------------------------------------------------
1509 SUBROUTINE wdm5init(den0,denr,dens,cl,cpv,ccn0,allowed_to_read)
1510 !-------------------------------------------------------------------
1512 !-------------------------------------------------------------------
1513 !.... constants which may not be tunable
1514 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
1515 LOGICAL, INTENT(IN) :: allowed_to_read
1520 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1521 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1530 bvtr2o5 = 2.5+.5*bvtr
1531 bvtr3o5 = 3.5+.5*bvtr
1532 g1pbr = rgmma(bvtr1)
1533 g2pbr = rgmma(bvtr2)
1534 g3pbr = rgmma(bvtr3)
1535 g4pbr = rgmma(bvtr4) ! 17.837825
1536 g5pbr = rgmma(bvtr5)
1537 g7pbr = rgmma(bvtr7)
1538 g5pbro2 = rgmma(bvtr2o5)
1539 g7pbro2 = rgmma(bvtr3o5)
1540 pvtr = avtr*g5pbr/24.
1543 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1545 precr2 = 2.*pi*.31*avtr**.5*g7pbro2
1546 pidn0r = pi*denr*n0r
1548 xmmax = (dimax/dicon)**2
1549 roqimax = 2.08e22*dimax**8
1555 g1pbs = rgmma(bvts1) !.8875
1556 g3pbs = rgmma(bvts3)
1557 g4pbs = rgmma(bvts4) ! 12.0786
1558 g5pbso2 = rgmma(bvts2)
1559 pvts = avts*g4pbs/6.
1560 pacrs = pi*n0s*avts*g3pbs*.25
1562 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1563 pidn0s = pi*dens*n0s
1564 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1566 rslopecmax = 1./lamdacmax
1567 rslopermax = 1./lamdarmax
1568 rslopesmax = 1./lamdasmax
1569 rsloperbmax = rslopermax ** bvtr
1570 rslopesbmax = rslopesmax ** bvts
1571 rslopec2max = rslopecmax * rslopecmax
1572 rsloper2max = rslopermax * rslopermax
1573 rslopes2max = rslopesmax * rslopesmax
1574 rslopec3max = rslopec2max * rslopecmax
1575 rsloper3max = rsloper2max * rslopermax
1576 rslopes3max = rslopes2max * rslopesmax
1578 END SUBROUTINE wdm5init
1579 !------------------------------------------------------------------------------
1580 subroutine slope_wdm5(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1581 vt,vtn,its,ite,kts,kte)
1583 INTEGER :: its,ite, jts,jte, kts,kte
1584 REAL, DIMENSION( its:ite , kts:kte,2) :: &
1591 REAL, DIMENSION( its:ite , kts:kte) :: &
1597 REAL, PARAMETER :: t0c = 273.15
1598 REAL, DIMENSION( its:ite , kts:kte ) :: &
1600 REAL :: lamdar, lamdas, x, y, z, supcol
1602 !----------------------------------------------------------------
1603 ! size distributions: (x=mixing ratio, y=air density):
1604 ! valid for mixing ratio > 1.e-9 kg/kg.
1606 ! Optimizatin : A**B => exp(log(A)*(B))
1607 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1608 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1613 !---------------------------------------------------------------
1614 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1615 !---------------------------------------------------------------
1616 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1617 if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
1618 rslope(i,k,1) = rslopermax
1619 rslopeb(i,k,1) = rsloperbmax
1620 rslope2(i,k,1) = rsloper2max
1621 rslope3(i,k,1) = rsloper3max
1623 rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
1624 rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1625 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1626 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1628 if(qrs(i,k,2).le.qcrmin) then
1629 rslope(i,k,2) = rslopesmax
1630 rslopeb(i,k,2) = rslopesbmax
1631 rslope2(i,k,2) = rslopes2max
1632 rslope3(i,k,2) = rslopes3max
1634 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1635 rslopeb(i,k,2) = rslope(i,k,2)**bvts
1636 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1637 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1639 vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1640 vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1641 vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
1642 if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1643 if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1644 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1647 END subroutine slope_wdm5
1648 !-----------------------------------------------------------------------------
1649 subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1650 vt,vtn,its,ite,kts,kte)
1652 INTEGER :: its,ite, jts,jte, kts,kte
1653 REAL, DIMENSION( its:ite , kts:kte) :: &
1665 REAL, PARAMETER :: t0c = 273.15
1666 REAL, DIMENSION( its:ite , kts:kte ) :: &
1668 REAL :: lamdar, x, y, z, supcol
1670 !----------------------------------------------------------------
1671 ! size distributions: (x=mixing ratio, y=air density):
1672 ! valid for mixing ratio > 1.e-9 kg/kg.
1673 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1677 if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
1678 rslope(i,k) = rslopermax
1679 rslopeb(i,k) = rsloperbmax
1680 rslope2(i,k) = rsloper2max
1681 rslope3(i,k) = rsloper3max
1683 rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
1684 rslopeb(i,k) = rslope(i,k)**bvtr
1685 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1686 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1688 vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1689 vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
1690 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1691 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1694 END subroutine slope_rain
1695 !------------------------------------------------------------------------------
1696 subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1699 INTEGER :: its,ite, jts,jte, kts,kte
1700 REAL, DIMENSION( its:ite , kts:kte) :: &
1710 REAL, PARAMETER :: t0c = 273.15
1711 REAL, DIMENSION( its:ite , kts:kte ) :: &
1713 REAL :: lamdas, x, y, z, supcol
1715 !----------------------------------------------------------------
1716 ! size distributions: (x=mixing ratio, y=air density):
1717 ! valid for mixing ratio > 1.e-9 kg/kg.
1718 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1723 !---------------------------------------------------------------
1724 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1725 !---------------------------------------------------------------
1726 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1727 if(qrs(i,k).le.qcrmin)then
1728 rslope(i,k) = rslopesmax
1729 rslopeb(i,k) = rslopesbmax
1730 rslope2(i,k) = rslopes2max
1731 rslope3(i,k) = rslopes3max
1733 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1734 rslopeb(i,k) = rslope(i,k)**bvts
1735 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1736 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1738 vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1739 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1742 END subroutine slope_snow
1743 !-------------------------------------------------------------------
1744 SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1745 !-------------------------------------------------------------------
1747 ! for non-iteration semi-Lagrangain forward advection for cloud
1748 ! with mass conservation and positive definite advection
1749 ! 2nd order interpolation with monotonic piecewise linear method
1750 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1752 ! dzl depth of model layer in meter
1753 ! wwl terminal velocity at model layer m/s
1754 ! rql cloud density*mixing ration
1755 ! precip precipitation
1757 ! id kind of precip: 0 test case; 1 raindrop 2: snow
1758 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1759 ! 0 : use departure wind for advection
1760 ! 1 : use mean wind for advection
1761 ! > 1 : use mean wind after iter-1 iterations
1763 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1764 ! implemented by song-you hong
1769 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1770 real denl(im,km),denfacl(im,km),tkl(im,km)
1772 integer i,k,n,m,kk,kb,kt,iter
1773 real tl,tl2,qql,dql,qqd
1775 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1776 real allold, allnew, zz, dzamin, cflmax, decfl
1777 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1778 real den(km), denfac(km), tk(km)
1779 real wi(km+1), zi(km+1), za(km+1)
1780 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1781 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1786 ! -----------------------------------
1791 denfac(:) = denfacl(i,:)
1793 ! skip for no precipitation for all layers
1796 allold = allold + qq(k)
1798 if(allold.le.0.0) then
1802 ! compute interface values
1805 zi(k+1) = zi(k)+dz(k)
1808 ! save departure wind
1812 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1813 ! 2nd order interpolation to get wi
1817 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1819 ! 3rd order interpolation to get wi
1823 wi(2) = 0.5*(ww(2)+ww(1))
1825 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1827 wi(km) = 0.5*(ww(km)+ww(km-1))
1830 ! terminate of top of raingroup
1832 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1838 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1839 if( decfl .gt. con1 ) then
1840 wi(k) = wi(k+1) - con1*dz(k)/dt
1843 ! compute arrival point
1845 za(k) = zi(k) - wi(k)*dt
1849 dza(k) = za(k+1)-za(k)
1851 dza(km+1) = zi(km+1) - za(km+1)
1853 ! computer deformation at arrival point
1855 qa(k) = qq(k)*dz(k)/dza(k)
1856 qr(k) = qa(k)/den(k)
1859 ! call maxmin(km,1,qa,' arrival points ')
1861 ! compute arrival terminal velocity, and estimate mean terminal velocity
1862 ! then back to use mean terminal velocity
1863 if( n.le.iter ) then
1866 ! call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1868 call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1870 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1873 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1875 ! mean wind is average of departure and new arrival winds
1876 ww(k) = 0.5* ( wd(k)+wa(k) )
1883 ! estimate values at arrival cell interface with monotone
1885 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1886 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1887 if( dip*dim.le.0.0 ) then
1891 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1892 qmi(k)=2.0*qa(k)-qpi(k)
1893 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1904 ! interpolation to regular point
1912 if( zi(k).ge.za(km+1) ) then
1915 find_kb : do kk=kb,km
1916 if( zi(k).le.za(kk+1) ) then
1923 find_kt : do kk=kt,km
1924 if( zi(k+1).le.za(kk) ) then
1932 ! compute q with piecewise constant method
1934 tl=(zi(k)-za(kb))/dza(kb)
1935 th=(zi(k+1)-za(kb))/dza(kb)
1938 qqd=0.5*(qpi(kb)-qmi(kb))
1939 qqh=qqd*th2+qmi(kb)*th
1940 qql=qqd*tl2+qmi(kb)*tl
1941 qn(k) = (qqh-qql)/(th-tl)
1942 else if( kt.gt.kb ) then
1943 tl=(zi(k)-za(kb))/dza(kb)
1945 qqd=0.5*(qpi(kb)-qmi(kb))
1946 qql=qqd*tl2+qmi(kb)*tl
1948 zsum = (1.-tl)*dza(kb)
1950 if( kt-kb.gt.1 ) then
1952 zsum = zsum + dza(m)
1953 qsum = qsum + qa(m) * dza(m)
1956 th=(zi(k+1)-za(kt))/dza(kt)
1958 qqd=0.5*(qpi(kt)-qmi(kt))
1959 dqh=qqd*th2+qmi(kt)*th
1960 zsum = zsum + th*dza(kt)
1961 qsum = qsum + dqh*dza(kt)
1970 sum_precip: do k=1,km
1971 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1972 precip(i) = precip(i) + qa(k)*dza(k)
1974 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1975 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1981 ! replace the new values
1984 ! ----------------------------------
1987 END SUBROUTINE nislfv_rain_plm
1988 END MODULE module_mp_wdm5