#undef DEBUG
#include <misc.h>


module cldcond 4,4

!---------------------------------------------------------------------------------
! Purpose:
!
! Provides the CAM interface to the prognostic cloud water and ice parameterization
!
! Author: Byron Boville  Sept 04, 2002
!
!---------------------------------------------------------------------------------

  use shr_kind_mod,  only: r8=>shr_kind_r8
  use ppgrid,        only: pcols, pver
  use physconst,     only: gravit, latvap, latice
  use abortutils,    only: endrun

  implicit none
  private
  save

  public :: cldcond_register, cldcond_init_cnst, cldcond_implements_cnst
  public :: cldcond_init, cldcond_tend
  public :: cldcond_zmconv_detrain
  public :: cldcond_sediment

! Private module data

  integer, parameter :: ncnst=2                      ! number of constituents
  character(len=8), dimension(ncnst), parameter :: & ! constituent names
     cnst_names = (/'CLDLIQ', 'CLDICE'/)

  integer :: &
       ixcldice,     &! cloud ice water index
       ixcldliq,     &! cloud liquid water index
       qcwat_idx,    &! qcwat index in physics buffer
       tcwat_idx,    &! tcwat index in physics buffer
       cld_idx,      &! cld index in physics buffer
       lcwat_idx      ! lcwat index in physics buffer

contains

!===============================================================================

  subroutine cldcond_sediment(state, ptend, dtime, & 1,12
       cloud, icefrac, landfrac, ocnfrac, prec, snow, landm, snowh)
!
! Interface to sedimentation of cloud liquid and ice particles
!
! NOTE: initial implementation by B.A. Boville from earlier code by P.J. Rasch
!
! B. A. Boville, Sept 20, 2002
!
!-----------------------------------------------------------------------
    use physics_types,    only: physics_state, physics_ptend
    use pkg_cld_sediment, only: cld_sediment_vel, cld_sediment_tend
    use history,          only: outfld

    implicit none

! Arguments
    type(physics_state), intent(in)    :: state   ! state variables
    type(physics_ptend), intent(inout) :: ptend   ! package tendencies

    real(r8), intent(in)  :: cloud(pcols,pver)    ! cloud fraction
    real(r8), intent(in)  :: icefrac (pcols)      ! sea ice fraction (fraction)
    real(r8), intent(in)  :: landfrac(pcols)      ! land fraction (fraction)
    real(r8), intent(in)  :: ocnfrac (pcols)      ! ocean fraction (fraction)
    real(r8), intent(in)  :: dtime                ! timestep

    real(r8), intent(out) :: prec(pcols)          ! surface flux of total cloud water
    real(r8), intent(out) :: snow(pcols)          ! surface flux of cloud ice
    real(r8), intent(in) :: landm(pcols)          ! land fraction ramped over water
    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)

! Local variables
    integer  :: i,k                               ! loop indexes
    integer  :: ncol                              ! number of atmospheric columns in chunk
    integer  :: lchnk                             ! chunk index
    real(r8) :: rain(pcols)                       ! surface flux of cloud liquid
    real(r8) :: pvliq(pcols,pver+1)               ! vertical velocity of cloud liquid drops (Pa/s)
    real(r8) :: pvice(pcols,pver+1)               ! vertical velocity of cloud ice particles (Pa/s)
!-----------------------------------------------------------------------
    ncol = state%ncol

    ptend%name         = 'pcwsediment'
    ptend%ls           = .TRUE.
    ptend%lq(1)        = .TRUE.
    ptend%lq(ixcldice) = .TRUE.
    ptend%lq(ixcldliq) = .TRUE.

    call cld_sediment_vel (ncol,                                    &
         icefrac, landfrac, ocnfrac, state%pmid, state%pdel, state%t, &
         cloud, state%q(:,:,ixcldliq), state%q(:,:,ixcldice), pvliq, pvice, landm, snowh)

    call cld_sediment_tend (ncol, dtime ,                                       &
         state%pint, state%pmid           , state%pdel            , state%t     , &
         cloud  ,    state%q(:,:,ixcldliq), state%q(:,:,ixcldice) , pvliq, pvice, &
         ptend%q(:,:,ixcldliq), ptend%q(:,:,ixcldice), ptend%q(:,:,1), ptend%s  , &
         rain   , snow   )

! convert rain and snow from kg/m2 to m/s
    snow(:ncol) = snow(:ncol)/1000.
    rain(:ncol) = rain(:ncol)/1000.
! compute total precip (m/s)
    prec(:ncol) = rain(:ncol) + snow(:ncol)

