Divergence_TwoPointVector2D_t Subroutine

  subroutine Divergence_TwoPointVector2D_t(this,df)
    !! Computes the split-form (two-point) divergence of a 2-D vector field
    !! in the reference element (computational coordinates).
    !!
    !! The split-form divergence at node (i,j) is
    !!
    !!   (nabla.F)_{i,j} = 2 sum_n [ D_{n,i} F^1(n,i,j) + D_{n,j} F^2(n,i,j) ]
    !!
    !! where D is the standard derivative matrix (dMatrix) and
    !! F^idir(n,i,j,...) = interior(n,i,j,iEl,iVar,idir) stores the two-point
    !! flux between nodes i (or j) and n in the idir-th coordinate direction.
    !!
    !! The interior array is assumed to hold contravariant (J-scaled) two-point
    !! fluxes.  To obtain the physical divergence, divide the result by the
    !! element Jacobian J(i,j,iEl).
    implicit none
    class(TwoPointVector2D_t),intent(in) :: this
    real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar)
    ! Local
    integer :: i,j,nn,iEl,iVar
    real(prec) :: dfLoc

    do concurrent(i=1:this%N+1,j=1:this%N+1,iEl=1:this%nElem,iVar=1:this%nVar)

      dfLoc = 0.0_prec
      do nn = 1,this%N+1
        dfLoc = dfLoc+ &
                this%interp%dSplitMatrix(nn,i)*this%interior(nn,i,j,iEl,iVar,1)+ &
                this%interp%dSplitMatrix(nn,j)*this%interior(nn,i,j,iEl,iVar,2)
      enddo
      df(i,j,iEl,iVar) = 2.0_prec*dfLoc

    enddo

  endsubroutine Divergence_TwoPointVector2D_t