coriolis_adv_compute_tendencies_sadourny_energy Subroutine

   pure subroutine coriolis_adv_compute_tendencies_sadourny_energy(grid, metrics, this, ms, &
                                                                   u, v, h, use_state_fluxes)
      !! Faithful MOM6 SADOURNY75_ENERGY (Sadourny 1975 energy-conserving)
      !! per-layer Coriolis + horizontal-advection tendency.  This is the
      !! TRANSPORT form: the absolute-vorticity flux is the potential
      !! vorticity `q = (f + ζ)/h_at_corner` times the layer MASS TRANSPORT
      !! (`vh`/`uh`), so the discrete Coriolis term produces zero net domain
      !! kinetic energy (energy-conserving).  The default enstrophy form
      !! (`_sadourny`, `(f+ζ)·v`) only matches this under uniform thickness.
      !!
      !!   CAu(i,j) = 0.25·( q_N·(vh_NW+vh_NE) + q_S·(vh_SW+vh_SE) )·idxCu − ∂x KE
      !!   CAv(i,j) = −0.25·( q_E·(uh_NE+uh_SE) + q_W·(uh_NW+uh_SW) )·idyCv − ∂y KE
      !!
      !! with q_N/q_S the north/south corner PVs of the u-face (q_E/q_W the
      !! east/west corner PVs of the v-face) and vh/uh the per-face mass
      !! transports.  Reduction: uniform thickness ⇒ each q·(Σvh) collapses
      !! to `(f+ζ)·v`, bit-identical to the enstrophy form (regression test).
      !!
      !! Passes 1-4 (ζ at corners → PV `q = (f+ζ)/h_corner` with the
      !! wet-area-weighted corner thickness → mass transports → centre KE)
      !! mirror `coriolis_adv_compute_tendencies_hk` exactly; only the final
      !! corner→face stencil differs (Sadourny 2-corner vs HK 12-point).
      !! KEEP THE PREP PASSES IN SYNC with `_hk` (shared-prep extraction is
      !! tracked as a follow-up cleanup).
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(coriolis_adv_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in) :: u(grid%nx_total + 1, grid%ny_total, ms%nz_ml)
         !! Face-velocity / thickness source arrays (outer-shim; the
         !! dispatcher forwards either the prognostic components or the
         !! `u_av` time-mean family under `split_scheme = "pred_corr"`).
      real(wp), intent(in) :: v(grid%nx_total, grid%ny_total + 1, ms%nz_ml)
      real(wp), intent(in) :: h(grid%nx_total, grid%ny_total, ms%nz_ml)
      logical, intent(in), optional :: use_state_fluxes
         !! Mass-consistent CorAdCalc (MOM6 parity): fill the transport
         !! buffers from `ms%mass_flux_*_layer` — the continuity solve's
         !! renormalised uh/vh (same `u·h_face·dy_cu` m³/s convention,
         !! same shape, physical walls already zeroed) — instead of
         !! recomputing from the u/h source arrays.  In the pred_corr
         !! CORRECTOR those are the PREDICTOR chain's fluxes, i.e. exactly
         !! the transport field that produced the `u_av` evaluation state,
         !! so the q·vh product is energy-consistent on rim columns where
         !! the renorm/wall-zero and the naive `u·h_face` recompute
         !! disagree.  Absent / `.false.` ⇒ bit-identical recompute path.

      integer :: i, j, k, nx, ny, nz, nu, nv
      real(wp) :: zeta_corner, h_corner, h_face
      real(wp) :: aSW, aSE, aNW, aNE, hm_num, hm_den
      integer :: iw, ie, js, jn
      real(wp) :: q_S, q_N, q_W, q_E, ke_grad_x, ke_grad_y
      real(wp) :: ns
      logical :: usf
      real(wp) :: pv_part, av_a, av_b, fv1, fv2, fv3, fv4
      logical :: do_bound, use_mom6_ch
      real(wp), parameter :: H_MIN_PV = CORIOLIS_H_MIN_PV
         !! Floor for the corner-h divide; vanishing-layer columns get
         !! q ≈ (f+ζ)/H_MIN_PV (large but finite) paired with mass_flux ≈ 0
         !! at the same column so the product decays toward zero.  (MOM6
         !! instead floors the denominator, Area_q/(hArea_q+vol_neglect),
         !! so its q→0 as h→0; the forms differ only in the vanishing-layer
         !! limit, immaterial for non-vanishing envelopes.)

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nu = size(u, 1)
      nv = size(v, 2)
      ns = merge(1.0_wp, 0.0_wp, this%no_slip)
      usf = .false.
      if (present(use_state_fluxes)) usf = use_state_fluxes
      ! BOUND_CORIOLIS host flag (read once; the clamp branch is untaken when
      ! off ⇒ bit-identical).  Energy scheme only (config-guaranteed).
      do_bound = this%bound_coriolis
      ! corner_h host flag: MOM6 area-weighted PV corner thickness.  Read once;
      ! the mom6 branch is untaken (and the div-then-cap path bit-identical to
      ! pre-knob) when off.  Energy scheme only (config-guaranteed).
      use_mom6_ch = this%corner_h_variant == CORNER_H_MOM6_AREA

      ! ---- Pass 1: relative vorticity at corners (circulation/area) ----
      ! Inline twin of `rdb_rvc_zeta_corner` (shared_module_utilities/
      ! rdb_rel_vort_corner.inc, read by the `vorticity_z` diag): keep in step.
      ! Slip factor masks the rel-vort at land corners (C1); Pass 2 reads
      ! this back and adds the UNMASKED planetary f.
      do concurrent(k=1:nz, j=2:ny, i=2:nx)
         this%q_corner%data(i, j, k) = &
            ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
            ((v(i, j, k)*metrics%dyCv(i, j) - &
              v(i - 1, j, k)*metrics%dyCv(i - 1, j)) - &
             (u(i, j, k)*metrics%dxCu(i, j) - &
              u(i, j - 1, k)*metrics%dxCu(i, j - 1)))* &
            metrics%iareaBu(i, j)
      end do
      do concurrent(k=1:nz, j=1:ny + 1)
         this%q_corner%data(1, j, k) = 0.0_wp
         this%q_corner%data(nx + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx + 1)
         this%q_corner%data(i, 1, k) = 0.0_wp
         this%q_corner%data(i, ny + 1, k) = 0.0_wp
      end do

      ! ---- Pass 2: PV q = (f + ζ) / h_at_corner ----
      ! h_at_corner is the wet-area-weighted 4-cell mean (MOM6 `Area_h =
      ! mask2dT·areaT`; land cells contribute zero area + thickness).  q =
      ! abs_vort·Area_q/hArea_q = abs_vort/h_corner.  All-wet ⇒ plain mean.
      ! `hm_num` ≡ MOM6 `hArea_q` (Σ area·h), `hm_den` ≡ MOM6 `Area_q` (Σ area)
      ! — the two constructions share these EXACTLY; they differ only in the
      ! vanishing-thickness guard (cell_mean caps h_corner at H_MIN_PV;
      ! mom6_area's PV_VOL_NEGLECT is pure 1/0 armor, matching MOM6).
      do concurrent(k=1:nz, j=1:ny + 1, i=1:nx + 1) &
         local(zeta_corner, h_corner, iw, ie, js, jn, &
               aSW, aSE, aNW, aNE, hm_num, hm_den)
         zeta_corner = this%q_corner%data(i, j, k)
         iw = max(1, i - 1)
         ie = min(nx, i)
         js = max(1, j - 1)
         jn = min(ny, j)
         aSW = metrics%wet_T(iw, js)*metrics%areaT(iw, js)
         aSE = metrics%wet_T(ie, js)*metrics%areaT(ie, js)
         aNW = metrics%wet_T(iw, jn)*metrics%areaT(iw, jn)
         aNE = metrics%wet_T(ie, jn)*metrics%areaT(ie, jn)
         hm_num = aSW*h(iw, js, k) + aSE*h(ie, js, k) + &
                  aNW*h(iw, jn, k) + aNE*h(ie, jn, k)
         hm_den = aSW + aSE + aNW + aNE
         if (use_mom6_ch) then
            ! MOM6 area form: q = abs_vort·Area_q/(hArea_q + vol_neglect).
            ! No thickness cap — vol_neglect is pure 1/0 armor.
            this%q_corner%data(i, j, k) = (this%f_corner(i, j) + zeta_corner)* &
                                          hm_den/(hm_num + PV_VOL_NEGLECT)
         else
            h_corner = hm_num/max(hm_den, H_DIV_EPS)
            h_corner = max(h_corner, H_MIN_PV)
            this%q_corner%data(i, j, k) = (this%f_corner(i, j) + zeta_corner)/h_corner
         end if
      end do

      ! ---- Pass 3a/3b: face transports uh / vh ----
      ! Two sources, same convention (u·h_face·dy_cu, m³/s):
      !   recompute (default) — self-contained `u·h_face` from the u/h
      !     source arrays; order-independent of continuity.
      !   state fluxes (`use_state_fluxes`) — copy the continuity
      !     solve's renormalised `ms%mass_flux_*_layer` (MOM6
      !     mass-consistent CorAdCalc; pred_corr-corrector stage only,
      !     where they still hold the predictor chain's fluxes — the
      !     transport field that produced the `u_av` evaluation state).
      if (usf) then
         do concurrent(k=1:nz, j=1:ny, i=1:nu)
            this%mass_flux_u%data(i, j, k) = ms%mass_flux_x_layer(i, j, k)
         end do
         do concurrent(k=1:nz, j=1:nv, i=1:nx)
            this%mass_flux_v%data(i, j, k) = ms%mass_flux_y_layer(i, j, k)
         end do
      else
         ! Pass 3a: u-face transport uh = u·h·dy_cu
         do concurrent(k=1:nz, j=1:ny, i=1:nu) local(h_face)
            if (i == 1) then
               h_face = h(1, j, k)
            else if (i == nu) then
               h_face = h(nx, j, k)
            else
               h_face = 0.5_wp*(h(i - 1, j, k) + h(i, j, k))
            end if
            this%mass_flux_u%data(i, j, k) = u(i, j, k)*h_face*metrics%dy_cu(i, j)
         end do

         ! Pass 3b: v-face transport vh = v·h·dx_cv
         do concurrent(k=1:nz, j=1:nv, i=1:nx) local(h_face)
            if (j == 1) then
               h_face = h(i, 1, k)
            else if (j == nv) then
               h_face = h(i, ny, k)
            else
               h_face = 0.5_wp*(h(i, j - 1, k) + h(i, j, k))
            end if
            this%mass_flux_v%data(i, j, k) = v(i, j, k)*h_face*metrics%dx_cv(i, j)
         end do

         ! ---- Porous barriers (Adcroft 2013) ----
         ! INSIDE the `else` only: the `usf` branch above copies
         ! continuity's mass fluxes, which are ALREADY narrowed, so
         ! applying the fraction again would square it.
         if (metrics%use_porous) then
            call porous_narrow_3d(nu, ny, nz, metrics%por_face_area_u, &
                                  this%mass_flux_u%data)
            call porous_narrow_3d(nx, nv, nz, metrics%por_face_area_v, &
                                  this%mass_flux_v%data)
         end if

         ! ---- z-level closed faces ----
         ! INSIDE the same `else` and for the same reason: the `usf`
         ! branch copies continuity's fluxes, which the mask already
         ! closed.
         if (metrics%use_closed_faces) then
            call porous_narrow_3d(nu, ny, nz, metrics%open_u, &
                                  this%mass_flux_u%data)
            call porous_narrow_3d(nx, nv, nz, metrics%open_v, &
                                  this%mass_flux_v%data)
         end if
      end if

      ! ---- Pass 4: KE at cell centres (area-weighted) ----
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         this%ke_centre%data(i, j, k) = 0.25_wp*metrics%iareaT(i, j)*( &
                                        metrics%areaCu(i, j)*u(i, j, k)**2 + &
                                        metrics%areaCu(i + 1, j)*u(i + 1, j, k)**2 + &
                                        metrics%areaCv(i, j)*v(i, j, k)**2 + &
                                        metrics%areaCv(i, j + 1)*v(i, j + 1, k)**2)
      end do

      ! ---- Pass 5: u-face energy tendency (2-corner Sadourny, q·vh) ----
      ! q·vh sum is a transport-weighted PV flux (m³/s); idxCu closes it to
      ! a per-length acceleration (= /dx on uniform).
      do concurrent(k=1:nz, j=1:ny, i=2:nx) &
         local(q_S, q_N, ke_grad_x, pv_part, av_a, av_b, fv1, fv2, fv3, fv4)
         q_S = this%q_corner%data(i, j, k)
         q_N = this%q_corner%data(i, j + 1, k)
         ke_grad_x = (this%ke_centre%data(i, j, k) - &
                      this%ke_centre%data(i - 1, j, k))*metrics%idxCu(i, j)
         pv_part = 0.25_wp*( &
                   q_N*(this%mass_flux_v%data(i - 1, j + 1, k) + &
                        this%mass_flux_v%data(i, j + 1, k)) &
                   + q_S*(this%mass_flux_v%data(i - 1, j, k) + &
                          this%mass_flux_v%data(i, j, k)))* &
                   metrics%idxCu(i, j)
         if (do_bound) then
            ! BOUND_CORIOLIS (MOM6): clamp the PV flux into the range
            ! of the four neighbour (f+ζ)·v velocity-form estimates — north
            ! corner (i,j+1) × v(i-1/i,j+1); south corner (i,j) × v(i-1/i,j) —
            ! BEFORE subtracting the KE gradient.  abs_vort = q·h_corner.
            av_b = corner_abs_vort(i, j + 1, k, nx, ny, nz, q_N, h, &
                                   metrics%wet_T, metrics%areaT, use_mom6_ch)
            av_a = corner_abs_vort(i, j, k, nx, ny, nz, q_S, h, &
                                   metrics%wet_T, metrics%areaT, use_mom6_ch)
            fv1 = av_b*v(i - 1, j + 1, k)
            fv2 = av_b*v(i, j + 1, k)
            fv3 = av_a*v(i - 1, j, k)
            fv4 = av_a*v(i, j, k)
            pv_part = min(pv_part, max(max(fv1, fv2), max(fv3, fv4)))
            pv_part = max(pv_part, min(min(fv1, fv2), min(fv3, fv4)))
         end if
         this%pv_flux_x%data(i, j, k) = pv_part - ke_grad_x
      end do
      do concurrent(k=1:nz, j=1:ny)
         this%pv_flux_x%data(1, j, k) = 0.0_wp
         this%pv_flux_x%data(nx + 1, j, k) = 0.0_wp
      end do

      ! ---- Pass 6: v-face energy tendency (2-corner Sadourny, −q·uh) ----
      do concurrent(k=1:nz, j=2:ny, i=1:nx) &
         local(q_W, q_E, ke_grad_y, pv_part, av_a, av_b, fv1, fv2, fv3, fv4)
         q_W = this%q_corner%data(i, j, k)
         q_E = this%q_corner%data(i + 1, j, k)
         ke_grad_y = (this%ke_centre%data(i, j, k) - &
                      this%ke_centre%data(i, j - 1, k))*metrics%idyCv(i, j)
         pv_part = -0.25_wp*( &
                   q_E*(this%mass_flux_u%data(i + 1, j, k) + &
                        this%mass_flux_u%data(i + 1, j - 1, k)) &
                   + q_W*(this%mass_flux_u%data(i, j, k) + &
                          this%mass_flux_u%data(i, j - 1, k)))* &
                   metrics%idyCv(i, j)
         if (do_bound) then
            ! BOUND_CORIOLIS (MOM6): clamp into the four neighbour
            ! −(f+ζ)·u estimates — east corner (i+1,j) × u(i+1,j/j-1); west
            ! corner (i,j) × u(i,j/j-1) — BEFORE subtracting the KE gradient.
            av_b = corner_abs_vort(i + 1, j, k, nx, ny, nz, q_E, h, &
                                   metrics%wet_T, metrics%areaT, use_mom6_ch)
            av_a = corner_abs_vort(i, j, k, nx, ny, nz, q_W, h, &
                                   metrics%wet_T, metrics%areaT, use_mom6_ch)
            fv1 = -av_b*u(i + 1, j, k)
            fv2 = -av_b*u(i + 1, j - 1, k)
            fv3 = -av_a*u(i, j, k)
            fv4 = -av_a*u(i, j - 1, k)
            pv_part = min(pv_part, max(max(fv1, fv2), max(fv3, fv4)))
            pv_part = max(pv_part, min(min(fv1, fv2), min(fv3, fv4)))
         end if
         this%pv_flux_y%data(i, j, k) = pv_part - ke_grad_y
      end do
      do concurrent(k=1:nz, i=1:nx)
         this%pv_flux_y%data(i, 1, k) = 0.0_wp
         this%pv_flux_y%data(i, ny + 1, k) = 0.0_wp
      end do
   end subroutine coriolis_adv_compute_tendencies_sadourny_energy