mom_ice_shelf Module Reference

Detailed Description

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.

section_ICE_SHELF

This module implements the thermodynamic aspects of ocean/ice-shelf inter-actions using the MOM framework and coding style.

Derived from code by Chris Little, early 2010.

The ice-sheet dynamics subroutines do the following: initialize_shelf_mass - Initializes the ice shelf mass distribution.

• Initializes h_shelf, h_mask, area_shelf_h
• CURRENTLY: initializes mass_shelf as well, but this is unnecessary, as mass_shelf is initialized based on h_shelf and density_ice immediately afterwards. Possibly subroutine should be renamed update_shelf_mass - updates ice shelf mass via netCDF file USER_update_shelf_mass (TODO). solo_step_ice_shelf - called only in ice-only mode. shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. currently mass_shelf is updated immediately after ice_shelf_advect in fully dynamic mode.

NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement, for subroutines in the velocity solve, and for thickness boundary conditions (this last one may be removed). in other words, interfering with its updates will have implications you might not expect.

Overall issues: Many variables need better documentation and units and the subgrid on which they are discretized.

ICE_SHELF equations

The three fundamental equations are: Heat flux

\[ \qquad \rho_w C_{pw} \gamma_T (T_w - T_b) = \rho_i \dot{m} L_f \]

Salt flux

\[ \qquad \rho_w \gamma_s (S_w - S_b) = \rho_i \dot{m} S_b \]

Freezing temperature

\[ \qquad T_b = a S_b + b + c P \]

where ....

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.

Data Types
type	ice_shelf_cs
	Control structure that contains ice shelf parameters and diagnostics handles. More...

Functions/Subroutines
subroutine, public	shelf_calc_flux (sfc_state, fluxes, Time, time_step, CS, forces)
	Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to use just two equations). See ICE_SHELF equations. More...
subroutine	change_thickness_using_melt (ISS, G, US, time_step, fluxes, density_ice, debug)
	Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting. More...
subroutine, public	add_shelf_forces (G, US, CS, forces, do_shelf_area)
	This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on the ice state in ice_shelf_CS. More...
subroutine	add_shelf_pressure (G, US, CS, fluxes)
	This subroutine adds the ice shelf pressure to the fluxes type. More...
subroutine, public	add_shelf_flux (G, US, CS, sfc_state, fluxes)
	Updates surface fluxes that are influenced by sub-ice-shelf melting. More...
subroutine, public	initialize_ice_shelf (param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, solo_ice_sheet_in)
	Initializes shelf model data, parameters and diagnostics. More...
subroutine	initialize_shelf_mass (G, param_file, CS, ISS, new_sim)
	Initializes shelf mass based on three options (file, zero and user) More...
subroutine	update_shelf_mass (G, US, CS, ISS, Time)
	Updates the ice shelf mass using data from a file. More...
subroutine, public	ice_shelf_save_restart (CS, Time, directory, time_stamped, filename_suffix)
	Save the ice shelf restart file. More...
subroutine, public	ice_shelf_end (CS)
	Deallocates all memory associated with this module. More...
subroutine, public	solo_step_ice_shelf (CS, time_interval, nsteps, Time, min_time_step_in)
	This routine is for stepping a stand-alone ice shelf model without an ocean. More...

Variables
integer	id_clock_shelf
	CPU Clock for the ice shelf code.
integer	id_clock_pass
	CPU Clock for group pass calls.

Function/Subroutine Documentation

add_shelf_flux()

subroutine, public mom_ice_shelf::add_shelf_flux	(	type(ocean_grid_type), intent(inout)	G,
		type(unit_scale_type), intent(in)	US,
		type(ice_shelf_cs), pointer	CS,
		type(surface), intent(inout)	sfc_state,
		type(forcing), intent(inout)	fluxes
	)

Updates surface fluxes that are influenced by sub-ice-shelf melting.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[in,out]	sfc_state	Surface ocean state
[in,out]	fluxes	A structure of surface fluxes that may be used/updated.

Definition at line 885 of file MOM_ice_shelf.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    pointer       :: cs   !< This module's control structure.
  type(surface),         intent(inout) :: sfc_state !< Surface ocean state
  type(forcing),         intent(inout) :: fluxes  !< A structure of surface fluxes that may be used/updated.

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

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

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

  iss => cs%ISS


  call add_shelf_pressure(g, us, cs, fluxes)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  endif ! constant_sea_level

