#include "Includes.h"

module PlaneSampling
   use SMConstants
   use HexMeshClass
   use MonitorDefinitions
   use PhysicsStorage
   use VariableConversion
   use MPI_Process_Info
   use FluidData
   use FileReadingUtilities, only: getRealArrayFromString, GetRealValue, getCharArrayFromString
   use NodalStorageClass   , only: NodalStorage
#ifdef _HAS_MPI_
   use mpi
#endif
   implicit none
   
   private
   public   PlaneSampling_t
!
!  *******************************
!  PlaneSamplings class definition
!  *******************************
!

   type PlaneSampling_t
      logical      , allocatable      :: active(:)
      integer      , allocatable      :: rank (:)
      integer                         :: ID
      integer                         :: nVariables
      integer                         :: interval
      integer                         :: bufferSize
      integer                         :: bufferLine
      integer                         :: intervalCount
      integer                         :: N(2)
      integer                         :: nNodes
      integer      , allocatable      :: eID (:)
      real(kind=RP), allocatable      :: lxi(:,:) , leta(:,:), lzeta(:,:)
      real(kind=RP), allocatable      :: values(:,:,:)
      real(kind=RP), allocatable      :: x(:,:)
      real(kind=RP), allocatable      :: xi(:,:)
      logical                         :: disturbanceData =.false.
      character(len=STR_LEN_MONITORS),  allocatable :: fileName (:)
      character(len=STR_LEN_MONITORS)               :: planeName
      character(len=STR_LEN_MONITORS)               :: fileInput
      character(len=STR_LEN_MONITORS),  allocatable :: variable (:)
      contains
         procedure   :: Initialization => Plane_Initialization
         procedure   :: Update         => Plane_Update
		 procedure   :: UpdateInterp   => Plane_UpdateLagrangeInterp
         procedure   :: WriteToFile    => Plane_WriteToFile
         procedure   :: LookInOtherPartitions => Plane_LookInOtherPartitions
         procedure   :: destruct       => Plane_Destruct
         procedure   :: copy           => Plane_Assign
         generic     :: assignment(=)  => copy
   end type PlaneSampling_t

   contains

      subroutine Plane_Initialization(self, mesh, ID, solution_file, FirstCall)
         use Headers
         use ParamfileRegions
         use MPI_Process_Info
         use Utilities, only: toLower
         implicit none
         class(PlaneSampling_t)  :: self
         class(HexMesh)          :: mesh
         integer                 :: ID
         character(len=*)        :: solution_file
         logical, intent(in)     :: FirstCall
!
!        ---------------
!        Local variables
!        ---------------
!
         integer                          :: i, j, k, l, fid, nNodes,inputType, ierr
         logical                          :: firstOutDomain
         real(kind=RP)                    :: point_1(NDIM,1),point_2(NDIM,1),point_3(NDIM,1)
         real(kind=RP), allocatable       :: pStart(:,:), pEnd(:,:)
         real(kind=RP)                    :: x(NDIM), xi(NDIM)
         real(kind=RP)                    :: delta(NDIM,2), invNodes(2)
         real(kind=RP), allocatable       :: deltaArray1(:,:)
         character(len=STR_LEN_MONITORS)  :: in_label
         character(len=STR_LEN_MONITORS)  :: fileName
         character(len=STR_LEN_MONITORS)  :: paramFile
         character(len=STR_LEN_MONITORS)  :: interval
         character(len=STR_LEN_MONITORS)  :: point1_char,point2_char,point3_char
         character(len=STR_LEN_MONITORS)  :: N12,N23, inputType_char
         character(len=STR_LEN_MONITORS)  :: variables
         character(len=STR_LEN_MONITORS)  :: fileFormat
         character(len=STR_LEN_MONITORS)  :: writeInterval
         character(len=STR_LEN_MONITORS)  :: distData
         
         firstOutDomain = .FALSE.
         
         if (FirstCall) then
!
!           Get monitor ID, assign zero to bufferLine and intervalCount
!           --------------
            self % ID = ID
            self % bufferLine = 0
            self % intervalCount = 0
