MOM_ice_shelf.F90

!> Implements the thermodynamic aspects of ocean / ice-shelf interactions,
!!  along with a crude placeholder for a later implementation of full
!!  ice shelf dynamics, all using the MOM framework and coding style.
module mom_ice_shelf

! This file is part of MOM6. See LICENSE.md for the license.

use mom_constants, only : hlf
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, diag_ctrl, time_type
use mom_diag_mediator, only : enable_averages, 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_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
use mom_dyn_horgrid, only : rescale_dyn_horgrid_bathymetry
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_grid_initialize, only : set_grid_metrics
use mom_fixed_initialization, only : mom_initialize_topography
use mom_fixed_initialization, only : mom_initialize_rotation
use user_initialization, only : user_initialize_topography
use mom_io, only : field_exists, file_exists, mom_read_data, write_version_number
use mom_io, only : slasher, fieldtype
use mom_io, only : write_field, close_file, single_file, multiple
use mom_restart, only : register_restart_field, query_initialized, save_restart
use mom_restart, only : restart_init, restore_state, mom_restart_cs
use mom_time_manager, only : time_type, time_type_to_real, real_to_time, operator(>), operator(-)
use mom_transcribe_grid, only : copy_dyngrid_to_mom_grid, copy_mom_grid_to_dyngrid
use mom_unit_scaling, only : unit_scale_type, unit_scaling_init, fix_restart_unit_scaling
use mom_variables, only : surface
use mom_forcing_type, only : forcing, allocate_forcing_type, mom_forcing_chksum
use mom_forcing_type, only : mech_forcing, allocate_mech_forcing, mom_mech_forcing_chksum
use mom_forcing_type, only : copy_common_forcing_fields
use mom_get_input, only : directories, get_mom_input
use mom_eos, only : calculate_density, calculate_density_derivs, calculate_tfreeze, eos_domain
use mom_eos, only : eos_type, eos_init
use mom_ice_shelf_dynamics, only : ice_shelf_dyn_cs, update_ice_shelf
use mom_ice_shelf_dynamics, only : register_ice_shelf_dyn_restarts, initialize_ice_shelf_dyn
use mom_ice_shelf_dynamics, only : ice_shelf_min_thickness_calve
use mom_ice_shelf_dynamics, only : ice_time_step_cfl, ice_shelf_dyn_end
use mom_ice_shelf_initialize, only : initialize_ice_thickness
!MJH use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary
use mom_ice_shelf_state, only : ice_shelf_state, ice_shelf_state_end, ice_shelf_state_init
use user_shelf_init, only : user_initialize_shelf_mass, user_update_shelf_mass
use user_shelf_init, only : user_ice_shelf_cs
use mom_coms, only : reproducing_sum
use mom_spatial_means, only : global_area_integral
use mom_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum
use time_interp_external_mod, only : init_external_field, time_interp_external
use time_interp_external_mod, only : time_interp_external_init
implicit none ; private

#include <MOM_memory.h>
#ifdef SYMMETRIC_LAND_ICE
#  define GRID_SYM_ .true.
#else
#  define GRID_SYM_ .false.
#endif

public shelf_calc_flux, add_shelf_flux, initialize_ice_shelf, ice_shelf_end
public ice_shelf_save_restart, solo_step_ice_shelf, add_shelf_forces

! 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.

!> Control structure that contains ice shelf parameters and diagnostics handles
type, public :: ice_shelf_cs ; private
  ! Parameters
  type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control
                                          !! structure for the ice shelves
  type(ocean_grid_type) :: grid           !< Grid for the ice-shelf model
  type(unit_scale_type), pointer :: &
    us => null()       !< A structure containing various unit conversion factors
  !type(dyn_horgrid_type), pointer :: dG  !< Dynamic grid for the ice-shelf model
  type(ocean_grid_type), pointer :: ocn_grid => null() !< A pointer to the ocean model grid
                                          !! The rest is private
  real ::   flux_factor = 1.0             !< A factor that can be used to turn off ice shelf
                                          !! melting (flux_factor = 0) [nondim].
  character(len=128) :: restart_output_dir = ' ' !< The directory in which to write restart files
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                          !! the ice-shelf state
  type(ice_shelf_dyn_cs), pointer :: dcs => null() !< The control structure for the ice-shelf dynamics.

  real, pointer, dimension(:,:) :: &
    utide   => null()  !< An unresolved tidal velocity [L T-1 ~> m s-1]

  real :: ustar_bg     !< A minimum value for ustar under ice shelves [Z T-1 ~> m s-1].
  real :: cdrag        !< drag coefficient under ice shelves [nondim].
  real :: g_earth      !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
  real :: cp           !< The heat capacity of sea water [Q degC-1 ~> J kg-1 degC-1].
  real :: rho_ocn      !< A reference ocean density [R ~> kg m-3].
  real :: cp_ice       !< The heat capacity of fresh ice [Q degC-1 ~> J kg-1 degC-1].
  real :: gamma_t      !< The (fixed) turbulent exchange velocity in the
                       !< 2-equation formulation [Z T-1 ~> m s-1].
  real :: salin_ice    !< The salinity of shelf ice [ppt].
  real :: temp_ice     !< The core temperature of shelf ice [degC].
  real :: kv_ice       !< The viscosity of ice [L4 Z-2 T-1 ~> m2 s-1].
  real :: density_ice  !< A typical density of ice [R ~> kg m-3].
  real :: kv_molec     !< The molecular kinematic viscosity of sea water [Z2 T-1 ~> m2 s-1].
  real :: kd_molec_salt!< The molecular diffusivity of salt [Z2 T-1 ~> m2 s-1].
  real :: kd_molec_temp!< The molecular diffusivity of heat [Z2 T-1 ~> m2 s-1].
  real :: lat_fusion   !< The latent heat of fusion [Q ~> J kg-1].
  real :: gamma_t_3eq  !<  Nondimensional heat-transfer coefficient, used in the 3Eq. formulation
  real :: gamma_s_3eq  !<  Nondimensional salt-transfer coefficient, used in the 3Eq. formulation
                       !<  This number should be specified by the user.
  real :: col_mass_melt_threshold !< An ocean column mass below the iceshelf below which melting
                       !! does not occur [R Z ~> kg m-2]
  logical :: mass_from_file !< Read the ice shelf mass from a file every dt

  !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!!

  real :: time_step    !< this is the shortest timestep that the ice shelf sees, and
                       !! is equal to the forcing timestep (it is passed in when the shelf
                       !! is initialized - so need to reorganize MOM driver.
                       !! it will be the prognistic timestep ... maybe.

  logical :: solo_ice_sheet !< whether the ice model is running without being
                            !! coupled to the ocean
  logical :: gl_regularize  !< whether to regularize the floatation condition
                            !! at the grounding line a la Goldberg Holland Schoof 2009
  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)
  logical :: calve_to_mask  !< If true, calve any ice that passes outside of a masked area
  real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
  real :: t0                !< temperature at ocean surface in the restoring region [degC]
  real :: s0                !< Salinity at ocean surface in the restoring region [ppt].
  real :: input_flux        !< Ice volume flux at an upstream open boundary [m3 s-1].
  real :: input_thickness   !< Ice thickness at an upstream open boundary [m].

  type(time_type) :: time                !< The component's time.
  type(eos_type), pointer :: eqn_of_state => null() !< Type that indicates the
                                         !! equation of state to use.
  logical :: active_shelf_dynamics       !< True if the ice shelf mass changes as a result
                                         !! the dynamic ice-shelf model.
  logical :: override_shelf_movement     !< If true, user code specifies the shelf movement
                                         !! instead of using the dynamic ice-shelf mode.
  logical :: isthermo                    !< True if the ice shelf can exchange heat and
                                         !! mass with the underlying ocean.
  logical :: threeeq                     !< If true, the 3 equation consistency equations are
                                         !! used to calculate the flux at the ocean-ice
                                         !! interface.
  logical :: insulator                   !< If true, ice shelf is a perfect insulator
  logical :: const_gamma                 !< If true, gamma_T is specified by the user.
  logical :: constant_sea_level          !< if true, apply an evaporative, heat and salt
                                         !! fluxes. It will avoid large increase in sea level.
  real    :: min_ocean_mass_float        !< The minimum ocean mass per unit area before the ice
                                         !! shelf is considered to float when constant_sea_level
                                         !! is used [R Z ~> kg m-2]
  real    :: cutoff_depth                !< Depth above which melt is set to zero (>= 0) [Z ~> m].
  logical :: find_salt_root              !< If true, if true find Sbdry using a quadratic eq.
  real    :: tfr_0_0                     !< The freezing point at 0 pressure and 0 salinity [degC]
  real    :: dtfr_ds                     !< Partial derivative of freezing temperature with salinity [degC ppt-1]
  real    :: dtfr_dp                     !< Partial derivative of freezing temperature with
                                         !! pressure [degC T2 R-1 L-2 ~> degC Pa-1]
  !>@{ Diagnostic handles
  integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, &
             id_tfreeze = -1, id_tfl_shelf = -1, &
             id_thermal_driving = -1, id_haline_driving = -1, &
             id_u_ml = -1, id_v_ml = -1, id_sbdry = -1, &
             id_h_shelf = -1, id_h_mask = -1, &
             id_surf_elev = -1, id_bathym = -1, &
             id_area_shelf_h = -1, &
             id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1
  !>@}

  integer :: id_read_mass !< An integer handle used in time interpolation of
                          !! the ice shelf mass read from a file
  integer :: id_read_area !< An integer handle used in time interpolation of
                          !! the ice shelf mass read from a file

  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to control diagnostic output.
  type(user_ice_shelf_cs), pointer :: user_cs => null() !< A pointer to the control structure for
                                  !! user-supplied modifications to the ice shelf code.

  logical :: debug                !< If true, write verbose checksums for debugging purposes
                                  !! and use reproducible sums
end type ice_shelf_cs

integer :: id_clock_shelf !< CPU Clock for the ice shelf code
integer :: id_clock_pass !< CPU Clock for group pass calls

contains

