remap_column_pqm Subroutine

   pure subroutine remap_column_pqm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)
      !$acc routine seq
      !! Piecewise-quartic (PQM_IH4IH3) conservative remap (White & Adcroft 2008).
      !! Implicit-h4 edge VALUES + implicit-h3 edge SLOPES (each a
      !! diagonally-dominant tridiagonal solve with one-sided 4-cell boundary
      !! closure), per-cell quartic, W&A monotonicity limiter, conservative
      !! quartic overlap integral. Cuts diapycnal mixing per remap vs PPM/PPM_H4
      !! (Ilicak et al. 2012).
      !! Per cell k, xi in [0,1]: q_hat(xi) = a + b*xi + c*xi^2 + d*xi^3 + e*xi^4.
      !! Boundary cells reconstruct as PCM. nz < 5 falls back to REMAP_PPM
      !! (W&A boundary closure needs >= 4 cells).
      integer, intent(in) :: nz
      real(wp), intent(in) :: dz_old(nz)
         !! Old layer thicknesses (must sum to same total as dz_new)
      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)

      logical, intent(in), optional :: bnd_extrap
         !! Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
      logical, intent(in), optional :: nonunif
         !! Non-uniform-grid weights.  The implicit-h4/h3 stencils below are
         !! ALREADY thickness-weighted, so this only reaches the `nz < 5`
         !! PPM fallback; passed through for consistency.

      logical :: be, nu
      real(wp) :: d_bnd
      real(wp) :: z_old(0:NZ_STACK_MAX), z_new(0:NZ_STACK_MAX)
      ! Edge values / slopes, two per cell: index 1 = left, 2 = right.
      real(wp) :: ev_l(NZ_STACK_MAX), ev_r(NZ_STACK_MAX)
      real(wp) :: es_l(NZ_STACK_MAX), es_r(NZ_STACK_MAX)
      ! Per-cell quartic coefficients a..e.
      real(wp) :: pa(NZ_STACK_MAX), pb(NZ_STACK_MAX), pc(NZ_STACK_MAX)
      real(wp) :: pd(NZ_STACK_MAX), pe(NZ_STACK_MAX)
      ! Tridiagonal workspace (N+1 edges).
      real(wp) :: tri_l(NZ_STACK_MAX + 1), tri_c(NZ_STACK_MAX + 1)
      real(wp) :: tri_u(NZ_STACK_MAX + 1), tri_b(NZ_STACK_MAX + 1)
      real(wp) :: tri_x(NZ_STACK_MAX + 1)
      real(wp) :: dz4(4), u4(4), csys(4)
      real(wp) :: h0, h1, i_h2, alpha, beta, abmix, aco, bco
      real(wp) :: i_h, h0h1, h0_2, h1_2, h0_3, h1_3, i_d
      real(wp) :: z_lo, z_hi, overlap, integral, xi_lo, xi_hi
      real(wp) :: h_c, u0_l, u0_r, u1_l, u1_r, u_l, u_c, u_r, h_l, h_r
      real(wp) :: sigma_l, sigma_c, sigma_r, slope, slope_x_h
      real(wp) :: u0_avg
      real(wp) :: a, b, cco, dco, eco, alpha1, alpha2, alpha3
      real(wp) :: rho, sqrt_rho, x1, x2, grad1, grad2
      integer :: k, ko, ko_start, np1, inflexion_l, inflexion_r

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

      ! Trivial / degenerate cases — fall back to lower-order safe paths.
      if (nz == 1) then
         q_new(1) = q_old(1)
         return
      end if
      if (nz < 5) then
         call remap_column_ppm(nz, dz_old, dz_new, q_old, q_new, be, nu)
         return
      end if

      np1 = nz + 1

      ! 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: implicit-h4 edge VALUES (tridiagonal, N+1 edges) ----
      ! Interior edges i+1 (between cells i and i+1) — roundoff-safe stencil
      ! al*x(i)+x(i+1)+be*x(i+2) with diagonal OFFSET tri_c = 2*abmix.
      do k = 1, nz - 1
         h0 = max(dz_old(k), H_NEGLECT)
         h1 = max(dz_old(k + 1), H_NEGLECT)
         if (abs(h0) < H_RATIO_FLOOR*abs(h1)) h0 = H_RATIO_FLOOR*h1
         if (abs(h1) < H_RATIO_FLOOR*abs(h0)) h1 = H_RATIO_FLOOR*h0
         i_h2 = 1.0_wp/((h0 + h1)**2)
         alpha = (h1*h1)*i_h2
         beta = (h0*h0)*i_h2
         abmix = (h0*h1)*i_h2
         aco = 2.0_wp*alpha*(alpha + 2.0_wp*beta + 3.0_wp*abmix)
         bco = 2.0_wp*beta*(beta + 2.0_wp*alpha + 3.0_wp*abmix)
         tri_l(k + 1) = alpha
         tri_c(k + 1) = 2.0_wp*abmix
         tri_u(k + 1) = beta
         tri_b(k + 1) = aco*q_old(k) + bco*q_old(k + 1)
      end do
      ! Top boundary (edge 1): exact 4-cell one-sided fit, value = csys(1).
      do k = 1, 4
         dz4(k) = max(H_NEGLECT, dz_old(k))
         u4(k) = q_old(k)
      end do
      call pqm_end_value_h4(dz4, u4, csys)
      tri_b(1) = csys(1)
      tri_c(1) = 1.0_wp
      tri_u(1) = 0.0_wp
      ! Bottom boundary (edge N+1): layers REVERSED, value = csys(1).
      do k = 1, 4
         dz4(k) = max(H_NEGLECT, dz_old(nz + 1 - k))
         u4(k) = q_old(nz + 1 - k)
      end do
      call pqm_end_value_h4(dz4, u4, csys)
      tri_b(np1) = csys(1)
      tri_c(np1) = 1.0_wp
      tri_l(np1) = 0.0_wp

      call pqm_solve_diag_dominant(np1, tri_l, tri_c, tri_u, tri_b, tri_x)

      ! Scatter the N+1 edge values into per-cell (left,right) pairs.
      ev_l(1) = tri_x(1)
      do k = 2, nz
         ev_l(k) = tri_x(k)
         ev_r(k - 1) = tri_x(k)
      end do
      ev_r(nz) = tri_x(np1)

      ! ---- Step 2: implicit-h3 edge SLOPES (tridiagonal, N+1 edges) ----
      ! Nondimensionalised stencil; diagonal OFFSET tri_c (NOT 1.0 — the
      ! diagonal-dominant solver adds tri_l+tri_u to form the pivot).
      do k = 1, nz - 1
         h0 = max(dz_old(k), H_NEGLECT)
         h1 = max(dz_old(k + 1), H_NEGLECT)
         i_h = 1.0_wp/(h0 + h1)
         h0 = h0*i_h
         h1 = h1*i_h
         h0h1 = h0*h1
         h0_2 = h0*h0
         h1_2 = h1*h1
         h0_3 = h0_2*h0
         h1_3 = h1_2*h1
         i_d = 1.0_wp/(4.0_wp*h0h1*(h0 + h1) + h1_3 + h0_3)
         tri_l(k + 1) = (h1*((h0_2 + h0h1) - h1_2))*i_d
         tri_c(k + 1) = 2.0_wp*((h0_2 + h1_2)*(h0 + h1))*i_d
         tri_u(k + 1) = (h0*((h1_2 + h0h1) - h0_2))*i_d
         tri_b(k + 1) = 12.0_wp*(h0h1*i_d)*((q_old(k + 1) - q_old(k))*i_h)
      end do
      ! Top boundary slope = csys(2) of the 4-cell fit.
      do k = 1, 4
         dz4(k) = max(H_NEGLECT, dz_old(k))
         u4(k) = q_old(k)
      end do
      call pqm_end_value_h4(dz4, u4, csys)
      tri_b(1) = csys(2)
      tri_c(1) = 1.0_wp
      tri_u(1) = 0.0_wp
      ! Bottom boundary slope = -csys(2) (layers reversed → sign flip).
      do k = 1, 4
         dz4(k) = max(H_NEGLECT, dz_old(nz + 1 - k))
         u4(k) = q_old(nz + 1 - k)
      end do
      call pqm_end_value_h4(dz4, u4, csys)
      tri_b(np1) = -csys(2)
      tri_c(np1) = 1.0_wp
      tri_l(np1) = 0.0_wp

      call pqm_solve_diag_dominant(np1, tri_l, tri_c, tri_u, tri_b, tri_x)

      es_l(1) = tri_x(1)
      do k = 2, nz
         es_l(k) = tri_x(k)
         es_r(k - 1) = tri_x(k)
      end do
      es_r(nz) = tri_x(np1)

      ! ---- Step 3: PQM limiter (White & Adcroft 2008) ----
      ! 3a. bound_edge_values: van-Leer edge limiting + neighbour-mean clamp.
      do k = 1, nz
         u_l = q_old(max(1, k - 1))
         u_c = q_old(k)
         u_r = q_old(min(k + 1, nz))
         h_l = dz_old(max(1, k - 1))
         h_c = dz_old(k)
         h_r = dz_old(min(k + 1, nz))
         slope_x_h = 0.0_wp
         if (((h_l + h_r) + 2.0_wp*h_c) > 0.0_wp) then
            sigma_l = (u_c - u_l)
            sigma_c = (u_r - u_l)*(h_c/((h_l + h_r) + 2.0_wp*h_c))
            sigma_r = (u_r - u_c)
            if ((sigma_l*sigma_r) > 0.0_wp) then
               slope_x_h = sign(min(abs(sigma_l), abs(sigma_c), abs(sigma_r)), sigma_c)
            end if
         end if
         if ((u_l - ev_l(k))*(ev_l(k) - u_c) < 0.0_wp) then
            ev_l(k) = u_c - sign(min(abs(slope_x_h), abs(ev_l(k) - u_c)), slope_x_h)
         end if
         if ((u_r - ev_r(k))*(ev_r(k) - u_c) < 0.0_wp) then
            ev_r(k) = u_c + sign(min(abs(slope_x_h), abs(ev_r(k) - u_c)), slope_x_h)
         end if
         ev_l(k) = max(min(ev_l(k), max(u_l, u_c)), min(u_l, u_c))
         ev_r(k) = max(min(ev_r(k), max(u_r, u_c)), min(u_r, u_c))
      end do

      ! 3b. check_discontinuous_edge_values: average non-monotonic collocated
      ! edges.  Sweep low→high; ev_km1_r holds the (possibly updated) right
      ! edge of cell k so the pair update stays consistent.
      do k = 1, nz - 1
         if ((ev_l(k + 1) - ev_r(k))*(q_old(k + 1) - q_old(k)) < 0.0_wp) then
            u0_avg = 0.5_wp*(ev_r(k) + ev_l(k + 1))
            u0_avg = max(min(u0_avg, max(q_old(k), q_old(k + 1))), &
                         min(q_old(k), q_old(k + 1)))
            ev_r(k) = u0_avg
            ev_l(k + 1) = u0_avg
         end if
      end do

      ! 3c. interior cells: PLM-slope consistency, extremum flatten, quartic
      ! curvature / inflexion test, collapse + post-collapse resets.
      do k = 2, nz - 1
         inflexion_l = 0
         inflexion_r = 0
         u0_l = ev_l(k)
         u0_r = ev_r(k)
         u1_l = es_l(k)
         u1_r = es_r(k)
         h_l = dz_old(k - 1)
         h_c = dz_old(k)
         h_r = dz_old(k + 1)
         u_l = q_old(k - 1)
         u_c = q_old(k)
         u_r = q_old(k + 1)

         sigma_l = 2.0_wp*(u_c - u_l)/(h_c + H_NEGLECT)
         sigma_c = 2.0_wp*(u_r - u_l)/(h_l + 2.0_wp*h_c + h_r + H_NEGLECT)
         sigma_r = 2.0_wp*(u_r - u_c)/(h_c + H_NEGLECT)
         if ((sigma_l*sigma_r) > 0.0_wp) then
            slope = sign(min(abs(sigma_l), abs(sigma_c), abs(sigma_r)), sigma_c)
         else
            slope = 0.0_wp
         end if

         if (u1_l*slope <= 0.0_wp) u1_l = slope
         if (u1_r*slope <= 0.0_wp) u1_r = slope

         if ((u0_r - u_c)*(u_c - u0_l) <= 0.0_wp) then
            u0_l = u_c
            u0_r = u_c
            u1_l = 0.0_wp
            u1_r = 0.0_wp
            inflexion_l = -1
            inflexion_r = -1
         end if

         if ((inflexion_l == 0) .and. (inflexion_r == 0)) then
            a = u0_l
            b = h_c*u1_l
            cco = 30.0_wp*u_c - 12.0_wp*u0_r - 18.0_wp*u0_l + 1.5_wp*h_c*(u1_r - 3.0_wp*u1_l)
            dco = -60.0_wp*u_c + h_c*(6.0_wp*u1_l - 4.0_wp*u1_r) + 28.0_wp*u0_r + 32.0_wp*u0_l
            eco = 30.0_wp*u_c + 2.5_wp*h_c*(u1_r - u1_l) - 15.0_wp*(u0_l + u0_r)

            alpha1 = 6.0_wp*eco
            alpha2 = 3.0_wp*dco
            alpha3 = cco
            rho = alpha2*alpha2 - 4.0_wp*alpha1*alpha3

            if ((alpha1 /= 0.0_wp) .and. (rho >= 0.0_wp)) then
               sqrt_rho = sqrt(rho)
               x1 = 0.5_wp*(-alpha2 - sqrt_rho)/alpha1
               x2 = 0.5_wp*(-alpha2 + sqrt_rho)/alpha1
               if ((x1 >= 0.0_wp) .and. (x1 <= 1.0_wp) .and. &
                   (x2 >= 0.0_wp) .and. (x2 <= 1.0_wp)) then
                  grad1 = 4.0_wp*eco*(x1**3) + 3.0_wp*dco*(x1**2) + 2.0_wp*cco*x1 + b
                  grad2 = 4.0_wp*eco*(x2**3) + 3.0_wp*dco*(x2**2) + 2.0_wp*cco*x2 + b
                  if ((grad1*slope < 0.0_wp) .or. (grad2*slope < 0.0_wp)) then
                     if (abs(sigma_l) < abs(sigma_r)) then
                        inflexion_l = 1
                     else
                        inflexion_r = 1
                     end if
                  end if
               else if ((x1 >= 0.0_wp) .and. (x1 <= 1.0_wp)) then
                  grad1 = 4.0_wp*eco*(x1**3) + 3.0_wp*dco*(x1**2) + 2.0_wp*cco*x1 + b
                  if (grad1*slope < 0.0_wp) then
                     if (abs(sigma_l) < abs(sigma_r)) then
                        inflexion_l = 1
                     else
                        inflexion_r = 1
                     end if
                  end if
               else if ((x2 >= 0.0_wp) .and. (x2 <= 1.0_wp)) then
                  grad2 = 4.0_wp*eco*(x2**3) + 3.0_wp*dco*(x2**2) + 2.0_wp*cco*x2 + b
                  if (grad2*slope < 0.0_wp) then
                     if (abs(sigma_l) < abs(sigma_r)) then
                        inflexion_l = 1
                     else
                        inflexion_r = 1
                     end if
                  end if
               end if
            end if

            if ((alpha1 == 0.0_wp) .and. (alpha2 /= 0.0_wp)) then
               x1 = -alpha3/alpha2
               if ((x1 >= 0.0_wp) .and. (x1 <= 1.0_wp)) then
                  grad1 = 4.0_wp*eco*(x1**3) + 3.0_wp*dco*(x1**2) + 2.0_wp*cco*x1 + b
                  if (grad1*slope < 0.0_wp) then
                     if (abs(sigma_l) < abs(sigma_r)) then
                        inflexion_l = 1
                     else
                        inflexion_r = 1
                     end if
                  end if
               end if
            end if
         end if

         if (inflexion_l == 1) then
            ! Collapse both inflexion points onto the LEFT edge.
            u1_l = (10.0_wp*u_c - 2.0_wp*u0_r - 8.0_wp*u0_l)/(3.0_wp*h_c + H_NEGLECT)
            u1_r = (-10.0_wp*u_c + 6.0_wp*u0_r + 4.0_wp*u0_l)/(h_c + H_NEGLECT)
            if (u1_l*slope < 0.0_wp) then
               u1_l = 0.0_wp
               u0_r = 5.0_wp*u_c - 4.0_wp*u0_l
               u1_r = 20.0_wp*(u_c - u0_l)/(h_c + H_NEGLECT)
            else if (u1_r*slope < 0.0_wp) then
               u1_r = 0.0_wp
               u0_l = (5.0_wp*u_c - 3.0_wp*u0_r)/2.0_wp
               u1_l = 10.0_wp*(-u_c + u0_r)/(3.0_wp*h_c + H_NEGLECT)
            end if
         else if (inflexion_r == 1) then
            ! Collapse both inflexion points onto the RIGHT edge.
            u1_r = (-10.0_wp*u_c + 8.0_wp*u0_r + 2.0_wp*u0_l)/(3.0_wp*h_c + H_NEGLECT)
            u1_l = (10.0_wp*u_c - 4.0_wp*u0_r - 6.0_wp*u0_l)/(h_c + H_NEGLECT)
            if (u1_l*slope < 0.0_wp) then
               u1_l = 0.0_wp
               u0_r = (5.0_wp*u_c - 3.0_wp*u0_l)/2.0_wp
               u1_r = 10.0_wp*(u_c - u0_l)/(3.0_wp*h_c + H_NEGLECT)
            else if (u1_r*slope < 0.0_wp) then
               u1_r = 0.0_wp
               u0_l = 5.0_wp*u_c - 4.0_wp*u0_r
               u1_l = 20.0_wp*(-u_c + u0_r)/(h_c + H_NEGLECT)
            end if
         end if

         ev_l(k) = u0_l
         ev_r(k) = u0_r
         es_l(k) = u1_l
         es_r(k) = u1_r
      end do

      ! Boundary cells: PCM (constant reconstruction).
      ev_l(1) = q_old(1)
      ev_r(1) = q_old(1)
      es_l(1) = 0.0_wp
      es_r(1) = 0.0_wp
      ev_l(nz) = q_old(nz)
      ev_r(nz) = q_old(nz)
      es_l(nz) = 0.0_wp
      es_r(nz) = 0.0_wp
      ! Opt-in boundary closure: the symmetric one-sided edge pair plus the
      ! matching constant edge slope 2d/h.  Substituting those into Step 4
      ! gives pc = pd = pe = 0 identically, so the boundary cell carries
      ! the exact straight line (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)
         ev_l(1) = q_old(1) - d_bnd
         ev_r(1) = q_old(1) + d_bnd
         es_l(1) = 2.0_wp*d_bnd/max(dz_old(1), H_NEGLECT)
         es_r(1) = es_l(1)
         call boundary_half_jump(dz_old(nz), dz_old(nz - 1), &
                                 q_old(nz) - q_old(nz - 1), d_bnd)
         ev_l(nz) = q_old(nz) - d_bnd
         ev_r(nz) = q_old(nz) + d_bnd
         es_l(nz) = 2.0_wp*d_bnd/max(dz_old(nz), H_NEGLECT)
         es_r(nz) = es_l(nz)
      end if

      ! ---- Step 4: per-cell quartic coefficients (xi in [0,1]) ----
      do k = 1, nz
         h_c = dz_old(k)
         u0_l = ev_l(k)
         u0_r = ev_r(k)
         u1_l = es_l(k)
         u1_r = es_r(k)
         u_c = q_old(k)
         pa(k) = u0_l
         pb(k) = h_c*u1_l
         pc(k) = 30.0_wp*u_c - 12.0_wp*u0_r - 18.0_wp*u0_l + 1.5_wp*h_c*(u1_r - 3.0_wp*u1_l)
         pd(k) = -60.0_wp*u_c + h_c*(6.0_wp*u1_l - 4.0_wp*u1_r) + 28.0_wp*u0_r + 32.0_wp*u0_l
         pe(k) = 30.0_wp*u_c + 2.5_wp*h_c*(u1_r - u1_l) - 15.0_wp*(u0_l + u0_r)
      end do

      ! ---- Step 5: conservative quartic overlap integration ----
      ! For old cell ko, the mean of the quartic over [xi_lo,xi_hi] is the
      ! analytic average_value_ppoly form:
      !   a + b*<xi> + c*<xi^2> + d*<xi^3> + e*<xi^4>
      ! with the symmetric power-mean expressions over the sub-interval.
      ! Multiplying the mean by the overlap thickness gives the contribution.
      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)
               integral = integral + overlap*( &
                          pa(ko) &
                          + pb(ko)*0.5_wp*(xi_lo + xi_hi) &
                          + pc(ko)*(1.0_wp/3.0_wp)*(xi_lo*xi_lo + xi_hi*xi_hi + xi_lo*xi_hi) &
                          + pd(ko)*0.25_wp*((xi_lo*xi_lo + xi_hi*xi_hi)*(xi_lo + xi_hi)) &
                          + pe(ko)*0.2_wp*((xi_hi**3 + xi_lo**3)*(xi_lo + xi_hi) &
                                           + xi_lo*xi_lo*xi_hi*xi_hi))
            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_pqm