diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90
index 9d04feb559..d0f802cf9a 100644
--- a/src/biogeophys/BalanceCheckMod.F90
+++ b/src/biogeophys/BalanceCheckMod.F90
@@ -22,8 +22,6 @@ module BalanceCheckMod
   use WaterDiagnosticBulkType     , only : waterdiagnosticbulk_type
   use WaterBalanceType     , only : waterbalance_type
   use WaterfluxType      , only : waterflux_type
-  use IrrigationMod      , only : irrigation_type
-  use GlacierSurfaceMassBalanceMod, only : glacier_smb_type
   use TotalWaterAndHeatMod, only : ComputeWaterMassNonLake, ComputeWaterMassLake
   use GridcellType       , only : grc                
   use LandunitType       , only : lun                
@@ -101,7 +99,7 @@ end subroutine BeginWaterBalance
    !-----------------------------------------------------------------------
    subroutine BalanceCheck( bounds, &
         atm2lnd_inst, solarabs_inst, waterflux_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, &
-        irrigation_inst, glacier_smb_inst, energyflux_inst, canopystate_inst)
+        energyflux_inst, canopystate_inst)
      !
      ! !DESCRIPTION:
      ! This subroutine accumulates the numerical truncation errors of the water
@@ -133,8 +131,6 @@ subroutine BalanceCheck( bounds, &
      type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst
      type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst
      type(waterbalance_type) , intent(inout) :: waterbalance_inst
-     type(irrigation_type) , intent(in)    :: irrigation_inst
-     type(glacier_smb_type), intent(in)    :: glacier_smb_inst
      type(energyflux_type) , intent(inout) :: energyflux_inst
      type(canopystate_type), intent(inout) :: canopystate_inst
      !
@@ -200,9 +196,9 @@ subroutine BalanceCheck( bounds, &
           snow_sources            =>    waterflux_inst%snow_sources_col         , & ! Output: [real(r8) (:)   ]  snow sources (mm H2O /s)  
           snow_sinks              =>    waterflux_inst%snow_sinks_col           , & ! Output: [real(r8) (:)   ]  snow sinks (mm H2O /s)    
 
-          qflx_irrig              =>    irrigation_inst%qflx_irrig_col          , & ! Input:  [real(r8) (:)   ]  irrigation flux (mm H2O /s)             
+          qflx_irrig              =>    waterflux_inst%qflx_irrig_col          , & ! Input:  [real(r8) (:)   ]  irrigation flux (mm H2O /s)             
 
-          qflx_glcice_dyn_water_flux => glacier_smb_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)]  water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
+          qflx_glcice_dyn_water_flux => waterflux_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)]  water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
 
           eflx_lwrad_out          =>    energyflux_inst%eflx_lwrad_out_patch    , & ! Input:  [real(r8) (:)   ]  emitted infrared (longwave) radiation (W/m**2)
           eflx_lwrad_net          =>    energyflux_inst%eflx_lwrad_net_patch    , & ! Input:  [real(r8) (:)   ]  net infrared (longwave) rad (W/m**2) [+ = to atm]
diff --git a/src/biogeophys/CanopyHydrologyMod.F90 b/src/biogeophys/CanopyHydrologyMod.F90
index a542cf0707..da37e72bc7 100644
--- a/src/biogeophys/CanopyHydrologyMod.F90
+++ b/src/biogeophys/CanopyHydrologyMod.F90
@@ -24,7 +24,6 @@ module CanopyHydrologyMod
   use WaterfluxType   , only : waterflux_type
   use WaterStateBulkType  , only : waterstatebulk_type
   use WaterDiagnosticBulkType  , only : waterdiagnosticbulk_type
-  use IrrigationMod   , only : irrigation_type
   use ColumnType      , only : col                
   use PatchType       , only : patch                
   !
@@ -143,7 +142,7 @@ end subroutine CanopyHydrology_readnl
    subroutine CanopyHydrology(bounds, &
         num_nolakec, filter_nolakec, num_nolakep, filter_nolakep, &
         atm2lnd_inst, canopystate_inst, temperature_inst, &
-        aerosol_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterflux_inst, irrigation_inst)
+        aerosol_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterflux_inst)
      !
      ! !DESCRIPTION:
      ! Calculation of
@@ -178,7 +177,6 @@ subroutine CanopyHydrology(bounds, &
      type(waterstatebulk_type)  , intent(inout) :: waterstatebulk_inst
      type(waterdiagnosticbulk_type)  , intent(inout) :: waterdiagnosticbulk_inst
      type(waterflux_type)   , intent(inout) :: waterflux_inst
-     type(irrigation_type)  , intent(in)    :: irrigation_inst
      !
      ! !LOCAL VARIABLES:
      integer  :: f                                            ! filter index
@@ -266,8 +264,7 @@ subroutine CanopyHydrology(bounds, &
           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_irrig           => waterflux_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]    
           )
diff --git a/src/biogeophys/GlacierSurfaceMassBalanceMod.F90 b/src/biogeophys/GlacierSurfaceMassBalanceMod.F90
index e4fa704775..a16877b3d1 100644
--- a/src/biogeophys/GlacierSurfaceMassBalanceMod.F90
+++ b/src/biogeophys/GlacierSurfaceMassBalanceMod.F90
@@ -8,9 +8,8 @@ module GlacierSurfaceMassBalanceMod
 #include "shr_assert.h"
   use shr_kind_mod   , only : r8 => shr_kind_r8
   use shr_log_mod    , only : errMsg => shr_log_errMsg
-  use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
   use decompMod      , only : bounds_type
-  use clm_varcon     , only : spval, secspday
+  use clm_varcon     , only : secspday
   use clm_varpar     , only : nlevgrnd
   use clm_varctl     , only : glc_snow_persistence_max_days
   use clm_time_manager, only : get_step_size
@@ -29,20 +28,6 @@ module GlacierSurfaceMassBalanceMod
   type, public :: glacier_smb_type
      private
 
-     ! ------------------------------------------------------------------------
-     ! Public data
-     ! ------------------------------------------------------------------------
-
-     real(r8), pointer, public :: qflx_glcice_col(:) ! col net flux of new glacial ice (growth - melt) (mm H2O/s), passed to GLC; only valid inside the do_smb_c filter
-     real(r8), pointer, public :: qflx_glcice_dyn_water_flux_col(:) ! col water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system); valid for all columns
-
-     ! ------------------------------------------------------------------------
-     ! Private data
-     ! ------------------------------------------------------------------------
-
-     real(r8), pointer :: qflx_glcice_frz_col (:)   ! col ice growth (positive definite) (mm H2O/s); only valid inside the do_smb_c filter
-     real(r8), pointer :: qflx_glcice_melt_col(:)   ! col ice melt (positive definite) (mm H2O/s); only valid inside the do_smb_c filter
-
    contains
 
      ! ------------------------------------------------------------------------
@@ -57,14 +42,6 @@ module GlacierSurfaceMassBalanceMod
      procedure, public :: ComputeSurfaceMassBalance ! compute fluxes other than ice melt
      procedure, public :: AdjustRunoffTerms    ! adjust liquid and ice runoff fluxes due to glacier fluxes
 
-     ! ------------------------------------------------------------------------
-     ! Private routines
-     ! ------------------------------------------------------------------------
-
-     procedure, private :: InitAllocate
-     procedure, private :: InitHistory
-     procedure, private :: InitCold
-
   end type glacier_smb_type
 
   character(len=*), parameter, private :: sourcefile = &
@@ -82,86 +59,16 @@ subroutine Init(this, bounds)
     type(bounds_type), intent(in) :: bounds
     !-----------------------------------------------------------------------
 
-    call this%InitAllocate(bounds)
-    call this%InitHistory(bounds)
-    call this%InitCold(bounds)
+    ! Nothing to do for now
   end subroutine Init
 
