pure subroutine ocean_lateral_mix_compute_leith_biharm(grid, metrics, this, ms)
!! Public only for the unit-test suite; ignore in production code.
!! Populate `nu4_face_x`/`nu4_face_y` (m⁴/s) with the 2-D Leith
!! biharmonic viscosity
!! A_4(face) = clamp(C_lb · grid_sp⁶ · inv_PI6 · |∇²ζ|,
!! nu4_bg, nu4_max)
!! where ζ is C-grid corner relative vorticity, `∇²ζ` its 5-point
!! corner Laplacian, `grid_sp⁶ = grid_sp_h2³`, and `inv_PI6 = (1/π)⁶`.
!! Per-face |∇²ζ| is the mean of the two adjacent corner Laplacians.
!! Wall faces get `nu4_bg`. Leith (1968); Griffies & Hallberg (2000).
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
integer :: i, j, k, nx, ny, nz
real(wp) :: c_leith_bi_local, nu4_bg_local, nu4_max_local, ns, inv_pi6
real(wp) :: dx2, dy2, grid_sp_h2, grid_sp6, leith_bi_scale
real(wp) :: del2_a, del2_b, del2_face, A_raw
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 (see Leith/Smag
! siblings) — `do concurrent` bodies see plain real(wp).
nu4_bg_local = this%nu4_bg
nu4_max_local = this%nu4_max
c_leith_bi_local = this%c_leith_bi
ns = merge(1.0_wp, 0.0_wp, this%no_slip) ! free-slip(0)/no-slip(1) selector (C1)
inv_pi6 = (1.0_wp/PI)**6
! ---- Pass 1: relative vorticity at SW corners, per layer ----
! Identical circulation/area form to compute_leith (consistent
! with the Coriolis kernel's converted zeta). Outer-most corners
! stay zero (closed-wall convention) so the corner Laplacian below
! sees a finite neighbourhood at the first interior corners.
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_4 at u-faces (i-1/2, j) ----
! ∇²ζ at the two adjacent corners Bu(i,j) (SW) and Bu(i,j+1) (NW)
! averaged onto the face; |∇²ζ| scaled by C_lb·grid_sp⁶·inv_PI6.
! Corner Laplacian needs j∈[2,ny-1] (j-1/j+1 in range for the NW
! corner at j+1); wall-adjacent rows are filled by row-copy below.
do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) &
local(dx2, dy2, grid_sp_h2, grid_sp6, leith_bi_scale, &
del2_a, del2_b, del2_face, A_raw)
dx2 = metrics%dx2h(i, j)
dy2 = metrics%dy2h(i, j)
grid_sp_h2 = (2.0_wp*dx2*dy2)/(dx2 + dy2)
grid_sp6 = grid_sp_h2*grid_sp_h2*grid_sp_h2
leith_bi_scale = c_leith_bi_local*grid_sp6*inv_pi6
! ∇²ζ at SW corner Bu(i,j): 5-point corner Laplacian.
del2_a = (this%vort_corner%data(i + 1, j, k) - &
2.0_wp*this%vort_corner%data(i, j, k) + &
this%vort_corner%data(i - 1, j, k))/metrics%dx2q(i, j) + &
(this%vort_corner%data(i, j + 1, k) - &
2.0_wp*this%vort_corner%data(i, j, k) + &
this%vort_corner%data(i, j - 1, k))/metrics%dy2q(i, j)
! ∇²ζ at NW corner Bu(i,j+1).
del2_b = (this%vort_corner%data(i + 1, j + 1, k) - &
2.0_wp*this%vort_corner%data(i, j + 1, k) + &
this%vort_corner%data(i - 1, j + 1, k))/metrics%dx2q(i, j + 1) + &
(this%vort_corner%data(i, j + 2, k) - &
2.0_wp*this%vort_corner%data(i, j + 1, k) + &
this%vort_corner%data(i, j, k))/metrics%dy2q(i, j + 1)
del2_face = 0.5_wp*(del2_a + del2_b)
A_raw = leith_bi_scale*abs(del2_face)
this%nu4_face_x(i, j, k) = min(nu4_max_local, max(nu4_bg_local, A_raw))
end do
do concurrent(k=1:nz, i=2:nx)
this%nu4_face_x(i, 1, k) = this%nu4_face_x(i, 2, k)
this%nu4_face_x(i, ny, k) = this%nu4_face_x(i, ny - 1, k)
end do
do concurrent(k=1:nz, j=1:ny)
this%nu4_face_x(1, j, k) = nu4_bg_local
this%nu4_face_x(nx + 1, j, k) = nu4_bg_local
end do
! ---- Pass 2b: A_4 at v-faces (i, j-1/2) ----
! Adjacent corners Bu(i,j) (SW) and Bu(i+1,j) (SE) averaged onto
! the face. i∈[2,nx-1] keeps i-1/i+1 in range for the SE corner.
do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) &
local(dx2, dy2, grid_sp_h2, grid_sp6, leith_bi_scale, &
del2_a, del2_b, del2_face, A_raw)
dx2 = metrics%dx2h(i, j)
dy2 = metrics%dy2h(i, j)
grid_sp_h2 = (2.0_wp*dx2*dy2)/(dx2 + dy2)
grid_sp6 = grid_sp_h2*grid_sp_h2*grid_sp_h2
leith_bi_scale = c_leith_bi_local*grid_sp6*inv_pi6
! ∇²ζ at SW corner Bu(i,j).
del2_a = (this%vort_corner%data(i + 1, j, k) - &
2.0_wp*this%vort_corner%data(i, j, k) + &
this%vort_corner%data(i - 1, j, k))/metrics%dx2q(i, j) + &
(this%vort_corner%data(i, j + 1, k) - &
2.0_wp*this%vort_corner%data(i, j, k) + &
this%vort_corner%data(i, j - 1, k))/metrics%dy2q(i, j)
! ∇²ζ at SE corner Bu(i+1,j).
del2_b = (this%vort_corner%data(i + 2, j, k) - &
2.0_wp*this%vort_corner%data(i + 1, j, k) + &
this%vort_corner%data(i, j, k))/metrics%dx2q(i + 1, j) + &
(this%vort_corner%data(i + 1, j + 1, k) - &
2.0_wp*this%vort_corner%data(i + 1, j, k) + &
this%vort_corner%data(i + 1, j - 1, k))/metrics%dy2q(i + 1, j)
del2_face = 0.5_wp*(del2_a + del2_b)
A_raw = leith_bi_scale*abs(del2_face)
this%nu4_face_y(i, j, k) = min(nu4_max_local, max(nu4_bg_local, A_raw))
end do
do concurrent(k=1:nz, j=2:ny)
this%nu4_face_y(1, j, k) = this%nu4_face_y(2, j, k)
this%nu4_face_y(nx, j, k) = this%nu4_face_y(nx - 1, j, k)
end do
do concurrent(k=1:nz, i=1:nx)
this%nu4_face_y(i, 1, k) = nu4_bg_local
this%nu4_face_y(i, ny + 1, k) = nu4_bg_local
end do
end subroutine ocean_lateral_mix_compute_leith_biharm