      module module_ras
      USE MACHINE , ONLY : kind_phys
      use physcons, grav => con_g, cp => con_cp, alhl => con_hvap       &
     &,             alhf => con_hfus, rgas => con_rd, rkap => con_rocp  &
     &,             nu => con_FVirt
      implicit none
      SAVE
!
      integer, parameter :: nrcmax=32 ! Maximum # of random clouds per 1200s

      real (kind=kind_phys), parameter :: delt_c=1800.0/3600.0          &
!     Adjustment time scales in hrs for deep and shallow clouds
!    &,                                   adjts_d=3.0, adjts_s=0.5
!    &,                                   adjts_d=2.5, adjts_s=0.5
     &,                                   adjts_d=2.0, adjts_s=0.5
!
      logical,               parameter :: fix_ncld_hr=.true.
!
      real(kind=kind_phys), parameter  :: ZERO=0.0, HALF=0.5,  ONE=1.0, &
     &                                    TWO=2.0,  FOUR_P2=4.E2,       &
     &                                    FOUR=4.,  ONE_M10=1.E-10,     &
     &                                    ONE_M6=1.E-6, ONE_M5=1.E-5,   &
     &                                    ONE_M2=1.E-2, ONE_M1=1.E-1
!
      real(kind=kind_phys), parameter  :: cmb2pa = 100.0  ! Conversion from Mb to Pa
!
      real(kind=kind_phys), parameter  ::                               &
     &           ONEBG   = ONE / GRAV,    GRAVCON = cmb2pa * ONEBG      &
     &,          GRAVFAC = GRAV / CMB2PA, ELOCP   = ALHL / CP           &
     &,          ELFOCP  = (ALHL+ALHF) / CP                             &
     &,          RKAPI   = ONE / RKAP,    RKPP1I  = ONE / (ONE+RKAP)    &
     &,          CMPOR   = CMB2PA / RGAS                                &
     &,          zfac    = 0.28888889E-4 * ONEBG
!
!     logical, parameter :: advcld=.true., advups=.true., advtvd=.false.
      logical, parameter :: advcld=.true., advups=.false., advtvd=.true.
!     logical, parameter :: advcld=.true., advups=.false.,advtvd=.false.
!
      real(kind=kind_phys)  RHMAX,  qudfac, QUAD_LAM, RHRAM, TESTMB,    &
     &                      TSTMBI, HCRITD, DD_DP,    RKNOB,  AFC, EKNOB&
     &,                     shalfac,HCRITS, HPERT_FAC

!     PARAMETER (DD_DP=1000.0, RKNOB=1.0, EKNOB=1.0)   ! No downdraft!
      PARAMETER (DD_DP=500.0,  RKNOB=1.0, EKNOB=1.0)
!     PARAMETER (DD_DP=500.0,  RKNOB=2.0, EKNOB=1.0)
!
      PARAMETER (RHMAX=1.0   )   ! MAX RELATIVE HUMIDITY
      PARAMETER (QUAD_LAM=1.0)   ! MASK FOR QUADRATIC LAMBDA
!     PARAMETER (RHRAM=0.15)     ! PBL RELATIVE HUMIDITY RAMP
      PARAMETER (RHRAM=0.05)     ! PBL RELATIVE HUMIDITY RAMP
      PARAMETER (HCRITD=4000.0)  ! Critical Moist Static Energy
      PARAMETER (HCRITS=2000.0)  ! Critical Moist Static Energy

!     parameter (hpert_fac=1.01) ! Perturbation on hbl when ctei=.true.
!     parameter (hpert_fac=1.005)! Perturbation on hbl when ctei=.true.
!     parameter (hpert_fac=1.00) ! Perturbation on hbl when ctei=.true.
!     parameter (qudfac=quad_lam*half, shalfac=1.0)
!     parameter (qudfac=quad_lam*half, shalfac=2.0)
      parameter (qudfac=quad_lam*half, shalfac=3.0)
!     parameter (qudfac=quad_lam*0.25)    ! Yogesh's
      parameter (testmb=0.1, tstmbi=one/testmb)
!
      real(kind=kind_phys) ALMIN1, ALMIN2, ALMAX
      real(kind=kind_phys) facdt
!
!     PARAMETER (ALMIN1=0.00E-6, ALMIN2=2.50E-5, ALMAX=1.0E-2)
      PARAMETER (ALMIN1=0.00E-6, ALMIN2=4.00E-5, ALMAX=1.0E-2)
!     PARAMETER (ALMIN1=0.00E-6, ALMIN2=4.00E-5, ALMAX=1.0E-2)
!cnt  PARAMETER (ALMIN1=0.00E-6, ALMIN2=2.50E-5, ALMAX=5.0E-3)
!
!     real(kind=kind_phys), parameter :: BLDMAX = 200.0
      real(kind=kind_phys), parameter :: BLDMAX = 300.0
!!    real(kind=kind_phys), parameter :: BLDMAX = 350.0
!
      real(kind=kind_phys) C0, C0I, QI0, QW0, c00, c00i, dlq_fac
      PARAMETER (QI0=1.0E-5, QW0=1.0E-5)
!     PARAMETER (QI0=1.0E-4, QW0=1.0E-5) ! 20050509
!     PARAMETER (QI0=1.0E-5, QW0=1.0E-6)
!!!   PARAMETER (C0I=1.0E-3)
      PARAMETER (C00I=1.0E-3)
!     parameter (c0=1.0e-3)
!     parameter (c0=1.5e-3)
!!!   parameter (c0=2.0e-3)
      parameter (c00=2.0e-3)
!
      real(kind=kind_phys) TF, TCR, TCRF, TCL
!     parameter (TF=130.16, TCR=160.16, TCRF=1.0/(TCR-TF),TCL=2.0)
!     parameter (TF=230.16, TCR=260.16, TCRF=1.0/(TCR-TF))
      parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF),TCL=2.0)
!
!     For Tilting Angle Specification
!
      real(kind=kind_phys) REFP(6), REFR(6), TLAC(8), PLAC(8), TLBPL(7) &
     &,                    drdp(5), VTP
!
      DATA PLAC/100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0/
      DATA TLAC/ 35.0,  25.0,  20.0,  17.5,  15.0,  12.5,  10.0,  7.5/
      DATA REFP/500.0, 300.0, 250.0, 200.0, 150.0, 100.0/
      DATA REFR/ 1.0,   2.0,  3.0,   4.0,   6.0,   8.0/
!
      real(kind=kind_phys) AC(16), AD(16)
!
      integer, parameter :: nqrp=500001
      real(kind=kind_phys) C1XQRP, C2XQRP, TBQRP(NQRP), TBQRA(NQRP)     &
     &,                    TBQRB(NQRP)
!
      integer, parameter :: nvtp=10001
      real(kind=kind_phys) C1XVTP, C2XVTP, TBVTP(NVTP)
!
      contains
!
      subroutine set_ras_afc(dt)
      implicit none
      real(kind=kind_phys) DT
!     AFC = -(1.04E-4*DT)*(3600./DT)**0.578
      AFC = -(1.01097E-4*DT)*(3600./DT)**0.57777778
      end subroutine set_ras_afc

      subroutine ras_init(levs,  me)
!
      Implicit none
!
      integer levs
!
      real(kind=kind_phys) actp,   facm, tem,  actop, tem1, tem2
      integer              i, l, me
      PARAMETER (ACTP=1.7,   FACM=1.00)
!
      real(kind=kind_phys) PH(15),    A(15)
!
      DATA PH/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0    &
     &,       550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/
!
       DATA A/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677           &
     &,       0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664            &
     &,       0.0553, 0.0445, 0.0633/
!
      logical first
      data first/.true./
!
      if (first) then
!                                   set critical workfunction arrays
        ACTOP = ACTP*FACM
        DO L=1,15
          A(L) = A(L)*FACM
        ENDDO
        DO L=2,15
          TEM   = 1.0 / (PH(L) - PH(L-1))
          AC(L) = (PH(L)*A(L-1) - PH(L-1)*A(L)) * TEM
          AD(L) = (A(L) - A(L-1)) * TEM
        ENDDO
        AC(1)  = ACTOP
        AC(16) = A(15)
        AD(1)  = 0.0
        AD(16) = 0.0
!
        CALL SETQRP
        CALL SETVTP
!
        do i=1,7
          tlbpl(i) = (tlac(i)-tlac(i+1)) / (plac(i)-plac(i+1))
        enddo
        do i=1,5
          drdp(i)  = (REFR(i+1)-REFR(i)) / (REFP(i+1)-REFP(i))
        enddo
!
        VTP    = 36.34*SQRT(1.2)* (0.001)**0.1364
!
        if (me .eq. 0) print *,' NO DOWNDRAFT FOR CLOUD TYPES'          &
     &,        ' DETRAINING WITHIN THE BOTTOM ',DD_DP,' hPa LAYERS'
!
        first = .false.
      endif
!
      end subroutine ras_init
      end module module_ras
!
      module module_rascnv
!
      USE MACHINE , ONLY : kind_phys
      implicit none
      SAVE
!
      logical REVAP, CUMFRC
      LOGICAL WRKFUN, CALKBL, CRTFUN, UPDRET, BOTOP, vsmooth

      real(kind=kind_phys), parameter :: frac=0.5,    crtmsf=0.0        &
     &,                                  rhfacs=0.70, rhfacl=0.70       &
     &,                                  face=5.0,    delx=10000.0      &
     &,                                  ddfac=face*delx*0.001          &
     &,                                  max_neg_bouy=0.15
!    &,                                  max_neg_bouy=0.25

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!    real(kind=kind_phys) FRAC, CRTMSF, MAX_NEG_BOUY, rhfacs, rhfacl   &
!!   &,                    FACE, DELX,   DDFAC
!     parameter (frac=0.1, crtmsf=0.0)
!     parameter (frac=0.25, crtmsf=0.0)
!!    parameter (frac=0.5, crtmsf=0.0)
!     PARAMETER (MAX_NEG_BOUY=0.15, REVAP=.true., CUMFRC=.false.)
!     PARAMETER (MAX_NEG_BOUY=0.15, REVAP=.true., CUMFRC=.true.)
!     PARAMETER (MAX_NEG_BOUY=0.10, REVAP=.true., CUMFRC=.true.)
!     PARAMETER (MAX_NEG_BOUY=0.20, REVAP=.true., CUMFRC=.true.)
!!    PARAMETER (MAX_NEG_BOUY=0.25, REVAP=.true., CUMFRC=.true.)
!     PARAMETER (MAX_NEG_BOUY=0.30, REVAP=.true., CUMFRC=.true.)
!!    PARAMETER (MAX_NEG_BOUY=0.05, REVAP=.true., CUMFRC=.true.)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      PARAMETER (                   REVAP=.true., CUMFRC=.true.)
      PARAMETER (WRKFUN = .FALSE.,  UPDRET = .FALSE., vsmooth=.false.)
!     PARAMETER (CRTFUN = .TRUE.,   CALKBL = .false., BOTOP=.true.)
      PARAMETER (CRTFUN = .TRUE.,   CALKBL = .true., BOTOP=.true.)
!
!!    parameter (rhfacs=0.70, rhfacl=0.70)
!     parameter (rhfacs=0.75, rhfacl=0.75)
!     parameter (rhfacs=0.85, rhfacl=0.85)
!     parameter (rhfacs=0.80, rhfacl=0.80)   ! August 26, 2008
!     parameter (rhfacs=0.80, rhfacl=0.85)
!!    PARAMETER (FACE=5.0, DELX=10000.0, DDFAC=FACE*DELX*0.001)
!
!     real (kind=kind_phys), parameter :: pgftop=0.7, pgfbot=0.3        &
!     real (kind=kind_phys), parameter :: pgftop=0.75, pgfbot=0.35      &
!    For pressure gradient force in momentum mixing
!     real (kind=kind_phys), parameter :: pgftop=0.80, pgfbot=0.30      &
!    No pressure gradient force in momentum mixing
      real (kind=kind_phys), parameter :: pgftop=0.0, pgfbot=0.0        &
!     real (kind=kind_phys), parameter :: pgftop=0.55, pgfbot=0.55      &
     &,          pgfgrad=(pgfbot-pgftop)*0.001
!
      end module module_rascnv
!
!
      subroutine rascnv(IM,    IX,     k,      dt,    dtf,  rannum      &
     &,                 tin,   qin,    uin,    vin,   ccin,  trac       &
     &,                 prsi,  prsl,   prsik,  prslk, phil,  phii       &
     &,                 KPBL,  CDRAG,  RAINC,  kbot,  ktop,  kcnv       &
     &,                 DDVEL, FLIPV,  facmb,  me,    garea, lmh, ccwfac&
     &,                 nrcm,  rhc,    ud_mf, dd_mf,  det_mf, dlqfac    &
     &,                 lprnt, ipr, kdt, fscav)
!    &,                 lprnt, ipr, kdt, fscav, ctei_r, ctei_rm)
!
!*********************************************************************
!*********************************************************************
!************         Relaxed Arakawa-Schubert      ******************
!************             Parameterization          ******************
!************          Plug Compatible Driver       ******************
!************               23 May 2002             ******************
!************                                       ******************
!************               Developed By            ******************
!************                                       ******************
!************             Shrinivas Moorthi         ******************
!************                                       ******************
!************                  EMC/NCEP             ******************
!*********************************************************************
!*********************************************************************
!
!
      USE MACHINE , ONLY : kind_phys
      use module_ras, DPD => DD_DP
      use module_rascnv
      Implicit none
!
      LOGICAL FLIPV, lprnt
!
!      input
!
      Integer IM, IX, k, ncrnd, me, trac, ipr, nrcm
      Integer kbot(im), ktop(im), kcnv(im), KPBL(im), lmh(im), kdt
!
      real(kind=kind_phys) tin(ix,k),     qin(ix,k),  uin(ix,k)         &
     &,                    vin(ix,k),     prsi(ix,k+1)                  &
     &,                    prsik(ix,k+1), prsl(ix,k), prslk(ix,k+1)     &
     &,                    phil(ix,k),    phii(ix,k+1),ccwfac(im)       &
     &,                    ccin(ix,k,trac+2)                            &
!    &,                    prsik(ix,k+1), clt(ix,k)                     &
     &,                    RAINC(im),     CDRAG(im),  DDVEL(im)         &
     &,                    rannum(ix,nrcm),dlqfac                       &
     &,                    ud_mf(im,k), dd_mf(im,k), det_mf(im,k)
      real(kind=kind_phys) DT, facmb, garea(im), dtf, rhc(im,k)         &
!
!     Added for aerosol scavenging for GOCART
!
      real(kind=kind_phys), intent(in) :: fscav(trac)

!    &,                                   ctei_r(im), ctei_rm
!
!     locals
!
!     real(kind=kind_phys) RAIN,     toi(k), qoi(k),  uvi(k,trac+2)     &
      real(kind=kind_phys) RAIN,     toi(k), qoi(k)                     &
     &,                    TCU(k),   QCU(k), PCU(k),   clw(k), cli(k)   &
     &,                    QII(k),   QLI(k), PRS(k+1), PSJ(k+1)         &
     &,                    phi_l(k), phi_h(k+1)                         &
     &,                    wfnc,     flx(k+1), FLXD(K+1)
      real(kind=kind_phys) daylen,pfac,tla,pl,clwmin
      integer icm,irnd,ib

      PARAMETER (ICM=100, DAYLEN=86400.0, PFAC=1.0/450.0,clwmin=1.0e-10)
      Integer  IC(ICM)
!
      real(kind=kind_phys), allocatable ::  ALFINT(:,:), uvi(:,:)
     &,                                     trcfac(:,:), rcu(:,:)
      real(kind=kind_phys)            ALFINQ(K),   PRSM(K),  PSJM(K)    &
     &,                    alfind(K), rhc_l(k), dtvd(2,4)
!    &,                    DPI(K),    psjp(k+1)              
      real(kind=kind_phys) CFAC, TEM, dpi, sgc, ccwf, tem1, tem2
!
      Integer              KCR,  KFX, NCMX, NC,  KTEM, I,   L, lm1      &
     &,                    ntrc, ia,  ll,   km1, kp1,  ipt, lv, KBL, n  &
     &,                    lmhij, KRMIN, KRMAX, KFMAX, kblmx
      real(kind=kind_phys) sgcs(k,im)
!
      LOGICAL  DNDRFT, lprint
!     LOGICAL  DNDRFT, lprint, ctei
!
!  Scavenging related parameters
!
      real                fscav_(trac+2)  ! Fraction scavenged per km
!
      fscav_ = 0.0                        ! By default no scavenging
      if (trac > 0) then
        do i=1,trac
          fscav_(i) = fscav(i)
        enddo
      endif
!
      km1    = k - 1
      kp1    = k + 1
!
      dlq_fac = dlqfac
      tem     = 1.0 + dlq_fac
      c0      = c00  * tem
      c0i     = c00i * tem
!
      ntrc = trac
      IF (CUMFRC) THEN
        ntrc = ntrc + 2
      ENDIF
      if (ntrc > 0) then
        if (.not. allocated(trcfac)) allocate (trcfac(k,ntrc))
        if (.not. allocated(uvi)) allocate (uvi(k,ntrc))
        if (.not. allocated(rcu)) allocate (rcu(k,ntrc))
        do n=1, ntrc
          do l=1,k
            trcfac(l,n) = 1.0         !  For other tracers
          enddo
        enddo
      endif
!
      if (.not. allocated(alfint)) allocate(alfint(k,ntrc+4))
!
      call set_ras_afc(dt)
!
      DO IPT=1,IM

        ccwf = 0.5
        if (ccwfac(ipt) >= 0.0) ccwf = ccwfac(ipt)

!
!       ctei = .false.
!       if (ctei_r(ipt) > ctei_rm) ctei = .true.
!

        do l=1,k
          ud_mf(ipt,l)  = 0.0
          dd_mf(ipt,l)  = 0.0
          det_mf(ipt,l) = 0.0
        enddo
!
!     Compute NCRND  : here LMH is the number of layers above the
!                      bottom surface.  For sigma coordinate LMH=K.
!                      if flipv is true, then input variables are from bottom
!                      to top while RAS goes top to bottom
!
        LMHIJ = LMH(ipt)
        if (flipv) then
           ll  = kp1 - LMHIJ
           tem = 1.0 / prsi(ipt,ll)
        else
           ll  = LMHIJ
           tem = 1.0 / prsi(ipt,ll+1)
        endif
        KRMIN = 1
        KRMAX = km1
        KFMAX = KRMAX
        kblmx = 1
        DO L=1,LMHIJ-1
          ll = l
          if (flipv) ll = kp1 -l ! Input variables are bottom to top!
          SGC = prsl(ipt,ll) * tem
          sgcs(l,ipt) = sgc
          IF (SGC .LE. 0.050) KRMIN = L
!         IF (SGC .LE. 0.700) KRMAX = L
!         IF (SGC .LE. 0.800) KRMAX = L
          IF (SGC .LE. 0.760) KRMAX = L
!         IF (SGC .LE. 0.930) KFMAX = L
          IF (SGC .LE. 0.970) KFMAX = L    ! Commented on 20060202
!         IF (SGC .LE. 0.700) kblmx = L    ! Commented on 20101015
          IF (SGC .LE. 0.600) kblmx = L    ! 
!         IF (SGC .LE. 0.650) kblmx = L    ! Commented on 20060202
        ENDDO
        krmin = max(krmin,2)

!     if (lprnt .and. ipt .eq. ipr) print *,' krmin=',krmin,' krmax=',
!    &krmax,' kfmax=',kfmax,' lmhij=',lmhij,' tem=',tem
!
        if (fix_ncld_hr) then
!!!       NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.50001
          NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.10001
!!        NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/600) + 0.50001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/360) + 0.50001
!    &                                         + 0.50001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * min(1.0,DTF/360) + 0.1
          facdt = delt_c / dt
        else
          NCRND = min(nrcmax, (KRMAX-KRMIN+1))
          facdt = 1.0 / 3600.0
        endif
        NCRND   = min(nrcm,max(NCRND, 1))
!
        KCR     = MIN(LMHIJ,KRMAX)
        KTEM    = MIN(LMHIJ,KFMAX)
        KFX     = KTEM - KCR

!     if(lprnt)print*,' enter RASCNV k=',k,' ktem=',ktem,' LMHIJ='
!    &,                 LMHIJ
!    &,               ' krmax=',krmax,' kfmax=',kfmax
!    &,               ' kcr=',kcr, ' cdrag=',cdrag(ipr)
 
        IF (KFX .GT. 0) THEN
           IF (BOTOP) THEN
              DO NC=1,KFX
                IC(NC) = KTEM + 1 - NC
              ENDDO
           ELSE
              DO NC=KFX,1,-1
               IC(NC) = KTEM + 1 - NC
              ENDDO
           ENDIF
        ENDIF
!
        NCMX  = KFX + NCRND
        IF (NCRND .GT. 0) THEN
          DO I=1,NCRND
            IRND = (RANNUM(ipt,I)-0.0005)*(KCR-KRMIN+1)
            IC(KFX+I) = IRND + KRMIN
          ENDDO
        ENDIF
!
!     ia = 1
!
!     print *,' in rascnv: k=',k,'lat=',lat,' lprnt=',lprnt
!     if (lprnt) then
!        if (me .eq. 0) then
!        print *,' tin',(tin(ia,l),l=k,1,-1)
!        print *,' qin',(qin(ia,l),l=k,1,-1)
!     endif
!
!
        lprint = lprnt .and. ipt .eq. ipr
!       lprint = lprnt
        do l=1,k
          ll = l
          if (flipv) ll = kp1 -l ! Input variables are bottom to top!
          CLW(l)     = 0.0       ! Assumes initial value of Cloud water
          CLI(l)     = 0.0       ! Assumes initial value of Cloud ice
                                 ! to be zero i.e. no environmental condensate!!!
          QII(l)     = 0.0
          QLI(l)     = 0.0
!                          Initialize heating, drying, cloudiness etc.
          tcu(l)     = 0.0
          qcu(l)     = 0.0
          pcu(l)     = 0.0
          flx(l)     = 0.0
          flxd(l)    = 0.0
          rcu(l,1)   = 0.0
          rcu(l,2)   = 0.0
!                          Transfer input prognostic data into local variable
          toi(l)     = tin(ipt,ll)
          qoi(l)     = qin(ipt,ll)
!
          if (ntrc > trac) then               ! CUMFRC is true 
            uvi(l,trac+1) = uin(ipt,ll)
            uvi(l,trac+2) = vin(ipt,ll)
          endif
!
          if (trac > 0) then                  ! tracers such as O3, dust etc
            do n=1,trac
              uvi(l,n) = ccin(ipt,ll,n+2)
              if (abs(uvi(l,n)) < 1.0e-20) uvi(l,n) = 0.0
            enddo
          endif
!
        enddo
        flx(k+1)  = 0.0
        flxd(k+1) = 0.0
