subroutine varmix_compute_impl(nx, ny, nz, p, alpha, resoln_khth, &
resoln_khtr, interp_res, do_visbeck, s2max, &
khth, khtr, khth_cff, khtr_cff, khth_min, &
khth_max, khtr_min, khtr_max, &
f2_dx2_u, f2_dx2_v, beta_dx2_u, beta_dx2_v, &
l2_u, l2_v, cg1, h_layer, slope_x, slope_y, &
n2_u, n2_v, res_fn_u, res_fn_v, sn_u, sn_v, &
khth_u, khth_v, khtr_u, khtr_v)
!! Flat-impl VarMix kernel (explicit-shape; NVHPC descriptor-walk-free).
!! Three phases: (1) Res_fn at faces (cg1 interpolated to faces or the
!! centre-Res_fn averaged, per `interp_res`), (2) Eady SN at u/v faces
!! (thickness-weighted column reductions with the orthogonal slope folded
!! into S^2, scalar accumulators — SN is final, no SN_v combine), (3) the
!! assembly (Visbeck addend, Res_fn scale, clamp) into the KhTh/KhTr base.
integer, intent(in) :: nx, ny, nz, p
real(wp), intent(in) :: alpha, s2max
logical, intent(in) :: resoln_khth, resoln_khtr, interp_res, do_visbeck
real(wp), intent(in) :: khth, khtr, khth_cff, khtr_cff
real(wp), intent(in) :: khth_min, khth_max, khtr_min, khtr_max
real(wp), intent(in) :: f2_dx2_u(nx + 1, ny)
real(wp), intent(in) :: f2_dx2_v(nx, ny + 1)
real(wp), intent(in) :: beta_dx2_u(nx + 1, ny)
real(wp), intent(in) :: beta_dx2_v(nx, ny + 1)
real(wp), intent(in) :: l2_u(nx + 1, ny)
real(wp), intent(in) :: l2_v(nx, ny + 1)
real(wp), intent(in) :: cg1(nx, ny)
real(wp), intent(in) :: h_layer(nx, ny, nz)
real(wp), intent(in) :: slope_x(nx + 1, ny, nz + 1)
real(wp), intent(in) :: slope_y(nx, ny + 1, nz + 1)
real(wp), intent(in) :: n2_u(nx + 1, ny, nz + 1)
real(wp), intent(in) :: n2_v(nx, ny + 1, nz + 1)
real(wp), intent(out) :: res_fn_u(nx + 1, ny)
real(wp), intent(out) :: res_fn_v(nx, ny + 1)
real(wp), intent(out) :: sn_u(nx + 1, ny)
real(wp), intent(out) :: sn_v(nx, ny + 1)
real(wp), intent(out) :: khth_u(nx + 1, ny)
real(wp), intent(out) :: khth_v(nx, ny + 1)
real(wp), intent(out) :: khtr_u(nx + 1, ny)
real(wp), intent(out) :: khtr_v(nx, ny + 1)
integer :: i, j
real(wp) :: cg1u, cg1v
! ---- 1. Resolution function at faces. ----
do concurrent(j=1:ny, i=1:nx + 1) local(cg1u)
res_fn_u(i, j) = 0.0_wp
if (i >= 2 .and. i <= nx) then
if (interp_res) then
! Centre Res_fn (own + west) then 2-pt average.
res_fn_u(i, j) = 0.5_wp*( &
varmix_res_fn(f2_dx2_u(i, j), beta_dx2_u(i, j), &
cg1(i - 1, j), alpha, p) + &
varmix_res_fn(f2_dx2_u(i, j), beta_dx2_u(i, j), &
cg1(i, j), alpha, p))
else
cg1u = 0.5_wp*(cg1(i - 1, j) + cg1(i, j))
res_fn_u(i, j) = varmix_res_fn(f2_dx2_u(i, j), beta_dx2_u(i, j), &
cg1u, alpha, p)
end if
end if
end do
do concurrent(j=1:ny + 1, i=1:nx) local(cg1v)
res_fn_v(i, j) = 0.0_wp
if (j >= 2 .and. j <= ny) then
if (interp_res) then
res_fn_v(i, j) = 0.5_wp*( &
varmix_res_fn(f2_dx2_v(i, j), beta_dx2_v(i, j), &
cg1(i, j - 1), alpha, p) + &
varmix_res_fn(f2_dx2_v(i, j), beta_dx2_v(i, j), &
cg1(i, j), alpha, p))
else
cg1v = 0.5_wp*(cg1(i, j - 1) + cg1(i, j))
res_fn_v(i, j) = varmix_res_fn(f2_dx2_v(i, j), beta_dx2_v(i, j), &
cg1v, alpha, p)
end if
end if
end do
! ---- 2. Eady SN (thickness-weighted) at u/v faces. These store the
! FINAL SN: S^2 = slope_x^2 + (h-weighted 4-corner slope_y^2) at
! u (mirror at v), so the orthogonal slope is already included —
! SN_u/SN_v are complete after the thickness-weighted sum (MOM6
! calc_Visbeck_coeffs_old; no separate SN_v combine). ----
call varmix_sn_u(nx, ny, nz, do_visbeck, s2max, h_layer, slope_x, &
slope_y, n2_u, sn_u)
call varmix_sn_v(nx, ny, nz, do_visbeck, s2max, h_layer, slope_x, &
slope_y, n2_v, sn_v)
! ---- 3. Assembly (Visbeck addend, Res_fn scale, clamp) into the pre-CFL
! base KhTh/KhTr face fields, consuming the final SN_u/SN_v.
! Assembly: SN_u/SN_v ALREADY carry the orthogonal slope (h4-weighted
! into S^2 in varmix_sn_*, per MOM6 calc_Visbeck_coeffs_old), so they
! are the FINAL Eady growth rate — NO extra 4-corner SN_v combine (that
! belongs to the separate calc_Eady_growth_rate_2D path and here would
! double-count the orthogonal slope).
do concurrent(j=1:ny, i=1:nx + 1)
khth_u(i, j) = varmix_assemble(khth, khth_cff, l2_u(i, j), sn_u(i, j), &
res_fn_u(i, j), resoln_khth, khth_min, &
khth_max, do_visbeck)
khtr_u(i, j) = varmix_assemble(khtr, khtr_cff, l2_u(i, j), sn_u(i, j), &
res_fn_u(i, j), resoln_khtr, khtr_min, &
khtr_max, do_visbeck)
end do
do concurrent(j=1:ny + 1, i=1:nx)
khth_v(i, j) = varmix_assemble(khth, khth_cff, l2_v(i, j), sn_v(i, j), &
res_fn_v(i, j), resoln_khth, khth_min, &
khth_max, do_visbeck)
khtr_v(i, j) = varmix_assemble(khtr, khtr_cff, l2_v(i, j), sn_v(i, j), &
res_fn_v(i, j), resoln_khtr, khtr_min, &
khtr_max, do_visbeck)
end do
end subroutine varmix_compute_impl