! Copyright (c) 2016-2018,  Los Alamos National Security, LLC (LANS)
! and the University Corporation for Atmospheric Research (UCAR).
!
! Unless noted otherwise source code is licensed under the BSD license.
! Additional copyright and license information can be found in the LICENSE file
! distributed with this code, or at http://mpas-dev.github.io/license.html
!

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
!  li_SGH_driver
!
!> \brief MPAS land ice SGH primary routines
!> \author Matt Hoffman
!> \date   27 June 2016
!> \details
!>  This module contains the main driver routines for
!>  for subglacial hydro.
!
!-----------------------------------------------------------------------

module li_subglacial_hydro

   use mpas_derived_types
   use mpas_pool_routines
   use mpas_dmpar
   use mpas_timekeeping
   use mpas_timer
   use mpas_log

   use li_setup
   use li_constants
   use li_mask

   implicit none
   private

   !--------------------------------------------------------------------
   !
   ! Public parameters
   !
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   !
   ! Public member functions
   !
   !--------------------------------------------------------------------

   public :: li_SGH_init, &
             li_SGH_solve, &
             li_SGH_finalize

   !--------------------------------------------------------------------
   !
   ! Private module variables
   !
   !--------------------------------------------------------------------


!***********************************************************************
   contains


!***********************************************************************
!
!  routine li_SGH_init
!
!> \brief   Initialize SGH
!> \author  Matt Hoffman
!> \date    27 June 2016
!> \details
!>  This routine initializes the subglacial hydro model
!-----------------------------------------------------------------------
   subroutine li_SGH_init(domain, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain  !< Input/Output: domain object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      logical, pointer :: config_SGH
      logical, pointer :: config_ocean_connection_N
      type (block_type), pointer :: block
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: velocityPool
      real (kind=RKIND), pointer :: deltatSGH
      real (kind=RKIND), dimension(:), pointer :: waterThickness
      real (kind=RKIND), dimension(:), pointer :: tillWaterThickness
      real (kind=RKIND), dimension(:), pointer :: waterPressure
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      integer, dimension(:), pointer :: cellMask
      real (kind=RKIND), pointer :: tillMax
      real (kind=RKIND), pointer :: rhoi, rhoo
      real (kind=RKIND), pointer :: config_sea_level
      integer, pointer :: config_num_halos
      integer :: err_tmp


      err = 0
      err_tmp = 0

      call mpas_pool_get_config(liConfigs, 'config_SGH', config_SGH)
      if (.not. config_SGH) then
         call mpas_pool_get_config(liConfigs, 'config_ocean_connection_N', config_ocean_connection_N)
         if (config_ocean_connection_N) then
            call ocean_connection_N(domain)
         endif
         ! If SGH is not active, skip everything
         return
      endif


      call mpas_timer_start("hydro init")

      call mpas_log_write('Beginning subglacial hydro init.')

      ! Check number of halos
      call mpas_pool_get_config(liConfigs, 'config_num_halos', config_num_halos)
      if (config_num_halos < 1) then
         call mpas_log_write("Subglacial hydrology requires config_num_halos >= 1", MPAS_LOG_ERR)
         err = ior(err, 1)
      endif

      call mpas_pool_get_config(liConfigs, 'config_SGH_till_max', tillMax)
      call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
      call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)
      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)

      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)

         call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)

         ! Until init is done properly, make this tiny.  It will be updated at the end of the first subcycle.
         ! TODO: Set time step appropriately on first subcycle of init
         deltatSGH = 1.0e-4_RKIND ! in seconds

         ! Mask needs to be initialized for pressure calcs to be correct
         call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
         err = ior(err, err_tmp)

         call calc_hydro_mask(domain)

         ! remove invalid values - not necessary on restart, but shouldn't hurt
         call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
         waterThickness = max(0.0_RKIND, waterThickness)

         call mpas_pool_get_array(hydroPool, 'tillWaterThickness', tillWaterThickness)
         tillWaterThickness = max(0.0_RKIND, tillWaterThickness)
         tillWaterThickness = min(tillMax, tillWaterThickness)

         call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
         call mpas_pool_get_array(geometryPool, 'thickness', thickness)
         waterPressure = max(0.0_RKIND, waterPressure)
         waterPressure = min(waterPressure, rhoi * gravity * thickness)
         ! set pressure correctly under floating ice and open ocean
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
         call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
         where ( (li_mask_is_floating_ice(cellMask)) .or. &
                 ( (.not. li_mask_is_ice(cellMask)) .and. (bedTopography < config_sea_level) ) )
            waterPressure = rhoo * gravity * (config_sea_level - bedTopography)
         end where

         ! Initialize diagnostic pressure variables
         call calc_pressure_diag_vars(block, err_tmp)
         err = ior(err, err_tmp)

         block => block % next
      end do

      ! === error check
      if (err > 0) then
          call mpas_log_write("An error has occurred in li_SGH_init.", MPAS_LOG_ERR)
      endif

      call mpas_timer_stop("hydro init")

   !--------------------------------------------------------------------
   end subroutine li_SGH_init