!
        if (ccin(ipt,1,2) .le. -999.0) then  ! input ice/water are together 
          do l=1,k
            ll = l
            if (flipv) ll = kp1 -l ! Input variables are bottom to top!
              tem = ccin(ipt,ll,1)                                      &
     &            * MAX(ZERO, MIN(ONE, (TCR-toi(L))*TCRF))
              ccin(ipt,ll,2) = ccin(ipt,ll,1) - tem
              ccin(ipt,ll,1) = tem
          enddo
        endif
        if (advcld) then
          do l=1,k
            ll = l
            if (flipv) ll = kp1 -l ! Input variables are bottom to top!
            QII(L) = ccin(ipt,ll,1)
            QLI(L) = ccin(ipt,ll,2)
          enddo
        endif
!
        KBL  = KPBL(ipt)
        if (flipv) KBL  = MAX(MIN(k, kp1-KPBL(ipt)), k/2)
        rain = 0.0
!
        DO L=1,kp1
          ll = l
          if (flipv) ll = kp1 + 1 - l      ! Input variables are bottom to top!
          PRS(LL)   = prsi(ipt, L) * facmb ! facmb is for conversion to MB
          PSJ(LL)   = prsik(ipt,L)
          phi_h(LL) = phii(ipt,L)
        ENDDO
!
        DO L=1,k
          ll = l
          if (flipv) ll = kp1 - l          ! Input variables are bottom to top!
          PRSM(LL)  = prsl(ipt, L) * facmb ! facmb is for conversion to MB
          PSJM(LL)  = prslk(ipt,L)
          phi_l(LL) = phil(ipt,L)
          rhc_l(LL) = rhc(ipt,L)
        ENDDO
!
!     if (lprnt .and. ipt == ipr) print *,' phi_h=',phi_h(:)
!     if(lprint) print *,' PRS=',PRS
!     if(lprint) print *,' PRSM=',PRSM
!     if (lprint) then
!        print *,' qns=',qns(ia),' qoi=',qn0(ia,k),'qin=',qin(ia,1)
!        if (me .eq. 0) then
!        print *,' toi',(tn0(ia,l),l=1,k)
!        print *,' qoi',(qn0(ia,l),l=1,k),' kbl=',kbl
!     endif
!
!
!       do l=k,kctop(1),-1
!!        DPI(L)  = 1.0 / (PRS(L+1) - PRS(L))
!       enddo
!
!     print *,' ipt=',ipt

        if (advups) then               ! For first order upstream for updraft
          alfint(:,:) = 1.0
        elseif (advtvd) then           ! TVD flux limiter scheme for updraft
          alfint(:,:) = 1.0
          l   = krmin
          lm1 = l - 1
          dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)        &
     &              + alhl*(qoi(l)-qoi(lm1))
          dtvd(1,2) = qoi(l) - qoi(lm1)
          dtvd(1,3) = qli(l) - qli(lm1)
          dtvd(1,4) = qii(l) - qii(lm1)
          do l=krmin+1,k
            lm1 = l - 1

!     print *,' toi=',toi(l),toi(lm1),' phi_l=',phi_l(l),phi_l(lm1)
!    &,' qoi=',qoi(l),qoi(lm1),' cp=',cp,' alhl=',alhl

            dtvd(2,1)   = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)    &
     &                  + alhl*(qoi(l)-qoi(lm1))

!     print *,' l=',l,' dtvd=',dtvd(:,1)

            if (abs(dtvd(2,1)) > 1.0e-10) then
              tem1        = dtvd(1,1) / dtvd(2,1)
              tem2        = abs(tem1)
              alfint(l,1) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for h
            endif

!     print *,' alfint=',alfint(l,1),' l=',l,' ipt=',ipt

            dtvd(1,1)   = dtvd(2,1)
!
            dtvd(2,2)   = qoi(l) - qoi(lm1)

!     print *,' l=',l,' dtvd2=',dtvd(:,2)

            if (abs(dtvd(2,2)) > 1.0e-10) then
              tem1        = dtvd(1,2) / dtvd(2,2)
              tem2        = abs(tem1)
              alfint(l,2) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for q
            endif
            dtvd(1,2)   = dtvd(2,2)
!
            dtvd(2,3)   = qli(l) - qli(lm1)

!     print *,' l=',l,' dtvd3=',dtvd(:,3)

            if (abs(dtvd(2,3)) > 1.0e-10) then
              tem1        = dtvd(1,3) / dtvd(2,3)
              tem2        = abs(tem1)
              alfint(l,3) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for ql
            endif
            dtvd(1,3)   = dtvd(2,3)
!
            dtvd(2,4)   = qii(l) - qii(lm1)

!     print *,' l=',l,' dtvd4=',dtvd(:,4)

            if (abs(dtvd(2,4)) > 1.0e-10) then
              tem1        = dtvd(1,4) / dtvd(2,4)
              tem2        = abs(tem1)
              alfint(l,4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for qi
            endif
            dtvd(1,4)   = dtvd(2,4)
          enddo
!
          if (ntrc > 0) then
            do n=1,ntrc
              l = krmin
              dtvd(1,1)   = uvi(l,n) - uvi(l-1,n)
              do l=krmin+1,k
                dtvd(2,1)     = uvi(l,n) - uvi(l-1,n)

!     print *,' l=',l,' dtvdn=',dtvd(:,1),' n=',n,' l=',l

                if (abs(dtvd(2,1)) > 1.0e-10) then
                  tem1          = dtvd(1,1) / dtvd(2,1)
                  tem2          = abs(tem1)
                  alfint(l,n+4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for tracers
                endif
                dtvd(1,1)     = dtvd(2,1)
              enddo
            enddo
          endif
        else
          alfint(:,:) = 0.5              ! For second order scheme
        endif
        alfind(:)   = 0.5
!
!     print *,' after alfint for ipt=',ipt

! Resolution dependent press grad correction momentum mixing

        if (CUMFRC) then
          do l=krmin,k
            tem = 1.0 - max(pgfbot, min(pgftop, pgftop+pgfgrad*prsm(l)))
            trcfac(l,trac+1) = tem
            trcfac(l,trac+2) = tem
          enddo
        endif
!
!       lprint = lprnt .and. ipt .eq. ipr

!     if (lprint) then
!       print *,' trcfac=',trcfac(krmin:k,1+trac)
!       print *,' alfint=',alfint(krmin:k,1)
!       print *,' alfinq=',alfint(krmin:k,2)
!       print *,' alfini=',alfint(krmin:k,4)
!       print *,' alfinu=',alfint(krmin:k,5)
!     endif
!
        if (calkbl) kbl = k
        DO NC=1,NCMX
!
          IB = IC(NC)
          if (ib > kbl) cycle

!         lprint = lprnt .and. ipt .eq. ipr
!         lprint = lprnt .and. ipt .eq. ipr .and. ib .eq. 41
!
          DNDRFT = DPD > 0.0
!
!         if (lprint) print *,' calling cloud type ib=',ib,' kbl=',kbl
!    *,   ' kpbl=',kpbl,' alfint=',alfint,' frac=',frac
!    *,   ' ntrc=',ntrc,' ipt=',ipt
!
!****************************************************************************
!         if (advtvd) then           ! TVD flux limiter scheme for updraft
!           l   = ib
!           lm1 = l - 1
!           dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
!    &                + alhl*(qoi(l)-qoi(lm1))
!           dtvd(1,2) = qoi(l) - qoi(lm1)
!           dtvd(1,3) = qli(l) - qli(lm1)
!           dtvd(1,4) = qii(l) - qii(lm1)
!           do l=ib+1,k
!             lm1 = l - 1
!             dtvd(2,1)   = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
!    &                    + alhl*(qoi(l)-qoi(lm1))
!             if (abs(dtvd(2,1)) > 1.0e-10) then
!               tem1        = dtvd(1,1) / dtvd(2,1)
!               tem2        = abs(tem1)
!               alfint(l,1) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for h
!             endif
!             dtvd(1,1)   = dtvd(2,1)
!
!             dtvd(2,2)   = qoi(l) - qoi(lm1)
!             if (abs(dtvd(2,2)) > 1.0e-10) then
!               tem1        = dtvd(1,2) / dtvd(2,2)
!               tem2        = abs(tem1)
!               alfint(l,2) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for q
!             endif
!             dtvd(1,2)   = dtvd(2,2)
!
!             dtvd(2,3)   = qli(l) - qli(lm1)
!             if (abs(dtvd(2,3)) > 1.0e-10) then
!               tem1        = dtvd(1,3) / dtvd(2,3)
!               tem2        = abs(tem1)
!               alfint(l,3) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for ql
!             endif
!             dtvd(1,3)   = dtvd(2,3)
!
!             dtvd(2,4)   = qii(l) - qii(lm1)
!             if (abs(dtvd(2,4)) > 1.0e-10) then
!               tem1        = dtvd(1,4) / dtvd(2,4)
!               tem2        = abs(tem1)
!               alfint(l,4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for qi
!             endif
!             dtvd(1,4)   = dtvd(2,4)
!           enddo
!
!           if (ntrc > 0) then
!             do n=1,ntrc
!               l = ib
!               dtvd(1,1)   = uvi(l,n) - uvi(l-1,n)
!               do l=ib+1,k
!                 dtvd(2,1)     = uvi(l,n) - uvi(l-1,n)
!                 if (abs(dtvd(2,1)) > 1.0e-10) then
!                   tem1        = dtvd(1,1) / dtvd(2,1)
!                   tem2          = abs(tem1)
!                   alfint(l,n+4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for tracers
!                 endif
!                 dtvd(1,1)     = dtvd(2,1)
!               enddo
!             enddo
!           endif
!         endif
!****************************************************************************
!
!     if (lprint) then
!       ia = ipt
!       print *,' toi=',(toi(ia,l),l=1,K)
!       print *,' qoi=',(qoi(ia,l),l=1,K),' kbl=',kbl
!       print *,' toi=',(toi(l),l=1,K)
!       print *,' qoi=',(qoi(l),l=1,K),' kbl=',kbl
!       print *,' prs=',(prs(l),l=1,K)
!     endif
!
          WFNC = 0.0
          do L=IB,K+1
            FLX(L) = 0.0
            FLXD(L)= 0.0
          enddo
!
!     if(lprint)then
!       print *,' CALLING CLOUD TYPE IB= ', IB,' DT=',DT,' K=',K
!    &,   'ipt=',ipt
!       print *,' TOI=',(TOI(L),L=IB,K)
!       print *,' QOI=',(QOI(L),L=IB,K)
!       print *,' qliin=',qli
!       print *,' qiiin=',qii
!     endif
!
          TLA = -10.0
!
          CALL CLOUD(lmhij, IB, ntrc, kblmx                             &
     &,              FRAC,  MAX_NEG_BOUY, vsmooth                       &
     &,              REVAP, WRKFUN, CALKBL, CRTFUN, DNDRFT, lprint      &
     &,              DT, KDT, TLA, DPD                                  &
     &,              ALFINT, rhfacl, rhfacs, garea(ipt)                 &
     &,              ccwf,   CDRAG(ipt), trcfac                         & 
     &,              alfind, rhc_l, phi_l, phi_h, PRS, PRSM,sgcs(1,ipt) &
     &,              TOI, QOI, UVI, QLI, QII, KBL, DDVEL(ipt)           &
     &,              TCU, QCU, RCU, PCU, FLX, FLXD, RAIN, WFNC, fscav_  &
     &               )
!    &,              ctei)

!     if (lprint) then
!       print *,' rain=',rain,' ipt=',ipt
!       print *,' after calling CLOUD TYPE IB= ', IB                    &
!    &,' rain=',rain,' prskd=',prs(ib),' qli=',qli(ib),' qii=',qii(ib)
!       print *,' phi_h=',phi_h(K-5:K+1)
!       print *,' TOI=',(TOI(L),L=1,K),' me=',me,' ib=',ib
!       print *,' QOI=',(QOI(L),L=1,K)
!       print *,' qliou=',qli
!       print *,' qiiou=',qii
!     endif
!
          do L=IB,K
            ll = l
            if (flipv) ll  = kp1 -l    ! Input variables are bottom to top!
            ud_mf(ipt,ll)  = ud_mf(ipt,ll)  + flx(l+1)
            dd_mf(ipt,ll)  = dd_mf(ipt,ll)  + flxd(l+1)
          enddo
          ll = ib
          if (flipv) ll  = kp1 - ib
          det_mf(ipt,ll) = det_mf(ipt,ll) + flx(ib)
! 
!
!   Warining!!!!
!   ------------
!   By doing the following, CLOUD does not contain environmental
!   condensate!
!
          if (.not. advcld) then
            do l=1,K
              clw(l ) = clw(l) + QLI(L)
              cli(l ) = cli(l) + QII(L)
              QLI(L)  = 0.0
              QII(L)  = 0.0
            enddo
          endif
!
        ENDDO                        ! End of the NC loop!
!
        RAINC(ipt) = rain * 0.001    ! Output rain is in meters

!     if (lprint) then
!       print*,' convective precip=',rain*86400/dt,' mm/day'
!    1,               ' ipt=',ipt
!        print *,' toi',(tn0(imax,l),l=1,k)
!        print *,' qoi',(qn0(imax,l),l=1,k)
!     endif
!
        do l=1,k
          ll = l
          if (flipv) ll  = kp1 - l
          tin(ipt,ll)    = toi(l)                   ! Temperature
          qin(ipt,ll)    = qoi(l)                   ! Specific humidity
          uin(ipt,ll)    = uvi(l,trac+1)            ! U momentum
          vin(ipt,ll)    = uvi(l,trac+2)            ! V momentum
          if (trac > 0) then
            do n=1,trac
              ccin(ipt,ll,n+2) = uvi(l,n)           ! Tracers
            enddo
          endif
        enddo
        if (advcld) then
          do l=1,k
            ll = l
            if (flipv) ll  = kp1 - l
            ccin(ipt,ll,1) = qii(l)          ! Cloud ice
            ccin(ipt,ll,2) = qli(l)          ! Cloud water
          enddo
        else
          do l=1,k
            ll = l
            if (flipv) ll  = kp1 - l
            ccin(ipt,ll,1) = ccin(ipt,ll,1) + cli(l)
            ccin(ipt,ll,2) = ccin(ipt,ll,2) + clw(l)
          enddo
        endif
!
        ktop(ipt) = kp1
        kbot(ipt) = 0

        kcnv(ipt) = 0

        do l=lmhij-1,1,-1
          if (sgcs(l,ipt) < 0.93 .and. tcu(l) .ne. 0.0) then
!         if (sgcs(l,ipt) < 0.90 .and. tcu(l) .ne. 0.0) then
!         if (sgcs(l,ipt) < 0.85 .and. tcu(l) .ne. 0.0) then
!         if (sgcs(l,ipt) < 0.70 .and. tcu(l) .ne. 0.0) then
!         if (sgcs(l,ipt) < 0.60 .and. tcu(l) .ne. 0.0) then
!         if (tcu(l) .ne. 0.0) then
             kcnv(ipt) = 1
          endif
!  New test for convective clouds ! added in 08/21/96
          if (clw(l)+cli(l) .gt. 0.0 .OR.                               &
     &        qli(l)+qii(l) .gt. clwmin) ktop(ipt) = l
        enddo
        do l=1,km1
          if (clw(l)+cli(l) .gt. 0.0 .OR.                               &
     &        qli(l)+qii(l) .gt. clwmin) kbot(ipt) = l
        enddo
        if (flipv) then
          ktop(ipt) = kp1 - ktop(ipt)
          kbot(ipt) = kp1 - kbot(ipt)
        endif
!
!     if (lprint) then
!        print *,' tin',(tin(ia,l),l=k,1,-1)
!        print *,' qin',(qin(ia,l),l=k,1,-1)
!     endif
!
!     Velocity scale from the downdraft!
!
        DDVEL(ipt) = DDVEL(ipt) * DDFAC * GRAV / (prs(K+1)-prs(k))
!
      ENDDO                            ! End of the IPT Loop!

      deallocate (alfint,uvi,trcfac,rcu)
!
      RETURN
      END
      SUBROUTINE CRTWRK(PL, CCWF, ACR)
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : ac, ad
      Implicit none
!
      real(kind=kind_phys) PL, CCWF, ACR
      INTEGER IWK
!
      IWK = PL * 0.02 - 0.999999999
      IWK = MAX(1, MIN(IWK,16))
      ACR = (AC(IWK) + PL * AD(IWK)) * CCWF
!
      RETURN
      END
      SUBROUTINE CLOUD(                                                 &
     &                  K, KD, NTRC, KBLMX                              &
     &,                 FRACBL, MAX_NEG_BOUY, vsmooth                   &
     &,                 REVAP, WRKFUN, CALKBL, CRTFUN, DNDRFT, lprnt    &
     &,                 DT, KDT, TLA, DPD                               &
     &,                 ALFINT, RHFACL, RHFACS, garea, ccwf, cd, trcfac &
     &,                 alfind, rhc_ls, phil, phih, prs, prsm, sgcs     &
     &,                 TOI, QOI, ROI,  QLI, QII, KPBL, DSFC            &
     &,                 TCU, QCU, RCU, PCU, FLX, FLXD, CUP, WFNC,fscav_ &
     &                  )
!    &,                 ctei)

!
!***********************************************************************
!******************** Relaxed  Arakawa-Schubert ************************
!****************** Plug Compatible Scalar Version *********************
!************************ SUBROUTINE CLOUD  ****************************
!************************  October 2004     ****************************
!********************  VERSION 2.0  (modified) *************************
!************* Shrinivas.Moorthi@noaa.gov (301) 763 8000(X7233) ********
!***********************************************************************
!*Reference:
!-----------
!     NOAA Technical Report NWS/NCEP 99-01:
!     Documentation of Version 2 of Relaxed-Arakawa-Schubert
!     Cumulus Parameterization with Convective Downdrafts, June 1999.
!     by S. Moorthi and M. J. Suarez.
!
!***********************************************************************
!
!===>    UPDATES CLOUD TENDENCIES DUE TO A SINGLE CLOUD
!===>    DETRAINING AT LEVEL KD.
!
!***********************************************************************
!
!===>  TOI(K)     INOUT   TEMPERATURE             KELVIN
!===>  QOI(K)     INOUT   SPECIFIC HUMIDITY       NON-DIMENSIONAL
!===>  ROI(K,NTRC)INOUT   TRACER                  ARBITRARY
!===>  QLI(K)     INOUT   LIQUID WATER            NON-DIMENSIONAL
!===>  QII(K)     INOUT   ICE                     NON-DIMENSIONAL

!===>  PRS(K+1)   INPUT   PRESSURE @ EDGES        MB
!===>  PRSM(K)    INPUT   PRESSURE @ LAYERS       MB
!===>  SGCS(K)    INPUT   Local sigma
!===>  PHIH(K+1)  INPUT   GEOPOTENTIAL @ EDGES  IN MKS units
!===>  PHIL(K)    INPUT   GEOPOTENTIAL @ LAYERS IN MKS units
!===>  PRJ(K+1)   INPUT   (P/P0)^KAPPA  @ EDGES   NON-DIMENSIONAL
!===>  PRJM(K)    INPUT   (P/P0)^KAPPA  @ LAYERS  NON-DIMENSIONAL

!===>  K          INPUT   THE RISE & THE INDEX OF THE SUBCLOUD LAYER
!===>  KD         INPUT   DETRAINMENT LEVEL ( 1<= KD < K )          
!===>  NTRC       INPUT   NUMBER OF TRACERS. MAY BE ZERO.
!===>  kblmx      INPUT   highest level the pbl can take
!===>  DNDRFT     INPUT   LOGICAL .TRUE. OR .FALSE.
!===>  DPD        INPUT   Minumum Cloud Depth for DOWNDRFAT Computation hPa
!
!===>  TCU(K  )   UPDATE  TEMPERATURE TENDENCY       DEG
!===>  QCU(K  )   UPDATE  WATER VAPOR TENDENCY       (G/G)
!===>  RCU(K,NTRC)UPDATE  TRACER TENDENCIES          ND
!===>  PCU(K-1)   UPDATE  PRECIP @ BASE OF LAYER     KG/M^2
!===>  FLX(K  )   UPDATE  MASS FLUX @ TOP OF LAYER   KG/M^2
!===>  CUP        UPDATE  PRECIPITATION AT THE SURFACE KG/M^2
!
      USE MACHINE , ONLY : kind_phys
      use module_ras
      IMPLICIT NONE
!
!  INPUT ARGUMENTS

!     LOGICAL REVAP, DNDRFT, WRKFUN, CALKBL, CRTFUN, CALCUP, ctei
      LOGICAL REVAP, DNDRFT, WRKFUN, CALKBL, CRTFUN, CALCUP
      logical vsmooth, lprnt
      INTEGER K, KD, NTRC, kblmx


      real(kind=kind_phys) TOI(K),    QOI(K ),  PRS(K+1),  PRSM(K)      &
     &,                    QLI(K),    QII(K),   PHIH(K+1), PHIL(K)      &
     &,                    ROI(K,NTRC), SGCS(K)
      real(kind=kind_phys) CD,        UFN,     DSFC
      INTEGER KPBL,   KBL,     KB1, kdt

      real(kind=kind_phys) FRACBL, MAX_NEG_BOUY,                        &
     &                     ALFINT(K,NTRC+4), RHFACL, RHFACS, garea, ccwf
      real(kind=kind_phys) DPD, alfind(k), rhc_ls(k)
      real(kind=kind_phys) trcfac(k,NTRC)
 
!  UPDATE ARGUMENTS

      real(kind=kind_phys) TCU(K),   QCU(K),    RCU(K,NTRC)             &
     &,                    TCD(K),   QCD(K),    PCU(K)                  &
     &,                    FLX(K+1), FLXD(K+1), CUP

!  TEMPORARY WORK SPACE

      real(kind=kind_phys) HOL(KD:K),  QOL(KD:K),   GAF(KD:K+1)         &
     &,                    HST(KD:K),  QST(KD:K),   TOL(KD:K)           &
     &,                    GMH(KD:K),  GMS(KD:K+1), GAM(KD:K+1)         &
     &,                    AKT(KD:K),  AKC(KD:K),   BKC(KD:K)           &
     &,                    LTL(KD:K),  RNN(KD:K),   FCO(KD:K)           &
     &,                                PRI(KD:K)                        &
     &,                    QIL(KD:K),  QLL(KD:K)                        &
     &,                    ZET(KD:K),  XI(KD:K),    RNS(KD:K)           &
     &,                    Q0U(KD:K),  Q0D(KD:K),   vtf(KD:K)           &
     &,                    DLB(KD:K+1),DLT(KD:K+1), ETA(KD:K+1)         &
     &,                    PRL(KD:K+1)                                  &
     &,                    CIL(KD:K),  CLL(KD:K),   ETAI(KD:K)          &
     &,                    dlq(kd:k)

      real(kind=kind_phys) ALM,   DET,    HCC,  CLP                     &
     &,                    HSU,   HSD,    QTL,  QTV                     &
     &,                    AKM,   WFN,    HOS,  QOS                     &
     &,                    AMB,   TX1,    TX2,  TX3                     &
     &,                    TX4,   TX5,    QIS,  QLS                     &
     &,                    HBL,   QBL,    RBL(NTRC)                     &
     &,                    QLB,   QIB,    PRIS                          &
     &,                    WFNC,  TX6,    ACR                           &
     &,                    TX7,   TX8,    TX9,  RHC                     &
     &,                    hstkd, qstkd,  ltlkd, q0ukd, q0dkd, dlbkd    &
     &,                    qtp, qw00, qi00, qrbkd                       &
     &,                    hstold, rel_fac, prism

      real(kind=kind_phys) fscav_(ntrc)

      real(kind=kind_phys) wrk1(kd:k), wrk2(kd:k)


      LOGICAL ep_wfn, cnvflg

      LOGICAL LOWEST, SKPDD

      real(kind=kind_phys) TL, PL, QL, QS, DQS, ST1, SGN, TAU,          &
     &                     QTVP, HB, QB, TB, QQQ,                       &
     &                     HCCP, DS, DH, AMBMAX, X00, EPP, QTLP,        &
     &                     DPI, DPHIB, DPHIT, DEL_ETA, DETP,            &
     &                     TEM, TEM1, TEM2, TEM3, TEM4,                 &
     &                     ST2, ST3, ST4, ST5,                          &
     &                     ERRMIN, ERRMI2, ERRH, ERRW, ERRE, TEM5,      &
     &                     TEM6, HBD, QBD, st1s, shal_fac, hmax, hmin,  &
     &                     dhdpmn, dhdp(kd:k)
      parameter (ERRMIN=0.0001, ERRMI2=0.1*ERRMIN)
      INTEGER I, L,  N,  KD1, II                                        &
     &,       KP1, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1, kbls, kmxh
     &,       kblh, kblm, kblpmn, kmax, kmaxm1, kmaxp1, klcl, kmin, kmxb

      real avt, avq, avr, avh
!

!     real(kind=kind_phys), parameter :: rainmin=1.0e-9
      real(kind=kind_phys), parameter :: rainmin=1.0e-8
      real(kind=kind_phys), parameter :: oneopt9=1.0/0.09
      real(kind=kind_phys), parameter :: oneopt4=1.0/0.04

      real(kind=kind_phys) CLFRAC, DT, clf, clvfr

      real(kind=kind_phys) ACTEVAP,AREARAT,DELTAQ,MASS,MASSINV,POTEVAP  &
     &,                    TEQ,QSTEQ,DQDT,QEQ
!
!     Temporary workspace and parameters needed for downdraft
!
      real(kind=kind_phys) TLA, GMF
!
      real(kind=kind_phys) BUY(KD:K+1), QRB(KD:K),   QRT(KD:K)          &
     &,                    ETD(KD:K+1), HOD(KD:K+1), QOD(KD:K+1)        &
     &,                    GHD(KD:K),   GSD(KD:K),   EVP(KD:K)          &
     &,                    ETZ(KD:K),   CLDFR(KD:K), ETZI(KD:K-1)       &
     &,                    TRAIN, DOF, CLDFRD                           &
     &,                    FAC, RSUM1, RSUM2, RSUM3, dpneg, hcrit
      INTEGER IDH, lcon
      LOGICAL DDFT, UPDRET
!
!  Scavenging related parameters
!
      real                delzkm          ! layer thickness in km
      real                fnoscav         ! fraction of tracer *not* scavenged
!
!***********************************************************************
!
      do l=1,K
        tcd(L) = 0.0
        qcd(L) = 0.0
      enddo
!
      KP1     = K  + 1
      KM1     = K  - 1
      KD1     = KD + 1
!
!     if (lprnt) then
!       print *,' IN CLOUD for KD=',kd
!       print *,' prs=',prs(Kd:K+1)
!       print *,' phil=',phil(KD:K)
!       print *,' phih=',phih(1:K+1),' kdt=',kdt
!       print *,' phih=',phih(KD:K+1)
!       print *,' toi=',toi
!       print *,' qoi=',qoi
!     endif
!
      CLDFRD   = 0.0
      DOF      = 0.0
      PRL(KP1) = PRS(KP1)
!
      DO L=KD,K
        RNN(L) = 0.0
        ZET(L) = 0.0
        XI(L)  = 0.0
!
        TOL(L) = TOI(L)
        QOL(L) = QOI(L)
        PRL(L) = PRS(L)
        BUY(L) = 0.0
        CLL(L) = QLI(L)
        CIL(L) = QII(L)
      ENDDO
!
      if (vsmooth) then
        do l=kd,k
          wrk1(l) = tol(l)
          wrk2(l) = qol(l)
        enddo
        do l=kd1,km1
          tol(l) = 0.25*wrk1(l-1) + 0.5*wrk1(l) + 0.25*wrk1(l+1)
          qol(l) = 0.25*wrk2(l-1) + 0.5*wrk2(l) + 0.25*wrk2(l+1)
        enddo
      endif
!
      DO L=KD, K
        DPI    = ONE / (PRL(L+1) - PRL(L))
        PRI(L) = GRAVFAC * DPI
!
        PL     = PRSM(L)
        TL     = TOL(L)

        AKT(L) = (PRL(L+1) - PL) * DPI
!
        CALL QSATCN(TL, PL, QS, DQS)
!       CALL QSATCN(TL, PL, QS, DQS,lprnt)
!
        QST(L) = QS
        GAM(L) = DQS * ELOCP
        ST1    = ONE + GAM(L)
        GAF(L) = (ONE/ALHL) * GAM(L)/ST1
 
        QL     = MAX(MIN(QS*RHMAX,QOL(L)), ONE_M10)
        QOL(L) = QL
 
        TEM    = CP * TL
        LTL(L) = TEM * ST1 / (ONE+NU*(QST(L)+TL*DQS))
        vtf(L) = 1.0 + NU * QL
        ETA(L) = ONE / (LTL(L) * VTF(L))

        HOL(L) = TEM + QL * ALHL
        HST(L) = TEM + QS * ALHL
!
      ENDDO
!
      ETA(K+1) = ZERO
      GMS(K)   = ZERO
!
      AKT(KD)  = HALF
      GMS(KD)  = ZERO
!
      CLP      = ZERO
!
      GAM(K+1) = GAM(K)
      GAF(K+1) = GAF(K)
!
      DO L=K,KD1,-1
        DPHIB  = PHIL(L) - PHIH(L+1)
        DPHIT  = PHIH(L) - PHIL(L)
!
        DLB(L) = DPHIB * ETA(L)
        DLT(L) = DPHIT * ETA(L)
!
        QRB(L) = DPHIB
        QRT(L) = DPHIT
!
        ETA(L) = ETA(L+1) + DPHIB

        HOL(L) = HOL(L) + ETA(L)
        hstold = hst(l)
        HST(L) = HST(L) + ETA(L)
!
        ETA(L) = ETA(L) + DPHIT
      ENDDO
!
!     For the cloud top layer
!
      L = KD

      DPHIB  = PHIL(L) - PHIH(L+1)
!
      DLB(L) = DPHIB * ETA(L)
!
      QRB(L) = DPHIB
      QRT(L) = DPHIB
!
      ETA(L) = ETA(L+1) + DPHIB

      HOL(L) = HOL(L) + ETA(L)
      HST(L) = HST(L) + ETA(L)
!
!     if (kd .eq. 12) then
!       if (lprnt) then
!         print *,' IN CLOUD for KD=',KD,' K=',K
!         print *,' l=',l,' hol=',hol(l),' hst=',hst(l)
!         print *,' TOL=',tol
!         print *,' qol=',qol
!         print *,' hol=',hol
!         print *,' hst=',hst
!       endif
!     endif
!
!     To determine KBL internally -- If KBL is defined externally
!     the following two loop should be skipped
!
!     if (lprnt) print *,' calkbl=',calkbl

      hcrit = hcritd
      if (sgcs(kd) > 0.65) hcrit = hcrits
      IF (CALKBL) THEN
         KTEM = MAX(KD+1, KBLMX)
         hmin = hol(k)
         kmin = k
         do l=km1,kd,-1
           if (hmin > hol(l)) then
             hmin = hol(l)
             kmin = l
           endif
         enddo
         if (kmin == k) return
         hmax = hol(k)
         kmax = k
         do l=km1,ktem,-1
           if (hmax < hol(l)) then
             hmax = hol(l)
             kmax = l
           endif
         enddo
         kmxb = kmax
         if (kmax < kmin) then
           kmax = k
           kmxb = k
           hmax = hol(kmax)
         elseif (kmax < k) then
           do l=kmax+1,k
             if (abs(hol(kmax)-hol(l)) > 0.5 * hcrit) then
               kmxb = l - 1
               exit
             endif
           enddo
         endif
         kmaxm1 = kmax - 1
         kmaxp1 = kmax + 1
         kblpmn = kmax
!
         dhdp(kmax:k) = 0.0
         dhdpmn = dhdp(kmax)
         do l=kmaxm1,ktem,-1
           dhdp(l) = (HOL(L)-HOL(L+1)) / (PRL(L+2)-PRL(L))
           if (dhdp(l) < dhdpmn) then
             dhdpmn = dhdp(l)
             kblpmn = l + 1
           elseif (dhdp(l) > 0.0 .and. l <= kmin) then
             exit
           endif
         enddo
         kbl = kmax
         if (kblpmn < kmax) then
           do l=kblpmn,kmaxm1
             if (hmax-hol(l) < 0.5*hcrit) then
               kbl = l
               exit
             endif
           enddo
         endif
       
!     if(lprnt) print *,' kbl=',kbl,' kbls=',kbls,' kmax=',kmax
!
         klcl = kd1
         if (kmax > kd1) then
           do l=kmaxm1,kd1,-1
             if (hmax > hst(l)) then
               klcl = l+1
               exit
             endif
           enddo
         endif
!        if(lprnt) print *,' klcl=',klcl,' ii=',ii
!        if (klcl == kd .or. klcl < ktem) return

!        This is to handle mid-level convection from quasi-uniform h

         if (kmax < kmxb) then
           kmax   = max(kd1, min(kmxb,k))
           kmaxm1 = kmax - 1
           kmaxp1 = kmax + 1
         endif


!        if (prl(Kmaxp1) - prl(klcl) > 250.0 ) return

         ii  = max(kbl,kd1)
         kbl = max(klcl,kd1)
         tem = min(50.0,max(10.0,(prl(kmaxp1)-prl(kd))*0.10))
         if (prl(kmaxp1) - prl(ii) > tem .and. ii > kbl) kbl = ii

!        if(lprnt) print *,' kbl2=',kbl,' ii=',ii

         if (kbl .ne. ii) then
           if (PRL(kmaxp1)-PRL(KBL) > bldmax) kbl = max(kbl,ii)
         endif
         if (kbl < ii) then
           if (hol(ii)-hol(ii-1) > 0.5*hcrit) kbl = ii
         endif

!        if (prl(kbl) - prl(klcl) > 300.0 ) return
         if (prl(kbl) - prl(klcl) > 250.0 ) return
!
         KBL  = min(kmax, MAX(KBL,KBLMX))
!        kbl  = min(kblh,kbl)
!!!
!        tem1 = max(prl(k+1)-prl(k),                                    &
!    &                     min((prl(kbl) - prl(kd))*0.05, 10.0))
!!   &                     min((prl(kbl) - prl(kd))*0.05, 20.0))
!!   &                     min((prl(kbl) - prl(kd))*0.05, 30.0))
!        if (prl(k+1)-prl(kbl) .lt. tem1) then
!          KTEM = MAX(KD+1, KBLMX)
!          do l=k,KTEM,-1
!            tem = prl(k+1) - prl(l)
!            if (tem .gt. tem1) then
!              kbl = min(kbl,l)
!              exit
!            endif
!          enddo
!        endif
!        if (kbl == kblmx .and. kmax >= k-1) kbl = k - 1
!!!

         KPBL = KBL