!
!           Search for the parameters in the case file
!           ------------------------------------------
            write(in_label , '(A,I0)') "#define plane sampling " , self % ID
            
            call get_command_argument(1, paramFile)
            call readValueInRegion(trim(paramFile), "name"           , self % planeName  , in_label, "#end" )
            call readValueInRegion(trim(paramFile), "variables"      , variables         , in_label, "#end" )
            call readValueInRegion(trim(paramFile), "node input type", inputType_char    , in_label, "#end" )
			if (len(trim(inputType_char))>0) then
				inputType  = GetRealValue(inputType_char)
			else
				inputType  = 0
			end if 
            if (inputType.ne.0) then
                call readValueInRegion(trim(paramFile), "input file"           , self % fileInput  , in_label, "#end" )
            else
                call readValueInRegion(trim(paramFile), "point 1 [x,y,z]", point1_char       , in_label, "#end" )
                call readValueInRegion(trim(paramFile), "point 2 [x,y,z]", point2_char       , in_label, "#end" )
                call readValueInRegion(trim(paramFile), "point 3 [x,y,z]", point3_char       , in_label, "#end" )
                call readValueInRegion(trim(paramFile), "n12"            , N12               , in_label, "#end" )
                call readValueInRegion(trim(paramFile), "n23"            , N23               , in_label, "#end" )
            end if 
            call readValueInRegion(trim(paramFile), "sampling interval", interval        , in_label, "#end" )
            call readValueInRegion(trim(paramFile), "write interval", writeInterval      , in_label, "#end" )
            call readValueInRegion(trim(paramFile), "data type", distData, in_label, "#end" )
            
!
!           Get the variables, points, N discretization, and interval
!           ----------------------------------------------
            call getCharArrayFromString(variables, STR_LEN_MONITORS, self % variable)
            self % nVariables = size(self % variable)
            if (inputType.eq.0) then
                point_1(:,1) = getRealArrayFromString(point1_char)
                point_2(:,1) = getRealArrayFromString(point2_char)
                point_3(:,1) = getRealArrayFromString(point3_char)
                self % N(1)  = GetRealValue(N12)
                self % N(2)  = GetRealValue(N23)
            else 
!               Open file input for nodes
                open(newunit=fid, file= trim(self % fileInput),status="old",action="read", iostat=ierr)
                if (ierr /= 0) then
                    write(*,*) "ERROR-opening input file for plane sampling: ", trim(self % fileInput)
                    stop
                end if
                read(fid, *, iostat=ierr) self % N(1), self % N(2)
            end if 
            self % interval = GetRealValue(interval)
            
            
            interval=TRIM(ADJUSTL(interval))
            writeInterval=TRIM(ADJUSTL(writeInterval))

!
!           Failsafe if the number of node is less than 2
!           ---------------------------------------------
            IF ((self % N(1).LT.2)) THEN
                self % N(1)=2
            END IF
            IF ((self % N(2).LT.2)) THEN
                self % N(2)=2
            END IF
            
            self % nNodes = self % N(1) * self % N(2)
            nNodes = self % nNodes
!
!           Get the max. number of timestep in the buffer file before being written
!           -----------------------------------------------------------------------         
            IF (LEN(TRIM(writeInterval)) .EQ. 0) THEN
                self % bufferSize = 1;
            ELSE
                self % bufferSize = GetRealValue(writeInterval)
!
!               Failsafe to prevent too many data being written at one time
!               -----------------------------------------------------------             
                IF (nNodes * self % bufferSize .GT. 500000) THEN   
                    self % bufferSize = 1;
                END IF
            END IF
!
!           Allocate Variables
!           ------------------          
            ALLOCATE(self % x (NDIM, nNodes), self % xi (NDIM, nNodes), self % active (nNodes) &
                    , self % eID (nNodes), self % rank (nNodes), self % fileName (self % nVariables) &
                    , self % values(nNodes+1,self % bufferSize,self % nVariables), deltaArray1(NDIM,self % N(1))  &
                    , pStart(NDIM, self % N(1)))
!
!           Allocate Lagrange interpolants - Assumed identical polynomial for all elements
!           -----------------------------------------------------------------------------
            ALLOCATE(self % lxi( 0 : mesh % elements (1) % Nxyz(1), 1:nNodes), &
                        self % leta( 0 : mesh % elements (1) % Nxyz(2), 1:nNodes), &
                        self % lzeta( 0 : mesh % elements (1) % Nxyz(3), 1:nNodes))
