remap_column_ppm Subroutine

   pure subroutine remap_column_ppm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)
      !$acc routine seq
      !! Piecewise-parabolic (Colella & Woodward 1984) remap.
      !! Per old layer k, xi in [0,1]:
      !!   q_hat(xi) = q_L + xi*(q_R - q_L + q6*(1 - xi)), q6 = 6*q_bar - 3*(q_L+q_R)
      !! Edge values: 4th-order interp + CW monotonicity limiting; boundary
      !! layers fall back to PLM-quality edges, or — under `bnd_extrap` —
      !! to the linear-exact one-sided pair (`boundary_half_jump`), which
      !! zeroes `q6` there so the boundary cell carries a straight line.
      !! `nonunif` swaps the `(7/12, -1/12)` edge estimate — an EQUAL-
      !! thickness specialisation — for CW84 (1.6) on the true stencil
      !! thicknesses, and the `1|2` / `(nz-1)|nz` edges for the
      !! thickness-weighted two-cell value; see `ppm_edge_nonuniform`.
      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)
      logical, intent(in), optional :: bnd_extrap
         !! Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
      logical, intent(in), optional :: nonunif
         !! Non-uniform-grid edge weights (CW84 1.6-1.8).

      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) :: dq_cw(NZ_STACK_MAX)
      real(wp) :: z_lo, z_hi, overlap, integral
      real(wp) :: xi_lo, xi_hi
      real(wp) :: edge, dq, dq_l, dq_r, q_min, q_max
      real(wp) :: d_bnd
      logical :: be, nu
      integer :: k, ko, ko_start

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

      ! Trivial cases
      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: Compute unlimited edge values via 4th-order interp ----
      ! Interior edges (between layers k and k+1) use the 4-cell stencil.
      ! Store in q_R(k) = right edge of layer k = left edge of layer k+1.

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

      if (nu) then
         ! ---- Non-uniform weights: CW84 (1.6) on the true thicknesses ----
         ! Limited per-cell jumps (1.7)+(1.8) feed the (1.6) correction; they
         ! exist only where a centred triple does, k = 2 .. nz-1.
         do k = 2, nz - 1
            call ppm_jump_nonuniform(dz_old(k - 1), dz_old(k), dz_old(k + 1), &
                                     q_old(k - 1), q_old(k), q_old(k + 1), dq_cw(k))
         end do
         ! The two edges CW84 (1.6) has no stencil for: thickness-weighted
         ! two-cell interpolation, which is linear-exact (the non-uniform
         ! generalisation of the 0.5 average the uniform path uses there).
         call ppm_edge_two_cell(dz_old(1), dz_old(2), q_old(1), q_old(2), edge)
         q_R(1) = edge
         q_L(2) = edge
         call ppm_edge_two_cell(dz_old(nz - 1), dz_old(nz), q_old(nz - 1), q_old(nz), edge)
         q_R(nz - 1) = edge
         q_L(nz) = edge
         ! Interior edges with the full four-cell stencil.
         do k = 2, nz - 2
            call ppm_edge_nonuniform(dz_old(k - 1), dz_old(k), dz_old(k + 1), dz_old(k + 2), &
                                     q_old(k), q_old(k + 1), dq_cw(k), dq_cw(k + 1), edge)
            q_R(k) = edge
            q_L(k + 1) = edge
         end do
      else
         ! Layer 1 right edge = layer 2 left edge: use 3-cell stencil (one-sided)
         q_R(1) = 0.5_wp*(q_old(1) + q_old(2))
         ! Layer nz left edge = layer nz-1 right edge: use 3-cell stencil
         q_L(nz) = 0.5_wp*(q_old(nz - 1) + q_old(nz))

         ! Interior edges: 4th-order Colella-Woodward interpolation
         ! For uniform layers this gives (7/12)(q_k + q_{k+1}) - (1/12)(q_{k-1} + q_{k+2})
         ! For non-uniform layers, use the simpler weighted average
         do k = 2, nz - 1
            edge = 0.5_wp*(q_old(k) + q_old(k + 1))
            if (k >= 2 .and. k + 1 <= nz) then
               ! Add 4th-order correction when stencil is available
               dq_l = q_old(k) - q_old(k - 1)
               dq_r = q_old(k + 1) - q_old(k)
               if (k - 1 >= 1 .and. k + 2 <= nz) then
                  edge = (7.0_wp/12.0_wp)*(q_old(k) + q_old(k + 1)) &
                         - (1.0_wp/12.0_wp)*(q_old(k - 1) + q_old(k + 2))
               end if
            end if
            q_R(k) = edge
            q_L(k + 1) = edge
         end do

         ! Layer 2 left edge (if nz >= 3, was set above; otherwise use average)
         if (nz >= 3) then
            q_L(2) = q_R(1)
         end if
      end if

      ! ---- Step 2: Colella-Woodward monotonicity limiting ----
      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 = 6*q_bar - 3*(q_L + q_R)
            ! The parabola has an extremum inside [0,1] if q6*(q_R - q_L) < 0
            ! and the extremum value exceeds the local bounds.
            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) ----
      ! The CW limiter above bounds every edge by the cell means it can
      ! see, and at k=1 / k=nz that is a ONE-SIDED bound, so the default
      ! closure collapses those two cells to PCM.  With extrapolation on,
      ! the symmetric one-sided pair replaces it and
      ! q6 = 6q - 3(q_L + q_R) = 0, so the boundary cell carries the exact
      ! straight line whenever q(z) is linear.  Deliberately written AFTER
      ! the limiter: the one-sided clip is precisely what has to be
      ! bypassed here.
      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 ----
      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