hvisc_compute_face_impl Subroutine

   pure subroutine hvisc_compute_face_impl(u_face, v_face, ah_face_x, ah_face_y, &
                                           du_visc, dv_visc, &
                                           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 × spatially-varying viscosity.
      !! Explicit-shape dummies so NVHPC stdpar emits a device kernel
      !! without per-launch descriptor walks.  See
      !! `metric_lap_u`/`metric_lap_v` for the curvilinear FV form.
      !! `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(in)    :: ah_face_x(nx + 1, ny, nz), ah_face_y(nx, ny + 1, nz)
      real(wp), intent(inout) :: du_visc(nx + 1, ny, nz), dv_visc(nx, ny + 1, nz)
      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`).  ABSENT (the default
         !! path) => the interior loops below are textually the ones this
         !! routine has always run => bit-identical.
         !!
         !! PRESENT => the closed faces become FREE-SLIP walls, which is
         !! what a z-level partial step is (Adcroft, Hill & Marshall
         !! 1997).  Two things happen, and both are needed:
         !!
         !! 1. each neighbour difference is multiplied by the NEIGHBOUR
         !!    face's open flag, so a closed neighbour -- whose velocity
         !!    `mask_layer_velocities` has zeroed -- contributes nothing.
         !!    Without it the zero reads as a Dirichlet-0 boundary, i.e.
         !!    NO-SLIP: the wall would exert `nu*u/dx^2` of drag on the
         !!    live face beside it every step, which is the opposite of
         !!    the free-slip a partial step is supposed to be;
         !! 2. the whole tendency is multiplied by the face's OWN open
         !!    flag, so a closed face gets exactly zero viscous tendency
         !!    (it is a wall; `mask_layer_velocities` would zero it
         !!    anyway, but leaving a tendency there would make
         !!    `ke_diss` -- which MEKE reads -- account for work done on
         !!    water that is not there).
      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

      ! u-face Laplacian interior + zero boundaries (curvilinear FV form)
      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 = ah_face_x(i, j, k)
            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 = ah_face_x(i, j, k)
            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

      ! v-face Laplacian interior + zero boundaries (curvilinear FV form)
      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 = ah_face_y(i, j, k)
            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 = ah_face_y(i, j, k)
            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_face_impl