!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.3 (r3163) - 09/25/2009 09:04
!
MODULE A_MODULE_MP_MKESSLER
  IMPLICIT NONE

CONTAINS
!  Differentiation of kessler in reverse (adjoint) mode:
!   gradient, with respect to input variables: qc p t qr qv rainnc
!                rainncv rho pii
!   of linear combination of output variables: qc p t qr qv rainnc
!                rainncv rho pii
!
  SUBROUTINE A_MKESSLER(t, tb, qv, qvb, qc, qcb, qr, qrb, rho, rhob, p, &
&    pb, pii, piib, dt_in, z, xlv, cp, ep2, svp1, svp2, svp3, svpt0, &
&    rhowater, dz8w, rainnc, rainncb, rainncv, rainncvb, ids, ide, jds, &
&    jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts&
&    , kte)
    
    USE module_mp_mkessler, ONLY : SMALLSTEP, RFALL, AUTOCA, SATADJ

    IMPLICIT NONE
!----------------------------------------------------------------
! Restructered from WRF Kessler Warm rain process
! H.L. Wang Aug. 1 2009
!----------------------------------------------------------------
    REAL, PARAMETER :: c1=.001
    REAL, PARAMETER :: c2=.001
    REAL, PARAMETER :: c3=2.2
    REAL, PARAMETER :: c4=.875
!----------------------------------------------------------------
    INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
&    jme, kms, kme, its, ite, jts, jte, kts, kte
    REAL, INTENT(IN) :: xlv, cp
    REAL, INTENT(IN) :: ep2, svp1, svp2, svp3, svpt0
    REAL, INTENT(IN) :: rhowater
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t, qv, &
&    qc, qr
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: tb, qvb, qcb, qrb
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rho, p, &
&    pii, dz8w
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: rhob, pb, piib
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: z
    REAL, INTENT(IN) :: dt_in
    REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: rainnc, rainncv
    REAL, DIMENSION(ims:ime, jms:jme) :: rainncb, rainncvb
! local variables
    REAL :: qrprod, ern, gam, rcgs, rcgsi
    REAL, DIMENSION(its:ite, kts:kte, jts:jte) :: prod
    REAL, DIMENSION(kts:kte) :: vt, prodk, vtden, rdzk, rhok, piik, &
&    factor, rdzw
    REAL, DIMENSION(kts:kte) :: rhokb, piikb
    INTEGER :: i, j, k
    INTEGER :: nfall, n, nfall_new
    REAL :: qrr, pressure, temp, es, qvs, dz, dt
    REAL :: f5, dtfall, rdz, product
    REAL :: vtmax, crmax, factorn
    REAL :: qcr, factorr, ppt
    REAL, PARAMETER :: max_cr_sedimentation=0.75
!----------------------------------------------------------------
    INTEGER :: imax, kmax
! whl
    REAL, DIMENSION(kts:kte) :: qv1d, qc1d, qr1d, t1d, p1d
    REAL, DIMENSION(kts:kte) :: qv1db, qc1db, qr1db, t1db, p1db
    REAL :: dtleft, rainncv0, max_cr
    REAL :: rainncv0b
    INTEGER :: kvts, kvte, kn
    INTEGER :: ad_to
    dt = dt_in
!  print*,'begin'
!   print*,its,ite,jts,jte
!   print*,ims,ime,jms,jme
!   print*,ids,ide,jds,jde
    rdzk = 0.0
    rdzw = 0.0
    DO j=jts,jte
      DO i=its,ite
        DO k=1,kte-1
          rdzk(k) = 1./(z(i, k+1, j)-z(i, k, j))
        END DO
        rdzk(kte) = 1./(z(i, kte, j)-z(i, kte-1, j))
      END DO
    END DO
    DO j=jts,jte
      DO i=its,ite
        DO k=1,kte
          qv1d(k) = qv(i, k, j)
          qc1d(k) = qc(i, k, j)
          qr1d(k) = qr(i, k, j)
          t1d(k) = t(i, k, j)
          CALL PUSHREAL8(p1d(k))
          p1d(k) = p(i, k, j)
          CALL PUSHREAL8(rhok(k))
          rhok(k) = rho(i, k, j)
          CALL PUSHREAL8(piik(k))
          piik(k) = pii(i, k, j)
          rdzw(k) = 1./dz8w(i, k, j)
        END DO
        CALL PUSHINTEGER4(kvts)