-  !-----------------------------------------------------------------------
-  subroutine InitAllocate(this, bounds)
-    class(glacier_smb_type), intent(inout) :: this
-    type(bounds_type), intent(in) :: bounds
-
-    integer :: begc, endc
-    !-----------------------------------------------------------------------
-
-    begc = bounds%begc; endc = bounds%endc
-
-    allocate(this%qflx_glcice_col               (begc:endc)) ; this%qflx_glcice_col                (:) = nan
-    allocate(this%qflx_glcice_dyn_water_flux_col(begc:endc)) ; this%qflx_glcice_dyn_water_flux_col (:) = nan
-    allocate(this%qflx_glcice_frz_col           (begc:endc)) ; this%qflx_glcice_frz_col            (:) = nan
-    allocate(this%qflx_glcice_melt_col          (begc:endc)) ; this%qflx_glcice_melt_col           (:) = nan
-  end subroutine InitAllocate
-
-  !-----------------------------------------------------------------------
-  subroutine InitHistory(this, bounds)
-    !
-    ! !USES:
-    use histFileMod , only : hist_addfld1d
-    !
-    ! !ARGUMENTS:
-    class(glacier_smb_type), intent(inout) :: this
-    type(bounds_type), intent(in) :: bounds
-    !
-    ! !LOCAL VARIABLES:
-    integer :: begc, endc
-    !-----------------------------------------------------------------------
-
-    begc = bounds%begc; endc = bounds%endc
-
-    this%qflx_glcice_col(begc:endc) = spval
-    call hist_addfld1d (fname='QICE',  units='mm/s',  &
-         avgflag='A', long_name='ice growth/melt', &
-         ptr_col=this%qflx_glcice_col, l2g_scale_type='ice')
-
-    this%qflx_glcice_frz_col(begc:endc) = spval
-    call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
-         avgflag='A', long_name='ice growth', &
-         ptr_col=this%qflx_glcice_frz_col, l2g_scale_type='ice')
-
-    this%qflx_glcice_melt_col(begc:endc) = spval
-    call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
-         avgflag='A', long_name='ice melt', &
-         ptr_col=this%qflx_glcice_melt_col, l2g_scale_type='ice')
-
-  end subroutine InitHistory
-
-  !-----------------------------------------------------------------------
-  subroutine InitCold(this, bounds)
-    class(glacier_smb_type), intent(inout) :: this
-    type(bounds_type), intent(in) :: bounds
-
-    integer :: c
-    !-----------------------------------------------------------------------
-
-    ! Initialize qflx_glcice_dyn_water_flux_col to 0 for all columns because we want this
-    ! flux to remain 0 for columns where is is never set, including non-glacier columns.
-    !
-    ! Other fluxes intentionally remain unset (spval) outside the do_smb filter, so that
-    ! they are flagged as missing value outside that filter.
-    do c = bounds%begc, bounds%endc
-       this%qflx_glcice_dyn_water_flux_col(c) = 0._r8
-    end do
-
-  end subroutine InitCold
-
   ! ========================================================================
   ! Science routines
   ! ========================================================================
 
   !-----------------------------------------------------------------------
   subroutine HandleIceMelt(this, bounds, num_do_smb_c, filter_do_smb_c, &
-       waterstatebulk_inst)
+       waterstatebulk_inst, waterflux_inst)
     !
     ! !DESCRIPTION:
     ! Compute ice melt in glacier columns, and convert liquid back to ice
@@ -184,11 +91,12 @@ subroutine HandleIceMelt(this, bounds, num_do_smb_c, filter_do_smb_c, &
     ! may be best to keep this separate.
     !
     ! !ARGUMENTS:
-    class(glacier_smb_type), intent(inout) :: this
+    class(glacier_smb_type), intent(in) :: this
     type(bounds_type), intent(in) :: bounds
     integer, intent(in) :: num_do_smb_c        ! number of column points in filter_do_smb_c
     integer, intent(in) :: filter_do_smb_c(:)  ! column filter for points where SMB is calculated
     type(waterstatebulk_type), intent(inout) :: waterstatebulk_inst
+    type(waterflux_type), intent(inout) :: waterflux_inst
     !
     ! !LOCAL VARIABLES:
     integer :: j
@@ -199,9 +107,9 @@ subroutine HandleIceMelt(this, bounds, num_do_smb_c, filter_do_smb_c, &
     !-----------------------------------------------------------------------
 
     associate( &
-         qflx_glcice_melt => this%qflx_glcice_melt_col , & ! Output: [real(r8) (:)   ] ice melt (positive definite) (mm H2O/s)
          h2osoi_liq       => waterstatebulk_inst%h2osoi_liq_col      , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2)
-         h2osoi_ice       => waterstatebulk_inst%h2osoi_ice_col        & ! Output: [real(r8) (:,:) ] ice lens (kg/m2)
+         h2osoi_ice       => waterstatebulk_inst%h2osoi_ice_col      , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2)
+         qflx_glcice_melt => waterflux_inst%qflx_glcice_melt_col       & ! Output: [real(r8) (:)   ] ice melt (positive definite) (mm H2O/s)
          )
 
     dtime = get_step_size()
@@ -243,14 +151,13 @@ subroutine ComputeSurfaceMassBalance(this, bounds, num_allc, filter_allc, &
     ! !DESCRIPTION:
     ! Compute glacier fluxes other than ice melt.
     !
-    ! This sets the public fields qflx_glcice_col and qflx_glcice_dyn_water_flux_col to
-    ! their final values.
+    ! This sets qflx_glcice_col and qflx_glcice_dyn_water_flux_col.
     !
     ! Should be called after HandleIceMelt, and after waterflux_inst%qflx_snwcp_ice_col is
     ! computed
     !
     ! !ARGUMENTS:
-    class(glacier_smb_type), intent(inout) :: this
+    class(glacier_smb_type), intent(in) :: this
     type(bounds_type), intent(in) :: bounds
     integer, intent(in) :: num_allc        ! number of column points in filter_allc
     integer, intent(in) :: filter_allc(:)  ! column filter for all points
@@ -258,7 +165,7 @@ subroutine ComputeSurfaceMassBalance(this, bounds, num_allc, filter_allc, &
     integer, intent(in) :: filter_do_smb_c(:)  ! column filter for points where SMB is calculated
     type(glc2lnd_type), intent(in) :: glc2lnd_inst
     type(waterstatebulk_type), intent(in) :: waterstatebulk_inst
-    type(waterflux_type), intent(in) :: waterflux_inst
+    type(waterflux_type), intent(inout) :: waterflux_inst
     !
     ! !LOCAL VARIABLES:
     integer :: fc, c, l, g
@@ -267,13 +174,13 @@ subroutine ComputeSurfaceMassBalance(this, bounds, num_allc, filter_allc, &
     !-----------------------------------------------------------------------
 
     associate( &
-         qflx_glcice                => this%qflx_glcice_col                    , & ! Output: [real(r8) (:)] net flux of new glacial ice (growth - melt) (mm H2O/s)
-         qflx_glcice_frz            => this%qflx_glcice_frz_col                , & ! Output: [real(r8) (:)] ice growth (positive definite) (mm H2O/s)
-         qflx_glcice_dyn_water_flux => this%qflx_glcice_dyn_water_flux_col     , & ! Output: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
-         qflx_glcice_melt           => this%qflx_glcice_melt_col               , & ! Input: [real(r8) (:)] ice melt (positive definite) (mm H2O/s)
          glc_dyn_runoff_routing     => glc2lnd_inst%glc_dyn_runoff_routing_grc , & ! Input:  [real(r8) (:)]  whether we're doing runoff routing appropriate for having a dynamic icesheet
-         snow_persistence           => waterstatebulk_inst%snow_persistence_col    , & ! Input: [real(r8) (:)]  counter for length of time snow-covered
-         qflx_snwcp_ice             => waterflux_inst%qflx_snwcp_ice_col         & ! Input: [real(r8) (:)]  excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]
+         snow_persistence           => waterstatebulk_inst%snow_persistence_col, & ! Input: [real(r8) (:)]  counter for length of time snow-covered
+         qflx_snwcp_ice             => waterflux_inst%qflx_snwcp_ice_col       , & ! Input: [real(r8) (:)]  excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]
+         qflx_glcice_melt           => waterflux_inst%qflx_glcice_melt_col     , & ! Input: [real(r8) (:)] ice melt (positive definite) (mm H2O/s)
+         qflx_glcice                => waterflux_inst%qflx_glcice_col          , & ! Output: [real(r8) (:)] net flux of new glacial ice (growth - melt) (mm H2O/s)
+         qflx_glcice_frz            => waterflux_inst%qflx_glcice_frz_col      , & ! Output: [real(r8) (:)] ice growth (positive definite) (mm H2O/s)
+         qflx_glcice_dyn_water_flux => waterflux_inst%qflx_glcice_dyn_water_flux_col & ! Output: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
          )
 
     ! NOTE(wjs, 2016-06-29) The following initialization is done in case the columns
@@ -340,7 +247,7 @@ end subroutine ComputeSurfaceMassBalance
 
   !-----------------------------------------------------------------------
   subroutine AdjustRunoffTerms(this, bounds, num_do_smb_c, filter_do_smb_c, &
-       glc2lnd_inst, qflx_qrgwl, qflx_ice_runoff_snwcp)
+       waterflux_inst, glc2lnd_inst, qflx_qrgwl, qflx_ice_runoff_snwcp)
     !
     ! !DESCRIPTION:
     ! Adjust liquid and ice runoff fluxes due to glacier fluxes
