wrftladj/module_mp_mkessler.F

   1 MODULE MODULE_MP_MKESSLER
   2   IMPLICIT NONE
   3
   4 CONTAINS
   5 !
   6   SUBROUTINE MKESSLER(t, qv, qc, qr, rho, p, pii, dt_in, z, xlv, cp, ep2&
   7 &    , svp1, svp2, svp3, svpt0, rhowater, dz8w, rainnc, rainncv, ids, ide&
   8 &    , jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
   9 &    jte, kts, kte)
  10     IMPLICIT NONE
  11 !----------------------------------------------------------------
  12 ! Restructered from WRF Kessler Warm rain process
  13 ! H.L. Wang Aug. 1 2009
  14 !----------------------------------------------------------------
  15     REAL, PARAMETER :: c1=.001
  16     REAL, PARAMETER :: c2=.001
  17     REAL, PARAMETER :: c3=2.2
  18     REAL, PARAMETER :: c4=.875
  19 !----------------------------------------------------------------
  20     INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
  21 &    jme, kms, kme, its, ite, jts, jte, kts, kte
  22     REAL, INTENT(IN) :: xlv, cp
  23     REAL, INTENT(IN) :: ep2, svp1, svp2, svp3, svpt0
  24     REAL, INTENT(IN) :: rhowater
  25     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t, qv, &
  26 &    qc, qr
  27     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rho, p, &
  28 &    pii, dz8w
  29     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: z
  30     REAL, INTENT(IN) :: dt_in
  31     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: rainnc, rainncv
  32 ! local variables
  33     REAL :: qrprod, ern, gam, rcgs, rcgsi
  34     REAL, DIMENSION(its:ite, kts:kte, jts:jte) :: prod
  35     REAL, DIMENSION(kts:kte) :: vt, prodk, vtden, rdzk, rhok, piik, &
  36 &    factor, rdzw
  37     INTEGER :: i, j, k
  38     INTEGER :: nfall, n, nfall_new
  39     REAL :: qrr, pressure, temp, es, qvs, dz, dt
  40     REAL :: f5, dtfall, rdz, product
  41     REAL :: vtmax, crmax, factorn
  42     REAL :: qcr, factorr, ppt
  43     REAL, PARAMETER :: max_cr_sedimentation=0.75
  44 !----------------------------------------------------------------
  45     INTEGER :: imax, kmax
  46 ! whl
  47     REAL, DIMENSION(kts:kte) :: qv1d, qc1d, qr1d, t1d, p1d
  48     REAL :: dtleft, rainncv0, max_cr
  49     INTEGER :: kvts, kvte, kn
  50     dt = dt_in
  51 !  print*,'begin'
  52 !   print*,its,ite,jts,jte
  53 !   print*,ims,ime,jms,jme
  54 !   print*,ids,ide,jds,jde
  55     f5 = svp2*(svpt0-svp3)*xlv/cp
  56 !   print*,its,ite,jts,jte
  57 !   print*,ims,ime,jms,jme
  58 !   print*,ids,ide,jds,jde
  59     rdzk = 0.0
  60     rdzw = 0.0
  61     DO j=jts,jte
  62       DO i=its,ite
  63         DO k=1,kte-1
  64           rdzk(k) = 1./(z(i, k+1, j)-z(i, k, j))
  65         END DO
  66         rdzk(kte) = 1./(z(i, kte, j)-z(i, kte-1, j))
  67       END DO
  68     END DO
  69     DO j=jts,jte
  70       DO i=its,ite
  71         DO k=1,kte
  72           qv1d(k) = qv(i, k, j)
  73           qc1d(k) = qc(i, k, j)
  74           qr1d(k) = qr(i, k, j)
  75           t1d(k) = t(i, k, j)
  76           p1d(k) = p(i, k, j)
  77           rhok(k) = rho(i, k, j)
  78           piik(k) = pii(i, k, j)
  79           rdzw(k) = 1./dz8w(i, k, j)
  80         END DO
  81 !   print*,i,j
  82         kvts = kts
  83         kvte = kte
  84         max_cr = max_cr_sedimentation
  85         dtleft = dt
  86         CALL SMALLSTEP(qr1d, rdzk, rdzw, rhok, max_cr, dtleft, nfall, &
  87 &                 kvts, kvte)
  88         dtleft = dt/nfall
  89         rainncv0 = 0.0
  90         rainncv(i, j) = 0.0
  91         DO kn=1,nfall
  92           CALL RFALL(qr1d, rdzk, rdzw, rhok, rainncv0, rhowater, max_cr&
  93 &               , dtleft, kvts, kvte)
  94           rainncv(i, j) = rainncv(i, j) + rainncv0
  95         END DO
  96 !    print*,rainncv0
  97 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
  98 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
  99         rainnc(i, j) = rainnc(i, j) + rainncv(i, j)
 100 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
 101         CALL AUTOCA(qc1d, qr1d, kvts, kvte, c1, c2, c3, c4, dt)
 102 !satadj(qv,qc,qr, tmp, pii,rho,  kvts,kvte,xlv, cp,EP2,SVP1,SVP2,SVP3,SVPT0)
 103         CALL SATADJ(qv1d, qc1d, qr1d, t1d, p1d, piik, rhok, kvts, kvte, &
 104 &              xlv, dt, cp, ep2, svp1, svp2, svp3, svpt0)
 105         DO k=1,kte
 106           qv(i, k, j) = qv1d(k)
 107           qc(i, k, j) = qc1d(k)
 108           qr(i, k, j) = qr1d(k)
 109           t(i, k, j) = t1d(k)
 110         END DO
 111       END DO
 112     END DO
 113 ! print*,rainncv
 114     RETURN
 115   END SUBROUTINE MKESSLER
 116
 117   SUBROUTINE SMALLSTEP(prodk, rdzk, rdzw, rhok, max_cr, dtleft, nstep, &
 118 &    kvts, kvte)
 119     IMPLICIT NONE
 120     INTEGER :: nstep, k, kvts, kvte
 121     REAL, DIMENSION(kvts:kvte) :: vtden, vt, prodk, factor, rdzk, rdzw, &
 122 &    rhok
 123     REAL :: max_cr, ppt, dtleft, crmax, qrr
 124     REAL :: arg1
 125     INTRINSIC AMAX1
 126 !    INTRINSIC NINT
 127     INTRINSIC SQRT
 128     INTRINSIC NINT
 129     crmax = 0.0
 130     DO k=kvts,kvte-1
 131       qrr = prodk(k)*0.001*rhok(k)
 132       arg1 = rhok(1)/rhok(k)
 133       vtden(k) = SQRT(arg1)
 134       IF (qrr/(0.001*rhok(k)) .GE. 1d-5) THEN
 135         vt(k) = 36.34*qrr**0.1364*vtden(k)
 136       ELSE
 137         vt(k) = 0.0
 138       END IF
 139       IF (vt(k)*dtleft*rdzw(k) .LT. crmax) THEN
 140         crmax = crmax
 141       ELSE
 142         crmax = vt(k)*dtleft*rdzw(k)
 143       END IF
 144     END DO
 145 !    nstep = NINT(0.5 + crmax/max_cr)
 146     nstep = NINT(0.5 + crmax/0.75)
 147   END SUBROUTINE SMALLSTEP
 148
 149   SUBROUTINE RFALL(prodk, rdzk, rdzw, rhok, rainncv0, rhowat, max_cr, &
 150 &    dtfall, kvts, kvte)
 151     IMPLICIT NONE
 152     INTEGER :: k, kvts, kvte
 153     REAL, DIMENSION(kvts:kvte) :: vtden, vt, prodk, factor, rdzk, rdzw, &
 154 &    rhok
 155     REAL :: rainncv0, rhowat, max_cr, ppt, dtleft
 156     REAL :: qrr, dtfall
 157     REAL :: arg1
 158     INTRINSIC SQRT
 159     DO k=kvts,kvte
 160       IF (prodk(k) .LT. 0) prodk(k) = 0.0
 161     END DO
 162     DO k=kvts,kvte
 163       qrr = prodk(k)*0.001*rhok(k)
 164       arg1 = rhok(1)/rhok(k)
 165       vtden(k) = SQRT(arg1)
 166       IF (qrr/( 0.001*rhok(k) ) .GE. 1d-5) THEN
 167         vt(k) = 36.34*qrr**0.1364*vtden(k)
 168       ELSE
 169         vt(k) = 0.0
 170       END IF
 171     END DO
 172     DO k=kvts,kvte-1
 173       factor(k) = dtfall*rdzk(k)/rhok(k)
 174     END DO
 175     factor(kvte) = dtfall*rdzk(kvte)
 176     ppt = 0.
 177     k = 1
 178     ppt = rhok(k)*prodk(k)*vt(k)*dtfall/rhowat
 179 !mm
 180     rainncv0 = ppt*1000.
 181 !      print*,rainncv0
 182 !------------------------------------------------------------------------------
 183 ! Time split loop, Fallout done with flux upstream
 184 !------------------------------------------------------------------------------
 185     DO k=kvts,kvte-1
 186       prodk(k) = prodk(k) - factor(k)*(rhok(k)*prodk(k)*vt(k)-rhok(k+1)*&
 187 &        prodk(k+1)*vt(k+1))
 188     END DO
 189     k = kvte
 190     prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
 191     DO k=kvts,kvte
 192       IF (prodk(k) .LT. 0) prodk(k) = 0.0
 193     END DO
 194   END SUBROUTINE RFALL
 195
 196   SUBROUTINE AUTOCA(qc1d, qr1d, kvts, kvte, c1, c2, c3, c4, dt)
 197     IMPLICIT NONE
 198 !     print*,k,qrprod
 199     INTEGER :: kvts, kvte, k
 200     REAL, DIMENSION(kvts:kvte) :: qc1d, qr1d
 201     REAL :: c1, c2, c3, c4
 202     REAL :: qrrc, dt, factorn, qrprod, qrprod2
 203     REAL :: pwr1
 204     qrrc = 1.0e-5
 205     DO k=kvts,kvte
 206       IF (qr1d(k) .LT. 0.0) qr1d(k) = 0.0
 207       IF (qc1d(k) .LT. 0.0) qc1d(k) = 0.0
 208       IF (qr1d(k) .GE. qrrc) THEN
 209         pwr1 = qr1d(k)**c4
 210         factorn = 1.0/(1.+c3*dt*pwr1)
 211       ELSE
 212         factorn = 1.0
 213       END IF
 214       qrprod = qc1d(k)*(1.0-factorn)
 215       qrprod2 = 0.0
 216       IF (qc1d(k) - c2 .GT. 0) THEN
 217         qrprod2 = factorn*c1*dt*(qc1d(k)-c2)
 218         IF (qrprod2 .GT. qc1d(k) - c2) qrprod2 = qc1d(k) - c2
 219       END IF
 220 !        print*,k,qrprod2
 221       qrprod = qrprod + qrprod2
 222       IF (qc1d(k) - qrprod .GT. 0) THEN
 223         qc1d(k) = qc1d(k) - qrprod
 224         qr1d(k) = qr1d(k) + qrprod
 225       ELSE
 226         qc1d(k) = 0.0
 227         qrprod = qc1d(k)
 228         qr1d(k) = qr1d(k) + qrprod
 229       END IF
 230     END DO
 231   END SUBROUTINE AUTOCA
 232
 233   SUBROUTINE SATADJ(qv, qc, qr, tmp, p1d, pii, rhok, kvts, kvte, xlv, dt&
 234 &    , cp, ep2, svp1, svp2, svp3, svpt0)
 235     IMPLICIT NONE
 236     INTEGER :: kvts, kvte, k
 237     REAL, DIMENSION(kvts:kvte) :: qv, qc, qr, tmp, p1d, pii, rhok
 238     REAL, DIMENSION(kvts:kvte) :: rcgs, pressure, temp, es, qvs
 239     REAL, DIMENSION(kvts:kvte) :: ern, qv2cl, rn2qv
 240 ! local var
 241     REAL :: svp1, svp2, svp3, svpt0, ep2, xlv, cp, dt, f5
 242     REAL :: ernmax, product
 243     REAL :: arg1
 244     INTRINSIC EXP
 245     f5 = svp2*(svpt0-svp3)*xlv/cp
 246     DO k=kvts,kvte
 247 !constant
 248       rcgs(k) = 0.001*rhok(k)
 249       pressure(k) = p1d(k)
 250       temp(k) = pii(k)*tmp(k)
 251       arg1 = svp2*(temp(k)-svpt0)/(temp(k)-svp3)
 252       es(k) = 1000.*svp1*EXP(arg1)
 253       qvs(k) = ep2*es(k)/(pressure(k)-es(k))
 254       IF (qr(k) .LT. 0) qr(k) = 0.0
 255       IF (qv(k) .LT. 0) qv(k) = 0.0
 256       IF (qc(k) .LT. 0) qc(k) = 0.0
 257     END DO
 258     DO k=kvts,kvte
 259 !not related to time; maximum transform qv to cl (sat) or cl to qv (sub sat)
 260       qv2cl(k) = (qv(k)-qvs(k))/(1.+pressure(k)/(pressure(k)-es(k))*qvs(&
 261 &        k)*f5/(temp(k)-svp3)**2)
 262 ! sub sat rain evaperate
 263       rn2qv(k) = 0.0
 264       ern(k) = 0.0
 265       IF (qvs(k) .GT. qv(k)) THEN
 266         IF (qr(k) .GE. 1d-5) THEN
 267           rn2qv(k) = dt*((1.6+124.9*(rcgs(k)*qr(k))**.2046)*(rcgs(k)*qr(&
 268 &            k))**.525/(2.55e8/(pressure(k)*qvs(k))+5.4e5))*((qvs(k)-qv(k&
 269 &            ))/(rcgs(k)*qvs(k)))
 270         ELSE
 271           rn2qv(k) = 0.0
 272         END IF
 273         IF (rn2qv(k) .GT. qr(k)) rn2qv(k) = qr(k)
 274         ernmax = 0.0
 275         IF (-qv2cl(k) - qc(k) .GT. 0.0) ernmax = -qv2cl(k) - qc(k)
 276 !        ern(k)  = amin1(rn2qv(k), ernmax)
 277         ern(k) = rn2qv(k)
 278         IF (rn2qv(k) .GT. ernmax) ern(k) = ernmax
 279       END IF
 280 ! Update all variables
 281 !       product = amax1(qv2cl(k),-qc(k))
 282       product = qv2cl(k)
 283       IF (qv2cl(k) .LT. -qc(k)) product = -qc(k)
 284 !       qv(k) = amax1(qv(k) - product + ern(k),0.)
 285       qv(k) = qv(k) - product + ern(k)
 286       IF (qv(k) .LT. 0) qv(k) = 0.0
 287       qc(k) = qc(k) + product
 288       qr(k) = qr(k) - ern(k)
 289       temp(k) = temp(k) + xlv/cp*(product-ern(k))
 290       tmp(k) = temp(k)/pii(k)
 291     END DO
 292   END SUBROUTINE SATADJ
 293
 294 END MODULE MODULE_MP_MKESSLER