!
!           Nodes spacing
!           -------------           
            invNodes(1)=1.0_RP/(self % N(1)-1)
            invNodes(2)=1.0_RP/(self % N(2)-1)
            if (inputType.eq.0) then 
                delta(1,1)=(point_2(1,1)-point_1(1,1))*invNodes(1)
                delta(2,1)=(point_2(2,1)-point_1(2,1))*invNodes(1)
                delta(3,1)=(point_2(3,1)-point_1(3,1))*invNodes(1)
                delta(1,2)=(point_3(1,1)-point_2(1,1))*invNodes(2)
                delta(2,2)=(point_3(2,1)-point_2(2,1))*invNodes(2)
                delta(3,2)=(point_3(3,1)-point_2(3,1))*invNodes(2)
!
!           Create spacing Array for Nodes in direction 1
!           ---------------------------------------------               
                DO l=1,self % N(1)
                    deltaArray1(1,l)=(delta(1,1))*(l-1)
                    deltaArray1(2,l)=(delta(2,1))*(l-1)
                    deltaArray1(3,l)=(delta(3,1))*(l-1)
                END DO
!
!           Get Nodes location of the plane
!           -------------------------------             
                DO l=1,self % N(2) 
                    DO k=1, self % N(1)
                        DO j=1,3
                            pStart(j,k) = (point_1(j,1)+delta(j,2)*(l-1))
                            self % x(j,(l-1)*self % N(1)+k) = pStart(j,k) + deltaArray1(j,k)
                        END DO
                    END DO 
                END DO
            else  
                ALLOCATE(pEnd(NDIM, self % N(1)))
!
!           Read first points location from input file              
!           ----------------------------------      
                DO l=1,self % N(1)
                    read(fid, *, iostat=ierr) pStart(1,l), pStart(2,l), pStart(3,l)
                    if (ierr /= 0) exit
                END DO 
!
!           Read last points location from input file               
!           ----------------------------------      
                DO l=1,self % N(1)
                    read(fid, *, iostat=ierr) pEnd(1,l), pEnd(2,l), pEnd(3,l)
                    if (ierr /= 0) exit
                END DO 
                close(fid)
!
!           Create spacing Array for Nodes in sweep direction
!           -------------------------------------------------               
                DO l=1,self % N(1)
                    deltaArray1(1,l)=(pEnd(1,l)-pStart(1,l))*invNodes(2)
                    deltaArray1(2,l)=(pEnd(2,l)-pStart(2,l))*invNodes(2)
                    deltaArray1(3,l)=(pEnd(3,l)-pStart(3,l))*invNodes(2)
                END DO
!
!           Get Nodes location of the plane
!           -------------------------------             
                DO l=1,self % N(2) 
                    DO j=1,NDIM
                        DO k=1, self % N(1)
                            self % x(j,(l-1)*self % N(1)+k) = pStart(j,k) + deltaArray1(j,k)*(l-1)
                        END DO
                    END DO 
                END DO       
                DEALLOCATE(pEnd)				
            end if 

            
            DEALLOCATE (deltaArray1,pStart)

!
!           Find node location in the element ID
!           ------------------------------------                
            DO i=1,self % nNodes
        
                x(1)=self % x(1,i)
                x(2)=self % x(2,i)
                x(3)=self % x(3,i)
!
!               Find the requested point in the mesh
!               ------------------------------------
                self % active (i) = mesh % FindPointWithCoords(x, self % eID(i), xi)
                
                self % xi(1,i) = xi(1)
                self % xi(2,i) = xi(2)
                self % xi(3,i) = xi(3)
!
!               Check whether the Plane is located in other partition
!               -----------------------------------------------------
                call self % LookInOtherPartitions (i)
!
!               Disable the nodes if the point is not found
!               -------------------------------------------
                IF ( .not. self % active (i) ) then
                    IF ( .not. firstOutDomain) then
                        IF ( MPI_Process % isRoot ) then
                            firstOutDomain = .TRUE.
                        END IF
                    END IF
                END IF
!
!               Get the Lagrange interpolants
!               -----------------------------               
                if ( (MPI_Process % rank .ne. self % rank (i)).OR.( .not. self % active (i) ) ) then
                    self % lxi (:,i)= 0.0_RP
                    self % leta (:,i)= 0.0_RP
                    self % lzeta (:,i)= 0.0_RP
                ELSE
                    associate(e => mesh % elements(self % eID(i)))
                    associate( spAxi   => NodalStorage(e % Nxyz(1)), &
                                spAeta  => NodalStorage(e % Nxyz(2)), &
                                spAzeta => NodalStorage(e % Nxyz(3)) )
                    self % lxi  (:,i)  = spAxi % lj(self % xi(1,i))
                    self % leta  (:,i) = spAeta % lj(self % xi(2,i))
                    self % lzeta  (:,i)= spAzeta % lj(self % xi(3,i))
                    end associate
                    end associate
                END IF
                
            END DO
            
         end if
