MOM_regridding.F90

Go to the documentation of this file.

!> Generates vertical grids as part of the ALE algorithm
module mom_regridding

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

use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser,   only : param_file_type, get_param, log_param
use mom_io,            only : file_exists, field_exists, field_size, mom_read_data
use mom_io,            only : vardesc, var_desc, fieldtype, single_file
use mom_io,            only : create_file, write_field, close_file, slasher
use mom_variables,     only : ocean_grid_type, thermo_var_ptrs
use mom_verticalgrid,  only : verticalgrid_type
use mom_eos,           only : eos_type, calculate_density
use mom_string_functions,only : uppercase, extractword, extract_integer, extract_real

use mom_remapping, only : remapping_cs
use regrid_consts, only : state_dependent, coordinateunits
use regrid_consts, only : coordinatemode, default_coordinate_mode
use regrid_consts, only : regridding_layer, regridding_zstar
use regrid_consts, only : regridding_rho, regridding_sigma
use regrid_consts, only : regridding_arbitrary, regridding_sigma_shelf_zstar
use regrid_consts, only : regridding_hycom1, regridding_slight, regridding_adaptive
use regrid_interp, only : interp_cs_type, set_interp_scheme, set_interp_extrap

use coord_zlike,  only : init_coord_zlike, zlike_cs, set_zlike_params, build_zstar_column, end_coord_zlike
use coord_sigma,  only : init_coord_sigma, sigma_cs, set_sigma_params, build_sigma_column, end_coord_sigma
use coord_rho,    only : init_coord_rho, rho_cs, set_rho_params, build_rho_column, end_coord_rho
use coord_rho,    only : old_inflate_layers_1d
use coord_hycom,  only : init_coord_hycom, hycom_cs, set_hycom_params, build_hycom1_column, end_coord_hycom
use coord_slight, only : init_coord_slight, slight_cs, set_slight_params, build_slight_column, end_coord_slight
use coord_adapt,  only : init_coord_adapt, adapt_cs, set_adapt_params, build_adapt_column, end_coord_adapt

use netcdf ! Used by check_grid_def()

implicit none ; private

#include <MOM_memory.h>

!> Regridding control structure
type, public :: regridding_cs
  private

  !> This array is set by function setCoordinateResolution()
  !! It contains the "resolution" or delta coordinate of the target
  !! coorindate. It has the units of the target coordiante, e.g.
  !! meters for z*, non-dimensional for sigma, etc.
  real, dimension(:), allocatable :: coordinateresolution

  !> This array is set by function set_target_densities()
  !! This array is the nominal coordinate of interfaces and is the
  !! running sum of coordinateResolution. i.e.
  !!  target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
  !! It is only used in "rho" mode.
  real, dimension(:), allocatable :: target_density
  logical :: target_density_set = .false.

  !> This array is set by function set_regrid_max_depths()
  !! It specifies the maximum depth that every interface is allowed to take, in H.
  real, dimension(:), allocatable :: max_interface_depths

  !> This array is set by function set_regrid_max_thickness()
  !! It specifies the maximum depth that every interface is allowed to take, in H.
  real, dimension(:), allocatable :: max_layer_thickness

  integer :: nk !< Number of layers/levels

  !> Indicates which grid to use in the vertical (z*, sigma, target interface
  !! densities)
  integer :: regridding_scheme

  !> Interpolation control structure
  type(interp_cs_type) :: interp_cs

  !> Minimum thickness allowed when building the new grid through regridding
  real :: min_thickness

  !> Reference pressure for potential density calculations (Pa)
  real :: ref_pressure = 2.e7

  !> Weight given to old coordinate when blending between new and old grids (nondim)
  !! Used only below depth_of_time_filter_shallow, with a cubic variation
  !! from zero to full effect between depth_of_time_filter_shallow and
  !! depth_of_time_filter_deep.
  real :: old_grid_weight = 0.

  !> Depth above which no time-filtering of grid is applied (H units)
  real :: depth_of_time_filter_shallow = 0.

  !> Depth below which time-filtering of grid is applied at full effect (H units)
  real :: depth_of_time_filter_deep = 0.

  !> Fraction (between 0 and 1) of compressibility to add to potential density
  !! profiles when interpolating for target grid positions. (nondim)
  real :: compressibility_fraction = 0.

  !> If true, each interface is given a maximum depth based on a rescaling of
  !! the indexing of coordinateResolution.
  logical :: set_maximum_depths = .false.

  !> A scaling factor (> 1) of the rate at which the coordinateResolution list
  !! is traversed to set the minimum depth of interfaces.
  real :: max_depth_index_scale = 2.0

  !> If true, integrate for interface positions from the top downward.
  !! If false, integrate from the bottom upward, as does the rest of the model.
  logical :: integrate_downward_for_e = .true.

  type(zlike_cs),  pointer :: zlike_cs  => null()
  type(sigma_cs),  pointer :: sigma_cs  => null()
  type(rho_cs),    pointer :: rho_cs    => null()
  type(hycom_cs),  pointer :: hycom_cs  => null()
  type(slight_cs), pointer :: slight_cs => null()
  type(adapt_cs),  pointer :: adapt_cs  => null()

end type

! The following routines are visible to the outside world
public initialize_regridding, end_regridding, regridding_main
public inflate_vanished_layers_old, check_remapping_grid, check_grid_column
public adjust_interface_motion
public set_regrid_params, get_regrid_size
public uniformresolution, setcoordinateresolution
public build_rho_column
public set_target_densities_from_gv, set_target_densities
public set_regrid_max_depths, set_regrid_max_thickness
public getcoordinateresolution, getcoordinateinterfaces
public getcoordinateunits, getcoordinateshortname, getstaticthickness
public default_coordinate_mode
public get_zlike_cs, get_sigma_cs, get_rho_cs

!> Documentation for coordinate options
character(len=*), parameter, public :: regriddingcoordinatemodedoc = &
                 " LAYER - Isopycnal or stacked shallow water layers\n"//&
                 " ZSTAR, Z* - stetched geopotential z*\n"//&
                 " SIGMA_SHELF_ZSTAR - stetched geopotential z* ignoring shelf\n"//&
                 " SIGMA - terrain following coordinates\n"//&
                 " RHO   - continuous isopycnal\n"//&
                 " HYCOM1 - HyCOM-like hybrid coordinate\n"//&
                 " SLIGHT - stretched coordinates above continuous isopycnal\n"//&
                 " ADAPTIVE - optimize for smooth neutral density surfaces"

! Documentation for regridding interpolation schemes
character(len=*), parameter, public :: regriddinginterpschemedoc = &
                 " P1M_H2     (2nd-order accurate)\n"//&
                 " P1M_H4     (2nd-order accurate)\n"//&
                 " P1M_IH4    (2nd-order accurate)\n"//&
                 " PLM        (2nd-order accurate)\n"//&
                 " PPM_H4     (3rd-order accurate)\n"//&
                 " PPM_IH4    (3rd-order accurate)\n"//&
                 " P3M_IH4IH3 (4th-order accurate)\n"//&
                 " P3M_IH6IH5 (4th-order accurate)\n"//&
                 " PQM_IH4IH3 (4th-order accurate)\n"//&
                 " PQM_IH6IH5 (5th-order accurate)"
character(len=*), parameter, public :: regriddingdefaultinterpscheme = "P1M_H2"
logical, parameter, public :: regriddingdefaultboundaryextrapolation = .false.
real, parameter, public :: regriddingdefaultminthickness = 1.e-3

#undef __DO_SAFETY_CHECKS__

contains