!   print*,i,j
        kvts = kts
        CALL PUSHINTEGER4(kvte)
        kvte = kte
        CALL PUSHREAL8(dtleft)
        dtleft = dt
        CALL SMALLSTEP(qr1d, rdzk, rdzw, rhok, max_cr, dtleft, nfall, &
&                 kvts, kvte)
        dtleft = dt/nfall
        DO kn=1,nfall
          CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
          CALL RFALL(qr1d, rdzk, rdzw, rhok, rainncv0, rhowater, max_cr&
&               , dtleft, kvts, kvte)
        END DO
        CALL PUSHINTEGER4(kn - 1)
        CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
        CALL PUSHREAL8ARRAY(qc1d, kte - kts + 1)
!    print*,rainncv0
!autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
!autoca(qc1d,qr1d, kvts,kvte,c1,c2,c3,c4,dt )
        CALL AUTOCA(qc1d, qr1d, kvts, kvte, c1, c2, c3, c4, dt)
        CALL PUSHREAL8ARRAY(t1d, kte - kts + 1)
        CALL PUSHREAL8ARRAY(qr1d, kte - kts + 1)
        CALL PUSHREAL8ARRAY(qc1d, kte - kts + 1)
        CALL PUSHREAL8ARRAY(qv1d, kte - kts + 1)
!satadj(qv,qc,qr, tmp, pii,rho,  kvts,kvte,xlv, cp,EP2,SVP1,SVP2,SVP3,SVPT0)
        CALL SATADJ(qv1d, qc1d, qr1d, t1d, p1d, piik, rhok, kvts, kvte, &
&              xlv, dt, cp, ep2, svp1, svp2, svp3, svpt0)
!      END DO
!    END DO
    qv1db = 0.0
    qc1db = 0.0
    p1db = 0.0
    qr1db = 0.0
    piikb = 0.0
    rhokb = 0.0
    t1db = 0.0
