PLAG_RFLU_ModFindCells.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Suite of routines for particle tracking on Eulerian grid.
!
! Description: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: PLAG_RFLU_ModFindCells.F90,v 1.16 2008/12/06 08:44:35 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

MODULE plag_rflu_modfindcells
  
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moddatatypes
  USE modglobal,     ONLY: t_global
  USE modgrid,       ONLY: t_grid
  USE modpartlag,    ONLY: t_plag, t_plag_input
  USE modbndpatch,   ONLY: t_patch
  USE moddatastruct, ONLY: t_region
  USE modborder,     ONLY: t_border
  USE moderror
  USE modmpi
  USE plag_modparameters
    
  USE modtools, ONLY: floatequal
  
  USE rflu_modcellfaceedgeinfo, ONLY: rflu_getglobalcellkind, & 
                                      rflu_getfacekind  
  USE rflu_modgeometrytools
  USE rflu_modincelltest, ONLY: rflu_ict_testincell, & 
                                rflu_ict_testincellfancy, & 
                                rflu_ict_testincelllohner
  USE rflu_modoctree
  USE rflu_modrelatedpatches, ONLY: rflu_relp_transformvector

  USE plag_modinterfaces, ONLY: plag_reflectparticledata 
  USE plag_modsurfstats, ONLY: plag_gathersurfstats

  USE rflu_modmpi, ONLY: rflu_mpi_recreatebufferipclsend  

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: plag_rflu_computedisttot,    &
            plag_rflu_findcellsbrute,    &
            plag_rflu_findcellsbrutemod, &
            plag_rflu_findcellslohner,   & 
            plag_rflu_findcellsoct,      &
            plag_rflu_findcellsoctmod,   &
            plag_rflu_findcellstrajfast, &
            plag_rflu_findcellstrajsafe
        
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
   
  CHARACTER(CHRLEN) :: RCSIdentString = & 
    '$RCSfile: PLAG_RFLU_ModFindCells.F90,v $ $Revision: 1.16 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS











! ******************************************************************************
!
! Purpose: Compute total distance travelled.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE plag_rflu_computedisttot(pRegion)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  CHARACTER(CHRLEN) :: errorstring
  INTEGER ::ipcl
  REAL(RFREAL) :: xlocnew,xlocold,xtraj,ylocnew,ylocold,ytraj,zlocnew, &
                  zlocold,ztraj
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_ComputeDistTot',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls
 
! ==============================================================================  
!   Compute distance travelled and trajectory
! ==============================================================================  
  
    xlocnew = pplag%cv(cv_plag_xpos,ipcl)
    ylocnew = pplag%cv(cv_plag_ypos,ipcl)
    zlocnew = pplag%cv(cv_plag_zpos,ipcl)  

    xlocold = pplag%cvOld(cv_plag_xpos,ipcl)
    ylocold = pplag%cvOld(cv_plag_ypos,ipcl)
    zlocold = pplag%cvOld(cv_plag_zpos,ipcl)
        
    xtraj = xlocnew - xlocold
    ytraj = ylocnew - ylocold      
    ztraj = zlocnew - zlocold

    pplag%arv(arv_plag_distot,ipcl)  = sqrt( xtraj*xtraj + ytraj*ytraj &
                                           + ztraj*ztraj )
          
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_computedisttot









! ******************************************************************************
!
! Purpose: Determine cells which contain particles using brute force approach.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsbrute(pRegion)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: foundflag
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: errorflag,icg,ipcl
  REAL(RFREAL) :: xloc,yloc,zloc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsBrute',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls  

! ==============================================================================  
!   Get particle position
! ==============================================================================  

    xloc = pplag%cv(cv_plag_xpos,ipcl)
    yloc = pplag%cv(cv_plag_ypos,ipcl)
    zloc = pplag%cv(cv_plag_zpos,ipcl)

! ==============================================================================  
!   Loop over actual cells
! ==============================================================================  
  
    foundflag = .false.
  
    cellloop: DO icg = 1,pgrid%nCells         
      IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,icg) .EQV. .true. ) THEN
        pplag%aiv(aiv_plag_regini,ipcl) = pregion%iRegionGlobal
        pplag%aiv(aiv_plag_icells,ipcl) = icg

        foundflag = .true.

        EXIT cellloop
      END IF ! RFLU_ICT_TestInCell        
    END DO cellloop
      
! ==============================================================================  
!   Check whether particle was found
! ==============================================================================  
            
    IF ( foundflag .EQV. .false. ) THEN
      WRITE(errorstring,'(I6)') ipcl 
      CALL errorstop(global,err_plag_pcl_not_found,__line__,trim(errorstring))
    END IF ! foundFlag
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsbrute







! ******************************************************************************
!
! Purpose: Kernel for brute force search.
!
! Description: None.
!
! Input: 
!   pRegion             Pointer to region
!   xLoc,yLoc,zLoc      Location
!
! Output:
!   icgOut              Cell containing location
!
! Notes: 
!   1. icgOut is set to CRAZY_VALUE_INT if no cell contained location.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsbrutekernel(pRegion,xLoc,yLoc,zLoc,icgOut)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER, INTENT(OUT) :: icgout
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  INTEGER :: errorflag,icg
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsBruteKernel',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid

! ******************************************************************************
! Loop over cells
! ******************************************************************************

  icgout = crazy_value_int

  cellloop: DO icg = 1,pgrid%nCells         
    IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,icg) .EQV. .true. ) THEN
      icgout = icg

      EXIT cellloop
    END IF ! RFLU_ICT_TestInCell        
  END DO cellloop
      
! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsbrutekernel







! ******************************************************************************
!
! Purpose: Determine cells which contain particles using modified brute 
!  force approach.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: 
!   1. This routine is suitable to determine the cell containing a particle 
!      after the latters position has been updated. 
!   2. This routine cannot be used in parallel runs and for particles which 
!      leave the domain or bounce from solid walls.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsbrutemod(pRegion)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: errorflag,icg,ipcl
  REAL(RFREAL) :: xloc,yloc,zloc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsBruteMod',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls  

! ==============================================================================  
!   Get particle cell and position
! ==============================================================================  

    icg  = pplag%aivOld(aiv_plag_icells,ipcl) ! NOTE OLD cell index

    xloc = pplag%cv(cv_plag_xpos,ipcl) ! NOTE already NEW particle position
    yloc = pplag%cv(cv_plag_ypos,ipcl)
    zloc = pplag%cv(cv_plag_zpos,ipcl)

