fill_rho_layer_impl Subroutine

   pure subroutine fill_rho_layer_impl(h_layer, rho_layer, buf)
      !! `buf = rho_layer` on live layers, IEEE NaN on vanished ones.
      !!
      !! The EOS deliberately evaluates a vanished layer at the reference
      !! `T_ref`/`S_ref`, so `rho_layer` holds `rho_0` there (README
      !! consumer table: a filler must not perturb the PGF's density
      !! column).  That is a plausible-looking density, not a measurement,
      !! so the diagnostic reports it as missing — the same convention
      !! `fill_tracer_impl` applies to the T and S it is computed from.
      ! assumed-shape-ok: diag fill — fires once per output frame (cadence-bounded).
      real(wp), intent(in)    :: h_layer(:, :, :)
      real(wp), intent(in)    :: rho_layer(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      real(wp), intent(inout) :: buf(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      integer :: i, j, k, nx, ny, nz
      real(wp) :: qnan
      nx = min(size(buf, 1), size(rho_layer, 1), size(h_layer, 1))
      ny = min(size(buf, 2), size(rho_layer, 2), size(h_layer, 2))
      nz = min(size(buf, 3), size(rho_layer, 3), size(h_layer, 3))
      ! Sentinel computed once on the HOST (the `fill_tracer_impl` pattern).
      qnan = ieee_value(0.0_wp, ieee_quiet_nan)
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         ! vanished-ok: diagnostics report a vanished layer as NaN, not
         ! the EOS's reference-density substitution.
         if (rdb_vl_is_live(h_layer(i, j, k))) then
            buf(i, j, k) = rho_layer(i, j, k)
         else
            buf(i, j, k) = qnan
         end if
      end do
   end subroutine fill_rho_layer_impl