#undef OLDLIQSED

module pkg_cld_sediment 1,6

!---------------------------------------------------------------------------------
! Purpose:
!
! Contains routines to compute tendencies from sedimentation of cloud liquid and 
! ice particles
!
! Author: Byron Boville  Sept 19, 2002 from code by P. J. Rasch
!
!---------------------------------------------------------------------------------

  use shr_kind_mod,  only: r8=>shr_kind_r8
  use ppgrid,        only: pcols, pver, pverp
  use physconst,     only: gravit, latvap, latice, rair, rhoh2o
  use cldwat,        only: icritc
  use pkg_cldoptics, only: reitab, reltab
  use abortutils,    only: endrun

  private
  public :: cld_sediment_vel, cld_sediment_tend


  real (r8), parameter :: vland  = 1.5            ! liquid fall velocity over land  (cm/s)
  real (r8), parameter :: vocean = 2.8            ! liquid fall velocity over ocean (cm/s)
  real (r8), parameter :: mxsedfac   = 0.99       ! maximum sedimentation flux factor

  logical,   parameter :: stokes = .true.         ! use Stokes velocity instead of McFarquhar and Heymsfield

! parameter for modified McFarquhar and Heymsfield
  real (r8), parameter :: vice_small = 1.         ! ice fall velocity for small concentration (cm/s)

! parameters for Stokes velocity
  real (r8), parameter :: eta =  1.7e-5           ! viscosity of air (kg m / s)
  real (r8), parameter :: r40 =  40.              !  40 micron radius
  real (r8), parameter :: r400= 400.              ! 400 micron radius
  real (r8), parameter :: v40 = (2./9.) * rhoh2o * gravit/eta * r40**2 * 1.e-12  
                                                  ! Stokes fall velocity of 40 micron sphere (m/s)
  real (r8), parameter :: v400= 1.00              ! fall velocity of 400 micron sphere (m/s)
  real (r8), parameter :: vslope = (v400 - v40)/(r400 -r40) ! linear slope for large particles m/s/micron


contains

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

  subroutine cld_sediment_vel (ncol,                               & 1,2
       icefrac , landfrac, ocnfrac , pmid    , pdel    , t       , &
       cloud   , cldliq  , cldice  , pvliq   , pvice   , landm, snowh)

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

! Compute gravitational sedimentation velocities for cloud liquid water
! and ice, based on Lawrence and Crutzen (1998).

! LIQUID

! The fall velocities assume that droplets have a gamma distribution
! with effective radii for land and ocean as assessed by Han et al.;
! see Lawrence and Crutzen (1998) for a derivation.

! ICE

! The fall velocities are based on data from McFarquhar and Heymsfield
! or on Stokes terminal velocity for spheres and the effective radius.

! NEED TO BE CAREFUL - VELOCITIES SHOULD BE AT THE *LOWER* INTERFACE
! (THAT IS, FOR K+1), FLUXES ARE ALSO AT THE LOWER INTERFACE (K+1), 
! BUT MIXING RATIOS ARE AT THE MIDPOINTS (K)...

! NOTE THAT PVEL IS ON PVERP (INTERFACES), WHEREAS VFALL IS FOR THE CELL
! AVERAGES (I.E., MIDPOINTS); ASSUME THE FALL VELOCITY APPLICABLE TO THE 
! LOWER INTERFACE (K+1) IS THE SAME AS THAT APPLICABLE FOR THE CELL (V(K))

!-----------------------------------------------------------------------
!     MATCH-MPIC version 2.0, Author: mgl, March 1998
! adapted by P. J. Rasch
!            B. A. Boville, September 19, 2002
!            P. J. Rasch    May 22, 2003 (added stokes flow calc for liquid
!                                         drops based on effect radii)
!-----------------------------------------------------------------------

    implicit none

! Arguments
    integer, intent(in) :: ncol                     ! number of colums to process

    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)  :: pmid  (pcols,pver)     ! pressure of midpoint levels (Pa)
    real(r8), intent(in)  :: pdel  (pcols,pver)     ! pressure diff across layer (Pa)
    real(r8), intent(in)  :: cloud (pcols,pver)     ! cloud fraction (fraction)
    real(r8), intent(in)  :: t     (pcols,pver)     ! temperature (K)
    real(r8), intent(in)  :: cldliq(pcols,pver)     ! cloud water, liquid (kg/kg)
    real(r8), intent(in)  :: cldice(pcols,pver)     ! cloud water, ice    (kg/kg)
    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)

    real(r8), intent(out) :: pvliq (pcols,pverp)    ! vertical velocity of cloud liquid drops (Pa/s)
    real(r8), intent(out) :: pvice (pcols,pverp)    ! vertical velocity of cloud ice particles (Pa/s)
    real(r8), intent(in) :: landm(pcols)            ! land fraction ramped over water
! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))

! Local variables
    real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
    real (r8) :: vfall                              ! settling velocity of cloud particles (m/s)
    real (r8) :: icice                              ! in cloud ice water content (kg/kg)
    real (r8) :: iciwc                              ! in cloud ice water content in g/m3
    real (r8) :: icefac
    real (r8) :: logiwc

    real (r8) :: rei(pcols,pver)                    ! effective radius of ice particles (microns)
    real (r8) :: rel(pcols,pver)                    ! effective radius of liq particles (microns)
    real(r8)  pvliq2 (pcols,pverp)    ! vertical velocity of cloud liquid drops (Pa/s)

    integer i,k

    real (r8) :: lbound, ac, bc, cc

!-----------------------------------------------------------------------
!--------------------- liquid fall velocity ----------------------------
!-----------------------------------------------------------------------

    ! compute air density
    rho(:ncol,:) = pmid(:ncol,:) / (rair * t(:ncol,:))

    pvliq(:ncol,:) = 0.

    ! get effective radius of liquid drop
    call reltab(ncol, t, landfrac, landm, icefrac, rel, snowh)

    do k = 1,pver
       do i = 1,ncol
          if (cloud(i,k) > 0. .and. cldliq(i,k) > 0.) then

#ifdef OLDLIQSED
! oldway
             ! merge the liquid fall velocities for land and ocean (cm/s)
             ! SHOULD ALSO ACCOUNT FOR ICEFRAC
             vfall = vland*landfrac(i) + vocean*(1.-landfrac(i))
!!$          vfall = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)

             ! convert the fall speed to pressure units, but do not apply the traditional
             ! negative convention for pvel.
             pvliq(i,k+1) = vfall     &
                  * 0.01                 & ! cm to meters
                  * rho(i,k)*gravit        ! meters/sec to pascals/sec
#else

! newway
             if (rel(i,k) < 40. ) then
                vfall = 2./9. * rhoh2o * gravit * rel(i,k)**2 / eta  * 1.e-12  ! micons^2 -> m^2
             else
                vfall = v40 + vslope * (rel(i,k)-r40)      ! linear above 40 microns
             end if
             ! convert the fall speed to pressure units
             ! but do not apply the traditional
             ! negative convention for pvel.