!    DO j=jte,jts,-1
!      DO i=ite,its,-1
        DO k=kte,1,-1
          t1db(k) = t1db(k) + tb(i, k, j)
          tb(i, k, j) = 0.0
          qr1db(k) = qr1db(k) + qrb(i, k, j)
          qrb(i, k, j) = 0.0
          qc1db(k) = qc1db(k) + qcb(i, k, j)
          qcb(i, k, j) = 0.0
          qv1db(k) = qv1db(k) + qvb(i, k, j)
          qvb(i, k, j) = 0.0
        END DO
        CALL POPREAL8ARRAY(qv1d, kte - kts + 1)
        CALL POPREAL8ARRAY(qc1d, kte - kts + 1)
        CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
        CALL POPREAL8ARRAY(t1d, kte - kts + 1)
        CALL SATADJ_B(qv1d, qv1db, qc1d, qc1db, qr1d, qr1db, t1d, t1db, &
&                p1d, p1db, piik, piikb, rhok, rhokb, kvts, kvte, xlv, dt&
&                , cp, ep2, svp1, svp2, svp3, svpt0)
        CALL POPREAL8ARRAY(qc1d, kte - kts + 1)
        CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
        CALL AUTOCA_B(qc1d, qc1db, qr1d, qr1db, kvts, kvte, c1, c2, c3, &
&                c4, dt)
        rainncvb(i, j) = rainncvb(i, j) + rainncb(i, j)
        CALL POPINTEGER4(ad_to)
        DO kn=ad_to,1,-1
          rainncv0b = rainncvb(i, j)
          CALL POPREAL8ARRAY(qr1d, kte - kts + 1)
          CALL RFALL_B(qr1d, qr1db, rdzk, rdzw, rhok, rhokb, rainncv0, &
&                 rainncv0b, rhowater, max_cr, dtleft, kvts, kvte)
        END DO
        rainncvb(i, j) = 0.0
        CALL POPREAL8(dtleft)
        CALL POPINTEGER4(kvte)
        CALL POPINTEGER4(kvts)
        DO k=kte,1,-1
          CALL POPREAL8(piik(k))
          piib(i, k, j) = piib(i, k, j) + piikb(k)
          piikb(k) = 0.0
          CALL POPREAL8(rhok(k))
          rhob(i, k, j) = rhob(i, k, j) + rhokb(k)
          rhokb(k) = 0.0
          CALL POPREAL8(p1d(k))
          pb(i, k, j) = pb(i, k, j) + p1db(k)
          p1db(k) = 0.0
          tb(i, k, j) = tb(i, k, j) + t1db(k)
          t1db(k) = 0.0
          qrb(i, k, j) = qrb(i, k, j) + qr1db(k)
          qr1db(k) = 0.0
          qcb(i, k, j) = qcb(i, k, j) + qc1db(k)
          qc1db(k) = 0.0
          qvb(i, k, j) = qvb(i, k, j) + qv1db(k)
          qv1db(k) = 0.0
        END DO
      END DO
    END DO
  END SUBROUTINE A_MKESSLER
!
!  Differentiation of rfall in reverse (adjoint) mode:
!   gradient, with respect to input variables: prodk rhok
!   of linear combination of output variables: prodk rainncv0 rhok
  SUBROUTINE RFALL_B(prodk, prodkb, rdzk, rdzw, rhok, rhokb, rainncv0, &
&    rainncv0b, rhowat, max_cr, dtfall, kvts, kvte)
    IMPLICIT NONE
    INTEGER :: k, kvts, kvte
    REAL, DIMENSION(kvts:kvte) :: vtden, vt, prodk, factor, rdzk, rdzw, &
&    rhok
    REAL, DIMENSION(kvts:kvte) :: vtdenb, vtb, prodkb, factorb, rhokb
    REAL :: rainncv0, rhowat, max_cr, ppt, dtleft
    REAL :: rainncv0b, pptb
    REAL :: qrr, dtfall
    REAL :: qrrb
    REAL :: arg1
    REAL :: arg1b
    INTRINSIC SQRT
    INTEGER :: branch
    REAL :: temp0
    REAL :: tempb1
    REAL :: tempb0
    REAL :: temp0b
    REAL :: tempb
    REAL :: temp
    DO k=kvts,kvte
      IF (prodk(k) .LT. 0) THEN
        prodk(k) = 0.0
        CALL PUSHINTEGER4(2)
      ELSE
        CALL PUSHINTEGER4(1)
      END IF
    END DO
    DO k=kvts,kvte
      CALL PUSHREAL8(qrr)
      qrr = prodk(k)*0.001*rhok(k)
      CALL PUSHREAL8(arg1)
      arg1 = rhok(1)/rhok(k)
      vtden(k) = SQRT(arg1)
      IF (qrr/(0.001*rhok(k)) .GE. 1d-5) THEN
        vt(k) = 36.34*qrr**0.1364*vtden(k)
        CALL PUSHINTEGER4(1)
      ELSE
        vt(k) = 0.0
        CALL PUSHINTEGER4(2)
      END IF
    END DO
!     pause
    DO k=kvts,kvte-1
      factor(k) = dtfall*rdzk(k)/rhok(k)
    END DO
    factor(kvte) = dtfall*rdzk(kvte)
    k = 1