add_shelf_forces()

subroutine, public mom_ice_shelf::add_shelf_forces	(	type(ocean_grid_type), intent(inout)	G,
		type(unit_scale_type), intent(in)	US,
		type(ice_shelf_cs), pointer	CS,
		type(mech_forcing), intent(inout)	forces,
		logical, intent(in), optional	do_shelf_area
	)

This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on the ice state in ice_shelf_CS.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[in,out]	forces	A structure with the driving mechanical forces
[in]	do_shelf_area	If true find the shelf-covered areas.

Definition at line 776 of file MOM_ice_shelf.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    pointer       :: cs   !< This module's control structure.
  type(mech_forcing),    intent(inout) :: forces !< A structure with the driving mechanical forces
  logical, optional,     intent(in)    :: do_shelf_area !< If true find the shelf-covered areas.

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

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

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

  iss => cs%ISS

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

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

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

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

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

add_shelf_pressure()

subroutine mom_ice_shelf::add_shelf_pressure ( type(ocean_grid_type), intent(inout)  G, type(unit_scale_type), intent(in)  US, type(ice_shelf_cs), intent(in)  CS, type(forcing), intent(inout)  fluxes  )

private

This subroutine adds the ice shelf pressure to the fluxes type.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
[in]	cs	This module's control structure.
[in,out]	fluxes	A structure of surface fluxes that may be updated.

Definition at line 856 of file MOM_ice_shelf.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    intent(in)    :: cs   !< This module's control structure.
  type(forcing),         intent(inout) :: fluxes  !< A structure of surface fluxes that may be updated.

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

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

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

change_thickness_using_melt()

subroutine mom_ice_shelf::change_thickness_using_melt ( type(ice_shelf_state), intent(inout)  ISS, type(ocean_grid_type), intent(inout)  G, type(unit_scale_type), intent(in)  US, real, intent(in)  time_step, type(forcing), intent(inout)  fluxes, real, intent(in)  density_ice, logical, intent(in), optional  debug  )

private

Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting.

Parameters

[in,out]	g	The ocean's grid structure.
[in,out]	iss	A structure with elements that describe the ice-shelf state
[in]	us	A dimensional unit scaling type
[in]	time_step	The time step for this update [T ~> s].
[in,out]	fluxes	structure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in]	density_ice	The density of ice-shelf ice [R ~> kg m-3].
[in]	debug	If present and true, write chksums

Definition at line 728 of file MOM_ice_shelf.F90.

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

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

  i_rho_ice = 1.0 / density_ice


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

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

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

ice_shelf_end()

subroutine, public mom_ice_shelf::ice_shelf_end

(

type(ice_shelf_cs), pointer

CS

)

Deallocates all memory associated with this module.

Parameters

cs	A pointer to the ice shelf control structure

Definition at line 1790 of file MOM_ice_shelf.F90.

  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)

ice_shelf_save_restart()

subroutine, public mom_ice_shelf::ice_shelf_save_restart	(	type(ice_shelf_cs), pointer	CS,
		type(time_type), intent(in)	Time,
		character(len=*), intent(in), optional	directory,
		logical, intent(in), optional	time_stamped,
		character(len=*), intent(in), optional	filename_suffix
	)

Save the ice shelf restart file.

Parameters

	cs	ice shelf control structure
[in]	time	model time at this call
[in]	directory	An optional directory into which to write these restart files.
[in]	time_stamped	f true, the restart file names include a unique time stamp. The default is false.
[in]	filename_suffix	An optional suffix (e.g., a time-stamp) to append to the restart file names.

Definition at line 1767 of file MOM_ice_shelf.F90.

  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)

initialize_ice_shelf()

subroutine, public mom_ice_shelf::initialize_ice_shelf	(	type(param_file_type), intent(in)	param_file,
		type(ocean_grid_type), pointer	ocn_grid,
		type(time_type), intent(inout)	Time,
		type(ice_shelf_cs), pointer	CS,
		type(diag_ctrl), intent(in), target	diag,
		type(mech_forcing), intent(inout), optional	forces,
		type(forcing), intent(inout), optional	fluxes,
		type(time_type), intent(in), optional	Time_in,
		logical, intent(in), optional	solo_ice_sheet_in
	)

Initializes shelf model data, parameters and diagnostics.

Parameters

[in]	param_file	A structure to parse for run-time parameters
	ocn_grid	The calling ocean model's horizontal grid structure
