pure subroutine ocean_lateral_mix_compute_smag_ah(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 biharmonic
!! Smagorinsky viscosity
!! A_4(face) = clamp(C_b · L⁴ · |D|, nu4_bg, nu4_max)
!! where `L² = 2·dx²·dy²/(dx²+dy²)` (harmonic mean of dx²,dy²) and
!! `|D|` is the strain-rate magnitude from `compute_smag`. Wall
!! faces get `nu4_bg`. `SMAG_BI_CONST` ≈ 0.015–0.06.
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) :: smag_bi_const_local, smag_bi_scale, nu4_bg_local, nu4_max_local
real(wp) :: dx2, dy2, grid_sp_h2
real(wp) :: D_T_W, D_T_E, D_T_S, D_T_N, D_T_face
real(wp) :: D_S_S, D_S_N, D_S_W, D_S_E, D_S_face
real(wp) :: strain_mag, A_raw, ns
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
ns = merge(1.0_wp, 0.0_wp, this%no_slip) ! free-slip(0)/no-slip(1) selector (C1)
! Biharmonic ν₄ is deliberately NOT resolution-scaled: the
! Hallberg (2013) resolution function suppresses only the
! scale-non-selective Laplacian, whereas the ∝k⁴ biharmonic already
! spares the resolved (large) scales and needs no suppression.
! Per-cell `C_b · (grid_sp_h2)^2` (MOM6 `Biharm_const_xx`), with
! grid_sp_h2 = 2·dx2h·dy2h/(dx2h+dy2h) the harmonic mean of the
! per-cell dxT²/dyT² (design §2; = the uniform value on square
! metrics, bit-reducing). Strain inverses per stagger as in
! `compute_smag`.
nu4_bg_local = this%nu4_bg
nu4_max_local = this%nu4_max
smag_bi_const_local = this%smag_bi_const
! ---- u-face viscosity (i-1/2, j) ----
do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) &
local(D_T_W, D_T_E, D_S_S, D_S_N, D_T_face, D_S_face, &
strain_mag, A_raw, dx2, dy2, grid_sp_h2, smag_bi_scale)
dx2 = metrics%dx2h(i, j)
dy2 = metrics%dy2h(i, j)
grid_sp_h2 = (2.0_wp*dx2*dy2)/(dx2 + dy2)
smag_bi_scale = smag_bi_const_local*(grid_sp_h2*grid_sp_h2)
D_T_W = (ms%u_face_x_layer(i, j, k) - ms%u_face_x_layer(i - 1, j, k))*metrics%idxT(i - 1, j) &
- (ms%v_face_y_layer(i - 1, j + 1, k) - ms%v_face_y_layer(i - 1, j, k))*metrics%idyT(i - 1, j)
D_T_E = (ms%u_face_x_layer(i + 1, j, k) - ms%u_face_x_layer(i, j, k))*metrics%idxT(i, j) &
- (ms%v_face_y_layer(i, j + 1, k) - ms%v_face_y_layer(i, j, k))*metrics%idyT(i, j)
D_T_face = 0.5_wp*(D_T_W + D_T_E)
! Corner D_S: ratio-bundle form (design §2). Bu(i,j) = SW corner
! of T(i,j); v at Cv(i,j)/Cv(i-1,j), u at Cu(i,j)/Cu(i,j-1).
! Corner shear strain sh_xy masked by slip factor (C1): free-slip
! ×wet_q / no-slip ×(2-wet_q). Bit-identical for all-wet.
D_S_S = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
(metrics%dy_dxBu(i, j)* &
(ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j) - &
ms%v_face_y_layer(i - 1, j, k)*metrics%idyCv(i - 1, j)) + &
metrics%dx_dyBu(i, j)* &
(ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j) - &
ms%u_face_x_layer(i, j - 1, k)*metrics%idxCu(i, j - 1)))
D_S_N = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j + 1) + 2.0_wp*ns)* &
(metrics%dy_dxBu(i, j + 1)* &
(ms%v_face_y_layer(i, j + 1, k)*metrics%idyCv(i, j + 1) - &
ms%v_face_y_layer(i - 1, j + 1, k)*metrics%idyCv(i - 1, j + 1)) + &
metrics%dx_dyBu(i, j + 1)* &
(ms%u_face_x_layer(i, j + 1, k)*metrics%idxCu(i, j + 1) - &
ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j)))
D_S_face = 0.5_wp*(D_S_S + D_S_N)
strain_mag = sqrt(D_T_face*D_T_face + D_S_face*D_S_face)
A_raw = smag_bi_scale*strain_mag
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
! ---- v-face viscosity (i, j-1/2) ----
do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) &
local(D_T_S, D_T_N, D_S_W, D_S_E, D_T_face, D_S_face, &
strain_mag, A_raw, dx2, dy2, grid_sp_h2, smag_bi_scale)
dx2 = metrics%dx2h(i, j)
dy2 = metrics%dy2h(i, j)
grid_sp_h2 = (2.0_wp*dx2*dy2)/(dx2 + dy2)
smag_bi_scale = smag_bi_const_local*(grid_sp_h2*grid_sp_h2)
D_T_S = (ms%u_face_x_layer(i + 1, j - 1, k) - ms%u_face_x_layer(i, j - 1, k))*metrics%idxT(i, j - 1) &
- (ms%v_face_y_layer(i, j, k) - ms%v_face_y_layer(i, j - 1, k))*metrics%idyT(i, j - 1)
D_T_N = (ms%u_face_x_layer(i + 1, j, k) - ms%u_face_x_layer(i, j, k))*metrics%idxT(i, j) &
- (ms%v_face_y_layer(i, j + 1, k) - ms%v_face_y_layer(i, j, k))*metrics%idyT(i, j)
D_T_face = 0.5_wp*(D_T_S + D_T_N)
! Corner D_S: ratio-bundle form (design §2). Bu(i,j) = SW corner
! of T(i,j); v at Cv(i,j)/Cv(i-1,j), u at Cu(i,j)/Cu(i,j-1).
! Corner shear strain sh_xy masked by slip factor (C1).
D_S_W = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
(metrics%dy_dxBu(i, j)* &
(ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j) - &
ms%v_face_y_layer(i - 1, j, k)*metrics%idyCv(i - 1, j)) + &
metrics%dx_dyBu(i, j)* &
(ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j) - &
ms%u_face_x_layer(i, j - 1, k)*metrics%idxCu(i, j - 1)))
D_S_E = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i + 1, j) + 2.0_wp*ns)* &
(metrics%dy_dxBu(i + 1, j)* &
(ms%v_face_y_layer(i + 1, j, k)*metrics%idyCv(i + 1, j) - &
ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j)) + &
metrics%dx_dyBu(i + 1, j)* &
(ms%u_face_x_layer(i + 1, j, k)*metrics%idxCu(i + 1, j) - &
ms%u_face_x_layer(i + 1, j - 1, k)*metrics%idxCu(i + 1, j - 1)))
D_S_face = 0.5_wp*(D_S_W + D_S_E)
strain_mag = sqrt(D_T_face*D_T_face + D_S_face*D_S_face)
A_raw = smag_bi_scale*strain_mag
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_smag_ah