!     if(lprnt)print*,' 1st kbl=',kbl,' kblmx=',kblmx,' kd=',kd
!     if(lprnt)print*,' tx3=',tx3,' tx1=',tx1,' tem=',tem
!    1,               ' hcrit=',hcrit

      ELSE
         KBL  = KPBL
!     if(lprnt)print*,' 2nd kbl=',kbl
      ENDIF

!     if(lprnt)print*,' after CALKBL l=',l,' hol=',hol(l)
!    1,               ' hst=',hst(l)
!
      KBL      = min(kmax,MAX(KBL,KD+2))
      KB1      = KBL - 1
!!
!     if (lprnt) print *,' kbl=',kbl,' prlkbl=',prl(kbl),prl(k+1)

      if(PRL(Kmaxp1)-PRL(KBL) .gt. bldmax .or. kb1 .le. kd) then
        return
      endif
!
!     if (lprnt) print *,' kbl=',kbl
!     write(0,*)' kbl=',kbl,' kmax=',kmax,' kmaxp1=',kmaxp1,' k=',k
!
      PRIS     = ONE / (PRL(K+1)-PRL(KBL))
      PRISM    = ONE / (PRL(Kmaxp1)-PRL(KBL))
      TX1      = ETA(KBL)
!
      GMS(KBL) = 0.0
      XI(KBL)  = 0.0
      ZET(KBL) = 0.0
!
      shal_fac = 1.0
!     if (prl(kbl)-prl(kd) < 300.0 .and. kmax == k) shal_fac = shalfac
      if (prl(kbl)-prl(kd) < 350.0 .and. kmax == k) shal_fac = shalfac
      DO L=Kmax,KD,-1
        IF (L >= KBL) THEN
          ETA(L) = (PRL(Kmaxp1)-PRL(L)) * PRISM
        ELSE
          ZET(L) = (ETA(L) - TX1) * ONEBG
          XI(L)  =  ZET(L) * ZET(L) * (QUDFAC*shal_fac)
          ETA(L) =  ZET(L) - ZET(L+1)
          GMS(L) =  XI(L)  - XI(L+1)
        ENDIF
!       if (lprnt) print *,' l=',l,' eta=',eta(l),' kbl=',kbl
      ENDDO
      if (kmax < k) then
        do l=kmaxp1,kp1
          eta(l) = 0.0
        enddo
      endif
!
      HBL = HOL(Kmax) * ETA(Kmax)
      QBL = QOL(Kmax) * ETA(Kmax)
      QLB = CLL(Kmax) * ETA(Kmax)
      QIB = CIL(Kmax) * ETA(Kmax)
      TX1 = QST(Kmax) * ETA(Kmax)
!
      DO L=Kmaxm1,KBL,-1
         TEM = ETA(L) - ETA(L+1)
         HBL = HBL + HOL(L) * TEM
         QBL = QBL + QOL(L) * TEM
         QLB = QLB + CLL(L) * TEM
         QIB = QIB + CIL(L) * TEM
         TX1 = TX1 + QST(L) * TEM
      ENDDO

!     if (ctei .and. sgcs(kd) > 0.65) then
!        hbl = hbl * hpert_fac
!        qbl = qbl * hpert_fac
!     endif

!     if (lprnt) print *,' hbl=',hbl,' qbl=',qbl
!                                   Find Min value of HOL in TX2
      TX2 = HOL(KD)
      IDH = KD1
      DO L=KD1,KB1
        IF (HOL(L) .LT. TX2) THEN
           TX2 = HOL(L)
           IDH = L             ! Level of minimum moist static energy!
        ENDIF
      ENDDO
      IDH = 1
      IDH = MAX(KD1, IDH)
!
      TEM1 = HBL - HOL(KD)
      TEM  = HBL - HST(KD1) - LTL(KD1) * NU *(QOL(KD1)-QST(KD1))
      LOWEST = KD .EQ. KB1

      lcon = kd
      do l=kb1,kd1,-1
        if (hbl >= hst(l)) then
          lcon = l
          exit
        endif
      enddo
!
      if (lcon == kd .or. kbl <= kd .or. prl(kbl)-prsm(lcon) > 150.0)   &
     &                                    return
!
      TX1    = RHFACS - QBL / TX1       !     Average RH

      cnvflg = (TEM .GT. ZERO .OR. (LOWEST .AND. TEM1 .GE. ZERO))       &
     &         .AND. (TX1 .LT. RHRAM)

!     if(lprnt) print *,' cnvflg=',cnvflg,' tem=',tem,' tem1=',tem1
!    &,' tx1=',tx1,' rhram=',rhram,' kbl=',kbl,' kd=',kd,' lowest='
!    &,lowest,' rhfacs=',rhfacs,' ltl=',ltl(kd1),' qol=',qol(kd1)
!    &,' qst=',qst(kd1),' hst=',hst(kd1),' nu=',nu
!     if(lprnt .and. (.not. cnvflg)) print *,' tx1=',tx1,' rhfacs='
!    &,rhfacs, ' tem=',tem,' hst=',hst(kd1)

      IF (.NOT. cnvflg) RETURN
!
      RHC    = MAX(ZERO, MIN(ONE, EXP(-20.0*TX1) ))
!
      if (ntrc > 0) then
        DO N=1,NTRC
          RBL(N) = ROI(Kmax,N) * ETA(Kmax)
        ENDDO
        DO N=1,NTRC
          DO L=KmaxM1,KBL,-1
            RBL(N) = RBL(N) + ROI(L,N)*(ETA(L)-ETA(L+1))
          ENDDO
        ENDDO
      endif
!
      TX4    = 0.0
      TX5    = 0.0
!
      TX3      = QST(KBL) - GAF(KBL) * HST(KBL)
      DO L=KBL,K
        QIL(L) = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
      ENDDO
!
      DO L=KB1,KD1,-1
        TEM      = QST(L) - GAF(L) * HST(L)
        TEM1     = (TX3 + TEM) * 0.5
        ST2      = (GAF(L)+GAF(L+1)) * 0.5
!
        FCO(L+1) =            TEM1 + ST2 * HBL

!     if(lprnt) print *,' fco=',fco(l+1),' tem1=',tem1,' st2=',st2
!    &,' hbl=',hbl,' tx3=',tx3,' tem=',tem,' gaf=',gaf(l),' l=',l

        RNN(L+1) = ZET(L+1) * TEM1 + ST2 * TX4
        GMH(L+1) = XI(L+1)  * TEM1 + ST2 * TX5
!
        TX3      = TEM
        TX4      = TX4 + ETA(L) * HOL(L)
        TX5      = TX5 + GMS(L) * HOL(L)
!
        QIL(L)   = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
        QLL(L+1) = (0.5*ALHF) * ST2 * (QIL(L)+QIL(L+1)) + ONE
      ENDDO
!
!     FOR THE CLOUD TOP -- L=KD
!
      L = KD
!
      TEM      = QST(L) - GAF(L) * HST(L)
      TEM1     = (TX3 + TEM) * 0.5
      ST2      = (GAF(L)+GAF(L+1)) * 0.5
!
      FCO(L+1) =            TEM1 + ST2 * HBL
      RNN(L+1) = ZET(L+1) * TEM1 + ST2 * TX4
      GMH(L+1) = XI(L+1)  * TEM1 + ST2 * TX5
!
      FCO(L)   = TEM + GAF(L) * HBL
      RNN(L)   = TEM * ZET(L) + (TX4 + ETA(L)*HOL(L)) * GAF(L)
      GMH(L)   = TEM * XI(L)  + (TX5 + GMS(L)*HOL(L)) * GAF(L)
!
!   Replace FCO for the Bottom
!
      FCO(KBL) = QBL
      RNN(KBL) = 0.0
      GMH(KBL) = 0.0
!
      QIL(KD)  =  MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(KD))*TCRF))
      QLL(KD1) = (0.5*ALHF) * ST2 * (QIL(KD) + QIL(KD1)) + ONE
      QLL(KD ) = ALHF * GAF(KD) * QIL(KD) + ONE
!
!     if (lprnt) then
!       print *,' fco=',fco(kd:kbl)
!       print *,' qil=',qil(kd:kbl)
!       print *,' qll=',qll(kd:kbl)
!     endif
!
      st1  = qil(kd)
      st2  = c0i * st1
      tem  = c0  * (1.0-st1)
      tem2 = st2*qi0 + tem*qw0
