evp_u_momentum_impl Subroutine

   pure subroutine evp_u_momentum_impl(idxCu, idyCu, dy2h, dx2q, iareaCu, mask_u, mi_u, mi_v, &
                                       q, str_d, str_t, str_s, uo, vo, tau_ax, ui, vi, &
                                       fxoc, m_neglect, i_cdrhodt, cdrho, dt, &
                                       nx_phys, ny_phys, nghost, nx, ny, a_u, a_face_on)
      !! u-momentum (:1172-1231, requirement 1: fxic_now carries the FULL
      !! str_t force term). Loop over u-faces `ng+1..ng+nxp+1` x
      !! `ng+1..ng+nyp` — each iteration writes only its own face.
      !!
      !! PR 62 (`a_face_on`): weights BOTH the wind (`tau_ax`) AND the
      !! ice-ocean drag (`drag_u`) by the face ice concentration `a_u`, in
      !! the momentum balance ONLY — `fxoc` stays unweighted (per unit ice
      !! area) so `ice_ocean_stress_flux_impl`'s `a_u*fxoc` on the coupler
      !! side is the ocean's share (see the module docstring D7 + the
      !! `ice_ocean_stress_flux` F5 caveat). Weighting the wind alone would
      !! convert today's leak (zero at steady free drift) into a permanent
      !! one — do not "simplify" this to a single weighted term.
      !!
      !! The `a_fac > 0.0` branch is a MANDATORY 0/0 guard, not defensive
      !! tidiness: at an ice-free face `mi_u = 0`, so a naive
      !! `a_fac*drag_u` collapses the denominator to `m_neglect` alone
      !! against a generally-nonzero `dt*fxic_now`, producing `O(1e30)` on
      !! the first substep. The `else` branch (`uio_c = 0` => `ui = uo`) is
      !! the SIS2 limit (`set_wind_stresses_C`'s `else WindStr_x_Cu = 0.0`)
      !! reached without the division hazard. Do NOT floor `a_fac` instead
      !! of branching — a `max(a_fac, eps)` floor reintroduces a
      !! (much smaller but nonzero) ghost-drift artefact.
      !!
      !! The drag PREDICTOR (`b_vel0`/`uio_pred`, below) is deliberately NOT
      !! folded by `a_u`: conservation depends only on `drag_u`'s use in the
      !! `uio_c`/`fxoc` pair (§3.3 of the PR-62 plan), not on the predictor's
      !! accuracy, and `drag_u`'s own `max(uio_init**2, ...)` converges to
      !! the exact quadratic drag as the substep loop converges regardless.
      integer, intent(in) :: nx_phys, ny_phys, nghost, nx, ny
      real(wp), intent(in) :: idxCu(nx + 1, ny), idyCu(nx + 1, ny)
      real(wp), intent(in) :: dy2h(nx, ny), dx2q(nx + 1, ny + 1)
      real(wp), intent(in) :: iareaCu(nx + 1, ny)
      real(wp), intent(in) :: mask_u(nx + 1, ny)
      real(wp), intent(in) :: mi_u(nx + 1, ny), mi_v(nx, ny + 1)
      real(wp), intent(in) :: q(nx + 1, ny + 1)
      real(wp), intent(in) :: str_d(nx, ny), str_t(nx, ny), str_s(nx + 1, ny + 1)
      real(wp), intent(in) :: uo(nx + 1, ny), vo(nx, ny + 1)
      real(wp), intent(in) :: tau_ax(nx + 1, ny)
      real(wp), intent(inout) :: ui(nx + 1, ny)
      real(wp), intent(in) :: vi(nx, ny + 1)
      real(wp), intent(inout) :: fxoc(nx + 1, ny)
      real(wp), intent(in) :: m_neglect, i_cdrhodt, cdrho, dt
      real(wp), intent(in) :: a_u(nx + 1, ny)
         !! PR 62: face ice concentration. Valid ONLY when `a_face_on`.
      logical, intent(in) :: a_face_on
      integer :: i, j, i_lo, i_hi, j_lo, j_hi
      real(wp) :: cor, f2dt_u, i1_f2dt2_u
      real(wp) :: azon, bzon, czon, dzon
      real(wp) :: fxic_now, v2_at_u, uio_init
      real(wp) :: m_uio_explicit, b_vel0, uio_pred, drag_u, uio_c
      real(wp) :: a_fac, tau_eff, drag_eff

      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(cor, f2dt_u, i1_f2dt2_u, azon, bzon, czon, dzon, fxic_now, v2_at_u, &
               uio_init, m_uio_explicit, b_vel0, uio_pred, drag_u, uio_c, &
               a_fac, tau_eff, drag_eff)
         azon = 0.25_wp*mi_v(i, j + 1)*q(i, j + 1)
         bzon = 0.25_wp*mi_v(i - 1, j + 1)*q(i, j + 1)
         czon = 0.25_wp*mi_v(i - 1, j)*q(i, j)
         dzon = 0.25_wp*mi_v(i, j)*q(i, j)

         cor = 0.25_wp*(q(i, j + 1)*(mi_v(i, j + 1)*vi(i, j + 1) + mi_v(i - 1, j + 1)*vi(i - 1, j + 1)) + &
                        q(i, j)*(mi_v(i - 1, j)*vi(i - 1, j) + mi_v(i, j)*vi(i, j)))

         f2dt_u = dt*4.0_wp*((azon**2 + czon**2) + (bzon**2 + dzon**2))
         i1_f2dt2_u = 1.0_wp/(1.0_wp + dt*f2dt_u)

         fxic_now = idxCu(i, j)*(str_d(i, j) - str_d(i - 1, j)) + &
                    (idyCu(i, j)*(dy2h(i, j)*str_t(i, j) - dy2h(i - 1, j)*str_t(i - 1, j)) + &
                     idxCu(i, j)*(dx2q(i, j + 1)*str_s(i, j + 1) - dx2q(i, j)*str_s(i, j)))* &
                    iareaCu(i, j)

         v2_at_u = 0.25_wp*(((vi(i - 1, j + 1) - vo(i - 1, j + 1))**2 + &
                             (vi(i, j) - vo(i, j))**2) + &
                            ((vi(i, j + 1) - vo(i, j + 1))**2 + &
                             (vi(i - 1, j) - vo(i - 1, j))**2))

         uio_init = ui(i, j) - uo(i, j)

         ! TWO FULLY SEPARATE ARMS, not `a_fac = 1.0_wp` feeding one shared
         ! expression: the off arm below is TEXTUALLY UNCHANGED from the
         ! pre-PR kernel, predictor included.  A shared `tau_eff`/`drag_eff`
         ! computed via `a_fac = 1.0_wp` is mathematically exact (1.0*x==x
         ! in IEEE) but NVHPC's GPU codegen does not guarantee identical FMA
         ! contraction/rounding across ~5000 chained substeps for two
         ! syntactically different expression trees that merely evaluate to
         ! the same VALUE — confirmed empirically (`a_face_full_cover_
         ! bitident` drifted ~1e-14 rel under exactly that construction
         ! before this fix; the `7d283fe8` GPU FMA precedent, CLAUDE.md
         ! Gotchas).  Keep the off arm untouched, full stop.
         if (a_face_on) then
            a_fac = a_u(i, j)
            tau_eff = a_fac*tau_ax(i, j)

            drag_u = 0.0_wp
            if (mask_u(i, j) > 0.0_wp) then
               m_uio_explicit = uio_init*mi_u(i, j) + dt*(cor*mi_u(i, j) + (fxic_now + tau_eff))
               b_vel0 = mi_u(i, j)*i_cdrhodt + (sqrt(uio_init**2 + v2_at_u) - abs(uio_init))
               if (b_vel0**2 > EVP_DRAG_LINEARIZE_THRESHOLD*i_cdrhodt*abs(m_uio_explicit)) then
                  if (b_vel0 /= 0.0_wp) then
                     uio_pred = m_uio_explicit*i_cdrhodt/b_vel0
                  else
                     uio_pred = 0.0_wp
                  end if
               else
                  uio_pred = 0.5_wp*(sqrt(b_vel0**2 + 4.0_wp*i_cdrhodt*abs(m_uio_explicit)) - b_vel0)
               end if
               drag_u = cdrho*sqrt(max(uio_init**2, uio_pred**2) + v2_at_u)
            end if

            drag_eff = a_fac*drag_u
            if (a_fac > 0.0_wp) then
               uio_c = mask_u(i, j)*(mi_u(i, j)*((ui(i, j) + dt*cor)*i1_f2dt2_u - uo(i, j)) + &
                                     dt*(fxic_now + tau_eff))/ &
                       (mi_u(i, j) + m_neglect + dt*drag_eff)
            else
               uio_c = 0.0_wp     ! no ice at either neighbour: no ice momentum here
            end if
         else
            drag_u = 0.0_wp
            if (mask_u(i, j) > 0.0_wp) then
               m_uio_explicit = uio_init*mi_u(i, j) + dt*(cor*mi_u(i, j) + (fxic_now + tau_ax(i, j)))
               b_vel0 = mi_u(i, j)*i_cdrhodt + (sqrt(uio_init**2 + v2_at_u) - abs(uio_init))
               if (b_vel0**2 > EVP_DRAG_LINEARIZE_THRESHOLD*i_cdrhodt*abs(m_uio_explicit)) then
                  if (b_vel0 /= 0.0_wp) then
                     uio_pred = m_uio_explicit*i_cdrhodt/b_vel0
                  else
                     uio_pred = 0.0_wp
                  end if
               else
                  uio_pred = 0.5_wp*(sqrt(b_vel0**2 + 4.0_wp*i_cdrhodt*abs(m_uio_explicit)) - b_vel0)
               end if
               drag_u = cdrho*sqrt(max(uio_init**2, uio_pred**2) + v2_at_u)
            end if

            uio_c = mask_u(i, j)*(mi_u(i, j)*((ui(i, j) + dt*cor)*i1_f2dt2_u - uo(i, j)) + &
                                  dt*(fxic_now + tau_ax(i, j)))/ &
                    (mi_u(i, j) + m_neglect + dt*drag_u)
         end if

         ui(i, j) = (uio_c + uo(i, j))*mask_u(i, j)
         fxoc(i, j) = fxoc(i, j) + drag_u*uio_c          ! UNWEIGHTED: the coupler applies a_u
      end do
   end subroutine evp_u_momentum_impl