remap_layer_to_density_impl Subroutine

   pure subroutine remap_layer_to_density_impl(h_layer, hT, hS, eos, rho_ref_p, &
                                               rho_tgt, layer_buf, output_buf, &
                                               is_extensive, method)
      !! Per-column density-space remap, `do concurrent` over (j, i).
      !! Source column TOP-DOWN + layer potential density (EOS at
      !! `rho_ref_p`), invert profile to interface depths
      !! (`invert_density_targets`), then donor-cell remap.  Cells outside
      !! the column density range read 0.  Intensive/extensive as the z-remap.
      !!
      !! **Vanished source layers** carry zero weight exactly as in
      !! `remap_layer_to_vcoord_impl` (`dz = 0`, `q = 0`, value never read),
      !! and take the density of the nearest LIVE layer (the one above,
      !! else the first one below) rather than an EOS evaluation of the
      !! filler: the EOS substitutes reference T/S on a vanished layer, and
      !! the PPM edge between a live layer and a zero-thickness neighbour IS
      !! that neighbour's density (`invert_density_targets`), so a `rho_0`
      !! filler would kink the profile the targets are inverted against.
      !! A column with no vanished layer is bit-identical; a column with no
      !! live layer keeps the legacy `eos(0, 0)` fill.
      ! assumed-shape-ok: diag fill — fires once per output frame (cadence-bounded);
      ! size() used to derive loop bounds from the actual buffer dimensions.
      real(wp), intent(in)    :: h_layer(:, :, :)
      real(wp), intent(in)    :: hT(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      real(wp), intent(in)    :: hS(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      type(eos_t), intent(in) :: eos
      real(wp), intent(in)    :: rho_ref_p
      real(wp), intent(in)    :: rho_tgt(:)  ! 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
      integer :: i, j, k, m, nx, ny, nz, n_bin, 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)
      ! z_iface holds n_bin+2 entries (`invert_density_targets` writes
      ! z_new(1:n_int+2) with n_int = n_bin, and :798 reads z_iface(n_bin+2)),
      ! so it needs NZ_STACK_MAX+2 — at NZ_STACK_MAX+1 a density diagnostic
      ! with n_rho_out == NZ_STACK_MAX (which the output-level guard permits)
      ! wrote and read one element past the end.
      real(wp) :: rhoc(NZ_STACK_MAX), z_iface(NZ_STACK_MAX + 2)
      real(wp) :: rho_tgt_c(NZ_STACK_MAX)
      real(wp) :: hh, tt, ss, rho_live
      logical :: seen_live
      nx = size(layer_buf, 1)
      ny = size(layer_buf, 2)
      nz = size(layer_buf, 3)
      n_bin = size(rho_tgt)      ! number of density bins (output cells)
      n = max(nz, n_bin)
      do concurrent(j=1:ny, i=1:nx) &
         local(dz_old, dz_new, q_old, q_new, rhoc, z_iface, rho_tgt_c, &
               k, m, hh, tt, ss, rho_live, seen_live)
         ! --- source column TOP-DOWN + layer potential density ---
         ! Vanished layers: zero weight, density of the nearest live layer
         ! (see the docstring).
         seen_live = .false.
         rho_live = 0.0_wp
         do k = 1, nz
            hh = h_layer(i, j, nz - k + 1)
            if (rdb_vl_is_live(hh)) then
               dz_old(k) = hh
               q_old(k) = layer_buf(i, j, nz - k + 1)
               tt = rdb_vl_conc(hT(i, j, nz - k + 1), hh)
               ss = rdb_vl_conc(hS(i, j, nz - k + 1), hh)
               rhoc(k) = eos_density_point(eos, tt, ss, rho_ref_p)
               if (.not. seen_live) then
                  ! Back-fill the vanished run above the first live layer.
                  do m = 1, k - 1
                     rhoc(m) = rhoc(k)
                  end do
               end if
               seen_live = .true.
               rho_live = rhoc(k)
            else
               dz_old(k) = 0.0_wp
               q_old(k) = 0.0_wp
               if (seen_live) then
                  rhoc(k) = rho_live
               else
                  rhoc(k) = eos_density_point(eos, 0.0_wp, 0.0_wp, rho_ref_p)
               end if
            end if
         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)
               else
                  q_old(k) = 0.0_wp
               end if
            end do
         end if
         do k = nz + 1, n
            dz_old(k) = 0.0_wp
            q_old(k) = 0.0_wp
         end do
         ! --- invert instantaneous density to target-interface depths ---
         ! n_bin output cells; the LAST cell absorbs to the bed so
         ! Σ dz_new == H exactly (extensive conservation, densest mass kept).
         ! Copy assumed-shape rho_tgt(:) into a contiguous fixed-size local
         ! first: passing the assumed-shape actual to the explicit-shape
         ! dummy makes flang emit host copy helpers absent in the AMD device
         ! runtime → offload link fails.  Element reads off the descriptor OK.
         do m = 1, n_bin
            rho_tgt_c(m) = rho_tgt(m)
         end do
         call invert_density_targets(nz, dz_old, rhoc, n_bin, rho_tgt_c, z_iface)
         do m = 1, n_bin - 1
            dz_new(m) = z_iface(m + 1) - z_iface(m)
            if (dz_new(m) < 0.0_wp) dz_new(m) = 0.0_wp
         end do
         dz_new(n_bin) = z_iface(n_bin + 2) - z_iface(n_bin)
         if (dz_new(n_bin) < 0.0_wp) dz_new(n_bin) = 0.0_wp
         do m = n_bin + 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, n_bin
            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_density_impl