MOM_ALE.F90

!> This module contains the main regridding routines.
!!
!! Regridding comprises two steps:
!! 1. Interpolation and creation of a new grid based on target interface
!!    densities (or any other criterion).
!! 2. Remapping of quantities between old grid and new grid.
!!
!! Original module written by Laurent White, 2008.06.09
module mom_ale

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

use mom_debugging,        only : check_column_integrals
use mom_diag_mediator,    only : register_diag_field, post_data, diag_ctrl
use mom_diag_mediator,    only : time_type, diag_update_remap_grids
use mom_diag_vkernels,    only : interpolate_column, reintegrate_column
use mom_domains,          only : create_group_pass, do_group_pass, group_pass_type
use mom_eos,              only : calculate_density
use mom_domains,          only : create_group_pass, do_group_pass, group_pass_type
use mom_error_handler,    only : mom_error, fatal, warning
use mom_error_handler,    only : calltree_showquery
use mom_error_handler,    only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,      only : get_param, param_file_type, log_param
use mom_io,               only : vardesc, var_desc, fieldtype, single_file
use mom_io,               only : create_file, write_field, close_file
use mom_interface_heights,only : find_eta
use mom_open_boundary,    only : ocean_obc_type, obc_direction_e, obc_direction_w
use mom_open_boundary,    only : obc_direction_n, obc_direction_s
use mom_regridding,       only : initialize_regridding, regridding_main, end_regridding
use mom_regridding,       only : uniformresolution
use mom_regridding,       only : inflate_vanished_layers_old
use mom_regridding,       only : set_target_densities_from_gv, set_target_densities
use mom_regridding,       only : regriddingcoordinatemodedoc, default_coordinate_mode
use mom_regridding,       only : regriddinginterpschemedoc, regriddingdefaultinterpscheme
use mom_regridding,       only : regriddingdefaultboundaryextrapolation
use mom_regridding,       only : regriddingdefaultminthickness
use mom_regridding,       only : regridding_cs, set_regrid_params
use mom_regridding,       only : getcoordinateinterfaces, getcoordinateresolution
use mom_regridding,       only : getcoordinateunits, getcoordinateshortname
use mom_regridding,       only : getstaticthickness
use mom_remapping,        only : initialize_remapping, end_remapping
use mom_remapping,        only : remapping_core_h, remapping_core_w
use mom_remapping,        only : remappingschemesdoc, remappingdefaultscheme
use mom_remapping,        only : remapping_cs, dzfromh1h2
use mom_string_functions, only : uppercase, extractword, extract_integer
use mom_tracer_registry,  only : tracer_registry_type, tracer_type, mom_tracer_chkinv
use mom_unit_scaling,     only : unit_scale_type
use mom_variables,        only : ocean_grid_type, thermo_var_ptrs
use mom_verticalgrid,     only : get_thickness_units, verticalgrid_type

!use regrid_consts,       only : coordinateMode, DEFAULT_COORDINATE_MODE
use regrid_consts,        only : coordinateunits, coordinatemode, state_dependent
use regrid_edge_values,   only : edge_values_implicit_h4
use plm_functions,        only : plm_reconstruction, plm_boundary_extrapolation
use plm_functions,        only : plm_extrapolate_slope, plm_monotonized_slope, plm_slope_wa
use ppm_functions,        only : ppm_reconstruction, ppm_boundary_extrapolation

implicit none ; private
#include <MOM_memory.h>


!> ALE control structure
type, public :: ale_cs ; private
  logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
                                    !! method. If False, uses the new method that
                                    !! remaps between grids described by h.

  real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid
                            !! and the target (new) grid [T ~> s]

  type(regridding_cs) :: regridcs !< Regridding parameters and work arrays
  type(remapping_cs)  :: remapcs  !< Remapping parameters and work arrays

  integer :: nk             !< Used only for queries, not directly by this module

  logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.

  logical :: answers_2018   !< If true, use the order of arithmetic and expressions for remapping
                            !! that recover the answers from the end of 2018.  Otherwise, use more
                            !! robust and accurate forms of mathematically equivalent expressions.

  logical :: show_call_tree !< For debugging

  ! for diagnostics
  type(diag_ctrl), pointer           :: diag                          !< structure to regulate output
  integer, dimension(:), allocatable :: id_tracer_remap_tendency      !< diagnostic id
  integer, dimension(:), allocatable :: id_htracer_remap_tendency     !< diagnostic id
  integer, dimension(:), allocatable :: id_htracer_remap_tendency_2d  !< diagnostic id
  logical, dimension(:), allocatable :: do_tendency_diag              !< flag for doing diagnostics
  integer                            :: id_dzregrid = -1              !< diagnostic id

  ! diagnostic for fields prior to applying ALE remapping
  integer :: id_u_preale = -1 !< diagnostic id for zonal velocity before ALE.
  integer :: id_v_preale = -1 !< diagnostic id for meridional velocity before ALE.
  integer :: id_h_preale = -1 !< diagnostic id for layer thicknesses before ALE.
  integer :: id_t_preale = -1 !< diagnostic id for temperatures before ALE.
  integer :: id_s_preale = -1 !< diagnostic id for salinities before ALE.
  integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE.
  integer :: id_vert_remap_h = -1      !< diagnostic id for layer thicknesses used for remapping
  integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE

end type

! Publicly available functions
public ale_init
public ale_end
public ale_main
public ale_main_offline
public ale_offline_inputs
public ale_offline_tracer_final
public ale_build_grid
public ale_regrid_accelerated
public ale_remap_scalar
public ale_plm_edge_values
public ts_plm_edge_values
public ts_ppm_edge_values
public adjustgridforintegrity
public ale_initregridding
public ale_getcoordinate
public ale_getcoordinateunits
public ale_writecoordinatefile
public ale_updateverticalgridtype
public ale_initthicknesstocoord
public ale_update_regrid_weights
public ale_remap_init_conds
public ale_register_diags

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

contains