!mm
    CALL PUSHINTEGER4(k)
!      print*,rainncv0
!------------------------------------------------------------------------------
! Time split loop, Fallout done with flux upstream
!------------------------------------------------------------------------------
    DO k=kvts,kvte-1
      CALL PUSHREAL8(prodk(k))
      prodk(k) = prodk(k) - factor(k)*(rhok(k)*prodk(k)*vt(k)-rhok(k+1)*&
&        prodk(k+1)*vt(k+1))
    END DO
    k = kvte
    CALL PUSHREAL8(prodk(k))
    prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
    CALL PUSHINTEGER4(k)
    DO k=kvts,kvte
      IF (prodk(k) .LT. 0) THEN
        CALL PUSHINTEGER4(2)
      ELSE
        CALL PUSHINTEGER4(1)
      END IF
    END DO
    DO k=kvte,kvts,-1
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 2) prodkb(k) = 0.0
    END DO
    CALL POPINTEGER4(k)
    vtb = 0.0
    factorb = 0.0
    CALL POPREAL8(prodk(k))
    factorb(k) = -(vt(k)*prodk(k)*prodkb(k))
    vtb(k) = -(factor(k)*prodk(k)*prodkb(k))
    prodkb(k) = (1.0-vt(k)*factor(k))*prodkb(k)
    DO k=kvte-1,kvts,-1
      CALL POPREAL8(prodk(k))
      temp = rhok(k+1)*prodk(k+1)
      temp0 = rhok(k)*prodk(k)
      temp0b = -(factor(k)*prodkb(k))
      tempb1 = -(vt(k+1)*temp0b)
      factorb(k) = factorb(k) - (temp0*vt(k)-temp*vt(k+1))*prodkb(k)
      rhokb(k) = rhokb(k) + vt(k)*prodk(k)*temp0b
      vtb(k) = vtb(k) + temp0*temp0b
      rhokb(k+1) = rhokb(k+1) + prodk(k+1)*tempb1
      prodkb(k+1) = prodkb(k+1) + rhok(k+1)*tempb1
      vtb(k+1) = vtb(k+1) - temp*temp0b
      prodkb(k) = vt(k)*rhok(k)*temp0b + prodkb(k)
    END DO
    CALL POPINTEGER4(k)
    pptb = 1000.*rainncv0b
    tempb0 = dtfall*rhok(k)*pptb/rhowat
    rhokb(k) = rhokb(k) + dtfall*prodk(k)*vt(k)*pptb/rhowat
    prodkb(k) = prodkb(k) + vt(k)*tempb0
    vtb(k) = vtb(k) + prodk(k)*tempb0
    factorb(kvte) = 0.0
    DO k=kvte-1,kvts,-1
      rhokb(k) = rhokb(k) - dtfall*rdzk(k)*factorb(k)/rhok(k)**2
      factorb(k) = 0.0
    END DO
    vtdenb = 0.0
    DO k=kvte,kvts,-1
      CALL POPINTEGER4(branch)
      IF (branch .LT. 2) THEN
        qrrb = 36.34*vtden(k)*0.1364*qrr**(-0.8636)*vtb(k)
        vtdenb(k) = vtdenb(k) + 36.34*qrr**0.1364*vtb(k)
        vtb(k) = 0.0
      ELSE
        vtb(k) = 0.0
        qrrb = 0.0
      END IF
      IF (arg1 .EQ. 0.0) THEN
        arg1b = 0.0
      ELSE
        arg1b = vtdenb(k)/(2.0*SQRT(arg1))
      END IF
      vtdenb(k) = 0.0
      CALL POPREAL8(arg1)
      tempb = arg1b/rhok(k)
      rhokb(1) = rhokb(1) + tempb
      rhokb(k) = rhokb(k) + 0.001*prodk(k)*qrrb - rhok(1)*tempb/rhok(k)
      CALL POPREAL8(qrr)
      prodkb(k) = prodkb(k) + 0.001*rhok(k)*qrrb
    END DO
    DO k=kvte,kvts,-1
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 2) prodkb(k) = 0.0
    END DO
  END SUBROUTINE RFALL_B

