Divergence_Vector3D_t Subroutine

  subroutine Divergence_Vector3D_t(this,df)
    implicit none
    class(Vector3D_t),intent(in) :: this
    real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar)
    ! Local
    integer    :: i,j,k,ii,iel,ivar
    real(prec) :: dfLoc

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar)

      dfLoc = 0.0_prec
      do ii = 1,this%N+1
        dfLoc = dfLoc+this%interp%dMatrix(ii,i)*this%interior(ii,j,k,iel,ivar,1)
      enddo
      dF(i,j,k,iel,ivar) = dfLoc

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar)

      dfLoc = 0.0_prec
      do ii = 1,this%N+1
        dfLoc = dfLoc+this%interp%dMatrix(ii,j)*this%interior(i,ii,k,iel,ivar,2)
      enddo
      dF(i,j,k,iel,ivar) = dF(i,j,k,iel,ivar)+dfLoc

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar)

      dfLoc = 0.0_prec
      do ii = 1,this%N+1
        dfLoc = dfLoc+this%interp%dMatrix(ii,k)*this%interior(i,j,ii,iel,ivar,3)
      enddo
      dF(i,j,k,iel,ivar) = dF(i,j,k,iel,ivar)+dfLoc

    enddo

  endsubroutine Divergence_Vector3D_t