MOM_state_initialization.F90

Go to the documentation of this file.

!> Initialize state variables, u, v, h, T and S.
module mom_state_initialization

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

use mom_debugging, only : hchksum, qchksum, uvchksum
use mom_coms, only : max_across_pes, min_across_pes, reproducing_sum
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only :  clock_routine, clock_loop
use mom_domains, only : pass_var, pass_vector, sum_across_pes, broadcast
use mom_domains, only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
use mom_eos, only : find_depth_of_pressure_in_cell
use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser, only : get_param, read_param, log_param, param_file_type
use mom_file_parser, only : log_version
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type, ispointincell
use mom_interface_heights, only : find_eta
use mom_io, only : close_file, fieldtype, file_exists
use mom_io, only : open_file, read_data, read_axis_data, single_file, multiple
use mom_io, only : slasher, vardesc, write_field
use mom_io, only : east_face, north_face
use mom_open_boundary, only : ocean_obc_type, open_boundary_init
use mom_open_boundary, only : obc_none, obc_simple
use mom_open_boundary, only : open_boundary_query, set_tracer_data
!use MOM_open_boundary, only : set_3D_OBC_data
use mom_grid_initialize, only : initialize_masks, set_grid_metrics
use mom_restart, only : restore_state, mom_restart_cs
use mom_sponge, only : set_up_sponge_field, set_up_sponge_ml_density
use mom_sponge, only : initialize_sponge, sponge_cs
use mom_ale_sponge, only : set_up_ale_sponge_field, initialize_ale_sponge
use mom_ale_sponge, only : ale_sponge_cs
use mom_string_functions, only : uppercase, lowercase
use mom_time_manager, only : time_type, set_time
use mom_tracer_registry, only : add_tracer_obc_values, tracer_registry_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : setverticalgridaxes, verticalgrid_type
use mom_ale, only : pressure_gradient_plm
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
use mom_eos, only : int_specific_vol_dp, convert_temp_salt_for_teos10
use user_initialization, only : user_initialize_thickness, user_initialize_velocity
use user_initialization, only : user_init_temperature_salinity
use user_initialization, only : user_set_obc_data
use user_initialization, only : user_initialize_sponges
use dome_initialization, only : dome_initialize_thickness
use dome_initialization, only : dome_set_obc_data
use dome_initialization, only : dome_initialize_sponges
use isomip_initialization, only : isomip_initialize_thickness
use isomip_initialization, only : isomip_initialize_sponges
use isomip_initialization, only : isomip_initialize_temperature_salinity
use baroclinic_zone_initialization, only : baroclinic_zone_init_temperature_salinity
use benchmark_initialization, only : benchmark_initialize_thickness
use benchmark_initialization, only : benchmark_init_temperature_salinity
use circle_obcs_initialization, only : circle_obcs_initialize_thickness
use lock_exchange_initialization, only : lock_exchange_initialize_thickness
use external_gwave_initialization, only : external_gwave_initialize_thickness
use dome2d_initialization, only : dome2d_initialize_thickness
use dome2d_initialization, only : dome2d_initialize_temperature_salinity
use dome2d_initialization, only : dome2d_initialize_sponges
use adjustment_initialization, only : adjustment_initialize_thickness
use adjustment_initialization, only : adjustment_initialize_temperature_salinity
use sloshing_initialization, only : sloshing_initialize_thickness
use sloshing_initialization, only : sloshing_initialize_temperature_salinity
use seamount_initialization, only : seamount_initialize_thickness
use seamount_initialization, only : seamount_initialize_temperature_salinity
use phillips_initialization, only : phillips_initialize_thickness
use phillips_initialization, only : phillips_initialize_velocity
use phillips_initialization, only : phillips_initialize_sponges
use rossby_front_2d_initialization, only : rossby_front_initialize_thickness
use rossby_front_2d_initialization, only : rossby_front_initialize_temperature_salinity
use rossby_front_2d_initialization, only : rossby_front_initialize_velocity
use scm_idealized_hurricane, only : scm_idealized_hurricane_ts_init
use scm_cvmix_tests, only: scm_cvmix_tests_ts_init
use supercritical_initialization, only : supercritical_set_obc_data
use soliton_initialization, only : soliton_initialize_velocity
use soliton_initialization, only : soliton_initialize_thickness
use bfb_initialization, only : bfb_initialize_sponges_southonly
use dense_water_initialization, only : dense_water_initialize_ts
use dense_water_initialization, only : dense_water_initialize_sponges

use midas_vertmap, only : find_interfaces, tracer_z_init
use midas_vertmap, only : determine_temperature

use mom_ale, only : ale_initregridding, ale_cs, ale_initthicknesstocoord
use mom_ale, only : ale_remap_scalar, ale_build_grid, ale_regrid_accelerated
use mom_regridding, only : regridding_cs, set_regrid_params, getcoordinateresolution
use mom_remapping, only : remapping_cs, initialize_remapping
use mom_remapping, only : remapping_core_h
use mom_tracer_initialization_from_z, only : horiz_interp_and_extrap_tracer

implicit none ; private

#include <MOM_memory.h>

public mom_initialize_state

character(len=40)  :: mdl = "MOM_state_initialization" ! This module's name.

contains

! -----------------------------------------------------------------------------
subroutine mom_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, &
                                restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
                                ALE_sponge_CSp, OBC, Time_in)
  type(ocean_grid_type),      intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),    intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(out)   :: u    !< The zonal velocity that is being initialized,
                                                    !! in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(out)   :: v    !< The meridional velocity that is being
                                                    !! initialized, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                              intent(out)   :: h    !< Layer thicknesses, in H (usually m or
                                                    !! kg m-2)
  type(thermo_var_ptrs),      intent(inout) :: tv   !< A structure pointing to various thermodynamic
                                                    !! variables
  type(time_type),            intent(inout) :: Time !< Time at the start of the run segment.
  type(param_file_type),      intent(in)    :: PF   !< A structure indicating the open file to parse
                                                    !! for model parameter values.
  type(directories),          intent(in)    :: dirs !< A structure containing several relevant
                                                    !! directory paths.
  type(mom_restart_cs),       pointer       :: restart_CS !< A pointer to the restart control
                                                    !! structure.
  type(ale_cs),               pointer       :: ALE_CSp
  type(tracer_registry_type), pointer       :: tracer_Reg
  type(sponge_cs),            pointer       :: sponge_CSp
  type(ale_sponge_cs),        pointer       :: ALE_sponge_CSp
  type(ocean_obc_type),       pointer       :: OBC
  type(time_type), optional,  intent(in)    :: Time_in !< Time at the start of the run segment.
                                                     !! Time_in overrides any value set for Time.
! Arguments: u  - Zonal velocity, in m s-1.
!  (out)     v  - Meridional velocity, in m s-1.
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (out)     Time    - Time at the start of the run segment.
!  (inout)   G       - The ocean's grid structure.
!  (in)      GV      - The ocean's vertical grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.
!  (in)      dirs    - A structure containing several relevant directory paths.
!  (inout)   restart_CS - A pointer to the restart control structure.
!  (inout)   CS      - A structure of pointers to be exchanged with MOM.F90.
!  (in)      Time_in - Time at the start of the run segment. Time_in overrides
!                      any value set for Time.

  character(len=200) :: filename   ! The name of an input file.
  character(len=200) :: filename2  ! The name of an input files.
  character(len=200) :: inputdir   ! The directory where NetCDF input files are.
  character(len=200) :: config
  logical :: from_Z_file, useALE
  logical :: new_sim
  integer :: write_geom
  logical :: use_temperature, use_sponge
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  logical :: depress_sfc ! If true, remove the mass that would be displaced
                         ! by a large surface pressure by squeezing the column.
  logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced
                         ! by a large surface pressure, such as with an ice sheet.
  logical :: regrid_accelerate
  integer :: regrid_iterations
  logical :: Analytic_FV_PGF, obsol_test
  logical :: convert
  logical :: just_read  ! If true, only read the parameters because this
                        ! is a run from a restart file; this option
                        ! allows the use of Fatal unused parameters.
  type(eos_type), pointer :: eos => null()
  logical :: debug    ! indicates whether to write debugging output
! This include declares and sets the variable "version".
#include "version_variable.h"
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  call calltree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
  call log_version(pf, mdl, version, "")
  call get_param(pf, mdl, "DEBUG", debug, default=.false.)

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

  just_read = .not.new_sim

  call get_param(pf, mdl, "INPUTDIR", inputdir, &
         "The directory in which input files are found.", default=".")
  inputdir = slasher(inputdir)

  use_temperature = ASSOCIATED(tv%T)
  useale = associated(ale_csp)
  use_eos = associated(tv%eqn_of_state)
  if (use_eos) eos => tv%eqn_of_state

!====================================================================
!    Initialize temporally evolving fields, either as initial
!  conditions or by reading them from a restart (or saves) file.
!====================================================================

  if (new_sim) then
    call mom_mesg("Run initialized internally.", 3)

    if (present(time_in)) time = time_in
    ! Otherwise leave Time at its input value.

    ! This initialization should not be needed. Certainly restricting it
    ! to the computational domain helps detect possible uninitialized
    ! data in halos which should be covered by the pass_var(h) later.
    !do k = 1, nz; do j = js, je; do i = is, ie
    !  h(i,j,k) = 0.
    !enddo
  endif

  ! The remaining initialization calls are done, regardless of whether the
  ! fields are actually initialized here (if just_read=.false.) or whether it
  ! is just to make sure that all valid parameters are read to enable the
  ! detection of unused parameters.
  call get_param(pf, mdl, "INIT_LAYERS_FROM_Z_FILE", from_z_file, &
             "If true, intialize the layer thicknesses, temperatures, \n"//&
             "and salnities from a Z-space file on a latitude- \n"//&
             "longitude grid.", default=.false., do_not_log=just_read)

  if (from_z_file) then
