ocean_lateral_mix_compute_leith Subroutine

   pure subroutine ocean_lateral_mix_compute_leith(grid, metrics, this, ms, &
                                                   res_fn_u, res_fn_v)
      !! Public only for the unit-test suite; ignore in production code.
      !! Populate `ah_face_x`/`ah_face_y` (m^2/s) with the Leith viscosity
      !!     A_h(face) = max(ah_bg, min(ah_max, (C_L · dx)^3 · |∇ζ|))
      !! where ζ is relative vorticity at C-grid corners (pass 1) and |∇ζ|
      !! the 2D gradient magnitude at each face (pass 2).  Wall faces get
      !! the background viscosity.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_lateral_mix_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in), optional :: res_fn_u(grid%nx_total + 1, grid%ny_total)
         !! VarMix resolution function at u-faces (nondim, [0,1]).  When
         !! present AND `this%resoln_scaled_visc`, scales `A_h` before clamp.
      real(wp), intent(in), optional :: res_fn_v(grid%nx_total, grid%ny_total + 1)
         !! VarMix resolution function at v-faces.

      integer :: i, j, k, nx, ny, nz
      real(wp) :: c_leith_local, ah_bg_local, ah_max_local
      real(wp) :: dzeta_dx, dzeta_dy, grad_mag, A_raw, leith_scale
      real(wp) :: ns
      logical :: do_resoln

      if (.not. this%is_init) return
      if (.not. allocated(ms%u_face_x_layer)) return
      if (.not. allocated(ms%v_face_y_layer)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      ! Hoist scalar fields off `this` to locals — `do concurrent`
      ! bodies see plain real(wp) instead of a derived-type deref.
      ah_bg_local = this%ah_bg
      ah_max_local = this%ah_max
      c_leith_local = this%c_leith
      ns = merge(1.0_wp, 0.0_wp, this%no_slip)   ! free-slip(0)/no-slip(1) selector (C1)
      ! Resolution-function scaling active only when the knob is on AND the
      ! VarMix face fields were supplied (Gap 1, Hallberg 2013); gated INSIDE
      ! the face loops so the optional is referenced only when present.
      do_resoln = this%resoln_scaled_visc .and. present(res_fn_u) .and. &
                  present(res_fn_v)

      ! Leith dimensionful prefactor: (C_L · L_grid)^3 with the grid
      ! scale `L_grid = sqrt(dxT·dyT)` evaluated per cell (design §2;
      ! = dx on uniform square metrics, so bit-reducing).  |∇ζ| has
      ! units 1/(m·s); A ~ L³·|∇ζ| is the right order for mesoscale
      ! closures.  The per-face scale below picks the adjacent T cell.

      ! ---- Pass 1: relative vorticity at SW corners, per layer ----
      ! Circulation/area form (consistent with the Coriolis kernel's
      ! converted zeta): ζ = (Δ(v·dyCv) − Δ(u·dxCu))·iareaBu.  Reduces
      ! to (Δv)/dx − (Δu)/dy on uniform square metrics.  Outer-most
      ! corners stay zero (closed-wall convention).
      do concurrent(k=1:nz, j=2:ny, i=2:nx)
         this%vort_corner%data(i, j, k) = &
            ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
            ((ms%v_face_y_layer(i, j, k)*metrics%dyCv(i, j) - &
              ms%v_face_y_layer(i - 1, j, k)*metrics%dyCv(i - 1, j)) - &
             (ms%u_face_x_layer(i, j, k)*metrics%dxCu(i, j) - &
              ms%u_face_x_layer(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%vort_corner%data(1, j, k) = 0.0_wp
         this%vort_corner%data(nx + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx + 1)
         this%vort_corner%data(i, 1, k) = 0.0_wp
         this%vort_corner%data(i, ny + 1, k) = 0.0_wp
      end do

      ! ---- Pass 2a: A_h at u-faces (i-1/2, j) ----
      ! Adjacent corners: (i, j) at (i-1/2, j-1/2) and (i, j+1) at
      ! (i-1/2, j+1/2).  Across-face corners (one cell west/east):
      ! (i-1, j), (i-1, j+1), (i+1, j), (i+1, j+1).
      !
      ! ∂ζ/∂y at u-face: (ζ(i, j+1) - ζ(i, j))·idyCu  (along the face).
      ! ∂ζ/∂x at u-face: idxCu·[(ζ_E_S + ζ_E_N) - (ζ_W_S + ζ_W_N)]/4
      ! leith_scale = (C_L·sqrt(dxT·dyT))³ at the adjacent T cell
      ! (design §2; = (C_L·dx)³ on uniform square metrics).
      do concurrent(k=1:nz, j=1:ny, i=2:nx) &
         local(dzeta_dx, dzeta_dy, grad_mag, A_raw, leith_scale)
         leith_scale = (c_leith_local*sqrt(metrics%dxT(i, j)*metrics%dyT(i, j)))**3
         dzeta_dy = (this%vort_corner%data(i, j + 1, k) - &
                     this%vort_corner%data(i, j, k))*metrics%idyCu(i, j)
         dzeta_dx = 0.25_wp*((this%vort_corner%data(i + 1, j, k) - &
                              this%vort_corner%data(i - 1, j, k)) + &
                             (this%vort_corner%data(i + 1, j + 1, k) - &
                              this%vort_corner%data(i - 1, j + 1, k)))*metrics%idxCu(i, j)
         grad_mag = sqrt(dzeta_dx*dzeta_dx + dzeta_dy*dzeta_dy)
         ! Resolution scaling applied BEFORE the clamp, as one assignment to
         ! the `local()` var `A_raw` per `do_resoln` branch (assigning a
         ! `local()` var once on each path; a conditional REASSIGN of a
         ! `do concurrent local()` var miscompiles on gfortran 15.1).
         if (do_resoln) then
            A_raw = leith_scale*grad_mag*res_fn_u(i, j)
         else
            A_raw = leith_scale*grad_mag
         end if
         this%ah_face_x(i, j, k) = min(ah_max_local, max(ah_bg_local, A_raw))
      end do
      ! Wall faces and i=1 edge: use background viscosity.
      do concurrent(k=1:nz, j=1:ny)
         this%ah_face_x(1, j, k) = ah_bg_local
         this%ah_face_x(nx + 1, j, k) = ah_bg_local
      end do

      ! ---- Pass 2b: A_h at v-faces (i, j-1/2) ----
      do concurrent(k=1:nz, j=2:ny, i=1:nx) &
         local(dzeta_dx, dzeta_dy, grad_mag, A_raw, leith_scale)
         leith_scale = (c_leith_local*sqrt(metrics%dxT(i, j)*metrics%dyT(i, j)))**3
         dzeta_dx = (this%vort_corner%data(i + 1, j, k) - &
                     this%vort_corner%data(i, j, k))*metrics%idxCv(i, j)
         dzeta_dy = 0.25_wp*((this%vort_corner%data(i, j + 1, k) - &
                              this%vort_corner%data(i, j - 1, k)) + &
                             (this%vort_corner%data(i + 1, j + 1, k) - &
                              this%vort_corner%data(i + 1, j - 1, k)))*metrics%idyCv(i, j)
         grad_mag = sqrt(dzeta_dx*dzeta_dx + dzeta_dy*dzeta_dy)
         if (do_resoln) then
            A_raw = leith_scale*grad_mag*res_fn_v(i, j)
         else
            A_raw = leith_scale*grad_mag
         end if
         this%ah_face_y(i, j, k) = min(ah_max_local, max(ah_bg_local, A_raw))
      end do
      ! Wall faces: background viscosity.
      do concurrent(k=1:nz, i=1:nx)
         this%ah_face_y(i, 1, k) = ah_bg_local
         this%ah_face_y(i, ny + 1, k) = ah_bg_local
      end do
   end subroutine ocean_lateral_mix_compute_leith