plm_edges_layer Subroutine

   pure subroutine plm_edges_layer(k, nz, h_dn, h_c, h_up, q_dn, q_c, q_up, q_t, q_b)
      !$acc routine seq
      !! PLM top/bottom edge values of ONE layer `k` of a layer-mean field
      !! `q`, via a two-stage h-weighted van-Leer slope (White, Adcroft &
      !! Hallberg 2009 §2).  Returns the SHALLOWER edge in `q_t` (toward
      !! k+1) and the DEEPER edge in `q_b` (toward k-1), bottom-up.
      !! Boundary layers (k=1, k=nz) -> `boundary_edges_linear`, the
      !! linear-exact one-sided pair.
      !!
      !! Every layer's edges depend on its two neighbours only, so the
      !! caller runs this one thread per CELL (the FV-MOM6 reconstruct
      !! kernel's Pass 0 is a 3-D `do concurrent`); the neighbour arguments
      !! of a boundary layer that has none 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_dn, h_c, h_up
         !! Thicknesses (m) of layers k-1 (deeper), k, k+1 (shallower).
      real(wp), intent(in)  :: q_dn, q_c, q_up
         !! Layer means of layers k-1, k, k+1.
      real(wp), intent(out) :: q_t
         !! Top (shallower) edge value of layer k.
      real(wp), intent(out) :: q_b
         !! Bottom (deeper) edge value of layer k.

      real(wp) :: slp, sig_c, sig_l, sig_r, slp_max, e_t, e_b, q_lo, q_hi

      ! Single-layer column: PCM is the only option (no neighbour).
      if (nz <= 1) then
         q_t = q_c
         q_b = q_c
         return
      end if
      if (k == 1) then
         call boundary_edges_linear(h_c, h_up, q_c, q_up - q_c, q_t, q_b)
         return
      end if
      if (k == nz) then
         call boundary_edges_linear(h_c, h_dn, q_c, q_c - q_dn, q_t, q_b)
         return
      end if

      ! ---- Stage 1: h-weighted limited central slope ----
      ! sig_c is the change ACROSS the layer measured deeper->shallower:
      ! positive sig_c means q increases toward the surface (k+1).
      ! h-weighted central slope (van-Leer / White-Adcroft-Hallberg):
      !   sig_c = (q(k+1)-q(k-1)) * h(k) / (h(k-1)+2 h(k)+h(k+1))  * 2
      ! then limited to 2*min(|q(k)-q_deeper|,|q_shallower-q(k)|), zeroed
      ! at extrema.
      sig_l = q_c - q_dn     ! deeper one-sided (toward k-1)
      sig_r = q_up - q_c     ! shallower one-sided (toward k+1)
      if (sig_l*sig_r <= 0.0_wp) then
         slp = 0.0_wp        ! local extremum -> flatten
      else
         sig_c = 2.0_wp*(q_up - q_dn)*h_c/(h_dn + 2.0_wp*h_c + h_up)
         slp_max = 2.0_wp*min(abs(sig_l), abs(sig_r))
         slp = sign(min(abs(sig_c), slp_max), sig_c)
      end if

      ! ---- Stage 2: monotonized edges bounded against neighbour means ----
      ! Clamp each edge between the cell mean and the adjacent cell mean
      ! (White, Adcroft & Hallberg 2009 §2 monotonization — prevents the
      ! reconstructed edge from over/undershooting the neighbour mean,
      ! which would manufacture a density inversion under the EOS).
      e_t = q_c + 0.5_wp*slp   ! shallower edge (toward k+1)
      e_b = q_c - 0.5_wp*slp   ! deeper edge (toward k-1)
      q_lo = min(q_c, q_up)
      q_hi = max(q_c, q_up)
      q_t = max(q_lo, min(q_hi, e_t))
      q_lo = min(q_c, q_dn)
      q_hi = max(q_c, q_dn)
      q_b = max(q_lo, min(q_hi, e_b))
   end subroutine plm_edges_layer