MOM_ALE.F90

Go to the documentation of this file.

!> 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, time_type
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_regridding,       only : initialize_regridding, regridding_main, end_regridding
use mom_regridding,       only : uniformresolution
use mom_regridding,       only : inflate_vanished_layers_old, setcoordinateresolution
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 : set_regrid_max_depths
use mom_regridding,       only : set_regrid_max_thickness
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, mom_tracer_chkinv
use mom_variables,        only : ocean_grid_type, thermo_var_ptrs
use mom_verticalgrid,     only : verticalgrid_type

use regrid_defs,          only : pressure_reconstruction_plm
!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 ppm_functions,        only : ppm_reconstruction, ppm_boundary_extrapolation
use p1m_functions,        only : p1m_interpolation,  p1m_boundary_extrapolation
use p3m_functions,        only : p3m_interpolation,  p3m_boundary_extrapolation


implicit none ; private
#include <MOM_memory.h>


!> ALE control structure
type, public :: ale_cs
  private

  logical :: boundary_extrapolation_for_pressure !< Indicate whether high-order boundary
                                                 !! extrapolation should be used within boundary cells

  logical :: reconstructforpressure = .false.    !< Indicates whether integrals for FV
                                                 !! pressure gradient calculation will
                                                 !! use reconstruction of T/S.
                                                 !! By default, it is true if regridding
                                                 !! has been initialized, otherwise false.

  integer :: pressurereconstructionscheme        !< Form of the reconstruction of T/S
                                                 !! for FV pressure gradient calculation.
                                                 !! By default, it is =1 (PLM)

  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. (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
  integer :: degree_linear=1 !< Degree of linear piecewise polynomial
  integer :: degree_parab=2  !< Degree of parabolic piecewise polynomial

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

  logical :: show_call_tree !< For debugging
  real    :: c_p            !< seawater heat capacity (J/(kg deg C))

  ! 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

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 pressure_gradient_plm
public pressure_gradient_ppm
public usepressurereconstruction
public pressurereconstructionscheme
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

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, max_depth, CS)
  type(param_file_type),   intent(in) :: param_file !< Parameter file
  type(verticalgrid_type), intent(in) :: GV         !< Ocean vertical grid structure
  real,                    intent(in) :: max_depth  !< The maximum depth of the ocean, in 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                         :: check_reconstruction
  logical                         :: check_remapping
  logical                         :: force_bounds_in_subcell
  logical                         :: local_logical

  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")

  ! --- BOUNDARY EXTRAPOLATION --
  ! This sets whether high-order (rather than PCM) reconstruction schemes
  ! should be used within boundary cells
  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", &
                 cs%boundary_extrapolation_for_pressure, &
                 "When defined, the reconstruction is extrapolated\n"//&
                 "within boundary cells rather than assume PCM for the.\n"//&
                 "calculation of pressure. e.g. if PPM is used, a\n"//&
                 "PPM reconstruction will also be used within\n"//&
                 "boundary cells.", default=.true.)

  ! --- PRESSURE GRADIENT CALCULATION ---
  call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", &
                 cs%reconstructForPressure , &
                 "If True, use vertical reconstruction of T/S within\n"//&
                 "the integrals of teh FV pressure gradient calculation.\n"//&
                 "If False, use the constant-by-layer algorithm.\n"//&
                 "By default, this is True when using ALE and False otherwise.", &
                 default=.true. )

  call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", &
                 cs%pressureReconstructionScheme, &
                 "Type of vertical reconstruction of T/S to use in integrals\n"//&
                 "within the FV pressure gradient calculation."//&
                 " 1: PLM reconstruction.\n"//&
                 " 2: PPM reconstruction.", default=pressure_reconstruction_plm)

  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\n"//&
                 "remapping u and v. If false, uses the new method that remaps\n"//&
                 "between grids described by an old and new thickness.", &
                 default=.true.)

  ! Initialize and configure regridding
  call ale_initregridding( gv, 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\n"//&
                 "for vertical remapping for all variables.\n"//&
                 "It can be one of the following schemes:\n"//&
                 trim(remappingschemesdoc), default=remappingdefaultscheme)
  call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
                 "If true, cell-by-cell reconstructions are checked for\n"//&
                 "consistency and if non-monotonicty or an inconsistency is\n"//&
                 "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\n"//&
                 "conservation and new extrema and if an inconsistency is\n"//&
                 "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\n"//&
                 "are forced to be bounded, which might not be the case due to\n"//&
                 "round off.", default=.false.)
  call initialize_remapping( cs%remapCS, string, &
                             boundary_extrapolation=.false., &
                             check_reconstruction=check_reconstruction, &
                             check_remapping=check_remapping, &
                             force_bounds_in_subcell=force_bounds_in_subcell)

  call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &
                 "If true, applies regridding and remapping immediately after\n"//&
                 "initialization so that the state is ALE consistent. This is a\n"//&
                 "legacy step and should not be needed if the initialization is\n"//&
                 "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\n"//&
                 "and the target (new) grid. A short time-scale favors the target\n"//&
                 "grid (0. or anything less than DT_THERM) has no memory of the old\n"//&
                 "grid. A very long time-scale makes the model more Lagrangian.", &
                 units="s", default=0.)
  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\n"//&
                 "final grid exactly matches the target (new) grid.", units="m", default=0.)
  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\n"//&
                 "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and\n"//&
                 "REGRID_FILTER_SHALLOW_DEPTH the filter wieghts adopt a cubic profile.", &
                 units="m", default=0.)
  call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth*gv%m_to_H, &
                                      depth_of_time_filter_deep=filter_deep_depth*gv%m_to_H)
  call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
                 "If true, the regridding ntegrates upwards from the bottom for\n"//&
                 "interface positions, much as the main model does. If false\n"//&
                 "regridding integrates downward, consistant with the remapping\n"//&
                 "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, diag, C_p, Reg, CS)
  type(time_type),target,     intent(in)  :: Time  !< Time structure
  type(ocean_grid_type),      intent(in)  :: G     !< Grid structure
  type(diag_ctrl), target,    intent(in)  :: diag  !< Diagnostics control structure
  real,                       intent(in)  :: C_p   !< seawater heat capacity (J/(kg deg C))
  type(tracer_registry_type), pointer     :: Reg   !< Tracer registry
  type(ale_cs), pointer                   :: CS    !< Module control structure

  integer :: m, ntr, nsize

  if (associated(reg)) then
    ntr = reg%ntr
  else
    ntr = 0
  endif
  nsize = max(1,ntr)

  cs%diag => diag
  cs%C_p  = c_p

  allocate(cs%id_tracer_remap_tendency(nsize))
  allocate(cs%id_Htracer_remap_tendency(nsize))
  allocate(cs%id_Htracer_remap_tendency_2d(nsize))
  allocate(cs%do_tendency_diag(nsize))
  cs%do_tendency_diag(:)             = .false.
  cs%id_tracer_remap_tendency(:)     = -1
  cs%id_Htracer_remap_tendency(:)    = -1
  cs%id_Htracer_remap_tendency_2d(:) = -1

  cs%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,time, &
      'Change in interface height due to ALE regridding', 'meter')

  if(ntr > 0) then

    do m=1,ntr
      if(trim(reg%Tr(m)%name) == 'T') then

        cs%id_tracer_remap_tendency(m) = register_diag_field('ocean_model',                &
        trim(reg%Tr(m)%name)//'_tendency_vert_remap', diag%axesTL, time,                   &
        'Tendency from vertical remapping for tracer concentration '//trim(reg%Tr(m)%name),&
        'degC/s')

        cs%id_Htracer_remap_tendency(m) = register_diag_field('ocean_model',&
        trim(reg%Tr(m)%name)//'h_tendency_vert_remap', diag%axesTL, time,   &
        'Tendency from vertical remapping for heat',                        &
         'W/m2',v_extensive=.true.)

        cs%id_Htracer_remap_tendency_2d(m) = register_diag_field('ocean_model',&
        trim(reg%Tr(m)%name)//'h_tendency_vert_remap_2d', diag%axesT1, time,   &
        'Vertical sum of tendency from vertical remapping for heat',           &
         'W/m2')

      else

        cs%id_tracer_remap_tendency(m) = register_diag_field('ocean_model',                &
        trim(reg%Tr(m)%name)//'_tendency_vert_remap', diag%axesTL, time,                   &
        'Tendency from vertical remapping for tracer concentration '//trim(reg%Tr(m)%name),&
        'tracer conc / sec')

        cs%id_Htracer_remap_tendency(m) = register_diag_field('ocean_model',         &
        trim(reg%Tr(m)%name)//'h_tendency_vert_remap', diag%axesTL, time,            &
        'Tendency from vertical remapping for tracer content '//trim(reg%Tr(m)%name),&
        'kg m-2 s-1',v_extensive=.true.)

        cs%id_Htracer_remap_tendency_2d(m) = register_diag_field('ocean_model',                      &
        trim(reg%Tr(m)%name)//'h_tendency_vert_remap_2d', diag%axesT1, time,                         &
        'Vertical sum of tendency from vertical remapping for tracer content '//trim(reg%Tr(m)%name),&
        'kg m-2 s-1')

      endif

      if(cs%id_tracer_remap_tendency(m)     > 0) cs%do_tendency_diag(m) = .true.
      if(cs%id_Htracer_remap_tendency(m)    > 0) cs%do_tendency_diag(m) = .true.
      if(cs%id_Htracer_remap_tendency_2d(m) > 0) cs%do_tendency_diag(m) = .true.

    enddo   ! m loop over tracers

  endif ! ntr > 0

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 (m or Pa)
  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, h, u, v, tv, Reg, CS, dt, frac_shelf_h)
  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 last time step (m or Pa)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u   !< Zonal velocity field (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v   !< Meridional velocity field (m/s)
  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
  real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main()
  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)) :: h_new ! New 3D grid obtained after last time step (m or Pa)
  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")

  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)")

  ! 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, -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(none) shared(isc,iec,jsc,jec,nk,h,h_new,CS)
  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, 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 last time step (m or Pa)
  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
  real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main()
  ! 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 (m or Pa)
  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, 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(none) shared(isc,iec,jsc,jec,nk,h,h_new,CS)
  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)
  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
  ! 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, 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)
  call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%S, h_new, tv%S)

  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)
  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
                                                                   !! last time step (m or Pa)
  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 (m or Pa)
  type(tracer_registry_type),                 pointer       :: Reg !< Tracer registry structure
  type(ale_cs),                               pointer       :: CS  !< Regridding parameters and options
  ! 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, 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 units)
  real,                                      intent(in) :: threshold !< Value below which to flag issues (H units)
  ! 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 (m or Pa)
  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_orig, tv, n, h_new, u, v)
  type(ale_cs),                               intent(in)    :: 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(in)    :: h_orig !< Original thicknesses
  type(thermo_var_ptrs),                      intent(inout) :: tv     !< Thermo vars (T/S/EOS)
  integer,                                    intent(in)    :: n      !< Number of times to regrid
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(out)   :: h_new  !< Thicknesses after regridding
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u      !< Zonal velocity
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v      !< Meridional velocity

  ! 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 ! A working copy of layer thickesses
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: T, S ! 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, 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
  h(:,:,:) = h_orig(:,:,:)

  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, tv_local, h_new, dzinterface)
    dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)

    ! remap from original grid onto new grid
    ! we need to use remapping_core because there isn't a tracer registry set up in
    ! the state initialization routine
    do j = g%jsc,g%jec ; do i = g%isc,g%iec
      call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h_new(i,j,:), tv_local%S(i,j,:))
      call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h_new(i,j,:), tv_local%T(i,j,:))
    enddo ; enddo


    h(:,:,:) = h_new(:,:,:)
  enddo

  ! save the final temp/salt
  tv%S(:,:,:) = s(:,:,:)
  tv%T(:,:,:) = t(:,:,:)

  ! remap velocities
  call remap_all_state_vars(cs%remapCS, cs, g, gv, h_orig, h_new, null(), dzinttotal, u, v)
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, 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 (m or Pa)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),        intent(in)    :: h_new         !< Thickness of destination grid (m or Pa)
  type(tracer_registry_type),                       pointer       :: Reg           !< Tracer registry structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),optional, intent(in)    :: dxInterface  !< Change in interface position (Hm or Pa)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), optional, intent(inout) :: u          !< Zonal velocity component (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, intent(inout) :: v          !< Meridional velocity component (m/s)
  logical,                                    optional, intent(in)    :: debug      !< If true, show the call tree
  real,                                       optional, intent(in)    :: dt         !< time step for diagnostics
  ! 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, ppt2mks
  real, dimension(GV%ke)                      :: h2
  logical                                     :: show_call_tree

  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

  nz      = gv%ke
  ppt2mks = 0.001

  if (associated(reg)) then
    ntr = reg%ntr
  else
    ntr = 0
  endif

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

  ! Remap tracer
!$OMP parallel default(none) shared(G,GV,h_old,h_new,dxInterface,CS_remapping,nz,Reg,u,v,ntr,show_call_tree, &
!$OMP                               dt,CS_ALE,work_conc,work_cont,work_2d,Idt,ppt2mks) &
!$OMP                       private(h1,h2,dx,u_column)
  if (ntr>0) then
    if (show_call_tree) call calltree_waypoint("remapping tracers (remap_all_state_vars)")
!$OMP do
    do m=1,ntr ! For each tracer

      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, reg%Tr(m)%t(i,j,:), nz, h2, u_column)

            ! Intermediate steps for tendency of tracer concentration and tracer content.
            ! Note: do not merge the two if-tests, since do_tendency_diag(:) is not
            ! allocated during the time=0 initialization call to this routine.
            if(present(dt)) then
              if(cs_ale%do_tendency_diag(m)) then
                do k=1,gv%ke
                  work_conc(i,j,k) = (u_column(k)    - reg%Tr(m)%t(i,j,k)      ) * idt
                  work_cont(i,j,k) = (u_column(k)*h2(k) - reg%Tr(m)%t(i,j,k)*h1(k)) * idt * gv%H_to_kg_m2
                enddo
              endif
            endif

            ! update tracer concentration
            reg%Tr(m)%t(i,j,:) = u_column(:)

          endif

        enddo ! i
      enddo ! j


      ! tendency diagnostics.
      ! Note: do not merge the two if-tests if(present(dt)) and
      ! if(CS_ALE%do_tendency_diag(m)).  The reason is that
      ! do_tendency_diag(:) is not allocated when this routine is called
      ! during initialization (time=0). So need to keep the if-tests split.
      if(present(dt)) then
        if(cs_ale%do_tendency_diag(m)) then

          if(cs_ale%id_tracer_remap_tendency(m) > 0) then
            call post_data(cs_ale%id_tracer_remap_tendency(m), work_conc, cs_ale%diag)
          endif

          if (cs_ale%id_Htracer_remap_tendency(m) > 0 .or. cs_ale%id_Htracer_remap_tendency_2d(m) > 0) then
            if(trim(reg%Tr(m)%name) == 'T') then
              do k=1,gv%ke
                do j = g%jsc,g%jec
                  do i = g%isc,g%iec
                    work_cont(i,j,k) = work_cont(i,j,k) * cs_ale%C_p
                  enddo
                enddo
              enddo
            elseif(trim(reg%Tr(m)%name) == 'S') then
              do k=1,gv%ke
                do j = g%jsc,g%jec
                  do i = g%isc,g%iec
                    work_cont(i,j,k) = work_cont(i,j,k) * ppt2mks
                  enddo
                enddo
              enddo
            endif
          endif

          if (cs_ale%id_Htracer_remap_tendency(m) > 0) then
            call post_data(cs_ale%id_Htracer_remap_tendency(m), work_cont, cs_ale%diag)
          endif
          if (cs_ale%id_Htracer_remap_tendency_2d(m) > 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(cs_ale%id_Htracer_remap_tendency_2d(m), work_2d, cs_ale%diag)
          endif

        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 do
    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
          call remapping_core_h(cs_remapping, nz, h1, u(i,j,:), nz, h2, u_column)
          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 do
    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
          call remapping_core_h(cs_remapping, nz, h1, v(i,j,:), nz, h2, u_column)
          v(i,j,:) = u_column(:)
        endif
      enddo
    enddo
  endif
!$OMP end parallel

  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 )
  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 (m or Pa)
  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 (m or Pa)
  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".
  ! Local variables
  integer :: i, j, k, n_points
  real :: dx(gv%ke+1)
  logical :: ignore_vanished_layers, use_remapping_core_w

  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