@@ -355,6 +262,7 @@ subroutine AdjustRunoffTerms(this, bounds, num_do_smb_c, filter_do_smb_c, &
     type(bounds_type), intent(in) :: bounds
     integer, intent(in) :: num_do_smb_c        ! number of column points in filter_do_smb_c
     integer, intent(in) :: filter_do_smb_c(:)  ! column filter for points where SMB is calculated
+    type(waterflux_type), intent(in) :: waterflux_inst
     type(glc2lnd_type), intent(in) :: glc2lnd_inst
     real(r8), intent(inout) :: qflx_qrgwl( bounds%begc: ) ! col qflx_surf at glaciers, wetlands, lakes
     real(r8), intent(inout) :: qflx_ice_runoff_snwcp( bounds%begc: ) ! col solid runoff from snow capping (mm H2O /s)
@@ -369,8 +277,8 @@ subroutine AdjustRunoffTerms(this, bounds, num_do_smb_c, filter_do_smb_c, &
     SHR_ASSERT_ALL((ubound(qflx_ice_runoff_snwcp) == (/bounds%endc/)), errMsg(sourcefile, __LINE__))
 
     associate( &
-         qflx_glcice_frz        => this%qflx_glcice_frz_col                   , & ! Input: [real(r8) (:)] ice growth (positive definite) (mm H2O/s)
-         qflx_glcice_melt       => this%qflx_glcice_melt_col                  , & ! Input: [real(r8) (:)] ice melt (positive definite) (mm H2O/s)
+         qflx_glcice_frz        => waterflux_inst%qflx_glcice_frz_col         , & ! Input: [real(r8) (:)] ice growth (positive definite) (mm H2O/s)
+         qflx_glcice_melt       => waterflux_inst%qflx_glcice_melt_col        , & ! Input: [real(r8) (:)] ice melt (positive definite) (mm H2O/s)
          glc_dyn_runoff_routing =>    glc2lnd_inst%glc_dyn_runoff_routing_grc   & ! Input: [real(r8) (:)]  gridcell fraction coupled to dynamic ice sheet
          )
 
diff --git a/src/biogeophys/HydrologyDrainageMod.F90 b/src/biogeophys/HydrologyDrainageMod.F90
index 04c47f35b6..ffb0b8da61 100644
--- a/src/biogeophys/HydrologyDrainageMod.F90
+++ b/src/biogeophys/HydrologyDrainageMod.F90
@@ -18,7 +18,6 @@ module HydrologyDrainageMod
   use WaterStateBulkType    , only : waterstatebulk_type
   use WaterDiagnosticBulkType    , only : waterdiagnosticbulk_type
   use WaterBalanceType    , only : waterbalance_type
-  use IrrigationMod     , only : irrigation_type
   use GlacierSurfaceMassBalanceMod, only : glacier_smb_type
   use TotalWaterAndHeatMod, only : ComputeWaterMassNonLake
   use LandunitType      , only : lun                
@@ -42,7 +41,7 @@ subroutine HydrologyDrainage(bounds,               &
        num_do_smb_c, filter_do_smb_c,                &
        atm2lnd_inst, glc2lnd_inst, temperature_inst, &
        soilhydrology_inst, soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, &
-       irrigation_inst, glacier_smb_inst)
+       glacier_smb_inst)
     !
     ! !DESCRIPTION:
     ! Calculates soil/snow hydrology with drainage (subsurface runoff)
@@ -76,7 +75,6 @@ subroutine HydrologyDrainage(bounds,               &
     type(waterdiagnosticbulk_type)    , intent(inout) :: waterdiagnosticbulk_inst
     type(waterbalance_type)    , intent(inout) :: waterbalance_inst
     type(waterflux_type)     , intent(inout) :: waterflux_inst
-    type(irrigation_type)    , intent(in)    :: irrigation_inst
     type(glacier_smb_type)   , intent(in)    :: glacier_smb_inst
     !
     ! !LOCAL VARIABLES:
@@ -117,7 +115,7 @@ subroutine HydrologyDrainage(bounds,               &
          qflx_runoff_r      => waterflux_inst%qflx_runoff_r_col      , & ! Rural total runoff 
                                                                          ! (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s)
          qflx_ice_runoff_snwcp => waterflux_inst%qflx_ice_runoff_snwcp_col, &  ! solid runoff from snow capping (mm H2O /s)
-         qflx_irrig         => irrigation_inst%qflx_irrig_col          & ! irrigation flux (mm H2O /s)   
+         qflx_irrig         => waterflux_inst%qflx_irrig_col          & ! irrigation flux (mm H2O /s)   
          )
 
       ! Determine time step and step size
@@ -199,6 +197,7 @@ subroutine HydrologyDrainage(bounds,               &
       ! qflx_qrgwl and qflx_ice_runoff_snwcp, but before the ues of qflx_qrgwl in
       ! qflx_runoff.
       call glacier_smb_inst%AdjustRunoffTerms(bounds, num_do_smb_c, filter_do_smb_c, &
+           waterflux_inst, &
            glc2lnd_inst, &
            qflx_qrgwl = qflx_qrgwl(bounds%begc:bounds%endc), &
            qflx_ice_runoff_snwcp = qflx_ice_runoff_snwcp(bounds%begc:bounds%endc))
diff --git a/src/biogeophys/HydrologyNoDrainageMod.F90 b/src/biogeophys/HydrologyNoDrainageMod.F90
index 1658ea1935..64f4a75933 100644
--- a/src/biogeophys/HydrologyNoDrainageMod.F90
+++ b/src/biogeophys/HydrologyNoDrainageMod.F90
@@ -193,7 +193,7 @@ subroutine HydrologyNoDrainage(bounds, &
            soilhydrology_inst, soilstate_inst, waterflux_inst)
 
       call SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, &
-           waterflux_inst, saturated_excess_runoff_inst, waterdiagnosticbulk_inst)
+           waterflux_inst, waterdiagnosticbulk_inst)
 
       call infiltration_excess_runoff_inst%InfiltrationExcessRunoff( &
            bounds, num_hydrologyc, filter_hydrologyc, &
@@ -201,7 +201,7 @@ subroutine HydrologyNoDrainage(bounds, &
            waterdiagnosticbulk_inst)
 
       call RouteInfiltrationExcess(bounds, num_hydrologyc, filter_hydrologyc, &
-           waterflux_inst, infiltration_excess_runoff_inst, soilhydrology_inst)
+           waterflux_inst, soilhydrology_inst)
 
       call UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, &
            infiltration_excess_runoff_inst, &
@@ -213,7 +213,7 @@ subroutine HydrologyNoDrainage(bounds, &
 
       call TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
            num_urbanc, filter_urbanc, &
-           waterflux_inst, soilhydrology_inst, saturated_excess_runoff_inst, &
+           waterflux_inst, soilhydrology_inst, &
            waterstatebulk_inst)
 
       call UpdateUrbanPonding(bounds, num_urbanc, filter_urbanc, &
diff --git a/src/biogeophys/InfiltrationExcessRunoffMod.F90 b/src/biogeophys/InfiltrationExcessRunoffMod.F90
index 79feb8f477..a73d0d366c 100644
--- a/src/biogeophys/InfiltrationExcessRunoffMod.F90
+++ b/src/biogeophys/InfiltrationExcessRunoffMod.F90
@@ -30,14 +30,7 @@ module InfiltrationExcessRunoffMod
      private
      ! Public data members
      ! Note: these should be treated as read-only by other modules
-
-     ! These are valid within the hydrology filter.
-     !
-     ! Both of these give averages over the entire column. However, qinmax is implicitly
-     ! 0 over the fraction of the column given by fsat, and qflx_infl_excess is
-     ! implicitly 0 over both fsat and frac_h2osfc.
      real(r8), pointer, public :: qinmax_col(:)  ! maximum infiltration rate (mm H2O /s)
-     real(r8), pointer, public :: qflx_infl_excess_col(:)  ! infiltration excess runoff (mm H2O /s)
 
      ! Private data members
      integer :: qinmax_method
@@ -121,7 +114,6 @@ subroutine InitAllocate(this, bounds)
 
     begc = bounds%begc; endc= bounds%endc
     allocate(this%qinmax_col          (begc:endc)); this%qinmax_col          (:) = nan
-    allocate(this%qflx_infl_excess_col(begc:endc)); this%qflx_infl_excess_col(:) = nan
 
   end subroutine InitAllocate
 
@@ -183,8 +175,13 @@ subroutine InfiltrationExcessRunoff(this, bounds, num_hydrologyc, filter_hydrolo
     ! !DESCRIPTION:
     ! Calculate surface runoff due to infiltration excess
     !
+    ! Sets waterflux_inst%qflx_infl_excess_col and this%qinmax_col. These are valid within
+    ! the hydrology filter. Both of these give averages over the entire column. However,
+    ! qinmax is implicitly 0 over the fraction of the column given by fsat, and
+    ! qflx_infl_excess is implicitly 0 over both fsat and frac_h2osfc.
+    !
     ! !ARGUMENTS:
-    class(infiltration_excess_runoff_type) , intent(inout) :: this
+    class(infiltration_excess_runoff_type) , intent(in)    :: this
     type(bounds_type)                      , intent(in)    :: bounds               
     integer                                , intent(in)    :: num_hydrologyc       ! number of column soil points in column filter
     integer                                , intent(in)    :: filter_hydrologyc(:) ! column filter for soil points
@@ -203,10 +200,10 @@ subroutine InfiltrationExcessRunoff(this, bounds, num_hydrologyc, filter_hydrolo
 
     associate( &
          qinmax           =>    this%qinmax_col                     , & ! Output: [real(r8) (:)   ]  maximum infiltration rate (mm H2O /s)
-         qflx_infl_excess =>    this%qflx_infl_excess_col           , & ! Output: [real(r8) (:)   ]  infiltration excess runoff (mm H2O /s)
 
          fsat             =>   saturated_excess_runoff_inst%fsat_col, & ! Input:  [real(r8) (:)   ]  fractional area with water table at surface       
 
+         qflx_infl_excess =>    waterflux_inst%qflx_infl_excess_col , & ! Output: [real(r8) (:)   ]  infiltration excess runoff (mm H2O /s)
          qflx_in_soil     =>    waterflux_inst%qflx_in_soil_col     , & ! Input:  [real(r8) (:)   ]  surface input to soil (mm/s)
 
          frac_h2osfc      =>    waterdiagnosticbulk_inst%frac_h2osfc_col & ! Input:  [real(r8) (:)   ]  fraction of ground covered by surface water (0 to 1)
diff --git a/src/biogeophys/IrrigationMod.F90 b/src/biogeophys/IrrigationMod.F90
index cccb8dbc6b..056fe7ddbd 100644
--- a/src/biogeophys/IrrigationMod.F90
+++ b/src/biogeophys/IrrigationMod.F90
@@ -10,13 +10,9 @@ module IrrigationMod
   !     needed for the next call to ApplyIrrigation. This should be called once per
   !     timestep.
   ! 
-  !   - Call ApplyIrrigation in order to calculate qflx_irrig. This should be called
-  !     exactly once per time step, before the first time qflx_irrig is needed by other
-  !     parts of the code. It is acceptable for this to be called earlier in the timestep
-  !     than CalcIrrigationNeeded.
-  !
-  !   - Access the timestep's irrigation flux via qflx_irrig_patch or
-  !     qflx_irrig_col. These should be treated as read-only.
+  !   - Call ApplyIrrigation in order to calculate qflx_irrig_patch and qflx_irrig_col. It
+  !     is acceptable for this to be called earlier in the timestep than
+  !     CalcIrrigationNeeded.
   !
   ! Design notes:
   !
@@ -55,6 +51,7 @@ module IrrigationMod
   use clm_varpar       , only : nlevsoi, nlevgrnd
   use clm_time_manager , only : get_step_size
   use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
+  use WaterfluxType    , only : waterflux_type
   use GridcellType     , only : grc                
   use ColumnType       , only : col                
   use PatchType        , only : patch                
@@ -120,10 +117,6 @@ module IrrigationMod
 
   type, public :: irrigation_type
      private
-     ! Public data members
-     ! Note: these should be treated as read-only by other modules
-     real(r8), pointer, public :: qflx_irrig_patch(:) ! patch irrigation flux (mm H2O/s)
-     real(r8), pointer, public :: qflx_irrig_col  (:) ! col irrigation flux (mm H2O/s)
 
      ! Private data members; set in initialization:
      type(irrigation_params_type) :: params
@@ -506,9 +499,7 @@ subroutine IrrigationInitAllocate(this, bounds)
     begp = bounds%begp; endp= bounds%endp
     begc = bounds%begc; endc= bounds%endc
 
-    allocate(this%qflx_irrig_patch         (begp:endp))          ; this%qflx_irrig_patch         (:)   = nan
     allocate(this%qflx_irrig_demand_patch  (begp:endp))          ; this%qflx_irrig_demand_patch  (:)   = nan
-    allocate(this%qflx_irrig_col           (begc:endc))          ; this%qflx_irrig_col           (:)   = nan
     allocate(this%relsat_wilting_point_col (begc:endc,nlevsoi)) ; this%relsat_wilting_point_col (:,:) = nan
     allocate(this%relsat_target_col        (begc:endc,nlevsoi)) ; this%relsat_target_col        (:,:) = nan
     allocate(this%irrig_rate_patch         (begp:endp))          ; this%irrig_rate_patch         (:)   = nan
@@ -538,11 +529,6 @@ subroutine IrrigationInitHistory(this, bounds)
 
     begp = bounds%begp; endp= bounds%endp
 
-    this%qflx_irrig_patch(begp:endp) = spval
-    call hist_addfld1d (fname='QIRRIG', units='mm/s', &
-         avgflag='A', long_name='water added through irrigation', &
-         ptr_patch=this%qflx_irrig_patch)
-
     this%qflx_irrig_demand_patch(begp:endp) = spval
     call hist_addfld1d (fname='QIRRIG_DEMAND', units='mm/s', &
          avgflag='A', long_name='irrigation demand', &
@@ -733,9 +719,7 @@ subroutine IrrigationClean(this)
     character(len=*), parameter :: subname = 'Clean'
     !-----------------------------------------------------------------------
     
-    deallocate(this%qflx_irrig_patch)
     deallocate(this%qflx_irrig_demand_patch)
-    deallocate(this%qflx_irrig_col)
     deallocate(this%relsat_wilting_point_col)
     deallocate(this%relsat_target_col)
     deallocate(this%irrig_rate_patch)
@@ -750,19 +734,21 @@ end subroutine IrrigationClean
   ! ========================================================================
   
   !-----------------------------------------------------------------------
-  subroutine ApplyIrrigation(this, bounds)
+  subroutine ApplyIrrigation(this, bounds, waterflux_inst)
     !
     ! !DESCRIPTION:
     ! Apply the irrigation computed by CalcIrrigationNeeded to qflx_irrig.
     !
-    ! Should be called once, AND ONLY ONCE, per time step. After this is called, you may
-    ! access qflx_irrig_patch or qflx_irrig_col.
+    ! Sets waterflux_inst%qflx_irrig_patch and waterflux_inst%qflx_irrig_col.
+    !
+    ! Should be called once, AND ONLY ONCE, per time step.
     !
     ! !USES:
     !
     ! !ARGUMENTS:
     class(irrigation_type) , intent(inout) :: this
     type(bounds_type)      , intent(in)    :: bounds
+    type(waterflux_type)   , intent(inout) :: waterflux_inst
     !
     ! !LOCAL VARIABLES:
     integer :: p  ! patch index
@@ -772,28 +758,35 @@ subroutine ApplyIrrigation(this, bounds)
 
     !-----------------------------------------------------------------------
 
-    ! This should be called exactly once per time step, so that this counter decrease
+    ! This should be called exactly once per time step, so that the counter decrease
     ! works correctly.
 
+    associate( &
+         qflx_irrig_patch => waterflux_inst%qflx_irrig_patch , & ! Output: [real(r8) (:)] patch irrigation flux (mm H2O/s)
+         qflx_irrig_col   => waterflux_inst%qflx_irrig_col     & ! Output: [real(r8) (:)] col irrigation flux (mm H2O/s)
+         )
+
     do p = bounds%begp, bounds%endp
        g = patch%gridcell(p)
 
        if (this%n_irrig_steps_left_patch(p) > 0) then
-          this%qflx_irrig_patch(p)         = this%irrig_rate_patch(p)
+          qflx_irrig_patch(p)              = this%irrig_rate_patch(p)
           this%qflx_irrig_demand_patch(p)  = this%irrig_rate_demand_patch(p)
           this%n_irrig_steps_left_patch(p) = this%n_irrig_steps_left_patch(p) - 1
        else
-          this%qflx_irrig_patch(p)        = 0._r8
+          qflx_irrig_patch(p)             = 0._r8
           this%qflx_irrig_demand_patch(p) = 0._r8
        end if
 
     end do
 
     call p2c (bounds = bounds, &
-         parr = this%qflx_irrig_patch(bounds%begp:bounds%endp), &
-         carr = this%qflx_irrig_col(bounds%begc:bounds%endc), &
+         parr = qflx_irrig_patch(bounds%begp:bounds%endp), &
+         carr = qflx_irrig_col(bounds%begc:bounds%endc), &
          p2c_scale_type = 'unity')
 
+    end associate
+
   end subroutine ApplyIrrigation
 
 
diff --git a/src/biogeophys/SaturatedExcessRunoffMod.F90 b/src/biogeophys/SaturatedExcessRunoffMod.F90
index 51f1c8369d..38bfaa8805 100644
--- a/src/biogeophys/SaturatedExcessRunoffMod.F90
+++ b/src/biogeophys/SaturatedExcessRunoffMod.F90
@@ -29,7 +29,6 @@ module SaturatedExcessRunoffMod
      private
      ! Public data members
      ! Note: these should be treated as read-only by other modules
-     real(r8), pointer, public :: qflx_sat_excess_surf_col(:) ! surface runoff due to saturated surface (mm H2O /s)
      real(r8), pointer, public :: fsat_col(:) ! fractional area with water table at surface
 
      ! Private data members
@@ -106,7 +105,6 @@ subroutine InitAllocate(this, bounds)
 
     begc = bounds%begc; endc= bounds%endc
 
-    allocate(this%qflx_sat_excess_surf_col(begc:endc)) ; this%qflx_sat_excess_surf_col(:) = nan
     allocate(this%fsat_col(begc:endc))                 ; this%fsat_col(:)                 = nan
     allocate(this%fcov_col(begc:endc))                 ; this%fcov_col(:)                 = nan   
 
@@ -180,6 +178,8 @@ subroutine SaturatedExcessRunoff (this, bounds, num_hydrologyc, filter_hydrology
     ! !DESCRIPTION:
     ! Calculate surface runoff due to saturated surface
     !
+    ! Sets this%fsat_col and waterflux_inst%qflx_sat_excess_surf_col
+    !
     ! !ARGUMENTS:
     class(saturated_excess_runoff_type), intent(inout) :: this
     type(bounds_type)        , intent(in)    :: bounds               
@@ -199,10 +199,10 @@ subroutine SaturatedExcessRunoff (this, bounds, num_hydrologyc, filter_hydrology
     associate(                                                        & 
          fcov                   =>    this%fcov_col                          , & ! Output: [real(r8) (:)   ]  fractional impermeable area
          fsat                   =>    this%fsat_col                          , & ! Output: [real(r8) (:)   ]  fractional area with water table at surface
-         qflx_sat_excess_surf   =>    this%qflx_sat_excess_surf_col          , & ! Output: [real(r8) (:)   ]  surface runoff due to saturated surface (mm H2O /s)
 
          snl                    =>    col%snl                                , & ! Input:  [integer  (:)   ]  minus number of snow layers
 
+         qflx_sat_excess_surf   =>    waterflux_inst%qflx_sat_excess_surf_col, & ! Output: [real(r8) (:)   ]  surface runoff due to saturated surface (mm H2O /s)
          qflx_floodc            =>    waterflux_inst%qflx_floodc_col         , & ! Input:  [real(r8) (:)   ]  column flux of flood water from RTM
          qflx_rain_plus_snomelt => waterflux_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:)   ] rain plus snow melt falling on the soil (mm/s)
 
diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90
index ad1fb34356..76ce4b078e 100644
--- a/src/biogeophys/SoilHydrologyMod.F90
+++ b/src/biogeophys/SoilHydrologyMod.F90
@@ -22,7 +22,6 @@ module SoilHydrologyMod
   use InfiltrationExcessRunoffMod, only : infiltration_excess_runoff_type
   use SoilHydrologyType , only : soilhydrology_type  
   use SoilStateType     , only : soilstate_type
-  use SaturatedExcessRunoffMod, only : saturated_excess_runoff_type
   use WaterfluxType     , only : waterflux_type
   use WaterStateBulkType    , only : waterstatebulk_type
   use WaterDiagnosticBulkType    , only : waterdiagnosticbulk_type
@@ -189,7 +188,7 @@ end subroutine SetSoilWaterFractions
 
    !-----------------------------------------------------------------------
    subroutine SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, &
-        waterflux_inst, saturated_excess_runoff_inst, waterdiagnosticbulk_inst)
+        waterflux_inst, waterdiagnosticbulk_inst)
      !
      ! !DESCRIPTION:
      ! Set various input fluxes of water
@@ -199,7 +198,6 @@ subroutine SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, &
      integer                    , intent(in)    :: num_hydrologyc       ! number of column soil points in column filter
      integer                    , intent(in)    :: filter_hydrologyc(:) ! column filter for soil points
      type(waterflux_type)       , intent(inout) :: waterflux_inst
-     type(saturated_excess_runoff_type) , intent(in) :: saturated_excess_runoff_inst
      type(waterdiagnosticbulk_type)     , intent(in) :: waterdiagnosticbulk_inst
      !
      ! !LOCAL VARIABLES:
@@ -222,8 +220,7 @@ subroutine SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, &
           qflx_ev_soil            =>    waterflux_inst%qflx_ev_soil_col           , & ! Input:  [real(r8) (:)]  evaporation flux from soil (W/m**2) [+ to atm]
           qflx_evap_grnd          =>    waterflux_inst%qflx_evap_grnd_col         , & ! Input:  [real(r8) (:)]  ground surface evaporation rate (mm H2O/s) [+]
           qflx_ev_h2osfc          =>    waterflux_inst%qflx_ev_h2osfc_col         , & ! Input:  [real(r8) (:)]  evaporation flux from h2osfc (W/m**2) [+ to atm]
-
-          qflx_sat_excess_surf    => saturated_excess_runoff_inst%qflx_sat_excess_surf_col, & ! Input:  [real(r8) (:)   ]  surface runoff due to saturated surface (mm H2O /s)
+          qflx_sat_excess_surf    =>    waterflux_inst%qflx_sat_excess_surf_col   , & ! Input:  [real(r8) (:)]  surface runoff due to saturated surface (mm H2O /s)
 
           frac_sno                =>    waterdiagnosticbulk_inst%frac_sno_eff_col          , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)
           frac_h2osfc             =>    waterdiagnosticbulk_inst%frac_h2osfc_col          & ! Input:  [real(r8) (:)   ]  fraction of ground covered by surface water (0 to 1)
@@ -262,7 +259,7 @@ end subroutine SetQflxInputs
 
    !-----------------------------------------------------------------------
    subroutine RouteInfiltrationExcess(bounds, num_hydrologyc, filter_hydrologyc, &
-        waterflux_inst, infiltration_excess_runoff_inst, soilhydrology_inst)
+        waterflux_inst, soilhydrology_inst)
      !
      ! !DESCRIPTION:
      ! Route the infiltration excess runoff flux
@@ -274,7 +271,6 @@ subroutine RouteInfiltrationExcess(bounds, num_hydrologyc, filter_hydrologyc, &
      integer                  , intent(in)    :: num_hydrologyc       ! number of column soil points in column filter
      integer                  , intent(in)    :: filter_hydrologyc(:) ! column filter for soil points
      type(waterflux_type)     , intent(inout) :: waterflux_inst
-     type(infiltration_excess_runoff_type), intent(in) :: infiltration_excess_runoff_inst
      type(soilhydrology_type) , intent(in)    :: soilhydrology_inst
      !
      ! !LOCAL VARIABLES:
@@ -289,8 +285,7 @@ subroutine RouteInfiltrationExcess(bounds, num_hydrologyc, filter_hydrologyc, &
           qflx_infl_excess_surf   => waterflux_inst%qflx_infl_excess_surf_col             , & ! Output: [real(r8) (:)   ] surface runoff due to infiltration excess (mm H2O /s)
           qflx_in_soil            => waterflux_inst%qflx_in_soil_col                      , & ! Input:  [real(r8) (:)   ] surface input to soil (mm/s)
           qflx_top_soil_to_h2osfc => waterflux_inst%qflx_top_soil_to_h2osfc_col           , & ! Input:  [real(r8) (:)   ] portion of qflx_top_soil going to h2osfc, minus evaporation (mm/s)
-
-          qflx_infl_excess        => infiltration_excess_runoff_inst%qflx_infl_excess_col , & ! Input:  [real(r8) (:)   ] infiltration excess runoff (mm H2O /s)
+          qflx_infl_excess        => waterflux_inst%qflx_infl_excess_col                  , & ! Input:  [real(r8) (:)   ] infiltration excess runoff (mm H2O /s)
 
           h2osfcflag              => soilhydrology_inst%h2osfcflag                          & ! Input:  integer
           )
@@ -574,7 +569,7 @@ end subroutine QflxH2osfcDrain
    !-----------------------------------------------------------------------
    subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
         num_urbanc, filter_urbanc, &
-        waterflux_inst, soilhydrology_inst, saturated_excess_runoff_inst, &
+        waterflux_inst, soilhydrology_inst, &
         waterstatebulk_inst)
      !
      ! !DESCRIPTION:
@@ -593,7 +588,6 @@ subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
      integer              , intent(in)    :: filter_urbanc(:)     ! column filter for urban points
      type(waterflux_type) , intent(inout) :: waterflux_inst
      type(soilhydrology_type) , intent(inout) :: soilhydrology_inst
-     type(saturated_excess_runoff_type), intent(in) :: saturated_excess_runoff_inst
      type(waterstatebulk_type), intent(in)    :: waterstatebulk_inst
      !
      ! !LOCAL VARIABLES:
@@ -612,11 +606,10 @@ subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
           qflx_rain_plus_snomelt => waterflux_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:)   ] rain plus snow melt falling on the soil (mm/s)
           qflx_evap_grnd   =>    waterflux_inst%qflx_evap_grnd_col   , & ! Input:  [real(r8) (:)   ]  ground surface evaporation rate (mm H2O/s) [+]    
           qflx_floodc      =>    waterflux_inst%qflx_floodc_col      , & ! Input:  [real(r8) (:)   ]  column flux of flood water from RTM               
+          qflx_sat_excess_surf => waterflux_inst%qflx_sat_excess_surf_col , & ! Input:  [real(r8) (:)   ]  surface runoff due to saturated surface (mm H2O /s)
 
           xs_urban         =>    soilhydrology_inst%xs_urban_col     , & ! Output: [real(r8) (:)   ]  excess soil water above urban ponding limit
 
-          qflx_sat_excess_surf => saturated_excess_runoff_inst%qflx_sat_excess_surf_col , & ! Input:  [real(r8) (:)   ]  surface runoff due to saturated surface (mm H2O /s)
-
           h2osoi_liq       =>    waterstatebulk_inst%h2osoi_liq_col        & ! Input:  [real(r8) (:,:) ]  liquid water (kg/m2)                            
           )
 
