ocean_lateral_mix_compute_vel_scale Subroutine

   pure subroutine ocean_lateral_mix_compute_vel_scale(grid, metrics, this, ms, seed_bg)
      !! Public only for the unit-test suite; ignore in production code.
      !! Live velocity-scale viscosity (MOM6 `KH_VEL_SCALE`, Kh = U·Δ):
      !! per face `A_vel = kh_vel_scale_live · L_grid · |u_face|`
      !! (`L_grid = sqrt(dxT·dyT)`), `max`-combined into `ah_face_*` so it
      !! floors — never reduces — the active closure.  `seed_bg = .true.`
      !! first fills every face with `ah_bg` (used when no closure ran);
      !! `.false.` only raises faces where `A_vel` exceeds the closure.
      !! No-op when `kh_vel_scale_live <= 0`.
      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
      logical, intent(in) :: seed_bg

      integer :: i, j, k, nx, ny, nz
      real(wp) :: vel_scale_local, ah_bg_local, ah_max_local
      real(wp) :: a_vel, l_grid

      if (.not. this%is_init) return
      if (this%kh_vel_scale_live <= 0.0_wp) 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
      vel_scale_local = this%kh_vel_scale_live
      ah_bg_local = this%ah_bg
      ah_max_local = this%ah_max

      ! ---- u-faces (i-1/2, j): A_vel from |u_face_x| ----
      ! Interior faces i=2:nx (the closures' range); L_grid taken at
      ! the west-adjacent T cell (i-1).  Cap and floor as the closures
      ! do.  Wall faces (i=1, i=nx+1) get the background under seeding.
      do concurrent(k=1:nz, j=1:ny, i=2:nx) local(a_vel, l_grid)
         l_grid = sqrt(metrics%dxT(i - 1, j)*metrics%dyT(i - 1, j))
         a_vel = min(ah_max_local, vel_scale_local*l_grid*abs(ms%u_face_x_layer(i, j, k)))
         if (seed_bg) then
            this%ah_face_x(i, j, k) = max(ah_bg_local, a_vel)
         else
            this%ah_face_x(i, j, k) = max(this%ah_face_x(i, j, k), a_vel)
         end if
      end do
      do concurrent(k=1:nz, j=1:ny)
         if (seed_bg) then
            this%ah_face_x(1, j, k) = ah_bg_local
            this%ah_face_x(nx + 1, j, k) = ah_bg_local
         end if
      end do

      ! ---- v-faces (i, j-1/2): A_vel from |v_face_y| ----
      ! Interior faces j=2:ny; L_grid at the south-adjacent T cell.
      do concurrent(k=1:nz, j=2:ny, i=1:nx) local(a_vel, l_grid)
         l_grid = sqrt(metrics%dxT(i, j - 1)*metrics%dyT(i, j - 1))
         a_vel = min(ah_max_local, vel_scale_local*l_grid*abs(ms%v_face_y_layer(i, j, k)))
         if (seed_bg) then
            this%ah_face_y(i, j, k) = max(ah_bg_local, a_vel)
         else
            this%ah_face_y(i, j, k) = max(this%ah_face_y(i, j, k), a_vel)
         end if
      end do
      do concurrent(k=1:nz, i=1:nx)
         if (seed_bg) then
            this%ah_face_y(i, 1, k) = ah_bg_local
            this%ah_face_y(i, ny + 1, k) = ah_bg_local
         end if
      end do
   end subroutine ocean_lateral_mix_compute_vel_scale