MOM_ALE_sponge.F90

!> This module contains the routines used to apply sponge layers when using
!! the ALE mode.
!!
!! Applying sponges requires the following:
!! 1. initialize_ALE_sponge
!! 2. set_up_ALE_sponge_field (tracers) and set_up_ALE_sponge_vel_field (vel)
!! 3. apply_ALE_sponge
!! 4. init_ALE_sponge_diags (not being used for now)
!! 5. ALE_sponge_end (not being used for now)

module mom_ale_sponge


! This file is part of MOM6. See LICENSE.md for the license.
use mom_array_transform, only: rotate_array
use mom_coms, only : sum_across_pes
use mom_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field
use mom_diag_mediator, only : diag_ctrl
use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_horizontal_regridding, only : horiz_interp_and_extrap_tracer
use mom_spatial_means, only : global_i_mean
use mom_time_manager, only : time_type, init_external_field, get_external_field_size, time_interp_external_init
use mom_remapping, only : remapping_cs, remapping_core_h, initialize_remapping
use mom_unit_scaling, only : unit_scale_type
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

!> Store the reference profile at h points for a variable
interface set_up_ale_sponge_field
  module procedure set_up_ale_sponge_field_fixed
  module procedure set_up_ale_sponge_field_varying
end interface

!> This subroutine stores the reference profile at u and v points for a vector
interface set_up_ale_sponge_vel_field
  module procedure set_up_ale_sponge_vel_field_fixed
  module procedure set_up_ale_sponge_vel_field_varying
end interface

!> Ddetermine the number of points which are within sponges in this computational domain.
!!
!! Only points that have positive values of Iresttime and which mask2dT indicates are ocean
!! points are included in the sponges.  It also stores the target interface heights.
interface initialize_ale_sponge
  module procedure initialize_ale_sponge_fixed
  module procedure initialize_ale_sponge_varying
end interface

!  Publicly available functions
public set_up_ale_sponge_field, set_up_ale_sponge_vel_field
public get_ale_sponge_thicknesses, get_ale_sponge_nz_data
public initialize_ale_sponge, apply_ale_sponge, ale_sponge_end, init_ale_sponge_diags
public rotate_ale_sponge, update_ale_sponge_field

! 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.

!> A structure for creating arrays of pointers to 3D arrays with extra gridding information
type :: p3d
  integer :: id !< id for FMS external time interpolator
  integer :: nz_data !< The number of vertical levels in the input field.
  integer :: num_tlevs !< The number of time records contained in the file
  real, dimension(:,:,:), pointer :: mask_in => null() !< pointer to the data mask.
  real, dimension(:,:,:), pointer :: p => null() !< pointer to the data.
  real, dimension(:,:,:), pointer :: h => null() !< pointer to the data grid.
end type p3d

!> A structure for creating arrays of pointers to 2D arrays with extra gridding information
type :: p2d
  integer :: id !< id for FMS external time interpolator
  integer :: nz_data !< The number of vertical levels in the input field
  integer :: num_tlevs !< The number of time records contained in the file
  real, dimension(:,:), pointer :: mask_in => null()!< pointer to the data mask.
  real, dimension(:,:), pointer :: p => null() !< pointer the data.
  real, dimension(:,:), pointer :: h => null() !< pointer the data grid.
end type p2d

!> ALE sponge control structure
type, public :: ale_sponge_cs ; private
  integer :: nz        !< The total number of layers.
  integer :: nz_data   !< The total number of arbritary layers (used by older code).
  integer :: isc       !< The starting i-index of the computational domain at h.
  integer :: iec       !< The ending i-index of the computational domain at h.
  integer :: jsc       !< The starting j-index of the computational domain at h.
  integer :: jec       !< The ending j-index of the computational domain at h.
  integer :: iscb      !< The starting I-index of the computational domain at u/v.
  integer :: iecb      !< The ending I-index of the computational domain at u/v.
  integer :: jscb      !< The starting J-index of the computational domain at u/v.
  integer :: jecb      !< The ending J-index of the computational domain at h.
  integer :: isd       !< The starting i-index of the data domain at h.
  integer :: ied       !< The ending i-index of the data domain at h.
  integer :: jsd       !< The starting j-index of the data domain at h.
  integer :: jed       !< The ending j-index of the data domain at h.
  integer :: num_col   !< The number of sponge tracer points within the computational domain.
  integer :: num_col_u !< The number of sponge u-points within the computational domain.
  integer :: num_col_v !< The number of sponge v-points within the computational domain.
  integer :: fldno = 0 !< The number of fields which have already been
                       !! registered by calls to set_up_sponge_field
  logical :: sponge_uv !< Control whether u and v are included in sponge
  integer, pointer :: col_i(:) => null()   !< Array of the i-indicies of each tracer columns being damped.
  integer, pointer :: col_j(:) => null()   !< Array of the j-indicies of each tracer columns being damped.
  integer, pointer :: col_i_u(:) => null() !< Array of the i-indicies of each u-columns being damped.
  integer, pointer :: col_j_u(:) => null() !< Array of the j-indicies of each u-columns being damped.
  integer, pointer :: col_i_v(:) => null() !< Array of the i-indicies of each v-columns being damped.
  integer, pointer :: col_j_v(:) => null() !< Array of the j-indicies of each v-columns being damped.

  real, pointer :: iresttime_col(:)   => null() !< The inverse restoring time of each tracer column [T-1 ~> s-1].
  real, pointer :: iresttime_col_u(:) => null() !< The inverse restoring time of each u-column [T-1 ~> s-1].
  real, pointer :: iresttime_col_v(:) => null() !< The inverse restoring time of each v-column [T-1 ~> s-1].

  type(p3d) :: var(max_fields_)      !< Pointers to the fields that are being damped.
  type(p2d) :: ref_val(max_fields_) !< The values to which the fields are damped.
  type(p2d) :: ref_val_u  !< The values to which the u-velocities are damped.
  type(p2d) :: ref_val_v  !< The values to which the v-velocities are damped.
  type(p3d) :: var_u  !< Pointer to the u velocities. that are being damped.
  type(p3d) :: var_v  !< Pointer to the v velocities. that are being damped.
  type(p2d) :: ref_h  !< Grid on which reference data is provided (older code).
  type(p2d) :: ref_hu !< u-point grid on which reference data is provided (older code).
  type(p2d) :: ref_hv !< v-point grid on which reference data is provided (older code).

  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                   !! timing of diagnostic output.

  type(remapping_cs) :: remap_cs   !< Remapping parameters and work arrays
  logical :: remap_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.
  logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding
                                   !! that recovers the answers from the end of 2018.  Otherwise, use
                                   !! rotationally symmetric forms of the same expressions.

  logical :: time_varying_sponges  !< True if using newer sponge code
  logical :: spongedataongrid !< True if the sponge data are on the model horizontal grid
