ocean_horizontal_viscosity_compute_tendencies_on Subroutine

   subroutine ocean_horizontal_viscosity_compute_tendencies_on(grid, metrics, this, ms, u, v, h, lateral_mix, dt)
      !! Fill `du_visc` and `dv_visc` with `nu_h * Laplacian` of the
      !! face velocities, per layer.  Closed-wall faces (i=1, i=nx+1
      !! for u; j=1, j=ny+1 for v) get zero tendency.  Interior y-
      !! boundary rows on u (j=1, j=ny) and interior x-boundary
      !! columns on v (i=1, i=nx) also get zero — equivalent to a
      !! free-slip wall condition on the tangential velocity.
      !!
      !! Loop order: (k, j, i) — `j` outermost-but-one for NVHPC GPU
      !! coalescing on the innermost-array-dimension `i`.
      !!
      !! Curvilinear (design §2): the face-velocity Laplacian is the
      !! finite-volume divergence of the velocity gradient over the
      !! face control volume — x-flux differenced across T-points
      !! (`dy_dxT` ratio), y-flux across Bu corners (`dx_dyBu`),
      !! normalised by `iareaCu`; the v-face mirror uses `dx_dyT` /
      !! `dy_dxBu` / `iareaCv`.  On uniform SQUARE metrics every ratio
      !! is 1 and `iareaCu = 1/(dx·dy)`, so the form collapses to the
      !! decoupled `Δ²u·(1/dx²) + Δ²u·(1/dy²)` to round-off (the only
      !! departure is FP reassociation of the 3-point grouping).
      !!
      !! Closure dispatch: when `lateral_mix` is present and its
      !! `closure` field is non-default (i.e. `/= LMIX_NONE`), the
      !! kernel reads per-face viscosity from
      !! `lateral_mix%ah_face_x` / `ah_face_y` instead of the scalar
      !! `nu_h`.  Caller must invoke
      !! `ocean_lateral_mix_compute_leith` (or equivalent) first to
      !! populate those fields.  Omitting the argument or leaving
      !! `closure = LMIX_NONE` preserves the scalar-`nu_h` behaviour
      !! bit-identically — existing call sites stay valid.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_horizontal_viscosity_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      ! assumed-shape-ok: forwarded straight to explicit-shape impl dummies,
      ! never indexed here (outer-shim source-selection layer, SPEC S3).
      real(wp), intent(in) :: u(:, :, :), v(:, :, :), h(:, :, :)
         !! Velocity/thickness source arrays — the prognostic components on
         !! the historical path, the `u_av` time-mean family under
         !! `split_scheme = "pred_corr"` (MOM6 evaluates horizontal_viscosity
         !! on `u_av`/`h_av`).
      type(ocean_lateral_mix_t), intent(in), optional :: lateral_mix
      real(wp), intent(in), optional :: dt
         !! Outer (or RK2-stage) time step.  Required when
         !! `this%stress_tensor` is true (drives the per-cell CFL
         !! viscosity limiter) and whenever a biharmonic add-on is
         !! active (drives its own per-face CFL clamp) — the velocity-
         !! Laplacian dispatch itself does not consume it, but the
         !! biharmonic block reached from every dispatch arm does.
         !! Omitting it defaults the biharmonic clamp to `dt_local = 1`
         !! (a huge, effectively-inactive bound); the caller is
         !! responsible for CFL in that case.

      integer :: nx, ny, nz
      real(wp) :: nu_h, nu_4, dt_local, slip_ns
      logical :: use_face_visc

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nu_h = this%nu_h
      nu_4 = this%nu_4
      ! Coastal slip selector for the velocity-Laplacian corner fluxes
      ! (0 = free-slip, 1 = no-slip); hoisted so no kernel reads `this`.
      slip_ns = merge(1.0_wp, 0.0_wp, this%no_slip)

      ! dt for the per-cell CFL clamps (harmonic BOUND_KH on the
      ! stress-tensor path, biharmonic on the add-on below).  Defaults
      ! to 1 so the bound is huge when no dt was supplied — caller
      ! respects CFL.  Single hoisted site: both the stress branch and
      ! the biharmonic add-on below now read it.
      dt_local = 1.0_wp
      if (present(dt)) dt_local = dt

      ! Flow-aware HARMONIC face viscosity (`ah_face_*`) is read for the
      ! Laplacian closures (Leith / Smagorinsky) and for the live
      ! velocity-scale floor (`kh_vel_scale_live > 0`), which also
      ! populates `ah_face_*`.  The biharmonic closures
      ! (`LMIX_LEITH_BIHARM`, `smag_ah_active`) fill `nu4_face_*` instead
      ! and engage via the biharmonic add-on below; they must NOT pull the
      ! (unfilled) `ah_face_*` into the Laplacian — the scalar `nu_h`
      ! Laplacian stays in force underneath them.
      use_face_visc = .false.
      if (present(lateral_mix)) then
         if (lateral_mix%is_init .and. &
             (lateral_mix%closure == LMIX_LEITH .or. &
              lateral_mix%closure == LMIX_SMAGORINSKY .or. &
              lateral_mix%kh_vel_scale_live > 0.0_wp)) then
            use_face_visc = .true.
         end if
      end if

      ! ---- MOM6-faithful thickness-weighted stress-divergence path ----
      ! Replaces the velocity Laplacian with a single momentum-
      ! conserving stress assembly; per-cell CFL limiter + wet_*
      ! coast-masking baked in.  Composes with the biharmonic add-on
      ! below exactly as the Laplacian arms do (MOM6 `BIHARMONIC` "may
      ! be used with `LAPLACIAN`") — the harmonic part is momentum-
      ! conserving/coast-masked, the velocity-form biharmonic part is
      ! not (a documented fidelity divergence; a stress-level
      ! biharmonic is the follow-on).
      if (this%stress_tensor) then
         ! Step A: average the per-face harmonic A onto T-cells (ah_t)
         ! and Bu corners (ah_q), reading the flow-aware lateral-mix
         ! fields when active, else the scalar nu_h.  Then CFL-clamp
         ! each per its discrete stencil bound.
         if (use_face_visc) then
            call hvisc_avg_A_face( &
               lateral_mix%ah_face_x, lateral_mix%ah_face_y, &
               this%ah_t%data, this%ah_q%data, nx, ny, nz)
         else
            call hvisc_fill_A_scalar(this%ah_t%data, this%ah_q%data, nu_h, nx, ny, nz)
         end if

         ! Step A2 (anisotropic, Smith & McWilliams 2003): ADD the
         ! direction-tensor coefficients to the isotropic A before the
         ! CFL clamp (MOM6 adds anisotropy into Kh prior to BOUND_KH).
         ! Tension (T-cell): Kh += kh_aniso·(1−n1n2²);
         ! shear (Bu corner): Kh += kh_aniso·n1n2².  No-op when
         ! kh_aniso ≤ 0 ⇒ bit-identical isotropic path.
         if (this%kh_aniso > 0.0_wp) then
            call hvisc_add_aniso_coef( &
               this%ah_t%data, this%ah_q%data, &
               this%kh_aniso, this%aniso_n1n2, nx, ny, nz)
         end if

         call hvisc_clamp_A( &
            this%ah_t%data, this%ah_q%data, this%bound_coef, dt_local, &
            metrics%idxT, metrics%idyT, metrics%idxCu, metrics%idyCu, &
            metrics%idxCv, metrics%idyCv, metrics%iareaCu, metrics%iareaCv, &
            nx, ny, nz)

         ! Step B: stress assembly + thickness-weighted divergence.
         ! The anisotropic CROSS terms (kh_aniso·n1n2·(n1²−n2²)·strain)
         ! are folded into str_xx / str_xy inside the assembly when
         ! kh_aniso > 0; they vanish for the default grid-i direction
         ! (n1n2 = 0).
         call hvisc_compute_stress( &
            u, v, h, &
            this%str_xx%data, this%str_xy%data, &
            this%ah_t%data, this%ah_q%data, &
            this%du_visc%data, this%dv_visc%data, &
            merge(1.0_wp, 0.0_wp, this%no_slip), &
            this%kh_aniso, this%aniso_n1n2, this%aniso_n1n1_m_n2n2, &
            metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, &
            metrics%dx2h, metrics%dy2h, metrics%dx2q, metrics%dy2q, &
            metrics%dy_dxT, metrics%dx_dyT, metrics%dy_dxBu, metrics%dx_dyBu, &
            metrics%iareaCu, metrics%iareaCv, &
            metrics%wet_u, metrics%wet_v, metrics%wet_q, &
            nx, ny, nz)
         ! Outer shim: dereference the multilayer + lateral-mix + this
         ! components on host, then dispatch to a flat-impl with explicit-
         ! shape array dummies.  Without this layering NVHPC's `do
         ! concurrent` body sees `this%du_visc%data`, `lateral_mix%ah_face_x`,
         ! `u` as derived-type deep derefs and emits per-
         ! iteration descriptor-walk memcpys.
      else if (use_face_visc) then
         ! z-level closed faces: the mask actuals are ABSENT on the
         ! default path (the `(1,1,1)` placeholder must never reach an
         ! explicit-shape dummy), so the call is written twice rather
         ! than the argument once.  Cold dispatcher code.
         if (metrics%use_closed_faces) then
            call hvisc_compute_face_impl( &
               u, v, &
               lateral_mix%ah_face_x, lateral_mix%ah_face_y, &
               this%du_visc%data, this%dv_visc%data, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               this%bound_kh, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, &
               metrics%wet_q, slip_ns, &
               nx, ny, nz, metrics%open_u, metrics%open_v)
         else
            call hvisc_compute_face_impl( &
               u, v, &
               lateral_mix%ah_face_x, lateral_mix%ah_face_y, &
               this%du_visc%data, this%dv_visc%data, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               this%bound_kh, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, &
               metrics%wet_q, slip_ns, &
               nx, ny, nz)
         end if
      else
         if (metrics%use_closed_faces) then
            call hvisc_compute_scalar_impl( &
               u, v, &
               this%du_visc%data, this%dv_visc%data, &
               nu_h, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               this%bound_kh, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, &
               metrics%wet_q, slip_ns, &
               nx, ny, nz, metrics%open_u, metrics%open_v)
         else
            call hvisc_compute_scalar_impl( &
               u, v, &
               this%du_visc%data, this%dv_visc%data, &
               nu_h, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               this%bound_kh, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, &
               metrics%wet_q, slip_ns, &
               nx, ny, nz)
         end if
      end if

      ! Biharmonic add-on.  Independent of the harmonic dispatch above
      ! (stress-tensor, face, or scalar) — sums into the same tendency
      ! buffer whichever of the three wrote it.  When `nu_4 = 0` the
      ! kernel is a no-op (early return).  Scale-selective damping for
      ! stratified closed-basin runs: damps grid-scale modes
      ! O(ν_4·k⁴) without spilling into basin-scale modes the way a
      ! large Laplacian ν_h does.  Required for production Tasman-class
      ! runs where the dyn-core has a parametric instability that the
      ! standard Laplacian operator cannot reach.
      ! Flow-aware biharmonic ν₄ (`nu4_face_*`) is filled either by the
      ! strain-rate Smagorinsky_AH (`smag_ah_active`) or by the 2-D Leith
      ! biharmonic closure (`LMIX_LEITH_BIHARM`); both engage the same
      ! per-face biharmonic apply.  `dt_local` (hoisted above) drives the
      ! per-cell biharmonic CFL clamp.
      if (present(lateral_mix)) then
         if (lateral_mix%is_init .and. &
             (lateral_mix%smag_ah_active .or. &
              lateral_mix%closure == LMIX_LEITH_BIHARM)) then
            ! z-level closed faces: written twice, as for the harmonic
            ! arms above (the `(1,1,1)` placeholder masks must never reach
            ! the explicit-shape dummies).
            if (metrics%use_closed_faces) then
               call hvisc_compute_biharmonic_face_impl( &
                  u, v, &
                  this%lap_u%data, this%lap_v%data, &
                  this%du_visc%data, this%dv_visc%data, &
                  lateral_mix%nu4_face_x, lateral_mix%nu4_face_y, &
                  this%bound_coef, dt_local, &
                  metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, metrics%wet_q, &
                  metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
                  metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
                  nx, ny, nz, metrics%open_u, metrics%open_v)
            else
               call hvisc_compute_biharmonic_face_impl( &
                  u, v, &
                  this%lap_u%data, this%lap_v%data, &
                  this%du_visc%data, this%dv_visc%data, &
                  lateral_mix%nu4_face_x, lateral_mix%nu4_face_y, &
                  this%bound_coef, dt_local, &
                  metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, metrics%wet_q, &
                  metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
                  metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
                  nx, ny, nz)
            end if
            return
         end if
      end if
      if (nu_4 > 0.0_wp) then
         if (metrics%use_closed_faces) then
            call hvisc_compute_biharmonic_impl( &
               u, v, &
               this%lap_u%data, this%lap_v%data, &
               this%du_visc%data, this%dv_visc%data, &
               nu_4, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, metrics%wet_q, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               nx, ny, nz, metrics%open_u, metrics%open_v)
         else
            call hvisc_compute_biharmonic_impl( &
               u, v, &
               this%lap_u%data, this%lap_v%data, &
               this%du_visc%data, this%dv_visc%data, &
               nu_4, this%bound_coef, dt_local, &
               metrics%idxCu, metrics%idyCu, metrics%idxCv, metrics%idyCv, metrics%wet_q, &
               metrics%dy_dxT, metrics%dx_dyBu, metrics%iareaCu, &
               metrics%dx_dyT, metrics%dy_dxBu, metrics%iareaCv, &
               nx, ny, nz)
         end if
      end if
   end subroutine ocean_horizontal_viscosity_compute_tendencies_on