+Add mask_edges argument to interpolate_column

  Added a new logical argument to interpolate_column to specify whether the
interpolated interface values outside of massless layers should be masked to
zero.  Also refactored the code in interpolate_column to separate out the
determination of the interface position from the interpolation and the masking
to facilitate the extension of this code to use higher order interpolation in
planned subsequent changes.  All answers are bitwise identical, although there
is a new mandatory argument for a public interface.

src/ALE/MOM_ALE.F90

@@ -554,7 +554,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
       if (check_column_integrals(nk, h_src, nk, h_dest)) then
         call MOM_error(FATAL, "ALE_offline_inputs: Kd interpolation columns do not match")
       endif
-      call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), Kd(i,j,:))
+      call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), Kd(i,j,:), .true.)
     endif
   enddo ; enddo
 

src/ALE/MOM_remapping.F90

@@ -840,78 +840,72 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth
 end subroutine remap_via_sub_cells
 
 !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest
-subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest)
+subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, mask_edges)
   integer,                  intent(in)    :: nsrc   !< Number of source cells
   real, dimension(nsrc),    intent(in)    :: h_src  !< Thickness of source cells [H]
   real, dimension(nsrc+1),  intent(in)    :: u_src  !< Values at source cell interfaces [A]
   integer,                  intent(in)    :: ndest  !< Number of destination cells
   real, dimension(ndest),   intent(in)    :: h_dest !< Thickness of destination cells [H]
   real, dimension(ndest+1), intent(inout) :: u_dest !< Interpolated value at destination cell interfaces [A]
+  logical,                  intent(in)    :: mask_edges !< If true, mask the values outside of massless
+                                                    !! layers at the top and bottom of the column.
 
   ! Local variables
-  real :: x_dest ! Relative position of target interface [H]
-  real :: dh ! Source cell thickness [H]
-  real :: u1, u2 ! Values to interpolate between [A]
-  real :: weight_a, weight_b ! Weights for interpolation [nondim]
-  integer :: k_src, k_dest ! Index of cell in src and dest columns
-  logical :: still_vanished ! Used for figuring out what to mask as missing
-
-  ! Initial values for the loop
-  still_vanished = .true.
+  real :: x_dest            ! Relative position of target interface [H]
+  real :: dh                ! Source cell thickness [H]
+  real :: frac_pos(ndest+1) ! Fractional position of the destination interface
+                            ! within the source layer [nondim], 0 <= frac_pos <= 1.
+  integer :: k_src(ndest+1) ! Source grid layer index of destination interface, 1 <= k_src <= ndest.
+  integer :: ks, k_dest     ! Index of cell in src and dest columns
 
   ! The following forces the "do while" loop to do one cycle that will set u1, u2, dh.
-  k_src = 0
+  ks = 0
   dh = 0.
   x_dest = 0.
 
-  do k_dest=1, ndest+1
-    do while (dh<=x_dest .and. k_src<nsrc)
+  ! Find the layer index and fractional position of the interfaces of the target
+  ! grid on the source grid.
+  do k_dest=1,ndest+1
+    do while (dh<=x_dest .and. ks<nsrc)
       ! Move positions pointers forward until the interval 0 .. dh spans x_dest.
       x_dest = x_dest - dh
-      k_src = k_src + 1
-      dh = h_src(k_src) ! Thickness of layer k_src
-
-      ! Values that will be used for the interpolation.
-      u1 = u_src(k_src) ! Value on left of source cell
-      u2 = u_src(k_src+1) ! Value on right of source cell
+      ks = ks + 1
+      dh = h_src(ks) ! Thickness of layer ks
     enddo
+    k_src(k_dest) = ks
 
     if (dh>0.) then
-      weight_b = max(0., min(1., x_dest / dh)) ! Weight of u2
+      frac_pos(k_dest) = max(0., min(1., x_dest / dh)) ! Weight of u2
     else  ! For a vanished source layer we need to do something reasonable...
-      weight_b = 0.5
+      frac_pos(k_dest) = 0.5
     endif
-    weight_a = 1.0 - weight_b   ! Weight of u1
-    ! Linear interpolation between u1 and u2
-    u_dest(k_dest) = weight_a * u1 + weight_b * u2
 
-    ! Mask vanished layers at the surface which would be under an ice-shelf.
-    ! TODO: Need to figure out what to do for an isopycnal coordinate diagnostic that could
-    !       also have vanished layers at the surface.
-    if (k_dest<=ndest) then
+    if (k_dest <= ndest) then
       x_dest = x_dest + h_dest(k_dest) ! Position of interface k_dest+1
-      if (still_vanished .and. h_dest(k_dest)==0.) then
-        ! When the layer k_dest is vanished and all layers above are also vanished, the k_dest
-        ! interface value should be missing.
-        u_dest(k_dest) = 0.0
-      else
-        still_vanished = .false.
-      endif
     endif
+  enddo
 
+  do k_dest=1,ndest+1
+    ! Linear interpolation between surrounding edge values.
+    ks = k_src(k_dest)
+    u_dest(k_dest) = (1.0 - frac_pos(k_dest)) * u_src(ks) + frac_pos(k_dest) * u_src(ks+1)
   enddo
 
-  ! Mask vanished layers on topography
-  still_vanished = .true.
-  do k_dest=ndest, 1, -1
-    if (still_vanished .and. h_dest(k_dest)==0.) then
-      ! When the layer k_dest is vanished and all layers below are also vanished, the k_dest+1
-      ! interface value should be missing.
+  if (mask_edges) then
+    ! Mask vanished layers at the surface which would be under an ice-shelf.
+    ! When the layer k_dest is vanished and all layers above are also vanished,
+    ! the k_dest interface value should be missing.
+    do k_dest=1,ndest
+      if (h_dest(k_dest) > 0.) exit
+      u_dest(k_dest) = 0.0
+    enddo
+
+    ! Mask interfaces below vanished layers at the bottom
+    do k_dest=ndest,1,-1
+      if (h_dest(k_dest) > 0.) exit
       u_dest(k_dest+1) = 0.0
-    else
-      exit
-    endif
-  enddo
+    enddo
+  endif
 
 end subroutine interpolate_column
 
@@ -1717,7 +1711,7 @@ logical function test_interp(verbose, msg, nsrc, h_src, u_src, ndest, h_dest, u_
   real :: error
 
   ! Interpolate from src to dest
-  call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest)
+  call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, .true.)
 
   test_interp = .false.
   do k=1,ndest+1

src/framework/MOM_diag_remap.F90

@@ -608,7 +608,7 @@ subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, sta
         h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
         h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
         call interpolate_column(nz_src, h_src, field(I1,j,:), &
-                                nz_dest, h_dest, interpolated_field(I1,j,:))
+                                nz_dest, h_dest, interpolated_field(I1,j,:), .true.)
       enddo
     enddo
   elseif (staggered_in_y .and. .not. staggered_in_x) then
@@ -623,7 +623,7 @@ subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, sta
         h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
         h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
         call interpolate_column(nz_src, h_src, field(i,J1,:), &
-                                nz_dest, h_dest, interpolated_field(i,J1,:))
+                                nz_dest, h_dest, interpolated_field(i,J1,:), .true.)
       enddo
     enddo
   elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
@@ -636,7 +636,7 @@ subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, sta
         h_src(:) = h(i,j,:)
         h_dest(:) = remap_cs%h(i,j,:)
         call interpolate_column(nz_src, h_src, field(i,j,:), &
-                                nz_dest, h_dest, interpolated_field(i,j,:))
+                                nz_dest, h_dest, interpolated_field(i,j,:), .true.)
       enddo
     enddo
   else