end type ale_sponge_cs

contains

!> This subroutine determines the number of points which are within sponges in this computational
!! domain.  Only points that have positive values of Iresttime and which mask2dT indicates are ocean
!! points are included in the sponges.  It also stores the target interface heights.
subroutine initialize_ale_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_data)

  type(ocean_grid_type),            intent(in) :: G !< The ocean's grid structure.
  integer,                          intent(in) :: nz_data !< The total number of sponge input layers.
  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Iresttime !< The inverse of the restoring time [T-1 ~> s-1].
  type(param_file_type),            intent(in) :: param_file !< A structure indicating the open file
                                                             !! to parse for model parameter values.
  type(ale_sponge_cs),              pointer    :: CS !< A pointer that is set to point to the control
                                                     !! structure for this module (in/out).
  real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The thicknesses of the sponge
                                                     !! input layers [H ~> m or kg m-2].


! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_sponge"  ! This module's name.
  logical :: use_sponge
  real, allocatable, dimension(:,:,:) :: data_hu !< thickness at u points [H ~> m or kg m-2]
  real, allocatable, dimension(:,:,:) :: data_hv !< thickness at v points [H ~> m or kg m-2]
  real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1]
  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1]
  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
  logical :: default_2018_answers
  integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
  character(len=10)  :: remapScheme
  if (associated(cs)) then
    call mom_error(warning, "initialize_ALE_sponge_fixed called with an associated "// &
                            "control structure.")
    return
  endif

! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. "//&
                 "The exact location and properties of those sponges are "//&
                 "specified from MOM_initialization.F90.", default=.false.)

  if (.not.use_sponge) return

  allocate(cs)

  call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
                 "Apply sponges in u and v, in addition to tracers.", &
                 default=.false.)

  call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
                 "This sets the reconstruction scheme used "//&
                 " for vertical remapping for all variables.", &
                 default="PLM", do_not_log=.true.)

  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
                 "When defined, a proper high-order reconstruction "//&
                 "scheme is used within boundary cells rather "//&
                 "than PCM. E.g., if PPM is used for remapping, a "//&
                 "PPM reconstruction will also be used within boundary cells.", &
                 default=.false., do_not_log=.true.)
  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%remap_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 get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", cs%hor_regrid_answers_2018, &
                 "If true, use the order of arithmetic for horizonal regridding that recovers "//&
                 "the answers from the end of 2018.  Otherwise, use rotationally symmetric "//&
                 "forms of the same expressions.", default=default_2018_answers)

  cs%time_varying_sponges = .false.
  cs%nz = g%ke
  cs%isc = g%isc ; cs%iec = g%iec ; cs%jsc = g%jsc ; cs%jec = g%jec
  cs%isd = g%isd ; cs%ied = g%ied ; cs%jsd = g%jsd ; cs%jed = g%jed
  cs%iscB = g%iscB ; cs%iecB = g%iecB; cs%jscB = g%jscB ; cs%jecB = g%jecB

  ! number of columns to be restored
  cs%num_col = 0 ; cs%fldno = 0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
      cs%num_col = cs%num_col + 1
  enddo ; enddo

  if (cs%num_col > 0) then
    allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
    allocate(cs%col_i(cs%num_col))         ; cs%col_i = 0
    allocate(cs%col_j(cs%num_col))         ; cs%col_j = 0
    ! pass indices, restoring time to the CS structure
    col = 1
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) then
        cs%col_i(col) = i ; cs%col_j(col) = j
        cs%Iresttime_col(col) = iresttime(i,j)
        col = col +1
      endif
    enddo ; enddo
    ! same for total number of arbritary layers and correspondent data
    cs%nz_data = nz_data
    allocate(cs%Ref_h%p(cs%nz_data,cs%num_col))
    do col=1,cs%num_col ; do k=1,cs%nz_data
      cs%Ref_h%p(k,col) = data_h(cs%col_i(col),cs%col_j(col),k)
    enddo ; enddo
  endif

  total_sponge_cols = cs%num_col
  call sum_across_pes(total_sponge_cols)

! Call the constructor for remapping control structure
  call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation, &
                            answers_2018=cs%remap_answers_2018)

  call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
                 "The total number of columns where sponges are applied at h points.", like_default=.true.)

  if (cs%sponge_uv) then

    allocate(data_hu(g%isdB:g%iedB,g%jsd:g%jed,nz_data)) ; data_hu(:,:,:) = 0.0
    allocate(data_hv(g%isd:g%ied,g%jsdB:g%jedB,nz_data)) ; data_hv(:,:,:) = 0.0
    allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed)) ; iresttime_u(:,:) = 0.0
    allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB)) ; iresttime_v(:,:) = 0.0

    ! u points
    cs%num_col_u = 0 ; !CS%fldno_u = 0
    do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
       data_hu(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:))
       iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
       if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) &
          cs%num_col_u = cs%num_col_u + 1
    enddo ; enddo

    if (cs%num_col_u > 0) then

       allocate(cs%Iresttime_col_u(cs%num_col_u)) ; cs%Iresttime_col_u(:) = 0.0
       allocate(cs%col_i_u(cs%num_col_u))         ; cs%col_i_u(:) = 0
       allocate(cs%col_j_u(cs%num_col_u))         ; cs%col_j_u(:) = 0

       ! pass indices, restoring time to the CS structure
       col = 1
       do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
         if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) then
           cs%col_i_u(col) = i ; cs%col_j_u(col) = j
           cs%Iresttime_col_u(col) = iresttime_u(i,j)
           col = col + 1
         endif
       enddo ; enddo

       ! same for total number of arbritary layers and correspondent data

       allocate(cs%Ref_hu%p(cs%nz_data,cs%num_col_u))
       do col=1,cs%num_col_u ; do k=1,cs%nz_data
         cs%Ref_hu%p(k,col) = data_hu(cs%col_i_u(col),cs%col_j_u(col),k)
       enddo ; enddo
    endif
    total_sponge_cols_u = cs%num_col_u
    call sum_across_pes(total_sponge_cols_u)
    call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
                "The total number of columns where sponges are applied at u points.", like_default=.true.)

    ! v points
    cs%num_col_v = 0 ; !CS%fldno_v = 0
    do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec
      data_hv(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i,j+1,:))
      iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
      if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) &
        cs%num_col_v = cs%num_col_v + 1
    enddo ; enddo

    if (cs%num_col_v > 0) then

      allocate(cs%Iresttime_col_v(cs%num_col_v)) ; cs%Iresttime_col_v = 0.0
      allocate(cs%col_i_v(cs%num_col_v))         ; cs%col_i_v = 0
      allocate(cs%col_j_v(cs%num_col_v))         ; cs%col_j_v = 0

      ! pass indices, restoring time to the CS structure
      col = 1
      do j=cs%jscB,cs%jecB ; do i=cs%isc,cs%iec
        if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) then
          cs%col_i_v(col) = i ; cs%col_j_v(col) = j
          cs%Iresttime_col_v(col) = iresttime_v(i,j)
          col = col + 1
        endif
      enddo ; enddo

      ! same for total number of arbritary layers and correspondent data
      allocate(cs%Ref_hv%p(cs%nz_data,cs%num_col_v))
      do col=1,cs%num_col_v ; do k=1,cs%nz_data
        cs%Ref_hv%p(k,col) = data_hv(cs%col_i_v(col),cs%col_j_v(col),k)
      enddo ; enddo
    endif
    total_sponge_cols_v = cs%num_col_v
    call sum_across_pes(total_sponge_cols_v)
    call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
                 "The total number of columns where sponges are applied at v points.", like_default=.true.)
  endif

