Update version info for release v4.6.1 (#2122)
[WRF.git] / phys / module_mp_kessler.F
blob7dd2e36894f8dd4ab308f8b4cb6f53a003c86bf0
1 !WRF:MODEL_LAYER:PHYSICS
4 MODULE module_mp_kessler
6 CONTAINS
7 !----------------------------------------------------------------
8    SUBROUTINE kessler( t, qv, qc, qr, rho, pii                  &
9                       ,dt_in, z, xlv, cp                        &
10                       ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater        &
11                       ,dz8w                                     &
12                       ,RAINNC, RAINNCV                          &
13                       ,ids,ide, jds,jde, kds,kde                & ! domain dims
14                       ,ims,ime, jms,jme, kms,kme                & ! memory dims
15                       ,its,ite, jts,jte, kts,kte                & ! tile   dims
16                                                                 )
17 !----------------------------------------------------------------
18    IMPLICIT NONE
19 !----------------------------------------------------------------
20    !  taken from the COMMAS code - WCS 10 May 1999.
21    !  converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
22 !----------------------------------------------------------------
23    REAL    , PARAMETER ::  c1 = .001 
24    REAL    , PARAMETER ::  c2 = .001 
25    REAL    , PARAMETER ::  c3 = 2.2 
26    REAL    , PARAMETER ::  c4 = .875 
27    REAL    , PARAMETER ::  fudge = 1.0 
28    REAL    , PARAMETER ::  mxfall = 10.0 
29 !----------------------------------------------------------------
30    INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde, &
31                                      ims,ime, jms,jme, kms,kme, &
32                                      its,ite, jts,jte, kts,kte
33    REAL   ,      INTENT(IN   )    :: xlv, cp
34    REAL   ,      INTENT(IN   )    :: EP2,SVP1,SVP2,SVP3,SVPT0
35    REAL   ,      INTENT(IN   )    :: rhowater
37    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
38          INTENT(INOUT) ::                                       &
39                                                             t , &
40                                                             qv, &
41                                                             qc, &
42                                                             qr
44    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
45          INTENT(IN   ) ::                                       &
46                                                            rho, &
47                                                            pii, &
48                                                           dz8w 
50    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
51          INTENT(IN   ) ::                                    z
53    REAL, INTENT(IN   ) :: dt_in
55    REAL, DIMENSION( ims:ime , jms:jme ),                        &
56          INTENT(INOUT) ::                               RAINNC, &
57                                                        RAINNCV
61    ! local variables
63    REAL :: qrprod, ern, gam, rcgs, rcgsi
64    REAL, DIMENSION( its:ite , kts:kte, jts:jte ) ::     prod
65    REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
66    INTEGER :: i,j,k
67    INTEGER :: nfall, n, nfall_new
68    REAL    :: qrr, pressure, temp, es, qvs, dz, dt
69    REAL    :: f5, dtfall, rdz, product
70    REAL    :: max_heating, max_condense, max_rain, maxqrp
71    REAL    :: vtmax, ernmax, crmax, factorn, time_sediment
72    REAL    :: qcr, factorr, ppt
73    REAL, PARAMETER :: max_cr_sedimentation = 0.75
74 !----------------------------------------------------------------
76    INTEGER :: imax, kmax
78     dt = dt_in
80 !   f5 = 237.3 * 17.27 * 2.5e6 / cp 
81     f5 = svp2*(svpt0-svp3)*xlv/cp
82     ernmax = 0.
83     maxqrp = -100.
85 !------------------------------------------------------------------------------
86 ! parameters for the time split terminal advection
87 !------------------------------------------------------------------------------
89       max_heating = 0.
90       max_condense = 0.
91       max_rain = 0.
93 !-----------------------------------------------------------------------------
94 ! outer J loop for entire microphysics, outer i loop for sedimentation
95 !-----------------------------------------------------------------------------
97   microphysics_outer_j_loop: DO j = jts, jte
99   sedimentation_outer_i_loop: DO i = its,ite
101 !   vtmax = 0.
102    crmax = 0.
105 !------------------------------------------------------------------------------
106 ! Terminal velocity calculation and advection, set up coefficients and
107 ! compute stable timestep
108 !------------------------------------------------------------------------------
110    DO k = 1, kte
111      prodk(k)   = qr(i,k,j)
112      rhok(k) = rho(i,k,j)
113      qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k))
114      vtden(k) = sqrt(rhok(1)/rhok(k))
115      vt(k) = 36.34*(qrr**0.1364) * vtden(k)
116 !     vtmax = amax1(vt(k), vtmax)
117      rdzw(k) = 1./dz8w(i,k,j)
118      crmax = amax1(vt(k)*dt*rdzw(k),crmax)
119    ENDDO
120    DO k = 1, kte-1
121      rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
122    ENDDO
123    rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))
125    nfall = max(1,nint(0.5+crmax/max_cr_sedimentation))  ! courant number for big timestep.
126    dtfall = dt / float(nfall)                           ! splitting so courant number for sedimentation
127    time_sediment = dt                                   ! is stable
129 !------------------------------------------------------------------------------
130 ! Terminal velocity calculation and advection
131 ! Do a time split loop on this for stability.
132 !------------------------------------------------------------------------------
134    column_sedimentation: DO WHILE ( nfall > 0 )
136    time_sediment = time_sediment - dtfall
137    DO k = 1, kte-1
138      factor(k) = dtfall*rdzk(k)/rhok(k)
139    ENDDO
140    factor(kte) = dtfall*rdzk(kte)
142    ppt=0.
144       k = 1
145       ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater
146       RAINNCV(i,j)=ppt*1000.
147       RAINNC(i,j)=RAINNC(i,j)+ppt*1000.  ! unit = mm
149 !------------------------------------------------------------------------------
150 ! Time split loop, Fallout done with flux upstream
151 !------------------------------------------------------------------------------
153       DO k = kts, kte-1
154         prodk(k) = prodk(k) - factor(k)           &
155                   * (rhok(k)*prodk(k)*vt(k)       &
156                     -rhok(k+1)*prodk(k+1)*vt(k+1))
157       ENDDO
159       k = kte
160       prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
162 !------------------------------------------------------------------------------
163 ! compute new sedimentation velocity, and check/recompute new 
164 ! sedimentation timestep if this isn't the last split step.
165 !------------------------------------------------------------------------------
167       IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep
169         nfall = nfall - 1
170         crmax = 0.
171         DO k = kts, kte 
172           qrr = amax1(0.,prodk(k)*0.001*rhok(k))
173           vt(k) = 36.34*(qrr**0.1364) * vtden(k)
174 !          vtmax = amax1(vt(k), vtmax)
175           crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax)
176         ENDDO
178         nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
179         if (nfall_new /= nfall ) then
180           nfall = nfall_new
181           dtfall = time_sediment/nfall
182         end if
184       ELSE  ! this was the last timestep
186         DO k=kts,kte
187           prod(i,k,j) = prodk(k)
188         ENDDO
189         nfall = 0  ! exit condition for sedimentation loop
191       END IF
193    ENDDO column_sedimentation
195    ENDDO sedimentation_outer_i_loop
197 !------------------------------------------------------------------------------
198 ! Production of rain and deletion of qc
199 ! Production of qc from supersaturation
200 ! Evaporation of QR
201 !------------------------------------------------------------------------------
203      DO k = kts, kte
204      DO i = its, ite
205        factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4)
206        qrprod = qc(i,k,j) * (1.0 - factorn)           &
207              + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.)      
208        rcgs = 0.001*rho(i,k,j)
210        qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
211        qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
212        qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)
214        temp      = pii(i,k,j)*t(i,k,j)
215        pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
216        gam = 2.5e+06/(1004.*pii(i,k,j))
217 !      qvs       = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure
218        es        = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3))
219        qvs       = ep2*es/(pressure-es)
220 !      prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2)
221        prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2)
222        ern  = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046)   &
223           *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs)       &
224           +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)),             &
225           amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))
227 ! Update all variables
229        product = amax1(prod(i,k,j),-qc(i,k,j))
230        t (i,k,j) = t(i,k,j) + gam*(product - ern)
231        qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
232        qc(i,k,j) =       qc(i,k,j) + product
233        qr(i,k,j) = qr(i,k,j) - ern
235      ENDDO
236      ENDDO
238   ENDDO  microphysics_outer_j_loop
240   RETURN
242   END SUBROUTINE kessler
244 END MODULE module_mp_kessler