!     Initialize thickness and T/S from z-coordinate data in a file.
    if (.NOT.use_temperature) call mom_error(fatal,"MOM_initialize_state : "//&
       "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")

    call mom_temp_salt_initialize_from_z(h, tv, g, gv, pf, just_read_params=just_read)

  else
!     Initialize thickness, h.
    call get_param(pf, mdl, "THICKNESS_CONFIG", config, &
             "A string that determines how the initial layer \n"//&
             "thicknesses are specified for a new run: \n"//&
             " \t file - read interface heights from the file specified \n"//&
             " \t thickness_file - read thicknesses from the file specified \n"//&
             " \t\t by (THICKNESS_FILE).\n"//&
             " \t coord - determined by ALE coordinate.\n"//&
             " \t uniform - uniform thickness layers evenly distributed \n"//&
             " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
             " \t DOME - use a slope and channel configuration for the \n"//&
             " \t\t DOME sill-overflow test case. \n"//&
             " \t ISOMIP - use a configuration for the \n"//&
             " \t\t ISOMIP test case. \n"//&
             " \t benchmark - use the benchmark test case thicknesses. \n"//&
             " \t search - search a density profile for the interface \n"//&
             " \t\t densities. This is not yet implemented. \n"//&
             " \t circle_obcs - the circle_obcs test case is used. \n"//&
             " \t DOME2D - 2D version of DOME initialization. \n"//&
             " \t adjustment2d - TBD AJA. \n"//&
             " \t sloshing - TBD AJA. \n"//&
             " \t seamount - TBD AJA. \n"//&
             " \t soliton - Equatorial Rossby soliton. \n"//&
             " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
             " \t USER - call a user modified routine.", &
             fail_if_missing=new_sim, do_not_log=just_read)
    select case (trim(config))
       case ("file"); call initialize_thickness_from_file(h, g, gv, pf, .false., just_read_params=just_read)
       case ("thickness_file"); call initialize_thickness_from_file(h, g, gv, pf, .true., just_read_params=just_read)
       case ("coord")
         if (new_sim .and. useale) then
           call ale_initthicknesstocoord( ale_csp, g, gv, h )
         elseif (new_sim) then
           call mom_error(fatal, "MOM_initialize_state: USE_REGRIDDING must be True "//&
                                 "for THICKNESS_CONFIG of 'coord'")
         endif
       case ("uniform"); call initialize_thickness_uniform(h, g, gv, pf, &
                                  just_read_params=just_read)
       case ("DOME"); call dome_initialize_thickness(h, g, gv, pf, &
                               just_read_params=just_read)
       case ("ISOMIP"); call isomip_initialize_thickness(h, g, gv, pf, tv, &
                                 just_read_params=just_read)
       case ("benchmark"); call benchmark_initialize_thickness(h, g, gv, pf, &
                                    tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
       case ("search"); call initialize_thickness_search
       case ("circle_obcs"); call circle_obcs_initialize_thickness(h, g, gv, pf, &
                                      just_read_params=just_read)
       case ("lock_exchange"); call lock_exchange_initialize_thickness(h, g, gv, &
                                        pf, just_read_params=just_read)
       case ("external_gwave"); call external_gwave_initialize_thickness(h, g, &
                                         pf, just_read_params=just_read)
       case ("DOME2D"); call dome2d_initialize_thickness(h, g, gv, pf, &
                                 just_read_params=just_read)
       case ("adjustment2d"); call adjustment_initialize_thickness(h, g, gv, &
                                       pf, just_read_params=just_read)
       case ("sloshing"); call sloshing_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("seamount"); call seamount_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("soliton"); call soliton_initialize_thickness(h, g)
       case ("phillips"); call phillips_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("rossby_front"); call rossby_front_initialize_thickness(h, g, gv, &
                                       pf, just_read_params=just_read)
       case ("USER"); call user_initialize_thickness(h, g, pf, tv%T, &
                               just_read_params=just_read)
       case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
            "Unrecognized layer thickness configuration "//trim(config))
    end select

!     Initialize temperature and salinity (T and S).
    if ( use_temperature ) then
      call get_param(pf, mdl, "TS_CONFIG", config, &
             "A string that determines how the initial tempertures \n"//&
             "and salinities are specified for a new run: \n"//&
             " \t file - read velocities from the file specified \n"//&
             " \t\t by (TS_FILE). \n"//&
             " \t fit - find the temperatures that are consistent with \n"//&
             " \t\t the layer densities and salinity S_REF. \n"//&
             " \t TS_profile - use temperature and salinity profiles \n"//&
             " \t\t (read from TS_FILE) to set layer densities. \n"//&
             " \t benchmark - use the benchmark test case T & S. \n"//&
             " \t linear - linear in logical layer space. \n"//&
             " \t DOME2D - 2D DOME initialization. \n"//&
             " \t ISOMIP - ISOMIP initialization. \n"//&
             " \t adjustment2d - TBD AJA. \n"//&
             " \t sloshing - TBD AJA. \n"//&
             " \t seamount - TBD AJA. \n"//&
             " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
             " \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//&
             " \t SCM_CVmix_tests - used in the SCM CVmix tests.\n"//&
             " \t USER - call a user modified routine.", &
             fail_if_missing=new_sim, do_not_log=just_read)
!              " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
      select case (trim(config))
        case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, g, gv, pf, &
                               eos, tv%P_Ref, just_read_params=just_read)
        case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, g, &
                                pf, just_read_params=just_read)
        case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
                                     g, gv, pf, eos, tv%P_Ref, just_read_params=just_read)
        case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, &
                                       g, pf, just_read_params=just_read)
        case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, g, pf, &
                                  just_read_params=just_read)
        case ("DOME2D"); call dome2d_initialize_temperature_salinity ( tv%T, &
                                  tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("ISOMIP"); call isomip_initialize_temperature_salinity ( tv%T, &
                                  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
        case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
                                        tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
                                           tv%S, h, g, pf, just_read_params=just_read)
        case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
                                    tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
                                    tv%S, h, g, gv, pf, eos, just_read_params=just_read)
        case ("rossby_front"); call rossby_front_initialize_temperature_salinity ( tv%T, &
                                        tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("SCM_ideal_hurr"); call scm_idealized_hurricane_ts_init ( tv%T, &
                                          tv%S, h, g, gv, pf, just_read_params=just_read)
        case ("SCM_CVmix_tests"); call scm_cvmix_tests_ts_init (tv%T, &
                                           tv%S, h, g, gv, pf, just_read_params=just_read)
        case ("dense"); call dense_water_initialize_ts(g, gv, pf, eos, tv%T, tv%S, &
                                 h, just_read_params=just_read)
        case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, g, pf, eos, &
                                just_read_params=just_read)
        case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
               "Unrecognized Temp & salt configuration "//trim(config))
      end select
    endif
  endif  ! not from_Z_file.

  ! The thicknesses in halo points might be needed to initialize the velocities.
  if (new_sim) call pass_var(h, g%Domain)

!   Initialize velocity components, u and v
  call get_param(pf, mdl, "VELOCITY_CONFIG", config, &
       "A string that determines how the initial velocities \n"//&
       "are specified for a new run: \n"//&
       " \t file - read velocities from the file specified \n"//&
       " \t\t by (VELOCITY_FILE). \n"//&
       " \t zero - the fluid is initially at rest. \n"//&
       " \t uniform - the flow is uniform (determined by\n"//&
       " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
       " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
       " \t soliton - Equatorial Rossby soliton.\n"//&
       " \t USER - call a user modified routine.", default="zero", &
       do_not_log=just_read)
  select case (trim(config))
     case ("file"); call initialize_velocity_from_file(u, v, g, pf, &
                             just_read_params=just_read)
     case ("zero"); call initialize_velocity_zero(u, v, g, pf, &
                             just_read_params=just_read)
     case ("uniform"); call initialize_velocity_uniform(u, v, g, pf, &
                                just_read_params=just_read)
     case ("circular"); call initialize_velocity_circular(u, v, g, pf, &
                                 just_read_params=just_read)
     case ("phillips"); call phillips_initialize_velocity(u, v, g, gv, pf, &
                                 just_read_params=just_read)
     case ("rossby_front"); call rossby_front_initialize_velocity(u, v, h, &
                                     g, gv, pf, just_read_params=just_read)
     case ("soliton"); call soliton_initialize_velocity(u, v, h, g)
     case ("USER"); call user_initialize_velocity(u, v, g, pf, &
                             just_read_params=just_read)
     case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
          "Unrecognized velocity configuration "//trim(config))
  end select

  if (new_sim) call pass_vector(u, v, g%Domain)
  if (debug .and. new_sim) then
    call uvchksum("MOM_initialize_state [uv]", u, v, g%HI, haloshift=1)
  endif

!   Optionally convert the thicknesses from m to kg m-2.  This is particularly
! useful in a non-Boussinesq model.
  call get_param(pf, mdl, "CONVERT_THICKNESS_UNITS", convert, &
               "If true,  convert the thickness initial conditions from \n"//&
               "units of m to kg m-2 or vice versa, depending on whether \n"//&
               "BOUSSINESQ is defined. This does not apply if a restart \n"//&
               "file is read.", default=.false., do_not_log=just_read)
  if (new_sim) then
    if (convert .and. .not. gv%Boussinesq) then
      ! Convert h from m to kg m-2 then to thickness units (H)
      call convert_thickness(h, g, gv, tv)
    elseif (gv%Boussinesq) then
      ! Convert h from m to thickness units (H)
      do k = 1, nz; do j = js, je; do i = is, ie
        h(i,j,k) = h(i,j,k)*gv%m_to_H
      enddo ; enddo ; enddo
    else
      do k = 1, nz; do j = js, je; do i = is, ie
        h(i,j,k) = h(i,j,k)*gv%kg_m2_to_H
      enddo ; enddo ; enddo
    endif
  endif

