diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90
index 9f8b2336f9..47ebf30bfb 100644
--- a/src/ALE/MOM_ALE.F90
+++ b/src/ALE/MOM_ALE.F90
@@ -129,6 +129,8 @@ module MOM_ALE
 public ALE_remap_scalar
 public ALE_remap_tracers
 public ALE_remap_velocities
+public ALE_remap_interface_vals
+public ALE_remap_vertex_vals
 public ALE_PLM_edge_values
 public TS_PLM_edge_values
 public TS_PPM_edge_values
@@ -974,6 +976,89 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old, h_new, u, v, OBC, dzInterface,
 
 end subroutine ALE_remap_velocities
 
+!> Interpolate to find an updated array of values at interfaces after remapping.
+subroutine ALE_remap_interface_vals(CS, G, GV, h_old, h_new, int_val)
+  type(ALE_CS),                              intent(in)    :: CS       !< ALE control structure
+  type(ocean_grid_type),                     intent(in)    :: G        !< Ocean grid structure
+  type(verticalGrid_type),                   intent(in)    :: GV       !< Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_old    !< Thickness of source grid
+                                                                       !! [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_new    !< Thickness of destination grid
+                                                                       !! [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                                             intent(inout) :: int_val  !< The interface values to interpolate [A]
+
+  real :: val_src(GV%ke)  ! A column of interface values on the source grid [A]
+  real :: val_tgt(GV%ke)  ! A column of interface values on the target grid [A]
+  real :: h_src(GV%ke)    ! A column of source grid layer thicknesses [H ~> m or kg m-2]
+  real :: h_tgt(GV%ke)    ! A column of target grid layer thicknesses [H ~> m or kg m-2]
+  real :: h_neglect, h_neglect_edge  ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
+  logical :: show_call_tree
+  integer :: i, j, k, nz
+
+  do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j)>0.) then
+    do k=1,nz
+      h_src(k) = h_old(i,j,k)
+      h_tgt(k) = h_new(i,j,k)
+    enddo
+
+    do K=1,nz+1
+      val_src(K) = int_val(i,j,K)
+    enddo
+
+    call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
+
+    do K=1,nz+1
+      int_val(i,j,K) = val_tgt(K)
+    enddo
+  endif ; enddo ; enddo
+
+end subroutine ALE_remap_interface_vals
+
+!> Interpolate to find an updated array of values at vertices of tracer cells after remapping.
+subroutine ALE_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val)
+  type(ALE_CS),                              intent(in)    :: CS       !< ALE control structure
+  type(ocean_grid_type),                     intent(in)    :: G        !< Ocean grid structure
+  type(verticalGrid_type),                   intent(in)    :: GV       !< Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_old    !< Thickness of source grid
+                                                                       !! [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_new    !< Thickness of destination grid
+                                                                       !! [H ~> m or kg m-2]
+  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
+                                             intent(inout) :: vert_val  !< The interface values to interpolate [A]
+
+  real :: val_src(GV%ke)  ! A column of interface values on the source grid [A]
+  real :: val_tgt(GV%ke)  ! A column of interface values on the target grid [A]
+  real :: h_src(GV%ke)    ! A column of source grid layer thicknesses [H ~> m or kg m-2]
+  real :: h_tgt(GV%ke)    ! A column of target grid layer thicknesses [H ~> m or kg m-2]
+  real :: I_mask_sum      ! The inverse of the tracer point masks surrounding a corner [nondim]
+  real :: h_neglect, h_neglect_edge  ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
+  logical :: show_call_tree
+  integer :: i, j, k, nz
+
+  do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB
+    if ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) > 0.0 ) then
+      I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)))
+
+    do k=1,nz
+      h_src(k) = ((G%mask2dT(i,j) * h_old(i,j,k) + G%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + &
+                  (G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum
+      h_tgt(k) = ((G%mask2dT(i,j) * h_new(i,j,k) + G%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + &
+                  (G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum
+    enddo
+
+    do K=1,nz+1
+      val_src(K) = vert_val(I,J,K)
+    enddo
+
+    call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
+
+    do K=1,nz+1
+      vert_val(I,J,K) = val_tgt(K)
+    enddo
+  endif ; enddo ; enddo
+
+end subroutine ALE_remap_vertex_vals
 
 !> Mask out thicknesses to 0 when their running sum exceeds a specified value.
 subroutine apply_partial_cell_mask(h1, h_mask)