!             pvliq2(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
             pvliq(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
#endif
          end if
       end do
    end do

!-----------------------------------------------------------------------
!--------------------- ice fall velocity -------------------------------
!--------------------- stokes terminal velocity < 40 microns -----------
!-----------------------------------------------------------------------
    pvice(:ncol,:) = 0.

! stokes terminal velocity
    if (stokes) then

! get effective radius
       call reitab(ncol, t, rei)

       do k = 1,pver
          do i = 1,ncol
             if (cloud(i,k) > 0. .and. cldice(i,k) > 0.) then
                if (rei(i,k) < 40. ) then
                   vfall = 2./9. * rhoh2o * gravit * rei(i,k)**2 / eta  * 1.e-12  ! micons^2 -> m^2
                else
                   vfall = v40 + vslope * (rei(i,k)-r40)      ! linear above 40 microns
                end if
!!$                print *, t(i,k), rei(i,k), vfall*100.
! convert the fall speed to pressure units, but do not apply the traditional
! negative convention for pvel.
                pvice(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
             end if
          end do
       end do

       return
    end if

!-----------------------------------------------------------------------
!--------------------- ice fall velocity -------------------------------
!--------------------- McFarquhar and Heymsfield > icritc --------------
!-----------------------------------------------------------------------

! lower bound for iciwc

    cc = 128.64 
    bc = 53.242
    ac = 5.4795
    lbound = (-bc + sqrt(bc*bc-4*ac*cc))/(2*ac)
    lbound = 10.**lbound

    do k = 1,pver
       do i = 1,ncol
          if (cloud(i,k) > 0. .and. cldice(i,k) > 0.) then

! compute the in-cloud ice concentration (kg/kg)
             icice = cldice(i,k) / cloud(i,k)

! compute the ice water content in g/m3
             iciwc = icice * rho(i,k) * 1.e3

! set the fall velocity (cm/s) to depend on the ice water content in g/m3,
             if (iciwc > lbound) then ! need this because of log10
                logiwc = log10(iciwc)
!          Median - 
                vfall = 128.64 + 53.242*logiwc + 5.4795*logiwc**2
!          Average - 
!!$             vfall = 122.63 + 44.111*logiwc + 4.2144*logiwc**2
             else
                vfall = 0.
             end if

! set ice velocity to 1 cm/s if ice mixing ratio < icritc, ramp to value
! calculated above at 2*icritc
             if (icice <= icritc) then
                vfall = vice_small
             else if(icice < 2*icritc) then
                icefac = (icice-icritc)/icritc
                vfall = vice_small * (1.-icefac) + vfall * icefac
             end if

! bound the terminal velocity of ice particles at high concentration
             vfall = min(100.0_r8, vfall)

! convert the fall speed to pressure units, but do not apply the traditional
! negative convention for pvel.
             pvice(i,k+1) = vfall     &
                  * 0.01                 & ! cm to meters
                  * rho(i,k)*gravit        ! meters/sec to pascals/sec
          end if
       end do
    end do

    return
  end subroutine cld_sediment_vel


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

  subroutine cld_sediment_tend (ncol, dtime  ,               & 1,2
       pint   , pmid   , pdel   , t      ,                   &
       cloud  , cldliq , cldice , pvliq  , pvice  ,          &
       liqtend, icetend, wvtend , htend  , sfliq  , sfice   )

!----------------------------------------------------------------------
!     Apply Cloud Particle Gravitational Sedimentation to Condensate
!----------------------------------------------------------------------

    implicit none

! Arguments
    integer,  intent(in)  :: ncol                      ! number of colums to process

    real(r8), intent(in)  :: dtime                     ! time step
    real(r8), intent(in)  :: pint  (pcols,pverp)       ! interfaces pressure (Pa)
    real(r8), intent(in)  :: pmid  (pcols,pver)        ! midpoint pressures (Pa)
    real(r8), intent(in)  :: pdel  (pcols,pver)        ! pressure diff across layer (Pa)
    real(r8), intent(in)  :: cloud (pcols,pver)        ! cloud fraction (fraction)
    real(r8), intent(in)  :: t     (pcols,pver)        ! temperature (K)
    real(r8), intent(in)  :: cldliq(pcols,pver)        ! cloud liquid water (kg/kg)
    real(r8), intent(in)  :: cldice(pcols,pver)        ! cloud ice water    (kg/kg)
    real(r8), intent(in)  :: pvliq (pcols,pverp)       ! vertical velocity of liquid drops  (Pa/s)
    real(r8), intent(in)  :: pvice (pcols,pverp)       ! vertical velocity of ice particles (Pa/s)
! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))

    real(r8), intent(out) :: liqtend(pcols,pver)       ! liquid condensate tend
    real(r8), intent(out) :: icetend(pcols,pver)       ! ice condensate tend
    real(r8), intent(out) :: wvtend (pcols,pver)       ! water vapor tend
    real(r8), intent(out) :: htend  (pcols,pver)       ! heating rate
    real(r8), intent(out) :: sfliq  (pcols)            ! surface flux of liquid (rain, kg/m/s)
    real(r8), intent(out) :: sfice  (pcols)            ! surface flux of ice    (snow, kg/m/s)

! Local variables
    real(r8) :: fxliq(pcols,pverp)                     ! fluxes at the interfaces, liquid (positive = down)
    real(r8) :: fxice(pcols,pverp)                     ! fluxes at the interfaces, ice    (positive = down)
    real(r8) :: cldab(pcols)                           ! cloud in layer above
    real(r8) :: evapliq                                ! evaporation of cloud liquid into environment
    real(r8) :: evapice                                ! evaporation of cloud ice into environment
    real(r8) :: cldovrl                                ! cloud overlap factor

    integer :: i,k
!----------------------------------------------------------------------

! initialize variables
    fxliq  (:ncol,:) = 0. ! flux at interfaces (liquid)
    fxice  (:ncol,:) = 0. ! flux at interfaces (ice)
    liqtend(:ncol,:) = 0. ! condensate tend (liquid)
    icetend(:ncol,:) = 0. ! condensate tend (ice)
    wvtend(:ncol,:)  = 0. ! environmental moistening
    htend(:ncol,:)   = 0. ! evaporative cooling
    sfliq(:ncol)     = 0. ! condensate sedimentation flux out bot of column (liquid)
    sfice(:ncol)     = 0. ! condensate sedimentation flux out bot of column (ice)

! fluxes at interior points
    call getflx(ncol, pint, cldliq, pvliq, dtime, fxliq)
    call getflx(ncol, pint, cldice, pvice, dtime, fxice)

! calculate fluxes at boundaries
    do i = 1,ncol
       fxliq(i,1) = 0
       fxice(i,1) = 0
! surface flux by upstream scheme
       fxliq(i,pverp) = cldliq(i,pver) * pvliq(i,pverp) * dtime
       fxice(i,pverp) = cldice(i,pver) * pvice(i,pverp) * dtime
    end do

! filter out any negative fluxes from the getflx routine
! (typical fluxes are of order > 1e-3 when clouds are present)
    do k = 2,pver
       fxliq(:ncol,k) = max(0._r8, fxliq(:ncol,k))
       fxice(:ncol,k) = max(0._r8, fxice(:ncol,k))
    end do

! Limit the flux out of the bottom of each cell to the water content in each phase.
! Apply mxsedfac to prevent generating very small negative cloud water/ice
! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
! ***Should we include the flux in the top, to allow for thin surface layers?
! ***Requires simple treatment of cloud overlap, already included below.
    do k = 1,pver
       do i = 1,ncol
          fxliq(i,k+1) = min( fxliq(i,k+1), mxsedfac * cldliq(i,k) * pdel(i,k) )
          fxice(i,k+1) = min( fxice(i,k+1), mxsedfac * cldice(i,k) * pdel(i,k) )
!!$        fxliq(i,k+1) = min( fxliq(i,k+1), cldliq(i,k) * pdel(i,k) + fxliq(i,k))
!!$        fxice(i,k+1) = min( fxice(i,k+1), cldice(i,k) * pdel(i,k) + fxice(i,k))
!!$        fxliq(i,k+1) = min( fxliq(i,k+1), cloud(i,k) * cldliq(i,k) * pdel(i,k) )
!!$        fxice(i,k+1) = min( fxice(i,k+1), cloud(i,k) * cldice(i,k) * pdel(i,k) )
       end do
    end do

! Now calculate the tendencies assuming that condensate evaporates when
! it falls into environment, and does not when it falls into cloud.
! All flux out of the layer comes from the cloudy part.
! Assume maximum overlap for stratiform clouds
!  if cloud above < cloud,  all water falls into cloud below
!  if cloud above > cloud,  water split between cloud  and environment
    do k = 1,pver
       cldab(:ncol) = 0.
       do i = 1,ncol
! cloud overlap cloud factor
          cldovrl  = min( cloud(i,k) / (cldab(i)+.0001), 1._r8 )
          cldab(i) = cloud(i,k)
! evaporation into environment cause moistening and cooling
          evapliq = fxliq(i,k) * (1.-cldovrl) / (dtime * pdel(i,k))  ! into env (kg/kg/s)
          evapice = fxice(i,k) * (1.-cldovrl) / (dtime * pdel(i,k))  ! into env (kg/kg/s)
          wvtend(i,k) = evapliq + evapice                          ! evaporation into environment (kg/kg/s)
          htend (i,k) = -latvap*evapliq -(latvap+latice)*evapice   ! evaporation (W/kg)
! net flux into cloud changes cloud liquid/ice (all flux is out of cloud)
          liqtend(i,k)  = (fxliq(i,k)*cldovrl - fxliq(i,k+1)) / (dtime * pdel(i,k))
          icetend(i,k)  = (fxice(i,k)*cldovrl - fxice(i,k+1)) / (dtime * pdel(i,k))
       end do
    end do

! convert flux out the bottom to mass units Pa -> kg/m2/s
    sfliq(:ncol) = fxliq(:ncol,pverp) / (dtime*gravit)
    sfice(:ncol) = fxice(:ncol,pverp) / (dtime*gravit)

    return
  end subroutine cld_sediment_tend

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

  subroutine getflx(ncol, xw, phi, vel, deltat, flux) 3,6

!.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
!....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
!....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
!.........phi1......phi2.......phi3.....phi4.......phi5.......


    implicit none

    integer, intent(in) :: ncol                      ! number of colums to process

    integer i
    integer k

    real (r8), intent(in) :: vel(pcols,pverp)
    real (r8) flux(pcols,pverp)
    real (r8) xw(pcols,pverp)
    real (r8) psi(pcols,pverp)
    real (r8), intent(in) :: phi(pcols,pverp-1)
    real (r8) fdot(pcols,pverp)
    real (r8) xx(pcols)
    real (r8) fxdot(pcols)
    real (r8) fxdd(pcols)

    real (r8) psistar(pcols)
    real (r8) deltat

    real (r8) xxk(pcols,pver)

    do i = 1,ncol
!        integral of phi
       psi(i,1) = 0.
!        fluxes at boundaries
       flux(i,1) = 0
       flux(i,pverp) = 0.
    end do

!     integral function
    do k = 2,pverp
       do i = 1,ncol
          psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
       end do
    end do


!     calculate the derivatives for the interpolating polynomial
    call cfdotmc_pro (ncol, xw, psi, fdot)

!  NEW WAY
!     calculate fluxes at interior pts
    do k = 2,pver
       do i = 1,ncol
          xxk(i,k) = xw(i,k)-vel(i,k)*deltat
       end do
    end do
    do k = 2,pver
       call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
       do i = 1,ncol
          flux(i,k) = (psi(i,k)-psistar(i))
       end do
    end do


    return
  end subroutine getflx



!##############################################################################


  subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar) 2,1


    implicit none

