hvisc_compute_biharmonic_face_impl Subroutine

   pure subroutine hvisc_compute_biharmonic_face_impl(u_face, v_face, lap_u, lap_v, &
                                                      du_visc, dv_visc, &
                                                      nu4_face_x, nu4_face_y, &
                                                      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)
      !! Flow-aware biharmonic friction (MOM6 SMAGORINSKY_AH analogue).
      !! Identical to `hvisc_compute_biharmonic_impl` except Pass 2
      !! multiplies the second Laplacian by the per-face viscosity
      !! `nu4_face_x/y` instead of the scalar `nu_4`.  The face fields
      !! are filled upstream by `ocean_lateral_mix_compute_smag_ah`,
      !! which sets them to `C_b · L⁴ · |D|` clamped to
      !! `[nu4_bg, nu4_max]`.  Pass 2 additionally clamps each face
      !! coefficient to the per-face explicit-biharmonic CFL ceiling
      !! (`hvisc_nu4_cfl_bound`, scaled by `bound_coef`) on top of the
      !! static `nu4_max` floor/ceiling applied upstream — so a strain
      !! spike on a fine cell can never violate the local CFL bound.
      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)    :: nu4_face_x(nx + 1, ny, nz), nu4_face_y(nx, ny + 1, nz)
      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(nu4_face_x(i, j, k), &
                        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(nu4_face_y(i, j, k), &
                        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/v faces (metric FV) ----
      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: -nu4_face · ∇²(lap_*) added into du_visc / dv_visc ----
      ! Each per-face nu4 is clamped to the local CFL ceiling; below it
      ! `min(...) == nu4_face` ⇒ 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(nu4_face_x(i, j, k), &
                     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(nu4_face_y(i, j, k), &
                     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_face_impl