diff --git a/src/biogeophys/WaterfluxType.F90 b/src/biogeophys/WaterfluxType.F90
index f2217d36a3..084c3d9ac3 100644
--- a/src/biogeophys/WaterfluxType.F90
+++ b/src/biogeophys/WaterfluxType.F90
@@ -29,6 +29,8 @@ module WaterfluxType
      real(r8), pointer :: qflx_rain_grnd_col       (:)   ! col rain on ground after interception (mm H2O/s) [+]
      real(r8), pointer :: qflx_snow_grnd_patch     (:)   ! patch snow on ground after interception (mm H2O/s) [+]
      real(r8), pointer :: qflx_snow_grnd_col       (:)   ! col snow on ground after interception (mm H2O/s) [+]
+     real(r8), pointer :: qflx_irrig_patch         (:)   ! patch irrigation flux (mm H2O/s) [+]
+     real(r8), pointer :: qflx_irrig_col           (:)   ! col irrigation flux (mm H2O/s) [+]
      real(r8), pointer :: qflx_sub_snow_patch      (:)   ! patch sublimation rate from snow pack (mm H2O /s) [+]
      real(r8), pointer :: qflx_sub_snow_col        (:)   ! col sublimation rate from snow pack (mm H2O /s) [+]
      real(r8), pointer :: qflx_evap_soi_patch      (:)   ! patch soil evaporation (mm H2O/s) (+ = to atm)
