remap_column_ppm_h4 Subroutine

   pure subroutine remap_column_ppm_h4(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)
      !$acc routine seq
      !! PPM with non-uniform 4th-order (H4) edge values (White & Adcroft 2008).
      !! As `remap_column_ppm` but the interior edge estimate is the
      !! thickness-weighted exactly-4th-order stencil (reduces to PPM's
      !! (7/12,-1/12) on uniform layers), cutting spurious diapycnal mixing per
      !! remap (Ilicak et al. 2012). Limiter/reconstruction/integration are
      !! identical to PPM. H4 edge algebra is inlined (helper extraction costs
      !! 4-6% on this hot kernel). Boundary edges: outermost = PCM,
      !! second-from-boundary = non-uniform 3-cell (H3) quadratic; CW limiter
      !! clamps all edges to local monotone bounds.
      integer, intent(in) :: nz
      real(wp), intent(in) :: dz_old(nz)
         !! Old layer thicknesses
      real(wp), intent(in) :: dz_new(nz)
         !! New layer thicknesses
      real(wp), intent(in) :: q_old(nz)
         !! Old cell-average scalar values
      real(wp), intent(out) :: q_new(nz)
         !! New cell-average scalar values (conservative)

      real(wp) :: z_old(0:NZ_STACK_MAX), z_new(0:NZ_STACK_MAX)
      real(wp) :: q_L(NZ_STACK_MAX), q_R(NZ_STACK_MAX), q6(NZ_STACK_MAX)
      real(wp) :: z_lo, z_hi, overlap, integral
      real(wp) :: xi_lo, xi_hi
      logical, intent(in), optional :: bnd_extrap
         !! Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
      logical, intent(in), optional :: nonunif
         !! Non-uniform-grid weights.  The H4/H3 stencils below are ALREADY
         !! thickness-weighted, so this only reaches the `nz == 2` PLM
         !! fallback; passed through for consistency.

      real(wp) :: dq, dq_l, dq_r, q_min, q_max
      real(wp) :: d_bnd
      logical :: be, nu
      real(wp) :: h0, h1, h2, h3, hf, h_sum
      real(wp) :: h01, h12, h23, h012, h123, h0123
      real(wp) :: f1, f2, f3, et1, et2, et3
      real(wp) :: m11, m12, m13, m21, m22, m23, m31, m32, m33, det
      real(wp) :: z1, z2, z3, ca, cb, cc, edge_val
      integer :: k, ko, ko_start

      be = .false.
      if (present(bnd_extrap)) be = bnd_extrap
      nu = .false.
      if (present(nonunif)) nu = nonunif

      ! Trivial cases (identical to PPM)
      if (nz == 1) then
         q_new(1) = q_old(1)
         return
      end if
      if (nz == 2) then
         ! With only 2 layers, PPM reduces to PLM
         call remap_column_plm(nz, dz_old, dz_new, q_old, q_new, be, nu)
         return
      end if

      ! Build interface positions
      z_old(0) = 0.0_wp
      z_new(0) = 0.0_wp
      do k = 1, nz
         z_old(k) = z_old(k - 1) + dz_old(k)
         z_new(k) = z_new(k - 1) + dz_new(k)
      end do

      ! ---- Step 1 (PPM_H4): non-uniform 4th-order edge values ----
      ! Store edges in q_R(k) = right edge of layer k = q_L(k+1).

      ! Outermost edges: PCM (layer 1 left, layer nz right)
      q_L(1) = q_old(1)
      q_R(nz) = q_old(nz)

      ! Second-from-boundary edges: non-uniform 3-cell H3 quadratic.
      ! Bed side: edge between layer 1 and 2 — fit a quadratic through layers
      ! 1,2,3 (cell-average constraints) and evaluate at the 1|2 interface.
      ! z-coordinate with the left of layer 1 at 0:
      h1 = dz_old(1)
      h2 = dz_old(2)
      h3 = dz_old(3)
      ! H3 boundary singularity guard (mirror of the interior floor)
      if (h1 <= 0.0_wp .or. h2 <= 0.0_wp .or. h3 <= 0.0_wp) then
         hf = H_MIN_FRAC*max(H_NEGLECT, dz_old(1) + dz_old(2) + dz_old(3))
         h1 = max(h1, hf)
         h2 = max(h2, hf)
         h3 = max(h3, hf)
      end if
      z1 = h1
      z2 = h1 + h2
      z3 = h1 + h2 + h3
      ! Cell moments [int z^2 / h, int z / h, 1] over each layer:
      m11 = (z1*z1*z1)/(3.0_wp*h1)
      m12 = (z1*z1)/(2.0_wp*h1)
      m13 = 1.0_wp
      m21 = (z2*z2*z2 - z1*z1*z1)/(3.0_wp*h2)
      m22 = (z2*z2 - z1*z1)/(2.0_wp*h2)
      m23 = 1.0_wp
      m31 = (z3*z3*z3 - z2*z2*z2)/(3.0_wp*h3)
      m32 = (z3*z3 - z2*z2)/(2.0_wp*h3)
      m33 = 1.0_wp
      det = m11*(m22*m33 - m23*m32) - m12*(m21*m33 - m23*m31) + m13*(m21*m32 - m22*m31)
      ! Cramer's rule for the quadratic coefficients ca*z^2 + cb*z + cc:
      ca = (q_old(1)*(m22*m33 - m23*m32) - m12*(q_old(2)*m33 - m23*q_old(3)) &
            + m13*(q_old(2)*m32 - m22*q_old(3)))/det
      cb = (m11*(q_old(2)*m33 - m23*q_old(3)) - q_old(1)*(m21*m33 - m23*m31) &
            + m13*(m21*q_old(3) - q_old(2)*m31))/det
      cc = (m11*(m22*q_old(3) - q_old(2)*m32) - m12*(m21*q_old(3) - q_old(2)*m31) &
            + q_old(1)*(m21*m32 - m22*m31))/det
      edge_val = ca*z1*z1 + cb*z1 + cc   ! evaluate at 1|2 interface
      q_R(1) = edge_val
      q_L(2) = edge_val

      ! Surface side: edge between layers nz-1 and nz — quadratic through
      ! layers nz-2, nz-1, nz, evaluated at the (nz-1)|nz interface.
      h1 = dz_old(nz - 2)
      h2 = dz_old(nz - 1)
      h3 = dz_old(nz)
      ! H3 boundary singularity guard (mirror of the interior floor)
      if (h1 <= 0.0_wp .or. h2 <= 0.0_wp .or. h3 <= 0.0_wp) then
         hf = H_MIN_FRAC*max(H_NEGLECT, dz_old(nz - 2) + dz_old(nz - 1) + dz_old(nz))
         h1 = max(h1, hf)
         h2 = max(h2, hf)
         h3 = max(h3, hf)
      end if
      z1 = h1
      z2 = h1 + h2
      z3 = h1 + h2 + h3
      m11 = (z1*z1*z1)/(3.0_wp*h1)
      m12 = (z1*z1)/(2.0_wp*h1)
      m13 = 1.0_wp
      m21 = (z2*z2*z2 - z1*z1*z1)/(3.0_wp*h2)
      m22 = (z2*z2 - z1*z1)/(2.0_wp*h2)
      m23 = 1.0_wp
      m31 = (z3*z3*z3 - z2*z2*z2)/(3.0_wp*h3)
      m32 = (z3*z3 - z2*z2)/(2.0_wp*h3)
      m33 = 1.0_wp
      det = m11*(m22*m33 - m23*m32) - m12*(m21*m33 - m23*m31) + m13*(m21*m32 - m22*m31)
      ca = (q_old(nz - 2)*(m22*m33 - m23*m32) - m12*(q_old(nz - 1)*m33 - m23*q_old(nz)) &
            + m13*(q_old(nz - 1)*m32 - m22*q_old(nz)))/det
      cb = (m11*(q_old(nz - 1)*m33 - m23*q_old(nz)) - q_old(nz - 2)*(m21*m33 - m23*m31) &
            + m13*(m21*q_old(nz) - q_old(nz - 1)*m31))/det
      cc = (m11*(m22*q_old(nz) - q_old(nz - 1)*m32) - m12*(m21*q_old(nz) - q_old(nz - 1)*m31) &
            + q_old(nz - 2)*(m21*m32 - m22*m31))/det
      edge_val = ca*z2*z2 + cb*z2 + cc   ! evaluate at (nz-1)|nz interface
      q_R(nz - 1) = edge_val
      q_L(nz) = edge_val

      ! Interior edges: full non-uniform H4 stencil (4 cells k-1..k+2).
      ! q_R(k) is the edge between layer k and k+1; stencil thicknesses
      ! h0=dz(k-1), h1=dz(k), h2=dz(k+1), h3=dz(k+2).
      do k = 2, nz - 2
         h0 = dz_old(k - 1)
         h1 = dz_old(k)
         h2 = dz_old(k + 1)
         h3 = dz_old(k + 2)

         ! Conditional singularity guard: only floor when a consecutive
         ! thickness pair sums to ~0 (vanishing layers under ZSTAR_FULL).
         h_sum = h0 + h1 + h2 + h3
         if (h0 + h1 <= 0.0_wp .or. h1 + h2 <= 0.0_wp .or. h2 + h3 <= 0.0_wp) then
            hf = H_MIN_FRAC*max(H_NEGLECT, h_sum)
            h0 = max(h0, hf)
            h1 = max(h1, hf)
            h2 = max(h2, hf)
            h3 = max(h3, hf)
         end if

         h01 = h0 + h1
         h12 = h1 + h2
         h23 = h2 + h3
         h012 = h0 + h1 + h2
         h123 = h1 + h2 + h3
         h0123 = h0 + h1 + h2 + h3

         f1 = h01*h23/h12
         f2 = h2*q_old(k) + h1*q_old(k + 1)
         f3 = 1.0_wp/h012 + 1.0_wp/h123
         et1 = f1*f2*f3
         et2 = (h2*h23/(h012*h01))*((h0 + 2.0_wp*h1)*q_old(k) - h1*q_old(k - 1))
         et3 = (h1*h01/(h123*h23))*((2.0_wp*h2 + h3)*q_old(k + 1) - h2*q_old(k + 2))

         q_R(k) = (et1 + et2 + et3)/h0123
         q_L(k + 1) = q_R(k)
      end do

      ! ---- Step 2: Colella-Woodward monotonicity limiting (verbatim) ----
      do k = 1, nz
         q_min = q_old(k)
         q_max = q_old(k)
         if (k > 1) then
            q_min = min(q_min, q_old(k - 1))
            q_max = max(q_max, q_old(k - 1))
         end if
         if (k < nz) then
            q_min = min(q_min, q_old(k + 1))
            q_max = max(q_max, q_old(k + 1))
         end if

         ! Clip edges to local bounds
         q_L(k) = max(q_min, min(q_max, q_L(k)))
         q_R(k) = max(q_min, min(q_max, q_R(k)))

         ! CW monotonicity: if the cell is a local extremum, flatten
         dq = q_R(k) - q_L(k)
         dq_l = q_old(k) - q_L(k)
         dq_r = q_R(k) - q_old(k)
         if (dq_l*dq_r <= 0.0_wp) then
            ! Local extremum: flatten to PCM
            q_L(k) = q_old(k)
            q_R(k) = q_old(k)
         else
            ! Check if parabola overshoots
            q6(k) = 6.0_wp*q_old(k) - 3.0_wp*(q_L(k) + q_R(k))
            if (abs(q6(k)) > abs(dq)) then
               if (q6(k)*dq > 0.0_wp) then
                  ! Overshoot near left edge: adjust q_L
                  q_L(k) = 3.0_wp*q_old(k) - 2.0_wp*q_R(k)
               else
                  ! Overshoot near right edge: adjust q_R
                  q_R(k) = 3.0_wp*q_old(k) - 2.0_wp*q_L(k)
               end if
            end if
         end if

         ! Recompute q6 after limiting
         q6(k) = 6.0_wp*q_old(k) - 3.0_wp*(q_L(k) + q_R(k))
      end do

      ! ---- Step 2b: boundary-cell closure (opt-in; see remap_column_ppm) ----
      if (be) then
         call boundary_half_jump(dz_old(1), dz_old(2), q_old(2) - q_old(1), d_bnd)
         q_L(1) = q_old(1) - d_bnd
         q_R(1) = q_old(1) + d_bnd
         q6(1) = 0.0_wp
         call boundary_half_jump(dz_old(nz), dz_old(nz - 1), &
                                 q_old(nz) - q_old(nz - 1), d_bnd)
         q_L(nz) = q_old(nz) - d_bnd
         q_R(nz) = q_old(nz) + d_bnd
         q6(nz) = 0.0_wp
      end if

      ! ---- Step 3: Integrate parabolic reconstruction over new layers (verbatim) ----
      ko_start = 1
      do k = 1, nz
         if (dz_new(k) <= 0.0_wp) then
            q_new(k) = 0.0_wp
            cycle
         end if

         integral = 0.0_wp
         do ko = ko_start, nz
            z_lo = max(z_new(k - 1), z_old(ko - 1))
            z_hi = min(z_new(k), z_old(ko))
            overlap = z_hi - z_lo

            if (overlap <= 0.0_wp) then
               if (z_old(ko) > z_new(k)) exit
               cycle
            end if

            if (dz_old(ko) > 0.0_wp) then
               xi_lo = (z_lo - z_old(ko - 1))/dz_old(ko)
               xi_hi = (z_hi - z_old(ko - 1))/dz_old(ko)

               ! Parabolic integral
               integral = integral + dz_old(ko)*( &
                          (xi_hi - xi_lo)*q_L(ko) &
                          + 0.5_wp*(xi_hi*xi_hi - xi_lo*xi_lo)*(q_R(ko) - q_L(ko) + q6(ko)) &
                          - (xi_hi*xi_hi*xi_hi - xi_lo*xi_lo*xi_lo)*q6(ko)/3.0_wp)
            else
               integral = integral + q_old(ko)*overlap
            end if

            if (z_old(ko) <= z_new(k)) ko_start = ko
         end do

         q_new(k) = integral/dz_new(k)
      end do
   end subroutine remap_column_ppm_h4