! ==============================================================================  
!   Check if particle still in same cell 
! ==============================================================================  

    IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,icg) .EQV. .false. ) THEN    

! ------------------------------------------------------------------------------  
!     Find cell containing particle
! ------------------------------------------------------------------------------  
  
      CALL plag_rflu_findcellsbrutekernel(pregion,xloc,yloc,zloc,icg)
  
      IF ( icg /= crazy_value_int ) THEN 
        pplag%aiv(aiv_plag_icells,ipcl) = icg
      ELSE 

! TEMPORARY
        WRITE(*,*) 'timeCurrent,iPcl,xLocOld,yLocOld,zLocOld,xLoc,yLoc,zLoc = ',&
        global%currentTime,ipcl,pplag%cvOld(cv_plag_xpos:cv_plag_zpos,ipcl),&
        xloc,yloc,zloc 
! END TEMPORARY

        WRITE(errorstring,'(I6)') ipcl 
        CALL errorstop(global,err_plag_pcl_not_found,__line__,trim(errorstring))     
      END IF ! icg
            
! ==============================================================================  
!   Particle still in same cell, copy cell index
! ==============================================================================        
      
    ELSE   
      pplag%aiv(aiv_plag_icells,ipcl) = icg
    END IF ! RFLU_ICT_TestInCell
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsbrutemod







! ******************************************************************************
!
! Purpose: Determine cells which contain particles using Lohners approach.
!
! Description: Follow particle path by passing particle to face which failed
!   in-cell test by largest amount.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: 
!   1. See R. Lohner and J. Ambrosiano, A Vectorized Particle Tracer for 
!      Unstructured Grids, J. Comp. Phys., Vol. 91, pp. 22-31, 1990.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellslohner(pRegion)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: testincell
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: c1,c1k,c2,c2k,errorflag,icg,icgout,ifg,iloc,ipcl,loopcounter
  REAL(RFREAL) :: xloc,yloc,zloc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsLohner',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls
    loopcounter = 0 ! Reset loop counter 

    xloc = pplag%cv(cv_plag_xpos,ipcl)
    yloc = pplag%cv(cv_plag_ypos,ipcl)
    zloc = pplag%cv(cv_plag_zpos,ipcl)  
                   
    icg = pplag%aivOld(aiv_plag_icells,ipcl)

! ==============================================================================
!   Loop until distance travelled along trajectory consumed 
! ==============================================================================        
        
    infloop: DO     
      loopcounter = loopcounter + 1  
        
! ------------------------------------------------------------------------------
!     Find appropriate intersection and associated face  
! ------------------------------------------------------------------------------
                    
      CALL rflu_ict_testincelllohner(pregion,xloc,yloc,zloc,icg,testincell, &
                                     iloc,ifg)                                                        
                                   
! ------------------------------------------------------------------------------
!     Check whether have remaining total distance 
! ------------------------------------------------------------------------------

! --- In-cell test successful --------------------------------------------------

      IF ( testincell .EQV. .true. ) THEN 
        pplag%aiv(aiv_plag_icells,ipcl) = icg

        EXIT infloop
        
! --- In-cell test failed ------------------------------------------------------
        
      ELSE 

! ----- Intersect interior face

        IF ( iloc == 0 ) THEN           
          c1 = pgrid%f2c(1,ifg)
          c2 = pgrid%f2c(2,ifg)
          
          c1k = rflu_getglobalcellkind(global,pgrid,c1)
          c2k = rflu_getglobalcellkind(global,pgrid,c2)

          SELECT CASE ( rflu_getfacekind(global,c1k,c2k) ) 
            CASE ( face_kind_aa ) ! Actual-actual face
              IF ( c1 == icg ) THEN 
                icg = c2
              ELSE 
                icg = c1
              END IF ! c1 
! TO DO
!            CASE ( FACE_KIND_AV ) ! Actual-virtual face
! Count particle intersections for each buffer           
! END TO DO
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! RFLU_GetFaceKind                    
                  
! ----- Intersect boundary face          
          
        ELSE ! Boundary face
          ppatch => pregion%patches(iloc)

          SELECT CASE ( ppatch%bcType )             
            CASE default 
              CALL errorstop(global,err_reached_default,__line__)           
          END SELECT ! pPatch%bcType
        END IF ! iLoc 
      END IF ! testInCell
      
! ------------------------------------------------------------------------------
!     Guard against infinite loop 
! ------------------------------------------------------------------------------  
  
      IF ( loopcounter >= limit_infinite_loop ) THEN 
        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! loopCounter
                    
    END DO infloop     
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellslohner







! ******************************************************************************
!
! Purpose: Determine cells which contain particles using Octree approach.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: 
!   1. It is not enough to simply get closest cell from Octree, because 
!      proximity to cell centroid does not guarantee that the specified 
!      location is actually inside the cell.
!   2. This routine is suitable to determine the cells which contain an 
!      initial distribution of particles, but should not be used to track 
!      moving particles, because it will be inefficient on account of it
!      not using knowledge of past position. It will also overwrite entry 
!      containing initial region in particle data structure. 
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsoct(pRegion)

  USE rflu_modgeometrytools, ONLY: rflu_testinboundbox

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: foundflag
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: ccsize,errorflag,icg,ipcl
  INTEGER, DIMENSION(:), ALLOCATABLE :: cc
  REAL(RFREAL) :: delfrac,xdel,xloc,xmax,xmin,ydel,yloc,ymax,ymin,zdel,zloc, & 
                  zmax,zmin
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsOct',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

  delfrac = 0.01_rfreal
  
  ccsize = min(100,pgrid%nCells) ! Must be larger than unity 

  ALLOCATE(cc(ccsize),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'cc')
  END IF ! global%error

! ******************************************************************************  
! Build bounding box
! ******************************************************************************

  xmin = minval(pgrid%xyz(xcoord,1:pgrid%nVert))
  xmax = maxval(pgrid%xyz(xcoord,1:pgrid%nVert))
  ymin = minval(pgrid%xyz(ycoord,1:pgrid%nVert))
  ymax = maxval(pgrid%xyz(ycoord,1:pgrid%nVert))
  zmin = minval(pgrid%xyz(zcoord,1:pgrid%nVert))
  zmax = maxval(pgrid%xyz(zcoord,1:pgrid%nVert))

  xdel = xmax - xmin 
  ydel = ymax - ymin
  zdel = zmax - zmin

  xmin = xmin - delfrac*xdel
  xmax = xmax + delfrac*xdel 
  ymin = ymin - delfrac*ydel
  ymax = ymax + delfrac*ydel 
  zmin = zmin - delfrac*zdel
  zmax = zmax + delfrac*zdel  

