eos_wright_pgf_column_sweep_impl Subroutine

   pure subroutine eos_wright_pgf_column_sweep_impl(h_layer, hS_layer, hT_layer, &
                                                    rho_layer_seed, p_top, p_edge_out, &
                                                    rho_insitu_out, gravity, rho_0, &
                                                    nx, ny, nz)
      !! FV-Wright PGF column sweep.  Top-down per-column traversal
      !! that simultaneously produces the hydrostatic pressure stack
      !! `p_edge_out` and the in-situ density `rho_insitu_out` at each
      !! layer centre.
      !!
      !! For each layer k from the surface (k=nz) down to the bed
      !! (k=1) we use a single Picard step:
      !!
      !!   1. Seed the half-layer pressure with the potential density
      !!      `rho_layer_seed` (= ms%rho_layer, the EOS output at the
      !!      uniform `eos%p_ref`):
      !!         p_centre_seed = p_top + p_above + 0.5 * g * rho_seed * h
      !!   2. Re-evaluate Wright at the seed pressure:
      !!         rho_insitu(k) = ρ(T_k, S_k, p_centre_seed)
      !!   3. Accumulate the bottom edge of layer k using the in-situ
      !!      density:
      !!         p_edge(k) = p_above + g * rho_insitu(k) * h
      !!
      !! This is the ONE genuinely IN-SITU pressure the EOS core builds:
      !! `p_centre_seed` is a true per-layer hydrostatic pressure, so a
      !! per-column `p_top(i,j)` belongs in it (`&ocean_psurf_nml in_eos`).
      !! It is NOT the potential-density trap that keeps `eos%p_ref` a
      !! scalar: `rho_insitu` is consumed only (a) as the integrand of
      !! THIS column's `p_edge` stack and (b) as `rho_face`, the
      !! two-point AVERAGE coefficient multiplying `Δz_centre` in the
      !! PGF's z-correction (`ocean_pressure_force_compute` Pass 2/3).
      !! Neither differences it along a layer, so a horizontally varying
      !! `p_top` cannot manufacture an along-layer density gradient here —
      !! and the density really IS higher under a thicker draft, so the
      !! `rho_face` coefficient becomes MORE correct, not less.
      !!
      !! **`p_top` reaches the EOS ARGUMENT only.**  `p_edge_out` stays an
      !! anomaly stack seeded at `p_edge_out(nz+1) = 0` exactly as before,
      !! so the PGF top boundary condition THIS kernel feeds (FV_WRIGHT's
      !! `p_edge`) is untouched.  Consequence, stated plainly: under a
      !! SLOPING load the along-layer difference
      !! `p_centre(i) − p_centre(i−1)` omits `Δp_top`.  That term is
      !! depth-uniform and is already carried by the barotropic
      !! `eta_forcing` seam as `−(1/ρ₀)∇p_surf`, so the momentum is not
      !! missing it.
      !!
      !! **Amended (P5.0).**  The original wording here said adding the
      !! load to a PGF top BC "would DOUBLE-COUNT".  That is the
      !! conservative statement, and it is stronger than the truth.  A
      !! depth-uniform `p_top` in the top BC perturbs EVERY layer's `PFu`
      !! by the SAME `−(1/ρ₀)∇p_top`, and the split solver replaces the
      !! depth mean of the layer PGF with the barotropic solution
      !! (`F_bt_u_fast = F_bt_u − ⟨PFu⟩_h`), so the uniform piece cancels
      !! identically and the seam keeps sole ownership of the barotropic
      !! response — the two are ORTHOGONAL, not additive.  That is what
      !! `&ocean_pgf_nml p_top_in_bc` does for FV_MOM6 (theorem in
      !! `compute_fv_mom6_impl`'s docstring).  It is NOT done here:
      !! FV_WRIGHT's `p_edge` seed is a separate follow-up, and this
      !! kernel's contract remains "EOS argument only".  What moves in
      !! this kernel is the COMPRESSIBILITY: `rho_insitu` is
      !! evaluated at the pressure the water actually sits at, which is
      !! the ~4-5 kg/m^3 systematic error an ice-shelf load introduces.
      !! Bit-identical when `p_top` is the zero array it ships as
      !! (`p_top + p_above` is `p_above` exactly under IEEE-754).
      !!
      !! This is one Picard iteration of the implicit
      !!   p_centre(k) = p_above + 0.5*g*ρ(T, S, p_centre(k))*h.
      !! For ocean conditions (Δρ along path << ρ_0) one iteration is
      !! within ~1e-5 of the converged value.
      !!
      !! `rho_layer_seed` must be `ms%rho_layer` from
      !! `ocean_eos_compute` with `EOS_VARIANT_WRIGHT_97` —
      !! the seed is the *full* nonlinear ρ at `eos%p_ref`, not a
      !! Boussinesq constant; this avoids a second Picard iteration in
      !! 99% of cases.  The seed enters only a HALF-LAYER increment, so
      !! a seed offset `Δρ` costs `0.5·g·Δρ·h` of pressure (≈ 5 kPa out
      !! of 1e7 Pa for `Δρ = 5`, i.e. ~5e-4 relative); referencing
      !! `p_ref` near the working pressure makes the seed better still.
      !!
      !! Vanishing-layer fallback: if `h_layer(k) <= 0` the Wright eval
      !! is skipped and `rho_insitu(k) = rho_0` — matches the existing
      !! `eos_wright_impl` defensive branch.
      !!
      !! Loop order: outer `do concurrent (j, i)` for GPU parallelism;
      !! inner serial k loop for the column recurrence (same shape as
      !! the existing PGF Pass 1 + vdiff column kernels).
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: hS_layer(nx, ny, nz)
      real(wp), intent(in) :: hT_layer(nx, ny, nz)
      real(wp), intent(in) :: rho_layer_seed(nx, ny, nz)
      real(wp), intent(in) :: p_top(nx, ny)
         !! Top-of-column pressure (Pa, `>= 0`) the EOS argument is
         !! measured down from.  Does NOT enter `p_edge_out`.
      real(wp), intent(out) :: p_edge_out(nx, ny, nz + 1)
      real(wp), intent(out) :: rho_insitu_out(nx, ny, nz)
      real(wp), intent(in) :: gravity, rho_0

      integer :: i, j, k
      real(wp) :: p_above, p_centre_seed, p_top_ij, inv_h, S_k, T_k, T_sq, T_cu
      real(wp) :: alpha_0, p_0, lambda, p_plus_p0, denom, rho_k

      do concurrent(j=1:ny, i=1:nx) &
         local(k, p_above, p_centre_seed, p_top_ij, inv_h, S_k, T_k, T_sq, T_cu, &
               alpha_0, p_0, lambda, p_plus_p0, denom, rho_k)
         p_edge_out(i, j, nz + 1) = 0.0_wp
         p_above = 0.0_wp
         p_top_ij = p_top(i, j)
         do k = nz, 1, -1
            if (h_layer(i, j, k) > H_VANISHED) then
               inv_h = 1.0_wp/h_layer(i, j, k)
               S_k = hS_layer(i, j, k)*inv_h
               T_k = hT_layer(i, j, k)*inv_h
               T_sq = T_k*T_k
               T_cu = T_sq*T_k

               p_centre_seed = p_top_ij + p_above + &
                               0.5_wp*gravity*rho_layer_seed(i, j, k)*h_layer(i, j, k)

               alpha_0 = WRIGHT_A0 + WRIGHT_A1*T_k + WRIGHT_A2*S_k
               p_0 = WRIGHT_B0 + WRIGHT_B1*T_k + WRIGHT_B2*T_sq + WRIGHT_B3*T_cu + &
                     WRIGHT_B4*S_k + WRIGHT_B5*S_k*T_k
               lambda = WRIGHT_C0 + WRIGHT_C1*T_k + WRIGHT_C2*T_sq + WRIGHT_C3*T_cu + &
                        WRIGHT_C4*S_k + WRIGHT_C5*S_k*T_k

               p_plus_p0 = p_centre_seed + p_0
               denom = lambda + alpha_0*p_plus_p0
               rho_k = p_plus_p0/denom
            else
               rho_k = rho_0
            end if
            rho_insitu_out(i, j, k) = rho_k
            p_above = p_above + gravity*rho_k*h_layer(i, j, k)
            p_edge_out(i, j, k) = p_above
         end do
      end do
   end subroutine eos_wright_pgf_column_sweep_impl