!
!        Check Variables, Create Files, and Write Header Files
!        -----------------------------------------------------
         do i=1,self % nVariables
!
!           Prepare the file in which the Plane is exported
!           -----------------------------------------------
            write( self % fileName (i) , '(A,A,A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % planeName) ,"_", trim(self % variable(i)) &
                                    , "_plane_" , self % ID,".sampling"  
!
!           Check the variable
!           ------------------
            call tolower(self % variable (i))

            select case ( trim(self % variable (i)) )
#ifdef NAVIERSTOKES
            case ("density")
            case ("pressure")
            case ("ptotal")
            case ("velocity")
            case ("viscosity")
            case ("u")
            case ("v")
            case ("w")
            case ("mach")
            case ("k")
            case ("omegax")
            case ("omegay")
            case ("omegaz")
            case ("q1")
            case ("q2")
            case ("q3")
            case ("q4")
            case ("q5")
            case default
                if ( MPI_Process % isRoot ) then 
                   print*, 'Plane variable "',trim(self % variable(i)),'" not implemented.'
                   print*, "Options available are:"
                   print*, "   * density"
                   print*, "   * pressure"
                   print*, "   * velocity"
                   print*, "   * viscosity"
                   print*, "   * u"
                   print*, "   * v"
                   print*, "   * w"
                   print*, "   * Mach"
                   print*, "   * K"
                   print*, "   * q1"
                   print*, "   * q2"
                   print*, "   * q3"
                   print*, "   * q4"
                   print*, "   * q5"
                end if 
#endif
#ifdef INCNS
            case default
               print*, "Planes are not implemented for the incompressible NSE"
#endif
#ifdef MULTIPHASE
            case ("density")
            case ("pressure")
			case ("static-pressure")
            case ("concentration")
            case ("velocity")
            case ("u")
            case ("v")
            case ("w")
            case ("omegax")
            case ("omegay")
            case ("omegaz")
            case ("q1")
            case ("q2")
            case ("q3")
            case ("q4")
            case ("q5")
            case default
                if ( MPI_Process % isRoot ) then 
                   print*, 'Plane variable "',trim(self % variable(i)),'" not implemented.'
                   print*, "Options available are:"
                   print*, "   * density"
                   print*, "   * pressure"
				   print*, "   * static-pressure"
                   print*, "   * velocity"
                   print*, "   * concentration"
                   print*, "   * u"
                   print*, "   * v"
                   print*, "   * w"
                   print*, "   * omegax"
				   print*, "   * omegay"
				   print*, "   * omegaz"
                   print*, "   * q1"
                   print*, "   * q2"
                   print*, "   * q3"
                   print*, "   * q4"
                   print*, "   * q5"
                end if 
#endif
		end select

        if ( .not. MPI_Process % isRoot ) return


!
!        Create file
!        -----------
         if (FirstCall) then
            if (MPI_Process % isRoot ) then 
                open ( newunit = fID , file = trim(self % fileName (i)) , action = "write" , access = "stream" , status = "replace", position='append' )
!
!           Write the file headers
!           ----------------------
                write( fID) self % ID
                write( fID) self % N(1)
                write( fID) self % N(2)
                write( fID) self % nNodes
                write( fID ) self % interval
#if defined(NAVIERSTOKES) && (!(INCNS))
                write( fID ) refValues % rho
                write( fID ) refValues % V
                write( fID ) refValues % p
                write( fID ) refValues % T
                write( fID ) refValues % mu
                write( fID ) refValues % AoATheta
                write( fID ) refValues % AoAPhi
#elif defined(MULTIPHASE)
                write( fID ) refValues % rho
                write( fID ) refValues % V
                write( fID ) refValues % p
                write( fID ) thermodynamics % rho(2)
                write( fID ) thermodynamics % mu(1)
                write( fID ) thermodynamics % mu(2)
                write( fID ) refValues % g0
#endif   
                write( fID ) transpose(self % x)
                close ( fID )
            end if
         end if
   
        end do 
!
!           File Format 
!           -----------------------------------------------
            write( fileFormat , '(A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % planeName) ,"_'variable'_plane_", self % ID, ".sampling"
!
!        Write Information 
!        -----------------------------------------------         
         write(STD_OUT,'(/)')
         call SubSection_Header("Plane Samplings")
            write(STD_OUT,'(30X,A,A27,I4)') "->" , "Plane ID: " , self % ID
            write(STD_OUT,'(30X,A,A27,A128)') "->" , "Variables: " , variables
			if (inputType.eq.0) then
			    write(STD_OUT,'(30X,A,A27,A36)') "->" , "Input Type: " , "Specified 3 points location"
				write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x1 [m]: ","[", &
													   point_1(1,1), ", ", point_1(2,1), ", ", point_1(3,1), "]"
				write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x2 [m]: ","[", &
													   point_2(1,1), ", ", point_2(2,1), ", ", point_2(3,1), "]"
				write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x3 [m]: ","[", &
													   point_3(1,1), ", ", point_3(2,1), ", ", point_3(3,1), "]"
	        else
			    write(STD_OUT,'(30X,A,A27,A27)') "->" , "Input Type: " , "Specified Input File"
				write(STD_OUT,'(30X,A,A27,A27)') "->" , "Input File: " , self % fileInput
				write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "First Point [m]: ","[", &
													   self % x(1,1), ", ", self % x(2,1), ", ", self % x(3,1), "]"
				write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Last Point [m]: ","[", &
													   self % x(1,self % nNodes), ", ", self % x(2,self % nNodes), ", ", self % x(3,self % nNodes), "]"
			end if 
            write(STD_OUT,'(30X,A,A27,A,I4,A,I4,A)') "->" , "Nodes [N12,N23]: ","[", &
                                                   self % N(1), ", ", self % N(2), "]"
            write(STD_OUT,'(30X,A,A27,I4)') "->" , "Samplings Interval: ", self % interval
            write(STD_OUT,'(30X,A,A27,A128)') "->" , "Filename: ", fileFormat
            if (firstOutDomain) then    
                write(STD_OUT,'(30X,A,A27,A61)') "->" ,"Note: ","Node/Nodes are located out of domain. A NaN will be assigned"
            end if
    
      end subroutine Plane_Initialization
!
!     Update Lagrange Interpolants After pAdaptation
!     -----------------------------------------------------
      subroutine Plane_UpdateLagrangeInterp(self, mesh)
         use Physics
         use MPI_Process_Info
         use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN
         implicit none
         class(PlaneSampling_t)           :: self
         class(HexMesh)                   :: mesh
	  
!
!        ---------------
!        Local variables
!        ---------------
!
         integer                          :: i
!
!           Find node location in the element ID
!           ------------------------------------                
            DO i=1,self % nNodes
!
!               Get the Lagrange interpolants
!               -----------------------------               
                if ( (MPI_Process % rank .ne. self % rank (i)).OR.( .not. self % active (i) ) ) then
                    self % lxi (:,i)= 0.0_RP
                    self % leta (:,i)= 0.0_RP
                    self % lzeta (:,i)= 0.0_RP
                ELSE
                    associate(e => mesh % elements(self % eID(i)))
                    associate( spAxi   => NodalStorage(e % Nxyz(1)), &
                                spAeta  => NodalStorage(e % Nxyz(2)), &
                                spAzeta => NodalStorage(e % Nxyz(3)) )
                    self % lxi  (:,i)  = spAxi % lj(self % xi(1,i))
                    self % leta  (:,i) = spAeta % lj(self % xi(2,i))
                    self % lzeta  (:,i)= spAzeta % lj(self % xi(3,i))
                    end associate
                    end associate
                END IF
                
            END DO
	  end subroutine Plane_UpdateLagrangeInterp

      subroutine Plane_Update(self, mesh, bufferPosition, t)
         use Physics
         use MPI_Process_Info
         use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN
         implicit none
         class(PlaneSampling_t)           :: self
         class(HexMesh)                   :: mesh
         integer                          :: bufferPosition
         real(kind=RP)                    :: t
         
!
!        ---------------
!        Local variables
!        ---------------
!
         integer        :: i, j, k, l, m, ierr
         real(kind=RP)  :: value, kappa
		 real(kind=RP)  :: c, sqrtRho
         real(kind=RP)  :: NaN
         real(kind=RP)  :: Sym, Asym
         real(kind=RP)  :: invRhoBase, invRhoF
         real(kind=RP), allocatable   :: var(:,:,:)   
                              
         if (self % intervalCount .EQ. 0 ) then
            self % bufferLine = self % bufferLine + 1

            DO m=1,self % nVariables
            
              self % values(1, bufferPosition, m) = t
              
              DO l=1, self % nNodes 
!
!               Assign NaN if the node located outside of domain
!               ------------------------------------------------              
                IF ( .not. self % active (l) ) THEN  
                    self % values(l, bufferPosition,m) = IEEE_VALUE(NaN, IEEE_QUIET_NAN)
                ELSE
                    if ( MPI_Process % rank .eq. self % rank (l) ) then
!
!                   Update the Node
!                   ----------------
                    associate( e => mesh % elements(self % eID(l)) )
					allocate (var(0:e % Nxyz(1),0:e % Nxyz(2),0:e % Nxyz(3)  ))	 
                    associate( Q => e % storage % Q,  S => e % storage )
   
#ifdef NAVIERSTOKES
            select case (trim(self % variable(m)))
            case("density")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHO,i,j,k)
               end do            ; end do             ; end do
            case("pressure")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Pressure(Q(:,i,j,k))
               end do            ; end do             ; end do
               
            case("viscosity")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  call get_laminar_mu_kappa(Q(:,i,j,k),var(i,j,k),kappa)
               end do            ; end do             ; end do
               
            case("ptotal")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k))     ! Vabs**2
                  var(i,j,k) =  var(i,j,k) / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * var(i,j,k)) ) ! Mach ^2
                  var(i,j,k) = Pressure(Q(:,i,j,k))*(1.0_RP+0.5_RP*(thermodynamics % gamma-1.0_RP)*var(i,j,k))**( thermodynamics % gamma/(thermodynamics % gamma-1.0_RP))
               end do            ; end do             ; end do
   
            case("velocity")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                  var(i,j,k) = sqrt(POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k)
               end do         ; end do         ; end do
            
            case("omegax")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_y(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k)) &
                                - (1/Q(IRHO,i,j,k) * S % U_z(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k))
               end do            ; end do             ; end do
               
            case("omegay")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_z(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k)) &
                                - (1/Q(IRHO,i,j,k) * S % U_x(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k))
               end do            ; end do             ; end do
               
            case("omegaz")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_x(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k)) &
                                - (1/Q(IRHO,i,j,k) * S % U_y(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k))
               end do            ; end do             ; end do
   
            case("u")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOU,i,j,k) / Q(IRHO,i,j,k)
               end do            ; end do             ; end do
   
            case("v")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOV,i,j,k) / Q(IRHO,i,j,k)
               end do            ; end do             ; end do
   
            case("w")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOW,i,j,k) / Q(IRHO,i,j,k)
               end do            ; end do             ; end do
   
            case("mach")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                  var(i,j,k) = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k))     ! Vabs**2
                  var(i,j,k) = sqrt( var(i,j,k) / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * var(i,j,k)) ) )
               end do         ; end do         ; end do
      
            case("k")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                  var(i,j,k) = 0.5_RP * (POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k)
               end do         ; end do         ; end do
               
            case("q1")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHO,i,j,k)
               end do            ; end do             ; end do
               
            case("q2")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOU,i,j,k)
               end do            ; end do             ; end do
               
            case("q3")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOV,i,j,k)
               end do            ; end do             ; end do
               
            case("q4")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOW,i,j,k)
               end do            ; end do             ; end do
               
            case("q5")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IRHOE,i,j,k)
               end do            ; end do             ; end do
               
            end select
