source: trunk/WRF.COMMON/WRFV3/phys/module_mp_wsm3.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size

File size: 39.1 KB
Line 
1#if ( RWORDSIZE == 4 )
2#  define VREC vsrec
3#  define VSQRT vssqrt
4#else
5#  define VREC vrec
6#  define VSQRT vsqrt
7#endif
8
9MODULE module_mp_wsm3
10!
11!
12   REAL, PARAMETER, PRIVATE :: dtcldcr     = 120.
13   REAL, PARAMETER, PRIVATE :: n0r = 8.e6
14   REAL, PARAMETER, PRIVATE :: avtr = 841.9
15   REAL, PARAMETER, PRIVATE :: bvtr = 0.8
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
21   REAL, PARAMETER, PRIVATE :: bvts = .41
22   REAL, PARAMETER, PRIVATE :: n0smax =  1.e11 ! t=-90C unlimited
23   REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
24   REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
25   REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
26   REAL, PARAMETER, PRIVATE :: betai = .6
27   REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
28   REAL, PARAMETER, PRIVATE :: dicon = 11.9
29   REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
30   REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
31   REAL, PARAMETER, PRIVATE :: n0s = 2.e6             ! temperature dependent n0s
32   REAL, PARAMETER, PRIVATE :: alpha = .12        ! .122 exponen factor for n0s
33   REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
34   REAL, SAVE ::                                     &
35             qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
36             g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr,   &
37             precr1,precr2,xm0,xmmax,roqimax,bvts1,  &
38             bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs,    &
39             g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
40             pidn0s,xlv1,                            &
41             rslopermax,rslopesmax,rslopegmax,       &
42             rsloperbmax,rslopesbmax,rslopegbmax,    &
43             rsloper2max,rslopes2max,rslopeg2max,    &
44             rsloper3max,rslopes3max,rslopeg3max
45!
46! Specifies code-inlining of fpvs function in WSM32D below. JM 20040507
47!
48CONTAINS
49!===================================================================
50!
51  SUBROUTINE wsm3(th, q, qci, qrs                                  &
52                   , w, den, pii, p, delz                          &
53                   , delt,g, cpd, cpv, rd, rv, t0c                 &
54                   , ep1, ep2, qmin                                &
55                   , XLS, XLV0, XLF0, den0, denr                   &
56                   , cliq,cice,psat                                &
57                   , rain, rainncv                                 &
58                   ,snow, snowncv                                  &
59                   ,sr                                             &
60                   , ids,ide, jds,jde, kds,kde                     &
61                   , ims,ime, jms,jme, kms,kme                     &
62                   , its,ite, jts,jte, kts,kte                     &
63                                                                   )
64!-------------------------------------------------------------------
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
106  REAL, INTENT(IN   ) ::                                    delt, &
107                                                               g, &
108                                                              rd, &
109                                                              rv, &
110                                                             t0c, &
111                                                            den0, &
112                                                             cpd, &
113                                                             cpv, &
114                                                             ep1, &
115                                                             ep2, &
116                                                            qmin, &
117                                                             XLS, &
118                                                            XLV0, &
119                                                            XLF0, &
120                                                            cliq, &
121                                                            cice, &
122                                                            psat, &
123                                                            denr
124  REAL, DIMENSION( ims:ime , jms:jme ),                           &
125        INTENT(INOUT) ::                                    rain, &
126                                                         rainncv, &
127                                                              sr
128
129  REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                &
130        INTENT(INOUT) ::                                    snow, &
131                                                         snowncv
132
133! LOCAL VAR
134  REAL, DIMENSION( its:ite , kts:kte ) ::   t
135  INTEGER ::               i,j,k
136!-------------------------------------------------------------------
137      DO j=jts,jte
138         DO k=kts,kte
139         DO i=its,ite
140            t(i,k)=th(i,k,j)*pii(i,k,j)
141         ENDDO
142         ENDDO
143         CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j)               &
144                    ,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j)   &
145                    ,p(ims,kms,j), delz(ims,kms,j)                 &
146                    ,delt,g, cpd, cpv, rd, rv, t0c                 &
147                    ,ep1, ep2, qmin                                &
148                    ,XLS, XLV0, XLF0, den0, denr                   &
149                    ,cliq,cice,psat                                &
150                    ,j                                             &
151                    ,rain(ims,j), rainncv(ims,j)                   &
152                    ,sr(ims,j)                                     &
153                    ,ids,ide, jds,jde, kds,kde                     &
154                    ,ims,ime, jms,jme, kms,kme                     &
155                    ,its,ite, jts,jte, kts,kte                     &
156                    ,snow(ims,j),snowncv(ims,j)                    &
157                                                                   )
158         DO K=kts,kte
159         DO I=its,ite
160            th(i,k,j)=t(i,k)/pii(i,k,j)
161         ENDDO
162         ENDDO
163      ENDDO
164  END SUBROUTINE wsm3
165!===================================================================
166!
167  SUBROUTINE wsm32D(t, q, qci, qrs,w, den, p, delz                &
168                   ,delt,g, cpd, cpv, rd, rv, t0c                 &
169                   ,ep1, ep2, qmin                                &
170                   ,XLS, XLV0, XLF0, den0, denr                   &
171                   ,cliq,cice,psat                                &
172                   ,lat                                           &
173                   ,rain, rainncv                                 &
174                   ,sr                                            &
175                   ,ids,ide, jds,jde, kds,kde                     &
176                   ,ims,ime, jms,jme, kms,kme                     &
177                   ,its,ite, jts,jte, kts,kte                     &
178                   ,snow,snowncv                                  &
179                                                                  )
180!-------------------------------------------------------------------
181  IMPLICIT NONE
182!-------------------------------------------------------------------
183  INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
184                                      ims,ime, jms,jme, kms,kme , &
185                                      its,ite, jts,jte, kts,kte,  &
186                                      lat
187  REAL, DIMENSION( its:ite , kts:kte ),                           &
188        INTENT(INOUT) ::                                          &
189                                                               t
190  REAL, DIMENSION( ims:ime , kms:kme ),                           &
191        INTENT(INOUT) ::                                          &
192                                                               q, &
193                                                             qci, &
194                                                             qrs
195  REAL, DIMENSION( ims:ime , kms:kme ),                           &
196        INTENT(IN   ) ::                                       w, &
197                                                             den, &
198                                                               p, &
199                                                            delz
200  REAL, INTENT(IN   ) ::                                    delt, &
201                                                               g, &
202                                                             cpd, &
203                                                             cpv, &
204                                                             t0c, &
205                                                            den0, &
206                                                              rd, &
207                                                              rv, &
208                                                             ep1, &
209                                                             ep2, &
210                                                            qmin, &
211                                                             XLS, &
212                                                            XLV0, &
213                                                            XLF0, &
214                                                            cliq, &
215                                                            cice, &
216                                                            psat, &
217                                                            denr
218  REAL, DIMENSION( ims:ime ),                                     &
219        INTENT(INOUT) ::                                    rain, &
220                                                         rainncv, &
221                                                              sr
222
223  REAL, DIMENSION( ims:ime ),     OPTIONAL,                       &
224        INTENT(INOUT) ::                                    snow, &
225                                                         snowncv
226! LOCAL VAR
227  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
228        rh, qs, denfac, rslope, rslope2, rslope3, rslopeb,        &
229        pgen, paut, pacr, pisd, pres, pcon, fall, falk,           &
230        xl, cpm, work1, work2, xni, qs0, n0sfac
231  INTEGER, DIMENSION( its:ite ) ::  kwork1, kwork2
232  REAL, DIMENSION( its:ite , kts:kte ) ::                         &
233              falkc, work1c, work2c, fallc
234! variables for optimization
235  REAL, DIMENSION( its:ite )           :: tvec1
236  INTEGER, DIMENSION( its:ite ) :: mstep, numdt
237  LOGICAL, DIMENSION( its:ite ) :: flgcld
238  REAL  ::  pi,                                                   &
239            cpmcal, xlcal, lamdar, lamdas, diffus,                &
240            viscos, xka, venfac, conden, diffac,                  &
241            x, y, z, a, b, c, d, e,                               &
242            fallsum, fallsum_qsi, vt2i,vt2s,acrfac,               &     
243            qdt, pvt, qik, delq, facq, qrsci, frzmlt,             &
244            snomlt, hold, holdrs, facqci, supcol, coeres,         &
245            supsat, dtcld, xmi, qciik, delqci, eacrs, satdt,      &
246            qimax, diameter, xni0, roqi0, supice
247  REAL  :: holdc, holdci
248  INTEGER :: i, j, k, mstepmax,                                   &
249            iprt, latd, lond, loop, loops, ifsat, kk, n
250! Temporaries used for inlining fpvs function
251  REAL  :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
252!
253!=================================================================
254!   compute internal functions
255!
256      cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
257      xlcal(x) = xlv0-xlv1*(x-t0c)
258!----------------------------------------------------------------
259!     size distributions: (x=mixing ratio, y=air density):
260!     valid for mixing ratio > 1.e-9 kg/kg.
261!
262! Optimizatin : A**B => exp(log(A)*(B))
263      lamdar(x,y)=   sqrt(sqrt(pidn0r/(x*y)))      ! (pidn0r/(x*y))**.25
264      lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))    ! (pidn0s*z/(x*y))**.25
265!
266!----------------------------------------------------------------
267!     diffus: diffusion coefficient of the water vapor
268!     viscos: kinematic viscosity(m2s-1)
269!
270      diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y        ! 8.794e-5*x**1.81/y
271      viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y  ! 1.496e-6*x**1.5/(x+120.)/y
272      xka(x,y) = 1.414e3*viscos(x,y)*y
273      diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
274!      venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333)       &
275!                      /viscos(b,c)**(.5)*(den0/c)**0.25
276      venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333)))    &
277                     /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
278      conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
279!
280      pi = 4. * atan(1.)
281!
282!----------------------------------------------------------------
283!     paddint 0 for negative values generated by dynamics
284!
285      do k = kts, kte
286        do i = its, ite
287          qci(i,k) = max(qci(i,k),0.0)
288          qrs(i,k) = max(qrs(i,k),0.0)
289        enddo
290      enddo
291!
292!----------------------------------------------------------------
293!     latent heat for phase changes and heat capacity. neglect the
294!     changes during microphysical process calculation
295!     emanuel(1994)
296!
297      do k = kts, kte
298        do i = its, ite
299          cpm(i,k) = cpmcal(q(i,k))
300          xl(i,k) = xlcal(t(i,k))
301        enddo
302      enddo
303!
304!----------------------------------------------------------------
305!     compute the minor time steps.
306!
307      loops = max(nint(delt/dtcldcr),1)
308      dtcld = delt/loops
309      if(delt.le.dtcldcr) dtcld = delt
310!
311      do loop = 1,loops
312!
313!----------------------------------------------------------------
314!     initialize the large scale variables
315!
316      do i = its, ite
317        mstep(i) = 1
318        flgcld(i) = .true.
319      enddo
320!
321!     do k = kts, kte
322!       do i = its, ite
323!         denfac(i,k) = sqrt(den0/den(i,k))
324!       enddo
325!     enddo
326      do k = kts, kte
327        CALL VREC( tvec1(its), den(its,k), ite-its+1)
328        do i = its, ite
329          tvec1(i) = tvec1(i)*den0
330        enddo
331        CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
332      enddo
333!
334! Inline expansion for fpvs
335!         qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
336!         qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
337      cvap = cpv
338      hvap=xlv0
339      hsub=xls
340      ttp=t0c+0.01
341      dldt=cvap-cliq
342      xa=-dldt/rv
343      xb=xa+hvap/(rv*ttp)
344      dldti=cvap-cice
345      xai=-dldti/rv
346      xbi=xai+hsub/(rv*ttp)
347      do k = kts, kte
348        do i = its, ite
349!         tr=ttp/t(i,k)
350!         if(t(i,k).lt.ttp) then
351!           qs(i,k) =psat*(tr**xai)*exp(xbi*(1.-tr))
352!         else
353!           qs(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
354!         endif
355!         qs0(i,k)  =psat*(tr**xa)*exp(xb*(1.-tr))
356          tr=ttp/t(i,k)
357          if(t(i,k).lt.ttp) then
358            qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
359          else
360            qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
361          endif
362          qs0(i,k)  =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
363          qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
364          qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
365          qs(i,k) = max(qs(i,k),qmin)
366          rh(i,k) = max(q(i,k) / qs(i,k),qmin)
367        enddo
368      enddo
369!
370!----------------------------------------------------------------
371!     initialize the variables for microphysical physics
372!
373!
374      do k = kts, kte
375        do i = its, ite
376          pres(i,k) = 0.
377          paut(i,k) = 0.
378          pacr(i,k) = 0.
379          pgen(i,k) = 0.
380          pisd(i,k) = 0.
381          pcon(i,k) = 0.
382          fall(i,k) = 0.
383          falk(i,k) = 0.
384          fallc(i,k) = 0.
385          falkc(i,k) = 0.
386          xni(i,k) = 1.e3
387        enddo
388      enddo
389!
390!----------------------------------------------------------------
391!     compute the fallout term:
392!     first, vertical terminal velosity for minor loops
393!---------------------------------------------------------------
394! n0s: Intercept parameter for snow [m-4] [HDC 6]
395!---------------------------------------------------------------
396      do k = kts, kte
397        do i = its, ite
398          supcol = t0c-t(i,k)
399          n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
400          if(t(i,k).ge.t0c) then
401            if(qrs(i,k).le.qcrmin)then
402              rslope(i,k) = rslopermax
403              rslopeb(i,k) = rsloperbmax
404              rslope2(i,k) = rsloper2max
405              rslope3(i,k) = rsloper3max
406            else
407              rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
408!             rslopeb(i,k) = rslope(i,k)**bvtr
409              rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
410              rslope2(i,k) = rslope(i,k)*rslope(i,k)
411              rslope3(i,k) = rslope2(i,k)*rslope(i,k)
412            endif
413          else
414            if(qrs(i,k).le.qcrmin)then
415              rslope(i,k) = rslopesmax
416              rslopeb(i,k) = rslopesbmax
417              rslope2(i,k) = rslopes2max
418              rslope3(i,k) = rslopes3max
419            else
420              rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
421!             rslopeb(i,k) = rslope(i,k)**bvts
422              rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
423              rslope2(i,k) = rslope(i,k)*rslope(i,k)
424              rslope3(i,k) = rslope2(i,k)*rslope(i,k)
425            endif
426          endif
427!-------------------------------------------------------------
428! Ni: ice crystal number concentraiton   [HDC 5c]
429!-------------------------------------------------------------
430!         xni(i,k) = min(max(5.38e7*(den(i,k)                           &
431!                   *max(qci(i,k),qmin))**0.75,1.e3),1.e6)
432          xni(i,k) = min(max(5.38e7*exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
433        enddo
434      enddo
435!
436      mstepmax = 1
437      numdt = 1
438      do k = kte, kts, -1
439        do i = its, ite
440          if(t(i,k).lt.t0c) then
441            pvt = pvts
442          else
443            pvt = pvtr
444          endif
445          work1(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
446          work2(i,k) = work1(i,k)/delz(i,k)
447          numdt(i) = max(nint(work2(i,k)*dtcld+.5),1)
448          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
449        enddo
450      enddo
451      do i = its, ite
452        if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
453      enddo
454!
455      do n = 1, mstepmax
456        k = kte
457        do i = its, ite
458          if(n.le.mstep(i)) then
459            falk(i,k) = den(i,k)*qrs(i,k)*work2(i,k)/mstep(i)
460            hold = falk(i,k)
461            fall(i,k) = fall(i,k)+falk(i,k)
462            holdrs = qrs(i,k)
463            qrs(i,k) = max(qrs(i,k)-falk(i,k)*dtcld/den(i,k),0.)
464          endif
465        enddo
466        do k = kte-1, kts, -1
467          do i = its, ite
468            if(n.le.mstep(i)) then
469              falk(i,k) = den(i,k)*qrs(i,k)*work2(i,k)/mstep(i)
470              hold = falk(i,k)
471              fall(i,k) = fall(i,k)+falk(i,k)
472              holdrs = qrs(i,k)
473              qrs(i,k) = max(qrs(i,k)-(falk(i,k)                        &
474                        -falk(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
475            endif
476          enddo
477        enddo
478      enddo
479!---------------------------------------------------------------
480! Vice [ms-1] : fallout of ice crystal [HDC 5a]
481!---------------------------------------------------------------
482      mstepmax = 1
483      mstep = 1
484      numdt = 1
485      do k = kte, kts, -1
486        do i = its, ite
487          if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
488            xmi = den(i,k)*qci(i,k)/xni(i,k)
489!           diameter  = dicon * sqrt(xmi)
490!           work1c(i,k) = 1.49e4*diameter**1.31
491            diameter  = max(dicon * sqrt(xmi), 1.e-25)
492            work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
493          else
494            work1c(i,k) = 0.
495          endif
496          if(qci(i,k).le.0.) then
497            work2c(i,k) = 0.
498          else
499            work2c(i,k) = work1c(i,k)/delz(i,k)
500          endif
501          numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
502          if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
503        enddo
504      enddo
505      do i = its, ite
506        if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
507      enddo
508!
509      do n = 1, mstepmax
510        k = kte
511        do i = its, ite
512          if (n.le.mstep(i)) then
513            falkc(i,k) = den(i,k)*qci(i,k)*work2c(i,k)/mstep(i)
514            holdc = falkc(i,k)
515            fallc(i,k) = fallc(i,k)+falkc(i,k)
516            holdci = qci(i,k)
517            qci(i,k) = max(qci(i,k)-falkc(i,k)*dtcld/den(i,k),0.)
518          endif
519        enddo
520        do k = kte-1, kts, -1
521          do i = its, ite
522            if (n.le.mstep(i)) then
523              falkc(i,k) = den(i,k)*qci(i,k)*work2c(i,k)/mstep(i)
524              holdc = falkc(i,k)
525              fallc(i,k) = fallc(i,k)+falkc(i,k)
526              holdci = qci(i,k)
527              qci(i,k) = max(qci(i,k)-(falkc(i,k)                       &
528                        -falkc(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
529            endif
530          enddo
531        enddo
532      enddo
533!
534!----------------------------------------------------------------
535!     compute the freezing/melting term. [D89 B16-B17]
536!     freezing occurs one layer above the melting level
537!
538      do i = its, ite
539        mstep(i) = 0
540      enddo
541      do k = kts, kte
542        do i = its, ite
543          if(t(i,k).ge.t0c) then
544            mstep(i) = k
545          endif
546        enddo
547      enddo
548!
549      do i = its, ite
550        kwork2(i) = mstep(i)
551        kwork1(i) = mstep(i)
552        if(mstep(i).ne.0) then
553          if (w(i,mstep(i)).gt.0.) then
554            kwork1(i) = mstep(i) + 1
555          endif
556        endif
557      enddo
558!
559      do i = its, ite
560        k  = kwork1(i)
561        kk = kwork2(i)
562        if(k*kk.ge.1) then
563          qrsci = qrs(i,k) + qci(i,k)
564          if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
565            frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld),    &
566                     qrsci/dtcld)
567            snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld),    &
568                     qrs(i,k)/dtcld)
569            if(k.eq.kk) then
570              t(i,k) = t(i,k) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld
571            else
572              t(i,k) = t(i,k) - xlf0/cpm(i,k)*frzmlt*dtcld
573              t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)*snomlt*dtcld
574            endif
575          endif
576        endif
577      enddo
578!
579!----------------------------------------------------------------
580!      rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
581!
582      do i = its, ite
583        fallsum = fall(i,1)
584        fallsum_qsi = 0.
585        if((t0c-t(i,1)).gt.0) then
586        fallsum = fallsum+fallc(i,1)
587        fallsum_qsi = fall(i,1)+fallc(i,1)
588        endif
589        rainncv(i) = 0.
590        if(fallsum.gt.0.) then
591          rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
592          rain(i) = fallsum*delz(i,1)/denr*dtcld*1000.                &
593                  + rain(i)
594        endif
595        IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
596        snowncv(i) = 0.
597        if(fallsum_qsi.gt.0.) then
598          snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
599          snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
600        endif
601        ENDIF
602        sr(i) = 0.
603        if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000./(rainncv(i)+1.e-12)
604      enddo
605!
606!----------------------------------------------------------------
607!     rsloper: reverse of the slope parameter of the rain(m)
608!     xka:    thermal conductivity of air(jm-1s-1k-1)
609!     work1:  the thermodynamic term in the denominator associated with
610!             heat conduction and vapor diffusion
611!             (ry88, y93, h85)
612!     work2: parameter associated with the ventilation effects(y93)
613!
614      do k = kts, kte
615        do i = its, ite
616          if(t(i,k).ge.t0c) then
617            if(qrs(i,k).le.qcrmin)then
618              rslope(i,k) = rslopermax
619              rslopeb(i,k) = rsloperbmax
620              rslope2(i,k) = rsloper2max
621              rslope3(i,k) = rsloper3max
622            else
623              rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
624              rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
625              rslope2(i,k) = rslope(i,k)*rslope(i,k)
626              rslope3(i,k) = rslope2(i,k)*rslope(i,k)
627            endif
628          else
629            if(qrs(i,k).le.qcrmin)then
630              rslope(i,k) = rslopesmax
631              rslopeb(i,k) = rslopesbmax
632              rslope2(i,k) = rslopes2max
633              rslope3(i,k) = rslopes3max
634            else
635              rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
636              rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
637              rslope2(i,k) = rslope(i,k)*rslope(i,k)
638              rslope3(i,k) = rslope2(i,k)*rslope(i,k)
639            endif
640          endif
641        enddo
642      enddo
643!
644      do k = kts, kte
645        do i = its, ite
646          if(t(i,k).ge.t0c) then
647            work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
648          else
649            work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
650          endif
651          work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
652        enddo
653      enddo
654!
655      do k = kts, kte
656        do i = its, ite
657          supsat = max(q(i,k),qmin)-qs(i,k)
658          satdt = supsat/dtcld
659          if(t(i,k).ge.t0c) then
660!
661!===============================================================
662!
663! warm rain processes
664!
665! - follows the processes in RH83 and LFO except for autoconcersion
666!
667!===============================================================
668!---------------------------------------------------------------
669! praut: auto conversion rate from cloud to rain [HDC 16]
670!        (C->R)
671!---------------------------------------------------------------
672            if(qci(i,k).gt.qc0) then
673!             paut(i,k) = qck1*qci(i,k)**(7./3.)
674              paut(i,k) = qck1*exp(log(qci(i,k))*((7./3.)))
675              paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
676            endif
677!---------------------------------------------------------------
678! pracw: accretion of cloud water by rain [HL A40] [D89 B15]
679!        (C->R)
680!---------------------------------------------------------------
681            if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
682                pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k)          &
683                     *qci(i,k)*denfac(i,k),qci(i,k)/dtcld)
684            endif
685!---------------------------------------------------------------
686! prevp: evaporation/condensation rate of rain [HDC 14]
687!        (V->R or R->V)
688!---------------------------------------------------------------
689            if(qrs(i,k).gt.0.) then
690                coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
691                pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k)            &
692                         +precr2*work2(i,k)*coeres)/work1(i,k)
693              if(pres(i,k).lt.0.) then
694                pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
695                pres(i,k) = max(pres(i,k),satdt/2)
696              else
697                pres(i,k) = min(pres(i,k),satdt/2)
698              endif
699            endif
700          else
701!
702!===============================================================
703!
704! cold rain processes
705!
706! - follows the revised ice microphysics processes in HDC
707! - the processes same as in RH83 and LFO behave
708!   following ice crystal hapits defined in HDC, inclduing
709!   intercept parameter for snow (n0s), ice crystal number
710!   concentration (ni), ice nuclei number concentration
711!   (n0i), ice diameter (d)
712!
713!===============================================================
714!
715            supcol = t0c-t(i,k)
716            ifsat = 0
717!-------------------------------------------------------------
718! Ni: ice crystal number concentraiton   [HDC 5c]
719!-------------------------------------------------------------
720!           xni(i,k) = min(max(5.38e7*(den(i,k)                         &
721!                     *max(qci(i,k),qmin))**0.75,1.e3),1.e6)
722            xni(i,k) = min(max(5.38e7*exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
723            eacrs = exp(0.07*(-supcol))
724!
725            if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
726              xmi = den(i,k)*qci(i,k)/xni(i,k)
727              diameter  = min(dicon * sqrt(xmi),dimax)
728              vt2i = 1.49e4*diameter**1.31
729              vt2s = pvts*rslopeb(i,k)*denfac(i,k)
730!-------------------------------------------------------------
731! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
732!        (T<T0: I->R)
733!-------------------------------------------------------------
734              acrfac = 2.*rslope3(i,k)+2.*diameter*rslope2(i,k)          &
735                      +diameter**2*rslope(i,k)
736              pacr(i,k) = min(pi*qci(i,k)*eacrs*n0s*n0sfac(i,k)          &
737                       *abs(vt2s-vt2i)*acrfac/4.,qci(i,k)/dtcld)
738            endif
739!-------------------------------------------------------------
740! pidep: Deposition/Sublimation rate of ice [HDC 9]
741!       (T<T0: V->I or I->V)
742!-------------------------------------------------------------
743            if(qci(i,k).gt.0.) then
744              xmi = den(i,k)*qci(i,k)/xni(i,k)
745              diameter = dicon * sqrt(xmi)
746              pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.)/work1(i,k)
747              if(pisd(i,k).lt.0.) then
748                pisd(i,k) = max(pisd(i,k),satdt/2)
749                pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
750              else
751                pisd(i,k) = min(pisd(i,k),satdt/2)
752              endif
753              if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
754            endif
755!-------------------------------------------------------------
756! psdep: deposition/sublimation rate of snow [HDC 14]
757!        (V->S or S->V)
758!-------------------------------------------------------------
759            if(qrs(i,k).gt.0..and.ifsat.ne.1) then
760              coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
761              pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k)   &
762                        +precs2*work2(i,k)*coeres)/work1(i,k)
763              supice = satdt-pisd(i,k)
764              if(pres(i,k).lt.0.) then
765                pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
766                pres(i,k) = max(max(pres(i,k),satdt/2),supice)
767              else
768                pres(i,k) = min(min(pres(i,k),satdt/2),supice)
769              endif
770              if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
771            endif
772!-------------------------------------------------------------
773! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
774!       (T<T0: V->I)
775!-------------------------------------------------------------
776            if(supsat.gt.0.and.ifsat.ne.1) then
777              supice = satdt-pisd(i,k)-pres(i,k)
778              xni0 = 1.e3*exp(0.1*supcol)
779!             roqi0 = 4.92e-11*xni0**1.33
780              roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
781              pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld)
782              pgen(i,k) = min(min(pgen(i,k),satdt),supice)
783            endif
784!-------------------------------------------------------------
785! psaut: conversion(aggregation) of ice to snow [HDC 12]
786!       (T<T0: I->S)
787!-------------------------------------------------------------
788            if(qci(i,k).gt.0.) then
789              qimax = roqimax/den(i,k)
790              paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld)
791            endif
792          endif
793        enddo
794      enddo
795!
796!----------------------------------------------------------------
797!     check mass conservation of generation terms and feedback to the
798!     large scale
799!
800      do k = kts, kte
801        do i = its, ite
802          qciik = max(qmin,qci(i,k))
803          delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
804          if(delqci.ge.qciik) then
805            facqci = qciik/delqci
806            paut(i,k) = paut(i,k)*facqci
807            pacr(i,k) = pacr(i,k)*facqci
808            pgen(i,k) = pgen(i,k)*facqci
809            pisd(i,k) = pisd(i,k)*facqci
810          endif
811          qik = max(qmin,q(i,k))
812          delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
813          if(delq.ge.qik) then
814            facq = qik/delq
815            pres(i,k) = pres(i,k)*facq
816            pgen(i,k) = pgen(i,k)*facq
817            pisd(i,k) = pisd(i,k)*facq
818          endif
819          work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
820          q(i,k) = q(i,k)+work2(i,k)*dtcld
821          qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k)     &
822                   -pisd(i,k))*dtcld,0.)
823          qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k)               &
824                   +pres(i,k))*dtcld,0.)
825          if(t(i,k).lt.t0c) then
826            t(i,k) = t(i,k)-xls*work2(i,k)/cpm(i,k)*dtcld
827          else
828            t(i,k) = t(i,k)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld
829          endif
830        enddo
831      enddo
832!
833      cvap = cpv
834      hvap = xlv0
835      hsub = xls
836      ttp=t0c+0.01
837      dldt=cvap-cliq
838      xa=-dldt/rv
839      xb=xa+hvap/(rv*ttp)
840      dldti=cvap-cice
841      xai=-dldti/rv
842      xbi=xai+hsub/(rv*ttp)
843      do k = kts, kte
844        do i = its, ite
845          tr=ttp/t(i,k)
846!         qs(i,k)=psat*(tr**xa)*exp(xb*(1.-tr))
847          qs(i,k)=psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
848          qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
849          qs(i,k) = max(qs(i,k),qmin)
850          denfac(i,k) = sqrt(den0/den(i,k))
851        enddo
852      enddo
853!
854!----------------------------------------------------------------
855!  pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
856!     if there exists additional water vapor condensated/if
857!     evaporation of cloud water is not enough to remove subsaturation
858!
859      do k = kts, kte
860        do i = its, ite
861          work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
862          work2(i,k) = qci(i,k)+work1(i,k)
863          pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
864          if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c)      &
865            pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
866          q(i,k) = q(i,k)-pcon(i,k)*dtcld
867          qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
868          t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
869        enddo
870      enddo
871!
872!----------------------------------------------------------------
873!     padding for small values
874!
875      do k = kts, kte
876        do i = its, ite
877          if(qci(i,k).le.qmin) qci(i,k) = 0.0
878        enddo
879      enddo
880!
881      enddo                  ! big loops
882  END SUBROUTINE wsm32D
883! ...................................................................
884      REAL FUNCTION rgmma(x)
885!-------------------------------------------------------------------
886  IMPLICIT NONE
887!-------------------------------------------------------------------
888!     rgmma function:  use infinite product form
889      REAL :: euler
890      PARAMETER (euler=0.577215664901532)
891      REAL :: x, y
892      INTEGER :: i
893      if(x.eq.1.)then
894        rgmma=0.
895          else
896        rgmma=x*exp(euler*x)
897        do i=1,10000
898          y=float(i)
899          rgmma=rgmma*(1.000+x/y)*exp(-x/y)
900        enddo
901        rgmma=1./rgmma
902      endif
903      END FUNCTION rgmma
904!
905!--------------------------------------------------------------------------
906      REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
907!--------------------------------------------------------------------------
908      IMPLICIT NONE
909!--------------------------------------------------------------------------
910      REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,   &
911           xai,xbi,ttp,tr
912      INTEGER ice
913! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
914      ttp=t0c+0.01
915      dldt=cvap-cliq
916      xa=-dldt/rv
917      xb=xa+hvap/(rv*ttp)
918      dldti=cvap-cice
919      xai=-dldti/rv
920      xbi=xai+hsub/(rv*ttp)
921      tr=ttp/t
922      if(t.lt.ttp.and.ice.eq.1) then
923        fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
924      else
925        fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
926      endif
927! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
928      END FUNCTION fpvs
929!-------------------------------------------------------------------
930  SUBROUTINE wsm3init(den0,denr,dens,cl,cpv,allowed_to_read)
931!-------------------------------------------------------------------
932  IMPLICIT NONE
933!-------------------------------------------------------------------
934!.... constants which may not be tunable
935   REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
936   LOGICAL, INTENT(IN) :: allowed_to_read
937   REAL :: pi
938!
939   pi = 4.*atan(1.)
940   xlv1 = cl-cpv
941!
942   qc0  = 4./3.*pi*denr*r0**3*xncr/den0  ! 0.419e-3 -- .61e-3
943   qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
944!
945   bvtr1 = 1.+bvtr
946   bvtr2 = 2.5+.5*bvtr
947   bvtr3 = 3.+bvtr
948   bvtr4 = 4.+bvtr
949   g1pbr = rgmma(bvtr1)
950   g3pbr = rgmma(bvtr3)
951   g4pbr = rgmma(bvtr4)            ! 17.837825
952   g5pbro2 = rgmma(bvtr2)          ! 1.8273
953   pvtr = avtr*g4pbr/6.
954   eacrr = 1.0
955   pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
956   precr1 = 2.*pi*n0r*.78
957   precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
958   xm0  = (di0/dicon)**2
959   xmmax = (dimax/dicon)**2
960   roqimax = 2.08e22*dimax**8
961!
962   bvts1 = 1.+bvts
963   bvts2 = 2.5+.5*bvts
964   bvts3 = 3.+bvts
965   bvts4 = 4.+bvts
966   g1pbs = rgmma(bvts1)    !.8875
967   g3pbs = rgmma(bvts3)
968   g4pbs = rgmma(bvts4)    ! 12.0786
969   g5pbso2 = rgmma(bvts2)
970   pvts = avts*g4pbs/6.
971   pacrs = pi*n0s*avts*g3pbs*.25
972   precs1 = 4.*n0s*.65
973   precs2 = 4.*n0s*.44*avts**.5*g5pbso2
974   pidn0r =  pi*denr*n0r
975   pidn0s =  pi*dens*n0s
976!
977   rslopermax = 1./lamdarmax
978   rslopesmax = 1./lamdasmax
979   rsloperbmax = rslopermax ** bvtr
980   rslopesbmax = rslopesmax ** bvts
981   rsloper2max = rslopermax * rslopermax
982   rslopes2max = rslopesmax * rslopesmax
983   rsloper3max = rsloper2max * rslopermax
984   rslopes3max = rslopes2max * rslopesmax
985!
986  END SUBROUTINE wsm3init
987END MODULE module_mp_wsm3
Note: See TracBrowser for help on using the repository browser.