hvisc_compute_stress Subroutine

   pure subroutine hvisc_compute_stress(u_face, v_face, h_layer, &
                                        str_xx, str_xy, ah_t, ah_q, &
                                        du_visc, dv_visc, ns, &
                                        kh_aniso, n1n2, n1n1_m_n2n2, &
                                        idxCu, idyCu, idxCv, idyCv, &
                                        dx2h, dy2h, dx2q, dy2q, &
                                        dy_dxT, dx_dyT, dy_dxBu, dx_dyBu, &
                                        iareaCu, iareaCv, &
                                        wet_u, wet_v, wet_q, nx, ny, nz)
      !! MOM6 thickness-weighted stress-divergence operator.  Three
      !! phases: (1) tension `str_xx` at T-cells, (2) shear `str_xy` at
      !! Bu corners, (3) the divergence `(1/(h_u+h_neglect))·∂str`.
      !! `wet_u/wet_v/wet_q` mask the stress so no momentum is diffused
      !! across a coastline; `h_neglect = H_VANISHED` floors the
      !! velocity-point thickness.  See module header for the form +
      !! the all-wet uniform-h Laplacian reduction.
      !!
      !! Anisotropic cross terms (Smith & McWilliams 2003, §2): when
      !! `kh_aniso > 0` the tension stress at T-cells gains
      !! `+kh_aniso·n1n2·(n1²−n2²)·sh_xy_at_T·h_T` and the shear stress
      !! at corners gains `+kh_aniso·n1n2·(n1²−n2²)·sh_xx_at_q·h_q`,
      !! where the off-component strain is 4-point averaged onto the
      !! stress location.  (Sign is `+` here because Roundabout's stress
      !! convention `str = +A·strain·h` reduces the divergence to
      !! `+A·∇²u`, opposite to MOM6's `str = −Kh·sh`.)  For the default
      !! grid-i direction `n1n2 = 0` ⇒ the cross terms vanish and the
      !! assembly is bit-identical to the isotropic path.
      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)    :: h_layer(nx, ny, nz)
      real(wp), intent(inout) :: str_xx(nx, ny, nz), str_xy(nx + 1, ny + 1, nz)
      real(wp), intent(in)    :: ah_t(nx, ny, nz), ah_q(nx + 1, ny + 1, nz)
      real(wp), intent(inout) :: du_visc(nx + 1, ny, nz), dv_visc(nx, ny + 1, nz)
      real(wp), intent(in)    :: ns
      real(wp), intent(in)    :: kh_aniso, n1n2, n1n1_m_n2n2
      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)    :: dx2h(nx, ny), dy2h(nx, ny)
      real(wp), intent(in)    :: dx2q(nx + 1, ny + 1), dy2q(nx + 1, ny + 1)
      real(wp), intent(in)    :: dy_dxT(nx, ny), dx_dyT(nx, ny)
      real(wp), intent(in)    :: dy_dxBu(nx + 1, ny + 1), dx_dyBu(nx + 1, ny + 1)
      real(wp), intent(in)    :: iareaCu(nx + 1, ny), iareaCv(nx, ny + 1)
      real(wp), intent(in)    :: wet_u(nx + 1, ny), wet_v(nx, ny + 1)
      real(wp), intent(in)    :: wet_q(nx + 1, ny + 1)
      integer :: i, j, k
      real(wp) :: dudx, dvdy, dvdx, dudy, h_q, slip
      real(wp) :: aniso_cross, shxy_at_t, shxx_at_q
      real(wp), parameter :: H_NEGLECT = H_VANISHED

      aniso_cross = kh_aniso*n1n2*n1n1_m_n2n2

      ! ---- Phase 1: tension stress at T-cells (i,j) ----
      ! sh_xx = du/dx − dv/dy.  Masking the velocity gradients by the
      ! face wet masks keeps a land neighbour from contributing strain
      ! (MOM6 masks via reduction_xx + pre-masked metrics; the wet_u/
      ! wet_v product is the equivalent here).  All-wet ⇒ ×1.
      do concurrent(k=1:nz, j=1:ny, i=1:nx) local(dudx, dvdy)
         dudx = dy_dxT(i, j)*(idyCu(i + 1, j)*wet_u(i + 1, j)*u_face(i + 1, j, k) - &
                              idyCu(i, j)*wet_u(i, j)*u_face(i, j, k))
         dvdy = dx_dyT(i, j)*(idxCv(i, j + 1)*wet_v(i, j + 1)*v_face(i, j + 1, k) - &
                              idxCv(i, j)*wet_v(i, j)*v_face(i, j, k))
         str_xx(i, j, k) = ah_t(i, j, k)*(dudx - dvdy)*h_layer(i, j, k)
      end do

      ! ---- Phase 2: shear stress at Bu corners (i,j) = SW of cell (i,j) ----
      ! sh_xy = dv/dx + du/dy.  Corner thickness h_q = 4-pt mean of the
      ! surrounding cells (no halo issues — only interior corners are
      ! read by the divergence).  Slip factor: free-slip ×wet_q,
      ! no-slip ×(2−wet_q); all-wet ⇒ ×1.
      do concurrent(k=1:nz, j=2:ny, i=2:nx) local(dvdx, dudy, h_q, slip)
         dvdx = dy_dxBu(i, j)*(idyCv(i, j)*v_face(i, j, k) - idyCv(i - 1, j)*v_face(i - 1, j, k))
         dudy = dx_dyBu(i, j)*(idxCu(i, j)*u_face(i, j, k) - idxCu(i, j - 1)*u_face(i, j - 1, k))
         h_q = 0.25_wp*((h_layer(i - 1, j - 1, k) + h_layer(i, j, k)) + &
                        (h_layer(i - 1, j, k) + h_layer(i, j - 1, k)))
         slip = (1.0_wp - 2.0_wp*ns)*wet_q(i, j) + 2.0_wp*ns
         str_xy(i, j, k) = ah_q(i, j, k)*(dvdx + dudy)*h_q*slip
      end do
      do concurrent(k=1:nz, j=1:ny + 1)
         str_xy(1, j, k) = 0.0_wp
         str_xy(nx + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx + 1)
         str_xy(i, 1, k) = 0.0_wp
         str_xy(i, ny + 1, k) = 0.0_wp
      end do

      ! ---- Phase 2.5: anisotropic cross terms (Smith & McWilliams 2003) ----
      ! Folded ONLY when kh_aniso·n1n2·(n1²−n2²) /= 0 — i.e. an
      ! off-grid-axis anisotropy direction.  The default grid-i
      ! direction (n1n2 = 0) skips both passes ⇒ bit-identical.  The
      ! cross terms couple the tension stress to the shear strain and
      ! vice versa, using the RAW velocity-gradient strains (not the
      ! already-formed stresses), recomputed inline:
      !   sh_xx(T)  = du/dx − dv/dy   (Phase-1 strain)
      !   sh_xy(Bu) = dv/dx + du/dy   (Phase-2 strain)
      ! str_xx(T) += aniso_cross · ⟨sh_xy⟩_corners→T · h_T
      ! str_xy(Bu)+= aniso_cross · ⟨sh_xx⟩_cells→Bu  · h_q · slip
      if (aniso_cross /= 0.0_wp) then
         ! str_xx gains the shear-strain contribution: average sh_xy
         ! from the 4 corners around T-cell (i,j): (i,j),(i+1,j),
         ! (i,j+1),(i+1,j+1).  Interior corners are valid for i=1:nx-1,
         ! j=1:ny-1; the boundary T-rows keep the isotropic value (the
         ! divergence weights them through the metric stencil).
         do concurrent(k=1:nz, j=2:ny - 1, i=2:nx - 1) local(shxy_at_t)
            shxy_at_t = 0.25_wp*( &
                        (raw_sh_xy(u_face, v_face, idxCu, idyCv, dy_dxBu, dx_dyBu, i, j, k, nx, ny, nz) + &
                         raw_sh_xy(u_face, v_face, idxCu, idyCv, dy_dxBu, dx_dyBu, i + 1, j + 1, k, nx, ny, nz)) + &
                        (raw_sh_xy(u_face, v_face, idxCu, idyCv, dy_dxBu, dx_dyBu, i + 1, j, k, nx, ny, nz) + &
                         raw_sh_xy(u_face, v_face, idxCu, idyCv, dy_dxBu, dx_dyBu, i, j + 1, k, nx, ny, nz)))
            str_xx(i, j, k) = str_xx(i, j, k) + aniso_cross*shxy_at_t*h_layer(i, j, k)
         end do
         ! str_xy gains the tension-strain contribution: average sh_xx
         ! from the 4 T-cells around corner (i,j): (i-1,j-1),(i,j-1),
         ! (i-1,j),(i,j).  Slip + h_q reuse Phase-2 forms.
         do concurrent(k=1:nz, j=2:ny, i=2:nx) local(shxx_at_q, h_q, slip)
            shxx_at_q = 0.25_wp*( &
                        (raw_sh_xx(u_face, v_face, idxCu, idyCu, idxCv, dy_dxT, dx_dyT, i - 1, j - 1, k, nx, ny, nz) + &
                         raw_sh_xx(u_face, v_face, idxCu, idyCu, idxCv, dy_dxT, dx_dyT, i, j, k, nx, ny, nz)) + &
                        (raw_sh_xx(u_face, v_face, idxCu, idyCu, idxCv, dy_dxT, dx_dyT, i, j - 1, k, nx, ny, nz) + &
                         raw_sh_xx(u_face, v_face, idxCu, idyCu, idxCv, dy_dxT, dx_dyT, i - 1, j, k, nx, ny, nz)))
            h_q = 0.25_wp*((h_layer(i - 1, j - 1, k) + h_layer(i, j, k)) + &
                           (h_layer(i - 1, j, k) + h_layer(i, j - 1, k)))
            slip = (1.0_wp - 2.0_wp*ns)*wet_q(i, j) + 2.0_wp*ns
            str_xy(i, j, k) = str_xy(i, j, k) + aniso_cross*shxx_at_q*h_q*slip
         end do
      end if

      ! ---- Phase 3a: u-face divergence ----
      ! diffu(i,j) = iareaCu·( idyCu·(dy2h·str_xx(i,j) − dy2h·str_xx(i-1,j))
      !                      + idxCu·(dx2q·str_xy(i,j+1) − dx2q·str_xy(i,j)) )
      !            / (h_u + h_neglect).   (signs verified to reduce to
      ! +A·∇²u on uniform-h all-wet; see module header.)
      do concurrent(k=1:nz, j=2:ny - 1, i=2:nx)
         du_visc(i, j, k) = wet_u(i, j)*iareaCu(i, j)* &
                            (idyCu(i, j)*(dy2h(i, j)*str_xx(i, j, k) - dy2h(i - 1, j)*str_xx(i - 1, j, k)) + &
                             idxCu(i, j)*(dx2q(i, j + 1)*str_xy(i, j + 1, k) - dx2q(i, j)*str_xy(i, j, k)))/ &
                            (0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k)) + H_NEGLECT)
      end do
      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

      ! ---- Phase 3b: v-face divergence ----
      ! diffv(i,j) = iareaCv·( idxCv·(−dx2h·str_xx(i,j) + dx2h·str_xx(i,j-1))
      !                      + idyCv·(dy2q·str_xy(i+1,j) − dy2q·str_xy(i,j)) )
      !            / (h_v + h_neglect).  Tension enters with the
      ! opposite sign on the v-face (MOM6 diffv).
      do concurrent(k=1:nz, j=2:ny, i=2:nx - 1)
         dv_visc(i, j, k) = wet_v(i, j)*iareaCv(i, j)* &
                            (idyCv(i, j)*(dy2q(i + 1, j)*str_xy(i + 1, j, k) - dy2q(i, j)*str_xy(i, j, k)) - &
                             idxCv(i, j)*(dx2h(i, j)*str_xx(i, j, k) - dx2h(i, j - 1)*str_xx(i, j - 1, k)))/ &
                            (0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k)) + H_NEGLECT)
      end do
      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_stress