! record history variables
    lchnk = state%lchnk
    call outfld('DQSED'   ,ptend%q(:,:,1)       , pcols,lchnk)
    call outfld('DISED'   ,ptend%q(:,:,ixcldice), pcols,lchnk)
    call outfld('DLSED'   ,ptend%q(:,:,ixcldliq), pcols,lchnk)
    call outfld('HSED'    ,ptend%s              , pcols,lchnk)
    call outfld('PRECSED' ,prec                 , pcols,lchnk)
    call outfld('SNOWSED' ,snow                 , pcols,lchnk)
    call outfld('RAINSED' ,rain                 , pcols,lchnk)

    return
  end subroutine cldcond_sediment

!===============================================================================

  subroutine cldcond_zmconv_detrain(dlf, cld, state, ptend) 1,3
!
! Partition the detrained condensed water from the ZM convection scheme.
!
! The ZM scheme does not have an ice phase, so the detrained water is paritioned
! soley between cloud liquid water and the environment. The ice/liquid partitioning 
! happens in cldcond_tend.
!
! NOTE: initial implementation by B.A. Boville just moves code here from TPHYSBC.
!
! B. A. Boville, Sept 09, 2002
!
!-----------------------------------------------------------------------
    use physics_types, only: physics_state, physics_ptend
    use history,       only: outfld

    implicit none

! Arguments
    type(physics_state), intent(in  )  :: state   ! state variables
    type(physics_ptend), intent(inout) :: ptend   ! package tendencies

    real(r8), intent(in) :: dlf(pcols,pver)       ! detrained water from ZM
    real(r8), intent(in) :: cld(pcols,pver)       ! cloud fraction

! Local variables
    integer :: i,k                            ! loop indexes
!-----------------------------------------------------------------------

    ptend%name         = 'pcwdetrain'
!!$    ptend%ls           = .TRUE.
!!$    ptend%lq(1)        = .TRUE.
!!$    ptend%lq(ixcldice) = .TRUE.
    ptend%lq(ixcldliq) = .TRUE.
!
! Put all of the detraining cloud water from convection into the large scale cloud.
! It all goes in liquid for the moment.
    do k = 1,pver
       do i = 1,state%ncol
!!$          ptend%q(i,k,1)        = dlf(i,k) * (1.-cld(i,k))
!!$          ptend%s(i,k)          =-dlf(i,k) * (1.-cld(i,k))*latvap
!!$          ptend%q(i,k,ixcldice) = 0.
!!$          ptend%q(i,k,ixcldliq) = dlf(i,k) * cld(i,k)
          ptend%q(i,k,ixcldliq) = dlf(i,k)
       end do
    end do
    call outfld('ZMDLF' ,dlf  , pcols,state%lchnk)

    return
  end subroutine cldcond_zmconv_detrain

!===============================================================================

  subroutine cldcond_register() 1,11
!
! Register the constituents (cloud liquid and cloud ice) and the fields
! in the physics buffer.
! 
!-----------------------------------------------------------------------
    use constituents, only: cnst_add, advected, nonadvec
    use physconst,    only: mwdry, cpair
    use phys_buffer,  only: pbuf_times, pbuf_add

    implicit none

!    logical, parameter :: cldw_adv=.false.  ! true => cloud water is treated as advected tracer
    logical, parameter :: cldw_adv=.true.  ! true => cloud water is treated as advected tracer

    integer flag
!-----------------------------------------------------------------------

! Register cloud water and determine index (either advected or non-adv).
    if (cldw_adv) then
       flag = advected
    else
       flag = nonadvec
    endif
    call cnst_add(cnst_names(1), flag, mwdry, cpair, 0._r8, ixcldliq, &
         longname='Grid box averaged liquid condensate amount')
    call cnst_add(cnst_names(2), flag, mwdry, cpair, 0._r8, ixcldice, &
         longname='Grid box averaged ice condensate amount')

! Request physics buffer space for fields that persist across timesteps.
    call pbuf_add('QCWAT', 'global', 1,pver,pbuf_times, qcwat_idx)
    call pbuf_add('TCWAT', 'global', 1,pver,pbuf_times, tcwat_idx)
    call pbuf_add('CLD',   'global', 1,pver,pbuf_times, cld_idx)
    call pbuf_add('LCWAT', 'global', 1,pver,pbuf_times, lcwat_idx)

    call pbuf_add('QINI' , 'physpkg', 1,pver,      1, flag)
    call pbuf_add('TINI' , 'physpkg', 1,pver,      1, flag)

  end subroutine cldcond_register