!> Calculates fluxes between the ocean and ice-shelf using the three-equations
!! formulation (optional to use just two equations).
!! See \ref section_ICE_SHELF_equations
subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces)
  type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                !! describe the surface state of the ocean.  The
                                                !! intent is only inout to allow for halo updates.
  type(forcing),         intent(inout) :: fluxes !< structure containing pointers to any possible
                                                !! thermodynamic or mass-flux forcing fields.
  type(time_type),       intent(in)    :: time  !< Start time of the fluxes.
  real,                  intent(in)    :: time_step !< Length of time over which these fluxes
                                                !! will be applied [s].
  type(ice_shelf_cs),    pointer       :: cs    !< A pointer to the control structure returned
                                                !! by a previous call to initialize_ice_shelf.
  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces

  ! Local variables
  type(ocean_grid_type), pointer :: g => null()  !< The grid structure used by the ice shelf.
  type(unit_scale_type), pointer :: us => null() !< Pointer to a structure containing
                                                 !! various unit conversion factors
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                                 !! the ice-shelf state

  real, dimension(SZI_(CS%grid)) :: &
    rhoml, &   !< Ocean mixed layer density [R ~> kg m-3].
    dr0_dt, &  !< Partial derivative of the mixed layer density
               !< with temperature [R degC-1 ~> kg m-3 degC-1].
    dr0_ds, &  !< Partial derivative of the mixed layer density
               !< with salinity [R ppt-1 ~> kg m-3 ppt-1].
    p_int      !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa].

  real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
    exch_vel_t, &  !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1]
    exch_vel_s     !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1]

  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
    mass_flux  !< Total mass flux of freshwater across the ice-ocean interface. [R Z L2 T-1 ~> kg/s]
  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
    haline_driving !< (SSS - S_boundary) ice-ocean
               !! interface, positive for melting and negative for freezing.
               !! This is computed as part of the ISOMIP diagnostics.
  real, parameter :: vk    = 0.40 !< Von Karman's constant - dimensionless
  real :: zeta_n = 0.052 !> The fraction of the boundary layer over which the
               !! viscosity is linearly increasing [nondim]. (Was 1/8. Why?)
  real, parameter :: rc    = 0.20     ! critical flux Richardson number.
  real :: i_zeta_n !< The inverse of ZETA_N [nondim].
  real :: i_lf     !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1].
  real :: i_vk     !< The inverse of the Von Karman constant [nondim].
  real :: pr, sc   !< The Prandtl number and Schmidt number [nondim].

  ! 3 equations formulation variables
  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
    sbdry     !< Salinities in the ocean at the interface with the ice shelf [ppt].
  real :: sbdry_it
  real :: sbdry1, sbdry2
  real :: s_a, s_b, s_c  ! Variables used to find salt roots
  real :: ds_it    !< The interface salinity change during an iteration [ppt].
  real :: hbl_neut !< The neutral boundary layer thickness [Z ~> m].
  real :: hbl_neut_h_molec !< The ratio of the neutral boundary layer thickness
                   !! to the molecular boundary layer thickness [nondim].
  real :: wt_flux !< The downward vertical flux of heat just inside the ocean [degC Z T-1 ~> degC m s-1].
  real :: wb_flux !< The downward vertical flux of buoyancy just inside the ocean [Z2 T-3 ~> m2 s-3].
  real :: db_ds   !< The derivative of buoyancy with salinity [Z T-2 ppt-1 ~> m s-2 ppt-1].
  real :: db_dt   !< The derivative of buoyancy with temperature [Z T-2 degC-1 ~> m s-2 degC-1].
  real :: i_n_star ! [nondim]
  real :: n_star_term ! A term in the expression for nstar [T3 Z-2 ~> s3 m-2]
  real :: absf     ! The absolute value of the Coriolis parameter [T-1 ~> s-1]
  real :: dins_dwb !< The partial derivative of I_n_star with wB_flux, in [T3 Z-2 ~> s3 m-2]
  real :: dt_ustar ! The difference between the the freezing point and the ocean boundary layer
                   ! temperature times the friction velocity [degC Z T-1 ~> degC m s-1]
  real :: ds_ustar ! The difference between the salinity at the ice-ocean interface and the ocean
                   ! boundary layer salinity times the friction velocity [ppt Z T-1 ~> ppt m s-1]
  real :: ustar_h  ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1]
  real :: gam_turb ! [nondim]
  real :: gam_mol_t, gam_mol_s ! Relative coefficients of molecular diffusivites [nondim]
  real :: rhocp     ! A typical ocean density times the heat capacity of water [Q R ~> J m-3]
  real :: ln_neut
  real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1]
  real :: sb_min, sb_max
  real :: ds_min, ds_max
  ! Variables used in iterating for wB_flux.
  real :: wb_flux_new, ddwb_dwb_in
  real :: i_gam_t, i_gam_s
  real :: dg_dwb   ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2]
  real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2].
  real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2]
  real :: asu1, asu2   ! Ocean areas covered by ice shelves at neighboring u-
  real :: asv1, asv2   ! and v-points [L2 ~> m2].
  real :: i_au, i_av   ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2]
  real :: irho0        ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1]
  logical :: sb_min_set, sb_max_set
  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
  logical :: coupled_gl     ! If true, the grouding line position is determined based on
                            ! coupled ice-ocean dynamics.

  real, parameter :: c2_3 = 2.0/3.0
  character(len=160) :: mesg  ! The text of an error message
  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
  integer :: i, j, is, ie, js, je, ied, jed, it1, it3

  if (.not. associated(cs)) call mom_error(fatal, "shelf_calc_flux: "// &
       "initialize_ice_shelf must be called before shelf_calc_flux.")
  call cpu_clock_begin(id_clock_shelf)

  g => cs%grid ; us => cs%US
  iss => cs%ISS

  ! useful parameters
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; ied = g%ied ; jed = g%jed
  i_zeta_n = 1.0 / zeta_n
  i_lf = 1.0 / cs%Lat_fusion
  sc = cs%kv_molec/cs%kd_molec_salt
  pr = cs%kv_molec/cs%kd_molec_temp
  i_vk = 1.0/vk
  rhocp = cs%Rho_ocn * cs%Cp

  !first calculate molecular component
  gam_mol_t = 12.5 * (pr**c2_3) - 6.0
  gam_mol_s = 12.5 * (sc**c2_3) - 6.0

  ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS)
  ! these fields are already set to zero during initialization
  ! However, they seem to be changed somewhere and, for diagnostic
  ! reasons, it is better to set them to zero again.
  exch_vel_t(:,:) = 0.0 ; exch_vel_s(:,:) = 0.0
  iss%tflux_shelf(:,:) = 0.0 ; iss%water_flux(:,:) = 0.0
  iss%salt_flux(:,:) = 0.0 ; iss%tflux_ocn(:,:) = 0.0 ; iss%tfreeze(:,:) = 0.0
  ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed.
  haline_driving(:,:) = 0.0
  sbdry(:,:) = sfc_state%sss(:,:)

  !update time
  cs%Time = time

  if (cs%override_shelf_movement) then
    cs%time_step = time_step
    ! update shelf mass
    if (cs%mass_from_file) call update_shelf_mass(g, us, cs, iss, time)
  endif

  if (cs%debug) then
    call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", g%HI, haloshift=0)
    call hchksum(sfc_state%sst, "sst before apply melting", g%HI, haloshift=0)
    call hchksum(sfc_state%sss, "sss before apply melting", g%HI, haloshift=0)
    call hchksum(sfc_state%u, "u_ml before apply melting", g%HI, haloshift=0, scale=us%L_T_to_m_s)
    call hchksum(sfc_state%v, "v_ml before apply melting", g%HI, haloshift=0, scale=us%L_T_to_m_s)
    call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", g%HI, haloshift=0, &
                 scale=us%RZ_to_kg_m2)
  endif

  ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible.
  if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
    call pass_vector(sfc_state%taux_shelf, sfc_state%tauy_shelf, g%domain, to_all, cgrid_ne)
  endif
  irho0 = us%Z_to_L / cs%Rho_ocn
  do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then
    taux2 = 0.0 ; tauy2 = 0.0 ; u2_av = 0.0 ; v2_av = 0.0
    asu1 = (iss%area_shelf_h(i-1,j) + iss%area_shelf_h(i,j))
    asu2 = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j))
    asv1 = (iss%area_shelf_h(i,j-1) + iss%area_shelf_h(i,j))
    asv2 = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1))
    i_au = 0.0 ; if (asu1 + asu2 > 0.0) i_au = 1.0 / (asu1 + asu2)
    i_av = 0.0 ; if (asv1 + asv2 > 0.0) i_av = 1.0 / (asv1 + asv2)
    if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
      taux2 = (asu1 * sfc_state%taux_shelf(i-1,j)**2 + asu2 * sfc_state%taux_shelf(i,j)**2  ) * i_au
      tauy2 = (asv1 * sfc_state%tauy_shelf(i,j-1)**2 + asv2 * sfc_state%tauy_shelf(i,j)**2  ) * i_av
    endif
    u2_av = (asu1 * sfc_state%u(i-1,j)**2 + asu2 * sfc_state%u(i,j)**2) * i_au
    v2_av = (asv1 * sfc_state%v(i,j-1)**2 + asu2 * sfc_state%v(i,j)**2) * i_av

    if (taux2 + tauy2 > 0.0) then
      fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%L_to_Z * &
          sqrt(irho0 * sqrt(taux2 + tauy2) + cs%cdrag*cs%utide(i,j)**2))
    else   ! Take care of the cases when taux_shelf is not set or not allocated.
      fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%L_TO_Z * &
          sqrt(cs%cdrag*((u2_av + v2_av) + cs%utide(i,j)**2)))
    endif
  else ! There is no shelf here.
    fluxes%ustar_shelf(i,j) = 0.0
  endif ; enddo ; enddo

  eosdom(:) = eos_domain(g%HI)
  do j=js,je
    ! Find the pressure at the ice-ocean interface, averaged only over the
    ! part of the cell covered by ice shelf.
    do i=is,ie ; p_int(i) = cs%g_Earth * iss%mass_shelf(i,j) ; enddo

    ! Calculate insitu densities and expansion coefficients
    call calculate_density(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, rhoml(:), &
                                 cs%eqn_of_state, eosdom)
    call calculate_density_derivs(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, dr0_dt, dr0_ds, &
                                 cs%eqn_of_state, eosdom)

    do i=is,ie
      if ((sfc_state%ocean_mass(i,j) > cs%col_mass_melt_threshold) .and. &
          (iss%area_shelf_h(i,j) > 0.0) .and. cs%isthermo) then

        if (cs%threeeq) then
          !   Iteratively determine a self-consistent set of fluxes, with the ocean
          ! salinity just below the ice-shelf as the variable that is being
          ! iterated for.

          ustar_h = fluxes%ustar_shelf(i,j)

          ! Estimate the neutral ocean boundary layer thickness as the minimum of the
          ! reported ocean mixed layer thickness and the neutral Ekman depth.
          absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
                                 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
          if (absf*sfc_state%Hml(i,j) <= vk*ustar_h) then ; hbl_neut = sfc_state%Hml(i,j)
          else ; hbl_neut = (vk*ustar_h) / absf ; endif
          hbl_neut_h_molec = zeta_n * ((hbl_neut * ustar_h) / (5.0 * cs%kv_molec))

          ! Determine the mixed layer buoyancy flux, wB_flux.
          db_ds = (us%L_to_Z**2*cs%g_Earth / rhoml(i)) * dr0_ds(i)
          db_dt = (us%L_to_Z**2*cs%g_Earth / rhoml(i)) * dr0_dt(i)
          ln_neut = 0.0 ; if (hbl_neut_h_molec > 1.0) ln_neut = log(hbl_neut_h_molec)

          if (cs%find_salt_root) then
            ! Solve for the skin salinity using the linearized liquidus parameters and
            ! balancing the turbulent fresh water flux in the near-boundary layer with
            ! the net fresh water or salt added by melting:
            ! (Cp/Lat_fusion)*Gamma_T_3Eq*(TFr_skin-T_ocn) = Gamma_S_3Eq*(S_skin-S_ocn)/S_skin

            ! S_a is always < 0.0 with a realistic expression for the freezing point.
            s_a = cs%dTFr_dS * cs%Gamma_T_3EQ * cs%Cp
            s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%TFr_0_0 + cs%dTFr_dp*p_int(i) - sfc_state%sst(i,j)) - &
                  cs%Lat_fusion * cs%Gamma_S_3EQ    ! S_b Can take either sign, but is usually negative.
            s_c = cs%Lat_fusion * cs%Gamma_S_3EQ * sfc_state%sss(i,j) ! Always >= 0

            if (s_c == 0.0) then  ! The solution for fresh water.
              sbdry(i,j) = 0.0
            elseif (s_a < 0.0) then ! This is the usual ocean case
              if (s_b < 0.0) then ! This is almost always the case
                sbdry(i,j) = 2.0*s_c / (-s_b + sqrt(s_b*s_b - 4.*s_a*s_c))
              else
                sbdry(i,j) = (s_b + sqrt(s_b*s_b - 4.*s_a*s_c)) / (-2.*s_a)
              endif
            elseif ((s_a == 0.0) .and. (s_b < 0.0)) then ! It should be the case that S_b < 0.
              sbdry(i,j) = -s_c / s_b
            else
              call mom_error(fatal, "Impossible conditions found in 3-equation skin salinity calculation.")
            endif

            ! Safety check
            if (sbdry(i,j) < 0.) then
              write(mesg,*) 'sfc_state%sss(i,j) = ',sfc_state%sss(i,j), 'S_a, S_b, S_c', s_a, s_b, s_c
              call mom_error(warning, mesg, .true.)
              write(mesg,*) 'I,J,Sbdry1,Sbdry2',i,j,sbdry1,sbdry2
              call mom_error(warning, mesg, .true.)
              call mom_error(fatal, "shelf_calc_flux: Negative salinity (Sbdry).")
            endif
          else
            ! Guess sss as the iteration starting point for the boundary salinity.
            sbdry(i,j) = sfc_state%sss(i,j) ; sb_max_set = .false.
            sb_min_set = .false.
          endif !find_salt_root

          do it1 = 1,20
            ! Determine the potential temperature at the ice-ocean interface.
            call calculate_tfreeze(sbdry(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state, &
                                   pres_scale=us%RL2_T2_to_Pa)

            dt_ustar = (iss%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h
            ds_ustar = (sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h

            ! First, determine the buoyancy flux assuming no effects of stability
            ! on the turbulence.  Following H & J '99, this limit also applies
            ! when the buoyancy flux is destabilizing.

            if (cs%const_gamma) then ! if using a constant gamma_T
              ! note the different form, here I_Gam_T is NOT 1/Gam_T!
              i_gam_t = cs%Gamma_T_3EQ
              i_gam_s = cs%Gamma_S_3EQ
            else
              gam_turb = i_vk * (ln_neut + (0.5 * i_zeta_n - 1.0))
              i_gam_t = 1.0 / (gam_mol_t + gam_turb)
              i_gam_s = 1.0 / (gam_mol_s + gam_turb)
            endif

            wt_flux = dt_ustar * i_gam_t
            wb_flux = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux

            if (wb_flux < 0.0) then
              ! The buoyancy flux is stabilizing and will reduce the tubulent
              ! fluxes, and iteration is required.
              n_star_term = (zeta_n/rc) * (hbl_neut * vk) / (ustar_h)**3
              do it3 = 1,30
               ! n_star <= 1.0 is the ratio of working boundary layer thickness
               ! to the neutral thickness.
               ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL

                i_n_star = sqrt(1.0 - n_star_term * wb_flux)
                dins_dwb = 0.5 * n_star_term / i_n_star
                if (hbl_neut_h_molec > i_n_star**2) then
                  gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + &
                                    (0.5*i_zeta_n*i_n_star - 1.0))
                  dg_dwb =  i_vk * ( -2.0 / i_n_star + (0.5 * i_zeta_n)) * dins_dwb
                else
                  !   The layer dominated by molecular viscosity is smaller than
                  ! the assumed boundary layer.  This should be rare!
                  gam_turb = i_vk * (0.5 * i_zeta_n*i_n_star - 1.0)
                  dg_dwb = i_vk * (0.5 * i_zeta_n) * dins_dwb
                endif

                if (cs%const_gamma) then ! if using a constant gamma_T
                  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
                  i_gam_t = cs%Gamma_T_3EQ
                  i_gam_s = cs%Gamma_S_3EQ
                else
                  i_gam_t = 1.0 / (gam_mol_t + gam_turb)
                  i_gam_s = 1.0 / (gam_mol_s + gam_turb)
                endif

                wt_flux = dt_ustar * i_gam_t
                wb_flux_new = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux

                ! Find the root where wB_flux_new = wB_flux.  Make the 1.0e-4 below into a parameter?
                if (abs(wb_flux_new - wb_flux) < 1.0e-4*(abs(wb_flux_new) + abs(wb_flux))) exit

                ddwb_dwb_in = dg_dwb * (db_ds * (ds_ustar * i_gam_s**2) + &
                                        db_dt * (dt_ustar * i_gam_t**2)) - 1.0
                ! This is Newton's method without any bounds.  Should bounds be needed?
                wb_flux_new = wb_flux - (wb_flux_new - wb_flux) / ddwb_dwb_in
              enddo !it3
            endif

            iss%tflux_ocn(i,j)  = rhocp * wt_flux
            exch_vel_t(i,j) = ustar_h * i_gam_t
            exch_vel_s(i,j) = ustar_h * i_gam_s

            ! Calculate the heat flux inside the ice shelf.
            ! Vertical adv/diff as in H+J 1999, eqns (26) & approx from (31).
            !   Q_ice = density_ice * CS%Cp_ice * K_ice * dT/dz (at interface)
            ! vertical adv/diff as in H+J 1999, eqs (31) & (26)...
            !   dT/dz ~= min( (lprec/(density_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 )
            ! If this approximation is not made, iterations are required... See H+J Fig 3.

            if (iss%tflux_ocn(i,j) >= 0.0) then
              ! Freezing occurs due to downward ocean heat flux, so zero iout ce heat flux.
              iss%water_flux(i,j) = -i_lf * iss%tflux_ocn(i,j)
              iss%tflux_shelf(i,j) = 0.0
            else
              if (cs%insulator) then
                !no conduction/perfect insulator
                iss%tflux_shelf(i,j) = 0.0
                iss%water_flux(i,j) = i_lf * (iss%tflux_shelf(i,j) - iss%tflux_ocn(i,j))

              else
                ! With melting, from H&J 1999, eqs (31) & (26)...
                !   Q_ice ~= Cp_ice * (CS%Temp_Ice-T_freeze) * lprec
                !   RhoLF*lprec = Q_ice - ISS%tflux_ocn(i,j)
                !   lprec = -(ISS%tflux_ocn(i,j)) / (CS%Lat_fusion + Cp_ice * (T_freeze-CS%Temp_Ice))
                iss%water_flux(i,j) = -iss%tflux_ocn(i,j) / &
                     (cs%Lat_fusion + cs%Cp_ice * (iss%tfreeze(i,j) - cs%Temp_Ice))

                iss%tflux_shelf(i,j) = iss%tflux_ocn(i,j) + cs%Lat_fusion*iss%water_flux(i,j)
              endif

            endif
            !other options: dTi/dz linear through shelf, with draft in [Z ~> m], KTI in [Z2 T-1 ~> m2 s-1]
            !    dTi_dz = (CS%Temp_Ice - ISS%tfreeze(i,j)) / draft(i,j)
            !    ISS%tflux_shelf(i,j) = Rho_Ice * CS%Cp_ice * KTI * dTi_dz


            if (cs%find_salt_root) then
              exit ! no need to do interaction, so exit loop
            else

              mass_exch = exch_vel_s(i,j) * cs%Rho_ocn
              sbdry_it = (sfc_state%sss(i,j) * mass_exch + cs%Salin_ice * iss%water_flux(i,j)) / &
                         (mass_exch + iss%water_flux(i,j))
              ds_it = sbdry_it - sbdry(i,j)
              if (abs(ds_it) < 1.0e-4*(0.5*(sfc_state%sss(i,j) + sbdry(i,j) + 1.0e-10))) exit


              if (ds_it < 0.0) then ! Sbdry is now the upper bound.
                if (sb_max_set .and. (sbdry(i,j) > sb_max)) &
                  call mom_error(fatal,"shelf_calc_flux: Irregular iteration for Sbdry (max).")
                sb_max = sbdry(i,j) ; ds_max = ds_it ; sb_max_set = .true.
              else ! Sbdry is now the lower bound.
                if (sb_min_set .and. (sbdry(i,j) < sb_min)) &
                   call mom_error(fatal, "shelf_calc_flux: Irregular iteration for Sbdry (min).")
                 sb_min = sbdry(i,j) ; ds_min = ds_it ; sb_min_set = .true.
              endif ! dS_it < 0.0

              if (sb_min_set .and. sb_max_set) then
                ! Use the false position method for the next iteration.
                sbdry(i,j) = sb_min + (sb_max-sb_min) * (ds_min / (ds_min - ds_max))
              else
                sbdry(i,j) = sbdry_it
              endif ! Sb_min_set

              sbdry(i,j) = sbdry_it
            endif ! CS%find_salt_root

          enddo !it1
          ! Check for non-convergence and/or non-boundedness?

        else
          !   In the 2-equation form, the mixed layer turbulent exchange velocity
          ! is specified and large enough that the ocean salinity at the interface
          ! is about the same as the boundary layer salinity.

          call calculate_tfreeze(sfc_state%sss(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state, &
                                 pres_scale=us%RL2_T2_to_Pa)

          exch_vel_t(i,j) = cs%gamma_t
          iss%tflux_ocn(i,j) = rhocp * exch_vel_t(i,j) * (iss%tfreeze(i,j) - sfc_state%sst(i,j))
          iss%tflux_shelf(i,j) = 0.0
          iss%water_flux(i,j) = -i_lf * iss%tflux_ocn(i,j)
          sbdry(i,j) = 0.0
        endif
      elseif (iss%area_shelf_h(i,j) > 0.0) then ! This is an ice-sheet, not a floating shelf.
        iss%tflux_ocn(i,j) = 0.0
      else ! There is no ice shelf or sheet here.
        iss%tflux_ocn(i,j) = 0.0
      endif

!      haline_driving(i,j) = sfc_state%sss(i,j) - Sbdry(i,j)

    enddo ! i-loop
  enddo ! j-loop


  do j=js,je ; do i=is,ie
    ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1]
    fluxes%iceshelf_melt(i,j) = iss%water_flux(i,j) * cs%flux_factor

    if ((sfc_state%ocean_mass(i,j) > cs%col_mass_melt_threshold) .and. &
        (iss%area_shelf_h(i,j) > 0.0) .and.  (cs%isthermo)) then

      ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth).
      ! This is needed for the ISOMIP test case.
      if (iss%mass_shelf(i,j) < cs%Rho_ocn*cs%cutoff_depth) then
        iss%water_flux(i,j) = 0.0
        fluxes%iceshelf_melt(i,j) = 0.0
      endif
      ! Compute haline driving, which is one of the diags. used in ISOMIP
      haline_driving(i,j) = (iss%water_flux(i,j) * sbdry(i,j)) / (cs%Rho_ocn * exch_vel_s(i,j))

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!!
      !1)Check if haline_driving computed above is consistent with
      ! haline_driving = sfc_state%sss - Sbdry
      !if (fluxes%iceshelf_melt(i,j) /= 0.0) then
      !   if (haline_driving(i,j) /= (sfc_state%sss(i,j) - Sbdry(i,j))) then
      !     write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',haline_driving(i,j), &
      !                   (sfc_state%sss(i,j) - Sbdry(i,j))
      !     call MOM_error(FATAL, &
      !            "shelf_calc_flux: Inconsistency in melt and haline_driving"//trim(mesg))
      !   endif
      !endif

      ! 2) check if |melt| > 0 when ustar_shelf = 0.
      ! this should never happen
      if ((abs(fluxes%iceshelf_melt(i,j))>0.0) .and. (fluxes%ustar_shelf(i,j) == 0.0)) then
        write(mesg,*) "|melt| = ",fluxes%iceshelf_melt(i,j)," > 0 and ustar_shelf = 0. at i,j", i, j
        call mom_error(fatal, "shelf_calc_flux: "//trim(mesg))
      endif
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
    elseif (iss%area_shelf_h(i,j) > 0.0) then
      ! This is grounded ice, that could be modified to melt if a geothermal heat flux were used.
      haline_driving(i,j) = 0.0
      iss%water_flux(i,j) = 0.0
      fluxes%iceshelf_melt(i,j) = 0.0
    endif ! area_shelf_h

    ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags.
    mass_flux(i,j) = iss%water_flux(i,j) * iss%area_shelf_h(i,j)
  enddo ; enddo ! i- and j-loops

  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
    call cpu_clock_begin(id_clock_pass)
    call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
    call pass_var(iss%mass_shelf, g%domain)
    call cpu_clock_end(id_clock_pass)
  endif

  ! Melting has been computed, now is time to update thickness and mass
  if ( cs%override_shelf_movement .and. (.not.cs%mass_from_file)) then
    call change_thickness_using_melt(iss, g, us, us%s_to_T*time_step, fluxes, cs%density_ice, cs%debug)

    if (cs%debug) then
      call hchksum(iss%h_shelf, "h_shelf after change thickness using melt", g%HI, haloshift=0, scale=us%Z_to_m)
      call hchksum(iss%mass_shelf, "mass_shelf after change thickness using melt", g%HI, haloshift=0, &
                   scale=us%RZ_to_kg_m2)
    endif
  endif

  if (cs%debug) call mom_forcing_chksum("Before add shelf flux", fluxes, g, cs%US, haloshift=0)

  call add_shelf_flux(g, us, cs, sfc_state, fluxes)

  ! now the thermodynamic data is passed on... time to update the ice dynamic quantities

  if (cs%active_shelf_dynamics) then
    update_ice_vel = .false.
    coupled_gl = (cs%GL_couple .and. .not.cs%solo_ice_sheet)

    ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well..
    ! when we decide on how to do it
    call update_ice_shelf(cs%dCS, iss, g, us, us%s_to_T*time_step, time, &
                          sfc_state%ocean_mass, coupled_gl)

  endif

  call enable_averaging(time_step,time,cs%diag)
  if (cs%id_shelf_mass > 0) call post_data(cs%id_shelf_mass, iss%mass_shelf, cs%diag)
  if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
  if (cs%id_ustar_shelf > 0) call post_data(cs%id_ustar_shelf, fluxes%ustar_shelf, cs%diag)
  if (cs%id_melt > 0) call post_data(cs%id_melt, fluxes%iceshelf_melt, cs%diag)
  if (cs%id_thermal_driving > 0) call post_data(cs%id_thermal_driving, (sfc_state%sst-iss%tfreeze), cs%diag)
  if (cs%id_Sbdry > 0) call post_data(cs%id_Sbdry, sbdry, cs%diag)
  if (cs%id_haline_driving > 0) call post_data(cs%id_haline_driving, haline_driving, cs%diag)
  if (cs%id_mass_flux > 0) call post_data(cs%id_mass_flux, mass_flux, cs%diag)
  if (cs%id_u_ml > 0) call post_data(cs%id_u_ml, sfc_state%u, cs%diag)
  if (cs%id_v_ml > 0) call post_data(cs%id_v_ml, sfc_state%v, cs%diag)
  if (cs%id_tfreeze > 0) call post_data(cs%id_tfreeze, iss%tfreeze, cs%diag)
  if (cs%id_tfl_shelf > 0) call post_data(cs%id_tfl_shelf, iss%tflux_shelf, cs%diag)
  if (cs%id_exch_vel_t > 0) call post_data(cs%id_exch_vel_t, exch_vel_t, cs%diag)
  if (cs%id_exch_vel_s > 0) call post_data(cs%id_exch_vel_s, exch_vel_s, cs%diag)
  if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf, iss%h_shelf, cs%diag)
  if (cs%id_h_mask > 0) call post_data(cs%id_h_mask,iss%hmask,cs%diag)
  call disable_averaging(cs%diag)

  if (present(forces)) then
    call add_shelf_forces(g, us, cs, forces, do_shelf_area=(cs%active_shelf_dynamics .or. &
                                                        cs%override_shelf_movement))
  endif

  call cpu_clock_end(id_clock_shelf)

  if (cs%debug) call mom_forcing_chksum("End of shelf calc flux", fluxes, g, cs%US, haloshift=0)

end subroutine shelf_calc_flux

!> Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting
subroutine change_thickness_using_melt(ISS, G, US, time_step, fluxes, density_ice, debug)
  type(ocean_grid_type), intent(inout) :: G  !< The ocean's grid structure.
  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
                                              !! the ice-shelf state
  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
  real,                  intent(in)    :: time_step !< The time step for this update [T ~> s].
  type(forcing),         intent(inout) :: fluxes !< structure containing pointers to any possible
                                                 !! thermodynamic or mass-flux forcing fields.
  real,                  intent(in)    :: density_ice !< The density of ice-shelf ice [R ~> kg m-3].
  logical,     optional, intent(in)    :: debug !< If present and true, write chksums

  ! locals
  real :: I_rho_ice ! Ice specific volume [R-1 ~> m3 kg-1]
  integer :: i, j

  i_rho_ice = 1.0 / density_ice


  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
      ! first, zero out fluxes applied during previous time step
      if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
      if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
      if (associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
      if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0

      if (iss%water_flux(i,j) * time_step / density_ice < iss%h_shelf(i,j)) then
        iss%h_shelf(i,j) = iss%h_shelf(i,j) - iss%water_flux(i,j) * time_step / density_ice
      else
        ! the ice is about to melt away, so set thickness, area, and mask to zero
        ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
        iss%h_shelf(i,j) = 0.0
        iss%hmask(i,j) = 0.0
        iss%area_shelf_h(i,j) = 0.0
      endif
      iss%mass_shelf(i,j) = iss%h_shelf(i,j) * density_ice
    endif
  enddo ; enddo

  call pass_var(iss%area_shelf_h, g%domain)
  call pass_var(iss%h_shelf, g%domain)
  call pass_var(iss%hmask, g%domain)
  call pass_var(iss%mass_shelf, g%domain)

end subroutine change_thickness_using_melt

!> This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on
!! the ice state in ice_shelf_CS.
subroutine add_shelf_forces(G, US, CS, forces, do_shelf_area)
  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    pointer       :: cs   !< This module's control structure.
  type(mech_forcing),    intent(inout) :: forces !< A structure with the driving mechanical forces
  logical, optional,     intent(in)    :: do_shelf_area !< If true find the shelf-covered areas.

  real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 T-1 R-1 Z-2 ~> m5 kg-1 s-1].
  real :: press_ice  ! The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa].
  logical :: find_area ! If true find the shelf areas at u & v points.
  type(ice_shelf_state), pointer :: iss => null() ! A structure with elements that describe
                                          ! the ice-shelf state

  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed

  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
      (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
    call mom_error(fatal,"add_shelf_forces: Incompatible ocean and ice shelf grids.")

  iss => cs%ISS

  find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area

  if (find_area) then
    ! The frac_shelf is set over the widest possible area. Could it be smaller?
    do j=jsd,jed ; do i=isd,ied-1
      forces%frac_shelf_u(i,j) = 0.0
      if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%areaCu(I,j) > 0.0)) &
        forces%frac_shelf_u(i,j) = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j)) / &
                                   (g%areaT(i,j) + g%areaT(i+1,j))
    enddo ; enddo
    do j=jsd,jed-1 ; do i=isd,ied
      forces%frac_shelf_v(i,j) = 0.0
      if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%areaCv(i,J) > 0.0)) &
        forces%frac_shelf_v(i,j) = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1)) / &
                                   (g%areaT(i,j) + g%areaT(i,j+1))
    enddo ; enddo
    call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain, to_all, cgrid_ne)
  endif

  do j=js,je ; do i=is,ie
    press_ice = (iss%area_shelf_h(i,j) * g%IareaT(i,j)) * (cs%g_Earth * iss%mass_shelf(i,j))
    if (associated(forces%p_surf)) then
      if (.not.forces%accumulate_p_surf) forces%p_surf(i,j) = 0.0
      forces%p_surf(i,j) = forces%p_surf(i,j) + press_ice
    endif
    if (associated(forces%p_surf_full)) then
      if (.not.forces%accumulate_p_surf) forces%p_surf_full(i,j) = 0.0
      forces%p_surf_full(i,j) = forces%p_surf_full(i,j) + press_ice
    endif
  enddo ; enddo

  ! For various reasons, forces%rigidity_ice_[uv] is always updated here. Note
  ! that it may have been zeroed out where IOB is translated to forces and
  ! contributions from icebergs and the sea-ice pack added subsequently.
  !### THE RIGIDITY SHOULD ALSO INCORPORATE AREAL-COVERAGE INFORMATION.
  kv_rho_ice = cs%kv_ice / cs%density_ice
  do j=js,je ; do i=is-1,ie
    if (.not.forces%accumulate_rigidity) forces%rigidity_ice_u(i,j) = 0.0
    forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
            kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i+1,j))
  enddo ; enddo
  do j=js-1,je ; do i=is,ie
    if (.not.forces%accumulate_rigidity) forces%rigidity_ice_v(i,j) = 0.0
    forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
            kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i,j+1))
  enddo ; enddo

  if (cs%debug) then
    call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, forces%rigidity_ice_v, &
                  g%HI, symmetric=.true., scale=us%L_to_m**3*us%L_to_Z*us%s_to_T)
    call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, forces%frac_shelf_v, &
                  g%HI, symmetric=.true.)
  endif