! ******************************************************************************  
! Build Octree with cell centroids: NOTE only for actual cells
! ******************************************************************************

  CALL rflu_createoctree(global,pgrid%nCells)

  CALL rflu_buildoctree(pgrid%cofg(xcoord,1:pgrid%nCells), & 
                        pgrid%cofg(ycoord,1:pgrid%nCells), & 
                        pgrid%cofg(zcoord,1:pgrid%nCells), & 
                        xmin,xmax,ymin,ymax,zmin,zmax)

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls

! ==============================================================================  
!   Test particle location against bounding box of partition
! ==============================================================================  
  
    xloc = pplag%cv(cv_plag_xpos,ipcl)
    yloc = pplag%cv(cv_plag_ypos,ipcl)
    zloc = pplag%cv(cv_plag_zpos,ipcl)
  
    IF ( rflu_testinboundbox(global,xloc,yloc,zloc,xmin,xmax,ymin,ymax, & 
                             zmin,zmax) .EQV. .true. ) THEN  
                               
! ------------------------------------------------------------------------------  
!     Query octree to get closest cells  
! ------------------------------------------------------------------------------ 
    
      CALL rflu_queryoctree(xloc,yloc,zloc,ccsize,cc)
        
! ------------------------------------------------------------------------------
!     Test cells obtained from Octree if they contain specified location
! ------------------------------------------------------------------------------        
           
      foundflag = .false.     
           
      cellloop: DO icg = 1,ccsize           
        IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,cc(icg)) .EQV. .true. ) THEN
          pplag%aiv(aiv_plag_regini,ipcl) = pregion%iRegionGlobal
          pplag%aiv(aiv_plag_icells,ipcl) = cc(icg)

          foundflag = .true.

          EXIT cellloop
        END IF ! RFLU_ICT_TestInCell        
      END DO cellloop
      
! ------------------------------------------------------------------------------
!     Check whether particle was found
! ------------------------------------------------------------------------------        
            
      IF ( foundflag .EQV. .false. ) THEN
        WRITE(errorstring,'(I6)') ipcl 
        CALL errorstop(global,err_plag_pcl_not_found,__line__, &
                       trim(errorstring))
      END IF ! foundFlag

! ==============================================================================  
!   Particle located outside bounding box of partition
! ==============================================================================  

    ELSE   
      WRITE(errorstring,'(I6)') ipcl 
      CALL errorstop(global,err_plag_pcl_not_found,__line__,trim(errorstring))      
    END IF ! RFLU_TestInBoundNox 
  END DO ! iPcl

! ******************************************************************************
! Destroy Octree and deallocate temporary memory
! ******************************************************************************

  CALL rflu_destroyoctree(global)

  DEALLOCATE(cc,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'cc')
  END IF ! global%error

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsoct








! ******************************************************************************
!
! Purpose: Kernel for Octree search
!
! Description: None.
!
! Input: 
!   pRegion             Pointer to region
!   xLoc,yLoc,zLoc      Location
!
! Output:
!   icgOut              Cell containing location
!
! Notes: 
!   1. It is not enough to simply get closest cell from Octree, because 
!      proximity to cell centroid does not guarantee that the specified 
!      location is actually inside the cell.
!   2. This routine is suitable to determine the cell containing a particle 
!      after the latters position has been updated. 
!   3. This routine cannot be used in parallel runs and for particles which 
!      leave the domain or bounce from solid walls.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsoctkernel(pRegion,xLoc,yLoc,zLoc,xMin,xMax,&
                                        ymin,ymax,zmin,zmax,icgout)

  USE rflu_modgeometrytools, ONLY: rflu_testinboundbox
  
  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER, INTENT(INOUT)   :: icgout
  REAL(RFREAL), INTENT(IN) :: xloc,xmax,xmin,yloc,ymax,ymin,zloc,zmax,zmin
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: foundflag
  INTEGER :: ccsize,errorflag,icgloop,ipcl
  INTEGER, DIMENSION(:), ALLOCATABLE :: cc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsOctKernel',&
  'PLAG_RFLU_ModFindCells.F90')
  
! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid

!  ccSize = MIN(100,pGrid%nCells) ! Must be larger than unity 

! TEMPORARY
  ccsize = min(50,pgrid%nCells) ! Must be larger than unity 
! END TEMPORARY

  ALLOCATE(cc(ccsize),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'cc')
  END IF ! global%error

! ******************************************************************************
! Test particle location against bounding box of partition
! ******************************************************************************

  IF ( rflu_testinboundbox(global,xloc,yloc,zloc,xmin,xmax,ymin,ymax, & 
                             zmin,zmax) .EQV. .true. ) THEN  

! ============================================================================== 
!  Query octree to get closest cells  
! ============================================================================== 
    
   CALL rflu_queryoctree(xloc,yloc,zloc,ccsize,cc)

! ============================================================================== 
!  Test cells obtained from Octree if they contain specified location
! ============================================================================== 
! ------------------------------------------------------------------------------ 
           
   foundflag = .false.     

   cellloop: DO icgloop = 1,ccsize           
      IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,cc(icgloop)) .EQV. .true. ) THEN
        icgout = cc(icgloop)

        foundflag = .true.

        EXIT cellloop
      END IF ! RFLU_ICT_TestInCell        
    END DO cellloop

! ============================================================================== 
!   Check whether particle was found 
! ============================================================================== 
            
    IF ( foundflag .EQV. .false. ) THEN
      icgout = crazy_value_int
    END IF ! foundFlag

! ******************************************************************************
! Particle located outside bounding box of partition
! ******************************************************************************

  ELSE   
    icgout = crazy_value_int    
  END IF ! RFLU_TestInBoundNox

! ******************************************************************************
! Deallocate temporary memory
! ******************************************************************************

  DEALLOCATE(cc,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'cc')
  END IF ! global%error
    
! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsoctkernel








! ******************************************************************************
!
! Purpose: Determine cells which contain particles using modified Octree 
!   approach.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!
! Output: None.
!
! Notes: 
!   1. It is not enough to simply get closest cell from Octree, because 
!      proximity to cell centroid does not guarantee that the specified 
!      location is actually inside the cell.
!   2. This routine is suitable to determine the cell containing a particle 
!      after the latters position has been updated. 
!   3. This routine cannot be used in parallel runs and for particles which 
!      leave the domain or bounce from solid walls.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellsoctmod(pRegion)

  USE rflu_modgeometrytools, ONLY: rflu_testinboundbox

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: foundflag
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: ccsize,errorflag,icg,ipcl
  INTEGER, DIMENSION(:), ALLOCATABLE :: cc
  REAL(RFREAL) :: delfrac,xdel,xloc,xmax,xmin,ydel,yloc,ymax,ymin,zdel,zloc, & 
                  zmax,zmin
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsOctMod',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

  delfrac = 0.01_rfreal
  
!  ccSize = MIN(100,pGrid%nCells) ! Must be larger than unity 

! TEMPORARY
  ccsize = min(50,pgrid%nCells) ! Must be larger than unity 
! END TEMPORARY

  ALLOCATE(cc(ccsize),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'cc')
  END IF ! global%error

! ******************************************************************************  
! Build bounding box
! ******************************************************************************

  xmin = minval(pgrid%xyz(xcoord,1:pgrid%nVert))
  xmax = maxval(pgrid%xyz(xcoord,1:pgrid%nVert))
  ymin = minval(pgrid%xyz(ycoord,1:pgrid%nVert))
  ymax = maxval(pgrid%xyz(ycoord,1:pgrid%nVert))
  zmin = minval(pgrid%xyz(zcoord,1:pgrid%nVert))
  zmax = maxval(pgrid%xyz(zcoord,1:pgrid%nVert))

  xdel = xmax - xmin 
  ydel = ymax - ymin
  zdel = zmax - zmin

  xmin = xmin - delfrac*xdel
  xmax = xmax + delfrac*xdel 
  ymin = ymin - delfrac*ydel
  ymax = ymax + delfrac*ydel 
  zmin = zmin - delfrac*zdel
  zmax = zmax + delfrac*zdel  

! ******************************************************************************  
! Build Octree with cell centroids: NOTE only for actual cells
! ******************************************************************************

  CALL rflu_createoctree(global,pgrid%nCells)

  CALL rflu_buildoctree(pgrid%cofg(xcoord,1:pgrid%nCells), & 
                        pgrid%cofg(ycoord,1:pgrid%nCells), & 
                        pgrid%cofg(zcoord,1:pgrid%nCells), & 
                        xmin,xmax,ymin,ymax,zmin,zmax)

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = 1,pplag%nPcls

! ==============================================================================  
!   Get particle cell and position
! ==============================================================================  

    icg  = pplag%aivOld(aiv_plag_icells,ipcl) ! NOTE OLD cell index
  
    xloc = pplag%cv(cv_plag_xpos,ipcl) ! NOTE already NEW particle position
    yloc = pplag%cv(cv_plag_ypos,ipcl)
    zloc = pplag%cv(cv_plag_zpos,ipcl)

! ==============================================================================  
!   Check if particle still in same cell 
! ==============================================================================  
    
    IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,icg) .EQV. .false. ) THEN    
 
! ------------------------------------------------------------------------------  
!     Test particle location against bounding box of partition
! ------------------------------------------------------------------------------  
  
      IF ( rflu_testinboundbox(global,xloc,yloc,zloc,xmin,xmax,ymin,ymax, & 
                               zmin,zmax) .EQV. .true. ) THEN  
                               
! ----- Query octree to get closest cells -------------------------------------- 
    
        CALL rflu_queryoctree(xloc,yloc,zloc,ccsize,cc)
        
! ----- Test cells obtained from Octree if they contain specified location -----
           
        foundflag = .false.     

        cellloop: DO icg = 1,ccsize           
          IF ( rflu_ict_testincell(pregion,xloc,yloc,zloc,cc(icg)) .EQV. .true. ) THEN
            pplag%aiv(aiv_plag_icells,ipcl) = cc(icg)

            foundflag = .true.

            EXIT cellloop
          END IF ! RFLU_ICT_TestInCell        
        END DO cellloop
      
! ----- Check whether particle was found ---------------------------------------
            
        IF ( foundflag .EQV. .false. ) THEN
          WRITE(errorstring,'(I6)') ipcl 
          CALL errorstop(global,err_plag_pcl_not_found,__line__, &
                         trim(errorstring))
        END IF ! foundFlag

! ------------------------------------------------------------------------------  
!     Particle located outside bounding box of partition
! ------------------------------------------------------------------------------  

      ELSE   

! TEMPORARY
        WRITE(*,*) 'timeCurrent,iPcl,xLocOld,yLocOld,zLocOld,xLoc,yLoc,zLoc = ',&
        global%currentTime,ipcl,pplag%cvOld(cv_plag_xpos:cv_plag_zpos,ipcl),&
        xloc,yloc,zloc
! END TEMPORARY

        WRITE(errorstring,'(I6)') ipcl 
        CALL errorstop(global,err_plag_pcl_not_found,__line__,trim(errorstring))      
      END IF ! RFLU_TestInBoundNox

! ==============================================================================  
!   Particle still in same cell, copy cell index 
! ==============================================================================  
    
    ELSE   
      pplag%aiv(aiv_plag_icells,ipcl) = icg
    END IF ! RFLU_ICT_TestInCell 
  END DO ! iPcl

! ******************************************************************************
! Destroy Octree and deallocate temporary memory
! ******************************************************************************

  CALL rflu_destroyoctree(global)

  DEALLOCATE(cc,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'cc')
  END IF ! global%error

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellsoctmod







! ******************************************************************************
!
! Purpose: Determine cells which contain particles by following trajectory 
!   using fast algorithm.
!
! Description: Follow particle path by intersecting particle trajectory with
!   faces and taking appropriate action determined by type of intersected face.
!
! Input: 
!  pRegion      Pointer to region
!  iPclBeg      Beginning particle index
!  iPclEnd      Ending particle index
!
! Output: None.
!
! Notes:
!   1. Need to properly handle particle deletion from outflow or motion
!      between domain during an intermediate RK-stage as its final 
!      position at (n+1) step might not be outside the region.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellstrajfast(pRegion,iPclBeg,iPclEnd)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER :: ipclbeg,ipclend
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: incellcheckflag
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: c1,c1k,c2,c2k,errorflag,iborder,icg,ifg,ifgout,ifl,iloc,ilocout, & 
             ipatch,ipcl,loopcounter
  REAL(RFREAL) :: dist,disttot,disttotcutoff,eps,fnx,fny,fnz,imagtraj,theta, &
                  xloc,xlocnew,xlocold,xtraj,yloc,ylocnew,ylocold,ytraj,zloc, &
                  zlocnew,zlocold,ztraj
  TYPE(t_border), POINTER :: pborder
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsTrajFast',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

  eps = epsilon(1.0_rfreal)
  disttotcutoff = 10*epsilon(1.0_rfreal) ! NOTE must be less than tolerICT

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = ipclbeg,ipclend
    loopcounter = 0 ! Reset loop counter 

! ==============================================================================  
!   Set distance travelled and compute trajectory. NOTE that cannot compute 
!   trajectory (unit) vector from distTot because distTot diminishes as the 
!   particle travels along the trajectory but xLocNew and xLocOld (as well as
!   y and z components) do not change.
! ==============================================================================  

    disttot = pplag%arv(arv_plag_distot,ipcl)

    IF ( disttot < disttotcutoff ) THEN  
      cycle 
    END IF ! distTot

    xlocnew = pplag%cv(cv_plag_xpos,ipcl)
    ylocnew = pplag%cv(cv_plag_ypos,ipcl)
    zlocnew = pplag%cv(cv_plag_zpos,ipcl)  

    xlocold = pplag%cvOld(cv_plag_xpos,ipcl)
    ylocold = pplag%cvOld(cv_plag_ypos,ipcl)
    zlocold = pplag%cvOld(cv_plag_zpos,ipcl)
        
    xtraj = xlocnew - xlocold
    ytraj = ylocnew - ylocold     
    ztraj = zlocnew - zlocold
             
    imagtraj = 1.0_rfreal/(sqrt(xtraj*xtraj + ytraj*ytraj + ztraj*ztraj) + eps)
    
    xtraj = imagtraj*xtraj
    ytraj = imagtraj*ytraj
    ztraj = imagtraj*ztraj    
    
! ------------------------------------------------------------------------------
!   Set variables for trajectory loop
! ------------------------------------------------------------------------------    
    
    xloc = xlocnew - disttot*xtraj
    yloc = ylocnew - disttot*ytraj
    zloc = zlocnew - disttot*ztraj

    icg = pplag%aivOld(aiv_plag_icells,ipcl)

! ==============================================================================
!   Loop until distance travelled along trajectory consumed 
! ==============================================================================        
        
    trajloop: DO     
      loopcounter = loopcounter + 1  
        
! ------------------------------------------------------------------------------
!     Find appropriate intersection and associated face  
! ------------------------------------------------------------------------------
                    
      CALL rflu_computelinecellxsectfast(pregion,xloc,yloc,zloc,xtraj,ytraj, &
                                         ztraj,icg,dist,iloc,ifg)

! ------------------------------------------------------------------------------
!     Update total distance travelled                             
! ------------------------------------------------------------------------------
                                   
      disttot = disttot - dist

! ------------------------------------------------------------------------------
!     Check whether have remaining total distance 
! ------------------------------------------------------------------------------

! --- No distance remaining: Found cell containing new location ----------------

      IF ( disttot <= 0.0_rfreal ) THEN 
        pplag%aiv(aiv_plag_icells,ipcl) = icg

        EXIT trajloop
        
! --- Distance remaining: Keep searching ---------------------------------------        
        
      ELSE 

! ----- Trajectory intersects interior face

        IF ( iloc == 0 ) THEN           
          c1 = pgrid%f2c(1,ifg)
          c2 = pgrid%f2c(2,ifg)
          
          c1k = rflu_getglobalcellkind(global,pgrid,c1)
          c2k = rflu_getglobalcellkind(global,pgrid,c2)

          SELECT CASE ( rflu_getfacekind(global,c1k,c2k) ) 
            CASE ( face_kind_aa ) ! Actual-actual face
              IF ( c1 == icg ) THEN 
                icg = c2
              ELSE 
                icg = c1
              END IF ! c1                                    
            CASE ( face_kind_av ) ! Actual-virtual face
              ifl = ifg - pgrid%nFaces + pgrid%nFacesAV

              iborder = pgrid%avf2b(1,ifl)

              pborder => pgrid%borders(iborder)
              
              pborder%nPclsSend = pborder%nPclsSend + 1

              IF ( pborder%nPclsSend > pborder%nPclsSendMax ) &
                CALL rflu_mpi_recreatebufferipclsend(pregion,pborder)
          
              pborder%iPclSend(1,pborder%nPclsSend) = ipcl
              pborder%iPclSend(2,pborder%nPclsSend) = pgrid%avf2b(2,ifl)            

              pplag%aiv(aiv_plag_status,ipcl) = plag_status_comm
              pplag%arv(arv_plag_distot,ipcl) = disttot

              IF ( c1 == icg ) THEN 
                icg = c2
              ELSE 
                icg = c1
              END IF ! c1
              
              pplag%aiv(aiv_plag_icells,ipcl) = icg
              
              ipatch = pgrid%avf2p(ifl)
              
              IF ( ipatch /= crazy_value_int ) THEN ! Transform coordinates
                ppatch => pregion%patches(ipatch)
                
                IF ( ppatch%bcType == bc_periodic ) THEN 
                  CALL rflu_relp_transformvector(pregion,ppatch,xlocold, &
                                                 ylocold,zlocold)
                  CALL rflu_relp_transformvector(pregion,ppatch,xlocnew, &
                                                 ylocnew,zlocnew)                                                
                ELSE
                  CALL errorstop(global,err_reached_default,__line__)
                END IF ! pPatch%bcType
              END IF ! iPatch         
              
              EXIT trajloop                  
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! RFLU_GetFaceKind                    
                  