!> This routine is typically called (from initialize_MOM in file MOM.F90)
!! before the main time integration loop to initialize the regridding stuff.
!! We read the MOM_input file to register the values of different
!! regridding/remapping parameters.
subroutine ale_init( param_file, GV, US, max_depth, CS)
  type(param_file_type),   intent(in) :: param_file !< Parameter file
  type(verticalgrid_type), intent(in) :: gv         !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: us         !< A dimensional unit scaling type
  real,                    intent(in) :: max_depth  !< The maximum depth of the ocean [Z ~> m].
  type(ale_cs),            pointer    :: cs         !< Module control structure

  ! Local variables
  real, dimension(:), allocatable :: dz
  character(len=40)               :: mdl = "MOM_ALE" ! This module's name.
  character(len=80)               :: string ! Temporary strings
  real                            :: filter_shallow_depth, filter_deep_depth
  logical       :: default_2018_answers
  logical                         :: check_reconstruction
  logical                         :: check_remapping
  logical                         :: force_bounds_in_subcell
  logical                         :: local_logical
  logical                         :: remap_boundary_extrap

  if (associated(cs)) then
    call mom_error(warning, "ALE_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  cs%show_call_tree = calltree_showquery()
  if (cs%show_call_tree) call calltree_enter("ALE_init(), MOM_ALE.F90")

  call get_param(param_file, mdl, "REMAP_UV_USING_OLD_ALG", cs%remap_uv_using_old_alg, &
                 "If true, uses the old remapping-via-a-delta-z method for "//&
                 "remapping u and v. If false, uses the new method that remaps "//&
                 "between grids described by an old and new thickness.", &
                 default=.false.)

  ! Initialize and configure regridding
  call ale_initregridding(gv, us, max_depth, param_file, mdl, cs%regridCS)

  ! Initialize and configure remapping
  call get_param(param_file, mdl, "REMAPPING_SCHEME", string, &
                 "This sets the reconstruction scheme used "//&
                 "for vertical remapping for all variables. "//&
                 "It can be one of the following schemes: "//&
                 trim(remappingschemesdoc), default=remappingdefaultscheme)
  call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
                 "If true, cell-by-cell reconstructions are checked for "//&
                 "consistency and if non-monotonicity or an inconsistency is "//&
                 "detected then a FATAL error is issued.", default=.false.)
  call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &
                 "If true, the results of remapping are checked for "//&
                 "conservation and new extrema and if an inconsistency is "//&
                 "detected then a FATAL error is issued.", default=.false.)
  call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
                 "If true, the values on the intermediate grid used for remapping "//&
                 "are forced to be bounded, which might not be the case due to "//&
                 "round off.", default=.false.)
  call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
                 "If true, values at the interfaces of boundary cells are "//&
                 "extrapolated instead of piecewise constant", default=.false.)
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
  call initialize_remapping( cs%remapCS, string, &
                             boundary_extrapolation=remap_boundary_extrap, &
                             check_reconstruction=check_reconstruction, &
                             check_remapping=check_remapping, &
                             force_bounds_in_subcell=force_bounds_in_subcell, &
                             answers_2018=cs%answers_2018)

  call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &
                 "If true, applies regridding and remapping immediately after "//&
                 "initialization so that the state is ALE consistent. This is a "//&
                 "legacy step and should not be needed if the initialization is "//&
                 "consistent with the coordinate mode.", default=.true.)

  call get_param(param_file, mdl, "REGRID_TIME_SCALE", cs%regrid_time_scale, &
                 "The time-scale used in blending between the current (old) grid "//&
                 "and the target (new) grid. A short time-scale favors the target "//&
                 "grid (0. or anything less than DT_THERM) has no memory of the old "//&
                 "grid. A very long time-scale makes the model more Lagrangian.", &
                 units="s", default=0., scale=us%s_to_T)
  call get_param(param_file, mdl, "REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &
                 "The depth above which no time-filtering is applied. Above this depth "//&
                 "final grid exactly matches the target (new) grid.", &
                 units="m", default=0., scale=gv%m_to_H)
  call get_param(param_file, mdl, "REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &
                 "The depth below which full time-filtering is applied with time-scale "//&
                 "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&
                 "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
                 units="m", default=0., scale=gv%m_to_H)
  call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
                         depth_of_time_filter_deep=filter_deep_depth)
  call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
                 "If true, the regridding ntegrates upwards from the bottom for "//&
                 "interface positions, much as the main model does. If false "//&
                 "regridding integrates downward, consistant with the remapping "//&
                 "code.", default=.true., do_not_log=.true.)
  call set_regrid_params(cs%regridCS, integrate_downward_for_e=.not.local_logical)

  ! Keep a record of values for subsequent queries
  cs%nk = gv%ke

  if (cs%show_call_tree) call calltree_leave("ALE_init()")
end subroutine ale_init

!> Initialize diagnostics for the ALE module.
subroutine ale_register_diags(Time, G, GV, US, diag, CS)
  type(time_type),target,     intent(in)  :: time  !< Time structure
  type(ocean_grid_type),      intent(in)  :: g     !< Grid structure
  type(unit_scale_type),      intent(in)  :: us    !< A dimensional unit scaling type
  type(verticalgrid_type),    intent(in)  :: gv    !< Ocean vertical grid structure
  type(diag_ctrl), target,    intent(in)  :: diag  !< Diagnostics control structure
  type(ale_cs), pointer                   :: cs    !< Module control structure

  cs%diag => diag

  ! These diagnostics of the state variables before ALE are useful for
  ! debugging the ALE code.
  cs%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, time, &
      'Zonal velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, time, &
      'Meridional velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, time, &
      'Layer Thickness before remapping', get_thickness_units(gv), &
      conversion=gv%H_to_MKS, v_extensive=.true.)
  cs%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, time, &
      'Temperature before remapping', 'degC')
  cs%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, time, &
      'Salinity before remapping', 'PSU')
  cs%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, time, &
      'Interface Heights before remapping', 'm', conversion=us%Z_to_m)

  cs%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,time, &
      'Change in interface height due to ALE regridding', 'm', &
      conversion=gv%H_to_m)
  cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', &
      diag%axestl, time, 'layer thicknesses after ALE regridding and remapping', 'm', &
      conversion=gv%H_to_m, v_extensive=.true.)
  cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, &
      'Layer thicknesses tendency due to ALE regridding and remapping', 'm', &
      conversion=gv%H_to_m*us%s_to_T, v_extensive = .true.)

end subroutine ale_register_diags

!> Crudely adjust (initial) grid for integrity.
!! This routine is typically called (from initialize_MOM in file MOM.F90)
!! before the main time integration loop to initialize the regridding stuff.
!! We read the MOM_input file to register the values of different
!! regridding/remapping parameters.
subroutine adjustgridforintegrity( CS, G, GV, h )
  type(ale_cs),                              pointer       :: cs  !< Regridding parameters and options
  type(ocean_grid_type),                     intent(in)    :: g   !< Ocean grid informations
  type(verticalgrid_type),                   intent(in)    :: gv  !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h   !< Current 3D grid thickness that
                                                                  !! are to be adjusted [H ~> m or kg-2]
  call inflate_vanished_layers_old( cs%regridCS, g, gv, h(:,:,:) )

end subroutine adjustgridforintegrity


