wrftladj/module_mp_mkessler_ad.F

   1 !        Generated by TAPENADE     (INRIA, Tropics team)
   2 !  Tapenade 3.3 (r3163) - 09/25/2009 09:04
   3 !
   4 MODULE A_MODULE_MP_MKESSLER
   5   IMPLICIT NONE
   6
   7 CONTAINS
   8 !  Differentiation of kessler in reverse (adjoint) mode:
   9 !   gradient, with respect to input variables: qc p t qr qv rainnc
  10 !                rainncv rho pii
  11 !   of linear combination of output variables: qc p t qr qv rainnc
  12 !                rainncv rho pii
  13 !
  14   SUBROUTINE A_MKESSLER(t, tb, qv, qvb, qc, qcb, qr, qrb, rho, rhob, p, &
  15 &    pb, pii, piib, dt_in, z, xlv, cp, ep2, svp1, svp2, svp3, svpt0, &
  16 &    rhowater, dz8w, rainnc, rainncb, rainncv, rainncvb, ids, ide, jds, &
  17 &    jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts&
  18 &    , kte)
  19
  20     USE module_mp_mkessler, ONLY : SMALLSTEP, RFALL, AUTOCA, SATADJ
  21
  22     IMPLICIT NONE
  23 !----------------------------------------------------------------
  24 ! Restructered from WRF Kessler Warm rain process
  25 ! H.L. Wang Aug. 1 2009
  26 !----------------------------------------------------------------
  27     REAL, PARAMETER :: c1=.001
  28     REAL, PARAMETER :: c2=.001
  29     REAL, PARAMETER :: c3=2.2
  30     REAL, PARAMETER :: c4=.875
  31 !----------------------------------------------------------------
  32     INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
  33 &    jme, kms, kme, its, ite, jts, jte, kts, kte
  34     REAL, INTENT(IN) :: xlv, cp
  35     REAL, INTENT(IN) :: ep2, svp1, svp2, svp3, svpt0
  36     REAL, INTENT(IN) :: rhowater
  37     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t, qv, &
  38 &    qc, qr
  39     REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: tb, qvb, qcb, qrb
  40     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rho, p, &
  41 &    pii, dz8w
  42     REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: rhob, pb, piib
  43     REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: z
  44     REAL, INTENT(IN) :: dt_in
  45     REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: rainnc, rainncv
  46     REAL, DIMENSION(ims:ime, jms:jme) :: rainncb, rainncvb
  47 ! local variables
  48     REAL :: qrprod, ern, gam, rcgs, rcgsi
  49     REAL, DIMENSION(its:ite, kts:kte, jts:jte) :: prod
  50     REAL, DIMENSION(kts:kte) :: vt, prodk, vtden, rdzk, rhok, piik, &
  51 &    factor, rdzw
  52     REAL, DIMENSION(kts:kte) :: rhokb, piikb
  53     INTEGER :: i, j, k
  54     INTEGER :: nfall, n, nfall_new
  55     REAL :: qrr, pressure, temp, es, qvs, dz, dt
  56     REAL :: f5, dtfall, rdz, product
  57     REAL :: vtmax, crmax, factorn
  58     REAL :: qcr, factorr, ppt
  59     REAL, PARAMETER :: max_cr_sedimentation=0.75
  60 !----------------------------------------------------------------
  61     INTEGER :: imax, kmax
  62 ! whl
  63     REAL, DIMENSION(kts:kte) :: qv1d, qc1d, qr1d, t1d, p1d
  64     REAL, DIMENSION(kts:kte) :: qv1db, qc1db, qr1db, t1db, p1db
  65     REAL :: dtleft, rainncv0, max_cr
  66     REAL :: rainncv0b
  67     INTEGER :: kvts, kvte, kn
  68     INTEGER :: ad_to
  69     dt = dt_in
  70 !  print*,'begin'
  71 !   print*,its,ite,jts,jte
  72 !   print*,ims,ime,jms,jme
  73 !   print*,ids,ide,jds,jde
  74     rdzk = 0.0
  75     rdzw = 0.0
  76     DO j=jts,jte
  77       DO i=its,ite
  78         DO k=1,kte-1
  79           rdzk(k) = 1./(z(i, k+1, j)-z(i, k, j))
  80         END DO
  81         rdzk(kte) = 1./(z(i, kte, j)-z(i, kte-1, j))
  82       END DO
  83     END DO
  84     DO j=jts,jte
  85       DO i=its,ite
  86         DO k=1,kte
  87           qv1d(k) = qv(i, k, j)
  88           qc1d(k) = qc(i, k, j)
  89           qr1d(k) = qr(i, k, j)
  90           t1d(k) = t(i, k, j)
  91           CALL PUSHREAL8(p1d(k))
  92           p1d(k) = p(i, k, j)
  93           CALL PUSHREAL8(rhok(k))
  94           rhok(k) = rho(i, k, j)
  95           CALL PUSHREAL8(piik(k))
  96           piik(k) = pii(i, k, j)
  97           rdzw(k) = 1./dz8w(i, k, j)
  98         END DO
  99         CALL PUSHINTEGER4(kvts)
 100 !   print*,i,j
 101         kvts = kts
 102         CALL PUSHINTEGER4(kvte)
 103         kvte = kte
 104         CALL PUSHREAL8(dtleft)
 105         dtleft = dt
 106         CALL SMALLSTEP(qr1d, rdzk, rdzw, rhok, max_cr, dtleft, nfall, &
 107 &                 kvts, kvte)
 108         dtleft = dt/nfall
 109         DO kn=1,nfall
 110           CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
 111           CALL RFALL(qr1d, rdzk, rdzw, rhok, rainncv0, rhowater, max_cr&
 112 &               , dtleft, kvts, kvte)
 113         END DO
 114         CALL PUSHINTEGER4(kn - 1)
 115         CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
 116         CALL PUSHREAL8ARRAY(qc1d, kte - kts + 1)
 117 !    print*,rainncv0
 118 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
 119 !autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
 120         CALL AUTOCA(qc1d, qr1d, kvts, kvte, c1, c2, c3, c4, dt)
 121         CALL PUSHREAL8ARRAY(t1d, kte - kts + 1)
 122         CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
 123         CALL PUSHREAL8ARRAY(qc1d, kte - kts + 1)
 124         CALL PUSHREAL8ARRAY(qv1d, kte - kts + 1)
 125 !satadj(qv,qc,qr, tmp, pii,rho,  kvts,kvte,xlv, cp,EP2,SVP1,SVP2,SVP3,SVPT0)
 126         CALL SATADJ(qv1d, qc1d, qr1d, t1d, p1d, piik, rhok, kvts, kvte, &
 127 &              xlv, dt, cp, ep2, svp1, svp2, svp3, svpt0)
 128 !      END DO
 129 !    END DO
 130     qv1db = 0.0
 131     qc1db = 0.0
 132     p1db = 0.0
 133     qr1db = 0.0
 134     piikb = 0.0
 135     rhokb = 0.0
 136     t1db = 0.0
 137 !    DO j=jte,jts,-1
 138 !      DO i=ite,its,-1
 139         DO k=kte,1,-1
 140           t1db(k) = t1db(k) + tb(i, k, j)
 141           tb(i, k, j) = 0.0
 142           qr1db(k) = qr1db(k) + qrb(i, k, j)
 143           qrb(i, k, j) = 0.0
 144           qc1db(k) = qc1db(k) + qcb(i, k, j)
 145           qcb(i, k, j) = 0.0
 146           qv1db(k) = qv1db(k) + qvb(i, k, j)
 147           qvb(i, k, j) = 0.0
 148         END DO
 149         CALL POPREAL8ARRAY(qv1d, kte - kts + 1)
 150         CALL POPREAL8ARRAY(qc1d, kte - kts + 1)
 151         CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
 152         CALL POPREAL8ARRAY(t1d, kte - kts + 1)
 153         CALL SATADJ_B(qv1d, qv1db, qc1d, qc1db, qr1d, qr1db, t1d, t1db, &
 154 &                p1d, p1db, piik, piikb, rhok, rhokb, kvts, kvte, xlv, dt&
 155 &                , cp, ep2, svp1, svp2, svp3, svpt0)
 156         CALL POPREAL8ARRAY(qc1d, kte - kts + 1)
 157         CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
 158         CALL AUTOCA_B(qc1d, qc1db, qr1d, qr1db, kvts, kvte, c1, c2, c3, &
 159 &                c4, dt)
 160         rainncvb(i, j) = rainncvb(i, j) + rainncb(i, j)
 161         CALL POPINTEGER4(ad_to)
 162         DO kn=ad_to,1,-1
 163           rainncv0b = rainncvb(i, j)
 164           CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
 165           CALL RFALL_B(qr1d, qr1db, rdzk, rdzw, rhok, rhokb, rainncv0, &
 166 &                 rainncv0b, rhowater, max_cr, dtleft, kvts, kvte)
 167         END DO
 168         rainncvb(i, j) = 0.0
 169         CALL POPREAL8(dtleft)
 170         CALL POPINTEGER4(kvte)
 171         CALL POPINTEGER4(kvts)
 172         DO k=kte,1,-1
 173           CALL POPREAL8(piik(k))
 174           piib(i, k, j) = piib(i, k, j) + piikb(k)
 175           piikb(k) = 0.0
 176           CALL POPREAL8(rhok(k))
 177           rhob(i, k, j) = rhob(i, k, j) + rhokb(k)
 178           rhokb(k) = 0.0
 179           CALL POPREAL8(p1d(k))
 180           pb(i, k, j) = pb(i, k, j) + p1db(k)
 181           p1db(k) = 0.0
 182           tb(i, k, j) = tb(i, k, j) + t1db(k)
 183           t1db(k) = 0.0
 184           qrb(i, k, j) = qrb(i, k, j) + qr1db(k)
 185           qr1db(k) = 0.0
 186           qcb(i, k, j) = qcb(i, k, j) + qc1db(k)
 187           qc1db(k) = 0.0
 188           qvb(i, k, j) = qvb(i, k, j) + qv1db(k)
 189           qv1db(k) = 0.0
 190         END DO
 191       END DO
 192     END DO
 193   END SUBROUTINE A_MKESSLER
 194 !
 195 !  Differentiation of rfall in reverse (adjoint) mode:
 196 !   gradient, with respect to input variables: prodk rhok
 197 !   of linear combination of output variables: prodk rainncv0 rhok
 198   SUBROUTINE RFALL_B(prodk, prodkb, rdzk, rdzw, rhok, rhokb, rainncv0, &
 199 &    rainncv0b, rhowat, max_cr, dtfall, kvts, kvte)
 200     IMPLICIT NONE
 201     INTEGER :: k, kvts, kvte
 202     REAL, DIMENSION(kvts:kvte) :: vtden, vt, prodk, factor, rdzk, rdzw, &
 203 &    rhok
 204     REAL, DIMENSION(kvts:kvte) :: vtdenb, vtb, prodkb, factorb, rhokb
 205     REAL :: rainncv0, rhowat, max_cr, ppt, dtleft
 206     REAL :: rainncv0b, pptb
 207     REAL :: qrr, dtfall
 208     REAL :: qrrb
 209     REAL :: arg1
 210     REAL :: arg1b
 211     INTRINSIC SQRT
 212     INTEGER :: branch
 213     REAL :: temp0
 214     REAL :: tempb1
 215     REAL :: tempb0
 216     REAL :: temp0b
 217     REAL :: tempb
 218     REAL :: temp
 219     DO k=kvts,kvte
 220       IF (prodk(k) .LT. 0) THEN
 221         prodk(k) = 0.0
 222         CALL PUSHINTEGER4(2)
 223       ELSE
 224         CALL PUSHINTEGER4(1)
 225       END IF
 226     END DO
 227     DO k=kvts,kvte
 228       CALL PUSHREAL8(qrr)
 229       qrr = prodk(k)*0.001*rhok(k)
 230       CALL PUSHREAL8(arg1)
 231       arg1 = rhok(1)/rhok(k)
 232       vtden(k) = SQRT(arg1)
 233       IF (qrr/(0.001*rhok(k)) .GE. 1d-5) THEN
 234         vt(k) = 36.34*qrr**0.1364*vtden(k)
 235         CALL PUSHINTEGER4(1)
 236       ELSE
 237         vt(k) = 0.0
 238         CALL PUSHINTEGER4(2)
 239       END IF
 240     END DO
 241 !     pause
 242     DO k=kvts,kvte-1
 243       factor(k) = dtfall*rdzk(k)/rhok(k)
 244     END DO
 245     factor(kvte) = dtfall*rdzk(kvte)
 246     k = 1
 247 !mm
 248     CALL PUSHINTEGER4(k)
 249 !      print*,rainncv0
 250 !------------------------------------------------------------------------------
 251 ! Time split loop, Fallout done with flux upstream
 252 !------------------------------------------------------------------------------
 253     DO k=kvts,kvte-1
 254       CALL PUSHREAL8(prodk(k))
 255       prodk(k) = prodk(k) - factor(k)*(rhok(k)*prodk(k)*vt(k)-rhok(k+1)*&
 256 &        prodk(k+1)*vt(k+1))
 257     END DO
 258     k = kvte
 259     CALL PUSHREAL8(prodk(k))
 260     prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
 261     CALL PUSHINTEGER4(k)
 262     DO k=kvts,kvte
 263       IF (prodk(k) .LT. 0) THEN
 264         CALL PUSHINTEGER4(2)
 265       ELSE
 266         CALL PUSHINTEGER4(1)
 267       END IF
 268     END DO
 269     DO k=kvte,kvts,-1
 270       CALL POPINTEGER4(branch)
 271       IF (.NOT.branch .LT. 2) prodkb(k) = 0.0
 272     END DO
 273     CALL POPINTEGER4(k)
 274     vtb = 0.0
 275     factorb = 0.0
 276     CALL POPREAL8(prodk(k))
 277     factorb(k) = -(vt(k)*prodk(k)*prodkb(k))
 278     vtb(k) = -(factor(k)*prodk(k)*prodkb(k))
 279     prodkb(k) = (1.0-vt(k)*factor(k))*prodkb(k)
 280     DO k=kvte-1,kvts,-1
 281       CALL POPREAL8(prodk(k))
 282       temp = rhok(k+1)*prodk(k+1)
 283       temp0 = rhok(k)*prodk(k)
 284       temp0b = -(factor(k)*prodkb(k))
 285       tempb1 = -(vt(k+1)*temp0b)
 286       factorb(k) = factorb(k) - (temp0*vt(k)-temp*vt(k+1))*prodkb(k)
 287       rhokb(k) = rhokb(k) + vt(k)*prodk(k)*temp0b
 288       vtb(k) = vtb(k) + temp0*temp0b
 289       rhokb(k+1) = rhokb(k+1) + prodk(k+1)*tempb1
 290       prodkb(k+1) = prodkb(k+1) + rhok(k+1)*tempb1
 291       vtb(k+1) = vtb(k+1) - temp*temp0b
 292       prodkb(k) = vt(k)*rhok(k)*temp0b + prodkb(k)
 293     END DO
 294     CALL POPINTEGER4(k)
 295     pptb = 1000.*rainncv0b
 296     tempb0 = dtfall*rhok(k)*pptb/rhowat
 297     rhokb(k) = rhokb(k) + dtfall*prodk(k)*vt(k)*pptb/rhowat
 298     prodkb(k) = prodkb(k) + vt(k)*tempb0
 299     vtb(k) = vtb(k) + prodk(k)*tempb0
 300     factorb(kvte) = 0.0
 301     DO k=kvte-1,kvts,-1
 302       rhokb(k) = rhokb(k) - dtfall*rdzk(k)*factorb(k)/rhok(k)**2
 303       factorb(k) = 0.0
 304     END DO
 305     vtdenb = 0.0
 306     DO k=kvte,kvts,-1
 307       CALL POPINTEGER4(branch)
 308       IF (branch .LT. 2) THEN
 309         qrrb = 36.34*vtden(k)*0.1364*qrr**(-0.8636)*vtb(k)
 310         vtdenb(k) = vtdenb(k) + 36.34*qrr**0.1364*vtb(k)
 311         vtb(k) = 0.0
 312       ELSE
 313         vtb(k) = 0.0
 314         qrrb = 0.0
 315       END IF
 316       IF (arg1 .EQ. 0.0) THEN
 317         arg1b = 0.0
 318       ELSE
 319         arg1b = vtdenb(k)/(2.0*SQRT(arg1))
 320       END IF
 321       vtdenb(k) = 0.0
 322       CALL POPREAL8(arg1)
 323       tempb = arg1b/rhok(k)
 324       rhokb(1) = rhokb(1) + tempb
 325       rhokb(k) = rhokb(k) + 0.001*prodk(k)*qrrb - rhok(1)*tempb/rhok(k)
 326       CALL POPREAL8(qrr)
 327       prodkb(k) = prodkb(k) + 0.001*rhok(k)*qrrb
 328     END DO
 329     DO k=kvte,kvts,-1
 330       CALL POPINTEGER4(branch)
 331       IF (.NOT.branch .LT. 2) prodkb(k) = 0.0
 332     END DO
 333   END SUBROUTINE RFALL_B
 334
 335 !  Differentiation of autoca in reverse (adjoint) mode:
 336 !   gradient, with respect to input variables: qc1d qr1d
 337 !   of linear combination of output variables: qc1d qr1d
 338   SUBROUTINE AUTOCA_B(qc1d, qc1db, qr1d, qr1db, kvts, kvte, c1, c2, c3, &
 339 &    c4, dt)
 340     IMPLICIT NONE
 341 !     print*,k,qrprod
 342     INTEGER :: kvts, kvte, k
 343     REAL, DIMENSION(kvts:kvte) :: qc1d, qr1d
 344     REAL, DIMENSION(kvts:kvte) :: qc1db, qr1db
 345     REAL :: c1, c2, c3, c4
 346     REAL :: qrrc, dt, factorn, qrprod, qrprod2
 347     REAL :: factornb, qrprodb, qrprod2b
 348     REAL :: pwr1
 349     REAL :: pwr1b
 350     INTEGER :: branch
 351     REAL :: temp0b
 352     REAL :: temp
 353     qrrc = 1.0e-5
 354     DO k=kvts,kvte
 355       IF (qr1d(k) .LT. 0.0) THEN
 356         qr1d(k) = 0.0
 357         CALL PUSHINTEGER4(1)
 358       ELSE
 359         CALL PUSHINTEGER4(0)
 360       END IF
 361       IF (qc1d(k) .LT. 0.0) THEN
 362         qc1d(k) = 0.0
 363         CALL PUSHINTEGER4(1)
 364       ELSE
 365         CALL PUSHINTEGER4(0)
 366       END IF
 367       IF (qr1d(k) .GE. qrrc) THEN
 368         CALL PUSHREAL8(pwr1)
 369         pwr1 = qr1d(k)**c4
 370         CALL PUSHREAL8(factorn)
 371         factorn = 1.0/(1.+c3*dt*pwr1)
 372         CALL PUSHINTEGER4(0)
 373       ELSE
 374         CALL PUSHREAL8(factorn)
 375         factorn = 1.0
 376         CALL PUSHINTEGER4(1)
 377       END IF
 378       qrprod = qc1d(k)*(1.0-factorn)
 379       qrprod2 = 0.0
 380       IF (qc1d(k) - c2 .GT. 0) THEN
 381         qrprod2 = factorn*c1*dt*(qc1d(k)-c2)
 382         IF (qrprod2 .GT. qc1d(k) - c2) THEN
 383           qrprod2 = qc1d(k) - c2
 384           CALL PUSHINTEGER4(2)
 385         ELSE
 386           CALL PUSHINTEGER4(1)
 387         END IF
 388       ELSE
 389         CALL PUSHINTEGER4(0)
 390       END IF
 391 !        print*,k,qrprod2
 392       qrprod = qrprod + qrprod2
 393       IF (qc1d(k) - qrprod .GT. 0) THEN
 394         CALL PUSHINTEGER4(1)
 395       ELSE
 396         CALL PUSHINTEGER4(2)
 397       END IF
 398     END DO
 399     DO k=kvte,kvts,-1
 400       CALL POPINTEGER4(branch)
 401       IF (branch .LT. 2) THEN
 402         qrprodb = qr1db(k) - qc1db(k)
 403       ELSE
 404         qrprodb = qr1db(k)
 405         qc1db(k) = 0.0
 406         qrprodb = 0.0
 407       END IF
 408       qrprod2b = qrprodb
 409       CALL POPINTEGER4(branch)
 410       IF (branch .LT. 2) THEN
 411         IF (branch .LT. 1) THEN
 412           factornb = 0.0
 413           GOTO 100
 414         END IF
 415       ELSE
 416         qc1db(k) = qc1db(k) + qrprod2b
 417         qrprod2b = 0.0
 418       END IF
 419       temp0b = c1*dt*qrprod2b
 420       factornb = (qc1d(k)-c2)*temp0b
 421       qc1db(k) = qc1db(k) + factorn*temp0b
 422  100  qc1db(k) = qc1db(k) + (1.0-factorn)*qrprodb
 423       factornb = factornb - qc1d(k)*qrprodb
 424       CALL POPINTEGER4(branch)
 425       IF (branch .LT. 1) THEN
 426         CALL POPREAL8(factorn)
 427         temp = c3*dt*pwr1 + 1.
 428         pwr1b = -(c3*dt*factornb/temp**2)
 429         CALL POPREAL8(pwr1)
 430         IF (.NOT.(qr1d(k) .LE. 0.0 .AND. (c4 .EQ. 0.0 .OR. c4 .NE. INT(&
 431 &            c4)))) qr1db(k) = qr1db(k) + c4*qr1d(k)**(c4-1)*pwr1b
 432       ELSE
 433         CALL POPREAL8(factorn)
 434       END IF
 435       CALL POPINTEGER4(branch)
 436       IF (.NOT.branch .LT. 1) qc1db(k) = 0.0
 437       CALL POPINTEGER4(branch)
 438       IF (.NOT.branch .LT. 1) qr1db(k) = 0.0
 439     END DO
 440   END SUBROUTINE AUTOCA_B
 441
 442 !  Differentiation of satadj in reverse (adjoint) mode:
 443 !   gradient, with respect to input variables: qc qr qv p1d rhok
 444 !                tmp pii
 445 !   of linear combination of output variables: qc qr qv p1d rhok
 446 !                tmp pii
 447   SUBROUTINE SATADJ_B(qv, qvb, qc, qcb, qr, qrb, tmp, tmpb, p1d, p1db, &
 448 &    pii, piib, rhok, rhokb, kvts, kvte, xlv, dt, cp, ep2, svp1, svp2, &
 449 &    svp3, svpt0)
 450     IMPLICIT NONE
 451     INTEGER :: kvts, kvte, k
 452     REAL, DIMENSION(kvts:kvte) :: qv, qc, qr, tmp, p1d, pii, rhok
 453     REAL, DIMENSION(kvts:kvte) :: qvb, qcb, qrb, tmpb, p1db, piib, rhokb
 454     REAL, DIMENSION(kvts:kvte) :: rcgs, pressure, temp, es, qvs
 455     REAL, DIMENSION(kvts:kvte) :: rcgsb, pressureb, tempb, esb, qvsb
 456     REAL, DIMENSION(kvts:kvte) :: ern, qv2cl, rn2qv
 457     REAL, DIMENSION(kvts:kvte) :: ernb, qv2clb, rn2qvb
 458 ! local var
 459     REAL :: svp1, svp2, svp3, svpt0, ep2, xlv, cp, dt, f5
 460     REAL :: ernmax, product
 461     REAL :: ernmaxb, productb
 462     REAL :: arg1
 463     REAL :: arg1b
 464     INTRINSIC EXP
 465     INTEGER :: branch
 466     REAL :: temp3
 467     REAL :: temp2
 468     REAL :: temp1
 469     REAL :: temp0
 470     REAL :: temp13b
 471     REAL :: temp7b
 472     REAL :: temp13b0
 473     REAL :: temp0b
 474     REAL :: temp6b
 475     REAL :: temp12
 476     REAL :: temp11
 477     REAL :: temp10
 478     REAL :: temp9b
 479     REAL :: temp0b3
 480     REAL :: temp0b2
 481     REAL :: temp0b1
 482     REAL :: temp0b0
 483     REAL :: temp2b
 484     REAL :: temp5b
 485     REAL :: temp8b
 486     REAL :: temp1b
 487     REAL :: temp9
 488     REAL :: temp8
 489     REAL :: temp7
 490     REAL :: temp6
 491     REAL :: temp4b
 492     REAL :: temp5
 493     REAL :: temp4
 494     f5 = svp2*(svpt0-svp3)*xlv/cp
 495     DO k=kvts,kvte
 496 !constant
 497       rcgs(k) = 0.001*rhok(k)
 498       pressure(k) = p1d(k)
 499       temp(k) = pii(k)*tmp(k)
 500       CALL PUSHREAL8(arg1)
 501       arg1 = svp2*(temp(k)-svpt0)/(temp(k)-svp3)
 502       es(k) = 1000.*svp1*EXP(arg1)
 503       qvs(k) = ep2*es(k)/(pressure(k)-es(k))
 504       IF (qr(k) .LT. 0) THEN
 505         qr(k) = 0.0
 506         CALL PUSHINTEGER4(1)
 507       ELSE
 508         CALL PUSHINTEGER4(0)
 509       END IF
 510       IF (qv(k) .LT. 0) THEN
 511         qv(k) = 0.0
 512         CALL PUSHINTEGER4(1)
 513       ELSE
 514         CALL PUSHINTEGER4(0)
 515       END IF
 516       IF (qc(k) .LT. 0) THEN
 517         qc(k) = 0.0
 518         CALL PUSHINTEGER4(2)
 519       ELSE
 520         CALL PUSHINTEGER4(1)
 521       END IF
 522     END DO
 523     DO k=kvts,kvte
 524 !not related to time; maximum transform qv to cl (sat) or cl to qv (sub sat)
 525       qv2cl(k) = (qv(k)-qvs(k))/(1.+pressure(k)/(pressure(k)-es(k))*qvs(&
 526 &        k)*f5/(temp(k)-svp3)**2)
 527 ! sub sat rain evaperate
 528       ern(k) = 0.0
 529       IF (qvs(k) .GT. qv(k)) THEN
 530         IF (qr(k) .GE. 1d-5) THEN
 531           rn2qv(k) = dt*((1.6+124.9*(rcgs(k)*qr(k))**.2046)*(rcgs(k)*qr(&
 532 &            k))**.525/(2.55e8/(pressure(k)*qvs(k))+5.4e5))*((qvs(k)-qv(k&
 533 &            ))/(rcgs(k)*qvs(k)))
 534           CALL PUSHINTEGER4(0)
 535         ELSE
 536           rn2qv(k) = 0.0
 537           CALL PUSHINTEGER4(1)
 538         END IF
 539         IF (rn2qv(k) .GT. qr(k)) THEN
 540           rn2qv(k) = qr(k)
 541           CALL PUSHINTEGER4(1)
 542         ELSE
 543           CALL PUSHINTEGER4(0)
 544         END IF
 545         ernmax = 0.0
 546         IF (-qv2cl(k) - qc(k) .GT. 0.0) THEN
 547           ernmax = -qv2cl(k) - qc(k)
 548           CALL PUSHINTEGER4(1)
 549         ELSE
 550           CALL PUSHINTEGER4(0)
 551         END IF
 552 !        ern(k)  = amin1(rn2qv(k), ernmax)
 553         ern(k) = rn2qv(k)
 554         IF (rn2qv(k) .GT. ernmax) THEN
 555           ern(k) = ernmax
 556           CALL PUSHINTEGER4(2)
 557         ELSE
 558           CALL PUSHINTEGER4(1)
 559         END IF
 560       ELSE
 561         CALL PUSHINTEGER4(0)
 562       END IF
 563 ! Update all variables
 564 !       product = amax1(qv2cl(k),-qc(k))
 565       product = qv2cl(k)
 566       IF (qv2cl(k) .LT. -qc(k)) THEN
 567         product = -qc(k)
 568         CALL PUSHINTEGER4(1)
 569       ELSE
 570         CALL PUSHINTEGER4(0)
 571       END IF
 572       CALL PUSHREAL8(qv(k))
 573 !       qv(k) = amax1(qv(k) - product + ern(k),0.)
 574       qv(k) = qv(k) - product + ern(k)
 575       IF (qv(k) .LT. 0) THEN
 576         CALL PUSHINTEGER4(1)
 577       ELSE
 578         CALL PUSHINTEGER4(0)
 579       END IF
 580       CALL PUSHREAL8(temp(k))
 581       temp(k) = temp(k) + xlv/cp*(product-ern(k))
 582     END DO
 583     ernb = 0.0
 584     tempb = 0.0
 585     rcgsb = 0.0
 586     pressureb = 0.0
 587     qv2clb = 0.0
 588     esb = 0.0
 589     qvsb = 0.0
 590     rn2qvb = 0.0
 591     DO k=kvte,kvts,-1
 592       temp13b = tmpb(k)/pii(k)
 593       tempb(k) = tempb(k) + temp13b
 594       piib(k) = piib(k) - temp(k)*temp13b/pii(k)
 595       tmpb(k) = 0.0
 596       CALL POPREAL8(temp(k))
 597       temp13b0 = xlv*tempb(k)/cp
 598       productb = qcb(k) + temp13b0
 599       ernb(k) = ernb(k) - qrb(k) - temp13b0
 600       CALL POPINTEGER4(branch)
 601       IF (.NOT.branch .LT. 1) qvb(k) = 0.0
 602       CALL POPREAL8(qv(k))
 603       productb = productb - qvb(k)
 604       ernb(k) = ernb(k) + qvb(k)
 605       CALL POPINTEGER4(branch)
 606       IF (.NOT.branch .LT. 1) THEN
 607         qcb(k) = qcb(k) - productb
 608         productb = 0.0
 609       END IF
 610       qv2clb(k) = qv2clb(k) + productb
 611       CALL POPINTEGER4(branch)
 612       IF (branch .LT. 2) THEN
 613         IF (branch .LT. 1) THEN
 614           GOTO 100
 615         ELSE
 616           ernmaxb = 0.0
 617         END IF
 618       ELSE
 619         ernmaxb = ernb(k)
 620         ernb(k) = 0.0
 621       END IF
 622       rn2qvb(k) = rn2qvb(k) + ernb(k)
 623       ernb(k) = 0.0
 624       CALL POPINTEGER4(branch)
 625       IF (.NOT.branch .LT. 1) THEN
 626         qv2clb(k) = qv2clb(k) - ernmaxb
 627         qcb(k) = qcb(k) - ernmaxb
 628       END IF
 629       CALL POPINTEGER4(branch)
 630       IF (.NOT.branch .LT. 1) THEN
 631         qrb(k) = qrb(k) + rn2qvb(k)
 632         rn2qvb(k) = 0.0
 633       END IF
 634       CALL POPINTEGER4(branch)
 635       IF (branch .LT. 1) THEN
 636         temp12 = rcgs(k)*qvs(k)
 637         temp4 = pressure(k)*qvs(k)
 638         temp11 = 2.55e8/temp4
 639         temp5 = (temp11+5.4e5)*temp12
 640         temp6 = rcgs(k)*qr(k)
 641         temp10 = 124.9*temp6**.2046 + 1.6
 642         temp7 = temp10/temp5
 643         temp8 = rcgs(k)*qr(k)
 644         temp9 = temp8**.525
 645         temp9b = dt*temp7*rn2qvb(k)
 646         temp8b = (qvs(k)-qv(k))*.525*temp8**(-0.475)*temp9b
 647         temp7b = dt*temp9*(qvs(k)-qv(k))*rn2qvb(k)/temp5
 648         temp6b = 124.9*.2046*temp6**(-0.7954)*temp7b
 649         temp5b = -(temp7*temp7b)
 650         temp4b = -(temp12*temp11*temp5b/temp4)
 651         rcgsb(k) = rcgsb(k) + (temp11+5.4e5)*qvs(k)*temp5b + qr(k)*&
 652 &          temp6b + qr(k)*temp8b
 653         qrb(k) = qrb(k) + rcgs(k)*temp6b + rcgs(k)*temp8b
 654         qvsb(k) = qvsb(k) + (temp11+5.4e5)*rcgs(k)*temp5b + pressure(k)*&
 655 &          temp4b + temp9*temp9b
 656         qvb(k) = qvb(k) - temp9*temp9b
 657         pressureb(k) = pressureb(k) + qvs(k)*temp4b
 658         rn2qvb(k) = 0.0
 659       ELSE
 660         rn2qvb(k) = 0.0
 661       END IF
 662  100  ernb(k) = 0.0
 663       rn2qvb(k) = 0.0
 664       temp3 = (temp(k)-svp3)**2
 665       temp0 = (pressure(k)-es(k))*temp3
 666       temp2 = pressure(k)*qvs(k)
 667       temp1 = temp2/temp0
 668       temp2b = qv2clb(k)/(f5*temp1+1.)
 669       temp1b = -((qv(k)-qvs(k))*f5*temp2b/((f5*temp1+1.)*temp0))
 670       temp0b2 = -(temp1*temp1b)
 671       temp0b3 = temp3*temp0b2
 672       qvb(k) = qvb(k) + temp2b
 673       qvsb(k) = qvsb(k) + pressure(k)*temp1b - temp2b
 674       pressureb(k) = pressureb(k) + temp0b3 + qvs(k)*temp1b
 675       esb(k) = esb(k) - temp0b3
 676       tempb(k) = tempb(k) + (pressure(k)-es(k))*2*(temp(k)-svp3)*temp0b2
 677       qv2clb(k) = 0.0
 678     END DO
 679     DO k=kvte,kvts,-1
 680       CALL POPINTEGER4(branch)
 681       IF (.NOT.branch .LT. 2) qcb(k) = 0.0
 682       CALL POPINTEGER4(branch)
 683       IF (.NOT.branch .LT. 1) qvb(k) = 0.0
 684       CALL POPINTEGER4(branch)
 685       IF (.NOT.branch .LT. 1) qrb(k) = 0.0
 686       temp0b = ep2*qvsb(k)/(pressure(k)-es(k))
 687       temp0b0 = -(es(k)*temp0b/(pressure(k)-es(k)))
 688       esb(k) = esb(k) + temp0b - temp0b0
 689       pressureb(k) = pressureb(k) + temp0b0
 690       qvsb(k) = 0.0
 691       arg1b = svp1*1000.*EXP(arg1)*esb(k)
 692       esb(k) = 0.0
 693       CALL POPREAL8(arg1)
 694       temp0b1 = svp2*arg1b/(temp(k)-svp3)
 695       tempb(k) = tempb(k) + (1.0-(temp(k)-svpt0)/(temp(k)-svp3))*temp0b1
 696       piib(k) = piib(k) + tmp(k)*tempb(k)
 697       tmpb(k) = tmpb(k) + pii(k)*tempb(k)
 698       tempb(k) = 0.0
 699       p1db(k) = p1db(k) + pressureb(k)
 700       pressureb(k) = 0.0
 701       rhokb(k) = rhokb(k) + 0.001*rcgsb(k)
 702       rcgsb(k) = 0.0
 703     END DO
 704   END SUBROUTINE SATADJ_B
 705
 706 END MODULE A_MODULE_MP_MKESSLER