MOM_ice_shelf_dynamics.F90

!> Implements a crude placeholder for a later implementation of full
!! ice shelf dynamics.
module mom_ice_shelf_dynamics
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_component, clock_routine
use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid
use mom_diag_mediator, only : diag_ctrl, time_type, enable_averaging, disable_averaging
use mom_domains, only : mom_domains_init, clone_mom_domain
use mom_domains, only : pass_var, pass_vector, to_all, cgrid_ne, bgrid_ne, corner
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
use mom_grid, only : mom_grid_init, ocean_grid_type
use mom_io, only : file_exists, slasher, mom_read_data
use mom_restart, only : register_restart_field, query_initialized
use mom_restart, only : mom_restart_cs
use mom_time_manager, only : time_type, set_time
use mom_unit_scaling, only : unit_scale_type, unit_scaling_init
!MJH use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary
use mom_ice_shelf_state, only : ice_shelf_state
use mom_coms, only : reproducing_sum, sum_across_pes, max_across_pes, min_across_pes
use mom_checksums, only : hchksum, qchksum
 
implicit none ; private
 
#include <MOM_memory.h>
 
public register_ice_shelf_dyn_restarts, initialize_ice_shelf_dyn, update_ice_shelf
public ice_time_step_cfl, ice_shelf_dyn_end
public shelf_advance_front, ice_shelf_min_thickness_calve, calve_to_mask
 
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
 
!> The control structure for the ice shelf dynamics.
type, public :: ice_shelf_dyn_cs ; private
  real, pointer, dimension(:,:) :: u_shelf => null() !< the zonal (?) velocity of the ice shelf/sheet
                                       !! on q-points (B grid) [m s-1]??
  real, pointer, dimension(:,:) :: v_shelf => null() !< the meridional velocity of the ice shelf/sheet
                                       !! on q-points (B grid) [m s-1]??
 
  real, pointer, dimension(:,:) :: u_face_mask => null() !< mask for velocity boundary conditions on the C-grid
                                       !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER,
                                       !! not vertices. Will represent boundary conditions on computational boundary
                                       !! (or permanent boundary between fast-moving and near-stagnant ice
                                       !! FOR NOW: 1=interior bdry, 0=no-flow boundary, 2=stress bdry condition,
                                       !! 3=inhomogeneous dirichlet boundary, 4=flux boundary: at these faces a flux
                                       !! will be specified which will override velocities; a homogeneous velocity
                                       !! condition will be specified (this seems to give the solver less difficulty)
  real, pointer, dimension(:,:) :: v_face_mask => null()  !< A mask for velocity boundary conditions on the C-grid
                                       !! v-face, with valued defined similarly to u_face_mask.
  real, pointer, dimension(:,:) :: u_face_mask_bdry => null() !< A duplicate copy of u_face_mask?
  real, pointer, dimension(:,:) :: v_face_mask_bdry => null() !< A duplicate copy of v_face_mask?
  real, pointer, dimension(:,:) :: u_flux_bdry_val => null() !< The ice volume flux into the cell through open boundary
                                       !! u-faces (where u_face_mask=4) [Z m2 s-1 ~> m3 s-1]??
  real, pointer, dimension(:,:) :: v_flux_bdry_val => null() !< The ice volume flux into the cell through open boundary
                                       !! v-faces (where v_face_mask=4) [Z m2 s-1 ~> m3 s-1]??
   ! needed where u_face_mask is equal to 4, similary for v_face_mask
  real, pointer, dimension(:,:) :: umask => null()      !< u-mask on the actual degrees of freedom (B grid)
                                       !! 1=normal node, 3=inhomogeneous boundary node,
                                       !!  0 - no flow node (will also get ice-free nodes)
  real, pointer, dimension(:,:) :: vmask => null()      !< v-mask on the actual degrees of freedom (B grid)
                                       !! 1=normal node, 3=inhomogeneous boundary node,
                                       !!  0 - no flow node (will also get ice-free nodes)
  real, pointer, dimension(:,:) :: calve_mask => null() !< a mask to prevent the ice shelf front from
                                          !! advancing past its initial position (but it may retreat)
  real, pointer, dimension(:,:) :: t_shelf => null() !< Veritcally integrated temperature in the ice shelf/stream,
                                                     !! on corner-points (B grid) [degC]
  real, pointer, dimension(:,:) :: tmask => null()   !< A mask on tracer points that is 1 where there is ice.
  real, pointer, dimension(:,:) :: ice_visc => null()   !< Glen's law ice viscosity, perhaps in [m].
  real, pointer, dimension(:,:) :: thickness_bdry_val => null() !< The ice thickness at an inflowing boundary [Z ~> m].
  real, pointer, dimension(:,:) :: u_bdry_val => null() !< The zonal ice velocity at inflowing boundaries [m s-1]??
  real, pointer, dimension(:,:) :: v_bdry_val => null() !< The meridional ice velocity at inflowing boundaries [m s-1]??
  real, pointer, dimension(:,:) :: h_bdry_val => null() !< The ice thickness at inflowing boundaries [m].
  real, pointer, dimension(:,:) :: t_bdry_val => null() !< The ice temperature at inflowing boundaries [degC].
 
  real, pointer, dimension(:,:) :: taub_beta_eff => null() !< nonlinear part of "linearized" basal stress.
                !!  The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011
 
  real, pointer, dimension(:,:) :: od_rt => null()         !< A running total for calculating OD_av.
  real, pointer, dimension(:,:) :: float_frac_rt => null() !< A running total for calculating float_frac.
  real, pointer, dimension(:,:) :: od_av => null()         !< The time average open ocean depth [Z ~> m].
  real, pointer, dimension(:,:) :: float_frac => null()   !< Fraction of the time a cell is "exposed", i.e. the column
                               !! thickness is below a threshold.
                       !### [if float_frac = 1 ==> grounded; obviously counterintuitive; might fix]
  integer :: od_rt_counter = 0 !< A counter of the number of contributions to OD_rt.
 
  real :: velocity_update_time_step !< The time interval over which to update the ice shelf velocity
                    !! using the nonlinear elliptic equation, or 0 to update every timestep [s].
                    ! DNGoldberg thinks this should be done no more often than about once a day
                    ! (maybe longer) because it will depend on ocean values  that are averaged over
                    ! this time interval, and solving for the equiliabrated flow will begin to lose
                    ! meaning if it is done too frequently.
  real :: elapsed_velocity_time  !< The elapsed time since the ice velocies were last udated [s].
 
  real :: g_earth      !< The gravitational acceleration [m s-2].
  real :: density_ice  !< A typical density of ice [kg m-3].
 
  logical :: gl_regularize  !< Specifies whether to regularize the floatation condition
                            !! at the grounding line as in Goldberg Holland Schoof 2009
  integer :: n_sub_regularize
                            !< partition of cell over which to integrate for
                            !! interpolated grounding line the (rectangular) is
                            !! divided into nxn equally-sized rectangles, over which
                            !!  basal contribution is integrated (iterative quadrature)
  logical :: gl_couple      !< whether to let the floatation condition be
                            !! determined by ocean column thickness means update_OD_ffrac
                            !! will be called (note: GL_regularize and GL_couple
                            !! should be exclusive)
 
  real    :: cfl_factor     !< A factor used to limit subcycled advective timestep in uncoupled runs
                            !! i.e. dt <= CFL_factor * min(dx / u)
 
  real :: a_glen_isothermal !< Ice viscosity parameter in Glen's Lawa, [Pa-1/3 year].
  real :: n_glen            !< Nonlinearity exponent in Glen's Law
  real :: eps_glen_min      !< Min. strain rate to avoid infinite Glen's law viscosity, [year-1].
  real :: c_basal_friction  !< Ceofficient in sliding law tau_b = C u^(n_basal_friction), in
                            !!  units="Pa (m-a)-(n_basal_friction)
  real :: n_basal_friction  !< Exponent in sliding law tau_b = C u^(m_slide)
  real :: density_ocean_avg !< This does not affect ocean circulation or thermodynamics.
                            !! It is used to estimate the gravitational driving force at the
                            !! shelf front (until we think of a better way to do it,
                            !! but any difference will be negligible).
  real :: thresh_float_col_depth !< The water column depth over which the shelf if considered to be floating
  logical :: moving_shelf_front  !< Specify whether to advance shelf front (and calve).
  logical :: calve_to_mask       !< If true, calve off the ice shelf when it passes the edge of a mask.
  real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
 
  real :: cg_tolerance !< The tolerance in the CG solver, relative to initial residual, that
                       !! deterimnes when to stop the conguage gradient iterations.
  real :: nonlinear_tolerance !< The fractional nonlinear tolerance, relative to the initial error,
                              !! that sets when to stop the iterative velocity solver
  integer :: cg_max_iterations !< The maximum number of iterations that can be used in the CG solver
  integer :: nonlin_solve_err_mode  !< 1: exit vel solve based on nonlin residual
                    !! 2: exit based on "fixed point" metric (|u - u_last| / |u| < tol where | | is infty-norm
  logical :: use_reproducing_sums !< Use reproducing global sums.
 
  ! ids for outputting intermediate thickness in advection subroutine (debugging)
  !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1
 
  logical :: debug                !< If true, write verbose checksums for debugging purposes
                                  !! and use reproducible sums
  logical :: module_is_initialized = .false. !< True if this module has been initialized.
 
  !>@{ Diagnostic handles
  integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, &
             id_float_frac = -1, id_col_thick = -1, id_od_av = -1, &
             id_u_mask = -1, id_v_mask = -1, id_t_mask = -1
  !!@}
  ! ids for outputting intermediate thickness in advection subroutine (debugging)
  !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1
 
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to control diagnostic output.
 
end type ice_shelf_dyn_cs
 
contains
 
!> used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia)
function slope_limiter(num, denom)
  real, intent(in)    :: num   !< The numerator of the ratio used in the Van Leer slope limiter
  real, intent(in)    :: denom !< The denominator of the ratio used in the Van Leer slope limiter
  real :: slope_limiter
  real :: r
 
  if (denom == 0) then
    slope_limiter = 0
  elseif (num*denom <= 0) then
    slope_limiter = 0
  else
    r = num/denom
    slope_limiter = (r+abs(r))/(1+abs(r))
  endif
 
end function slope_limiter
 
!> Calculate area of quadrilateral.
function quad_area (X, Y)
  real, dimension(4), intent(in) :: x !< The x-positions of the vertices of the quadrilateral.
  real, dimension(4), intent(in) :: y !< The y-positions of the vertices of the quadrilateral.
  real :: quad_area, p2, q2, a2, c2, b2, d2
 
! X and Y must be passed in the form
    !  3 - 4
    !  |   |
    !  1 - 2
 
  p2 = (x(4)-x(1))**2 + (y(4)-y(1))**2 ; q2 = (x(3)-x(2))**2 + (y(3)-y(2))**2
  a2 = (x(3)-x(4))**2 + (y(3)-y(4))**2 ; c2 = (x(1)-x(2))**2 + (y(1)-y(2))**2
  b2 = (x(2)-x(4))**2 + (y(2)-y(4))**2 ; d2 = (x(3)-x(1))**2 + (y(3)-y(1))**2
  quad_area = .25 * sqrt(4*p2*q2-(b2+d2-a2-c2)**2)
 
end function quad_area
 
!> This subroutine is used to register any fields related to the ice shelf
!! dynamics that should be written to or read from the restart file.
subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS)
  type(ocean_grid_type),  intent(inout) :: g    !< The grid type describing the ice shelf grid.
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(ice_shelf_dyn_cs), pointer       :: cs !< A pointer to the ice shelf dynamics control structure
  type(mom_restart_cs),   pointer       :: restart_cs !< A pointer to the restart control structure.
 
  logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
  character(len=40)  :: mdl = "MOM_ice_shelf_dyn"  ! This module's name.
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  if (associated(cs)) then
    call mom_error(fatal, "MOM_ice_shelf_dyn.F90, register_ice_shelf_dyn_restarts: "// &
                          "called with an associated control structure.")
    return
  endif
  allocate(cs)
 
  override_shelf_movement = .false. ; active_shelf_dynamics = .false.
  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
                 "If true, the ice sheet mass can evolve with time.", &
                 default=.false., do_not_log=.true.)
  if (shelf_mass_is_dynamic) then
    call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
                 "If true, user provided code specifies the ice-shelf "//&
                 "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
    active_shelf_dynamics = .not.override_shelf_movement
  endif
 
  if (active_shelf_dynamics) then
    allocate( cs%u_shelf(isdb:iedb,jsdb:jedb) ) ; cs%u_shelf(:,:) = 0.0
    allocate( cs%v_shelf(isdb:iedb,jsdb:jedb) ) ; cs%v_shelf(:,:) = 0.0
    allocate( cs%t_shelf(isd:ied,jsd:jed) )   ; cs%t_shelf(:,:) = -10.0
    allocate( cs%ice_visc(isd:ied,jsd:jed) )    ; cs%ice_visc(:,:) = 0.0
    allocate( cs%taub_beta_eff(isd:ied,jsd:jed) ) ; cs%taub_beta_eff(:,:) = 0.0
    allocate( cs%OD_av(isd:ied,jsd:jed) )       ; cs%OD_av(:,:) = 0.0
    allocate( cs%float_frac(isd:ied,jsd:jed) )  ; cs%float_frac(:,:) = 0.0
 
    ! additional restarts for ice shelf state
    call register_restart_field(cs%u_shelf, "u_shelf", .false., restart_cs, &
                                "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu')
    call register_restart_field(cs%v_shelf, "v_shelf", .false., restart_cs, &
                                "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu')
    call register_restart_field(cs%t_shelf, "t_shelf", .true., restart_cs, &
                                "ice sheet/shelf vertically averaged temperature", "deg C")
    call register_restart_field(cs%OD_av, "OD_av", .true., restart_cs, &
                                "Average open ocean depth in a cell","m")
    call register_restart_field(cs%float_frac, "float_frac", .true., restart_cs, &
                                "fractional degree of grounding", "nondim")
    call register_restart_field(cs%ice_visc, "viscosity", .true., restart_cs, &
                                "Glens law ice viscosity", "m (seems wrong)")
    call register_restart_field(cs%taub_beta_eff, "tau_b_beta", .true., restart_cs, &
                                "Coefficient of basal traction", "m (seems wrong)")
  endif
 
end subroutine register_ice_shelf_dyn_restarts
 
!> Initializes shelf model data, parameters and diagnostics
subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_sim, solo_ice_sheet_in)
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(ocean_grid_type),   pointer       :: ocn_grid   !< The calling ocean model's horizontal grid structure
  type(time_type),         intent(inout) :: time !< The clock that that will indicate the model time
  type(ice_shelf_state),   intent(in)    :: iss  !< A structure with elements that describe
                                                 !! the ice-shelf state
  type(ice_shelf_dyn_cs),  pointer       :: cs   !< A pointer to the ice shelf dynamics control structure
  type(ocean_grid_type),   intent(inout) :: g    !< The grid type describing the ice shelf grid.
  type(unit_scale_type),   intent(in)    :: us   !< Pointer to a structure containing unit conversion factors
  type(diag_ctrl), target, intent(in)    :: diag !< A structure that is used to regulate the diagnostic output.
  logical,                 intent(in)    :: new_sim !< If true this is a new simulation, otherwise
                                                 !! has been started from a restart file.
  logical,       optional, intent(in)    :: solo_ice_sheet_in !< If present, this indicates whether
                                                 !! a solo ice-sheet driver.
 
  ! Local variables
  real    :: z_rescale  ! A rescaling factor for heights from the representation in
                        ! a reastart fole to the internal representation in this run.
  !This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=200) :: config
  character(len=200) :: ic_file,filename,inputdir
  character(len=40)  :: var_name
  character(len=40)  :: mdl = "MOM_ice_shelf_dyn"  ! This module's name.
  logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
  logical :: debug
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, isdq, iedq, jsdq, jedq, iters
 
  if (.not.associated(cs)) then
    call mom_error(fatal, "MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn: "// &
                          "called with an associated control structure.")
    return
  endif
  if (cs%module_is_initialized) then
    call mom_error(warning, "MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn was "//&
             "called with a control structure that has already been initialized.")
  endif
  cs%module_is_initialized = .true.
 
  cs%diag => diag ! ; CS%Time => Time
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
  call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
                 "If true, write verbose debugging messages for the ice shelf.", &
                 default=debug)
  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
                 "If true, the ice sheet mass can evolve with time.", &
                 default=.false.)
  override_shelf_movement = .false. ; active_shelf_dynamics = .false.
  if (shelf_mass_is_dynamic) then
    call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
                 "If true, user provided code specifies the ice-shelf "//&
                 "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
    active_shelf_dynamics = .not.override_shelf_movement
 
    call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
                 "If true, regularize the floatation condition at the "//&
                 "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
    call get_param(param_file, mdl, "GROUNDING_LINE_INTERP_SUBGRID_N", cs%n_sub_regularize, &
                 "The number of sub-partitions of each cell over which to "//&
                 "integrate for the interpolated grounding line. Each cell "//&
                 "is divided into NxN equally-sized rectangles, over which the "//&
                 "basal contribution is integrated by iterative quadrature.", &
                 default=0)
    call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
                 "If true, let the floatation condition be determined by "//&
                 "ocean column thickness. This means that update_OD_ffrac "//&
                 "will be called.  GL_REGULARIZE and GL_COUPLE are exclusive.", &
                 default=.false., do_not_log=cs%GL_regularize)
    if (cs%GL_regularize) cs%GL_couple = .false.
    if (cs%GL_regularize .and. (cs%n_sub_regularize == 0)) call mom_error (fatal, &
      "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used")
    call get_param(param_file, mdl, "ICE_SHELF_CFL_FACTOR", cs%CFL_factor, &
                 "A factor used to limit timestep as CFL_FACTOR * min (\Delta x / u). "//&
                 "This is only used with an ice-only model.", default=0.25)
  endif
  call get_param(param_file, mdl, "RHO_0", cs%density_ocean_avg, &
                 "avg ocean density used in floatation cond", &
                 units="kg m-3", default=1035.)
  if (active_shelf_dynamics) then
    call get_param(param_file, mdl, "ICE_VELOCITY_TIMESTEP", cs%velocity_update_time_step, &
                 "seconds between ice velocity calcs", units="s", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80)
 
    call get_param(param_file, mdl, "A_GLEN_ISOTHERM", cs%A_glen_isothermal, &
                 "Ice viscosity parameter in Glen's Law", &
                 units="Pa -1/3 a", default=9.461e-18)
    call get_param(param_file, mdl, "GLEN_EXPONENT", cs%n_glen, &
                 "nonlinearity exponent in Glen's Law", &
                  units="none", default=3.)
    call get_param(param_file, mdl, "MIN_STRAIN_RATE_GLEN", cs%eps_glen_min, &
                 "min. strain rate to avoid infinite Glen's law viscosity", &
                 units="a-1", default=1.e-12)
    call get_param(param_file, mdl, "BASAL_FRICTION_COEFF", cs%C_basal_friction, &
                 "ceofficient in sliding law \tau_b = C u^(n_basal_friction)", &
                 units="Pa (m-a)-(n_basal_friction)", fail_if_missing=.true.)
    call get_param(param_file, mdl, "BASAL_FRICTION_EXP", cs%n_basal_friction, &
                 "exponent in sliding law \tau_b = C u^(m_slide)", &
                 units="none", fail_if_missing=.true.)
    call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
                 "A typical density of ice.", units="kg m-3", default=917.0)
    call get_param(param_file, mdl, "CONJUGATE_GRADIENT_TOLERANCE", cs%cg_tolerance, &
                "tolerance in CG solver, relative to initial residual", default=1.e-6)
    call get_param(param_file, mdl, "ICE_NONLINEAR_TOLERANCE", cs%nonlinear_tolerance, &
                "nonlin tolerance in iterative velocity solve",default=1.e-6)
    call get_param(param_file, mdl, "CONJUGATE_GRADIENT_MAXIT", cs%cg_max_iterations, &
                "max iteratiions in CG solver", default=2000)
    call get_param(param_file, mdl, "THRESH_FLOAT_COL_DEPTH", cs%thresh_float_col_depth, &
                "min ocean thickness to consider ice *floating*; "//&
                "will only be important with use of tides", &
                units="m", default=1.e-3, scale=us%m_to_Z)
    call get_param(param_file, mdl, "NONLIN_SOLVE_ERR_MODE", cs%nonlin_solve_err_mode, &
                "Choose whether nonlin error in vel solve is based on nonlinear "//&
                "residual (1) or relative change since last iteration (2)", default=1)
    call get_param(param_file, mdl, "SHELF_DYN_REPRODUCING_SUMS", cs%use_reproducing_sums, &
                "If true, use the reproducing extended-fixed-point sums in "//&
                "the ice shelf dynamics solvers.", default=.true.)
 
    call get_param(param_file, mdl, "SHELF_MOVING_FRONT", cs%moving_shelf_front, &
                 "Specify whether to advance shelf front (and calve).", &
                 default=.true.)
    call get_param(param_file, mdl, "CALVE_TO_MASK", cs%calve_to_mask, &
                 "If true, do not allow an ice shelf where prohibited by a mask.", &
                 default=.false.)
  endif
  call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", &
                 cs%min_thickness_simple_calve, &
                 "Min thickness rule for the VERY simple calving law",&
                 units="m", default=0.0, scale=us%m_to_Z)
 
  ! Allocate memory in the ice shelf dynamics control structure that was not
  ! previously allocated for registration for restarts.
  ! OVS vertically integrated Temperature
 
  if (active_shelf_dynamics) then
    ! DNG
    allocate( cs%u_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%u_bdry_val(:,:) = 0.0
    allocate( cs%v_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%v_bdry_val(:,:) = 0.0
    allocate( cs%t_bdry_val(isd:ied,jsd:jed) )   ; cs%t_bdry_val(:,:) = -15.0
    allocate( cs%h_bdry_val(isd:ied,jsd:jed) ) ; cs%h_bdry_val(:,:) = 0.0
    allocate( cs%thickness_bdry_val(isd:ied,jsd:jed) ) ; cs%thickness_bdry_val(:,:) = 0.0
    allocate( cs%u_face_mask(isdq:iedq,jsd:jed) ) ; cs%u_face_mask(:,:) = 0.0
    allocate( cs%v_face_mask(isd:ied,jsdq:jedq) ) ; cs%v_face_mask(:,:) = 0.0
    allocate( cs%u_face_mask_bdry(isdq:iedq,jsd:jed) ) ; cs%u_face_mask_bdry(:,:) = -2.0
    allocate( cs%v_face_mask_bdry(isd:ied,jsdq:jedq) ) ; cs%v_face_mask_bdry(:,:) = -2.0
    allocate( cs%u_flux_bdry_val(isdq:iedq,jsd:jed) ) ; cs%u_flux_bdry_val(:,:) = 0.0
    allocate( cs%v_flux_bdry_val(isd:ied,jsdq:jedq) ) ; cs%v_flux_bdry_val(:,:) = 0.0
    allocate( cs%umask(isdq:iedq,jsdq:jedq) ) ; cs%umask(:,:) = -1.0
    allocate( cs%vmask(isdq:iedq,jsdq:jedq) ) ; cs%vmask(:,:) = -1.0
    allocate( cs%tmask(isdq:iedq,jsdq:jedq) ) ; cs%tmask(:,:) = -1.0
 
    cs%OD_rt_counter = 0
    allocate( cs%OD_rt(isd:ied,jsd:jed) ) ; cs%OD_rt(:,:) = 0.0
    allocate( cs%float_frac_rt(isd:ied,jsd:jed) ) ; cs%float_frac_rt(:,:) = 0.0
 
    if (cs%calve_to_mask) then
      allocate( cs%calve_mask(isd:ied,jsd:jed) ) ; cs%calve_mask(:,:) = 0.0
    endif
 
    cs%elapsed_velocity_time = 0.0
 
    call update_velocity_masks(cs, g, iss%hmask, cs%umask, cs%vmask, cs%u_face_mask, cs%v_face_mask)
  endif
 
  ! Take additional initialization steps, for example of dependent variables.
  if (active_shelf_dynamics .and. .not.new_sim) then
    if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) then
      z_rescale = us%m_to_Z / us%m_to_Z_restart
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        cs%OD_av(i,j) = z_rescale * cs%OD_av(i,j)
      enddo ; enddo
    endif
 
    ! this is unfortunately necessary; if grid is not symmetric the boundary values
    !  of u and v are otherwise not set till the end of the first linear solve, and so
    !  viscosity is not calculated correctly.
    ! This has to occur after init_boundary_values or some of the arrays on the
    ! right hand side have not been set up yet.
    if (.not. g%symmetric) then
      do j=g%jsd,g%jed ; do i=g%isd,g%ied
        if (((i+g%idg_offset) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3)) then
          cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
          cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
        endif
        if (((j+g%jdg_offset) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3)) then
          cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
          cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
        endif
      enddo ; enddo
    endif
 
    call pass_var(cs%OD_av,g%domain)
    call pass_var(cs%float_frac,g%domain)
    call pass_var(cs%ice_visc,g%domain)
    call pass_var(cs%taub_beta_eff,g%domain)
    call pass_vector(cs%u_shelf, cs%v_shelf, g%domain, to_all, bgrid_ne)
  endif
 
  if (active_shelf_dynamics) then
    ! If we are calving to a mask, i.e. if a mask exists where a shelf cannot, read the mask from a file.
    if (cs%calve_to_mask) then
      call mom_mesg("  MOM_ice_shelf.F90, initialize_ice_shelf: reading calving_mask")
 
      call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
      inputdir = slasher(inputdir)
      call get_param(param_file, mdl, "CALVING_MASK_FILE", ic_file, &
                   "The file with a mask for where calving might occur.", &
                   default="ice_shelf_h.nc")
      call get_param(param_file, mdl, "CALVING_MASK_VARNAME", var_name, &
                   "The variable to use in masking calving.", &
                   default="area_shelf_h")
 
      filename = trim(inputdir)//trim(ic_file)
      call log_param(param_file, mdl, "INPUTDIR/CALVING_MASK_FILE", filename)
      if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
         " calving mask file: Unable to open "//trim(filename))
 
      call mom_read_data(filename,trim(var_name),cs%calve_mask,g%Domain)
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        if (cs%calve_mask(i,j) > 0.0) cs%calve_mask(i,j) = 1.0
      enddo ; enddo
      call pass_var(cs%calve_mask,g%domain)
    endif
 
!    call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim)
 
    if (new_sim) then
      call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.")
      call update_od_ffrac_uncoupled(cs, g, iss%h_shelf(:,:))
      call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, time)
 
      if (cs%id_u_shelf > 0) call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
      if (cs%id_v_shelf > 0) call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
    endif
 
  ! Register diagnostics.
    cs%id_u_shelf = register_diag_field('ocean_model','u_shelf',cs%diag%axesCu1, time, &
       'x-velocity of ice', 'm yr-1')
    cs%id_v_shelf = register_diag_field('ocean_model','v_shelf',cs%diag%axesCv1, time, &
       'y-velocity of ice', 'm yr-1')
    cs%id_u_mask = register_diag_field('ocean_model','u_mask',cs%diag%axesCu1, time, &
       'mask for u-nodes', 'none')
    cs%id_v_mask = register_diag_field('ocean_model','v_mask',cs%diag%axesCv1, time, &
       'mask for v-nodes', 'none')