end subroutine initialize_ale_sponge_fixed

!> Return the number of layers in the data with a fixed ALE sponge, or 0 if there are
!! no sponge columns on this PE.
function get_ale_sponge_nz_data(CS)
  type(ale_sponge_cs),   pointer       :: cs !< A pointer that is set to point to the control
                                             !! structure for the ALE_sponge module.
  integer :: get_ale_sponge_nz_data  !< The number of layers in the fixed sponge data.

  if (associated(cs)) then
    get_ale_sponge_nz_data = cs%nz_data
  else
    get_ale_sponge_nz_data = 0
  endif
end function get_ale_sponge_nz_data

!> Return the thicknesses used for the data with a fixed ALE sponge
subroutine get_ale_sponge_thicknesses(G, data_h, sponge_mask, CS)
  type(ocean_grid_type), intent(in)    :: g !< The ocean's grid structure (in).
  real, allocatable, dimension(:,:,:), &
                         intent(inout) :: data_h !< The thicknesses of the sponge input layers [H ~> m or kg m-2].
  logical, dimension(SZI_(G),SZJ_(G)), &
                         intent(out)   :: sponge_mask !< A logical mask that is true where
                                                 !! sponges are being applied.
  type(ale_sponge_cs),   pointer       :: cs !< A pointer that is set to point to the control
                                             !! structure for the ALE_sponge module.
  integer :: c, i, j, k

  if (allocated(data_h)) call mom_error(fatal, &
    "get_ALE_sponge_thicknesses called with an allocated data_h.")

  if (.not.associated(cs)) then
    ! There are no sponge points on this PE.
    allocate(data_h(g%isd:g%ied,g%jsd:g%jed,1)) ; data_h(:,:,:) = -1.0
    sponge_mask(:,:) = .false.
    return
  endif

  allocate(data_h(g%isd:g%ied,g%jsd:g%jed,cs%nz_data)) ; data_h(:,:,:) = -1.0
  sponge_mask(:,:) = .false.

  do c=1,cs%num_col
    i = cs%col_i(c) ; j = cs%col_j(c)
    sponge_mask(i,j) = .true.
    do k=1,cs%nz_data
      data_h(i,j,k) = cs%Ref_h%p(k,c)
    enddo
  enddo

end subroutine get_ale_sponge_thicknesses

!> This subroutine determines the number of points which are to be restoref in the computational
!! domain.  Only points that have positive values of Iresttime and which mask2dT indicates are ocean
!! points are included in the sponges.
subroutine initialize_ale_sponge_varying(Iresttime, G, param_file, CS)

  type(ocean_grid_type),            intent(in) :: G !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Iresttime !< The inverse of the restoring time [T-1 ~> s-1].
  type(param_file_type),            intent(in) :: param_file !< A structure indicating the open file to parse
                                                             !! for model parameter values.
  type(ale_sponge_cs),              pointer    :: CS !< A pointer that is set to point to the control
                                                     !! structure for this module (in/out).

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_sponge"  ! This module's name.
  logical :: use_sponge
    real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1]
  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1]
  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
  logical :: default_2018_answers
  logical :: spongeDataOngrid = .false.
  integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
  character(len=10)  :: remapScheme

  if (associated(cs)) then
    call mom_error(warning, "initialize_ALE_sponge_varying called with an associated "// &
                            "control structure.")
    return
  endif
! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. "//&
                 "The exact location and properties of those sponges are "//&
                 "specified from MOM_initialization.F90.", default=.false.)
  if (.not.use_sponge) return
  allocate(cs)
  call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
                 "Apply sponges in u and v, in addition to tracers.", &
                 default=.false.)
  call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
                 "This sets the reconstruction scheme used "//&
                 " for vertical remapping for all variables.", &
                 default="PLM", do_not_log=.true.)
  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
                 "When defined, a proper high-order reconstruction "//&
                 "scheme is used within boundary cells rather "//&
                 "than PCM. E.g., if PPM is used for remapping, a "//&
                 "PPM reconstruction will also be used within boundary cells.", &
                 default=.false., do_not_log=.true.)
  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%remap_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 get_param(param_file, mdl, "SPONGE_DATA_ONGRID", cs%spongeDataOngrid, &
                 "When defined, the incoming sponge data are "//&
                 "assumed to be on the model grid " , &
                 default=.false.)
  cs%time_varying_sponges = .true.
  cs%nz = g%ke
  cs%isc = g%isc ; cs%iec = g%iec ; cs%jsc = g%jsc ; cs%jec = g%jec
  cs%isd = g%isd ; cs%ied = g%ied ; cs%jsd = g%jsd ; cs%jed = g%jed
  cs%iscB = g%iscB ; cs%iecB = g%iecB; cs%jscB = g%jscB ; cs%jecB = g%jecB

  ! number of columns to be restored
  cs%num_col = 0 ; cs%fldno = 0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
      cs%num_col = cs%num_col + 1
  enddo ; enddo
  if (cs%num_col > 0) then
    allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
    allocate(cs%col_i(cs%num_col))         ; cs%col_i = 0
    allocate(cs%col_j(cs%num_col))         ; cs%col_j = 0
    ! pass indices, restoring time to the CS structure
    col = 1
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) then
        cs%col_i(col) = i ; cs%col_j(col) = j
        cs%Iresttime_col(col) = iresttime(i,j)
        col = col +1
      endif
    enddo ; enddo
  endif
  total_sponge_cols = cs%num_col
  call sum_across_pes(total_sponge_cols)