!  Remove the mass that would be displaced by an ice shelf or inverse barometer.
  call get_param(pf, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
               "If true,  depress the initial surface to avoid huge \n"//&
               "tsunamis when a large surface pressure is applied.", &
               default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
               "If true, cuts way the top of the column for initial conditions\n"//&
               "at the depth where the hydrostatic presure matches the imposed\n"//&
               "surface pressure which is read from file.", default=.false., &
               do_not_log=just_read)
  if (depress_sfc .and. trim_ic_for_p_surf) call mom_error(fatal, "MOM_initialize_state: "//&
           "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
  if (depress_sfc) call depress_surface(h, g, gv, pf, tv, just_read_params=just_read)
  if (trim_ic_for_p_surf) call trim_for_ice(pf, g, gv, ale_csp, tv, h, just_read_params=just_read)

  ! Perhaps we want to run the regridding coordinate generator for multiple
  ! iterations here so the initial grid is consistent with the coordinate
  if (useale) then
    call get_param(pf, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, &
         "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding\n"//&
         "algorithm to push the initial grid to be consistent with the initial\n"//&
         "condition. Useful only for state-based and iterative coordinates.", &
         default=.false., do_not_log=just_read)
    if (regrid_accelerate) then
      call get_param(pf, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
           "The number of regridding iterations to perform to generate\n"//&
           "an initial grid that is consistent with the initial conditions.", &
           default=1, do_not_log=just_read)

      if (new_sim) &
        call ale_regrid_accelerated(ale_csp, g, gv, h, tv, regrid_iterations, h, u, v)
    endif
  endif
  ! This is the end of the block of code that might have initialized fields
  ! internally at the start of a new run.

  if (.not.new_sim) then ! This block restores the state from a restart file.
    !    This line calls a subroutine that reads the initial conditions  !
    !  from a previously generated file.                                 !
    call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
                       g, restart_cs)
    if (present(time_in)) time = time_in
  endif

  if ( use_temperature ) then
    call pass_var(tv%T, g%Domain, complete=.false.)
    call pass_var(tv%S, g%Domain, complete=.false.)
  endif
  call pass_var(h, g%Domain)

  if (debug) then
    call hchksum(h, "MOM_initialize_state: h ", g%HI, haloshift=1, scale=gv%H_to_m)
    if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", g%HI, haloshift=1)
    if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", g%HI, haloshift=1)
  endif

  call get_param(pf, mdl, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. \n"//&
                 "The exact location and properties of those sponges are \n"//&
                 "specified via SPONGE_CONFIG.", default=.false.)
  if ( use_sponge ) then
    call get_param(pf, mdl, "SPONGE_CONFIG", config, &
                 "A string that sets how the sponges are configured: \n"//&
                 " \t file - read sponge properties from the file \n"//&
                 " \t\t specified by (SPONGE_FILE).\n"//&
                 " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
                 " \t DOME - use a slope and channel configuration for the \n"//&
                 " \t\t DOME sill-overflow test case. \n"//&
                 " \t BFB - Sponge at the southern boundary of the domain\n"//&
                 " \t\t for buoyancy-forced basin case.\n"//&
                 " \t USER - call a user modified routine.", default="file")
    select case (trim(config))
      case ("DOME"); call dome_initialize_sponges(g, gv, tv, pf, sponge_csp)
      case ("DOME2D"); call dome2d_initialize_sponges(g, gv, tv, pf, useale, &
                                                      sponge_csp, ale_sponge_csp)
      case ("ISOMIP"); call isomip_initialize_sponges(g, gv, tv, pf, useale, &
                                                     sponge_csp, ale_sponge_csp)
      case ("USER"); call user_initialize_sponges(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("BFB"); call bfb_initialize_sponges_southonly(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("phillips"); call phillips_initialize_sponges(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("dense"); call dense_water_initialize_sponges(g, gv, tv, pf, useale, &
           sponge_csp, ale_sponge_csp)
      case ("file"); call initialize_sponges_file(g, gv, use_temperature, tv, &
                                               pf, sponge_csp)
      case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
             "Unrecognized sponge configuration "//trim(config))
    end select
  endif

  ! Reads OBC parameters not pertaining to the location of the boundaries
  call open_boundary_init(g, pf, obc)

  ! This controls user code for setting open boundary data
  if (associated(obc)) then
    call get_param(pf, mdl, "OBC_USER_CONFIG", config, &
                 "A string that sets how the user code is invoked to set open\n"//&
                 " boundary data: \n"//&
                 "   DOME - specified inflow on northern boundary\n"//&
                 "   tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
                 "   supercritical - now only needed here for the allocations\n"//&
                 "   Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
                 "   shelfwave - Flather with shelf wave forcing on western boundary\n"//&
                 "   USER - user specified", default="none")
    if (trim(config) == "DOME") then
      call dome_set_obc_data(obc, tv, g, gv, pf, tracer_reg)
    elseif (lowercase(trim(config)) == "supercritical") then
      call supercritical_set_obc_data(obc, g, pf)
    elseif (trim(config) == "tidal_bay") then
      obc%update_OBC = .true.
    elseif (trim(config) == "Kelvin") then
      obc%update_OBC = .true.
    elseif (trim(config) == "shelfwave") then
      obc%update_OBC = .true.
    elseif (trim(config) == "USER") then
      call user_set_obc_data(obc, tv, g, pf, tracer_reg)
    elseif (.not. trim(config) == "none") then
      call mom_error(fatal, "The open boundary conditions specified by "//&
              "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.")
    endif
    if (open_boundary_query(obc, apply_open_obc=.true.)) then
      call set_tracer_data(obc, tv, h, g, pf, tracer_reg)
    endif
  endif
! if (open_boundary_query(OBC, apply_nudged_OBC=.true.)) then
!   call set_3D_OBC_data(OBC, tv, h, G, PF, tracer_Reg)
! endif
  ! Still need a way to specify the boundary values
  if (debug.and.associated(obc)) then
    call hchksum(g%mask2dT, 'MOM_initialize_state: mask2dT ', g%HI)
    call uvchksum('MOM_initialize_state: mask2dC[uv]', g%mask2dCu,  &
                  g%mask2dCv, g%HI)
    call qchksum(g%mask2dBu, 'MOM_initialize_state: mask2dBu ', g%HI)
  endif

  call calltree_leave('MOM_initialize_state()')

end subroutine mom_initialize_state
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!>  This subroutine reads the layer thicknesses or interface heights from a file.
subroutine initialize_thickness_from_file(h, G, GV, param_file, file_has_thickness, just_read_params)
  type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized, in m.
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,                 intent(in)  :: file_has_thickness !< If true, this file contains layer
                                                      !! thicknesses; otherwise it contains
                                                      !! interface heights.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      file_has_thickness - If true, this file contains thicknesses;
!                                 otherwise it contains interface heights.

!  This subroutine reads the layer thicknesses from file.
  real :: eta(szi_(g),szj_(g),szk_(g)+1)
  integer :: inconsistent = 0
  real :: dilate     ! The amount by which each layer is dilated to agree
                     ! with the bottom depth and free surface height, nondim.
  logical :: correct_thickness
  logical :: just_read    ! If true, just read parameters but set nothing.  character(len=20) :: verticalCoordinate
  character(len=40)  :: mdl = "initialize_thickness_from_file" ! This subroutine's name.
  character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) &
    call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".", do_not_log=just_read)
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "THICKNESS_FILE", thickness_file, &
                 "The name of the thickness file.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)

  filename = trim(inputdir)//trim(thickness_file)
  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/THICKNESS_FILE", filename)

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

  if (file_has_thickness) then
    if (just_read) return ! All run-time parameters have been read, so return.
    call read_data(filename,"h",h(:,:,:),domain=g%Domain%mpp_domain)
  else
    call get_param(param_file, mdl, "ADJUST_THICKNESS", correct_thickness, &
                 "If true, all mass below the bottom removed if the \n"//&
                 "topography is shallower than the thickness input file \n"//&
                 "would indicate.", default=.false., do_not_log=just_read)
    if (just_read) return ! All run-time parameters have been read, so return.

    call read_data(filename,"eta",eta(:,:,:),domain=g%Domain%mpp_domain)

    if (correct_thickness) then
      call adjustetatofitbathymetry(g, gv, eta, h)
    else
      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
        if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) then
          eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
          h(i,j,k) = gv%Angstrom_z
        else
          h(i,j,k) = eta(i,j,k) - eta(i,j,k+1)
        endif
      enddo ; enddo ; enddo

      do j=js,je ; do i=is,ie
        if (abs(eta(i,j,nz+1) + g%bathyT(i,j)) > 1.0) &
          inconsistent = inconsistent + 1
      enddo ; enddo
      call sum_across_pes(inconsistent)

      if ((inconsistent > 0) .and. (is_root_pe())) then
        write(mesg,'("Thickness initial conditions are inconsistent ",'// &
                 '"with topography in ",I8," places.")') inconsistent
        call mom_error(warning, mesg)
      endif
    endif

  endif
  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_thickness_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