!    CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1, Time, &
!       'ice surf elev', 'm')
    cs%id_float_frac = register_diag_field('ocean_model','ice_float_frac',cs%diag%axesT1, time, &
       'fraction of cell that is floating (sort of)', 'none')
    cs%id_col_thick = register_diag_field('ocean_model','col_thick',cs%diag%axesT1, time, &
       'ocean column thickness passed to ice model', 'm', conversion=us%Z_to_m)
    cs%id_OD_av = register_diag_field('ocean_model','OD_av',cs%diag%axesT1, time, &
       'intermediate ocean column thickness passed to ice model', 'm', conversion=us%Z_to_m)
    !CS%id_h_after_uflux = register_diag_field('ocean_model','h_after_uflux',CS%diag%axesh1, Time, &
    !   'thickness after u flux ', 'none')
    !CS%id_h_after_vflux = register_diag_field('ocean_model','h_after_vflux',CS%diag%axesh1, Time, &
    !   'thickness after v flux ', 'none')
    !CS%id_h_after_adv = register_diag_field('ocean_model','h_after_adv',CS%diag%axesh1, Time, &
    !   'thickness after front adv ', 'none')
 
!!! OVS vertically integrated temperature
    cs%id_t_shelf = register_diag_field('ocean_model','t_shelf',cs%diag%axesT1, time, &
       'T of ice', 'oC')
    cs%id_t_mask = register_diag_field('ocean_model','tmask',cs%diag%axesT1, time, &
       'mask for T-nodes', 'none')
  endif
 
end subroutine initialize_ice_shelf_dyn
 
 
subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time)
  type(ice_shelf_dyn_cs), intent(inout) :: CS  !< A pointer to the ice shelf control structure
  type(ice_shelf_state),  intent(in)    :: ISS !< A structure with elements that describe
                                               !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: G   !< The grid structure used by the ice shelf.
  type(unit_scale_type),  intent(in)    :: US  !< Pointer to a structure containing unit conversion factors
  type(time_type),        intent(in)    :: Time !< The current model time
 
  integer         :: i, j, iters, isd, ied, jsd, jed
  real            :: rhoi_rhow, OD
  type(time_type) :: dummy_time
 
  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
  dummy_time = set_time(0,0)
  isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
 
  do j=jsd,jed
    do i=isd,ied
      od = g%bathyT(i,j) - rhoi_rhow * iss%h_shelf(i,j)
      if (od >= 0) then
    ! ice thickness does not take up whole ocean column -> floating
        cs%OD_av(i,j) = od
        cs%float_frac(i,j) = 0.
      else
        cs%OD_av(i,j) = 0.
        cs%float_frac(i,j) = 1.
      endif
    enddo
  enddo
 
  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, dummy_time)
 
end subroutine initialize_diagnostic_fields
 
!> This function returns the global maximum timestep that can be taken based on the current
!! ice velocities.  Because it involves finding a global minimum, it can be suprisingly expensive.
function ice_time_step_cfl(CS, ISS, G)
  type(ice_shelf_dyn_cs), intent(inout) :: cs  !< The ice shelf dynamics control structure
  type(ice_shelf_state),  intent(inout) :: iss !< A structure with elements that describe
                                               !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: g   !< The grid structure used by the ice shelf.
  real :: ice_time_step_cfl !< The maximum permitted timestep based on the ice velocities [s].
 
  real :: ratio, min_ratio
  real :: local_u_max, local_v_max
  integer :: i, j
 
  min_ratio = 1.0e16 ! This is just an arbitrary large value.
  do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (iss%hmask(i,j) == 1.0) then
    local_u_max = max(abs(cs%u_shelf(i,j)), abs(cs%u_shelf(i+1,j+1)), &
                      abs(cs%u_shelf(i+1,j)), abs(cs%u_shelf(i,j+1)))
    local_v_max = max(abs(cs%v_shelf(i,j)), abs(cs%v_shelf(i+1,j+1)), &
                      abs(cs%v_shelf(i+1,j)), abs(cs%v_shelf(i,j+1)))
 
    ratio = min(g%areaT(i,j) / (local_u_max+1.0e-12), g%areaT(i,j) / (local_v_max+1.0e-12))
    min_ratio = min(min_ratio, ratio)
  endif ; enddo ; enddo ! i- and j- loops
 
  call min_across_pes(min_ratio)
 
  ! solved velocities are in m/yr; we want time_step_int in seconds
  ice_time_step_cfl = cs%CFL_factor * min_ratio * (365*86400)
 
end function ice_time_step_cfl
 
!> This subroutine updates the ice shelf velocities, mass, stresses and properties due to the
!! ice shelf dynamics.
subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled_grounding, must_update_vel)
  type(ice_shelf_dyn_cs), intent(inout) :: cs !< The ice shelf dynamics control structure
  type(ice_shelf_state),  intent(inout) :: iss !< A structure with elements that describe
                                              !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: g  !< The grid structure used by the ice shelf.
  type(unit_scale_type),  intent(in)    :: us !< Pointer to a structure containing unit conversion factors
  real,                   intent(in)    :: time_step !< time step [s]
  type(time_type),        intent(in)    :: time !< The current model time
  real, dimension(SZDI_(G),SZDJ_(G)), &
                optional, intent(in)    :: ocean_mass !< If present this is the mass per unit area
                                              !! of the ocean [kg m-2].
  logical,      optional, intent(in)    :: coupled_grounding !< If true, the grounding line is
                                              !! determined by coupled ice-ocean dynamics
  logical,      optional, intent(in)    :: must_update_vel !< Always update the ice velocities if true.
 
  integer :: iters
  logical :: update_ice_vel, coupled_gl
 
  update_ice_vel = .false.
  if (present(must_update_vel)) update_ice_vel = must_update_vel
 
  coupled_gl = .false.
  if (present(ocean_mass) .and. present(coupled_grounding)) coupled_gl = coupled_grounding
 
  call ice_shelf_advect(cs, iss, g, time_step, time)
  cs%elapsed_velocity_time = cs%elapsed_velocity_time + time_step
  if (cs%elapsed_velocity_time >= cs%velocity_update_time_step) update_ice_vel = .true.
 
  if (coupled_gl) then
    call update_od_ffrac(cs, g, us, ocean_mass, update_ice_vel)
  elseif (update_ice_vel) then
    call update_od_ffrac_uncoupled(cs, g, iss%h_shelf(:,:))
  endif
 
  if (update_ice_vel) then
    call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, time)
  endif
 
  call ice_shelf_temp(cs, iss, g, us, time_step, iss%water_flux, time)
 
  if (update_ice_vel) then
    call enable_averaging(cs%elapsed_velocity_time, time, cs%diag)
    if (cs%id_col_thick > 0) call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
    if (cs%id_u_shelf > 0) call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
    if (cs%id_v_shelf > 0) call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
    if (cs%id_t_shelf > 0) call post_data(cs%id_t_shelf,cs%t_shelf,cs%diag)
    if (cs%id_float_frac > 0) call post_data(cs%id_float_frac,cs%float_frac,cs%diag)
    if (cs%id_OD_av >0) call post_data(cs%id_OD_av, cs%OD_av,cs%diag)
 
    if (cs%id_u_mask > 0) call post_data(cs%id_u_mask,cs%umask,cs%diag)
    if (cs%id_v_mask > 0) call post_data(cs%id_v_mask,cs%vmask,cs%diag)
    if (cs%id_t_mask > 0) call post_data(cs%id_t_mask,cs%tmask,cs%diag)
 
    call disable_averaging(cs%diag)
 
    cs%elapsed_velocity_time = 0.0
  endif
 
end subroutine update_ice_shelf
 
!> This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once.
!! Additionally, it will update the volume of ice in partially-filled cells, and update
!! hmask accordingly
subroutine ice_shelf_advect(CS, ISS, G, time_step, Time)
  type(ice_shelf_dyn_cs), intent(inout) :: CS !< The ice shelf dynamics control structure
  type(ice_shelf_state),  intent(inout) :: ISS !< A structure with elements that describe
                                               !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: G  !< The grid structure used by the ice shelf.
  real,                   intent(in)    :: time_step !< time step [s]
  type(time_type),        intent(in)    :: Time !< The current model time
 
 