!===============================================================================


  function cldcond_implements_cnst(name)
!----------------------------------------------------------------------- 
! 
! Purpose: return true if specified constituent is implemented by this package
! 
! Author: B. Eaton
! 
!-----------------------------------------------------------------------
     implicit none
!-----------------------------Arguments---------------------------------

     character(len=*), intent(in) :: name   ! constituent name
     logical :: cldcond_implements_cnst     ! return value
!---------------------------Local workspace-----------------------------
     integer :: m
!-----------------------------------------------------------------------

     cldcond_implements_cnst = .false.
     do m = 1, ncnst
        if (name == cnst_names(m)) then
           cldcond_implements_cnst = .true.
           return
        end if
     end do
  end function cldcond_implements_cnst

!===============================================================================

  subroutine cldcond_init_cnst(name, q) 1,1
!
! Initialize the cloud water mixing ratios (liquid and ice), if they are
! not read from the initial file
! 
!-----------------------------------------------------------------------
    use pmgrid,        only: plon, plev, plat

    implicit none

! Arguments
    character(len=*), intent(in)  :: name                ! constituent name
    real(r8),         intent(out) :: q(plon,plev,plat)   ! mass mixing ratio
!-----------------------------------------------------------------------

    if ( name == 'CLDLIQ' ) then
       q = 0.0
       return
    else if ( name == 'CLDICE' ) then
       q = 0.0
       return
    end if

  end subroutine cldcond_init_cnst

!===============================================================================

  subroutine cldcond_init() 1,30
!
! Initialize the cloud water parameterization
! 
!-----------------------------------------------------------------------
    use cldwat,        only: inimc
    use history,       only: addfld, add_default, phys_decomp
    use physconst,     only: tmelt, rh2o, rhodair

    implicit none
!-----------------------------------------------------------------------

! initialization routine for prognostic cloud water
    call inimc (tmelt, rhodair/1000.0, gravit, rh2o )

! register history variables
    call addfld ('FWAUT   ','fraction',pver, 'A','Relative importance of liquid autoconversion' ,phys_decomp)
    call addfld ('FSAUT   ','fraction',pver, 'A','Relative importance of ice autoconversion'    ,phys_decomp)
    call addfld ('FRACW   ','fraction',pver, 'A','Relative  importance of rain accreting liquid',phys_decomp)
    call addfld ('FSACW   ','fraction',pver, 'A','Relative  importance of snow accreting liquid',phys_decomp)
    call addfld ('FSACI   ','fraction',pver, 'A','Relative  importance of snow accreting ice'   ,phys_decomp)
    call addfld ('CME     ','kg/kg/s ',pver, 'A','Rate of cond-evap within the cloud'           ,phys_decomp)
    call addfld ('ZMDLF   ','kg/kg/s ',pver, 'A','Detrained liquid water from ZM convection'    ,phys_decomp)
    call addfld ('PRODPREC','kg/kg/s ',pver, 'A','Rate of conversion of condensate to precip'   ,phys_decomp)
    call addfld ('EVAPPREC','kg/kg/s ',pver, 'A','Rate of evaporation of falling precip'        ,phys_decomp)
    call addfld ('EVAPSNOW','kg/kg/s ',pver, 'A','Rate of evaporation of falling snow'          ,phys_decomp)
    call addfld ('HPROGCLD','W/kg'    ,pver, 'A','Heating from prognostic clouds'               ,phys_decomp)
    call addfld ('HEVAP   ','W/kg'    ,pver, 'A','Heating from evaporation of falling precip'   ,phys_decomp)
    call addfld ('HMELT   ','W/kg'    ,pver, 'A','Heating from snow melt'                       ,phys_decomp)
    call addfld ('HREPART ','W/kg'    ,pver, 'A','Heating from cloud ice/liquid repartitioning' ,phys_decomp)
    call addfld ('FICE    ','fraction',pver, 'A','Fractional ice content within cloud'          ,phys_decomp)
    call addfld ('ICWMR   ','kg/kg   ',pver, 'A','Prognostic in-cloud water mixing ratio'       ,phys_decomp)
    call addfld ('ICIMR   ','kg/kg   ',pver, 'A','Prognostic in-cloud ice mixing ratio'         ,phys_decomp)
    call addfld ('PCSNOW  ','m/s'     ,1   , 'A','Snow fall from prognostic clouds'             ,phys_decomp)

    call addfld ('DQSED   ','kg/kg/s' ,pver, 'A','Water vapor tendency from cloud sedimentation',phys_decomp)
    call addfld ('DLSED   ','kg/kg/s' ,pver, 'A','Cloud liquid tendency from sedimentation'     ,phys_decomp)
    call addfld ('DISED   ','kg/kg/s' ,pver, 'A','Cloud ice tendency from sedimentation'        ,phys_decomp)
    call addfld ('HSED    ','W/kg'    ,pver, 'A','Heating from cloud sediment evaporation'      ,phys_decomp)
    call addfld ('SNOWSED ','m/s'     ,1   , 'A','Snow from cloud ice sedimentation'            ,phys_decomp)
    call addfld ('RAINSED ','m/s'     ,1   , 'A','Rain from cloud liquid sedimentation'         ,phys_decomp)
    call addfld ('PRECSED ','m/s'     ,1   , 'A','Precipitation from cloud sedimentation'       ,phys_decomp)

    call add_default ('FICE    ', 1, ' ')