end subroutine add_shelf_forces

!> This subroutine adds the ice shelf pressure to the fluxes type.
subroutine add_shelf_pressure(G, US, CS, fluxes)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
  type(ice_shelf_CS),    intent(in)    :: CS   !< This module's control structure.
  type(forcing),         intent(inout) :: fluxes  !< A structure of surface fluxes that may be updated.

  real :: press_ice       !< The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa].
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
      (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
    call mom_error(fatal,"add_shelf_pressure: Incompatible ocean and ice shelf grids.")

  do j=js,je ; do i=is,ie
    press_ice = (cs%ISS%area_shelf_h(i,j) * g%IareaT(i,j)) * (cs%g_Earth * cs%ISS%mass_shelf(i,j))
    if (associated(fluxes%p_surf)) then
      if (.not.fluxes%accumulate_p_surf) fluxes%p_surf(i,j) = 0.0
      fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + press_ice
    endif
    if (associated(fluxes%p_surf_full)) then
      if (.not.fluxes%accumulate_p_surf) fluxes%p_surf_full(i,j) = 0.0
      fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + press_ice
    endif
  enddo ; enddo

end subroutine add_shelf_pressure

!> Updates surface fluxes that are influenced by sub-ice-shelf melting
subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    pointer       :: cs   !< This module's control structure.
  type(surface),         intent(inout) :: sfc_state !< Surface ocean state
  type(forcing),         intent(inout) :: fluxes  !< A structure of surface fluxes that may be used/updated.

  ! local variables
  real :: frac_shelf       !< The fractional area covered by the ice shelf [nondim].
  real :: frac_open        !< The fractional area of the ocean that is not covered by the ice shelf [nondim].
  real :: delta_mass_shelf !< Change in ice shelf mass over one time step [R Z m2 T-1 ~> kg s-1]
  real :: balancing_flux   !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1]
  real :: balancing_area   !< total area where the balancing flux is applied [m2]
  type(time_type) :: dtime !< The time step as a time_type
  type(time_type) :: time0 !< The previous time (Time-dt)
  real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac  !< Fraction of the cel1 where the mass flux
                          !! balancing the net melt flux occurs, 0 to 1 [nondim]
  real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
                          !! at at previous time (Time-dt) [R Z ~> kg m-2]
  real, dimension(SZDI_(G),SZDJ_(G)) :: delta_float_mass   !< The change in the floating mass between
                          !! the two timesteps at (Time) and (Time-dt) [R Z ~> kg m-2].
  real, dimension(SZDI_(G),SZDJ_(G))  :: last_h_shelf !< Ice shelf thickness [Z ~> m]
                          !! at at previous time (Time-dt)
  real, dimension(SZDI_(G),SZDJ_(G))  :: last_hmask !< Ice shelf mask [nondim]
                          !! at at previous time (Time-dt)
  real, dimension(SZDI_(G),SZDJ_(G))  :: last_area_shelf_h !< Ice shelf area [L2 ~> m2]
                          !! at at previous time (Time-dt)
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                          !! the ice-shelf state

  character(len=160) :: mesg  ! The text of an error message
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed

  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
      (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
    call mom_error(fatal,"add_shelf_flux: Incompatible ocean and ice shelf grids.")

  iss => cs%ISS


  call add_shelf_pressure(g, us, cs, fluxes)

  ! Determine ustar and the square magnitude of the velocity in the
  ! bottom boundary layer. Together these give the TKE source and
  ! vertical decay scale.

  if (cs%debug) then
    if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
      call uvchksum("tau[xy]_shelf", sfc_state%taux_shelf, sfc_state%tauy_shelf, &
                    g%HI, haloshift=0, scale=us%RZ_T_to_kg_m2s*us%L_T_to_m_s)
    endif
  endif

  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
    do j=jsd,jed ; do i=isd,ied
      if (g%areaT(i,j) > 0.0) &
        fluxes%frac_shelf_h(i,j) = min(1.0, iss%area_shelf_h(i,j) * g%IareaT(i,j))
    enddo ; enddo
  endif

  if (cs%debug) then
    call mom_forcing_chksum("Before adding shelf fluxes", fluxes, g, cs%US, haloshift=0)
  endif

  do j=js,je ; do i=is,ie ; if (iss%area_shelf_h(i,j) > 0.0) then
    ! Replace fluxes intercepted by the ice shelf with fluxes from the ice shelf
    frac_shelf = min(1.0, iss%area_shelf_h(i,j) * g%IareaT(i,j))
    frac_open = max(0.0, 1.0 - frac_shelf)

    if (associated(fluxes%sw)) fluxes%sw(i,j) = frac_open * fluxes%sw(i,j)
    if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = frac_open * fluxes%sw_vis_dir(i,j)
    if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = frac_open * fluxes%sw_vis_dif(i,j)
    if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = frac_open * fluxes%sw_nir_dir(i,j)
    if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = frac_open * fluxes%sw_nir_dif(i,j)
    if (associated(fluxes%lw)) fluxes%lw(i,j) = frac_open * fluxes%lw(i,j)
    if (associated(fluxes%latent)) fluxes%latent(i,j) = frac_open * fluxes%latent(i,j)
    if (associated(fluxes%evap)) fluxes%evap(i,j) = frac_open * fluxes%evap(i,j)
    if (associated(fluxes%lprec)) then
      if (iss%water_flux(i,j) > 0.0) then
        fluxes%lprec(i,j) =  frac_shelf*iss%water_flux(i,j)*cs%flux_factor + frac_open * fluxes%lprec(i,j)
      else
        fluxes%lprec(i,j) = frac_open * fluxes%lprec(i,j)
        fluxes%evap(i,j) = fluxes%evap(i,j) + frac_shelf*iss%water_flux(i,j)*cs%flux_factor
      endif
    endif

    if (associated(fluxes%sens)) &
      fluxes%sens(i,j) = frac_shelf*iss%tflux_ocn(i,j)*cs%flux_factor + frac_open * fluxes%sens(i,j)
    ! The salt flux should be mostly from sea ice, so perhaps none should be intercepted and this should be changed.
    if (associated(fluxes%salt_flux)) &
      fluxes%salt_flux(i,j) = frac_shelf * iss%salt_flux(i,j)*cs%flux_factor + frac_open * fluxes%salt_flux(i,j)
  endif ; enddo ; enddo

  if (cs%debug) then
    call hchksum(iss%water_flux, "water_flux add shelf fluxes", g%HI, haloshift=0, scale=us%RZ_T_to_kg_m2s)
    call hchksum(iss%tflux_ocn, "tflux_ocn add shelf fluxes", g%HI, haloshift=0, scale=us%QRZ_T_to_W_m2)
    call mom_forcing_chksum("After adding shelf fluxes", fluxes, g, cs%US, haloshift=0)
  endif

  ! Keep sea level constant by removing mass via a balancing flux that might be applied
  ! in the open ocean or the sponge region (via virtual precip, vprec). Apply additional
  ! salt/heat fluxes so that the resultant surface buoyancy forcing is ~ 0.
  ! This is needed for some of the ISOMIP+ experiments.

  if (cs%constant_sea_level) then
    if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je))
    if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je))
    fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0

    ! take into account changes in mass (or thickness) when imposing ice shelf mass
    if (cs%override_shelf_movement .and. cs%mass_from_file) then
      dtime = real_to_time(cs%time_step)

      ! Compute changes in mass after at least one full time step
      if (cs%Time > dtime) then
        time0 = cs%Time - dtime
        do j=js,je ; do i=is,ie
          last_hmask(i,j) = iss%hmask(i,j) ; last_area_shelf_h(i,j) = iss%area_shelf_h(i,j)
        enddo ; enddo
        call time_interp_external(cs%id_read_mass, time0, last_mass_shelf)
        do j=js,je ; do i=is,ie
        ! This should only be done if time_interp_external did an update.
          last_mass_shelf(i,j) = us%kg_m3_to_R*us%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp
          last_h_shelf(i,j) = last_mass_shelf(i,j) / cs%density_ice
        enddo ; enddo

        ! apply calving
        if (cs%min_thickness_simple_calve > 0.0) then
          call ice_shelf_min_thickness_calve(g, last_h_shelf, last_area_shelf_h, last_hmask, &
                                       cs%min_thickness_simple_calve, halo=0)
          ! convert to mass again
          do j=js,je ; do i=is,ie
            last_mass_shelf(i,j) = last_h_shelf(i,j) * cs%density_ice
          enddo ; enddo
        endif

        ! get total ice shelf mass at (Time-dt) and (Time), in kg
        do j=js,je ; do i=is,ie
          ! Just consider the change in the mass of the floating shelf.
          if ((sfc_state%ocean_mass(i,j) > cs%min_ocean_mass_float) .and. &
              (iss%area_shelf_h(i,j) > 0.0)) then
            delta_float_mass(i,j) = iss%mass_shelf(i,j) - last_mass_shelf(i,j)
          else
            delta_float_mass(i,j) = 0.0
          endif
        enddo ; enddo
        delta_mass_shelf = us%kg_m2s_to_RZ_T*(global_area_integral(delta_float_mass, g, scale=us%RZ_to_kg_m2, &
                                                                   area=iss%area_shelf_h) / cs%time_step)
      else! first time step
        delta_mass_shelf = 0.0
      endif
    else ! ice shelf mass does not change
      delta_mass_shelf = 0.0
    endif

    ! average total melt flux over sponge area
    do j=js,je ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .AND. (iss%area_shelf_h(i,j) * g%IareaT(i,j) < 1.0)) then
         ! Uncomment this for some ISOMIP cases:
         !  .AND. (G%geoLonT(i,j) >= 790.0) .AND. (G%geoLonT(i,j) <= 800.0)) then
        bal_frac(i,j) = max(1.0 - iss%area_shelf_h(i,j) * g%IareaT(i,j), 0.0)
      else
        bal_frac(i,j) = 0.0
      endif
    enddo ; enddo

    balancing_area = global_area_integral(bal_frac, g)
    if (balancing_area > 0.0) then
      balancing_flux = ( us%kg_m2s_to_RZ_T*global_area_integral(iss%water_flux, g, scale=us%RZ_T_to_kg_m2s, &
                                                                area=iss%area_shelf_h) + &
                         delta_mass_shelf ) / balancing_area
    else
      balancing_flux = 0.0
    endif

    ! apply fluxes
    do j=js,je ; do i=is,ie
      if (bal_frac(i,j) > 0.0) then
        ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1]
        fluxes%vprec(i,j) = -balancing_flux
        fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0 ! [ Q R Z T-1 ~> W /m^2 ]
        fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3 ! [kgSalt/kg R Z T-1 ~> kgSalt m-2 s-1]
      endif
    enddo ; enddo

    if (cs%debug) then
      write(mesg,*) 'Balancing flux (kg/(m^2 s)), dt = ', balancing_flux*us%RZ_T_to_kg_m2s, cs%time_step
      call mom_mesg(mesg)
      call mom_forcing_chksum("After constant sea level", fluxes, g, cs%US, haloshift=0)
    endif

  endif ! constant_sea_level