!***********************************************************************
!
!  routine li_SGH_solve
!
!> \brief   Solve and update SGH for current time step
!> \author  Matt Hoffman
!> \date    27 June 2016
!> \details
!>  This routine solves and updates the subglacial hydro model
!>  for the current ISM time step.
!-----------------------------------------------------------------------
   subroutine li_SGH_solve(domain, err)
      use mpas_vector_reconstruction
      use li_diagnostic_vars

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain  !< Input/Output: domain object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      logical, pointer :: config_SGH
      logical, pointer :: config_ocean_connection_N
      logical, pointer :: config_SGH_chnl_active
      character (len=StrKIND), pointer :: config_SGH_basal_melt
      type (block_type), pointer :: block
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: thermalPool
      type (mpas_pool_type), pointer :: velocityPool
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:,:), pointer :: temperature, flowParamA
      real (kind=RKIND), dimension(:), pointer :: Wtill, WtillOld
      real (kind=RKIND), dimension(:), pointer :: basalMeltInput
      real (kind=RKIND), dimension(:), pointer :: groundedBasalMassBal
      real (kind=RKIND), dimension(:), pointer :: basalHeatFlux
      real (kind=RKIND), dimension(:), pointer :: basalFrictionFlux
      real (kind=RKIND), dimension(:), pointer :: externalWaterInput
      real (kind=RKIND), dimension(:), pointer :: divergence
      real (kind=RKIND), dimension(:), pointer :: waterFlux
      real (kind=RKIND), dimension(:), pointer :: waterThickness
      real (kind=RKIND), dimension(:), pointer :: waterThicknessOld
      real (kind=RKIND), dimension(:), pointer :: waterThicknessTendency
      real (kind=RKIND), dimension(:), pointer :: divergenceChannel
      real (kind=RKIND), dimension(:), pointer :: channelAreaChangeCell
      real (kind=RKIND), dimension(:), pointer :: dvEdge
      real (kind=RKIND), dimension(:), pointer :: areaCell
      real (kind=RKIND), dimension(:), pointer :: waterVelocity
      real (kind=RKIND), dimension(:), pointer :: waterVelocityCellX
      real (kind=RKIND), dimension(:), pointer :: waterVelocityCellY
      real (kind=RKIND), dimension(:), allocatable :: cellJunk
      integer, dimension(:), pointer :: nEdgesOnCell
      integer, dimension(:,:), pointer :: edgesOnCell
      integer, dimension(:,:), pointer :: cellsOnEdge
      integer, dimension(:,:), pointer :: edgeSignOnCell
      integer, dimension(:), pointer :: cellMask
      real (kind=RKIND), pointer :: deltatSGH
      real (kind=RKIND), pointer :: masterDeltat
      real (kind=RKIND), pointer :: Cd
      real (kind=RKIND), pointer :: tillMax
      integer, pointer :: nCellsSolve
      integer, pointer :: nCells
      integer :: iCell, iEdge, iEdgeOnCell
      real (kind=RKIND) :: timeLeft ! subcycling time remaining until master dt is reached
      integer :: numSubCycles ! number of subcycles
      integer :: err_tmp



      err = 0
      err_tmp = 0

      call mpas_pool_get_config(liConfigs, 'config_SGH', config_SGH)
      if (.not. config_SGH) then
         call mpas_pool_get_config(liConfigs, 'config_ocean_connection_N', config_ocean_connection_N)
         if (config_ocean_connection_N) then
            call ocean_connection_N(domain)
         endif
         ! If SGH is not active, skip everything
         return
      endif

      call calc_hydro_mask(domain)

      call mpas_log_write('Beginning subglacial hydro solve.')
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_active', config_SGH_chnl_active)
      call mpas_pool_get_config(liConfigs, 'config_SGH_till_drainage', Cd)
      call mpas_pool_get_config(liConfigs, 'config_SGH_till_max', tillMax)
      call mpas_pool_get_config(liConfigs, 'config_SGH_basal_melt', config_SGH_basal_melt)

      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'thermal', thermalPool)
         call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_array(thermalPool, 'temperature', temperature)
         call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
         call mpas_pool_get_array(geometryPool, 'thickness', thickness)

         call li_calculate_flowParamA(meshPool, temperature, thickness, flowParamA, err_tmp)
         err = ior(err, err_tmp)

         block => block % next
      end do

      ! initialize while loop
      call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)  ! can get from any block
      call mpas_pool_get_array(meshPool, 'deltat', masterDeltat)
      timeLeft = masterDeltat  ! in seconds
      numSubCycles = 0
      ! =============
      ! =============
      ! =============
      ! subcycle hydro model until master dt is reached
      ! =============
      ! =============
      ! =============
      timecycling: do while (timeLeft > 0.0_RKIND)
      numSubCycles = numSubCycles + 1


      ! =============
      ! Calculate time-varying forcing, as needed
      ! =============
      block => domain % blocklist
      do while (associated(block))

         ! Decide where the basal melt input comes from
         ! Either prescribed, from the temperature model, or calculated here
         if (trim(config_SGH_basal_melt) == 'file') then
             ! do nothing
         elseif (trim(config_SGH_basal_melt) == 'thermal') then
             call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
             call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
             call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
             call mpas_pool_get_array(geometryPool,'groundedBasalMassBal',groundedBasalMassBal)
             basalMeltInput = -1.0_RKIND * groundedBasalMassBal ! TODO: Ensure positive flux?
         elseif (trim(config_SGH_basal_melt) == 'basal_heat') then
             call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
             call mpas_pool_get_subpool(block % structs, 'thermal', thermalPool)
             call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
             call mpas_pool_get_array(thermalPool, 'basalHeatFlux', basalHeatFlux)
             call mpas_pool_get_array(thermalPool, 'basalFrictionFlux', basalFrictionFlux)
             basalMeltInput = (basalHeatFlux + basalFrictionFlux) / latent_heat_ice
         else
            call mpas_log_write("Unknown value provided for config_SGH_basal_melt: " // trim(config_SGH_basal_melt), MPAS_LOG_ERR)
            err = ior(err, 1)
         endif

         ! SHMIP forcing will override basalMeltInput.
         call shmip_timevarying_forcing(block, err_tmp)
         err = ior(err, err_tmp)

         block => block % next
      end do

      ! Perform halo update, if needed
      if ((trim(config_SGH_basal_melt) == 'thermal') .or. (trim(config_SGH_basal_melt) == 'basal_heat')) then
         call mpas_timer_start("halo updates")
         ! (groundedBasalMassBal will not have had a halo update,
         !  because a halo update for it is unneeded except for the SGH model.)
         call mpas_dmpar_field_halo_exch(domain, 'basalMeltInput')
         call mpas_timer_stop("halo updates")
      endif


      ! =============
      ! Update till water layer thickness
      ! =============
      block => domain % blocklist
      do while (associated(block))

         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)

         call mpas_pool_get_array(hydroPool, 'tillWaterThickness', Wtill)
         call mpas_pool_get_array(hydroPool, 'tillWaterThicknessOld', WtillOld)
         call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
         call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
         call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

         WtillOld = Wtill
         Wtill = Wtill + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - Cd)
         Wtill = Wtill * li_mask_is_grounded_ice_int(cellMask)  ! zero Wtill in non-grounded locations
         Wtill = min(Wtill, tillmax)
         Wtill = max(0.0_RKIND, Wtill)

         block => block % next
      end do


      ! =============
      ! Calculate edge quantities and advective fluxes
      ! =============
      block => domain % blocklist
      do while (associated(block))

         call calc_edge_quantities(block, err_tmp)
         err = ior(err, err_tmp)

         block => block % next
      end do
      ! Update halos on edge quantities
      call mpas_timer_start("halo updates")
      call mpas_dmpar_field_halo_exch(domain, 'waterFlux')
      ! intermediate fields will be out of date, but will be correct in output files
      ! waterVelocity needs to be updated in order to get waterVelocityCellX/Y correct
      call mpas_dmpar_field_halo_exch(domain, 'waterVelocity')
      call mpas_timer_stop("halo updates")

      ! Calculate reconstructed velocities on cell centers for viz
      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
         call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity)
         call mpas_pool_get_array(hydroPool, 'waterVelocityCellX', waterVelocityCellX)
         call mpas_pool_get_array(hydroPool, 'waterVelocityCellY', waterVelocityCellY)
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
         allocate(cellJunk(nCells+1))
         call mpas_reconstruct(meshPool, waterVelocity, waterVelocityCellX, waterVelocityCellY, &
              cellJunk, cellJunk, cellJunk)
         deallocate(cellJunk)
         waterVelocityCellX = waterVelocityCellX * li_mask_is_grounded_ice_int(cellMask)  ! zero in non-grounded locations
         waterVelocityCellY = waterVelocityCellY * li_mask_is_grounded_ice_int(cellMask)  ! zero in non-grounded locations

         block => block % next
      end do


      ! =============
      ! Update Channel fields
      ! =============
      if (config_SGH_chnl_active) then
         block => domain % blocklist
         do while (associated(block))

            call update_channel(block, err_tmp)
            err = ior(err, err_tmp)

            block => block % next
         end do
         ! Update halos on channel
         call mpas_timer_start("halo updates")
         call mpas_dmpar_field_halo_exch(domain, 'channelChangeRate')
         call mpas_dmpar_field_halo_exch(domain, 'channelDischarge')
         call mpas_timer_stop("halo updates")
      endif


      ! =============
      ! Calculate adaptive time step
      ! =============
      call check_timestep(domain, timeLeft, numSubCycles, err_tmp)
      err = ior(err, err_tmp)
      if (err > 0) then
         exit timecycling
      endif


      ! =============
      ! Compute flux divergence
      ! =============
      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_array(hydroPool, 'divergence', divergence)
         call mpas_pool_get_array(hydroPool, 'waterFlux', waterFlux)
         call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
         call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
         call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
         call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
         call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
         call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
         call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

         divergence(:) = 0.0_RKIND  ! zero div before starting
         ! loop over locally owned cells
         do iCell = 1, nCellsSolve
            if (li_mask_is_grounded_ice(cellMask(iCell))) then  ! can skip for non-grounded ice
               ! compute fluxes for each edge of the cell
               do iEdgeOnCell = 1, nEdgesOnCell(iCell)
                  iEdge = edgesOnCell(iEdgeOnCell, iCell)
                  ! add on advective & diffusive fluxes
                  divergence(iCell) = divergence(iCell) - waterFlux(iEdge) * dvEdge(iEdge) * edgeSignOnCell(iEdgeOnCell, iCell)
               end do ! edges
            end if
         end do ! cells
         divergence(1:nCellsSolve) = divergence(1:nCellsSolve) / areaCell(1:nCellsSolve)

         block => block % next
      end do
      ! Update halos on divergence
      call mpas_timer_start("halo updates")
      call mpas_dmpar_field_halo_exch(domain, 'divergence')
      call mpas_timer_stop("halo updates")


      ! =============
      ! Update channel area now that we have time step
      ! =============
      if (config_SGH_chnl_active) then
         block => domain % blocklist
         do while (associated(block))

            call evolve_channel(block, err_tmp)
            err = ior(err, err_tmp)

            block => block % next
         end do
         call mpas_timer_start("halo updates")
         call mpas_dmpar_field_halo_exch(domain, 'divergenceChannel')
         call mpas_dmpar_field_halo_exch(domain, 'channelAreaChangeCell')
         call mpas_dmpar_field_halo_exch(domain, 'channelArea')
         call mpas_timer_stop("halo updates")
      endif


      ! =============
      ! Calculate pressure field
      ! =============
      block => domain % blocklist
      do while (associated(block))

         call calc_pressure(block, err_tmp)
         err = ior(err, err_tmp)

         block => block % next
      end do


      ! =============
      ! Update water layer thickness
      ! =============
      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
         call mpas_pool_get_array(hydroPool, 'waterThicknessOld', waterThicknessOld)
         call mpas_pool_get_array(hydroPool, 'waterThicknessTendency', waterThicknessTendency)
         call mpas_pool_get_array(hydroPool, 'tillWaterThickness', Wtill)
         call mpas_pool_get_array(hydroPool, 'tillWaterThicknessOld', WtillOld)
         call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
         call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
         call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
         call mpas_pool_get_array(hydroPool, 'divergence', divergence)
         call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
         call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

         waterThicknessOld = waterThickness
         waterThickness = waterThicknessOld + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - divergence  &
             - divergenceChannel - channelAreaChangeCell  &
             - (Wtill - WtillOld) / deltatSGH)
         waterThickness = waterThickness * li_mask_is_grounded_ice_int(cellMask)  ! zero in non-grounded locations
         waterThickness = max(0.0_RKIND, waterThickness)
         divergence = divergence * li_mask_is_grounded_ice_int(cellMask)  ! zero in non-grounded locations for more convenient viz
         waterThicknessTendency = (waterThickness - waterThicknessOld) / deltatSGH

         block => block % next
      end do


      ! =============
      ! =============
      ! =============
      end do timecycling ! while timeLeft>0
      ! =============
      ! =============
      ! =============

      call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))


      ! === error check
      if (err > 0) then
          call mpas_log_write("An error has occurred in li_SGH_solve.", MPAS_LOG_ERR)
      endif

   !--------------------------------------------------------------------
   end subroutine li_SGH_solve