!
      DO L=KD,KB1
         tx2    = akt(l) * eta(l)
         tx1    = tx2 * tem2
         q0u(l) = tx1
         FCO(L) = FCO(L+1) - FCO(L) + tx1
         RNN(L) = RNN(L+1) - RNN(L)                                     &
     &          + ETA(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*zet(l)
         GMH(L) = GMH(L+1) - GMH(L)                                     &
     &          + GMS(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*xi(l)
!
         tem1   = (1.0-akt(l)) * eta(l)

!     if(lprnt) print *,' qll=',qll(l),' st2=',st2,' tem=',tem
!    &,' tx2=',tx2,' akt=',akt(l),' eta=',eta(l)

         AKT(L) = QLL(L)   + (st2 + tem) * tx2

!     if(lprnt) print *,' akt==',akt(l),' l==',l

         AKC(L) = 1.0 / AKT(L)
!
         st1    = 0.5 * (qil(l)+qil(l+1))
         st2    = c0i * st1
         tem    = c0  * (1.0-st1)
         tem2   = st2*qi0 + tem*qw0
!
         BKC(L) = QLL(L+1) - (st2 + tem) * tem1
!
         tx1    = tem1*tem2
         q0d(l) = tx1
         FCO(L) = FCO(L) + tx1
         RNN(L) = RNN(L) + tx1*zet(l+1)
         GMH(L) = GMH(L) + tx1*xi(l+1)
      ENDDO

!     if(lprnt) print *,' akt=',akt(kd:kb1)
!     if(lprnt) print *,' akc=',akc(kd:kb1)

      qw00 = qw0
      qi00 = qi0
      ii = 0
  777 continue
!
!     if (lprnt) print *,' after 777 ii=',ii,' ep_wfn=',ep_wfn
!
      ep_wfn = .false.
      RNN(KBL) = 0.0
      TX3      = bkc(kb1) * (QIB + QLB)
      TX4      = 0.0
      TX5      = 0.0
      DO L=KB1,KD1,-1
        TEM    = BKC(L-1)       * AKC(L)
!     if (lprnt) print *,' tx3=',tx3,' fco=',fco(l),' akc=',akc(l)
!    &,' bkc=',bkc(l-1), ' l=',l
        TX3    = (TX3 + FCO(L)) * TEM
        TX4    = (TX4 + RNN(L)) * TEM
        TX5    = (TX5 + GMH(L)) * TEM
      ENDDO
      IF (KD .LT. KB1) THEN
         HSD   = HST(KD1) + LTL(KD1) *  NU *(QOL(KD1)-QST(KD1))
      ELSE
         HSD   = HBL
      ENDIF
!
!     if (lprnt) print *,' tx3=',tx3,' fco=',fco(kd),' akc=',akc(kd)

      TX3 = (TX3 + FCO(KD)) * AKC(KD)
      TX4 = (TX4 + RNN(KD)) * AKC(KD)
      TX5 = (TX5 + GMH(KD)) * AKC(KD)
      ALM = ALHF*QIL(KD) - LTL(KD) * VTF(KD)
!
      HSU = HST(KD) + LTL(KD) * NU * (QOL(KD)-QST(KD))

!     if (lprnt) print *,' hsu=',hsu,' hst=',hst(kd),
!    &' ltl=',ltl(kd),' qol=',qol(kd),' qst=',qst(kd)
!
!===> VERTICAL INTEGRALS NEEDED TO COMPUTE THE ENTRAINMENT PARAMETER
!
      TX1 = ALM * TX4
      TX2 = ALM * TX5

      DO L=KD,KB1
        TAU = HOL(L) - HSU
        TX1 = TX1 + TAU * ETA(L)
        TX2 = TX2 + TAU * GMS(L)
      ENDDO
!
!     MODIFY HSU TO INCLUDE CLOUD LIQUID WATER AND ICE TERMS
!
!     if (lprnt) print *,' hsu=',hsu,' alm=',alm,' tx3=',tx3

      HSU    = HSU - ALM * TX3
!
      CLP    = ZERO
      ALM    = -100.0
      HOS    = HOL(KD)
      QOS    = QOL(KD)
      QIS    = CIL(KD)
      QLS    = CLL(KD)
      cnvflg = HBL .GT. HSU .and. abs(tx1) .gt. 1.0e-4

!     if (lprnt) print *,' ii=',ii,' cnvflg=',cnvflg,' hsu=',hsu
!    &,' hbl=',hbl,' tx1=',tx1,' hsd=',hsd


!***********************************************************************


       ST1  = HALF*(HSU + HSD)
       IF (cnvflg) THEN
!
!  STANDARD CASE:
!   CLOUD CAN BE NEUTRALLY BOUYANT AT MIDDLE OF LEVEL KD W/ +VE LAMBDA.
!   EPP < .25 IS REQUIRED TO HAVE REAL ROOTS.
!
       clp = 1.0
       st2 = hbl - hsu

!     if(lprnt) print *,' tx2=',tx2,' tx1=',tx1,' st2=',st2
!
       if (tx2 .eq. 0.0) then
         alm = - st2 / tx1
         if (alm .gt. almax) alm = -100.0
       else
         x00 = tx2 + tx2
         epp = tx1 * tx1 - (x00+x00)*st2
         if (epp .gt. 0.0) then
           x00  = 1.0 / x00
           tem  = sqrt(epp)
           tem1 = (-tx1-tem)*x00
           tem2 = (-tx1+tem)*x00
           if (tem1 .gt. almax) tem1 = -100.0
           if (tem2 .gt. almax) tem2 = -100.0
           alm  = max(tem1,tem2)

!     if (lprnt) print *,' tem1=',tem1,' tem2=',tem2,' alm=',alm
!    &,' tx1=',tx1,' tem=',tem,' epp=',epp,' x00=',x00,' st2=',st2

         endif
       endif

!     if (lprnt) print *,' almF=',alm,' ii=',ii,' qw00=',qw00
!    &,' qi00=',qi00
!
!  CLIP CASE:
!   NON-ENTRAINIG CLOUD DETRAINS IN LOWER HALF OF TOP LAYER.
!   NO CLOUDS ARE ALLOWED TO DETRAIN BELOW THE TOP LAYER.
!
       ELSEIF ( (HBL .LE. HSU) .AND.                                    &
     &          (HBL .GT. ST1   )     ) THEN
         ALM = ZERO
!        CLP = (HBL-ST1) / (HSU-ST1)    ! commented on Jan 16, 2010
       ENDIF
!
      cnvflg = .TRUE.
      IF (ALMIN1 .GT. 0.0) THEN
        IF (ALM .GE. ALMIN1) cnvflg = .FALSE.
      ELSE
        LOWEST   = KD .EQ. KB1
        IF ( (ALM .GT. ZERO) .OR.                                       &
     &      (.NOT. LOWEST .AND. ALM .EQ. ZERO) ) cnvflg = .FALSE.
      ENDIF
!
!===>  IF NO SOUNDING MEETS SECOND CONDITION, RETURN
!
      IF (cnvflg) THEN
         IF (ii .gt. 0 .or. (qw00 .eq. 0.0 .and. qi00 .eq. 0.0)) RETURN
         CLP = 1.0
         ep_wfn = .true.
         GO TO 888
      ENDIF
!
!     if (lprnt) print *,' hstkd=',hst(kd),' qstkd=',qst(kd)
!    &,' ii=',ii,' clp=',clp

      st1s = ONE
      IF(CLP.GT.ZERO .AND. CLP.LT.ONE) THEN
        ST1     = HALF*(ONE+CLP)
        ST2     = ONE - ST1
        st1s    = st1
        hstkd   = hst(kd)
        qstkd   = qst(kd)
        ltlkd   = ltl(kd)
        q0ukd   = q0u(kd)
        q0dkd   = q0d(kd)
        dlbkd   = dlb(kd)
        qrbkd   = qrb(kd)
!
        HST(KD) = HST(KD)*ST1 + HST(KD1)*ST2
        HOS     = HOL(KD)*ST1 + HOL(KD1)*ST2
        QST(KD) = QST(KD)*ST1 + QST(KD1)*ST2
        QOS     = QOL(KD)*ST1 + QOL(KD1)*ST2
        QLS     = CLL(KD)*ST1 + CLL(KD1)*ST2
        QIS     = CIL(KD)*ST1 + CIL(KD1)*ST2
        LTL(KD) = LTL(KD)*ST1 + LTL(KD1)*ST2
!
        DLB(KD) = DLB(KD)*CLP
        qrb(KD) = qrb(KD)*CLP
        ETA(KD) = ETA(KD)*CLP
        GMS(KD) = GMS(KD)*CLP
        Q0U(KD) = Q0U(KD)*CLP
        Q0D(KD) = Q0D(KD)*CLP
      ENDIF
!
!
!***********************************************************************
!
!    Critical workfunction is included in this version
!
      ACR = 0.0
      TEM = PRL(KD1) - (PRL(KD1)-PRL(KD)) * CLP * HALF
      tx1 = PRL(KBL) - TEM
      tx2 = min(900.0,max(tx1,100.0))
      tem1    = log(tx2*0.01) / log(10.0)
      if ( kdt == 1 ) then
        rel_fac = (dt * facdt)  / (tem1*12.0 + (1-tem1)*3.0)
      else
        rel_fac = (dt * facdt) / (tem1*adjts_d + (1-tem1)*adjts_s)
      endif
!
!     rel_fac = max(zero, min(one,rel_fac))
      rel_fac = max(zero, min(half,rel_fac))
      
      IF (CRTFUN) THEN
        CALL CRTWRK(TEM, CCWF, ST1)
        ACR = TX1 * ST1
      ENDIF
!
!===>  NORMALIZED MASSFLUX
!
!  ETA IS THE THICKNESS COMING IN AND THE MASS FLUX GOING OUT.
!  GMS IS THE THICKNESS OF THE SQUARE; IT IS LATER REUSED FOR GAMMA_S
!
!     ETA(K) = ONE

      DO L=KB1,KD,-1
        ETA(L)  = ETA(L+1) + ALM * (ETA(L) + ALM * GMS(L))
      ENDDO
      DO L=KD,KBL
        ETAI(L) = 1.0 / ETA(L)
      ENDDO

!     if (lprnt) print *,' eta=',eta,' ii=',ii,' alm=',alm
!
!===>  CLOUD WORKFUNCTION
!
      WFN    = ZERO
      AKM    = ZERO
      DET    = ZERO
      HCC    = HBL
      cnvflg = .FALSE.
      QTL    = QST(KB1) - GAF(KB1)*HST(KB1)
      TX1    = HBL
!
      qtv    = qbl
      det    = qlb + qib
!
      tx2    = 0.0
      dpneg  = 0.0
!
      DO L=KB1,KD1,-1
         DEL_ETA = ETA(L) - ETA(L+1)
         HCCP = HCC + DEL_ETA*HOL(L)
!
         QTLP = QST(L-1) - GAF(L-1)*HST(L-1)
         QTVP = 0.5 * ((QTLP+QTL)*ETA(L)                                &
     &              + (GAF(L)+GAF(L-1))*HCCP)
         ST1  = ETA(L)*Q0U(L) + ETA(L+1)*Q0D(L)
         DETP = (BKC(L)*DET - (QTVP-QTV)                                &
     &        + DEL_ETA*(QOL(L)+CLL(L)+CIL(L)) + ST1)  * AKC(L)

!     if(lprnt) print *,' detp=',detp,' bkc=',bkc(l),' det=',det
!     if (lprnt .and. kd .eq. 15) 
!    &          print *,' detp=',detp,' bkc=',bkc(l),' det=',det
!    &,' qtvp=',qtvp,' qtv=',qtv,' del_eta=',del_eta,' qol='
!    &,qol(l),' st1=',st1,' akc=',akc(l)
!
         TEM1   = AKT(L)   - QLL(L)
         TEM2   = QLL(L+1) - BKC(L)
         RNS(L) = TEM1*DETP  + TEM2*DET - ST1

         qtp    = 0.5 * (qil(L)+qil(L-1))
         tem2   = min(qtp*(detp-eta(l)*qw00),                           &
     &               (1.0-qtp)*(detp-eta(l)*qi00))
         st1    = min(tx2,tem2)
         tx2    = tem2
!
         IF (rns(l) .lt. zero .or. st1 .lt. zero) ep_wfn = .TRUE.
         IF (DETP <= ZERO) cnvflg = .TRUE.

         ST1  = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))


         TEM2 = HCCP   + DETP   * QTP * ALHF
!
!     if(lprnt) print *,' hst=',hst(l),' ltl=',ltl(l),' nu=',nu
!     if (lprnt .and. kd .eq. 15) 
!    &          print *,' hst=',hst(l),' ltl=',ltl(l),' nu=',nu
!    &,' qst=',qst(l),' qol=',qol(l),' hccp=',hccp,' detp=',detp
!    *,' qtp=',qtp,' alhf=',alhf,' vtf=',vtf(l)

         ST2  = LTL(L) * VTF(L)
         TEM5 = CLL(L) + CIL(L)
         TEM3 = (TX1  - ETA(L+1)*ST1 - ST2*(DET-TEM5*eta(l+1))) * DLB(L)
         TEM4 = (TEM2 - ETA(L  )*ST1 - ST2*(DETP-TEM5*eta(l)))  * DLT(L)
!
!     if (lprnt) then
!     if (lprnt .and. kd .eq. 12) then 
!       print *,' tem3=',tem3,' tx1=',tx1,' st1=',st1,' eta1=',eta(l+1)
!    &, ' st2=',st2,' det=',det,' tem5=',tem5,' dlb=',dlb(l)
!       print *,' tem4=',tem4,' tem2=',tem2,' detp=',detp
!    &, ' eta=',eta(l),' dlt=',dlt(l),' rns=',rns(l),' l=',l
!       print *,' bt1=',tem3/(eta(l+1)*qrb(l))
!    &,         ' bt2=',tem4/(eta(l)*qrt(l))
!      endif

         ST1  = TEM3 + TEM4

!     if (lprnt) print *,' wfn=',wfn,' st1=',st1,' l=',l,' ep_wfn=',
!    &ep_wfn,' akm=',akm

         WFN = WFN + ST1       
         AKM = AKM - min(ST1,ZERO)

!     if (lprnt) print *,' wfn=',wfn,' akm=',akm

         if (st1 .lt. zero .and. wfn .lt. zero) then
           dpneg = dpneg + prl(l+1) - prl(l)
         endif

         BUY(L) = 0.5 * (tem3/(eta(l+1)*qrb(l)) + tem4/(eta(l)*qrt(l)))
!
         HCC = HCCP
         DET = DETP
         QTL = QTLP
         QTV = QTVP
         TX1 = TEM2

      ENDDO

      DEL_ETA = ETA(KD) - ETA(KD1)
      HCCP    = HCC + DEL_ETA*HOS
!
      QTLP    = QST(KD) - GAF(KD)*HST(KD)
      QTVP    = QTLP*ETA(KD) + GAF(KD)*HCCP
      ST1     = ETA(KD)*Q0U(KD) + ETA(KD1)*Q0D(KD)
      DETP    = (BKC(KD)*DET - (QTVP-QTV)                               &
     &        + DEL_ETA*(QOS+QLS+QIS) + ST1) * AKC(KD)
!
      TEM1    = AKT(KD)  - QLL(KD)
      TEM2    = QLL(KD1) - BKC(KD)
      RNS(KD) = TEM1*DETP  + TEM2*DET - ST1
!
      IF (rns(kd) .lt. zero) ep_wfn = .TRUE.
      IF (DETP <= ZERO) cnvflg = .TRUE.
!
  888 continue

!     if (lprnt) print *,' ep_wfn=',ep_wfn,' ii=',ii,' rns=',rns(kd)
!    &,' clp=',clp,' hst(kd)=',hst(kd)

      if (ep_wfn) then
        IF ((qw00 .eq. 0.0 .and. qi00 .eq. 0.0)) RETURN
        if (ii .eq. 0) then
          ii  = 1
          if (clp .gt. 0.0 .and. clp .lt. 1.0) then
            hst(kd) = hstkd
            qst(kd) = qstkd
            ltl(kd) = ltlkd
            q0u(kd) = q0ukd
            q0d(kd) = q0dkd
            dlb(kd) = dlbkd
            qrb(kd) = qrbkd
          endif
          do l=kd,kb1
            FCO(L) = FCO(L) - q0u(l) - q0d(l)
            RNN(L) = RNN(L) - q0u(l)*zet(l) - q0d(l)*zet(l+1)
            GMH(L) = GMH(L) - q0u(l)*xi(l)  - q0d(l)*zet(l+1)
            ETA(L) = ZET(L) - ZET(L+1)
            GMS(L) = XI(L)  - XI(L+1)
            Q0U(L) = 0.0
            Q0D(L) = 0.0
          ENDDO
          qw00 = 0.0
          qi00 = 0.0

!     if (lprnt) print *,' returning to 777 : ii=',ii,' qw00=',qw00,qi00
!    &,' clp=',clp,' hst(kd)=',hst(kd)

          go to 777
        else
          cnvflg = .true.
        endif
      endif
!
!
!     ST1 = 0.5 * (HST(KD)  - LTL(KD)*NU*(QST(KD)-QOS)
!    &          +  HST(KD1) - LTL(KD1)*NU*(QST(KD1)-QOL(KD1)))
!
      ST1 = HST(KD)  - LTL(KD)*NU*(QST(KD)-QOS)
      ST2 = LTL(KD)  * VTF(KD)
      TEM5 = (QLS + QIS) * eta(kd1)
      ST1  = HALF * (TX1-ETA(KD1)*ST1-ST2*(DET-TEM5))*DLB(KD)
!
!     if (lprnt) print *,' st1=',st1,' st2=',st2,' ltl=',ltl(kd)
!    *,ltl(kd1),' qos=',qos,qol(kd1)

      WFN = WFN + ST1
      AKM = AKM - min(ST1,ZERO)   ! Commented on 08/26/02 - does not include top
!

      BUY(KD) = ST1 / (ETA(KD1)*qrb(kd))
!
!     if (lprnt) print *,' wfn=',wfn,' akm=',akm,' st1=',st1
!    &,' dpneg=',dpneg

      DET = DETP
      HCC = HCCP
      AKM = AKM / WFN


!***********************************************************************
!
!     If only to calculate workfunction save it and return
!
      IF (WRKFUN) THEN
        IF (WFN .GE. 0.0) WFNC = WFN
        RETURN
      ELSEIF (.NOT. CRTFUN) THEN
        ACR = WFNC
      ENDIF
!
!===>  THIRD CHECK BASED ON CLOUD WORKFUNCTION
!
      CALCUP = .FALSE.

      TEM  =  max(0.05, MIN(CD*200.0, MAX_NEG_BOUY))
      IF (WFN .GT. ACR .AND.  (.NOT. cnvflg)                             &
!    & .and. dpneg .lt. 100.0  .AND. AKM .LE. TEM) THEN
     & .and. dpneg .lt. 150.0  .AND. AKM .LE. TEM) THEN
!    & .and. dpneg .lt. 200.0  .AND. AKM .LE. TEM) THEN
!
        CALCUP = .TRUE.
      ENDIF

!     if (lprnt) print *,' calcup=',calcup,' akm=',akm,' tem=',tem
!    *,' cnvflg=',cnvflg,' clp=',clp,' rhc=',rhc,' cd=',cd,' acr=',acr
!
!===>  IF NO SOUNDING MEETS THIRD CONDITION, RETURN
!
      IF (.NOT. CALCUP) RETURN
!
! This is for not LL - 20050601
      IF (ALMIN2 .NE. 0.0) THEN
        IF (ALMIN1 .NE. ALMIN2) ST1 = 1.0 / max(ONE_M10,(ALMIN2-ALMIN1))
        IF (ALM .LT. ALMIN2) THEN
           CLP = CLP * max(0.0, min(1.0,(0.3 + 0.7*(ALM-ALMIN1)*ST1)))
!          CLP = CLP * max(0.0, min(1.0,(0.2 + 0.8*(ALM-ALMIN1)*ST1)))
!          CLP = CLP * max(0.0, min(1.0,(0.1 + 0.9*(ALM-ALMIN1)*ST1)))
        ENDIF
      ENDIF
!
!     if (lprnt) print *,' clp=',clp
!
      CLP = CLP * RHC
      dlq = 0.0
      tem = 1.0 / (1.0 + dlq_fac)
      do l=kd,kb1
        rnn(l) = rns(l) * tem
        dlq(l) = rns(l) * tem * dlq_fac
      enddo
      DO L=KBL,K 
        RNN(L) = 0.0 
      ENDDO
!     if (lprnt) print *,' rnn=',rnn
!
!     If downdraft is to be invoked, do preliminary check to see
!     if enough rain is available and then call DDRFT.
!
      DDFT = .FALSE.
      IF (DNDRFT) THEN
!
        TRAIN = 0.0
        IF (CLP .GT. 0.0) THEN
          DO L=KD,KB1
            TRAIN = TRAIN + RNN(L)
          ENDDO
        ENDIF

        PL = (PRL(KD1) + PRL(KD))*HALF
        TEM = PRL(K+1)*(1.0-DPD*0.001)
        IF (TRAIN .GT. 1.0E-4 .AND. PL .LE. TEM) DDFT  = .TRUE.
!
      ENDIF
!
!     if (lprnt) then
!       print *,' BEFORE CALLING DDRFT KD=',kd,' DDFT=',DDFT
!    &,                  ' PL=',PL,' TRAIN=',TRAIN
!       print *,' buy=',(buy(l),l=kd,kb1)
!     endif

      IF (DDFT) THEN ! Downdraft scheme based on (Cheng and Arakawa, 1997)
        CALL DDRFT(                                                     &
     &              K, KD                                               &
     &,             TLA, ALFIND                                         &
     &,             TOL, QOL, HOL,   PRL, QST, HST, GAM, GAF            &
!    &,             TOL, QOL, HOL,   PRL, QST, HST, GAM, GAF, HBL, QBL  &
     &,             QRB, QRT, BUY,   KBL, IDH, ETA, RNN, ETAI           &
     &,             ALM, WFN, TRAIN, DDFT                               &
     &,             ETD, HOD, QOD,   EVP, DOF, CLDFR, ETZ               &
     &,             GMS, GSD, GHD,   lprnt)               

      ENDIF
!
!  No Downdraft case (including case with no downdraft solution)
!  ---------------------------------------------------------
!
      IF (.NOT. DDFT) THEN
        DO L=KD,K+1
          ETD(L) = 0.0
          HOD(L) = 0.0
          QOD(L) = 0.0
        ENDDO
        DO L=KD,K
          EVP(L) = 0.0
          ETZ(L) = 0.0
        ENDDO

      ENDIF

!     if (lprnt) print *,' hod=',hod
!     if (lprnt) print *,' etd=',etd
!
!
!===> CALCULATE GAMMAS  i.e. TENDENCIES PER UNIT CLOUD BASE MASSFLUX
!           Includes downdraft terms!

      avh = 0.0

!
!     Fraction of detrained condensate evaporated
!
!     tem1 = max(ZERO, min(HALF, (prl(kd)-FOUR_P2)*ONE_M2))
!     tem1 = max(ZERO, min(HALF, (prl(kd)-300.0)*0.005))
      tem1 = 0.0
!     tem1 = 1.0
!     if (kd1 .eq. kbl) tem1 = 0.0
!
      tem2    = 1.0 - tem1
      TEM = DET * QIL(KD)


      st1 = (HCC+ALHF*TEM-ETA(KD)*HST(KD)) / (1.0+gam(KD))
      DS  = ETA(KD1) * (HOS- HOL(KD)) - ALHL*(QOS - QOL(KD))
      DH  = ETA(KD1) * (HOS- HOL(KD))


      GMS(KD) = (DS + st1 - tem1*det*alhl-tem*alhf) * PRI(KD)
      GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOS + DH)


!     if (lprnt) print *,' gmhkd=',gmh(kd),' gmskd=',gms(kd)
!    &,' det=',det,' tem=',tem,' tem1=',tem1,' tem2=',tem2
!
!      TENDENCY FOR SUSPENDED ENVIRONMENTAL ICE AND/OR LIQUID WATER
!
      QLL(KD) = (tem2*(DET-TEM) + ETA(KD1)*(QLS-CLL(KD))                &
     &        + (1.0-QIL(KD))*dlq(kd) - ETA(KD)*QLS ) * PRI(KD)

      QIL(KD) =     (tem2*TEM + ETA(KD1)*(QIS-CIL(KD))                  &
     &        + QIL(KD)*dlq(kd) - ETA(KD)*QIS ) * PRI(KD)
!
      GHD(KD) = 0.0
      GSD(KD) = 0.0
!
      DO L=KD1,K
         ST1 = ONE - ALFINT(L,1)
         ST2 = ONE - ALFINT(L,2)
         ST3 = ONE - ALFINT(L,3)
         ST4 = ONE - ALFINT(L,4)
         ST5 = ONE - ALFIND(L)
         HB       = ALFINT(L,1)*HOL(L-1) + ST1*HOL(L)
         QB       = ALFINT(L,2)*QOL(L-1) + ST2*QOL(L)

         TEM      = ALFINT(L,4)*CIL(L-1) + ST4*CIL(L)
         TEM2     = ALFINT(L,3)*CLL(L-1) + ST3*CLL(L)
 
         TEM1     = ETA(L) * (TEM - CIL(L))
         TEM3     = ETA(L) * (TEM2 - CLL(L))

         HBD      = ALFIND(L)*HOL(L-1) + ST5*HOL(L)
         QBD      = ALFIND(L)*QOL(L-1) + ST5*QOL(L)

         TEM5     = ETD(L) * (HOD(L) - HBD)
         TEM6     = ETD(L) * (QOD(L) - QBD)
!
         DH       = ETA(L) * (HB - HOL(L)) + TEM5
         DS       = DH - ALHL * (ETA(L) * (QB - QOL(L)) + TEM6)

         GMH(L)   = DH * PRI(L)
         GMS(L)   = DS * PRI(L)