end subroutine add_shelf_flux


!> Initializes shelf model data, parameters and diagnostics
subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, 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_cs),           pointer       :: cs   !< A pointer to the ice shelf control structure
  type(diag_ctrl),    target,   intent(in)    :: diag !< A structure that is used to regulate the diagnostic output.
  type(forcing),      optional, intent(inout) :: fluxes !< A structure containing pointers to any possible
                                                   !! thermodynamic or mass-flux forcing fields.
  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
  type(time_type),    optional, intent(in)    :: time_in !< The time at initialization.
  logical,            optional, intent(in)    :: solo_ice_sheet_in !< If present, this indicates whether
                                                   !! a solo ice-sheet driver.

  type(ocean_grid_type), pointer :: g  => null(), og  => null() ! Pointers to grids for convenience.
  type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
                                                 ! various unit conversion factors
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                          !! the ice-shelf state
  type(directories)  :: dirs
  type(dyn_horgrid_type), pointer :: dg => null()
  real    :: z_rescale  ! A rescaling factor for heights from the representation in
                        ! a restart file to the internal representation in this run.
  real    :: rz_rescale ! A rescaling factor for mass loads from the representation in
                        ! a restart file to the internal representation in this run.
  real    :: l_rescale  ! A rescaling factor for horizontal lengths from the representation in
                        ! a restart file to the internal representation in this run.
  real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic.
  real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered
                        ! to be floating when CONST_SEA_LEVEL = True [Z ~> m].
  real :: cdrag, drag_bg_vel
  logical :: new_sim, save_ic, var_force
  !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)  :: mdl = "MOM_ice_shelf"  ! This module's name.
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, isdq, iedq, jsdq, jedq
  integer :: wd_halos(2)
  logical :: read_tideamp, shelf_mass_is_dynamic, debug
  character(len=240) :: tideamp_file
  real    :: utide  ! A tidal velocity [L T-1 ~> m s-1]
  real    :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting
                       ! does not occur [Z ~> m]
  if (associated(cs)) then
    call mom_error(fatal, "MOM_ice_shelf.F90, initialize_ice_shelf: "// &
                          "called with an associated control structure.")
    return
  endif
  allocate(cs)

  !   Go through all of the infrastructure initialization calls, since this is
  ! being treated as an independent component that just happens to use the
  ! MOM's grid and infrastructure.
  call get_mom_input(dirs=dirs)

  ! Determining the internal unit scaling factors for this run.
  call unit_scaling_init(param_file, cs%US)

  ! Set up the ice-shelf domain and grid
  wd_halos(:)=0
  call mom_domains_init(cs%grid%domain, param_file, min_halo=wd_halos, symmetric=grid_sym_)
  ! call diag_mediator_init(CS%grid,param_file,CS%diag)
  ! this needs to be fixed - will probably break when not using coupled driver 0
  call mom_grid_init(cs%grid, param_file, cs%US)

  call create_dyn_horgrid(dg, cs%grid%HI)
  call clone_mom_domain(cs%grid%Domain, dg%Domain)

  call set_grid_metrics(dg, param_file, cs%US)
  ! call set_diag_mediator_grid(CS%grid, CS%diag)

  ! The ocean grid possibly uses different symmetry.
  if (associated(ocn_grid)) then ; cs%ocn_grid => ocn_grid
  else ; cs%ocn_grid => cs%grid ; endif

  ! Convenience pointers
  g => cs%grid
  og => cs%ocn_grid
  us => cs%US

  if (is_root_pe()) then
    write(0,*) 'OG: ', og%isd, og%isc, og%iec, og%ied, og%jsd, og%jsc, og%jsd, og%jed
    write(0,*) 'IG: ', g%isd, g%isc, g%iec, g%ied, g%jsd, g%jsc, g%jsd, g%jed
  endif

  cs%diag => diag

  ! Are we being called from the solo ice-sheet driver? When called by the ocean
  ! model solo_ice_sheet_in is not preset.
  cs%solo_ice_sheet = .false.
  if (present(solo_ice_sheet_in)) cs%solo_ice_sheet = solo_ice_sheet_in

  if (present(time_in)) time = time_in

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
  isdq = g%IsdB ; iedq = g%IedB ; jsdq = g%JsdB ; jedq = g%JedB

  cs%override_shelf_movement = .false. ; cs%active_shelf_dynamics = .false.

  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.)
  if (shelf_mass_is_dynamic) then
    call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
                 "If true, user provided code specifies the ice-shelf "//&
                 "movement instead of the dynamic ice model.", default=.false.)
    cs%active_shelf_dynamics = .not.cs%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_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.
  endif

  call get_param(param_file, mdl, "SHELF_THERMO", cs%isthermo, &
                 "If true, use a thermodynamically interactive ice shelf.", &
                 default=.false.)
  call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%Lat_fusion, &
                 "The latent heat of fusion.", units="J/kg", default=hlf, scale=us%J_kg_to_Q)
  call get_param(param_file, mdl, "SHELF_THREE_EQN", cs%threeeq, &
                 "If true, use the three equation expression of "//&
                 "consistency to calculate the fluxes at the ice-ocean "//&
                 "interface.", default=.true.)
  call get_param(param_file, mdl, "SHELF_INSULATOR", cs%insulator, &
                 "If true, the ice shelf is a perfect insulatior "//&
                 "(no conduction).", default=.false.)
  call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", cs%cutoff_depth, &
                 "Depth above which the melt is set to zero (it must be >= 0) "//&
                 "Default value won't affect the solution.", units="m", default=0.0, scale=us%m_to_Z)
  if (cs%cutoff_depth < 0.) &
    call mom_error(warning,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.")

  call get_param(param_file, mdl, "CONST_SEA_LEVEL", cs%constant_sea_level, &
                 "If true, apply evaporative, heat and salt fluxes in "//&
                 "the sponge region. This will avoid a large increase "//&
                 "in sea level. This option is needed for some of the "//&
                 "ISOMIP+ experiments (Ocean3 and Ocean4). "//&
                 "IMPORTANT: it is not currently possible to do "//&
                 "prefect restarts using this flag.", default=.false.)
  call get_param(param_file, mdl, "MIN_OCEAN_FLOAT_THICK", dz_ocean_min_float, &
                 "The minimum ocean thickness above which the ice shelf is considered to be "//&
                 "floating when CONST_SEA_LEVEL = True.", &
                 default=0.1, units="m", scale=us%m_to_Z, do_not_log=.not.cs%constant_sea_level)

  call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", cs%S0, &
                 "Surface salinity in the restoring region.", &
                default=33.8, units='ppt', do_not_log=.true.)

  call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", cs%T0, &
                "Surface temperature in the restoring region.", &
                default=-1.9, units='degC', do_not_log=.true.)

  call get_param(param_file, mdl, "SHELF_3EQ_GAMMA", cs%const_gamma, &
                 "If true, user specifies a constant nondimensional heat-transfer coefficient "//&
                 "(GAMMA_T_3EQ), from which the default salt-transfer coefficient is set "//&
                 "as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
  if (cs%threeeq) then
    call get_param(param_file, mdl, "SHELF_S_ROOT", cs%find_salt_root, &
                 "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//&
                 "is computed from a quadratic equation. Otherwise, the previous "//&
                 "interactive method to estimate Sbdry is used.", default=.false.)
  else
    call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", cs%gamma_t, &
                 "If SHELF_THREE_EQN is false, this the fixed turbulent "//&
                 "exchange velocity at the ice-ocean interface.", &
                 units="m s-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
  endif
  if (cs%const_gamma .or. cs%find_salt_root) then
    call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
                 "Nondimensional heat-transfer coefficient.", &
                  units="nondim", default=2.2e-2)
    call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_S", cs%Gamma_S_3EQ, &
                 "Nondimensional salt-transfer coefficient.", &
                 default=cs%Gamma_T_3EQ/35.0, units="nondim")
  endif

  call get_param(param_file, mdl, "ICE_SHELF_MASS_FROM_FILE", &
                 cs%mass_from_file, "Read the mass of the "//&
                 "ice shelf (every time step) from a file.", default=.false.)

  if (cs%find_salt_root) then ! read liquidus coeffs.
     call get_param(param_file, mdl, "TFREEZE_S0_P0", cs%TFr_0_0, &
                 "this is the freezing potential temperature at "//&
                 "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.)
    call get_param(param_file, mdl, "DTFREEZE_DS", cs%dTFr_dS, &
                 "this is the derivative of the freezing potential temperature with salinity.", &
                 units="degC psu-1", default=-0.054, do_not_log=.true.)
    call get_param(param_file, mdl, "DTFREEZE_DP", cs%dTFr_dp, &
                 "this is the derivative of the freezing potential temperature with pressure.", &
                 units="degC Pa-1", default=0.0, scale=us%RL2_T2_to_Pa, do_not_log=.true.)
  endif

  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80, scale=us%m_s_to_L_T**2*us%Z_to_m)
  call get_param(param_file, mdl, "C_P", cs%Cp, &
                 "The heat capacity of sea water, approximated as a constant. "//&
                 "The default value is from the TEOS-10 definition of conservative temperature.", &
                 units="J kg-1 K-1", default=3991.86795711963, scale=us%J_kg_to_Q)
  call get_param(param_file, mdl, "RHO_0", cs%Rho_ocn, &
                 "The mean ocean density used with BOUSSINESQ true to "//&
                 "calculate accelerations and the mass for conservation "//&
                 "properties, or with BOUSSINSEQ false to convert some "//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
  call get_param(param_file, mdl, "C_P_ICE", cs%Cp_ice, &
                 "The heat capacity of ice.", units="J kg-1 K-1", scale=us%J_kg_to_Q, &
                 default=2.10e3)
  if (cs%constant_sea_level) cs%min_ocean_mass_float = dz_ocean_min_float*cs%Rho_ocn

  call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", cs%flux_factor, &
                 "Non-dimensional factor applied to shelf thermodynamic "//&
                 "fluxes.", units="none", default=1.0)

  call get_param(param_file, mdl, "KV_ICE", cs%kv_ice, &
                 "The viscosity of the ice.", &
                 units="m2 s-1", default=1.0e10, scale=us%Z_to_L**2*us%m_to_L**2*us%T_to_s)
  call get_param(param_file, mdl, "KV_MOLECULAR", cs%kv_molec, &
                 "The molecular kinimatic viscosity of sea water at the "//&
                 "freezing temperature.", units="m2 s-1", default=1.95e-6, scale=us%m2_s_to_Z2_T)
  call get_param(param_file, mdl, "ICE_SHELF_SALINITY", cs%Salin_ice, &
                 "The salinity of the ice inside the ice shelf.", units="psu", &
                 default=0.0)
  call get_param(param_file, mdl, "ICE_SHELF_TEMPERATURE", cs%Temp_ice, &
                 "The temperature at the center of the ice shelf.", &
                 units = "degC", default=-15.0)
  call get_param(param_file, mdl, "KD_SALT_MOLECULAR", cs%kd_molec_salt, &
                 "The molecular diffusivity of salt in sea water at the "//&
                 "freezing point.", units="m2 s-1", default=8.02e-10, scale=us%m2_s_to_Z2_T)
  call get_param(param_file, mdl, "KD_TEMP_MOLECULAR", cs%kd_molec_temp, &
                 "The molecular diffusivity of heat in sea water at the "//&
                 "freezing point.", units="m2 s-1", default=1.41e-7, scale=us%m2_s_to_Z2_T)
  call get_param(param_file, mdl, "DT_FORCING", cs%time_step, &
                 "The time step for changing forcing, coupling with other "//&
                 "components, or potentially writing certain diagnostics. "//&
                 "The default value is given by DT.", units="s", default=0.0)

  call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", col_thick_melt_thresh, &
                 "The minimum ocean column thickness where melting is allowed.", &
                 units="m", scale=us%m_to_Z, default=0.0)
  cs%col_mass_melt_threshold =  cs%Rho_ocn * col_thick_melt_thresh

  call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
                 "If true, read a file (given by TIDEAMP_FILE) containing "//&
                 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)

  call safe_alloc_ptr(cs%utide,isd,ied,jsd,jed)   ; cs%utide(:,:) = 0.0

  if (read_tideamp) then
    call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
                 "The path to the file containing the spatially varying "//&
                 "tidal amplitudes.", &
                 default="tideamp.nc")
    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
    inputdir = slasher(inputdir)
    tideamp_file = trim(inputdir) // trim(tideamp_file)
    call mom_read_data(tideamp_file, 'tideamp', cs%utide, g%domain, timelevel=1, scale=us%m_s_to_L_T)
  else
    call get_param(param_file, mdl, "UTIDE", utide, &
                 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
                 units="m s-1", default=0.0 , scale=us%m_s_to_L_T)
    cs%utide(:,:) = utide
  endif

  call eos_init(param_file, cs%eqn_of_state)

  !! new parameters that need to be in MOM_input

  if (cs%active_shelf_dynamics) then

    call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
                 "A typical density of ice.", units="kg m-3", default=917.0, scale=us%kg_m3_to_R)

    call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", cs%input_flux, &
                 "volume flux at upstream boundary", units="m2 s-1", default=0.)
    call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", cs%input_thickness, &
                 "flux thickness at upstream boundary", units="m", default=1000.)
  else
    ! This is here because of inconsistent defaults.  I don't know why.  RWH
    call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
                 "A typical density of ice.", units="kg m-3", default=900.0, scale=us%kg_m3_to_R)
  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)

  call get_param(param_file, mdl, "USTAR_SHELF_BG", cs%ustar_bg, &
                 "The minimum value of ustar under ice shelves.", &
                 units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
  call get_param(param_file, mdl, "CDRAG_SHELF", cdrag, &
       "CDRAG is the drag coefficient relating the magnitude of "//&
       "the velocity field to the surface stress.", units="nondim", &
       default=0.003)
  cs%cdrag = cdrag
  if (cs%ustar_bg <= 0.0) then
    call get_param(param_file, mdl, "DRAG_BG_VEL_SHELF", drag_bg_vel, &
                 "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
                 "LINEAR_DRAG) or an unresolved  velocity that is "//&
                 "combined with the resolved velocity to estimate the "//&
                 "velocity magnitude.", units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
    if (cs%cdrag*drag_bg_vel > 0.0) cs%ustar_bg = sqrt(cs%cdrag)*drag_bg_vel
  endif

  ! Allocate and initialize state variables to default values
  call ice_shelf_state_init(cs%ISS, cs%grid)
  iss => cs%ISS

  ! Allocate the arrays for passing ice-shelf data through the forcing type.
  if (.not. cs%solo_ice_sheet) then
    call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.")
     ! GMM: the following assures that water/heat fluxes are just allocated
     ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
     ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
    if (present(fluxes)) &
      call allocate_forcing_type(cs%ocn_grid, fluxes, ustar=.true., shelf=.true., &
                                 press=.true., water=cs%isthermo, heat=cs%isthermo)
    if (present(forces)) &
      call allocate_mech_forcing(cs%ocn_grid, forces, ustar=.true., shelf=.true., press=.true.)
  else
    call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
    if (present(fluxes)) &
      call allocate_forcing_type(g, fluxes, ustar=.true., shelf=.true., press=.true.)
    if (present(forces)) &
      call allocate_mech_forcing(g, forces, ustar=.true., shelf=.true., press=.true.)
  endif

  ! Set up the bottom depth, G%D either analytically or from file
  call mom_initialize_topography(dg%bathyT, g%max_depth, dg, param_file)
  call rescale_dyn_horgrid_bathymetry(dg, us%Z_to_m)

  ! Set up the Coriolis parameter, G%f, usually analytically.
  call mom_initialize_rotation(dg%CoriolisBu, dg, param_file, us)
  ! This copies grid elements, including bathyT and CoriolisBu from dG to CS%grid.
  call copy_dyngrid_to_mom_grid(dg, cs%grid, us)

  call destroy_dyn_horgrid(dg)

  ! Set up the restarts.
  call restart_init(param_file, cs%restart_CSp, "Shelf.res")
  call register_restart_field(iss%mass_shelf, "shelf_mass", .true., cs%restart_CSp, &
                              "Ice shelf mass", "kg m-2")
  call register_restart_field(iss%area_shelf_h, "shelf_area", .true., cs%restart_CSp, &
                              "Ice shelf area in cell", "m2")
  call register_restart_field(iss%h_shelf, "h_shelf", .true., cs%restart_CSp, &
                              "ice sheet/shelf thickness", "m")
  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., cs%restart_CSp, &
                              "Height unit conversion factor", "Z meter-1")
  call register_restart_field(us%m_to_L_restart, "m_to_L", .false., cs%restart_CSp, &
                              "Length unit conversion factor", "L meter-1")
  call register_restart_field(us%kg_m3_to_R_restart, "kg_m3_to_R", .false., cs%restart_CSp, &
                              "Density unit conversion factor", "R m3 kg-1")
  if (cs%active_shelf_dynamics) then
    call register_restart_field(iss%hmask, "h_mask", .true., cs%restart_CSp, &
                                "ice sheet/shelf thickness mask" ,"none")
  endif

  if (cs%active_shelf_dynamics) then
    ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics
    call register_ice_shelf_dyn_restarts(g, param_file, cs%dCS, cs%restart_CSp)
  endif

  !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file
  !if (.not. CS%solo_ice_sheet) then
  !  call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, &
  !                              "Friction velocity under ice shelves", "m s-1")
  !endif

  cs%restart_output_dir = dirs%restart_output_dir

  new_sim = .false.
  if ((dirs%input_filename(1:1) == 'n') .and. &
      (len_trim(dirs%input_filename) == 1)) new_sim = .true.

  if (cs%override_shelf_movement .and. cs%mass_from_file) then

    ! initialize the ids for reading shelf mass from a netCDF
    call initialize_shelf_mass(g, param_file, cs, iss)

    if (new_sim) then
      ! new simulation, initialize ice thickness as in the static case
      call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)

    ! next make sure mass is consistent with thickness
      do j=g%jsd,g%jed ; do i=g%isd,g%ied
        if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
          iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%density_ice
        endif
      enddo ; enddo

      if (cs%min_thickness_simple_calve > 0.0) &
        call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
                                           cs%min_thickness_simple_calve)
    endif
  endif

  if (cs%active_shelf_dynamics) then
    ! the only reason to initialize boundary conds is if the shelf is dynamic - MJH

    ! call initialize_ice_shelf_boundary ( CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
    !                                      CS%u_flux_bdry_val, CS%v_flux_bdry_val, &
    !                                      CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, &
    !                                      ISS%hmask, G, param_file)

  endif

  if (new_sim .and. (.not. (cs%override_shelf_movement .and. cs%mass_from_file))) then

    ! This model is initialized internally or from a file.
    call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)

    ! next make sure mass is consistent with thickness
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
        iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%density_ice
      endif
    enddo ; enddo

  ! else ! Previous block for new_sim=.T., this block restores the state.
  elseif (.not.new_sim) then
    ! This line calls a subroutine that reads the initial conditions from a restart file.
    call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.")
    call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
                       g, cs%restart_CSp)

    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
        iss%h_shelf(i,j) = z_rescale * iss%h_shelf(i,j)
      enddo ; enddo
    endif

    if ((us%m_to_Z_restart*us%kg_m3_to_R_restart /= 0.0) .and. &
        (us%m_to_Z*us%kg_m3_to_R /= us%m_to_Z_restart * us%kg_m3_to_R_restart)) then
      rz_rescale = us%m_to_Z*us%kg_m3_to_R / (us%m_to_Z_restart * us%kg_m3_to_R_restart)
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        iss%mass_shelf(i,j) = rz_rescale * iss%mass_shelf(i,j)
      enddo ; enddo
    endif

    if ((us%m_to_L_restart /= 0.0) .and. (us%m_to_L_restart /= us%m_to_L)) then
      l_rescale = us%m_to_L / us%m_to_L_restart
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        iss%area_shelf_h(i,j) = l_rescale**2 * iss%area_shelf_h(i,j)
      enddo ; enddo
    endif

  endif ! .not. new_sim

  cs%Time = time

  call cpu_clock_begin(id_clock_pass)
  call pass_var(iss%area_shelf_h, g%domain)
  call pass_var(iss%h_shelf, g%domain)
  call pass_var(iss%mass_shelf, g%domain)
  call pass_var(iss%hmask, g%domain)
  call pass_var(g%bathyT, g%domain)
  call cpu_clock_end(id_clock_pass)

  do j=jsd,jed ; do i=isd,ied
    if (iss%area_shelf_h(i,j) > g%areaT(i,j)) then
      call mom_error(warning,"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.")
      iss%area_shelf_h(i,j) = g%areaT(i,j)
    endif
  enddo ; enddo
  if (present(fluxes)) then ; do j=jsd,jed ; do i=isd,ied
    if (g%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = iss%area_shelf_h(i,j) / g%areaT(i,j)
  enddo ; enddo ; endif

  if (cs%debug) then
    call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", g%HI, haloshift=0)
  endif

  if (present(forces)) &
    call add_shelf_forces(g, us, cs, forces, do_shelf_area=.not.cs%solo_ice_sheet)

  if (present(fluxes)) call add_shelf_pressure(g, us, cs, fluxes)

  if (cs%active_shelf_dynamics .and. .not.cs%isthermo) then
    iss%water_flux(:,:) = 0.0
  endif

  if (shelf_mass_is_dynamic) &
    call initialize_ice_shelf_dyn(param_file, time, iss, cs%dCS, g, us, diag, new_sim, solo_ice_sheet_in)

  call fix_restart_unit_scaling(us)

  call get_param(param_file, mdl, "SAVE_INITIAL_CONDS", save_ic, &
                 "If true, save the ice shelf initial conditions.", &
                 default=.false.)
  if (save_ic) call get_param(param_file, mdl, "SHELF_IC_OUTPUT_FILE", ic_file,&
                 "The name-root of the output file for the ice shelf "//&
                 "initial conditions.", default="MOM_Shelf_IC")

  if (save_ic .and. .not.((dirs%input_filename(1:1) == 'r') .and. &
                          (len_trim(dirs%input_filename) == 1))) then
    call save_restart(dirs%output_directory, cs%Time, g, &
                      cs%restart_CSp, filename=ic_file)
  endif


  cs%id_area_shelf_h = register_diag_field('ocean_model', 'area_shelf_h', cs%diag%axesT1, cs%Time, &
     'Ice Shelf Area in cell', 'meter-2', conversion=us%L_to_m**2)
  cs%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', cs%diag%axesT1, cs%Time, &
     'mass of shelf', 'kg/m^2', conversion=us%RZ_to_kg_m2)
  cs%id_h_shelf = register_diag_field('ocean_model', 'h_shelf', cs%diag%axesT1, cs%Time, &
       'ice shelf thickness', 'm', conversion=us%Z_to_m)
  cs%id_mass_flux = register_diag_field('ocean_model', 'mass_flux', cs%diag%axesT1,&
     cs%Time, 'Total mass flux of freshwater across the ice-ocean interface.', &
     'kg/s', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2)

  if (cs%const_gamma) then ! use ISOMIP+ eq. with rho_fw = 1000. kg m-3
    meltrate_conversion = 86400.0*365.0*us%Z_to_m*us%s_to_T / (1000.0*us%kg_m3_to_R)
  else ! use original eq.
    meltrate_conversion = 86400.0*365.0*us%Z_to_m*us%s_to_T / cs%density_ice
  endif
  cs%id_melt = register_diag_field('ocean_model', 'melt', cs%diag%axesT1, cs%Time, &
     'Ice Shelf Melt Rate', 'm yr-1', conversion= meltrate_conversion)
  cs%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', cs%diag%axesT1, cs%Time, &
     'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', 'Celsius')
  cs%id_haline_driving = register_diag_field('ocean_model', 'haline_driving', cs%diag%axesT1, cs%Time, &
     'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'psu')
  cs%id_Sbdry = register_diag_field('ocean_model', 'sbdry', cs%diag%axesT1, cs%Time, &
     'salinity at the ice-ocean interface.', 'psu')
  cs%id_u_ml = register_diag_field('ocean_model', 'u_ml', cs%diag%axesCu1, cs%Time, &
     'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_v_ml = register_diag_field('ocean_model', 'v_ml', cs%diag%axesCv1, cs%Time, &
     'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', cs%diag%axesT1, cs%Time, &
     'Sub-shelf salinity exchange velocity', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
  cs%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', cs%diag%axesT1, cs%Time, &
     'Sub-shelf thermal exchange velocity', 'm s-1' , conversion=us%Z_to_m*us%s_to_T)
  cs%id_tfreeze = register_diag_field('ocean_model', 'tfreeze', cs%diag%axesT1, cs%Time, &
     'In Situ Freezing point at ice shelf interface', 'degC')
  cs%id_tfl_shelf = register_diag_field('ocean_model', 'tflux_shelf', cs%diag%axesT1, cs%Time, &
     'Heat conduction into ice shelf', 'W m-2', conversion=-us%QRZ_T_to_W_m2)
  cs%id_ustar_shelf = register_diag_field('ocean_model', 'ustar_shelf', cs%diag%axesT1, cs%Time, &
     'Fric vel under shelf', 'm/s', conversion=us%Z_to_m*us%s_to_T)
  if (cs%active_shelf_dynamics) then
    cs%id_h_mask = register_diag_field('ocean_model', 'h_mask', cs%diag%axesT1, cs%Time, &
       'ice shelf thickness mask', 'none')
  endif

  id_clock_shelf = cpu_clock_id('Ice shelf', grain=clock_component)
  id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=clock_routine)

