evp_v_momentum_impl Subroutine

   pure subroutine evp_v_momentum_impl(idyCv, idxCv, dx2h, dy2q, iareaCv, mask_v, mi_v, mi_u, &
                                       q, str_d, str_t, str_s, uo, vo, tau_ay, u_tmp, vi, &
                                       fyoc, m_neglect, i_cdrhodt, cdrho, dt, &
                                       nx_phys, ny_phys, nghost, nx, ny, a_v, a_face_on)
      !! v-momentum (:1257-1334, mirror of u). D4: reads `u_tmp` (the
      !! PRE-update u), never the just-updated `ui`. **Minus** on the
      !! str_t divergence term (:1263-1267).
      !!
      !! PR 62 (`a_face_on`): exact mirror of `evp_u_momentum_impl`'s
      !! weighting — see that kernel's docstring for the full rationale
      !! (both terms weighted, `fyoc` unweighted, the `a_fac > 0.0` 0/0
      !! guard, and the untouched drag predictor).
      integer, intent(in) :: nx_phys, ny_phys, nghost, nx, ny
      real(wp), intent(in) :: idyCv(nx, ny + 1), idxCv(nx, ny + 1)
      real(wp), intent(in) :: dx2h(nx, ny), dy2q(nx + 1, ny + 1)
      real(wp), intent(in) :: iareaCv(nx, ny + 1)
      real(wp), intent(in) :: mask_v(nx, ny + 1)
      real(wp), intent(in) :: mi_v(nx, ny + 1), mi_u(nx + 1, ny)
      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_ay(nx, ny + 1)
      real(wp), intent(in) :: u_tmp(nx + 1, ny)
      real(wp), intent(inout) :: vi(nx, ny + 1)
      real(wp), intent(inout) :: fyoc(nx, ny + 1)
      real(wp), intent(in) :: m_neglect, i_cdrhodt, cdrho, dt
      real(wp), intent(in) :: a_v(nx, ny + 1)
         !! 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_v, i1_f2dt2_v
      real(wp) :: amer, bmer, cmer, dmer
      real(wp) :: fyic_now, u2_at_v, vio_init
      real(wp) :: m_vio_explicit, b_vel0, vio_pred, drag_v, vio_c
      real(wp) :: a_fac, tau_eff, drag_eff

      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(cor, f2dt_v, i1_f2dt2_v, amer, bmer, cmer, dmer, fyic_now, u2_at_v, &
               vio_init, m_vio_explicit, b_vel0, vio_pred, drag_v, vio_c, &
               a_fac, tau_eff, drag_eff)
         amer = 0.25_wp*mi_u(i, j - 1)*q(i, j)
         bmer = 0.25_wp*mi_u(i + 1, j - 1)*q(i + 1, j)
         cmer = 0.25_wp*mi_u(i + 1, j)*q(i + 1, j)
         dmer = 0.25_wp*mi_u(i, j)*q(i, j)

         cor = -0.25_wp*(q(i, j)*(mi_u(i, j - 1)*u_tmp(i, j - 1) + mi_u(i, j)*u_tmp(i, j)) + &
                         q(i + 1, j)*(mi_u(i + 1, j - 1)*u_tmp(i + 1, j - 1) + &
                                      mi_u(i + 1, j)*u_tmp(i + 1, j)))

         f2dt_v = dt*4.0_wp*((amer**2 + cmer**2) + (bmer**2 + dmer**2))
         i1_f2dt2_v = 1.0_wp/(1.0_wp + dt*f2dt_v)

         fyic_now = idyCv(i, j)*(str_d(i, j) - str_d(i, j - 1)) + &
                    (-idxCv(i, j)*(dx2h(i, j)*str_t(i, j) - dx2h(i, j - 1)*str_t(i, j - 1)) + &
                     idyCv(i, j)*(dy2q(i + 1, j)*str_s(i + 1, j) - dy2q(i, j)*str_s(i, j)))* &
                    iareaCv(i, j)

         u2_at_v = 0.25_wp*(((u_tmp(i + 1, j - 1) - uo(i + 1, j - 1))**2 + &
                             (u_tmp(i, j) - uo(i, j))**2) + &
                            ((u_tmp(i + 1, j) - uo(i + 1, j))**2 + &
                             (u_tmp(i, j - 1) - uo(i, j - 1))**2))

         vio_init = vi(i, j) - vo(i, j)

         ! TWO FULLY SEPARATE ARMS -- see evp_u_momentum_impl's comment at
         ! the mirror site for why the off arm (predictor included) must
         ! stay textually untouched rather than routed through a shared
         ! `a_fac = 1.0_wp` expression.
         if (a_face_on) then
            a_fac = a_v(i, j)
            tau_eff = a_fac*tau_ay(i, j)

            drag_v = 0.0_wp
            if (mask_v(i, j) > 0.0_wp) then
               m_vio_explicit = vio_init*mi_v(i, j) + dt*(cor*mi_v(i, j) + (fyic_now + tau_eff))
               b_vel0 = mi_v(i, j)*i_cdrhodt + (sqrt(vio_init**2 + u2_at_v) - abs(vio_init))
               if (b_vel0**2 > EVP_DRAG_LINEARIZE_THRESHOLD*i_cdrhodt*abs(m_vio_explicit)) then
                  if (b_vel0 /= 0.0_wp) then
                     vio_pred = m_vio_explicit*i_cdrhodt/b_vel0
                  else
                     vio_pred = 0.0_wp
                  end if
               else
                  vio_pred = 0.5_wp*(sqrt(b_vel0**2 + 4.0_wp*i_cdrhodt*abs(m_vio_explicit)) - b_vel0)
               end if
               drag_v = cdrho*sqrt(max(vio_init**2, vio_pred**2) + u2_at_v)
            end if

            drag_eff = a_fac*drag_v
            if (a_fac > 0.0_wp) then
               vio_c = mask_v(i, j)*(mi_v(i, j)*((vi(i, j) + dt*cor)*i1_f2dt2_v - vo(i, j)) + &
                                     dt*(fyic_now + tau_eff))/ &
                       (mi_v(i, j) + m_neglect + dt*drag_eff)
            else
               vio_c = 0.0_wp     ! no ice at either neighbour: no ice momentum here
            end if
         else
            drag_v = 0.0_wp
            if (mask_v(i, j) > 0.0_wp) then
               m_vio_explicit = vio_init*mi_v(i, j) + dt*(cor*mi_v(i, j) + (fyic_now + tau_ay(i, j)))
               b_vel0 = mi_v(i, j)*i_cdrhodt + (sqrt(vio_init**2 + u2_at_v) - abs(vio_init))
               if (b_vel0**2 > EVP_DRAG_LINEARIZE_THRESHOLD*i_cdrhodt*abs(m_vio_explicit)) then
                  if (b_vel0 /= 0.0_wp) then
                     vio_pred = m_vio_explicit*i_cdrhodt/b_vel0
                  else
                     vio_pred = 0.0_wp
                  end if
               else
                  vio_pred = 0.5_wp*(sqrt(b_vel0**2 + 4.0_wp*i_cdrhodt*abs(m_vio_explicit)) - b_vel0)
               end if
               drag_v = cdrho*sqrt(max(vio_init**2, vio_pred**2) + u2_at_v)
            end if

            vio_c = mask_v(i, j)*(mi_v(i, j)*((vi(i, j) + dt*cor)*i1_f2dt2_v - vo(i, j)) + &
                                  dt*(fyic_now + tau_ay(i, j)))/ &
                    (mi_v(i, j) + m_neglect + dt*drag_v)
         end if

         vi(i, j) = (vio_c + vo(i, j))*mask_v(i, j)
         fyoc(i, j) = fyoc(i, j) + drag_v*vio_c          ! UNWEIGHTED: the coupler applies a_v
      end do
   end subroutine evp_v_momentum_impl