!> End of regridding (memory deallocation).
!! This routine is typically called (from MOM_end in file MOM.F90)
!! after the main time integration loop to deallocate the regridding stuff.
subroutine ale_end(CS)
  type(ale_cs), pointer :: cs  !< module control structure

  ! Deallocate memory used for the regridding
  call end_remapping( cs%remapCS )
  call end_regridding( cs%regridCS )

  deallocate(cs)

end subroutine ale_end

!> Takes care of (1) building a new grid and (2) remapping all variables 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.
subroutine ale_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
  type(ocean_grid_type),                      intent(in)    :: g   !< Ocean grid informations
  type(verticalgrid_type),                    intent(in)    :: gv  !< Ocean vertical grid structure
  type(unit_scale_type),                      intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: h   !< Current 3D grid obtained after the
                                                                   !! last time step [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u   !< Zonal velocity field [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v   !< Meridional velocity field [L T-1 ~> m s-1]
  type(thermo_var_ptrs),                      intent(inout) :: tv  !< Thermodynamic variable structure
  type(tracer_registry_type),                 pointer       :: reg !< Tracer registry structure
  type(ale_cs),                               pointer       :: cs  !< Regridding parameters and options
  type(ocean_obc_type),                       pointer       :: obc !< Open boundary structure
  real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main [T ~> s]
  real, dimension(:,:),             optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg-2]
  integer :: nk, i, j, k, isc, iec, jsc, jec
  logical :: ice_shelf

  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec

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

  if (cs%show_call_tree) call calltree_enter("ALE_main(), MOM_ALE.F90")

  ! These diagnostics of the state before ALE is applied are mostly used for debugging.
  if (cs%id_u_preale > 0) call post_data(cs%id_u_preale, u,    cs%diag)
  if (cs%id_v_preale > 0) call post_data(cs%id_v_preale, v,    cs%diag)
  if (cs%id_h_preale > 0) call post_data(cs%id_h_preale, h,    cs%diag)
  if (cs%id_T_preale > 0) call post_data(cs%id_T_preale, tv%T, cs%diag)
  if (cs%id_S_preale > 0) call post_data(cs%id_S_preale, tv%S, cs%diag)
  if (cs%id_e_preale > 0) then
    call find_eta(h, tv, g, gv, us, eta_preale)
    call post_data(cs%id_e_preale, eta_preale, cs%diag)
  endif

  if (present(dt)) then
    call ale_update_regrid_weights( dt, cs )
  endif
  dzregrid(:,:,:) = 0.0

  ! Build new grid. The new grid is stored in h_new. The old grid is h.
  ! Both are needed for the subsequent remapping of variables.
  if (ice_shelf) then
     call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, frac_shelf_h)
  else
     call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid)
  endif

  call check_grid( g, gv, h, 0. )

  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_main)")

  ! The presence of dt is used for expediency to distinguish whether ALE_main is being called during init
  ! or in the main loop. Tendency diagnostics in remap_all_state_vars also rely on this logic.
  if (present(dt)) then
    call diag_update_remap_grids(cs%diag)
  endif
  ! Remap all variables from old grid h onto new grid h_new
  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, -dzregrid, &
                             u, v, cs%show_call_tree, dt )

  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_main)")

  ! Override old grid with new one. The new grid 'h_new' is built in
  ! one of the 'build_...' routines above.
  !$OMP parallel do default(shared)
  do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
    h(i,j,k) = h_new(i,j,k)
  enddo ; enddo ; enddo

  if (cs%show_call_tree) call calltree_leave("ALE_main()")

  if (cs%id_dzRegrid>0 .and. present(dt)) call post_data(cs%id_dzRegrid, dzregrid, cs%diag)


end subroutine ale_main

!> Takes care of (1) building a new grid and (2) remapping all variables 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.
subroutine ale_main_offline( G, GV, h, tv, Reg, CS, OBC, dt)
  type(ocean_grid_type),                      intent(in)    :: g   !< Ocean grid informations
  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 [H ~> m or kg-2]
  type(thermo_var_ptrs),                      intent(inout) :: tv  !< Thermodynamic variable structure
  type(tracer_registry_type),                 pointer       :: reg !< Tracer registry structure
  type(ale_cs),                               pointer       :: cs  !< Regridding parameters and options
  type(ocean_obc_type),                       pointer       :: obc !< Open boundary structure
  real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main [T ~> s]
  ! Local variables
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg-2]
  integer :: nk, i, j, k, isc, iec, jsc, jec

  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec

  if (cs%show_call_tree) call calltree_enter("ALE_main_offline(), MOM_ALE.F90")

  if (present(dt)) then
    call ale_update_regrid_weights( dt, cs )
  endif
  dzregrid(:,:,:) = 0.0

  ! Build new grid. The new grid is stored in h_new. The old grid is h.
  ! Both are needed for the subsequent remapping of variables.
  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid )

  call check_grid( g, gv, h, 0. )

  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_main)")

  ! Remap all variables from old grid h onto new grid h_new

  call remap_all_state_vars(cs%remapCS, cs, g, gv, h, h_new, reg, obc, &
                            debug=cs%show_call_tree, dt=dt )

  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_main)")

  ! Override old grid with new one. The new grid 'h_new' is built in
  ! one of the 'build_...' routines above.
  !$OMP parallel do default(shared)
  do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
    h(i,j,k) = h_new(i,j,k)
  enddo ; enddo ; enddo

  if (cs%show_call_tree) call calltree_leave("ALE_main()")
  if (cs%id_dzRegrid>0 .and. present(dt)) call post_data(cs%id_dzRegrid, dzregrid, cs%diag)

end subroutine ale_main_offline

!> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have
!! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This
!! routine builds a grid on the runtime specified vertical coordinate
subroutine ale_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
  type(ale_cs),                                 pointer       :: cs    !< Regridding parameters and options
  type(ocean_grid_type),                        intent(in   ) :: g     !< Ocean grid informations
  type(verticalgrid_type),                      intent(in   ) :: gv    !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(inout) :: h     !< Layer thicknesses
  type(thermo_var_ptrs),                        intent(inout) :: tv    !< Thermodynamic variable structure
  type(tracer_registry_type),                   pointer       :: reg   !< Tracer registry structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),   intent(inout) :: uhtr  !< Zonal mass fluxes
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),   intent(inout) :: vhtr  !< Meridional mass fluxes
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),  intent(inout) :: kd    !< Input diffusivites
  logical,                                      intent(in   ) :: debug !< If true, then turn checksums
  type(ocean_obc_type),                         pointer       :: obc   !< Open boundary structure
  ! Local variables
  integer :: nk, i, j, k, isc, iec, jsc, jec
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV))   :: h_new    ! Layer thicknesses after regridding
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
  real, dimension(SZK_(GV)) :: h_src
  real, dimension(SZK_(GV)) :: h_dest, uh_dest
  real, dimension(SZK_(GV)) :: temp_vec

  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
  dzregrid(:,:,:) = 0.0
  h_new(:,:,:) = 0.0

  if (debug) call mom_tracer_chkinv("Before ALE_offline_inputs", g, h, reg%Tr, reg%ntr)

  ! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored
  ! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective
  ! adjustment right now is not used because it is unclear what to do with vanished layers
  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, conv_adjust = .false. )
  call check_grid( g, gv, h_new, 0. )
  if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_offline_inputs)")

  ! Remap all variables from old grid h onto new grid h_new
  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, debug=cs%show_call_tree )
  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_inputs)")

  ! Reintegrate mass transports from Zstar to the offline vertical coordinate
  do j=jsc,jec ; do i=g%iscB,g%iecB
    if (g%mask2dCu(i,j)>0.) then
      h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
      h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
      call reintegrate_column(nk, h_src, uhtr(i,j,:), nk, h_dest, 0., temp_vec)
      uhtr(i,j,:) = temp_vec
    endif
  enddo ; enddo
  do j=g%jscB,g%jecB ; do i=isc,iec
    if (g%mask2dCv(i,j)>0.) then
      h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
      h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
      call reintegrate_column(nk, h_src, vhtr(i,j,:), nk, h_dest, 0., temp_vec)
      vhtr(i,j,:) = temp_vec
    endif
  enddo ; enddo

  do j = jsc,jec ; do i=isc,iec
    if (g%mask2dT(i,j)>0.) then
      if (check_column_integrals(nk, h_src, nk, h_dest)) then
        call mom_error(fatal, "ALE_offline_inputs: Kd interpolation columns do not match")
      endif
      call interpolate_column(nk, h(i,j,:), kd(i,j,:), nk, h_new(i,j,:), 0., kd(i,j,:))
    endif
  enddo ; enddo

  call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%T, h_new, tv%T, answers_2018=cs%answers_2018)
  call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%S, h_new, tv%S, answers_2018=cs%answers_2018)

  if (debug) call mom_tracer_chkinv("After ALE_offline_inputs", g, h_new, reg%Tr, reg%ntr)

  ! Copy over the new layer thicknesses
  do k = 1,nk  ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
      h(i,j,k) = h_new(i,j,k)
  enddo ; enddo ; enddo

  if (cs%show_call_tree) call calltree_leave("ALE_offline_inputs()")
end subroutine ale_offline_inputs


!> Remaps all tracers from h onto h_target. This is intended to be called when tracers
!! are done offline. In the case where transports don't quite conserve, we still want to
!! make sure that layer thicknesses offline do not drift too far away from the online model
subroutine ale_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS, OBC)
  type(ocean_grid_type),                      intent(in)    :: g   !< Ocean grid informations
  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 [H ~> m or kg-2]
  type(thermo_var_ptrs),                      intent(inout) :: tv  !< Thermodynamic variable structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: h_target !< Current 3D grid obtained after
                                                                        !! last time step  [H ~> m or kg-2]
  type(tracer_registry_type),                 pointer       :: reg !< Tracer registry structure
  type(ale_cs),                               pointer       :: cs  !< Regridding parameters and options
  type(ocean_obc_type),                       pointer       :: obc !< Open boundary structure
  ! Local variables

  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid !< The change in grid interface positions
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV))   :: h_new    !< Regridded target thicknesses
  integer :: nk, i, j, k, isc, iec, jsc, jec

  nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec

  if (cs%show_call_tree) call calltree_enter("ALE_offline_tracer_final(), MOM_ALE.F90")
  ! Need to make sure that h_target is consistent with the current offline ALE confiuration
  call regridding_main( cs%remapCS, cs%regridCS, g, gv, h_target, tv, h_new, dzregrid )
  call check_grid( g, gv, h_target, 0. )


  if (cs%show_call_tree) call calltree_waypoint("Source and target grids checked (ALE_offline_tracer_final)")

  ! Remap all variables from old grid h onto new grid h_new

  call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, obc, debug=cs%show_call_tree )

  if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_offline_tracer_final)")

  ! Override old grid with new one. The new grid 'h_new' is built in
  ! one of the 'build_...' routines above.
  !$OMP parallel do default(shared)
  do k = 1,nk
    do j = jsc-1,jec+1 ; do i = isc-1,iec+1
      h(i,j,k) = h_new(i,j,k)
    enddo ; enddo
  enddo
  if (cs%show_call_tree) call calltree_leave("ALE_offline_tracer_final()")
end subroutine ale_offline_tracer_final

!> Check grid for negative thicknesses
subroutine check_grid( G, GV, h, threshold )
  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  !< Current 3D grid obtained after the
                                                              !! last time step [H ~> m or kg m-2]
  real,                                      intent(in) :: threshold !< Value below which to flag issues,
                                                              !! [H ~> m or kg m-2]
  ! Local variables
  integer :: i, j

  do j = g%jsc,g%jec ; do i = g%isc,g%iec
    if (g%mask2dT(i,j)>0.) then
      if (minval(h(i,j,:)) < threshold) then
        write(0,*) 'check_grid: i,j=',i,j,'h(i,j,:)=',h(i,j,:)
        if (threshold <= 0.) then
          call mom_error(fatal,"MOM_ALE, check_grid: negative thickness encountered.")
        else
          call mom_error(fatal,"MOM_ALE, check_grid: too tiny thickness encountered.")
        endif
      endif
    endif
  enddo ; enddo


end subroutine check_grid

!> Generates new grid
subroutine ale_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
  type(ocean_grid_type),                   intent(in)    :: g        !< Ocean grid structure
  type(verticalgrid_type),                 intent(in)    :: gv       !< Ocean vertical grid structure
  type(regridding_cs),                     intent(in)    :: regridcs !< Regridding parameters and options
  type(remapping_cs),                      intent(in)    :: remapcs  !< Remapping parameters and options
  type(thermo_var_ptrs),                   intent(inout) :: tv       !< Thermodynamical variable structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h     !< Current 3D grid obtained after the
                                                                     !! last time step [H ~> m or kg-2]
  logical,                       optional, intent(in)    :: debug    !< If true, show the call tree
  real, dimension(:,:),          optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage
  ! Local variables
  integer :: nk, i, j, k
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! The new grid thicknesses
  logical :: show_call_tree, use_ice_shelf

  show_call_tree = .false.
  if (present(debug)) show_call_tree = debug
  if (show_call_tree) call calltree_enter("ALE_build_grid(), MOM_ALE.F90")
  use_ice_shelf = .false.
  if (present(frac_shelf_h)) then
    if (associated(frac_shelf_h)) use_ice_shelf = .true.
  endif

  ! Build new grid. The new grid is stored in h_new. The old grid is h.
  ! Both are needed for the subsequent remapping of variables.
  if (use_ice_shelf) then
     call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid, frac_shelf_h )
  else
     call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid )
  endif

  ! Override old grid with new one. The new grid 'h_new' is built in
  ! one of the 'build_...' routines above.