#endif
#ifdef MULTIPHASE
            select case (trim(self % variable(m)))
            case("static-pressure")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                  var(i,j,k) = Q(IMP,i,j,k) + Q(IMC,i,j,k)*e % storage % mu(1,i,j,k) - 12.0_RP*multiphase%sigma*multiphase%invEps*(POW2(Q(IMC,i,j,k)*(1.0_RP-Q(IMC,i,j,k)))) &
                               - 0.25_RP*3.0_RP*multiphase % sigma * multiphase % eps * (POW2(e % storage % c_x(1,i,j,k))+POW2(e % storage % c_y(1,i,j,k))+POW2(e % storage % c_z(1,i,j,k)))
               end do         ; end do         ; end do
			   
            case("concentration")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
               end do            ; end do             ; end do
			   
            case("density")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  var(i,j,k) = c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)
               end do            ; end do             ; end do
			   
            case("pressure")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(5,i,j,k) 
               end do            ; end do             ; end do
   
            case("velocity")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = sqrt(POW2(Q(2,i,j,k)) + POW2(Q(3,i,j,k)) + POW2(Q(4,i,j,k)))/sqrtRho
               end do         ; end do         ; end do
            
            case("omegax")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = (1/sqrtRho * S % U_y(IMSQRHOW,i,j,k) - Q(IMSQRHOW,i,j,k)/(sqrtRho**2) * S % U_y(IMC,i,j,k)) &
                                - (1/sqrtRho * S % U_z(IMSQRHOV,i,j,k) - Q(IMSQRHOV,i,j,k)/(sqrtRho**2) * S % U_z(IMC,i,j,k))
               end do            ; end do             ; end do
               
            case("omegay")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = (1/sqrtRho * S % U_z(IMSQRHOU,i,j,k) - Q(IMSQRHOU,i,j,k)/(sqrtRho**2) * S % U_z(IMC,i,j,k)) &
                                - (1/sqrtRho * S % U_x(IMSQRHOW,i,j,k) - Q(IMSQRHOW,i,j,k)/(sqrtRho**2) * S % U_x(IMC,i,j,k))
               end do            ; end do             ; end do
               
            case("omegaz")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = (1/sqrtRho * S % U_x(IMSQRHOV,i,j,k) - Q(IMSQRHOV,i,j,k)/(sqrtRho**2) * S % U_x(IMC,i,j,k)) &
                                - (1/sqrtRho * S % U_y(IMSQRHOU,i,j,k) - Q(IMSQRHOU,i,j,k)/(sqrtRho**2) * S % U_y(IMC,i,j,k))
               end do            ; end do             ; end do
   
            case("u")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = Q(IMSQRHOU,i,j,k) / sqrtRho
               end do            ; end do             ; end do
   
            case("v")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = Q(IMSQRHOV,i,j,k) / sqrtRho
               end do            ; end do             ; end do
   
            case("w")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
			      c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP)
                  sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2))
                  var(i,j,k) = Q(IMSQRHOW,i,j,k) / sqrtRho
               end do            ; end do             ; end do
   
            case("q1")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IMC,i,j,k)
               end do            ; end do             ; end do
               
            case("q2")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IMSQRHOU,i,j,k)
               end do            ; end do             ; end do
               
            case("q3")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IMSQRHOV,i,j,k)
               end do            ; end do             ; end do
               
            case("q4")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IMSQRHOW,i,j,k)
               end do            ; end do             ; end do
               
            case("q5")
               do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1) 
                  var(i,j,k) = Q(IMP,i,j,k)
               end do            ; end do             ; end do
			   
            end select