@@ -49,6 +51,11 @@ module WaterfluxType
      real(r8), pointer :: qflx_snwcp_ice_col       (:)   ! col excess solid h2o due to snow capping (outgoing) (mm H2O /s)
      real(r8), pointer :: qflx_snwcp_discarded_liq_col(:) ! col excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s)
      real(r8), pointer :: qflx_snwcp_discarded_ice_col(:) ! col excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s)
+     real(r8), pointer :: qflx_glcice_col(:)              ! col net flux of new glacial ice (growth - melt) (mm H2O/s), passed to GLC; only valid inside the do_smb_c filter
+     real(r8), pointer :: qflx_glcice_dyn_water_flux_col(:) ! col water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system); valid for all columns
+     real(r8), pointer :: qflx_glcice_frz_col (:)         ! col ice growth (positive definite) (mm H2O/s); only valid inside the do_smb_c filter
+     real(r8), pointer :: qflx_glcice_melt_col(:)         ! col ice melt (positive definite) (mm H2O/s); only valid inside the do_smb_c filter
+
 
      real(r8), pointer :: qflx_tran_veg_patch      (:)   ! patch vegetation transpiration (mm H2O/s) (+ = to atm)
      real(r8), pointer :: qflx_tran_veg_col        (:)   ! col vegetation transpiration (mm H2O/s) (+ = to atm)
