ocean_slopes_pass_x Subroutine

   pure subroutine ocean_slopes_pass_x(nx, ny, nz, eos, rho0, min_dz, &
                                       h_layer, t_fill, s_fill, e_int, &
                                       idxCu, wet_u, slope_x, n2_u)
      !! u-face slope + N² pass.  Interface `K` (interior 2..nz) straddles
      !! layer `k=K` (above, surface side) and `k=K-1` (below, bed side).
      !! Bed (K=1) + surface (K=nz+1) are forced to zero.  The u-face at
      !! (i,j) sits between cells (i-1,j) and (i,j); pairs columns
      !! `iw=i-1` (west) and `i` (east), so loop `i=2:nx`.
      integer, intent(in) :: nx, ny, nz
      type(eos_t), intent(in) :: eos
      real(wp), intent(in) :: rho0, min_dz
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: t_fill(nx, ny, nz)
      real(wp), intent(in) :: s_fill(nx, ny, nz)
      real(wp), intent(in) :: e_int(nx, ny, nz + 1)
      real(wp), intent(in) :: idxCu(nx + 1, ny)
      real(wp), intent(in) :: wet_u(nx + 1, ny)
      real(wp), intent(out) :: slope_x(nx + 1, ny, nz + 1)
      real(wp), intent(out) :: n2_u(nx + 1, ny, nz + 1)

      integer :: i, j, k, iw, ka, kb
      real(wp) :: pres_u, t_u, s_u, rho_u, dsv_dt, dsv_ds, drdt, drds
      real(wp) :: drdiA, drdiB, drdkL, drdkR
      real(wp) :: hg2A, hg2B, hg2L, hg2R, haA, haB, haL, haR
      real(wp) :: dzaL, dzaR, wtA, wtB, wtL, wtR
      real(wp) :: drdx, drdz, mag2, slope, presL, presR
      real(wp) :: g_rho0, mask

      g_rho0 = GRAVITY/rho0

      ! Bed + surface interfaces: zero everywhere.
      do concurrent(j=1:ny, i=1:nx + 1)
         slope_x(i, j, 1) = 0.0_wp
         slope_x(i, j, nz + 1) = 0.0_wp
         n2_u(i, j, 1) = 0.0_wp
         n2_u(i, j, nz + 1) = 0.0_wp
      end do
      ! Wall faces (i=1, i=nx+1): zero at all interfaces (no interior pair).
      do concurrent(k=1:nz + 1, j=1:ny)
         slope_x(1, j, k) = 0.0_wp
         slope_x(nx + 1, j, k) = 0.0_wp
         n2_u(1, j, k) = 0.0_wp
         n2_u(nx + 1, j, k) = 0.0_wp
      end do

      ! Interior interfaces K = 2..nz, interior u-faces i = 2..nx.
      do concurrent(k=2:nz, j=1:ny, i=2:nx) &
         local(iw, ka, kb, pres_u, t_u, s_u, rho_u, dsv_dt, dsv_ds, &
               drdt, drds, drdiA, drdiB, drdkL, drdkR, &
               hg2A, hg2B, hg2L, hg2R, haA, haB, haL, haR, &
               dzaL, dzaR, wtA, wtB, wtL, wtR, drdx, drdz, &
               mag2, slope, presL, presR, mask)
         iw = i - 1
         ka = k          ! layer ABOVE the interface (surface side)
         kb = k - 1      ! layer BELOW the interface (bed side)

         ! Interface pressure: accumulate from the surface (k=nz) down to
         ! the layer above this interface.  Surface-relative hydrostatic
         ! pressure at the interface = g·ρ₀·Σ_{above} h.
         presL = pressure_above_x(nx, ny, nz, h_layer, iw, j, ka, rho0)
         presR = pressure_above_x(nx, ny, nz, h_layer, i, j, ka, rho0)
         pres_u = 0.5_wp*(presL + presR)

         ! 4-point interface T/S (two columns × two adjacent layers).
         t_u = 0.25_wp*((t_fill(iw, j, ka) + t_fill(i, j, ka)) + &
                        (t_fill(iw, j, kb) + t_fill(i, j, kb)))
         s_u = 0.25_wp*((s_fill(iw, j, ka) + s_fill(i, j, ka)) + &
                        (s_fill(iw, j, kb) + s_fill(i, j, kb)))

         ! Locally-referenced density derivatives: drho_dX = -ρ²·dSV/dX.
         call eos_density_specvol_derivs(eos, t_u, s_u, pres_u, rho_u, dsv_dt, dsv_ds)
         drdt = -(rho_u*rho_u)*dsv_dt
         drds = -(rho_u*rho_u)*dsv_ds

         ! Along-layer horizontal ρ-gradients, above (A=ka) / below (B=kb).
         drdiA = drdt*(t_fill(i, j, ka) - t_fill(iw, j, ka)) + &
                 drds*(s_fill(i, j, ka) - s_fill(iw, j, ka))
         drdiB = drdt*(t_fill(i, j, kb) - t_fill(iw, j, kb)) + &
                 drds*(s_fill(i, j, kb) - s_fill(iw, j, kb))

         ! Vertical ρ-difference (below - above): drho_dX·(X[kb]-X[ka]).
         ! For stable stratification (lighter water above) this gives
         ! drdk>0 ⇒ drdz>0 ⇒ N²>0 (punch-list sign fix #2).
         drdkL = drdt*(t_fill(iw, j, kb) - t_fill(iw, j, ka)) + &
                 drds*(s_fill(iw, j, kb) - s_fill(iw, j, ka))
         drdkR = drdt*(t_fill(i, j, kb) - t_fill(i, j, ka)) + &
                 drds*(s_fill(i, j, kb) - s_fill(i, j, ka))

         ! Harmonic-mean thickness weights.
         hg2A = h_layer(iw, j, ka)*h_layer(i, j, ka) + H_DIV_EPS*H_DIV_EPS
         hg2B = h_layer(iw, j, kb)*h_layer(i, j, kb) + H_DIV_EPS*H_DIV_EPS
         hg2L = h_layer(iw, j, ka)*h_layer(iw, j, kb) + H_DIV_EPS*H_DIV_EPS
         hg2R = h_layer(i, j, ka)*h_layer(i, j, kb) + H_DIV_EPS*H_DIV_EPS
         haA = 0.5_wp*(h_layer(iw, j, ka) + h_layer(i, j, ka)) + H_DIV_EPS
         haB = 0.5_wp*(h_layer(iw, j, kb) + h_layer(i, j, kb)) + H_DIV_EPS
         haL = 0.5_wp*(h_layer(iw, j, ka) + h_layer(iw, j, kb)) + H_DIV_EPS
         haR = 0.5_wp*(h_layer(i, j, ka) + h_layer(i, j, kb)) + H_DIV_EPS
         ! Vertical centre spacing across the interface (floored).
         dzaL = max(haL, min_dz)
         dzaR = max(haR, min_dz)
         wtA = hg2A*haB
         wtB = hg2B*haA
         wtL = hg2L*(haR*dzaR)
         wtR = hg2R*(haL*dzaL)

         drdz = ((wtL*drdkL) + (wtR*drdkR))/((dzaL*wtL) + (dzaR*wtR))

         ! Interface-tilt rotation term + metric scaling.  `e_int` is
         ! geopotential (bed datum −D), so `e_W − e_E` is the real tilt of
         ! the interface, never the bathymetry step.
         drdx = ((wtA*drdiA + wtB*drdiB)/(wtA + wtB) - &
                 drdz*(e_int(iw, j, k) - e_int(i, j, k)))*idxCu(i, j)

         mag2 = drdx*drdx + drdz*drdz
         if (mag2 > 0.0_wp) then
            slope = drdx/sqrt(mag2)
         else
            slope = 0.0_wp
         end if

         mask = wet_u(i, j)
         slope_x(i, j, k) = slope*mask
         n2_u(i, j, k) = g_rho0*drdz*mask
      end do
   end subroutine ocean_slopes_pass_x