VectorVolumeIntegral Function

      function VectorVolumeIntegral(mesh, integralType, num_of_vars) result(val)
!
!        -----------------------------------------------------------
!           This function computes vector integrals, that is, those
!           in the form:
!                 val = \int \vec{v} dx
!           Implemented integrals are:
!              * VELOCITY
!              * MOMENTUM
!        -----------------------------------------------------------
!
         implicit none
         class(HexMesh),      intent(in)  :: mesh
         integer,             intent(in)  :: integralType
         integer,             intent(in)  :: num_of_vars
         real(kind=RP)                    :: val(num_of_vars)
!
!        ---------------
!        Local variables
!        ---------------
!
         logical       :: fiveVars
         integer       :: eID, ierr
         real(kind=RP) :: localVal(num_of_vars)
         real(kind=RP) :: valAux(num_of_vars)
         real(kind=RP) :: val1, val2, val3, val4, val5

         if (num_of_vars == 5) then    ! Ugly hack.. But only way to make it work with ifort....
            fiveVars = .TRUE.
         else
            fiveVars = .FALSE.
         end if

!
!        Initialization
!        --------------
         val1 = 0.0_RP
         val2 = 0.0_RP
         val3 = 0.0_RP
         val4 = 0.0_RP
         val5 = 0.0_RP
!
!        Loop the mesh
!        -------------
!$omp parallel do reduction(+:val1,val2,val3,val4,val5) private(valAux) schedule(guided)
         do eID = 1, mesh % no_of_elements
!
!           Compute the integral
!           --------------------
            valAux = VectorVolumeIntegral_Local(mesh % elements(eID), integralType  , num_of_vars  )

            val1 = val1 + valAux(1)
            val2 = val2 + valAux(2)
            val3 = val3 + valAux(3)

            if (fiveVars) then
               val4 = val4 + valAux(4)
               val5 = val5 + valAux(5)
            end if
         end do
!$omp end parallel do

         val(1:3) = [val1, val2, val3]

         if (fiveVars) val(4:5) = [val4, val5]

#ifdef _HAS_MPI_
         localVal = val
         call mpi_allreduce(localVal, val, num_of_vars, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
#endif

      end function VectorVolumeIntegral