remap_layer_to_vcoord_impl Subroutine

   pure subroutine remap_layer_to_vcoord_impl(coord_type, h_layer, levels, layer_buf, &
                                              output_buf, is_extensive, method, &
                                              mask_vanished, missing)
      !! Per-column conservative layer→output-coordinate remap as a
      !! `do concurrent` over (j, i).  Builds the source column TOP-DOWN
      !! (work index 1 = surface = state k=nz) and the target cells from the
      !! `levels` interface positions (implicit 0 surface), clips both to the
      !! column total, then runs the donor-cell overlap integral
      !! (`remap_column`, scheme `method`) on a common `n = max(nz, nz_out)`
      !! padded partition.  Intensive remaps the value directly; extensive
      !! divides in / multiplies out by thickness so the column integral
      !! redistributes (sum preserved).
      !!
      !! `coord_type` selects how `levels(m)` maps to a target interface
      !! depth (a per-column scale hoisted out of the inner loop):
      !!   * `Z_FIXED`: `levels` are absolute depths (m, positive-down) —
      !!     scale 1 (the bit-identical legacy fixed-z path).
      !!   * `SIGMA`:   `levels` are cumulative fractions (0..1) — depth =
      !!     `levels(m) * col_h` (terrain-following).
      !!   * `ZSTAR`:   `levels` are reference depths (deepest = H_ref) —
      !!     depth = `levels(m) * col_h / H_ref` (SSH-tracking: when
      !!     col_h == H_ref the grid is the reference grid).  Identical to
      !!     SIGMA under uniform levels; differs only with a non-uniform
      !!     (e.g. fine-near-surface) reference.
      !!
      !! **Vanished source layers carry ZERO WEIGHT, not a value.**  A
      !! source layer that is not live (`rdb_vl_is_live`) enters the
      !! overlap integral with `dz = 0` and `q = 0`, and its thickness is
      !! left out of the column total the targets are clipped to.  Its
      !! `layer_buf` value is never read: for a concentration it is the
      !! NaN missing-data sentinel (`fill_tracer_impl`), which the
      !! donor-cell reconstruction would otherwise smear into every target
      !! cell of the column; for content it is `h·c_live` by I1′, a copy of
      !! the donor layer the remap already counts.  A column with
      !! no vanished layer is bit-identical.
      !!
      !! Fixed-size `NZ_STACK_MAX` stack locals via `local(...)` — automatic
      !! arrays sized from a dummy crash NVHPC stdpar device codegen.
      ! assumed-shape-ok: diag fill — fires once per output frame (cadence-bounded);
      ! size() used to derive loop bounds from the actual buffer dimensions.
      integer, intent(in)    :: coord_type
      real(wp), intent(in)    :: h_layer(:, :, :)
      real(wp), intent(in)    :: levels(:)  ! assumed-shape-ok: diag fill — cadence-bounded
      real(wp), intent(in)    :: layer_buf(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      real(wp), intent(inout) :: output_buf(:, :, :)
      logical, intent(in)    :: is_extensive
      integer, intent(in)    :: method
      logical, intent(in)    :: mask_vanished
      real(wp), intent(in)    :: missing
      integer :: i, j, k, m, nx, ny, nz, nz_out, n
      real(wp) :: dz_old(NZ_STACK_MAX), dz_new(NZ_STACK_MAX)
      real(wp) :: q_old(NZ_STACK_MAX), q_new(NZ_STACK_MAX)
      real(wp) :: col_h, zf_prev, zf, dz, lvl_scale, href
      nx = size(layer_buf, 1)
      ny = size(layer_buf, 2)
      nz = size(layer_buf, 3)
      nz_out = size(levels)
      n = max(nz, nz_out)
      href = levels(nz_out)   ! ZSTAR reference total depth (deepest interface)
      do concurrent(j=1:ny, i=1:nx) &
         local(dz_old, dz_new, q_old, q_new, k, m, col_h, zf_prev, zf, dz, lvl_scale)
         ! --- source column TOP-DOWN: work index k = state index nz-k+1 ---
         ! A vanished layer gets zero weight (see the docstring).
         col_h = 0.0_wp
         do k = 1, nz
            dz_old(k) = h_layer(i, j, nz - k + 1)
            if (rdb_vl_is_live(dz_old(k))) then
               q_old(k) = layer_buf(i, j, nz - k + 1)
            else
               dz_old(k) = 0.0_wp
               q_old(k) = 0.0_wp
            end if
            col_h = col_h + dz_old(k)
         end do
         if (is_extensive) then
            do k = 1, nz
               if (dz_old(k) > 1.0e-12_wp) then
                  q_old(k) = q_old(k)/dz_old(k)   ! integral -> concentration
               else
                  q_old(k) = 0.0_wp
               end if
            end do
         end if
         ! pad the source to n with zero-thickness layers at the deep end
         do k = nz + 1, n
            dz_old(k) = 0.0_wp
            q_old(k) = 0.0_wp
         end do
         ! --- per-column level->depth scale (loop-invariant, hoisted) ---
         select case (coord_type)
         case (DIAG_VGRID_SIGMA)
            lvl_scale = col_h
         case (DIAG_VGRID_ZSTAR)
            if (href > 1.0e-12_wp) then
               lvl_scale = col_h/href
            else
               lvl_scale = 1.0_wp
            end if
         case default   ! DIAG_VGRID_Z_FIXED: levels are absolute depths
            lvl_scale = 1.0_wp
         end select
         ! --- target cells from scaled interfaces (implicit 0 surface),
         !     clipped to the column total H so the totals match exactly ---
         zf_prev = 0.0_wp
         do m = 1, nz_out
            zf = min(max(levels(m)*lvl_scale, 0.0_wp), col_h)
            dz = zf - zf_prev
            if (dz < 0.0_wp) dz = 0.0_wp
            dz_new(m) = dz
            zf_prev = zf
         end do
         do m = nz_out + 1, n
            dz_new(m) = 0.0_wp
         end do
         call remap_column(method, n, dz_old, dz_new, q_old, q_new)
         do m = 1, nz_out
            if (mask_vanished .and. dz_new(m) <= VANISHED_TARGET_FLOOR) then
               ! Target cell overlaps no water (below-bottom / pinched-out):
               ! emit the missing sentinel instead of a misleading 0.
               output_buf(i, j, m) = missing
            else if (is_extensive) then
               output_buf(i, j, m) = q_new(m)*dz_new(m)
            else
               output_buf(i, j, m) = q_new(m)
            end if
         end do
      end do
   end subroutine remap_layer_to_vcoord_impl