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_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
! GMM - Planned extension:  Support for time varying sponge targets.
 
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
 
! 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 [s-1].
  real, pointer :: iresttime_col_u(:) => null() !< The inverse restoring time of each u-column [s-1].
  real, pointer :: iresttime_col_v(:) => null() !< The inverse restoring time of each v-column [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 :: new_sponges  !< True if using newer sponge code
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 [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 [s-1]
  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [s-1]
  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
  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_sponge 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.)
 
  cs%new_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)
 
  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.")
 
  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.")
 
    ! 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.")
  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 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.
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 [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 [s-1]
  real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [s-1]
  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
  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_sponge 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.)
 
  cs%new_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)
 
  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.")
 
  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.")
 
     ! 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.")
  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 arbritary 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_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, 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
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                   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 :: hsrc ! Source thicknesses [Z ~> m].
  real, dimension(:), allocatable :: tmpT1d
  real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m].
  type(remapping_cs) :: remapCS ! Remapping parameters and work arrays
 
  if (.not.associated(cs)) return
 
  ! Call this in case it was not previously done.
  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_sponge." )') cs%fldno
    call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
  endif
 
  ! 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(cs%fldno)%id = init_external_field(filename, fieldname)
  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 !< each individual sponge field is assumed to reside on a different grid
  cs%Ref_val(cs%fldno)%num_tlevs = fld_sz(4)
 
  allocate( sp_val(isd:ied,jsd:jed, nz_data) )
  allocate( mask_z(isd:ied,jsd:jed, nz_data) )
 
  ! initializes the current reference profile array
  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
 
  ! 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(CS%fldno)%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)
 
! Do not think halo updates are needed (mjh)
!  call pass_var(sp_val,G%Domain)
!  call pass_var(mask_z,G%Domain)
 
! Done with horizontal interpolation.
! Now remap to model coordinates
! First we reserve a work space for reconstructions of the source data
  allocate( hsrc(nz_data) )
  allocate( tmpt1d(nz_data) )
 
  do col=1,cs%num_col
    ! 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(col),cs%col_j(col),k) == 1.0) then
        zbottomofcell = -min( z_edges_in(k+1), g%bathyT(cs%col_i(col),cs%col_j(col)) )
!        tmpT1d(k) = sp_val(CS%col_i(col),CS%col_j(col),k)
      elseif (k>1) then
        zbottomofcell = -g%bathyT(cs%col_i(col),cs%col_j(col))
!        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(col),cs%col_j(col)) )
    cs%Ref_val(cs%fldno)%h(1:nz_data,col) = 0.
    cs%Ref_val(cs%fldno)%p(1:nz_data,col) = -1.e24
    cs%Ref_val(cs%fldno)%h(1:nz_data,col) = gv%Z_to_H*hsrc(1:nz_data)
!    CS%Ref_val(CS%fldno)%p(1:nz_data,col) = tmpT1d(1:nz_data)
  enddo
 
  cs%var(cs%fldno)%p => f_ptr
 
  deallocate( hsrc )
  deallocate( tmpt1d )
  deallocate(sp_val, mask_z)
 
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)
 
!!! 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)
 
  ! 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_(G)), &
                             intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2] (in)
  real,                      intent(in)    :: dt !< The amount of time covered by this call [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_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 :: idt                                   ! 1.0/dt [s-1].
  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
  integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz, nz_data
  real, allocatable, dimension(:), target :: z_in, z_edges_in
  real :: missing_value
  real :: h_neglect, h_neglect_edge
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  if (.not.associated(cs)) return
 
  if (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%new_sponges) then
    if (.not. present(time)) &
         call mom_error(fatal,"apply_ALE_sponge: No time information provided")
 
! Interpolate new grid in time-space
    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(cs%fldno)%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)
 
!      call pass_var(sp_val,G%Domain)
!      call pass_var(mask_z,G%Domain)
 
 
      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)
        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)
    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%new_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%new_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)
 
!      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)
 
!      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%new_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%new_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%new_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
 
! 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_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(cs%fldno)%p)) deallocate(cs%Ref_val(cs%fldno)%p)
  enddo
 
  deallocate(cs)
 
end subroutine ale_sponge_end
end module mom_ale_sponge