remap_column_plm Subroutine

   pure subroutine remap_column_plm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)
      !$acc routine seq
      !! Piecewise-linear (minmod-limited) remap. Monotone (no new extrema).
      !! Per old layer k: q_hat(xi) = q(k) + slope(k)*(2*xi - 1), xi in [0,1],
      !! slope(k) = 0.5*minmod(q(k+1)-q(k), q(k)-q(k-1)).
      !! `bnd_extrap` (absent/.false. = default) closes the boundary cells
      !! with `boundary_half_jump` instead of the PCM flatten — see
      !! `remap_column`.  `nonunif` (absent/.false. = default) replaces the
      !! minmod half-difference — which assumes EQUAL source thicknesses —
      !! with the thickness-weighted CW84 (1.7)/(1.8) slope; see
      !! `plm_slope_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 slope weights (CW84 1.7/1.8).

      real(wp) :: z_old(0:NZ_STACK_MAX), z_new(0:NZ_STACK_MAX)
      real(wp) :: slope(NZ_STACK_MAX)
      real(wp) :: z_lo, z_hi, overlap, integral
      real(wp) :: xi_lo, xi_hi, dq_l, dq_r
      logical :: nu
      integer :: k, ko, ko_start

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

      ! Single-layer case: identity remap
      if (nz == 1) then
         q_new(1) = q_old(1)
         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

      ! Compute minmod-limited slopes
      ! slope(k) = half the limited difference across the layer
      slope(1) = 0.0_wp
      if (nu) then
         do k = 2, nz - 1
            call plm_slope_nonuniform(dz_old(k - 1), dz_old(k), dz_old(k + 1), &
                                      q_old(k - 1), q_old(k), q_old(k + 1), slope(k))
         end do
      else
         do k = 2, nz - 1
            dq_l = q_old(k) - q_old(k - 1)
            dq_r = q_old(k + 1) - q_old(k)
            if (dq_l*dq_r > 0.0_wp) then
               slope(k) = 0.5_wp*sign(min(abs(dq_l), abs(dq_r)), dq_l)
            else
               slope(k) = 0.0_wp
            end if
         end do
      end if
      slope(nz) = 0.0_wp
      ! Boundary cells: PCM flatten by default; the linear-exact one-sided
      ! half-jump when boundary extrapolation is requested.
      if (present(bnd_extrap)) then
         if (bnd_extrap) then
            call boundary_half_jump(dz_old(1), dz_old(2), q_old(2) - q_old(1), slope(1))
            call boundary_half_jump(dz_old(nz), dz_old(nz - 1), &
                                    q_old(nz) - q_old(nz - 1), slope(nz))
         end if
      end if

      ! Sweep: for each new layer, integrate PLM from old 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
               ! Normalised coordinates within old layer ko
               xi_lo = (z_lo - z_old(ko - 1))/dz_old(ko)
               xi_hi = (z_hi - z_old(ko - 1))/dz_old(ko)

               ! Integral of q_hat(xi) = q + slope*(2*xi - 1) over [xi_lo, xi_hi]
               ! = (xi_hi - xi_lo) * (q + slope*(xi_hi + xi_lo - 1))
               ! scaled to physical space: * dz_old(ko)
               integral = integral + overlap* &
                          (q_old(ko) + slope(ko)*(xi_lo + xi_hi - 1.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_plm