!$OMP parallel do default(none) shared(G,h,h_new)
  do j = g%jsc,g%jec ; do i = g%isc,g%iec
    if (g%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:)
  enddo ; enddo

  if (show_call_tree) call calltree_leave("ALE_build_grid()")
end subroutine ale_build_grid

!> For a state-based coordinate, accelerate the process of regridding by
!! repeatedly applying the grid calculation algorithm
subroutine ale_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial)
  type(ale_cs),            pointer       :: cs     !< ALE control structure
  type(ocean_grid_type),   intent(inout) :: g      !< Ocean grid
  type(verticalgrid_type), intent(in)    :: gv     !< Vertical grid
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: h      !< Original thicknesses [H ~> m or kg-2]
  type(thermo_var_ptrs),   intent(inout) :: tv     !< Thermo vars (T/S/EOS)
  integer,                 intent(in)    :: n      !< Number of times to regrid
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: u      !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                           intent(inout) :: v      !< Meridional velocity [L T-1 ~> m s-1]
  type(ocean_obc_type),    pointer       :: obc    !< Open boundary structure
  type(tracer_registry_type), &
                 optional, pointer       :: reg    !< Tracer registry to remap onto new grid
  real,          optional, intent(in)    :: dt     !< Model timestep to provide a timescale for regridding [T ~> s]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                 optional, intent(inout) :: dzregrid !< Final change in interface positions
  logical,       optional, intent(in)    :: initial !< Whether we're being called from an initialization
                                                    !! routine (and expect diagnostics to work)

  ! Local variables
  integer :: i, j, k, nz
  type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt
  type(group_pass_type) :: pass_t_s_h ! group pass if the coordinate has a stencil
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV))         :: h_loc, h_orig ! A working copy of layer thicknesses
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: t, s ! local temporary state
  ! we have to keep track of the total dzInterface if for some reason
  ! we're using the old remapping algorithm for u/v
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinterface, dzinttotal

  nz = gv%ke

  ! initial total interface displacement due to successive regridding
  dzinttotal(:,:,:) = 0.

  call create_group_pass(pass_t_s_h, t, g%domain)
  call create_group_pass(pass_t_s_h, s, g%domain)
  call create_group_pass(pass_t_s_h, h_loc, g%domain)

  ! copy original temp/salt and set our local tv_pointers to them
  tv_local = tv
  t(:,:,:) = tv%T(:,:,:)
  s(:,:,:) = tv%S(:,:,:)
  tv_local%T => t
  tv_local%S => s

  ! get local copy of thickness and save original state for remapping
  h_loc(:,:,:) = h(:,:,:)
  h_orig(:,:,:) = h(:,:,:)

  ! Apply timescale to regridding (for e.g. filtered_grid_motion)
  if (present(dt)) &
    call ale_update_regrid_weights(dt, cs)

  do k = 1, n
    call do_group_pass(pass_t_s_h, g%domain)

    ! generate new grid
    call regridding_main(cs%remapCS, cs%regridCS, g, gv, h_loc, tv_local, h, dzinterface)
    dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)

    ! remap from original grid onto new grid
    do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
      call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
      call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
    enddo ; enddo

    ! starting grid for next iteration
    h_loc(:,:,:) = h(:,:,:)
  enddo

  ! remap all state variables (including those that weren't needed for regridding)
  call remap_all_state_vars(cs%remapCS, cs, g, gv, h_orig, h, reg, obc, dzinttotal, u, v)

  ! save total dzregrid for diags if needed?
  if (present(dzregrid)) dzregrid(:,:,:) = dzinttotal(:,:,:)
end subroutine ale_regrid_accelerated

!> This routine takes care of remapping all variable between the old and the
!! new grids. When velocity components need to be remapped, thicknesses at
!! velocity points are taken to be arithmetic averages of tracer thicknesses.
!! This routine is called during initialization of the model at time=0, to
!! remap initiali conditions to the model grid.  It is also called during a
!! time step to update the state.
subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, OBC, &
                                dxInterface, u, v, debug, dt)
  type(remapping_CS),                        intent(in)    :: CS_remapping !< Remapping control structure
  type(ALE_CS),                              intent(in)    :: CS_ALE       !< ALE 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_old        !< Thickness of source grid
                                                                           !! [H ~> m or kg-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_new        !< Thickness of destination grid
                                                                           !! [H ~> m or kg-2]
  type(tracer_registry_type),                pointer       :: Reg          !< Tracer registry structure
  type(ocean_OBC_type),                      pointer       :: OBC          !< Open boundary structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                                   optional, intent(in)    :: dxInterface  !< Change in interface position
                                                                           !! [H ~> m or kg-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                                   optional, intent(inout) :: u      !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                                   optional, intent(inout) :: v      !< Meridional velocity [L T-1 ~> m s-1]
  logical,                         optional, intent(in)    :: debug  !< If true, show the call tree
  real,                            optional, intent(in)    :: dt     !< time step for diagnostics [T ~> s]
  ! Local variables
  integer                                     :: i, j, k, m
  integer                                     :: nz, ntr
  real, dimension(GV%ke+1)                    :: dx
  real, dimension(GV%ke)                      :: h1, u_column
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc
  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont
  real, dimension(SZI_(G), SZJ_(G))           :: work_2d
  real                                        :: Idt ! The inverse of the timestep [T-1 ~> s-1]
  real                                        :: ppt2mks
  real, dimension(GV%ke)                      :: h2
  real :: h_neglect, h_neglect_edge
  logical                                     :: show_call_tree
  type(tracer_type), pointer                  :: Tr => null()

  show_call_tree = .false.
  if (present(debug)) show_call_tree = debug
  if (show_call_tree) call calltree_enter("remap_all_state_vars(), MOM_ALE.F90")

  ! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dxInterface. Otherwise,
  ! u and v can be remapped without dxInterface
  if ( .not. present(dxinterface) .and. (cs_ale%remap_uv_using_old_alg .and. (present(u) .or. present(v))) ) then
    call mom_error(fatal, "remap_all_state_vars: dxInterface must be present if using old algorithm "// &
                          "and u/v are to be remapped")
  endif

  if (.not.cs_ale%answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif

  nz      = gv%ke
  ppt2mks = 0.001

  ntr = 0 ; if (associated(reg)) ntr = reg%ntr

  if (present(dt)) then
    idt = 1.0/dt
    work_conc(:,:,:) = 0.0
    work_cont(:,:,:) = 0.0
  endif

  ! Remap tracer
  if (ntr>0) then
    if (show_call_tree) call calltree_waypoint("remapping tracers (remap_all_state_vars)")
    !$OMP parallel do default(shared) private(h1,h2,u_column,Tr)
    do m=1,ntr ! For each tracer
      tr => reg%Tr(m)
      do j = g%jsc,g%jec ; do i = g%isc,g%iec ; if (g%mask2dT(i,j)>0.) then
        ! Build the start and final grids
        h1(:) = h_old(i,j,:)
        h2(:) = h_new(i,j,:)
        call remapping_core_h(cs_remapping, nz, h1, tr%t(i,j,:), nz, h2, &
                              u_column, h_neglect, h_neglect_edge)

        ! Intermediate steps for tendency of tracer concentration and tracer content.
        if (present(dt)) then
          if (tr%id_remap_conc > 0) then
            do k=1,gv%ke
              work_conc(i,j,k) = (u_column(k) - tr%t(i,j,k)) * idt
            enddo
          endif
          if (tr%id_remap_cont > 0 .or. tr%id_remap_cont_2d > 0) then
            do k=1,gv%ke
              work_cont(i,j,k) = (u_column(k)*h2(k) - tr%t(i,j,k)*h1(k)) * idt
            enddo
          endif
        endif
        ! update tracer concentration
        tr%t(i,j,:) = u_column(:)
      endif ; enddo ; enddo

      ! tendency diagnostics.
      if (present(dt)) then
        if (tr%id_remap_conc > 0) then
          call post_data(tr%id_remap_conc, work_conc, cs_ale%diag)
        endif
        if (tr%id_remap_cont > 0) then
          call post_data(tr%id_remap_cont, work_cont, cs_ale%diag)
        endif
        if (tr%id_remap_cont_2d > 0) then
          do j = g%jsc,g%jec ; do i = g%isc,g%iec
            work_2d(i,j) = 0.0
            do k = 1,gv%ke
              work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)
            enddo
          enddo ; enddo
          call post_data(tr%id_remap_cont_2d, work_2d, cs_ale%diag)
        endif
      endif
    enddo ! m=1,ntr

  endif   ! endif for ntr > 0

  if (show_call_tree) call calltree_waypoint("tracers remapped (remap_all_state_vars)")

  ! Remap u velocity component
  if ( present(u) ) then
    !$OMP parallel do default(shared) private(h1,h2,dx,u_column)
    do j = g%jsc,g%jec ; do i = g%iscB,g%iecB ; if (g%mask2dCu(i,j)>0.) then
      ! Build the start and final grids
      h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i+1,j,:) )
      if (cs_ale%remap_uv_using_old_alg) then
        dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i+1,j,:) )
        do k = 1, nz
          h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
        enddo
      else
        h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
      endif
      if (associated(obc)) then
        if (obc%segnum_u(i,j) .ne. 0) then
          if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
            h1(:) = h_old(i,j,:)
            h2(:) = h_new(i,j,:)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
            h1(:) = h_old(i+1,j,:)
            h2(:) = h_new(i+1,j,:)
          endif
        endif
      endif
      call remapping_core_h(cs_remapping, nz, h1, u(i,j,:), nz, h2, &
                            u_column, h_neglect, h_neglect_edge)
      u(i,j,:) = u_column(:)
    endif ; enddo ; enddo
  endif

  if (show_call_tree) call calltree_waypoint("u remapped (remap_all_state_vars)")

  ! Remap v velocity component
  if ( present(v) ) then
    !$OMP parallel do default(shared) private(h1,h2,dx,u_column)
    do j = g%jscB,g%jecB ; do i = g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then
      ! Build the start and final grids
      h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i,j+1,:) )
      if (cs_ale%remap_uv_using_old_alg) then
        dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i,j+1,:) )
        do k = 1, nz
          h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
        enddo
      else
        h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
      endif
      if (associated(obc)) then
        if (obc%segnum_v(i,j) .ne. 0) then
          if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
            h1(:) = h_old(i,j,:)
            h2(:) = h_new(i,j,:)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
            h1(:) = h_old(i,j+1,:)
            h2(:) = h_new(i,j+1,:)
          endif
        endif
      endif
      call remapping_core_h(cs_remapping, nz, h1, v(i,j,:), nz, h2, &
                            u_column, h_neglect, h_neglect_edge)
      v(i,j,:) = u_column(:)
    endif ; enddo ; enddo
  endif

  if (cs_ale%id_vert_remap_h > 0) call post_data(cs_ale%id_vert_remap_h, h_old, cs_ale%diag)
  if ((cs_ale%id_vert_remap_h_tendency > 0) .and. present(dt)) then
    do k = 1, nz ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
      work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs_ale%id_vert_remap_h_tendency, work_cont, cs_ale%diag)
  endif
  if (show_call_tree) call calltree_waypoint("v remapped (remap_all_state_vars)")
  if (show_call_tree) call calltree_leave("remap_all_state_vars()")