! input
    integer ncol                      ! number of colums to process

    real (r8) x(pcols, pverp)
    real (r8) f(pcols, pverp)
    real (r8) fdot(pcols, pverp)
    real (r8) xin(pcols)

! output
    real (r8) fxdot(pcols)
    real (r8) fxdd(pcols)
    real (r8) psistar(pcols)

    integer i
    integer k
    integer intz(pcols)
    real (r8) dx
    real (r8) s
    real (r8) c2
    real (r8) c3
    real (r8) xx
    real (r8) xinf
    real (r8) psi1, psi2, psi3, psim
    real (r8) cfint
    real (r8) cfnew
    real (r8) xins(pcols)

!     the minmod function 
    real (r8) a, b, c
    real (r8) minmod
    real (r8) medan
    logical found_error

    minmod(a,b) = 0.5*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
    medan(a,b,c) = a + minmod(b-a,c-a)

    do i = 1,ncol
       xins(i) = medan(x(i,1), xin(i), x(i,pverp))
       intz(i) = 0
    end do

! first find the interval 
    do k =  1,pverp-1
       do i = 1,ncol
          if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0) then
             intz(i) = k
          endif
       end do
    end do

    found_error=.false.
    do i = 1,ncol
       if (intz(i).eq.0.) found_error=.true.
    end do
    if(found_error) then
       do i = 1,ncol
          if (intz(i).eq.0.) then
             write (6,*) ' interval was not found for col i ', i
             call endrun('CFINT2')
          endif
       end do
    endif