!     if (lprnt) print *,' gmh=',gmh(l),' gms=',gms(l)
!    &,' dh=',dh,' ds=',ds,' qb=',qb,' qol=',qol(l),' eta=',eta(l)
!    &,' hb=',hb,' hol=',hol(l),' l=',l,' hod=',hod(l)
!    &,' etd=',etd(l),' qod=',qod(l),' tem5=',tem5,' tem6=',tem6
!
         GHD(L)   = TEM5 * PRI(L)
         GSD(L)   = (TEM5 - ALHL * TEM6) * PRI(L)
!
         QLL(L)   = (TEM3 + (1.0-QIL(L))*dlq(l)) * PRI(L)
         QIL(L)   = (TEM1 + QIL(L)*dlq(l)) * PRI(L)

         TEM1     = ETA(L) * (CIL(L-1) - TEM)
         TEM3     = ETA(L) * (CLL(L-1) - TEM2)

         DH       = ETA(L) * (HOL(L-1) - HB) - TEM5
         DS       = DH - ALHL * ETA(L) * (QOL(L-1) - QB)                &
     &                 + ALHL * (TEM6 - EVP(L-1))

         GMH(L-1) = GMH(L-1) + DH * PRI(L-1)
         GMS(L-1) = GMS(L-1) + DS * PRI(L-1)
!
!     if (lprnt) print *,' gmh1=',gmh(l-1),' gms1=',gms(l-1)
!    &,' dh=',dh,' ds=',ds,' qb=',qb,' qol=',qol(l-1)
!    &,' hb=',hb,' hol=',hol(l-1),' evp=',evp(l-1)
!
         GHD(L-1) = GHD(L-1) - TEM5 * PRI(L-1)
         GSD(L-1) = GSD(L-1) - (TEM5-ALHL*(TEM6-EVP(L-1))) * PRI(L-1)

         QIL(L-1) = QIL(L-1) + TEM1 * PRI(L-1)
         QLL(L-1) = QLL(L-1) + TEM3 * PRI(L-1)


!     if (lprnt) print *,' gmh=',gmh(l),' gms=',gms(l)
!    &,' dh=',dh,' ds=',ds,' qb=',qb,' qol=',qol(l),' eta=',eta(l)
!    &,' hb=',hb,' hol=',hol(l),' l=',l
!
        avh = avh + gmh(l-1)*(prs(l)-prs(l-1))

      ENDDO
!
      HBD    = HOL(K)
      QBD    = QOL(K)
      TEM5   =  ETD(K+1) * (HOD(K+1) - HBD)
      TEM6   =  ETD(K+1) * (QOD(K+1) - QBD)
      DH     = - TEM5
      DS     = DH  + ALHL * TEM6
      TEM1   = DH * PRI(K)
      TEM2   = (DS - ALHL * EVP(K)) * PRI(K)
      GMH(K) = GMH(K) + TEM1
      GMS(K) = GMS(K) + TEM2
      GHD(K) = GHD(K) + TEM1
      GSD(K) = GSD(K) + TEM2

!     if (lprnt) print *,' gmhk=',gmh(k),' gmsk=',gms(k)
!    &,' tem1=',tem1,' tem2=',tem2,' dh=',dh,' ds=',ds
!
      avh    = avh + gmh(K)*(prs(KP1)-prs(K))
!
      tem4   = - GRAVFAC * pris
      TX1    = DH * tem4
      TX2    = DS * tem4
!
      DO L=KBL,K
        GMH(L) = GMH(L) + TX1
        GMS(L) = GMS(L) + TX2
        GHD(L) = GHD(L) + TX1
        GSD(L) = GSD(L) + TX2
!
        avh = avh + tx1*(prs(l+1)-prs(l))
      ENDDO

!
!     if (lprnt) then
!        print *,' gmh=',gmh
!        print *,' gms=',gms(KD:K)
!     endif
!
!***********************************************************************
!***********************************************************************

!===>  KERNEL (AKM) CALCULATION BEGINS

!===>  MODIFY SOUNDING WITH UNIT MASS FLUX
!
      DO L=KD,K

         TEM1   = GMH(L)
         TEM2   = GMS(L)
         HOL(L) = HOL(L) +  TEM1*TESTMB
         QOL(L) = QOL(L) + (TEM1-TEM2)  * (TESTMB/ALHL)
         HST(L) = HST(L) +  TEM2*(ONE+GAM(L))*TESTMB
         QST(L) = QST(L) +  TEM2*GAM(L)*(TESTMB/ALHL)
         CLL(L) = CLL(L) + QLL(L) * TESTMB
         CIL(L) = CIL(L) + QIL(L) * TESTMB
      ENDDO
!
      if (alm .gt. 0.0) then
        HOS     = HOS + GMH(KD)  * TESTMB
        QOS     = QOS + (GMH(KD)-GMS(KD)) * (TESTMB/ALHL)
        QLS     = QLS + QLL(KD) * TESTMB
        QIS     = QIS + QIL(KD) * TESTMB
      else
        st2 = 1.0 - st1s
        HOS     = HOS + (st1s*GMH(KD)+st2*GMH(KD1))  * TESTMB
        QOS     = QOS + (st1s * (GMH(KD)-GMS(KD))                       &
     &                +  st2  * (GMH(KD1)-GMS(KD1))) * (TESTMB/ALHL)
        HST(kd) = HST(kd) + (st1s*GMS(kd)*(ONE+GAM(kd))                 &
     &                    +  st2*gms(kd1)*(ONE+GAM(kd1))) * TESTMB
        QST(kd) = QST(kd) + (st1s*GMS(kd)*GAM(kd)                       &
     &                    +  st2*gms(kd1)*gam(kd1)) * (TESTMB/ALHL)

        QLS     = QLS + (st1s*QLL(KD)+st2*QLL(KD1)) * TESTMB
        QIS     = QIS + (st1s*QIL(KD)+st2*QIL(KD1)) * TESTMB
      endif

!
      TEM = PRL(Kmaxp1) - PRL(Kmax)
      HBL = HOL(Kmax) * TEM
      QBL = QOL(Kmax) * TEM
      QLB = CLL(Kmax) * TEM
      QIB = CIL(Kmax) * TEM
      DO L=KmaxM1,KBL,-1
        TEM = PRL(L+1) - PRL(L)
        HBL = HBL + HOL(L) * TEM
        QBL = QBL + QOL(L) * TEM
        QLB = QLB + CLL(L) * TEM
        QIB = QIB + CIL(L) * TEM
      ENDDO
      HBL = HBL * PRISM
      QBL = QBL * PRISM
      QLB = QLB * PRISM
      QIB = QIB * PRISM

!     if (ctei .and. sgcs(kd) > 0.65) then
!        hbl = hbl * hpert_fac
!        qbl = qbl * hpert_fac
!     endif

!     if (lprnt) print *,' hbla=',hbl,' qbla=',qbl

!***********************************************************************

!===>  CLOUD WORKFUNCTION FOR MODIFIED SOUNDING, THEN KERNEL (AKM)
!
      AKM = ZERO
      TX1 = ZERO
      QTL = QST(KB1) - GAF(KB1)*HST(KB1)
      QTV = QBL
      HCC = HBL
      TX2 = HCC
      TX4 = (ALHF*0.5)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(KB1))*TCRF))
!
      qtv = qbl
      tx1 = qib + qlb
!

      DO L=KB1,KD1,-1
         DEL_ETA = ETA(L) - ETA(L+1)
         HCCP = HCC + DEL_ETA*HOL(L)
!
         QTLP = QST(L-1) - GAF(L-1)*HST(L-1)
         QTVP = 0.5 * ((QTLP+QTL)*ETA(L) + (GAF(L)+GAF(L-1))*HCCP)

         DETP = (BKC(L)*TX1 - (QTVP-QTV)                                &
     &        +  DEL_ETA*(QOL(L)+CLL(L)+CIL(L))                         &
     &        +  ETA(L)*Q0U(L) + ETA(L+1)*Q0D(L)) * AKC(L)
         IF (DETP .LE. ZERO) cnvflg = .TRUE.

         ST1  = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))

         TEM2 = (ALHF*0.5)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(L-1))*TCRF))
         TEM1 = HCCP + DETP * (TEM2+TX4)

         ST2  = LTL(L) * VTF(L)
         TEM5 = CLL(L) + CIL(L)
         AKM  = AKM +                                                   &
     &     (  (TX2  -ETA(L+1)*ST1-ST2*(TX1-TEM5*eta(l+1))) * DLB(L)     &
     &      + (TEM1 -ETA(L  )*ST1-ST2*(DETP-TEM5*eta(l)))  * DLT(L) )
!
         HCC  = HCCP
         TX1  = DETP
         TX2  = TEM1
         QTL  = QTLP
         QTV  = QTVP
         TX4  = TEM2
      ENDDO
!
      if (cnvflg) return
!
!  Eventhough we ignore the change in lambda, we still assume
!  that the cLoud-top contribution is zero; as though we still
!  had non-bouyancy there.
!
!
      ST1 = HST(KD)  - LTL(KD)*NU*(QST(KD)-QOS)
      ST2 = LTL(KD)  * VTF(KD)
      TEM5 = (QLS + QIS) * eta(kd1)
      AKM  = AKM + HALF * (TX2-ETA(KD1)*ST1-ST2*(TX1-TEM5)) * DLB(KD)
!
      AKM = (AKM - WFN) * (ONE/TESTMB)


!***********************************************************************

!===>   MASS FLUX

      tem2 = rel_fac
!
      AMB = - (WFN-ACR) / AKM
!
!     if(lprnt) print *,' wfn=',wfn,' acr=',acr,' akm=',akm             &
!    &,' amb=',amb,' KD=',kd,' cldfrd=',cldfrd,' tem2=',tem2            &
!    &,' rel_fac=',rel_fac,' prskd=',prs(kd)

!===>   RELAXATION AND CLIPPING FACTORS
!
      AMB = AMB * CLP * tem2

!!!   if (DDFT) AMB = MIN(AMB, ONE/CLDFRD)
       
!===>   SUB-CLOUD LAYER DEPTH LIMIT ON MASS FLUX

      AMBMAX = (PRL(KMAXP1)-PRL(KBL))*(FRACBL*GRAVCON)
      AMB    = MAX(MIN(AMB, AMBMAX),ZERO)


!     if(lprnt) print *,' AMB=',amb,' clp=',clp,' ambmax=',ambmax
!***********************************************************************
!*************************RESULTS***************************************
!***********************************************************************

!===>  PRECIPITATION AND CLW DETRAINMENT
!
      avt = 0.0
      avq = 0.0
      avr = dof

!
      DSFC = DSFC + AMB * ETD(K) * (1.0/DT)
!
!     DO L=KBL,KD,-1
      DO L=K,KD,-1
          PCU(L) = PCU(L) + AMB*RNN(L)      !  (A40)
          avr    = avr + rnn(l)
!     if(lprnt) print *,' avr=',avr,' rnn=',rnn(l),' l=',l
      ENDDO
!
!===> TEMPARATURE AND Q CHANGE AND CLOUD MASS FLUX DUE TO CLOUD TYPE KD
!
      TX1  = AMB * (ONE/CP)
      TX2  = AMB * (ONE/ALHL)
      DO L=KD,K
        ST1    = GMS(L)*TX1
        TOI(L) = TOI(L) + ST1
        TCU(L) = TCU(L) + ST1
        TCD(L) = TCD(L) + GSD(L) * TX1
!
        st1 = st1 - (alhl/cp) * (QIL(L) + QLL(L)) * AMB

        avt = avt + st1 * (prs(l+1)-prs(l))

        FLX(L)  = FLX(L)  + ETA(L)*AMB
        FLXD(L) = FLXD(L) + ETD(L)*AMB
!
        QII(L)  = QII(L) + QIL(L) * AMB
        TEM     = 0.0

        QLI(L)  = QLI(L) + QLL(L) * AMB + TEM

        ST1     = (GMH(L)-GMS(L)) * TX2

        QOI(L)  = QOI(L) + ST1
        QCU(L)  = QCU(L) + ST1
        QCD(L)  = QCD(L) + (GHD(L)-GSD(L)) * TX2
!
        avq = avq + (st1+(QLL(L)+QIL(L))*amb) * (prs(l+1)-prs(l))
!       avq = avq + st1 * (prs(l+1)-prs(l))
!       avr = avr + (QLL(L) + QIL(L)*(1+alhf/alhl))
!       avr = avr + (QLL(L) + QIL(L))
!    *                  * (prs(l+1)-prs(l)) * gravcon

!     if(lprnt) print *,' avr=',avr,' qll=',qll(l),' l=',l
!    &,' qil=',qil(l)

      ENDDO
      avr = avr * amb
!
!      Correction for negative condensate!
!     if (advcld) then
!       do l=kd,k
!         if (qli(l) .lt. 0.0) then
!           qoi(l) = qoi(l) + qli(l)
!           toi(l) = toi(l) - (alhl/cp) * qli(l)
!           qli(l) = 0.0
!         endif
!         if (qii(l) .lt. 0.0) then
!           qoi(l) = qoi(l) + qii(l)
!           toi(l) = toi(l) - ((alhl+alhf)/cp) * qii(l)
!           qii(l) = 0.0
!         endif
!       enddo
!     endif

!
!
!     if (lprnt) then
!       print *,' For KD=',KD
!       avt = avt * cp * 100.0*86400.0 / (alhl*DT*grav)
!       avq = avq *  100.0*86400.0 / (DT*grav)
!       avr = avr * 86400.0 / DT
!       print *,' avt=',avt,' avq=',avq,' avr=',avr,' avh='
!    *   ,avh,' alm=',alm,' DDFT=',DDFT,' KD=',KD
!    &,' TOIK-',toi(k),' TOIK-1=',toi(k-1),' TOIK-2=',toi(k-2)
!        if (kd .eq. 12 .and. .not. ddft) stop
!       if (avh .gt. 0.1 .or. abs(avt+avq) .gt. 1.0e-5 .or.
!    &      abs(avt-avr) .gt. 1.0e-5 .or. abs(avr+avq) .gt. 1.0e-5) stop
!
!     if (lprnt) then
!       print *,' For KD=',KD
!       print *,' TCU=',(tcu(l),l=kd,k)
!       print *,' QCU=',(Qcu(l),l=kd,k)
!     endif
!
      TX1 = 0.0
      TX2 = 0.0
!
      IF (REVAP) THEN !     REEVAPORATION OF FALLING CONVECTIVE RAIN
!
       tem = 0.0
       do l=kd,kbl
         IF (L .lt. IDH .or. (.not. DDFT)) THEN
           tem = tem + amb * rnn(l)
         endif
       enddo
       tem = tem + amb * dof
       tem = tem * (3600.0/dt)
!!!!   tem1 = max(1.0, min(100.0,sqrt((5.0E10/max(garea,one)))))
!      tem1 = max(1.0, min(100.0,(7.5E10/max(garea,one))))
!      tem1 = max(1.0, min(100.0,(5.0E10/max(garea,one))))
!      tem1 = max(1.0, min(100.0,(4.0E10/max(garea,one))))
!!     tem1 = sqrt(max(1.0, min(100.0,(4.0E10/max(garea,one))))) ! 20100902
       tem1 = sqrt(max(1.0, min(100.0,(6.25E10/max(garea,one))))) ! 20110530

!      if (lprnt) print *,' clfr0=',clf(tem),' tem=',tem,' tem1=',tem1

!      clfrac = max(ZERO, min(ONE, rknob*clf(tem)*tem1))
!      clfrac = max(ZERO, min(0.25, rknob*clf(tem)*tem1))
       clfrac = max(ZERO, min(half, rknob*clf(tem)*tem1))

!      if (lprnt) then
!        print *,' cldfrd=',cldfrd,' amb=',amb,' clfrac=',clfrac
!        print *,' tx3=',tx3,' etakd=',eta(kd),' pri=',pri(kd)
!        print *,' RNN=',RNN(kd:k)
!      endif
!
!cnt   DO L=KD,K
       DO L=KD,KBL         ! Testing on 20070926
!                                                 for L=KD,K
         IF (L .GE. IDH .AND. DDFT) THEN
           TX2    = TX2 + AMB * RNN(L)
           CLDFRD = MIN(AMB*CLDFR(L), clfrac)
         ELSE
           TX1 = TX1 + AMB * RNN(L)
         ENDIF
         tx4 = zfac * phil(l)
         tx4 = (one - tx4 * (one - half*tx4)) * afc
!
         IF (TX1 .GT. 0. .OR. TX2 .GT. 0.0) THEN
          TEQ     = TOI(L)
          QEQ     = QOI(L)
          PL      = 0.5 * (PRL(L+1)+PRL(L))

          ST1     = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
          ST2     = ST1*ELFOCP + (1.0-ST1)*ELOCP

          CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
!         CALL QSATCN ( TEQ,PL,QSTEQ,DQDT,.false.)
!
          DELTAQ = 0.5 * (QSTEQ*rhc_ls(l)-QEQ) / (1.+ST2*DQDT)
!
          QEQ    = QEQ + DELTAQ
          TEQ    = TEQ - DELTAQ*ST2
!
          TEM1   = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
          TEM2   = TEM1*ELFOCP + (1.0-TEM1)*ELOCP

          CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
!         CALL QSATCN ( TEQ,PL,QSTEQ,DQDT,.false.)
!
          DELTAQ = (QSTEQ*rhc_ls(l)-QEQ) / (1.+TEM2*DQDT)
!
          QEQ    = QEQ + DELTAQ
          TEQ    = TEQ - DELTAQ*TEM2

          IF (QEQ .GT. QOI(L)) THEN
            POTEVAP = (QEQ-QOI(L))*(PRL(L+1)-PRL(L))*GRAVCON

            tem4    = 0.0
            if (tx1 .gt. 0.0)                                           &
     &      TEM4    = POTEVAP * (1. - EXP( tx4*TX1**0.57777778 ) )
!    &      TEM4    = POTEVAP * (1. - EXP( AFC*tx4*SQRT(TX1) ) )
            ACTEVAP = MIN(TX1, TEM4*CLFRAC)

!     if(lprnt) print *,' L=',L,' actevap=',actevap,' tem4=',tem4,
!    &' clfrac='
!    &,clfrac,' potevap=',potevap,'efac=',AFC*SQRT(TX1*TEM3)
!    &,' tx1=',tx1

            if (tx1 .lt. rainmin*dt) actevap = min(tx1, potevap)
!
            tem4    = 0.0
            if (tx2 .gt. 0.0)                                           &
     &      TEM4    = POTEVAP * (1. - EXP( tx4*TX2**0.57777778 ) )
!    &      TEM4    = POTEVAP * (1. - EXP( AFC*tx4*SQRT(TX2) ) )
            TEM4    = min(MIN(TX2, TEM4*CLDFRD), potevap-actevap)
            if (tx2 .lt. rainmin*dt) tem4 = min(tx2, potevap-actevap)
!
            TX1     = TX1 - ACTEVAP
            TX2     = TX2 - TEM4
            ST1     = (ACTEVAP+TEM4) * PRI(L)
            QOI(L)  = QOI(L) + ST1
            QCU(L)  = QCU(L) + ST1
!

            ST1     = ST1 * ELOCP
            TOI(L)  = TOI(L) - ST1 
            TCU(L)  = TCU(L) - ST1
          ENDIF
         ENDIF
       ENDDO
!
       CUP = CUP + TX1 + TX2 + DOF * AMB
      ELSE
       DO L=KD,K
         TX1 = TX1 + AMB * RNN(L)
       ENDDO
       CUP = CUP + TX1 + DOF * AMB
      ENDIF

!     if (lprnt) print *,' tx1=',tx1,' tx2=',tx2,' dof=',dof
!    &,' cup=',cup*86400/dt,' amb=',amb
!    &,' amb=',amb,' cup=',cup,' clfrac=',clfrac,' cldfrd=',cldfrd
!    &,' ddft=',ddft,' kd=',kd,' kbl=',kbl,' k=',k
!
!    Convective transport (mixing) of passive tracers
!
      if (NTRC > 0) then
        do l=kd,k-1
          if (etz(l) .ne. zero) etzi(l) = one / etz(l)
        enddo
        DO N=1,NTRC        ! Tracer loop ; first two are u and v

          DO L=KD,K
            HOL(L) = ROI(L,N)
          ENDDO
!
          HCC     = RBL(N)
          HOD(KD) = HOL(KD)
!      Compute downdraft properties for the tracer
          DO L=KD1,K
            ST1 = ONE - ALFIND(L)
            HB  = ALFIND(L)  * HOL(L-1) + ST1 * HOL(L)
            IF (ETZ(L-1) .NE. ZERO) THEN
              TEM = ETZI(L-1)
              IF (ETD(L)  > ETD(L-1)) THEN
                HOD(L) = (ETD(L-1)*(HOD(L-1)-HOL(L-1))                  &
     &                 +  ETD(L)  *(HOL(L-1)-HB) +  ETZ(L-1)*HB) * TEM
              ELSE
                HOD(L) = (ETD(L-1)*(HOD(L-1)-HB) + ETZ(L-1)*HB) * TEM
              ENDIF
            ELSE
              HOD(L) = HB
            ENDIF
          ENDDO
             
          DO L=KB1,KD,-1
            HCC = HCC + (ETA(L)-ETA(L+1))*HOL(L)
          ENDDO
!
!         Scavenging -- fscav   - fraction scavenged [km-1]
!                       delz    - distance from the entrainment to detrainment layer [km]
!                       fnoscav - the fraction not scavenged
!                                 following Liu et al. [JGR,2001] Eq 1

          if (FSCAV_(N) > 0.0) then
            DELZKM = ( PHIL(KD) - PHIH(KD1) ) *(onebg*0.001)
            FNOSCAV = exp(- FSCAV_(N) * DELZKM)
          else
            FNOSCAV = 1.0
          endif

          GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOL(KD)) * trcfac(kd,n)      &
     &                                              * FNOSCAV
          DO L=KD1,K
            if (FSCAV_(N) > 0.0) then
              DELZKM = ( PHIL(KD) - PHIH(L+1) ) *(onebg*0.001)
              FNOSCAV = exp(- FSCAV_(N) * DELZKM)
            endif
            ST1      = ONE - ALFINT(L,N+4)
            ST2      = ONE - ALFIND(L)
            HB       = ALFINT(L,N+4) * HOL(L-1) + ST1 * HOL(L)
            HBD      = ALFIND(L) * HOL(L-1) + ST2 * HOL(L)
            TEM5     = ETD(L)    * (HOD(L) - HBD)
            DH       = ETA(L)    * (HB - HOL(L))   * FNOSCAV + TEM5
            GMH(L  ) = DH * PRI(L) * trcfac(l,n)
            DH       = ETA(L)    * (HOL(L-1) - HB) * FNOSCAV - TEM5
            GMH(L-1) = GMH(L-1)  + DH * PRI(L-1) * trcfac(l,n)
          ENDDO
!
          DO L=KD,K
            ST1      = GMH(L)*AMB
            ROI(L,N) = HOL(L)   + ST1
            RCU(L,N) = RCU(L,N) + ST1
          ENDDO
        ENDDO                             ! Tracer loop NTRC
      endif