! 3/8/11 DNG
!
!    This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once.
!    ADDITIONALLY, it will update the volume of ice in partially-filled cells, and update
!        hmask accordingly
!
!    The flux overflows are included here. That is because they will be used to advect 3D scalars
!    into partial cells
 
  !
  ! flux_enter: this is to capture flow into partially covered cells; it gives the mass flux into a given
  ! cell across its boundaries.
  ! ###Perhaps flux_enter should be changed into u-face and v-face
  ! ###fluxes, which can then be used in halo updates, etc.
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !  THESE ARE NOT CONSISTENT ==> FIND OUT WHAT YOU IMPLEMENTED
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   o--- (4) ---o
  !   |           |
  !  (1)         (2)
  !   |           |
  !   o--- (3) ---o
  !
 
  real, dimension(SZDI_(G),SZDJ_(G))   :: h_after_uflux, h_after_vflux ! Ice thicknesses [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
  integer                           :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
  real                              :: rho, spy
 
  rho = cs%density_ice
  spy = 365 * 86400 ! seconds per year; is there a global constant for this?  No - it is dependent upon a calendar.
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
  flux_enter(:,:,:) = 0.0
 
  h_after_uflux(:,:) = 0.0
  h_after_vflux(:,:) = 0.0
  ! call MOM_mesg("MOM_ice_shelf.F90: ice_shelf_advect called")
 
  do j=jsd,jed ; do i=isd,ied ; if (cs%thickness_bdry_val(i,j) /= 0.0) then
    iss%h_shelf(i,j) = cs%thickness_bdry_val(i,j)
  endif ; enddo ; enddo
 
  call ice_shelf_advect_thickness_x(cs, g, time_step/spy, iss%hmask, iss%h_shelf, h_after_uflux, flux_enter)
 
!  call enable_averaging(time_step,Time,CS%diag)
 ! call pass_var(h_after_uflux, G%domain)
!  if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
!  call disable_averaging(CS%diag)
 
  call ice_shelf_advect_thickness_y(cs, g, time_step/spy, iss%hmask, h_after_uflux, h_after_vflux, flux_enter)
 
!  call enable_averaging(time_step,Time,CS%diag)
!  call pass_var(h_after_vflux, G%domain)
!  if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
!  call disable_averaging(CS%diag)
 
  do j=jsd,jed
    do i=isd,ied
      if (iss%hmask(i,j) == 1) iss%h_shelf(i,j) = h_after_vflux(i,j)
    enddo
  enddo
 
  if (cs%moving_shelf_front) then
    call shelf_advance_front(cs, iss, g, flux_enter)
    if (cs%min_thickness_simple_calve > 0.0) then
      call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
                                         cs%min_thickness_simple_calve)
    endif
    if (cs%calve_to_mask) then
      call calve_to_mask(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, cs%calve_mask)
    endif
  endif
 
  !call enable_averaging(time_step,Time,CS%diag)
  !if (CS%id_h_after_adv > 0) call post_data(CS%id_h_after_adv, ISS%h_shelf, CS%diag)
  !call disable_averaging(CS%diag)
 
  !call change_thickness_using_melt(ISS, G, time_step, fluxes, CS%density_ice)
 
  call update_velocity_masks(cs, g, iss%hmask, cs%umask, cs%vmask, cs%u_face_mask, cs%v_face_mask)
 
end subroutine ice_shelf_advect
 
subroutine ice_shelf_solve_outer(CS, ISS, G, US, u, v, iters, time)
  type(ice_shelf_dyn_cs), intent(inout) :: CS !< The ice shelf dynamics control structure
  type(ice_shelf_state),  intent(in)    :: ISS !< A structure with elements that describe
                                               !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: G  !< The grid structure used by the ice shelf.
  type(unit_scale_type),  intent(in)    :: US !< Pointer to a structure containing unit conversion factors
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: u  !< The zonal ice shelf velocity at vertices [m year-1]
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: v  !< The meridional ice shelf velocity at vertices [m year-1]
  integer,                intent(out)   :: iters !< The number of iterations used in the solver.
  type(time_type),        intent(in)    :: Time !< The current model time
 
  real, dimension(SZDIB_(G),SZDJB_(G)) :: TAUDX, TAUDY, u_prev_iterate, v_prev_iterate, &
                        u_bdry_cont, v_bdry_cont, Au, Av, err_u, err_v, &
                        u_last, v_last
  real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! An array indicating where the ice
                                                ! shelf is floating: 0 if floating, 1 if not.
  character(len=160) :: mesg  ! The text of an error message
  integer :: conv_flag, i, j, k,l, iter
  integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
  real                     :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi_rhow
  real, pointer, dimension(:,:,:,:) :: Phi => null()
  real, pointer, dimension(:,:,:,:,:,:) :: Phisub => null()
  real, dimension(8,4)       :: Phi_temp
  real, dimension(2,2)       :: X,Y
  character(2)                :: iternum
  character(2)                :: numproc
 
  ! for GL interpolation - need to make this a readable parameter
  nsub = cs%n_sub_regularize
 
  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
 
  taudx(:,:) = 0.0 ; taudy(:,:) = 0.0
  u_bdry_cont(:,:) = 0.0 ; v_bdry_cont(:,:) = 0.0
  au(:,:) = 0.0 ; av(:,:) = 0.0
 
  ! need to make these conditional on GL interpolation
  float_cond(:,:) = 0.0 ; h_node(:,:)=0
  allocate(phisub(nsub,nsub,2,2,2,2)) ; phisub = 0.0
 
  isumstart = g%isc
  ! Include the edge if tile is at the western bdry;  Should add a test to avoid this if reentrant.
  if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
 
  jsumstart = g%jsc
  ! Include the edge if tile is at the southern bdry;  Should add a test to avoid this if reentrant.
  if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
 
  call calc_shelf_driving_stress(cs, iss, g, us, taudx, taudy, cs%OD_av)
 
  ! this is to determine which cells contain the grounding line,
  !  the criterion being that the cell is ice-covered, with some nodes
  !  floating and some grounded
  ! floatation condition is estimated by assuming topography is cellwise constant
  !  and H is bilinear in a cell; floating where rho_i/rho_w * H_node + D is nonpositive
 
  ! need to make this conditional on GL interp
 
  if (cs%GL_regularize) then
 
    call interpolate_h_to_b(g, iss%h_shelf, iss%hmask, h_node)
 
    do j=g%jsc,g%jec
      do i=g%isc,g%iec
        nodefloat = 0
        do k=0,1
          do l=0,1
            if ((iss%hmask(i,j) == 1) .and. &
              (rhoi_rhow * h_node(i-1+k,j-1+l) - g%bathyT(i,j) <= 0)) then
              nodefloat = nodefloat + 1
            endif
          enddo
        enddo
        if ((nodefloat > 0) .and. (nodefloat < 4)) then
          float_cond(i,j) = 1.0
          cs%float_frac(i,j) = 1.0
        endif
      enddo
    enddo
 
    call pass_var(float_cond, g%Domain)
 
    call bilinear_shape_functions_subgrid(phisub, nsub)
 
  endif
 
  ! make above conditional
 
  u_prev_iterate(:,:) = u(:,:)
  v_prev_iterate(:,:) = v(:,:)
 
  ! must prepare phi
  allocate(phi(isd:ied,jsd:jed,1:8,1:4)) ; phi(:,:,:,:) = 0.0
 
  do j=jsd,jed ; do i=isd,ied
    if (((i > isd) .and. (j > jsd))) then
      x(:,:) = g%geoLonBu(i-1:i,j-1:j)*1000
      y(:,:) = g%geoLatBu(i-1:i,j-1:j)*1000
    else
      x(2,:) = g%geoLonBu(i,j)*1000
      x(1,:) = g%geoLonBu(i,j)*1000-g%dxT(i,j)
      y(:,2) = g%geoLatBu(i,j)*1000
      y(:,1) = g%geoLatBu(i,j)*1000-g%dyT(i,j)
    endif
 
    call bilinear_shape_functions(x, y, phi_temp, area)
    phi(i,j,:,:) = phi_temp
  enddo ; enddo
 
  call calc_shelf_visc(cs, iss, g, us, u, v)
 
  call pass_var(cs%ice_visc, g%domain)
  call pass_var(cs%taub_beta_eff, g%domain)
 
  ! makes sure basal stress is only applied when it is supposed to be
 
  do j=g%jsd,g%jed ; do i=g%isd,g%ied
    cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
  enddo ; enddo
 
  call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
                                      cs%taub_beta_eff, float_cond, &
                                      rhoi_rhow, u_bdry_cont, v_bdry_cont)
 
  au(:,:) = 0.0 ; av(:,:) = 0.0
 
  call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
            cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%areaT, &
            g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
 
  err_init = 0 ; err_tempu = 0; err_tempv = 0
  do j=jsumstart,g%jecB
    do i=isumstart,g%iecB
      if (cs%umask(i,j) == 1) then
        err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
      endif
      if (cs%vmask(i,j) == 1) then
        err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
      endif
      if (err_tempv >= err_init) then
        err_init = err_tempv
      endif
    enddo
  enddo
 
  call max_across_pes(err_init)
 
  write(mesg,*) "ice_shelf_solve_outer: INITIAL nonlinear residual = ",err_init
  call mom_mesg(mesg, 5)
 
  u_last(:,:) = u(:,:) ; v_last(:,:) = v(:,:)
 
  !! begin loop
 
  do iter=1,100
 
    call ice_shelf_solve_inner(cs, iss, g, u, v, taudx, taudy, h_node, float_cond, &
                               iss%hmask, conv_flag, iters, time, phi, phisub)
 
    if (cs%debug) then
      call qchksum(u, "u shelf", g%HI, haloshift=2)
      call qchksum(v, "v shelf", g%HI, haloshift=2)
    endif
 
    write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations"
    call mom_mesg(mesg, 5)
 
    call calc_shelf_visc(cs, iss, g, us, u, v)
    call pass_var(cs%ice_visc, g%domain)
    call pass_var(cs%taub_beta_eff, g%domain)
 
    ! makes sure basal stress is only applied when it is supposed to be
 
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
    enddo ; enddo
 
    u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0
 
    call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
                                        cs%taub_beta_eff, float_cond, &
                                        rhoi_rhow, u_bdry_cont, v_bdry_cont)
 
    au(:,:) = 0 ; av(:,:) = 0
 
    call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
            cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%areaT, &
            g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
 
    err_max = 0
 
      if (cs%nonlin_solve_err_mode == 1) then
 
      do j=jsumstart,g%jecB
        do i=isumstart,g%iecB
          if (cs%umask(i,j) == 1) then
            err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
          endif
          if (cs%vmask(i,j) == 1) then
            err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
          endif
          if (err_tempv >= err_max) then
            err_max = err_tempv
          endif
        enddo
      enddo
 
      call max_across_pes(err_max)
 
    elseif (cs%nonlin_solve_err_mode == 2) then
 
      max_vel = 0 ; tempu = 0 ; tempv = 0
 
      do j=jsumstart,g%jecB
        do i=isumstart,g%iecB
          if (cs%umask(i,j) == 1) then
            err_tempu = abs(u_last(i,j)-u(i,j))
            tempu = u(i,j)
          endif
          if (cs%vmask(i,j) == 1) then
            err_tempv = max(abs(v_last(i,j)- v(i,j)), err_tempu)
            tempv = sqrt(v(i,j)**2+tempu**2)
          endif
          if (err_tempv >= err_max) then
            err_max = err_tempv
          endif
          if  (tempv >= max_vel) then
            max_vel = tempv
          endif
        enddo
      enddo
 
      u_last(:,:) = u(:,:)
      v_last(:,:) = v(:,:)
 
      call max_across_pes(max_vel)
      call max_across_pes(err_max)
      err_init = max_vel
 
    endif
 
    write(mesg,*) "ice_shelf_solve_outer: nonlinear residual = ",err_max/err_init
    call mom_mesg(mesg, 5)
 
    if (err_max <= cs%nonlinear_tolerance * err_init) then
      write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations"
      call mom_mesg(mesg, 5)
      exit
    endif
 
  enddo
 
  deallocate(phi)
  deallocate(phisub)
 
end subroutine ice_shelf_solve_outer
 
subroutine ice_shelf_solve_inner(CS, ISS, G, u, v, taudx, taudy, H_node, float_cond, &
                                 hmask, conv_flag, iters, time, Phi, Phisub)
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state),  intent(in)    :: ISS !< A structure with elements that describe
                                           !! the ice-shelf state
  type(ocean_grid_type),  intent(inout) :: G  !< The grid structure used by the ice shelf.
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: u  !< The zonal ice shelf velocity at vertices [m year-1]
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: v  !< The meridional ice shelf velocity at vertices [m year-1]
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: taudx !< The x-direction driving stress, in ???
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: taudy  !< The y-direction driving stress, in ???
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: H_node !< The ice shelf thickness at nodal (corner)
                                             !! points [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: float_cond !< An array indicating where the ice
                                                !! shelf is floating: 0 if floating, 1 if not.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  integer,                intent(out)   :: conv_flag !< A flag indicating whether (1) or not (0) the
                                           !! iterations have converged to the specified tolerence
  integer,                intent(out)   :: iters !< The number of iterations used in the solver.
  type(time_type),        intent(in)    :: Time !< The current model time
  real, dimension(SZDI_(G),SZDJ_(G),8,4), &
                          intent(in)    :: Phi !< The gradients of bilinear basis elements at Gaussian
                                             !! quadrature points surrounding the cell verticies.
  real, dimension(:,:,:,:,:,:), &
                          intent(in)    :: Phisub !< Quadrature structure weights at subgridscale
                                            !! locations for finite element calculations
! one linear solve (nonlinear iteration) of the solution for velocity
 
! in this subroutine:
!    boundary contributions are added to taud to get the RHS
!    diagonal of matrix is found (for Jacobi precondition)
!    CG iteration is carried out for max. iterations or until convergence
 
! assumed - u, v, taud, visc, beta_eff are valid on the halo
 
  real, dimension(SZDIB_(G),SZDJB_(G)) ::  &
                        Ru, Rv, Zu, Zv, DIAGu, DIAGv, RHSu, RHSv, &
                        ubd, vbd, Au, Av, Du, Dv, &
                        Zu_old, Zv_old, Ru_old, Rv_old, &
                        sum_vec, sum_vec_2
  integer                 :: iter, i, j, isd, ied, jsd, jed, &
                        isc, iec, jsc, jec, is, js, ie, je, isumstart, jsumstart, &
                        isdq, iedq, jsdq, jedq, iscq, iecq, jscq, jecq, nx_halo, ny_halo
  real                     :: tol, beta_k, alpha_k, area, dot_p1, dot_p2, resid0, cg_halo, dot_p1a, dot_p2a
  character(2)                       :: gridsize
 
  real, dimension(8,4)              :: Phi_temp
  real, dimension(2,2)              :: X,Y
 
  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
  ny_halo = g%domain%njhalo ; nx_halo = g%domain%nihalo
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
 
  zu(:,:) = 0 ; zv(:,:) = 0 ; diagu(:,:) = 0 ; diagv(:,:) = 0
  ru(:,:) = 0 ; rv(:,:) = 0 ; au(:,:) = 0 ; av(:,:) = 0
  du(:,:) = 0 ; dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0
  dot_p1 = 0 ; dot_p2 = 0
 
  isumstart = g%isc
  ! Include the edge if tile is at the western bdry;  Should add a test to avoid this if reentrant.
  if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
 
  jsumstart = g%jsc
  ! Include the edge if tile is at the southern bdry;  Should add a test to avoid this if reentrant.
  if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
 
  call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
                                      cs%taub_beta_eff, float_cond, &
                                      cs%density_ice/cs%density_ocean_avg, ubd, vbd)
 
  rhsu(:,:) = taudx(:,:) - ubd(:,:)
  rhsv(:,:) = taudy(:,:) - vbd(:,:)
 
 
  call pass_vector(rhsu, rhsv, g%domain, to_all, bgrid_ne)
 
  call matrix_diagonal(cs, g, float_cond, h_node, cs%ice_visc, &
           cs%taub_beta_eff, hmask, &
           cs%density_ice/cs%density_ocean_avg, phisub, diagu, diagv)