!  These fields removed 10/30/2003 per CRB decision.
!    call add_default ('CME     ', 1, ' ')
!    call add_default ('ZMDLF   ', 1, ' ')
!    call add_default ('PRODPREC', 1, ' ')
!    call add_default ('EVAPPREC', 1, ' ')
!    call add_default ('EVAPSNOW', 1, ' ')
!    call add_default ('HPROGCLD', 1, ' ')
!    call add_default ('HEVAP   ', 1, ' ')
!    call add_default ('HMELT   ', 1, ' ')
!    call add_default ('HREPART ', 1, ' ')
!    call add_default ('ICWMR   ', 1, ' ')
!    call add_default ('ICIMR   ', 1, ' ')
!    call add_default ('PCSNOW  ', 1, ' ')
!    call add_default ('DQSED   ', 1, ' ')
!    call add_default ('DLSED   ', 1, ' ')
!    call add_default ('DISED   ', 1, ' ')
!    call add_default ('HSED    ', 1, ' ')
!    call add_default ('SNOWSED ', 1, ' ')
!    call add_default ('RAINSED ', 1, ' ')
!    call add_default ('PRECSED ', 1, ' ')

    return
  end subroutine cldcond_init

!===============================================================================
!!$  subroutine cldcond_tend(state, ptend, dt, pbuf)

  subroutine cldcond_tend(state, ptend, dt, & 1,28
       tcwato, qcwato, lcwato, precip, snow, icefrac, rhdfda, rhu00, cldn, evapprec, prodprec, cme, snowh)
!
! Compute the tendency of cloud water mixing ratios (liquid and ice) from
! microphysical processes and sedimenation of cloud particles.
! 
! Author: Byron Boville  Sept 04, 2002
!  modified pjr: march 2003 to repartition snow/rain
!
!-----------------------------------------------------------------------

    use physics_types, only: physics_state, physics_ptend
!!$    use phys_buffer,   only: pbuf_size_max, pbuf_fld
    use history,       only: outfld
    use cldwat,        only: pcond, cldwat_fice
    use physconst,     only: tmelt

    implicit none

! Arguments
    type(physics_state), intent(inout) :: state          ! state variables
    type(physics_ptend), intent(inout) :: ptend          ! package tendencies
    real(r8),            intent(in)    :: dt             ! timestep
!!$    type(pbuf_fld), intent(inout), dimension(pbuf_size_max) :: pbuf  ! physics buffer
    real(r8), intent(in)  :: tcwato(pcols,pver)        !cloud water old temperature
    real(r8), intent(in)  :: qcwato(pcols,pver)        ! cloud water old q
    real(r8), intent(in)  :: lcwato(pcols,pver)        ! cloud liquid water old q
    real(r8), intent(in)  :: icefrac(pcols)            ! sea ice fraction (fraction)
    real(r8), intent(in)  :: rhdfda(pcols,pver)        ! dRh/dcloud, old 
    real(r8), intent(in)  :: rhu00 (pcols,pver)        ! Rh threshold for cloud, old
    real(r8), intent(in)  :: cldn(pcols,pver)          !new cloud fraction
    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)

    real(r8), intent(out) :: precip(pcols)             ! sfc flux of precip (m/s)
    real(r8), intent(out) :: snow (pcols)              ! sfc flux of snow   (m/s)
    real(r8), intent(out) :: evapprec(pcols,pver)          ! local evaporation of precipitation
    real(r8), intent(out) :: prodprec(pcols,pver)          ! local production of precipitation
    real(r8), intent(out) :: cme     (pcols,pver)          ! local condensation - evaporation of cloud water