!> Initialization and configures a regridding control structure based on customizable run-time parameters
subroutine initialize_regridding(CS, GV, max_depth, param_file, mod, coord_mode, param_prefix, param_suffix)
  type(regridding_cs),        intent(inout) :: CS !< Regridding control structure
  type(verticalgrid_type),    intent(in)    :: GV         !< Ocean vertical grid structure
  real,                       intent(in)    :: max_depth  !< The maximum depth of the ocean, in m.
  type(param_file_type),      intent(in)    :: param_file !< Parameter file
  character(len=*),           intent(in)    :: mod        !< Name of calling module.
  character(len=*),           intent(in)    :: coord_mode !< Coordinate mode
  character(len=*),           intent(in)    :: param_prefix !< String to prefix to parameter names.
                                                            !! If empty, causes main model parameters to be used.
  character(len=*),           intent(in)    :: param_suffix !< String to append to parameter names.
  ! Local variables
  integer :: ke ! Number of levels
  character(len=80)  :: string, string2, varName ! Temporary strings
  character(len=40)  :: coord_units, param_name, coord_res_param ! Temporary strings
  character(len=200) :: inputdir, fileName
  character(len=320) :: message ! Temporary strings
  character(len=12) :: expected_units ! Temporary strings
  logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters
  logical :: coord_is_state_dependent, ierr
  real :: filt_len, strat_tol, index_scale, tmpReal
  real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
  real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
  integer :: nz_fixed_sfc, k, nzf(4)
  real, dimension(:), allocatable :: dz     ! Resolution (thickness) in units of coordinate
  real, dimension(:), allocatable :: h_max  ! Maximum layer thicknesses, in m.
  real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths, in m.
  real, dimension(:), allocatable :: z_max  ! Maximum interface depths, in m.
  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode
  ! Thicknesses that give level centers corresponding to table 2 of WOA09
  real, dimension(40) :: woa09_dz = (/ 5.,  10.,  10.,  15.,  22.5, 25., 25.,  25.,  &
                                      37.5, 50.,  50.,  75., 100., 100., 100., 100., &
                                     100., 100., 100., 100., 100., 100., 100., 175., &
                                     250., 375., 500., 500., 500., 500., 500., 500., &
                                     500., 500., 500., 500., 500., 500., 500., 500. /)

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

  main_parameters=.false.
  if (len_trim(param_prefix)==0) main_parameters=.true.
  if (main_parameters .and. len_trim(param_suffix)>0) call mom_error(fatal,trim(mod)//', initialize_regridding: '// &
              'Suffix provided without prefix for parameter names!')

  cs%nk = 0
  cs%regridding_scheme = coordinatemode(coord_mode)
  coord_is_state_dependent = state_dependent(coord_mode)

  if (main_parameters) then
    ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
    call get_param(param_file, mod, "REGRIDDING_COORDINATE_UNITS", coord_units, &
                 "Units of the regridding coordinuate.",&
                 default=coordinateunits(coord_mode))
  else
    coord_units=coordinateunits(coord_mode)
  endif

  if (coord_is_state_dependent) then
    if (main_parameters) then
      param_name = "INTERPOLATION_SCHEME"
      string2 = regriddingdefaultinterpscheme
    else
      param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
      string2 = 'PPM_H4' ! Default for diagnostics
    endif
    call get_param(param_file, mod, "INTERPOLATION_SCHEME", string, &
                 "This sets the interpolation scheme to use to\n"//&
                 "determine the new grid. These parameters are\n"//&
                 "only relevant when REGRIDDING_COORDINATE_MODE is\n"//&
                 "set to a function of state. Otherwise, it is not\n"//&
                 "used. It can be one of the following schemes:\n"//&
                 trim(regriddinginterpschemedoc), default=trim(string2))
    call set_regrid_params(cs, interp_scheme=string)
  endif

  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mod, "BOUNDARY_EXTRAPOLATION", tmplogical, &
                 "When defined, a proper high-order reconstruction\n"//&
                 "scheme is used within boundary cells rather\n"//&
                 "than PCM. E.g., if PPM is used for remapping, a\n"//&
                 "PPM reconstruction will also be used within\n"//&
                 "boundary cells.", default=regriddingdefaultboundaryextrapolation)
    call set_regrid_params(cs, boundary_extrapolation=tmplogical)
  else
    call set_regrid_params(cs, boundary_extrapolation=.false.)
  endif

  ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
  if (main_parameters) then
    param_name = "ALE_COORDINATE_CONFIG"
    coord_res_param = "ALE_RESOLUTION"
    string2 = 'UNIFORM'
  else
    param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
    coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
    string2 = 'UNIFORM'
    if (max_depth>3000.) string2='WOA09' ! For convenience
  endif
  call get_param(param_file, mod, param_name, string, &
                 "Determines how to specify the coordinate\n"//&
                 "resolution. Valid options are:\n"//&
                 " PARAM       - use the vector-parameter "//trim(coord_res_param)//"\n"//&
                 " UNIFORM[:N] - uniformly distributed\n"//&
                 " FILE:string - read from a file. The string specifies\n"//&
                 "               the filename and variable name, separated\n"//&
                 "               by a comma or space, e.g. FILE:lev.nc,dz\n"//&
                 "               or FILE:lev.nc,interfaces=zw\n"//&
                 " WOA09[:N]   - the WOA09 vertical grid (approximately)\n"//&
                 " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
                 " HYBRID:string - read from a file. The string specifies\n"//&
                 "               the filename and two variable names, separated\n"//&
                 "               by a comma or space, for sigma-2 and dz. e.g.\n"//&
                 "               HYBRID:vgrid.nc,sigma2,dz",&
                 default=trim(string2))
  message = "The distribution of vertical resolution for the target\n"//&
            "grid used for Eulerian-like coordinates. For example,\n"//&
            "in z-coordinate mode, the parameter is a list of level\n"//&
            "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
            "is of non-dimensional fractions of the water column."
  if (index(trim(string),'UNIFORM')==1) then
    if (len_trim(string)==7) then
      ke = gv%ke ! Use model nk by default
      tmpreal = max_depth
    elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
      ! Format is "UNIFORM:N" or "UNIFORM:N,dz"
      ke = extract_integer(string(9:len_trim(string)),'',1)
      tmpreal = extract_real(string(9:len_trim(string)),',',2,missing_value=max_depth)
    else
      call mom_error(fatal,trim(mod)//', initialize_regridding: '// &
          'Unable to interpret "'//trim(string)//'".')
    endif
    allocate(dz(ke))
    if (ke==1) then
      dz(:) = uniformresolution(ke, coord_mode, tmpreal, gv%Rlay(1), gv%Rlay(1))
    else
      dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
                   gv%Rlay(1)+0.5*(gv%Rlay(1)-gv%Rlay(2)), &
                   gv%Rlay(ke)+0.5*(gv%Rlay(ke)-gv%Rlay(ke-1)) )
    endif
    if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, &
                   trim(message), units=trim(coord_units))
  elseif (trim(string)=='PARAM') then
    ! Read coordinate resolution (main model = ALE_RESOLUTION)
    ke = gv%ke ! Use model nk by default
    allocate(dz(ke))
    call get_param(param_file, mod, coord_res_param, dz, &
                   trim(message), units=trim(coord_units), fail_if_missing=.true.)
  elseif (index(trim(string),'FILE:')==1) then
    ! FILE:filename,var_name is assumed to be reading level thickness variables
    ! FILE:filename,interfaces=var_name reads positions
    if (string(6:6)=='.' .or. string(6:6)=='/') then
      ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
      filename = trim( extractword(trim(string(6:80)), 1) )
    else
      ! Otherwise assume we should look for the file in INPUTDIR
      filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
    endif
    if (.not. file_exists(filename)) call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
            "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")

    varname = trim( extractword(trim(string(6:)), 2) )
    if (len_trim(varname)==0) then
      if (field_exists(filename,'dz')) then; varname = 'dz'
      elseif (field_exists(filename,'dsigma')) then; varname = 'dsigma'
      elseif (field_exists(filename,'ztest')) then; varname = 'ztest'
      else ;  call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
                    "Coordinate variable not specified and none could be guessed.")
      endif
    endif
    ! This check fails when the variable is a dimension variable! -AJA
   !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// &
   !             "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
    if (cs%regridding_scheme == regridding_sigma) then
      expected_units = 'nondim'
    elseif (cs%regridding_scheme == regridding_rho) then
      expected_units = 'kg m-3'
    else
      expected_units = 'meters'
    endif
    if (index(trim(varname),'interfaces=')==1) then
      varname=trim(varname(12:))
      call check_grid_def(filename, varname, expected_units, message, ierr)
      if (ierr) call mom_error(fatal,trim(mod)//", initialize_regridding: "//&
                  "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message))
      call field_size(trim(filename), trim(varname), nzf)
      ke = nzf(1)-1
      if (cs%regridding_scheme == regridding_rho) then
        allocate(rho_target(ke+1))
        call mom_read_data(trim(filename), trim(varname), rho_target)
      else
        allocate(dz(ke))
        allocate(z_max(ke+1))
        call mom_read_data(trim(filename), trim(varname), z_max)
        dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
        deallocate(z_max)
      endif
    else
      ! Assume reading resolution
      call field_size(trim(filename), trim(varname), nzf)
      ke = nzf(1)
      allocate(dz(ke))
      call mom_read_data(trim(filename), trim(varname), dz)
    endif
    if (main_parameters .and. ke/=gv%ke) then
      call mom_error(fatal,trim(mod)//', initialize_regridding: '// &
                 'Mismatch in number of model levels and "'//trim(string)//'".')
    endif
    if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
  elseif (index(trim(string),'FNC1:')==1) then
    ke = gv%ke; allocate(dz(ke))
    call dz_function1( trim(string(6:)), dz )
    if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
  elseif (index(trim(string),'HYBRID:')==1) then
    ke = gv%ke; allocate(dz(ke))
    ! The following assumes the FILE: syntax of above but without "FILE:" in the string
    allocate(rho_target(ke+1))
    filename = trim( extractword(trim(string(8:)), 1) )
    if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
    if (.not. file_exists(filename)) call mom_error(fatal,trim(mod)//", initialize_regridding: HYBRID "// &
      "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
    varname = trim( extractword(trim(string(8:)), 2) )
    if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mod)//", initialize_regridding: HYBRID "// &
      "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
    call mom_read_data(trim(filename), trim(varname), rho_target)
    varname = trim( extractword(trim(string(8:)), 3) )
    if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
      call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
    else ! Read dz from file
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mod)//", initialize_regridding: HYBRID "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      call mom_read_data(trim(filename), trim(varname), dz)
    endif
    if (main_parameters) then
      call log_param(param_file, mod, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
      call log_param(param_file, mod, "!TARGET_DENSITIES", rho_target, &
               'HYBRID target densities for itnerfaces', units=coordinateunits(coord_mode))
    endif
  elseif (index(trim(string),'WOA09')==1) then
    if (len_trim(string)==5) then
      tmpreal = 0. ; ke = 0
      do while (tmpreal<max_depth)
        ke = ke + 1
        tmpreal = tmpreal + woa09_dz(ke)
      enddo
    elseif (index(trim(string),'WOA09:')==1) then
      if (len_trim(string)==6) call mom_error(fatal,trim(mod)//', initialize_regridding: '// &
                 'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
      ke = extract_integer(string(7:len_trim(string)),'',1)
    endif
    if (ke>40 .or. ke<1) call mom_error(fatal,trim(mod)//', initialize_regridding: '// &
                 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
    allocate(dz(ke))
    dz(1:ke) = woa09_dz(1:ke)
  else
    call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
      "Unrecognized coordinate configuraiton"//trim(string))
  endif

  if (main_parameters) then
    ! This is a work around to apparently needed to work with the from_Z initialization...  ???
    if (coordinatemode(coord_mode) == regridding_zstar .or. &
        coordinatemode(coord_mode) == regridding_hycom1 .or. &
        coordinatemode(coord_mode) == regridding_slight .or. &
        coordinatemode(coord_mode) == regridding_adaptive) then
      ! Adjust target grid to be consistent with max_depth
      tmpreal = sum( dz(:) )
      if (tmpreal < max_depth) then
        dz(ke) = dz(ke) + ( max_depth - tmpreal )
      elseif (tmpreal > max_depth) then
        if ( dz(ke) + ( max_depth - tmpreal ) > 0. ) then
          dz(ke) = dz(ke) + ( max_depth - tmpreal )
        else
          call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
            "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
        endif
      endif
    endif
  endif

  cs%nk=ke

  ! Target resolution (for fixed coordinates)
  allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
  if (state_dependent(cs%regridding_scheme)) then
    ! Target values
    allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30
  endif

  if (allocated(dz)) call setcoordinateresolution(dz, cs)

  if (allocated(rho_target)) then
    call set_target_densities(cs, rho_target)
    deallocate(rho_target)

  ! \todo This line looks like it would overwrite the target densities set just above?
  elseif (coordinatemode(coord_mode) == regridding_rho) then
    call set_target_densities_from_gv(gv, cs)
    call log_param(param_file, mod, "!TARGET_DENSITIES", cs%target_density, &
             'RHO target densities for interfaces', units=coordinateunits(coord_mode))
  endif

  ! initialise coordinate-specific control structure
  call initcoord(cs, coord_mode)

  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mod, "REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
                 "When interpolating potential density profiles we can add\n"//&
                 "some artificial compressibility solely to make homogenous\n"//&
                 "regions appear stratified.", default=0.)
    call set_regrid_params(cs, compress_fraction=tmpreal)
  endif

  if (main_parameters) then
    call get_param(param_file, mod, "MIN_THICKNESS", tmpreal, &
                 "When regridding, this is the minimum layer\n"//&
                 "thickness allowed.", units="m",&
                 default=regriddingdefaultminthickness )
    call set_regrid_params(cs, min_thickness=tmpreal)
  else
    call set_regrid_params(cs, min_thickness=0.)
  endif

  if (coordinatemode(coord_mode) == regridding_slight) then
    ! Set SLight-specific regridding parameters.
    call get_param(param_file, mod, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
                 "The nominal thickness of fixed thickness near-surface\n"//&
                 "layers with the SLight coordinate.", units="m", default=1.0)
    call get_param(param_file, mod, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
                 "The number of fixed-depth surface layers with the SLight\n"//&
                 "coordinate.", units="nondimensional", default=2)
    call get_param(param_file, mod, "SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
                 "The thickness of the surface region over which to average\n"//&
                 "when calculating the density to use to define the interior\n"//&
                 "with the SLight coordinate.", units="m", default=1.0)
    call get_param(param_file, mod, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
                 "The number of layers to offset the surface density when\n"//&
                 "defining where the interior ocean starts with SLight.", &
                 units="nondimensional", default=2.0)
    call get_param(param_file, mod, "SLIGHT_FIX_HALOCLINES", fix_haloclines, &
                 "If true, identify regions above the reference pressure\n"//&
                 "where the reference pressure systematically underestimates\n"//&
                 "the stratification and use this in the definition of the\n"//&
                 "interior with the SLight coordinate.", default=.false.)

    call set_regrid_params(cs, dz_min_surface=dz_fixed_sfc, &
                nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
                nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
    if (fix_haloclines) then
      ! Set additional parameters related to SLIGHT_FIX_HALOCLINES.
      call get_param(param_file, mod, "HALOCLINE_FILTER_LENGTH", filt_len, &
                 "A length scale over which to smooth the temperature and\n"//&
                 "salinity before identifying erroneously unstable haloclines.", &
                 units="m", default=2.0)
      call get_param(param_file, mod, "HALOCLINE_STRAT_TOL", strat_tol, &
                 "A tolerance for the ratio of the stratification of the\n"//&
                 "apparent coordinate stratification to the actual value\n"//&
                 "that is used to identify erroneously unstable haloclines.\n"//&
                 "This ratio is 1 when they are equal, and sensible values \n"//&
                 "are between 0 and 0.5.", units="nondimensional", default=0.2)
      call set_regrid_params(cs, halocline_filt_len=filt_len, &
                             halocline_strat_tol=strat_tol)
    endif

  endif

  if (coordinatemode(coord_mode) == regridding_adaptive) then
    call get_param(param_file, mod, "ADAPT_TIME_RATIO", adapttimeratio, &
         "Ratio of ALE timestep to grid timescale.", units="s", default=1e-1)
    call get_param(param_file, mod, "ADAPT_ZOOM_DEPTH", adaptzoom, &
         "Depth of near-surface zooming region.", units="m", default=200.0)
    call get_param(param_file, mod, "ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
         "Coefficient of near-surface zooming diffusivity.", &
         units="nondim", default=0.2)
    call get_param(param_file, mod, "ADAPT_BUOY_COEFF", adaptbuoycoeff, &
         "Coefficient of buoyancy diffusivity.", &
         units="nondim", default=0.8)
    call get_param(param_file, mod, "ADAPT_ALPHA", adaptalpha, &
         "Scaling on optimization tendency.", &
         units="nondim", default=1.0)
    call get_param(param_file, mod, "ADAPT_DO_MIN_DEPTH", tmplogical, &
         "If true, make a HyCOM-like mixed layer by preventing interfaces\n"//&
         "from being shallower than the depths specified by the regridding coordinate.", &
         default=.false.)

    call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
         adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
         adaptdomin=tmplogical)
  endif

  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mod, "MAXIMUM_INT_DEPTH_CONFIG", string, &
                 "Determines how to specify the maximum interface depths.\n"//&
                 "Valid options are:\n"//&
                 " NONE        - there are no maximum interface depths\n"//&
                 " PARAM       - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
                 " FILE:string - read from a file. The string specifies\n"//&
                 "               the filename and variable name, separated\n"//&
                 "               by a comma or space, e.g. FILE:lev.nc,Z\n"//&
                 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
                 default='NONE')
    message = "The list of maximum depths for each interface."
    allocate(z_max(ke+1))
    allocate(dz_max(ke))
    if ( trim(string) == "NONE") then
      ! Do nothing.
    elseif ( trim(string) ==  "PARAM") then
      call get_param(param_file, mod, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
                   trim(message), units="m", fail_if_missing=.true.)
      call set_regrid_max_depths(cs, z_max, gv%m_to_H)
    elseif (index(trim(string),'FILE:')==1) then
      if (string(6:6)=='.' .or. string(6:6)=='/') then
        ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
        filename = trim( extractword(trim(string(6:80)), 1) )
      else
        ! Otherwise assume we should look for the file in INPUTDIR
        filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
      endif
      if (.not. file_exists(filename)) call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")

      do_sum = .false.
      varname = trim( extractword(trim(string(6:)), 2) )
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      if (len_trim(varname)==0) then
        if (field_exists(filename,'z_max')) then; varname = 'z_max'
        elseif (field_exists(filename,'dz')) then; varname = 'dz' ; do_sum = .true.
        elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max' ; do_sum = .true.
        else ; call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
          "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
        endif
      endif
      if (do_sum) then
        call mom_read_data(trim(filename), trim(varname), dz_max)
        z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
      else
        call mom_read_data(trim(filename), trim(varname), z_max)
      endif
      call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_depths(cs, z_max, gv%m_to_H)
    elseif (index(trim(string),'FNC1:')==1) then
      call dz_function1( trim(string(6:)), dz_max )
      if ((coordinatemode(coord_mode) == regridding_slight) .and. &
          (dz_fixed_sfc > 0.0)) then
        do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo
      endif
      z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
      call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_depths(cs, z_max, gv%m_to_H)
    else
      call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
    endif
    deallocate(z_max)
    deallocate(dz_max)

    ! Optionally specify maximum thicknesses for each layer, enforced by moving
    ! the interface below a layer downward.
    call get_param(param_file, mod, "MAX_LAYER_THICKNESS_CONFIG", string, &
                   "Determines how to specify the maximum layer thicknesses.\n"//&
                   "Valid options are:\n"//&
                   " NONE        - there are no maximum layer thicknesses\n"//&
                   " PARAM       - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
                   " FILE:string - read from a file. The string specifies\n"//&
                   "               the filename and variable name, separated\n"//&
                   "               by a comma or space, e.g. FILE:lev.nc,Z\n"//&
                   " FNC1:string - FNC1:dz_min,H_total,power,precision",&
                   default='NONE')
    message = "The list of maximum thickness for each layer."
    allocate(h_max(ke))
    if ( trim(string) == "NONE") then
      ! Do nothing.
    elseif ( trim(string) ==  "PARAM") then
      call get_param(param_file, mod, "MAX_LAYER_THICKNESS", h_max, &
                   trim(message), units="m", fail_if_missing=.true.)
      call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
    elseif (index(trim(string),'FILE:')==1) then
      if (string(6:6)=='.' .or. string(6:6)=='/') then
        ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
        filename = trim( extractword(trim(string(6:80)), 1) )
      else
        ! Otherwise assume we should look for the file in INPUTDIR
        filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
      endif
      if (.not. file_exists(filename)) call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")

      varname = trim( extractword(trim(string(6:)), 2) )
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      if (len_trim(varname)==0) then
        if (field_exists(filename,'h_max')) then; varname = 'h_max'
        elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max'
        else ; call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
          "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
        endif
      endif
      call mom_read_data(trim(filename), trim(varname), h_max)
      call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
    elseif (index(trim(string),'FNC1:')==1) then
      call dz_function1( trim(string(6:)), h_max )
      call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
    else
      call mom_error(fatal,trim(mod)//", initialize_regridding: "// &
        "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
    endif
    deallocate(h_max)
  endif

  if (allocated(dz)) deallocate(dz)
end subroutine initialize_regridding

!> Do some basic checks on the vertical grid definition file, variable
subroutine check_grid_def(filename, varname, expected_units, msg, ierr)
  character(len=*), intent(in)    :: filename !< File name
  character(len=*), intent(in)    :: varname !< Variable name
  character(len=*), intent(in)    :: expected_units !< Expected units of variable
  character(len=*), intent(inout) :: msg !< Message to use for errors
  logical,          intent(out)   :: ierr !< True if an error occurs
  ! Local variables
  character (len=200) :: units, long_name
  integer :: ncid, status, intid, vid
  integer :: i

  ierr = .false.
  status = nf90_open(trim(filename), nf90_nowrite, ncid);
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'File not found: '//trim(filename)
    return
  endif

  status = nf90_inq_varid(ncid, trim(varname), vid)
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'Var not found: '//trim(varname)
    return
  endif

  status = nf90_get_att(ncid, vid, "units", units)
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'Attribute not found: units'
    return
  endif
  ! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.
  ! This loop replaces any null characters with a space so that the following check between
  ! the read units and the expected units will pass
  do i=1,len_trim(units)
    if (units(i:i) == char(0)) units(i:i) = " "
  enddo

  if (trim(units) /= trim(expected_units)) then
    if (trim(expected_units) == "meters") then
      if (trim(units) /= "m") then
        ierr = .true.
      endif
    else
      ierr = .true.
    endif
  endif

  if (ierr) then
    msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)
  endif

end subroutine check_grid_def

!> Deallocation of regridding memory
subroutine end_regridding(CS)
  type(regridding_cs), intent(inout) :: CS !< Regridding control structure

  if (associated(cs%zlike_CS))  call end_coord_zlike(cs%zlike_CS)
  if (associated(cs%sigma_CS))  call end_coord_sigma(cs%sigma_CS)
  if (associated(cs%rho_CS))    call end_coord_rho(cs%rho_CS)
  if (associated(cs%hycom_CS))  call end_coord_hycom(cs%hycom_CS)
  if (associated(cs%slight_CS)) call end_coord_slight(cs%slight_CS)
  if (associated(cs%adapt_CS))  call end_coord_adapt(cs%adapt_CS)

  deallocate( cs%coordinateResolution )
  if (allocated(cs%target_density)) deallocate( cs%target_density )
  if (allocated(cs%max_interface_depths) ) deallocate( cs%max_interface_depths )
  if (allocated(cs%max_layer_thickness) ) deallocate( cs%max_layer_thickness )

end subroutine end_regridding

!------------------------------------------------------------------------------
! Dispatching regridding routine: regridding & remapping
!------------------------------------------------------------------------------
subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
!------------------------------------------------------------------------------
! This routine takes care of (1) building a new grid and (2) remapping between
! the old grid and the new grid. The creation of the new grid can be based
! on z coordinates, target interface densities, sigma coordinates or any
! arbitrary coordinate system.
!   The MOM6 interface positions are always calculated from the bottom up by
! accumulating the layer thicknesses starting at z=-G%bathyT.  z increases
! upwards (decreasing k-index).
!   The new grid is defined by the change in position of those interfaces in z
!       dzInterface = zNew - zOld.
!   Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
!       hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
! IMPORTANT NOTE:
!   This is the converse of the sign convention used in the remapping code!
!------------------------------------------------------------------------------

  ! Arguments
  type(remapping_cs),                         intent(in)    :: remapCS !< Remapping parameters and options
  type(regridding_cs),                        intent(in)    :: CS     !< Regridding control structure
  type(ocean_grid_type),                      intent(in)    :: G      !< Ocean grid structure
  type(verticalgrid_type),                    intent(in)    :: GV     !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h      !< Current 3D grid obtained after the last time step
  type(thermo_var_ptrs),                      intent(inout) :: tv     !< Thermodynamical variables (T, S, ...)
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h_new  !< New 3D grid consistent with target coordinate
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in position of each interface
  real, dimension(:,:),                   optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
  logical,                          optional, intent(in   ) :: conv_adjust ! If true, do convective adjustment
  ! Local variables
  real :: trickGnuCompiler
  logical :: use_ice_shelf
  logical :: do_convective_adjustment

  do_convective_adjustment = .true.
  if (present(conv_adjust)) do_convective_adjustment = conv_adjust

  use_ice_shelf = .false.
  if (present(frac_shelf_h)) then
    if (associated(frac_shelf_h)) use_ice_shelf = .true.
  endif

  select case ( cs%regridding_scheme )

    case ( regridding_zstar )
      if (use_ice_shelf) then
        call build_zstar_grid( cs, g, gv, h, dzinterface, frac_shelf_h )
      else
        call build_zstar_grid( cs, g, gv, h, dzinterface )
      endif
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_sigma_shelf_zstar)
      call build_zstar_grid( cs, g, gv, h, dzinterface )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_sigma )
      call build_sigma_grid( cs, g, gv, h, dzinterface )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_rho )
      if (do_convective_adjustment) call convective_adjustment(g, gv, h, tv)
      call build_rho_grid( g, gv, h, tv, dzinterface, remapcs, cs )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_arbitrary )
      call build_grid_arbitrary( g, gv, h, dzinterface, trickgnucompiler, cs )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_hycom1 )
      call build_grid_hycom1( g, gv, h, tv, dzinterface, cs )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_slight )
      call build_grid_slight( g, gv, h, tv, dzinterface, cs )
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case ( regridding_adaptive )
      call build_grid_adaptive(g, gv, h, tv, dzinterface, remapcs, cs)
      call calc_h_new_by_dz(g, gv, h, dzinterface, h_new)

    case default
      call mom_error(fatal,'MOM_regridding, regridding_main: '//&
                     'Unknown regridding scheme selected!')

  end select ! type of grid

