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_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_component, clock_routine
use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid
use mom_diag_mediator, only : diag_ctrl, time_type, enable_averaging, disable_averaging
use mom_domains, only : mom_domains_init, clone_mom_domain
use mom_domains, only : pass_var, pass_vector, to_all, cgrid_ne, bgrid_ne, corner
use mom_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, time_type_to_real, real_to_time
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
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, sum_across_pes
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
use time_manager_mod, only : print_time
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_time_step, 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()  !< tidal velocity [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 [m s-2]
  real :: cp           !< The heat capacity of sea water [J kg-1 degC-1].
  real :: rho0         !< A reference ocean density [kg m-3].
  real :: cp_ice       !< The heat capacity of fresh ice [J kg-1 degC-1].
  real :: gamma_t      !< The (fixed) turbulent exchange velocity in the
                       !< 2-equation formulation [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 [m2 s-1].
  real :: density_ice  !< A typical density of ice [kg m-3].
  real :: rho_ice      !< Nominal ice density [kg m-2 Z-1 ~> kg m-3].
  real :: kv_molec     !< The molecular kinematic viscosity of sea water [m2 s-1].
  real :: kd_molec_salt!< The molecular diffusivity of salt [m2 s-1].
  real :: kd_molec_temp!< The molecular diffusivity of heat [m2 s-1].
  real :: lat_fusion   !< The latent heat of fusion [J kg-1].
  real :: gamma_t_3eq  !<  Nondimensional heat-transfer coefficient, used in the 3Eq. formulation
                       !<  This number should be specified by the user.
  real :: col_thick_melt_threshold !< if the mixed layer is below this threshold, melt rate
  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)
  real :: density_ocean_avg !< this does not affect ocean circulation OR thermodynamics
                            !! it is to estimate the gravitational driving force at the
                            !! shelf front (until we think of a better way to do it,
                            !! but any difference will be negligible)
  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 :: find_salt_root              !< If true, if true find Sbdry using a quadratic eq.
  logical :: constant_sea_level          !< if true, apply an evaporative, heat and salt
                                         !! fluxes. It will avoid large increase in sea level.
  real    :: cutoff_depth                !< depth above which melt is set to zero (>= 0).
  real    :: lambda1                     !< liquidus coeff., Needed if find_salt_root = true
  real    :: lambda2                     !< liquidus coeff., Needed if find_salt_root = true
  real    :: lambda3                     !< liquidus coeff., Needed if find_salt_root = true
  !>@{ 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(state, fluxes, Time, time_step, CS, forces)
  type(surface),         intent(inout) :: state !< structure containing fields that
                                                !!describe the surface state of the ocean
  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
 
  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 [kg m-3].
    dr0_dt, &  !< Partial derivative of the mixed layer density
               !< with temperature [kg m-3 degC-1].
    dr0_ds, &  !< Partial derivative of the mixed layer density
               !< with salinity [kg m-3 ppt-1].
    p_int      !< The pressure at the ice-ocean interface [Pa].
 
  real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
    exch_vel_t, &  !< Sub-shelf thermal exchange velocity [m s-1]
    exch_vel_s     !< Sub-shelf salt exchange velocity [m s-1]
 
  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
    mass_flux  !< total mass flux of freshwater across
  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. (Was 1/8. Why?)
  real, parameter :: rc    = 0.20     ! critical flux Richardson number.
  real :: i_zeta_n !< The inverse of ZETA_N.
  real :: lf, i_lf !< Latent Heat of fusion [J kg-1] and its inverse.
  real :: i_vk     !< The inverse of VK.
  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, s_a, s_b, s_c  ! use to find salt roots
  real :: ds_it    !< The interface salinity change during an iteration [ppt].
  real :: hbl_neut !< The neutral boundary layer thickness [m].
  real :: hbl_neut_h_molec !< The ratio of the neutral boundary layer thickness
                   !! to the molecular boundary layer thickness [nondim].
  !### THESE ARE CURRENTLY POSITIVE UPWARD.
  real :: wt_flux !< The vertical flux of heat just inside the ocean [degC m s-1].
  real :: wb_flux !< The vertical flux of heat just inside the ocean [m2 s-3].
  real :: db_ds   !< The derivative of buoyancy with salinity [m s-2 ppt-1].
  real :: db_dt   !< The derivative of buoyancy with temperature [m s-2 degC-1].
  real :: i_n_star, n_star_term, absf
  real :: dins_dwb !< The partial derivative of I_n_star with wB_flux, in ???.
  real :: dt_ustar, ds_ustar
  real :: ustar_h
  real :: gam_turb
  real :: gam_mol_t, gam_mol_s
  real :: rhocp
  real :: i_rholf
  real :: ln_neut
  real :: mass_exch
  real :: sb_min, sb_max
  real :: ds_min, ds_max
  ! Variables used in iterating for wB_flux.
  real :: wb_flux_new, dwb, ddwb_dwb_in
  real :: i_gam_t, i_gam_s, dg_dwb, idens
  real :: u_at_h, v_at_h, isqrt2
  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 :: i, j, is, ie, js, je, ied, jed, it1, it3
  real, parameter :: rho_fw = 1000.0 ! fresh water density
 
  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
  lf = cs%Lat_fusion
  i_rholf = 1.0/(cs%Rho0*lf)
  i_lf = 1.0 / lf
  sc = cs%kv_molec/cs%kd_molec_salt
  pr = cs%kv_molec/cs%kd_molec_temp
  i_vk = 1.0/vk
  rhocp = cs%Rho0 * cs%Cp
  isqrt2 = 1.0/sqrt(2.0)
 
  !first calculate molecular component
  gam_mol_t = 12.5 * (pr**c2_3) - 6
  gam_mol_s = 12.5 * (sc**c2_3) - 6
 
  idens = 1.0/cs%density_ocean_avg
 
  ! 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(:,:) = 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, 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(state%sst, "sst before apply melting", g%HI, haloshift=0)
    call hchksum(state%sss, "sss before apply melting", g%HI, haloshift=0)
    call hchksum(state%u, "u_ml before apply melting", g%HI, haloshift=0)
    call hchksum(state%v, "v_ml before apply melting", g%HI, haloshift=0)
    call hchksum(state%ocean_mass, "ocean_mass before apply melting", g%HI, haloshift=0)
  endif
 
  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(state%sst(:,j),state%sss(:,j), p_int, &
             rhoml(:), is, ie-is+1, cs%eqn_of_state)
    call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, &
             dr0_dt, dr0_ds, is, ie-is+1, cs%eqn_of_state)
 
    do i=is,ie
      ! set ustar_shelf to zero. This is necessary if shelf_mass_is_dynamic
      ! but it won't make a difference otherwise.
      fluxes%ustar_shelf(i,j)= 0.0
 
      ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells
      !       propose instead to allow where Hml > [some threshold]
 
      if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
          (iss%area_shelf_h(i,j) > 0.0) .and. &
          (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) 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.
          ! ### SHOULD I SET USTAR_SHELF YET?
 
          u_at_h = state%u(i,j)
          v_at_h = state%v(i,j)
 
          !### I think that CS%utide**1 should be CS%utide**2
          fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%m_to_Z*us%T_to_s * &
              sqrt(cs%cdrag*((u_at_h**2 + v_at_h**2) + cs%utide(i,j)**1)))
 
          ustar_h = us%Z_to_m*us%s_to_T*fluxes%ustar_shelf(i,j)
 
          if (associated(state%taux_shelf) .and. associated(state%tauy_shelf)) then
            state%taux_shelf(i,j) = ustar_h*ustar_h*cs%Rho0*isqrt2
            state%tauy_shelf(i,j) = state%taux_shelf(i,j)
          endif
 
          ! 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*us%s_to_T*((abs(us%s_to_T*g%CoriolisBu(i,j)) + abs(us%s_to_T*g%CoriolisBu(i-1,j-1))) + &
                                 (abs(us%s_to_T*g%CoriolisBu(i,j-1)) + abs(us%s_to_T*g%CoriolisBu(i-1,j))))
          if (absf*state%Hml(i,j) <= vk*ustar_h) then ; hbl_neut = 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 = (cs%g_Earth / rhoml(i)) * dr0_ds(i)
          db_dt = (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
            ! read liquidus parameters
 
            s_a = cs%lambda1 * cs%Gamma_T_3EQ * cs%Cp
!            S_b = -CS%Gamma_T_3EQ*(CS%lambda2-CS%lambda3*p_int(i)-state%sst(i,j)) &
!               -LF*CS%Gamma_T_3EQ/35.0
 
            s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%lambda2+cs%lambda3*p_int(i)- &
                  state%sst(i,j))-lf*cs%Gamma_T_3EQ/35.0
            s_c = lf*(cs%Gamma_T_3EQ/35.0)*state%sss(i,j)
 
            !### Depending on the sign of S_b, one of these will be inaccurate!
            sbdry1 = (-s_b + sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
            sbdry2 = (-s_b - sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
            sbdry(i,j) = max(sbdry1, sbdry2)
            ! Safety check
            if (sbdry(i,j) < 0.) then
              write(mesg,*) 'state%sss(i,j) = ',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) = 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)
 
            dt_ustar = (state%sst(i,j) - iss%tfreeze(i,j)) * ustar_h
            ds_ustar = (state%sss(i,j) - sbdry(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_T_3EQ/35.
            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_T_3EQ/35.
                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 dwB = 0.0
                dwb = wb_flux_new - wb_flux
                if (abs(wb_flux_new - wb_flux) < &
                    1e-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 - dwb / 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 = rho_ice * CS%CP_Ice * K_ice * dT/dz (at interface)
    !vertical adv/diff as in H+J 199, eqs (31) & (26)...
    !  dT/dz ~= min( (lprec/(rho_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, so zero ice 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)) / (LF + cp_ice * (T_freeze-CS%Temp_Ice))
                iss%water_flux(i,j) = iss%tflux_ocn(i,j) / &
                      (lf + cs%CP_Ice * (iss%tfreeze(i,j) - cs%Temp_Ice))
 
                iss%tflux_shelf(i,j) = iss%tflux_ocn(i,j) - lf*iss%water_flux(i,j)
              endif
 
            endif
            !other options: dTi/dz linear through shelf
            !    dTi_dz = (CS%Temp_Ice - ISS%tfreeze(i,j))/G%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%Rho0
              sbdry_it = (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) < 1e-4*(0.5*(state%sss(i,j) + sbdry(i,j) + 1.e-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(state%sss(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state)
 
          exch_vel_t(i,j) = cs%gamma_t
          iss%tflux_ocn(i,j) = rhocp * exch_vel_t(i,j) * (state%sst(i,j) - iss%tfreeze(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
      else !not shelf
        iss%tflux_ocn(i,j) = 0.0
      endif
 
!      haline_driving(:,:) = state%sss(i,j) - Sbdry(i,j)
 
    enddo ! i-loop
  enddo ! j-loop
 
  ! ISS%water_flux = net liquid water into the ocean ( kg/(m^2 s) )
  ! We want melt in m/year
  if (cs%const_gamma) then ! use ISOMIP+ eq. with rho_fw
    fluxes%iceshelf_melt = iss%water_flux  * (86400.0*365.0/rho_fw) * cs%flux_factor
  else ! use original eq.
    fluxes%iceshelf_melt = iss%water_flux  * (86400.0*365.0/cs%density_ice) * cs%flux_factor
  endif
 
  do j=js,je ; do i=is,ie
    if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
        (iss%area_shelf_h(i,j) > 0.0) .and. &
        (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) then
 
      ! Set melt to zero above a cutoff pressure
      ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip
      ! test case.
      if ((cs%g_Earth * iss%mass_shelf(i,j)) < cs%Rho0*cs%cutoff_depth* &
         cs%g_Earth) 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%Rho0 * exch_vel_s(i,j))
 
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!!
      !1)Check if haline_driving computed above is consistent with
      ! haline_driving = state%sss - Sbdry
      !if (fluxes%iceshelf_melt(i,j) /= 0.0) then
      !   if (haline_driving(i,j) /= (state%sss(i,j) - Sbdry(i,j))) then
      !     write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',haline_driving(i,j), &
      !                   (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
    endif ! area_shelf_h
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
  enddo ; enddo ! i- and j-loops
 
  ! mass flux [kg s-1], part of ISOMIP diags.
  mass_flux(:,:) = 0.0
  mass_flux(:,:) = iss%water_flux(:,:) * iss%area_shelf_h(:,:)
 
  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, time_step, fluxes, cs%rho_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)
    endif
  endif
 
  if (cs%debug) call mom_forcing_chksum("Before add shelf flux", fluxes, g, cs%US, haloshift=0)
 
  call add_shelf_flux(g, cs, 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, time_step, time, 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, (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, state%u, cs%diag)
  if (cs%id_v_ml > 0) call post_data(cs%id_v_ml, 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, 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, time_step, fluxes, rho_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
  real,                  intent(in)    :: time_step !< The time step for this update [s].
  type(forcing),         intent(inout) :: fluxes !< structure containing pointers to any possible
                                                 !! thermodynamic or mass-flux forcing fields.
  real,                  intent(in)    :: rho_ice !< The density of ice-shelf ice [kg m-2 Z-1 ~> kg m-3].
  logical,     optional, intent(in)    :: debug !< If present and true, write chksums
 
  ! locals
  real :: I_rho_ice
  integer :: i, j
 
  i_rho_ice = 1.0 / rho_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) / rho_ice * time_step < iss%h_shelf(i,j)) then
        iss%h_shelf(i,j) = iss%h_shelf(i,j) - iss%water_flux(i,j) / rho_ice * time_step
      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
    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)
 
  !### combine this with the loops above.
  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)*rho_ice
    endif
  enddo ; enddo
 
  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, CS, forces, do_shelf_area)
  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  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 [m5 kg-1 s-1].
  real :: press_ice  ! The pressure of the ice shelf per unit area of ocean (not ice) [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
 
  !### Consider working over a smaller array range.
  do j=jsd,jed ; do i=isd,ied
    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.)
    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, CS, fluxes)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
  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) [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, CS, state, fluxes)
  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(ice_shelf_cs),    pointer       :: cs   !< This module's control structure.
  type(surface),         intent(inout) :: state!< Surface ocean state
  type(forcing),         intent(inout) :: fluxes  !< A structure of surface fluxes that may be used/updated.
 
  ! local variables
  real :: irho0           !< The inverse of the mean density [m3 kg-1].
  real :: frac_area       !< The fractional area covered by the ice shelf [nondim].
  real :: shelf_mass0     !< Total ice shelf mass at previous time (Time-dt).
  real :: shelf_mass1     !< Total ice shelf mass at current time (Time).
  real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1]
  real :: taux2, tauy2    !< The squared surface stresses [Pa].
  real :: press_ice       !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
  real :: asu1, asu2      !< Ocean areas covered by ice shelves at neighboring u-
  real :: asv1, asv2      !< and v-points [m2].
  real :: fraz            !< refreezing rate [kg m-2 s-1]
  real :: mean_melt_flux  !< spatial mean melt flux [kg s-1] or [kg m-2 s-1] at various points in the code.
  real :: sponge_area     !< total area of sponge region [m2]
  real :: t0              !< The previous time (Time-dt) [s].
  type(time_type) :: time0!< The previous time (Time-dt)
  real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
                          !! at at previous time (Time-dt) [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
                          !! at at previous time (Time-dt)
  real, dimension(SZDI_(G),SZDJ_(G))  :: last_area_shelf_h !< Ice shelf area [m2]
                          !! at at previous time (Time-dt)
  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
                                          !! the ice-shelf state
 
  real :: kv_rho_ice ! The viscosity of ice divided by its density [m5 kg-1 s-1]
  real, parameter :: rho_fw = 1000.0 ! fresh water density
  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, 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 (associated(state%taux_shelf) .and. associated(state%tauy_shelf)) then
      call uvchksum("tau[xy]_shelf", state%taux_shelf, state%tauy_shelf, &
                    g%HI, haloshift=0)
    endif
  endif
 
  if (associated(state%taux_shelf) .and. associated(state%tauy_shelf)) then
    call pass_vector(state%taux_shelf, state%tauy_shelf, g%domain, to_all, cgrid_ne)
  endif
  ! GMM: melting is computed using ustar_shelf (and not ustar), which has already
  ! been passed, I so believe we do not need to update fluxes%ustar.
!  Irho0 = 1.0 / CS%Rho0
!  do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then
    ! ### THIS SHOULD BE AN AREA WEIGHTED AVERAGE OF THE ustar_shelf POINTS.
    ! taux2 = 0.0 ; tauy2 = 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))
    ! if ((asu1 + asu2 > 0.0) .and. associated(state%taux_shelf)) &
    !   taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + &
    !            asu2 * state%taux_shelf(I,j)**2  ) / (asu1 + asu2)
    ! if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) &
    !   tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + &
    !            asv2 * state%tauy_shelf(i,J)**2  ) / (asv1 + asv2)
 
    ! fluxes%ustar(i,j) = MAX(CS%ustar_bg, US%m_to_Z*US%T_to_s*sqrt(Irho0 * sqrt(taux2 + tauy2)))