[in,out]	time	The clock that that will indicate the model time
	cs	A pointer to the ice shelf control structure
[in]	diag	A structure that is used to regulate the diagnostic output.
[in,out]	fluxes	A structure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in,out]	forces	A structure with the driving mechanical forces
[in]	time_in	The time at initialization.
[in]	solo_ice_sheet_in	If present, this indicates whether a solo ice-sheet driver.

Definition at line 1084 of file MOM_ice_shelf.F90.

  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(ocean_grid_type),        pointer       :: ocn_grid   !< The calling ocean model's horizontal grid structure
  type(time_type),              intent(inout) :: time !< The clock that that will indicate the model time
  type(ice_shelf_cs),           pointer       :: cs   !< A pointer to the ice shelf control structure
  type(diag_ctrl),    target,   intent(in)    :: diag !< A structure that is used to regulate the diagnostic output.
  type(forcing),      optional, intent(inout) :: fluxes !< A structure containing pointers to any possible
                                                   !! thermodynamic or mass-flux forcing fields.
  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
  type(time_type),    optional, intent(in)    :: time_in !< The time at initialization.
  logical,            optional, intent(in)    :: solo_ice_sheet_in !< If present, this indicates whether
                                                   !! a solo ice-sheet driver.

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

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

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

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

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

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

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

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

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

  cs%diag => diag

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

  if (present(time_in)) time = time_in

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

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

  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
  call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
                 "If true, write verbose debugging messages for the ice shelf.", &
                 default=debug)
  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
                 "If true, the ice sheet mass can evolve with time.", &
                 default=.false.)
  if (shelf_mass_is_dynamic) then
    call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
                 "If true, user provided code specifies the ice-shelf "//&
                 "movement instead of the dynamic ice model.", default=.false.)
    cs%active_shelf_dynamics = .not.cs%override_shelf_movement
    call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
                 "If true, regularize the floatation condition at the "//&
                 "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
    call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
                 "If true, let the floatation condition be determined by "//&
                 "ocean column thickness. This means that update_OD_ffrac "//&
                 "will be called.  GL_REGULARIZE and GL_COUPLE are exclusive.", &
                 default=.false., do_not_log=cs%GL_regularize)
    if (cs%GL_regularize) cs%GL_couple = .false.
  endif

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  call eos_init(param_file, cs%eqn_of_state)

  !! new parameters that need to be in MOM_input

  if (cs%active_shelf_dynamics) then

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

    call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", cs%input_flux, &
                 "volume flux at upstream boundary", units="m2 s-1", default=0.)
    call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", cs%input_thickness, &
                 "flux thickness at upstream boundary", units="m", default=1000.)
  else
    ! This is here because of inconsistent defaults.  I don't know why.  RWH
    call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
                 "A typical density of ice.", units="kg m-3", default=900.0, scale=us%kg_m3_to_R)
  endif
  call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", &
                cs%min_thickness_simple_calve, &
                 "Min thickness rule for the very simple calving law",&
                 units="m", default=0.0, scale=us%m_to_Z)

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

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

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

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

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

  call destroy_dyn_horgrid(dg)

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

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

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

  cs%restart_output_dir = dirs%restart_output_dir

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

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

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

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

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

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

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

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

  endif

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

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

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

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

    if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) then
      z_rescale = us%m_to_Z / us%m_to_Z_restart
      do j=g%jsc,g%jec ; do i=g%isc,g%iec
        iss%h_shelf(i,j) = z_rescale * iss%h_shelf(i,j)
      enddo ; enddo
    endif

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

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

  endif ! .not. new_sim

  cs%Time = time

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

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

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

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

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

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

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

  call fix_restart_unit_scaling(us)

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

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


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

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

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

initialize_shelf_mass()

subroutine mom_ice_shelf::initialize_shelf_mass ( type(ocean_grid_type), intent(in)  G, type(param_file_type), intent(in)  param_file, type(ice_shelf_cs), pointer  CS, type(ice_shelf_state), intent(inout)  ISS, logical, intent(in), optional  new_sim  )

private

Initializes shelf mass based on three options (file, zero and user)

Parameters

[in]	g	The ocean's grid structure.
[in]	param_file	A structure to parse for run-time parameters
	cs	A pointer to the ice shelf control structure
[in,out]	iss	The ice shelf state type that is being updated
[in]	new_sim	If present and false, this run is being restarted

Definition at line 1643 of file MOM_ice_shelf.F90.

  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

