ppm_edges_layer Subroutine

   pure subroutine ppm_edges_layer(k, nz, h_m2, h_m1, h_c, h_p1, h_p2, &
                                   s_m2, s_m1, s_c, s_p1, s_p2, &
                                   t_m2, t_m1, t_c, t_p1, t_p2, &
                                   s_top, s_bot, t_top, t_bot)
      !$acc routine seq
      !! PPM top/bottom edge values of salinity and temperature in ONE layer
      !! `k`: the implicit-h4 interface estimates (`ppm_interface_values`)
      !! + the Colella & Woodward (1984) limiter (`ppm_limit_edges`).
      !! Returns the SHALLOWER edges in `*_top`, the DEEPER in `*_bot`
      !! (bottom-up).
      !!
      !! NOTE the PGF integrand on top of these edges is PARABOLIC (see
      !! `boole_dpa_intz_layer`): q6 = 3*(2*q_mean - (q_t + q_b)) is the
      !! in-layer curvature — why PPM differs from PLM at the density
      !! integral even at identical edge values.
      !!
      !! Boundary layers (k=1, k=nz) -> `boundary_edges_linear`.  The pair
      !! is symmetric about the layer mean, so `q6 = 3*(2*q - (q_t+q_b))`
      !! is identically zero there: the boundary layer carries a straight
      !! line, which is the exact profile whenever `q(z)` is linear.
      !!
      !! A layer's edges depend on the five layers `k-2 .. k+2` only, so the
      !! caller runs this one thread per CELL (the FV-MOM6 reconstruct
      !! kernel's Pass 0 is a 3-D `do concurrent`); arguments of layers
      !! outside the column are not referenced.
      integer, intent(in) :: k
         !! Layer index (1 = bed, nz = surface).
      integer, intent(in) :: nz
         !! Number of layers in the column.
      real(wp), intent(in)  :: h_m2, h_m1, h_c, h_p1, h_p2
         !! Thicknesses (m) of layers k-2 .. k+2.
      real(wp), intent(in)  :: s_m2, s_m1, s_c, s_p1, s_p2
         !! Salinity layer means of layers k-2 .. k+2.
      real(wp), intent(in)  :: t_m2, t_m1, t_c, t_p1, t_p2
         !! Temperature layer means of layers k-2 .. k+2.
      real(wp), intent(out) :: s_top, s_bot, t_top, t_bot
         !! Top (shallower) / bottom (deeper) edge values of layer k.

      real(wp) :: sl, sr, tl, tr

      if (nz <= 1) then
         s_top = s_c
         s_bot = s_c
         t_top = t_c
         t_bot = t_c
         return
      end if
      ! Boundary layers: linear-exact one-sided pair (q6 == 0 there).
      if (k == 1) then
         call boundary_edges_linear(h_c, h_p1, s_c, s_p1 - s_c, s_top, s_bot)
         call boundary_edges_linear(h_c, h_p1, t_c, t_p1 - t_c, t_top, t_bot)
         return
      end if
      if (k == nz) then
         call boundary_edges_linear(h_c, h_m1, s_c, s_c - s_m1, s_top, s_bot)
         call boundary_edges_linear(h_c, h_m1, t_c, t_c - t_m1, t_top, t_bot)
         return
      end if

      call ppm_interface_values(k - 1, nz, h_m2, h_m1, h_c, h_p1, s_m2, s_m1, s_c, s_p1, &
                                t_m2, t_m1, t_c, t_p1, sl, tl)
      call ppm_interface_values(k, nz, h_m1, h_c, h_p1, h_p2, s_m1, s_c, s_p1, s_p2, &
                                t_m1, t_c, t_p1, t_p2, sr, tr)
      call ppm_limit_edges(s_m1, s_c, s_p1, sl, sr, s_top, s_bot)
      call ppm_limit_edges(t_m1, t_c, t_p1, tl, tr, t_top, t_bot)
   end subroutine ppm_edges_layer