!$OMP parallel default(none) shared(CS,G,GV,h_src,s_src,h_dst,s_dst &
!$OMP                               ,ignore_vanished_layers, use_remapping_core_w, nk_src ) &
!$OMP                        firstprivate(n_points,dx)
!$OMP do
  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,:))
        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,:))
        endif
      else
        s_dst(i,j,:) = 0.
      endif
    enddo
  enddo
!$OMP end parallel

end subroutine ale_remap_scalar


!> Use plm reconstruction for pressure gradient (determine edge values)
!! By using a PLM (limited piecewise linear method) reconstruction, this
!! routine determines the edge values for the salinity and temperature
!! within each layer. These edge values are returned and are used to compute
!! the pressure gradient (by computing the densities).
subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h )
  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

  ! Local variables
  integer :: i, j, k
  real    :: hTmp(gv%ke)
  real    :: tmp(gv%ke)
  real, dimension(CS%nk,2)                  :: ppoly_linear_E            !Edge value of polynomial
  real, dimension(CS%nk,CS%degree_linear+1) :: ppoly_linear_coefficients !Coefficients of polynomial

  ! NOTE: the variables 'CS%grid_generic' and 'CS%ppoly_linear' are declared at
  ! the module level.

  ! Determine reconstruction within each column
!$OMP parallel do default(none) shared(G,GV,h,tv,CS,S_t,S_b,T_t,T_b)                  &
!$OMP                          private(hTmp,ppoly_linear_E,ppoly_linear_coefficients,tmp)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
      ! Build current grid
      htmp(:) = h(i,j,:)*gv%H_to_m
      tmp(:) = tv%S(i,j,:)
      ! Reconstruct salinity profile
      ppoly_linear_e = 0.0
      ppoly_linear_coefficients = 0.0
      call plm_reconstruction( gv%ke, htmp, tmp, ppoly_linear_e, ppoly_linear_coefficients )
      if (cs%boundary_extrapolation_for_pressure) call &
        plm_boundary_extrapolation( gv%ke, htmp, tmp, ppoly_linear_e, ppoly_linear_coefficients )

      do k = 1,gv%ke
        s_t(i,j,k) = ppoly_linear_e(k,1)
        s_b(i,j,k) = ppoly_linear_e(k,2)
      end do

      ! Reconstruct temperature profile
      ppoly_linear_e = 0.0
      ppoly_linear_coefficients = 0.0
      tmp(:) = tv%T(i,j,:)
      call plm_reconstruction( gv%ke, htmp, tmp, ppoly_linear_e, ppoly_linear_coefficients )
      if (cs%boundary_extrapolation_for_pressure) call &
        plm_boundary_extrapolation( gv%ke, htmp, tmp, ppoly_linear_e, ppoly_linear_coefficients )

      do k = 1,gv%ke
        t_t(i,j,k) = ppoly_linear_e(k,1)
        t_b(i,j,k) = ppoly_linear_e(k,2)
      end do

    end do
  end do

