hvisc_compute_scalar_impl Subroutine

   pure subroutine hvisc_compute_scalar_impl(u_face, v_face, du_visc, dv_visc, &
                                             nu_h, &
                                             dy_dxT, dx_dyBu, iareaCu, &
                                             dx_dyT, dy_dxBu, iareaCv, &
                                             bound_kh, bound_coef, dt, &
                                             idxCu, idyCu, idxCv, idyCv, wet_q, ns, &
                                             nx, ny, nz, open_u, open_v)
      !! Per-face metric Laplacian × scalar viscosity.  Used when no
      !! lateral-mix closure is active — falls back to constant
      !! `nu_h`.  When `nu_h = 0` the kernel still zeros all interior +
      !! boundary cells so the apply step sees a defined state.
      !! `bound_kh` engages the per-face harmonic CFL ceiling
      !! (`hvisc_kh_cfl_bound`); `.false.` ⇒ bit-identical.
      integer, intent(in)    :: nx, ny, nz
      real(wp), intent(in)    :: u_face(nx + 1, ny, nz), v_face(nx, ny + 1, nz)
      real(wp), intent(inout) :: du_visc(nx + 1, ny, nz), dv_visc(nx, ny + 1, nz)
      real(wp), intent(in)    :: nu_h
      real(wp), intent(in)    :: dy_dxT(nx, ny), dx_dyBu(nx + 1, ny + 1), iareaCu(nx + 1, ny)
      real(wp), intent(in)    :: dx_dyT(nx, ny), dy_dxBu(nx + 1, ny + 1), iareaCv(nx, ny + 1)
      logical, intent(in)     :: bound_kh
      real(wp), intent(in)    :: bound_coef, dt
      real(wp), intent(in)    :: idxCu(nx + 1, ny), idyCu(nx + 1, ny)
      real(wp), intent(in)    :: idxCv(nx, ny + 1), idyCv(nx, ny + 1)
      real(wp), intent(in)    :: wet_q(nx + 1, ny + 1)
         !! Bu-corner wet mask (`metrics%wet_q`).  Each corner (shear) flux of
         !! the velocity Laplacian is scaled by the slip factor
         !! `(1 - 2·ns)·wet_q + 2·ns` -- the same C1 factor the Smagorinsky
         !! strain, the Coriolis corner vorticity and the stress-tensor path
         !! use (MOM6 `sh_xy = mask2dBu·(dvdx+dudy)` free-slip, `(2-mask2dBu)`
         !! no-slip).  Free-slip: a land corner carries NO shear flux, so the
         !! zero stored at a land face never acts as a Dirichlet-0 wall.
         !! All-wet corner => factor 1 => bit-identical.
      real(wp), intent(in)    :: ns
         !! 1 = no-slip (`&ocean_hvisc_nml no_slip`), 0 = free-slip.
      real(wp), intent(in), optional :: open_u(nx + 1, ny, nz)
         !! Per-layer 0/1 u-face open mask
         !! (`&vcoord_nml zfixed_closed_faces`) — the FREE-SLIP closure
         !! of a z-level partial step.  See `hvisc_compute_face_impl`'s
         !! `open_u` docstring for the full argument; ABSENT (the default
         !! path) ⇒ the loops below are textually unchanged ⇒
         !! bit-identical.
      real(wp), intent(in), optional :: open_v(nx, ny + 1, nz)
         !! v-face twin.  Present iff `open_u` is.
      integer :: i, j, k
      real(wp) :: lap_u, lap_v, nu_eff, idt, sa, sb

      idt = 0.0_wp
      if (dt > 0.0_wp) idt = 1.0_wp/dt
      sa = 1.0_wp - 2.0_wp*ns
      sb = 2.0_wp*ns

      if (nu_h <= 0.0_wp) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            du_visc(i, j, k) = 0.0_wp
         end do
         do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
            dv_visc(i, j, k) = 0.0_wp
         end do
         return
      end if

      if (present(open_u)) then
         do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) local(lap_u, nu_eff)
            lap_u = iareaCu(i, j)*( &
                    (dy_dxT(i, j)*(u_face(i + 1, j, k) - u_face(i, j, k))*open_u(i + 1, j, k) - &
                     dy_dxT(i - 1, j)*(u_face(i, j, k) - u_face(i - 1, j, k))*open_u(i - 1, j, k)) + &
                    (dx_dyBu(i, j + 1)*(sa*wet_q(i, j + 1) + sb)*(u_face(i, j + 1, k) - u_face(i, j, k))*open_u(i, j + 1, k) - &
                     dx_dyBu(i, j)*(sa*wet_q(i, j) + sb)*(u_face(i, j, k) - u_face(i, j - 1, k))*open_u(i, j - 1, k)))
            nu_eff = nu_h
            if (bound_kh) then
               nu_eff = min(nu_eff, hvisc_kh_cfl_bound(idxCu(i, j), idyCu(i, j), bound_coef, idt))
            end if
            du_visc(i, j, k) = nu_eff*lap_u*open_u(i, j, k)
         end do
      else
         do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) local(lap_u, nu_eff)
            lap_u = iareaCu(i, j)*( &
                    (dy_dxT(i, j)*(u_face(i + 1, j, k) - u_face(i, j, k)) - &
                     dy_dxT(i - 1, j)*(u_face(i, j, k) - u_face(i - 1, j, k))) + &
                    (dx_dyBu(i, j + 1)*(sa*wet_q(i, j + 1) + sb)*(u_face(i, j + 1, k) - u_face(i, j, k)) - &
                     dx_dyBu(i, j)*(sa*wet_q(i, j) + sb)*(u_face(i, j, k) - u_face(i, j - 1, k))))
            nu_eff = nu_h
            if (bound_kh) then
               nu_eff = min(nu_eff, hvisc_kh_cfl_bound(idxCu(i, j), idyCu(i, j), bound_coef, idt))
            end if
            du_visc(i, j, k) = nu_eff*lap_u
         end do
      end if
      do concurrent(k=1:nz, j=1:ny)
         du_visc(1, j, k) = 0.0_wp
         du_visc(nx + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx + 1)
         du_visc(i, 1, k) = 0.0_wp
         du_visc(i, ny, k) = 0.0_wp
      end do

      if (present(open_v)) then
         do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) local(lap_v, nu_eff)
            lap_v = iareaCv(i, j)*( &
                    (dx_dyT(i, j)*(v_face(i, j + 1, k) - v_face(i, j, k))*open_v(i, j + 1, k) - &
                     dx_dyT(i, j - 1)*(v_face(i, j, k) - v_face(i, j - 1, k))*open_v(i, j - 1, k)) + &
                    (dy_dxBu(i + 1, j)*(sa*wet_q(i + 1, j) + sb)*(v_face(i + 1, j, k) - v_face(i, j, k))*open_v(i + 1, j, k) - &
                     dy_dxBu(i, j)*(sa*wet_q(i, j) + sb)*(v_face(i, j, k) - v_face(i - 1, j, k))*open_v(i - 1, j, k)))
            nu_eff = nu_h
            if (bound_kh) then
               nu_eff = min(nu_eff, hvisc_kh_cfl_bound(idxCv(i, j), idyCv(i, j), bound_coef, idt))
            end if
            dv_visc(i, j, k) = nu_eff*lap_v*open_v(i, j, k)
         end do
      else
         do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) local(lap_v, nu_eff)
            lap_v = iareaCv(i, j)*( &
                    (dx_dyT(i, j)*(v_face(i, j + 1, k) - v_face(i, j, k)) - &
                     dx_dyT(i, j - 1)*(v_face(i, j, k) - v_face(i, j - 1, k))) + &
                    (dy_dxBu(i + 1, j)*(sa*wet_q(i + 1, j) + sb)*(v_face(i + 1, j, k) - v_face(i, j, k)) - &
                     dy_dxBu(i, j)*(sa*wet_q(i, j) + sb)*(v_face(i, j, k) - v_face(i - 1, j, k))))
            nu_eff = nu_h
            if (bound_kh) then
               nu_eff = min(nu_eff, hvisc_kh_cfl_bound(idxCv(i, j), idyCv(i, j), bound_coef, idt))
            end if
            dv_visc(i, j, k) = nu_eff*lap_v
         end do
      end if
      do concurrent(k=1:nz, i=1:nx)
         dv_visc(i, 1, k) = 0.0_wp
         dv_visc(i, ny + 1, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, j=1:ny + 1)
         dv_visc(1, j, k) = 0.0_wp
         dv_visc(nx, j, k) = 0.0_wp
      end do
   end subroutine hvisc_compute_scalar_impl