end subroutine remap_all_state_vars


!> Remaps a single scalar between grids described by thicknesses h_src and h_dst.
!! h_dst must be dimensioned as a model array with GV%ke layers while h_src can
!! have an arbitrary number of layers specified by nk_src.
subroutine ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap, answers_2018 )
  type(remapping_cs),                      intent(in)    :: cs        !< Remapping control structure
  type(ocean_grid_type),                   intent(in)    :: g         !< Ocean grid structure
  type(verticalgrid_type),                 intent(in)    :: gv        !< Ocean vertical grid structure
  integer,                                 intent(in)    :: nk_src    !< Number of levels on source grid
  real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in)    :: h_src     !< Level thickness of source grid
                                                                      !! [H ~> m or kg-2]
  real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in)    :: s_src     !< Scalar on source grid
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in)   :: h_dst     !< Level thickness of destination grid
                                                                      !! [H ~> m or kg-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst    !< Scalar on destination grid
  logical, optional,                       intent(in)    :: all_cells !< If false, only reconstruct for
                                                                      !! non-vanished cells. Use all vanished
                                                                      !! layers otherwise (default).
  logical, optional,                       intent(in)    :: old_remap !< If true, use the old "remapping_core_w"
                                                                      !! method, otherwise use "remapping_core_h".
  logical,                       optional, intent(in)    :: answers_2018 !< If true, use the order of arithmetic
                                                                      !! and expressions that recover the answers for
                                                                      !! remapping from the end of 2018. Otherwise,
                                                                      !! use more robust forms of the same expressions.
  ! Local variables
  integer :: i, j, k, n_points
  real :: dx(gv%ke+1)
  real :: h_neglect, h_neglect_edge
  logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap

  ignore_vanished_layers = .false.
  if (present(all_cells)) ignore_vanished_layers = .not. all_cells
  use_remapping_core_w = .false.
  if (present(old_remap)) use_remapping_core_w = old_remap
  n_points = nk_src
  use_2018_remap = .true. ; if (present(answers_2018)) use_2018_remap = answers_2018

  if (.not.use_2018_remap) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif

  !$OMP parallel do default(shared) firstprivate(n_points,dx)
  do j = g%jsc,g%jec ; do i = g%isc,g%iec
    if (g%mask2dT(i,j) > 0.) then
      if (ignore_vanished_layers) then
        n_points = 0
        do k = 1, nk_src
          if (h_src(i,j,k)>0.) n_points = n_points + 1
        enddo
        s_dst(i,j,:) = 0.
      endif
      if (use_remapping_core_w) then
        call dzfromh1h2( n_points, h_src(i,j,1:n_points), gv%ke, h_dst(i,j,:), dx )
        call remapping_core_w(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
                              gv%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge)
      else
        call remapping_core_h(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
                              gv%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge)
      endif
    else
      s_dst(i,j,:) = 0.
    endif
  enddo ; enddo

end subroutine ale_remap_scalar


!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction
!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
subroutine ts_plm_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
  type(ocean_grid_type),   intent(in)    :: g    !< ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< Ocean vertical grid structure
  type(ale_cs),            intent(inout) :: cs   !< module control structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: s_t  !< Salinity at the top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: s_b  !< Salinity at the bottom edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: t_t  !< Temperature at the top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: t_b  !< Temperature at the bottom edge of each layer
  type(thermo_var_ptrs),   intent(in)    :: tv   !< thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: h    !< layer thickness [H ~> m or kg m-2]
  logical,                 intent(in)    :: bdry_extrap !< If true, use high-order boundary
                                                 !! extrapolation within boundary cells

  call ale_plm_edge_values( cs, g, gv, h, tv%S, bdry_extrap, s_t, s_b )
  call ale_plm_edge_values( cs, g, gv, h, tv%T, bdry_extrap, t_t, t_b )

end subroutine ts_plm_edge_values

!> Calculate edge values (top and bottom of layer) 3d scalar array.
!! Boundary reconstructions are PCM unless bdry_extrap is true.
subroutine ale_plm_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b )
  type(ale_cs),            intent(in)    :: cs   !< module 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 thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: q    !< 3d scalar array
  logical,                 intent(in)    :: bdry_extrap !< If true, use high-order boundary
                                                 !! extrapolation within boundary cells
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: q_t  !< Scalar at the top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: q_b  !< Scalar at the bottom edge of each layer
  ! Local variables
  integer :: i, j, k
  real :: slp(gv%ke)
  real :: mslp
  real :: h_neglect

  if (.not.cs%answers_2018) then
    h_neglect = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30
  endif

  !$OMP parallel do default(shared) private(slp,mslp)
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    slp(1) = 0.
    do k = 2, gv%ke-1
      slp(k) = plm_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, q(i,j,k-1), q(i,j,k), q(i,j,k+1))
    enddo
    slp(gv%ke) = 0.

    do k = 2, gv%ke-1
      mslp = plm_monotonized_slope(q(i,j,k-1), q(i,j,k), q(i,j,k+1), slp(k-1), slp(k), slp(k+1))
      q_t(i,j,k) = q(i,j,k) - 0.5 * mslp
      q_b(i,j,k) = q(i,j,k) + 0.5 * mslp
    enddo
    if (bdry_extrap) then
      mslp = - plm_extrapolate_slope(h(i,j,2), h(i,j,1), h_neglect, q(i,j,2), q(i,j,1))
      q_t(i,j,1) = q(i,j,1) - 0.5 * mslp
      q_b(i,j,1) = q(i,j,1) + 0.5 * mslp
      mslp = plm_extrapolate_slope(h(i,j,gv%ke-1), h(i,j,gv%ke), h_neglect, q(i,j,gv%ke-1), q(i,j,gv%ke))
      q_t(i,j,gv%ke) = q(i,j,gv%ke) - 0.5 * mslp
      q_b(i,j,gv%ke) = q(i,j,gv%ke) + 0.5 * mslp
    else
      q_t(i,j,1) = q(i,j,1)
      q_b(i,j,1) = q(i,j,1)
      q_t(i,j,gv%ke) = q(i,j,gv%ke)
      q_b(i,j,gv%ke) = q(i,j,gv%ke)
    endif

  enddo ; enddo