!    DIAGu(:,:) = 1 ; DIAGv(:,:) = 1
 
  call pass_vector(diagu, diagv, g%domain, to_all, bgrid_ne)
 
  call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, hmask, &
          h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
          g%areaT, g, isc-1, iec+1, jsc-1, jec+1, cs%density_ice/cs%density_ocean_avg)
 
  call pass_vector(au, av, g%domain, to_all, bgrid_ne)
 
  ru(:,:) = rhsu(:,:) - au(:,:) ; rv(:,:) = rhsv(:,:) - av(:,:)
 
  if (.not. cs%use_reproducing_sums) then
 
    do j=jsumstart,jecq
      do i=isumstart,iecq
        if (cs%umask(i,j) == 1) dot_p1 = dot_p1 + ru(i,j)**2
        if (cs%vmask(i,j) == 1) dot_p1 = dot_p1 + rv(i,j)**2
      enddo
    enddo
 
    call sum_across_pes(dot_p1)
 
  else
 
    sum_vec(:,:) = 0.0
 
    do j=jsumstart,jecq
      do i=isumstart,iecq
        if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
        if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
      enddo
    enddo
 
    dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                       jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
 
  endif
 
  resid0 = sqrt(dot_p1)
 
  do j=jsdq,jedq
    do i=isdq,iedq
      if (cs%umask(i,j) == 1) zu(i,j) = ru(i,j) / diagu(i,j)
      if (cs%vmask(i,j) == 1) zv(i,j) = rv(i,j) / diagv(i,j)
    enddo
  enddo
 
  du(:,:) = zu(:,:) ; dv(:,:) = zv(:,:)
 
  cg_halo = 3
  conv_flag = 0
 
  !!!!!!!!!!!!!!!!!!
  !!              !!
  !! MAIN CG LOOP !!
  !!              !!
  !!!!!!!!!!!!!!!!!!
 
 
 
  ! initially, c-grid data is valid up to 3 halo nodes out
 
  do iter = 1,cs%cg_max_iterations
 
    ! assume asymmetry
    ! thus we can never assume that any arrays are legit more than 3 vertices past
    ! the computational domain - this is their state in the initial iteration
 
 
    is = isc - cg_halo ; ie = iecq + cg_halo
    js = jscq - cg_halo ; je = jecq + cg_halo
 
    au(:,:) = 0 ; av(:,:) = 0
 
    call cg_action(au, av, du, dv, phi, phisub, cs%umask, cs%vmask, hmask, &
          h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
          g%areaT, g, is, ie, js, je, cs%density_ice/cs%density_ocean_avg)
 
    ! Au, Av valid region moves in by 1
 
    if ( .not. cs%use_reproducing_sums) then
 
 
      ! alpha_k = (Z \dot R) / (D \dot AD}
      dot_p1 = 0 ; dot_p2 = 0
      do j=jsumstart,jecq
        do i=isumstart,iecq
          if (cs%umask(i,j) == 1) then
            dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
            dot_p2 = dot_p2 + du(i,j)*au(i,j)
          endif
          if (cs%vmask(i,j) == 1) then
              dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
              dot_p2 = dot_p2 + dv(i,j)*av(i,j)
          endif
        enddo
      enddo
      call sum_across_pes(dot_p1) ; call sum_across_pes(dot_p2)
    else
 
      sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
 
      do j=jscq,jecq
        do i=iscq,iecq
          if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
          if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + zv(i,j) * rv(i,j)
 
          if (cs%umask(i,j) == 1) sum_vec_2(i,j) = du(i,j) * au(i,j)
          if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + dv(i,j) * av(i,j)
        enddo
      enddo
 
      dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                         jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
 
      dot_p2 = reproducing_sum( sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                           jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
    endif
 
    alpha_k = dot_p1/dot_p2
 
    do j=jsd,jed
      do i=isd,ied
        if (cs%umask(i,j) == 1) u(i,j) = u(i,j) + alpha_k * du(i,j)
        if (cs%vmask(i,j) == 1) v(i,j) = v(i,j) + alpha_k * dv(i,j)
      enddo
    enddo
 
    do j=jsd,jed
      do i=isd,ied
        if (cs%umask(i,j) == 1) then
          ru_old(i,j) = ru(i,j) ; zu_old(i,j) = zu(i,j)
        endif
        if (cs%vmask(i,j) == 1) then
          rv_old(i,j) = rv(i,j) ; zv_old(i,j) = zv(i,j)
        endif
      enddo
    enddo
 
!    Ru(:,:) = Ru(:,:) - alpha_k * Au(:,:)
!    Rv(:,:) = Rv(:,:) - alpha_k * Av(:,:)
 
    do j=jsd,jed
      do i=isd,ied
        if (cs%umask(i,j) == 1) ru(i,j) = ru(i,j) - alpha_k * au(i,j)
        if (cs%vmask(i,j) == 1) rv(i,j) = rv(i,j) - alpha_k * av(i,j)
      enddo
    enddo
 
 
    do j=jsdq,jedq
      do i=isdq,iedq
        if (cs%umask(i,j) == 1) then
          zu(i,j) = ru(i,j) / diagu(i,j)
        endif
        if (cs%vmask(i,j) == 1) then
          zv(i,j) = rv(i,j) / diagv(i,j)
        endif
      enddo
    enddo
 
    ! R,u,v,Z valid region moves in by 1
 
    if (.not. cs%use_reproducing_sums) then
 
    ! beta_k = (Z \dot R) / (Zold \dot Rold}
      dot_p1 = 0 ; dot_p2 = 0
      do j=jsumstart,jecq
        do i=isumstart,iecq
          if (cs%umask(i,j) == 1) then
            dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
            dot_p2 = dot_p2 + zu_old(i,j)*ru_old(i,j)
          endif
          if (cs%vmask(i,j) == 1) then
            dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
            dot_p2 = dot_p2 + zv_old(i,j)*rv_old(i,j)
          endif
        enddo
      enddo
      call sum_across_pes(dot_p1) ; call sum_across_pes(dot_p2)
 
 
    else
 
      sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
 
      do j=jsumstart,jecq
        do i=isumstart,iecq
          if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
          if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + &
                                                zv(i,j) * rv(i,j)
 
          if (cs%umask(i,j) == 1) sum_vec_2(i,j) = zu_old(i,j) * ru_old(i,j)
          if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + &
                                                zv_old(i,j) * rv_old(i,j)
        enddo
      enddo
 
 
      dot_p1 = reproducing_sum(sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                        jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
 
      dot_p2 = reproducing_sum(sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                          jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
 
    endif
 
    beta_k = dot_p1/dot_p2
 
 
!    Du(:,:) = Zu(:,:) + beta_k * Du(:,:)
!    Dv(:,:) = Zv(:,:) + beta_k * Dv(:,:)
 
    do j=jsd,jed
      do i=isd,ied
        if (cs%umask(i,j) == 1) du(i,j) = zu(i,j) + beta_k * du(i,j)
        if (cs%vmask(i,j) == 1) dv(i,j) = zv(i,j) + beta_k * dv(i,j)
      enddo
    enddo
 
   ! D valid region moves in by 1
 
    dot_p1 = 0
 
    if (.not. cs%use_reproducing_sums) then
 
      do j=jsumstart,jecq
        do i=isumstart,iecq
          if (cs%umask(i,j) == 1) then
            dot_p1 = dot_p1 + ru(i,j)**2
          endif
          if (cs%vmask(i,j) == 1) then
            dot_p1 = dot_p1 + rv(i,j)**2
          endif
        enddo
      enddo
      call sum_across_pes(dot_p1)
 
    else
 
      sum_vec(:,:) = 0.0
 
      do j=jsumstart,jecq
        do i=isumstart,iecq
          if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
          if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
        enddo
      enddo
 
      dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
                                         jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
    endif
 
    dot_p1 = sqrt(dot_p1)
 
    if (dot_p1 <= cs%cg_tolerance * resid0) then
      iters = iter
      conv_flag = 1
      exit
    endif
 
    cg_halo = cg_halo - 1
 
    if (cg_halo == 0) then
      ! pass vectors
      call pass_vector(du, dv, g%domain, to_all, bgrid_ne)
      call pass_vector(u, v, g%domain, to_all, bgrid_ne)
      call pass_vector(ru, rv, g%domain, to_all, bgrid_ne)
      cg_halo = 3
    endif
 
  enddo ! end of CG loop
 
  do j=jsdq,jedq
    do i=isdq,iedq
      if (cs%umask(i,j) == 3) then
        u(i,j) = cs%u_bdry_val(i,j)
      elseif (cs%umask(i,j) == 0) then
        u(i,j) = 0
      endif
 
      if (cs%vmask(i,j) == 3) then
        v(i,j) = cs%v_bdry_val(i,j)
      elseif (cs%vmask(i,j) == 0) then
        v(i,j) = 0
      endif
    enddo
  enddo
 
  call pass_vector(u,v, g%domain, to_all, bgrid_ne)
 
  if (conv_flag == 0) then
    iters = cs%cg_max_iterations
  endif
 
end subroutine ice_shelf_solve_inner
 
subroutine ice_shelf_advect_thickness_x(CS, G, time_step, hmask, h0, h_after_uflux, flux_enter)
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  real,                   intent(in)    :: time_step !< The time step for this update [s].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: h0 !< The initial ice shelf thicknesses [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
                                              !! the zonal mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G),4), &
                          intent(inout) :: flux_enter !< The ice volume flux into the cell
                                              !! through the 4 cell boundaries [Z m2 ~> m3].
 
  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
 
  ! if there is an input bdry condition, the thickness there will be set in initialization
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !        o--- (4) ---o
  !        |           |
  !       (1)         (2)
  !        |           |
  !        o--- (3) ---o
  !
 
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
  integer :: i_off, j_off
  logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
  real, dimension(-2:2) :: stencil ! Thicknesses [Z ~> m].
  real :: u_face, &  ! positive if out
      flux_diff_cell, phi, dxh, dyh, dxdyh
  character (len=1)        :: debug_str
 
  is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  do j=jsd+1,jed-1
    if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
        ((j+j_off) >= g%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries
 
      stencil(:) = -1.
!     if (i+i_off == G%domain%nihalo+G%domain%nihalo)
      do i=is,ie
 
        if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
             ((i+i_off) >= g%domain%nihalo+1)) then
 
          if (i+i_off == g%domain%nihalo+1) then
            at_west_bdry=.true.
          else
            at_west_bdry=.false.
          endif
 
          if (i+i_off == g%domain%niglobal+g%domain%nihalo) then
            at_east_bdry=.true.
          else
            at_east_bdry=.false.
          endif
 
          if (hmask(i,j) == 1) then
 
            dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
 
            h_after_uflux(i,j) = h0(i,j)
 
            stencil(:) = h0(i-2:i+2,j)  ! fine as long has nx_halo >= 2
 
            flux_diff_cell = 0
 
            ! 1ST DO LEFT FACE
 
            if (cs%u_face_mask(i-1,j) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) / dxdyh
 
            else
 
              ! get u-velocity at center of left face
              u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
 
              if (u_face > 0) then !flux is into cell - we need info from h(i-2), h(i-1) if available
 
              ! i may not cover all the cases.. but i cover the realistic ones
 
                if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a
                              ! thickness bdry condition, and the stencil contains it
                  stencil(-1) = cs%thickness_bdry_val(i-1,j)
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
 
                elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then  ! h(i-2) and h(i-1) are valid
                  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
                           (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
 
                else                            ! h(i-1) is valid
                                    ! (o.w. flux would most likely be out of cell)
                                    !  but h(i-2) is not
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(-1)
 
                endif
 
              elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
                if (hmask(i-1,j) * hmask(i+1,j) == 1) then         ! h(i-1) and h(i+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
                             (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
 
                else
                  flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
 
                  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
                    flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
                  endif
                endif
              endif
            endif
 
            ! NEXT DO RIGHT FACE
 
            ! get u-velocity at center of right face
 
            if (cs%u_face_mask(i+1,j) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) / dxdyh
 
            else
 
              u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
 
              if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available
 
                if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a
                              ! thickness bdry condition, and the stencil contains it
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
 
                elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then  ! h(i+2) and h(i+1) are valid
 
                  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
                      (stencil(1) - phi * (stencil(1)-stencil(0))/2)
 
                else                            ! h(i+1) is valid
                                            ! (o.w. flux would most likely be out of cell)
                                            !  but h(i+2) is not
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(1)
 
                endif
 
              elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
 
                if (hmask(i-1,j) * hmask(i+1,j) == 1) then         ! h(i-1) and h(i+1) are both valid
 
                  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(1))/2)
 
                else                            ! h(i+1) is valid
                                            ! (o.w. flux would most likely be out of cell)
                                            !  but h(i+2) is not
 
                  flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
 
                  if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
                    flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step  * stencil(0)
                  endif
 
                endif
 
              endif
 
              h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
 
            endif
 
          elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
 
            if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then
              u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
              flux_enter(i,j,1) = abs(u_face) * g%dyT(i,j) * time_step * cs%thickness_bdry_val(i-1,j)
            elseif (cs%u_face_mask(i-1,j) == 4.) then
              flux_enter(i,j,1) = g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)
            endif
 
            if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then
              u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
              flux_enter(i,j,2) = abs(u_face) * g%dyT(i,j) * time_step * cs%thickness_bdry_val(i+1,j)
            elseif (cs%u_face_mask(i+1,j) == 4.) then
              flux_enter(i,j,2) = g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j)
            endif
 
            if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then
              ! this is solely for the purposes of keeping the mask consistent while advancing
              ! the front without having to call pass_var - if cell is empty and cell to left
              ! is ice-covered then this cell will become partly covered
 
              hmask(i,j) = 2
            elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then
              ! this is solely for the purposes of keeping the mask consistent while advancing
              ! the front without having to call pass_var - if cell is empty and cell to left
              ! is ice-covered then this cell will become partly covered
 
              hmask(i,j) = 2
 
            endif
 
          endif
 
        endif
 
      enddo ! i loop
 
    endif
 
  enddo ! j loop
 
end subroutine ice_shelf_advect_thickness_x
 
subroutine ice_shelf_advect_thickness_y(CS, G, time_step, hmask, h_after_uflux, h_after_vflux, flux_enter)
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  real,                   intent(in)    :: time_step !< The time step for this update [s].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: h_after_uflux !< The ice shelf thicknesses after
                                              !! the zonal mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
                                              !! the meridional mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G),4), &
                          intent(inout) :: flux_enter !< The ice volume flux into the cell
                                              !! through the 4 cell boundaries [Z m2 ~> m3].
 
  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
 
  ! if there is an input bdry condition, the thickness there will be set in initialization
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !        o--- (4) ---o
  !        |           |
  !       (1)         (2)
  !        |           |
  !        o--- (3) ---o
  !
 
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
  integer :: i_off, j_off
  logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
  real, dimension(-2:2) :: stencil  ! Thicknesses [Z ~> m].
  real :: v_face, &  ! positive if out
      flux_diff_cell, phi, dxh, dyh, dxdyh
  character(len=1)        :: debug_str
 
  is = g%isc ; ie = g%iec ; js = g%jsc-1 ; je = g%jec+1 ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  do i=isd+2,ied-2
    if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
       ((i+i_off) >= g%domain%nihalo+1)) then  ! based on Mehmet's code - only if btw east & west boundaries
 
      stencil(:) = -1
 
      do j=js,je
 
        if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
             ((j+j_off) >= g%domain%njhalo+1)) then
 
          if (j+j_off == g%domain%njhalo+1) then
            at_south_bdry=.true.
          else
            at_south_bdry=.false.
          endif
 
          if (j+j_off == g%domain%njglobal+g%domain%njhalo) then
            at_north_bdry=.true.
          else
            at_north_bdry=.false.
          endif
 
          if (hmask(i,j) == 1) then
            dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
            h_after_vflux(i,j) = h_after_uflux(i,j)
 
            stencil(:) = h_after_uflux(i,j-2:j+2)  ! fine as long has ny_halo >= 2
            flux_diff_cell = 0
 
            ! 1ST DO south FACE
 
            if (cs%v_face_mask(i,j-1) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) / dxdyh
 
            else
 
              ! get u-velocity at center of left face
              v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
 
              if (v_face > 0) then !flux is into cell - we need info from h(j-2), h(j-1) if available
 
                ! i may not cover all the cases.. but i cover the realistic ones
 
                if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a
                                            ! thickness bdry condition, and the stencil contains it
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
 
                elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then  ! h(j-2) and h(j-1) are valid
 
                  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
 
                else     ! h(j-1) is valid
                         ! (o.w. flux would most likely be out of cell)
                         !  but h(j-2) is not
                  flux_diff_cell = flux_diff_cell + abs(v_face) * (dxh * time_step / dxdyh) * stencil(-1)
                endif
 
              elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
 
                if (hmask(i,j-1) * hmask(i,j+1) == 1) then  ! h(j-1) and h(j+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
                else
                  flux_diff_cell = flux_diff_cell - abs(v_face) * (dxh * time_step / dxdyh) * stencil(0)
 
                  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
                    flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
                  endif
 
                endif
 
              endif
 
            endif
 
            ! NEXT DO north FACE
 
            if (cs%v_face_mask(i,j+1) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) / dxdyh
 
            else
 
            ! get u-velocity at center of right face
              v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
 
              if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available
 
                if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a
                                            ! thickness bdry condition, and the stencil contains it
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
                elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then  ! h(j+2) and h(j+1) are valid
                  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(1) - phi * (stencil(1)-stencil(0))/2)
                else     ! h(j+1) is valid
                         ! (o.w. flux would most likely be out of cell)
                         !  but h(j+2) is not
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
                endif
 
              elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
 
                if (hmask(i,j-1) * hmask(i,j+1) == 1) then         ! h(j-1) and h(j+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(1))/2)
                else   ! h(j+1) is valid
                       ! (o.w. flux would most likely be out of cell)
                       !  but h(j+2) is not
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
                  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
                    flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
                  endif
                endif
 
              endif
 
            endif
 
            h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
 
          elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
 
            if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then
              v_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i,j-1))
              flux_enter(i,j,3) = abs(v_face) * g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j-1)
            elseif (cs%v_face_mask(i,j-1) == 4.) then
              flux_enter(i,j,3) = g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)
            endif
 
            if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then
              v_face = 0.5 * (cs%u_shelf(i-1,j) + cs%u_shelf(i,j))
              flux_enter(i,j,4) = abs(v_face) * g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j+1)
            elseif (cs%v_face_mask(i,j+1) == 4.) then
              flux_enter(i,j,4) = g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)
            endif
 
            if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then
                ! this is solely for the purposes of keeping the mask consistent while advancing
                ! the front without having to call pass_var - if cell is empty and cell to left
                ! is ice-covered then this cell will become partly covered
              hmask(i,j) = 2
            elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then
                ! this is solely for the purposes of keeping the mask consistent while advancing
                ! the front without having to call pass_var - if cell is empty and cell to left
                ! is ice-covered then this cell will become partly covered
              hmask(i,j) = 2
            endif
 
          endif
        endif
      enddo ! j loop
    endif
  enddo ! i loop
 