! Call the constructor for remapping control structure
  call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation, &
                            answers_2018=cs%remap_answers_2018)
  call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
                 "The total number of columns where sponges are applied at h points.", like_default=.true.)
  if (cs%sponge_uv) then
    allocate(iresttime_u(g%isdB:g%iedB,g%jsd:g%jed)) ; iresttime_u(:,:) = 0.0
    allocate(iresttime_v(g%isd:g%ied,g%jsdB:g%jedB)) ; iresttime_v(:,:) = 0.0
    ! u points
    cs%num_col_u = 0 ; !CS%fldno_u = 0
    do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB
      iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
      if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) &
        cs%num_col_u = cs%num_col_u + 1
    enddo ; enddo
    if (cs%num_col_u > 0) then
      allocate(cs%Iresttime_col_u(cs%num_col_u)) ; cs%Iresttime_col_u = 0.0
      allocate(cs%col_i_u(cs%num_col_u))         ; cs%col_i_u = 0
      allocate(cs%col_j_u(cs%num_col_u))         ; cs%col_j_u = 0
      ! pass indices, restoring time to the CS structure
      col = 1
      do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
        if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) then
          cs%col_i_u(col) = i ; cs%col_j_u(col) = j
          cs%Iresttime_col_u(col) = iresttime_u(i,j)
          col = col +1
        endif
      enddo ; enddo
      ! same for total number of arbritary layers and correspondent data
    endif
    total_sponge_cols_u = cs%num_col_u
    call sum_across_pes(total_sponge_cols_u)
    call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
                "The total number of columns where sponges are applied at u points.", like_default=.true.)
    ! v points
    cs%num_col_v = 0 ; !CS%fldno_v = 0
    do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec
      iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
      if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) &
        cs%num_col_v = cs%num_col_v + 1
    enddo ; enddo
    if (cs%num_col_v > 0) then
      allocate(cs%Iresttime_col_v(cs%num_col_v)) ; cs%Iresttime_col_v = 0.0
      allocate(cs%col_i_v(cs%num_col_v))         ; cs%col_i_v = 0
      allocate(cs%col_j_v(cs%num_col_v))         ; cs%col_j_v = 0
      ! pass indices, restoring time to the CS structure
      col = 1
      do j=cs%jscB,cs%jecB ; do i=cs%isc,cs%iec
        if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) then
          cs%col_i_v(col) = i ; cs%col_j_v(col) = j
          cs%Iresttime_col_v(col) = iresttime_v(i,j)
          col = col +1
        endif
      enddo ; enddo
    endif
    total_sponge_cols_v = cs%num_col_v
    call sum_across_pes(total_sponge_cols_v)
    call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
                "The total number of columns where sponges are applied at v points.", like_default=.true.)
  endif

end subroutine initialize_ale_sponge_varying

!> Initialize diagnostics for the ALE_sponge module.
! GMM: this routine is not being used for now.
subroutine init_ale_sponge_diags(Time, G, diag, CS)
  type(time_type), target, intent(in)    :: time !< The current model time
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(ale_sponge_cs),     pointer       :: cs   !< ALE sponge control structure

  if (.not.associated(cs)) return

  cs%diag => diag

end subroutine init_ale_sponge_diags