end subroutine pressure_gradient_plm


!> Use ppm reconstruction for pressure gradient (determine edge values)
!> By using a PPM (limited piecewise linear method) reconstruction, this
!> routine determines the edge values for the salinity and temperature
!> within each layer. These edge values are returned and are used to compute
!> the pressure gradient (by computing the densities).
subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h )
  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 top edge of each layer
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: S_b  !< Salinity at 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   !< ocean thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h    !< layer thickness

  ! Local variables
  integer :: i, j, k
  real    :: hTmp(gv%ke)
  real    :: tmp(gv%ke)
  real, dimension(CS%nk,2) :: &
      ppoly_parab_E            !Edge value of polynomial
  real, dimension(CS%nk,CS%degree_parab+1) :: &
      ppoly_parab_coefficients !Coefficients of polynomial


  ! NOTE: the variables 'CS%grid_generic' and 'CS%ppoly_parab' are declared at
  ! the module level.

  ! Determine reconstruction within each column
!$OMP parallel do default(none) shared(G,GV,h,tv,CS,S_t,S_b,T_t,T_b) &
!$OMP                          private(hTmp,tmp,ppoly_parab_E,ppoly_parab_coefficients)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1

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

      ! Reconstruct salinity profile
      ppoly_parab_e = 0.0
      ppoly_parab_coefficients = 0.0
      call edge_values_implicit_h4( gv%ke, htmp, tmp, ppoly_parab_e )
      call ppm_reconstruction( gv%ke, htmp, tmp, ppoly_parab_e, ppoly_parab_coefficients )
      if (cs%boundary_extrapolation_for_pressure) call &
        ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppoly_parab_e, ppoly_parab_coefficients )

      do k = 1,gv%ke
        s_t(i,j,k) = ppoly_parab_e(k,1)
        s_b(i,j,k) = ppoly_parab_e(k,2)
      end do

      ! Reconstruct temperature profile
      ppoly_parab_e = 0.0
      ppoly_parab_coefficients = 0.0
      tmp(:) = tv%T(i,j,:)
      call edge_values_implicit_h4( gv%ke, htmp, tmp, ppoly_parab_e )
      call ppm_reconstruction( gv%ke, htmp, tmp, ppoly_parab_e, ppoly_parab_coefficients )
      if (cs%boundary_extrapolation_for_pressure) call &
        ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppoly_parab_e, ppoly_parab_coefficients )

      do k = 1,gv%ke
        t_t(i,j,k) = ppoly_parab_e(k,1)
        t_b(i,j,k) = ppoly_parab_e(k,2)
      end do

    end do
  end do

