evp_truncate_velocity_impl Subroutine

   pure subroutine evp_truncate_velocity_impl(areaT, dy_cu, dx_cv, ui, vi, cfl_trunc, dt_tr, &
                                              backoff, nghost, nx_phys, ny_phys, nx, ny)
      !! PR 36: the shared CFL-clip algebra -- the transport-CFL bound on
      !! the ice velocity (SIS2 `SIS_dyn_cgrid.F90:839-870`, the in-loop
      !! half at `:1338-1361`, the final half at `:1443-1500`; this
      !! routine is the counting-free, caller-chosen-backoff form both
      !! reuse; `evp_truncate_final_impl` below wraps it with the 0.95
      !! back-off and the `mi > m_neglect` count).
      !!
      !! `u_max(face) = +cfl_trunc*areaT(donor for u>0)/(dt_tr*dy_cu(face))`,
      !! `u_min(face) = -cfl_trunc*areaT(donor for u<0)/(dt_tr*dy_cu(face))`
      !! -- "the flux out of a cell in one slow step cannot exceed
      !! cfl_trunc of its volume". The donor asymmetry is load-bearing:
      !! `+u` at rdb u-face `(i,j)` (the WEST face of cell `(i,j)`,
      !! module docstring §1) drains the WEST cell `(i-1,j)`; `-u` drains
      !! the EAST cell `(i,j)`. v-mirror: `+v` drains the SOUTH cell
      !! `(i,j-1)`, `-v` drains the NORTH cell `(i,j)`.
      !!
      !! `dy_cu`/`dx_cv` (NOT the unmasked `dyCu`/`dxCv`) are the
      !! topography-aware OPEN face widths -- zero at a closed/land face,
      !! which is why the bound is guarded `> 0.0`: a closed face gets
      !! `u_hi = u_lo = 0` (forces `ui = 0` there), not a finite spurious
      !! bound from dividing by a nonzero length at land.
      !!
      !! Loop ranges are copied VERBATIM from `evp_u_momentum_impl` (u) and
      !! `evp_v_momentum_impl` (v) -- physical faces only. Over that range
      !! `i-1 >= nghost >= 1` (resp. `j-1 >= nghost >= 1`) always, so no
      !! array-edge branch is needed (unlike `evp_mi_face_impl`, which
      !! loops the full `1:nx+1`/`1:ny+1` and does need one).
      integer, intent(in) :: nghost, nx_phys, ny_phys, nx, ny
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: dy_cu(nx + 1, ny), dx_cv(nx, ny + 1)
      real(wp), intent(inout) :: ui(nx + 1, ny), vi(nx, ny + 1)
      real(wp), intent(in) :: cfl_trunc, dt_tr, backoff
      integer :: i, j, i_lo, i_hi, j_lo, j_hi
      real(wp) :: u_hi, u_lo, v_hi, v_lo, loc_scale

      i_lo = nghost + 1
      i_hi = nghost + nx_phys + 1
      j_lo = nghost + 1
      j_hi = nghost + ny_phys
      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) local(u_hi, u_lo, loc_scale)
         u_hi = 0.0_wp
         u_lo = 0.0_wp
         if (dy_cu(i, j) > 0.0_wp) then
            loc_scale = cfl_trunc/(dt_tr*dy_cu(i, j))
            u_hi = backoff*loc_scale*areaT(i - 1, j)
            u_lo = -backoff*loc_scale*areaT(i, j)
         end if
         if (ui(i, j) > u_hi) then
            ui(i, j) = u_hi
         else if (ui(i, j) < u_lo) then
            ui(i, j) = u_lo
         end if
      end do

      i_lo = nghost + 1
      i_hi = nghost + nx_phys
      j_lo = nghost + 1
      j_hi = nghost + ny_phys + 1
      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) local(v_hi, v_lo, loc_scale)
         v_hi = 0.0_wp
         v_lo = 0.0_wp
         if (dx_cv(i, j) > 0.0_wp) then
            loc_scale = cfl_trunc/(dt_tr*dx_cv(i, j))
            v_hi = backoff*loc_scale*areaT(i, j - 1)
            v_lo = -backoff*loc_scale*areaT(i, j)
         end if
         if (vi(i, j) > v_hi) then
            vi(i, j) = v_hi
         else if (vi(i, j) < v_lo) then
            vi(i, j) = v_lo
         end if
      end do
   end subroutine evp_truncate_velocity_impl