!***********************************************************************
!
!  routine li_SGH_finalize
!
!> \brief   Finalize SGH
!> \author  Matt Hoffman
!> \date    27 June 2016
!> \details
!>  This routine finalizes the subglacial hydro model
!-----------------------------------------------------------------------
   subroutine li_SGH_finalize(domain, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain  !< Input/Output: domain object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      logical, pointer :: config_SGH
      type (block_type), pointer :: block
      integer :: err_tmp

      err = 0
      err_tmp = 0

      call mpas_pool_get_config(liConfigs, 'config_SGH', config_SGH)
      if (.not. config_SGH) then
         return
      endif

      block => domain % blocklist
      do while (associated(block))

         block => block % next
      end do

      ! === error check
      if (err > 0) then
          call mpas_log_write("An error has occurred in li_SGH_finalize.", MPAS_LOG_ERR)
      endif

   !--------------------------------------------------------------------
   end subroutine li_SGH_finalize



   !--------------------------------------------------------------------
   !--------------------------------------------------------------------
   ! Local routines
   !--------------------------------------------------------------------
   !--------------------------------------------------------------------



!***********************************************************************
!
!  routine calc_edge_quantities
!
!> \brief   Calculate SGH fields on edges
!> \author  Matt Hoffman
!> \date    27 June 2016
!> \details
!>  This routine calculates needed SGH fields on edges
!-----------------------------------------------------------------------
   subroutine calc_edge_quantities(block, err)

      use mpas_geometry_utils, only: mpas_cells_to_points_using_baryweights
      use mpas_vector_operations, only: mpas_tangential_vector_1d
      use li_setup, only: li_cells_to_vertices_1dfield_using_kiteAreas

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: hydroPool
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBase
      real (kind=RKIND), dimension(:), pointer :: hydropotential
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBaseVertex
      real (kind=RKIND), dimension(:), pointer :: waterPressure
      real (kind=RKIND), dimension(:), pointer :: waterThicknessEdge
      real (kind=RKIND), dimension(:), pointer :: waterThicknessEdgeUpwind
      real (kind=RKIND), dimension(:), pointer :: waterThickness
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBaseSlopeNormal
      real (kind=RKIND), dimension(:), pointer :: hydropotentialSlopeNormal
      real (kind=RKIND), dimension(:), pointer :: waterPressureSlopeNormal
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBaseSlopeTangent
      real (kind=RKIND), dimension(:), pointer :: gradMagPhiEdge
      real (kind=RKIND), dimension(:), pointer :: effectiveConducEdge
      real (kind=RKIND), dimension(:), pointer :: diffusivity
      real (kind=RKIND), dimension(:), pointer :: dcEdge
      real (kind=RKIND), dimension(:), pointer :: dvEdge
      real (kind=RKIND), dimension(:), pointer :: waterVelocity
      real (kind=RKIND), dimension(:), pointer :: waterFlux
      real (kind=RKIND), dimension(:), pointer :: waterFluxAdvec
      real (kind=RKIND), dimension(:), pointer :: waterFluxDiffu
      integer, dimension(:), pointer :: waterFluxMask
      integer, dimension(:), pointer :: hydroMarineMarginMask
      integer, dimension(:,:), pointer :: edgeSignOnCell
      integer, dimension(:), pointer :: cellMask
      integer, dimension(:), pointer :: edgeMask
      integer, dimension(:,:), pointer :: cellsOnEdge
      integer, dimension(:,:), pointer :: verticesOnEdge
      integer, dimension(:,:), pointer :: baryCellsOnVertex
      real (kind=RKIND), dimension(:,:), pointer :: baryWeightsOnVertex
      real (kind=RKIND), pointer :: alpha, beta
      real (kind=RKIND), pointer :: conduc_coeff
      real (kind=RKIND), pointer :: conduc_coeff_drowned
      real (kind=RKIND), pointer :: bedRoughMax
      real (kind=RKIND) :: conduc_coeff_wtd
      character (len=StrKIND), pointer :: config_SGH_tangent_slope_calculation
      real (kind=RKIND), pointer :: config_sea_level
      real (kind=RKIND), pointer :: rhoo
      integer, pointer :: nEdges
      integer, pointer :: nCells
      integer, pointer :: nVertices
      integer :: iEdge, cell1, cell2
      real (kind=RKIND) :: velSign
      integer :: numGroundedCells
      integer :: err_tmp

      err = 0
      err_tmp = 0


      ! Get pools things
      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)

      call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
      call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
      call mpas_pool_get_dimension(meshPool, 'nVertices', nVertices)

      call mpas_pool_get_config(liConfigs, 'config_SGH_alpha', alpha)
      call mpas_pool_get_config(liConfigs, 'config_SGH_beta', beta)
      call mpas_pool_get_config(liConfigs, 'config_SGH_conduc_coeff', conduc_coeff)
      call mpas_pool_get_config(liConfigs, 'config_SGH_conduc_coeff_drowned', conduc_coeff_drowned)
      call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness_max', bedRoughMax)
      call mpas_pool_get_config(liConfigs, 'config_SGH_tangent_slope_calculation', config_SGH_tangent_slope_calculation)
      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
      call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)

      call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
      call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
      call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase)
      call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential)
      call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
      call mpas_pool_get_array(geometryPool, 'thickness', thickness)
      call mpas_pool_get_array(hydroPool, 'waterThicknessEdge', waterThicknessEdge)
      call mpas_pool_get_array(hydroPool, 'waterThicknessEdgeUpwind', waterThicknessEdgeUpwind)
      call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal)
      call mpas_pool_get_array(hydroPool, 'hydropotentialSlopeNormal', hydropotentialSlopeNormal)
      call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal)
      call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
      call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
      call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
      call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeTangent', hydropotentialBaseSlopeTangent)
      call mpas_pool_get_array(hydroPool, 'gradMagPhiEdge', gradMagPhiEdge)
      call mpas_pool_get_array(hydroPool, 'effectiveConducEdge', effectiveConducEdge)
      call mpas_pool_get_array(hydroPool, 'diffusivity', diffusivity)
      call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity)
      call mpas_pool_get_array(hydroPool, 'waterFlux', waterFlux)
      call mpas_pool_get_array(hydroPool, 'waterFluxAdvec', waterFluxAdvec)
      call mpas_pool_get_array(hydroPool, 'waterFluxDiffu', waterFluxDiffu)
      call mpas_pool_get_array(hydroPool, 'waterFluxMask', waterFluxMask)
      call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
      call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)


      do iEdge = 1, nEdges
         cell1 = cellsOnEdge(1, iEdge)
         cell2 = cellsOnEdge(2, iEdge)

         !waterThicknessEdge(iEdge) = 0.5_RKIND * ( waterThickness(cell1) + waterThickness(cell2) )
         ! This version ignores the thickness where there is no grounded ice (one-sided average at margin)
         numGroundedCells = li_mask_is_grounded_ice_int(cellMask(cell1)) + li_mask_is_grounded_ice_int(cellMask(cell2))
         if (numGroundedCells > 0) then
            ! Assuming here that waterThickness has been zeroed in non-grounded locations
            waterThicknessEdge(iEdge) = ( waterThickness(cell1) + waterThickness(cell2) ) / real(numGroundedCells)
         else
            waterThicknessEdge(iEdge) = 0.0_RKIND
         endif

         hydropotentialBaseSlopeNormal(iEdge) = (hydropotentialBase(cell2) - hydropotentialBase(cell1)) / dcEdge(iEdge)
         hydropotentialSlopeNormal(iEdge) = (hydropotential(cell2) - hydropotential(cell1)) / dcEdge(iEdge)
         waterPressureSlopeNormal(iEdge) = (waterPressure(cell2) - waterPressure(cell1)) / dcEdge(iEdge)
      end do

      ! At terrestrial margin, ignore the downslope bed topography gradient.  Including it can lead to unrealistically large
      ! hydropotential gradients and unstable channel growth.
      ! We also want to do this at marine margins because otherwise the offshore topography can create a barrier to flow,
      ! but that is unrealistic.
      ! So for all boundaries of the hydro system, the hydropotential at the margin should be determined by the geometry
      ! at the edge of the cell in a 1-sided sense
      do iEdge = 1, nEdges
         if ( (li_mask_is_margin(edgeMask(iEdge)) .and. li_mask_is_grounded_ice(edgeMask(iEdge))) .or. &
             (hydroMarineMarginMask(iEdge)==1)) then
            cell1 = cellsOnEdge(1, iEdge)
            cell2 = cellsOnEdge(2, iEdge)
            if (li_mask_is_grounded_ice(cellMask(cell1))) then ! cell2 is the icefree cell - replace phi there with cell1 Phig
               hydropotentialBaseSlopeNormal(iEdge) = (rho_water * gravity * bedTopography(cell1) + &
                  max(rhoo * gravity * (config_sea_level - bedTopography(cell1)), 0.0_RKIND) - &
                  hydropotentialBase(cell1)) / dcEdge(iEdge)
               hydropotentialSlopeNormal(iEdge) = (rho_water * gravity * bedTopography(cell1) + &
                  max(rhoo * gravity * (config_sea_level - bedTopography(cell1)), 0.0_RKIND) - &
                  hydropotential(cell1)) / dcEdge(iEdge)
            else ! cell1 is the icefree cell - replace phi there with cell2 Phig
               hydropotentialBaseSlopeNormal(iEdge) = (hydropotentialBase(cell2) - &
                  ( rho_water * gravity * bedTopography(cell2) + &
                    max(rhoo * gravity * (config_sea_level - bedTopography(cell2)), 0.0_RKIND) ) ) / dcEdge(iEdge)
               hydropotentialSlopeNormal(iEdge) = (hydropotential(cell2) - &
                  ( rho_water * gravity * bedTopography(cell2) + &
                    max(rhoo * gravity * (config_sea_level - bedTopography(cell2)), 0.0_RKIND) ) ) / dcEdge(iEdge)
            endif ! which cell is icefree
         endif ! if edge of grounded ice
      end do

      ! Disallow flow from ocean to glacier, or land to glacier,
      ! which can occur under some circumstances
      ! For ocean this is invalid because ocean water has a different density!
      ! For land this would only happen if there is a supply of water, which is not currently handled.
      ! Do this by simply zeroing the hydropotential gradient in those cases.
      ! (Do this step only after the other hydropotential special cases are treated above.)
      do iEdge = 1, nEdges
         ! Find edges along GL or margin to check for 'backwards' flow
         if ((hydroMarineMarginMask(iEdge)==1) .or. &
             li_mask_is_margin(edgeMask(iEdge)) ) then
            ! Now check if flow is backwards
            cell1 = cellsOnEdge(1, iEdge)
            cell2 = cellsOnEdge(2, iEdge)
            if (hydropotentialBaseSlopeNormal(iEdge) > 0.0_RKIND) then
               ! flow is from cell2 to cell1
               if (.not. li_mask_is_grounded_ice(cellMask(cell2))) then
                  hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND
                  hydropotentialSlopeNormal(iEdge) = 0.0_RKIND
               endif
            elseif (hydropotentialBaseSlopeNormal(iEdge) < 0.0_RKIND) then
               ! flow is from cell1 to cell2
               if (.not. li_mask_is_grounded_ice(cellMask(cell1))) then
                  hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND
                  hydropotentialSlopeNormal(iEdge) = 0.0_RKIND
               endif
            endif
         endif ! GL edge or grounded margin
      end do


      ! Calculate tangent slope of hydropotentialBase - three possible methods to consider

      ! Calculate hydropotentialBaseVertex if needed
      call mpas_pool_get_array(hydroPool, 'hydropotentialBaseVertex', hydropotentialBaseVertex)
         ! < this array could be protected by logic if desired
      select case (trim(config_SGH_tangent_slope_calculation))
      case ('from_vertex_barycentric')
         call mpas_pool_get_array(meshPool, 'baryCellsOnVertex', baryCellsOnVertex)
         call mpas_pool_get_array(meshPool, 'baryWeightsOnVertex', baryWeightsOnVertex)
         call mpas_cells_to_points_using_baryweights(meshPool, baryCellsOnVertex(:, 1:nVertices), &
              baryWeightsOnVertex(:, 1:nVertices), hydropotentialBase, hydropotentialBaseVertex(1:nVertices), err_tmp)
              err = ior(err, err_tmp)
      case ('from_vertex_barycentric_kiteareas')
         call li_cells_to_vertices_1dfield_using_kiteAreas(meshPool, hydropotentialBase, hydropotentialBaseVertex)
      end select

      ! Now perform tangent slope calculation based on method chosen
      select case (trim(config_SGH_tangent_slope_calculation))
      case ('from_vertex_barycentric', 'from_vertex_barycentric_kiteareas')
         call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
         call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
         call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
         do iEdge = 1, nEdges
            ! Only calculate slope for edges that have ice on at least one side.
            if ( li_mask_is_ice(edgeMask(iEdge)) ) then
               hydropotentialBaseSlopeTangent(iEdge) = ( hydropotentialBaseVertex(verticesOnEdge(1,iEdge)) -  &
                     hydropotentialBaseVertex(verticesOnEdge(2,iEdge)) ) / dvEdge(iEdge)
            else
               hydropotentialBaseSlopeTangent(iEdge) = 0.0_RKIND
            endif
         end do  ! edges
      case ('from_normal_slope')
         call mpas_tangential_vector_1d(hydropotentialBaseSlopeNormal, meshPool, &
                  includeHalo=.false., tangentialVector=hydropotentialBaseSlopeTangent)
      case default
         call mpas_log_write('Invalid value for config_SGH_tangent_slope_calculation.', MPAS_LOG_ERR)
         err = 1
      end select

      ! calculate magnitude of gradient of Phi
      gradMagPhiEdge = sqrt(hydropotentialBaseSlopeNormal**2 + hydropotentialBaseSlopeTangent**2)

      ! calculate effective conductivity on edges
      if (conduc_coeff_drowned > 0.0_RKIND) then
         ! Use a thickness weighted conductivity coeff. when water thickness exceeds bump height
         do iEdge = 1, nEdges
            conduc_coeff_wtd = (conduc_coeff * min(waterThicknessEdge(iEdge), bedRoughMax) + &
                             conduc_coeff_drowned * max(waterThicknessEdge(iEdge) - bedRoughMax, 0.0_RKIND)) / &
                             (waterThicknessEdge(iEdge) + 1.0e-16_RKIND)  ! Regularization only applies where value doesn't matter
            effectiveConducEdge(iEdge) = conduc_coeff_wtd * waterThicknessEdge(iEdge)**(alpha-1.0_RKIND) * &
               (gradMagPhiEdge(iEdge)+1.0e-30_RKIND)**(beta - 2.0_RKIND)   ! small value used for regularization
         end do
      else
         ! Just use a single conductivity coeff.
         effectiveConducEdge(:) = conduc_coeff * waterThicknessEdge(:)**(alpha-1.0_RKIND) *&
            (gradMagPhiEdge(:)+1.0e-30_RKIND)**(beta - 2.0_RKIND)   ! small value used for regularization
      endif

      ! calculate diffusivity on edges
      diffusivity(:) = rho_water * gravity * effectiveConducEdge(:) * waterThicknessEdge(:)

      do iEdge = 1, nEdges
         cell1 = cellsOnEdge(1, iEdge)
         cell2 = cellsOnEdge(2, iEdge)
         waterVelocity(iEdge) = -1.0_RKIND * effectiveConducEdge(iEdge) * hydropotentialBaseSlopeNormal(iEdge)
         velSign = sign(1.0_RKIND, waterVelocity(iEdge))
         waterThicknessEdgeUpwind(iEdge) = max(velSign * waterThickness(cell1),   &
                     velSign * (-1.0_RKIND) * waterThickness(cell2))

         ! advective flux
         waterFluxAdvec(iEdge) = waterVelocity(iEdge) * waterThicknessEdgeUpwind(iEdge)

         ! diffusive flux
         ! Note: At the GL, the water thickness for one cell will be 0, meaning a large gradient.  However, the
         ! diffusivity uses a one sided water thickness.  It's unclear what really happens to lakes at the GL/margin.
         waterFluxDiffu(iEdge) = -1.0_RKIND * diffusivity(iEdge) * (waterThickness(cell2) - waterThickness(cell1)) &
               / dcEdge(iEdge)
         !endif
      end do
      where (waterFluxMask == 2)
         waterFluxAdvec = 0.0_RKIND
         waterFluxDiffu = 0.0_RKIND
         waterVelocity = 0.0_RKIND
      end where
      where (waterThicknessEdgeUpwind == 0.0_RKIND)  ! if no water to give, should have no flux!
         waterFluxAdvec = 0.0_RKIND ! Should already be 0
         waterFluxDiffu = 0.0_RKIND ! Might not be 0
         waterVelocity = 0.0_RKIND
      end where

      waterFlux(:) = waterFluxAdvec(:) + waterFluxDiffu(:)

   !--------------------------------------------------------------------
   end subroutine calc_edge_quantities