#endif 
   
                    value = 0.0_RP
                    do k = 0, e % Nxyz(3)    ; do j = 0, e % Nxyz(2)  ; do i = 0, e % Nxyz(1)
                        value = value + var(i,j,k) * self % lxi(i,l) * self % leta(j,l) * self % lzeta(k,l)
                    end do               ; end do             ; end do
   
                    self % values(l+1, bufferPosition, m) = value
   
                    end associate
					deallocate(var)
                    end associate
#ifdef _HAS_MPI_            
                    if ( MPI_Process % doMPIAction ) then

!                   Share the result with the rest of the processes
!                   -----------------------------------------------         
                        call mpi_bcast(value, 1, MPI_DOUBLE, self % rank (l), MPI_COMM_WORLD, ierr)
                    end if
#endif
                    else

!                   Receive the result from the rank that contains the Plane
!                   --------------------------------------------------------
#ifdef _HAS_MPI_
                        if ( MPI_Process % doMPIAction ) then
                            call mpi_bcast(self % values(l+1, bufferPosition, m), 1, MPI_DOUBLE, self % rank (l), MPI_COMM_WORLD, ierr)
                        end if
#endif
                    end if
                END IF  
              END DO
            END DO
         end if
         self % intervalCount = self % intervalCount + 1
        
         
         if (self % intervalCount .EQ. self % interval) then
            self % intervalCount = 0 
         end if 

      end subroutine Plane_Update

      subroutine Plane_WriteToFile ( self, no_of_lines)