!  Differentiation of autoca in reverse (adjoint) mode:
!   gradient, with respect to input variables: qc1d qr1d
!   of linear combination of output variables: qc1d qr1d
  SUBROUTINE AUTOCA_B(qc1d, qc1db, qr1d, qr1db, kvts, kvte, c1, c2, c3, &
&    c4, dt)
    IMPLICIT NONE
!     print*,k,qrprod
    INTEGER :: kvts, kvte, k
    REAL, DIMENSION(kvts:kvte) :: qc1d, qr1d
    REAL, DIMENSION(kvts:kvte) :: qc1db, qr1db
    REAL :: c1, c2, c3, c4
    REAL :: qrrc, dt, factorn, qrprod, qrprod2
    REAL :: factornb, qrprodb, qrprod2b
    REAL :: pwr1
    REAL :: pwr1b
    INTEGER :: branch
    REAL :: temp0b
    REAL :: temp
    qrrc = 1.0e-5
    DO k=kvts,kvte
      IF (qr1d(k) .LT. 0.0) THEN
        qr1d(k) = 0.0
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      IF (qc1d(k) .LT. 0.0) THEN
        qc1d(k) = 0.0
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      IF (qr1d(k) .GE. qrrc) THEN
        CALL PUSHREAL8(pwr1)
        pwr1 = qr1d(k)**c4
        CALL PUSHREAL8(factorn)
        factorn = 1.0/(1.+c3*dt*pwr1)
        CALL PUSHINTEGER4(0)
      ELSE
        CALL PUSHREAL8(factorn)
        factorn = 1.0
        CALL PUSHINTEGER4(1)
      END IF
      qrprod = qc1d(k)*(1.0-factorn)
      qrprod2 = 0.0
      IF (qc1d(k) - c2 .GT. 0) THEN
        qrprod2 = factorn*c1*dt*(qc1d(k)-c2)
        IF (qrprod2 .GT. qc1d(k) - c2) THEN
          qrprod2 = qc1d(k) - c2
          CALL PUSHINTEGER4(2)
        ELSE
          CALL PUSHINTEGER4(1)
        END IF
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
!        print*,k,qrprod2
      qrprod = qrprod + qrprod2
      IF (qc1d(k) - qrprod .GT. 0) THEN
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(2)
      END IF
    END DO
    DO k=kvte,kvts,-1
      CALL POPINTEGER4(branch)
      IF (branch .LT. 2) THEN
        qrprodb = qr1db(k) - qc1db(k)
      ELSE
        qrprodb = qr1db(k)
        qc1db(k) = 0.0
        qrprodb = 0.0
      END IF
      qrprod2b = qrprodb
      CALL POPINTEGER4(branch)
      IF (branch .LT. 2) THEN
        IF (branch .LT. 1) THEN
          factornb = 0.0
          GOTO 100
        END IF
      ELSE
        qc1db(k) = qc1db(k) + qrprod2b
        qrprod2b = 0.0
      END IF
      temp0b = c1*dt*qrprod2b
      factornb = (qc1d(k)-c2)*temp0b
      qc1db(k) = qc1db(k) + factorn*temp0b
 100  qc1db(k) = qc1db(k) + (1.0-factorn)*qrprodb
      factornb = factornb - qc1d(k)*qrprodb
      CALL POPINTEGER4(branch)
      IF (branch .LT. 1) THEN
        CALL POPREAL8(factorn)
        temp = c3*dt*pwr1 + 1.
        pwr1b = -(c3*dt*factornb/temp**2)
        CALL POPREAL8(pwr1)
        IF (.NOT.(qr1d(k) .LE. 0.0 .AND. (c4 .EQ. 0.0 .OR. c4 .NE. INT(&
&            c4)))) qr1db(k) = qr1db(k) + c4*qr1d(k)**(c4-1)*pwr1b
      ELSE
        CALL POPREAL8(factorn)
      END IF
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) qc1db(k) = 0.0
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) qr1db(k) = 0.0
    END DO
  END SUBROUTINE AUTOCA_B

