module CanopyHydrologyMod

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Calculation of
  ! (1) water storage of intercepted precipitation
  ! (2) direct throughfall and canopy drainage of precipitation
  ! (3) the fraction of foliage covered by water and the fraction
  !     of foliage that is dry and transpiring.
  ! (4) snow layer initialization if the snow accumulation exceeds 10 mm.
  !
  ! !USES:
  use shr_kind_mod    , only : r8 => shr_kind_r8
  use shr_log_mod     , only : errMsg => shr_log_errMsg
  use shr_sys_mod     , only : shr_sys_flush
  use decompMod       , only : bounds_type
  use abortutils      , only : endrun
  use clm_varctl      , only : iulog
  use LandunitType    , only : lun                
  use atm2lndType     , only : atm2lnd_type
  use AerosolMod      , only : aerosol_type
  use CanopyStateType , only : canopystate_type
  use TemperatureType , only : temperature_type
  use WaterfluxType   , only : waterflux_type
  use WaterstateType  , only : waterstate_type
  use IrrigationMod   , only : irrigation_type
  use ColumnType      , only : col                
  use PatchType       , only : patch                
  !
  ! !PUBLIC TYPES:
  implicit none
  save
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  public :: CanopyHydrology_readnl ! Read namelist
  public :: CanopyHydrology        ! Run
  public :: IsSnowvegFlagOff       ! Returns true if snowveg_flag is OFF
  public :: IsSnowvegFlagOn        ! Returns true if snowveg_flag is ON
  public :: IsSnowvegFlagOnRad     ! Returns true if snowveg_flag is ON_RAD
  !
  ! !PRIVATE MEMBER FUNCTIONS:
  private :: FracWet    ! Determine fraction of vegetated surface that is wet
  private :: FracH2oSfc ! Determine fraction of land surfaces which are submerged  
  !
  ! !PRIVATE DATA MEMBERS:
  integer :: oldfflag=0  ! use old fsno parameterization (N&Y07) 
  real(r8) :: interception_fraction ! Fraction of intercepted precipitation
  real(r8) :: maximum_leaf_wetted_fraction ! Maximum fraction of leaf that may be wet
  logical, private :: use_clm5_fpi    = .false. ! use clm5 fpi equation
  ! Snow in vegetation canopy namelist options.
  logical, private :: snowveg_off    = .false.  ! snowveg_flag = 'OFF'
  logical, private :: snowveg_on     = .false.  ! snowveg_flag = 'ON'
  logical, private :: snowveg_onrad  = .true.   ! snowveg_flag = 'ON_RAD'
  ! for now, all mods on by default:
  character(len= 10), public           :: snowveg_flag = 'ON_RAD'

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

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