end subroutine ale_plm_edge_values

!> Calculate edge values (top and bottom of layer) for T and S consistent with a PPM reconstruction
!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
subroutine ts_ppm_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
  type(ocean_grid_type),   intent(in)    :: g    !< ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< Ocean vertical grid structure
  type(ale_cs),            intent(inout) :: cs   !< module control structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: s_t  !< Salinity at the top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: s_b  !< Salinity at the bottom edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: t_t  !< Temperature at the top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: t_b  !< Temperature at the bottom edge of each layer
  type(thermo_var_ptrs),   intent(in)    :: tv   !< thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: h    !< layer thicknesses [H ~> m or kg m-2]
  logical,                 intent(in)    :: bdry_extrap !< If true, use high-order boundary
                                                 !! extrapolation within boundary cells

  ! Local variables
  integer :: i, j, k
  real    :: htmp(gv%ke) ! A 1-d copy of h [H ~> m or kg m-2]
  real    :: tmp(gv%ke)  ! A 1-d copy of a column of temperature [degC] or salinity [ppt]
  real, dimension(CS%nk,2) :: &
      ppol_e            ! Edge value of polynomial in [degC] or [ppt]
  real, dimension(CS%nk,3) :: &
      ppol_coefs        ! Coefficients of polynomial, all in [degC] or [ppt]
  real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]

  if (.not.cs%answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif

  ! Determine reconstruction within each column
  !$OMP parallel do default(shared) private(hTmp,tmp,ppol_E,ppol_coefs)
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

    ! Build current grid
    htmp(:) = h(i,j,:)
    tmp(:) = tv%S(i,j,:)

    ! Reconstruct salinity profile
    ppol_e(:,:) = 0.0
    ppol_coefs(:,:) = 0.0
    call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=h_neglect_edge, &
                                  answers_2018=cs%answers_2018 )
    call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect, &
                                  answers_2018=cs%answers_2018 )
    if (bdry_extrap) &
      call ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )

    do k = 1,gv%ke
      s_t(i,j,k) = ppol_e(k,1)
      s_b(i,j,k) = ppol_e(k,2)
    enddo

    ! Reconstruct temperature profile
    ppol_e(:,:) = 0.0
    ppol_coefs(:,:) = 0.0
    tmp(:) = tv%T(i,j,:)
    if (cs%answers_2018) then
      call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=1.0e-10*gv%m_to_H, &
                                  answers_2018=cs%answers_2018 )
    else
      call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=gv%H_subroundoff, &
                                  answers_2018=cs%answers_2018 )
    endif
    call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect, &
                                  answers_2018=cs%answers_2018 )
    if (bdry_extrap) &
      call ppm_boundary_extrapolation(gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )

    do k = 1,gv%ke
      t_t(i,j,k) = ppol_e(k,1)
      t_b(i,j,k) = ppol_e(k,2)
    enddo

  enddo ; enddo

