compute_gtot_faces Subroutine

public pure subroutine compute_gtot_faces(grid, bt_work, ms, metrics)

Face-centred depth-weighted column averages of pbce (gtot_E/W/N/S). Wall cells fall back to pbce(:,:,nz). By construction Σ_k h_face(k)·(pbce(k) − gtot_face) = 0 per column, making the bc-PGF Δu correction depth-mean zero – UNDER THE SAME h_face/weight derive_bt_from_layers and apply_bt_correction’s folds use.

Root cause of the fix(bt) frhat-matrix mass leak (budget c055-class cells, eulerian_z + correction_bc_pgf, 2026-10-06): this routine hand-rolled h_face = 0.5*(h_L+h_R) – the plain arithmetic mean, UNCONDITIONALLY – while every other depth mean in this module (derive_bt_from_layers’s bt_ubt, the BT forcing’s F_bt_u/_fast, apply_bt_correction’s own folds) already dispatches on &ocean_bt_nml frhat_scheme via frhat_h_face_step, gated to metrics%use_closed_faces like every other call site (see derive_bt_from_layers’s matching comment). gtot_face’s “depth-mean zero” property holds for ANY weight consistent with ITSELF (true by construction of a weighted mean), but apply_bt_correction’s do_bc_pgf fold adds Δu_bc directly, LAYER BY LAYER, to u_face_x_layer – so what cancels the barotropic mode’s own mass transport is Δu_bc’s depth mean UNDER THE WEIGHT THE FAST LOOP ACTUALLY INTEGRATED, not under some independent, always-arithmetic weight. Fixing ONLY this routine does not, by itself, close the c055 leak (isolated by bisection: derive_bt_from_layers’s own bt_ubt, gated above, is the dominant term – HYBRID’s shelf-test/harmonic suppression assumes the z_fixed/zstar FILLER convention, where a layer with h -> 0 genuinely carries no water; on every other family (hycom, eulerian_z, rho, sigma, lagrangian) a thin layer is REAL water merely scaled by a shallow column, and HYBRID suppresses it anyway, breaking exactly this mass consistency) – but this routine is its own latent instance of the same inconsistency class and would bite the moment correction_bc_pgf composes with a closed-faces z_fixed/zstar run (currently refused at configure, closed_faces_bc_pgf), so it is fixed here too, independently.

frhat_h_face_step is hand-inlined into all four face loops below rather than called – the SAME NVHPC -stdpar=gpu call-site miscompile apply_bt_correction’s visc_rem-weighted fold already documents (see that routine’s comment for the original compute-sanitizer evidence): the plain call version of this routine crashed with an “Invalid global read” inside this kernel’s launch (compute-sanitizer memcheck, 2026-10-06, gtot_faces_hybrid_gated_on_closed_faces / gtot_depth_mean_zero_invariant), reproduced with ms/bt_work unmapped (this routine’s normal calling convention – see apply_bt_correction’s matching comment) and not a bounds bug (i/j neighbours are always in range by construction). The inline copy is textually identical to rdb_frhat_face.inc’s body.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(barotropic_workstate_t), intent(inout) :: bt_work
type(multilayer_state_t), intent(in) :: ms
type(ocean_metrics_t), intent(in) :: metrics

REQUIRED, like every other frhat call site’s metrics – see derive_bt_from_layers’s docstring for why an optional here is the defect class, not a defect.

Calls

proc~~compute_gtot_faces~~CallsGraph proc~compute_gtot_faces compute_gtot_faces local local proc~compute_gtot_faces->local

Called by

proc~~compute_gtot_faces~~CalledByGraph proc~compute_gtot_faces compute_gtot_faces proc~run_stage_split run_stage_split proc~run_stage_split->proc~compute_gtot_faces proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: d_shallow_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: e_cur_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: e_prev
real(kind=wp), private :: h_arith_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: h_face
real(kind=wp), private :: h_harm_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: h_sum
real(kind=wp), private :: hl_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: hr_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.
real(kind=wp), private :: href_l
real(kind=wp), private :: href_r
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: nx
integer, private :: ny
integer, private :: nz
real(kind=wp), private :: p_sum
integer, private :: scheme
real(kind=wp), private :: wt_arith_loc

