ppm_interface_values Subroutine

   pure subroutine ppm_interface_values(m, nz, h0, h1, h2, h3, s0, s1, s2, s3, &
                                        t0, t1, t2, t3, edge_s, edge_t)
      !$acc routine seq
      !! The PPM estimates of salinity and temperature at the interface
      !! between layer `m` (deeper) and `m+1` (shallower), `1 <= m <= nz-1`,
      !! from the four-layer stencil `m-1 .. m+2` (thicknesses `h0..h3`,
      !! layer means `s0..s3`, `t0..t3`).  Interior interfaces
      !! (`2 <= m <= nz-2`) take the implicit-h4 estimate in its explicit
      !! form (White & Adcroft 2008 non-uniform stencil, exactly 4th-order
      !! on non-uniform layers; matches `remap_column_ppm_h4` Step 1 to
      !! round-off).  The two near-boundary interfaces (`m = 1`,
      !! `m = nz-1`) take the thickness-weighted (h2) estimate — the
      !! linear-exact value at the shared face of two piecewise-linear
      !! cells, `(q_m h_{m+1} + q_{m+1} h_m)/(h_m + h_{m+1})` (a plain mean
      !! biases the thick interior layer's edge on non-uniform
      !! thicknesses); layer `m-1` (resp. `m+2`) is then not referenced.
      !!
      !! Both tracers at once: the stencil's thickness factors (five of the
      !! six divisions) are the same for S and T, so they are formed once.
      !! Each tracer's value is the single-tracer expression, operation for
      !! operation.
      integer, intent(in) :: m
         !! Interface index: between layers m and m+1.
      integer, intent(in) :: nz
         !! Number of layers in the column.
      real(wp), intent(in) :: h0, h1, h2, h3
         !! Thicknesses of layers m-1, m, m+1, m+2 (m).
      real(wp), intent(in) :: s0, s1, s2, s3
         !! Salinity layer means of layers m-1 .. m+2.
      real(wp), intent(in) :: t0, t1, t2, t3
         !! Temperature layer means of layers m-1 .. m+2.
      real(wp), intent(out) :: edge_s, edge_t
         !! Interface salinity / temperature.

      real(wp) :: g0, g1, g2, g3, hf, h_sum
      real(wp) :: h01, h12, h23, h012, h123, h0123
      real(wp) :: f1, f3, w2, w3
      real(wp), parameter :: H_NEGLECT = 1.0e-30_wp
      real(wp), parameter :: H_MIN_FRAC = 1.0e-5_wp

      if (m == 1 .or. m == nz - 1) then
         edge_s = (s1*h2 + s2*h1)/(h1 + h2)
         edge_t = (t1*h2 + t2*h1)/(h1 + h2)
         return
      end if
      g0 = h0
      g1 = h1
      g2 = h2
      g3 = h3
      h_sum = g0 + g1 + g2 + g3
      if (g0 + g1 <= 0.0_wp .or. g1 + g2 <= 0.0_wp .or. g2 + g3 <= 0.0_wp) then
         hf = H_MIN_FRAC*max(H_NEGLECT, h_sum)
         g0 = max(g0, hf)
         g1 = max(g1, hf)
         g2 = max(g2, hf)
         g3 = max(g3, hf)
      end if
      h01 = g0 + g1
      h12 = g1 + g2
      h23 = g2 + g3
      h012 = g0 + g1 + g2
      h123 = g1 + g2 + g3
      h0123 = g0 + g1 + g2 + g3
      f1 = h01*h23/h12
      f3 = 1.0_wp/h012 + 1.0_wp/h123
      w2 = g2*h23/(h012*h01)
      w3 = g1*h01/(h123*h23)
      edge_s = (f1*(g2*s1 + g1*s2)*f3 + w2*((g0 + 2.0_wp*g1)*s1 - g1*s0) &
                + w3*((2.0_wp*g2 + g3)*s2 - g2*s3))/h0123
      edge_t = (f1*(g2*t1 + g1*t2)*f3 + w2*((g0 + 2.0_wp*g1)*t1 - g1*t0) &
                + w3*((2.0_wp*g2 + g3)*t2 - g2*t3))/h0123
   end subroutine ppm_interface_values