! now interpolate
    do i = 1,ncol
       k = intz(i)
       dx = (x(i,k+1)-x(i,k))
       s = (f(i,k+1)-f(i,k))/dx
       c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
       c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
       xx = (xins(i)-x(i,k))
       fxdot(i) =  (3*c3*xx + 2*c2)*xx + fdot(i,k)
       fxdd(i) = 6*c3*xx + 2*c2
       cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)

!        limit the interpolant
       psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
       if (k.eq.1) then
          psi2 = f(i,1)
       else
          psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
       endif
       if (k+1.eq.pverp) then
          psi3 = f(i,pverp)
       else
          psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
       endif
       psim = medan(psi1, psi2, psi3)
       cfnew = medan(cfint, psi1, psim)
       if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36)  .gt..03) then
!     CHANGE THIS BACK LATER!!!
!     $        .gt..1) then


!     UNCOMMENT THIS LATER!!!
!            write (6,*) ' cfint2 limiting important ', cfint, cfnew


       endif
       psistar(i) = cfnew
    end do

    return
  end subroutine cfint2



!##############################################################################


  subroutine cfdotmc_pro (ncol, x, f, fdot) 2

!     prototype version; eventually replace with final SPITFIRE scheme

!     calculate the derivative for the interpolating polynomial
!     multi column version


    implicit none

! input
    integer ncol                      ! number of colums to process

    real (r8) x(pcols, pverp)
    real (r8) f(pcols, pverp)
! output
    real (r8) fdot(pcols, pverp)          ! derivative at nodes