#ifdef __DO_SAFETY_CHECKS__
  call check_remapping_grid(g, gv, h, dzinterface,'in regridding_main')
#endif

end subroutine regridding_main

!> Calculates h_new from h + delta_k dzInterface
subroutine calc_h_new_by_dz(G, GV, h, dzInterface, h_new)
  type(ocean_grid_type),                      intent(in)    :: G !< Grid structure
  type(verticalgrid_type),                    intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h !< Old layer thicknesses (m)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: dzInterface !< Change in interface positions (m)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(inout) :: h_new !< New layer thicknesses (m)
  ! Local variables
  integer :: i, j, k

!$OMP parallel do default(none) shared(G,GV,h,dzInterface,h_new)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
      if (g%mask2dT(i,j)>0.) then
        do k=1,gv%ke
          h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
        enddo
      else
        h_new(i,j,:) = h(i,j,:)
      endif
    enddo
  enddo

end subroutine calc_h_new_by_dz

!> Check that the total thickness of two grids match
subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
  type(ocean_grid_type),                       intent(in) :: G !< Grid structure
  type(verticalgrid_type),                     intent(in) :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in) :: h !< Layer thicknesses (m)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions (m)
  character(len=*),                            intent(in) :: msg !< Message to append to errors
  ! Local variables
  integer :: i, j

