coriolis_adv_compute_tendencies_hk Subroutine

   pure subroutine coriolis_adv_compute_tendencies_hk(grid, metrics, this, ms, u, v, h)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Per-layer PV-conserving Coriolis + horizontal-advection tendency
      !! in the Arakawa-Hsu (1990) form ("HK correction").  The wider
      !! 3-corner PV stencil at each face suppresses the spurious
      !! Hollingsworth-Källén instability that biases the simpler
      !! Sadourny 2-corner form at eddy-resolving resolutions.
      !!
      !! Algorithm — 6 passes per call:
      !!   1. Pass 1: relative vorticity ζ at interior corners (same
      !!      stencil as the Sadourny multilayer kernel) into `q_corner`.
      !!      Wall corners get ζ = 0 (free-slip BC).
      !!   2. Pass 2: per-mass PV `q = (f + ζ) / h_at_corner` rewritten
      !!      into `q_corner`.  `h_at_corner` is the AREA-WEIGHTED 4-cell
      !!      mean (each T thickness weighted by its `areaT`, normalised
      !!      by the summed areas), with `min/max` clamps so wall corners
      !!      collapse onto the available cells.  Reduces to the plain
      !!      4-cell mean on uniform Cartesian (equal areas).
      !!   3. Pass 3: per-face mass fluxes — `mass_flux_u = u_face_x ·
      !!      h_at_u_face` (face-averaged thickness) and `mass_flux_v`
      !!      symmetrically.  Wall faces fall through to the single
      !!      available cell (velocities are zero there anyway).
      !!   4. Pass 4: KE at cell centres — identical to Sadourny.
      !!   5. Pass 5: u-face tendency
      !!        CAu(i,j) = a · mass_flux_v(i,   j+1)   (NE)
      !!                 + b · mass_flux_v(i-1, j+1)   (NW)
      !!                 + c · mass_flux_v(i-1, j)     (SW)
      !!                 + d · mass_flux_v(i,   j)     (SE)
      !!                 - grad_KE_x
      !!      where each coefficient combines the face's two end
      !!      corners + one diagonal corner:
      !!        a = (q(i,j+1) + q(i+1,j+1) + q(i,j))   / 12
      !!        b = (q(i,j+1) + q(i-1,j+1) + q(i,j))   / 12
      !!        c = (q(i,j+1) + q(i-1,j)   + q(i,j))   / 12
      !!        d = (q(i,j+1) + q(i+1,j)   + q(i,j))   / 12
      !!   6. Pass 6: v-face tendency — symmetric construction; overall
      !!      minus sign on the q-stencil sum (Coriolis on v is `-f·u`):
      !!        CAv(i,j) = -[a' · mass_flux_u(i+1, j)
      !!                   + b' · mass_flux_u(i,   j)
      !!                   + c' · mass_flux_u(i,   j-1)
      !!                   + d' · mass_flux_u(i+1, j-1)]
      !!                   - grad_KE_y
      !!
      !! Reduction property (uniform h, uniform v): each coefficient
      !! evaluates to `q/4`, so the sum of 4 mass-flux terms is `q · vh`
      !! and the kernel collapses to the Sadourny `(f+ζ)·v` form
      !! bit-identically.  This is the basis for the regression test.
      !!
      !! Under `&vcoord_nml zfixed_closed_faces` (`metrics%use_closed_faces`)
      !! Passes 5/6 run a PAIR-FLOORED twin: every PV in a pair coefficient
      !! is evaluated with a corner thickness of at least half the larger of
      !! the pair's two face thicknesses (`hk_pair_coef`).  Without it the
      !! "cross" pairs — a corner PV times a transport whose far cell lies
      !! outside that corner — carry an unbounded `h_face/h_corner`, which a
      !! z-level staircase (a live partial bottom cell as thin as
      !! `H_VANISHED` against a full-depth neighbour) turns into a runaway
      !! (1-degree Southern Ocean: NaN at step 11).  The floor keeps the
      !! pair coefficients symmetric, so the energy-conserving antisymmetry
      !! is kept, and is inactive wherever no cell outweighs the other three
      !! of its corner.  The same twin runs with the knob OFF on the
      !! coordinates that lay static bed fillers (`hk_pair_floor`): there
      !! the open step face pairs a live cell with a `1e-4 m` filler, and
      !! without the floor the matrix staircase leaves the remap a negative
      !! thickness at step 2 (compat matrix `vcoord=zstar` x
      !! `coriolis=sadourny_hk`).  Every other coordinate ⇒ the original
      !! passes, bit-identical.
      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)

      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, q_NE, q_NW, q_SE, q_SW
      real(wp) :: a_NE, b_NW, c_SW, d_SE
      real(wp) :: ke_grad_x, ke_grad_y
      real(wp) :: h_S, h_N, h_W, h_E, h_NE, h_NW, h_SE, h_SW
         !! Closed-face branch only: corner thicknesses of the stencil.
      real(wp) :: h_u, h_v, hv_NE, hv_NW, hv_SW, hv_SE, hu_NE, hu_NW, hu_SW, hu_SE
         !! Closed-face branch only: face thicknesses of the stencil.
      real(wp) :: ns
      real(wp), parameter :: C1_12 = 1.0_wp/12.0_wp
      real(wp), parameter :: H_MIN_PV = 1.0e-12_wp
         !! Floor for the corner-h divide; vanishing-layer columns
         !! get q ≈ (f+ζ)/H_MIN_PV which is large but finite — paired
         !! with `mass_flux ≈ 0` at the same column so the product
         !! decays cleanly toward zero rather than blowing up.

      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)

      ! ---- 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 as `zeta_corner` 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 AREA-WEIGHTED mean of the four surrounding
      ! T cells, with min/max wall clamps — wall corners collapse onto
      ! the available cells (single cell at the four grid corners;
      ! two-cell mean along an edge).  Area weighting (each T thickness
      ! weighted by its own areaT, normalised by the summed areas) is the
      ! conservative corner thickness on curvilinear grids where
      ! areaT varies cell-to-cell; on uniform Cartesian every areaT is
      ! equal so it reduces to the plain 4-cell mean (same value).
      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)
         ! Mask-weighted corner area (C1, MOM6 `Area_h = mask2dT·areaT`):
         ! a blocked T-column contributes zero area + zero thickness, so
         ! `h_corner` is the wet-column mean only.  All-wet ⇒ `wet_T≡1` ⇒
         ! plain area-weighted mean (bit-identical).  `hm_den` floored at
         ! H_DIV_EPS for a fully-land corner (zeta=0 there, mass fluxes 0).
         !
         ! This is the MOM6 SADOURNY75_ENERGY (transport/energy-form) land
         ! treatment, NOT a divergence from it.  MOM6 masks corner area
         ! UNCONDITIONALLY at land (`Area_h = mask2dT·areaT`, zero on a land
         ! T-column) then forms `q = abs_vort·Area_q/(hArea_q + vol_neglect)`
         ! with `Area_q = Σ Area_h` over the four corners — algebraically
         ! `abs_vort·Σ(wet·areaT) / Σ(wet·areaT·h)`, IDENTICAL to our
         ! `(f+ζ)/h_corner` with `h_corner = Σ(wet·areaT·h)/Σ(wet·areaT)`.
         ! MOM6's *additional* `Area_h` area-mirroring across a boundary is
         ! OBC-SEGMENT-ONLY (inside `if (associated(OBC))`), not a land-coast
         ! rule — it does not apply here.  No HK-specific corner liberty.
         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
         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 do

      ! ---- 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) ----
      ! The PV/advection transports must carry the SAME narrowed face
      ! width continuity uses, or the two mass-flux definitions disagree.
      ! Host-side gate => no kernel launch and no textual change to the
      ! passes above when the knob is off (byte-identical).
      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 (&vcoord_nml zfixed_closed_faces) ----
      ! The SAME argument as the porous pass above, and a SEPARATE factor:
      ! the transport-Coriolis PV flux must see the same per-layer walls
      ! continuity does, or a closed layer contributes transport the mass
      ! budget never moved.  See the composition rule on
      ! `ocean_metrics_t%open_v`.
      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

      ! ---- 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

      if (metrics%use_closed_faces .or. this%hk_pair_floor) then
         ! ---- Passes 5f/6f: the same stencil, PAIR-FLOORED PV ----
         ! `&vcoord_nml zfixed_closed_faces`, or an OPEN staircase of static
         ! bed fillers (`hk_pair_floor`: z_fixed / zstar / zstar_full);
         ! otherwise the ELSE branch, the original Passes 5/6, textually
         ! untouched ⇒ bit-identical.  Inline,
         ! not a call: a host-gated call handing the tendency buffers to another
         ! procedure pessimises every loop of this routine on nvfortran.
         ! Each pair coefficient sums three corner PVs; the "cross" ones meet a
         ! transport whose far cell lies outside the corner, so `h_face/h_corner`
         ! is unbounded there (a thin live partial cell next to a full level).
         ! `hk_pair_coef` re-evaluates each PV at a corner thickness of at least
         ! half the pair's larger face thickness — the bound `sadourny_energy`
         ! has by construction — and the floor belongs to the PAIR, so the
         ! u- and v-tendencies share the coefficient (HK energy antisymmetry).
         ! See the routine docstring and `hk_pair_coef`.
         ! ---- Pass 5f: u-face ----
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(q_S, q_N, q_NE, q_NW, q_SE, q_SW, &
                  a_NE, b_NW, c_SW, d_SE, ke_grad_x, &
                  h_S, h_N, h_NE, h_NW, h_SE, h_SW, h_u, &
                  hv_NE, hv_NW, hv_SW, hv_SE)
            q_S = this%q_corner%data(i, j, k)
            q_N = this%q_corner%data(i, j + 1, k)
            q_NE = this%q_corner%data(i + 1, j + 1, k)
            q_NW = this%q_corner%data(i - 1, j + 1, k)
            q_SE = this%q_corner%data(i + 1, j, k)
            q_SW = this%q_corner%data(i - 1, j, k)
            h_S = hk_corner_h(i, j, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_N = hk_corner_h(i, j + 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_NE = hk_corner_h(i + 1, j + 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_NW = hk_corner_h(i - 1, j + 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_SE = hk_corner_h(i + 1, j, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_SW = hk_corner_h(i - 1, j, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            ! Face thicknesses exactly as Pass 3 builds the transports
            ! (one-sided at the array edge: 0.5·(a+a) == a).
            h_u = 0.5_wp*(h(i - 1, j, k) + h(i, j, k))
            hv_NE = 0.5_wp*(h(i, j, k) + h(i, min(ny, j + 1), k))
            hv_NW = 0.5_wp*(h(i - 1, j, k) + h(i - 1, min(ny, j + 1), k))
            hv_SW = 0.5_wp*(h(i - 1, max(1, j - 1), k) + h(i - 1, j, k))
            hv_SE = 0.5_wp*(h(i, max(1, j - 1), k) + h(i, j, k))
            a_NE = hk_pair_coef(q_N, h_N, q_NE, h_NE, q_S, h_S, max(h_u, hv_NE))
            b_NW = hk_pair_coef(q_N, h_N, q_NW, h_NW, q_S, h_S, max(h_u, hv_NW))
            c_SW = hk_pair_coef(q_N, h_N, q_SW, h_SW, q_S, h_S, max(h_u, hv_SW))
            d_SE = hk_pair_coef(q_N, h_N, q_SE, h_SE, q_S, h_S, max(h_u, hv_SE))
            ke_grad_x = (this%ke_centre%data(i, j, k) - &
                         this%ke_centre%data(i - 1, j, k))*metrics%idxCu(i, j)
            this%pv_flux_x%data(i, j, k) = &
               (a_NE*this%mass_flux_v%data(i, j + 1, k) + &
                b_NW*this%mass_flux_v%data(i - 1, j + 1, k) + &
                c_SW*this%mass_flux_v%data(i - 1, j, k) + &
                d_SE*this%mass_flux_v%data(i, j, k))*metrics%idxCu(i, j) - ke_grad_x
         end do
         ! ---- Pass 6f: v-face ----
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(q_W, q_E, q_NE, q_NW, q_SE, q_SW, &
                  a_NE, b_NW, c_SW, d_SE, ke_grad_y, &
                  h_W, h_E, h_NE, h_NW, h_SE, h_SW, h_v, &
                  hu_NE, hu_NW, hu_SW, hu_SE)
            q_W = this%q_corner%data(i, j, k)
            q_E = this%q_corner%data(i + 1, j, k)
            q_NE = this%q_corner%data(i + 1, j + 1, k)
            q_NW = this%q_corner%data(i, j + 1, k)
            q_SE = this%q_corner%data(i + 1, j - 1, k)
            q_SW = this%q_corner%data(i, j - 1, k)
            h_W = hk_corner_h(i, j, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_E = hk_corner_h(i + 1, j, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_NE = hk_corner_h(i + 1, j + 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_NW = hk_corner_h(i, j + 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_SE = hk_corner_h(i + 1, j - 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_SW = hk_corner_h(i, j - 1, k, nx, ny, nz, h, metrics%wet_T, metrics%areaT)
            h_v = 0.5_wp*(h(i, j - 1, k) + h(i, j, k))
            hu_NE = 0.5_wp*(h(i, j, k) + h(min(nx, i + 1), j, k))
            hu_NW = 0.5_wp*(h(max(1, i - 1), j, k) + h(i, j, k))
            hu_SW = 0.5_wp*(h(max(1, i - 1), j - 1, k) + h(i, j - 1, k))
            hu_SE = 0.5_wp*(h(i, j - 1, k) + h(min(nx, i + 1), j - 1, k))
            a_NE = hk_pair_coef(q_W, h_W, q_NE, h_NE, q_E, h_E, max(h_v, hu_NE))
            b_NW = hk_pair_coef(q_W, h_W, q_NW, h_NW, q_E, h_E, max(h_v, hu_NW))
            c_SW = hk_pair_coef(q_W, h_W, q_SW, h_SW, q_E, h_E, max(h_v, hu_SW))
            d_SE = hk_pair_coef(q_W, h_W, q_SE, h_SE, q_E, h_E, max(h_v, hu_SE))
            ke_grad_y = (this%ke_centre%data(i, j, k) - &
                         this%ke_centre%data(i, j - 1, k))*metrics%idyCv(i, j)
            this%pv_flux_y%data(i, j, k) = &
               -(a_NE*this%mass_flux_u%data(i + 1, j, k) + &
                 b_NW*this%mass_flux_u%data(i, j, k) + &
                 c_SW*this%mass_flux_u%data(i, j - 1, k) + &
                 d_SE*this%mass_flux_u%data(i + 1, j - 1, k))*metrics%idyCv(i, j) - ke_grad_y
         end do
      else
         ! ---- Pass 5: u-face HK tendency ----
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(q_S, q_N, q_NE, q_NW, q_SE, q_SW, &
                  a_NE, b_NW, c_SW, d_SE, ke_grad_x)
            q_S = this%q_corner%data(i, j, k)
            q_N = this%q_corner%data(i, j + 1, k)
            q_NE = this%q_corner%data(i + 1, j + 1, k)
            q_NW = this%q_corner%data(i - 1, j + 1, k)
            q_SE = this%q_corner%data(i + 1, j, k)
            q_SW = this%q_corner%data(i - 1, j, k)
            a_NE = (q_N + q_NE + q_S)*C1_12
            b_NW = (q_N + q_NW + q_S)*C1_12
            c_SW = (q_N + q_SW + q_S)*C1_12
            d_SE = (q_N + q_SE + q_S)*C1_12
            ke_grad_x = (this%ke_centre%data(i, j, k) - &
                         this%ke_centre%data(i - 1, j, k))*metrics%idxCu(i, j)
            ! q·vh sum is a transport-weighted PV flux (m³/s); the u-face
            ! IdxCu closes it to a per-length acceleration (= /dx on uniform).
            this%pv_flux_x%data(i, j, k) = &
               (a_NE*this%mass_flux_v%data(i, j + 1, k) + &
                b_NW*this%mass_flux_v%data(i - 1, j + 1, k) + &
                c_SW*this%mass_flux_v%data(i - 1, j, k) + &
                d_SE*this%mass_flux_v%data(i, j, k))*metrics%idxCu(i, j) - ke_grad_x
         end do

         ! ---- Pass 6: v-face HK tendency ----
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(q_W, q_E, q_NE, q_NW, q_SE, q_SW, &
                  a_NE, b_NW, c_SW, d_SE, ke_grad_y)
            q_W = this%q_corner%data(i, j, k)
            q_E = this%q_corner%data(i + 1, j, k)
            q_NE = this%q_corner%data(i + 1, j + 1, k)
            q_NW = this%q_corner%data(i, j + 1, k)
            q_SE = this%q_corner%data(i + 1, j - 1, k)
            q_SW = this%q_corner%data(i, j - 1, k)
            a_NE = (q_W + q_NE + q_E)*C1_12
            b_NW = (q_W + q_NW + q_E)*C1_12
            c_SW = (q_W + q_SW + q_E)*C1_12
            d_SE = (q_W + q_SE + q_E)*C1_12
            ke_grad_y = (this%ke_centre%data(i, j, k) - &
                         this%ke_centre%data(i, j - 1, k))*metrics%idyCv(i, j)
            this%pv_flux_y%data(i, j, k) = &
               -(a_NE*this%mass_flux_u%data(i + 1, j, k) + &
                 b_NW*this%mass_flux_u%data(i, j, k) + &
                 c_SW*this%mass_flux_u%data(i, j - 1, k) + &
                 d_SE*this%mass_flux_u%data(i + 1, j - 1, k))*metrics%idyCv(i, j) - ke_grad_y
         end do
      end if
      ! Array-edge faces carry no tendency (both branches).
      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
      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_hk