Hand-inlined frhat_h_face_step locals – see the docstring above for why this routine inlines rather than calls.

Source Code

   pure subroutine compute_gtot_faces(grid, bt_work, ms, metrics)
      !! Face-centred depth-weighted column averages of `pbce` (gtot_E/W/N/S).
      !! Wall cells fall back to `pbce(:,:,nz)`. By construction
      !! Σ_k h_face(k)·(pbce(k) − gtot_face) = 0 per column, making the bc-PGF
      !! Δu correction depth-mean zero -- UNDER THE SAME h_face/weight
      !! `derive_bt_from_layers` and `apply_bt_correction`'s folds use.
      !!
      !! Root cause of the fix(bt) frhat-matrix mass leak (budget
      !! c055-class cells, `eulerian_z` + `correction_bc_pgf`, 2026-10-06):
      !! this routine hand-rolled `h_face = 0.5*(h_L+h_R)` -- the plain
      !! arithmetic mean, UNCONDITIONALLY -- while every other depth mean
      !! in this module (`derive_bt_from_layers`'s `bt_ubt`, the BT
      !! forcing's `F_bt_u/_fast`, `apply_bt_correction`'s own folds)
      !! already dispatches on `&ocean_bt_nml frhat_scheme` via
      !! `frhat_h_face_step`, gated to `metrics%use_closed_faces` like
      !! every other call site (see `derive_bt_from_layers`'s matching
      !! comment). `gtot_face`'s "depth-mean zero" property holds for ANY
      !! weight consistent with ITSELF (true by construction of a
      !! weighted mean), but `apply_bt_correction`'s `do_bc_pgf` fold adds
      !! `Δu_bc` directly, LAYER BY LAYER, to `u_face_x_layer` -- so what
      !! cancels the barotropic mode's own mass transport is Δu_bc's depth
      !! mean UNDER THE WEIGHT THE FAST LOOP ACTUALLY INTEGRATED, not under
      !! some independent, always-arithmetic weight. Fixing ONLY this
      !! routine does not, by itself, close the c055 leak (isolated by
      !! bisection: `derive_bt_from_layers`'s own `bt_ubt`, gated above, is
      !! the dominant term -- HYBRID's shelf-test/harmonic suppression
      !! assumes the z_fixed/zstar FILLER convention, where a layer with
      !! h -> 0 genuinely carries no water; on every other family
      !! (hycom, eulerian_z, rho, sigma, lagrangian) a thin layer is REAL
      !! water merely scaled by a shallow column, and HYBRID suppresses it
      !! anyway, breaking exactly this mass consistency) -- but this
      !! routine is its own latent instance of the same inconsistency
      !! class and would bite the moment `correction_bc_pgf` composes with
      !! a closed-faces z_fixed/zstar run (currently refused at configure,
      !! `closed_faces_bc_pgf`), so it is fixed here too, independently.
      !!
      !! `frhat_h_face_step` is hand-inlined into all four face loops below
      !! rather than `call`ed -- the SAME NVHPC `-stdpar=gpu` call-site
      !! miscompile `apply_bt_correction`'s visc_rem-weighted fold already
      !! documents (see that routine's comment for the original
      !! compute-sanitizer evidence): the plain `call` version of this
      !! routine crashed with an "Invalid __global__ read" inside this
      !! kernel's launch (compute-sanitizer memcheck, 2026-10-06,
      !! `gtot_faces_hybrid_gated_on_closed_faces` /
      !! `gtot_depth_mean_zero_invariant`), reproduced with `ms`/`bt_work`
      !! unmapped (this routine's normal calling convention -- see
      !! `apply_bt_correction`'s matching comment) and not a bounds bug
      !! (`i`/`j` neighbours are always in range by construction). The
      !! inline copy is textually identical to `rdb_frhat_face.inc`'s body.
      type(hgrid_t), intent(in) :: grid
      type(barotropic_workstate_t), intent(inout) :: bt_work
      type(multilayer_state_t), intent(in) :: ms
      type(ocean_metrics_t), intent(in) :: metrics
         !! REQUIRED, like every other frhat call site's `metrics` --
         !! see `derive_bt_from_layers`'s docstring for why an optional
         !! here is the defect class, not a defect.

      integer :: i, j, k, nx, ny, nz, scheme
      real(wp) :: h_face, h_sum, p_sum, href_l, href_r, e_prev
      real(wp) :: hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc
         !! Hand-inlined `frhat_h_face_step` locals -- see the docstring
         !! above for why this routine inlines rather than calls.

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      ! Gated exactly like every other frhat call site in this module --
      ! see `derive_bt_from_layers`'s matching comment.
      scheme = merge(bt_work%frhat_scheme, FRHAT_ARITHMETIC, metrics%use_closed_faces)

      do concurrent(j=1:ny, i=1:nx - 1) &
         local(k, h_face, h_sum, p_sum, href_l, href_r, e_prev, &
               hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
         href_l = bt_work%bt_H_ref(i, j)
         href_r = bt_work%bt_H_ref(i + 1, j)
         e_prev = -0.5_wp*(href_l + href_r)
         h_sum = 0.0_wp
         p_sum = 0.0_wp
         do k = 1, nz
            ! `frhat_h_face_step` hand-inlined -- see this routine's
            ! docstring above (NVHPC call-site miscompile).
            hl_loc = ms%h_layer(i, j, k)
            hr_loc = ms%h_layer(i + 1, j, k)
            h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
            if (scheme == FRHAT_HYBRID) then
               d_shallow_loc = -min(href_l, href_r)
               e_cur_loc = e_prev + h_arith_loc
               ! vanished-ok: hand-inlined frhat_h_face_step copy (NVHPC call-site miscompile, see docstring above)
               if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                  h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
               else if (e_prev >= d_shallow_loc) then
                  h_face = h_arith_loc
               else
                  h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  if (e_cur_loc <= d_shallow_loc) then
                     h_face = h_harm_loc
                  else
                     wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                     h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                  end if
               end if
               e_prev = e_cur_loc
            else
               h_face = h_arith_loc
            end if
            h_sum = h_sum + h_face
            p_sum = p_sum + h_face*bt_work%pbce(i, j, k)
         end do
         if (h_sum > 0.0_wp) then
            bt_work%gtot_E(i, j) = p_sum/h_sum
         else
            bt_work%gtot_E(i, j) = bt_work%pbce(i, j, nz)
         end if
      end do
      do concurrent(j=1:ny)
         bt_work%gtot_E(nx, j) = bt_work%pbce(nx, j, nz)
      end do

      do concurrent(j=1:ny, i=2:nx) &
         local(k, h_face, h_sum, p_sum, href_l, href_r, e_prev, &
               hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
         href_l = bt_work%bt_H_ref(i - 1, j)
         href_r = bt_work%bt_H_ref(i, j)
         e_prev = -0.5_wp*(href_l + href_r)
         h_sum = 0.0_wp
         p_sum = 0.0_wp
         do k = 1, nz
            ! See the matching comment in the E-face loop above.
            hl_loc = ms%h_layer(i - 1, j, k)
            hr_loc = ms%h_layer(i, j, k)
            h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
            if (scheme == FRHAT_HYBRID) then
               d_shallow_loc = -min(href_l, href_r)
               e_cur_loc = e_prev + h_arith_loc
               ! vanished-ok: hand-inlined frhat_h_face_step copy, same reason as the E-face loop above
               if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                  h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
               else if (e_prev >= d_shallow_loc) then
                  h_face = h_arith_loc
               else
                  h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  if (e_cur_loc <= d_shallow_loc) then
                     h_face = h_harm_loc
                  else
                     wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                     h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                  end if
               end if
               e_prev = e_cur_loc
            else
               h_face = h_arith_loc
            end if
            h_sum = h_sum + h_face
            p_sum = p_sum + h_face*bt_work%pbce(i, j, k)
         end do
         if (h_sum > 0.0_wp) then
            bt_work%gtot_W(i, j) = p_sum/h_sum
         else
            bt_work%gtot_W(i, j) = bt_work%pbce(i, j, nz)
         end if
      end do
      do concurrent(j=1:ny)
         bt_work%gtot_W(1, j) = bt_work%pbce(1, j, nz)
      end do

      do concurrent(j=1:ny - 1, i=1:nx) &
         local(k, h_face, h_sum, p_sum, href_l, href_r, e_prev, &
               hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
         href_l = bt_work%bt_H_ref(i, j)
         href_r = bt_work%bt_H_ref(i, j + 1)
         e_prev = -0.5_wp*(href_l + href_r)
         h_sum = 0.0_wp
         p_sum = 0.0_wp
         do k = 1, nz
            ! See the matching comment in the E-face loop above.
            hl_loc = ms%h_layer(i, j, k)
            hr_loc = ms%h_layer(i, j + 1, k)
            h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
            if (scheme == FRHAT_HYBRID) then
               d_shallow_loc = -min(href_l, href_r)
               e_cur_loc = e_prev + h_arith_loc
               ! vanished-ok: hand-inlined frhat_h_face_step copy, same reason as the E-face loop above
               if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                  h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
               else if (e_prev >= d_shallow_loc) then
                  h_face = h_arith_loc
               else
                  h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  if (e_cur_loc <= d_shallow_loc) then
                     h_face = h_harm_loc
                  else
                     wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                     h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                  end if
               end if
               e_prev = e_cur_loc
            else
               h_face = h_arith_loc
            end if
            h_sum = h_sum + h_face
            p_sum = p_sum + h_face*bt_work%pbce(i, j, k)
         end do
         if (h_sum > 0.0_wp) then
            bt_work%gtot_N(i, j) = p_sum/h_sum
         else
            bt_work%gtot_N(i, j) = bt_work%pbce(i, j, nz)
         end if
      end do
      do concurrent(i=1:nx)
         bt_work%gtot_N(i, ny) = bt_work%pbce(i, ny, nz)
      end do

      do concurrent(j=2:ny, i=1:nx) &
         local(k, h_face, h_sum, p_sum, href_l, href_r, e_prev, &
               hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
         href_l = bt_work%bt_H_ref(i, j - 1)
         href_r = bt_work%bt_H_ref(i, j)
         e_prev = -0.5_wp*(href_l + href_r)
         h_sum = 0.0_wp
         p_sum = 0.0_wp
         do k = 1, nz
            ! See the matching comment in the E-face loop above.
            hl_loc = ms%h_layer(i, j - 1, k)
            hr_loc = ms%h_layer(i, j, k)
            h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
            if (scheme == FRHAT_HYBRID) then
               d_shallow_loc = -min(href_l, href_r)
               e_cur_loc = e_prev + h_arith_loc
               ! vanished-ok: hand-inlined frhat_h_face_step copy, same reason as the E-face loop above
               if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                  h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
               else if (e_prev >= d_shallow_loc) then
                  h_face = h_arith_loc
               else
                  h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  if (e_cur_loc <= d_shallow_loc) then
                     h_face = h_harm_loc
                  else
                     wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                     h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                  end if
               end if
               e_prev = e_cur_loc
            else
               h_face = h_arith_loc
            end if
            h_sum = h_sum + h_face
            p_sum = p_sum + h_face*bt_work%pbce(i, j, k)
         end do
         if (h_sum > 0.0_wp) then
            bt_work%gtot_S(i, j) = p_sum/h_sum
         else
            bt_work%gtot_S(i, j) = bt_work%pbce(i, j, nz)
         end if
      end do
      do concurrent(i=1:nx)
         bt_work%gtot_S(i, 1) = bt_work%pbce(i, 1, nz)
      end do
   end subroutine compute_gtot_faces