end subroutine ice_shelf_advect_thickness_y
 
subroutine shelf_advance_front(CS, ISS, G, flux_enter)
  type(ice_shelf_dyn_cs), intent(in)    :: cs !< A pointer to the ice shelf control structure
  type(ice_shelf_state),  intent(inout) :: iss !< A structure with elements that describe
                                           !! the ice-shelf state
  type(ocean_grid_type),  intent(in)    :: g  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G),4), &
                          intent(inout) :: flux_enter !< The ice volume flux into the cell
                                              !! through the 4 cell boundaries [Z m2 ~> m3].
 
  ! in this subroutine we go through the computational cells only and, if they are empty or partial cells,
  ! we find the reference thickness and update the shelf mass and partial area fraction and the hmask if necessary
 
  ! if any cells go from partial to complete, we then must set the thickness, update hmask accordingly,
  ! and divide the overflow across the adjacent EMPTY (not partly-covered) cells.
  ! (it is highly unlikely there will not be any; in which case this will need to be rethought.)
 
  ! most likely there will only be one "overflow". if not, though, a pass_var of all relevant variables
  ! is done; there will therefore be a loop which, in practice, will hopefully not have to go through
  ! many iterations
 
  ! when 3d advected scalars are introduced, they will be impacted by what is done here
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !        o--- (4) ---o
  !        |           |
  !       (1)         (2)
  !        |           |
  !        o--- (3) ---o
  !
 
  integer :: i, j, isc, iec, jsc, jec, n_flux, k, l, iter_count
  integer :: i_off, j_off
  integer :: iter_flag
 
  real :: h_reference, dxh, dyh, dxdyh, rho, partial_vol, tot_flux
  character(len=160) :: mesg  ! The text of an error message
  integer, dimension(4) :: mapi, mapj, new_partial
!   real, dimension(size(flux_enter,1),size(flux_enter,2),size(flux_enter,2)) :: flux_enter_replace
  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter_replace
 
  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
  i_off = g%idg_offset ; j_off = g%jdg_offset
  rho = cs%density_ice
  iter_count = 0 ; iter_flag = 1
 
 
  mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0
  mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0
 
  do while (iter_flag == 1)
 
    iter_flag = 0
 
    if (iter_count > 0) then
      flux_enter(:,:,:) = flux_enter_replace(:,:,:)
    endif
    flux_enter_replace(:,:,:) = 0.0
 
    iter_count = iter_count + 1
 
    ! if iter_count >= 3 then some halo updates need to be done...
 
    do j=jsc-1,jec+1
 
      if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
         ((j+j_off) >= g%domain%njhalo+1)) then
 
        do i=isc-1,iec+1
 
          if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
              ((i+i_off) >= g%domain%nihalo+1)) then
        ! first get reference thickness by averaging over cells that are fluxing into this cell
            n_flux = 0
            h_reference = 0.0
            tot_flux = 0.0
 
            do k=1,2
              if (flux_enter(i,j,k) > 0) then
                n_flux = n_flux + 1
                h_reference = h_reference + iss%h_shelf(i+2*k-3,j)
                tot_flux = tot_flux + flux_enter(i,j,k)
                flux_enter(i,j,k) = 0.0
              endif
            enddo
 
            do k=1,2
              if (flux_enter(i,j,k+2) > 0) then
                n_flux = n_flux + 1
                h_reference = h_reference + iss%h_shelf(i,j+2*k-3)
                tot_flux = tot_flux + flux_enter(i,j,k+2)
                flux_enter(i,j,k+2) = 0.0
              endif
            enddo
 
            if (n_flux > 0) then
              dxdyh = g%areaT(i,j)
              h_reference = h_reference / real(n_flux)
              partial_vol = iss%h_shelf(i,j) * iss%area_shelf_h(i,j) + tot_flux
 
              if ((partial_vol / dxdyh) == h_reference) then ! cell is exactly covered, no overflow
                iss%hmask(i,j) = 1
                iss%h_shelf(i,j) = h_reference
                iss%area_shelf_h(i,j) = dxdyh
              elseif ((partial_vol / dxdyh) < h_reference) then
                iss%hmask(i,j) = 2
        !         ISS%mass_shelf(i,j) = partial_vol * rho
                iss%area_shelf_h(i,j) = partial_vol / h_reference
                iss%h_shelf(i,j) = h_reference
              else
 
                iss%hmask(i,j) = 1
                iss%area_shelf_h(i,j) = dxdyh
                !h_temp(i,j) = h_reference
                partial_vol = partial_vol - h_reference * dxdyh
 
                iter_flag  = 1
 
                n_flux = 0 ; new_partial(:) = 0
 
                do k=1,2
                  if (cs%u_face_mask(i-2+k,j) == 2) then
                    n_flux = n_flux + 1
                  elseif (iss%hmask(i+2*k-3,j) == 0) then
                    n_flux = n_flux + 1
                    new_partial(k) = 1
                  endif
                enddo
                do k=1,2
                  if (cs%v_face_mask(i,j-2+k) == 2) then
                    n_flux = n_flux + 1
                  elseif (iss%hmask(i,j+2*k-3) == 0) then
                    n_flux = n_flux + 1
                    new_partial(k+2) = 1
                  endif
                enddo
 
                if (n_flux == 0) then ! there is nowhere to put the extra ice!
                  iss%h_shelf(i,j) = h_reference + partial_vol / dxdyh
                else
                  iss%h_shelf(i,j) = h_reference
 
                  do k=1,2
                    if (new_partial(k) == 1) &
                      flux_enter_replace(i+2*k-3,j,3-k) = partial_vol / real(n_flux)
                  enddo
                  do k=1,2 ! ### Combine these two loops?
                    if (new_partial(k+2) == 1) &
                      flux_enter_replace(i,j+2*k-3,5-k) = partial_vol / real(n_flux)
                  enddo
                endif
 
              endif ! Parital_vol test.
            endif ! n_flux gt 0 test.
 
          endif
        enddo ! j-loop
      endif
    enddo
 
  !  call max_across_PEs(iter_flag)
 
  enddo ! End of do while(iter_flag) loop
 
  call max_across_pes(iter_count)
 
  if (is_root_pe() .and. (iter_count > 1)) then
    write(mesg,*) "shelf_advance_front: ", iter_count, " max iterations"
    call mom_mesg(mesg, 5)
  endif
 
end subroutine shelf_advance_front
 
!> Apply a very simple calving law using a minimum thickness rule
subroutine ice_shelf_min_thickness_calve(G, h_shelf, area_shelf_h, hmask, thickness_calve)
  type(ocean_grid_type), intent(in)    :: g  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [m2].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real,                  intent(in)    :: thickness_calve !< The thickness at which to trigger calving [Z ~> m].
 
  integer                        :: i,j
 
  do j=g%jsd,g%jed
    do i=g%isd,g%ied
!      if ((h_shelf(i,j) < CS%thickness_calve) .and. (hmask(i,j) == 1) .and. &
!           (CS%float_frac(i,j) == 0.0)) then
      if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.)) then
        h_shelf(i,j) = 0.0
        area_shelf_h(i,j) = 0.0
        hmask(i,j) = 0.0
      endif
    enddo
  enddo
 
end subroutine ice_shelf_min_thickness_calve
 
subroutine calve_to_mask(G, h_shelf, area_shelf_h, hmask, calve_mask)
  type(ocean_grid_type), intent(in) :: g  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [m2].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(inout) :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: calve_mask !< A mask that indicates where the ice shelf
                                             !! can exist, and where it will calve.
 
  integer                        :: i,j
 
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if ((calve_mask(i,j) == 0.0) .and. (hmask(i,j) /= 0.0)) then
      h_shelf(i,j) = 0.0
      area_shelf_h(i,j) = 0.0
      hmask(i,j) = 0.0
    endif
  enddo ; enddo
 
end subroutine calve_to_mask
 
subroutine calc_shelf_driving_stress(CS, ISS, G, US, TAUD_X, TAUD_Y, OD)
  type(ice_shelf_dyn_cs), intent(in)   :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(in)    :: ISS !< A structure with elements that describe
                                             !! the ice-shelf state
  type(ocean_grid_type), intent(inout) :: G  !< The grid structure used by the ice shelf.
  type(unit_scale_type), intent(in)    :: US !< Pointer to a structure containing unit conversion factors
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: OD  !< ocean floor depth at tracer points [Z ~> m].
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(inout) :: TAUD_X  !< X-direction driving stress at q-points
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(inout) :: TAUD_Y  !< Y-direction driving stress at q-points
 
! driving stress!
 