end subroutine initialize_ice_shelf

!> Initializes shelf mass based on three options (file, zero and user)
subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim)

  type(ocean_grid_type), intent(in) :: G   !< The ocean's grid structure.
  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
  type(ice_shelf_CS),    pointer    :: CS !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
  logical,     optional, intent(in) :: new_sim !< If present and false, this run is being restarted

  integer :: i, j, is, ie, js, je
  logical :: read_shelf_area, new_sim_2
  character(len=240) :: config, inputdir, shelf_file, filename
  character(len=120) :: shelf_mass_var  ! The name of shelf mass in the file.
  character(len=120) :: shelf_area_var ! The name of shelf area in the file.
  character(len=40)  :: mdl = "MOM_ice_shelf"
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  new_sim_2 = .true. ; if (present(new_sim)) new_sim_2 = new_sim

  call get_param(param_file, mdl, "ICE_SHELF_CONFIG", config, &
                 "A string that specifies how the ice shelf is "//&
                 "initialized. Valid options include:\n"//&
                 " \tfile\t Read from a file.\n"//&
                 " \tzero\t Set shelf mass to 0 everywhere.\n"//&
                 " \tUSER\t Call USER_initialize_shelf_mass.\n", &
                 fail_if_missing=.true.)

  select case ( trim(config) )
    case ("file")

      call time_interp_external_init()

      call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
      inputdir = slasher(inputdir)

      call get_param(param_file, mdl, "SHELF_FILE", shelf_file, &
              "If DYNAMIC_SHELF_MASS = True, OVERRIDE_SHELF_MOVEMENT = True "//&
              "and ICE_SHELF_MASS_FROM_FILE = True, this is the file from "//&
              "which to read the shelf mass and area.", &
               default="shelf_mass.nc")
      call get_param(param_file, mdl, "SHELF_MASS_VAR", shelf_mass_var, &
                 "The variable in SHELF_FILE with the shelf mass.", &
                 default="shelf_mass")
      call get_param(param_file, mdl, "READ_SHELF_AREA", read_shelf_area, &
                 "If true, also read the area covered by ice-shelf from SHELF_FILE.", &
                 default=.false.)

      filename = trim(slasher(inputdir))//trim(shelf_file)
      call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename)

      cs%id_read_mass = init_external_field(filename, shelf_mass_var, &
                          domain=g%Domain%mpp_domain, verbose=cs%debug)

      if (read_shelf_area) then
         call get_param(param_file, mdl, "SHELF_AREA_VAR", shelf_area_var, &
                  "The variable in SHELF_FILE with the shelf area.", &
                  default="shelf_area")

         cs%id_read_area = init_external_field(filename,shelf_area_var, &
                             domain=g%Domain%mpp_domain)
      endif

      if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
           " initialize_shelf_mass: Unable to open "//trim(filename))

    case ("zero")
      do j=js,je ; do i=is,ie
        iss%mass_shelf(i,j) = 0.0
        iss%area_shelf_h(i,j) = 0.0
      enddo ; enddo

    case ("USER")
      call user_initialize_shelf_mass(iss%mass_shelf, iss%area_shelf_h, &
                   iss%h_shelf, iss%hmask, g, cs%US, cs%user_CS, param_file, new_sim_2)

    case default ;  call mom_error(fatal,"initialize_ice_shelf: "// &
      "Unrecognized ice shelf setup "//trim(config))
  end select

end subroutine initialize_shelf_mass

!> Updates the ice shelf mass using data from a file.
subroutine update_shelf_mass(G, US, CS, ISS, Time)
  type(ocean_grid_type), intent(inout) :: G   !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
  type(ice_shelf_CS),    intent(in)    :: CS  !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
  type(time_type),       intent(in)    :: Time !< The current model time

  ! local variables
  integer :: i, j, is, ie, js, je
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  call time_interp_external(cs%id_read_mass, time, iss%mass_shelf)
  ! This should only be done if time_interp_external did an update.
  do j=js,je ; do i=is,ie
    iss%mass_shelf(i,j) = us%kg_m3_to_R*us%m_to_Z * iss%mass_shelf(i,j) ! Rescale after time_interp
  enddo ; enddo

  do j=js,je ; do i=is,ie
    iss%area_shelf_h(i,j) = 0.0
    iss%hmask(i,j) = 0.
    if (iss%mass_shelf(i,j) > 0.0) then
      iss%area_shelf_h(i,j) = g%areaT(i,j)
      iss%h_shelf(i,j) = iss%mass_shelf(i,j) / cs%density_ice
      iss%hmask(i,j) = 1.
    endif
  enddo ; enddo

  !call USER_update_shelf_mass(ISS%mass_shelf, ISS%area_shelf_h, ISS%h_shelf, &
  !                            ISS%hmask, CS%grid, CS%user_CS, Time, .true.)

  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, halo=0)
  endif

  call pass_var(iss%area_shelf_h, g%domain)
  call pass_var(iss%h_shelf, g%domain)
  call pass_var(iss%hmask, g%domain)
  call pass_var(iss%mass_shelf, g%domain)