!***********************************************************************
!
!  routine check_timestep
!
!> \brief   Calculate SGH timesteps and check current timestep
!> \author  Matt Hoffman
!> \date    7 July 2016
!> \details
!>  This routine calculates the three timesteps associated with the
!>  SGH solver and compares them to the current model timestep.
!-----------------------------------------------------------------------
   subroutine check_timestep(domain, timeLeft, numSubCycles, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------
      integer, intent(in) :: numSubCycles  !< Input: number of subcycles taken so far

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain    !< Input/Output: domain object
      real (kind=RKIND), intent(inout) :: timeLeft  !< Input/Output: time remaining for subcycling (seconds)

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: hydroPool
      real (kind=RKIND), dimension(:), pointer :: waterVelocity
      real (kind=RKIND), dimension(:), pointer :: channelVelocity
      real (kind=RKIND), dimension(:), pointer :: diffusivity
      real (kind=RKIND), dimension(:), pointer :: channelDiffusivity
      real (kind=RKIND), dimension(:), pointer :: dcEdge
      integer, dimension(:), pointer :: indexToEdgeID
      real (kind=RKIND), pointer :: deltatSGH
      real (kind=RKIND), pointer :: deltatSGHadvec
      real (kind=RKIND), pointer :: deltatSGHdiffu
      real (kind=RKIND), pointer :: deltatSGHpressure
      real (kind=RKIND), pointer :: deltatSGHadvecChannel
      real (kind=RKIND), pointer :: deltatSGHdiffuChannel
      real (kind=RKIND), pointer :: porosity
      type (block_type), pointer :: block
      real (kind=RKIND), pointer :: deltat
      integer, pointer :: nEdgesSolve
      logical, pointer :: config_SGH_chnl_active, config_SGH_chnl_include_DCFL
      ! in the following variables, "Proc" indicates the value on the current processor,
      ! and "Block" indicates value on current block
      real (kind=RKIND) :: dtSGHadvecBlock, dtSGHadvecProc
      real (kind=RKIND) :: dtSGHdiffuBlock, dtSGHdiffuProc
      real (kind=RKIND) :: dtSGHpressureBlock, dtSGHpressureProc
      real (kind=RKIND) :: dtSGHadvecChanBlock, dtSGHadvecChanProc
      real (kind=RKIND) :: dtSGHdiffuChanBlock, dtSGHdiffuChanProc
      integer :: err_tmp
      real(kind=RKIND), parameter :: bigNumber = 1.0e36_RKIND
      real(kind=RKIND), pointer :: CFLfraction
      real(kind=RKIND), pointer :: maxDt
      real(kind=RKIND) :: proposedDt
      real(kind=RKIND) :: masterDt
      real(kind=RKIND), dimension(5) :: localMinValues, reducedMinValues
      integer, dimension(5) :: localLocValues, reducedLocValues
      integer :: iEdge
      real(kind=RKIND), save :: minDtSoFar
      integer, save :: numDtLess1
      logical :: printDtInfo


      err = 0
      err_tmp = 0

      printDtInfo = .false.
      if (numSubCycles == 1) then
         minDtSoFar = 1.0e36_RKIND
         numDtLess1 = 0
      endif

      call mpas_pool_get_config(liConfigs, 'config_SGH_englacial_porosity', porosity)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_active', config_SGH_chnl_active)

      dtSGHadvecProc = bigNumber
      dtSGHdiffuProc = bigNumber
      dtSGHpressureProc = bigNumber
      dtSGHadvecChanProc = bigNumber
      dtSGHdiffuChanProc = bigNumber

      ! ---
      ! Find local (block) limiting dt for various CFL conditions
      ! ---
      block => domain % blocklist
      do while (associated(block))

         ! Get pools things
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)

         call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)

         call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
         call mpas_pool_get_array(hydroPool, 'deltatSGHadvec', deltatSGHadvec)
         call mpas_pool_get_array(hydroPool, 'deltatSGHdiffu', deltatSGHdiffu)
         call mpas_pool_get_array(hydroPool, 'deltatSGHpressure', deltatSGHpressure)
         call mpas_pool_get_array(hydroPool, 'deltatSGHadvecChannel', deltatSGHadvecChannel)
         call mpas_pool_get_array(hydroPool, 'deltatSGHdiffuChannel', deltatSGHdiffuChannel)
         call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity)
         call mpas_pool_get_array(hydroPool, 'channelVelocity', channelVelocity)
         call mpas_pool_get_array(hydroPool, 'diffusivity', diffusivity)
         call mpas_pool_get_array(hydroPool, 'channelDiffusivity', channelDiffusivity)
         call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
         call mpas_pool_get_array(meshPool, 'indexToEdgeID', indexToEdgeID)

         ! Calculate advective CFL-limited time step
         dtSGHadvecBlock = 0.5_RKIND * minval(dcEdge(1:nEdgesSolve) / (abs(waterVelocity(1:nEdgesSolve)) &
            + 1.0e-12_RKIND))  ! regularize
         dtSGHadvecProc = min(dtSGHadvecProc, dtSGHadvecBlock)

         ! Calculate diffusive CFL-limited time step
         dtSGHdiffuBlock = 0.25_RKIND * minval(dcEdge(1:nEdgesSolve)**2 / (diffusivity(1:nEdgesSolve) + 1.0e-12_RKIND))
         dtSGHdiffuProc = min(dtSGHdiffuProc, dtSGHdiffuBlock)

         ! Calculate pressure limited time step
         dtSGHpressureBlock = 1.0_RKIND * minval(porosity * dcEdge(1:nEdgesSolve)**2 &
                   / (2.0_RKIND * diffusivity(1:nEdgesSolve) + 1.0e-12_RKIND))
         dtSGHpressureProc = min(dtSGHpressureProc, dtSGHpressureBlock)

         if (config_SGH_chnl_active) then
            ! Calculate channel advection limited time step
            dtSGHadvecChanBlock = 0.5_RKIND * minval(dcEdge(1:nEdgesSolve) / (abs(channelVelocity(1:nEdgesSolve)) &
               + 1.0e-12_RKIND))
            ! regularize
            dtSGHadvecChanProc = min(dtSGHadvecChanProc, dtSGHadvecChanBlock)
            ! Calculate channel diffusion limited time step
            dtSGHdiffuChanBlock = 0.25_RKIND * minval(dcEdge(1:nEdgesSolve)**2 / (channelDiffusivity(1:nEdgesSolve) + &
               1.0e-12_RKIND))
            dtSGHdiffuChanProc = min(dtSGHdiffuChanProc, dtSGHdiffuChanBlock)
         endif

         ! Master deltat is needed below, so grab it in this block loop
         call mpas_pool_get_array(meshPool, 'deltat', deltat)

         block => block % next
      end do


      ! ---
      ! reduce across procs
      ! ---
      localMinValues(1) = dtSGHadvecProc
      localMinValues(2) = dtSGHdiffuProc
      localMinValues(3) = dtSGHpressureProc
      if (config_SGH_chnl_active) then
         localMinValues(4) = dtSGHadvecChanProc
         localMinValues(5) = dtSGHdiffuChanProc
         call mpas_timer_start("global reduce")
         call mpas_dmpar_min_real_array(domain % dminfo, 5, localMinValues, reducedMinValues)
         call mpas_timer_stop("global reduce")
         deltatSGHadvecChannel = reducedMinValues(4)
         deltatSGHdiffuChannel = reducedMinValues(5)
      else
         call mpas_timer_start("global reduce")
         call mpas_dmpar_min_real_array(domain % dminfo, 3, localMinValues(1:3), reducedMinValues(1:3))
         call mpas_timer_stop("global reduce")
      endif
      deltatSGHadvec = reducedMinValues(1)
      deltatSGHdiffu = reducedMinValues(2)
      deltatSGHpressure = reducedMinValues(3)


      ! ---
      ! Find and set the new hydro subcycling dt
      ! ---
      call mpas_pool_get_config(liConfigs, 'config_SGH_adaptive_timestep_fraction', CFLfraction)
      call mpas_pool_get_config(liConfigs, 'config_SGH_max_adaptive_timestep', maxDt)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_include_DCFL', config_SGH_chnl_include_DCFL)

      ! Find smallest of 3 or 5 limiting time steps
      proposedDt = min(deltatSGHadvec, deltatSGHdiffu, deltatSGHpressure)
      if (config_SGH_chnl_active) then
         proposedDt = min(proposedDt, deltatSGHadvecChannel)
         if (config_SGH_chnl_include_DCFL) then
            proposedDt = min(proposedDt, deltatSGHdiffuChannel)
         endif
      endif
      proposedDt = proposedDt * CFLfraction

      if (proposedDt < 1.0_RKIND) then
         numDtLess1 = numDtLess1 + 1
      endif
      minDtSoFar = min(minDtSoFar, proposedDt)

      ! Write out dt stats every now and then
      if (mod(numSubCycles, 2000)==0) then
         call mpas_log_write("SGH subcycle #$i. Time left=$r s.", intArgs=(/numSubCycles/), realArgs=(/timeLeft/))
         printDtInfo=.true.
      endif

      ! Write out dt info on the final subcycle
      if (proposedDt >= timeLeft) then
         call mpas_log_write("SGH completed subcycling with $i timesteps.", intArgs=(/numSubCycles/))
         printDtInfo=.true.
      endif

      if (proposedDt < 1.0E-2_RKIND) then
         call mpas_log_write("CFL conditions are limiting subglacial hydrology timestep to <1e-2s ($r) on subcycle $i.", &
                 MPAS_LOG_WARN, realArgs=(/proposedDt/), intArgs=(/numSubCycles/))
         printDtInfo=.true.
      endif

      if (proposedDt < 1.0E-4_RKIND) then
         call mpas_log_write("CFL conditions are limiting subglacial hydrology timestep to <1e-4s.", MPAS_LOG_ERR)
         printDtInfo=.true.
         err = ior(err, 1)
      endif

      if (printDtInfo) then
         ! Determine location of limiting Dt
         ! NOTE: ignoring multiple blocks here.
         if (localMinValues(1) == reducedMinValues(1)) then ! limiting edge is on this proc
            localLocValues(1) = indexToEdgeID(minloc(dcEdge(1:nEdgesSolve) / (abs(waterVelocity(1:nEdgesSolve)) &
               + 1.0e-12_RKIND), dim=1))
         else
            localLocValues(1) = -1
         endif

         if (localMinValues(2) == reducedMinValues(2)) then ! limiting edge is on this proc
            localLocValues(2) = indexToEdgeID(minloc(dcEdge(1:nEdgesSolve)**2 / (diffusivity(1:nEdgesSolve) + 1.0e-12_RKIND), &
                    dim=1))
         else
            localLocValues(2) = -1
         endif

         if (localMinValues(3) == reducedMinValues(3)) then ! limiting edge is on this proc
            localLocValues(3) = indexToEdgeID(minloc(porosity * dcEdge(1:nEdgesSolve)**2 &
                      / (2.0_RKIND * diffusivity(1:nEdgesSolve) + 1.0e-12_RKIND), dim=1))
         else
            localLocValues(3) = -1
         endif

         if (config_SGH_chnl_active) then
            if (localMinValues(4) == reducedMinValues(4)) then ! limiting edge is on this proc
               localLocValues(4) = indexToEdgeID(minloc(dcEdge(1:nEdgesSolve) / (abs(channelVelocity(1:nEdgesSolve)) &
                  + 1.0e-12_RKIND), dim=1))
            else
               localLocValues(4) = -1
            endif

            if (localMinValues(5) == reducedMinValues(5)) then ! limiting edge is on this proc
               localLocValues(5) = indexToEdgeID(minloc(dcEdge(1:nEdgesSolve)**2 / (channelDiffusivity(1:nEdgesSolve) + &
                  1.0e-12_RKIND), dim=1))
            else
               localLocValues(5) = -1
            endif

            ! global reduce of locations
            call mpas_timer_start("global reduce")
            call mpas_dmpar_max_int_array(domain % dminfo, 5, localLocValues, reducedLocValues)
            call mpas_timer_stop("global reduce")
         else
            call mpas_timer_start("global reduce")
            call mpas_dmpar_max_int_array(domain % dminfo, 3, localLocValues, reducedLocValues)
            call mpas_timer_stop("global reduce")
         endif


         if (config_SGH_chnl_active) then
            call mpas_log_write("deltatSGH: used=$r, prev=$r, advec=$r, diffu=$r, pressure=$r, advecChannel=$r, diffuChannel=$r",&
              realArgs=(/proposedDt, deltatSGH, deltatSGHadvec, deltatSGHdiffu, deltatSGHpressure, deltatSGHadvecChannel, &
                         deltatSGHdiffuChannel/))
            call mpas_log_write("Limiting edge IDs: advec=$i, diffu=$i, pressure=$i, advecChannel=$i, diffuChannel=$i", &
               intArgs=reducedLocValues(1:5))
         else
            call mpas_log_write("deltatSGH: used=$r, previous=$r, advec=$r, diffu=$r, pressure=$r", &
              realArgs=(/proposedDt, deltatSGH, deltatSGHadvec, deltatSGHdiffu, deltatSGHpressure/))
            call mpas_log_write("Limiting edge IDs: advec=$i, diffu=$i, pressure=$i", &
               intArgs=reducedLocValues(1:3))
         endif
         call mpas_log_write("  min hydro dt used=$r; There were $i subcyles with dt<1s.", &
                 realArgs=(/minDtSoFar/), intArgs=(/numDtLess1/))
      endif

      ! Don't exceed the maximum allowed hydro time step
      proposedDt = min(proposedDt, maxDt)
      ! Don't let SGH time step exceed time left in the master model dt
      proposedDt = min(proposedDt, timeLeft)

      !call mpas_log_write('  Setting SGH time step to: $r seconds', realArgs=(/proposedDt/))

      timeLeft = timeLeft - proposedDt
      !call mpas_log_write("dt=$r, new TIMELEFT=$r", realArgs=(/proposedDt, timeLeft/))


      ! ---
      ! Assign new time step to all blocks
      ! ---
      block => domain % blocklist
      do while (associated(block))

         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
         deltatSGH = proposedDt

         block => block % next
      end do


      ! compare timesteps
      if (deltatSGH > deltatSGHadvec) then
         call mpas_log_write("deltatSGH > deltatSGHadvec  $r > $r", MPAS_LOG_WARN, realArgs=(/deltatSGH, deltatSGHadvec/))
      endif
      if (deltatSGH > deltatSGHdiffu) then
         call mpas_log_write("deltatSGH > deltatSGHdiffu  $r > $r", MPAS_LOG_WARN, realArgs=(/deltatSGH, deltatSGHdiffu/))
      endif
      if (deltatSGH > deltatSGHpressure) then
         call mpas_log_write("deltatSGH > deltatSGHpressure  $r > $r", MPAS_LOG_WARN, realArgs=(/deltatSGH, deltatSGHpressure/))
      endif
      if (config_SGH_chnl_active .and. (deltatSGH > deltatSGHadvecChannel)) then
         call mpas_log_write("deltatSGH > deltatSGHadvecChannel  $r > $r", MPAS_LOG_WARN, &
            realArgs=(/deltatSGH, deltatSGHadvecChannel/))
      endif
      if (config_SGH_chnl_active .and. config_SGH_chnl_include_DCFL .and. (deltatSGH > deltatSGHdiffuChannel)) then
         call mpas_log_write("deltatSGH > deltatSGHdiffuChannel  $r > $r", MPAS_LOG_WARN, &
            realArgs=(/deltatSGH, deltatSGHdiffuChannel/))
      endif


   !--------------------------------------------------------------------
   end subroutine check_timestep