!! If the bottom most interface is below the topography then the bottom-most
!! layers are contracted to GV%Angstrom_z.
!! If the bottom most interface is above the topography then the entire column
!! is dilated (expanded) to fill the void.
!!   @remark{There is a (hard-wired) "tolerance" parameter such that the
!! criteria for adjustment must equal or exceed 10cm.}
!!   @param[in]     G   Grid type
!!   @param[in,out] eta Interface heights
!!   @param[out]    h   Layer thicknesses
subroutine adjustetatofitbathymetry(G, GV, eta, h)
  type(ocean_grid_type),                          intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                        intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)+1),    intent(inout) :: eta
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)),      intent(inout) :: h    !< Layer thicknesses, in m
  ! Local variables
  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
  real, parameter :: hTolerance = 0.1 !<  Tolerance to exceed adjustment criteria (m)
  real :: hTmp, eTmp, dilate
  character(len=100) :: mesg

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  contractions = 0
  do j=js,je ; do i=is,ie
    if (-eta(i,j,nz+1) > g%bathyT(i,j) + htolerance) then
      eta(i,j,nz+1) = -g%bathyT(i,j)
      contractions = contractions + 1
    endif
  enddo ; enddo
  call sum_across_pes(contractions)
  if ((contractions > 0) .and. (is_root_pe())) then
    write(mesg,'("Thickness initial conditions were contracted ",'// &
               '"to fit topography in ",I8," places.")') contractions
    call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
  endif

  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    ! Collapse layers to thinnest possible if the thickness less than
    ! the thinnest possible (or negative).
    if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) then
      eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
      h(i,j,k) = gv%Angstrom_z
    else
      h(i,j,k) = eta(i,j,k) - eta(i,j,k+1)
    endif
  enddo ; enddo ; enddo

  dilations = 0
  do j=js,je ; do i=is,ie
    !   The whole column is dilated to accommodate deeper topography than
    ! the bathymetry would indicate.
    ! This should be...  if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. &
    if (-eta(i,j,nz+1) < g%bathyT(i,j) - htolerance) then
      dilations = dilations + 1
      if (eta(i,j,1) <= eta(i,j,nz+1)) then
        do k=1,nz ; h(i,j,k) = (eta(i,j,1)+g%bathyT(i,j)) / real(nz) ; enddo
      else
        dilate = (eta(i,j,1)+g%bathyT(i,j)) / (eta(i,j,1)-eta(i,j,nz+1))
        do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
      endif
      do k=nz, 2, -1; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k); enddo
    endif
  enddo ; enddo
  call sum_across_pes(dilations)
  if ((dilations > 0) .and. (is_root_pe())) then
    write(mesg,'("Thickness initial conditions were dilated ",'// &
               '"to fit topography in ",I8," places.")') dilations
    call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
  endif

end subroutine adjustetatofitbathymetry
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params)
  type(ocean_grid_type),   intent(in)  :: G           !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)  :: GV          !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized, in m.
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!  This subroutine initializes the layer thicknesses to be uniform.
  character(len=40)  :: mdl = "initialize_thickness_uniform" ! This subroutine's name.
  real :: e0(szk_(g)+1)   ! The resting interface heights, in m, usually !
                          ! negative because it is positive upward.      !
  real :: eta1D(szk_(g)+1)! Interface height relative to the sea surface !
                          ! positive upward, in m.                       !
  logical :: just_read    ! If true, just read parameters but set nothing.  character(len=20) :: verticalCoordinate
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (just_read) return ! This subroutine has no run-time parameters.

  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  if (g%max_depth<=0.) call mom_error(fatal,"initialize_thickness_uniform: "// &
      "MAXIMUM_DEPTH has a non-sensical value! Was it set?")

  do k=1,nz
    e0(k) = -g%max_depth * real(k-1) / real(nz)
  enddo

  do j=js,je ; do i=is,ie
!    This sets the initial thickness (in m) of the layers.  The      !
!  thicknesses are set to insure that: 1.  each layer is at least an !
!  Angstrom thick, and 2.  the interfaces are where they should be   !
!  based on the resting depths and interface height perturbations,   !
!  as long at this doesn't interfere with 1.                         !
    eta1d(nz+1) = -1.0*g%bathyT(i,j)
    do k=nz,1,-1
      eta1d(k) = e0(k)
      if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_z)) then
        eta1d(k) = eta1d(k+1) + gv%Angstrom_z
        h(i,j,k) = gv%Angstrom_z
      else
        h(i,j,k) = eta1d(k) - eta1d(k+1)
      endif
    enddo
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_thickness_uniform
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_thickness_search
! search density space for location of layers
  call mom_error(fatal,"  MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
end subroutine initialize_thickness_search
! -----------------------------------------------------------------------------

subroutine convert_thickness(h, G, GV, tv)
  type(ocean_grid_type),                  intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                                          intent(inout) :: h    !< Layer thicknesses, being
                                                                !! converted from m to H (m or kg
                                                                !! m-2)
  type(thermo_var_ptrs),                  intent(in)    :: tv   !< A structure pointing to various
                                                                !! thermodynamic variables
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    p_top, p_bot
  real :: dz_geo(szi_(g),szj_(g))      ! The change in geopotential height
                                       ! across a layer, in m2 s-2.
  real :: rho(szi_(g))
  real :: I_gEarth
  logical :: Boussinesq
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: itt, max_itt

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  max_itt = 10
  boussinesq = gv%Boussinesq
  i_gearth = 1.0 / gv%g_Earth

  if (boussinesq) then
    call mom_error(fatal,"Not yet converting thickness with Boussinesq approx.")
  else
    if (associated(tv%eqn_of_state)) then
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
      enddo ; enddo
      do k=1,nz
        do j=js,je
          do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
                                 is, ie-is+1, tv%eqn_of_state)
          do i=is,ie
            p_bot(i,j) = p_top(i,j) + gv%g_Earth * h(i,j,k) * rho(i)
          enddo
        enddo

        do itt=1,max_itt
          call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
                                   0.0, g%HI, tv%eqn_of_state, dz_geo)
          if (itt < max_itt) then ; do j=js,je
            call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
                                   is, ie-is+1, tv%eqn_of_state)
            ! Use Newton's method to correct the bottom value.
            !   The hydrostatic equation is linear to such a
            ! high degree that no bounds-checking is needed.
            do i=is,ie
              p_bot(i,j) = p_bot(i,j) + rho(i) * (gv%g_Earth*h(i,j,k) - dz_geo(i,j))
            enddo
          enddo ; endif
        enddo

        do j=js,je ; do i=is,ie
          h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * gv%kg_m2_to_H * i_gearth
        enddo ; enddo
      enddo
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        h(i,j,k) = h(i,j,k) * gv%Rlay(k) * gv%kg_m2_to_H
      enddo ; enddo ; enddo
    endif
  endif

end subroutine convert_thickness

subroutine depress_surface(h, G, GV, param_file, tv, just_read_params)
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    eta_sfc  ! The free surface height that the model should use, in m.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    eta  ! The free surface height that the model should use, in m.
  real :: dilate  ! A ratio by which layers are dilated, nondim.
  real :: scale_factor ! A scaling factor for the eta_sfc values that are read
                       ! in, which can be used to change units, for example.
  character(len=40)  :: mdl = "depress_surface" ! This subroutine's name.
  character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
  character(len=200) :: filename, eta_srf_var  ! Strings for file/path
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  ! Read the surface height (or pressure) from a file.

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
                 "The initial condition file for the surface height.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
                 "The initial condition variable for the surface height.",&
                 default="SSH", do_not_log=just_read)
  filename = trim(inputdir)//trim(eta_srf_file)
  if (.not.just_read) &
    call log_param(param_file,  mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)

  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
                 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into \n"//&
                 "units of m", units="variable", default=1.0, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call read_data(filename,eta_srf_var,eta_sfc,domain=g%Domain%mpp_domain)

  if (scale_factor /= 1.0) then ; do j=js,je ; do i=is,ie
    eta_sfc(i,j) = eta_sfc(i,j) * scale_factor
  enddo ; enddo ; endif

  ! Convert thicknesses to interface heights.
  call find_eta(h, tv, gv%g_Earth, g, gv, eta)

  do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
!    if (eta_sfc(i,j) < eta(i,j,nz+1)) then
      ! Issue a warning?
!    endif
    if (eta_sfc(i,j) > eta(i,j,1)) then
      ! Dilate the water column to agree, but only up to 10-fold.
      if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
        dilate = 10.0
        call mom_error(warning, "Free surface height dilation attempted "//&
               "to exceed 10-fold.", all_print=.true.)
      else
        dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
      endif
      do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
    elseif (eta(i,j,1) > eta_sfc(i,j)) then
      ! Remove any mass that is above the target free surface.
      do k=1,nz
        if (eta(i,j,k) <= eta_sfc(i,j)) exit
        if (eta(i,j,k+1) >= eta_sfc(i,j)) then
          h(i,j,k) = gv%Angstrom
        else
          h(i,j,k) = max(gv%Angstrom, h(i,j,k) * &
              (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
        endif
      enddo
    endif
  endif ; enddo ; enddo

end subroutine depress_surface

!> Adjust the layer thicknesses by cutting away the top of each model column at the depth
!! where the hydrostatic pressure matches an imposed surface pressure read from file.
subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h, just_read_params)
  type(param_file_type),   intent(in)    :: PF !< Parameter file structure
  type(ocean_grid_type),   intent(in)    :: G !< Ocean grid structure
  type(verticalgrid_type), intent(in)    :: GV !< Vertical grid structure
  type(ale_cs),            pointer       :: ALE_CSp !< ALE control structure
  type(thermo_var_ptrs),   intent(inout) :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(inout) :: h !< Layer thickness (H units, m or Pa)
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  ! Local variables
  character(len=200) :: mdl = "trim_for_ice"
  real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface (Pa)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_t, S_b, T_t, T_b ! Top and bottom edge values for reconstructions
                                                                 ! of salinity and temperature within each layer.
  character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
  real :: scale_factor, min_thickness
  integer :: i, j, k
  logical :: just_read    ! If true, just read parameters but set nothing.
  logical :: use_remapping
  type(remapping_cs), pointer :: remap_CS => null()

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(pf, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
                 "The initial condition file for the surface height.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(pf, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
                 "The initial condition variable for the surface height.", &
                 units="kg m-2", default="", do_not_log=just_read)
  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
  filename = trim(slasher(inputdir))//trim(p_surf_file)
  if (.not.just_read) call log_param(pf,  mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)

  call get_param(pf, mdl, "SURFACE_PRESSURE_SCALE", scale_factor, &
                 "A scaling factor to convert SURFACE_PRESSURE_VAR from\n"//&
                 "file SURFACE_PRESSURE_FILE into a surface pressure.", &
                 units="file dependent", default=1., do_not_log=just_read)
  call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
                 units='m', default=1.e-3, do_not_log=just_read)
  call get_param(pf, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
                 'When trimming the column, also remap T and S.', &
                 default=.false., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call read_data(filename, p_surf_var, p_surf, domain=g%Domain%mpp_domain)
  if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:)

  if (use_remapping) then
    allocate(remap_cs)
    call initialize_remapping(remap_cs, 'PLM', boundary_extrapolation=.true.)
  endif

  ! Find edge values of T and S used in reconstructions
  if ( associated(ale_csp) ) then ! This should only be associated if we are in ALE mode
