!
!//////////////////////////////////////////////////////
!
!This class represents the observer for the fW-H acoustic analogy, including the relations with several observers

#include "Includes.h"
Module  FWHObseverClass  !

   use SMConstants
   use FaceClass
   use Physics
   use PhysicsStorage
   use NodalStorageClass
   use FWHDefinitions, only: OB_BUFFER_SIZE, OB_BUFFER_SIZE_DEFAULT, STR_LEN_OBSERVER
   use ZoneClass
   use HexMeshClass
   use MPI_Process_Info
#ifdef _HAS_MPI_
   use mpi
#endif
   Implicit None

!
!  *****************************
!  Observer source pair class definition
!  class for the coupling of each pair of observer and source(face)
!  mainly acoustic geometrical relations and face link
!  *****************************
   type ObserverSourcePairClass
       real(kind=RP), dimension(:,:,:), allocatable        :: rVect
       real(kind=RP), dimension(:,:),   allocatable        :: r
       real(kind=RP), dimension(:,:),   allocatable        :: re
       real(kind=RP), dimension(:,:,:), allocatable        :: reUnitVect 
       real(kind=RP), dimension(:,:),   allocatable        :: reStar
       real(kind=RP), dimension(:,:,:), allocatable        :: reStarUnitVect 
       real(kind=RP)                                       :: tDelay
       integer                                             :: faceIDinMesh    ! ID of the source (face) at the Mesh array (linked list)
       integer                                             :: elementSide
       real(kind=RP)                                       :: normalCorrection
       real(kind=RP), dimension(:,:),   allocatable        :: Pacc ! temporal solution of acoustic pressure for each pair
       real(kind=RP), dimension(:),     allocatable        :: tInterp ! time array for interpolation

       contains

           procedure :: construct       => ObserverSourcePairConstruct
           procedure :: destruct        => ObserverSourcePairDestruct
           procedure :: allocPacc       => ObserverSourcePairAllocSolution
           procedure :: interpolateSolF => ObserverSourcePairInterpolateSolFirst
           procedure :: newUpdate       => ObserverSourcePairNewUpdate
           procedure :: interpolateSolS => ObserverSourcePairInterpolateSolSecond
           procedure :: updateOneStep   => ObserverSourcePairUpdateOneStep
           procedure :: FWHSurfaceIntegral

   end type ObserverSourcePairClass
!
!  *****************************
!  General observer class definition
!   (similar to a monitor, mostly surface monitor in many behaviours)
!  *****************************
!  
   type ObserverClass

       integer                                                         :: ID
       real(kind=RP), dimension(NDIM)                                  :: x        ! position of the observer at global coordinates
       integer                                                         :: numberOfFaces
       class(ObserverSourcePairClass), dimension(:), allocatable       :: sourcePair
       real(kind=RP), dimension(:,:), allocatable                      :: Pac      ! acoustic pressure, two components and the total (sum)
       real(kind=RP)                                                   :: tDelay
       real(kind=RP)                                                   :: tDelayMax
       logical                                                         :: active
       character(len=STR_LEN_OBSERVER)                                 :: observerName
       character(len=STR_LEN_OBSERVER)                                 :: fileName

       contains

           procedure :: construct      => ObserverConstruct
           procedure :: destruct       => ObserverDestruct
           procedure :: update         => ObserverUpdate
           procedure :: writeToFile    => ObserverWriteToFile
           procedure :: updateTdelay   => ObserverUpdateTdelay
           procedure :: interpolateSol => ObserverInterpolateSol
           procedure :: sumIntegrals   => ObserverSumIntegrals
           procedure :: updateOneStep  => ObserverUpdateOneStep

   end type ObserverClass

   contains

!/////////////////////////////////////////////////////////////////////////
!           OBSERVER CLASS PROCEDURES --------------------------
!/////////////////////////////////////////////////////////////////////////

   Subroutine ObserverConstruct(self, sourceZone, mesh, ID, solution_file, FirstCall, interpolate, totalNumberOfFaces, elementSide)

!        *****************************************************************************
!              This subroutine initializes the observer similar to a monitor. The following
!           data is obtained from the case file:
!              -> Name: The observer name (10 characters maximum)
!              -> x: The observer position
!        *****************************************************************************

       use ParamfileRegions
       use FileReadingUtilities, only: getRealArrayFromString
       implicit none

       class(ObserverClass)                                 :: self
       class(Zone_t), intent(in)                            :: sourceZone
       class(HexMesh), intent(in)                           :: mesh
       integer, intent(in)                                  :: ID, totalNumberOfFaces
       character(len=*), intent(in)                         :: solution_file
       logical, intent(in)                                  :: FirstCall, interpolate
       integer, dimension(:), intent(in)                    :: elementSide

       ! local variables
       character(len=STR_LEN_OBSERVER)  :: in_label
       character(len=STR_LEN_OBSERVER)  :: fileName
       character(len=STR_LEN_OBSERVER)  :: paramFile
       character(len=STR_LEN_OBSERVER)  :: coordinates
       integer                          :: fID
       integer                          :: MeshFaceID, zoneFaceID
       ! integer                          :: elementSide