!$OMP parallel do default(none) shared(G,GV,h,dzInterface,msg)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
      if (g%mask2dT(i,j)>0.) call check_grid_column( gv%ke, g%bathyT(i,j), h(i,j,:), dzinterface(i,j,:), msg )
    enddo
  enddo

end subroutine check_remapping_grid

!> Check that the total thickness of new and old grids are consistent
subroutine check_grid_column( nk, depth, h, dzInterface, msg )
  integer,               intent(in) :: nk !< Number of cells
  real,                  intent(in) :: depth !< Depth of bottom (m)
  real, dimension(nk),   intent(in) :: h !< Cell thicknesses (m)
  real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (m)
  character(len=*),      intent(in) :: msg !< Message to append to errors
  ! Local variables
  integer :: k
  real    :: eps, total_h_old, total_h_new, h_new, z_old, z_new

  eps =1. ; eps = epsilon(eps)

  ! Total thickness of grid h
  total_h_old = 0.
  do k = 1,nk
    total_h_old = total_h_old + h(k)
  enddo

  ! Integrate upwards for the interfaces consistent with the rest of MOM6
  z_old = - depth
  if (depth == 0.) z_old = - total_h_old
  total_h_new = 0.
  do k = nk,1,-1
    z_old = z_old + h(k) ! Old interface position above layer k
    z_new = z_old + dzinterface(k) ! New interface position based on dzInterface
    h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) ) ! New thickness
    if (h_new<0.) then
      write(0,*) 'k,h,hnew=',k,h(k),h_new
      write(0,*) 'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
      call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
        'Negative layer thickness implied by re-gridding, '//trim(msg))
    endif
    total_h_new = total_h_new + h_new

  enddo

  ! Conservation by implied h_new
  if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
    write(0,*) 'nk=',nk
    do k = 1,nk
      write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
    enddo
    write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
    write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
    call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
      'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
  endif

  ! Check that the top and bottom are intentionally moving
  if (dzinterface(1) /= 0.) call mom_error( fatal, &
    'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
  if (dzinterface(nk+1) /= 0.) call mom_error( fatal, &
    'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))

end subroutine check_grid_column

!> Returns the change in interface position motion after filtering and
!! assuming the top and bottom interfaces do not move.  The filtering is
!! a function of depth, and is applied as the integrated average filtering
!! over the trajectory of the interface.  By design, this code can not give
!! tangled interfaces provided that z_old and z_new are not already tangled.
subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
  type(regridding_cs),                           intent(in)    :: CS !< Regridding control structure
  integer,               intent(in)    :: nk !< Number of cells
  real, dimension(nk+1), intent(in)    :: z_old !< Old grid position (m)
  real, dimension(nk+1), intent(in)    :: z_new !< New grid position (m)
  real, dimension(nk+1), intent(inout) :: dz_g !< Change in interface positions (m)
  ! Local variables
  real :: sgn  ! The sign convention for downward.
  real :: dz_tgt, zr1
  real :: Aq, Bq, dz0, z0, F0
  real :: zs, zd, dzwt, Idzwt
  real :: wtd, Iwtd
  real :: Int_zs, Int_zd, dInt_zs_zd
! For debugging:
  real, dimension(nk+1) :: z_act
!  real, dimension(nk+1) :: ddz_g_s, ddz_g_d
  logical :: debug = .false.
  integer :: k

  if ((z_old(nk+1) - z_old(1)) * (z_new(nk+1) - z_new(1)) < 0.0) then
    call mom_error(fatal, "filtered_grid_motion: z_old and z_new use different sign conventions.")
  elseif ((z_old(nk+1) - z_old(1)) * (z_new(nk+1) - z_new(1)) == 0.0) then
    ! This is a massless column, so do nothing and return.
    do k=1,nk+1 ; dz_g(k) = 0.0 ; enddo ; return
  elseif ((z_old(nk+1) - z_old(1)) + (z_new(nk+1) - z_new(1)) > 0.0) then
    sgn = 1.0
  else
    sgn = -1.0
  endif

  if (debug) then
    do k=2,nk+1
      if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_new is tangled.")
      if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_old is tangled.")
    enddo
    ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
  endif

  zs = cs%depth_of_time_filter_shallow
  zd = cs%depth_of_time_filter_deep
  wtd = 1.0 - cs%old_grid_weight
  iwtd = 1.0 / wtd

  dzwt = (zd - zs)
  idzwt = 0.0 ; if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
  dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
  aq = 0.5*(iwtd - 1.0)

  dz_g(1) = 0.0
  do k = 2,nk
    ! zr1 is positive and increases with depth, and dz_tgt is positive downward.
    dz_tgt = sgn*(z_new(k) - z_old(k))
    zr1 = sgn*(z_old(k) - z_old(1))

    !   First, handle the two simple and common cases that do not pass through
    ! the adjustment rate transition zone.
    if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
      dz_g(k) = sgn * wtd * dz_tgt
    elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
      dz_g(k) = sgn * dz_tgt
    else
      ! Find the new value by inverting the equation
      !   integral(0 to dz_new) Iwt(z) dz = dz_tgt
      ! This is trivial where Iwt is a constant, and agrees with the two limits above.

      ! Take test values at the transition points to figure out which segment
      ! the new value will be found in.
      if (zr1 >= zd) then
        int_zd = iwtd*(zd - zr1)
        int_zs = int_zd - dint_zs_zd
      elseif (zr1 <= zs) then
        int_zs = (zs - zr1)
        int_zd = dint_zs_zd + (zs - zr1)
      else
!        Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs))
        int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
        int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
        ! It has been verified that  Int_zs = Int_zd - dInt_zs_zd to within roundoff.
      endif

      if (dz_tgt >= int_zd) then ! The new location is in the deep, slow region.
        dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
      elseif (dz_tgt <= int_zs) then ! The new location is in the shallow region.
        dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
      else  ! We need to solve a quadratic equation for z_new.
        ! For accuracy, do the integral from the starting depth or the nearest
        ! edge of the transition region.  The results with each choice are
        ! mathematically equivalent, but differ in roundoff, and this choice
        ! should minimize the likelihood of inadvertently overlapping interfaces.
        if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
        elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
        else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ; endif

        bq = (dzwt + 2.0*aq*(z0-zs))
        ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0
        ! Note that b>=0, and the two terms in the standard form cancel for the right root.
        dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))