!  endif ; enddo ; enddo
 
  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) = iss%area_shelf_h(i,j) * g%IareaT(i,j)
    enddo ; enddo
  endif
 
  do j=js,je ; do i=is,ie ; if (iss%area_shelf_h(i,j) > 0.0) then
    frac_area = fluxes%frac_shelf_h(i,j)  !### Should this be 1-frac_shelf_h?
    if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
    if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = 0.0
    if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = 0.0
    if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = 0.0
    if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = 0.0
    if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
    if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
    if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
    if (associated(fluxes%lprec)) then
      if (iss%water_flux(i,j) > 0.0) then
        fluxes%lprec(i,j) =  frac_area*iss%water_flux(i,j)*cs%flux_factor
      else
        fluxes%lprec(i,j) = 0.0
        fluxes%evap(i,j) = frac_area*iss%water_flux(i,j)*cs%flux_factor
      endif
    endif
 
    if (associated(fluxes%sens)) &
      fluxes%sens(i,j) = -frac_area*iss%tflux_ocn(i,j)*cs%flux_factor
    if (associated(fluxes%salt_flux)) &
      fluxes%salt_flux(i,j) = frac_area * iss%salt_flux(i,j)*cs%flux_factor
  endif ; enddo ; enddo
 
  ! keep sea level constant by removing mass in 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
    !### This code has lots of problems with hard coded constants and the use of
    !### of non-reproducing sums.  It needs to be refactored. -RWH
 
    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
 
    mean_melt_flux = 0.0; sponge_area = 0.0
    do j=js,je ; do i=is,ie
      frac_area = fluxes%frac_shelf_h(i,j)
      if (frac_area > 0.0) &
        mean_melt_flux = mean_melt_flux + (iss%water_flux(i,j)) * iss%area_shelf_h(i,j)
 
      !### These hard-coded limits need to be corrected.  They are inappropriate here.
      if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
        sponge_area = sponge_area + g%areaT(i,j)
      endif
    enddo ; enddo
 
    ! take into account changes in mass (or thickness) when imposing ice shelf mass
    if (cs%override_shelf_movement .and. cs%mass_from_file) then
      t0 = time_type_to_real(cs%Time) - cs%time_step
 
      ! just compute changes in mass after first time step
      if (t0>0.0) then
        time0 = real_to_time(t0)
        last_hmask(:,:) = iss%hmask(:,:) ; last_area_shelf_h(:,:) = iss%area_shelf_h(:,:)
        call time_interp_external(cs%id_read_mass, time0, last_mass_shelf)
        last_h_shelf(:,:) = last_mass_shelf(:,:) / cs%rho_ice
 
        ! 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)
          ! convert to mass again
          last_mass_shelf(:,:) = last_h_shelf(:,:) * cs%rho_ice
        endif
 
        shelf_mass0 = 0.0; shelf_mass1 = 0.0
        ! get total ice shelf mass at (Time-dt) and (Time), in kg
        do j=js,je ; do i=is,ie
          ! just floating shelf (0.1 is a threshold for min ocean thickness)
          if (((1.0/cs%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. &
              (iss%area_shelf_h(i,j) > 0.0)) then
            shelf_mass0 = shelf_mass0 + (last_mass_shelf(i,j) * iss%area_shelf_h(i,j))
            shelf_mass1 = shelf_mass1 + (iss%mass_shelf(i,j) * iss%area_shelf_h(i,j))
          endif
        enddo ; enddo
        call sum_across_pes(shelf_mass0); call sum_across_pes(shelf_mass1)
        delta_mass_shelf = (shelf_mass1 - shelf_mass0)/cs%time_step
!          delta_mass_shelf = (shelf_mass1 - shelf_mass0)* &
!                         (rho_fw/CS%density_ice)/CS%time_step
!        write(mesg,*)'delta_mass_shelf = ',delta_mass_shelf
!        call MOM_mesg(mesg,5)
      else! first time step
        delta_mass_shelf = 0.0
      endif
    else ! ice shelf mass does not change
      delta_mass_shelf = 0.0
    endif
 
    call sum_across_pes(mean_melt_flux)
    call sum_across_pes(sponge_area)
 
    ! average total melt flux over sponge area
    mean_melt_flux = (mean_melt_flux+delta_mass_shelf) / sponge_area !kg/(m^2 s)
 
    ! apply fluxes
    do j=js,je ; do i=is,ie
      ! Note the following is hard coded for ISOMIP
      if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
        fluxes%vprec(i,j) = -mean_melt_flux * cs%density_ice/1000. ! evap is negative
        fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0 ! W /m^2
        fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3 ! kg (salt)/(m^2 s)
      endif
    enddo ; enddo
 
    if (cs%debug) then
      write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, 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 reastart fole to the internal representation in this run.
  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
  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)
 
  call create_dyn_horgrid(dg, cs%grid%HI)
  call clone_mom_domain(cs%grid%Domain, dg%Domain)
 
  call set_grid_metrics(dg, param_file)
  ! 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%Time = time ! ### This might not be in the right place?
  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%Lat_fusion = 3.34e5
  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, "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.", default=0.0)
  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, "ISOMIP_S_SUR_SPONGE", &
                cs%S0, "Surface salinity in the resoring region.", &
                default=33.8, do_not_log=.true.)
 
  call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", &
                cs%T0, "Surface temperature in the resoring region.", &
                default=-1.9, 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 salt-transfer coefficient is then computed "//&
                 " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
  if (cs%const_gamma) call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
                 "Nondimensional heat-transfer coefficient.",default=2.2e-2, &
                  units="nondim.", fail_if_missing=.true.)
 
  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%threeeq) &
    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.)
  if (cs%find_salt_root) then ! read liquidus coeffs.
     call get_param(param_file, mdl, "TFREEZE_S0_P0",cs%lambda1, &
                 "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%lambda1, &
                 "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%lambda3, &
                 "this is the derivative of the freezing potential "//&
                 "temperature with pressure.", &
                 units="degC Pa-1", default=0.0, do_not_log=.true.)
 
  endif
 
  if (.not.cs%threeeq) &
    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", fail_if_missing=.true.)
 
  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80)
  call get_param(param_file, mdl, "C_P", cs%Cp, &
                 "The heat capacity of sea water.", units="J kg-1 K-1", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
                 "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) !### MAKE THIS A SEPARATE PARAMETER.
  call get_param(param_file, mdl, "C_P_ICE", cs%Cp_ice, &
                 "The heat capacity of ice.", units="J kg-1 K-1", &
                 default=2.10e3)
 
  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)
  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)
  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)
  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)
  call get_param(param_file, mdl, "RHO_0", cs%density_ocean_avg, &
                 "avg ocean density used in floatation cond", &
                 units="kg m-3", default=1035.)
  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", cs%col_thick_melt_threshold, &
                 "The minimum ML thickness where melting is allowed.", units="m", &
                 default=0.0)
 
  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)
  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)
    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)
 
    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)
  endif
  cs%rho_ice = cs%density_ice*us%Z_to_m
  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 sheves.", &
                 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)
 
  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")
  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  !### Consider adding an ice shelf dynamics switch.
    ! 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")
  !  call register_restart_field(fluxes%iceshelf_melt, "iceshelf_melt", .false., CS%restart_CSp, &
  !                              "Ice Shelf Melt Rate", "m year-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%rho_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%rho_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
 
  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, cs, forces, do_shelf_area=.not.cs%solo_ice_sheet)
 
  if (present(fluxes)) call add_shelf_pressure(g, 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')
  cs%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', cs%diag%axesT1, cs%Time, &
     'mass of shelf', 'kg/m^2')
  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')
  cs%id_melt = register_diag_field('ocean_model', 'melt', cs%diag%axesT1, cs%Time, &
     'Ice Shelf Melt Rate', 'm yr-1')
  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')
  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')
  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')
  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')
  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')
  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, CS, ISS, Time)
  type(ocean_grid_type), intent(inout) :: G   !< The ocean's grid structure.
  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)
 
  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%rho_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)
  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_time_step(CS, time_step, nsteps, Time, min_time_step_in)
  type(ice_shelf_cs), pointer    :: cs !< A pointer to the ice shelf control structure
  real,            intent(in)    :: time_step !< 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 [s].
 
  type(ocean_grid_type), pointer :: g => null()
  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
  integer :: is, iec, js, jec, i, j
  real :: time_step_remain
  real :: time_step_int, min_time_step
  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.
 
  g => cs%grid
  us => cs%US
  iss => cs%ISS
  is = g%isc ; iec = g%iec ; js = g%jsc ; jec = g%jec
 
  time_step_remain = time_step
  if (present (min_time_step_in)) then
    min_time_step = min_time_step_in
  else
    min_time_step = 1000.0 ! This is in seconds - at 1 km resolution it 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_time_step: "//mesg)
 
  do while (time_step_remain > 0.0)
    nsteps = nsteps+1
 
    ! If time_step is not too long, this is unnecessary.
    time_step_int = min(ice_time_step_cfl(cs%dCS, iss, g), time_step)
 
    write (mesg,*) "Ice model timestep = ", time_step_int, " seconds"
    if (time_step_int < min_time_step) then
      call mom_error(fatal, "MOM_ice_shelf:solo_time_step: abnormally small timestep "//mesg)
    else
      call mom_mesg("solo_time_step: "//mesg)
    endif
 
    if (time_step_int >= time_step_remain) then
      time_step_int = time_step_remain
      time_step_remain = 0.0
    else
      time_step_remain = time_step_remain - time_step_int
    endif
 
    ! 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_int > min_time_step) .or. (time_step_int >= time_step))
    coupled_gl = .false.
 
    call update_ice_shelf(cs%dCS, iss, g, us, time_step_int, time, must_update_vel=update_ice_vel)
 
    call enable_averaging(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_time_step
 
!> \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_time_step - 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