redi_build_column Subroutine

   pure subroutine redi_build_column(nz, kb, kt, h_col, thtr_col, shtr_col, eos, &
                                     Pint, Tint, Sint, dRdT, dRdS)
      !! Build one column's TOP-DOWN interface P/T/S + density derivs from the
      !! BOTTOM-UP native column, restricted to the face's open window
      !! `kb..kt` (`nk = kt-kb+1` layers; the whole column `1..nz` off the
      !! z-level closed-face path).  Layer T/S = the I1′ column read
      !! (`rdb_vl_column_conc`: `hTr/h` on a live layer, the donor's
      !! concentration on a vanished one); interface
      !! T/S = PPM edge reconstruction on the flipped window; interface P =
      !! surface-relative hydrostatic, seeded with the column ABOVE the
      !! window (fillers under an ice draft; nothing — exactly 0 — when
      !! `kt = nz`); dR/dT, dR/dS = -rho² dSV/dX at each interface.  Only
      !! the first `nk+1` entries of the outputs are written.
      !$acc routine seq
      integer, intent(in) :: nz
      integer, intent(in) :: kb, kt
         !! Native (bottom-up) open window, `1 <= kb <= kt <= nz`.
      real(wp), intent(in) :: h_col(nz), thtr_col(nz), shtr_col(nz)
      type(eos_t), intent(in) :: eos
      real(wp), intent(out) :: Pint(nz + 1), Tint(nz + 1), Sint(nz + 1)
      real(wp), intent(out) :: dRdT(nz + 1), dRdS(nz + 1)
      real(wp) :: htd(NZ_STACK_MAX), ttd(NZ_STACK_MAX), std(NZ_STACK_MAX)
      real(wp) :: hs(NZ_STACK_MAX), qs(NZ_STACK_MAX), cT(NZ_STACK_MAX), cS(NZ_STACK_MAX)
      real(wp) :: rho_i, dsv_dt, dsv_ds
      integer :: k, kf, nk

      nk = kt - kb + 1
      ! Layer T/S from the FULL native column by the I1′ rule
      ! (`rdb_vl_column_conc`): a vanished layer reads its donor's
      ! concentration, never `hTr/max(h, H_DIV_EPS)`.  The floored divide is
      ! `hTr/1e-20` on a layer the continuity step has taken to (or below)
      ! zero, which on an OPEN z-like step (closed faces off, so the window
      ! is the whole column and fillers are paired) read T ~ 1e17 and put
      ! 1e12-1e27 of content into the face (1-degree Southern Ocean, zstar,
      ! step 13).  A live layer reads `hTr/h`, the same divide as before.
      do k = 1, nz
         hs(k) = h_col(k)
         qs(k) = thtr_col(k)
      end do
      call rdb_vl_column_conc(nz, hs, qs, cT)
      do k = 1, nz
         qs(k) = shtr_col(k)
      end do
      call rdb_vl_column_conc(nz, hs, qs, cS)
      ! Flip the native bottom-up window -> top-down layer arrays.  Native
      ! layer k maps to top-down layer (kt+1-k).
      do k = kb, kt
         kf = kt + 1 - k
         htd(kf) = h_col(k)
         ttd(kf) = cT(k)
         std(kf) = cS(k)
      end do

      ! PPM interface edge values (top-down) for T and S.
      call redi_interface_scalar(nk, htd, ttd, Tint)
      call redi_interface_scalar(nk, htd, std, Sint)

      ! Top-down interface pressure: K=1 the window top — the surface (p=0)
      ! unless layers sit above the window — accumulate g*rho0*h downward.
      ! rho0 from the EOS reference.
      Pint(1) = 0.0_wp
      do k = nz, kt + 1, -1
         Pint(1) = Pint(1) + GRAVITY*eos%rho0*h_col(k)
      end do
      do k = 1, nk
         Pint(k + 1) = Pint(k) + GRAVITY*eos%rho0*htd(k)
      end do

      ! Interface density derivs (locally referenced at the interface P).
      do k = 1, nk + 1
         call eos_density_specvol_derivs(eos, Tint(k), Sint(k), Pint(k), rho_i, dsv_dt, dsv_ds)
         dRdT(k) = -(rho_i*rho_i)*dsv_dt
         dRdS(k) = -(rho_i*rho_i)*dsv_ds
      end do
   end subroutine redi_build_column