!   if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then
      call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h)
!   elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then
!     call pressure_gradient_ppm(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
!   endif
  else
!    call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
    do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
      t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
      s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    call cut_off_column_top(gv%ke, tv, gv%Rho0, gv%g_Earth, g%bathyT(i,j), min_thickness, &
               tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), &
               p_surf(i,j), h(i,j,:), remap_cs)
  enddo ; enddo

end subroutine trim_for_ice

!> Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic
!! pressure matches p_surf
subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, &
                              T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS)
  integer,               intent(in)    :: nk !< Number of layers
  type(thermo_var_ptrs), intent(in)    :: tv !< Thermodynamics structure
  real,                  intent(in)    :: Rho0 !< Reference density (kg/m3)
  real,                  intent(in)    :: G_earth !< Gravitational acceleration (m/s2)
  real,                  intent(in)    :: depth !< Depth of ocean column (m)
  real,                  intent(in)    :: min_thickness !< Smallest thickness allowed (m)
  real, dimension(nk),   intent(inout) :: T !< Layer mean temperature
  real, dimension(nk),   intent(in)    :: T_t !< Temperature at top of layer
  real, dimension(nk),   intent(in)    :: T_b !< Temperature at bottom of layer
  real, dimension(nk),   intent(inout) :: S !< Layer mean salinity
  real, dimension(nk),   intent(in)    :: S_t !< Salinity at top of layer
  real, dimension(nk),   intent(in)    :: S_b !< Salinity at bottom of layer
  real,                  intent(in)    :: p_surf !< Imposed pressure on ocean at surface (Pa)
  real, dimension(nk),   intent(inout) :: h !< Layer thickness (H units, m or Pa)
  type(remapping_cs),    pointer       :: remap_CS ! Remapping structure for remapping T and S, if associated
  ! Local variables
  real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions
  real, dimension(nk) :: h0, S0, T0, h1, S1, T1
  real :: P_t, P_b, z_out, e_top
  integer :: k

  ! Calculate original interface positions
  e(nk+1) = -depth
  do k=nk,1,-1
    e(k) = e(k+1) + h(k)
    h0(k) = h(nk+1-k) ! Keep a copy to use in remapping
  enddo

  p_t = 0.
  e_top = e(1)
  do k=1,nk
    call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
                                        p_t, p_surf, rho0, g_earth, tv%eqn_of_state, p_b, z_out)
    if (z_out>=e(k)) then
      ! Imposed pressure was less that pressure at top of cell
      exit
    elseif (z_out<=e(k+1)) then
      ! Imposed pressure was greater than pressure at bottom of cell
      e_top = e(k+1)
    else
      ! Imposed pressure was fell between pressures at top and bottom of cell
      e_top = z_out
      exit
    endif
    p_t = p_b
  enddo
  if (e_top<e(1)) then
    ! Clip layers from the top down, if at all
    do k=1,nk
      if (e(k)>e_top) then
        ! Original e(K) is too high
        e(k) = e_top
        e_top = e_top - min_thickness ! Next interface must be at least this deep
      endif
      ! This layer needs trimming
      h(k) = max( min_thickness, e(k) - e(k+1) )
      if (e(k)<e_top) exit ! No need to go further
    enddo
  endif

  ! Now we need to remap but remapping assumes the surface is at the
  ! same place in the two columns so we turn the column upside down.
  if (associated(remap_cs)) then
    do k=1,nk
      s0(k) = s(nk+1-k)
      t0(k) = t(nk+1-k)
      h1(k) = h(nk+1-k)
    enddo
    call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1)
    call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1)
    do k=1,nk
      s(k) = s1(nk+1-k)
      t(k) = t1(nk+1-k)
    enddo
  endif

end subroutine cut_off_column_top

! -----------------------------------------------------------------------------
subroutine initialize_velocity_from_file(u, v, G, param_file, just_read_params)
  type(ocean_grid_type),   intent(in)  :: G           !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine reads the initial velocity components from file
  character(len=40)  :: mdl = "initialize_velocity_from_file" ! This subroutine's name.
  character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path
  logical :: just_read    ! If true, just read parameters but set nothing.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "VELOCITY_FILE", velocity_file, &
                 "The name of the velocity initial condition file.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  if (just_read) return ! All run-time parameters have been read, so return.

  filename = trim(inputdir)//trim(velocity_file)
  call log_param(param_file, mdl, "INPUTDIR/VELOCITY_FILE", filename)

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

  !  Read the velocities from a netcdf file.
  call read_data(filename,"u",u(:,:,:),domain=g%Domain%mpp_domain,position=east_face)
  call read_data(filename,"v",v(:,:,:),domain=g%Domain%mpp_domain,position=north_face)

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_velocity_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_zero(u, v, G, param_file, just_read_params)
  type(ocean_grid_type),                   intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to zero
  character(len=200) :: mdl = "initialize_velocity_zero" ! This subroutine's name.
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = 0.0
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = 0.0
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_velocity_zero
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_uniform(u, v, G, param_file, just_read_params)
  type(ocean_grid_type),                   intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to uniform
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  real    :: initial_u_const, initial_v_const
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(param_file, mdl, "INITIAL_U_CONST", initial_u_const, &
                 "A initial uniform value for the zonal flow.", &
                 units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INITIAL_V_CONST", initial_v_const, &
                 "A initial uniform value for the meridional flow.", &
                 units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = initial_u_const
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = initial_v_const
  enddo ; enddo ; enddo

end subroutine initialize_velocity_uniform
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_circular(u, v, G, param_file, just_read_params)
  type(ocean_grid_type),                   intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to be circular with
! no flow at edges of domain and center.
  character(len=200) :: mdl = "initialize_velocity_circular"
  real :: circular_max_u
  real :: dpi, psi1, psi2
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(param_file, mdl, "CIRCULAR_MAX_U", circular_max_u, &
                 "The amplitude of zonal flow from which to scale the\n"// &
                 "circular stream function (m/s).", &
                 units="m s-1", default=0., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  dpi=acos(0.0)*2.0 ! pi

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    psi1 = my_psi(i,j)
    psi2 = my_psi(i,j-1)
    u(i,j,k) = (psi1-psi2)/g%dy_Cu(i,j)! *(circular_max_u*G%len_lon/(2.0*dpi))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    psi1 = my_psi(i,j)
    psi2 = my_psi(i-1,j)
    v(i,j,k) = (psi2-psi1)/g%dx_Cv(i,j)! *(circular_max_u*G%len_lon/(2.0*dpi))
  enddo ; enddo ; enddo

  contains

  real function my_psi(ig,jg) ! in-line function
    integer :: ig, jg
    real :: x, y, r
    x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon)/g%len_lon-1.0  ! -1<x<1
    y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat)/g%len_lat-1.0 ! -1<y<1
    r = sqrt( x**2 + y**2 ) ! Circulat stream fn nis fn of radius only
    r = min(1.0,r) ! Flatten stream function in corners of box
    my_psi = 0.5*(1.0 - cos(dpi*r))
    my_psi = my_psi * (circular_max_u*g%len_lon*1e3/dpi) ! len_lon is in km
  end function my_psi

end subroutine initialize_velocity_circular
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_from_file(T, S, G, param_file, just_read_params)
  type(ocean_grid_type),                  intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< The potential temperature that is being initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< The salinity that is being initialized.
  type(param_file_type),                  intent(in)  :: param_file !< A structure to parse for run-time parameters
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: filename, salt_filename ! Full paths to input files
  character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path
  character(len=40)  :: mdl = "initialize_temp_salt_from_file"
  character(len=64)  :: temp_var, salt_var ! Temperature and salinity names in files

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "TS_FILE", ts_file, &
                 "The initial condition file for temperature.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  filename = trim(inputdir)//trim(ts_file)
  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
  call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
                 "The initial condition variable for potential temperature.", &
                 default="PTEMP", do_not_log=just_read)
  call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
                 "The initial condition variable for salinity.", &
                 default="SALT", do_not_log=just_read)
  call get_param(param_file, mdl, "SALT_FILE", salt_file, &
                 "The initial condition file for salinity.", &
                 default=trim(ts_file), do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

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

! Read the temperatures and salinities from netcdf files.           !
  call read_data(filename, temp_var, t(:,:,:), domain=g%Domain%mpp_domain)

  salt_filename = trim(inputdir)//trim(salt_file)
  if (.not.file_exists(salt_filename, g%Domain)) call mom_error(fatal, &
     " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))

  call read_data(salt_filename, salt_var, s(:,:,:), domain=g%Domain%mpp_domain)

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_temp_salt_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_from_profile(T, S, G, param_file, just_read_params)
  type(ocean_grid_type),                  intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S
  type(param_file_type),                  intent(in)  :: param_file !< A structure to parse for run-time parameters
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  real, dimension(SZK_(G)) :: T0, S0
  integer :: i, j, k
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
  character(len=40)  :: mdl = "initialize_temp_salt_from_profile"

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "TS_FILE", ts_file, &
                 "The file with the reference profiles for temperature \n"//&
                 "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(ts_file)
  call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
  if (.not.file_exists(filename)) call mom_error(fatal, &
     " initialize_temp_salt_from_profile: Unable to open "//trim(filename))

! Read the temperatures and salinities from a netcdf file.           !
  call read_data(filename,"PTEMP",t0(:),domain=g%Domain%mpp_domain)
  call read_data(filename,"SALT", s0(:),domain=g%Domain%mpp_domain)

  do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
    t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_temp_salt_from_profile
! -----------------------------------------------------------------------------


! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref, just_read_params)
  type(ocean_grid_type),   intent(in)  :: G            !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)  :: GV           !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(out) :: T            !< The potential temperature that is being
                                                       !! initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(out) :: S            !< The salinity that is being initialized.
  type(param_file_type),   intent(in)  :: param_file   !< A structure to parse for run-time
                                                       !! parameters.
  type(eos_type),          pointer     :: eqn_of_state !< Integer that selects the equatio of state.
  real,                    intent(in)  :: P_Ref        !< The coordinate-density reference pressure
                                                       !! in Pa.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                       !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer that selects the equatio of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.
  real :: T0(szk_(g)), S0(szk_(g))
  real :: T_Ref         ! Reference Temperature
  real :: S_Ref         ! Reference Salinity
  real :: pres(szk_(g))      ! An array of the reference pressure in Pa.
  real :: drho_dT(szk_(g))   ! Derivative of density with temperature in kg m-3 K-1.                              !
  real :: drho_dS(szk_(g))   ! Derivative of density with salinity in kg m-3 PSU-1.                             !
  real :: rho_guess(szk_(g)) ! Potential density at T0 & S0 in kg m-3.
  logical :: fit_salin       ! If true, accept the prescribed temperature and fit the salinity.
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "initialize_temp_salt_fit" ! This subroutine's name.
  integer :: i, j, k, itt, nz
  nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "T_REF", t_ref, &
                 "A reference temperature used in initialization.", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_REF", s_ref, &
                 "A reference salinity used in initialization.", units="PSU", &
                 default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
                 "If true, accept the prescribed temperature and fit the \n"//&
                 "salinity; otherwise take salinity and fit temperature.", &
                 default=.false., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz
    pres(k) = p_ref ; s0(k) = s_ref
    t0(k) = t_ref
  enddo

  call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),eqn_of_state)
  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,eqn_of_state)

  if (fit_salin) then