@@ -72,6 +79,8 @@ module WaterfluxType
 
      real(r8), pointer :: qflx_adv_col             (:,:) ! col advective flux across different soil layer interfaces [mm H2O/s] [+ downward]
      real(r8), pointer :: qflx_rootsoi_col         (:,:) ! col root and soil water exchange [mm H2O/s] [+ into root]
+     real(r8), pointer :: qflx_sat_excess_surf_col (:)   ! col surface runoff due to saturated surface (mm H2O /s)
+     real(r8), pointer :: qflx_infl_excess_col     (:)   ! col infiltration excess runoff (mm H2O /s)
      real(r8), pointer :: qflx_infl_col            (:)   ! col infiltration (mm H2O /s)
      real(r8), pointer :: qflx_infl_excess_surf_col(:)   ! col surface runoff due to infiltration excess (mm H2O /s)
      real(r8), pointer :: qflx_h2osfc_surf_col     (:)   ! col surface water runoff (mm H2O /s)
@@ -175,6 +184,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qflx_prec_grnd_patch     (begp:endp))              ; this%qflx_prec_grnd_patch     (:)   = nan
     allocate(this%qflx_rain_grnd_patch     (begp:endp))              ; this%qflx_rain_grnd_patch     (:)   = nan
     allocate(this%qflx_snow_grnd_patch     (begp:endp))              ; this%qflx_snow_grnd_patch     (:)   = nan
+    allocate(this%qflx_irrig_patch         (begp:endp))              ; this%qflx_irrig_patch         (:)   = nan
     allocate(this%qflx_sub_snow_patch      (begp:endp))              ; this%qflx_sub_snow_patch      (:)   = 0.0_r8
     allocate(this%qflx_tran_veg_patch      (begp:endp))              ; this%qflx_tran_veg_patch      (:)   = nan
 
@@ -190,11 +200,16 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qflx_prec_grnd_col       (begc:endc))              ; this%qflx_prec_grnd_col       (:)   = nan
     allocate(this%qflx_rain_grnd_col       (begc:endc))              ; this%qflx_rain_grnd_col       (:)   = nan
     allocate(this%qflx_snow_grnd_col       (begc:endc))              ; this%qflx_snow_grnd_col       (:)   = nan
+    allocate(this%qflx_irrig_col           (begc:endc))              ; this%qflx_irrig_col           (:)   = nan
     allocate(this%qflx_sub_snow_col        (begc:endc))              ; this%qflx_sub_snow_col        (:)   = 0.0_r8
     allocate(this%qflx_snwcp_liq_col       (begc:endc))              ; this%qflx_snwcp_liq_col       (:)   = nan
     allocate(this%qflx_snwcp_ice_col       (begc:endc))              ; this%qflx_snwcp_ice_col       (:)   = nan
     allocate(this%qflx_snwcp_discarded_liq_col(begc:endc))           ; this%qflx_snwcp_discarded_liq_col(:) = nan
     allocate(this%qflx_snwcp_discarded_ice_col(begc:endc))           ; this%qflx_snwcp_discarded_ice_col(:) = nan
+    allocate(this%qflx_glcice_col          (begc:endc))              ; this%qflx_glcice_col          (:)   = nan
+    allocate(this%qflx_glcice_dyn_water_flux_col(begc:endc))         ; this%qflx_glcice_dyn_water_flux_col (:) = nan
+    allocate(this%qflx_glcice_frz_col      (begc:endc))              ; this%qflx_glcice_frz_col      (:)   = nan
+    allocate(this%qflx_glcice_melt_col     (begc:endc))              ; this%qflx_glcice_melt_col     (:)   = nan
     allocate(this%qflx_tran_veg_col        (begc:endc))              ; this%qflx_tran_veg_col        (:)   = nan
     allocate(this%qflx_evap_veg_col        (begc:endc))              ; this%qflx_evap_veg_col        (:)   = nan
     allocate(this%qflx_evap_can_col        (begc:endc))              ; this%qflx_evap_can_col        (:)   = nan
@@ -220,6 +235,8 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qflx_drain_vr_col      (begc:endc,1:nlevsoi))      ; this%qflx_drain_vr_col        (:,:) = nan
     allocate(this%qflx_adv_col             (begc:endc,0:nlevsoi))    ; this%qflx_adv_col             (:,:) = nan
     allocate(this%qflx_rootsoi_col         (begc:endc,1:nlevsoi))    ; this%qflx_rootsoi_col         (:,:) = nan
+    allocate(this%qflx_sat_excess_surf_col (begc:endc))              ; this%qflx_sat_excess_surf_col (:)   = nan
+    allocate(this%qflx_infl_excess_col     (begc:endc))              ; this%qflx_infl_excess_col     (:)   = nan
     allocate(this%qflx_infl_col            (begc:endc))              ; this%qflx_infl_col            (:)   = nan
     allocate(this%qflx_surf_col            (begc:endc))              ; this%qflx_surf_col            (:)   = nan
     allocate(this%qflx_drain_col           (begc:endc))              ; this%qflx_drain_col           (:)   = nan
@@ -460,6 +477,21 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='excess solid h2o due to snow capping not including correction for land use change', &
          ptr_col=this%qflx_snwcp_ice_col, c2l_scale_type='urbanf')
 
+    this%qflx_glcice_col(begc:endc) = spval
+    call hist_addfld1d (fname='QICE',  units='mm/s',  &
+         avgflag='A', long_name='ice growth/melt', &
+         ptr_col=this%qflx_glcice_col, l2g_scale_type='ice')
+
+    this%qflx_glcice_frz_col(begc:endc) = spval
+    call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
+         avgflag='A', long_name='ice growth', &
+         ptr_col=this%qflx_glcice_frz_col, l2g_scale_type='ice')
+
+    this%qflx_glcice_melt_col(begc:endc) = spval
+    call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
+         avgflag='A', long_name='ice melt', &
+         ptr_col=this%qflx_glcice_melt_col, l2g_scale_type='ice')
+
     this%qflx_rain_grnd_patch(begp:endp) = spval
     call hist_addfld1d (fname='QFLX_RAIN_GRND', units='mm H2O/s', &
          avgflag='A', long_name='rain on ground after interception', &
@@ -470,6 +502,11 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='snow on ground after interception', &
          ptr_patch=this%qflx_snow_grnd_patch, default='inactive', c2l_scale_type='urbanf')
 