end subroutine update_shelf_mass

!> Save the ice shelf restart file
subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix)
  type(ice_shelf_cs),         pointer    :: cs !< ice shelf control structure
  type(time_type),            intent(in) :: time !< model time at this call
  character(len=*), optional, intent(in) :: directory !< An optional directory into which to write
                                               !! these restart files.
  logical,          optional, intent(in) :: time_stamped !< f true, the restart file names include
                                               !! a unique time stamp.  The default is false.
  character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a
                                               !! time-stamp) to append to the restart file names.
  ! local variables
  type(ocean_grid_type), pointer :: g => null()
  character(len=200) :: restart_dir

  g => cs%grid

  if (present(directory)) then ; restart_dir = directory
  else ; restart_dir = cs%restart_output_dir ; endif

  call save_restart(restart_dir, time, cs%grid, cs%restart_CSp, time_stamped)

end subroutine ice_shelf_save_restart

!> Deallocates all memory associated with this module
subroutine ice_shelf_end(CS)
  type(ice_shelf_cs), pointer   :: cs !< A pointer to the ice shelf control structure

  if (.not.associated(cs)) return

  call ice_shelf_state_end(cs%ISS)

  if (cs%active_shelf_dynamics) call ice_shelf_dyn_end(cs%dCS)

  deallocate(cs)