!***********************************************************************
!
!  routine calc_pressure
!
!> \brief   Calculate SGH water pressure
!> \author  Matt Hoffman
!> \date    5 July 2016
!> \details
!>  This routine calculates SGH water pressure
!-----------------------------------------------------------------------
   subroutine calc_pressure(block, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: velocityPool
      real (kind=RKIND), dimension(:), pointer :: waterPressure
      real (kind=RKIND), dimension(:), pointer :: waterPressureOld
      real (kind=RKIND), dimension(:), pointer :: waterPressureTendency
      real (kind=RKIND), dimension(:), pointer :: waterThickness
      real (kind=RKIND), dimension(:), pointer :: effectivePressure
      real (kind=RKIND), dimension(:), pointer :: zeroOrderSum
      real (kind=RKIND), dimension(:), pointer :: closingRate
      real (kind=RKIND), dimension(:), pointer :: openingRate
      real (kind=RKIND), dimension(:), pointer :: basalMeltInput
      real (kind=RKIND), dimension(:), pointer :: externalWaterInput
      real (kind=RKIND), dimension(:), pointer :: Wtill, WtillOld
      real (kind=RKIND), dimension(:), pointer :: divergence
      real (kind=RKIND), dimension(:), pointer :: basalSpeed
      real (kind=RKIND), dimension(:,:), pointer :: flowParamA
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:), pointer :: divergenceChannel
      real (kind=RKIND), dimension(:), pointer :: channelAreaChangeCell
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      integer, dimension(:), pointer :: hydroMarineMarginMask
      integer, dimension(:), pointer :: cellMask
      integer, dimension(:), pointer :: nEdgesOnCell
      integer, dimension(:,:), pointer :: edgesOnCell
      real (kind=RKIND), pointer :: deltatSGH
      real (kind=RKIND), pointer :: bedRough, bedRoughMax
      real (kind=RKIND), pointer :: rhoi
      real (kind=RKIND), pointer :: creepCoeff
      real (kind=RKIND), pointer :: porosity
      integer, pointer :: nVertLevels
      integer, pointer :: nCells
      character (len=StrKIND), pointer :: config_SGH_pressure_calc
      real (kind=RKIND), pointer :: config_sea_level
      real (kind=RKIND), pointer :: rhoo
      integer :: err_tmp
      integer :: iCell
      integer :: iEdge
      logical :: onMarineMargin

      err = 0
      err_tmp = 0


      ! Get pools things
      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
      call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)

      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
      call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)

      call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
      call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness', bedRough)
      call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness_max', bedRoughMax)
      call mpas_pool_get_config(liConfigs, 'config_SGH_creep_coefficient', creepCoeff)
      call mpas_pool_get_config(liConfigs, 'config_SGH_englacial_porosity', porosity)
      call mpas_pool_get_config(liConfigs, 'config_SGH_pressure_calc', config_SGH_pressure_calc)

      call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
      call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
      call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
      call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)

      call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
      call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
      call mpas_pool_get_array(hydroPool, 'waterPressureOld', waterPressureOld)
      call mpas_pool_get_array(hydroPool, 'waterPressureTendency', waterPressureTendency)
      call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
      call mpas_pool_get_array(hydroPool, 'zeroOrderSum', zeroOrderSum)
      call mpas_pool_get_array(hydroPool, 'closingRate', closingRate)
      call mpas_pool_get_array(hydroPool, 'openingRate', openingRate)
      call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
      call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
      call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
      call mpas_pool_get_array(hydroPool, 'tillWaterThickness', Wtill)
      call mpas_pool_get_array(hydroPool, 'tillWaterThicknessOld', WtillOld)
      call mpas_pool_get_array(hydroPool, 'divergence', divergence)
      call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
      call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
      call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
      call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed)
      call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
      call mpas_pool_get_array(geometryPool, 'thickness', thickness)
      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
      call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)

      openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:))
      !openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:)) + &
      !   basalMeltInput(:) / rhoi  ! Hewitt 2011 opening
      openingRate = max(0.0_RKIND, openingRate)

      closingRate(:) = creepCoeff * flowParamA(nVertLevels, :) * &
              effectivePressure(:)**3 * waterThickness(:)