!
!      Get observer ID
!      --------------
       self % ID = ID
!
!      Search for the parameters in the case file
!      ------------------------------------------
       write(in_label , '(A,I0)') "#define acoustic observer " , self % ID

       call get_command_argument(1, paramFile)
       call readValueInRegion(trim ( paramFile), "name",   self % observerName, in_label, "# end" ) 
       call readValueInRegion(trim(paramFile), "position", coordinates        , in_label, "# end" )

!      Get the coordinates
!      -------------------
       self % x = getRealArrayFromString(coordinates)

!     Enable the observer
!     ------------------
      self % active = .true.
      allocate ( self % Pac(OB_BUFFER_SIZE,3) )

!     Get source information
!     ------------------
      self % numberOfFaces = sourceZone % no_of_faces

!     Construct each pair observer-source
!     ------------------
      allocate( self % sourcePair(self % numberOfFaces) )
!     Loop the zone to get faces
      do zoneFaceID = 1, self % numberOfFaces
!         Face global ID
          MeshFaceID = sourceZone % faces(zoneFaceID)
          call self % sourcePair(zoneFaceID) % construct(self % x, mesh % faces(MeshFaceID), MeshFaceID, FirstCall, elementSide(zoneFaceID))
      end do  

!     Allocate variables for interpolation
!     -------------------------------------------------
      if (interpolate) then
          do zoneFaceID = 1, self % numberOfFaces
              call self % sourcePair(zoneFaceID) % allocPacc(OB_BUFFER_SIZE)
          end do
      end if 

!     Set the average time delay of the observer
!     -------------------------------------------------
      call self % updateTdelay(totalNumberOfFaces)

!     Prepare the file in which the observer is exported
!     -------------------------------------------------
      write( self % fileName , '(A,A,A,A)') trim(solution_file) , "." , trim(self % observerName) , ".observer"
!
!     Create file
!     -----------
      if (FirstCall) then
         open ( newunit = fID , file = trim(self % fileName) , status = "unknown" , action = "write" ) 

!        Write the file headers
!        ----------------------
         write( fID , '(A20,A  )') "Observer name:      ", trim(self % observerName)
         write( fID , '(A20,ES24.10)') "x coordinate: ", self % x(1)
         write( fID , '(A20,ES24.10)') "y coordinate: ", self % x(2)
         write( fID , '(A20,ES24.10)') "z coordinate: ", self % x(3)

         write( fID , * )
         write( fID , '(A10,5(2X,A24))' ) "Iteration" , "Time" , "Observer_Time", "P'T", "P'L", "P'"

         close ( fID )
      end if

   End Subroutine ObserverConstruct

   Subroutine ObserverUpdate(self, mesh, isSolid, BufferPosition, interpolate)

!     *******************************************************************
!        This subroutine updates the observer acoustic pressure computing it from
!        the mesh storage. It is stored in the "bufferPosition" position of the 
!        buffer.
!     *******************************************************************
!
use VariableConversion, only: Pressure, PressureDot
      implicit none
      class (ObserverClass)                                :: self
      class (HexMesh), intent(in)                          :: mesh
      integer,intent(in), optional                         :: bufferPosition
      logical, intent(in)                                  :: isSolid, interpolate

      ! local variables
      real(kind=RP)                                        :: Pt, Pl  ! pressure of each pair
      real(kind=RP), dimension(3)                          :: localPacc, Pacc   ! temporal variable to store the sum of the pressure
      real(kind=RP), dimension(3)                          :: mInterp ! slope of interpolation
      real(kind=RP)                                        :: valx, valy, valz
      integer                                              :: zoneFaceID, meshFaceID,  ierr
      integer                                             :: storePosition

!     Initialization
!     --------------            
      if (present(bufferPosition)) self % Pac(bufferPosition,:) = 0.0_RP
      Pacc = 0.0_RP
      valx = 0.0_RP
      valy = 0.0_RP
      valz = 0.0_RP

!     Loop the pairs (equivalent to loop the zone) and get the values
!     ---------------------------------------
      interp_cond: if (interpolate) then
!        For this case only save the values of the solution of each pair, at the corresponding position
!        ---------------------------------------
!$omp parallel private(meshFaceID,storePosition,localPacc) shared(mesh,isSolid,interpolate,Pacc,NodalStorage,&
!$omp&                                                     self,bufferPosition)
!$omp do private(meshFaceID,storePosition,localPacc) schedule(runtime)
         do zoneFaceID = 1, self % numberOfFaces
!            Compute the integral
!            --------------------
             meshFaceID = self % sourcePair(zoneFaceID) % faceIDinMesh
             localPacc = self % sourcePair(zoneFaceID) % FWHSurfaceIntegral( mesh % faces(meshFaceID), isSolid )

             !save solution at bufferPosition or last position
             if (present(bufferPosition)) then
                 storePosition = bufferPosition
             else
                 storePosition = size(self % sourcePair(zoneFaceID) % Pacc, dim=1)
             end if
             self % sourcePair(zoneFaceID) % Pacc(storePosition,:) = localPacc
         end do  