!  Differentiation of satadj in reverse (adjoint) mode:
!   gradient, with respect to input variables: qc qr qv p1d rhok
!                tmp pii
!   of linear combination of output variables: qc qr qv p1d rhok
!                tmp pii
  SUBROUTINE SATADJ_B(qv, qvb, qc, qcb, qr, qrb, tmp, tmpb, p1d, p1db, &
&    pii, piib, rhok, rhokb, kvts, kvte, xlv, dt, cp, ep2, svp1, svp2, &
&    svp3, svpt0)
    IMPLICIT NONE
    INTEGER :: kvts, kvte, k
    REAL, DIMENSION(kvts:kvte) :: qv, qc, qr, tmp, p1d, pii, rhok
    REAL, DIMENSION(kvts:kvte) :: qvb, qcb, qrb, tmpb, p1db, piib, rhokb
    REAL, DIMENSION(kvts:kvte) :: rcgs, pressure, temp, es, qvs
    REAL, DIMENSION(kvts:kvte) :: rcgsb, pressureb, tempb, esb, qvsb
    REAL, DIMENSION(kvts:kvte) :: ern, qv2cl, rn2qv
    REAL, DIMENSION(kvts:kvte) :: ernb, qv2clb, rn2qvb
! local var
    REAL :: svp1, svp2, svp3, svpt0, ep2, xlv, cp, dt, f5
    REAL :: ernmax, product
    REAL :: ernmaxb, productb
    REAL :: arg1
    REAL :: arg1b
    INTRINSIC EXP
    INTEGER :: branch
    REAL :: temp3
    REAL :: temp2
    REAL :: temp1
    REAL :: temp0
    REAL :: temp13b
    REAL :: temp7b
    REAL :: temp13b0
    REAL :: temp0b
    REAL :: temp6b
    REAL :: temp12
    REAL :: temp11
    REAL :: temp10
    REAL :: temp9b
    REAL :: temp0b3
    REAL :: temp0b2
    REAL :: temp0b1
    REAL :: temp0b0
    REAL :: temp2b
    REAL :: temp5b
    REAL :: temp8b
    REAL :: temp1b
    REAL :: temp9
    REAL :: temp8
    REAL :: temp7
    REAL :: temp6
    REAL :: temp4b
    REAL :: temp5
    REAL :: temp4
    f5 = svp2*(svpt0-svp3)*xlv/cp
    DO k=kvts,kvte
!constant
      rcgs(k) = 0.001*rhok(k)
      pressure(k) = p1d(k)
      temp(k) = pii(k)*tmp(k)
      CALL PUSHREAL8(arg1)
      arg1 = svp2*(temp(k)-svpt0)/(temp(k)-svp3)
      es(k) = 1000.*svp1*EXP(arg1)
      qvs(k) = ep2*es(k)/(pressure(k)-es(k))
      IF (qr(k) .LT. 0) THEN
        qr(k) = 0.0
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      IF (qv(k) .LT. 0) THEN
        qv(k) = 0.0
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      IF (qc(k) .LT. 0) THEN
        qc(k) = 0.0
        CALL PUSHINTEGER4(2)
      ELSE
        CALL PUSHINTEGER4(1)
      END IF
    END DO
    DO k=kvts,kvte