!      closingRate(:) = waterThickness(:) * effectivePressure(:) / 1.0e13_RKIND
!          ! Hewitt 2011 creep closure form.  Denominator is ice viscosity

      zeroOrderSum = closingRate - openingRate + (basalMeltInput + externalWaterInput) / rho_water - &
         (Wtill - WtillOld) / deltatSGH

      waterPressureOld = waterPressure

      select case (trim(config_SGH_pressure_calc))
      case ('cavity')

         where (li_mask_is_floating_ice(cellMask))
            waterPressure = rhoi * gravity * thickness
         elsewhere (.not. li_mask_is_ice(cellMask))
            waterPressure = 0.0_RKIND
         elsewhere
            waterPressure = (zeroOrderSum - divergence - divergenceChannel - channelAreaChangeCell) * &
               rho_water * gravity * deltatSGH / porosity + waterPressureOld
         end where

      case ('overburden')
         where (li_mask_is_floating_ice(cellMask))
            waterPressure = rhoi * gravity * thickness
         elsewhere (.not. li_mask_is_ice(cellMask))
            waterPressure = 0.0_RKIND
         elsewhere
            waterPressure = rhoi * gravity * thickness
         end where

      case default
         call mpas_log_write("Invalid option specified for config_SGH_pressure_calc:" // config_SGH_pressure_calc, MPAS_LOG_ERR)
         err = ior(err, 1)
      end select

      waterPressure = max(0.0_RKIND, waterPressure)
      waterPressure = min(waterPressure, rhoi * gravity * thickness)

      do iCell = 1, nCells
        if ( li_mask_is_floating_ice(cellMask(iCell)) .or. &
             ((.not. li_mask_is_ice(cellMask(iCell))) .and. (bedTopography(iCell) < config_sea_level) ) ) then
           ! set pressure correctly under floating ice and open ocean
           waterPressure(iCell) = rhoo * gravity * (config_sea_level - bedTopography(iCell))
        else
           onMarineMargin = .false.
           do iEdge = 1, nEdgesOnCell(iCell)
              if (hydroMarineMarginMask(edgesOnCell(iEdge, iCell)) == 1) then
                 onMarineMargin = .true.
                 exit
              endif
           enddo
           if (onMarineMargin) then
              ! At marine margin, don't let water pressure fall below ocean pressure
              ! TODO: Not sure if this should include the water layer thickness term.  Leaving it off.
              if (waterPressure(iCell) < rhoo * gravity * (config_sea_level - bedTopography(iCell))) then
                  waterPressure(iCell) = rhoo * gravity * (config_sea_level - bedTopography(iCell))
              endif
           endif
        endif
      enddo

      waterPressureTendency = (waterPressure - waterPressureOld) / deltatSGH

      call calc_pressure_diag_vars(block, err_tmp)
      err = ior(err, err_tmp)

   !--------------------------------------------------------------------
   end subroutine calc_pressure


!***********************************************************************
!
!  routine calc_pressure_diag_vars
!
!> \brief   Calculate SGH diagnostic variables related to  pressure
!> \author  Matt Hoffman
!> \date    5 July 2016
!> \details
!>  This routine calculates variables related to water pressure
!-----------------------------------------------------------------------
   subroutine calc_pressure_diag_vars(block, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: hydroPool
      real (kind=RKIND), pointer :: rhoi, rhoo
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:), pointer :: waterPressure
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBase
      real (kind=RKIND), dimension(:), pointer :: waterThickness
      real (kind=RKIND), dimension(:), pointer :: hydropotential
      real (kind=RKIND), dimension(:), pointer :: effectivePressure
      integer, dimension(:), pointer :: cellMask
      real (kind=RKIND), pointer :: config_sea_level

      err = 0

      ! Get pools things
      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)

      call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
      call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)

      call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
      call mpas_pool_get_array(geometryPool, 'thickness', thickness)
      call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
      call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
      call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase)
      call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
      call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential)
      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

      effectivePressure = rhoi * gravity * thickness - waterPressure
         ! < this should evalute to 0 for floating ice if Pw set correctly there.
      where (.not. li_mask_is_grounded_ice(cellmask))
         effectivePressure = 0.0_RKIND  ! zero effective pressure where no ice to avoid confusion
      end where

      hydropotentialBase = rho_water * gravity * bedTopography + waterPressure
      ! This is still correct under ice shelves/open ocean because waterPressure has been set appropriately there already.
      ! Note this leads to a nonuniform hydropotential at sea level that is a function of the ocean depth.
      ! That is what we want because we use this as a boundary condition on the subglacial system,
      ! and we want the subglacial system to feel the pressure of the ocean column at its edge.

      ! hydropotential with water thickness
      hydropotential = hydropotentialBase + rho_water * gravity * waterThickness

   !--------------------------------------------------------------------
   end subroutine calc_pressure_diag_vars