!
!        *************************************************************
!              This subroutine writes the buffer to the file.
!        *************************************************************
!
         implicit none  
         class(PlaneSampling_t)  :: self
         integer                 :: no_of_lines
!
!        ---------------
!        Local variables
!        ---------------
!
         integer                    :: i, j, ierr
         integer                    :: fID
        
         if ( MPI_Process % isRoot ) then
         DO j=1, self % nVariables
            open( newunit = fID , file = trim ( self % fileName (j) ) , action = "write" , access = "stream" , status = "old", position='append' )
            do i = 1 , no_of_lines
                write( fID ) self % values(:,i,j)
            end do
            close ( fID )
         END DO
         end if
         if ( no_of_lines .ne. 0 ) self % values(:,1,:) = self % values(:,no_of_lines,:)
         
         self % bufferLine = 0
      
      end subroutine Plane_WriteToFile 

      subroutine Plane_LookInOtherPartitions(self, i)
         use MPI_Process_Info
         implicit none
         class(PlaneSampling_t)     :: self
         integer ,intent(in)        :: i
         integer                    :: allActives(MPI_Process % nProcs)
         integer                    :: active, ierr

         if ( MPI_Process % doMPIAction ) then
#ifdef _HAS_MPI_
!
!           Cast the logicals onto integers
!           -------------------------------
            if ( self % active (i) ) then
               active = 1
            else
               active = 0
            end if
            allActives=0