! Local variables
    integer :: lchnk                          ! chunk identifier
    integer :: ncol                           ! number of atmospheric columns in chunk
    integer :: i,k                            ! loop indexes
!!$    real(r8), pointer, dimension(:)   :: buffld1  ! physics buffer field1
!!$    real(r8), pointer, dimension(:,:) :: buffld2  ! physics buffer field2
    real(r8) :: rdt                          ! 1./dt
    real(r8) :: qtend (pcols,pver)            ! moisture tendencies
    real(r8) :: ttend (pcols,pver)            ! temperature tendencies
    real(r8) :: ltend (pcols,pver)            ! cloud liquid water tendencies
!    real(r8) :: cme     (pcols,pver)          ! local condensation - evaporation of cloud water
    real(r8) :: evapheat(pcols,pver)          ! heating rate due to evaporation of precip
!    real(r8) :: evapprec(pcols,pver)          ! local evaporation of precipitation
    real(r8) :: evapsnow(pcols,pver)          ! local evaporation of snow
    real(r8) :: prfzheat(pcols,pver)          ! heating rate due to freezing of precip (W/kg)
    real(r8) :: meltheat(pcols,pver)          ! heating rate due to phase change of precip
!    real(r8) :: prodprec(pcols,pver)          ! local production of precipitation
    real(r8) :: prodsnow(pcols,pver)          ! local production of snow
    real(r8) :: totcw   (pcols,pver)          ! total cloud water mixing ratio
    real(r8) :: fice    (pcols,pver)          ! Fractional ice content within cloud
    real(r8) :: fsnow   (pcols,pver)          ! Fractional snow production
    real(r8) :: repartht(pcols,pver)          ! heating rate due to phase repartition of input precip
    real(r8) :: icimr(pcols,pver)             ! in cloud ice mixing ratio
    real(r8) :: icwmr(pcols,pver)             ! in cloud water mixing ratio
    real(r8) fwaut(pcols,pver)              
    real(r8) fsaut(pcols,pver)              
    real(r8) fracw(pcols,pver)              
    real(r8) fsacw(pcols,pver)              
    real(r8) fsaci(pcols,pver)              
    real(r8) ice2pr(pcols,pver)   ! rate of conversion of ice to precip
    real(r8) liq2pr(pcols,pver)   ! rate of conversion of liquid to precip
    real(r8) liq2snow(pcols,pver)   ! rate of conversion of liquid to snow
    real(r8) hs1, qv1, ql1, qi1, qs1, qr1, fice2, pr1, w1, w2, w3, fliq, res(pcols,pver)
    real(r8) temp(pcols), w4, wl, wv, wi, wlf, wvf, wif, qif, qlf, qvf, res2

!-----------------------------------------------------------------------

! Set output flags
    ptend%name         = 'cldwat'
    ptend%ls           = .true.
    ptend%lq(1)        = .true.
    ptend%lq(ixcldice) = .true.
    ptend%lq(ixcldliq) = .true.

! Initialize chunk id and size
    lchnk = state%lchnk
    ncol  = state%ncol

! associate local pointers with fields in the physics buffer
!!$    buffld1 => pbuf(ixbuffld1)%fld_ptr(1,1:pcols,1,     lchnk,1)
!!$    buffld2 => pbuf(ixbuffld2)%fld_ptr(1,1:pcols,1:pver,lchnk,1)

! Define fractional amount of cloud condensate in ice phase
    call cldwat_fice(ncol, state%t, fice, fsnow)

! compute total cloud water
    totcw(:ncol,:pver) = state%q(:ncol,:pver,ixcldice) + state%q(:ncol,:pver,ixcldliq)

! save input cloud ice
    repartht(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)

! Repartition ice and liquid
    state%q(:ncol,:pver,ixcldice) = totcw(:ncol,:pver) * fice(:ncol,:pver)
    state%q(:ncol,:pver,ixcldliq) = totcw(:ncol,:pver) * (1.0_r8 - fice(:ncol,:pver))

! Determine heating from change in cloud ice
    repartht(:ncol,:pver) = latice/dt * (state%q(:ncol,:pver,ixcldice) - repartht(:ncol,:pver))

! calculate the tendencies for moisture, temperature and cloud fraction
    rdt = 1./dt
    qtend(:ncol,:pver) = (state%q(:ncol,:pver,1) - qcwato(:ncol,:pver))*rdt
    ttend(:ncol,:pver) = (state%t(:ncol,:pver)   - tcwato(:ncol,:pver))*rdt
    ltend(:ncol,:pver) = (totcw  (:ncol,:pver)   - lcwato(:ncol,:pver))*rdt