!***********************************************************************
!
!  routine update_channel
!
!> \brief   Calculate SGH channel area
!> \author  Matt Hoffman
!> \date    28 July 2016
!> \details
!>  This routine updates the channel area in the subglacial hydrology model.
!>  It uses the conduit space evolution equation.
!-----------------------------------------------------------------------
   subroutine update_channel(block, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
!!      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: meshPool
      type (mpas_pool_type), pointer :: velocityPool
      type (mpas_pool_type), pointer :: geometryPool
      real (kind=RKIND), pointer :: Kc
      real (kind=RKIND), pointer :: alpha_c
      real (kind=RKIND), pointer :: beta_c
      real (kind=RKIND), pointer :: creep_coeff
      real (kind=RKIND), pointer :: rhoi
      real (kind=RKIND), pointer :: config_SGH_incipient_channel_width
      logical, pointer :: config_SGH_include_pressure_melt

      real (kind=RKIND), dimension(:), pointer :: channelArea
      real (kind=RKIND), dimension(:), pointer :: channelMelt
      real (kind=RKIND), dimension(:), pointer :: channelPressureFreeze
      real (kind=RKIND), dimension(:), pointer :: channelDischarge
      real (kind=RKIND), dimension(:), pointer :: channelVelocity
      real (kind=RKIND), dimension(:), pointer :: gradMagPhiEdge
      real (kind=RKIND), dimension(:), pointer :: waterFlux
      real (kind=RKIND), dimension(:), pointer :: hydropotentialBaseSlopeNormal
      real (kind=RKIND), dimension(:), pointer :: waterPressureSlopeNormal
      real (kind=RKIND), dimension(:), pointer :: channelOpeningRate
      real (kind=RKIND), dimension(:), pointer :: channelClosingRate
      real (kind=RKIND), dimension(:), pointer :: channelChangeRate
      real (kind=RKIND), dimension(:), pointer :: flowParamAChannel
      real (kind=RKIND), dimension(:), pointer :: channelEffectivePressure
      real (kind=RKIND), dimension(:), pointer :: effectivePressure
      real (kind=RKIND), dimension(:), pointer :: channelDiffusivity
      integer, dimension(:), pointer :: waterFluxMask
      integer, dimension(:), pointer :: hydroMarineMarginMask
      integer, dimension(:), pointer :: edgeMask
      real (kind=RKIND), dimension(:,:), pointer :: flowParamA
      integer, dimension(:,:), pointer :: cellsOnEdge
      integer, pointer :: nVertLevels

      integer, pointer :: nEdgesSolve
      integer :: iEdge, cell1, cell2


      err = 0

      ! Get pools things
      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
!      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
      call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)

      call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_conduc_coeff', Kc)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_alpha', alpha_c)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_beta', beta_c)
      call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_creep_coefficient', creep_coeff)
      call mpas_pool_get_config(liConfigs, 'config_SGH_incipient_channel_width', config_SGH_incipient_channel_width)
      call mpas_pool_get_config(liConfigs, 'config_SGH_include_pressure_melt', config_SGH_include_pressure_melt)

      call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
      call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)

      call mpas_pool_get_array(hydroPool, 'channelArea', channelArea)
      call mpas_pool_get_array(hydroPool, 'channelMelt', channelMelt)
      call mpas_pool_get_array(hydroPool, 'channelPressureFreeze', channelPressureFreeze)
      call mpas_pool_get_array(hydroPool, 'channelDischarge', channelDischarge)
      call mpas_pool_get_array(hydroPool, 'channelVelocity', channelVelocity)
      call mpas_pool_get_array(hydroPool, 'gradMagPhiEdge', gradMagPhiEdge)
      call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal)
      call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal)
      call mpas_pool_get_array(hydroPool, 'waterFlux', waterFlux)
      call mpas_pool_get_array(hydroPool, 'channelOpeningRate', channelOpeningRate)
      call mpas_pool_get_array(hydroPool, 'channelClosingRate', channelClosingRate)
      call mpas_pool_get_array(hydroPool, 'channelChangeRate', channelChangeRate)
      call mpas_pool_get_array(hydroPool, 'flowParamAChannel', flowParamAChannel)
      call mpas_pool_get_array(hydroPool, 'channelEffectivePressure', channelEffectivePressure)
      call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
      call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
      call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
      call mpas_pool_get_array(hydroPool, 'waterFluxMask', waterFluxMask)
      call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
      call mpas_pool_get_array(hydroPool, 'channelDiffusivity', channelDiffusivity)
      call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)

      ! Calculate terms needed for opening (melt) rate

      where(gradMagPhiEdge < 0.01_RKIND)
         channelDischarge(:) = 0.0_RKIND
      elsewhere
         channelDischarge = -1.0_RKIND * Kc * channelArea**alpha_c * gradMagPhiEdge**(beta_c - 2.0_RKIND) * &
            hydropotentialBaseSlopeNormal
      end where

      where (waterFluxMask == 2)
         channelDischarge = 0.0_RKIND
         channelArea = 0.0_RKIND
      end where

      ! Note: an edge with only one grounded cell neighbor is called floating, so this logic retains channel vars
      ! on those edges to allow channel discharge across GL
      where (.not. ( (li_mask_is_grounded_ice(edgeMask)) .or. (hydroMarineMarginMask==1) ) )
         channelArea = 0.0_RKIND
         channelDischarge = 0.0_RKIND
      end where

      ! Disable channels from forming if there is no sheet flux
      ! TODO: Make a function of sheet dissipation threshold?
      where (abs(waterFlux) <= 1e-10_RKIND)
         channelArea = 0.0_RKIND
         channelDischarge = 0.0_RKIND
      end where

      channelVelocity = channelDischarge / (channelArea + 1.0e-12_RKIND)

      ! diffusivity used only to limit channel dt right now
      where(gradMagPhiEdge < 0.01_RKIND)
         channelDiffusivity = 0.0_RKIND
      elsewhere
         channelDiffusivity = abs(rho_water * gravity * channelArea *  &
            Kc * channelArea**(alpha_c - 1.0_RKIND) * gradMagPhiEdge**(beta_c - 2.0_RKIND))
      end where

      channelMelt = (abs(channelDischarge * hydropotentialBaseSlopeNormal) &  ! channel dissipation
                  +  abs(waterFlux * hydropotentialBaseSlopeNormal * config_SGH_incipient_channel_width) & !some sheet dissipation
                  ) / latent_heat_ice
      channelPressureFreeze = -1.0_RKIND * iceMeltingPointPressureDependence * cp_freshwater * rho_water * &
         (channelDischarge + waterFlux * config_SGH_incipient_channel_width) &
         * waterPressureSlopeNormal / latent_heat_ice

      if (config_SGH_include_pressure_melt) then
         channelOpeningRate = (channelMelt - channelPressureFreeze) / rhoi
      else
         channelOpeningRate = channelMelt / rhoi
      endif

      ! Calculate terms needed for closing (creep) rate
      ! Need cell center quantities on edges
      do iEdge = 1, nEdgesSolve
         cell1 = cellsOnEdge(1, iEdge)
         cell2 = cellsOnEdge(2, iEdge)

         ! Not sure if these ought to be upwind average, but using centered
         flowParamAChannel(iEdge) = 0.5_RKIND * (flowParamA(nVertLevels, cell1) + flowParamA(nVertLevels, cell2) )
         channelEffectivePressure(iEdge) = 0.5_RKIND * (effectivePressure(cell1) + effectivePressure(cell2))
      end do
      channelClosingRate(:) = creep_coeff * channelArea(:) * flowParamAChannel(:) * channelEffectivePressure(:)**3

      where (waterFluxMask == 2)
         channelOpeningRate = 0.0_RKIND
         channelClosingRate = 0.0_RKIND
      end where
      channelChangeRate = channelOpeningRate - channelClosingRate

   !--------------------------------------------------------------------
   end subroutine update_channel


!***********************************************************************
!
!  routine evolve_channel
!
!> \brief   Evolve SGH channel area, calculate changes on cells
!> \author  Matt Hoffman
!> \date    28 July 2016
!> \details
!>  This routine updates the channel area in the subglacial hydrology model.
!>  It uses the conduit space evolution equation.
!-----------------------------------------------------------------------
   subroutine evolve_channel(block, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: meshPool
      real (kind=RKIND), dimension(:), pointer :: channelArea
      real (kind=RKIND), dimension(:), pointer :: channelDischarge
      real (kind=RKIND), dimension(:), pointer :: channelChangeRate
      real (kind=RKIND), dimension(:), pointer :: divergenceChannel
      real (kind=RKIND), dimension(:), pointer :: channelAreaChangeCell
      real (kind=RKIND), pointer :: deltatSGH
      integer, dimension(:), pointer :: nEdgesOnCell
      integer, dimension(:,:), pointer :: edgesOnCell
      integer, dimension(:,:), pointer :: cellsOnEdge
      integer, dimension(:,:), pointer :: edgeSignOnCell
      real (kind=RKIND), dimension(:), pointer :: areaCell
      real (kind=RKIND), dimension(:), pointer :: dcEdge
      integer, pointer :: nCellsSolve
      integer :: iCell, iEdgeOnCell, iEdge

      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
      call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
      call mpas_pool_get_array(hydroPool, 'channelArea', channelArea)
      call mpas_pool_get_array(hydroPool, 'channelDischarge', channelDischarge)
      call mpas_pool_get_array(hydroPool, 'channelChangeRate', channelChangeRate)
      call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
      call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
      call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
      call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
      call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
      call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
      call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
      call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
      call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)

      ! Calculate flux divergence on cells and channel area change on cells
      divergenceChannel(:) = 0.0_RKIND  ! zero div before accumulating
      channelAreaChangeCell(:) = 0.0_RKIND  ! zero before accumulating
      ! loop over locally owned cells
      do iCell = 1, nCellsSolve
         ! TODO: could limit to grounded cells only
         ! compute fluxes for each edge of the cell
         do iEdgeOnCell = 1, nEdgesOnCell(iCell)
            iEdge = edgesOnCell(iEdgeOnCell, iCell)
            ! add on advective & diffusive fluxes
            divergenceChannel(iCell) = divergenceChannel(iCell) - channelDischarge(iEdge) * edgeSignOnCell(iEdgeOnCell, iCell)
            channelAreaChangeCell(iCell) = channelChangeRate(iEdge) * dcEdge(iEdge) * 0.5_RKIND
               ! < only half of channel is in this cell
         end do ! edges
      end do ! cells
      divergenceChannel(1:nCellsSolve) = divergenceChannel(1:nCellsSolve) / areaCell(1:nCellsSolve)
      channelAreaChangeCell(1:nCellsSolve) = channelAreaChangeCell(1:nCellsSolve) / areaCell(1:nCellsSolve)

      ! Now update channel area
      channelArea = channelChangeRate * deltatSGH + channelArea
      ! If sheet dissipation contributes to channel, there should be no need for a minimum channel size
      channelArea = max(1.0e-8_RKIND, channelArea)  ! make some tiny value when it goes negative

   !--------------------------------------------------------------------
   end subroutine evolve_channel