shelf_calc_flux()

subroutine, public mom_ice_shelf::shelf_calc_flux	(	type(surface), intent(inout)	sfc_state,
		type(forcing), intent(inout)	fluxes,
		type(time_type), intent(in)	Time,
		real, intent(in)	time_step,
		type(ice_shelf_cs), pointer	CS,
		type(mech_forcing), intent(inout), optional	forces
	)

Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to use just two equations). See ICE_SHELF equations.

Parameters

[in,out]	sfc_state	A structure containing fields that describe the surface state of the ocean. The intent is only inout to allow for halo updates.
[in,out]	fluxes	structure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in]	time	Start time of the fluxes.
[in]	time_step	Length of time over which these fluxes will be applied [s].
	cs	A pointer to the control structure returned by a previous call to initialize_ice_shelf.
[in,out]	forces	A structure with the driving mechanical forces

Definition at line 195 of file MOM_ice_shelf.F90.

  type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                !! describe the surface state of the ocean.  The
                                                !! intent is only inout to allow for halo updates.
  type(forcing),         intent(inout) :: fluxes !< structure containing pointers to any possible
                                                !! thermodynamic or mass-flux forcing fields.
  type(time_type),       intent(in)    :: time  !< Start time of the fluxes.
  real,                  intent(in)    :: time_step !< Length of time over which these fluxes
                                                !! will be applied [s].
  type(ice_shelf_cs),    pointer       :: cs    !< A pointer to the control structure returned
                                                !! by a previous call to initialize_ice_shelf.
  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces

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

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

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

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

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

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

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

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

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

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

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

  !update time
  cs%Time = time

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

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

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

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

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

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

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

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

          ustar_h = fluxes%ustar_shelf(i,j)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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


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

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

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

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

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

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

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

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

    enddo ! i-loop
  enddo ! j-loop


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

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

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

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

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

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

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

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

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

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

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

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

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

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

  endif

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

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

  call cpu_clock_end(id_clock_shelf)

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

solo_step_ice_shelf()

subroutine, public mom_ice_shelf::solo_step_ice_shelf	(	type(ice_shelf_cs), pointer	CS,
		type(time_type), intent(in)	time_interval,
		integer, intent(inout)	nsteps,
		type(time_type), intent(inout)	Time,
		real, intent(in), optional	min_time_step_in
	)

This routine is for stepping a stand-alone ice shelf model without an ocean.

Parameters

	cs	A pointer to the ice shelf control structure
[in]	time_interval	The time interval for this update [s].
[in,out]	nsteps	The running number of ice shelf steps.
[in,out]	time	The current model time
[in]	min_time_step_in	The minimum permitted time step [T ~> s].

Definition at line 1804 of file MOM_ice_shelf.F90.

  type(ice_shelf_cs), pointer    :: cs      !< A pointer to the ice shelf control structure
  type(time_type), intent(in)    :: time_interval !< The time interval for this update [s].
  integer,         intent(inout) :: nsteps  !< The running number of ice shelf steps.
  type(time_type), intent(inout) :: time    !< The current model time
  real,  optional, intent(in)    :: min_time_step_in !< The minimum permitted time step [T ~> s].

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

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

  remaining_time = us%s_to_T*time_type_to_real(time_interval)

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

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

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

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

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

    remaining_time = remaining_time - time_step

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

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

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

  enddo

update_shelf_mass()

subroutine mom_ice_shelf::update_shelf_mass ( type(ocean_grid_type), intent(inout)  G, type(unit_scale_type), intent(in)  US, type(ice_shelf_cs), intent(in)  CS, type(ice_shelf_state), intent(inout)  ISS, type(time_type), intent(in)  Time  )

private

Updates the ice shelf mass using data from a file.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
[in]	cs	A pointer to the ice shelf control structure
[in,out]	iss	The ice shelf state type that is being updated
[in]	time	The current model time

Definition at line 1724 of file MOM_ice_shelf.F90.

  type(ocean_grid_type), intent(inout) :: g   !< The ocean's grid structure.
  type(unit_scale_type), intent(in)    :: us   !< A dimensional unit scaling type
  type(ice_shelf_cs),    intent(in)    :: cs  !< A pointer to the ice shelf control structure
  type(ice_shelf_state), intent(inout) :: iss !< The ice shelf state type that is being updated
  type(time_type),       intent(in)    :: time !< The current model time

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

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

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

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

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

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