! ----- Trajectory intersects boundary face          
          
        ELSE ! Boundary face
          ppatch => pregion%patches(iloc)
          
          fnx = ppatch%fn(xcoord,ifg)
          fny = ppatch%fn(ycoord,ifg)
          fnz = ppatch%fn(zcoord,ifg) 
                      
          theta = xtraj*fnx + ytraj*fny + ztraj*fnz 
          
          SELECT CASE ( ppatch%bcType )             
            CASE ( bc_slipwall:bc_slipwall+bc_range )
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_noslipwall:bc_noslipwall+bc_range ) 
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_injection:bc_injection+bc_range ) 
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_outflow:bc_outflow+bc_range )
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag

              pplag%aiv(aiv_plag_status,ipcl) = plag_status_delete
              EXIT trajloop            
            CASE ( bc_farfield:bc_farfield+bc_range )
              pplag%aiv(aiv_plag_status,ipcl) = plag_status_delete
              EXIT trajloop 
            CASE ( bc_symmetry:bc_symmetry+bc_range )
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_virtual:bc_virtual+bc_range )
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE default 
              CALL errorstop(global,err_reached_default,__line__)           
          END SELECT ! pPatch%bcType

        END IF ! iLoc 
      END IF ! distTot
    
! ------------------------------------------------------------------------------
!     Guard against infinite loop 
! ------------------------------------------------------------------------------  
  
      IF ( loopcounter >= limit_infinite_loop ) THEN 
        WRITE(stdout,'(A,1X,A)') solver_name, &
             'Infinite loop encountered in particle cell search algorithm'
        WRITE(stdout,'(A,3X,A)') solver_name,'Module: Lagrangian Particle (PLAG).'

        WRITE(stdout,'(A,3X,A,1X,1PE12.5)') solver_name,'Current time:', &
                                            global%currentTime

        WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
                                         pregion%iRegionGlobal
        WRITE(stdout,'(A,6X,A,11(1X,A))') solver_name,'#', &
                                         '    iPcl     ', &
                                         '    idIni    ', &
                                         '    RegIni   ', &
                                         '    icg      ', &
                                         '  x-location ', &
                                         '  y-location ', &
                                         '  z-Location ', &
                                         '   Energy    ', &
                                         '   Diameter  '

        WRITE(stdout,'(A,4X,4(1X,I8),6(1X,E13.6))') solver_name,ipcl, &
                        pplag%aivOld(aiv_plag_pidini,ipcl),           &
                        pplag%aivOld(aiv_plag_regini,ipcl),           &
                        icg,xloc,yloc,zloc,                           &
                        pplag%cv(cv_plag_ener,ipcl),                  &
                        pplag%dv(dv_plag_diam,ipcl)

        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! loopCounter        
    END DO trajloop 

! ==============================================================================
!   Check that kept particles are indeed located in new cells
! ==============================================================================        

    IF ( (global%checkLevel == check_high) .AND. &
         (pplag%aiv(aiv_plag_status,ipcl) == plag_status_keep) ) THEN
      CALL rflu_ict_testincellfancy(pregion,xlocnew,ylocnew,zlocnew,icg, &
                                    incellcheckflag,ilocout,ifgout)

      IF ( incellcheckflag .EQV. .false. ) THEN 
        WRITE(stderr,'(A,1X,A,1X,I6)') solver_name,'Particle index:',ipcl
        WRITE(stderr,'(A,1X,A,1X,I6)') solver_name,'Cell which failed test:', & 
                                                   icg
        WRITE(stderr,'(A,1X,A,2(1X,I6))') solver_name, &
                                          'Face which failed test:', & 
                                          ilocout,ifgout                                                                 
        WRITE(stderr,'(A,1X,A,3(1X,E23.16))') solver_name, &
                                              'Particle old location:', &
                                              xlocold,ylocold,zlocold
        WRITE(stderr,'(A,1X,A,3(1X,E23.16))') solver_name, &
                                              'Particle new location:', &
                                              xlocnew,ylocnew,zlocnew

        WRITE(errorstring,'(I6)') ipcl
        CALL errorstop(global,err_plag_pcl_not_found,__line__, & 
                       trim(errorstring))           
      END IF ! inCellCheckFlag
    END IF ! global%checkLevel    
    
! ==============================================================================
!   Store new position. IMPORTANT: Need to store new position because it may 
!   have been reflected and old position because it may have been transformed.
! ==============================================================================

    pplag%cv(cv_plag_xpos,ipcl) = xlocnew
    pplag%cv(cv_plag_ypos,ipcl) = ylocnew
    pplag%cv(cv_plag_zpos,ipcl) = zlocnew
     
    pplag%cvOld(cv_plag_xpos,ipcl) = xlocold
    pplag%cvOld(cv_plag_ypos,ipcl) = ylocold
    pplag%cvOld(cv_plag_zpos,ipcl) = zlocold       
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellstrajfast







! ******************************************************************************
!
! Purpose: Determine cells which contain particles by following trajectory 
!   using safe algorithm.
!
! Description: Follow particle path by intersecting particle trajectory with
!   faces and taking appropriate action determined by type of intersected face.
!
! Input: 
!  pRegion      Pointer to region
!  iPclBeg      Beginning particle index
!  iPclEnd      Ending particle index
!
! Output: None.
!
! Notes:
!   1. Need to properly handle particle deletion from outflow or motion
!      between domain during an intermediate RK-stage as its final 
!      position at (n+1) step might not be outside the region.
!
! ******************************************************************************

SUBROUTINE plag_rflu_findcellstrajsafe(pRegion,iPclBeg,iPclEnd)

  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER :: ipclbeg,ipclend
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  LOGICAL :: incellcheckflag
  CHARACTER(CHRLEN) :: errorstring
  INTEGER :: c1,c1k,c2,c2k,errorflag,iborder,icg,ifg,ifgout,ifl,iloc, &
             ilocout,ipatch,ipcl,loopcounter
  REAL(RFREAL) :: dist,disttot,disttotcutoff,eps,fnx,fny,fnz,imagtraj,theta, &
                  xloc,xlocnew,xlocold,xtraj,yloc,ylocnew,ylocold,ytraj,zloc, &
                  zlocnew,zlocold,ztraj
  TYPE(t_border), POINTER :: pborder
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_plag), POINTER :: pplag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'PLAG_RFLU_FindCellsTrajSafe',&
  'PLAG_RFLU_ModFindCells.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pplag => pregion%plag

  eps = epsilon(1.0_rfreal)
  disttotcutoff = 10*eps ! NOTE must be less than tolerICT