!$omp end do
!$omp end parallel
     else interp_cond
!        For this case get the whole solution of the observer, adding all the pairs without saving
!        ---------------------------------------
!$omp parallel private(meshFaceID, localPacc) shared(mesh,isSolid,interpolate,Pacc,NodalStorage,&
!$omp&                                        self,valx,valy,valz)
!$omp do private(meshFaceID,localPacc) reduction(+:valx,valy,valz) schedule(runtime)
          do zoneFaceID = 1, self % numberOfFaces
!            Compute the integral
!            --------------------
             meshFaceID = self % sourcePair(zoneFaceID) % faceIDinMesh

             localPacc = self % sourcePair(zoneFaceID) % FWHSurfaceIntegral( mesh % faces(meshFaceID), isSolid )

             ! sum without interpolate: suppose little change of each tDelay
             valx = valx + localPacc(1)
             valy = valy + localPacc(2)
             valz = valz + localPacc(3)
         end do  
!$omp end do
!$omp end parallel

          Pacc = (/valx, valy, valz/)

#ifdef _HAS_MPI_
      localPacc = Pacc
      call mpi_allreduce(localPacc, Pacc, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
#endif

          self % Pac(bufferPosition,:) = Pacc
     end if interp_cond

   End Subroutine ObserverUpdate

   Subroutine ObserverUpdateTdelay(self, totalNumberOfFaces)

!     *******************************************************************
!        This subroutine updates the observer time delay. For static surfaces it 
!        doesn't need to be updated at every iteration.
!        Minimum time is used, for interpolation procedures
!     *******************************************************************
!
      use MPI_Process_Info
      implicit none
      class   (ObserverClass)                              :: self
      integer, intent(in)                                  :: totalNumberOfFaces
      
      ! local variables
      integer                                              :: i, ierr
      real(kind=RP)                                        :: t, tmax
      real(kind=RP), dimension(:), allocatable             :: alltDelay
      real(kind=RP), dimension(self % numberOfFaces)       :: tDelayArray
      integer, dimension(MPI_Process % nProcs)             :: no_of_faces_p, displs

      do i =1, self % numberOfFaces
        tDelayArray(i) = self % sourcePair(i) % tDelay
      end do

      if ( (MPI_Process % doMPIAction) ) then
#ifdef _HAS_MPI_
          call mpi_gather(self % numberOfFaces,1,MPI_INT,no_of_faces_p,1,MPI_INT,0,MPI_COMM_WORLD,ierr)

          if (MPI_Process % isRoot) then
              displs=0
              do i = 2, MPI_Process % nProcs 
                  displs(i) = displs(i-1) + no_of_faces_p(i-1)
              end do
          end if

          allocate(alltDelay(totalNumberOfFaces))
          call mpi_gatherv(tDelayArray, self % numberOfFaces, MPI_DOUBLE, &
                                     alltDelay, no_of_faces_p, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD, ierr)

          if (MPI_Process % isRoot) then
              t = minval(alltDelay)
              tmax = maxval(alltDelay)
          end if

          call mpi_Bcast(t, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD, ierr)
          call mpi_Bcast(tmax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD, ierr)
#endif
      else
          t = minval(tDelayArray)
          tmax = maxval(tDelayArray)
      end if

      self % tDelay = t
      self % tDelayMax = tmax

   End Subroutine ObserverUpdateTdelay

   Subroutine ObserverWriteToFile(self, iter, tsource, no_of_lines)
!
!     *************************************************************
!           This subroutine writes the buffer to the file.
!     *************************************************************
!
      implicit none  
      class(ObserverClass)                                 :: self
      integer, dimension(:)                                :: iter
      real(kind=RP), dimension(:)                          :: tsource
      integer                                              :: no_of_lines
!
!     ---------------
!     Local variables
!     ---------------
!
      integer                    :: i
      integer                    :: fID

      if ( MPI_Process % isRoot ) then

         open( newunit = fID , file = trim ( self % fileName ) , action = "write" , access = "append" , status = "old" )
      
         do i = 1 , no_of_lines
            write( fID , '(I10,5(2X,ES24.16))' ) iter(i) , tsource(i), tsource(i) + self % tDelay, self % Pac(i,:)
         end do
      
         close ( fID )
      end if
      
      if ( no_of_lines .ne. 0 ) self % Pac(1,:) = self % Pac(no_of_lines,:)
      
   End Subroutine ObserverWriteToFile

    Subroutine ObserverUpdateOneStep(self, mesh, BufferPosition, isSolid, tsource)

      implicit none

      class (ObserverClass)                                :: self
      class (HexMesh), intent(in)                          :: mesh
      integer,intent(in)                                   :: bufferPosition
      logical, intent(in)                                  :: isSolid
      real(kind=RP), intent(in)                            :: tsource

      ! local variables
      real(kind=RP)                                        :: tobserver
      integer                                              :: i
      integer, dimension(self%numberOfFaces)               :: nDiscard

      ! store the solution of each pair at the last position, by not giving the bufferPosition
      call self % update(mesh, isSolid, interpolate=.TRUE.)

      ! interpolate the solution of each pair at first position
      ! and save the time of each pair at its last position
      tobserver = tsource + self % tDelay
!$omp parallel shared(self)
!$omp do schedule(runtime)
      do i = 1, self % numberOfFaces
            if (self % sourcePair(i) % tDelay .eq. self % tDelay) cycle
            call self % sourcePair(i) % interpolateSolS(tobserver, tsource)
      end do 
!$omp end do
!$omp end parallel

      ! sum all the pair solution and save it at bufferPosition of the observer sol
      nDiscard = 0
      call self % sumIntegrals(nDiscard, 1, bufferPosition, bufferPosition)

      ! update the solution of each pair and its times for next iteration
!$omp parallel shared(self)
!$omp do schedule(runtime)
      do i = 1, self % numberOfFaces
            call self % sourcePair(i) % updateOneStep()
      end do 
!$omp end do
!$omp end parallel

    End Subroutine ObserverUpdateOneStep

   !interpolate the solution of all the pairs to get it at the mean observer time
   Subroutine ObserverInterpolateSol(self, tsource, no_of_lines)

      implicit none  

      class(ObserverClass)                                 :: self
      real(kind=RP), dimension(:), intent(in)              :: tsource
      integer, intent(in)                                  :: no_of_lines

      !local variables
      real(kind=RP), dimension(:), allocatable             :: tobserver
      integer                                              :: i, n, k, m
      integer, dimension(self % numberOfFaces)             :: nDiscard
      logical                                              :: sameDelay

      allocate(tobserver(no_of_lines))
      tobserver = tsource(1:no_of_lines) + self % tdelay

      ! get max tobserver that can be interpolated
      do k =1, no_of_lines
          if ( tobserver(k) .ge. (self % tDelayMax + tsource(1)) ) exit
      end do
      n = no_of_lines - k + 1 ! k is the min tobserver index

      safedeallocate(tobserver)
      allocate(tobserver(n))
      tobserver(1:n) = tsource(k:no_of_lines) + self % tdelay

!$omp parallel shared(self, nDiscard, n, no_of_lines, tobserver, tsource,k)
!$omp do schedule(runtime)
      do i = 1, self % numberOfFaces
          ! call interp of each pair that are not the minimum
          ! if (almostequal(self % sourcepair(i) % tdelay, self % tdelay)) then
          if (self % sourcepair(i) % tdelay .eq. self % tdelay) then
              nDiscard(i) = k-1
          else
              call self % sourcePair(i) % interpolateSolF(n, no_of_lines, tobserver, tsource(1:no_of_lines), nDiscard(i))
          end if
      end do
!$omp end do
!$omp end parallel

      ! set to 0 the first part of the solution, which cannot be interpolated
      ! in this case Pacc is written from 1:no_of_lines, which have a value of 0 at first positions, not need to change obs write proc
      self % pac(1:k-1,:) = 0.0_RP

      ! sum all values from k to no_of_lines
      call self % sumIntegrals(nDiscard, n, k, no_of_lines)

      ! update all the solution of the pair to save the future ones
      do i = 1, self % numberOfFaces
          ! sameDelay = almostequal(self % sourcepair(i) % tdelay, self % tdelay)
          sameDelay = self % sourcepair(i) % tdelay .eq. self % tdelay
          call self % sourcePair(i) % newUpdate(n, nDiscard(i), no_of_lines, tsource(1:no_of_lines), sameDelay)
      end do

   End Subroutine ObserverInterpolateSol

   ! sum all the interpolated solution of all pairs and save it at the observer solution
   Subroutine ObserverSumIntegrals(self, nDiscard, N, startIndex, no_of_lines)

      implicit none  

      class(observerclass)                                 :: self
      integer, intent(in)                                  :: no_of_lines, startIndex, N
      integer, dimension(self % numberOfFaces), intent(in) :: nDiscard

      ! local variables
      real(kind=RP), dimension(:,:), allocatable           :: localPacc, Pacc   ! temporal variable to store the sum of the pressure
      real(kind=RP), dimension(:), allocatable             :: valx, valy, valz
      integer                                              :: i, ierr

!     Initialization
!     --------------            
      ! 1:N must be equal to startIndex:no_of_lines
      allocate(Pacc(N,3), localPacc(N,3), valx(N), valy(N), valz(N))
      Pacc = 0.0_RP
      valx = 0.0_RP
      valy = 0.0_RP
      valz = 0.0_RP

!$omp parallel private(localPacc) shared(Pacc,nDiscard,N,self,valx,valy,valz)
!$omp do private(localPacc) reduction(+:valx,valy,valz) schedule(runtime)
      do i = 1, self % numberOfFaces
!        Get the array of interpolated values of each pair
         localPacc(1:N,:) = self % sourcePair(i) % Pacc(nDiscard(i)+1:nDiscard(i)+N,:)

         ! sum interpolated
         valx = valx + localPacc(:,1)
         valy = valy + localPacc(:,2)
         valz = valz + localPacc(:,3)

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

      Pacc(:,1) = valx(:)
      Pacc(:,2) = valy(:)
      Pacc(:,3) = valz(:)

#ifdef _HAS_MPI_
      localPacc = Pacc
      call mpi_allreduce(localPacc, Pacc, 3*N, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
#endif

      self % Pac(startIndex:no_of_lines,:) = Pacc(1:N,:)

   End Subroutine ObserverSumIntegrals

   Subroutine ObserverDestruct(self)

        implicit none
        class(ObserverClass), intent(inout)               :: self

       ! local variables
       integer                                            :: i
        safedeallocate (self % Pac)
        do i = 1, self % numberOfFaces
            call self % sourcePair(i) % destruct
        end do
        safedeallocate (self % sourcePair)

   End Subroutine ObserverDestruct

!/////////////////////////////////////////////////////////////////////////
!           OBSERVER SOURCE PAIR CLASS PROCEDURES --------------------------
!/////////////////////////////////////////////////////////////////////////

   Subroutine  ObserverSourcePairConstruct(self, x, f, fID, FirstCall, elementSide)

       ! use fluiddata
       use FWHDefinitions, only: rho0, P0, c0, U0, M0, fwGamma2
       implicit none

       class(ObserverSourcePairClass)                      :: self
       real(kind=RP), dimension(NDIM), intent(in)          :: x       ! observer position
       type(face), intent(in)                              :: f    ! source
       integer, intent(in)                                 :: fID, elementSide
       logical, intent(in)                                 :: FirstCall

       ! local variables
       integer                                             :: Nx,Ny
       integer                                             :: i, j
       real(kind=RP)                                       :: fwGammaInv

       self % faceIDinMesh = fID

       Nx = f % Nf(1)
       Ny = f % Nf(2)

       self % elementSide = elementSide

   select case (elementSide)
   case (1)
       self % normalCorrection = 1.0_RP
   case (2)
       self % normalCorrection = -1.0_RP
   end select

       allocate( self % r(0:Nx,0:Ny), self % re(0:Nx,0:Ny), self % reStar(0:Nx,0:Ny) )
       allocate( self % rVect(NDIM,0:Nx,0:Ny), self % reUnitVect(NDIM,0:Nx,0:Ny) ,self % reStarUnitVect(NDIM,0:Nx,0:Ny) )

       fwGammaInv = 1.0_RP / sqrt(fwGamma2)
       ! source position, for each node of the face
       associate (y => f % geom % x)
           do j= 0, Ny; do i = 0,Nx
               ! store geometrical acoustic relations for each node
               self % rVect(:,i,j) = x(:) - y(:,i,j)
               self % r(i,j) = norm2(self % rVect(:,i,j))
               self % reStar(i,j) = fwGammaInv*sqrt( self%r(i,j)**2 + fwGamma2*( dot_product(M0, self%rVect(:,i,j)) )**2 )
               self % reStarUnitVect(:,i,j) = ( self%rVect(:,i,j) + fwGamma2*dot_product(M0, self%rVect(:,i,j))*M0(:) ) / &
                                            (fwGamma2*self%reStar(i,j))
               self % re(i,j) = fwGamma2*( self%reStar(i,j) - dot_product(M0, self%rVect(:,i,j)) )
               self % reUnitVect(:,i,j) = fwGamma2*( self%reStarUnitVect(:,i,j) - M0(:) )
               self % tDelay = (sum(self%re))/real(size(self%re),RP) / c0
           end do; end do
       end associate

   End Subroutine ObserverSourcePairConstruct 

  elemental Subroutine ObserverSourcePairDestruct(self)

      Class(ObserverSourcePairClass), intent(inout)       :: self
 
      safedeallocate(self % rVect)
      safedeallocate(self % r)
      safedeallocate(self % re)
      safedeallocate(self % reUnitVect)
      safedeallocate(self % reStar)
      safedeallocate(self % reStarUnitVect)

  End Subroutine ObserverSourcePairDestruct 

  ! allocate time history solution for a posterior interpolation

  Subroutine ObserverSourcePairAllocSolution(self, buffer_size)

       class(ObserverSourcePairClass)                      :: self
       integer, intent(in)                                 :: buffer_size
       
       allocate(self % Pacc(buffer_size,3))

  End Subroutine ObserverSourcePairAllocSolution

  Subroutine ObserverSourcePairInterpolateSolFirst(self, N, M, tobserver, tsource, nd)

       implicit none

       class(ObserverSourcePairClass)                      :: self
       integer, intent(in)                                 :: N, M
       real(kind=RP), dimension(N), intent(in)             :: tobserver
       real(kind=RP), dimension(:), intent(in)             :: tsource
       integer, intent(out)                                :: nd

       ! local variables
       real(kind=RP), dimension(N,3)                       :: PaccInterp    ! solution of the pair interpolated
       real(kind=RP), dimension(M)                         :: tPair         ! time array of the panel
       integer                                             :: i, j, ii

       ! get the times of the pair for interpolation
       tPair = tsource + self % tDelay

       j = 1
       nd = 0
       do i = 2, M
           if (j .gt. N) exit
           ii = i - 1
           if (tPair(i) .lt. tobserver(1)) then
               nd = nd +1
               cycle
           end if 
           PaccInterp(j,:) = linearInterpolation(tobserver(j), tPair(ii), self%Pacc(ii,:), tPair(i), self%Pacc(i,:), 3)
           j = j + 1
       end do 

           !update solution of the pair with the interpolated one
           self % Pacc(nd+1:nd+N,:) = PaccInterp

  End Subroutine ObserverSourcePairInterpolateSolFirst

  Subroutine ObserverSourcePairNewUpdate(self, N, NDiscard, M, tsource, sameDelay)

       implicit none

       class(ObserverSourcePairClass)                      :: self
       integer, intent(in)                                 :: N, NDiscard, M
       ! real(kind=RP), dimension(N), intent(in)             :: tobserver
       real(kind=RP), dimension(:), intent(in)             :: tsource
       logical, intent(in)                                 :: sameDelay

       !local variables
       real(kind=RP), dimension(M)                         :: tPair         ! time array of the panel
       real(kind=RP), dimension(:,:), allocatable          :: PaccFuture    ! solution of the pair of future (have not been interpolainterpolated) values
       integer                                             :: Nfuture

       tPair = tsource + self % tDelay
       if (sameDelay) then
           !size = 1 for last value + 1(empty for next iter)
           Nfuture = 2
       else
           !size = M - interpolated values +1 + 1(empty for next iter)
           Nfuture = M - N -NDiscard + 1
       end if 

       ! save old results and kept last position empty
       allocate(PaccFuture(Nfuture-1,3))
       PaccFuture(:,:) = self % Pacc(M-Nfuture+2:M,:)

       safedeallocate(self % Pacc)
       allocate(self % Pacc(1:Nfuture,3))

       self % Pacc(1:Nfuture-1,:) = PaccFuture(:,:)

       if (.not. sameDelay) then
           allocate(self % tInterp(1:Nfuture))
           ! save old results and kept last position empty
           self % tInterp(1:Nfuture-1) = tPair(NDiscard+N+1:M)
       end if 

  End Subroutine ObserverSourcePairNewUpdate

  ! save the solution and times from position 2 to last as the first position, letting free the last one
  Subroutine ObserverSourcePairUpdateOneStep(self)

      implicit none
      class(ObserverSourcePairClass)                      :: self

      !local variables
      integer                                             :: M

      M = size(self%Pacc, dim=1)
      if (allocated(self%tInterp)) self % tInterp(1:M-1) = self % tInterp(2:M)
      self % Pacc(1:M-1,:) = self % Pacc(2:M,:)

  End Subroutine ObserverSourcePairUpdateOneStep

  Subroutine ObserverSourcePairInterpolateSolSecond(self, tobserver, tsource)

       implicit none

       class(ObserverSourcePairClass)                      :: self
       real(kind=RP),  intent(in)                          :: tobserver
       real(kind=RP),  intent(in)                          :: tsource

       ! local variables
       real(kind=RP), dimension(2)                         :: tPair         ! time array of the panel
       real(kind=RP), dimension(3)                         :: PaccInterp    ! solution of the pair interpolated

       ! save last time
       self % tInterp(size(self % tInterp)) = tsource + self % tDelay
       ! use 2 first values to interpolate
       tPair = self % tInterp(1:2)

       PaccInterp(:) = linearInterpolation(tobserver, tPair(1), self % Pacc(1,:), tPair(2), self % Pacc(2,:), 3)

       !update the first value of the solution of the pair with the interpolated one
       self % Pacc(1,:) = PaccInterp(:)

    End Subroutine ObserverSourcePairInterpolateSolSecond

   ! calculate the surface integrals of the FW-H analogy for stationary surfaces (permeable or impermeable) with a general flow
   ! direction of the medium
   ! the integrals are for a single face (pane in FWH terminology) for a single observer
!         TODO: check if is more efficient to store FWHvariables for each face instead of calculating it always
!               for many observers, its being recomputed as many as observers

   ! Function FWHSurfaceIntegral(self, f, isSolid, interpolate, bufferPosition) result(Pacc)
   Function FWHSurfaceIntegral(self, f, isSolid) result(Pacc)

       use FWHDefinitions, only: rho0, P0, c0, U0, M0
       use VariableConversion, only: Pressure, PressureDot
       use fluiddata, only: dimensionless
       implicit none

       class(ObserverSourcePairClass)                      :: self
       class(Face), intent(in)                             :: f
       logical, intent(in)                                 :: isSolid
       real(kind=RP),dimension(3)                          :: Pacc  ! acoustic pressure values
       ! logical, intent(in)                                 :: isSolid, interpolate
       ! integer, intent(in), optional                       :: bufferPosition

       ! local variables
       integer                                             :: i, j  ! face indexes
       real(kind=RP), dimension(NDIM)                      :: Qi,QiDot, n
       real(kind=RP), dimension(NDIM,NDIM)                 :: Lij, LijDot
       type(NodalStorage_t), pointer                       :: spAxi, spAeta
       real(kind=RP)                                       :: Pt, Pl
       real(kind=RP)                                       :: LR, MR, UmMr,LdotR, LM
       ! integer                                             :: storePosition

       ! Initialization
       Pt = 0.0_RP
       Pl = 0.0_RP
       spAxi  => NodalStorage(f % Nf(1))
       spAeta => NodalStorage(f % Nf(2))


       associate( Q => f % storage(self % elementSide) % Q )
           associate( Qdot => f % storage(self % elementSide) % Qdot )

    !           **********************************
    !           Computes the surface integral
    !              I = \int vec{f}·vec{n} * vec{g}·vec{r} dS
    !           **********************************
    !
                do j = 0, f % Nf(2) ;    do i = 0, f % Nf(1)

                   n = f % geom % normal(:,i,j) * self % normalCorrection
                   call calculateFWHVariables(Q(:,i,j), Qdot(:,i,j), isSolid, Qi, QiDot, Lij, LijDot)

                   LR = dot_product(matmul(Lij, n(:)), self%reUnitVect(:,i,j))
                   MR = dot_product(M0(:), self % reUnitVect(:,i,j))
                   UmMr = 1 - MR
                   LdotR = dot_product(matmul(LijDot, n(:)), self%reUnitVect(:,i,j))
                   LM = dot_product(matmul(Lij, n(:)), M0(:))

                   ! loading term integrals
                   Pl = Pl +  dot_product(matmul(LijDot,n(:)),self%reUnitVect(:,i,j)) / (self%reStar(i,j) * c0) * &
                             spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   Pl = Pl +  dot_product(matmul(Lij,n(:)),self%reStarUnitVect(:,i,j)) / (self%reStar(i,j)**2) * &
                             spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + LdotR / ( c0 * self % re(i,j) * (UmMr**2) ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + (LR - LM) / ( (self % re(i,j) * UmMr)**2 ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + (LR * (MR - (dimensionless % Mach**2))) / ( (self % re(i,j)**2) * (UmMr**3) ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)

                   ! thickness term integrals, only for permeable surfaces
                   if (.not. isSolid) then
                       Pt = Pt + (1 - dot_product(M0(:),self%reUnitVect(:,i,j))) * dot_product(QiDot(:),n(:)) / (self%reStar(i,j)) * &
                                 spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       Pt = Pt -  dot_product(U0(:),self%reStarUnitVect(:,i,j)) * dot_product(Qi(:),n(:)) / (self%reStar(i,j)**2) * &
                                 spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       ! Pt = Pt + dot_product(QiDot(:),n(:)) / (self%reStar(i,j) * (UmMr**2)) * &
                       !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       ! Pt = Pt + (dot_product(Qi(:), n(:)) * c0 * (MR - (dimensionless % Mach**2))) / ( (self % re(i,j)**2) * (UmMr**3) ) * &
                       !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)

                   end if  
                end do          ;    end do
           end associate
       end associate

       Pt = Pt / (4.0_RP * PI)
       Pl = Pl / (4.0_RP * PI)

      ! get total acoustic pressure as the sum of the two components (the quadrapol terms are being ignored)
       Pacc = (/Pt, Pl, Pt+Pl/)

   End Function FWHSurfaceIntegral 

!/////////////////////////////////////////////////////////////////////////
!           ZONE PROCEDURES --------------------------
!/////////////////////////////////////////////////////////////////////////

   Subroutine SourceProlongSolution(source_zone, mesh, eSides)

!     *******************************************************************
!        This subroutine prolong the solution from the mesh storage to the faces (source).
!         TODO: use openmp (commented)
!         TODO: use mpi (see surface integral)
!     *******************************************************************
!
      use ElementClass
      implicit none
      class (Zone_t), intent(in)                           :: source_zone
      class (HexMesh), intent(inout), target               :: mesh
      integer, dimension(:), intent(in)                    :: eSides

      ! local variables
      integer                                              :: zoneFaceID, meshFaceID, eID
      integer, dimension(6)                                :: meshFaceIDs
      class(Element), pointer                              :: elements(:)

!     *************************
!     Perform the interpolation
!     *************************
!
      elements => mesh % elements
!$omp parallel private(meshFaceID,eID,meshFaceIDs) shared(elements,mesh,NodalStorage)
!!$omp&                                        t)
!$omp single

!        Loop the zone to get faces and elements
!        ---------------------------------------
      do zoneFaceID = 1, source_zone % no_of_faces
          meshFaceID = source_zone % faces(zoneFaceID)

         eID = mesh % faces(meshFaceID) % elementIDs(eSides(zoneFaceID))
         meshFaceIDs = mesh % elements(eID) % faceIDs

!$omp task depend(inout:elements(eID))
         call elements(eID) % ProlongSolutionToFaces(NCONS,&
                                                            mesh % faces(meshFaceIDs(1)),&
                                                            mesh % faces(meshFaceIDs(2)),&
                                                            mesh % faces(meshFaceIDs(3)),&
                                                            mesh % faces(meshFaceIDs(4)),&
                                                            mesh % faces(meshFaceIDs(5)),&
                                                            mesh % faces(meshFaceIDs(6)),&
                                                             computeQdot = .TRUE.)

         ! if ( computeGradients ) then
         !    call elements(eID) % ProlongGradientsToFaces(NGRAD, mesh % faces(meshFaceIDs(1)),&
         !                                     mesh % faces(meshFaceIDs(2)),&
         !                                     mesh % faces(meshFaceIDs(3)),&
         !                                     mesh % faces(meshFaceIDs(4)),&
         !                                     mesh % faces(meshFaceIDs(5)),&
         !                                     mesh % faces(meshFaceIDs(6)) )
         ! end if
!$omp end task
      end do
!$omp end single
!$omp end parallel

   End Subroutine SourceProlongSolution
!
!/////////////////////////////////////////////////////////////////////////
!           AUXILIARY PROCEDURES --------------------------
!/////////////////////////////////////////////////////////////////////////

! get the interpolated value of an array
! ------------
  Function linearInterpolation(x, x1, y1, x2, y2, N) result(y)

      integer, intent(in)                           :: N
      real(kind=RP), intent(in)                     :: x1, x2, x
      real(kind=RP), dimension(N), intent(in)       :: y1, y2
      real(kind=RP), dimension(N)                   :: y

      y = y1 + (y2-y1)/(x2-x1) * (x - x1)

  End Function linearInterpolation

   Subroutine calculateFWHVariables(Q, Qdot, isSolid, Qi, QiDot, Lij, LijDot)

       use VariableConversion, only: Pressure, PressureDot
       use FWHDefinitions,     only: rho0, P0, c0, U0, M0
       use Utilities, only: AlmostEqual
       implicit none

       real(kind=RP), dimension(NCONS), intent(in)         :: Q        ! horses variables array
       real(kind=RP), dimension(NCONS), intent(in)         :: Qdot     ! horses time derivatives array
       logical, intent(in)                                 :: isSolid
       real(kind=RP), dimension(NDIM), intent(out)         :: Qi       ! fwh Qi array, related with the acoustic pressure thickness
       real(kind=RP), dimension(NDIM), intent(out)         :: Qidot
       real(kind=RP), dimension(NDIM,NDIM), intent(out)    :: Lij      ! fwh Lij tensor: related with the acoustic pressure loading
       real(kind=RP), dimension(NDIM,NDIM), intent(out)    :: LijDot

       !local variables
       real(kind=RP)                                       :: P, pDot
       real(kind=RP), dimension(NDIM,NDIM)                 :: Pij      ! fwh perturbation stress tensor
       ! real(kind=RP), dimension(NDIM:NDIM)                 :: tau
       integer                                             :: i, j, ii, jj

       P = Pressure(Q)
       pDot = PressureDot(Q,Qdot)

       Pij = 0.0_RP
       LijDot = 0.0_RP
       do i=1,NDIM
           Pij(i,i) = P - P0
           !pressure derivative of LijDot
           LijDot(i,i) = pDot
       end do

       !TODO use the stress tensor and the time derivative for Lij and LijDot respectively
       ! call getStressTensor(Q, U_x, U_y, U_z, tau)
       ! Pij = Pij - tau
       ! LijDot = LijDot - tauDot

       ! set values for solid (impermeable) surface
       Qi(:) = -rho0*U0(:)
       Qidot = 0.0_RP
       Lij = Pij
       ! Lij = 0.0_RP

       !calculate terms for permeable surface
       if (.not. isSolid) then
           Qi(:) = Qi(:) + Q(2:4)
           ! convert to complete velocity instead of perturbation velocity
           QiDot(:) = QiDot(:) + Qdot(2:4)

           do j = 1, NDIM
               jj = j + 1
               do i = 1, NDIM
                   ! one index is added since rhoV1 = Q(2), rhoV2 = Q(3) ...
                   ii = i + 1
                   Lij(i,j) = Lij(i,j) + (Q(ii) - Q(1)*U0(i))*(Q(jj)/Q(1))
                   LijDot(i,j) = LijDot(i,j) + ( Qdot(ii) - Q(ii)/Q(1)*Qdot(1) )/Q(1) * Q(jj) + &
                                               (Q(ii)/Q(1) - U0(i)) * Qdot(jj)
               end do  
           end do  
       end if

   End Subroutine calculateFWHVariables

End Module  FWHObseverClass 