+    this%qflx_irrig_patch(begp:endp) = spval
+    call hist_addfld1d (fname='QIRRIG', units='mm/s', &
+         avgflag='A', long_name='water added through irrigation', &
+         ptr_patch=this%qflx_irrig_patch)
+
     this%qflx_evap_grnd_patch(begp:endp) = spval
     call hist_addfld1d (fname='QFLX_EVAP_GRND', units='mm H2O/s', &
          avgflag='A', long_name='ground surface evaporation', &
@@ -695,6 +732,13 @@ subroutine InitCold(this, bounds)
     ! the sake of columns outside this filter
     this%qflx_ice_runoff_xs_col(bounds%begc:bounds%endc) = 0._r8
 
+    ! Initialize qflx_glcice_dyn_water_flux_col to 0 for all columns because we want this
+    ! flux to remain 0 for columns where is is never set, including non-glacier columns.
+    !
+    ! Other qflx_glcice fluxes intentionally remain unset (spval) outside the do_smb
+    ! filter, so that they are flagged as missing value outside that filter.
+    this%qflx_glcice_dyn_water_flux_col(bounds%begc:bounds%endc) = 0._r8
+
     this%AnnEt(bounds%begc:bounds%endc)                 = 0._r8
   
     ! needed for CNNLeaching 
diff --git a/src/biogeophys/test/Irrigation_test/test_irrigation.pf b/src/biogeophys/test/Irrigation_test/test_irrigation.pf
index 3c599aac02..99c5726b52 100644
--- a/src/biogeophys/test/Irrigation_test/test_irrigation.pf
+++ b/src/biogeophys/test/Irrigation_test/test_irrigation.pf
@@ -8,6 +8,7 @@ module test_irrigation
   use shr_kind_mod, only : r8 => shr_kind_r8
   use clm_varpar, only : nlevsoi, nlevgrnd
   use landunit_varcon, only : istsoil
+  use WaterfluxType, only : waterflux_type
   use PatchType     , only : patch
   use ColumnType    , only : col
   use GridcellType  , only : grc
@@ -26,6 +27,7 @@ module test_irrigation
      integer :: numf
      integer, allocatable :: filter(:)
      type(irrigation_type) :: irrigation
+     type(waterflux_type) :: waterflux
 
      ! Irrigation parameters
      type(irrigation_params_type) :: irrigation_params
@@ -84,6 +86,8 @@ contains
     class(TestIrrigation), intent(inout) :: this
 
     call this%irrigation%Clean()
+    deallocate(this%waterflux%qflx_irrig_patch)
+    deallocate(this%waterflux%qflx_irrig_col)
     call this%teardownEnvironment()
     call unittest_subgrid_teardown()
   end subroutine tearDown
@@ -178,6 +182,10 @@ contains
          irrig_river_volume_threshold = l_irrig_river_volume_threshold, &
          limit_irrigation_if_rof_enabled = limit_irrigation_if_rof_enabled)
 
+    ! Allocate fluxes output from irrigation routines
+    allocate(this%waterflux%qflx_irrig_patch(bounds%begp:bounds%endp))
+    allocate(this%waterflux%qflx_irrig_col(bounds%begc:bounds%endc))
+
     ! Set inputs to irrigation routines
     allocate(this%elai(bounds%begp:bounds%endp), source=10._r8)
     allocate(this%t_soisno(bounds%begc:bounds%endc, nlevgrnd), source=1000._r8)
@@ -333,7 +341,7 @@ contains
          volr = this%volr, &
          rof_prognostic = .true.)
     
-    call this%irrigation%ApplyIrrigation(bounds)
+    call this%irrigation%ApplyIrrigation(bounds, this%waterflux)
 
   end subroutine calculateAndApplyIrrigation
 
@@ -361,7 +369,7 @@ contains
     ! Check result
     call this%computeDeficits(deficits)
     expected = sum(deficits(bounds%begp,1:nlevsoi)) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
 
   end subroutine irrigation_flux_is_correct
 
@@ -378,7 +386,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
   end subroutine no_irrigation_for_wet_soil
 
   @Test
@@ -408,7 +416,7 @@ contains
     ! This first assertion makes sure the test has been set up reasonably - to give a net deficit
     @assertLessThan(0._r8, expected)
     ! Here is the main assertion:
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
   end subroutine surplus_offsets_deficit
 
   @Test
@@ -426,7 +434,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine no_irrigation_for_unirrigated_pfts
 
@@ -443,7 +451,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine no_irrigation_for_lai0
 
@@ -470,7 +478,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
   end subroutine no_irrigation_for_soil_moisture_above_threshold
 
   @Test
@@ -487,7 +495,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine no_irrigation_at_wrong_time
 
@@ -519,7 +527,7 @@ contains
 
     ! Check result
     expected = total_deficit / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
   end subroutine unlimited_irrigation_for_non_limiting_volr
 
   @Test
@@ -552,7 +560,7 @@ contains
     expected = ((this%volr(begg) * (1._r8 - irrig_river_volume_threshold)) &
                 / grc%area(begg) * m3_over_km2_to_mm) / &
                this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
   end subroutine limited_irrigation_for_limiting_volr
 
   @Test
@@ -572,7 +580,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertTrue(this%irrigation%qflx_irrig_patch(bounds%begp) > 0._r8)
+    @assertTrue(this%waterflux%qflx_irrig_patch(bounds%begp) > 0._r8)
 
   end subroutine irrigation_should_happen_for_big_longitude
 
@@ -597,7 +605,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
 
   end subroutine irrigation_continues_at_same_rate_for_multiple_time_steps
 
@@ -620,11 +628,11 @@ contains
        call this%calculateAndApplyIrrigation()
        this%time_prev = this%time_prev + dtime
     end do
-    @assertTrue(this%irrigation%qflx_irrig_patch(bounds%begp) > 0._r8)
+    @assertTrue(this%waterflux%qflx_irrig_patch(bounds%begp) > 0._r8)
 
     ! Ensure that irrigation flux goes to 0 in the following time step
     call this%calculateAndApplyIrrigation()
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine irrigation_continues_for_correct_number_of_time_steps
   
@@ -653,7 +661,7 @@ contains
        this%time_prev = this%time_prev + dtime
     end do
     ! The following assertion is mainly here to make sure the test is working as intended
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
     ! Now reset time, change soil moisture, and make sure that irrigation happens as expected
     this%time_prev = time_prev_orig
@@ -664,7 +672,7 @@ contains
     ! Make sure that the test has been set up reasonably - to give a net deficit
     @assertLessThan(0._r8, expected)
     ! Here's the main assertion:
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
 
   end subroutine irrigation_flux_is_correct_on_second_day
     
@@ -687,7 +695,7 @@ contains
     call this%computeDeficits(deficits)
     ! Now on to the real assertion
     expected = sum(deficits(bounds%begp,1:nlevirrig)) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
   end subroutine irrigation_excludes_deep_layers
 
   @Test
@@ -708,7 +716,7 @@ contains
     ! Check result
     call this%computeDeficits(deficits)
     expected = sum(deficits(bounds%begp,1:(nlevsoi-1))) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
   end subroutine irrigation_excludes_bedrock_layers
 
   @Test
@@ -724,7 +732,7 @@ contains
     call this%calculateAndApplyIrrigation()
 
     ! Check result
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine no_irrigation_for_frozen_soil
 
@@ -746,7 +754,7 @@ contains
     call this%computeDeficits(deficits)
     ! Only include deficit from top layer, since 2nd layer is frozen
     expected = deficits(bounds%begp, 1) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp))
 
   end subroutine no_irrigation_below_frozen_soil_layer
 
@@ -783,9 +791,9 @@ contains
     ! Check result
     call this%computeDeficits(deficits)
     expected1 = sum(deficits(bounds%begp,1:nlevsoi)) / this%irrigation_params%irrig_length
-    @assertEqual(expected1, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected1, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
     expected2 = sum(deficits(bounds%endp,1:nlevsoi)) / this%irrigation_params%irrig_length
-    @assertEqual(expected2, this%irrigation%qflx_irrig_patch(bounds%endp), tolerance=tol)
+    @assertEqual(expected2, this%waterflux%qflx_irrig_patch(bounds%endp), tolerance=tol)
     ! Make sure this test had some power, by ensuring that the two points differ:
     @assertTrue(expected1 /= expected2)
 
@@ -830,7 +838,7 @@ contains
     total_deficit = sum(deficits(bounds%begp+1, 1:nlevsoi))
     col_total_deficit = total_deficit * wt_irrig
     expected_col_flux = col_total_deficit / this%irrigation_params%irrig_length
-    @assertEqual(expected_col_flux, this%irrigation%qflx_irrig_col(bounds%begc), tolerance=tol)
+    @assertEqual(expected_col_flux, this%waterflux%qflx_irrig_col(bounds%begc), tolerance=tol)
   end subroutine multiple_patches_per_column
 
   @Test
@@ -860,10 +868,10 @@ contains
 
     ! Check result
     ! Irrigation happens within filter
-    @assertTrue(this%irrigation%qflx_irrig_patch(bounds%begp + 1) > 0._r8)
+    @assertTrue(this%waterflux%qflx_irrig_patch(bounds%begp + 1) > 0._r8)
     ! Irrigation does NOT happen outside filter
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%endp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%endp))
     
   end subroutine irrigation_only_happens_within_filter
 
@@ -904,13 +912,13 @@ contains
     ! Check result
     call this%computeDeficits(deficits)
     ! First patch should have no irrigation, because soil is all frozen
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp))
     ! Second patch should have irrigation just based on top layer, because 2nd layer is frozen
     expected = deficits(bounds%begp+1, 1) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp+1), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%begp+1), tolerance=tol)
     ! Third patch should have irrigation from all layers
     expected = sum(deficits(bounds%endp,1:nlevsoi)) / this%irrigation_params%irrig_length
-    @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%endp), tolerance=tol)
+    @assertEqual(expected, this%waterflux%qflx_irrig_patch(bounds%endp), tolerance=tol)
 
   end subroutine test_multiple_patches_with_different_frozen_soil
 
@@ -965,12 +973,12 @@ contains
 
     ! Check result
     expected1 = total_deficit1*volr_fraction / this%irrigation_params%irrig_length