contains

  !-----------------------------------------------------------------------
  subroutine CanopyHydrology_readnl( NLFilename )
    !
    ! !DESCRIPTION:
    ! Read the namelist for CanopyHydrology
    !
    ! !USES:
    use spmdMod       , only : masterproc, mpicom
    use fileutils     , only : getavu, relavu, opnfil
    use shr_nl_mod    , only : shr_nl_find_group_name
    use shr_mpi_mod   , only : shr_mpi_bcast
    !
    ! !ARGUMENTS:
    character(len=*), intent(IN) :: NLFilename ! Namelist filename
    !
    ! !LOCAL VARIABLES:
    integer :: ierr                 ! error code
    integer :: unitn                ! unit for namelist file
    character(len=32) :: subname = 'CanopyHydrology_readnl'  ! subroutine name
    !-----------------------------------------------------------------------

    namelist /clm_canopyhydrology_inparm/ &
         oldfflag, &
         interception_fraction, &
         maximum_leaf_wetted_fraction, &
         use_clm5_fpi, &
         snowveg_flag

    ! ----------------------------------------------------------------------
    ! Read namelist from standard input. 
    ! ----------------------------------------------------------------------

    if ( masterproc )then

       unitn = getavu()
       write(iulog,*) 'Read in clm_CanopyHydrology_inparm  namelist'
       call opnfil (NLFilename, unitn, 'F')
       call shr_nl_find_group_name(unitn, 'clm_canopyhydrology_inparm', status=ierr)
       if (ierr == 0) then
          read(unitn, clm_canopyhydrology_inparm, iostat=ierr)
          if (ierr /= 0) then
             call endrun(msg="ERROR reading clm_canopyhydrology_inparm namelist"//errmsg(sourcefile, __LINE__))
          end if
       else
          call endrun(msg="ERROR finding clm_canopyhydrology_inparm namelist"//errmsg(sourcefile, __LINE__))
       end if
       call relavu( unitn )

       snowveg_off   = IsSnowvegFlagOff()
       snowveg_on    = IsSnowvegFlagOn()
       snowveg_onrad = IsSnowvegFlagOnRad()
       write(iulog,*) 'snowveg_off = ',snowveg_off
       write(iulog,*) 'snowveg_on = ',snowveg_on
       write(iulog,*) 'snowveg_onrad = ',snowveg_onrad
       if (snowveg_off .or. snowveg_on .or. snowveg_onrad) then
          write(iulog,*) 'snowveg_flag = ',snowveg_flag
       else
          call endrun(msg="snowveg_flag is set incorrectly (not ON, ON_RAD, or OFF)"//errmsg(sourcefile, __LINE__))
       end if
    end if
    ! Broadcast namelist variables read in
    call shr_mpi_bcast(oldfflag, mpicom)
    call shr_mpi_bcast(interception_fraction, mpicom)
    call shr_mpi_bcast(maximum_leaf_wetted_fraction, mpicom)
    call shr_mpi_bcast(use_clm5_fpi, mpicom)
    call shr_mpi_bcast(snowveg_flag, mpicom)

    if (masterproc) then
       write(iulog,*) ' '
       write(iulog,*) 'canopyhydrology settings:'
       write(iulog,*) '  interception_fraction        = ',interception_fraction
       write(iulog,*) '  maximum_leaf_wetted_fraction = ',maximum_leaf_wetted_fraction
       write(iulog,*) '  use_clm5_fpi                 = ',use_clm5_fpi
    endif

   end subroutine CanopyHydrology_readnl

   !-----------------------------------------------------------------------
   subroutine CanopyHydrology(bounds, &
        num_nolakec, filter_nolakec, num_nolakep, filter_nolakep, &
        atm2lnd_inst, canopystate_inst, temperature_inst, &
        aerosol_inst, waterstate_inst, waterflux_inst, irrigation_inst)
     !
     ! !DESCRIPTION:
     ! Calculation of
     ! (1) water storage of intercepted precipitation
     ! (2) direct throughfall and canopy drainage of precipitation
     ! (3) the fraction of foliage covered by water and the fraction
     !     of foliage that is dry and transpiring.
     ! (4) snow layer initialization if the snow accumulation exceeds 10 mm.
     ! Note:  The evaporation loss is taken off after the calculation of leaf
     ! temperature in the subroutine clm\_leaftem.f90, not in this subroutine.
     !
     ! !USES:
     use clm_varcon         , only : hfus, denice, zlnd, rpi, spval, tfrz
     use column_varcon      , only : icol_roof, icol_sunwall, icol_shadewall
     use landunit_varcon    , only : istcrop, istice, istwet, istsoil, istice_mec 
     use clm_varctl         , only : subgridflag
     use clm_varpar         , only : nlevsoi,nlevsno
     use clm_time_manager   , only : get_step_size
     use subgridAveMod      , only : p2c
     use SnowHydrologyMod   , only : NewSnowBulkDensity
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds     
     integer                , intent(in)    :: num_nolakec          ! number of column non-lake points in column filter
     integer                , intent(in)    :: filter_nolakec(:)    ! column filter for non-lake points
     integer                , intent(in)    :: num_nolakep          ! number of pft non-lake points in pft filter
     integer                , intent(in)    :: filter_nolakep(:)    ! patch filter for non-lake points
     type(atm2lnd_type)     , intent(in)    :: atm2lnd_inst
     type(canopystate_type) , intent(in)    :: canopystate_inst
     type(temperature_type) , intent(inout) :: temperature_inst
     type(aerosol_type)     , intent(inout) :: aerosol_inst
     type(waterstate_type)  , intent(inout) :: waterstate_inst
     type(waterflux_type)   , intent(inout) :: waterflux_inst
     type(irrigation_type)  , intent(in)    :: irrigation_inst
     !
     ! !LOCAL VARIABLES:
     integer  :: f                                            ! filter index
     integer  :: pi                                           ! patch index
     integer  :: p                                            ! patch index
     integer  :: c                                            ! column index
     integer  :: l                                            ! landunit index
     integer  :: g                                            ! gridcell index
     integer  :: newnode                                      ! flag when new snow node is set, (1=yes, 0=no)
     real(r8) :: dtime                                        ! land model time step (sec)
     real(r8) :: h2ocanmx                                     ! maximum allowed water on canopy [mm]
     real(r8) :: snocanmx                                     ! maximum allowed snow on canopy [mm equiv. water]
     real(r8) :: liqcanmx                                     ! maximum allowed snowliquid water on canopy [mm equiv. water]
     real(r8) :: xsnorun                                      ! excess snow that exceeds the leaf capacity [mm/s]
     real(r8) :: xliqrun                                      ! excess liquid water 
     real(r8) :: qflx_snocanfall(bounds%begp:bounds%endp)     ! rate of excess canopy snow falling off canopy
     real(r8) :: qflx_liqcanfall(bounds%begp:bounds%endp)     ! rate of excess canopy liquid falling off canopy
     real(r8) :: fpi                                          ! coefficient of interception
     real(r8) :: fpisnow                                      ! coefficient of interception for snowfall
     real(r8) :: xrun                                         ! excess water that exceeds the leaf capacity [mm/s]
     real(r8) :: dz_snowf                                     ! layer thickness rate change due to precipitation [mm/s]
     real(r8) :: bifall(bounds%begc:bounds%endc)              ! bulk density of newly fallen dry snow [kg/m3]
     real(r8) :: fracsnow(bounds%begp:bounds%endp)            ! frac of precipitation that is snow
     real(r8) :: fracrain(bounds%begp:bounds%endp)            ! frac of precipitation that is rain
     real(r8) :: qflx_candrip(bounds%begp:bounds%endp)        ! rate of canopy runoff and snow falling off canopy [mm/s]
     real(r8) :: qflx_through_rain(bounds%begp:bounds%endp)   ! direct rain throughfall [mm/s]
     real(r8) :: qflx_through_snow(bounds%begp:bounds%endp)   ! direct snow throughfall [mm/s]
     real(r8) :: qflx_prec_grnd_snow(bounds%begp:bounds%endp) ! snow precipitation incident on ground [mm/s]
     real(r8) :: qflx_prec_grnd_rain(bounds%begp:bounds%endp) ! rain precipitation incident on ground [mm/s]
     real(r8) :: z_avg                                        ! grid cell average snow depth
     real(r8) :: rho_avg                                      ! avg density of snow column
     real(r8) :: temp_snow_depth,temp_intsnow                 ! temporary variables
     real(r8) :: fmelt
     real(r8) :: smr
     real(r8) :: delf_melt
     real(r8) :: fsno_new
     real(r8) :: accum_factor
     real(r8) :: newsnow(bounds%begc:bounds%endc)
     real(r8) :: snowmelt(bounds%begc:bounds%endc)
     integer  :: j
     !-----------------------------------------------------------------------

     associate(                                                             & 
          snl                  => col%snl                                 , & ! Input:  [integer  (:)   ]  number of snow layers                    
          n_melt               => col%n_melt                              , & ! Input:  [real(r8) (:)   ]  SCA shape parameter                     
          zi                   => col%zi                                  , & ! Output: [real(r8) (:,:) ]  interface level below a "z" level (m) 
          dz                   => col%dz                                  , & ! Output: [real(r8) (:,:) ]  layer depth (m)                       
          z                    => col%z                                   , & ! Output: [real(r8) (:,:) ]  layer thickness (m)                   

          forc_rain            => atm2lnd_inst%forc_rain_downscaled_col   , & ! Input:  [real(r8) (:)   ]  rain rate [mm/s]                        
          forc_snow            => atm2lnd_inst%forc_snow_downscaled_col   , & ! Input:  [real(r8) (:)   ]  snow rate [mm/s]                        
          forc_t               => atm2lnd_inst%forc_t_downscaled_col      , & ! Input:  [real(r8) (:)   ]  atmospheric temperature (Kelvin)        
          qflx_floodg          => atm2lnd_inst%forc_flood_grc             , & ! Input:  [real(r8) (:)   ]  gridcell flux of flood water from RTM   
          forc_wind            => atm2lnd_inst%forc_wind_grc              , & ! Input:  [real(r8) (:)   ]  atmospheric wind speed (m/s)
          dewmx                => canopystate_inst%dewmx_patch            , & ! Input:  [real(r8) (:)   ]  Maximum allowed dew [mm]                
          frac_veg_nosno       => canopystate_inst%frac_veg_nosno_patch   , & ! Input:  [integer  (:)   ]  fraction of vegetation not covered by snow (0 OR 1) [-]
          elai                 => canopystate_inst%elai_patch             , & ! Input:  [real(r8) (:)   ]  one-sided leaf area index with burying by snow
          esai                 => canopystate_inst%esai_patch             , & ! Input:  [real(r8) (:)   ]  one-sided stem area index with burying by snow

          t_grnd               => temperature_inst%t_grnd_col             , & ! Input:  [real(r8) (:)   ]  ground temperature (Kelvin)             
          t_soisno             => temperature_inst%t_soisno_col           , & ! Output: [real(r8) (:,:) ]  soil temperature (Kelvin)  

          h2ocan               => waterstate_inst%h2ocan_patch            , & ! Output: [real(r8) (:)   ]  total canopy water (mm H2O)             
          snocan               => waterstate_inst%snocan_patch            , & ! Output: [real(r8) (:)   ]  canopy snow (mm H2O)             
          liqcan               => waterstate_inst%liqcan_patch            , & ! Output: [real(r8) (:)   ]  canopy liquid (mm H2O)   
          snounload            => waterstate_inst%snounload_patch         , & ! Output: [real(r8) (:)   ]  canopy snow unloading (mm H2O)
          h2osfc               => waterstate_inst%h2osfc_col              , & ! Output: [real(r8) (:)   ]  surface water (mm)                      
          h2osno               => waterstate_inst%h2osno_col              , & ! Output: [real(r8) (:)   ]  snow water (mm H2O)                     
          snow_depth           => waterstate_inst%snow_depth_col          , & ! Output: [real(r8) (:)   ]  snow height (m)                         
          int_snow             => waterstate_inst%int_snow_col            , & ! Output: [real(r8) (:)   ]  integrated snowfall [mm]                
          frac_sno_eff         => waterstate_inst%frac_sno_eff_col        , & ! Output: [real(r8) (:)   ]  eff. fraction of ground covered by snow (0 to 1)
          frac_sno             => waterstate_inst%frac_sno_col            , & ! Output: [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)
          frac_h2osfc          => waterstate_inst%frac_h2osfc_col         , & ! Output: [real(r8) (:)   ]  fraction of ground covered by surface water (0 to 1)
          frac_iceold          => waterstate_inst%frac_iceold_col         , & ! Output: [real(r8) (:,:) ]  fraction of ice relative to the tot water
          h2osoi_ice           => waterstate_inst%h2osoi_ice_col          , & ! Output: [real(r8) (:,:) ]  ice lens (kg/m2)                      
          h2osoi_liq           => waterstate_inst%h2osoi_liq_col          , & ! Output: [real(r8) (:,:) ]  liquid water (kg/m2)                  
          swe_old              => waterstate_inst%swe_old_col             , & ! Output: [real(r8) (:,:) ]  snow water before update              

          qflx_floodc          => waterflux_inst%qflx_floodc_col          , & ! Output: [real(r8) (:)   ]  column flux of flood water from RTM     
          qflx_snow_drain       => waterflux_inst%qflx_snow_drain_col     , & ! Input: [real(r8) (:)   ]  drainage from snow pack from previous time step       
          qflx_snow_h2osfc     => waterflux_inst%qflx_snow_h2osfc_col     , & ! Output: [real(r8) (:)   ]  snow falling on surface water (mm/s)     
          qflx_snow_grnd_col   => waterflux_inst%qflx_snow_grnd_col       , & ! Output: [real(r8) (:)   ]  snow on ground after interception (mm H2O/s) [+]
          qflx_snow_grnd_patch => waterflux_inst%qflx_snow_grnd_patch     , & ! Output: [real(r8) (:)   ]  snow on ground after interception (mm H2O/s) [+]
          qflx_prec_intr       => waterflux_inst%qflx_prec_intr_patch     , & ! Output: [real(r8) (:)   ]  interception of precipitation [mm/s]    
          qflx_prec_grnd       => waterflux_inst%qflx_prec_grnd_patch     , & ! Output: [real(r8) (:)   ]  water onto ground including canopy runoff [kg/(m2 s)]
          qflx_rain_grnd       => waterflux_inst%qflx_rain_grnd_patch     , & ! Output: [real(r8) (:)   ]  rain on ground after interception (mm H2O/s) [+]

          qflx_irrig           => irrigation_inst%qflx_irrig_patch        , & ! Input: [real(r8) (:)    ]  irrigation amount (mm/s)           
          qflx_snowindunload   => waterflux_inst%qflx_snowindunload_patch , & ! Output: [real(r8) (:)   ]  canopy snow unloading from wind [mm/s]    
          qflx_snotempunload   => waterflux_inst%qflx_snotempunload_patch   & ! Output: [real(r8) (:)   ]  canopy snow unloading from temp. [mm/s]    
          )

       ! Compute time step
       
       dtime = get_step_size()

       ! Set status of snowveg_flag
       snowveg_on    = IsSnowvegFlagOn()
       snowveg_onrad = IsSnowvegFlagOnRad()

       ! Start patch loop

       do f = 1, num_nolakep
          p = filter_nolakep(f)
          g = patch%gridcell(p)
          l = patch%landunit(p)
          c = patch%column(p)

          ! Canopy interception and precipitation onto ground surface
          ! Add precipitation to leaf water

          if (lun%itype(l)==istsoil .or. lun%itype(l)==istwet .or. lun%urbpoi(l) .or. &
               lun%itype(l)==istcrop) then

             qflx_candrip(p) = 0._r8      ! rate of canopy runoff
             qflx_snocanfall(p) = 0._r8      ! rate of just snow canopy fall
             qflx_liqcanfall(p) = 0._r8
             qflx_snowindunload(p) = 0._r8
             qflx_snotempunload(p) = 0._r8
             snounload(p)=0._r8
             qflx_through_snow(p) = 0._r8 ! rain precipitation direct through canopy
             qflx_through_rain(p) = 0._r8 ! snow precipitation direct through canopy
             qflx_prec_intr(p) = 0._r8    ! total intercepted precipitation
             fracsnow(p) = 0._r8          ! fraction of input precip that is snow
             fracrain(p) = 0._r8          ! fraction of input precip that is rain


             if (col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall) then

                if (frac_veg_nosno(p) == 1 .and. (forc_rain(c) + forc_snow(c)) > 0._r8) then

                   ! determine fraction of input precipitation that is snow and rain
                   fracsnow(p) = forc_snow(c)/(forc_snow(c) + forc_rain(c))
                   fracrain(p) = forc_rain(c)/(forc_snow(c) + forc_rain(c))

                   ! The leaf water capacities for solid and liquid are different,
                   ! generally double for snow, but these are of somewhat less
                   ! significance for the water budget because of lower evap. rate at
                   ! lower temperature.  Hence, it is reasonable to assume that
                   ! vegetation storage of solid water is the same as liquid water.
                   h2ocanmx = dewmx(p) * (elai(p) + esai(p))
                   ! Coefficient of interception
                   if(use_clm5_fpi) then 
                      fpi = interception_fraction * tanh(elai(p) + esai(p))
                   else
                      fpi = 0.25_r8*(1._r8 - exp(-0.5_r8*(elai(p) + esai(p))))
                   endif

                   if (snowveg_on .or. snowveg_onrad) then
                      snocanmx = 60._r8*dewmx(p) * (elai(p) + esai(p))  ! 6*(LAI+SAI)
                      liqcanmx = h2ocanmx

                      fpisnow = (1._r8 - exp(-0.5_r8*(elai(p) + esai(p))))  ! max interception of 1 
                      ! Direct throughfall
                      qflx_through_snow(p) = forc_snow(c) * (1._r8-fpisnow)
                   else
                      ! Direct throughfall
                      qflx_through_snow(p) = forc_snow(c) * (1._r8-fpi)
                   end if

                   ! Direct throughfall
                   qflx_through_rain(p) = forc_rain(c) * (1._r8-fpi)

                   if (snowveg_on .or. snowveg_onrad) then
                      ! Intercepted precipitation [mm/s]
                      qflx_prec_intr(p) = forc_snow(c)*fpisnow + forc_rain(c)*fpi
                      ! storage of intercepted snowfall, rain, and dew
                      snocan(p) = max(0._r8, snocan(p) + dtime*forc_snow(c)*fpisnow)    
                      liqcan(p) = max(0._r8, liqcan(p) + dtime*forc_rain(c)*fpi)
                   else
                      ! Intercepted precipitation [mm/s]
                      qflx_prec_intr(p) = (forc_snow(c) + forc_rain(c)) * fpi
                   end if

                   ! Water storage of intercepted precipitation and dew
                   h2ocan(p) = max(0._r8, h2ocan(p) + dtime*qflx_prec_intr(p))

                   ! Initialize rate of canopy runoff and snow falling off canopy
                   qflx_candrip(p) = 0._r8
                   qflx_snocanfall(p) = 0._r8
                   qflx_liqcanfall(p) = 0._r8
                   qflx_snowindunload(p) = 0._r8
                   qflx_snotempunload(p) = 0._r8
                   snounload(p)=0._r8

                   if (snowveg_on .or. snowveg_onrad) then
                      if (forc_t(c) > tfrz) then ! Above freezing (Use t_veg?)
                         xliqrun = (liqcan(p) - liqcanmx)/dtime
                         if (xliqrun > 0._r8) then
                            qflx_liqcanfall(p) = xliqrun
                            liqcan(p) = liqcanmx
                         end if
                      else ! Below freezing
                         xsnorun = (snocan(p) - snocanmx)/dtime
                         if (xsnorun > 0._r8) then ! exceeds snow capacity
                            qflx_snocanfall(p) = xsnorun
                            snocan(p) = snocanmx
                         end if
                      end if
                      qflx_candrip(p) = qflx_snocanfall(p) + qflx_liqcanfall(p)
                      h2ocan(p) = snocan(p) + liqcan(p)
                   else 
                      ! Excess water that exceeds the leaf capacity
                      xrun = (h2ocan(p) - h2ocanmx)/dtime

                      ! Test on maximum dew on leaf
                      ! Note if xrun > 0 then h2ocan must be at least h2ocanmx
                      if (xrun > 0._r8) then
                         qflx_candrip(p) = xrun
                         h2ocan(p) = h2ocanmx
                      end if
                   end if
                end if
             end if

          else if (lun%itype(l)==istice .or. lun%itype(l)==istice_mec) then

             h2ocan(p)            = 0._r8
             qflx_candrip(p)      = 0._r8
             qflx_through_snow(p) = 0._r8
             qflx_through_rain(p) = 0._r8
             qflx_prec_intr(p)    = 0._r8
             fracsnow(p)          = 0._r8
             fracrain(p)          = 0._r8
             snocan(p)            = 0._r8
             liqcan(p)            = 0._r8
             qflx_snocanfall(p)    = 0._r8
             qflx_liqcanfall(p)    = 0._r8
             qflx_snowindunload(p) = 0._r8 
             qflx_snotempunload(p) = 0._r8 
             snounload(p)=0._r8

          end if

          ! Precipitation onto ground (kg/(m2 s))

          if (col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall) then
             if (frac_veg_nosno(p) == 0) then
                qflx_prec_grnd_snow(p) = forc_snow(c)
                qflx_prec_grnd_rain(p) = forc_rain(c)
             else
                if (snowveg_on .or. snowveg_onrad) then
                   qflx_snowindunload(p)=0._r8 
                   qflx_snotempunload(p)=0._r8 
                   snounload(p)=0._r8
                   if (snocan(p) > 0._r8) then
                      qflx_snotempunload(p) = max(0._r8,snocan(p)*(forc_t(c)-270.15_r8)/1.87e5_r8) 
                      qflx_snowindunload(p) = 0.5_r8*snocan(p)*forc_wind(g)/1.56e5_r8 
                      snounload(p) = (qflx_snowindunload(p)+qflx_snotempunload(p))*dtime ! total canopy unloading in timestep
                      if ( snounload(p) > snocan(p) ) then ! Limit unloading to snow in canopy
                         snounload(p) = snocan(p)
                         write(iulog,"(A,I2.2,A,ES13.4E2)") "snocan",p,": ",snocan(p)
                      end if
                      snocan(p) = snocan(p) - snounload(p)
                      h2ocan(p) = h2ocan(p) - snounload(p)
                   endif
                   qflx_prec_grnd_snow(p) = qflx_through_snow(p) + qflx_snocanfall(p)  + snounload(p)/dtime
                   qflx_prec_grnd_rain(p) = qflx_through_rain(p) + qflx_liqcanfall(p) 

                else
                   qflx_prec_grnd_snow(p) = qflx_through_snow(p) + (qflx_candrip(p) * fracsnow(p))
                   qflx_prec_grnd_rain(p) = qflx_through_rain(p) + (qflx_candrip(p) * fracrain(p))
                end if
             end if
             ! Urban sunwall and shadewall have no intercepted precipitation
          else
             qflx_prec_grnd_snow(p) = 0.
             qflx_prec_grnd_rain(p) = 0.
          end if

          ! Add irrigation water directly onto ground (bypassing canopy interception)
          ! Note that it's still possible that (some of) this irrigation water will runoff (as runoff is computed later)
          qflx_prec_grnd_rain(p) = qflx_prec_grnd_rain(p) + qflx_irrig(p)

          qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)

          qflx_snow_grnd_patch(p) = qflx_prec_grnd_snow(p)           ! ice onto ground (mm/s)
          qflx_rain_grnd(p)       = qflx_prec_grnd_rain(p)           ! liquid water onto ground (mm/s)

       end do ! (end patch loop)

       ! Determine the fraction of foliage covered by water and the
       ! fraction of foliage that is dry and transpiring.

       call FracWet(num_nolakep, filter_nolakep, &
            canopystate_inst, waterstate_inst)

       ! Update column level state variables for snow.

       call p2c(bounds, num_nolakec, filter_nolakec, &
            qflx_snow_grnd_patch(bounds%begp:bounds%endp), &
            qflx_snow_grnd_col(bounds%begc:bounds%endc))

       ! apply gridcell flood water flux to non-lake columns
       do f = 1, num_nolakec
          c = filter_nolakec(f)
          g = col%gridcell(c)
          if (col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall) then      
             qflx_floodc(c) = qflx_floodg(g)
          else
             qflx_floodc(c) = 0._r8
          endif
       enddo

       ! Determine snow height and snow water

       call NewSnowBulkDensity(bounds, num_nolakec, filter_nolakec, &
            atm2lnd_inst, bifall(bounds%begc:bounds%endc))

       do f = 1, num_nolakec
          c = filter_nolakec(f)
          l = col%landunit(c)
          g = col%gridcell(c)

          ! Use Alta relationship, Anderson(1976); LaChapelle(1961),
          ! U.S.Department of Agriculture Forest Service, Project F,
          ! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.

          qflx_snow_h2osfc(c) = 0._r8
          ! set temporary variables prior to updating
          temp_snow_depth=snow_depth(c)
          ! save initial snow content
          do j= -nlevsno+1,snl(c)
             swe_old(c,j) = 0.0_r8
          end do
          do j= snl(c)+1,0
             swe_old(c,j)=h2osoi_liq(c,j)+h2osoi_ice(c,j)
          enddo

          ! newsnow is all snow that doesn't fall on h2osfc
          newsnow(c) = (1._r8 - frac_h2osfc(c)) * qflx_snow_grnd_col(c) * dtime

          ! update int_snow
          int_snow(c) = max(int_snow(c),h2osno(c)) !h2osno could be larger due to frost

          ! snowmelt from previous time step * dtime
          snowmelt(c) = qflx_snow_drain(c) * dtime

          ! set shape factor for accumulation of snow
          accum_factor=0.1

          if (h2osno(c) > 0.0) then

             !======================  FSCA PARAMETERIZATIONS  ======================
             ! fsca parameterization based on *changes* in swe
             ! first compute change from melt during previous time step
             if(snowmelt(c) > 0._r8) then

                smr=min(1._r8,(h2osno(c))/(int_snow(c)))

                frac_sno(c) = 1. - (acos(min(1._r8,(2.*smr - 1._r8)))/rpi)**(n_melt(c))

             endif

             ! update fsca by new snow event, add to previous fsca
             if (newsnow(c) > 0._r8) then
                fsno_new = 1._r8 - (1._r8 - tanh(accum_factor*newsnow(c)))*(1._r8 - frac_sno(c))
                frac_sno(c) = fsno_new

                ! reset int_snow after accumulation events
                temp_intsnow= (h2osno(c) + newsnow(c)) &
                     / (0.5*(cos(rpi*(1._r8-max(frac_sno(c),1e-6_r8))**(1./n_melt(c)))+1._r8))
                int_snow(c) = min(1.e8_r8,temp_intsnow)
             endif

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

             ! for subgrid fluxes
             if (subgridflag ==1 .and. .not. lun%urbpoi(l)) then
                if (frac_sno(c) > 0._r8)then
                   snow_depth(c)=snow_depth(c) + newsnow(c)/(bifall(c) * frac_sno(c))
                else
                   snow_depth(c)=0._r8
                end if
             else
                ! for uniform snow cover
                snow_depth(c)=snow_depth(c)+newsnow(c)/bifall(c)
             endif

             ! use original fsca formulation (n&y 07)
             if (oldfflag == 1) then 
                ! snow cover fraction in Niu et al. 2007
                if(snow_depth(c) > 0.0_r8)  then
                   frac_sno(c) = tanh(snow_depth(c)/(2.5_r8*zlnd* &
                        (min(800._r8,(h2osno(c)+ newsnow(c))/snow_depth(c))/100._r8)**1._r8) )
                endif
                if(h2osno(c) < 1.0_r8)  then
                   frac_sno(c)=min(frac_sno(c),h2osno(c))
                endif
             endif

          else !h2osno == 0
             ! initialize frac_sno and snow_depth when no snow present initially
             if (newsnow(c) > 0._r8) then 
                z_avg = newsnow(c)/bifall(c)
                fmelt=newsnow(c)
                frac_sno(c) = tanh(accum_factor*newsnow(c))

                ! make int_snow consistent w/ new fsno, h2osno
                int_snow(c) = 0. !reset prior to adding newsnow below
                temp_intsnow= (h2osno(c) + newsnow(c)) &
                     / (0.5*(cos(rpi*(1._r8-max(frac_sno(c),1e-6_r8))**(1./n_melt(c)))+1._r8))
                int_snow(c) = min(1.e8_r8,temp_intsnow)

                ! update snow_depth and h2osno to be consistent with frac_sno, z_avg
                if (subgridflag ==1 .and. .not. lun%urbpoi(l)) then
                   snow_depth(c)=z_avg/frac_sno(c)
                else
                   snow_depth(c)=newsnow(c)/bifall(c)
                endif
                ! use n&y07 formulation
                if (oldfflag == 1) then 
                   ! snow cover fraction in Niu et al. 2007
                   if(snow_depth(c) > 0.0_r8)  then
                      frac_sno(c) = tanh(snow_depth(c)/(2.5_r8*zlnd* &
                           (min(800._r8,newsnow(c)/snow_depth(c))/100._r8)**1._r8) )
                   endif
                endif
             else
                z_avg = 0._r8
                snow_depth(c) = 0._r8
                frac_sno(c) = 0._r8
             endif
          endif ! end of h2osno > 0

          ! snow directly falling on surface water melts, increases h2osfc
          qflx_snow_h2osfc(c) = frac_h2osfc(c)*qflx_snow_grnd_col(c)

          ! update h2osno for new snow
          h2osno(c) = h2osno(c) + newsnow(c) 
          int_snow(c) = int_snow(c) + newsnow(c)

          ! update change in snow depth
          dz_snowf = (snow_depth(c) - temp_snow_depth) / dtime

          ! set frac_sno_eff variable
          if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then
             if (subgridflag ==1) then 
                frac_sno_eff(c) = frac_sno(c)
             else
                frac_sno_eff(c) = 1._r8
             endif
          else
             frac_sno_eff(c) = 1._r8
          endif

          if (lun%itype(l)==istwet .and. t_grnd(c)>tfrz) then
             h2osno(c)=0._r8
             snow_depth(c)=0._r8
          end if

          ! When the snow accumulation exceeds 10 mm, initialize snow layer
          ! Currently, the water temperature for the precipitation is simply set
          ! as the surface air temperature

          newnode = 0    ! flag for when snow node will be initialized
          if (snl(c) == 0 .and. frac_sno(c)*snow_depth(c) >= 0.01_r8) then
             newnode = 1
             snl(c) = -1
             dz(c,0) = snow_depth(c)                       ! meter
             z(c,0) = -0.5_r8*dz(c,0)
             zi(c,-1) = -dz(c,0)
             t_soisno(c,0) = min(tfrz, forc_t(c))      ! K
             h2osoi_ice(c,0) = h2osno(c)               ! kg/m2
             h2osoi_liq(c,0) = 0._r8                   ! kg/m2
             frac_iceold(c,0) = 1._r8

             ! intitialize SNICAR variables for fresh snow:
             call aerosol_inst%Reset(column=c)
             call waterstate_inst%Reset(column=c)
          end if

          ! The change of ice partial density of surface node due to precipitation.
          ! Only ice part of snowfall is added here, the liquid part will be added
          ! later.

          if (snl(c) < 0 .and. newnode == 0) then
             h2osoi_ice(c,snl(c)+1) = h2osoi_ice(c,snl(c)+1)+newsnow(c)
             dz(c,snl(c)+1) = dz(c,snl(c)+1)+dz_snowf*dtime
          end if

       end do

       ! update surface water fraction (this may modify frac_sno)
       call FracH2oSfc(bounds, num_nolakec, filter_nolakec, &
            waterstate_inst)

     end associate 

   end subroutine CanopyHydrology

   !-----------------------------------------------------------------------
   subroutine FracWet(numf, filter, canopystate_inst, waterstate_inst)
     !
     ! !DESCRIPTION:
     ! Determine fraction of vegetated surfaces which are wet and
     ! fraction of elai which is dry. The variable ``fwet'' is the
     ! fraction of all vegetation surfaces which are wet including
     ! stem area which contribute to evaporation. The variable ``fdry''
     ! is the fraction of elai which is dry because only leaves
     ! can transpire.  Adjusted for stem area which does not transpire.
     !
     ! ! USES:
     use clm_varcon         , only : tfrz
     ! !ARGUMENTS:
     integer                , intent(in)    :: numf                  ! number of filter non-lake points
     integer                , intent(in)    :: filter(numf)          ! patch filter for non-lake points
     type(canopystate_type) , intent(in)    :: canopystate_inst
     type(waterstate_type)  , intent(inout) :: waterstate_inst
     !
     ! !LOCAL VARIABLES:
     integer  :: fp,p             ! indices
     real(r8) :: vegt             ! lsai
     real(r8) :: dewmxi           ! inverse of maximum allowed dew [1/mm]
     !-----------------------------------------------------------------------

     associate(                                              & 
          frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input:  [integer (:)]  fraction of veg not covered by snow (0/1 now) [-]
          dewmx          => canopystate_inst%dewmx_patch          , & ! Input:  [real(r8) (:) ]  Maximum allowed dew [mm]                
          elai           => canopystate_inst%elai_patch           , & ! Input:  [real(r8) (:) ]  one-sided leaf area index with burying by snow
          esai           => canopystate_inst%esai_patch           , & ! Input:  [real(r8) (:) ]  one-sided stem area index with burying by snow

          h2ocan         => waterstate_inst%h2ocan_patch          , & ! Input:  [real(r8) (:) ]  total canopy water (mm H2O)             
          snocan         => waterstate_inst%snocan_patch          , & ! Output: [real(r8) (:)   ]  canopy snow (mm H2O)             
          liqcan         => waterstate_inst%liqcan_patch          , & ! Output: [real(r8) (:)   ]  canopy liquid (mm H2O)

          fwet           => waterstate_inst%fwet_patch            , & ! Output: [real(r8) (:) ]  fraction of canopy that is wet (0 to 1) 
          fcansno        => waterstate_inst%fcansno_patch         , & ! Output: [real(r8) (:) ]  fraction of canopy that is snow covered (0 to 1) 
          fdry           => waterstate_inst%fdry_patch              & ! Output: [real(r8) (:) ]  fraction of foliage that is green and dry [-] (new)
          )

       ! Set status of snowveg_flag
       snowveg_onrad = IsSnowvegFlagOnRad()

       do fp = 1,numf
          p = filter(fp)
          if (frac_veg_nosno(p) == 1) then
             if (h2ocan(p) > 0._r8) then
                vegt    = frac_veg_nosno(p)*(elai(p) + esai(p))
                dewmxi  = 1.0_r8/dewmx(p)
                fwet(p) = ((dewmxi/vegt)*h2ocan(p))**0.666666666666_r8
                fwet(p) = min (fwet(p),maximum_leaf_wetted_fraction)   ! Check for maximum limit of fwet
                if (snowveg_onrad) then
                   if (snocan(p) > 0._r8) then
                      dewmxi  = 1.0_r8/dewmx(p)
                      fcansno(p) = ((dewmxi/(vegt*6.0_r8*10.0_r8))*snocan(p))**0.15_r8 ! must match snocanmx 
                      fcansno(p) = min (fcansno(p),1.0_r8)
                   else
                      fcansno(p) = 0._r8
                   end if
                else
                   fcansno(p) = 0._r8
                end if
             else
                fwet(p) = 0._r8
                fcansno(p) = 0._r8
             end if
             fdry(p) = (1._r8-fwet(p))*elai(p)/(elai(p)+esai(p))
          else
             fwet(p) = 0._r8
             fdry(p) = 0._r8
          end if
       end do

     end associate
   end subroutine FracWet

   !-----------------------------------------------------------------------
   subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
        waterstate_inst, no_update)
     !
     ! !DESCRIPTION:
     ! Determine fraction of land surfaces which are submerged  
     ! based on surface microtopography and surface water storage.
     !
     ! !USES:
     use shr_const_mod   , only : shr_const_pi
     use shr_spfn_mod    , only : erf => shr_spfn_erf
     use landunit_varcon , only : istsoil, istcrop
     !
     ! !ARGUMENTS:
     type(bounds_type)     , intent(in)           :: bounds           
     integer               , intent(in)           :: num_h2osfc       ! number of column points in column filter
     integer               , intent(in)           :: filter_h2osfc(:) ! column filter 
     type(waterstate_type) , intent(inout)        :: waterstate_inst
     integer               , intent(in), optional :: no_update        ! flag to make calculation w/o updating variables
     !
     ! !LOCAL VARIABLES:
     integer :: c,f,l          ! indices
     real(r8):: d,fd,dfdd      ! temporary variable for frac_h2o iteration
     real(r8):: sigma          ! microtopography pdf sigma in mm
     real(r8):: min_h2osfc
     !-----------------------------------------------------------------------

     associate(                                              & 
          micro_sigma  => col%micro_sigma                  , & ! Input:  [real(r8) (:)   ] microtopography pdf sigma (m)                     
          h2osno       => waterstate_inst%h2osno_col       , & ! Input:  [real(r8) (:)   ] snow water (mm H2O)                               
          h2osoi_liq   => waterstate_inst%h2osoi_liq_col   , & ! Output: [real(r8) (:,:) ] liquid water (col,lyr) [kg/m2]                  
          h2osfc       => waterstate_inst%h2osfc_col       , & ! Output: [real(r8) (:)   ] surface water (mm)                                
          frac_sno     => waterstate_inst%frac_sno_col     , & ! Output: [real(r8) (:)   ] fraction of ground covered by snow (0 to 1)       
          frac_sno_eff => waterstate_inst%frac_sno_eff_col , & ! Output: [real(r8) (:)   ] eff. fraction of ground covered by snow (0 to 1)  
          frac_h2osfc  => waterstate_inst%frac_h2osfc_col,   & ! Output: [real(r8) (:)   ] col fractional area with surface water greater than zero 
          frac_h2osfc_nosnow  => waterstate_inst%frac_h2osfc_nosnow_col    & ! Output: [real(r8) (:)   ] col fractional area with surface water greater than zero (if no snow present)
          )

       ! arbitrary lower limit on h2osfc for safer numerics...
       min_h2osfc=1.e-8_r8

       do f = 1, num_h2osfc
          c = filter_h2osfc(f)
          l = col%landunit(c)

          ! h2osfc only calculated for soil vegetated land units
          if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then

             !  Use newton-raphson method to iteratively determine frac_h2osfc
             !  based on amount of surface water storage (h2osfc) and 
             !  microtopography variability (micro_sigma)

             if (h2osfc(c) > min_h2osfc) then
                ! a cutoff is needed for numerical reasons...(nonconvergence after 5 iterations)
                d=0.0

                sigma=1.0e3 * micro_sigma(c) ! convert to mm
                do l=1,10
                   fd = 0.5*d*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) &
                        +sigma/sqrt(2.0*shr_const_pi)*exp(-d**2/(2.0*sigma**2)) &
                        -h2osfc(c)
                   dfdd = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))

                   d = d - fd/dfdd
                enddo
                !--  update the submerged areal fraction using the new d value
                frac_h2osfc(c) = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0))))

             else
                frac_h2osfc(c) = 0._r8
                h2osoi_liq(c,1) = h2osoi_liq(c,1) + h2osfc(c)
                h2osfc(c)=0._r8
             endif

             frac_h2osfc_nosnow(c) = frac_h2osfc(c)


             if (.not. present(no_update)) then

                ! adjust fh2o, fsno when sum is greater than zero
                if (frac_sno(c) > (1._r8 - frac_h2osfc(c)) .and. h2osno(c) > 0) then

                   if (frac_h2osfc(c) > 0.01_r8) then             
                      frac_h2osfc(c) = max(1.0_r8 - frac_sno(c),0.01_r8)
                      frac_sno(c) = 1.0_r8 - frac_h2osfc(c)
                   else
                      frac_sno(c) = 1.0_r8 - frac_h2osfc(c)
                   endif
                   frac_sno_eff(c)=frac_sno(c)

                endif

             endif ! end of no_update construct

          else !if landunit not istsoil/istcrop, set frac_h2osfc to zero

             frac_h2osfc(c) = 0._r8

          endif

       end do

     end associate

   end subroutine FracH2OSfc

   !-----------------------------------------------------------------------
   !BOP
   !
   ! !IROUTINE: IsSnowvegFlagOff
   !
   ! !INTERFACE:
   !
   logical function IsSnowvegFlagOff( )
     !
     ! !DESCRIPTION:
     !
     ! Return true if snowveg_flag is OFF
     !
     ! !USES:
     implicit none
     !EOP
     !-----------------------------------------------------------------------

     IsSnowvegFlagOff = (trim(snowveg_flag) == 'OFF')

   end function IsSnowvegFlagOff

   !-----------------------------------------------------------------------
   !BOP
   !
   ! !IROUTINE: IsSnowvegFlagOn
   !
   ! !INTERFACE:
   !
   logical function IsSnowvegFlagOn( )
     !
     ! !DESCRIPTION:
     !
     ! Return true if snowveg_flag is ON
     !
     ! !USES:
     implicit none
     !EOP
     !-----------------------------------------------------------------------

     IsSnowvegFlagOn = (trim(snowveg_flag) == 'ON')

   end function IsSnowvegFlagOn

   !-----------------------------------------------------------------------
   !BOP
   !
   ! !IROUTINE: IsSnowvegFlagOnRad
   !
   ! !INTERFACE:
   !
   logical function IsSnowvegFlagOnRad( )
     !
     ! !DESCRIPTION:
     !
     ! Return true if snowveg_flag is ON_RAD
     !
     ! !USES:
     implicit none
     !EOP
     !-----------------------------------------------------------------------

     IsSnowvegFlagOnRad = (trim(snowveg_flag) == 'ON_RAD')

   end function IsSnowvegFlagOnRad

end module CanopyHydrologyMod