! A first guess of the layers' temperatures.
    do k=nz,1,-1
      s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
    enddo
! Refine the guesses for each layer.
    do itt=1,6
      call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
      call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
      do k=1,nz
        s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
      enddo
    enddo
  else
! A first guess of the layers' temperatures.
    do k=nz,1,-1
      t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
    enddo
    do itt=1,6
      call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
      call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
      do k=1,nz
        t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
      enddo
    enddo
  endif

  do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
    t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_temp_salt_fit
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_linear(T, S, G, param_file, just_read_params)
  type(ocean_grid_type),                     intent(in)  :: G          !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S
  type(param_file_type),                     intent(in)  :: param_file !< A structure to parse for
                                                                       !! run-time parameters
  logical,       optional,                   intent(in)  :: just_read_params !< If present and true,
                                                                       !! this call will only read
                                                                       !! parameters without
                                                                       !! changing h.

  ! This subroutine initializes linear profiles for T and S according to
  ! reference surface layer salinity and temperature and a specified range.
  ! Note that the linear distribution is set up with respect to the layer
  ! number, not the physical position).
  integer :: k;
  real  :: delta_S, delta_T
  real  :: S_top, T_top ! Reference salinity and temerature within surface layer
  real  :: S_range, T_range ! Range of salinities and temperatures over the vertical
  real  :: delta
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "initialize_temp_salt_linear" ! This subroutine's name.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
  call get_param(param_file, mdl, "T_TOP", t_top, &
                 "Initial temperature of the top surface.", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "T_RANGE", t_range, &
                 "Initial temperature difference (top-bottom).", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_TOP", s_top, &
                 "Initial salinity of the top surface.", &
                 units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_RANGE", s_range, &
                 "Initial salinity difference (top-bottom).", &
                 units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

! ! Prescribe salinity
! delta_S = S_range / ( G%ke - 1.0 );
! S(:,:,1) = S_top;
! do k = 2,G%ke
!   S(:,:,k) = S(:,:,k-1) + delta_S;
! end do
  do k = 1,g%ke
    s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(G%ke))
    t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(G%ke))
  end do

! ! Prescribe temperature
! delta_T = T_range / ( G%ke - 1.0 );
! T(:,:,1) = T_top;
! do k = 2,G%ke
!   T(:,:,k) = T(:,:,k-1) + delta_T;
! end do
! delta = 1;
! T(:,:,G%ke/2 - (delta-1):G%ke/2 + delta) = 1.0;

  call calltree_leave(trim(mdl)//'()')
end subroutine initialize_temp_salt_linear
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine sets the inverse restoration time (Idamp), and
!! the values towards which the interface heights and an arbitrary
!! number of tracers should be restored within each sponge. The
!!interface height is always subject to damping, and must always be
!! the first registered field.
subroutine initialize_sponges_file(G, GV, use_temperature, tv, param_file, CSp)
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure.
  logical,                 intent(in) :: use_temperature
  type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic
                                              !! variables.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
  type(sponge_cs),         pointer    :: CSp  !! A pointer that is set to point to the control
                                              !! structure for this module.
!   This subroutine sets the inverse restoration time (Idamp), and   !
! the values towards which the interface heights and an arbitrary    !
! number of tracers should be restored within each sponge. The       !
! interface height is always subject to damping, and must always be  !
! the first registered field.                                        !

! Arguments: from_file - .true. if the variables that are used here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read for all fields
!                       except the inverse damping rate.
!  (in)      damp_file - The name of the file from which to read the
!                        inverse damping rate.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      use_temperature - If true, T & S are state variables.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CSp - A pointer that is set to point to the control structure
!                  for this module

  real :: eta(szi_(g),szj_(g),szk_(g)+1) ! The target interface heights, in m.
  real, dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
    tmp, tmp2 ! A temporary array for tracers.
  real, dimension (SZI_(G),SZJ_(G)) :: &
    tmp_2d ! A temporary array for tracers.

  real :: Idamp(szi_(g),szj_(g))    ! The inverse damping rate, in s-1.
  real :: pres(szi_(g))     ! An array of the reference pressure, in Pa.

  integer :: i, j, k, is, ie, js, je, nz
  character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
  character(len=40) :: mdl = "initialize_sponges_file"
  character(len=200) :: damping_file, state_file  ! Strings for filenames
  character(len=200) :: filename, inputdir ! Strings for file/path and path.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  pres(:) = 0.0 ; eta(:,:,:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "SPONGE_DAMPING_FILE", damping_file, &
                 "The name of the file with the sponge damping rates.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
                 "The name of the file with the state to damp toward.", &
                 default=damping_file)

  call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "SPONGE_STATE_FILE.", default="PTEMP")
  call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "SPONGE_STATE_FILE.", default="SALT")
  call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
                 "The name of the interface height variable in \n"//&
                 "SPONGE_STATE_FILE.", default="ETA")
  call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", idamp_var, &
                 "The name of the inverse damping rate variable in \n"//&
                 "SPONGE_DAMPING_FILE.", default="IDAMP")

  filename = trim(inputdir)//trim(damping_file)
  call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
  if (.not.file_exists(filename, g%Domain)) &
    call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))


  call read_data(filename,"Idamp",idamp(:,:), domain=g%Domain%mpp_domain)

! Now register all of the fields which are damped in the sponge.     !
! By default, momentum is advected vertically within the sponge, but !
! momentum is typically not damped within the sponge.                !

  filename = trim(inputdir)//trim(state_file)
  call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_FILE", filename)
  if (.not.file_exists(filename, g%Domain)) &
    call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))


!  The first call to set_up_sponge_field is for the interface height.!
  call read_data(filename, eta_var, eta(:,:,:), domain=g%Domain%mpp_domain)

  do j=js,je ; do i=is,ie
    eta(i,j,nz+1) = -g%bathyT(i,j)
  enddo ; enddo
  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) &
      eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
  enddo ; enddo ; enddo
! Set the inverse damping rates so that the model will know where to !
! apply the sponges, along with the interface heights.               !
  call initialize_sponge(idamp, eta, g, param_file, csp)

!   Now register all of the tracer fields which are damped in the    !
! sponge. By default, momentum is advected vertically within the     !
! sponge, but momentum is typically not damped within the sponge.    !

  if ( gv%nkml>0 ) then
!   This call to set_up_sponge_ML_density registers the target values of the
! mixed layer density, which is used in determining which layers can be
! inflated without causing static instabilities.
    do i=is-1,ie ; pres(i) = tv%P_Ref ; enddo

    call read_data(filename, potemp_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call read_data(filename, salin_var, tmp2(:,:,:), domain=g%Domain%mpp_domain)

    do j=js,je
      call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), &
                             is, ie-is+1, tv%eqn_of_state)
    enddo

    call set_up_sponge_ml_density(tmp_2d, g, csp)
  endif

!  The remaining calls to set_up_sponge_field can be in any order.   !
  if ( use_temperature ) then
    call read_data(filename, potemp_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%T, g, nz, csp)
    call read_data(filename, salin_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%S, g, nz, csp)
  endif


