gm_column_x Subroutine

   pure subroutine gm_column_x(nx, ny, nz, i_smax2, i4dt, use_open, &
                               dy_cu, areaT, h_layer, bathy, slope_x, khth_u, open_u, uhD)
      !! u-face GM streamfunction + bolus-transport column recurrence.
      !! Interior u-face (i=2..nx) pairs columns iw=i-1 (west) and i (east).
      !! Bottom-up sweep: interior interfaces Kr=2 (bed-most) -> nz
      !! (surface-most), uhtot=0 at the bed; the surface BC (Sfn=0 at
      !! Kr=nz+1) is closed after the loop by `uhD(nz)=-uhtot`, giving
      !! `Sum_k uhD=0` exactly.  Interface Kr straddles ka=Kr (above) and
      !! kb=Kr-1 (below) and fills LAYER kb; the rsum bound keys on ka, the
      !! donor `h_frac` and per-layer cap on kb.
      !!
      !! OPEN COLUMN (`use_open`, `&vcoord_nml zfixed_closed_faces`): a
      !! layer is part of the face column only if `ok(k) = open_u(i,j,k)`
      !! and it is live on both sides.  A not-`ok` layer gets `h_avail = 0`
      !! and `uhD = 0`, leaving `uhtot` (the streamfunction at its top)
      !! unchanged, and the closure lands in the topmost `ok` layer
      !! `ktop` instead of `nz` — Sfn = 0 at the bottom AND the top of the
      !! open column (see the module docstring).  `use_open = .false.` ⇒
      !! `ok` all-true, `ktop = nz`: the full-column recurrence verbatim.
      !!
      !! DIVERGENCE (MOM6 nk_linear): MOM6's `thickness_diffuse_full`
      !! sets `nk_linear = max(GV%nkml, 1)`; in
      !! ALE mode `GV%nkml = 0`, so MOM6 always runs `nk_linear = 1` and it
      !! is NOT a namelist parameter — its top layer always takes a linear
      !! return-flow closure (comment: "Balance the deeper flow with a
      !! return flow uniformly distributed though the remaining
      !! near-surface layers"; MOM6 is top-down, k=1=surface, so
      !! `k <= nk_linear` selects the SURFACE region).  Roundabout is
      !! bottom-up (k=1=bed, CLAUDE.md); this column's `kb = k - 1` is the
      !! layer BELOW the interface (bed side), so `kb <= nk_linear` would
      !! select the BED-most layers — porting MOM6's `nk_linear=1` as
      !! written would put its SURFACE return-flow region at the SEABED.
      !! The dead `nk_linear` field + its `else` branch here were removed
      !! (PR-8) rather than wired: with the field permanently 0 (the only
      !! value ever set), `kb > nk_linear` was always true for every valid
      !! `kb >= 1`, so the limited path below ran unconditionally and the
      !! deleted branch never executed.  Adding a correct surface linear
      !! return-flow region is a live, unrecorded physics divergence from
      !! MOM6's default GM — it changes GM answers and needs its own PR +
      !! an analytical test (Psi -> 0 through the top layer), not a
      !! same-PR wire-up.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: i_smax2, i4dt
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: bathy(nx, ny)
      real(wp), intent(in) :: slope_x(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: khth_u(nx + 1, ny)
      logical, intent(in) :: use_open
      real(wp), intent(in) :: open_u(nx + 1, ny, nz)
      real(wp), intent(inout) :: uhD(nx + 1, ny, nz)

      integer :: i, j, k, iw, ka, kb, ktop
      real(wp) :: havL(NZ_STACK_MAX), havR(NZ_STACK_MAX)
      real(wp) :: rsumL(NZ_STACK_MAX + 1), rsumR(NZ_STACK_MAX + 1)
      real(wp) :: eL(NZ_STACK_MAX + 1), eR(NZ_STACK_MAX + 1)
      logical :: ok(NZ_STACK_MAX)
      real(wp) :: uhtot, slope, s2r, sfn_unlim, sfn_safe, sfn_est, sfn_in_h
      real(wp) :: h_frac_d, uhd_k, kh

      do concurrent(j=1:ny, i=2:nx) &
         local(k, iw, ka, kb, ktop, havL, havR, rsumL, rsumR, eL, eR, ok, uhtot, slope, &
               s2r, sfn_unlim, sfn_safe, sfn_est, sfn_in_h, h_frac_d, uhd_k, kh)
         iw = i - 1
         kh = khth_u(i, j)
         ! Interface heights (bed-up from z = -D) of the two columns, for
         ! the bottom-blocking limiter (`gm_block_below_bed`).
         eL(1) = -bathy(iw, j)
         eR(1) = -bathy(i, j)
         do k = 1, nz
            eL(k + 1) = eL(k) + h_layer(iw, j, k)
            eR(k + 1) = eR(k) + h_layer(i, j, k)
         end do

         ! The face's OPEN column (all-true without closed faces).
         do k = 1, nz
            ok(k) = .true.
         end do
         if (use_open) then
            do k = 1, nz
               ok(k) = open_u(i, j, k) > 0.5_wp .and. &
                       rdb_vl_is_live(h_layer(iw, j, k)) .and. &
                       rdb_vl_is_live(h_layer(i, j, k))
            end do
         end if
         ktop = 0
         do k = nz, 1, -1
            if (ok(k)) then
               ktop = k
               exit
            end if
         end do

         ! Per-layer availability + cumulative rsum from the SURFACE down:
         ! rsum*(k) = Sum_{k'=k}^{nz} h_avail(k') (mass above interface k).
         ! A layer outside the open column has nothing available here.
         rsumL(nz + 1) = 0.0_wp
         rsumR(nz + 1) = 0.0_wp
         do k = nz, 1, -1
            if (ok(k)) then
               havL(k) = max(i4dt*areaT(iw, j)*(h_layer(iw, j, k) - H_VANISHED), 0.0_wp)
               havR(k) = max(i4dt*areaT(i, j)*(h_layer(i, j, k) - H_VANISHED), 0.0_wp)
            else
               havL(k) = 0.0_wp
               havR(k) = 0.0_wp
            end if
            rsumL(k) = rsumL(k + 1) + havL(k)
            rsumR(k) = rsumR(k + 1) + havR(k)
         end do

         ! Sweep interior interfaces bed-most (Kr=2) -> surface-most (Kr=nz).
         uhtot = 0.0_wp
         uhD(i, j, nz) = 0.0_wp
         do k = 2, nz       ! k is the interface index Kr
            ka = k           ! layer above interface (surface side)
            kb = k - 1       ! layer below interface (bed side) -> uhD(kb)
            if (.not. ok(kb) .or. kb >= ktop) then
               ! Outside the open column, or its top layer (closed below):
               ! no transport, Sfn carried unchanged across it.
               uhD(i, j, kb) = 0.0_wp
               cycle
            end if
            slope = slope_x(i, j, k)
            ! NaN-safe: a non-finite slope must not reach the min/max
            ! limiters below (they launder NaN into a bound, CLAUDE.md).
            if (.not. ieee_is_finite(slope)) slope = 0.0_wp
            s2r = slope*slope*i_smax2
            sfn_unlim = gm_block_below_bed(-(kh*dy_cu(i, j))*slope, &
                                           eL(k), eL(k - 1), eR(1), eR(k), eR(k - 1), eL(1))
            if (uhtot <= 0.0_wp) then
               h_frac_d = gm_h_frac(havL(kb), rsumL(kb))
            else
               h_frac_d = gm_h_frac(havR(kb), rsumR(kb))
            end if
            sfn_safe = uhtot*(1.0_wp - h_frac_d)
            sfn_est = (sfn_unlim + s2r*sfn_safe)/(1.0_wp + s2r)
            ! Mass above interface Kr bounds the streamfunction.
            sfn_in_h = min(max(sfn_est, -rsumL(ka)), rsumR(ka))
            uhd_k = max(min(sfn_in_h - uhtot, havL(kb)), -havR(kb))
            uhD(i, j, kb) = uhd_k
            uhtot = uhtot + uhd_k
         end do
         ! Top BC (Sfn=0 above the topmost open layer, `nz` without closed
         ! faces): close the column so Sum=0.  No open layer ⇒ all zero.
         if (ktop > 0) uhD(i, j, ktop) = -uhtot
      end do
   end subroutine gm_column_x