! ******************************************************************************
! Loop over particles
! ******************************************************************************

  DO ipcl = ipclbeg,ipclend
    loopcounter = 0 ! Reset loop counter 

! ==============================================================================  
!   Set distance travelled and compute trajectory. NOTE that cannot compute 
!   trajectory (unit) vector from distTot because distTot diminishes as the 
!   particle travels along the trajectory but xLocNew and xLocOld (as well as
!   y and z components) do not change.
! ==============================================================================  

    disttot = pplag%arv(arv_plag_distot,ipcl)
    
    IF ( disttot < disttotcutoff ) THEN  
      cycle 
    END IF ! distTot

    xlocnew = pplag%cv(cv_plag_xpos,ipcl)
    ylocnew = pplag%cv(cv_plag_ypos,ipcl)
    zlocnew = pplag%cv(cv_plag_zpos,ipcl)  

    xlocold = pplag%cvOld(cv_plag_xpos,ipcl)
    ylocold = pplag%cvOld(cv_plag_ypos,ipcl)
    zlocold = pplag%cvOld(cv_plag_zpos,ipcl)
        
    xtraj = xlocnew - xlocold
    ytraj = ylocnew - ylocold      
    ztraj = zlocnew - zlocold

    imagtraj = 1.0_rfreal/(sqrt(xtraj*xtraj + ytraj*ytraj + ztraj*ztraj) + eps)
    
    xtraj = imagtraj*xtraj
    ytraj = imagtraj*ytraj
    ztraj = imagtraj*ztraj    
    
! ------------------------------------------------------------------------------
!   Set variables for trajectory loop
! ------------------------------------------------------------------------------    
    
    xloc = xlocnew - disttot*xtraj
    yloc = ylocnew - disttot*ytraj
    zloc = zlocnew - disttot*ztraj
    
    icg = pplag%aivOld(aiv_plag_icells,ipcl)

! ==============================================================================
!   Loop until distance travelled along trajectory consumed 
! ==============================================================================        
        
    trajloop: DO     
      loopcounter = loopcounter + 1  
        
! ------------------------------------------------------------------------------
!     Find appropriate intersection and associated face  
! ------------------------------------------------------------------------------
                    
      CALL rflu_computelinecellxsectsafe(pregion,xloc,yloc,zloc,xtraj,ytraj, &
                                         ztraj,icg,dist,iloc,ifg)  
                                                                                              
! ------------------------------------------------------------------------------
!     Update total distance travelled                             
! ------------------------------------------------------------------------------
                                   
      disttot = disttot - dist
      
! ------------------------------------------------------------------------------
!     Check whether have remaining total distance 
! ------------------------------------------------------------------------------

! --- No distance remaining: Found cell containing new location ----------------

      IF ( disttot <= 0.0_rfreal ) THEN       
        pplag%aiv(aiv_plag_icells,ipcl) = icg
        
        EXIT trajloop
        
! --- Distance remaining: Keep searching ---------------------------------------        
        
      ELSE 

! ----- Trajectory intersects interior face

        IF ( iloc == 0 ) THEN           
          c1 = pgrid%f2c(1,ifg)
          c2 = pgrid%f2c(2,ifg)
          
          c1k = rflu_getglobalcellkind(global,pgrid,c1)
          c2k = rflu_getglobalcellkind(global,pgrid,c2)

          SELECT CASE ( rflu_getfacekind(global,c1k,c2k) ) 
            CASE ( face_kind_aa ) ! Actual-actual face
              IF ( c1 == icg ) THEN 
                icg = c2
              ELSE 
                icg = c1
              END IF ! c1         
            CASE ( face_kind_av ) ! Actual-virtual face
              ifl = ifg - pgrid%nFaces + pgrid%nFacesAV

              iborder = pgrid%avf2b(1,ifl)

              pborder => pgrid%borders(iborder)
              
              pborder%nPclsSend = pborder%nPclsSend + 1

              IF ( pborder%nPclsSend > pborder%nPclsSendMax ) &
                CALL rflu_mpi_recreatebufferipclsend(pregion,pborder)
              
              pborder%iPclSend(1,pborder%nPclsSend) = ipcl
              pborder%iPclSend(2,pborder%nPclsSend) = pgrid%avf2b(2,ifl)            

              pplag%aiv(aiv_plag_status,ipcl) = plag_status_comm
              pplag%arv(arv_plag_distot,ipcl) = disttot
              
              IF ( c1 == icg ) THEN 
                icg = c2
              ELSE 
                icg = c1
              END IF ! c1
              
              pplag%aiv(aiv_plag_icells,ipcl) = icg  
              
              ipatch = pgrid%avf2p(ifl)
              
              IF ( ipatch /= crazy_value_int ) THEN ! Transform coordinates
                ppatch => pregion%patches(ipatch)
                
                IF ( ppatch%bcType == bc_periodic ) THEN 
                  CALL rflu_relp_transformvector(pregion,ppatch,xlocold, &
                                                 ylocold,zlocold)
                  CALL rflu_relp_transformvector(pregion,ppatch,xlocnew, &
                                                 ylocnew,zlocnew)                                                
                ELSE
                  CALL errorstop(global,err_reached_default,__line__)
                END IF ! pPatch%bcType
              END IF ! iPatch

              EXIT trajloop                             
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! RFLU_GetFaceKind                    
                  
! ----- Trajectory intersects boundary face          
          
        ELSE ! Boundary face
          ppatch => pregion%patches(iloc)
          
          fnx = ppatch%fn(xcoord,ifg)
          fny = ppatch%fn(ycoord,ifg)
          fnz = ppatch%fn(zcoord,ifg) 
                      
          theta = xtraj*fnx + ytraj*fny + ztraj*fnz 
          
          SELECT CASE ( ppatch%bcType )             
            CASE ( bc_slipwall:bc_slipwall+bc_range )
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_noslipwall:bc_noslipwall+bc_range ) 
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_injection:bc_injection+bc_range ) 
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_outflow:bc_outflow+bc_range )
              IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
                CALL plag_gathersurfstats(pregion,pplag,ppatch%statsPlag,&
                                          ifg,ipcl,theta)
              END IF ! pPatch%plotStatsFlag

              pplag%aiv(aiv_plag_status,ipcl) = plag_status_delete
              EXIT trajloop            
            CASE ( bc_farfield:bc_farfield+bc_range )
              pplag%aiv(aiv_plag_status,ipcl) = plag_status_delete
              EXIT trajloop 
            CASE ( bc_symmetry:bc_symmetry+bc_range )
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE ( bc_virtual:bc_virtual+bc_range )
              CALL plag_reflectparticledata(ppatch,pplag,ifg,ipcl,xlocold, &
                                            ylocold,zlocold,xlocnew,ylocnew, &
                                            zlocnew,xtraj,ytraj,ztraj)
            CASE default 
              CALL errorstop(global,err_reached_default,__line__)           
          END SELECT ! pPatch%bcType

        END IF ! iLoc 
      END IF ! distTot
    