!       if (debug) then
!         dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs))
!         ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
!         dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs))
!         ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
!
!         if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
!           call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
!         if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
!           call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
!       endif
      endif

    endif
  enddo
  dz_g(nk+1) = 0.0

  if (debug) then
    do k=1,nk+1 ; z_act(k) = z_old(k) + dz_g(k) ; enddo
    do k=2,nk+1
      if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_output is tangled.")
    enddo
  endif

end subroutine filtered_grid_motion

!> Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004).
!! z* is defined as
!!   z* = (z-eta)/(H+eta)*H  s.t. z*=0 when z=eta and z*=-H when z=-H .
subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)

  ! Arguments
  type(regridding_cs),                          intent(in)    :: CS !< Regridding control structure
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalgrid_type),                      intent(in)    :: GV !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in)    :: h  !< Layer thicknesses, in H
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H.
  real, dimension(:,:),               optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage.
  ! Local variables
  integer :: i, j, k
  integer :: nz
  real    :: nominalDepth, totalThickness, dh
  real, dimension(SZK_(GV)+1) :: zOld, zNew
  real :: minThickness
  logical :: ice_shelf

  nz = gv%ke
  minthickness = cs%min_thickness
  ice_shelf = .false.
  if (present(frac_shelf_h)) then
    if (associated(frac_shelf_h)) ice_shelf = .true.
  endif

!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, &
!$OMP                                  ice_shelf,minThickness) &
!$OMP                          private(nominalDepth,totalThickness, &
!$OMP                                  zNew,dh,zOld)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1

      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif

      ! Local depth (G%bathyT is positive)
      nominaldepth = g%bathyT(i,j)*gv%m_to_H

      ! Determine water column thickness
      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      end do

      zold(nz+1) = - nominaldepth
      do k = nz,1,-1
        zold(k) = zold(k+1) + h(i,j,k)
      enddo

      if (ice_shelf) then
        if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
           call build_zstar_column(cs%zlike_CS, nz, nominaldepth, totalthickness, znew, &
                                z_rigid_top = totalthickness-nominaldepth, &
                                eta_orig = zold(1))
        else
           call build_zstar_column(cs%zlike_CS, nz, nominaldepth, totalthickness, znew)
        endif
      else
        call build_zstar_column(cs%zlike_CS, nz, nominaldepth, totalthickness, znew)
      endif

      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )

#ifdef __DO_SAFETY_CHECKS__
      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',minthickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,nz
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
      endif
#endif

      call adjust_interface_motion( nz, cs%min_thickness, h(i,j,:), dzinterface(i,j,:) )

    end do
  end do

end subroutine build_zstar_grid

!------------------------------------------------------------------------------
! Build sigma grid
!------------------------------------------------------------------------------
subroutine build_sigma_grid( CS, G, GV, h, dzInterface )
!------------------------------------------------------------------------------
! This routine builds a grid based on terrain-following coordinates.
! The module parameter coordinateResolution(:) determines the resolution in
! sigma coordinate, dSigma(:). sigma-coordinates are defined by
!   sigma = (eta-z)/(H+eta)  s.t. sigma=0 at z=eta and sigma=1 at z=-H .
!------------------------------------------------------------------------------

  ! Arguments
  type(regridding_cs),                          intent(in)    :: CS !< Regridding control structure
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalgrid_type),                      intent(in)    :: GV !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in)    :: h  !< Layer thicknesses, in H
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H.

  ! Local variables
  integer :: i, j, k
  integer :: nz
  real    :: nominalDepth, totalThickness, dh
  real, dimension(SZK_(GV)+1) :: zOld, zNew

  nz = gv%ke

  do i = g%isc-1,g%iec+1
    do j = g%jsc-1,g%jec+1

      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif

      ! The rest of the model defines grids integrating up from the bottom
      nominaldepth = g%bathyT(i,j)*gv%m_to_H

      ! Determine water column height
      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      end do

      call build_sigma_column(cs%sigma_CS, nz, nominaldepth, totalthickness, znew)

      ! Calculate the final change in grid position after blending new and old grids
      zold(nz+1) =  -nominaldepth
      do k = nz,1,-1
        zold(k) = zold(k+1) + h(i, j, k)
      end do

      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )

#ifdef __DO_SAFETY_CHECKS__
      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',cs%min_thickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,nz
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
      endif
      dzinterface(i,j,1) = 0.
      dzinterface(i,j,nz+1) = 0.
#endif

    end do
  end do

end subroutine build_sigma_grid

!------------------------------------------------------------------------------
! Build grid based on target interface densities
!------------------------------------------------------------------------------
subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
!------------------------------------------------------------------------------
! This routine builds a new grid based on a given set of target interface
! densities (these target densities are computed by taking the mean value
! of given layer densities). The algorithn operates as follows within each
! column:
! 1. Given T & S within each layer, the layer densities are computed.
! 2. Based on these layer densities, a global density profile is reconstructed
!    (this profile is monotonically increasing and may be discontinuous)
! 3. The new grid interfaces are determined based on the target interface
!    densities.
! 4. T & S are remapped onto the new grid.
! 5. Return to step 1 until convergence or until the maximum number of
!    iterations is reached, whichever comes first.
!------------------------------------------------------------------------------

  ! Arguments
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalgrid_type),                      intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in)    :: h  !< Layer thicknesses, in H
  type(thermo_var_ptrs),                        intent(in)    :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H
  type(remapping_cs),                           intent(in)    :: remapCS !< The remapping control structure
  type(regridding_cs),                          intent(in)    :: CS !< Regridding control structure

  ! Local variables
  integer :: nz
  integer :: i, j, k
  real    :: nominalDepth, totalThickness
  real, dimension(SZK_(GV)+1) :: zOld, zNew
#ifdef __DO_SAFETY_CHECKS__
  real    :: dh
#endif

  nz = gv%ke

  if (.not.cs%target_density_set) call mom_error(fatal, "build_rho_grid: "//&
        "Target densities must be set before build_rho_grid is called.")

  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1

      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif


      ! Local depth (G%bathyT is positive)
      nominaldepth = g%bathyT(i,j)*gv%m_to_H

      call build_rho_column(cs%rho_CS, remapcs, nz, nominaldepth, h(i, j, :)*gv%H_to_m, &
                            tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew)

      if (cs%integrate_downward_for_e) then
        zold(1) = 0.
        do k = 1,nz
          zold(k+1) = zold(k) - h(i,j,k)
        enddo
      else
        ! The rest of the model defines grids integrating up from the bottom
        zold(nz+1) = - nominaldepth
        do k = nz,1,-1
          zold(k) = zold(k+1) + h(i,j,k)
        enddo
      endif

      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )

#ifdef __DO_SAFETY_CHECKS__
      do k = 2,nz
        if (znew(k) > zold(1)) then
          write(0,*) 'zOld=',zold
          write(0,*) 'zNew=',znew
          call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
               'interior interface above surface!' )
        endif
        if (znew(k) > znew(k-1)) then
          write(0,*) 'zOld=',zold
          write(0,*) 'zNew=',znew
          call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
               'interior interfaces cross!' )
        endif
      enddo

      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      enddo

      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',cs%min_thickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,nz
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_rho_grid: top surface has moved!!!' )
      endif
#endif

    end do  ! end loop on i
  end do  ! end loop on j

end subroutine build_rho_grid




!> Builds a simple HyCOM-like grid with the deepest location of potential
!! density interpolated from the column profile and a clipping of depth for
!! each interface to a fixed z* or p* grid.  This should probably be (optionally?)
!! changed to find the nearest location of the target density.
!! \remark { Based on Bleck, 2002: An oceanice general circulation model framed in
!! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88.
!! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 }
subroutine build_grid_hycom1( G, GV, h, tv, dzInterface, CS )
  type(ocean_grid_type),                       intent(in)    :: G  !< Grid structure
  type(verticalgrid_type),                     intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h  !< Existing model thickness, in H units
  type(thermo_var_ptrs),                       intent(in)    :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
  type(regridding_cs),                         intent(in)    :: CS !< Regridding control structure

  ! Local variables
  real, dimension(SZK_(GV)+1) :: z_col, z_col_new ! Interface positions relative to the surface in H units (m or kg m-2)
  real, dimension(SZK_(GV)+1) :: dz_col  ! The realized change in z_col in H units (m or kg m-2)
  real, dimension(SZK_(GV))   :: p_col   ! Layer pressure in Pa
  integer   :: i, j, k, nz
  real :: depth

  nz = gv%ke

  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_HyCOM1 : "//&
        "Target densities must be set before build_grid_HyCOM1 is called.")

  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j)>0.) then

      depth = g%bathyT(i,j) * gv%m_to_H

      z_col(1) = 0. ! Work downward rather than bottom up
      do k = 1, nz
        z_col(k+1) = z_col(k) + h(i,j,k) ! Work in units of h (m or Pa)
        p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
             ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
      enddo

      call build_hycom1_column(cs%hycom_CS, tv%eqn_of_state, nz, depth, &
                               h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new)

      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
      do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ; enddo

      ! This adjusts things robust to round-off errors
      call adjust_interface_motion( nz, cs%min_thickness, h(i,j,:), dzinterface(i,j,:) )

    else ! on land
      dzinterface(i,j,:) = 0.
    endif ! mask2dT
  enddo; enddo ! i,j

end subroutine build_grid_hycom1

subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS)
  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(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface
  type(remapping_cs),                          intent(in)    :: remapCS
  type(regridding_cs),                         intent(in)    :: CS

  ! local variables
  integer :: i, j, k, nz ! indices and dimension lengths
  ! temperature, salinity and pressure on interfaces
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
  ! current interface positions and after tendency term is applied
  ! positive downward
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
  real, dimension(SZK_(GV)+1) :: zNext

  nz = gv%ke

  ! position surface at z = 0.
  zint(:,:,1) = 0.

  ! work on interior interfaces
  do k = 2, nz ; do j = g%jsc-2,g%jec+2 ; do i = g%isc-2,g%iec+2
    tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
    sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
    zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1) ! zInt in [H]
  enddo ; enddo ; enddo

  ! top and bottom temp/salt interfaces are just the layer
  ! average values
  tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
  sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)

  ! set the bottom interface depth
  zint(:,:,nz+1)  = zint(:,:,nz) + h(:,:,nz)

  ! calculate horizontal density derivatives (alpha/beta)
  ! between cells in a 5-point stencil, columnwise
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j) < 0.5) then
      dzinterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
      cycle
    endif

    call build_adapt_column(cs%adapt_CS, g, gv, tv, i, j, zint, tint, sint, h, znext)

    call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
    ! convert from depth to z
    do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ; enddo
    call adjust_interface_motion(nz, cs%min_thickness, h(i,j,:), dzinterface(i,j,:))
  enddo ; enddo
end subroutine build_grid_adaptive

!> Builds a grid that tracks density interfaces for water that is denser than
!! the surface density plus an increment of some number of layers, and uses all
!! lighter layers uniformly above this location.  Note that this amounts to
!! interpolating to find the depth of an arbitrary (non-integer) interface index
!! which should make the results vary smoothly in space to the extent that the
!! surface density and interior stratification vary smoothly in space.  Over
!! shallow topography, this will tend to give a uniform sigma-like coordinate.
!! For sufficiently shallow water, a minimum grid spacing is used to avoid
!! certain instabilities.
subroutine build_grid_slight( G, GV, h, tv, dzInterface, CS )
  type(ocean_grid_type),                       intent(in)    :: G  !< Grid structure
  type(verticalgrid_type),                     intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h  !< Existing model thickness, in H units
  type(thermo_var_ptrs),                       intent(in)    :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
  type(regridding_cs),                         intent(in)    :: CS !< Regridding control structure

  real, dimension(SZK_(GV)+1) :: z_col, z_col_new ! Interface positions relative to the surface in H units (m or kg m-2)
  real, dimension(SZK_(GV)+1) :: dz_col  ! The realized change in z_col in H units (m or kg m-2)
  real, dimension(SZK_(GV))   :: p_col   ! Layer pressure in Pa
  real :: depth
  integer :: i, j, k, nz

  nz = gv%ke

  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_SLight : "//&
        "Target densities must be set before build_grid_SLight is called.")

  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j)>0.) then

      depth = g%bathyT(i,j) * gv%m_to_H
      z_col(1) = 0. ! Work downward rather than bottom up
      do k=1,nz
        z_col(k+1) = z_col(k) + h(i, j, k) ! Work in units of h (m or Pa)
        p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
                    ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
      enddo

      call build_slight_column(cs%slight_CS, tv%eqn_of_state, gv%H_to_Pa, gv%m_to_H, &
                          gv%H_subroundoff, nz, depth, &
                          h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new)

      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
      do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ; enddo
#ifdef __DO_SAFETY_CHECKS__
      if (dzinterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!'
      if (dzinterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!'
#endif

      ! This adjusts things robust to round-off errors
      call adjust_interface_motion( nz, cs%min_thickness, h(i,j,:), dzinterface(i,j,:) )

    else ! on land
      dzinterface(i,j,:) = 0.
    endif ! mask2dT
  enddo; enddo ! i,j

end subroutine build_grid_slight

!> Adjust dz_Interface to ensure non-negative future thicknesses
subroutine adjust_interface_motion( nk, min_thickness, h_old, dz_int )
  integer,               intent(in)    :: nk !< Number of layers
  real,                  intent(in)    :: min_thickness !< Minium allowed thickness of h (H units)
  real, dimension(nk),   intent(in)    :: h_old !< Minium allowed thickness of h (H units)
  real, dimension(nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h (H units)
  ! Local variables
  integer :: k
  real :: h_new, eps, h_total, h_err

  eps = 1. ; eps = epsilon(eps)

  h_total = 0. ; h_err = 0.
  do k = 1, nk
    h_total = h_total + h_old(k)
    h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new < -3.0*h_err) then
      write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
        'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
        'h_new=',h_new,'h_err=',h_err
      call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
                     'implied h<0 is larger than roundoff!')
    endif
  enddo
  do k = nk,2,-1
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<min_thickness) dz_int(k) = ( dz_int(k+1) - h_old(k) ) + min_thickness ! Implies next h_new = min_thickness
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<0.) dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<0.) then
      write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
        'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
        'h_new=',h_new
      stop 'Still did not work!'
      call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
                     'Repeated adjustment for roundoff h<0 failed!')
    endif
  enddo
 !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'

end subroutine adjust_interface_motion

!------------------------------------------------------------------------------
! Build arbitrary grid
!------------------------------------------------------------------------------
subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS )
!------------------------------------------------------------------------------
! This routine builds a grid based on arbitrary rules
!------------------------------------------------------------------------------

  ! Arguments
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalgrid_type),                      intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)),   intent(in)    :: h  !< Original ayer thicknesses, in H
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H
  real,                                         intent(inout) :: h_new !< New layer thicknesses, in H
  type(regridding_cs),                          intent(in)    :: CS !< Regridding control structure

  ! Local variables
  integer   :: i, j, k
  integer   :: nz
  real      :: z_inter(szk_(gv)+1)
  real      :: total_height
  real      :: delta_h
  real      :: max_depth
  real      :: min_thickness
  real      :: eta              ! local elevation
  real      :: local_depth
  real      :: x1, y1, x2, y2
  real      :: x, t

  nz = gv%ke

  max_depth = g%max_depth
  min_thickness = cs%min_thickness

  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1

      ! Local depth
      local_depth = g%bathyT(i,j)*gv%m_to_H

      ! Determine water column height
      total_height = 0.0
      do k = 1,nz
        total_height = total_height + h(i,j,k)
      end do

      eta = total_height - local_depth

      ! Compute new thicknesses based on stretched water column
      delta_h = (max_depth + eta) / nz

      ! Define interfaces
      z_inter(1) = eta
      do k = 1,nz
        z_inter(k+1) = z_inter(k) - delta_h
      end do

      ! Refine grid in the middle
      do k = 1,nz+1
        x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55

        x = - ( z_inter(k) - eta ) / max_depth

        if ( x <= x1 ) then
          t = y1*x/x1
        else if ( (x > x1 ) .and. ( x < x2 )) then
          t = y1 + (y2-y1) * (x-x1) / (x2-x1)
        else
          t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
        end if

        z_inter(k) = -t * max_depth + eta

      end do

      ! Modify interface heights to account for topography
      z_inter(nz+1) = - local_depth

      ! Modify interface heights to avoid layers of zero thicknesses
      do k = nz,1,-1
        if ( z_inter(k) < (z_inter(k+1) + min_thickness) ) then
          z_inter(k) = z_inter(k+1) + min_thickness
        end if
      end do

      ! Chnage in interface position
      x = 0. ! Left boundary at x=0
      dzinterface(i,j,1) = 0.
      do k = 2,nz
        x = x + h(i,j,k)
        dzinterface(i,j,k) = z_inter(k) - x
      end do
      dzinterface(i,j,nz+1) = 0.

    end do
  end do

