hvisc_biharm_lap_closed Subroutine

   pure subroutine hvisc_biharm_lap_closed(u_face, v_face, lap_u, lap_v, wet_q, &
                                           dy_dxT, dx_dyBu, iareaCu, dx_dyT, dy_dxBu, iareaCv, &
                                           nx, ny, nz, open_u, open_v)
      !! Pass 1 of BOTH velocity biharmonics (scalar `nu_4` and the
      !! flow-aware `nu4_face_*`) under `&vcoord_nml zfixed_closed_faces`:
      !! the intermediate Laplacian `lap_u`/`lap_v` with a closed
      !! face-layer treated as a FREE-SLIP wall.
      !!
      !! Rule (gated by `tests/test_ocean_hvisc_biharm_zfixed.F90`): every
      !! face-to-face difference `(f_b - f_a)` carries the weight
      !! `open(a)*open(b)` — the neighbour's flag on the difference, the
      !! face's own on the result.  This is TEXTUALLY the gating of the
      !! harmonic closed-face kernels (`hvisc_compute_scalar_impl` /
      !! `hvisc_compute_face_impl` with `open_u` present), with the
      !! biharmonic's fixed free-slip corner factor `wet_q`, so this pass
      !! IS the harmonic closed-face Laplacian; Pass 2 applies the same
      !! gating to `lap`.  Consequences:
      !!
      !! * a closed neighbour contributes NOTHING — the stencil behaves as
      !!   if the neighbour held the face's own value (Neumann / mirror),
      !!   the per-layer analogue of the `wet_q` free-slip land corner.
      !!   Ungated, the zero `mask_layer_velocities` stores there is read
      !!   as a Dirichlet value: no-slip drag at every staircase step, and
      !!   a `lap` that reads across the wall;
      !! * a closed face holds `lap = 0` and its Pass-2 tendency is
      !!   multiplied by its own flag, so it receives exactly zero;
      !! * every difference weight is symmetric in (a, b), so the gated
      !!   Laplacian is `A^-1 L` with `L` symmetric negative semi-
      !!   definite, and for a constant `nu_4` the biharmonic
      !!   `-nu_4 A^-1 L A^-1 L` dissipates: `dE/dt = -nu_4 (Lu)^T A^-1 (Lu)
      !!   <= 0` (A = face areas).
      !!
      !! The gate is applied to the NORMAL differences as well as the
      !! tangential ones, as in the harmonic kernel (a closed face's
      !! normal velocity is a wall value, not an interior sample).  The
      !! array-edge mirror rows below are unchanged from the ungated pass.
      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) :: lap_u(nx + 1, ny, nz), lap_v(nx, ny + 1, nz)
      real(wp), intent(in)    :: wet_q(nx + 1, ny + 1)
         !! Bu-corner wet mask (free-slip corner factor, as in the ungated pass).
      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)
      real(wp), intent(in)    :: open_u(nx + 1, ny, nz)
         !! Per-layer 0/1 u-face open mask (`metrics%open_u`).
      real(wp), intent(in)    :: open_v(nx, ny + 1, nz)
         !! Per-layer 0/1 v-face open mask (`metrics%open_v`).
      integer :: i, j, k
      real(wp) :: l_u, l_v

      do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) local(l_u)
         l_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)*wet_q(i, j + 1)*(u_face(i, j + 1, k) - u_face(i, j, k))*open_u(i, j + 1, k) - &
                dx_dyBu(i, j)*wet_q(i, j)*(u_face(i, j, k) - u_face(i, j - 1, k))*open_u(i, j - 1, k)))
         lap_u(i, j, k) = l_u*open_u(i, j, k)
      end do
      do concurrent(k=1:nz, j=1:ny)
         lap_u(1, j, k) = lap_u(2, j, k)
         lap_u(nx + 1, j, k) = lap_u(nx, j, k)
      end do
      do concurrent(k=1:nz, i=1:nx + 1)
         lap_u(i, 1, k) = lap_u(i, 2, k)
         lap_u(i, ny, k) = lap_u(i, ny - 1, k)
      end do

      do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) local(l_v)
         l_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)*wet_q(i + 1, j)*(v_face(i + 1, j, k) - v_face(i, j, k))*open_v(i + 1, j, k) - &
                dy_dxBu(i, j)*wet_q(i, j)*(v_face(i, j, k) - v_face(i - 1, j, k))*open_v(i - 1, j, k)))
         lap_v(i, j, k) = l_v*open_v(i, j, k)
      end do
      do concurrent(k=1:nz, i=1:nx)
         lap_v(i, 1, k) = lap_v(i, 2, k)
         lap_v(i, ny + 1, k) = lap_v(i, ny, k)
      end do
      do concurrent(k=1:nz, j=1:ny + 1)
         lap_v(1, j, k) = lap_v(2, j, k)
         lap_v(nx, j, k) = lap_v(nx - 1, j, k)
      end do
   end subroutine hvisc_biharm_lap_closed