!***********************************************************************
!
!  routine shmip_timevarying_forcing
!
!> \brief   Calculate time-varying forcings for SHMIP experiments
!> \author  Matt Hoffman
!> \date    18 January 2017
!> \details
!>  This routine calculates time-varying forcings for the Subglacial
!>  Hydrology Model Intercomparison Project's Experiments C and D.
!>  I chose to include these in model source code rather than through
!>  a forcing module because 1) this allows the SHMIP analytic forcing
!>  functions to be reproduced exactly and 2) with the short hydro model
!>  time steps, using a forcing i/o module would require a lot of file
!>  reading during each run.
!-----------------------------------------------------------------------
   subroutine shmip_timevarying_forcing(block, err)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (block_type), intent(inout) :: block    !< Input/Output: block object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------
      integer, intent(out) :: err !< Output: error flag

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      ! Pools pointers
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: meshPool
      character (len=StrKIND), pointer :: config_SGH_shmip_forcing
      real(kind=RKIND), pointer :: daysSinceStart
      real (kind=RKIND), dimension(:), pointer :: externalWaterInput
      real (kind=RKIND), dimension(:), pointer :: basalMeltInput
      real (kind=RKIND), dimension(:), pointer :: upperSurface
      real (kind=RKIND), dimension(:), pointer :: areaCell
      ! Test case specific variables
      real(kind=RKIND) :: ra ! relative amplitude for versions of C experiment
      real(kind=RKIND) :: temperature0  ! temperature at 0m elevation for D experiment
      real(kind=RKIND) :: DT ! temperature offset for versions of D experiment

      err = 0

      call mpas_pool_get_config(liConfigs, 'config_SGH_shmip_forcing', config_SGH_shmip_forcing)

      if (config_SGH_shmip_forcing=='none') then
         ! Default is to do nothing and skip rest of routine
         return
      endif

      call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
      call mpas_pool_get_array(meshPool, 'daysSinceStart', daysSinceStart)
      call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)

      select case (config_SGH_shmip_forcing(1:1))  ! Check on first character only
      case ('C')
         ! Note: Case C assumes that moulin locations are input via externalWaterInput
         ! and have a very small value that will not affect the sum of
         ! externalWaterInput+basalMeltInput.
         ! The actual forcing is stored in basalMeltInput.  Note this is the moulin
         ! and the basal contributions combined, unlike in the other tests
         ! where I keep those separate.

         call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
         call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)

         select case (config_SGH_shmip_forcing(2:2)) ! get experiment number
         case ('1')
            ra = 0.25_RKIND
         case ('2')
            ra = 0.5_RKIND
         case ('3')
            ra = 1.0_RKIND
         case ('4')
            ra = 2.0_RKIND
         case default
            call mpas_log_write("Unknown value for config_SGH_shmip_forcing:" // config_SGH_shmip_forcing, MPAS_LOG_ERR)
            err = ior(err, 1)
         end select

         ! basalMeltInput calculated in SHMIP m/s units here
         where (externalWaterInput > 0.0_RKIND)
            ! 0.9/areaCell is the value being input into each moulin in B5
            basalMeltInput = 0.9_RKIND / areaCell * (1.0_RKIND - ra * sin(2.0_RKIND * pii * daysSinceStart))
         elsewhere
            ! where the input file does not identify a moulin, set source term to 0
            basalMeltInput = 0.0_RKIND
         end where
         ! instructions say to zero negative values
         where (basalMeltInput < 0.0_RKIND)
            basalMeltInput = 0.0_RKIND
         end where
         ! Add in basal melting contribution specified by instructions AND convert from m/s to kg/m2/s
         basalMeltInput = (basalMeltInput + 7.93e-11_RKIND) * 1000.0_RKIND

      case ('D')
         select case (config_SGH_shmip_forcing(2:2)) ! get experiment number
         case ('1')
            DT = -4.0_RKIND
         case ('2')
            DT = -2.0_RKIND
         case ('3')
            DT = 0.0_RKIND
         case ('4')
            DT = 2.0_RKIND
         case ('5')
            DT = 4.0_RKIND
         case default
            call mpas_log_write("Unknown value for config_SGH_shmip_forcing:" // config_SGH_shmip_forcing, MPAS_LOG_ERR)
            err = ior(err, 1)
         end select

         call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
         temperature0 = -16.0_RKIND * cos(2.0_RKIND * pii / 365.0 * daysSinceStart) - 5.0_RKIND + DT
         externalWaterInput(:) = (upperSurface(:) * (-0.0075_RKIND) + temperature0) * (0.01_RKIND / 86400.0_RKIND) * 1000.0_RKIND
           ! < the 1000 factor converts from m/s to kg/m2/s
         where (externalWaterInput < 0.0_RKIND)
            externalWaterInput = 0.0_RKIND
         end where

      case default
         call mpas_log_write("Unknown value for config_SGH_shmip_forcing:" // config_SGH_shmip_forcing, MPAS_LOG_ERR)
         err = ior(err, 1)
      end select


   !--------------------------------------------------------------------
   end subroutine shmip_timevarying_forcing



!***********************************************************************
!
!  routine ocean_connection_N
!
!> \brief   Calculate effective pressure assuming perfect ocean connection
!> \author  Matt Hoffman
!> \date    02 June 2020
!> \details
!>  This routine calculates effective pressure assuming a perfect hydrologic
!>  connection with the ocean.  It can be used as a simple alternative to the
!>  full subglacial hydrology.
!-----------------------------------------------------------------------
   subroutine ocean_connection_N(domain)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain  !< Input/Output: domain object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      type (block_type), pointer :: block
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: geometryPool
      real (kind=RKIND), dimension(:), pointer :: thickness
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      real (kind=RKIND), dimension(:), pointer :: effectivePressure
      real (kind=RKIND), pointer :: rhoi, rhoo

      ! Calculate N assuming perfect ocean connection
      call mpas_log_write('Calculating N assuming perfect ocean connection.')

      call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
      call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)
      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_array(geometryPool, 'thickness', thickness)
         call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
         call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
         effectivePressure = rhoi * gravity * thickness - rhoi * gravity * max(0.0_RKIND,  -1.0_RKIND * rhoo/rhoi * bedTopography)
         effectivePressure = max(effectivePressure, 0.0_RKIND) ! This is just to zero out N in the open ocean to avoid confusion

         block => block % next
      end do

   !--------------------------------------------------------------------
   end subroutine ocean_connection_N


!***********************************************************************
!
!  routine calc_hydro_mask
!
!> \brief   Calculate the boundaries of the active hydrology domain
!> \author  Matt Hoffman
!> \date    24 October 2022
!> \details
!>  This routine calculates a mask of the boundaries of the active hydrology domain
!-----------------------------------------------------------------------
   subroutine calc_hydro_mask(domain)

      !-----------------------------------------------------------------
      ! input variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! input/output variables
      !-----------------------------------------------------------------
      type (domain_type), intent(inout) :: domain  !< Input/Output: domain object

      !-----------------------------------------------------------------
      ! output variables
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! local variables
      !-----------------------------------------------------------------
      type (block_type), pointer :: block
      type (mpas_pool_type), pointer :: hydroPool
      type (mpas_pool_type), pointer :: geometryPool
      type (mpas_pool_type), pointer :: meshPool
      real (kind=RKIND), dimension(:), pointer :: bedTopography
      integer, dimension(:), pointer :: hydroMarineMarginMask
      integer, dimension(:,:), pointer :: cellsOnEdge
      integer, dimension(:), pointer :: cellMask
      integer, pointer :: nEdgesSolve
      integer :: cell1, cell2, iEdge
      real (kind=RKIND), pointer :: config_sea_level

      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)

      block => domain % blocklist
      do while (associated(block))
         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
         call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
         call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
         call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
         call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)

         hydroMarineMarginMask(:) = 0
         do iEdge = 1, nEdgesSolve
            cell1 = cellsOnEdge(1, iEdge)
            cell2 = cellsOnEdge(2, iEdge)
            ! We are looking for edges with 1 edge grounded ice and the other edge floating ice or open ocean
            if ( (li_mask_is_grounded_ice(cellMask(cell1))) .and. &
                 (li_mask_is_floating_ice(cellMask(cell2)) .or. &
                 ((bedTopography(cell2) < config_sea_level) .and. (.not. li_mask_is_ice(cellMask(cell2)))) ) ) then
               hydroMarineMarginMask(iEdge) = 1
            elseif ( (li_mask_is_grounded_ice(cellMask(cell2))) .and. &
                     (li_mask_is_floating_ice(cellMask(cell1)) .or. &
                     ((bedTopography(cell1) < config_sea_level) .and. (.not. li_mask_is_ice(cellMask(cell1)))) ) ) then
               hydroMarineMarginMask(iEdge) = 1
            endif
         enddo

         block => block % next
      end do

      call mpas_timer_start("halo updates")
      call mpas_dmpar_field_halo_exch(domain, 'hydroMarineMarginMask')
      call mpas_timer_stop("halo updates")

   !--------------------------------------------------------------------
   end subroutine calc_hydro_mask

end module li_subglacial_hydro
