VectorSurfaceSampling Subroutine

      subroutine VectorSurfaceSampling(mesh, zoneID, integralType, monitorName, data_out)
		 use MonitorDefinitions 
#ifdef _HAS_MPI_
         use mpi
#endif
         implicit none
         class(HexMesh),      intent(inout), target :: mesh
         integer,             intent(in)    		:: zoneID
         integer,             intent(in)    		:: integralType
		 real(kind=RP), allocatable, intent(out)  	:: data_out(:)

!
!        ---------------
!        Local variables
!        ---------------
!
         integer  							:: zonefID, fID, eID, fIDs(6), ierr, Nx, Ny, fsID
         class(Element), pointer  			:: elements(:)
		 real(kind=RP) , allocatable        :: data_proc(:)
		 logical 				  			:: file_exists
		 character(len=STR_LEN_MONITORS) 	:: fileName
		 character(len=STR_LEN_MONITORS) 	:: monitorName
		 
		 !
!        Get the number of Order and the sampling interval
!        -------------------------------------------------
		 if (mesh % zones(zoneID) % no_of_faces .gt.0) then
			Nx = mesh % faces (mesh % zones(zoneID) % faces(1))%Nf(1)+1
			Ny = mesh % faces (mesh % zones(zoneID) % faces(1))%Nf(2)+1
		 else
			Nx=0
			Ny=0
		 end if 
		 
		 ALLOCATE (data_out(Nx*Ny*mesh % zones(zoneID) % no_of_faces), data_proc(Nx*Ny))
		 
!
!        *************************
!        Perform the interpolation
!        *************************
!
#if defined(NAVIERSTOKES) && (!(INCNS))
         elements => mesh % elements
!$omp parallel private(fID, eID, fIDs,data_proc) shared(elements,mesh,NodalStorage,zoneID,integralType,&
!$omp&                                        computeGradients,data_out,Nx)
!$omp single
         do zonefID = 1, mesh % zones(zoneID) % no_of_faces
            fID = mesh % zones(zoneID) % faces(zonefID)

            eID = mesh % faces(fID) % elementIDs(1)
            fIDs = mesh % elements(eID) % faceIDs

!$omp task depend(inout:elements(eID))
            call elements(eID) % ProlongSolutionToFaces(NCONS, mesh % faces(fIDs(1)),&
                                            mesh % faces(fIDs(2)),&
                                            mesh % faces(fIDs(3)),&
                                            mesh % faces(fIDs(4)),&
                                            mesh % faces(fIDs(5)),&
                                            mesh % faces(fIDs(6)) )
            if ( computeGradients ) then
               call elements(eID) % ProlongGradientsToFaces(NGRAD, mesh % faces(fIDs(1)),&
                                                mesh % faces(fIDs(2)),&
                                                mesh % faces(fIDs(3)),&
                                                mesh % faces(fIDs(4)),&
                                                mesh % faces(fIDs(5)),&
                                                mesh % faces(fIDs(6)) )
            end if
!$omp end task
         end do
!$omp end single 
		
!
!        Loop the zone to get faces and elements
!        ---------------------------------------
!$omp do  private(fID,data_proc) schedule(runtime)
         do zonefID = 1, mesh % zones(zoneID) % no_of_faces         
!
!           Face global ID
!           --------------
            fID = mesh % zones(zoneID) % faces(zonefID)
!
!           Compute the integral
!           --------------------
			
            CALL VectorSurfaceSampling_Face(mesh % faces(fID), integralType, fID, data_proc)

			data_out((zonefID-1)*Nx*Ny+1:zonefID*Nx*Ny)=data_proc

         end do
		 
!$omp end do 
!$omp end parallel

#endif

      end subroutine VectorSurfaceSampling