hvisc_compute_biharmonic_impl Subroutine

   pure subroutine hvisc_compute_biharmonic_impl(u_face, v_face, lap_u, lap_v, &
                                                 du_visc, dv_visc, &
                                                 nu_4, bound_coef, dt, &
                                                 idxCu, idyCu, idxCv, idyCv, wet_q, &
                                                 dy_dxT, dx_dyBu, iareaCu, &
                                                 dx_dyT, dy_dxBu, iareaCv, &
                                                 nx, ny, nz, open_u, open_v)
      !! Constant-coefficient biharmonic friction: applies
      !! `-ν₄ · ∇²(∇²u)` to the face velocities via two chained 5-point
      !! Laplacians.  Adds into the existing `du_visc` / `dv_visc`
      !! buffers (which Laplacian friction has already filled), so the
      !! caller can run with both `nu_h` and `nu_4` non-zero.
      !!
      !! Wall convention: the first Laplacian writes zero at every
      !! C-grid wall face (i=1, i=nx+1 on u; j=1, j=ny+1 on v) and at
      !! the interior tangential-boundary rows.  The second Laplacian
      !! reads `lap_u / lap_v` and sees those zeros — equivalent to a
      !! "no-stress on Δu" wall condition, which is the closure
      !! MOM6's `BIHARMONIC` block uses for closed boundaries.
      !!
      !! Sign convention: the operator is `du/dt = -ν₄·∇⁴u` so the
      !! sinusoid `u = e^{ikx}` has growth rate `-ν₄·k⁴` — damping for
      !! positive ν₄, scaling as `k⁴` rather than `k²` (Laplacian).
      !!
      !! Stability: forward-Euler biharmonic CFL is
      !! `ν₄ · dt · ((π/dx)² + (π/dy)²)² <= 2` (established project
      !! constant).  The constant scalar `nu_4` is clamped per face to
      !! this CFL ceiling (`hvisc_nu4_cfl_bound`, scaled by
      !! `bound_coef`) in Pass 2, so an over-large namelist `nu_4`
      !! cannot violate the local bound; below the ceiling the multiply
      !! is by `nu_4` bit-for-bit.
      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(inout) :: du_visc(nx + 1, ny, nz), dv_visc(nx, ny + 1, nz)
      real(wp), intent(in)    :: nu_4
      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`).  BOTH chained Laplacians mask
         !! their corner (shear) fluxes by it -- free-slip, MOM6
         !! `sh_xy = mask2dBu·(...)` for `Del2u` and `str_xy·mask2dBu` for the
         !! biharmonic stress.  Unmasked, the zero stored at a land face read
         !! as a Dirichlet-0 wall and the k^4 operator rang against it at every
         !! staircase step.  Always free-slip: MOM6 refuses NOSLIP with
         !! BIHARMONIC.  All-wet corner => factor 1 => bit-identical.
      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), 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 loops
         !! below are textually the ones this routine has always run =>
         !! bit-identical.  PRESENT => see `hvisc_biharm_lap_closed`: in
         !! BOTH chained Laplacians every face-to-face difference carries
         !! `open(a)*open(b)`, so a closed face-layer is a free-slip
         !! (Neumann/mirror) boundary of the stencil exactly as a land
         !! corner is under `wet_q`, and receives zero tendency.
      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) :: l_u, l_v, idt, nu4_u, nu4_v, sa, sb
      idt = 1.0_wp/dt
      sa = 1.0_wp
      sb = 0.0_wp

      if (present(open_u)) then
         ! z-level closed faces: Pass 1 is the gated Laplacian (closed
         ! face-layers hold lap = 0 and are never read as a value);
         ! Pass 2 gates every difference by the NEIGHBOUR's open flag
         ! and the tendency by the face's OWN.  Same per-face CFL clamp as below.
         call 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)
         do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) local(l_u, nu4_u)
            l_u = iareaCu(i, j)*( &
                  (dy_dxT(i, j)*(lap_u(i + 1, j, k) - lap_u(i, j, k))*open_u(i + 1, j, k) - &
                   dy_dxT(i - 1, j)*(lap_u(i, j, k) - lap_u(i - 1, j, k))*open_u(i - 1, j, k)) + &
                  (dx_dyBu(i, j + 1)*wet_q(i, j + 1)*(lap_u(i, j + 1, k) - lap_u(i, j, k))*open_u(i, j + 1, k) - &
                   dx_dyBu(i, j)*wet_q(i, j)*(lap_u(i, j, k) - lap_u(i, j - 1, k))*open_u(i, j - 1, k)))
            nu4_u = min(nu_4, hvisc_nu4_cfl_bound(idxCu(i, j), idyCu(i, j), bound_coef, idt))
            du_visc(i, j, k) = du_visc(i, j, k) - nu4_u*l_u*open_u(i, j, k)
         end do
         do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) local(l_v, nu4_v)
            l_v = iareaCv(i, j)*( &
                  (dx_dyT(i, j)*(lap_v(i, j + 1, k) - lap_v(i, j, k))*open_v(i, j + 1, k) - &
                   dx_dyT(i, j - 1)*(lap_v(i, j, k) - lap_v(i, j - 1, k))*open_v(i, j - 1, k)) + &
                  (dy_dxBu(i + 1, j)*wet_q(i + 1, j)*(lap_v(i + 1, j, k) - lap_v(i, j, k))*open_v(i + 1, j, k) - &
                   dy_dxBu(i, j)*wet_q(i, j)*(lap_v(i, j, k) - lap_v(i - 1, j, k))*open_v(i - 1, j, k)))
            nu4_v = min(nu_4, hvisc_nu4_cfl_bound(idxCv(i, j), idyCv(i, j), bound_coef, idt))
            dv_visc(i, j, k) = dv_visc(i, j, k) - nu4_v*l_v*open_v(i, j, k)
         end do
         return
      end if

      ! ---- Pass 1: lap_u, lap_v at interior u-faces and v-faces ----
      ! Wall BC on the intermediate Laplacian: MIRROR (Neumann),
      ! not zero.  Hard-zero on lap_u/lap_v creates a step
      ! discontinuity that Pass 2 reads as a high-k mode at the
      ! wall-adjacent interior face → injects spurious work into the
      ! boundary band.  Mirror BC (copy the adjacent interior value)
      ! preserves the "smooth interior field ⇒ zero biharmonic"
      ! property — a smooth lap_u extended by mirror gives ∇² = 0 at
      ! the boundary too.  This is the equivalent of MOM6's "stress
      ! free at the boundary" closure on biharmonic friction.
      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)) - &
                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))))
         lap_u(i, j, k) = l_u
      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)) - &
                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))))
         lap_v(i, j, k) = l_v
      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

      ! ---- Pass 2: -ν₄ · ∇²(lap_*) added into du_visc / dv_visc ----
      ! ν₄ is clamped per face to the local explicit-biharmonic CFL
      ! ceiling (bound_coef · 2/(dt·k⁴)).  Below the ceiling
      ! `min(nu_4, bound) == nu_4` ⇒ bit-identical to the unclamped add.
      do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) local(l_u, nu4_u)
         l_u = iareaCu(i, j)*( &
               (dy_dxT(i, j)*(lap_u(i + 1, j, k) - lap_u(i, j, k)) - &
                dy_dxT(i - 1, j)*(lap_u(i, j, k) - lap_u(i - 1, j, k))) + &
               (dx_dyBu(i, j + 1)*(sa*wet_q(i, j + 1) + sb)*(lap_u(i, j + 1, k) - lap_u(i, j, k)) - &
                dx_dyBu(i, j)*(sa*wet_q(i, j) + sb)*(lap_u(i, j, k) - lap_u(i, j - 1, k))))
         nu4_u = min(nu_4, hvisc_nu4_cfl_bound(idxCu(i, j), idyCu(i, j), bound_coef, idt))
         du_visc(i, j, k) = du_visc(i, j, k) - nu4_u*l_u
      end do

      do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) local(l_v, nu4_v)
         l_v = iareaCv(i, j)*( &
               (dx_dyT(i, j)*(lap_v(i, j + 1, k) - lap_v(i, j, k)) - &
                dx_dyT(i, j - 1)*(lap_v(i, j, k) - lap_v(i, j - 1, k))) + &
               (dy_dxBu(i + 1, j)*(sa*wet_q(i + 1, j) + sb)*(lap_v(i + 1, j, k) - lap_v(i, j, k)) - &
                dy_dxBu(i, j)*(sa*wet_q(i, j) + sb)*(lap_v(i, j, k) - lap_v(i - 1, j, k))))
         nu4_v = min(nu_4, hvisc_nu4_cfl_bound(idxCv(i, j), idyCv(i, j), bound_coef, idt))
         dv_visc(i, j, k) = dv_visc(i, j, k) - nu4_v*l_v
      end do
   end subroutine hvisc_compute_biharmonic_impl