!     if (lprnt) print *,' toio=',toi
!     if (lprnt) print *,' qoio=',qoi

      RETURN
      END

      SUBROUTINE DDRFT(                                                 &
     &                  K, KD                                           &
     &,                 TLA, ALFIND                                     &
     &,                 TOL, QOL, HOL, PRL, QST, HST, GAM, GAF          &
!    &,                 TOL, QOL, HOL, PRL, QST, HST, GAM, GAF, HBL, QBL&
     &,                 QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI         &
     &,                 ALM, WFN, TRAIN, DDFT                           &
     &,                 ETD, HOD, QOD, EVP, DOF, CLDFRD, WCB            &
     &,                 GMS, GSD, GHD,lprnt)                   

!
!***********************************************************************
!******************** Cumulus Downdraft Subroutine *********************
!****************** Based on Cheng and Arakawa (1997)  ****** **********
!************************ SUBROUTINE DDRFT  ****************************
!*************************  October 2004  ******************************
!***********************************************************************
!***********************************************************************
!************* Shrinivas.Moorthi@noaa.gov (301) 763 8000(X7233) ********
!***********************************************************************
!***********************************************************************
!23456789012345678901234567890123456789012345678901234567890123456789012
!
!===>  TOL(K)     INPUT   TEMPERATURE            KELVIN
!===>  QOL(K)     INPUT   SPECIFIC HUMIDITY      NON-DIMENSIONAL

!===>  PRL(K+1)   INPUT   PRESSURE @ EDGES       MB

!===>  K     INPUT   THE RISE & THE INDEX OF THE SUBCLOUD LAYER
!===>  KD    INPUT   DETRAINMENT LEVEL ( 1<= KD < K )          
!     
      USE MACHINE , ONLY : kind_phys
      use module_ras
      IMPLICIT NONE
!
!  INPUT ARGUMENTS
!
      INTEGER K, KD
      real(kind=kind_phys) ALFIND(K)
      INTEGER KBL, KB1

      LOGICAL SKPDD, SKPUP

      real(kind=kind_phys) HOL(KD:K),   QOL(KD:K),   GAF(KD:K+1)        &
     &,                    HST(KD:K),   QST(KD:K),   TOL(KD:K)          &
     &,                    BUY(KD:K+1), QRB(KD:K),   QRT(KD:K)          &
     &,                    GAM(KD:K+1), RNN(KD:K),   RNS(KD:K)                       &
     &,                    ETA(KD:K+1), PRL(KD:K+1), ETAI(KD:K)
!
!     real(kind=kind_phys)    HBL,     QBL,        PRIS                 &
!    &,                       TRAIN,   WFN,        ALM
      real(kind=kind_phys)    TRAIN,   WFN,        ALM
!
!     TEMPORARY WORK SPACE
!
      real(kind=kind_phys) GMS(KD:K+1)
      real(kind=kind_phys) TX1,    TX2,  TX3, TX4                       &
     &,                    TX5,    TX6,  TX7, TX8, TX9
      LOGICAL cnvflg

      real(kind=kind_phys) TL, PL, QL, QS, DQS, ST1,  HB, QB, TB        &
     &,                    QQQ, PICON, PIINV, DEL_ETA                   &
     &,                    TEM, TEM1, TEM2, TEM3, TEM4, ST2             &
     &,                    ERRMIN, ERRMI2, ERRH, ERRW, ERRE, TEM5       &
     &,                    TEM6, HBD, QBD
      INTEGER I, L,  N, IX, KD1, II                                     &
     &,       KP1, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1                 &
     &,       IP1, JJ, ntla

!
      integer, parameter :: NUMTLA=2
!     integer, parameter :: NUMTLA=4
      parameter (ERRMIN=0.0001, ERRMI2=0.1*ERRMIN)
!     parameter (ERRMIN=0.00001, ERRMI2=0.1*ERRMIN)
!
      real(kind=kind_phys) TLA,    STLA,  CTL2, CTL3
      real(kind=kind_phys) GMF,    PI,    ONPG, CTLA, VTRM, VTPEXP      &
     &,    RPART,  QRMIN,  AA1,    BB1,   CC1,  DD1                     &
     &,    WC2MIN, WCMIN,  WCBASE, F2,    F3,   F5, GMF1, GMF5          &
     &,    QRAF,   QRBF,   del_tla               
!    &,    sialf
!
      parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=0.0)
!     parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=1.0)
!     parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=0.5)
!     PARAMETER (AA1=1.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
!     PARAMETER (AA1=2.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
      PARAMETER (AA1=1.0, BB1=1.0, CC1=1.0, DD1=1.0, F3=CC1,  F5=1.0)
      parameter (QRMIN=1.0E-6, WC2MIN=0.01, GMF1=GMF/AA1, GMF5=GMF/F5)
!     parameter (QRMIN=1.0E-6, WC2MIN=1.00, GMF1=GMF/AA1, GMF5=GMF/F5)
!     parameter (sialf=0.5)
!
      PARAMETER (PI=3.1415926535897931, PIINV=1.0/PI)
      INTEGER ITR, ITRMU, ITRMD, KTPD, ITRMIN, ITRMND
!     PARAMETER (ITRMU=25, ITRMD=25, ITRMIN=7)
      PARAMETER (ITRMU=25, ITRMD=25, ITRMIN=12, ITRMND=12)
!     PARAMETER (ITRMU=25, ITRMD=25, ITRMIN=12)
!     PARAMETER (ITRMU=14, ITRMD=18, ITRMIN=7)
!     PARAMETER (ITRMU=10, ITRMD=10, ITRMIN=5)
      real(kind=kind_phys) QRP(KD:K+1), WVL(KD:K+1), AL2
      real(kind=kind_phys) WVLO(KD:K+1)
!
      real(kind=kind_phys) RNF(KD:K),   ETD(KD:K+1), WCB(KD:K)          &
     &,                    HOD(KD:K+1), QOD(KD:K+1), EVP(KD:K)          &
     &,                    ROR(KD:K+1), STLT(KD:K)                      &
     &,                    GHD(KD:K),   GSD(KD:K),   CLDFRD(KD:K)       &
     &,                    RNT,        RNB                              &
     &,                    ERRQ,       RNTP
      INTEGER IDW, IDH, IDN(K), idnm
      real(kind=kind_phys) ELM(K)
!     real(kind=kind_phys) EM(K*K), ELM(K)
      real(kind=kind_phys) EDZ, DDZ, CE, QHS, FAC, FACG, ASIN,          &
     &                     RSUM1, RSUM2, RSUM3, CEE
      LOGICAL DDFT, UPDRET, DDLGK
!
      real(kind=kind_phys) AA(KD:K,KD:K+1), QW(KD:K,KD:K)               &
     &,                    BUD(KD:K), VT(2), VRW(2), TRW(2)             &
     &,                    GQW(KD:K)                                    &
     &,                    QA(3),     WA(3),    DOF, DOFW               &
     &,                    QRPI(KD:K), QRPS(KD:K)
!    &,                    GQW(KD:K), WCB(KD:K)

!***********************************************************************

      real(kind=kind_phys) QRPF, VTPF
      logical lprnt

!     if(lprnt) print *,' K=',K,' KD=',KD,' In Downdrft'

      KD1    = KD + 1
      KP1    = K  + 1
      KM1    = K  - 1
      KB1    = KBL - 1
!
!     VTP    = 36.34*SQRT(1.2)* (0.001)**0.1364
      VTPEXP = -0.3636
!     PIINV  = 1.0 / PI
      PICON  = PI * ONEBG * 0.5
!
!     Compute Rain Water Budget of the Updraft (Cheng and Arakawa, 1997)
!
      CLDFRD = 0.0
      RNTP   = 0.0
      DOF    = 0.0
      ERRQ   = 10.0
      RNB    = 0.0
      RNT    = 0.0
      TX2    = PRL(KBL)
!
      TX1      = (PRL(KD) + PRL(KD1)) * 0.5
      ROR(KD)  = CMPOR*TX1 / (TOL(KD)*(1.0+NU*QOL(KD)))
!     GMS(KD)  = VTP * ROR(KD) ** VTPEXP
      GMS(KD)  = VTP * VTPF(ROR(KD))
!
      QRP(KD)  = QRMIN
!
      TEM      = TOL(K) * (1.0 + NU * QOL(K))
      ROR(K+1) = 0.5 * CMPOR * (PRL(K+1)+PRL(K)) / TEM
      GMS(K+1) = VTP * VTPF(ROR(K+1))
      QRP(K+1) = QRMIN
!
      kk = kbl
      DO L=KD1,K
        TEM = 0.5 * (TOL(L)+TOL(L-1))                                   &
     &      * (1.0 + (0.5*NU) * (QOL(L)+QOL(L-1)))
        ROR(L) = CMPOR * PRL(L) / TEM
!       GMS(L) = VTP * ROR(L) ** VTPEXP
        GMS(L) = VTP * VTPF(ROR(L))
        QRP(L) = QRMIN
        if (buy(l) .le. 0.0 .and. kk .eq. KBL) then
          kk = l
        endif
      ENDDO
      if (kk .ne. kbl) then
        do l=kk,kbl
          buy(l) = 0.9 * buy(l-1)
        enddo
      endif
!
      do l=kd,k
        qrpi(l) = buy(l)
      enddo
      do l=kd1,kb1
        buy(l) = 0.25 * (qrpi(l-1)+qrpi(l)+qrpi(l)+qrpi(l+1))
      enddo
      
!
!     CALL ANGRAD(TX1, ALM, STLA, CTL2, AL2, PI, TLA, TX2, WFN, TX3)
      tx1 = 1000.0 + tx1 - prl(k+1)
      CALL ANGRAD(TX1, ALM,  AL2, TLA, TX2, WFN, TX3)
!
!    Following Ucla approach for rain profile
!
      F2      = 2.0*BB1*ONEBG/(PI*0.2)
      WCMIN   = SQRT(WC2MIN)
      WCBASE  = WCMIN
!
!     del_tla = TLA * 0.2
!     del_tla = TLA * 0.25
      del_tla = TLA * 0.3
      TLA     = TLA - DEL_TLA
!
      DO L=KD,K
        RNF(L)   = 0.0
        RNS(L)   = 0.0
        WVL(L)   = 0.0
        STLT(L)  = 0.0
        GQW(L)   = 0.0
        QRP(L)   = QRMIN
        DO N=KD,K
          QW(N,L) = 0.0
        ENDDO
      ENDDO
!
!-----QW(N,L) = D(W(N)*W(N))/DQR(L)
!
      KK = KBL
      QW(KD,KD)  = -QRB(KD)  * GMF1
      GHD(KD)    = ETA(KD)   * ETA(KD)
      GQW(KD)    = QW(KD,KD) * GHD(KD)
      GSD(KD)    = ETAI(KD)  * ETAI(KD)
!
      GQW(KK)    = - QRB(KK-1) * (GMF1+GMF1)
!
      WCB(KK)    = WCBASE * WCBASE

      TX1        = WCB(KK)
      GSD(KK)    = 1.0
      GHD(KK)    = 1.0
!
      TEM        = GMF1 + GMF1
      DO L=KB1,KD1,-1
         GHD(L)  = ETA(L)  * ETA(L)
         GSD(L)  = ETAI(L) * ETAI(L)
         GQW(L)  = - GHD(L) * (QRB(L-1)+QRT(L)) * TEM
         QW(L,L) = - QRT(L) * TEM
!
         st1     = 0.5 * (eta(l) + eta(l+1))
         TX1     = TX1 + BUY(L) * TEM * (qrb(l)+qrt(l)) * st1 * st1
         WCB(L)  = TX1 * GSD(L)
      ENDDO
!
      TEM1        = (QRB(KD) + QRT(KD1) + QRT(KD1)) * GMF1
      GQW(KD1)    = - GHD(KD1) * TEM1
      QW(KD1,KD1) = - QRT(KD1) * TEM
      st1         = 0.5 * (eta(kd) + eta(kd1))
      WCB(KD)     = (TX1 + BUY(KD)*TEM*qrb(kd)*st1*st1) * GSD(KD)
!
      DO L=KD1,KBL
        DO N=KD,L-1
           QW(N,L) = GQW(L) * GSD(N)
        ENDDO
      ENDDO
      QW(KBL,KBL) = 0.0
!
      do ntla=1,numtla
!
!       if (errq .lt. 1.0 .or. tla .gt. 45.0) cycle
        if (errq .lt. 0.1 .or. tla .gt. 45.0) cycle
!
        tla = tla + del_tla
        STLA = SIN(TLA*PI/180.0)
        CTL2 = 1.0 - STLA * STLA
!
!       if (lprnt) print *,' tla=',tla,' al2=',al2,' ptop='
!    &,0.5*(prl(kd)+prl(kd1)),' ntla=',ntla,' f2=',f2,' stla=',stla
!       if (lprnt) print *,' buy=',(buy(l),l=kd,kbl)
!
        STLA = F2     * STLA * AL2
        CTL2 = DD1    * CTL2
        CTL3 = 0.1364 * CTL2
!
        DO L=KD,K
          RNF(L)   = 0.0
          WVL(L)   = 0.0
          STLT(L)  = 0.0
          QRP(L)   = QRMIN
        ENDDO
        WVL(KBL)   = WCBASE
        STLT(KBL)  = 1.0 / WCBASE
!
        DO L=KD,K+1
          DO N=KD,K
            AA(N,L) = 0.0
          ENDDO
        ENDDO
!
        SKPUP = .FALSE.
!
        DO ITR=1,ITRMU               ! Rain Profile Iteration starts!
          IF (.NOT. SKPUP) THEN
             wvlo = wvl
!
!-----CALCULATING THE VERTICAL VELOCITY
!
            TX1      = 0.0
            QRPI(KBL) = 1.0 / QRP(KBL)
            DO L=KB1,KD,-1
              TX1     = TX1    + QRP(L+1) * GQW(L+1)
              ST1     = WCB(L) + QW(L,L)  * QRP(L)                      &
     &                         + TX1      * GSD(L)
              if (st1 .gt. wc2min) then
!               WVL(L)  = SQRT(ST1)
                WVL(L)  = 0.5 * (SQRT(ST1) + WVL(L))
!               if (itr .eq. 1) wvl(l) = wvl(l) * 0.25
              else

!       if (lprnt)  print *,' l=',l,' st1=',st1,' wcb=',wcb(l),' qw='
!    &,qw(l,l),' qrp=',qrp(l),' tx1=',tx1,' gsd=',gsd(l),' ite=',itr

!               wvl(l) = 0.5*(wcmin+wvl(l))
                wvl(l) = 0.5*(wvl(l) + wvl(l+1))
                qrp(l) = 0.5*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l) &
     &                      + qrp(l))
!!              wvl(l) = 0.5 * (wvl(l) + wvl(l+1))
              endif
!             wvl(l)  = 0.5 * (wvl(l) + wvlo(l))
!             WVL(L)  = SQRT(MAX(ST1,WC2MIN))
              wvl(l)  = max(wvl(l), wcbase)
              STLT(L) = 1.0 / WVL(L)
              QRPI(L) = 1.0 / QRP(L)
            ENDDO
!
!       if (lprnt) then
!         print *,' ITR=',ITR,' ITRMU=',ITRMU
!         print *,' WVL=',(WVL(L),L=KD,KBL)
!         print *,' qrp=',(qrp(L),L=KD,KBL)
!         print *,' qrpi=',(qrpi(L),L=KD,KBL)
!         print *,' rnf=',(rnf(L),L=KD,KBL)
!       endif
!
!-----CALCULATING TRW, VRW AND OF
!
!           VT(1)   = GMS(KD) * QRP(KD)**0.1364
            VT(1)   = GMS(KD) * QRPF(QRP(KD))
            TRW(1)  = ETA(KD) * QRP(KD) * STLT(KD)
            TX6     = TRW(1) * VT(1)
            VRW(1)  = F3*WVL(KD) - CTL2*VT(1)
            BUD(KD) = STLA * TX6 * QRB(KD) * 0.5
            RNF(KD) = BUD(KD)
            DOF     = 1.1364 * BUD(KD) * QRPI(KD)
            DOFW    = -BUD(KD) * STLT(KD)
!
            RNT     = TRW(1) * VRW(1)
            TX2     = 0.0
            TX4     = 0.0
            RNB     = RNT
            TX1     = 0.5
            TX8     = 0.0
!
            IF (RNT .GE. 0.0) THEN
              TX3 = (RNT-CTL3*TX6) * QRPI(KD)
              TX5 = CTL2 * TX6 * STLT(KD)
            ELSE
              TX3 = 0.0
              TX5 = 0.0
              RNT = 0.0
              RNB = 0.0
            ENDIF
!
            DO L=KD1,KB1
              KTEM    = MAX(L-2, KD)
              LL      = L - 1
! 
!             VT(2)   = GMS(L) * QRP(L)**0.1364
              VT(2)   = GMS(L) * QRPF(QRP(L))
              TRW(2)  = ETA(L) * QRP(L) * STLT(L)
              VRW(2)  = F3*WVL(L) - CTL2*VT(2)
              QQQ     = STLA * TRW(2) * VT(2)
              ST1     = TX1  * QRB(LL)
              BUD(L)  = QQQ * (ST1 + QRT(L))
!
              QA(2)   = DOF
              WA(2)   = DOFW
              DOF     = 1.1364 * BUD(L) * QRPI(L)
              DOFW    = -BUD(L) * STLT(L)
!
              RNF(LL) = RNF(LL) + QQQ * ST1
              RNF(L)  =           QQQ * QRT(L)
!
              TEM3    = VRW(1) + VRW(2)
              TEM4    = TRW(1) + TRW(2)
!
              TX6     = .25 * TEM3 * TEM4
              TEM4    = TEM4 * CTL3
!
!-----BY QR ABOVE
!
!             TEM1    = .25*(TRW(1)*TEM3 - TEM4*VT(1))*TX7
              TEM1    = .25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
              ST1     = .25*(TRW(1)*(CTL2*VT(1)-VRW(2))                 &
     &                     * STLT(LL) + F3*TRW(2))
!-----BY QR BELOW
              TEM2    = .25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
              ST2     = .25*(TRW(2)*(CTL2*VT(2)-VRW(1))                 &
     &                     * STLT(L)  + F3*TRW(1))
!
!      From top to  the KBL-2 layer
!
              QA(1)   = TX2
              QA(2)   = QA(2) + TX3 - TEM1
              QA(3)   = -TEM2
!
              WA(1)   = TX4
              WA(2)   = WA(2) + TX5 - ST1
              WA(3)   = -ST2
!
              TX2     = TEM1
              TX3     = TEM2
              TX4     = ST1
              TX5     = ST2
!
              VT(1)   = VT(2)
              TRW(1)  = TRW(2)
              VRW(1)  = VRW(2)
!
              IF (WVL(KTEM) .EQ. WCMIN) WA(1) = 0.0
              IF (WVL(LL)   .EQ. WCMIN) WA(2) = 0.0
              IF (WVL(L)    .EQ. WCMIN) WA(3) = 0.0
              DO N=KTEM,KBL
                AA(LL,N) = (WA(1)*QW(KTEM,N) * STLT(KTEM)               &
     &                   +  WA(2)*QW(LL,N)   * STLT(LL)                 &
     &                   +  WA(3)*QW(L,N)    * STLT(L) ) * 0.5
              ENDDO
              AA(LL,KTEM) = AA(LL,KTEM) + QA(1)
              AA(LL,LL)   = AA(LL,LL)   + QA(2)
              AA(LL,L)    = AA(LL,L)    + QA(3)
              BUD(LL)     = (TX8 + RNN(LL)) * 0.5                       &
     &                      - RNB + TX6 - BUD(LL)
              AA(LL,KBL+1) = BUD(LL)
              RNB = TX6
              TX1 = 1.0
              TX8 = RNN(LL)
            ENDDO
            L  = KBL
            LL = L - 1
!           VT(2)   = GMS(L) * QRP(L)**0.1364
            VT(2)   = GMS(L) * QRPF(QRP(L))
            TRW(2)  = ETA(L) * QRP(L) * STLT(L)
            VRW(2)  = F3*WVL(L) - CTL2*VT(2)
            ST1     = STLA * TRW(2) * VT(2) * QRB(LL)
            BUD(L)  = ST1

            QA(2)   = DOF
            WA(2)   = DOFW
            DOF     = 1.1364 * BUD(L) * QRPI(L)
            DOFW    = -BUD(L) * STLT(L)
!
            RNF(LL) = RNF(LL) + ST1
!
            TEM3    = VRW(1) + VRW(2)
            TEM4    = TRW(1) + TRW(2)
!
            TX6     = .25 * TEM3 * TEM4
            TEM4    = TEM4 * CTL3
!
!-----BY QR ABOVE
!
            TEM1    = .25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
            ST1     = .25*(TRW(1)*(CTL2*VT(1)-VRW(2))                   &
     &                  * STLT(LL) + F3*TRW(2))
!-----BY QR BELOW
            TEM2    = .25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
            ST2     = .25*(TRW(2)*(CTL2*VT(2)-VRW(1))                   &
     &                   * STLT(L)  + F3*TRW(1))
!
!      For the layer next to the top of the boundary layer
!
            QA(1)   = TX2
            QA(2)   = QA(2) + TX3 - TEM1
            QA(3)   = -TEM2
!
            WA(1)   = TX4
            WA(2)   = WA(2) + TX5 - ST1
            WA(3)   = -ST2
!
            TX2     = TEM1
            TX3     = TEM2
            TX4     = ST1
            TX5     = ST2
!
            IDW     = MAX(L-2, KD)
!
            IF (WVL(IDW) .EQ. WCMIN) WA(1) = 0.0
            IF (WVL(LL)  .EQ. WCMIN) WA(2) = 0.0
            IF (WVL(L)   .EQ. WCMIN) WA(3) = 0.0
!
            KK = IDW
            DO N=KK,L
              AA(LL,N) = (WA(1)*QW(KK,N) * STLT(KK)                     &
     &                 +  WA(2)*QW(LL,N) * STLT(LL)                     &
     &                 +  WA(3)*QW(L,N)  * STLT(L) ) * 0.5

            ENDDO
!
            AA(LL,IDW) = AA(LL,IDW) + QA(1)
            AA(LL,LL)  = AA(LL,LL)  + QA(2)
            AA(LL,L)   = AA(LL,L)   + QA(3)
            BUD(LL)    = (TX8+RNN(LL)) * 0.5 - RNB + TX6 - BUD(LL)
!
            AA(LL,L+1) = BUD(LL)
!
            RNB        = TRW(2) * VRW(2)
!
!      For the top of the boundary layer
!
            IF (RNB .LT. 0.0) THEN
               KK    = KBL
               TEM   = VT(2) * TRW(2)
               QA(2) = (RNB - CTL3*TEM) * QRPI(KK)
               WA(2) = CTL2 * TEM * STLT(KK)
            ELSE
               RNB   = 0.0
               QA(2) = 0.0
               WA(2) = 0.0
            ENDIF