end subroutine initialize_sponges_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_velocity_depth_max(G)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! maximum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  enddo ; enddo
end subroutine set_velocity_depth_max
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine compute_global_grid_integrals(G)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
  ! Subroutine to pre-compute global integrals of grid quantities for
  ! later use in reporting diagnostics
  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
  integer :: i,j

  tmpforsumming(:,:) = 0.
  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    tmpforsumming(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
  enddo ; enddo
  g%areaT_global = reproducing_sum(tmpforsumming)
  g%IareaT_global = 1. / g%areaT_global
end subroutine compute_global_grid_integrals

! -----------------------------------------------------------------------------
subroutine set_velocity_depth_min(G)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! minimum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  enddo ; enddo
end subroutine set_velocity_depth_min
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine determines the isopycnal or other coordinate interfaces and
!! layer potential temperatures and salinities directly from a z-space file on
!! a latitude-longitude grid.
subroutine mom_temp_salt_initialize_from_z(h, tv, G, GV, PF, just_read_params)
! This subroutine was written by M. Harrison, with input from R. Hallberg & A. Adcroft.
!
! Arguments:
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (inout)   G       - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.

  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: h    !< Layer thicknesses being initialized, in m
  type(thermo_var_ptrs),   intent(inout) :: tv   !< A structure pointing to various thermodynamic
                                                 !! variables including temperature and salinity
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
  type(param_file_type),   intent(in)    :: PF   !< A structure indicating the open file
                                                 !! to parse for model parameter values.
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  character(len=200) :: filename   ! The name of an input file containing temperature
                                   ! and salinity in z-space; also used for  ice shelf area.
  character(len=200) :: shelf_file ! The name of an input file used for  ice shelf area.
  character(len=200) :: inputdir ! The directory where NetCDF input filesare.
  character(len=200) :: mesg, area_varname, ice_shelf_file

  type(eos_type), pointer :: eos => null()

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_initialize_layers_from_Z" ! This module's name.

  integer :: is, ie, js, je, nz ! compute domain indices
  integer :: isc,iec,jsc,jec    ! global compute domain indices
  integer :: isg, ieg, jsg, jeg ! global extent
  integer :: isd, ied, jsd, jed ! data domain indices

  integer :: i, j, k, ks, np, ni, nj
  integer :: idbg, jdbg
  integer :: nkml, nkbl         ! number of mixed and buffer layers

  integer :: kd, inconsistent
  real    :: PI_180             ! for conversion from degrees to radians

  real, dimension(:,:), pointer :: shelf_area
  real    :: min_depth
  real    :: dilate
  real    :: missing_value_temp, missing_value_salt
  logical :: correct_thickness
  character(len=40) :: potemp_var, salin_var
  character(len=8)  :: laynum

  integer, parameter :: niter=10   ! number of iterations for t/s adjustment to layer density
  logical :: just_read    ! If true, just read parameters but set nothing.
  logical            :: adjust_temperature = .true.  ! fit t/s to target densities
  real, parameter    :: missing_value = -1.e20
  real, parameter    :: temp_land_fill = 0.0, salt_land_fill = 35.0
  logical :: reentrant_x, tripolar_n,dbg
  logical :: debug = .false.  ! manually set this to true for verbose output

  !data arrays
  real, dimension(:), allocatable :: z_edges_in, z_in, Rb
  real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z
  real, dimension(:,:,:), allocatable :: rho_z
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi
  real, dimension(SZI_(G),SZJ_(G))  :: nlevs
  real, dimension(SZI_(G))   :: press


  ! Local variables for ALE remapping
  real, dimension(:), allocatable :: hTarget
  real, dimension(:,:), allocatable :: area_shelf_h
  real, dimension(:,:), allocatable, target  :: frac_shelf_h
  real, dimension(:,:,:), allocatable :: tmpT1dIn, tmpS1dIn, h1, tmp_mask_in
  real :: zTopOfCell, zBottomOfCell
  type(regridding_cs) :: regridCS ! Regridding parameters and work arrays
  type(remapping_cs) :: remapCS ! Remapping parameters and work arrays

  logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
  logical :: use_ice_shelf
  character(len=10) :: remappingScheme
  real :: tempAvg, saltAvg
  integer :: nPoints, ans
  integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE

  id_clock_routine = cpu_clock_id('(Initialize from Z)', grain=clock_routine)
  id_clock_ale = cpu_clock_id('(Initialize from Z) ALE', grain=clock_loop)

  call cpu_clock_begin(id_clock_routine)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg

  pi_180=atan(1.0)/45.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
  if (.not.just_read) call log_version(pf, mdl, version, "")

  inputdir = "." ;  call get_param(pf, mdl, "INPUTDIR", inputdir)
  inputdir = slasher(inputdir)

  eos => tv%eqn_of_state

! call mpp_get_compute_domain(G%domain%mpp_domain,isc,iec,jsc,jec)

  reentrant_x = .false. ; call get_param(pf, mdl, "REENTRANT_X", reentrant_x,default=.true.)
  tripolar_n = .false. ;  call get_param(pf, mdl, "TRIPOLAR_N", tripolar_n, default=.false.)
  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, default=0.0)

  call get_param(pf, mdl, "NKML",nkml,default=0)
  call get_param(pf, mdl, "NKBL",nkbl,default=0)

  call get_param(pf, mdl, "TEMP_SALT_Z_INIT_FILE",filename, &
                 "The name of the z-space input file used to initialize \n"//&
                 "the layer thicknesses, temperatures and salinities.", &
                 default="temp_salt_z.nc", do_not_log=just_read)
  filename = trim(inputdir)//trim(filename)
  call get_param(pf, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="ptemp", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
                 "If True, then horizontally homogenize the interpolated \n"//&
                 "initial conditions.", default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", usealeremapping, &
                 "If True, then remap straight to model coordinate from file.",&
                 default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remappingscheme, &
                 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING\n"//&
                 "is True.", default="PPM_IH4", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_GENERAL", remap_general, &
                 "If false, only initializes to z* coordinates.\n"//&
                 "If true, allows initialization directly to general coordinates.",&
                 default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
                 "If false, only reconstructs profiles for valid data points.\n"//&
                 "If true, inserts vanished layers below the valid data.",&
                 default=remap_general, do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
                 "If false, uses the preferred remapping algorithm for initialization.\n"//&
                 "If true, use an older, less robust algorithm for remapping.",&
                 default=.true., do_not_log=just_read)
  call get_param(pf, mdl, "ICE_SHELF", use_ice_shelf, default=.false.)
  if (use_ice_shelf) then
    call get_param(pf, mdl, "ICE_THICKNESS_FILE", ice_shelf_file, &
                 "The file from which the ice bathymetry and area are read.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
    shelf_file = trim(inputdir)//trim(ice_shelf_file)
    if (.not.just_read) call log_param(pf, mdl, "INPUTDIR/THICKNESS_FILE", shelf_file)
    call get_param(pf, mdl, "ICE_AREA_VARNAME", area_varname, &
                 "The name of the area variable in ICE_THICKNESS_FILE.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  endif
  if (.not.usealeremapping) then
    call get_param(pf, mdl, "ADJUST_THICKNESS", correct_thickness, &
                 "If true, all mass below the bottom removed if the \n"//&
                 "topography is shallower than the thickness input file \n"//&
                 "would indicate.", default=.false., do_not_log=just_read)

    call get_param(pf, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
                 "If true, all the interior layers are adjusted to \n"//&
                 "their target densities using mostly temperature \n"//&
                 "This approach can be problematic, particularly in the \n"//&
                 "high latitudes.", default=.true., do_not_log=just_read)
  endif
  if (just_read) then
    call cpu_clock_end(id_clock_routine)
    return ! All run-time parameters have been read, so return.
  endif

!   Read input grid coordinates for temperature and salinity field
!   in z-coordinate dataset. The file is REQUIRED to contain the
!   following:
!
!   dimension variables:
!            lon (degrees_E), lat (degrees_N), depth(meters)
!   variables:
!            ptemp(lon,lat,depth) : degC, potential temperature
!            salt (lon,lat,depth) : PSU, salinity
!
!   The first record will be read if there are multiple time levels.
!   The observation grid MUST tile the model grid. If the model grid extends
!   to the North/South Pole past the limits of the input data, they are extrapolated using the average
!   value at the northernmost/southernmost latitude.

  call horiz_interp_and_extrap_tracer(filename, potemp_var,1.0,1, &
       g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, tripolar_n, homogenize)

  call horiz_interp_and_extrap_tracer(filename, salin_var,1.0,1, &
       g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, tripolar_n, homogenize)

  kd = size(z_in,1)

  allocate(rho_z(isd:ied,jsd:jed,kd))
  allocate(area_shelf_h(isd:ied,jsd:jed))
  allocate(frac_shelf_h(isd:ied,jsd:jed))

  press(:)=tv%p_ref

  !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
  call convert_temp_salt_for_teos10(temp_z,salt_z, press, g, kd, mask_z, eos)

  do k=1,kd
    do j=js,je
      call calculate_density(temp_z(:,j,k),salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos)
    enddo
  enddo ! kd

  call pass_var(temp_z,g%Domain)
  call pass_var(salt_z,g%Domain)
  call pass_var(mask_z,g%Domain)
  call pass_var(rho_z,g%Domain)

  ! This is needed for building an ALE grid under ice shelves
  if (use_ice_shelf) then
    if (.not.file_exists(shelf_file, g%Domain)) call mom_error(fatal, &
      "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file))

    call read_data(shelf_file,trim(area_varname),area_shelf_h,domain=g%Domain%mpp_domain)

    ! initialize frac_shelf_h with zeros (open water everywhere)
    frac_shelf_h(:,:) = 0.0
    ! compute fractional ice shelf coverage of h
    do j=jsd,jed ; do i=isd,ied
      if (g%areaT(i,j) > 0.0) &
        frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
    enddo ; enddo
    ! pass to the pointer
    shelf_area => frac_shelf_h

  endif

! Done with horizontal interpolation.
! Now remap to model coordinates
  if (usealeremapping) then
    call cpu_clock_begin(id_clock_ale)
    ! The regridding tools (grid generation) are coded to work on model arrays of the same
    ! vertical shape. We need to re-write the regridding if the model has fewer layers
    ! than the data. -AJA
    if (kd>nz) call mom_error(fatal,"MOM_initialize_state, MOM_temp_salt_initialize_from_Z(): "//&
         "Data has more levels than the model - this has not been coded yet!")
    ! Build the source grid and copy data onto model-shaped arrays with vanished layers
    allocate( tmp_mask_in(isd:ied,jsd:jed,nz) ) ; tmp_mask_in(:,:,:) = 0.
    allocate( h1(isd:ied,jsd:jed,nz) ) ; h1(:,:,:) = 0.
    allocate( tmpt1din(isd:ied,jsd:jed,nz) ) ; tmpt1din(:,:,:) = 0.
    allocate( tmps1din(isd:ied,jsd:jed,nz) ) ; tmps1din(:,:,:) = 0.
    do j = js, je ; do i = is, ie
      if (g%mask2dT(i,j)>0.) then
        ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0
        tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
        do k = 1, nz
          if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then
            zbottomofcell = -min( z_edges_in(k+1), g%bathyT(i,j) )
            tmpt1din(i,j,k) = temp_z(i,j,k)
            tmps1din(i,j,k) = salt_z(i,j,k)
          elseif (k>1) then
            zbottomofcell = -g%bathyT(i,j)
            tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
            tmps1din(i,j,k) = tmps1din(i,j,k-1)
          else ! This next block should only ever be reached over land
            tmpt1din(i,j,k) = -99.9
            tmps1din(i,j,k) = -99.9
          endif
          h1(i,j,k) = ztopofcell - zbottomofcell
          if (h1(i,j,k)>0.) npoints = npoints + 1
          ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
        enddo
        h1(i,j,kd) = h1(i,j,kd) + ( ztopofcell + g%bathyT(i,j) ) ! In case data is deeper than model
      endif ! mask2dT
    enddo ; enddo
    deallocate( tmp_mask_in )
    call pass_var(h1, g%Domain)
    call pass_var(tmpt1din, g%Domain)
    call pass_var(tmps1din, g%Domain)

    ! Build the target grid (and set the model thickness to it)
    ! This call can be more general but is hard-coded for z* coordinates...  ????
    call ale_initregridding( gv, g%max_depth, pf, mdl, regridcs ) ! sets regridCS

    if (.not. remap_general) then
      ! This is the old way of initializing to z* coordinates only
      allocate( htarget(nz) )
      htarget = getcoordinateresolution( regridcs )
      do j = js, je ; do i = is, ie
        h(i,j,:) = 0.
        if (g%mask2dT(i,j)>0.) then
          ! Build the target grid combining hTarget and topography
          ztopofcell = 0. ; zbottomofcell = 0.
          do k = 1, nz
            zbottomofcell = max( ztopofcell - htarget(k), -g%bathyT(i,j) )
            h(i,j,k) = ztopofcell - zbottomofcell
            ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
          enddo
        else
          h(i,j,:) = 0.
        endif ! mask2dT
      enddo ; enddo
      call pass_var(h, g%Domain)
      deallocate( htarget )
    endif

    ! Now remap from source grid to target grid
    call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false. ) ! Reconstruction parameters
    if (remap_general) then
      call set_regrid_params( regridcs, min_thickness=0. )
      h(:,:,:) = h1(:,:,:) ; tv%T(:,:,:) = tmpt1din(:,:,:) ; tv%S(:,:,:) = tmps1din(:,:,:)
      do j = js, je ; do i = is, ie
        if (g%mask2dT(i,j)==0.) then ! Ensure there are no nonsense values on land
          h(i,j,:) = 0. ; tv%T(i,j,:) = 0. ; tv%S(i,j,:) = 0.
        endif
      enddo ; enddo
      call pass_var(h, g%Domain)    ! Regridding might eventually use spatial information and
      call pass_var(tv%T, g%Domain) ! thus needs to be up to date in the halo regions even though
      call pass_var(tv%S, g%Domain) ! ALE_build_grid() only updates h on the computational domain.

      if (use_ice_shelf) then
         call ale_build_grid( g, gv, regridcs, remapcs, h, tv, .true., shelf_area)
      else
         call ale_build_grid( g, gv, regridcs, remapcs, h, tv, .true. )
      endif
    endif
    call ale_remap_scalar( remapcs, g, gv, nz, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, old_remap=remap_old_alg )
    call ale_remap_scalar( remapcs, g, gv, nz, h1, tmps1din, h, tv%S, all_cells=remap_full_column, old_remap=remap_old_alg )
    deallocate( h1 )
    deallocate( tmpt1din )
    deallocate( tmps1din )

    call cpu_clock_end(id_clock_ale)

  else ! remap to isopycnal layer space

! next find interface positions using local arrays
! nlevs contains the number of valid data points in each column
    nlevs = sum(mask_z,dim=3)

! Rb contains the layer interface densities
    allocate(rb(nz+1))
    do k=2,nz ; rb(k)=0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
    rb(1) = 0.0 ;  rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)

    zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, rb, g%bathyT(is:ie,js:je), &
                         nlevs(is:ie,js:je), nkml, nkbl, min_depth)

    if (correct_thickness) then
      call adjustetatofitbathymetry(g, gv, zi, h)
    else
      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
        if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_z)) then
          zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_z
          h(i,j,k) = gv%Angstrom_z
        else
          h(i,j,k) = zi(i,j,k) - zi(i,j,k+1)
        endif
      enddo ; enddo ; enddo
      inconsistent=0
      do j=js,je ; do i=is,ie
        if (abs(zi(i,j,nz+1) + g%bathyT(i,j)) > 1.0) &
          inconsistent = inconsistent + 1
      enddo ; enddo
      call sum_across_pes(inconsistent)

      if ((inconsistent > 0) .and. (is_root_pe())) then
        write(mesg,'("Thickness initial conditions are inconsistent ",'// &
                 '"with topography in ",I5," places.")') inconsistent
        call mom_error(warning, mesg)
      endif
    endif

    tv%T(is:ie,js:je,:) = tracer_z_init(temp_z(is:ie,js:je,:),-1.0*z_edges_in,zi(is:ie,js:je,:), &
                                        nkml,nkbl,missing_value,g%mask2dT(is:ie,js:je),nz, &
                                        nlevs(is:ie,js:je),dbg,idbg,jdbg)
    tv%S(is:ie,js:je,:) = tracer_z_init(salt_z(is:ie,js:je,:),-1.0*z_edges_in,zi(is:ie,js:je,:), &
                                        nkml,nkbl,missing_value,g%mask2dT(is:ie,js:je),nz, &
                                        nlevs(is:ie,js:je))

    do k=1,nz
      npoints = 0 ; tempavg = 0. ; saltavg = 0.
      do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) >= 1.0) then
        npoints = npoints + 1
        tempavg = tempavg + tv%T(i,j,k)
        saltavg =saltavg + tv%S(i,j,k)
      endif ; enddo ; enddo

   ! Horizontally homogenize data to produce perfectly "flat" initial conditions
      if (homogenize) then
        call sum_across_pes(npoints)
        call sum_across_pes(tempavg)
        call sum_across_pes(saltavg)
        if (npoints>0) then
          tempavg = tempavg/real(npoints)
          saltavg = saltavg/real(npoints)
        endif
        tv%T(:,:,k) = tempavg
        tv%S(:,:,k) = saltavg
      endif
    enddo

  endif ! useALEremapping