-    @assertEqual(expected1, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol)
+    @assertEqual(expected1, this%waterflux%qflx_irrig_patch(bounds%begp), tolerance=tol)
     expected2 = total_deficit2*volr_fraction / this%irrigation_params%irrig_length
-    @assertEqual(expected2, this%irrigation%qflx_irrig_patch(bounds%begp+1), tolerance=tol)
+    @assertEqual(expected2, this%waterflux%qflx_irrig_patch(bounds%begp+1), tolerance=tol)
     ! Check of unirrigated patch isn't central to this test, but we might as well check
     ! it while we're at it:
-    @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp+2))
+    @assertEqual(0._r8, this%waterflux%qflx_irrig_patch(bounds%begp+2))
   end subroutine volr_limiting_with_multiple_columns
 
 end module test_irrigation
diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90
index 6266a13f1f..f43aad694d 100644
--- a/src/main/clm_driver.F90
+++ b/src/main/clm_driver.F90
@@ -424,7 +424,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
 
        ! Irrigation flux
        ! input is main channel storage
-       call irrigation_inst%ApplyIrrigation(bounds_clump)
+       call irrigation_inst%ApplyIrrigation(bounds_clump, waterflux_inst)
        call t_stopf('drvinit')
 
        ! ============================================================================
@@ -440,8 +440,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
             filter(nc)%num_nolakec, filter(nc)%nolakec, &
             filter(nc)%num_nolakep, filter(nc)%nolakep, &
             atm2lnd_inst, canopystate_inst, temperature_inst, &
-            aerosol_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterflux_inst, &
-            irrigation_inst)
+            aerosol_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterflux_inst)
        call t_stopf('canhydro')
 
        ! ============================================================================
@@ -653,7 +652,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        ! converted back to solid ice as soon as possible
        call glacier_smb_inst%HandleIceMelt(bounds_clump, &
             filter(nc)%num_do_smb_c, filter(nc)%do_smb_c, &
-            waterstatebulk_inst)
+            waterstatebulk_inst, waterflux_inst)
        call t_stopf('soiltemperature')
 
        ! ============================================================================
@@ -853,7 +852,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
             filter(nc)%num_do_smb_c, filter(nc)%do_smb_c,     &                
             atm2lnd_inst, glc2lnd_inst, temperature_inst,     &
             soilhydrology_inst, soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, &
-            irrigation_inst, glacier_smb_inst)
+            glacier_smb_inst)
 
        call t_stopf('hydro2_drainage')     
 
@@ -927,7 +926,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        call t_startf('balchk')
        call BalanceCheck(bounds_clump, &
             atm2lnd_inst, solarabs_inst, waterflux_inst, &
-            waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, irrigation_inst, glacier_smb_inst, &
+            waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, &
             energyflux_inst, canopystate_inst)
        call t_stopf('balchk')
 
@@ -1041,7 +1040,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
 
     call lnd2atm(bounds_proc,                                            &
          atm2lnd_inst, surfalb_inst, temperature_inst, frictionvel_inst, &
-         waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, irrigation_inst, energyflux_inst, &
+         waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, energyflux_inst, &
          solarabs_inst, drydepvel_inst,       &
          vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, &
          lnd2atm_inst, &
@@ -1059,7 +1058,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        call get_clump_bounds(nc, bounds_clump)
        call lnd2glc_inst%update_lnd2glc(bounds_clump,       &
             filter(nc)%num_do_smb_c, filter(nc)%do_smb_c,   &
-            temperature_inst, glacier_smb_inst, topo_inst,    &
+            temperature_inst, waterflux_inst, topo_inst,    &
             init=.false.)           
     end do
     !$OMP END PARALLEL DO
diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90
index 9ec58fd2bf..eddf5421c7 100644
--- a/src/main/clm_initializeMod.F90
+++ b/src/main/clm_initializeMod.F90
@@ -634,7 +634,7 @@ subroutine initialize2( )
        call t_startf('init_lnd2glc')
        call lnd2glc_inst%update_lnd2glc(bounds_clump,       &
             filter(nc)%num_do_smb_c, filter(nc)%do_smb_c,   &
-            temperature_inst, glacier_smb_inst, topo_inst, &
+            temperature_inst, waterflux_inst, topo_inst, &
             init=.true.)
        call t_stopf('init_lnd2glc')
     end do
diff --git a/src/main/lnd2atmMod.F90 b/src/main/lnd2atmMod.F90
index 3c52e3febb..65cf55a7b4 100644
--- a/src/main/lnd2atmMod.F90
+++ b/src/main/lnd2atmMod.F90
@@ -33,7 +33,6 @@ module lnd2atmMod
   use WaterStateBulkType       , only : waterstatebulk_type
   use WaterDiagnosticBulkType       , only : waterdiagnosticbulk_type
   use WaterBalanceType       , only : waterbalance_type
-  use IrrigationMod        , only : irrigation_type 
   use glcBehaviorMod       , only : glc_behavior_type
   use glc2lndMod           , only : glc2lnd_type
   use ColumnType           , only : col
@@ -125,7 +124,7 @@ end subroutine lnd2atm_minimal
   !------------------------------------------------------------------------
   subroutine lnd2atm(bounds, &
        atm2lnd_inst, surfalb_inst, temperature_inst, frictionvel_inst, &
-       waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, irrigation_inst, energyflux_inst, &
+       waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalance_inst, waterflux_inst, energyflux_inst, &
        solarabs_inst, drydepvel_inst,  &
        vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, &
        lnd2atm_inst, &
@@ -147,7 +146,6 @@ subroutine lnd2atm(bounds, &
     type(waterdiagnosticbulk_type)       , intent(inout) :: waterdiagnosticbulk_inst
     type(waterbalance_type)       , intent(inout) :: waterbalance_inst
     type(waterflux_type)        , intent(inout) :: waterflux_inst
-    type(irrigation_type)       , intent(in)    :: irrigation_inst
     type(energyflux_type)       , intent(in)    :: energyflux_inst
     type(solarabs_type)         , intent(in)    :: solarabs_inst
     type(drydepvel_type)        , intent(in)    :: drydepvel_inst
@@ -370,7 +368,7 @@ subroutine lnd2atm(bounds, &
          c2l_scale_type= 'urbanf', l2g_scale_type='unity' )
 
     call c2g( bounds, &
-         irrigation_inst%qflx_irrig_col (bounds%begc:bounds%endc), &
+         waterflux_inst%qflx_irrig_col (bounds%begc:bounds%endc), &
          lnd2atm_inst%qirrig_grc(bounds%begg:bounds%endg), &
          c2l_scale_type= 'urbanf', l2g_scale_type='unity' )
 
diff --git a/src/main/lnd2glcMod.F90 b/src/main/lnd2glcMod.F90
index 9de7eba3f3..422c8af14c 100644
--- a/src/main/lnd2glcMod.F90
+++ b/src/main/lnd2glcMod.F90
@@ -25,8 +25,8 @@ module lnd2glcMod
   use column_varcon   , only : col_itype_to_icemec_class
   use landunit_varcon , only : istice_mec, istsoil
   use abortutils      , only : endrun
-  use GlacierSurfaceMassBalanceMod, only : glacier_smb_type
   use TemperatureType , only : temperature_type
+  use WaterfluxType   , only : waterflux_type
   use LandunitType    , only : lun                
   use ColumnType      , only : col
   use TopoMod         , only : topo_type
@@ -146,7 +146,7 @@ end subroutine InitHistory
 
   !------------------------------------------------------------------------------
   subroutine update_lnd2glc(this, bounds, num_do_smb_c, filter_do_smb_c, &
-       temperature_inst, glacier_smb_inst, topo_inst, init)
+       temperature_inst, waterflux_inst, topo_inst, init)
     !
     ! !DESCRIPTION:
     ! Assign values to lnd2glc+
@@ -157,7 +157,7 @@ subroutine update_lnd2glc(this, bounds, num_do_smb_c, filter_do_smb_c, &
     integer                , intent(in)    :: num_do_smb_c       ! number of columns in filter_do_smb_c
     integer                , intent(in)    :: filter_do_smb_c(:) ! column filter: columns where smb calculations are performed
     type(temperature_type) , intent(in)    :: temperature_inst
-    type(glacier_smb_type) , intent(in)    :: glacier_smb_inst
+    type(waterflux_type)   , intent(in)    :: waterflux_inst
     type(topo_type)        , intent(in)    :: topo_inst
     logical                , intent(in)    :: init               ! if true=>only set a subset of fields
     !
@@ -221,7 +221,7 @@ subroutine update_lnd2glc(this, bounds, num_do_smb_c, filter_do_smb_c, &
       this%tsrf_grc(g,n) = temperature_inst%t_soisno_col(c,1)
       this%topo_grc(g,n) = topo_inst%topo_col(c)
       if (.not. init) then
-         this%qice_grc(g,n) = glacier_smb_inst%qflx_glcice_col(c) * flux_normalization
+         this%qice_grc(g,n) = waterflux_inst%qflx_glcice_col(c) * flux_normalization
 
          ! Check for bad values of qice
          if ( abs(this%qice_grc(g,n)) > 1.0_r8) then