ocean_lateral_mix_compute_smag Subroutine

   pure subroutine ocean_lateral_mix_compute_smag(grid, metrics, this, ms, &
                                                  res_fn_u, res_fn_v)
      !! Populate `ah_face_x`/`ah_face_y` (m^2/s) with the Smagorinsky
      !! Laplacian viscosity
      !!     A_h(face) = max(ah_bg, min(ah_max, (C_S · dx)^2 · |D|))
      !! where `|D| = sqrt(D_T^2 + D_S^2)` is the deformation-tensor
      !! magnitude — tension `D_T = ∂u/∂x − ∂v/∂y` (cell centred) and
      !! shear `D_S = ∂v/∂x + ∂u/∂y` (corner) — averaged onto the face.
      !! Wall faces get the background viscosity (wall-adjacent rows
      !! re-use the next interior row).  Smagorinsky (1963); C_S ≈ 0.15–0.2.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_lateral_mix_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in), optional :: res_fn_u(grid%nx_total + 1, grid%ny_total)
         !! VarMix resolution function at u-faces (nondim, [0,1]).  When
         !! present AND `this%resoln_scaled_visc`, scales `A_h` before clamp.
      real(wp), intent(in), optional :: res_fn_v(grid%nx_total, grid%ny_total + 1)
         !! VarMix resolution function at v-faces.

      integer :: i, j, k, nx, ny, nz
      real(wp) :: c_smag_local, smag_scale, ah_bg_local, ah_max_local
      real(wp) :: D_T_W, D_T_E, D_T_S, D_T_N, D_T_face
      real(wp) :: D_S_S, D_S_N, D_S_W, D_S_E, D_S_face
      real(wp) :: strain_mag, A_raw, ns
      logical :: do_resoln

      if (.not. this%is_init) return
      if (.not. allocated(ms%u_face_x_layer)) return
      if (.not. allocated(ms%v_face_y_layer)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      ! Hoist scalar reads off `this` to local variables (see Leith
      ! sibling) — guards against device-side descriptor walks.
      ah_bg_local = this%ah_bg
      ah_max_local = this%ah_max
      c_smag_local = this%c_smag
      ns = merge(1.0_wp, 0.0_wp, this%no_slip)   ! free-slip(0)/no-slip(1) selector (C1)
      ! Resolution-function scaling (Gap 1, Hallberg 2013) — see compute_leith.
      do_resoln = this%resoln_scaled_visc .and. present(res_fn_u) .and. &
                  present(res_fn_v)

      ! Strain components use per-stagger metric inverses (design §2):
      ! tension ∂u/∂x, ∂v/∂y on cell (idxT/idyT); shear ∂v/∂x, ∂u/∂y
      ! on the corner (idxBu/idyBu).  smag_scale = (C_S·sqrt(dxT·dyT))²
      ! per cell (= (C_S·dx)² on uniform square metrics, bit-reducing).

      ! ---- u-face viscosity (i-1/2, j) ----
      ! D_T at the face = mean of the two adjacent cell-centred values:
      !   D_T(i-1, j) and D_T(i, j).
      ! D_S at the face = mean of the two adjacent corner values along
      ! the face: D_S(i, j) at SW corner of (i, j) and D_S(i, j+1) at NW.
      !
      ! Curvilinear D_S at corner Bu(i,j) (= SW corner of T(i,j)):
      !   dvdx = dy_dxBu · (v(i,j)·idyCv(i,j)  - v(i-1,j)·idyCv(i-1,j))
      !   dudy = dx_dyBu · (u(i,j)·idxCu(i,j)  - u(i,j-1)·idxCu(i,j-1))
      !   D_S  = dvdx + dudy
      ! On uniform SQUARE metrics dy_dxBu=dx_dyBu=1 and idyCv=idxCu=1/dx,
      ! so this collapses to the old plain-difference form bit-for-bit.
      ! (design §2; mirrors MOM6 MOM_hor_visc shear-strain form)
      do concurrent(k=1:nz, j=2:ny - 1, i=2:nx) &
         local(D_T_W, D_T_E, D_S_S, D_S_N, D_T_face, D_S_face, &
               strain_mag, A_raw, smag_scale)
         smag_scale = (c_smag_local*sqrt(metrics%dxT(i, j)*metrics%dyT(i, j)))**2
         ! Cell-centred D_T at (i-1, j) and (i, j)
         D_T_W = (ms%u_face_x_layer(i, j, k) - ms%u_face_x_layer(i - 1, j, k))*metrics%idxT(i - 1, j) &
                 - (ms%v_face_y_layer(i - 1, j + 1, k) - ms%v_face_y_layer(i - 1, j, k))*metrics%idyT(i - 1, j)
         D_T_E = (ms%u_face_x_layer(i + 1, j, k) - ms%u_face_x_layer(i, j, k))*metrics%idxT(i, j) &
                 - (ms%v_face_y_layer(i, j + 1, k) - ms%v_face_y_layer(i, j, k))*metrics%idyT(i, j)
         D_T_face = 0.5_wp*(D_T_W + D_T_E)
         ! Corner D_S: ratio-bundle form (design §2).  Bu(i,j) = SW corner
         ! of T(i,j); v at Cv(i,j) / Cv(i-1,j), u at Cu(i,j) / Cu(i,j-1).
         ! Corner shear strain sh_xy masked by slip factor (C1): free-slip
         ! ×wet_q / no-slip ×(2-wet_q).  Bit-identical for all-wet.
         D_S_S = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
                 (metrics%dy_dxBu(i, j)* &
                  (ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j) - &
                   ms%v_face_y_layer(i - 1, j, k)*metrics%idyCv(i - 1, j)) + &
                  metrics%dx_dyBu(i, j)* &
                  (ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j) - &
                   ms%u_face_x_layer(i, j - 1, k)*metrics%idxCu(i, j - 1)))
         D_S_N = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j + 1) + 2.0_wp*ns)* &
                 (metrics%dy_dxBu(i, j + 1)* &
                  (ms%v_face_y_layer(i, j + 1, k)*metrics%idyCv(i, j + 1) - &
                   ms%v_face_y_layer(i - 1, j + 1, k)*metrics%idyCv(i - 1, j + 1)) + &
                  metrics%dx_dyBu(i, j + 1)* &
                  (ms%u_face_x_layer(i, j + 1, k)*metrics%idxCu(i, j + 1) - &
                   ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j)))
         D_S_face = 0.5_wp*(D_S_S + D_S_N)
         strain_mag = sqrt(D_T_face*D_T_face + D_S_face*D_S_face)
         if (do_resoln) then
            A_raw = smag_scale*strain_mag*res_fn_u(i, j)
         else
            A_raw = smag_scale*strain_mag
         end if
         this%ah_face_x(i, j, k) = min(ah_max_local, max(ah_bg_local, A_raw))
      end do
      ! j=1 and j=ny rows: re-use the j=2 / j=ny-1 values one row in.
      ! Avoids the j-1 / j+1 stencil walking into the wall.
      do concurrent(k=1:nz, i=2:nx)
         this%ah_face_x(i, 1, k) = this%ah_face_x(i, 2, k)
         this%ah_face_x(i, ny, k) = this%ah_face_x(i, ny - 1, k)
      end do
      ! Wall faces (i=1, i=nx+1): background.
      do concurrent(k=1:nz, j=1:ny)
         this%ah_face_x(1, j, k) = ah_bg_local
         this%ah_face_x(nx + 1, j, k) = ah_bg_local
      end do

      ! ---- v-face viscosity (i, j-1/2) ----
      ! Mirror of the u-face stencil.
      do concurrent(k=1:nz, j=2:ny, i=2:nx - 1) &
         local(D_T_S, D_T_N, D_S_W, D_S_E, D_T_face, D_S_face, &
               strain_mag, A_raw, smag_scale)
         smag_scale = (c_smag_local*sqrt(metrics%dxT(i, j)*metrics%dyT(i, j)))**2
         ! Cell-centred D_T at (i, j-1) and (i, j)
         D_T_S = (ms%u_face_x_layer(i + 1, j - 1, k) - ms%u_face_x_layer(i, j - 1, k))*metrics%idxT(i, j - 1) &
                 - (ms%v_face_y_layer(i, j, k) - ms%v_face_y_layer(i, j - 1, k))*metrics%idyT(i, j - 1)
         D_T_N = (ms%u_face_x_layer(i + 1, j, k) - ms%u_face_x_layer(i, j, k))*metrics%idxT(i, j) &
                 - (ms%v_face_y_layer(i, j + 1, k) - ms%v_face_y_layer(i, j, k))*metrics%idyT(i, j)
         D_T_face = 0.5_wp*(D_T_S + D_T_N)
         ! Corner D_S: ratio-bundle form (design §2).  Bu(i,j) = SW corner
         ! of T(i,j); v at Cv(i,j) / Cv(i-1,j), u at Cu(i,j) / Cu(i,j-1).
         ! Corner shear strain sh_xy masked by slip factor (C1).
         D_S_W = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
                 (metrics%dy_dxBu(i, j)* &
                  (ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j) - &
                   ms%v_face_y_layer(i - 1, j, k)*metrics%idyCv(i - 1, j)) + &
                  metrics%dx_dyBu(i, j)* &
                  (ms%u_face_x_layer(i, j, k)*metrics%idxCu(i, j) - &
                   ms%u_face_x_layer(i, j - 1, k)*metrics%idxCu(i, j - 1)))
         D_S_E = ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i + 1, j) + 2.0_wp*ns)* &
                 (metrics%dy_dxBu(i + 1, j)* &
                  (ms%v_face_y_layer(i + 1, j, k)*metrics%idyCv(i + 1, j) - &
                   ms%v_face_y_layer(i, j, k)*metrics%idyCv(i, j)) + &
                  metrics%dx_dyBu(i + 1, j)* &
                  (ms%u_face_x_layer(i + 1, j, k)*metrics%idxCu(i + 1, j) - &
                   ms%u_face_x_layer(i + 1, j - 1, k)*metrics%idxCu(i + 1, j - 1)))
         D_S_face = 0.5_wp*(D_S_W + D_S_E)
         strain_mag = sqrt(D_T_face*D_T_face + D_S_face*D_S_face)
         if (do_resoln) then
            A_raw = smag_scale*strain_mag*res_fn_v(i, j)
         else
            A_raw = smag_scale*strain_mag
         end if
         this%ah_face_y(i, j, k) = min(ah_max_local, max(ah_bg_local, A_raw))
      end do
      ! i=1, i=nx columns: re-use the i=2 / i=nx-1 values.
      do concurrent(k=1:nz, j=2:ny)
         this%ah_face_y(1, j, k) = this%ah_face_y(2, j, k)
         this%ah_face_y(nx, j, k) = this%ah_face_y(nx - 1, j, k)
      end do
      ! Wall faces (j=1, j=ny+1): background.
      do concurrent(k=1:nz, i=1:nx)
         this%ah_face_y(i, 1, k) = ah_bg_local
         this%ah_face_y(i, ny + 1, k) = ah_bg_local
      end do
   end subroutine ocean_lateral_mix_compute_smag