!
            QA(1) = TX2
            QA(2) = DOF + TX3 - QA(2)
            QA(3) = 0.0
!
            WA(1) = TX4
            WA(2) = DOFW + TX5 - WA(2)
            WA(3) = 0.0
!
            KK = KBL
            IF (WVL(KK-1) .EQ. WCMIN) WA(1) = 0.0
            IF (WVL(KK)   .EQ. WCMIN) WA(2) = 0.0
!
            DO II=1,2
               N = KK + II - 2
               AA(KK,N) = (WA(1)*QW(KK-1,N) * STLT(KK-1)                &
     &                  +  WA(2)*QW(KK,N)   * STLT(KK)) * 0.5
            ENDDO
            FAC = 0.5
            LL  = KBL
            L   = LL + 1
            LM1 = LL - 1
            AA(LL,LM1)  = AA(LL,LM1) + QA(1)
            AA(LL,LL)   = AA(LL,LL)  + QA(2)
            BUD(LL)     = 0.5*RNN(LM1) - TX6 + RNB - BUD(LL)
            AA(LL,LL+1) = BUD(LL)
!
!-----SOLVING THE BUDGET EQUATIONS FOR DQR
!
            DO L=KD1,KBL
              LM1  = L - 1
              cnvflg = ABS(AA(LM1,LM1)) .LT. ABS(AA(L,LM1))
              DO  N=LM1,KBL+1
                 IF (cnvflg) THEN
                    TX1       = AA(LM1,N)
                    AA(LM1,N) = AA(L,N)
                    AA(L,N)   = TX1
                 ENDIF
              ENDDO
              TX1 = AA(L,LM1) / AA(LM1,LM1)
              DO  N=L,KBL+1
                AA(L,N) = AA(L,N) - TX1 * AA(LM1,N)
              ENDDO
            ENDDO     
!
!-----BACK SUBSTITUTION AND CHECK IF THE SOLUTION CONVERGES
!
            KK = KBL
            KK1 = KK + 1
            AA(KK,KK1) = AA(KK,KK1) / AA(KK,KK)      !   Qr correction !
            TX2        = ABS(AA(KK,KK1)) * QRPI(KK)  !   Error Measure !
!     if (lprnt) print *,' tx2a=',tx2,' aa1=',aa(kk,kk1)
!    &,' qrpi=',qrpi(kk)
!
            KK = KBL + 1
            DO L=KB1,KD,-1
               LP1   = L + 1
               TX1  = 0.0
               DO N=LP1,KBL
                 TX1  = TX1 + AA(L,N) * AA(N,KK)
               ENDDO
               AA(L,KK) = (AA(L,KK) - TX1) / AA(L,L)       ! Qr correction !
               TX2      = MAX(TX2, ABS(AA(L,KK))*QRPI(L))  ! Error Measure !

!     if (lprnt) print *,' tx2b=',tx2,' aa1=',aa(l,kk)
!    &,' qrpi=',qrpi(l),' L=',L

            ENDDO
!
!           tem = 0.5
            if (tx2 .gt. 1.0 .and. abs(errq-tx2) .gt. 0.1) then
              tem = 0.5
!!          elseif (tx2 .lt. 0.1) then
!!            tem = 1.2
            else
              tem = 1.0
            endif
!
            DO L=KD,KBL
!              QRP(L) = MAX(QRP(L)+AA(L,KBL+1), QRMIN)
               QRP(L) = MAX(QRP(L)+AA(L,KBL+1)*tem, QRMIN)
            ENDDO
!
!       if (lprnt) print *,' itr=',itr,' tx2=',tx2

            IF (ITR .LT. ITRMIN) THEN
               TEM = ABS(ERRQ-TX2) 
               IF (TEM .GE. ERRMI2 .AND. TX2 .GE. ERRMIN) THEN 
                 ERRQ  = TX2                              ! Further iteration !
               ELSE 
                 SKPUP = .TRUE.                           ! Converges      !
                 ERRQ  = 0.0                              ! Rain profile exists!
!     if (lprnt) print *,' here1',' tem=',tem,' tx2=',tx2,' errmi2=',
!    *errmi2,' errmin=',errmin
               ENDIF 
            ELSE
               TEM = ERRQ - TX2
!              IF (TEM .LT. ZERO .AND. ERRQ .GT. 0.1) THEN
               IF (TEM .LT. ZERO .AND. ERRQ .GT. 0.5) THEN
!              IF (TEM .LT. ZERO .and.                                    &
!    &            (ntla .lt. numtla .or. ERRQ .gt. 0.5)) THEN
!     if (lprnt) print *,' tx2=',tx2,' errq=',errq,' tem=',tem
                 SKPUP = .TRUE.                           ! No convergence !
                 ERRQ = 10.0                              ! No rain profile!
!!!!           ELSEIF (ABS(TEM).LT.ERRMI2 .OR. TX2.LT.ERRMIN) THEN
               ELSEIF (TX2.LT.ERRMIN) THEN
                 SKPUP = .TRUE.                           ! Converges      !
                 ERRQ = 0.0                               ! Rain profile exists!
!     if (lprnt) print *,' here2'
               elseif (tem .lt. zero .and. errq .lt. 0.1) then
                 skpup = .true.
!                if (ntla .eq. numtla .or. tem .gt. -0.003) then
                   errq  = 0.0
!                else
!                  errq = 10.0
!                endif
               ELSE
                 ERRQ = TX2                               ! Further iteration !
!     if (lprnt) print *,' itr=',itr,' errq=',errq
!              if (itr .eq. itrmu .and. ERRQ .GT. ERRMIN*10             &
!    &            .and. ntla .eq. 1) ERRQ = 10.0 
               ENDIF
            ENDIF
!
!         if (lprnt) print *,' ERRQ=',ERRQ

          ENDIF                                           ! SKPUP  ENDIF!
!
        ENDDO                                          ! End of the ITR Loop!!
!
!     if(lprnt) then
!       print *,' QRP=',(QRP(L),L=KD,KBL)
!       print *,'RNF=',(RNF(L),L=KD,KBL),' RNT=',RNT,' RNB=',RNB
!    &,' errq=',errq
!     endif
!
        IF (ERRQ .LT. 0.1) THEN
          DDFT = .TRUE.
          RNB  = - RNB
   !      do l=kd1,kb1-1
   !        if (wvl(l)-wcbase .lt. 1.0E-9) ddft = .false.
   !      enddo
        ELSE
          DDFT = .FALSE.
        ENDIF
!
!     Caution !! Below is an adjustment to rain flux to maintain
!                conservation of precip!
!
        IF (DDFT) THEN
          TX1 = 0.0
          DO L=KD,KB1
            TX1 = TX1 + RNF(L)
          ENDDO
!     if (lprnt) print *,' tx1+rnt+rnb=',tx1+rnt+rnb, ' train=',train
          TX1 = TRAIN / (TX1+RNT+RNB)
          IF (ABS(TX1-1.0) .LT. 0.2) THEN
             RNT = MAX(RNT*TX1,ZERO)
             RNB = RNB * TX1
          ELSE
             DDFT = .FALSE.
             ERRQ = 10.0
          ENDIF
        ENDIF
      enddo                                          ! End of ntla loop
!
      DOF = 0.0
      IF (.NOT. DDFT) RETURN     ! Rain profile did not converge!
!

      DO L=KD,KB1
         RNF(L) = RNF(L) * TX1

      ENDDO
!     if (lprnt) print *,' TRAIN=',TRAIN
!     if (lprnt) print *,' RNF=',RNF
!
!     Adjustment is over
!
!     Downdraft
!
      DO L=KD,K
        WCB(L) = 0.0
      ENDDO
!
      SKPDD = .NOT. DDFT
!
      ERRQ  = 10.0
      IF (.NOT. SKPDD) THEN
!
!     Calculate Downdraft Properties
!

        KK = MAX(KB1,KD1)
        DO L=KK,K
          STLT(L) = STLT(L-1)
        ENDDO
        TEM1 = 1.0 / BB1
!
        DO L=KD,K
          IF (L .LE. KBL) THEN
            TEM     = STLA * TEM1
            STLT(L) = ETA(L) * STLT(L) * TEM / ROR(L)
          ELSE
            STLT(L) = 0.0
          ENDIF
        ENDDO
!       if (lprnt) print *,' STLT=',stlt

        rsum1 = 0.0
        rsum2 = 0.0

!
        IDN      = 99
        DO L=KD,K+1
          ETD(L)  = 0.0
          WVL(L)  = 0.0
!         QRP(L)  = 0.0
        ENDDO
        DO L=KD,K
          EVP(L)   = 0.0
          BUY(L)   = 0.0
          QRP(L+1) = 0.0
        ENDDO
        HOD(KD)  = HOL(KD)
        QOD(KD)  = QOL(KD)
        TX1      = 0.0                               ! sigma at the top
!!!     TX1      = STLT(KD)*QRB(KD)*ONE              ! sigma at the top
!       TX1      = MIN(STLT(KD)*QRB(KD)*ONE, ONE)    ! sigma at the top
!       TX1      = MIN(STLT(KD)*QRB(KD)*0.5, ONE)    ! sigma at the top
        RNTP     = 0.0
        TX5      = TX1
        QA(1)    = 0.0
!     if(lprnt) print *,' stlt=',stlt(kd),' qrb=',qrb(kd)
!    *,' tx1=',tx1,' ror=',ror(kd),' gms=',gms(kd),' rpart=',rpart
!    *,' rnt=',rnt
!
!       Here we assume RPART of detrained rain RNT goes to Pd
!
        IF (RNT .GT. 0.0) THEN
          if (TX1 .gt. 0.0) THEN
            QRP(KD) = (RPART*RNT / (ROR(KD)*TX1*GMS(KD)))               &
     &                                          ** (1.0/1.1364)
           else
             tx1 = RPART*RNT / (ROR(KD)*GMS(KD)*QRP(KD)**1.1364)
           endif
            RNTP    = (1.0 - RPART) * RNT
            BUY(KD) = - ROR(KD) * TX1 * QRP(KD)
        ELSE
          QRP(KD) = 0.0
        ENDIF
!
!     L-loop for the downdraft iteration from KD1 to K+1 (bottom surface)
!
!     BUD(KD) = ROR(KD)
      idnm = 1
      DO L=KD1,K+1

          QA(1) = 0.0
          ddlgk = idn(idnm) .eq. 99
          if (.not. ddlgk) cycle
          IF (L .LE. K) THEN
            ST1   = 1.0 - ALFIND(L)
            WA(1) = ALFIND(L)*HOL(L-1) + ST1*HOL(L)
            WA(2) = ALFIND(L)*QOL(L-1) + ST1*QOL(L)
            WA(3) = ALFIND(L)*TOL(L-1) + ST1*TOL(L)
            QA(2) = ALFIND(L)*HST(L-1) + ST1*HST(L)
            QA(3) = ALFIND(L)*QST(L-1) + ST1*QST(L)
          ELSE
            WA(1) = HOL(K)
            WA(2) = QOL(K)
            WA(3) = TOL(K)
            QA(2) = HST(K)
            QA(3) = QST(K)
          ENDIF
!
          FAC = 2.0
          IF (L .EQ. KD1) FAC = 1.0

          FACG    = FAC * 0.5 * GMF5     !  12/17/97
!
!         DDLGK   =  IDN(idnm) .EQ. 99
          BUD(KD) = ROR(L)

!         IF (DDLGK) THEN
            TX1    = TX5
            WVL(L) = MAX(WVL(L-1),ONE_M1)

            QRP(L) = MAX(QRP(L-1),QRP(L))
!
!           VT(1)  = GMS(L-1) * QRP(L-1) ** 0.1364
            VT(1)  = GMS(L-1) * QRPF(QRP(L-1))
            RNT    = ROR(L-1) * (WVL(L-1)+VT(1))*QRP(L-1)
!     if(lprnt) print *,' l=',l,' qa=',qa(1), ' tx1RNT=',RNT*tx1,
!    *' wvl=',wvl(l-1)
!    *,' qrp=',qrp(l-1),' tx5=',tx5,' tx1=',tx1,' rnt=',rnt

!

!           TEM    = MAX(ALM, 2.5E-4) * MAX(ETA(L), 1.0)
            TEM    = MAX(ALM,ONE_M6) * MAX(ETA(L), ONE)
!           TEM    = MAX(ALM, 1.0E-5) * MAX(ETA(L), 1.0)
            TRW(1) = PICON*TEM*(QRB(L-1)+QRT(L-1))
            TRW(2) = 1.0 / TRW(1)
!
            VRW(1) = 0.5 * (GAM(L-1) + GAM(L))
            VRW(2) = 1.0 / (VRW(1) + VRW(1))
!
            TX4    =  (QRT(L-1)+QRB(L-1))*(ONEBG*FAC*500.00*EKNOB)
!
            DOFW   = 1.0 / (WA(3) * (1.0 + NU*WA(2)))      !  1.0 / TVbar!
!
            ETD(L) = ETD(L-1)
            HOD(L) = HOD(L-1)
            QOD(L) = QOD(L-1)
!
            ERRQ   = 10.0

!
            IF (L .LE. KBL) THEN
              TX3 = STLT(L-1) * QRT(L-1) * (0.5*FAC)
              TX8 = STLT(L)   * QRB(L-1) * (0.5*FAC)
              TX9 = TX8 + TX3
            ELSE
              TX3 = 0.0
              TX8 = 0.0
              TX9 = 0.0
            ENDIF
!
            TEM  = WVL(L-1) + VT(1)
            IF (TEM .GT. 0.0) THEN
              TEM1 = 1.0 / (TEM*ROR(L-1))
              TX3 = VT(1) * TEM1 * ROR(L-1) * TX3
              TX6 = TX1 * TEM1
            ELSE
              TX6 = 1.0
            ENDIF
!         ENDIF
!
          IF (L .EQ. KD1) THEN
            IF (RNT .GT. 0.0) THEN
              TEM    = MAX(QRP(L-1),QRP(L))
              WVL(L) = TX1 * TEM * QRB(L-1)*(FACG*5.0)
            ENDIF
            WVL(L) = MAX(ONE_M2, WVL(L))
            TRW(1) = TRW(1) * 0.5
            TRW(2) = TRW(2) + TRW(2)
          ELSE
            IF (DDLGK) EVP(L-1) = EVP(L-2)
          ENDIF
!
!       No downdraft above level IDH
!

          IF (L .LT. IDH) THEN

            ETD(L)   = 0.0
            HOD(L)   = WA(1)
            QOD(L)   = WA(2)
            EVP(L-1) = 0.0
            WVL(L)   = 0.0
            QRP(L)   = 0.0
            BUY(L)   = 0.0
            TX5      = TX9
            ERRQ     = 0.0
            RNTP     = RNTP + RNT * TX1
            RNT      = 0.0
            WCB(L-1) = 0.0
          ENDIF
!         BUD(KD) = ROR(L)
!
!       Iteration loop for a given level L begins
!
!         if (lprnt) print *,' tx8=',tx8,' tx9=',tx9,' tx5=',tx5
!    &,                      ' tx1=',tx1
          DO ITR=1,ITRMD
!
!           cnvflg =  DDLGK .AND. (ERRQ .GT. ERRMIN)
            cnvflg =  ERRQ .GT. ERRMIN
            IF (cnvflg) THEN
!
!             VT(1)  = GMS(L) * QRP(L) ** 0.1364
              VT(1)  = GMS(L) * QRPF(QRP(L))
              TEM    =  WVL(L) + VT(1)
!
              IF (TEM .GT. 0.0) THEN
                ST1    = ROR(L) * TEM * QRP(L) + RNT
                IF (ST1 .NE. 0.0) ST1 = 2.0 * EVP(L-1) / ST1
                TEM1   = 1.0 / (TEM*ROR(L))
                TEM2   = VT(1) * TEM1 * ROR(L) * TX8
              ELSE
                TEM1   = 0.0
                TEM2   = TX8
                ST1    = 0.0
              ENDIF
!     if (lprnt) print *,' st1=',st1,' tem=',tem,' ror=',ror(l)
!    &,' qrp=',qrp(l),' rnt=',rnt,' ror1=',ror(l-1),' wvl=',wvl(l)
!    &,' wvl1=',wvl(l-1),' tem2=',tem2,' vt=',vt(1),' tx3=',tx3
!
              st2 = tx5
              TEM = ROR(L)*WVL(L) - ROR(L-1)*WVL(L-1)
              if (tem .gt. 0.0) then
                TX5 = (TX1 - ST1 + TEM2 + TX3)/(1.0+tem*tem1)
              else
                TX5 = TX1 - tem*tx6 - ST1 + TEM2 + TX3
              endif
              TX5   = MAX(TX5,ZERO)
              tx5 = 0.5 * (tx5 + st2)
!
!             qqq = 1.0 + tem * tem1 * (1.0 - sialf)
!
!             if (qqq .gt. 0.0) then
!               TX5   = (TX1 - sialf*tem*tx6 - ST1 + TEM2 + TX3) / qqq
!             else
!               TX5   = (TX1 - tem*tx6 - ST1 + TEM2 + TX3)
!             endif
!
!     if(lprnt) print *,' tx51=',tx5,' tx1=',tx1,' st1=',st1,' tem2='
!     if(tx5 .le. 0.0 .and. l .gt. kd+2)
!    * print *,' tx51=',tx5,' tx1=',tx1,' st1=',st1,' tem2='
!    *,tem2,' tx3=',tx3,' tem=',tem,' tem1=',tem1,' wvl=',wvl(l-1),
!    &wvl(l),' l=',l,' itr=',itr,' evp=',evp(l-1),' vt=',vt(1)
!    *,' qrp=',qrp(l),' rnt=',rnt,' kd=',kd
!     if (lprnt) print *,' etd=',etd(l),' wvl=',wvl(l)
!    &,' trw=',trw(1),trw(2),' ror=',ror(l),' wa=',wa


!
              TEM1   = ETD(L)
              ETD(L) = ROR(L) * TX5 * MAX(WVL(L),ZERO)
!
              if (etd(l) .gt. 0.0) etd(l) = 0.5 * (etd(l) + tem1)
!

              DEL_ETA = ETD(L) - ETD(L-1)

!               TEM       = DEL_ETA * TRW(2)
!               TEM2      = MAX(MIN(TEM, 1.0), -1.0)
!               IF (ABS(TEM) .GT. 1.0 .AND. ETD(L) .GT. 0.0 ) THEN
!                 DEL_ETA = TEM2 * TRW(1)
!                 ETD(L)  = ETD(L-1) + DEL_ETA
!               ENDIF
!               IF (WVL(L) .GT. 0.0) TX5 = ETD(L) / (ROR(L)*WVL(L))
!
                ERRE  = ETD(L) - TEM1
!
                tem  = max(abs(del_eta), trw(1))
                tem2 = del_eta / tem
                TEM1 = SQRT(MAX((tem+DEL_ETA)*(tem-DEL_ETA),ZERO))
!               TEM1 = SQRT(MAX((TRW(1)+DEL_ETA)*(TRW(1)-DEL_ETA),0.0))

                EDZ  = (0.5 + ASIN(TEM2)*PIINV)*DEL_ETA + TEM1*PIINV

              DDZ   = EDZ - DEL_ETA
              WCB(L-1) = ETD(L) + DDZ
!
              TEM1  = HOD(L)
              IF (DEL_ETA .GT. 0.0) THEN
                QQQ    = 1.0 / (ETD(L) + DDZ)
                HOD(L) = (ETD(L-1)*HOD(L-1) + DEL_ETA*HOL(L-1)          &
     &                                            + DDZ*WA(1)) * QQQ
                QOD(L) = (ETD(L-1)*QOD(L-1) + DEL_ETA*QOL(L-1)          &
     &                                            + DDZ*WA(2)) * QQQ
              ELSEif((ETD(L-1) + EDZ) .gt. 0.0) then
                QQQ    = 1.0 / (ETD(L-1) + EDZ)
                HOD(L) = (ETD(L-1)*HOD(L-1) + EDZ*WA(1)) * QQQ
                QOD(L) = (ETD(L-1)*QOD(L-1) + EDZ*WA(2)) * QQQ
              ENDIF
              ERRH  = HOD(L) - TEM1
              ERRQ  = ABS(ERRH/HOD(L))  + ABS(ERRE/MAX(ETD(L),ONE_M5))
!     if (lprnt) print *,' ERRQP=',errq,' errh=',errh,' hod=',hod(l)
!    &,' erre=',erre,' etd=',etd(l),' del_eta=',del_eta
              DOF   = DDZ
              VT(2) = QQQ

!
              DDZ  = DOF
              TEM4 = QOD(L)
              TEM1 = VRW(1)
!
              QHS  = QA(3) + 0.5 * (GAF(L-1)+GAF(L))                    &
     &                           * (HOD(L)-QA(2))
!
!                                           First iteration       !
!
              ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
              TEM2 = ROR(L) * QRP(L)
              CALL QRABF(TEM2,QRAF,QRBF)
              TEM6 = TX5 * (1.6 + 124.9 * QRAF) * QRBF * TX4