! Fill land values
  do k=1,nz ; do j=js,je ; do i=is,ie
    if (tv%T(i,j,k) == missing_value) then
      tv%T(i,j,k)=temp_land_fill
      tv%S(i,j,k)=salt_land_fill
    endif
  enddo ; enddo ; enddo

! Finally adjust to target density
  ks=max(0,nkml)+max(0,nkbl)+1

  if (adjust_temperature .and. .not. usealeremapping) then
    call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
            gv%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos)

  endif

  deallocate(z_in,z_edges_in,temp_z,salt_z,mask_z)

  call pass_var(h, g%Domain)
  call pass_var(tv%T, g%Domain)
  call pass_var(tv%S, g%Domain)

  call calltree_leave(trim(mdl)//'()')
  call cpu_clock_end(id_clock_routine)

end subroutine mom_temp_salt_initialize_from_z

!> Run simple unit tests
subroutine mom_state_init_tests(G, GV, tv)
  type(ocean_grid_type),     intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),   intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(thermo_var_ptrs),     intent(in)    :: tv   !< Thermodynamics structure.
  ! Local variables
  integer, parameter :: nk=5
  real, dimension(nk) :: T, T_t, T_b, S, S_t, S_b, rho, h, z
  real, dimension(nk+1) :: e
  integer :: k
  real :: P_tot, P_t, P_b, z_out
  type(remapping_cs), pointer :: remap_CS => null()

  do k = 1, nk
    h(k) = 100.
  enddo
  e(1) = 0.
  do k = 1, nk
    e(k+1) = e(k) - h(k)
  enddo
  p_tot = 0.
  do k = 1, nk
    z(k) = 0.5 * ( e(k) + e(k+1) )
    t_t(k) = 20.+(0./500.)*e(k)
    t(k)   = 20.+(0./500.)*z(k)
    t_b(k) = 20.+(0./500.)*e(k+1)
    s_t(k) = 35.-(0./500.)*e(k)
    s(k)   = 35.+(0./500.)*z(k)
    s_b(k) = 35.-(0./500.)*e(k+1)
    call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -gv%Rho0*gv%g_Earth*z(k), rho(k), tv%eqn_of_state)
    p_tot = p_tot + gv%g_Earth * rho(k) * h(k)
  enddo

  p_t = 0.
  do k = 1, nk
    call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
                                        p_t, 0.5*p_tot, gv%Rho0, gv%g_Earth, tv%eqn_of_state, p_b, z_out)
    write(0,*) k,p_t,p_b,0.5*p_tot,e(k),e(k+1),z_out
    p_t = p_b
  enddo
  write(0,*) p_b,p_tot

  write(0,*) ''
  write(0,*) ' ==================================================================== '
  write(0,*) ''
  write(0,*) h
  call cut_off_column_top(nk, tv, gv%Rho0, gv%g_Earth, -e(nk+1), gv%Angstrom, &
               t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)
  write(0,*) h

end subroutine mom_state_init_tests

end module mom_state_initialization