end subroutine ts_ppm_edge_values


!> Initializes regridding for the main ALE algorithm
subroutine ale_initregridding(GV, US, max_depth, param_file, mdl, regridCS)
  type(verticalgrid_type), intent(in)  :: gv         !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)  :: us         !< A dimensional unit scaling type
  real,                    intent(in)  :: max_depth  !< The maximum depth of the ocean [Z ~> m].
  type(param_file_type),   intent(in)  :: param_file !< parameter file
  character(len=*),        intent(in)  :: mdl        !< Name of calling module
  type(regridding_cs),     intent(out) :: regridcs   !< Regridding parameters and work arrays
  ! Local variables
  character(len=30) :: coord_mode

  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", coord_mode, &
                 "Coordinate mode for vertical regridding. "//&
                 "Choose among the following possibilities: "//&
                 trim(regriddingcoordinatemodedoc), &
                 default=default_coordinate_mode, fail_if_missing=.true.)

  call initialize_regridding(regridcs, gv, us, max_depth, param_file, mdl, coord_mode, '', '')

end subroutine ale_initregridding

!> Query the target coordinate interfaces positions
function ale_getcoordinate( CS )
  type(ale_cs), pointer    :: cs                  !< module control structure

  real, dimension(CS%nk+1) :: ale_getcoordinate
  ale_getcoordinate(:) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )

end function ale_getcoordinate


!> Query the target coordinate units
function ale_getcoordinateunits( CS )
  type(ale_cs), pointer :: cs   !< module control structure

  character(len=20)     :: ale_getcoordinateunits

  ale_getcoordinateunits = getcoordinateunits( cs%regridCS )

end function ale_getcoordinateunits


!> Returns true if initial conditions should be regridded and remapped
logical function ale_remap_init_conds( CS )
  type(ale_cs), pointer :: cs   !< module control structure

  ale_remap_init_conds = .false.
  if (associated(cs)) ale_remap_init_conds = cs%remap_after_initialization
end function ale_remap_init_conds

!> Updates the weights for time filtering the new grid generated in regridding
subroutine ale_update_regrid_weights( dt, CS )
  real,         intent(in) :: dt !< Time-step used between ALE calls [T ~> s]
  type(ale_cs), pointer    :: cs !< ALE control structure
  ! Local variables
  real :: w  ! An implicit weighting estimate.

  if (associated(cs)) then
    w = 0.0
    if (cs%regrid_time_scale > 0.0) then
      w = cs%regrid_time_scale / (cs%regrid_time_scale + dt)
    endif
    call set_regrid_params(cs%regridCS, old_grid_weight=w)
  endif

end subroutine ale_update_regrid_weights

!> Update the vertical grid type with ALE information.
!! This subroutine sets information in the verticalGrid_type to be
!! consistent with the use of ALE mode.
subroutine ale_updateverticalgridtype(CS, GV)
  type(ale_cs),            pointer :: cs  !< ALE control structure
  type(verticalgrid_type), pointer :: gv  !< vertical grid information

  integer :: nk

  nk = gv%ke
  gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
  gv%sLayer(1:nk) = 0.5*( gv%sInterface(1:nk) + gv%sInterface(2:nk+1) )
  gv%zAxisUnits = getcoordinateunits( cs%regridCS )
  gv%zAxisLongName = getcoordinateshortname( cs%regridCS )
  gv%direction = -1 ! Because of ferret in z* mode. Need method to set
                    ! as function of coordinae mode.

end subroutine ale_updateverticalgridtype


!> Write the vertical coordinate information into a file.
!! This subroutine writes out a file containing any available data related
!! to the vertical grid used by the MOM ocean model when in ALE mode.
subroutine ale_writecoordinatefile( CS, GV, directory )
  type(ale_cs),            pointer     :: cs         !< module control structure
  type(verticalgrid_type), intent(in)  :: gv         !< ocean vertical grid structure
  character(len=*),        intent(in)  :: directory  !< directory for writing grid info

  character(len=240) :: filepath
  type(vardesc)      :: vars(2)
  type(fieldtype)    :: fields(2)
  integer            :: unit
  real               :: ds(gv%ke), dsi(gv%ke+1)

  filepath    = trim(directory) // trim("Vertical_coordinate")
  ds(:)       = getcoordinateresolution( cs%regridCS, undo_scaling=.true. )
  dsi(1)      = 0.5*ds(1)
  dsi(2:gv%ke) = 0.5*( ds(1:gv%ke-1) + ds(2:gv%ke) )
  dsi(gv%ke+1) = 0.5*ds(gv%ke)

  vars(1) = var_desc('ds', getcoordinateunits( cs%regridCS ), &
                    'Layer Coordinate Thickness','1','L','1')
  vars(2) = var_desc('ds_interface', getcoordinateunits( cs%regridCS ), &
                    'Layer Center Coordinate Separation','1','i','1')

  call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
  call write_field(unit, fields(1), ds)
  call write_field(unit, fields(2), dsi)
  call close_file(unit)

end subroutine ale_writecoordinatefile


!> Set h to coordinate values for fixed coordinate systems
subroutine ale_initthicknesstocoord( CS, G, GV, h )
  type(ale_cs), intent(inout)                            :: cs  !< module control structure
  type(ocean_grid_type), intent(in)                      :: g   !< module grid structure
  type(verticalgrid_type), intent(in)                    :: gv  !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h   !< layer thickness [H ~> m or kg m-2]

  ! Local variables
  integer :: i, j, k

  do j = g%jsd,g%jed ; do i = g%isd,g%ied
    h(i,j,:) = gv%Z_to_H * getstaticthickness( cs%regridCS, 0., g%bathyT(i,j) )
  enddo ; enddo

end subroutine ale_initthicknesstocoord

end module mom_ale