end subroutine pressure_gradient_ppm


!> pressure reconstruction logical
logical function usepressurereconstruction(CS)
  type(ale_cs), pointer :: CS  !< control structure

  if (associated(cs)) then
    usepressurereconstruction=cs%reconstructForPressure
  else
    usepressurereconstruction=.false.
  endif

end function usepressurereconstruction


!> pressure reconstruction integer
integer function pressurereconstructionscheme(CS)
  type(ale_cs), pointer :: CS !< control structure

  if (associated(cs)) then
    pressurereconstructionscheme=cs%pressureReconstructionScheme
  else
    pressurereconstructionscheme=-1
  endif

end function pressurereconstructionscheme

!> Initializes regridding for the main ALE algorithm
subroutine ale_initregridding(GV, max_depth, param_file, mdl, regridCS)
  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)  :: 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.\n"//&
                 "Choose among the following possibilities:\n"//&
                 trim(regriddingcoordinatemodedoc), &
                 default=default_coordinate_mode, fail_if_missing=.true.)

  call initialize_regridding(regridcs, gv, 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 )

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
  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  ! module control structure
  type(verticalgrid_type), pointer :: GV  ! vertical grid information

  integer :: nk

  nk = gv%ke
  gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS )
  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 )
  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

  ! Local variables
  integer :: i, j, k

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

end subroutine ale_initthicknesstocoord

end module mom_ale