!not related to time; maximum transform qv to cl (sat) or cl to qv (sub sat)
      qv2cl(k) = (qv(k)-qvs(k))/(1.+pressure(k)/(pressure(k)-es(k))*qvs(&
&        k)*f5/(temp(k)-svp3)**2)
! sub sat rain evaperate
      ern(k) = 0.0
      IF (qvs(k) .GT. qv(k)) THEN
        IF (qr(k) .GE. 1d-5) THEN
          rn2qv(k) = dt*((1.6+124.9*(rcgs(k)*qr(k))**.2046)*(rcgs(k)*qr(&
&            k))**.525/(2.55e8/(pressure(k)*qvs(k))+5.4e5))*((qvs(k)-qv(k&
&            ))/(rcgs(k)*qvs(k)))
          CALL PUSHINTEGER4(0)
        ELSE
          rn2qv(k) = 0.0
          CALL PUSHINTEGER4(1)
        END IF
        IF (rn2qv(k) .GT. qr(k)) THEN
          rn2qv(k) = qr(k)
          CALL PUSHINTEGER4(1)
        ELSE
          CALL PUSHINTEGER4(0)
        END IF
        ernmax = 0.0
        IF (-qv2cl(k) - qc(k) .GT. 0.0) THEN
          ernmax = -qv2cl(k) - qc(k)
          CALL PUSHINTEGER4(1)
        ELSE
          CALL PUSHINTEGER4(0)
        END IF
!        ern(k)  = amin1(rn2qv(k), ernmax)
        ern(k) = rn2qv(k)
        IF (rn2qv(k) .GT. ernmax) THEN
          ern(k) = ernmax
          CALL PUSHINTEGER4(2)
        ELSE
          CALL PUSHINTEGER4(1)
        END IF
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
! Update all variables
!       product = amax1(qv2cl(k),-qc(k))
      product = qv2cl(k)
      IF (qv2cl(k) .LT. -qc(k)) THEN
        product = -qc(k)
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      CALL PUSHREAL8(qv(k))
!       qv(k) = amax1(qv(k) - product + ern(k),0.)
      qv(k) = qv(k) - product + ern(k)
      IF (qv(k) .LT. 0) THEN
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      CALL PUSHREAL8(temp(k))
      temp(k) = temp(k) + xlv/cp*(product-ern(k))
    END DO
    ernb = 0.0
    tempb = 0.0
    rcgsb = 0.0
    pressureb = 0.0
    qv2clb = 0.0
    esb = 0.0
    qvsb = 0.0
    rn2qvb = 0.0
    DO k=kvte,kvts,-1
      temp13b = tmpb(k)/pii(k)
      tempb(k) = tempb(k) + temp13b
      piib(k) = piib(k) - temp(k)*temp13b/pii(k)
      tmpb(k) = 0.0
      CALL POPREAL8(temp(k))
      temp13b0 = xlv*tempb(k)/cp
      productb = qcb(k) + temp13b0
      ernb(k) = ernb(k) - qrb(k) - temp13b0
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) qvb(k) = 0.0
      CALL POPREAL8(qv(k))
      productb = productb - qvb(k)
      ernb(k) = ernb(k) + qvb(k)
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) THEN
        qcb(k) = qcb(k) - productb
        productb = 0.0
      END IF
      qv2clb(k) = qv2clb(k) + productb
      CALL POPINTEGER4(branch)
      IF (branch .LT. 2) THEN
        IF (branch .LT. 1) THEN
          GOTO 100
        ELSE
          ernmaxb = 0.0
        END IF
      ELSE
        ernmaxb = ernb(k)
        ernb(k) = 0.0
      END IF
      rn2qvb(k) = rn2qvb(k) + ernb(k)
      ernb(k) = 0.0
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) THEN
        qv2clb(k) = qv2clb(k) - ernmaxb
        qcb(k) = qcb(k) - ernmaxb
      END IF
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) THEN
        qrb(k) = qrb(k) + rn2qvb(k)
        rn2qvb(k) = 0.0
      END IF
      CALL POPINTEGER4(branch)
      IF (branch .LT. 1) THEN
        temp12 = rcgs(k)*qvs(k)
        temp4 = pressure(k)*qvs(k)
        temp11 = 2.55e8/temp4
        temp5 = (temp11+5.4e5)*temp12
        temp6 = rcgs(k)*qr(k)
        temp10 = 124.9*temp6**.2046 + 1.6
        temp7 = temp10/temp5
        temp8 = rcgs(k)*qr(k)
        temp9 = temp8**.525
        temp9b = dt*temp7*rn2qvb(k)
        temp8b = (qvs(k)-qv(k))*.525*temp8**(-0.475)*temp9b
        temp7b = dt*temp9*(qvs(k)-qv(k))*rn2qvb(k)/temp5
        temp6b = 124.9*.2046*temp6**(-0.7954)*temp7b
        temp5b = -(temp7*temp7b)
        temp4b = -(temp12*temp11*temp5b/temp4)
        rcgsb(k) = rcgsb(k) + (temp11+5.4e5)*qvs(k)*temp5b + qr(k)*&