! assumed variable distribution
!     x1.......x2.......x3.......x4.......x5.......x6     1,pverp points
!     f1.......f2.......f3.......f4.......f5.......f6     1,pverp points
!     ...sh1.......sh2......sh3......sh4......sh5....     1,pver points
!     .........d2.......d3.......d4.......d5.........     2,pver points
!     .........s2.......s3.......s4.......s5.........     2,pver points
!     .............dh2......dh3......dh4.............     2,pver-1 points
!     .............eh2......eh3......eh4.............     2,pver-1 points
!     ..................e3.......e4..................     3,pver-1 points
!     .................ppl3......ppl4................     3,pver-1 points
!     .................ppr3......ppr4................     3,pver-1 points
!     .................t3........t4..................     3,pver-1 points
!     ................fdot3.....fdot4................     3,pver-1 points


! work variables


    integer i
    integer k

    real (r8) a                    ! work var
    real (r8) b                    ! work var
    real (r8) c                    ! work var
    real (r8) s(pcols,pverp)             ! first divided differences at nodes
    real (r8) sh(pcols,pverp)            ! first divided differences between nodes
    real (r8) d(pcols,pverp)             ! second divided differences at nodes
    real (r8) dh(pcols,pverp)            ! second divided differences between nodes
    real (r8) e(pcols,pverp)             ! third divided differences at nodes
    real (r8) eh(pcols,pverp)            ! third divided differences between nodes
    real (r8) pp                   ! p prime
    real (r8) ppl(pcols,pverp)           ! p prime on left
    real (r8) ppr(pcols,pverp)           ! p prime on right
    real (r8) qpl
    real (r8) qpr
    real (r8) ttt
    real (r8) t
    real (r8) tmin
    real (r8) tmax
    real (r8) delxh(pcols,pverp)


!     the minmod function 
    real (r8) minmod
    real (r8) medan
    minmod(a,b) = 0.5*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
    medan(a,b,c) = a + minmod(b-a,c-a)

    do k = 1,pver


!        first divided differences between nodes
       do i = 1, ncol
          delxh(i,k) = (x(i,k+1)-x(i,k))
          sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
       end do

!        first and second divided differences at nodes
       if (k.ge.2) then
          do i = 1,ncol
             d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
             s(i,k) = minmod(sh(i,k),sh(i,k-1))
          end do
       endif
    end do

!     second and third divided diffs between nodes
    do k = 2,pver-1
       do i = 1, ncol
          eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
          dh(i,k) = minmod(d(i,k),d(i,k+1))
       end do
    end do

!     treat the boundaries
    do i = 1,ncol
       e(i,2) = eh(i,2)
       e(i,pver) = eh(i,pver-1)
!        outside level
       fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1)  &
            - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
       fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
       fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver)  &
            + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
       fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
!        one in from boundary
       fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
       fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
       fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver)   &
            - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
       fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
    end do


    do k = 3,pver-1
       do i = 1,ncol
          e(i,k) = minmod(eh(i,k),eh(i,k-1))
       end do
    end do



    do k = 3,pver-1

       do i = 1,ncol

!           p prime at k-0.5
          ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)  
!           p prime at k+0.5
          ppr(i,k)=sh(i,k)   - dh(i,k)  *delxh(i,k)

          t = minmod(ppl(i,k),ppr(i,k))

!           derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
          pp = sh(i,k-1) + d(i,k)*delxh(i,k-1) 

!           quartic estimate of fdot
          fdot(i,k) = pp                            &
               - delxh(i,k-1)*delxh(i,k)            &
               *(  eh(i,k-1)*(x(i,k+2)-x(i,k  ))    &
               + eh(i,k  )*(x(i,k  )-x(i,k-2))      &
               )/(x(i,k+2)-x(i,k-2))

!           now limit it
          qpl = sh(i,k-1)       &
               + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
               d(i,k)  -e(i,k)*delxh(i,k))
          qpr = sh(i,k)         &
               + delxh(i,k  )*minmod(d(i,k)  +e(i,k)*delxh(i,k-1),        &
               d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))

          fdot(i,k) = medan(fdot(i,k), qpl, qpr)

          ttt = minmod(qpl, qpr)
          tmin = min(0._r8,3*s(i,k),1.5*t,ttt)
          tmax = max(0._r8,3*s(i,k),1.5*t,ttt)

          fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))

       end do

    end do

    return
  end subroutine cfdotmc_pro
end module pkg_cld_sediment