!
!           Gather all data from all processes
!           ----------------------------------
            call mpi_allgather(active, 1, MPI_INT, allActives, 1, MPI_INT, MPI_COMM_WORLD, ierr)
!
!           Check if any of them found the Plane
!           ------------------------------------
            if ( any(allActives .eq. 1) ) then
!
!              Assign the domain of the partition that contains the Plane
!              ----------------------------------------------------------
               self % active (i) = .true.
               self % rank (i) = maxloc(allActives, dim = 1) - 1

            else
!
!              Disable the Plane
!              -----------------
               self % active (i) = .false.
               self % rank (i)  = -1
               self % eID (i) = 1

            end if
#endif
         else
!
!           Without MPI select the rank 0 as default
!           ----------------------------------------
            self % rank (i)= 0

         end if

      end subroutine Plane_LookInOtherPartitions
      
      elemental subroutine Plane_Destruct (self)
         implicit none
         class(PlaneSampling_t), intent(inout) :: self
         
         safedeallocate (self % values)
         safedeallocate (self % lxi)
         safedeallocate (self % leta)
         safedeallocate (self % lzeta)
         safedeallocate (self % active)
         safedeallocate (self % rank)
         safedeallocate (self % eID)
         safedeallocate (self % x)
		 safedeallocate (self % xi)
         safedeallocate (self % fileName)
         safedeallocate (self % variable)

      end subroutine Plane_Destruct
      
      elemental subroutine Plane_Assign (to, from)
         implicit none
         class(PlaneSampling_t), intent(inout) :: to
         type(PlaneSampling_t) , intent(in) :: from
         
         safedeallocate ( to % active )
         allocate ( to % active ( size(from % active) ) )
         to % active = from % active
         
         safedeallocate ( to % rank )
         allocate ( to % rank ( size(from % rank) ) )
         to % rank = from % rank
         
         to % ID = from %  ID
         to % nVariables      = from % nVariables
         to % interval        = from % interval
         to % bufferSize      = from % bufferSize
         to % bufferLine      = from % bufferLine
         to % intervalCount   = from % intervalCount
         to % N               = from % N
         to % nNodes          = from % nNodes
         to % disturbanceData = from % disturbanceData
         
         safedeallocate ( to % eID )
         allocate ( to % eID ( size(from % eID) ) )
         to % eID = from % eID      

         safedeallocate ( to % lxi )
         allocate ( to % lxi ( size(from % lxi,1), size(from % lxi,2) ) )
         to % lxi = from % lxi
         
         safedeallocate ( to % leta )
         allocate ( to % leta ( size(from % leta,1), size(from % leta,2) ) )
         to % leta = from % leta
         
         safedeallocate ( to % lzeta )
         allocate ( to % lzeta ( size(from % lzeta,1), size(from % lzeta,2) ) )
         to % lzeta = from % lzeta
         
         safedeallocate ( to % values )
         allocate ( to % values( size(from % values,1),size(from % values,2),size(from % values,3) ) )
         to % values = from % values
         
         safedeallocate ( to % x )
         allocate ( to % x( size(from % x,1),size(from % x,2) ) )
         to % x = from % x       
         
         safedeallocate ( to % xi )
         allocate ( to % xi( size(from % xi,1),size(from % xi,2) ) )
         to % xi = from % xi    

         safedeallocate ( to % fileName )
         allocate ( to % fileName ( size(from % fileName) ) )
         to % fileName = from % fileName
         
         to % planeName = from % planeName

         safedeallocate ( to % variable )
         allocate ( to % variable ( size(from % variable) ) )
         to % variable = from % variable         

         
         
      end subroutine Plane_Assign
      
end module PlaneSampling