! call microphysics package to calculate tendencies
    call t_startf('pcond')
    call pcond (lchnk,   ncol, &
         state%t  , ttend      , state%q(1,1,1), qtend   , state%omega, &
         totcw    , state%pmid , state%pdel , cldn       , fice       , fsnow, &
         cme      , prodprec   , prodsnow, evapprec   , evapsnow   , evapheat   , prfzheat, &
         meltheat , precip     , snow       , dt         , fwaut      , &
         fsaut    , fracw      , fsacw      , fsaci      , ltend      , &
         rhdfda   , rhu00      , icefrac    , state%zi   , ice2pr, liq2pr, liq2snow, snowh)
    call t_stopf('pcond')

! make it interactive
    do k = 1,pver
       do i = 1,ncol
          ptend%s(i,k)          = cme(i,k) * (latvap + latice*fice(i,k)) &
               + evapheat(i,k) + prfzheat(i,k) + meltheat(i,k) + repartht(i,k)

          ptend%q(i,k,1)        =-cme(i,k) + evapprec(i,k)

          ptend%q(i,k,ixcldice) =cme(i,k)*fice(i,k)      - ice2pr(i,k)
          ptend%q(i,k,ixcldliq) =cme(i,k)*(1.-fice(i,k)) - liq2pr(i,k)
       end do
    end do

! sanity checks to be removed after 5 year run
    res(:ncol,:) = state%q(:ncol,:,1) + ptend%q(:ncol,:,1)*dt
    if (any(res(:ncol,:) < 0)) then
       do k = 1,pver
          do i = 1,ncol
             if (res(i,k).lt.0.) then
                write (6,*) ' predicted neg vapor ', i,k,lchnk, res(i,k)
                write (6,*) ' oldq, cme(i,k)*dt, evapprec(i,k)*dt ', state%q(i,k,1), &
                   -cme(i,k)*dt, evapprec(i,k)*dt, ptend%q(i,k,1)*dt
                call endrun('CLDCOND_TEND')
             endif
          end do
       end do
    endif
    res(:ncol,:) = state%q(:ncol,:,ixcldice) + ptend%q(:ncol,:,ixcldice)*dt
    if (any(res(:ncol,:) < 0)) then
       do k = 1,pver
          do i = 1,ncol
             if (res(i,k).lt.0.) then
                write (6,*) ' cldcond_tend: predicted neg ice ', i,k,lchnk, res(i,k), fice(i,k)
                write (6,*) ' oldice, cme(i,k)*dt*(fice(i,k)), ice2pr*dt, ptend_ice*dt', &
                   state%q(i,k,ixcldice), &
                   cme(i,k)*dt*(fice(i,k)), ice2pr(i,k)*dt, ptend%q(i,k,ixcldice)*dt
                write (6,*) ' oldcond, cme(i,k)*dt, prodprec*dt, ptend_cond*dt', &
                   state%q(i,k,ixcldice)+state%q(i,k,ixcldliq), &
                   cme(i,k)*dt, prodprec(i,k)*dt, &
                   (ptend%q(i,k,ixcldice)+ptend%q(i,k,ixcldliq))*dt
                call endrun('CLDCOND_TEND')
             endif
          end do
       end do
    endif
! end sanity check