! ! TAUD_X and TAUD_Y will hold driving stress in the x- and y- directions when done.
!    they will sit on the BGrid, and so their size depends on whether the grid is symmetric
!
! Since this is a finite element solve, they will actually have the form \int \phi_i rho g h \nabla s
!
! OD -this is important and we do not yet know where (in MOM) it will come from. It represents
!     "average" ocean depth -- and is needed to find surface elevation
!    (it is assumed that base_ice = bed + OD)
 
  real, dimension(SIZE(OD,1),SIZE(OD,2))  :: S, &     ! surface elevation [Z ~> m].
                            BASE     ! basal elevation of shelf/stream [Z ~> m].
 
 
  real      :: rho, rhow, sx, sy, neumann_val, dxh, dyh, dxdyh, grav
 
  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
  integer :: i_off, j_off
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
  isd = g%isd ; jsd = g%jsd
  iegq = g%iegB ; jegq = g%jegB
  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
  giec = g%domain%niglobal+g%domain%nihalo ; gjec = g%domain%njglobal+g%domain%njhalo
  is = iscq - 1; js = jscq - 1
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  rho = cs%density_ice
  rhow = cs%density_ocean_avg
  grav = us%Z_to_m**2 * cs%g_Earth
 
  ! prelim - go through and calculate S
 
  ! or is this faster?
  base(:,:) = -g%bathyT(:,:) + od(:,:)
  s(:,:) = base(:,:) + iss%h_shelf(:,:)
 
  do j=jsc-1,jec+1
    do i=isc-1,iec+1
      cnt = 0
      sx = 0
      sy = 0
      dxh = g%dxT(i,j)
      dyh = g%dyT(i,j)
      dxdyh = g%areaT(i,j)
 
      if (iss%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell
 
        ! calculate sx
        if ((i+i_off) == gisc) then ! at left computational bdry
          if (iss%hmask(i+1,j) == 1) then
            sx = (s(i+1,j)-s(i,j))/dxh
          else
            sx = 0
          endif
        elseif ((i+i_off) == giec) then ! at right computational bdry
          if (iss%hmask(i-1,j) == 1) then
            sx = (s(i,j)-s(i-1,j))/dxh
          else
            sx = 0
          endif
        else ! interior
          if (iss%hmask(i+1,j) == 1) then
            cnt = cnt+1
            sx = s(i+1,j)
          else
            sx = s(i,j)
          endif
          if (iss%hmask(i-1,j) == 1) then
            cnt = cnt+1
            sx = sx - s(i-1,j)
          else
            sx = sx - s(i,j)
          endif
          if (cnt == 0) then
            sx = 0
          else
            sx = sx / (cnt * dxh)
          endif
        endif
 
        cnt = 0
 
        ! calculate sy, similarly
        if ((j+j_off) == gjsc) then ! at south computational bdry
          if (iss%hmask(i,j+1) == 1) then
            sy = (s(i,j+1)-s(i,j))/dyh
          else
            sy = 0
          endif
        elseif ((j+j_off) == gjec) then ! at nprth computational bdry
          if (iss%hmask(i,j-1) == 1) then
            sy = (s(i,j)-s(i,j-1))/dyh
          else
            sy = 0
          endif
        else ! interior
          if (iss%hmask(i,j+1) == 1) then
            cnt = cnt+1
            sy = s(i,j+1)
          else
            sy = s(i,j)
          endif
          if (iss%hmask(i,j-1) == 1) then
            cnt = cnt+1
            sy = sy - s(i,j-1)
          else
            sy = sy - s(i,j)
          endif
          if (cnt == 0) then
            sy = 0
          else
            sy = sy / (cnt * dyh)
          endif
        endif
 
        ! SW vertex
        taud_x(i-1,j-1) = taud_x(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
        taud_y(i-1,j-1) = taud_y(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
 
        ! SE vertex
        taud_x(i,j-1) = taud_x(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
        taud_y(i,j-1) = taud_y(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
 
        ! NW vertex
        taud_x(i-1,j) = taud_x(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
        taud_y(i-1,j) = taud_y(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
 
        ! NE vertex
        taud_x(i,j) = taud_x(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
        taud_y(i,j) = taud_y(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
 
        if (cs%float_frac(i,j) == 1) then
          neumann_val = .5 * grav * (rho * iss%h_shelf(i,j)**2 - rhow * g%bathyT(i,j)**2)
        else
          neumann_val = .5 * grav * (1-rho/rhow) * rho * iss%h_shelf(i,j)**2
        endif
 
 
        if ((cs%u_face_mask(i-1,j) == 2) .OR. (iss%hmask(i-1,j) == 0) .OR. (iss%hmask(i-1,j) == 2) ) then
          ! left face of the cell is at a stress boundary
          ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated
          ! pressure on either side of the face
          ! on the ice side, it is rho g h^2 / 2
          ! on the ocean side, it is rhow g (delta OD)^2 / 2
          ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation
          !     is not above the base of the ice in the current cell
 
          ! note negative sign due to direction of normal vector
          taud_x(i-1,j-1) = taud_x(i-1,j-1) - .5 * dyh * neumann_val
          taud_x(i-1,j) = taud_x(i-1,j) - .5 * dyh * neumann_val
        endif
 
        if ((cs%u_face_mask(i,j) == 2) .OR. (iss%hmask(i+1,j) == 0) .OR. (iss%hmask(i+1,j) == 2) ) then
          ! right face of the cell is at a stress boundary
          taud_x(i,j-1) = taud_x(i,j-1) + .5 * dyh * neumann_val
          taud_x(i,j) = taud_x(i,j) + .5 * dyh * neumann_val
        endif
 
        if ((cs%v_face_mask(i,j-1) == 2) .OR. (iss%hmask(i,j-1) == 0) .OR. (iss%hmask(i,j-1) == 2) ) then
          ! south face of the cell is at a stress boundary
          taud_y(i-1,j-1) = taud_y(i-1,j-1) - .5 * dxh * neumann_val
          taud_y(i,j-1) = taud_y(i,j-1) - .5 * dxh * neumann_val
        endif
 
        if ((cs%v_face_mask(i,j) == 2) .OR. (iss%hmask(i,j+1) == 0) .OR. (iss%hmask(i,j+1) == 2) ) then
          ! north face of the cell is at a stress boundary
          taud_y(i-1,j) = taud_y(i-1,j) + .5 * dxh * neumann_val ! note negative sign due to direction of normal vector
          taud_y(i,j) = taud_y(i,j) + .5 * dxh * neumann_val
        endif
 
      endif
    enddo
  enddo
 
end subroutine calc_shelf_driving_stress
 
subroutine init_boundary_values(CS, G, time, hmask, input_flux, input_thick, new_sim)
  type(ice_shelf_dyn_cs),intent(inout) :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type), intent(inout) :: G  !< The grid structure used by the ice shelf.
  type(time_type),       intent(in)    :: Time !< The current model time
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real,                  intent(in)    :: input_flux !< The integrated inward ice thickness flux [Z m2 s-1 ~> m3 s-1]
  real,                  intent(in)    :: input_thick !< The ice thickness at boundaries [Z ~> m].
  logical,     optional, intent(in)    :: new_sim !< If present and false, this run is being restarted
 
! this will be a per-setup function. the boundary values of thickness and velocity
! (and possibly other variables) will be updated in this function
 
! FOR RESTARTING PURPOSES: if grid is not symmetric and the model is restarted, we will
!               need to update those velocity points not *technically* in any
!               computational domain -- if this function gets moves to another module,
!               DO NOT TAKE THE RESTARTING BIT WITH IT
  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
  integer :: gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
  integer :: i_off, j_off
  real :: A, n, ux, uy, vx, vy, eps_min, domain_width
 
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
!   iscq = G%iscq ; iecq = G%iecq ; jscq = G%jscq ; jecq = G%jecq
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
!   iegq = G%iegq ; jegq = G%jegq
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  domain_width = g%len_lat
 
  ! this loop results in some values being set twice but... eh.
 
  do j=jsd,jed
    do i=isd,ied
 
      if (hmask(i,j) == 3) then
        cs%thickness_bdry_val(i,j) = input_thick
      endif
 
      if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2)) then
        if ((i <= iec).and.(i >= isc)) then
          if (cs%u_face_mask(i-1,j) == 3) then
            cs%u_bdry_val(i-1,j-1) = (1 - ((g%geoLatBu(i-1,j-1) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
                  1.5 * input_flux / input_thick
            cs%u_bdry_val(i-1,j) = (1 - ((g%geoLatBu(i-1,j) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
                  1.5 * input_flux / input_thick
          endif
        endif
      endif
 
      if (.not.(new_sim)) then
        if (.not. g%symmetric) then
          if (((i+i_off) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3)) then
            cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
            cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
          endif
          if (((j+j_off) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3)) then
            cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
            cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
          endif
        endif
      endif
    enddo
  enddo
 
end subroutine init_boundary_values
 
 
subroutine cg_action(uret, vret, u, v, Phi, Phisub, umask, vmask, hmask, H_node, &
                nu, float_cond, bathyT, beta, dxdyh, G, is, ie, js, je, dens_ratio)
 
  type(ocean_grid_type), intent(in) :: G  !< The grid structure used by the ice shelf.
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                         intent(inout) :: uret !< The retarding stresses working at u-points.
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                         intent(inout) :: vret !< The retarding stresses working at v-points.
  real, dimension(SZDI_(G),SZDJ_(G),8,4), &
                         intent(in)   :: Phi !< The gradients of bilinear basis elements at Gaussian
                                             !! quadrature points surrounding the cell verticies.
  real, dimension(:,:,:,:,:,:), &
                         intent(in)    :: Phisub !< Quadrature structure weights at subgridscale
                                            !! locations for finite element calculations
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: u  !< The zonal ice shelf velocity at vertices [m year-1]
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: v  !< The meridional ice shelf velocity at vertices [m year-1]
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: umask !< A coded mask indicating the nature of the
                                             !! zonal flow at the corner point
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: vmask !< A coded mask indicating the nature of the
                                             !! meridional flow at the corner point
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: H_node !< The ice shelf thickness at nodal (corner)
                                             !! points [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: nu   !< A field related to the ice viscosity from Glen's
                                               !! flow law. The exact form and units depend on the
                                               !! basal law exponent.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: float_cond !< An array indicating where the ice
                                                !! shelf is floating: 0 if floating, 1 if not.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(in)    :: beta  !< A field related to the nonlinear part of the
                                                !! "linearized" basal stress.  The exact form and
                                                !! units depend on the basal law exponent.
                ! and/or whether flow is "hybridized"
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: dxdyh  !< The tracer cell area [m2]
  real,                  intent(in)    :: dens_ratio !< The density of ice divided by the density
                                                     !! of seawater, nondimensional
  integer,               intent(in)    :: is  !< The starting i-index to work on
  integer,               intent(in)    :: ie  !< The ending i-index to work on
  integer,               intent(in)    :: js  !< The starting j-index to work on
  integer,               intent(in)    :: je  !< The ending j-index to work on
 
! the linear action of the matrix on (u,v) with bilinear finite elements
! as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
! but this may change pursuant to conversations with others
!
! is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
!     in order to make less frequent halo updates
 
! the linear action of the matrix on (u,v) with bilinear finite elements
! Phi has the form
! Phi(i,j,k,q) - applies to cell i,j
 
    !  3 - 4
    !  |   |
    !  1 - 2
 
! Phi(i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
! Phi(i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
! Phi_k is equal to 1 at vertex k, and 0 at vertex l /= k, and bilinear
 
  real :: ux, vx, uy, vy, uq, vq, area, basel
  integer :: iq, jq, iphi, jphi, i, j, ilq, jlq
  real, dimension(2) :: xquad
  real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr,Ucontr
 
  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
 
  do j=js,je
    do i=is,ie ; if (hmask(i,j) == 1) then
!     dxh = G%dxh(i,j)
!     dyh = G%dyh(i,j)
!
!     X(:,:) = G%geoLonBu(i-1:i,j-1:j)
!     Y(:,:) = G%geoLatBu(i-1:i,j-1:j)
!
!     call bilinear_shape_functions (X, Y, Phi, area)
 
    ! X and Y must be passed in the form
        !  3 - 4
        !  |   |
        !  1 - 2
    ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
    ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
 
      area = dxdyh(i,j)
 
        ucontr=0
      do iq=1,2 ; do jq=1,2
 
 
        if (iq == 2) then
            ilq = 2
        else
            ilq = 1
        endif
 
        if (jq == 2) then
            jlq = 2
        else
            jlq = 1
        endif
 
        uq = u(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
        u(i,j-1) * xquad(iq) * xquad(3-jq) + &
        u(i-1,j) * xquad(3-iq) * xquad(jq) + &
        u(i,j) * xquad(iq) * xquad(jq)
 
        vq = v(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
        v(i,j-1) * xquad(iq) * xquad(3-jq) + &
        v(i-1,j) * xquad(3-iq) * xquad(jq) + &
        v(i,j) * xquad(iq) * xquad(jq)
 
        ux = u(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
        u(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
        u(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
        u(i,j) * phi(i,j,7,2*(jq-1)+iq)
 
        vx = v(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
        v(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
        v(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
        v(i,j) * phi(i,j,7,2*(jq-1)+iq)
 
        uy = u(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
        u(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
        u(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
        u(i,j) * phi(i,j,8,2*(jq-1)+iq)
 
        vy = v(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
        v(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
        v(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
        v(i,j) * phi(i,j,8,2*(jq-1)+iq)
 
        do iphi=1,2 ; do jphi=1,2
          if (umask(i-2+iphi,j-2+jphi) == 1) then
 
            uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
                .25 * area * nu(i,j) * ((4*ux+2*vy) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                                (uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
          endif
          if (vmask(i-2+iphi,j-2+jphi) == 1) then
 
            vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
                .25 * area * nu(i,j) * ((uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                                (4*vy+2*ux) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
          endif
 
          if (iq == iphi) then
            ilq = 2
          else
            ilq = 1
          endif
 
          if (jq == jphi) then
            jlq = 2
          else
            jlq = 1
          endif
 
          if (float_cond(i,j) == 0) then
 
            if (umask(i-2+iphi,j-2+jphi) == 1) then
 
              uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * area * uq * xquad(ilq) * xquad(jlq)
 
            endif
 
            if (vmask(i-2+iphi,j-2+jphi) == 1) then
 
              vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * area * vq * xquad(ilq) * xquad(jlq)
 
            endif
 
          endif
              ucontr(iphi,jphi) = ucontr(iphi,jphi) + .25 * area * uq * xquad(ilq) * xquad(jlq) * beta(i,j)
        enddo ; enddo
      enddo ; enddo
 
      if (float_cond(i,j) == 1) then
        usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = bathyt(i,j)
        ucell(:,:) = u(i-1:i,j-1:j) ; vcell(:,:) = v(i-1:i,j-1:j) ; hcell(:,:) = h_node(i-1:i,j-1:j)
        call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, area, basel, &
                                     dens_ratio, usubcontr, vsubcontr)
        do iphi=1,2 ; do jphi=1,2
          if (umask(i-2+iphi,j-2+jphi) == 1) then
            uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
          endif
          if (vmask(i-2+iphi,j-2+jphi) == 1) then
            vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
          endif
        enddo ; enddo
      endif
 
    endif
  enddo ; enddo
 
end subroutine cg_action
 
subroutine cg_action_subgrid_basal(Phisub, H, U, V, DXDYH, bathyT, dens_ratio, Ucontr, Vcontr)
  real, dimension(:,:,:,:,:,:), &
                        intent(in)    :: Phisub !< Quadrature structure weights at subgridscale
                                            !! locations for finite element calculations
  real, dimension(2,2), intent(in)    :: H  !< The ice shelf thickness at nodal (corner) points [Z ~> m].
  real, dimension(2,2), intent(in)    :: U  !< The zonal ice shelf velocity at vertices [m year-1]
  real, dimension(2,2), intent(in)    :: V  !< The meridional ice shelf velocity at vertices [m year-1]
  real,                 intent(in)    :: DXDYH  !< The tracer cell area [m2]
  real,                 intent(in)    :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
  real,                 intent(in)    :: dens_ratio !< The density of ice divided by the density
                                                    !! of seawater, nondimensional
  real, dimension(2,2), intent(inout) :: Ucontr !< A field related to the subgridscale contributions to
                                                !! the u-direction basal stress.
  real, dimension(2,2), intent(inout) :: Vcontr !< A field related to the subgridscale contributions to
                                                !! the v-direction basal stress.
 
  integer              :: nsub, i, j, k, l, qx, qy, m, n
  real                 :: subarea, hloc, uq, vq
 
  nsub = size(phisub,1)
  subarea = dxdyh / (nsub**2)
 
  do m=1,2
    do n=1,2
      do j=1,nsub
        do i=1,nsub
          do qx=1,2
            do qy = 1,2
 
              hloc = phisub(i,j,1,1,qx,qy)*h(1,1) + phisub(i,j,1,2,qx,qy)*h(1,2) + &
                     phisub(i,j,2,1,qx,qy)*h(2,1) + phisub(i,j,2,2,qx,qy)*h(2,2)
 
              if (dens_ratio * hloc - bathyt > 0) then
              !if (.true.) then
                uq = 0 ; vq = 0
                do k=1,2
                  do l=1,2
                    !Ucontr(m,n) = Ucontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * Phisub(i,j,k,l,qx,qy) * U(k,l)
                    !Vcontr(m,n) = Vcontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * Phisub(i,j,k,l,qx,qy) * V(k,l)
                    uq = uq + phisub(i,j,k,l,qx,qy) * u(k,l) ; vq = vq + phisub(i,j,k,l,qx,qy) * v(k,l)
                  enddo
                enddo
 
                ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * uq
                vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * vq
 
              endif
 
            enddo
          enddo
        enddo
      enddo
    enddo
  enddo
 
end subroutine cg_action_subgrid_basal
 
!> returns the diagonal entries of the matrix for a Jacobi preconditioning
subroutine matrix_diagonal(CS, G, float_cond, H_node, nu, beta, hmask, dens_ratio, &
                           Phisub, u_diagonal, v_diagonal)
 
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: float_cond !< An array indicating where the ice
                                                !! shelf is floating: 0 if floating, 1 if not.
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: H_node !< The ice shelf thickness at nodal
                                                 !! (corner) points [Z ~> m].
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: nu   !< A field related to the ice viscosity from Glen's
                                                !! flow law. The exact form and units depend on the
                                                !! basal law exponent.
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: beta !< A field related to the nonlinear part of the
                                                !! "linearized" basal stress.  The exact form and
                                                !! units depend on the basal law exponent
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real,                   intent(in)    :: dens_ratio !< The density of ice divided by the density
                                                     !! of seawater, nondimensional
  real, dimension(:,:,:,:,:,:), intent(in) :: Phisub !< Quadrature structure weights at subgridscale
                                            !! locations for finite element calculations
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: u_diagonal !< The diagonal elements of the u-velocity
                                            !! matrix from the left-hand side of the solver.
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: v_diagonal  !< The diagonal elements of the v-velocity
                                            !! matrix from the left-hand side of the solver.
 
 
! returns the diagonal entries of the matrix for a Jacobi preconditioning
 
  integer :: i, j, is, js, cnt, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq
  real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, area, uq, vq, basel
  real, dimension(8,4)  :: Phi
  real, dimension(4) :: X, Y
  real, dimension(2) :: xquad
  real, dimension(2,2) :: Hcell,Usubcontr,Vsubcontr
 
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
 
  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
 
! X and Y must be passed in the form
    !  3 - 4
    !  |   |
    !  1 - 2
! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
 
  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1) then
 
    dxh = g%dxT(i,j)
    dyh = g%dyT(i,j)
    dxdyh = g%areaT(i,j)
 
    x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
    x(3:4) = g%geoLonBu(i-1:i,j) *1000
    y(1:2) = g%geoLatBu(i-1:i,j-1) *1000
    y(3:4) = g%geoLatBu(i-1:i,j)*1000
 
    call bilinear_shape_functions(x, y, phi, area)
 
    ! X and Y must be passed in the form
        !  3 - 4
        !  |   |
        !  1 - 2
    ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
    ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
 
    do iq=1,2 ; do jq=1,2
 
      do iphi=1,2 ; do jphi=1,2
 
        if (iq == iphi) then
          ilq = 2
        else
          ilq = 1
        endif
 
        if (jq == jphi) then
          jlq = 2
        else
          jlq = 1
        endif
 
        if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
 
          ux = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
          uy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
          vx = 0.
          vy = 0.
 
          u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
              .25 * dxdyh * nu(i,j) * ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                              (uy+vy) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
 
          uq = xquad(ilq) * xquad(jlq)
 
          if (float_cond(i,j) == 0) then
            u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
          endif
 
        endif
 
        if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
 
          vx = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
          vy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
          ux = 0.
          uy = 0.
 
          v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
              .25 * dxdyh * nu(i,j) * ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                              (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
 
          vq = xquad(ilq) * xquad(jlq)
 
          if (float_cond(i,j) == 0) then
            v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
          endif
 
        endif
      enddo ; enddo
    enddo ; enddo
    if (float_cond(i,j) == 1) then
      usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
      hcell(:,:) = h_node(i-1:i,j-1:j)
      call cg_diagonal_subgrid_basal(phisub, hcell, dxdyh, basel, dens_ratio, usubcontr, vsubcontr)
      do iphi=1,2 ; do jphi=1,2
        if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
          u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
          v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
        endif
      enddo ; enddo
    endif
  endif ; enddo ; enddo
 
end subroutine matrix_diagonal
 
subroutine cg_diagonal_subgrid_basal (Phisub, H_node, DXDYH, bathyT, dens_ratio, Ucontr, Vcontr)
  real, dimension(:,:,:,:,:,:), &
                        intent(in) :: Phisub !< Quadrature structure weights at subgridscale
                                             !! locations for finite element calculations
  real, dimension(2,2), intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
                                             !! points [Z ~> m].
  real,              intent(in)    :: DXDYH  !< The tracer cell area [m2]
  real,              intent(in)    :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
  real,              intent(in)    :: dens_ratio !< The density of ice divided by the density
                                                 !! of seawater, nondimensional
  real, dimension(2,2), intent(inout) :: Ucontr !< A field related to the subgridscale contributions to
                                                !! the u-direction diagonal elements from basal stress.
  real, dimension(2,2), intent(inout) :: Vcontr !< A field related to the subgridscale contributions to
                                                !! the v-direction diagonal elements from basal stress.
 
  ! bathyT = cellwise-constant bed elevation
 
  integer              :: nsub, i, j, k, l, qx, qy, m, n
  real                 :: subarea, hloc
 
  nsub = size(phisub,1)
  subarea = dxdyh / (nsub**2)
 
  do m=1,2 ; do n=1,2 ; do j=1,nsub ; do i=1,nsub ; do qx=1,2 ; do qy = 1,2
 
    hloc = phisub(i,j,1,1,qx,qy)*h_node(1,1) + phisub(i,j,1,2,qx,qy)*h_node(1,2) + &
           phisub(i,j,2,1,qx,qy)*h_node(2,1) + phisub(i,j,2,2,qx,qy)*h_node(2,2)
 
    if (dens_ratio * hloc - bathyt > 0) then
      ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
      vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
    endif
 
  enddo ; enddo ; enddo ; enddo ; enddo ; enddo
 
end subroutine cg_diagonal_subgrid_basal
 
 
subroutine apply_boundary_values(CS, ISS, G, time, Phisub, H_node, nu, beta, float_cond, &
                                 dens_ratio, u_bdry_contr, v_bdry_contr)
 
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state),  intent(in)    :: ISS !< A structure with elements that describe
                                               !! the ice-shelf state
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  type(time_type),        intent(in)    :: Time !< The current model time
  real, dimension(:,:,:,:,:,:), &
                          intent(in)    :: Phisub !< Quadrature structure weights at subgridscale
                                            !! locations for finite element calculations
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: H_node !< The ice shelf thickness at nodal
                                                 !! (corner) points [Z ~> m].
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: nu   !< A field related to the ice viscosity from Glen's
                                                !! flow law. The exact form and units depend on the
                                                !! basal law exponent.
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(in)    :: beta !< A field related to the nonlinear part of the
                                                !! "linearized" basal stress.  The exact form and
                                                !! units depend on the basal law exponent
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: float_cond !< An array indicating where the ice
                                                !! shelf is floating: 0 if floating, 1 if not.
  real,                   intent(in)    :: dens_ratio !< The density of ice divided by the density
                                                     !! of seawater, nondimensional
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: u_bdry_contr !< Contributions to the zonal ice
                                                !! velocities due to the open boundaries
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                          intent(inout) :: v_bdry_contr !< Contributions to the zonal ice
                                                !! velocities due to the open boundaries
 
! this will be a per-setup function. the boundary values of thickness and velocity
! (and possibly other variables) will be updated in this function
 
  real, dimension(8,4)  :: Phi
  real, dimension(4) :: X, Y
  real, dimension(2) :: xquad
  integer :: i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq
  real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, uq, vq, area, basel
  real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr
 
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
 
  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
 
! X and Y must be passed in the form
    !  3 - 4
    !  |   |
    !  1 - 2
! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
 
  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (iss%hmask(i,j) == 1) then
 
    ! process this cell if any corners have umask set to non-dirichlet bdry.
    ! NOTE: vmask not considered, probably should be
 
    if ((cs%umask(i-1,j-1) == 3) .OR. (cs%umask(i,j-1) == 3) .OR. &
        (cs%umask(i-1,j) == 3) .OR. (cs%umask(i,j) == 3)) then
 
 
      dxh = g%dxT(i,j)
      dyh = g%dyT(i,j)
      dxdyh = g%areaT(i,j)
 
      x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
      x(3:4) = g%geoLonBu(i-1:i,j)*1000
      y(1:2) = g%geoLatBu(i-1:i,j-1)*1000
      y(3:4) = g%geoLatBu(i-1:i,j)*1000
 
      call bilinear_shape_functions(x, y, phi, area)
 
      ! X and Y must be passed in the form
          !  3 - 4
                 !  |   |
          !  1 - 2
      ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
      ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
 
 
 
      do iq=1,2 ; do jq=1,2
 
        uq = cs%u_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
              cs%u_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
             cs%u_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
             cs%u_bdry_val(i,j) * xquad(iq) * xquad(jq)
 
        vq = cs%v_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
              cs%v_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
             cs%v_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
             cs%v_bdry_val(i,j) * xquad(iq) * xquad(jq)
 
        ux = cs%u_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
              cs%u_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
             cs%u_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
             cs%u_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
 
        vx = cs%v_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
              cs%v_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
             cs%v_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
             cs%v_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
 
        uy = cs%u_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
              cs%u_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
             cs%u_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
             cs%u_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
 
        vy = cs%v_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
              cs%v_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
             cs%v_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
             cs%v_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
 
        do iphi=1,2 ; do jphi=1,2
 
          if (iq == iphi) then
            ilq = 2
          else
            ilq = 1
          endif
 
          if (jq == jphi) then
            jlq = 2
          else
            jlq = 1
          endif
 
          if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
 
 
            u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
            .25 * dxdyh * nu(i,j) * ( (4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                         (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
 
            if (float_cond(i,j) == 0) then
              u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
            endif
 
          endif
 
          if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
 
 
            v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
              .25 * dxdyh * nu(i,j) * ( (uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
                           (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
 
            if (float_cond(i,j) == 0) then
              v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
                .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
            endif
 
          endif
        enddo ; enddo
      enddo ; enddo
 
      if (float_cond(i,j) == 1) then
        usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
        ucell(:,:) = cs%u_bdry_val(i-1:i,j-1:j) ; vcell(:,:) = cs%v_bdry_val(i-1:i,j-1:j)
        hcell(:,:) = h_node(i-1:i,j-1:j)
        call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, dxdyh, basel, &
                                     dens_ratio, usubcontr, vsubcontr)
        do iphi=1,2 ; do jphi = 1,2
          if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
            u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
              usubcontr(iphi,jphi) * beta(i,j)
          endif
          if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
            v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
              vsubcontr(iphi,jphi) * beta(i,j)
          endif
        enddo ; enddo
      endif
    endif
  endif ; enddo ; enddo
 
end subroutine apply_boundary_values
 
!> Update depth integrated viscosity, based on horizontal strain rates, and also update the
!! nonlinear part of the basal traction.
subroutine calc_shelf_visc(CS, ISS, G, US, u, v)
  type(ice_shelf_dyn_cs),    intent(inout) :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(in)    :: ISS !< A structure with elements that describe
                                              !! the ice-shelf state
  type(ocean_grid_type), intent(in)    :: G  !< The grid structure used by the ice shelf.
  type(unit_scale_type), intent(in)    :: US !< Pointer to a structure containing unit conversion factors
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                         intent(inout) :: u !< The zonal ice shelf velocity [m year-1].
  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
                         intent(inout) :: v !< The meridional ice shelf velocity [m year-1].
 
! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve
! so there is an "upper" and "lower" bilinear viscosity
 
! also this subroutine updates the nonlinear part of the basal traction
 
! this may be subject to change later... to make it "hybrid"
 
  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
  real :: A, n, ux, uy, vx, vy, eps_min, umid, vmid, unorm, C_basal_friction, n_basal_friction, dxh, dyh, dxdyh
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
  iegq = g%iegB ; jegq = g%jegB
  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
  giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
  is = iscq - 1; js = jscq - 1
 
  a = cs%A_glen_isothermal ; n = cs%n_glen; eps_min = cs%eps_glen_min
  c_basal_friction = cs%C_basal_friction ; n_basal_friction = cs%n_basal_friction
 
  do j=jsd+1,jed-1
    do i=isd+1,ied-1
 
      dxh = g%dxT(i,j)
      dyh = g%dyT(i,j)
      dxdyh = g%areaT(i,j)
 
      if (iss%hmask(i,j) == 1) then
        ux = (u(i,j) + u(i,j-1) - u(i-1,j) - u(i-1,j-1)) / (2*dxh)
        vx = (v(i,j) + v(i,j-1) - v(i-1,j) - v(i-1,j-1)) / (2*dxh)
        uy = (u(i,j) - u(i,j-1) + u(i-1,j) - u(i-1,j-1)) / (2*dyh)
        vy = (v(i,j) - v(i,j-1) + v(i-1,j) - v(i-1,j-1)) / (2*dyh)
 
        cs%ice_visc(i,j) = .5 * a**(-1/n) * &
             (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * &
             us%Z_to_m*iss%h_shelf(i,j)
 
        umid = (u(i,j) + u(i,j-1) + u(i-1,j) + u(i-1,j-1))/4
        vmid = (v(i,j) + v(i,j-1) + v(i-1,j) + v(i-1,j-1))/4
        unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2)
        cs%taub_beta_eff(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
      endif
    enddo
  enddo
 
end subroutine calc_shelf_visc
 
subroutine update_od_ffrac(CS, G, US, ocean_mass, find_avg)
  type(ice_shelf_dyn_cs), intent(inout) :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(inout) :: G  !< The grid structure used by the ice shelf.
  type(unit_scale_type),  intent(in)    :: US !< Pointer to a structure containing unit conversion factors
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: ocean_mass !< The mass per unit area of the ocean [kg m-2].
  logical,                intent(in)    :: find_avg !< If true, find the average of OD and ffrac, and
                                              !! reset the underlying running sums to 0.
 
  integer :: isc, iec, jsc, jec, i, j
  real    :: I_rho_ocean
  real    :: I_counter
 
  i_rho_ocean = 1.0 / (us%Z_to_m*cs%density_ocean_avg)
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
 
  do j=jsc,jec ; do i=isc,iec
    cs%OD_rt(i,j) = cs%OD_rt(i,j) + ocean_mass(i,j)*i_rho_ocean
    if (ocean_mass(i,j)*i_rho_ocean > cs%thresh_float_col_depth) then
      cs%float_frac_rt(i,j) = cs%float_frac_rt(i,j) + 1.0
    endif
  enddo ; enddo
  cs%OD_rt_counter = cs%OD_rt_counter + 1
 
  if (find_avg) then
    i_counter = 1.0 / real(cs%OD_rt_counter)
    do j=jsc,jec ; do i=isc,iec
      cs%float_frac(i,j) = 1.0 - (cs%float_frac_rt(i,j) * i_counter)
      cs%OD_av(i,j) = cs%OD_rt(i,j) * i_counter
 
      cs%OD_rt(i,j) = 0.0 ; cs%float_frac_rt(i,j) = 0.0
    enddo ; enddo
 
    call pass_var(cs%float_frac, g%domain)
    call pass_var(cs%OD_av, g%domain)
  endif
 
end subroutine update_od_ffrac
 
subroutine update_od_ffrac_uncoupled(CS, G, h_shelf)
  type(ice_shelf_dyn_cs), intent(inout) :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: h_shelf !< the thickness of the ice shelf [Z ~> m].
 
  integer :: i, j, iters, isd, ied, jsd, jed
  real    :: rhoi_rhow, OD
 
  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
 
  do j=jsd,jed
    do i=isd,ied
      od = g%bathyT(i,j) - rhoi_rhow * h_shelf(i,j)
      if (od >= 0) then
    ! ice thickness does not take up whole ocean column -> floating
        cs%OD_av(i,j) = od
        cs%float_frac(i,j) = 0.
      else
        cs%OD_av(i,j) = 0.
        cs%float_frac(i,j) = 1.
      endif
    enddo
  enddo
 
end subroutine update_od_ffrac_uncoupled
 
!> This subroutine calculates the gradients of bilinear basis elements that
!! that are centered at the vertices of the cell. values are calculated at
!! points of gaussian quadrature.
subroutine bilinear_shape_functions (X, Y, Phi, area)
  real, dimension(4),   intent(in) :: X !< The x-positions of the vertices of the quadrilateral.
  real, dimension(4),   intent(in) :: Y !< The y-positions of the vertices of the quadrilateral.
  real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian
                                             !! quadrature points surrounding the cell verticies.
  real,                 intent(out) :: area !< The quadrilateral cell area [m2].
 
! X and Y must be passed in the form
    !  3 - 4
    !  |   |
    !  1 - 2
 
! this subroutine calculates the gradients of bilinear basis elements that
! that are centered at the vertices of the cell. values are calculated at
! points of gaussian quadrature. (in 1D: .5 * (1 +/- sqrt(1/3)) for [0,1])
!     (ordered in same way as vertices)
!
! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
! Phi_i is equal to 1 at vertex i, and 0 at vertex k /= i, and bilinear
!
! This should be a one-off; once per nonlinear solve? once per lifetime?
! ... will all cells have the same shape and dimension?
 
  real, dimension(4) :: xquad, yquad
  integer :: node, qpoint, xnode, xq, ynode, yq
  real :: a,b,c,d,e,f,xexp,yexp
 
  xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
  xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
 
  do qpoint=1,4
 
    a = -x(1)*(1-yquad(qpoint)) + x(2)*(1-yquad(qpoint)) - x(3)*yquad(qpoint) + x(4)*yquad(qpoint) ! d(x)/d(x*)
    b = -y(1)*(1-yquad(qpoint)) + y(2)*(1-yquad(qpoint)) - y(3)*yquad(qpoint) + y(4)*yquad(qpoint) ! d(y)/d(x*)
    c = -x(1)*(1-xquad(qpoint)) - x(2)*(xquad(qpoint)) + x(3)*(1-xquad(qpoint)) + x(4)*(xquad(qpoint)) ! d(x)/d(y*)
    d = -y(1)*(1-xquad(qpoint)) - y(2)*(xquad(qpoint)) + y(3)*(1-xquad(qpoint)) + y(4)*(xquad(qpoint)) ! d(y)/d(y*)
 
    do node=1,4
 
      xnode = 2-mod(node,2) ; ynode = ceiling(real(node)/2)
 
      if (ynode == 1) then
        yexp = 1-yquad(qpoint)
      else
        yexp = yquad(qpoint)
      endif
 
      if (1 == xnode) then
        xexp = 1-xquad(qpoint)
      else
        xexp = xquad(qpoint)
      endif
 
      phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c)
      phi(2*node,qpoint) = ( -c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c)
 
    enddo
  enddo
 
  area = quad_area(x, y)
 
end subroutine bilinear_shape_functions
 
 
subroutine bilinear_shape_functions_subgrid(Phisub, nsub)
  real, dimension(nsub,nsub,2,2,2,2), &
           intent(inout) :: Phisub !< Quadrature structure weights at subgridscale
                                   !! locations for finite element calculations
  integer, intent(in)    :: nsub   !< The nubmer of subgridscale quadrature locations in each direction
 
  ! this subroutine is a helper for interpolation of floatation condition
  ! for the purposes of evaluating the terms \int (u,v) \phi_i dx dy in a cell that is
  !     in partial floatation
  ! the array Phisub contains the values of \phi_i (where i is a node of the cell)
  !     at quad point j
  ! i think this general approach may not work for nonrectangular elements...
  !
 
  ! Phisub(i,j,k,l,q1,q2)
  !  i: subgrid index in x-direction
  !  j: subgrid index in y-direction
  !  k: basis function x-index
  !  l: basis function y-index
  !  q1: quad point x-index
  !  q2: quad point y-index
 
  ! e.g. k=1,l=1 => node 1
  !      q1=2,q2=1 => quad point 2
 
    !  3 - 4
    !  |   |
    !  1 - 2
 
  integer :: i, j, k, l, qx, qy, indx, indy
  real,dimension(2)    :: xquad
  real                 :: x0, y0, x, y, val, fracx
 
  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
  fracx = 1.0/real(nsub)
 
  do j=1,nsub
    do i=1,nsub
      x0 = (i-1) * fracx ; y0 = (j-1) * fracx
      do qx=1,2
        do qy=1,2
          x = x0 + fracx*xquad(qx)
          y = y0 + fracx*xquad(qy)
          do k=1,2
            do l=1,2
              val = 1.0
              if (k == 1) then
                val = val * (1.0-x)
              else
                val = val * x
              endif
              if (l == 1) then
                val = val * (1.0-y)
              else
                val = val * y
              endif
              phisub(i,j,k,l,qx,qy) = val
            enddo
          enddo
        enddo
      enddo
    enddo
  enddo
 
end subroutine bilinear_shape_functions_subgrid
 
 
subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face_mask)
  type(ice_shelf_dyn_cs),intent(in)    :: CS !< A pointer to the ice shelf dynamics control structure
  type(ocean_grid_type), intent(inout) :: G  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(out)   :: umask !< A coded mask indicating the nature of the
                                             !! zonal flow at the corner point
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(out)   :: vmask !< A coded mask indicating the nature of the
                                             !! meridional flow at the corner point
  real, dimension(SZDIB_(G),SZDJ_(G)), &
                         intent(out)   :: u_face_mask !< A coded mask for velocities at the C-grid u-face
  real, dimension(SZDI_(G),SZDJB_(G)), &
                         intent(out)   :: v_face_mask !< A coded mask for velocities at the C-grid v-face
  ! sets masks for velocity solve
  ! ignores the fact that their might be ice-free cells - this only considers the computational boundary
 
  ! !!!IMPORTANT!!! relies on thickness mask - assumed that this is called after hmask has been updated & halo-updated
 
  integer :: i, j, k, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
  integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec
  integer :: i_off, j_off
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
  i_off = g%idg_offset ; j_off = g%jdg_offset
  isd = g%isd ; jsd = g%jsd
  iegq = g%iegB ; jegq = g%jegB
  gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
  giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
 
  umask(:,:) = 0 ; vmask(:,:) = 0
  u_face_mask(:,:) = 0 ; v_face_mask(:,:) = 0
 
  if (g%symmetric) then
   is = isd ; js = jsd
  else
   is = isd+1 ; js = jsd+1
  endif
 
  do j=js,g%jed
    do i=is,g%ied
 
      if (hmask(i,j) == 1) then
 
        umask(i-1:i,j-1:j) = 1.
        vmask(i-1:i,j-1:j) = 1.
 
        do k=0,1
 
          select case (int(cs%u_face_mask_bdry(i-1+k,j)))
            case (3)
              umask(i-1+k,j-1:j)=3.
              vmask(i-1+k,j-1:j)=0.
              u_face_mask(i-1+k,j)=3.
            case (2)
              u_face_mask(i-1+k,j)=2.
            case (4)
              umask(i-1+k,j-1:j)=0.
              vmask(i-1+k,j-1:j)=0.
              u_face_mask(i-1+k,j)=4.
            case (0)
              umask(i-1+k,j-1:j)=0.
              vmask(i-1+k,j-1:j)=0.
              u_face_mask(i-1+k,j)=0.
            case (1)  ! stress free x-boundary
              umask(i-1+k,j-1:j)=0.
            case default
          end select
        enddo
 
        do k=0,1
 
          select case (int(cs%v_face_mask_bdry(i,j-1+k)))
            case (3)
              vmask(i-1:i,j-1+k)=3.
              umask(i-1:i,j-1+k)=0.
              v_face_mask(i,j-1+k)=3.
            case (2)
              v_face_mask(i,j-1+k)=2.
            case (4)
              umask(i-1:i,j-1+k)=0.
              vmask(i-1:i,j-1+k)=0.
              v_face_mask(i,j-1+k)=4.
            case (0)
              umask(i-1:i,j-1+k)=0.
              vmask(i-1:i,j-1+k)=0.
              u_face_mask(i,j-1+k)=0.
            case (1) ! stress free y-boundary
              vmask(i-1:i,j-1+k)=0.
            case default
          end select
        enddo
 
        !if (CS%u_face_mask_bdry(i-1,j) >= 0) then !left boundary
        !  u_face_mask(i-1,j) = CS%u_face_mask_bdry(i-1,j)
        !  umask(i-1,j-1:j) = 3.
        !  vmask(i-1,j-1:j) = 0.
        !endif
 
        !if (j_off+j == gjsc+1) then !bot boundary
        !  v_face_mask(i,j-1) = 0.
        !  umask (i-1:i,j-1) = 0.
        !  vmask (i-1:i,j-1) = 0.
        !elseif (j_off+j == gjec) then !top boundary
        !  v_face_mask(i,j) = 0.
        !  umask (i-1:i,j) = 0.
        !  vmask (i-1:i,j) = 0.
        !endif
 
        if (i < g%ied) then
          if ((hmask(i+1,j) == 0) &
              .OR. (hmask(i+1,j) == 2)) then
            !right boundary or adjacent to unfilled cell
            u_face_mask(i,j) = 2.
          endif
        endif
 
        if (i > g%isd) then
          if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
            !adjacent to unfilled cell
            u_face_mask(i-1,j) = 2.
          endif
        endif
 
        if (j > g%jsd) then
          if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
            !adjacent to unfilled cell
            v_face_mask(i,j-1) = 2.
          endif
        endif
 
        if (j < g%jed) then
          if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
            !adjacent to unfilled cell
            v_face_mask(i,j) = 2.
          endif
        endif
 
 
      endif
 
    enddo
  enddo
 
  ! note: if the grid is nonsymmetric, there is a part that will not be transferred with a halo update
  ! so this subroutine must update its own symmetric part of the halo
 
  call pass_vector(u_face_mask, v_face_mask, g%domain, to_all, cgrid_ne)
  call pass_vector(umask, vmask, g%domain, to_all, bgrid_ne)
 
end subroutine update_velocity_masks
 
!> Interpolate the ice shelf thickness from tracer point to nodal points,
!! subject to a mask.
subroutine interpolate_h_to_b(G, h_shelf, hmask, H_node)
  type(ocean_grid_type), intent(inout) :: G  !< The grid structure used by the ice shelf.
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDIB_(G),SZDJB_(G)), &
                         intent(inout) :: H_node !< The ice shelf thickness at nodal (corner)
                                             !! points [Z ~> m].
 
  integer :: i, j, isc, iec, jsc, jec, num_h, k, l
  real    :: summ
 
  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
 
  h_node(:,:) = 0.0
 
  ! H_node is node-centered; average over all cells that share that node
  ! if no (active) cells share the node then its value there is irrelevant
 
  do j=jsc-1,jec
    do i=isc-1,iec
      summ = 0.0
      num_h = 0
      do k=0,1
        do l=0,1
          if (hmask(i+k,j+l) == 1.0) then
            summ = summ + h_shelf(i+k,j+l)
            num_h = num_h + 1
          endif
        enddo
      enddo
      if (num_h > 0) then
        h_node(i,j) = summ / num_h
      endif
    enddo
  enddo
 
  call pass_var(h_node, g%domain, position=corner)
 
end subroutine interpolate_h_to_b
 
!> Deallocates all memory associated with the ice shelf dynamics module
subroutine ice_shelf_dyn_end(CS)
  type(ice_shelf_dyn_cs), pointer   :: cs !< A pointer to the ice shelf dynamics control structure
 
  if (.not.associated(cs)) return
 
  deallocate(cs%u_shelf, cs%v_shelf)
  deallocate(cs%t_shelf, cs%tmask)
  deallocate(cs%u_bdry_val, cs%v_bdry_val, cs%t_bdry_val)
  deallocate(cs%u_face_mask, cs%v_face_mask)
  deallocate(cs%umask, cs%vmask)
 
  deallocate(cs%ice_visc, cs%taub_beta_eff)
  deallocate(cs%OD_rt, cs%OD_av)
  deallocate(cs%float_frac, cs%float_frac_rt)
 
  deallocate(cs)
 
end subroutine ice_shelf_dyn_end
 
 
!> This subroutine updates the vertically averaged ice shelf temperature.
subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time)
  type(ice_shelf_dyn_cs), intent(inout) :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
                                           !! the ice-shelf state
  type(ocean_grid_type), intent(inout) :: G  !< The grid structure used by the ice shelf.
  type(unit_scale_type), intent(in)    :: US !< Pointer to a structure containing unit conversion factors
  real,                  intent(in) :: time_step !< The time step for this update [s].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                         intent(in) :: melt_rate !< basal melt rate [kg m-2 s-1]
  type(time_type),       intent(in) :: Time !< The current model time
 
! 5/23/12 OVS
!    This subroutine takes the velocity (on the Bgrid) and timesteps
!      (HT)_t = - div (uHT) + (adot Tsurf -bdot Tbot) once and then calculates T=HT/H
!
!    The flux overflows are included here. That is because they will be used to advect 3D scalars
!    into partial cells
 
  !
  ! flux_enter: this is to capture flow into partially covered cells; it gives the mass flux into a given
  ! cell across its boundaries.
  ! ###Perhaps flux_enter should be changed into u-face and v-face
  ! ###fluxes, which can then be used in halo updates, etc.
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !  THESE ARE NOT CONSISTENT ==> FIND OUT WHAT YOU IMPLEMENTED
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   o--- (4) ---o
  !   |           |
  !  (1)         (2)
  !   |           |
  !   o--- (3) ---o
  !
 
  real, dimension(SZDI_(G),SZDJ_(G))   :: th_after_uflux, th_after_vflux, TH
  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
  integer                           :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
  real                              :: rho, spy, t_bd, Tsurf, adot
 
  rho = cs%density_ice
  spy = 365 * 86400 ! seconds per year; is there a global constant for this?  No - it is dependent upon a calendar.
 
  adot = 0.1*us%m_to_Z/spy ! for now adot and Tsurf are defined here adot=surf acc 0.1m/yr, Tsurf=-20oC, vary them later
  tsurf = -20.0
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
  flux_enter(:,:,:) = 0.0
 
  th_after_uflux(:,:) = 0.0
  th_after_vflux(:,:) = 0.0
 
  do j=jsd,jed
    do i=isd,ied
      t_bd = cs%t_bdry_val(i,j)
!      if (ISS%hmask(i,j) > 1) then
      if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
        cs%t_shelf(i,j) = cs%t_bdry_val(i,j)
      endif
    enddo
  enddo
 
  do j=jsd,jed
    do i=isd,ied
      th(i,j) = cs%t_shelf(i,j)*iss%h_shelf(i,j)
    enddo
  enddo
 
 
!  call enable_averaging(time_step,Time,CS%diag)
!  call pass_var(h_after_uflux, G%domain)
!  call pass_var(h_after_vflux, G%domain)
!  if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
!  if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
!  call disable_averaging(CS%diag)
 
  call ice_shelf_advect_temp_x(cs, g, time_step/spy, iss%hmask, th, th_after_uflux, flux_enter)
  call ice_shelf_advect_temp_y(cs, g, time_step/spy, iss%hmask, th_after_uflux, th_after_vflux, flux_enter)
 
  do j=jsd,jed
    do i=isd,ied
!      if (ISS%hmask(i,j) == 1) then
      if (iss%h_shelf(i,j) > 0.0) then
        cs%t_shelf(i,j) = th_after_vflux(i,j)/(iss%h_shelf(i,j))
      else
        cs%t_shelf(i,j) = -10.0
      endif
    enddo
  enddo
 
  do j=jsd,jed
    do i=isd,ied
      t_bd = cs%t_bdry_val(i,j)
!      if (ISS%hmask(i,j) > 1) then
      if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
          cs%t_shelf(i,j) = t_bd
!          CS%t_shelf(i,j) = -15.0
      endif
    enddo
  enddo
 
  do j=jsc,jec
    do i=isc,iec
      if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
        if (iss%h_shelf(i,j) > 0.0) then
!         CS%t_shelf(i,j) = CS%t_shelf(i,j) + &
!             time_step*(adot*Tsurf - US%m_to_Z*melt_rate(i,j)*ISS%tfreeze(i,j))/(ISS%h_shelf(i,j))
          cs%t_shelf(i,j) = cs%t_shelf(i,j) + &
              time_step*(adot*tsurf - (3.0*us%m_to_Z/spy)*iss%tfreeze(i,j)) / iss%h_shelf(i,j)
        else
          ! the ice is about to melt away
          ! in this case set thickness, area, and mask to zero
          ! NOTE: not mass conservative
          ! should maybe scale salt & heat flux for this cell
 
          cs%t_shelf(i,j) = -10.0
          cs%tmask(i,j) = 0.0
        endif
      endif
    enddo
  enddo
 
  call pass_var(cs%t_shelf, g%domain)
  call pass_var(cs%tmask, g%domain)
 
  if (cs%debug) then
    call hchksum(cs%t_shelf, "temp after front", g%HI, haloshift=3)
  endif
 
end subroutine ice_shelf_temp
 
 
subroutine ice_shelf_advect_temp_x(CS, G, time_step, hmask, h0, h_after_uflux, flux_enter)
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(inout) :: G  !< The grid structure used by the ice shelf.
  real,                   intent(in)    :: time_step !< The time step for this update [s].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: h0 !< The initial ice shelf thicknesses [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
                                              !! the zonal mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G),4), &
                          intent(inout) :: flux_enter !< The integrated temperature flux into
                                              !! the cell through the 4 cell boundaries [degC Z m2 ~> degC m3]
 
  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
 
  ! if there is an input bdry condition, the thickness there will be set in initialization
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !        o--- (4) ---o
  !        |           |
  !       (1)         (2)
  !        |           |
  !        o--- (3) ---o
  !
 
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
  integer :: i_off, j_off
  logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
  real, dimension(-2:2) :: stencil
  real :: u_face, &  ! positive if out
      flux_diff_cell, phi, dxh, dyh, dxdyh
 
  character (len=1)        :: debug_str
 
 
  is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  do j=jsd+1,jed-1
    if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
        ((j+j_off) >= g%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries
 
      stencil(:) = -1
!     if (i+i_off == G%domain%nihalo+G%domain%nihalo)
      do i=is,ie
 
        if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
             ((i+i_off) >= g%domain%nihalo+1)) then
 
          if (i+i_off == g%domain%nihalo+1) then
            at_west_bdry=.true.
          else
            at_west_bdry=.false.
          endif
 
          if (i+i_off == g%domain%niglobal+g%domain%nihalo) then
            at_east_bdry=.true.
          else
            at_east_bdry=.false.
          endif
 
          if (hmask(i,j) == 1) then
 
            dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
 
            h_after_uflux(i,j) = h0(i,j)
 
            stencil(:) = h0(i-2:i+2,j)  ! fine as long has nx_halo >= 2
 
            flux_diff_cell = 0
 
            ! 1ST DO LEFT FACE
 
            if (cs%u_face_mask(i-1,j) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) * &
                               cs%t_bdry_val(i-1,j) / dxdyh
            else
 
              ! get u-velocity at center of left face
              u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
 
              if (u_face > 0) then !flux is into cell - we need info from h(i-2), h(i-1) if available
 
              ! i may not cover all the cases.. but i cover the realistic ones
 
                if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a
                              ! thickness bdry condition, and the stencil contains it
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
 
                elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then  ! h(i-2) and h(i-1) are valid
                  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
                           (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
 
                else                            ! h(i-1) is valid
                                    ! (o.w. flux would most likely be out of cell)
                                    !  but h(i-2) is not
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(-1)
 
                endif
 
              elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
                if (hmask(i-1,j) * hmask(i+1,j) == 1) then         ! h(i-1) and h(i+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
                             (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
 
                else
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
 
                  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
                    flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
                  endif
                endif
              endif
            endif
 
            ! NEXT DO RIGHT FACE
 
            ! get u-velocity at center of right face
 
            if (cs%u_face_mask(i+1,j) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) *&
                               cs%t_bdry_val(i+1,j)/ dxdyh
            else
 
              u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
 
              if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available
 
                if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a
                                            ! thickness bdry condition, and the stencil contains it
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
 
                elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then  ! h(i+2) and h(i+1) are valid
 
                  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
                      (stencil(1) - phi * (stencil(1)-stencil(0))/2)
 
                else                            ! h(i+1) is valid
                                            ! (o.w. flux would most likely be out of cell)
                                            !  but h(i+2) is not
 
                  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(1)
 
                endif
 
              elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
 
                if (hmask(i-1,j) * hmask(i+1,j) == 1) then         ! h(i-1) and h(i+1) are both valid
 
                  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(1))/2)
 
                else                            ! h(i+1) is valid
                                            ! (o.w. flux would most likely be out of cell)
                                            !  but h(i+2) is not
 
                  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
 
                  if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
 
                    flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step  * stencil(0)
                  endif
 
                endif
 
              endif
 
              h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
 
            endif
 
          elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
 
            if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then
              u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
              flux_enter(i,j,1) = abs(u_face) * g%dyT(i,j) * time_step * cs%t_bdry_val(i-1,j)* &
                                  cs%thickness_bdry_val(i+1,j)
            elseif (cs%u_face_mask(i-1,j) == 4.) then
              flux_enter(i,j,1) = g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)*cs%t_bdry_val(i-1,j)
            endif
 
            if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then
              u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
              flux_enter(i,j,2) = abs(u_face) * g%dyT(i,j) * time_step * cs%t_bdry_val(i+1,j)* &
                                  cs%thickness_bdry_val(i+1,j)
            elseif (cs%u_face_mask(i+1,j) == 4.) then
              flux_enter(i,j,2) = g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j) * cs%t_bdry_val(i+1,j)
            endif
 
!            if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then
              ! this is solely for the purposes of keeping the mask consistent while advancing
              ! the front without having to call pass_var - if cell is empty and cell to left
              ! is ice-covered then this cell will become partly covered
!              hmask(i,j) = 2
!            elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then
              ! this is solely for the purposes of keeping the mask consistent while advancing
              ! the front without having to call pass_var - if cell is empty and cell to left
              ! is ice-covered then this cell will become partly covered
!              hmask(i,j) = 2
 
!            endif
 
          endif
 
        endif
 
      enddo ! i loop
 
    endif
 
  enddo ! j loop
 
end subroutine ice_shelf_advect_temp_x
 
subroutine ice_shelf_advect_temp_y(CS, G, time_step, hmask, h_after_uflux, h_after_vflux, flux_enter)
  type(ice_shelf_dyn_cs), intent(in)    :: CS !< A pointer to the ice shelf control structure
  type(ocean_grid_type),  intent(in)    :: G  !< The grid structure used by the ice shelf.
  real,                   intent(in)    :: time_step !< The time step for this update [s].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: hmask !< A mask indicating which tracer points are
                                             !! partly or fully covered by an ice-shelf
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(in)    :: h_after_uflux !< The ice shelf thicknesses after
                                              !! the zonal mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G)), &
                          intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
                                              !! the meridional mass fluxes [Z ~> m].
  real, dimension(SZDI_(G),SZDJ_(G),4), &
                          intent(inout) :: flux_enter !< The integrated temperature flux into
                                              !! the cell through the 4 cell boundaries [degC Z m2 ~> degC m3]
 
  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
 
  ! if there is an input bdry condition, the thickness there will be set in initialization
 
  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
  !
  !   from left neighbor:   flux_enter(:,:,1)
  !   from right neighbor:  flux_enter(:,:,2)
  !   from bottom neighbor: flux_enter(:,:,3)
  !   from top neighbor:    flux_enter(:,:,4)
  !
  !        o--- (4) ---o
  !        |           |
  !       (1)         (2)
  !        |           |
  !        o--- (3) ---o
  !
 
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
  integer :: i_off, j_off
  logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
  real, dimension(-2:2) :: stencil
  real :: v_face, &  ! positive if out
      flux_diff_cell, phi, dxh, dyh, dxdyh
  character(len=1)        :: debug_str
 
  is = g%isc ; ie = g%iec ; js = g%jsc-1 ; je = g%jec+1 ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  i_off = g%idg_offset ; j_off = g%jdg_offset
 
  do i=isd+2,ied-2
    if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
       ((i+i_off) >= g%domain%nihalo+1)) then  ! based on mehmet's code - only if btw east & west boundaries
 
      stencil(:) = -1
 
      do j=js,je
 
        if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
             ((j+j_off) >= g%domain%njhalo+1)) then
 
          if (j+j_off == g%domain%njhalo+1) then
            at_south_bdry=.true.
          else
            at_south_bdry=.false.
          endif
          if (j+j_off == g%domain%njglobal+g%domain%njhalo) then
            at_north_bdry=.true.
          else
            at_north_bdry=.false.
          endif
 
          if (hmask(i,j) == 1) then
            dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
            h_after_vflux(i,j) = h_after_uflux(i,j)
 
            stencil(:) = h_after_uflux(i,j-2:j+2)  ! fine as long has ny_halo >= 2
            flux_diff_cell = 0
 
            ! 1ST DO south FACE
 
            if (cs%v_face_mask(i,j-1) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) * &
                                 cs%t_bdry_val(i,j-1)/ dxdyh
            else
 
              ! get u-velocity at center of left face
              v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
 
              if (v_face > 0) then !flux is into cell - we need info from h(j-2), h(j-1) if available
 
                ! i may not cover all the cases.. but i cover the realistic ones
 
                if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a
                                            ! thickness bdry condition, and the stencil contains it
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
 
                elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then  ! h(j-2) and h(j-1) are valid
 
                  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
 
                else     ! h(j-1) is valid
                         ! (o.w. flux would most likely be out of cell)
                         !  but h(j-2) is not
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(-1)
                endif
 
              elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
 
                if (hmask(i,j-1) * hmask(i,j+1) == 1) then  ! h(j-1) and h(j+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
                else
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
 
                  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
                    flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
                  endif
 
                endif
 
              endif
 
            endif
 
            ! NEXT DO north FACE
 
            if (cs%v_face_mask(i,j+1) == 4.) then
 
              flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) *&
                               cs%t_bdry_val(i,j+1)/ dxdyh
            else
 
            ! get u-velocity at center of right face
              v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
 
              if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available
 
                if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a
                                            ! thickness bdry condition, and the stencil contains it
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
                elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then  ! h(j+2) and h(j+1) are valid
                  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(1) - phi * (stencil(1)-stencil(0))/2)
                else     ! h(j+1) is valid
                         ! (o.w. flux would most likely be out of cell)
                         !  but h(j+2) is not
                  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
                endif
 
              elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
 
                if (hmask(i,j-1) * hmask(i,j+1) == 1) then         ! h(j-1) and h(j+1) are both valid
                  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
                      (stencil(0) - phi * (stencil(0)-stencil(1))/2)
                else   ! h(j+1) is valid
                       ! (o.w. flux would most likely be out of cell)
                       !  but h(j+2) is not
                  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
                  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
                    flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
                  endif
                endif
 
              endif
 
            endif
 
            h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
 
          elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
 
            if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then
              v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
              flux_enter(i,j,3) = abs(v_face) * g%dxT(i,j) * time_step * cs%t_bdry_val(i,j-1)* &
                                   cs%thickness_bdry_val(i,j-1)
            elseif (cs%v_face_mask(i,j-1) == 4.) then
              flux_enter(i,j,3) = g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)*cs%t_bdry_val(i,j-1)
            endif
 
            if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then
              v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
              flux_enter(i,j,4) = abs(v_face) * g%dxT(i,j) * time_step * cs%t_bdry_val(i,j+1)* &
                                   cs%thickness_bdry_val(i,j+1)
            elseif (cs%v_face_mask(i,j+1) == 4.) then
              flux_enter(i,j,4) = g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)*cs%t_bdry_val(i,j+1)
            endif
 
!            if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then
               ! this is solely for the purposes of keeping the mask consistent while advancing
               ! the front without having to call pass_var - if cell is empty and cell to left
               ! is ice-covered then this cell will become partly covered
 !             hmask(i,j) = 2
 !           elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then
                ! this is solely for the purposes of keeping the mask consistent while advancing the
                ! front without having to call pass_var - if cell is empty and cell to left is
                ! ice-covered then this cell will become partly covered
!              hmask(i,j) = 2
!            endif
 
          endif
        endif
      enddo ! j loop
    endif
  enddo ! i loop
 
end subroutine ice_shelf_advect_temp_y
 
end module mom_ice_shelf_dynamics