&          temp6b + qr(k)*temp8b
        qrb(k) = qrb(k) + rcgs(k)*temp6b + rcgs(k)*temp8b
        qvsb(k) = qvsb(k) + (temp11+5.4e5)*rcgs(k)*temp5b + pressure(k)*&
&          temp4b + temp9*temp9b
        qvb(k) = qvb(k) - temp9*temp9b
        pressureb(k) = pressureb(k) + qvs(k)*temp4b
        rn2qvb(k) = 0.0
      ELSE
        rn2qvb(k) = 0.0
      END IF
 100  ernb(k) = 0.0
      rn2qvb(k) = 0.0
      temp3 = (temp(k)-svp3)**2
      temp0 = (pressure(k)-es(k))*temp3
      temp2 = pressure(k)*qvs(k)
      temp1 = temp2/temp0
      temp2b = qv2clb(k)/(f5*temp1+1.)
      temp1b = -((qv(k)-qvs(k))*f5*temp2b/((f5*temp1+1.)*temp0))
      temp0b2 = -(temp1*temp1b)
      temp0b3 = temp3*temp0b2
      qvb(k) = qvb(k) + temp2b
      qvsb(k) = qvsb(k) + pressure(k)*temp1b - temp2b
      pressureb(k) = pressureb(k) + temp0b3 + qvs(k)*temp1b
      esb(k) = esb(k) - temp0b3
      tempb(k) = tempb(k) + (pressure(k)-es(k))*2*(temp(k)-svp3)*temp0b2
      qv2clb(k) = 0.0
    END DO
    DO k=kvte,kvts,-1
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 2) qcb(k) = 0.0
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) qvb(k) = 0.0
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) qrb(k) = 0.0
      temp0b = ep2*qvsb(k)/(pressure(k)-es(k))
      temp0b0 = -(es(k)*temp0b/(pressure(k)-es(k)))
      esb(k) = esb(k) + temp0b - temp0b0
      pressureb(k) = pressureb(k) + temp0b0
      qvsb(k) = 0.0
      arg1b = svp1*1000.*EXP(arg1)*esb(k)
      esb(k) = 0.0
      CALL POPREAL8(arg1)
      temp0b1 = svp2*arg1b/(temp(k)-svp3)
      tempb(k) = tempb(k) + (1.0-(temp(k)-svpt0)/(temp(k)-svp3))*temp0b1
      piib(k) = piib(k) + tmp(k)*tempb(k)
      tmpb(k) = tmpb(k) + pii(k)*tempb(k)
      tempb(k) = 0.0
      p1db(k) = p1db(k) + pressureb(k)
      pressureb(k) = 0.0
      rhokb(k) = rhokb(k) + 0.001*rcgsb(k)
      rcgsb(k) = 0.0
    END DO
  END SUBROUTINE SATADJ_B

END MODULE A_MODULE_MP_MKESSLER