#ifdef DEBUG
    do i = 1,ncol
       pr1 = 0
       hs1 = 0
       qv1 = 0
       ql1 = 0
       qi1 = 0
       qs1 = 0
       qr1 = 0
       w1 = 0
       wl = 0
       wv = 0
       wi = 0
       wlf = 0
       wvf = 0
       wif = 0
       do k = 1,pver

          if (lchnk.eq.248.and.i.eq.12) then

             write (6,*) 
             write (6,*) ' input state, t, q, l, i ', k, state%t(i,k), state%q(i,k,1), state%q(i,k,ixcldliq),  state%q(i,k,ixcldice)
             write (6,*) ' rain, snow, total from components before accumulation ', qr1, qs1, qr1+qs1
             write (6,*) ' total precip before accumulation                      ', k, pr1

             wv = wv + state%q(i,k,1       )*state%pdel(i,k)/gravit
             wl = wl + state%q(i,k,ixcldliq)*state%pdel(i,k)/gravit
             wi = wi + state%q(i,k,ixcldice)*state%pdel(i,k)/gravit

             qvf = state%q(i,k,1) + ptend%q(i,k,1)*dt
             qlf = state%q(i,k,ixcldliq) + ptend%q(i,k,ixcldliq)*dt
             qif = state%q(i,k,ixcldice) + ptend%q(i,k,ixcldice)*dt

             if (qvf.lt.0.) then
                write (6,*) ' qvf is negative *******', qvf
             endif
             if (qlf.lt.0.) then
                write (6,*) ' qlf is negative *******', qlf
             endif
             if (qif.lt.0.) then
                write (6,*) ' qif is negative *******', qif
             endif
             write (6,*) ' qvf, qlf, qif ', qvf, qlf, qif

             wvf = wvf + qvf*state%pdel(i,k)/gravit
             wlf = wlf + qlf*state%pdel(i,k)/gravit
             wif = wif + qif*state%pdel(i,k)/gravit

             hs1 = hs1 + ptend%s(i,k)*state%pdel(i,k)/gravit
             pr1 = pr1 + state%pdel(i,k)/gravit*(prodprec(i,k)-evapprec(i,k))
             qv1 = qv1 - (cme(i,k)-evapprec(i,k))*state%pdel(i,k)/gravit    ! vdot
             w1  = w1  + (cme(i,k)-prodprec(i,k))*state%pdel(i,k)/gravit    ! cdot
             qi1 = qi1 + ((cme(i,k))*fice(i,k)        -ice2pr(i,k) )*state%pdel(i,k)/gravit   ! idot
             ql1 = ql1 + ((cme(i,k))*(1._r8-fice(i,k))-liq2pr(i,k) )*state%pdel(i,k)/gravit   ! ldot

             qr1 = qr1 &
                  + ( liq2pr(i,k)-liq2snow(i,k)   &     ! production of rain
                   -(evapprec(i,k)-evapsnow(i,k)) &     ! rain evaporation
                    )*state%pdel(i,k)/gravit
             qs1 = qs1 &
                  + ( ice2pr(i,k) + liq2snow(i,k) &     ! production of snow.Note last term has phase change
                     -evapsnow(i,k)               &     ! snow evaporation
                    )*state%pdel(i,k)/gravit

             if (state%t(i,k).gt.tmelt) then
                qr1 = qr1 + qs1
                qs1 = 0.
             endif
             write (6,*) ' rain, snow, total after accumulation ', qr1, qs1, qr1+qs1
             write (6,*) ' total precip after accumulation      ', k, pr1
             write (6,*)
             write (6,*) ' layer prodprec, evapprec, pdel ', prodprec(i,k), evapprec(i,k), state%pdel(i,k)
             write (6,*) ' layer prodsnow, ice2pr+liq2snow ', prodsnow(i,k), ice2pr(i,k)+liq2snow(i,k)
             write (6,*) ' layer prodprec-prodsnow, liq2pr-liq2snow ', prodprec(i,k)-prodsnow(i,k), liq2pr(i,k)-liq2snow(i,k)
             write (6,*) ' layer evapsnow, evaprain ', k, evapsnow(i,k), evapprec(i,k)-evapsnow(i,k)
             write (6,*) ' layer ice2pr, liq2pr, liq2snow ', ice2pr(i,k), liq2pr(i,k), liq2snow(i,k)
             write (6,*) ' layer ice2pr+liq2pr, prodprec ', ice2pr(i,k)+liq2pr(i,k), prodprec(i,k)
             write (6,*)
             write (6,*) ' qv1 vapor removed from col after accum  (vdot)   ', k, qv1
             write (6,*) ' - (precip produced - vapor removed) after accum  ', k, -pr1-qv1
             write (6,*) ' condensate produce after accum                   ', k, w1
             write (6,*) ' liq+ice tends accum                              ', k, ql1+qi1
             write (6,*) ' change in total water after accum                ', k, qv1+ql1+qi1
             write (6,*) ' imbalance in colum after accum                   ', k, qs1+qr1+qv1+ql1+qi1
             write (6,*) ' fice at this lev ', fice(i,k)
             write (6,*)

             res2 = abs((qs1+qr1+qv1+ql1+qi1)/max(abs(qv1),abs(ql1),abs(qi1),abs(qs1),abs(qr1),1.e-36))
             write (6,*) ' relative residual in column method 1             ', k, res2

             write (6,*) ' relative residual in column method 2             ', k, abs((qs1+qr1+qv1+ql1+qi1)/max(abs(qv1+ql1+qi1),1.e-36))
!            if (abs((qs1+qr1+qv1+ql1+qi1)/(qs1+qr1+1.e-36)).gt.1.e-14) then
             if (res2.gt.1.e-14) then
                call endrun ('CLDCOND_TEND')
             endif

