wrftladj/module_mp_mkessler_tl.F

   1 !        Generated by TAPENADE     (INRIA, Tropics team)
   2 !  Tapenade 3.3 (r3163) - 09/25/2009 09:04
   3 !
   4 MODULE G_MODULE_MP_MKESSLER
   5   IMPLICIT NONE
   6
   7 CONTAINS
   8
   9   SUBROUTINE G_MKESSLER(t, td, qv, qvd, qc, qcd, qr, qrd, rho, rhod, p, &
  10 &    pd, pii, piid, dt_in, z, xlv, cp, ep2, svp1, svp2, svp3, svpt0, &
  11 &    rhowater, dz8w, rainnc, rainncd, rainncv, rainncvd, ids, ide, jds, &
  12 &    jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts&
  13 &    , kte)
  14
  15     USE module_mp_mkessler, ONLY : SMALLSTEP
  16
  17     IMPLICIT NONE
  18 !----------------------------------------------------------------
  19 ! Restructered from WRF Kessler Warm rain process
  20 ! H.L. Wang Aug. 1 2009
  21 !----------------------------------------------------------------
  22     REAL, PARAMETER :: c1=.001
  23     REAL, PARAMETER :: c2=.001
  24     REAL, PARAMETER :: c3=2.2
  25     REAL, PARAMETER :: c4=.875
  26 !----------------------------------------------------------------
  27     INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
  28 &    jme, kms, kme, its, ite, jts, jte, kts, kte
  29     REAL, INTENT(IN) :: xlv, cp
  30     REAL, INTENT(IN) :: ep2, svp1, svp2, svp3, svpt0
  31     REAL, INTENT(IN) :: rhowater
  32     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t, qv, &
  33 &    qc, qr
  34     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: td, qvd&
  35 &    , qcd, qrd
  36     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rho, p, &
  37 &    pii, dz8w
  38     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rhod, pd, &
  39 &    piid
  40     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: z
  41     REAL, INTENT(IN) :: dt_in
  42     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: rainnc, rainncv
  43     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: rainncd, &
  44 &    rainncvd
  45 ! local variables
  46     REAL :: qrprod, ern, gam, rcgs, rcgsi
  47     REAL, DIMENSION(its:ite, kts:kte, jts:jte) :: prod
  48     REAL, DIMENSION(kts:kte) :: vt, prodk, vtden, rdzk, rhok, piik, &
  49 &    factor, rdzw
  50     REAL, DIMENSION(kts:kte) :: rdzkd, rhokd, piikd
  51     INTEGER :: i, j, k
  52     INTEGER :: nfall, n, nfall_new
  53     REAL :: qrr, pressure, temp, es, qvs, dz, dt
  54     REAL :: f5, dtfall, rdz, product
  55     REAL :: vtmax, crmax, factorn
  56     REAL :: qcr, factorr, ppt
  57     REAL, PARAMETER :: max_cr_sedimentation=0.75
  58 !----------------------------------------------------------------
  59     INTEGER :: imax, kmax
  60 ! whl
  61     REAL, DIMENSION(kts:kte) :: qv1d, qc1d, qr1d, t1d, p1d
  62     REAL, DIMENSION(kts:kte) :: qv1dd, qc1dd, qr1dd, t1dd, p1dd
  63     REAL :: dtleft, rainncv0, max_cr
  64     REAL :: rainncv0d
  65     INTEGER :: kvts, kvte, kn
  66     dt = dt_in
  67     f5 = svp2*(svpt0-svp3)*xlv/cp
  68     rdzk = 0.0
  69     rdzw = 0.0
  70     DO j=jts,jte
  71       DO i=its,ite
  72         DO k=1,kte-1
  73           rdzkd(k) = 0.0
  74           rdzk(k) = 1./(z(i, k+1, j)-z(i, k, j))
  75         END DO
  76         rdzkd(kte) = 0.0
  77         rdzk(kte) = 1./(z(i, kte, j)-z(i, kte-1, j))
  78       END DO
  79     END DO
  80     qv1dd = 0.0
  81     qc1dd = 0.0
  82     p1dd = 0.0
  83     qr1dd = 0.0
  84     piikd = 0.0
  85     rhokd = 0.0
  86     t1dd = 0.0
  87     DO j=jts,jte
  88       DO i=its,ite
  89         DO k=1,kte
  90           qv1dd(k) = qvd(i, k, j)
  91           qv1d(k) = qv(i, k, j)
  92           qc1dd(k) = qcd(i, k, j)
  93           qc1d(k) = qc(i, k, j)
  94           qr1dd(k) = qrd(i, k, j)
  95           qr1d(k) = qr(i, k, j)
  96           t1dd(k) = td(i, k, j)
  97           t1d(k) = t(i, k, j)
  98           p1dd(k) = pd(i, k, j)
  99           p1d(k) = p(i, k, j)
 100           rhokd(k) = rhod(i, k, j)
 101           rhok(k) = rho(i, k, j)
 102           piikd(k) = piid(i, k, j)
 103           piik(k) = pii(i, k, j)
 104           rdzw(k) = 1./dz8w(i, k, j)
 105         END DO
 106 !   print*,i,j
 107         kvts = kts
 108         kvte = kte
 109         max_cr = max_cr_sedimentation
 110         dtleft = dt
 111         CALL SMALLSTEP(qr1d, rdzk, rdzw, rhok, max_cr, dtleft, nfall, &
 112 &                 kvts, kvte)
 113         dtleft = dt/nfall
 114         rainncv0 = 0.0
 115         rainncvd(i, j) = 0.0
 116         rainncv(i, j) = 0.0
 117         DO kn=1,nfall
 118           CALL RFALL_D(qr1d, qr1dd, rdzk, rdzw, rhok, rhokd, rainncv0, &
 119 &                 rainncv0d, rhowater, max_cr, dtleft, kvts, kvte)
 120           rainncvd(i, j) = rainncvd(i, j) + rainncv0d
 121           rainncv(i, j) = rainncv(i, j) + rainncv0
 122         END DO
 123 !    print*,rainncv0
 124 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
 125         rainncd(i, j) = rainncd(i, j) + rainncvd(i, j)
 126         rainnc(i, j) = rainnc(i, j) + rainncv(i, j)
 127 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
 128         CALL AUTOCA_D(qc1d, qc1dd, qr1d, qr1dd, kvts, kvte, c1, c2, c3, &
 129 &                c4, dt)
 130 !satadj(qv,qc,qr, tmp, pii,rho,  kvts,kvte,xlv, cp,EP2,SVP1,SVP2,SVP3,SVPT0)
 131         CALL SATADJ_D(qv1d, qv1dd, qc1d, qc1dd, qr1d, qr1dd, t1d, t1dd, &
 132 &                p1d, p1dd, piik, piikd, rhok, rhokd, kvts, kvte, xlv, dt&
 133 &                , cp, ep2, svp1, svp2, svp3, svpt0)
 134         DO k=1,kte
 135           qvd(i, k, j) = qv1dd(k)
 136           qv(i, k, j) = qv1d(k)
 137           qcd(i, k, j) = qc1dd(k)
 138           qc(i, k, j) = qc1d(k)
 139           qrd(i, k, j) = qr1dd(k)
 140           qr(i, k, j) = qr1d(k)
 141           td(i, k, j) = t1dd(k)
 142           t(i, k, j) = t1d(k)
 143         END DO
 144       END DO
 145     END DO
 146     RETURN
 147   END SUBROUTINE G_MKESSLER
 148 !
 149 !  Differentiation of rfall in forward (tangent) mode:
 150 !   variations  of output variables: prodk rainncv0
 151 !   with respect to input variables: prodk rhok
 152   SUBROUTINE RFALL_D(prodk, prodkd, rdzk, rdzw, rhok, rhokd, rainncv0, &
 153 &    rainncv0d, rhowat, max_cr, dtfall, kvts, kvte)
 154     IMPLICIT NONE
 155     INTEGER :: k, kvts, kvte
 156     REAL, DIMENSION(kvts:kvte) :: vtden, vt, prodk, factor, rdzk, rdzw, &
 157 &    rhok
 158     REAL, DIMENSION(kvts:kvte) :: vtdend, vtd, prodkd, factord, rhokd
 159     REAL :: rainncv0, rhowat, max_cr, ppt, dtleft
 160     REAL :: rainncv0d, pptd
 161     REAL :: qrr, dtfall
 162     REAL :: qrrd
 163     REAL :: arg1
 164     REAL :: arg1d
 165     INTRINSIC SQRT
 166     DO k=kvts,kvte
 167       IF (prodk(k) .LT. 0) THEN
 168         prodkd(k) = 0.0
 169         prodk(k) = 0.0
 170       END IF
 171     END DO
 172     vtd = 0.0
 173     vtdend = 0.0
 174     DO k=kvts,kvte
 175       qrrd = 0.001*(prodkd(k)*rhok(k)+prodk(k)*rhokd(k))
 176       qrr = prodk(k)*0.001*rhok(k)
 177       arg1d = (rhokd(1)*rhok(k)-rhok(1)*rhokd(k))/rhok(k)**2
 178       arg1 = rhok(1)/rhok(k)
 179       IF (arg1 .EQ. 0.0) THEN
 180         vtdend(k) = 0.0
 181       ELSE
 182         vtdend(k) = arg1d/(2.0*SQRT(arg1))
 183       END IF
 184       vtden(k) = SQRT(arg1)
 185       IF (qrr/(0.001*rhok(k)) .GE. 1d-5) THEN
 186         vtd(k) = 36.34*(0.1364*qrr**(-0.8636)*qrrd*vtden(k)+qrr**0.1364*&
 187 &          vtdend(k))
 188         vt(k) = 36.34*qrr**0.1364*vtden(k)
 189       ELSE
 190         vtd(k) = 0.0
 191         vt(k) = 0.0
 192       END IF
 193     END DO
 194     factord = 0.0
 195 !     pause
 196     DO k=kvts,kvte-1
 197       factord(k) = -(dtfall*rdzk(k)*rhokd(k)/rhok(k)**2)
 198       factor(k) = dtfall*rdzk(k)/rhok(k)
 199     END DO
 200     factord(kvte) = 0.0
 201     factor(kvte) = dtfall*rdzk(kvte)
 202     ppt = 0.
 203     k = 1
 204     pptd = dtfall*((rhokd(k)*prodk(k)+rhok(k)*prodkd(k))*vt(k)+rhok(k)*&
 205 &      prodk(k)*vtd(k))/rhowat
 206     ppt = rhok(k)*prodk(k)*vt(k)*dtfall/rhowat
 207 !mm
 208     rainncv0d = 1000.*pptd
 209     rainncv0 = ppt*1000.
 210 !      print*,rainncv0
 211 !------------------------------------------------------------------------------
 212 ! Time split loop, Fallout done with flux upstream
 213 !------------------------------------------------------------------------------
 214     DO k=kvts,kvte-1
 215       prodkd(k) = prodkd(k) - factord(k)*(rhok(k)*prodk(k)*vt(k)-rhok(k+&
 216 &        1)*prodk(k+1)*vt(k+1)) - factor(k)*((rhokd(k)*prodk(k)+rhok(k)*&
 217 &        prodkd(k))*vt(k)+rhok(k)*prodk(k)*vtd(k)-(rhokd(k+1)*prodk(k+1)+&
 218 &        rhok(k+1)*prodkd(k+1))*vt(k+1)-rhok(k+1)*prodk(k+1)*vtd(k+1))
 219       prodk(k) = prodk(k) - factor(k)*(rhok(k)*prodk(k)*vt(k)-rhok(k+1)*&
 220 &        prodk(k+1)*vt(k+1))
 221     END DO
 222     k = kvte
 223     prodkd(k) = prodkd(k) - (factord(k)*prodk(k)+factor(k)*prodkd(k))*vt&
 224 &      (k) - factor(k)*prodk(k)*vtd(k)
 225     prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
 226     DO k=kvts,kvte
 227       IF (prodk(k) .LT. 0) THEN
 228         prodkd(k) = 0.0
 229         prodk(k) = 0.0
 230       END IF
 231     END DO
 232   END SUBROUTINE RFALL_D
 233 !   with respect to input variables: qc1d qr1d
 234   SUBROUTINE AUTOCA_D(qc1d, qc1dd, qr1d, qr1dd, kvts, kvte, c1, c2, c3, &
 235 &    c4, dt)
 236     IMPLICIT NONE
 237 !     print*,k,qrprod
 238     INTEGER :: kvts, kvte, k
 239     REAL, DIMENSION(kvts:kvte) :: qc1d, qr1d
 240     REAL, DIMENSION(kvts:kvte) :: qc1dd, qr1dd
 241     REAL :: c1, c2, c3, c4
 242     REAL :: qrrc, dt, factorn, qrprod, qrprod2
 243     REAL :: factornd, qrprodd, qrprod2d
 244     REAL :: pwr1
 245     REAL :: pwr1d
 246     qrrc = 1.0e-5
 247     DO k=kvts,kvte
 248       IF (qr1d(k) .LT. 0.0) THEN
 249         qr1dd(k) = 0.0
 250         qr1d(k) = 0.0
 251       END IF
 252       IF (qc1d(k) .LT. 0.0) THEN
 253         qc1dd(k) = 0.0
 254         qc1d(k) = 0.0
 255       END IF
 256       IF (qr1d(k) .GE. qrrc) THEN
 257         IF (qr1d(k) .GT. 0.0 .OR. (qr1d(k) .LT. 0.0 .AND. c4 .EQ. INT(c4&
 258 &            ))) THEN
 259           pwr1d = c4*qr1d(k)**(c4-1)*qr1dd(k)
 260         ELSE IF (qr1d(k) .EQ. 0.0 .AND. c4 .EQ. 1.0) THEN
 261           pwr1d = qr1dd(k)
 262         ELSE
 263           pwr1d = 0.0
 264         END IF
 265         pwr1 = qr1d(k)**c4
 266         factornd = -(c3*dt*pwr1d/(1.+c3*dt*pwr1)**2)
 267         factorn = 1.0/(1.+c3*dt*pwr1)
 268       ELSE
 269         factorn = 1.0
 270         factornd = 0.0
 271       END IF
 272       qrprodd = qc1dd(k)*(1.0-factorn) - qc1d(k)*factornd
 273       qrprod = qc1d(k)*(1.0-factorn)
 274       qrprod2 = 0.0
 275       IF (qc1d(k) - c2 .GT. 0) THEN
 276         qrprod2d = c1*dt*(factornd*(qc1d(k)-c2)+factorn*qc1dd(k))
 277         qrprod2 = factorn*c1*dt*(qc1d(k)-c2)
 278         IF (qrprod2 .GT. qc1d(k) - c2) THEN
 279           qrprod2d = qc1dd(k)
 280           qrprod2 = qc1d(k) - c2
 281         END IF
 282       ELSE
 283         qrprod2d = 0.0
 284       END IF
 285 !        print*,k,qrprod2
 286       qrprodd = qrprodd + qrprod2d
 287       qrprod = qrprod + qrprod2
 288       IF (qc1d(k) - qrprod .GT. 0) THEN
 289         qc1dd(k) = qc1dd(k) - qrprodd
 290         qc1d(k) = qc1d(k) - qrprod
 291         qr1dd(k) = qr1dd(k) + qrprodd
 292         qr1d(k) = qr1d(k) + qrprod
 293       ELSE
 294         qc1dd(k) = 0.0
 295         qc1d(k) = 0.0
 296         qrprodd = qc1dd(k)
 297         qrprod = qc1d(k)
 298         qr1dd(k) = qr1dd(k) + qrprodd
 299         qr1d(k) = qr1d(k) + qrprod
 300       END IF
 301     END DO
 302   END SUBROUTINE AUTOCA_D
 303 !  Differentiation of satadj in forward (tangent) mode:
 304 !   variations  of output variables: qc qr qv tmp
 305 !   with respect to input variables: qc qr qv p1d rhok tmp pii
 306   SUBROUTINE SATADJ_D(qv, qvd, qc, qcd, qr, qrd, tmp, tmpd, p1d, p1dd, &
 307 &    pii, piid, rhok, rhokd, kvts, kvte, xlv, dt, cp, ep2, svp1, svp2, &
 308 &    svp3, svpt0)
 309     IMPLICIT NONE
 310     INTEGER :: kvts, kvte, k
 311     REAL, DIMENSION(kvts:kvte) :: qv, qc, qr, tmp, p1d, pii, rhok
 312     REAL, DIMENSION(kvts:kvte) :: qvd, qcd, qrd, tmpd, p1dd, piid, rhokd
 313     REAL, DIMENSION(kvts:kvte) :: rcgs, pressure, temp, es, qvs
 314     REAL, DIMENSION(kvts:kvte) :: rcgsd, pressured, tempd, esd, qvsd
 315     REAL, DIMENSION(kvts:kvte) :: ern, qv2cl, rn2qv
 316     REAL, DIMENSION(kvts:kvte) :: ernd, qv2cld, rn2qvd
 317 ! local var
 318     REAL :: svp1, svp2, svp3, svpt0, ep2, xlv, cp, dt, f5
 319     REAL :: ernmax, product
 320     REAL :: ernmaxd, productd
 321     REAL :: arg1
 322     REAL :: arg1d
 323     INTRINSIC EXP
 324     f5 = svp2*(svpt0-svp3)*xlv/cp
 325     tempd = 0.0
 326     rcgsd = 0.0
 327     pressured = 0.0
 328     esd = 0.0
 329     qvsd = 0.0
 330     DO k=kvts,kvte
 331 !constant
 332       rcgsd(k) = 0.001*rhokd(k)
 333       rcgs(k) = 0.001*rhok(k)
 334       pressured(k) = p1dd(k)
 335       pressure(k) = p1d(k)
 336       tempd(k) = piid(k)*tmp(k) + pii(k)*tmpd(k)
 337       temp(k) = pii(k)*tmp(k)
 338       arg1d = (svp2*tempd(k)*(temp(k)-svp3)-svp2*(temp(k)-svpt0)*tempd(k&
 339 &        ))/(temp(k)-svp3)**2
 340       arg1 = svp2*(temp(k)-svpt0)/(temp(k)-svp3)
 341       esd(k) = 1000.*svp1*arg1d*EXP(arg1)
 342       es(k) = 1000.*svp1*EXP(arg1)
 343       qvsd(k) = (ep2*esd(k)*(pressure(k)-es(k))-ep2*es(k)*(pressured(k)-&
 344 &        esd(k)))/(pressure(k)-es(k))**2
 345       qvs(k) = ep2*es(k)/(pressure(k)-es(k))
 346       IF (qr(k) .LT. 0) THEN
 347         qrd(k) = 0.0
 348         qr(k) = 0.0
 349       END IF
 350       IF (qv(k) .LT. 0) THEN
 351         qvd(k) = 0.0
 352         qv(k) = 0.0
 353       END IF
 354       IF (qc(k) .LT. 0) THEN
 355         qcd(k) = 0.0
 356         qc(k) = 0.0
 357       END IF
 358     END DO
 359     ernd = 0.0
 360     qv2cld = 0.0
 361     rn2qvd = 0.0
 362     DO k=kvts,kvte
 363 !not related to time; maximum transform qv to cl (sat) or cl to qv (sub sat)
 364       qv2cld(k) = ((qvd(k)-qvsd(k))*(1.+pressure(k)/(pressure(k)-es(k))*&
 365 &        qvs(k)*f5/(temp(k)-svp3)**2)-(qv(k)-qvs(k))*(f5*((pressured(k)*(&
 366 &        pressure(k)-es(k))-pressure(k)*(pressured(k)-esd(k)))*qvs(k)/(&
 367 &        pressure(k)-es(k))**2+pressure(k)*qvsd(k)/(pressure(k)-es(k)))*(&
 368 &        temp(k)-svp3)**2-pressure(k)*qvs(k)*f5*2*(temp(k)-svp3)*tempd(k)&
 369 &        /(pressure(k)-es(k)))/(temp(k)-svp3)**4)/(1.+pressure(k)/(&
 370 &        pressure(k)-es(k))*qvs(k)*f5/(temp(k)-svp3)**2)**2
 371       qv2cl(k) = (qv(k)-qvs(k))/(1.+pressure(k)/(pressure(k)-es(k))*qvs(&
 372 &        k)*f5/(temp(k)-svp3)**2)
 373 ! sub sat rain evaperate
 374       rn2qvd(k) = 0.0
 375       rn2qv(k) = 0.0
 376       ernd(k) = 0.0
 377       ern(k) = 0.0
 378       IF (qvs(k) .GT. qv(k)) THEN
 379         IF (qr(k) .GE. 1d-5) THEN
 380           rn2qvd(k) = dt*(((124.9*.2046*(rcgs(k)*qr(k))**(-0.7954)*(&
 381 &            rcgsd(k)*qr(k)+rcgs(k)*qrd(k))*(rcgs(k)*qr(k))**.525+(1.6+&
 382 &            124.9*(rcgs(k)*qr(k))**.2046)*.525*(rcgs(k)*qr(k))**(-0.475)&
 383 &            *(rcgsd(k)*qr(k)+rcgs(k)*qrd(k)))*(2.55e8/(pressure(k)*qvs(k&
 384 &            ))+5.4e5)+(1.6+124.9*(rcgs(k)*qr(k))**.2046)*(rcgs(k)*qr(k))&
 385 &            **.525*2.55e8*(pressured(k)*qvs(k)+pressure(k)*qvsd(k))/(&
 386 &            pressure(k)**2*qvs(k)**2))*(qvs(k)-qv(k))/((2.55e8/(pressure&
 387 &            (k)*qvs(k))+5.4e5)**2*rcgs(k)*qvs(k))+(1.6+124.9*(rcgs(k)*qr&
 388 &            (k))**.2046)*(rcgs(k)*qr(k))**.525*((qvsd(k)-qvd(k))*rcgs(k)&
 389 &            *qvs(k)-(qvs(k)-qv(k))*(rcgsd(k)*qvs(k)+rcgs(k)*qvsd(k)))/((&
 390 &            2.55e8/(pressure(k)*qvs(k))+5.4e5)*rcgs(k)**2*qvs(k)**2))
 391           rn2qv(k) = dt*((1.6+124.9*(rcgs(k)*qr(k))**.2046)*(rcgs(k)*qr(&
 392 &            k))**.525/(2.55e8/(pressure(k)*qvs(k))+5.4e5))*((qvs(k)-qv(k&
 393 &            ))/(rcgs(k)*qvs(k)))
 394         ELSE
 395           rn2qvd(k) = 0.0
 396           rn2qv(k) = 0.0
 397         END IF
 398         IF (rn2qv(k) .GT. qr(k)) THEN
 399           rn2qvd(k) = qrd(k)
 400           rn2qv(k) = qr(k)
 401         END IF
 402         ernmax = 0.0
 403         IF (-qv2cl(k) - qc(k) .GT. 0.0) THEN
 404           ernmaxd = -qv2cld(k) - qcd(k)
 405           ernmax = -qv2cl(k) - qc(k)
 406         ELSE
 407           ernmaxd = 0.0
 408         END IF
 409 !        ern(k)  = amin1(rn2qv(k), ernmax)
 410         ernd(k) = rn2qvd(k)
 411         ern(k) = rn2qv(k)
 412         IF (rn2qv(k) .GT. ernmax) THEN
 413           ernd(k) = ernmaxd
 414           ern(k) = ernmax
 415         END IF
 416       END IF
 417 ! Update all variables
 418 !       product = amax1(qv2cl(k),-qc(k))
 419       productd = qv2cld(k)
 420       product = qv2cl(k)
 421       IF (qv2cl(k) .LT. -qc(k)) THEN
 422         productd = -qcd(k)
 423         product = -qc(k)
 424       END IF
 425 !       qv(k) = amax1(qv(k) - product + ern(k),0.)
 426       qvd(k) = qvd(k) - productd + ernd(k)
 427       qv(k) = qv(k) - product + ern(k)
 428       IF (qv(k) .LT. 0) THEN
 429         qvd(k) = 0.0
 430         qv(k) = 0.0
 431       END IF
 432       qcd(k) = qcd(k) + productd
 433       qc(k) = qc(k) + product
 434       qrd(k) = qrd(k) - ernd(k)
 435       qr(k) = qr(k) - ern(k)
 436       tempd(k) = tempd(k) + xlv*(productd-ernd(k))/cp
 437       temp(k) = temp(k) + xlv/cp*(product-ern(k))
 438       tmpd(k) = (tempd(k)*pii(k)-temp(k)*piid(k))/pii(k)**2
 439       tmp(k) = temp(k)/pii(k)
 440     END DO
 441   END SUBROUTINE SATADJ_D
 442
 443 END MODULE G_MODULE_MP_MKESSLER