stop 'OOOOOOPS' ! For some reason the gnu compiler will not let me delete this
                ! routine????

end subroutine build_grid_arbitrary



!------------------------------------------------------------------------------
! Check grid integrity
!------------------------------------------------------------------------------
subroutine inflate_vanished_layers_old( CS, G, GV, h )
!------------------------------------------------------------------------------
! This routine is called when initializing the regridding options. The
! objective is to make sure all layers are at least as thick as the minimum
! thickness allowed for regridding purposes (this parameter is set in the
! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
! are inflated up to the minmum thickness.
!------------------------------------------------------------------------------

  ! Arguments
  type(regridding_cs),                    intent(in)    :: CS
  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(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)

  ! Local variables
  integer :: i, j, k
  real    :: hTmp(gv%ke)

  do i = g%isc-1,g%iec+1
    do j = g%jsc-1,g%jec+1

      ! Build grid for current column
      do k = 1,gv%ke
        htmp(k) = h(i,j,k)
      end do

      call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )

      ! Save modified grid
      do k = 1,gv%ke
        h(i,j,k) = htmp(k)
      end do

    end do
  end do

end subroutine inflate_vanished_layers_old

!------------------------------------------------------------------------------
! Convective adjustment by swapping layers
!------------------------------------------------------------------------------
subroutine convective_adjustment(G, GV, h, tv)
!------------------------------------------------------------------------------
! Check each water column to see whether it is stratified. If not, sort the
! layers by successive swappings of water masses (bubble sort algorithm)
!------------------------------------------------------------------------------

  ! Arguments
  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(inout) :: 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

  ! Local variables
  integer   :: i, j, k
  real      :: T0, T1       ! temperatures
  real      :: S0, S1       ! salinities
  real      :: r0, r1       ! densities
  real      :: h0, h1
  logical   :: stratified
  real, dimension(GV%ke) :: p_col, densities

  p_col(:) = 0.

  ! Loop on columns
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

    ! Compute densities within current water column
    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, &
                            densities, 1, gv%ke, tv%eqn_of_state )

    ! Repeat restratification until complete
    do
      stratified = .true.
      do k = 1,gv%ke-1
        ! Gather information of current and next cells
        t0 = tv%T(i,j,k)  ; t1 = tv%T(i,j,k+1)
        s0 = tv%S(i,j,k)  ; s1 = tv%S(i,j,k+1)
        r0 = densities(k) ; r1 = densities(k+1)
        h0 = h(i,j,k) ; h1 = h(i,j,k+1)
        ! If the density of the current cell is larger than the density
        ! below it, we swap the cells and recalculate the densitiies
        ! within the swapped cells
        if ( r0 > r1 ) then
          tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
          tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
          h(i,j,k)    = h1 ; h(i,j,k+1)    = h0
          ! Recompute densities at levels k and k+1
          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
                                  densities(k), tv%eqn_of_state )
          call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
                                  densities(k+1), tv%eqn_of_state )
          stratified = .false.
        end if
      enddo  ! k

      if ( stratified ) exit
    enddo

  enddo ; enddo  ! i & j

end subroutine convective_adjustment


!------------------------------------------------------------------------------
! Return uniform resolution vector based on coordiante mode
!------------------------------------------------------------------------------
function uniformresolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy)
!------------------------------------------------------------------------------
! Calculate a vector of uniform resolution in the units of the coordinate
!------------------------------------------------------------------------------
  ! Arguments
  integer,          intent(in) :: nk
  character(len=*), intent(in) :: coordMode
  real,             intent(in) :: maxDepth, rhoLight, rhoHeavy
  real                         :: uniformResolution(nk)

  ! Local variables
  integer :: scheme

  scheme = coordinatemode(coordmode)
  select case ( scheme )

    case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_sigma_shelf_zstar, &
           regridding_adaptive )
      uniformresolution(:) = maxdepth / real(nk)

    case ( regridding_rho )
      uniformresolution(:) = (rhoheavy - rholight) / real(nk)

    case ( regridding_sigma )
      uniformresolution(:) = 1. / real(nk)

    case default
      call mom_error(fatal, "MOM_regridding, uniformResolution: "//&
       "Unrecognized choice for coordinate mode ("//trim(coordmode)//").")

  end select ! type of grid

end function uniformresolution

subroutine initcoord(CS, coord_mode)
  type(regridding_cs), intent(inout) :: CS
  character(len=*),    intent(in)    :: coord_mode

  select case (coordinatemode(coord_mode))
  case (regridding_zstar)
    call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
  case (regridding_sigma_shelf_zstar)
    call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
  case (regridding_sigma)
    call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
  case (regridding_rho)
    call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS)
  case (regridding_hycom1)
    call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, cs%interp_CS)
  case (regridding_slight)
    call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS)
  case (regridding_adaptive)
    call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution)
  end select
end subroutine initcoord

!------------------------------------------------------------------------------
! Set the fixed resolution data
!------------------------------------------------------------------------------
subroutine setcoordinateresolution( dz, CS )
  real, dimension(:),  intent(in)    :: dz
  type(regridding_cs), intent(inout) :: CS

  if (size(dz)/=cs%nk) call mom_error( fatal, &
      'setCoordinateResolution: inconsistent number of levels' )

  cs%coordinateResolution(:) = dz(:)

end subroutine setcoordinateresolution

!> Set target densities based on the old Rlay variable
subroutine set_target_densities_from_gv( GV, CS )
  type(verticalgrid_type), intent(in)    :: GV !< Ocean vertical grid structure
  type(regridding_cs),     intent(inout) :: CS !< Regridding control structure
  ! Local variables
  integer :: k, nz

  nz = cs%nk
  cs%target_density(1)    = gv%Rlay(1)+0.5*(gv%Rlay(1)-gv%Rlay(2))
  cs%target_density(nz+1) = gv%Rlay(nz)+0.5*(gv%Rlay(nz)-gv%Rlay(nz-1))
  do k = 2,nz
    cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
  end do
  cs%target_density_set = .true.

end subroutine set_target_densities_from_gv

!> Set target densities based on vector of interface values
subroutine set_target_densities( CS, rho_int )
  type(regridding_cs),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: rho_int !< Interface densities

  cs%target_density(:) = rho_int(:)
  cs%target_density_set = .true.

end subroutine set_target_densities

!> Set maximum interface depths based on a vector of input values.
subroutine set_regrid_max_depths( CS, max_depths, units_to_H )
  type(regridding_cs),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: max_depths !< Maximum interface depths, in arbitrary units
  real, optional,           intent(in)    :: units_to_H !< A conversion factor for max_depths into H units
  ! Local variables
  real :: val_to_H
  integer :: K

  if (.not.allocated(cs%max_interface_depths)) allocate(cs%max_interface_depths(1:cs%nk+1))

  val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h
  if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h

  ! Check for sign reversals in the depths.
  if (max_depths(cs%nk+1) < max_depths(1)) then
    do k=1,cs%nk ; if (max_depths(k+1) > max_depths(k)) &
      call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths!")
    enddo
  else
    do k=1,cs%nk ; if (max_depths(k+1) < max_depths(k)) &
      call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths.")
    enddo
  endif

  do k=1,cs%nk+1
    cs%max_interface_depths(k) = val_to_h * max_depths(k)
  enddo

  ! set max depths for coordinate
  select case (cs%regridding_scheme)
  case (regridding_hycom1)
    call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
  case (regridding_slight)
    call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
  end select
end subroutine set_regrid_max_depths

!> Set maximum layer thicknesses based on a vector of input values.
subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
  type(regridding_cs),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: max_h !< Maximum interface depths, in arbitrary units
  real, optional,           intent(in)    :: units_to_H !< A conversion factor for max_h into H units
  ! Local variables
  real :: val_to_H
  integer :: K

  if (.not.allocated(cs%max_layer_thickness)) allocate(cs%max_layer_thickness(1:cs%nk))

  val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h

  do k=1,cs%nk
    cs%max_layer_thickness(k) = val_to_h * max_h(k)
  enddo

  ! set max thickness for coordinate
  select case (cs%regridding_scheme)
  case (regridding_hycom1)
    call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
  case (regridding_slight)
    call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
  end select
end subroutine set_regrid_max_thickness


!------------------------------------------------------------------------------
! Query the fixed resolution data
!------------------------------------------------------------------------------
function getcoordinateresolution( CS )
  type(regridding_cs), intent(in) :: CS
  real, dimension(CS%nk)          :: getCoordinateResolution

  getcoordinateresolution(:) = cs%coordinateResolution(:)

end function getcoordinateresolution