!             w3  = cme(i,k) * (latvap + latice*fice(i,k)) &
!               + evapheat(i,k) + prfzheat(i,k) + meltheat(i,k)

             res2 = qs1+qr1-pr1
             w4 = max(abs(qs1),abs(qr1),abs(pr1)) 
             if (w4.gt.0.)  then
                if (res2/w4.gt.1.e-14) then
                   write (6,*) ' imbalance in precips calculated two ways '
                   write (6,*) ' res2/w4, pr1, qr1, qs1, qr1+qs1 ', &
                        res2/w4, pr1, qr1, qs1, qr1+qs1
!                   call endrun()
                endif
             endif
             if (k.eq.pver) then
                write (6,*) ' pcond returned precip, rain and snow rates ', precip(i), precip(i)-snow(i), snow(i)
                write (6,*) ' I calculate ', pr1, qr1, qs1
!               call endrun
                write (6,*) ' byrons water check ', wv+wl+wi-pr1*dt, wvf+wlf+wif
             endif
          write (6,*)
          endif
       end do
    end do
#endif

#ifdef DEBUG
    if (.true.) then
    do i = 1,ncol
       if (snow(i).gt.0.01/8.64e4.and.state%t(i,pver).gt.273.16) then
          write (6,*) ' cldcond: snow, temp, ', i, lchnk, &
               snow(i), state%t(i,pver)
          write (6,*) ' t ', state%t(i,:)
          write (6,*) ' fsaut ', fsaut(i,:)
          write (6,*) ' fsacw ', fsacw(i,:)
          write (6,*) ' fsaci ', fsaci(i,:)
          write (6,*) ' meltheat ', meltheat(i,:)
          call endrun ('CLDCOND_TEND')
       endif

       if (snow(i)*8.64e4.lt.-1.e-5) then
          write (6,*) ' neg snow ', snow(i)*8.64e4
          write (6,*) ' cldcond: snow, temp, ', i, lchnk, &
               snow(i), state%t(i,pver)
          write (6,*) ' t ', state%t(i,:)
          write (6,*) ' fsaut ', fsaut(i,:)
          write (6,*) ' fsacw ', fsacw(i,:)
          write (6,*) ' fsaci ', fsaci(i,:)
          write (6,*) ' meltheat ', meltheat(i,:)
          call endrun ('CLDCOND_TEND')
       endif
    end do
    endif
#endif

! Compute in cloud ice and liquid mixing ratios
    do k=1,pver
       do i = 1,ncol
          icimr(i,k) = (state%q(i,k,ixcldice) + dt*ptend%q(i,k,ixcldice)) / max(0.01_r8,cldn(i,k))
          icwmr(i,k) = (state%q(i,k,ixcldliq) + dt*ptend%q(i,k,ixcldliq)) / max(0.01_r8,cldn(i,k))
       end do
    end do

! convert precipitation from kg/m2 to m/s
    snow  (:ncol) = snow  (:ncol)/1000.
    precip(:ncol) = precip(:ncol)/1000.

! record history variables
    call outfld('FWAUT',fwaut,    pcols,lchnk)
    call outfld('FSAUT',fsaut,    pcols,lchnk)
    call outfld('FRACW',fracw,    pcols,lchnk)
    call outfld('FSACW',fsacw,    pcols,lchnk)
    call outfld('FSACI',fsaci,    pcols,lchnk)
    call outfld('ICIMR',icimr,    pcols,lchnk)
    call outfld('ICWMR',icwmr,    pcols,lchnk)

    call outfld('PCSNOW'  ,snow    , pcols,lchnk)
    call outfld('FICE'    ,fice    , pcols,lchnk)
    call outfld('CME'     ,cme     , pcols,lchnk)
    call outfld('PRODPREC',prodprec, pcols,lchnk)
    call outfld('EVAPPREC',evapprec, pcols,lchnk)
    call outfld('EVAPSNOW',evapsnow, pcols,lchnk)
    call outfld('HPROGCLD',ptend%s , pcols,lchnk)
    call outfld('HEVAP   ',evapheat, pcols,lchnk)
    call outfld('HMELT'   ,meltheat, pcols,lchnk)
    call outfld('HREPART' ,repartht, pcols,lchnk)

! update boundary quantities
!!$    ptend%hflx_srf = 0.
!!$    ptend%hflx_top = 0.
!!$    ptend%cflx_srf = 0.
!!$    ptend%cflx_top = 0.
  end subroutine cldcond_tend
end module cldcond