end subroutine ice_shelf_end

!> This routine is for stepping a stand-alone ice shelf model without an ocean.
subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in)
  type(ice_shelf_cs), pointer    :: cs      !< A pointer to the ice shelf control structure
  type(time_type), intent(in)    :: time_interval !< The time interval for this update [s].
  integer,         intent(inout) :: nsteps  !< The running number of ice shelf steps.
  type(time_type), intent(inout) :: time    !< The current model time
  real,  optional, intent(in)    :: min_time_step_in !< The minimum permitted time step [T ~> s].

  type(ocean_grid_type), pointer :: g => null()  ! A pointer to the ocean's grid structure
  type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
                                                 ! various unit conversion factors
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                          !! the ice-shelf state
  real :: remaining_time    ! The remaining time in this call [T ~> s]
  real :: time_step         ! The internal time step during this call [T ~> s]
  real :: min_time_step     ! The minimal required timestep that would indicate a fatal problem [T ~> s]
  character(len=240) :: mesg
  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
  logical :: coupled_gl     ! If true the grouding line position is determined based on
                            ! coupled ice-ocean dynamics.
  integer :: is, iec, js, jec, i, j

  g => cs%grid
  us => cs%US
  iss => cs%ISS
  is = g%isc ; iec = g%iec ; js = g%jsc ; jec = g%jec

  remaining_time = us%s_to_T*time_type_to_real(time_interval)

  if (present (min_time_step_in)) then
    min_time_step = min_time_step_in
  else
    min_time_step = 1000.0*us%s_to_T ! At 1 km resolution this would imply ice is moving at ~1 meter per second
  endif

  write (mesg,*) "TIME in ice shelf call, yrs: ", time_type_to_real(time)/(365. * 86400.)
  call mom_mesg("solo_step_ice_shelf: "//mesg, 5)

  do while (remaining_time > 0.0)
    nsteps = nsteps+1

    ! If time_interval is not too long, this is unnecessary.
    time_step = min(ice_time_step_cfl(cs%dCS, iss, g), remaining_time)

    write (mesg,*) "Ice model timestep = ", us%T_to_s*time_step, " seconds"
    if ((time_step < min_time_step) .and. (time_step < remaining_time))  then
      call mom_error(fatal, "MOM_ice_shelf:solo_step_ice_shelf: abnormally small timestep "//mesg)
    else
      call mom_mesg("solo_step_ice_shelf: "//mesg, 5)
    endif

    remaining_time = remaining_time - time_step

    ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much.
    ! Do not update the velocities if the last step is very short.
    update_ice_vel = ((time_step > min_time_step) .or. (remaining_time > 0.0))
    coupled_gl = .false.

    call update_ice_shelf(cs%dCS, iss, g, us, time_step, time, must_update_vel=update_ice_vel)

    call enable_averages(time_step, time, cs%diag)
    if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
    if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf, iss%h_shelf, cs%diag)
    if (cs%id_h_mask > 0) call post_data(cs%id_h_mask, iss%hmask, cs%diag)
    call disable_averaging(cs%diag)

  enddo

end subroutine solo_step_ice_shelf

!> \namespace mom_ice_shelf
!!
!! \section section_ICE_SHELF
!!
!! This module implements the thermodynamic aspects of ocean/ice-shelf
!! inter-actions using the MOM framework and coding style.
!!
!! Derived from code by Chris Little, early 2010.
!!
!!   The ice-sheet dynamics subroutines do the following:
!!  initialize_shelf_mass - Initializes the ice shelf mass distribution.
!!      - Initializes h_shelf, h_mask, area_shelf_h
!!      - CURRENTLY: initializes mass_shelf as well, but this is unnecessary, as mass_shelf is initialized based on
!!             h_shelf and density_ice immediately afterwards. Possibly subroutine should be renamed
!!  update_shelf_mass - updates ice shelf mass via netCDF file
!!                      USER_update_shelf_mass (TODO).
!!    solo_step_ice_shelf - called only in ice-only mode.
!!    shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. currently mass_shelf is
!! updated immediately after ice_shelf_advect in fully dynamic mode.
!!
!!   NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement,
!! for subroutines in the velocity solve, and for thickness boundary conditions (this last one may be removed).
!! in other words, interfering with its updates will have implications you might not expect.
!!
!!  Overall issues: Many variables need better documentation and units and the
!!                  subgrid on which they are discretized.
!!
!! \subsection section_ICE_SHELF_equations ICE_SHELF equations
!!
!! The three fundamental equations are:
!! Heat flux
!! \f[ \qquad \rho_w  C_{pw} \gamma_T (T_w - T_b) = \rho_i  \dot{m}  L_f \f]
!! Salt flux
!! \f[  \qquad \rho_w \gamma_s (S_w - S_b) =  \rho_i \dot{m} S_b \f]
!! Freezing temperature
!! \f[  \qquad T_b = a S_b + b + c P \f]
!!
!! where ....
!!
!! \subsection section_ICE_SHELF_references References
!!
!! Asay-Davis, Xylar S., Stephen L. Cornford, Benjamin K. Galton-Fenzi, Rupert M. Gladstone, G. Hilmar Gudmundsson,
!! David M. Holland, Paul R. Holland, and Daniel F. Martin. Experimental design for three interrelated marine ice sheet
!! and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MISOMIP v. 1 (MISOMIP1).
!! Geoscientific Model Development 9, no. 7 (2016): 2471.
!!
!! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 1.
!!  Model description and behavior. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
!!
!! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 2.
!! Sensitivity to external forcings. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
!!
!! Holland, David M., and Adrian Jenkins. Modeling thermodynamic ice-ocean interactions at the base of an ice shelf.
!! Journal of Physical Oceanography 29.8 (1999): 1787-1800.

end module mom_ice_shelf