fill_vorticity_z_impl Subroutine

   pure subroutine fill_vorticity_z_impl(u_face, v_face, dyCv, dxCu, wet_q, wet_T, &
                                         iareaBu, no_slip, nz_ml, buf)
      ! assumed-shape-ok: diag fill — fires once per output frame (cadence-bounded);
      ! face-sized u_face/v_face and the Cv/Cu/Bu metrics carry the natural
      ! (nx[+1], ny[+1]) C-grid stagger; size() min-clips at call.
      !! Cell-centred relative vorticity that MIRRORS the dynamics' own
      !! corner ζ bit-for-bit — the shared `rdb_rvc_zeta_corner` body
      !! (`src/shared_module_utilities/rdb_rel_vort_corner.inc`), the same
      !! circulation-form formula `coriolis_adv_compute_tendencies` Pass 1
      !! evaluates inline for its `q_corner` scratch
      !! (`src/core/ocean/kernels/coriolis_adv/rdb_coriolis_adv.F90`), incl.
      !! the C1 slip factor: `no_slip=.false.` (default, free-slip) masks a
      !! land corner (`wet_q=0`) to zero rel-vort; `.true.` (no-slip) gives
      !! it the `2-wet_q` image-vorticity value instead.
      !!
      !! The 4 corners of each T-cell (SW=(i,j), SE=(i+1,j), NW=(i,j+1),
      !! NE=(i+1,j+1) in the corner-storage convention `q_corner(i,j)` sits
      !! at (i-1/2,j-1/2)) are averaged onto the cell centre with a plain
      !! `0.25*(sw+se+nw+ne)` — the SAME fixed-count average the dynamics
      !! itself uses to interpolate `q_corner` onto a u/v-face
      !! (`zeta_at_u = 0.5*(q_corner(i,j)+q_corner(i,j+1))`), just over 4
      !! corners instead of 2.  It is deliberately NOT re-weighted by
      !! `wet_q`: each corner's OWN value already carries the correct
      !! land/coast physics via the C1 factor above (free-slip -> exactly
      !! 0 at a land corner, so it dilutes a coastal cell's average toward
      !! 0 the same way the dynamics' own face-average does; no-slip -> the
      !! nonzero image-vorticity value, so the boundary condition actually
      !! reaches the diagnostic instead of being masked out a second time
      !! by a wet-only divisor). A land cell (`wet_T=0`) reports
      !! `DIAG_MISSING_VALUE`.
      real(wp), intent(in)    :: u_face(:, :, :), v_face(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      ! assumed-shape-ok: diag fill — cadence-bounded (once per output frame)
      real(wp), intent(in)    :: dyCv(:, :), dxCu(:, :), wet_q(:, :), wet_T(:, :)  ! assumed-shape-ok: diag fill
      real(wp), intent(in)    :: iareaBu(:, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      logical, intent(in)     :: no_slip
      integer, intent(in)    :: nz_ml
      real(wp), intent(inout) :: buf(:, :, :)  ! assumed-shape-ok: diag fill — cadence-bounded
      integer :: i, j, k, nx, ny, nz
      real(wp) :: ns, z_sw, z_se, z_nw, z_ne
      call zero3_impl(buf)
      nx = min(size(buf, 1), size(wet_T, 1))
      ny = min(size(buf, 2), size(wet_T, 2))
      nz = min(size(buf, 3), nz_ml)
      ns = merge(1.0_wp, 0.0_wp, no_slip)
      ! Interior cells only (i in [2, nx-1], j in [2, ny-1]), matching the
      ! predecessor stencil's rim convention: the ±1 corner reach needs
      ! i-1/i+1 and j-1/j+1 in bounds, so the outer rim stays at the zero
      ! seed exactly as before this fix (an orthogonal, pre-existing
      ! limitation of the diag's domain coverage, not part of the bug).
      do concurrent(k=1:nz, j=2:ny - 1, i=2:nx - 1) &
         local(z_sw, z_se, z_nw, z_ne)
         if (wet_T(i, j) > 0.5_wp) then
            z_sw = rdb_rvc_zeta_corner(v_face(i, j, k), v_face(i - 1, j, k), &
                                       dyCv(i, j), dyCv(i - 1, j), &
                                       u_face(i, j, k), u_face(i, j - 1, k), &
                                       dxCu(i, j), dxCu(i, j - 1), &
                                       iareaBu(i, j), wet_q(i, j), ns)
            z_se = rdb_rvc_zeta_corner(v_face(i + 1, j, k), v_face(i, j, k), &
                                       dyCv(i + 1, j), dyCv(i, j), &
                                       u_face(i + 1, j, k), u_face(i + 1, j - 1, k), &
                                       dxCu(i + 1, j), dxCu(i + 1, j - 1), &
                                       iareaBu(i + 1, j), wet_q(i + 1, j), ns)
            z_nw = rdb_rvc_zeta_corner(v_face(i, j + 1, k), v_face(i - 1, j + 1, k), &
                                       dyCv(i, j + 1), dyCv(i - 1, j + 1), &
                                       u_face(i, j + 1, k), u_face(i, j, k), &
                                       dxCu(i, j + 1), dxCu(i, j), &
                                       iareaBu(i, j + 1), wet_q(i, j + 1), ns)
            z_ne = rdb_rvc_zeta_corner(v_face(i + 1, j + 1, k), v_face(i, j + 1, k), &
                                       dyCv(i + 1, j + 1), dyCv(i, j + 1), &
                                       u_face(i + 1, j + 1, k), u_face(i + 1, j, k), &
                                       dxCu(i + 1, j + 1), dxCu(i + 1, j), &
                                       iareaBu(i + 1, j + 1), wet_q(i + 1, j + 1), ns)
            buf(i, j, k) = 0.25_wp*((z_sw + z_se) + (z_nw + z_ne))
         else
            buf(i, j, k) = DIAG_MISSING_VALUE
         end if
      end do
   end subroutine fill_vorticity_z_impl