!> Query the target coordinate interface positions
function getcoordinateinterfaces( CS )
  type(regridding_cs), intent(in) :: CS                      !< Regridding control structure
  real, dimension(CS%nk+1)        :: getCoordinateInterfaces !< Interface positions in target coordinate

  integer :: k

  ! When using a coordinate with target densities, we need to get the actual
  ! densities, rather than computing the interfaces based on resolution
  if (cs%regridding_scheme == regridding_rho) then
    if (.not. cs%target_density_set) &
      call mom_error(fatal, 'MOM_regridding, getCoordinateInterfaces: '//&
                            'target densities not set!')

    getcoordinateinterfaces(:) = cs%target_density(:)
  else
    getcoordinateinterfaces(1) = 0.
    do k = 1, cs%nk
      getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) &
                                    -cs%coordinateResolution(k)
    enddo
    ! The following line has an "abs()" to allow ferret users to reference
    ! data by index. It is a temporary work around...  :(  -AJA
    getcoordinateinterfaces(:) = abs( getcoordinateinterfaces(:) )
  end if

end function getcoordinateinterfaces

!------------------------------------------------------------------------------
! Query the target coordinate units
!------------------------------------------------------------------------------
function getcoordinateunits( CS )
  type(regridding_cs), intent(in) :: CS
  character(len=20)               :: getCoordinateUnits

  select case ( cs%regridding_scheme )
    case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
      getcoordinateunits = 'meter'
    case ( regridding_sigma_shelf_zstar )
      getcoordinateunits = 'meter/fraction'
    case ( regridding_sigma )
      getcoordinateunits = 'fraction'
    case ( regridding_rho )
      getcoordinateunits = 'kg/m3'
    case ( regridding_arbitrary )
      getcoordinateunits = 'unknown'
    case default
      call mom_error(fatal,'MOM_regridding, getCoordinateUnits: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid

end function getcoordinateunits

!------------------------------------------------------------------------------
! Query the short name of the coordinate
!------------------------------------------------------------------------------
function getcoordinateshortname( CS )
  type(regridding_cs), intent(in) :: CS
  character(len=20)               :: getCoordinateShortName

  select case ( cs%regridding_scheme )
    case ( regridding_zstar )
      !getCoordinateShortName = 'z*'
      ! The following line is a temporary work around...  :(  -AJA
      getcoordinateshortname = 'pseudo-depth, -z*'
    case ( regridding_sigma_shelf_zstar )
      getcoordinateshortname = 'pseudo-depth, -z*/sigma'
    case ( regridding_sigma )
      getcoordinateshortname = 'sigma'
    case ( regridding_rho )
      getcoordinateshortname = 'rho'
    case ( regridding_arbitrary )
      getcoordinateshortname = 'coordinate'
    case ( regridding_hycom1 )
      getcoordinateshortname = 'z-rho'
    case ( regridding_slight )
      getcoordinateshortname = 's-rho'
    case ( regridding_adaptive )
      getcoordinateshortname = 'adaptive'
    case default
      call mom_error(fatal,'MOM_regridding, getCoordinateShortName: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid

end function getcoordinateshortname

!> Can be used to set any of the parameters for MOM_regridding.
subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
             interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
             compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
             nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
             halocline_strat_tol, integrate_downward_for_e, &
             adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
  type(regridding_cs), intent(inout) :: CS !< Regridding control structure
  logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
  real,    optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the new grid (m)
  real,    optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid
  character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
  real,    optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic (H units)
  real,    optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic (H units)
  real,    optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density
  real,    optional, intent(in) :: dz_min_surface !< The fixed resolution in the topmost SLight_nkml_min layers (m)
  integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickess layers at the top of the model
  real,    optional, intent(in) :: Rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential density (m)
  real,    optional, intent(in) :: nlay_ML_to_interior !< Number of layers to offset the mixed layer density to find resolved stratification (nondim)
  logical, optional, intent(in) :: fix_haloclines !< Detect regions with much weaker stratification in the coordinate
  real,    optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for spuriously unstable water mass profiles (m)
  real,    optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic halocline region.
  logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward from the top.
  real, optional, intent(in) :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
  logical, optional, intent(in) :: adaptDoMin

  if (present(interp_scheme)) call set_interp_scheme(cs%interp_CS, interp_scheme)
  if (present(boundary_extrapolation)) call set_interp_extrap(cs%interp_CS, boundary_extrapolation)

  if (present(old_grid_weight)) then
    if (old_grid_weight<0. .or. old_grid_weight>1.) &
      call mom_error(fatal,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
    cs%old_grid_weight = old_grid_weight
  endif
  if (present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
  if (present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
  if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
    if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow) call mom_error(fatal,'MOM_regridding, '//&
                     'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
  endif

  if (present(min_thickness)) cs%min_thickness = min_thickness
  if (present(compress_fraction)) cs%compressibility_fraction = compress_fraction
  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e

  select case (cs%regridding_scheme)
  case (regridding_zstar)
    if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
  case (regridding_sigma_shelf_zstar)
    if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
  case (regridding_sigma)
    if (present(min_thickness)) call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
  case (regridding_rho)
    if (present(min_thickness)) call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
    if (present(integrate_downward_for_e)) call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
    if (associated(cs%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
  case (regridding_hycom1)
    if (associated(cs%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
  case (regridding_slight)
    if (present(min_thickness))       call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
    if (present(dz_min_surface))      call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
    if (present(nz_fixed_surface))    call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
    if (present(rho_ml_avg_depth))    call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
    if (present(nlay_ml_to_interior)) call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
    if (present(fix_haloclines))      call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
    if (present(halocline_filt_len))  call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
    if (present(halocline_strat_tol)) call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
    if (present(compress_fraction))   call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
    if (associated(cs%slight_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_slight_params(cs%slight_CS, interp_cs=cs%interp_CS)
  case (regridding_adaptive)
    if (present(adapttimeratio)) call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
    if (present(adaptzoom))      call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
    if (present(adaptzoomcoeff)) call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
    if (present(adaptbuoycoeff)) call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
    if (present(adaptalpha))     call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
    if (present(adaptdomin))     call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
  end select

end subroutine set_regrid_params

!> Returns the number of levels/layers in the regridding control structure
integer function get_regrid_size(CS)
  type(regridding_cs), intent(inout) :: CS !< Regridding control structure

  get_regrid_size = cs%nk

end function get_regrid_size

function get_zlike_cs(CS)
  type(regridding_cs), intent(in) :: CS
  type(zlike_cs) :: get_zlike_CS

  get_zlike_cs = cs%zlike_CS
end function get_zlike_cs

function get_sigma_cs(CS)
  type(regridding_cs), intent(in) :: CS
  type(sigma_cs) :: get_sigma_CS

  get_sigma_cs = cs%sigma_CS
end function get_sigma_cs

function get_rho_cs(CS)
  type(regridding_cs), intent(in) :: CS
  type(rho_cs) :: get_rho_CS

  get_rho_cs = cs%rho_CS
end function get_rho_cs

!------------------------------------------------------------------------------
! Return coordinate-derived thicknesses for fixed coordinate systems
!------------------------------------------------------------------------------
function getstaticthickness( CS, SSH, depth )
  type(regridding_cs), intent(in) :: CS
  real,                intent(in) :: SSH
  real,                intent(in) :: depth
  real, dimension(CS%nk)          :: getStaticThickness
  ! Local
  integer :: k
  real :: z, dz

  select case ( cs%regridding_scheme )
    case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
      if (depth>0.) then
        z = ssh
        do k = 1, cs%nk
          dz = cs%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
          dz = max(dz, 0.)                                     ! Avoid negative incase ssh=-depth
          dz = min(dz, depth - z)                              ! Clip if below topography
          z = z + dz                                           ! Bottom of layer
          getstaticthickness(k) = dz
        enddo
      else
        getstaticthickness(:) = 0. ! On land ...
      endif
    case ( regridding_sigma )
      getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
    case ( regridding_rho )
      getstaticthickness(:) = 0. ! Not applicable
    case ( regridding_arbitrary )
      getstaticthickness(:) = 0.  ! Not applicable
    case default
      call mom_error(fatal,'MOM_regridding, getStaticThickness: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid

end function getstaticthickness

!> Parses a string and generates a dz(:) profile that goes like k**power.
subroutine dz_function1( string, dz )
  character(len=*),   intent(in)    :: string !< String with list of parameters in form
                                              !! dz_min, H_total, power, precision
  real, dimension(:), intent(inout) :: dz     !< Profile of nominal thicknesses
  ! Local variables
  integer :: nk, k
  real    :: dz_min, power, prec, H_total

  nk = size(dz) ! Number of cells
  prec = -1024.
  read( string, *) dz_min, h_total, power, prec
  if (prec == -1024.) call mom_error(fatal,"dz_function1: "// &
          "Problem reading FNC1: string  ="//trim(string))
  ! Create profile of ( dz - dz_min )
  do k = 1, nk
    dz(k) = (real(k-1)/real(nk-1))**power
  enddo
  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min

end subroutine dz_function1

!> \namespace mom_regridding
!!
!! A vertical grid is defined solely by the cell thicknesses, \f$h\f$.
!! Most calculations in this module start with the coordinate at the bottom
!! of the column set to -depth, and use a increasing value of coordinate with
!! decreasing k. This is consistent with the rest of MOM6 that uses position,
!! \f$z\f$ which is a negative quantity for most of the ocean.
!!
!! A change in grid is define through a change in position of the interfaces:
!! \f[
!! z^n_{k+1/2} = z^{n-1}_{k+1/2} + \Delta z_{k+1/2}
!! \f]
!! with the positive upward coordinate convention
!! \f[
!! z_{k-1/2} = z_{k+1/2} + h_k
!! \f]
!! so that
!! \f[
!! h^n_k = h^{n-1}_k + ( \Delta z_{k-1/2} - \Delta z_{k+1/2} )
!! \f]
!!
!! Original date of creation: 2008.06.09 by L. White

end module mom_regridding