! ------------------------------------------------------------------------------
!     Guard against infinite loop 
! ------------------------------------------------------------------------------  

      IF ( loopcounter >= limit_infinite_loop ) THEN 
        WRITE(stdout,'(A,1X,A)') solver_name, & 
             'Infinite loop encountered in particle cell search algorithm'
        WRITE(stdout,'(A,3X,A)') solver_name,'Module: Lagrangian Particle (PLAG).'        
          
        WRITE(stdout,'(A,3X,A,1X,1PE12.5)') solver_name,'Current time:', &
                                            global%currentTime              
                                            
        WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', & 
                                         pregion%iRegionGlobal  
        WRITE(stdout,'(A,6X,A,11(1X,A))') solver_name,'#', &
                                         '    iPcl     ', &
                                         '    idIni    ', &
                                         '    RegIni   ', &
                                         '    icg      ', &
                                         '  x-location ', &
                                         '  y-location ', &
                                         '  z-Location ', &
                                         '   Energy    ', &
                                         '   Diameter  '       

        WRITE(stdout,'(A,4X,4(1X,I8),6(1X,E13.6))') solver_name,ipcl, & 
                        pplag%aivOld(aiv_plag_pidini,ipcl),           &
                        pplag%aivOld(aiv_plag_regini,ipcl),           &
                        icg,xloc,yloc,zloc,                           &
                        pplag%cv(cv_plag_ener,ipcl),                  &
                        pplag%dv(dv_plag_diam,ipcl)

        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! loopCounter        
    END DO trajloop 

! ==============================================================================
!   Check that kept particles are indeed located in new cells
! ==============================================================================        

    IF ( (global%checkLevel == check_high) .AND. &
         (pplag%aiv(aiv_plag_status,ipcl) == plag_status_keep) ) THEN
      CALL rflu_ict_testincellfancy(pregion,xlocnew,ylocnew,zlocnew,icg, &
                                    incellcheckflag,ilocout,ifgout)

      IF ( incellcheckflag .EQV. .false. ) THEN 
        WRITE(stderr,'(A,1X,A,1X,I6)') solver_name,'Particle index:',ipcl
        WRITE(stderr,'(A,1X,A,1X,I6)') solver_name,'Cell which failed test:', & 
                                                   icg
        WRITE(stderr,'(A,1X,A,2(1X,I6))') solver_name, &
                                          'Face which failed test:', & 
                                          ilocout,ifgout                                                                 
        WRITE(stderr,'(A,1X,A,3(1X,E23.16))') solver_name, &
                                              'Particle old location:', &
                                              xlocold,ylocold,zlocold
        WRITE(stderr,'(A,1X,A,3(1X,E23.16))') solver_name, &
                                              'Particle new location:', &
                                              xlocnew,ylocnew,zlocnew

        WRITE(errorstring,'(I6)') ipcl
        CALL errorstop(global,err_plag_pcl_not_found,__line__, & 
                       trim(errorstring))           
      END IF ! inCellCheckFlag
    END IF ! global%checkLevel    
    
! ==============================================================================
!   Store new position. IMPORTANT: Need to store new position because it may 
!   have been reflected and old position because it may have been transformed.
! ==============================================================================

    pplag%cv(cv_plag_xpos,ipcl) = xlocnew
    pplag%cv(cv_plag_ypos,ipcl) = ylocnew
    pplag%cv(cv_plag_zpos,ipcl) = zlocnew
     
    pplag%cvOld(cv_plag_xpos,ipcl) = xlocold
    pplag%cvOld(cv_plag_ypos,ipcl) = ylocold
    pplag%cvOld(cv_plag_zpos,ipcl) = zlocold                          
  END DO ! iPcl

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_rflu_findcellstrajsafe





! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE plag_rflu_modfindcells

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_RFLU_ModFindCells.F90,v $
! Revision 1.16  2008/12/06 08:44:35  mtcampbe
! Updated license.
!
! Revision 1.15  2008/11/19 22:17:48  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.14  2006/08/18 21:11:35  fnajjar
! Enabled periodicity, cosmetics
!
! Revision 1.13  2006/05/05 02:15:49  haselbac
! Bug fixes: Incorrect computation of traj, missing EXIT, wrong IF
!
! Revision 1.12  2006/05/02 17:46:23  fnajjar
! Allowed surface statistics to be gathered on active patches
!
! Revision 1.11  2006/04/07 15:19:24  haselbac
! Removed tabs
!
! Revision 1.10  2005/12/24 21:39:47  haselbac
! Adapted to changes in ICT
!
! Revision 1.9  2005/12/19 16:47:45  fnajjar
! Added verbosity around error trap for infinite loop counter
!
! Revision 1.8  2005/12/14 21:21:21  fnajjar
! Added call for dynamic allocation of iPclSend
!
! Revision 1.7  2005/12/02 20:07:50  fnajjar
! Added particle reflection for BC_VIRTUAL
!
! Revision 1.6  2005/11/12 00:34:08  fnajjar
! Bug fix for iMagTraj in Fast and Safe
!
! Revision 1.5  2005/09/20 15:46:35  fnajjar
! Fixed bug in TrajFast and added error trap for iPclSend overflow
!
! Revision 1.4  2005/07/18 20:47:41  fnajjar
! Aligned fast routine with safe adding proper construct for trajectory
!
! Revision 1.3  2005/05/18 22:17:50  fnajjar
! Added PLAG_RFLU_ComputeDistTot, changed comp of distTot and xyzLoc, 
! infrast for comm
!
! Revision 1.2  2005/05/02 21:59:38  haselbac
! Preparation for parallelization of FindCells routines, cosmetics
!
! Revision 1.1  2005/04/27 14:55:55  fnajjar
! Initial import
!
! ******************************************************************************