!> This subroutine stores the reference profile at h points for the variable
!! whose address is given by f_ptr.
subroutine set_up_ale_sponge_field_fixed(sp_val, G, f_ptr, CS)
  type(ocean_grid_type), intent(in) :: G  !< Grid structure
  type(ale_sponge_cs),   pointer    :: CS !< ALE sponge control structure (in/out).
  real, dimension(SZI_(G),SZJ_(G),CS%nz_data), &
                         intent(in) :: sp_val !< Field to be used in the sponge, it has arbitrary number of layers.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                 target, intent(in) :: f_ptr !< Pointer to the field to be damped

  integer :: j, k, col
  character(len=256) :: mesg ! String for error messages

  if (.not.associated(cs)) return

  cs%fldno = cs%fldno + 1

  if (cs%fldno > max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
           &the number of fields to be damped in the call to &
           &initialize_ALE_sponge." )') cs%fldno
    call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
  endif

  ! stores the reference profile
  allocate(cs%Ref_val(cs%fldno)%p(cs%nz_data,cs%num_col))
  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
  do col=1,cs%num_col
    do k=1,cs%nz_data
      cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
    enddo
  enddo

  cs%var(cs%fldno)%p => f_ptr

end subroutine set_up_ale_sponge_field_fixed

!> This subroutine stores the reference profile at h points for the variable
!! whose address is given by filename and fieldname.
subroutine set_up_ale_sponge_field_varying(filename, fieldname, Time, G, GV, US, f_ptr, CS)
  character(len=*),        intent(in) :: filename !< The name of the file with the
                                                  !! time varying field data
  character(len=*),        intent(in) :: fieldname !< The name of the field in the file
                                                  !! with the time varying field data
  type(time_type),         intent(in) :: Time  !< The current model time
  type(ocean_grid_type),   intent(in) :: G     !< Grid structure (in).
  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)), &
                   target, intent(in) :: f_ptr !< Pointer to the field to be damped (in).
  type(ale_sponge_cs),     pointer    :: CS    !< Sponge control structure (in/out).

  ! Local variables
  real, allocatable, dimension(:,:,:) :: sp_val !< Field to be used in the sponge
  real, allocatable, dimension(:,:,:) :: mask_z !< Field mask for the sponge data
  real, allocatable, dimension(:), target :: z_in, z_edges_in ! Heights [Z ~> m].
  real :: missing_value
  integer :: j, k, col
  integer :: isd,ied,jsd,jed
  integer :: nPoints
  integer, dimension(4) :: fld_sz
  integer :: nz_data !< the number of vertical levels in this input field
  character(len=256) :: mesg ! String for error messages
  ! Local variables for ALE remapping
  real, dimension(:), allocatable :: tmpT1d
  real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m].
  type(remapping_cs) :: remapCS ! Remapping parameters and work arrays

  if (.not.associated(cs)) return
  ! initialize time interpolator module
  call time_interp_external_init()
  isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed
  cs%fldno = cs%fldno + 1
  if (cs%fldno > max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
           &the number of fields to be damped in the call to &
           &initialize_ALE_sponge." )') cs%fldno
    call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
  endif
  ! get a unique time interp id for this field. If sponge data is ongrid, then setup
  ! to only read on the computational domain
  if (cs%spongeDataOngrid) then
    cs%Ref_val(cs%fldno)%id = init_external_field(filename, fieldname,domain=g%Domain%mpp_domain)
  else
    cs%Ref_val(cs%fldno)%id = init_external_field(filename, fieldname)
  endif
  fld_sz(1:4)=-1
  fld_sz = get_external_field_size(cs%Ref_val(cs%fldno)%id)
  nz_data = fld_sz(3)
  cs%Ref_val(cs%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid
  cs%Ref_val(cs%fldno)%num_tlevs = fld_sz(4)
  ! initializes the target profile array for this field
  ! for all columns which will be masked
  allocate(cs%Ref_val(cs%fldno)%p(nz_data,cs%num_col))
  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
  allocate( cs%Ref_val(cs%fldno)%h(nz_data,cs%num_col) )
  cs%Ref_val(cs%fldno)%h(:,:) = 0.0
  cs%var(cs%fldno)%p => f_ptr


end subroutine set_up_ale_sponge_field_varying

!> This subroutine stores the reference profile at u and v points for the variable
!! whose address is given by u_ptr and v_ptr.
subroutine set_up_ale_sponge_vel_field_fixed(u_val, v_val, G, u_ptr, v_ptr, CS)
  type(ocean_grid_type), intent(in) :: G  !< Grid structure (in).
  type(ale_sponge_cs),   pointer    :: CS !< Sponge structure (in/out).
  real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), &
                         intent(in) :: u_val !< u field to be used in the sponge, it has arbritary number of layers.
  real, dimension(SZI_(G),SZJB_(G),CS%nz_data), &
                         intent(in) :: v_val !< v field to be used in the sponge, it has arbritary number of layers.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped

  integer :: j, k, col
  character(len=256) :: mesg ! String for error messages

  if (.not.associated(cs)) return

  ! stores the reference profile
  allocate(cs%Ref_val_u%p(cs%nz_data,cs%num_col_u))
  cs%Ref_val_u%p(:,:) = 0.0
  do col=1,cs%num_col_u
    do k=1,cs%nz_data
      cs%Ref_val_u%p(k,col) = u_val(cs%col_i_u(col),cs%col_j_u(col),k)
    enddo
  enddo
  cs%var_u%p => u_ptr
  allocate(cs%Ref_val_v%p(cs%nz_data,cs%num_col_v))
  cs%Ref_val_v%p(:,:) = 0.0
  do col=1,cs%num_col_v
    do k=1,cs%nz_data
      cs%Ref_val_v%p(k,col) = v_val(cs%col_i_v(col),cs%col_j_v(col),k)
    enddo
  enddo
  cs%var_v%p => v_ptr

end subroutine set_up_ale_sponge_vel_field_fixed

!> This subroutine stores the reference profile at uand v points for the variable
!! whose address is given by u_ptr and v_ptr.
subroutine set_up_ale_sponge_vel_field_varying(filename_u, fieldname_u, filename_v, fieldname_v, &
                                               Time, G, US, CS, u_ptr, v_ptr)
  character(len=*), intent(in)    :: filename_u  !< File name for u field
  character(len=*), intent(in)    :: fieldname_u !< Name of u variable in file
  character(len=*), intent(in)    :: filename_v  !< File name for v field
  character(len=*), intent(in)    :: fieldname_v !< Name of v variable in file
  type(time_type),  intent(in)    :: Time        !< Model time
  type(ocean_grid_type), intent(inout) :: G      !< Ocean grid (in)
  type(unit_scale_type), intent(in)    :: US     !< A dimensional unit scaling type
  type(ale_sponge_cs), pointer    :: CS          !< Sponge structure (in/out).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped (in).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped (in).
  ! Local variables
  real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge.
  real, allocatable, dimension(:,:,:) :: mask_u !< U field mask for the sponge data.
  real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge.
  real, allocatable, dimension(:,:,:) :: mask_v !< V field mask for the sponge data.

  real, allocatable, dimension(:), target :: z_in, z_edges_in
  real :: missing_value

  integer :: j, k, col
  integer :: isd, ied, jsd, jed
  integer :: isdB, iedB, jsdB, jedB
  integer, dimension(4) :: fld_sz
  character(len=256) :: mesg ! String for error messages

  if (.not.associated(cs)) return

  isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed
  isdb = g%isdB; iedb = g%iedB; jsdb = g%jsdB; jedb = g%jedB
  ! get a unique id for this field which will allow us to return an array
  ! containing time-interpolated values from an external file corresponding
  ! to the current model date.
  cs%Ref_val_u%id = init_external_field(filename_u, fieldname_u)
  fld_sz(1:4)=-1
  fld_sz = get_external_field_size(cs%Ref_val_u%id)
  cs%Ref_val_u%nz_data = fld_sz(3)
  cs%Ref_val_u%num_tlevs = fld_sz(4)
  cs%Ref_val_v%id = init_external_field(filename_v, fieldname_v)
  fld_sz(1:4)=-1
  fld_sz = get_external_field_size(cs%Ref_val_v%id)
  cs%Ref_val_v%nz_data = fld_sz(3)
  cs%Ref_val_v%num_tlevs = fld_sz(4)
  allocate( u_val(isdb:iedb,jsd:jed, fld_sz(3)) )
  allocate( mask_u(isdb:iedb,jsd:jed, fld_sz(3)) )
  allocate( v_val(isd:ied,jsdb:jedb, fld_sz(3)) )
  allocate( mask_v(isd:ied,jsdb:jedb, fld_sz(3)) )
  ! Interpolate external file data to the model grid
  ! I am hard-wiring this call to assume that the input grid is zonally re-entrant
  ! In the future, this should be generalized using an interface to return the
  ! modulo attribute of the zonal axis (mjh).
  call horiz_interp_and_extrap_tracer(cs%Ref_val_u%id, time, 1.0, g, u_val, mask_u, z_in, z_edges_in, &
                                      missing_value, .true., .false., .false., m_to_z=us%m_to_Z, &
                                      answers_2018=cs%hor_regrid_answers_2018)
  !!! TODO: add a velocity interface! (mjh)
  ! Interpolate external file data to the model grid
  ! I am hard-wiring this call to assume that the input grid is zonally re-entrant
  ! In the future, this should be generalized using an interface to return the
  ! modulo attribute of the zonal axis (mjh).
  call horiz_interp_and_extrap_tracer(cs%Ref_val_v%id, time, 1.0, g, v_val, mask_v, z_in, z_edges_in, &
                                      missing_value, .true., .false., .false., m_to_z=us%m_to_Z, &
                                      answers_2018=cs%hor_regrid_answers_2018)
  ! stores the reference profile
  allocate(cs%Ref_val_u%p(fld_sz(3),cs%num_col_u))
  cs%Ref_val_u%p(:,:) = 0.0
  do col=1,cs%num_col_u
    do k=1,fld_sz(3)
      cs%Ref_val_u%p(k,col) = u_val(cs%col_i_u(col),cs%col_j_u(col),k)
    enddo
  enddo
  cs%var_u%p => u_ptr
  allocate(cs%Ref_val_v%p(fld_sz(3),cs%num_col_v))
  cs%Ref_val_v%p(:,:) = 0.0
  do col=1,cs%num_col_v
    do k=1,fld_sz(3)
      cs%Ref_val_v%p(k,col) = v_val(cs%col_i_v(col),cs%col_j_v(col),k)
    enddo
  enddo
  cs%var_v%p => v_ptr

end subroutine set_up_ale_sponge_vel_field_varying

!> This subroutine applies damping to the layers thicknesses, temp, salt and a variety of tracers
!! for every column where there is damping.
subroutine apply_ale_sponge(h, dt, G, GV, US, CS, Time)
  type(ocean_grid_type),     intent(inout) :: g  !< The ocean's grid structure (in).
  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  !< Layer thickness [H ~> m or kg m-2] (in)
  real,                      intent(in)    :: dt !< The amount of time covered by this call [T ~> s].
  type(ale_sponge_cs),       pointer       :: cs !< A pointer to the control structure for this module
                                                 !! that is set by a previous call to initialize_ALE_sponge (in).
  type(time_type), optional, intent(in)    :: time !< The current model date

  real :: damp                                  ! The timestep times the local damping coefficient [nondim].
  real :: i1pdamp                               ! I1pdamp is 1/(1 + damp). [nondim].
  real :: m_to_z                                ! A unit conversion factor from m to Z.
  real, allocatable, dimension(:) :: tmp_val2   ! data values on the original grid
  real, dimension(SZK_(G)) :: tmp_val1          ! data values remapped to model grid
  real :: hu(szib_(g), szj_(g), szk_(g))        ! A temporary array for h at u pts
  real :: hv(szi_(g), szjb_(g), szk_(g))        ! A temporary array for h at v pts
  real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields
  real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts
  real, dimension(:), allocatable :: hsrc       ! Source thicknesses [Z ~> m].
  ! Local variables for ALE remapping
  real, dimension(:), allocatable :: tmpt1d
  integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz, nz_data
  integer :: col, total_sponge_cols
  real, allocatable, dimension(:), target :: z_in, z_edges_in
  real :: missing_value
  real :: h_neglect, h_neglect_edge
  real :: ztopofcell, zbottomofcell ! Heights [Z ~> m].
  integer :: npoints

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  if (.not.associated(cs)) return

  if (.not.cs%remap_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

  if (cs%time_varying_sponges) then
    if (.not. present(time)) &
      call mom_error(fatal,"apply_ALE_sponge: No time information provided")
    do m=1,cs%fldno
       nz_data = cs%Ref_val(m)%nz_data
       allocate(sp_val(g%isd:g%ied,g%jsd:g%jed,1:nz_data))
       allocate(mask_z(g%isd:g%ied,g%jsd:g%jed,1:nz_data))
       sp_val(:,:,:) = 0.0
       mask_z(:,:,:) = 0.0
       call horiz_interp_and_extrap_tracer(cs%Ref_val(m)%id, time, 1.0, g, sp_val, mask_z, z_in, &
                      z_edges_in,  missing_value, .true., .false., .false., &
                      spongeongrid=cs%SpongeDataOngrid, m_to_z=us%m_to_Z, &
                      answers_2018=cs%hor_regrid_answers_2018)
       allocate( hsrc(nz_data) )
       allocate( tmpt1d(nz_data) )
       do c=1,cs%num_col
          i = cs%col_i(c) ; j = cs%col_j(c)
          cs%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
          ! Build the source grid
          ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0; hsrc(:) = 0.0; tmpt1d(:) = -99.9
          do k=1,nz_data
             if (mask_z(cs%col_i(c),cs%col_j(c),k) == 1.0) then
                zbottomofcell = -min( z_edges_in(k+1), g%bathyT(cs%col_i(c),cs%col_j(c)) )
                tmpt1d(k) = sp_val(cs%col_i(c),cs%col_j(c),k)
             elseif (k>1) then
                zbottomofcell = -g%bathyT(cs%col_i(c),cs%col_j(c))
                tmpt1d(k) = tmpt1d(k-1)
             else ! This next block should only ever be reached over land
                tmpt1d(k) = -99.9
             endif
             hsrc(k) = ztopofcell - zbottomofcell
             if (hsrc(k)>0.) npoints = npoints + 1
             ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
          enddo
          ! In case data is deeper than model
          hsrc(nz_data) = hsrc(nz_data) + ( ztopofcell + g%bathyT(cs%col_i(c),cs%col_j(c)) )
          cs%Ref_val(m)%h(1:nz_data,c) = gv%Z_to_H*hsrc(1:nz_data)
          cs%Ref_val(m)%p(1:nz_data,c) = tmpt1d(1:nz_data)
          do k=2,nz_data
             !          if (mask_z(i,j,k)==0.) &
             if (cs%Ref_val(m)%h(k,c) <= 0.001*gv%m_to_H) &
                  ! some confusion here about why the masks are not correct returning from horiz_interp
                  ! reverting to using a minimum thickness criteria
                  cs%Ref_val(m)%p(k,c) = cs%Ref_val(m)%p(k-1,c)
          enddo
      enddo
      deallocate(sp_val, mask_z, hsrc, tmpt1d)
    enddo
  else
    nz_data = cs%nz_data
  endif

  allocate(tmp_val2(nz_data))
  do m=1,cs%fldno
    do c=1,cs%num_col
! c is an index for the next 3 lines but a multiplier for the rest of the loop
! Therefore we use c as per C code and increment the index where necessary.
      i = cs%col_i(c) ; j = cs%col_j(c)
      damp = dt * cs%Iresttime_col(c)
      i1pdamp = 1.0 / (1.0 + damp)
      tmp_val2(1:nz_data) = cs%Ref_val(m)%p(1:nz_data,c)
      if (cs%time_varying_sponges) then
        call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val(m)%h(1:nz_data,c), tmp_val2, &
                              cs%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      else
        call remapping_core_h(cs%remap_cs,nz_data, cs%Ref_h%p(1:nz_data,c), tmp_val2, &
                              cs%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      endif
      !Backward Euler method
      cs%var(m)%p(i,j,1:cs%nz) = i1pdamp * (cs%var(m)%p(i,j,1:cs%nz) + tmp_val1 * damp)
    enddo
  enddo

  ! for debugging
  !c=CS%num_col
  !do m=1,CS%fldno
  !   write(*,*) 'APPLY SPONGE,m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1',&
  !               m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1
  !enddo

  if (cs%sponge_uv) then
    ! u points
    do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB; do k=1,nz
      hu(i,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
    enddo ; enddo ; enddo
    if (cs%time_varying_sponges) then
      if (.not. present(time)) &
         call mom_error(fatal,"apply_ALE_sponge: No time information provided")

      nz_data = cs%Ref_val_u%nz_data
      allocate(sp_val(g%isdB:g%iedB,g%jsd:g%jed,1:nz_data))
      allocate(mask_z(g%isdB:g%iedB,g%jsd:g%jed,1:nz_data))
      ! Interpolate from the external horizontal grid and in time
      call horiz_interp_and_extrap_tracer(cs%Ref_val_u%id, time, 1.0, g, sp_val, mask_z, z_in, &
                                          z_edges_in, missing_value, .true., .false., .false., &
                                          m_to_z=us%m_to_Z, answers_2018=cs%hor_regrid_answers_2018)
!     call pass_var(sp_val,G%Domain)
!     call pass_var(mask_z,G%Domain)
      do c=1,cs%num_col
        ! c is an index for the next 3 lines but a multiplier for the rest of the loop
        ! Therefore we use c as per C code and increment the index where necessary.
        i = cs%col_i(c) ; j = cs%col_j(c)
        cs%Ref_val_u%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
      enddo

      deallocate (sp_val, mask_z)

      nz_data = cs%Ref_val_v%nz_data
      allocate(sp_val(g%isd:g%ied,g%jsdB:g%jedB,1:nz_data))
      allocate(mask_z(g%isd:g%ied,g%jsdB:g%jedB,1:nz_data))
      ! Interpolate from the external horizontal grid and in time
      call horiz_interp_and_extrap_tracer(cs%Ref_val_v%id, time, 1.0, g, sp_val, mask_z, z_in, &
                                          z_edges_in, missing_value, .true., .false., .false., &
                                          m_to_z=us%m_to_Z, answers_2018=cs%hor_regrid_answers_2018)

!     call pass_var(sp_val,G%Domain)
!     call pass_var(mask_z,G%Domain)

      do c=1,cs%num_col
        ! c is an index for the next 3 lines but a multiplier for the rest of the loop
        ! Therefore we use c as per C code and increment the index where necessary.
        i = cs%col_i(c) ; j = cs%col_j(c)
        cs%Ref_val_v%p(1:nz_data,c) = sp_val(i,j,1:nz_data)
      enddo

      deallocate (sp_val, mask_z)

    else
      nz_data = cs%nz_data
    endif

    do c=1,cs%num_col_u
      i = cs%col_i_u(c) ; j = cs%col_j_u(c)
      damp = dt * cs%Iresttime_col_u(c)
      i1pdamp = 1.0 / (1.0 + damp)
      if (cs%time_varying_sponges) nz_data = cs%Ref_val(m)%nz_data
      tmp_val2(1:nz_data) = cs%Ref_val_u%p(1:nz_data,c)
      if (cs%time_varying_sponges) then
        call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_val_u%h(:,c), tmp_val2, &
                 cs%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      else
        call remapping_core_h(cs%remap_cs, nz_data, cs%Ref_hu%p(:,c), tmp_val2, &
                 cs%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      endif
      !Backward Euler method
      cs%var_u%p(i,j,:) = i1pdamp * (cs%var_u%p(i,j,:) + tmp_val1 * damp)
    enddo

    ! v points
    do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec; do k=1,nz
      hv(i,j,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
    enddo ; enddo ; enddo

    do c=1,cs%num_col_v
      i = cs%col_i_v(c) ; j = cs%col_j_v(c)
      damp = dt * cs%Iresttime_col_v(c)
      i1pdamp = 1.0 / (1.0 + damp)
      tmp_val2(1:nz_data) = cs%Ref_val_v%p(1:nz_data,c)
      if (cs%time_varying_sponges) then
        call remapping_core_h(cs%remap_cs, cs%nz_data, cs%Ref_val_v%h(:,c), tmp_val2, &
                 cs%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      else
        call remapping_core_h(cs%remap_cs, cs%nz_data, cs%Ref_hv%p(:,c), tmp_val2, &
                 cs%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge)
      endif
      !Backward Euler method
      cs%var_v%p(i,j,:) = i1pdamp * (cs%var_v%p(i,j,:) + tmp_val1 * damp)
    enddo

  endif

  deallocate(tmp_val2)

end subroutine apply_ale_sponge

!> Rotate the ALE sponge fields from the input to the model index map.
subroutine rotate_ale_sponge(sponge_in, G_in, sponge, G, turns, param_file)
  type(ale_sponge_cs),   intent(in) :: sponge_in !< The control structure for this module with the
                                                 !! original grid rotation
  type(ocean_grid_type), intent(in) :: g_in      !< The ocean's grid structure with the original rotation.
  type(ale_sponge_cs),   pointer    :: sponge    !< A pointer to the control that will be set up with
                                                 !! the new grid rotation
  type(ocean_grid_type), intent(in) :: g         !< The ocean's grid structure with the new rotation.
  integer,               intent(in) :: turns     !< The number of 90-degree turns between grids
  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
                                                 !! to parse for model parameter values.

  ! First part: Index construction
  !   1. Reconstruct Iresttime(:,:) from sponge_in
  !   2. rotate Iresttime(:,:)
  !   3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:)
  ! All the index adjustment should follow from the Iresttime rotation

  real, dimension(:,:), allocatable :: iresttime_in, iresttime
  real, dimension(:,:,:), allocatable :: data_h_in, data_h
  real, dimension(:,:,:), allocatable :: sp_val_in, sp_val
  real, dimension(:,:,:), pointer :: sp_ptr => null()
  integer :: c, c_i, c_j
  integer :: k, nz_data
  integer :: n
  logical :: fixed_sponge

  fixed_sponge = .not. sponge_in%time_varying_sponges
  ! NOTE: nz_data is only conditionally set when fixed_sponge is true.

  allocate(iresttime_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed))
  allocate(iresttime(g%isd:g%ied, g%jsd:g%jed))
  iresttime_in(:,:) = 0.0

  if (fixed_sponge) then
    nz_data = sponge_in%nz_data
    allocate(data_h_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz_data))
    allocate(data_h(g%isd:g%ied, g%jsd:g%jed, nz_data))
    data_h_in(:,:,:) = 0.
  endif

  ! Re-populate the 2D Iresttime and data_h arrays on the original grid
  do c=1,sponge_in%num_col
    c_i = sponge_in%col_i(c)
    c_j = sponge_in%col_j(c)
    iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c)
    if (fixed_sponge) then
      do k = 1, nz_data
        data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c)
      enddo
    endif
  enddo

  call rotate_array(iresttime_in, turns, iresttime)
  if (fixed_sponge) then
    call rotate_array(data_h_in, turns, data_h)
    call initialize_ale_sponge_fixed(iresttime, g, param_file, sponge, &
                                     data_h, nz_data)
  else
    call initialize_ale_sponge_varying(iresttime, g, param_file, sponge)
  endif

  deallocate(iresttime_in)
  deallocate(iresttime)
  if (fixed_sponge) then
    deallocate(data_h_in)
    deallocate(data_h)
  endif

  ! Second part: Provide rotated fields for which relaxation is applied

  sponge%fldno = sponge_in%fldno

  if (fixed_sponge) then
    allocate(sp_val_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz_data))
    allocate(sp_val(g%isd:g%ied, g%jsd:g%jed, nz_data))
  endif

  do n=1,sponge_in%fldno
    ! Assume that tracers are pointers and are remapped in other functions(?)
    sp_ptr => sponge_in%var(n)%p
    if (fixed_sponge) then
      sp_val_in(:,:,:) = 0.0
      do c = 1, sponge_in%num_col
        c_i = sponge_in%col_i(c)
        c_j = sponge_in%col_j(c)
        do k = 1, nz_data
          sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c)
        enddo
      enddo

      call rotate_array(sp_val_in, turns, sp_val)

      ! NOTE: This points sp_val with the unrotated field.  See note below.
      call set_up_ale_sponge_field(sp_val, g, sp_ptr, sponge)

      deallocate(sp_val_in)
    else
      ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying()
      ! (time_interp_external_init, init_external_field, etc), so we manually
      ! do a portion of this function below.
      sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id
      sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs

      nz_data = sponge_in%Ref_val(n)%nz_data
      sponge%Ref_val(n)%nz_data = nz_data

      allocate(sponge%Ref_val(n)%p(nz_data, sponge_in%num_col))
      allocate(sponge%Ref_val(n)%h(nz_data, sponge_in%num_col))
      sponge%Ref_val(n)%p(:,:) = 0.0
      sponge%Ref_val(n)%h(:,:) = 0.0

      ! TODO: There is currently no way to associate a generic field pointer to
      !   its rotated equivalent without introducing a new data structure which
      !   explicitly tracks the pairing.
      !
      !   As a temporary fix, we store the pointer to the unrotated field in
      !   the rotated sponge, and use this reference to replace the pointer
      !   to the rotated field update_ALE_sponge field.
      !
      !   This makes a lot of unverifiable assumptions, and should not be
      !   considered the final solution.
      sponge%var(n)%p => sp_ptr
    endif
  enddo

  ! TODO: var_u and var_v sponge dampling is not yet supported.
  if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) &
    call mom_error(fatal, "Rotation of ALE sponge velocities is not yet " &
      // "implemented.")

  ! Transfer any existing diag_CS reference pointer
  sponge%diag => sponge_in%diag

  ! NOTE: initialize_ALE_sponge_* resolves remap_cs
end subroutine rotate_ale_sponge


!> Scan the ALE sponge variables and replace a prescribed pointer to a new value.
! TODO: This function solely exists to replace field pointers in the sponge
!   after rotation.  This function is part of a temporary solution until
!   something more robust is developed.
subroutine update_ale_sponge_field(sponge, p_old, G, GV, p_new)
  type(ale_sponge_cs),     pointer    :: sponge !< A pointer to the control structure for this module
                                               !! that is set by a previous call to initialize_ALE_sponge.
  real, dimension(:,:,:), &
                   target, intent(in) :: p_old !< The previous array of target values
  type(ocean_grid_type),   intent(in) :: g     !< The updated ocean grid structure
  type(verticalgrid_type), intent(in) :: gv    !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                   target, intent(in) :: p_new !< The new array of target values

  integer :: n

  do n=1,sponge%fldno
    if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new
  enddo

end subroutine update_ale_sponge_field


! GMM: I could not find where sponge_end is being called, but I am keeping
!  ALE_sponge_end here so we can add that if needed.
!> This subroutine deallocates any memory associated with the ALE_sponge module.
subroutine ale_sponge_end(CS)
  type(ale_sponge_cs), pointer :: cs !< A pointer to the control structure that is
                                     !! set by a previous call to initialize_ALE_sponge.

  integer :: m

  if (.not.associated(cs)) return

  if (associated(cs%col_i)) deallocate(cs%col_i)
  if (associated(cs%col_i_u)) deallocate(cs%col_i_u)
  if (associated(cs%col_i_v)) deallocate(cs%col_i_v)
  if (associated(cs%col_j)) deallocate(cs%col_j)
  if (associated(cs%col_j_u)) deallocate(cs%col_j_u)
  if (associated(cs%col_j_v)) deallocate(cs%col_j_v)

  if (associated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
  if (associated(cs%Iresttime_col_u)) deallocate(cs%Iresttime_col_u)
  if (associated(cs%Iresttime_col_v)) deallocate(cs%Iresttime_col_v)

  do m=1,cs%fldno
    if (associated(cs%Ref_val(m)%p)) deallocate(cs%Ref_val(m)%p)
  enddo

  deallocate(cs)

end subroutine ale_sponge_end

end module mom_ale_sponge