!
              CE   = TEM6 * ST2 / ((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
!
              TEM2   = - ((1.0+TEM1)*(QHS+CE) + TEM1*QOD(L))
              TEM3   = (1.0 + TEM1) * QHS * (QOD(L)+CE)
              TEM    = MAX(TEM2*TEM2 - 4.0*TEM1*TEM3,ZERO)
              QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!

!
!                                            second iteration   !
!
              ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
              CE   = TEM6 * ST2 / ((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
!             CEE  = CE * (ETD(L)+DDZ)
!


              TEM2   = - ((1.0+TEM1)*(QHS+CE) + TEM1*tem4)
              TEM3   = (1.0 + TEM1) * QHS * (tem4+CE)
              TEM    = MAX(TEM2*TEM2 - 4.0*TEM1*TEM3,ZERO)
              QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!                                              Evaporation in Layer L-1
!

              EVP(L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)
!                                              Calculate Pd (L+1/2)
              QA(1)    = TX1*RNT + RNF(L-1) - EVP(L-1)
!
!     if(lprnt) print *,' etd=',etd(l),' tx5=',tx5,' rnt=',rnt
!    *,' rnf=',rnf(l-1),' evp=',evp(l-1),' itr=',itr,' L=',L

!
              if (qa(1) .gt. 0.0) then
              IF (ETD(L) .GT. 0.0) THEN
                TEM    = QA(1) / (ETD(L)+ROR(L)*TX5*VT(1))
                QRP(L) = MAX(TEM,ZERO)
              ELSEIF (TX5 .GT. 0.0) THEN
                QRP(L) = (MAX(ZERO,QA(1)/(ROR(L)*TX5*GMS(L))))           &
     &                                          ** (1.0/1.1364)
              ELSE
                QRP(L) = 0.0
              ENDIF
              else
                qrp(l) = 0.5 * qrp(l)
              endif
!                                              Compute Buoyancy
              TEM1   = WA(3)+(HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2)))         &
     &                                                  * (1.0/CP)
!             if (lprnt) print *,' tem1=',tem1,' wa3=',wa(3),' hod='
!    &,hod(l),' wa1=',wa(1),' qod=',qod(l),' wa2=',wa(2),' alhl=',alhl
!    &,' cmpor=',cmpor,' dofw=',dofw,' prl=',prl(l),' qrp=',qrp(l)
              TEM1   = TEM1 * (1.0 + NU*QOD(L))
              ROR(L) = CMPOR * PRL(L) / TEM1
              TEM1   = TEM1 * DOFW
!!!           TEM1   = TEM1 * (1.0 + NU*QOD(L)) * DOFW

              BUY(L) = (TEM1 - 1.0 - QRP(L)) * ROR(L) * TX5
!                                              Compute W (L+1/2)

              TEM1   = WVL(L)
!             IF (ETD(L) .GT. 0.0) THEN
              WVL(L) = VT(2) * (ETD(L-1)*WVL(L-1) - FACG                &
     &                 * (BUY(L-1)*QRT(L-1)+BUY(L)*QRB(L-1)))
!
!             if (lprnt) print *,' wvl=',wvl(l),'vt2=',vt(2),' buy1='
!    &,buy(l-1),' buy=',buy(l),' qrt1=',qrt(l-1),' qrb1=',qrb(l-1)
!    &,' etd1=',etd(l-1),' wvl1=',wvl(l-1)
!             ENDIF
!
              if (wvl(l) .lt. 0.0) then
!               WVL(L) = max(wvl(l), 0.1*tem1)
!               WVL(L) = 0.5*tem1
!               WVL(L) = 0.1*tem1
!               WVL(L) = 0.0
                WVL(L) = 1.0e-10
              else
                WVL(L) = 0.5*(WVL(L)+TEM1)
              endif

!
!             WVL(L) = max(0.5*(WVL(L)+TEM1), 0.0)

              ERRW   = WVL(L) - TEM1
!
              ERRQ   = ERRQ + ABS(ERRW/MAX(WVL(L),ONE_M5))

!     if (lprnt) print *,' errw=',errw,' wvl=',wvl(l)
!     if(lprnt .or. tx5 .eq. 0.0) then
!     if(tx5 .eq. 0.0 .and. l .gt. kbl) then
!        print *,' errq=',errq,' itr=',itr,' l=',l,' wvl=',wvl(l)
!    &,' tx5=',tx5,' idnm=',idnm,' etd1=',etd(l-1),' etd=',etd(l)
!    &,' kbl=',kbl
!     endif
!
!     if(lprnt) print *,' itr=',itr,' itrmnd=',itrmnd,' itrmd=',itrmd
!             IF (ITR .GE. MIN(ITRMIN,ITRMD/2)) THEN
              IF (ITR .GE. MIN(ITRMND,ITRMD/2)) THEN
!     if(lprnt) print *,' itr=',itr,' etd1=',etd(l-1),' errq=',errq
                IF (ETD(L-1) .EQ. 0.0 .AND. ERRQ .GT. 0.2) THEN
!     if(lprnt) print *,' bud=',bud(kd),' wa=',wa(1),wa(2)
                  ROR(L)   = BUD(KD)
                  ETD(L)   = 0.0
                  WVL(L)   = 0.0
                  ERRQ     = 0.0
                  HOD(L)   = WA(1)
                  QOD(L)   = WA(2)
!                 TX5      = TX1 + TX9
                  if (L .le. KBL) then
                    TX5      = TX9
                  else
                    TX5 = (STLT(KB1) * QRT(KB1)                         &
     &                  +  STLT(KBL) * QRB(KB1)) * (0.5*FAC)
                  endif

!     if(lprnt) print *,' tx1=',tx1,' rnt=',rnt,' rnf=',rnf(l-1)
!    *,' evp=',evp(l-1),' l=',l

                  EVP(L-1) = 0.0
                  TEM      = MAX(TX1*RNT+RNF(L-1),ZERO)
                  QA(1)    = TEM - EVP(L-1)
!                 IF (QA(1) .GT. 0.0) THEN

!     if(lprnt) print *,' ror=',ror(l),' tx5=',tx5,' tx1=',tx1
!    *,' tx9=',tx9,' gms=',gms(l),' qa=',qa(1)
!     if(lprnt) call mpi_quit(13)
!     if (tx5 .eq. 0.0 .or. gms(l) .eq. 0.0)
!     if (lprnt) 
!    *  print *,' Atx5=',tx5,' gms=',gms(l),' ror=',ror(l)
!    *,' L=',L,' QA=',QA(1),' tx1=',tx1,' tx9=',tx9
!    *,' kbl=',kbl,' etd1=',etd(l-1),' idnm=',idnm,' idn=',idn(idnm)
!    *,' errq=',errq

                  QRP(L)   = (QA(1) / (ROR(L)*TX5*GMS(L)))              &
     &                                            ** (1.0/1.1364)
!                 endif
                  BUY(L)   = - ROR(L) * TX5 * QRP(L)
                  WCB(L-1) = 0.0
                ENDIF
!
                DEL_ETA = ETD(L) - ETD(L-1)
                IF(DEL_ETA .LT. 0.0 .AND. ERRQ .GT. 0.1) THEN
                  ROR(L)   = BUD(KD)
                  ETD(L)   = 0.0
                  WVL(L)   = 0.0
!!!!!             TX5      = TX1 + TX9
                  CLDFRD(L-1) = TX5
!
                  DEL_ETA  = - ETD(L-1)
                  EDZ      = 0.0
                  DDZ      = -DEL_ETA
                  WCB(L-1) = DDZ
!
                  HOD(L)   = HOD(L-1)
                  QOD(L)   = QOD(L-1)
!
                  TEM4     = QOD(L)
                  TEM1     = VRW(1)
!
                  QHS      = QA(3) + 0.5 * (GAF(L-1)+GAF(L))            &
     &                                   * (HOD(L)-QA(2))

!
!                                           First iteration       !
!
                  ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                  TEM2 = ROR(L) * QRP(L-1)
                  CALL QRABF(TEM2,QRAF,QRBF)
                  TEM6 = TX5 * (1.6 + 124.9 * QRAF) * QRBF * TX4
!
                  CE   = TEM6*ST2/((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
!

                  TEM2   = - ((1.0+TEM1)*(QHS+CE) + TEM1*QOD(L))
                  TEM3   = (1.0 + TEM1) * QHS * (QOD(L)+CE)
                  TEM    = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
                  QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!
!                                            second iteration   !
!
                  ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                  CE   = TEM6*ST2/((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
!                 CEE  = CE * (ETD(L)+DDZ)
!


                  TEM2   = - ((1.0+TEM1)*(QHS+CE) + TEM1*tem4)
                  TEM3   = (1.0 + TEM1) * QHS * (tem4+CE)
                  TEM    = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
                  QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))

!                                              Evaporation in Layer L-1
!
                  EVP(L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)

!                                               Calculate Pd (L+1/2)
!                 RNN(L-1) = TX1*RNT + RNF(L-1) - EVP(L-1)

                  QA(1)    = TX1*RNT + RNF(L-1)
                  EVP(L-1) = min(EVP(L-1), QA(1))
                  QA(1)    = QA(1) - EVP(L-1)
                  qrp(l)   = 0.0

!
!     if (tx5 .eq. 0.0 .or. gms(l) .eq. 0.0)
!     if (lprnt)
!    *  print *,' Btx5=',tx5,' gms=',gms(l),' ror=',ror(l)
!    *,' L=',L,' QA=',QA(1),' tx1=',tx1,' tx9=',tx9
!    *,' kbl=',kbl,' etd1=',etd(l-1),' DEL_ETA=',DEL_ETA
!    &,' evp=',evp(l-1)
!
!                 IF (QA(1) .GT. 0.0) THEN
!!                  RNS(L-1) = QA(1)
!!!                 tx5      = tx9
!                   QRP(L) = (QA(1) / (ROR(L)*TX5*GMS(L)))              &
!    &                                         ** (1.0/1.1364)
!                 endif
!                 ERRQ   = 0.0
!                                              Compute Buoyancy
!                 TEM1   = WA(3)+(HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2)))     &
!    &                                                  * (1.0/CP)
!                 TEM1   = TEM1 * (1.0 + NU*QOD(L)) * DOFW
!                 BUY(L) = (TEM1 - 1.0 - QRP(L)) * ROR(L) * TX5
!
!                 IF (QA(1) .GT. 0.0) RNS(L) = QA(1)

                  IF (L .LE. K) THEN
                     RNS(L) = QA(1)
                     QA(1)  = 0.0
                  ENDIF
                  tx5      = tx9
                  ERRQ     = 0.0
                  QRP(L)   = 0.0
                  BUY(L)   = 0.0
!
                ENDIF
              ENDIF
            ENDIF
!
          ENDDO                ! End of the iteration loop  for a given L!
          IF (L .LE. K) THEN
            IF (ETD(L-1) .EQ. 0.0                                       &
     &         .AND. ERRQ .GT. 0.1 .and. l .le. kbl) THEN
!!!  &         .AND. ERRQ .GT. ERRMIN*10.0 .and. l .le. kbl) THEN
!    &         .AND. ERRQ .GT. ERRMIN*10.0) THEN
               ROR(L)   = BUD(KD)
               HOD(L)   = WA(1)
               QOD(L)   = WA(2)
               TX5      =       TX9     ! Does not make too much difference!
!              TX5      = TX1 + TX9
               EVP(L-1) = 0.0
!              EVP(L-1) = CEE * (1.0 - qod(l)/qa(3))
               QA(1)    = TX1*RNT + RNF(L-1)
               EVP(L-1) = min(EVP(L-1), QA(1))
               QA(1)    = QA(1) - EVP(L-1)

!              QRP(L)   = 0.0
!              if (tx5 .eq. 0.0 .or. gms(l) .eq. 0.0) then
!                print *,' Ctx5=',tx5,' gms=',gms(l),' ror=',ror(l)     &
!    &,          ' L=',L,' QA=',QA(1),' tx1=',tx1,' tx9=',tx9           &
!    &,          ' kbl=',kbl,' etd1=',etd(l-1),' DEL_ETA=',DEL_ETA
!              endif
!              IF (QA(1) .GT. 0.0) THEN

                 QRP(L) = (QA(1) / (ROR(L)*TX5*GMS(L)))                 &
     &                                         ** (1.0/1.1364)
!              ENDIF
               ETD(L)   = 0.0
               WVL(L)   = 0.0
               ST1      = 1.0 - ALFIND(L)

               ERRQ     = 0.0
               BUY(L)   = - ROR(L) * TX5 * QRP(L)
               WCB(L-1) = 0.0
            ENDIF
          ENDIF
!
          LL = MIN(IDN(idnm), K+1)
          IF (ERRQ .LT. 1.0 .AND. L .LE. LL) THEN
            IF (ETD(L-1) .GT. 0.0 .AND. ETD(L) .EQ. 0.0) THEN
             IDN(idnm) = L
             wvl(l)    = 0.0
             if (L .lt. KBL .or. tx5 .gt. 0.0) idnm  = idnm + 1
             errq      = 0.0
            ENDIF
            if (etd(l) .eq. 0.0 .and. l .gt. kbl) then
              idn(idnm) = l
              if (tx5 .gt. 0.0) idnm  = idnm + 1
            endif
          ENDIF

!         if (lprnt) then
!           print *,' ERRQ=',ERRQ,' IDN=',IDN(idnm),' idnm=',idnm
!           print *,' L=',L,' QRP=',QRP(L),' ETD=',ETD(L),' QA=',QA(1)
!    *,' evp=',evp(l-1),' rnf=',rnf(l-1)
!         endif

! 
!     If downdraft properties are not obtainable, (i.e.solution does
!      not converge) , no downdraft is assumed
!
!         IF (ERRQ .GT. ERRMIN*100.0 .AND. IDN(idnm) .EQ. 99)           &
          IF (ERRQ .GT. 0.1 .AND. IDN(idnm) .EQ. 99)                    &
     &                          DDFT = .FALSE.
!
!
          DOF = 0.0
          IF (.NOT. DDFT) RETURN
!
!         if (ddlgk .or. l .le. idn(idnm)) then
!           rsum2 = rsum2 + evp(l-1)
!           print *,' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' qa=',qa(1)&
!    &,   ' evp=',evp(l-1)
!         else
!           rsum1 = rsum1 + rnf(l-1)
!           print *,' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' rnf=',    &
!     &     rnf(l-1)
!         endif

        ENDDO                      ! End of the L Loop of downdraft !

        TX1 = 0.0

        DOF = QA(1)
!
!       print *,' dof=',dof,' rntp=',rntp,' rnb=',rnb
!       print *,' total=',(rsum1+dof+rntp+rnb)

      ENDIF                       ! SKPDD endif
!

      RNN(KD) = RNTP
      TX1     = EVP(KD)
      TX2     = RNTP + RNB + DOF

!     if (lprnt) print *,' tx2=',tx2
      II = IDH
      IF (II .GE. KD1+1) THEN
         RNN(KD)   = RNN(KD) + RNF(KD)
         TX2       = TX2 + RNF(KD)
         RNN(II-1) = 0.0
         TX1       = EVP(II-1)
      ENDIF
!     if (lprnt) print *,' tx2=',tx2,' idnm=',idnm,' idn=',idn(idnm)
      DO L=KD,K
        II = IDH

        IF (L .GT. KD1 .AND. L .LT. II) THEN
          RNN(L-1) = RNF(L-1)
          TX2      = TX2 + RNN(L-1)
        ELSEIF (L .GE. II .AND. L .LT. IDN(idnm)) THEN
          rnn(l)   = rns(l)
          tx2      = tx2 + rnn(l)
          TX1      = TX1 + EVP(L)
        ELSEIF (L .GE. IDN(idnm)) THEN
          ETD(L+1) = 0.0
          HOD(L+1) = 0.0
          QOD(L+1) = 0.0
          EVP(L)   = 0.0
          RNN(L)   = RNF(L) + RNS(L)
          TX2      = TX2    + RNN(L)
        ENDIF
!     if (lprnt) print *,' tx2=',tx2,' L=',L,' rnn=',rnn(l)
      ENDDO
!
!      For Downdraft case the rain is that falls thru the bottom

      L = KBL

      RNN(L)    = RNN(L) + RNB
      CLDFRD(L) = TX5

!
!     Caution !! Below is an adjustment to rain flux to maintain
!                conservation of precip!

!
!     if (lprnt) print *,' train=',train,' tx2=',tx2,' tx1=',tx1

      IF (TX1 .GT. 0.0) THEN
        TX1 = (TRAIN - TX2) / TX1
      ELSE
        TX1 = 0.0
      ENDIF

      DO L=KD,K
        EVP(L) = EVP(L) * TX1
      ENDDO
!
!***********************************************************************
!***********************************************************************

      RETURN
      END

      SUBROUTINE QSATCN(TT,P,Q,DQDT)
!     SUBROUTINE QSATCN(TT,P,Q,DQDT,lprnt)

      USE MACHINE , ONLY : kind_phys
      USE FUNCPHYS , ONLY : fpvs
      USE PHYSCONS, RV => con_RV, CVAP => con_CVAP, CLIQ => con_CLIQ    &
     &,             CSOL => con_CSOL, TTP => con_TTP, HVAP => con_HVAP  &
     &,             HFUS => con_HFUS, EPS => con_eps, EPSM1 => con_epsm1
      implicit none
!
      real(kind=kind_phys) TT, P, Q, DQDT
!
      real(kind=kind_phys) rvi, facw, faci, hsub, tmix, DEN
      real(kind=kind_phys) ZERO,ONE,ONE_M10
      PARAMETER (RVI=1.0/RV)
      PARAMETER (FACW=CVAP-CLIQ, FACI=CVAP-CSOL)
      PARAMETER (HSUB=HVAP+HFUS, tmix=TTP-20.0, DEN=1.0/(TTP-TMIX))
      PARAMETER (ZERO=0.,ONE=1.,ONE_M10=1.E-10)
!     logical lprnt
!
      real(kind=kind_phys) es, d, hlorv, W
!
!     es    = 10.0 * fpvs(tt)                ! fpvs is in centibars!
      es    = 0.01 * fpvs(tt)                ! fpvs is in Pascals!
      D     = 1.0 / max(p+epsm1*es,ONE_M10)
!
      q     = MIN(eps*es*D, ONE)
!
      W     = max(ZERO, min(ONE, (TT - TMIX)*DEN))
      hlorv = ( W      * (HVAP + FACW * (tt-ttp))                       &
     &       + (1.0-W) * (HSUB + FACI * (tt-ttp)) ) * RVI
      dqdt  = p * q * hlorv *  D / (tt*tt)
!
      return
      end

      SUBROUTINE ANGRAD( PRES, ALM,  AL2, TLA, PRB, WFN, UFN)
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : refp, refr, tlac, plac, tlbpl, drdp, almax
      implicit none

      real(kind=kind_phys) PRES                                         &
     &,                    ALM,  AL2,  TLA,  TEM, TEM1                  &
     &,                    PRB,  ACR,  WFN,  UFN
!
      integer i
!
      IF (TLA .LT. 0.0) THEN
          IF (PRES .LE. PLAC(1)) THEN
            TLA = TLAC(1)
          ELSEIF (PRES .LE. PLAC(2)) THEN
            TLA = TLAC(2) + (PRES-PLAC(2))*tlbpl(1)
          ELSEIF (PRES .LE. PLAC(3)) THEN
            TLA = TLAC(3) + (PRES-PLAC(3))*tlbpl(2)
          ELSEIF (PRES .LE. PLAC(4)) THEN
            TLA = TLAC(4) + (PRES-PLAC(4))*tlbpl(3)
          ELSEIF (PRES .LE. PLAC(5)) THEN
            TLA = TLAC(5) + (PRES-PLAC(5))*tlbpl(4)
          ELSEIF (PRES .LE. PLAC(6)) THEN
            TLA = TLAC(6) + (PRES-PLAC(6))*tlbpl(5)
          ELSEIF (PRES .LE. PLAC(7)) THEN
            TLA = TLAC(7) + (PRES-PLAC(7))*tlbpl(6)
          ELSEIF (PRES .LE. PLAC(8)) THEN
            TLA = TLAC(8) + (PRES-PLAC(8))*tlbpl(7)
          ELSE
            TLA = TLAC(8)
          ENDIF
      ENDIF
        IF (PRES .GE. REFP(1)) THEN
          TEM = REFR(1)
        ELSEIF (PRES .GE. REFP(2)) THEN
          TEM = REFR(1) + (PRES-REFP(1)) * drdp(1)
        ELSEIF (PRES .GE. REFP(3)) THEN
          TEM = REFR(2) + (PRES-REFP(2)) * drdp(2)
        ELSEIF (PRES .GE. REFP(4)) THEN
          TEM = REFR(3) + (PRES-REFP(3)) * drdp(3)
        ELSEIF (PRES .GE. REFP(5)) THEN
          TEM = REFR(4) + (PRES-REFP(4)) * drdp(4)
        ELSEIF (PRES .GE. REFP(6)) THEN
          TEM = REFR(5) + (PRES-REFP(5)) * drdp(5)
        ELSE
          TEM = REFR(6)
        ENDIF
!
        tem = 2.0E-4 / tem
        al2 = min(4.0*tem, max(alm, tem))
!
      RETURN
      END
      SUBROUTINE SETQRP
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
      implicit none

      real(kind=kind_phys) tem2,tem1,x,xinc,xmax,xmin
      integer jx
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     XMIN   = 1.0E-6
      XMIN   = 0.0
      XMAX   = 5.0
      XINC   = (XMAX-XMIN)/(NQRP-1)
      C2XQRP = 1.0/XINC
      C1XQRP = 1.0 - XMIN*C2XQRP
      TEM1   = 0.001 ** 0.2046
      TEM2   = 0.001 ** 0.525
      DO JX=1,NQRP
        X         = XMIN + (JX-1)*XINC
        TBQRP(JX) =        X ** 0.1364
        TBQRA(JX) = TEM1 * X ** 0.2046
        TBQRB(JX) = TEM2 * X ** 0.525
      ENDDO    
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
      FUNCTION QRPF(QRP)
!
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
      implicit none

      real(kind=kind_phys) QRP, QRPF, XJ, REAL_NQRP, ONE
      PARAMETER (ONE=1.0)
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NQRP = REAL(NQRP)
      XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
!     XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),FLOAT(NQRP))
      JX   = MIN(XJ,NQRP-ONE)
      QRPF = TBQRP(JX)  + (XJ-JX) * (TBQRP(JX+1)-TBQRP(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
      SUBROUTINE QRABF(QRP,QRAF,QRBF)
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
      implicit none
!
      real(kind=kind_phys) QRP, QRAF, QRBF, XJ, REAL_NQRP, ONE
      PARAMETER (ONE=1.0)
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NQRP=REAL(NQRP)
      XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
      JX   = MIN(XJ,NQRP-ONE)
      XJ   = XJ - JX
      QRAF = TBQRA(JX)  + XJ * (TBQRA(JX+1)-TBQRA(JX))
      QRBF = TBQRB(JX)  + XJ * (TBQRB(JX+1)-TBQRB(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
      SUBROUTINE SETVTP
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : NVTP,C1XVTP,C2XVTP,TBVTP
      implicit none

      real(kind=kind_phys) vtpexp,xinc,x,xmax,xmin
      integer jx
      PARAMETER(VTPEXP=-0.3636)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      XMIN   = 0.05
      XMAX   = 1.5
      XINC   = (XMAX-XMIN)/(NVTP-1)
      C2XVTP = 1.0/XINC
      C1XVTP = 1.0 - XMIN*C2XVTP
      DO JX=1,NVTP
        X         = XMIN + (JX-1)*XINC
        TBVTP(JX) =        X ** VTPEXP
      ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
      FUNCTION VTPF(ROR)
!
      USE MACHINE , ONLY : kind_phys
      use module_ras , only : NVTP,C1XVTP,C2XVTP,TBVTP
      implicit none
      real(kind=kind_phys) ROR, VTPF, XJ, REAL_NVTP, ONE
      PARAMETER (ONE=1.0)
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NVTP = REAL(NVTP)
      XJ   = MIN(MAX(C1XVTP+C2XVTP*ROR,ONE),REAL_NVTP)
      JX   = MIN(XJ,NVTP-ONE)
      VTPF = TBVTP(JX)  + (XJ-JX) * (TBVTP(JX+1)-TBVTP(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
      FUNCTION CLF(PRATE)
!
      USE MACHINE , ONLY : kind_phys
      implicit none
      real(kind=kind_phys) PRATE, CLF
!
      real (kind=kind_phys), parameter :: ccf1=0.30, ccf2=0.09          &
     &,                                   ccf3=0.04, ccf4=0.01          &
     &,                                   pr1=1.0,   pr2=5.0            &
     &,                                   pr3=20.0
!
      if (prate .lt. pr1) then
        clf = ccf1
      elseif (prate .lt. pr2) then
        clf = ccf2
      elseif (prate .lt. pr3) then
        clf = ccf3
      else
        clf = ccf4
      endif
!
      RETURN
      END
