RFLU_ModStencilsUtils.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 utility routines to construct stencils.
!
! Description: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: RFLU_ModStencilsUtils.F90,v 1.14 2008/12/06 08:44:24 mtcampbe Exp $
!
! Copyright: (c) 2003-2006 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modstencilsutils

  USE modglobal, ONLY: t_global 
  USE modparameters
  USE moddatatypes
  USE moderror
  USE modbndpatch, ONLY: t_patch
  USE moddatastruct, ONLY: t_region  
  USE modgrid, ONLY: t_grid
  USE modstencil, ONLY: t_stencil
  USE modsortsearch
  USE modmpi
  
  USE rflu_modcellfaceedgeinfo  
  
  IMPLICIT NONE

  PRIVATE
  PUBLIC :: rflu_addbfaces, & 
            rflu_addcelllayer, &  
            rflu_addcelllayer_1d, &
            rflu_addcelllayer_1d_g, &           
            rflu_addfacevertneighbs, &
            rflu_computestencilsize, &            
            rflu_computestencilweights, &             
            rflu_sortbfaces

       
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
   
  CHARACTER(CHRLEN) :: RCSIdentString = & 
    '$RCSfile: RFLU_ModStencilsUtils.F90,v $ $Revision: 1.14 $'        
             
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  
  
  


! ******************************************************************************
!
! Purpose: Extend basic stencil by adding boundary faces.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   nBFaceMembsMaxTemp  Maximum allowed number of boundary faces in stencil
!   nCellMembs          Number of cell members in stencil
!   cellMembs           List of cell members of stencil
!
! Output:
!   nBFaceMembs         Number of boundary faces in stencil
!   bFaceMembs          List of boundary faces in stencil
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_addbfaces(pRegion,nBFaceMembsMaxTemp,nCellMembs,cellMembs, &
                            nbfacemembs,bfacemembs)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: nbfacemembsmaxtemp,ncellmembs
    INTEGER, INTENT(OUT) :: nbfacemembs
    INTEGER, INTENT(IN) :: cellmembs(ncellmembs)
    INTEGER, INTENT(OUT) :: bfacemembs(2,nbfacemembsmaxtemp)                 
    TYPE(t_region), POINTER :: pregion      

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

    INTEGER :: errorflag,ifl,iloc,ipatch,isg,isl      
    TYPE(t_global), POINTER :: global 
    TYPE(t_grid), POINTER :: pgrid          
    TYPE(t_patch), POINTER :: ppatch      

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

    global => pregion%global

    pgrid => pregion%grid

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

    CALL registerfunction(global,'RFLU_AddBFaces',&
  'RFLU_ModStencilsUtils.F90')     

! ******************************************************************************
!   Build sorted bf2c lists (and sorting keys) so can be searched later on
! ******************************************************************************        

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( (ppatch%bcType == bc_noslipwall_hflux) .OR. & 
           (ppatch%bcType == bc_noslipwall_temp ) .OR. & 
           (ppatch%bcType == bc_slipwall        ) .OR. & 
           (ppatch%bcType == bc_injection       ) ) THEN 
        IF ( ppatch%nBFaces > 0 ) THEN 
          ALLOCATE(ppatch%bf2cSorted(ppatch%nBFaces),stat=errorflag)
          global%error = errorflag
          IF ( global%error /= err_none ) THEN 
            CALL errorstop(global,err_allocate,__line__,'pPatch%bf2cSorted')
          END IF ! global%error

          ALLOCATE(ppatch%bf2cSortedKeys(ppatch%nBFaces),stat=errorflag)
          global%error = errorflag
          IF ( global%error /= err_none ) THEN 
            CALL errorstop(global,err_allocate,__line__,'pPatch%bf2cSortedKeys')
          END IF ! global%error

          DO ifl = 1,ppatch%nBFaces
            ppatch%bf2cSorted(ifl) = ppatch%bf2c(ifl)
          END DO ! ifl

          CALL quicksortinteger(ppatch%bf2cSorted,ppatch%nBFaces)

          DO ifl = 1,ppatch%nBFaces
            CALL binarysearchinteger(ppatch%bf2cSorted,ppatch%nBFaces, &
                                     ppatch%bf2c(ifl),iloc)

            IF ( iloc /= element_not_found ) THEN 
              ppatch%bf2cSortedKeys(iloc) = ifl
            ELSE 
              CALL errorstop(global,err_binary_search,__line__)
            END IF ! iloc                                   
          END DO ! ifl 
        END IF ! pPatch%nBFaces       
      END IF ! pPatch%bcType        
    END DO ! iPatch 

! ******************************************************************************
!   Initialize list
! ******************************************************************************        

    nbfacemembs = 0

    DO isl = 1,nbfacemembsmaxtemp
      bfacemembs(1,isl) = crazy_value_int  
      bfacemembs(2,isl) = crazy_value_int                 
    END DO ! isl

! ******************************************************************************
!   Find boundary cells in stencil
! ******************************************************************************        

    islloop: DO isl = 1,ncellmembs
      isg = cellmembs(isl)

      DO ipatch = 1,pgrid%nPatches
        ppatch => pregion%patches(ipatch)

        IF ( (ppatch%bcType == bc_noslipwall_hflux) .OR. & 
             (ppatch%bcType == bc_noslipwall_temp ) .OR. & 
             (ppatch%bcType == bc_slipwall        ) .OR. & 
             (ppatch%bcType == bc_injection       ) ) THEN 
          IF ( ppatch%nBFaces > 0 ) THEN 
            CALL binarysearchinteger(ppatch%bf2cSorted,ppatch%nBFaces,isg, & 
                                     iloc)

            IF ( iloc /= element_not_found ) THEN
              IF ( nbfacemembs < nbfacemembsmaxtemp ) THEN  
                nbfacemembs = nbfacemembs + 1

                bfacemembs(1,nbfacemembs) = ipatch
                bfacemembs(2,nbfacemembs) = ppatch%bf2cSortedKeys(iloc)
              ELSE 
                EXIT islloop
              END IF ! nBFaceMembs
            END IF ! iloc
          END IF ! pPatch%nBFaces 
        END IF ! pPatch%bcType                             
      END DO ! iPatch                    
    END DO islloop          

! ******************************************************************************
!   Deallocate memory for sorted bf2c lists
! ******************************************************************************        

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( (ppatch%bcType == bc_noslipwall_hflux) .OR. & 
           (ppatch%bcType == bc_noslipwall_temp ) .OR. & 
           (ppatch%bcType == bc_slipwall        ) .OR. & 
           (ppatch%bcType == bc_injection       ) ) THEN
        IF ( ppatch%nBFaces > 0 ) THEN 
          DEALLOCATE(ppatch%bf2cSorted,stat=errorflag)
          global%error = errorflag
          IF ( global%error /= err_none ) THEN 
            CALL errorstop(global,err_deallocate,__line__,'pPatch%bf2cSorted')
          END IF ! global%error

          DEALLOCATE(ppatch%bf2cSortedKeys,stat=errorflag)
          global%error = errorflag
          IF ( global%error /= err_none ) THEN 
            CALL errorstop(global,err_deallocate,__line__, &
                           'pPatch%bf2cSortedKeys')
          END IF ! global%error 
        END IF ! pPatch%nBFaces
      END IF ! pPatch%bcType       
    END DO ! iPatch 

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

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_addbfaces
  
  
  





! ******************************************************************************
!
! Purpose: Extend basic stencil in only ONE direction
!
! Description: Loop over cells in the last layer of stencil and add cells
!  whose normal vector is aligned along a particular direction, i.e. x, y, z
!  Not just add vertex-neighbors because we only want to extend the stencil 
!  in one direction
!
! Input:
!   global              Pointer to global data
!   pGrid               Pointer to grid
!   stencilSizeMax      Maximum allowed size of stencil
!   ixg                 Index of cell for which cells are added
!   degr                Degree of stencil before addition of cell layer
!   x2csBeg             Beginning index for the last layer of cells in stencil
!   x2csEnd             Ending index for the last layer of cells in stencil
!   x2cs                Stencil before addition of cell layer
!   fnDir               Direction of stencil
!
! Output:
!   degr                Degree of stencil after addition of cell layer
!   x2cs                Stencil after addition of cell layer
!
! Notes: None.
!
! ******************************************************************************
 
  SUBROUTINE rflu_addcelllayer_1d(global,pGrid,stencilSizeMax,ixg,degr, &
                                  x2csbeg,x2csend,x2cs,fndir)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: fndir,ixg,stencilsizemax,x2csbeg,x2csend
    INTEGER, INTENT(INOUT) :: degr
    INTEGER, INTENT(INOUT) :: x2cs(stencilsizemax)
    TYPE(t_global), POINTER :: global
    TYPE(t_grid), POINTER :: pgrid

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

    INTEGER :: c1,c2,degrhelp,icg,icg2,icl,ifg,ifl,iloc,ipatch,isl,iv2c,ivg, &
               ivl,nfaces
    INTEGER :: x2cshelp(stencilsizemax),x2cssort(stencilsizemax)
    REAL(RFREAL) :: fn

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

    CALL registerfunction(global,'RFLU_AddCellLayer_1D',&
  'RFLU_ModStencilsUtils.F90')

! ******************************************************************************
!   Add layer of cells to existing stencil in 1D, store in separate array
! ******************************************************************************

! ==============================================================================      
!   Initialize
! ==============================================================================

    degrhelp = 0 ! renamed because 'degr' is already used

    DO isl = 1,stencilsizemax
      x2cshelp(isl) = 0 ! store new layer's members
      x2cssort(isl) = 0 ! duplicated array for x2cs
    END DO ! isl

    DO isl = 1,degr
      x2cssort(isl) = x2cs(isl)
    END DO ! isl

    CALL quicksortinteger(x2cssort(1:degr),degr)

! ==============================================================================
!   Loop over existing members in the last layer and find candidates for stencil 
! ==============================================================================

    outerloop: DO isl = x2csbeg,x2csend
      icg = x2cs(isl)
      icl = pgrid%cellGlob2Loc(2,icg)

      nfaces = SIZE(pgrid%hex2f,2)

      DO ifl = 1,nfaces
        ipatch = pgrid%hex2f(1,ifl,icl)
        ifg    = pgrid%hex2f(2,ifl,icl)
      
        IF ( ipatch == 0 ) THEN ! Interior face
          c1 = pgrid%f2c(1,ifg)
          c2 = pgrid%f2c(2,ifg)

          IF ( c1 /= ixg .AND. c2 /= ixg ) THEN
            IF ( c1 == icg ) THEN
              fn =  pgrid%fn(fndir,ifg)
            ELSE IF ( c2 == icg ) THEN
              fn = -pgrid%fn(fndir,ifg)
            ELSE ! defensive programming
              CALL errorstop(global,err_reached_default,__line__)
            END IF ! c1

            IF ( abs(fn) >= 0.999_rfreal ) THEN 
              IF ( c1 == icg ) THEN 
                icg2 = c2
              ELSE IF ( c2 == icg ) THEN 
                icg2 = c1
              ELSE ! defensive programming
                CALL errorstop(global,err_reached_default,__line__)
              END IF ! c1

              CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)
              
              IF ( iloc == element_not_found ) THEN 
                IF ( degrhelp > 0 ) THEN 
                  CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                           icg2,iloc)
                  IF ( iloc == element_not_found ) THEN
                    IF ( degrhelp < stencilsizemax ) THEN  
                      degrhelp = degrhelp + 1
                      x2cshelp(degrhelp) = icg2
                      IF ( degrhelp > 1 ) THEN
                        CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                      END IF ! degrHelp
                    ELSE
                      EXIT outerloop
                    END IF ! degrHelp
                  END IF ! iloc
                ELSE ! First member                
                  degrhelp = degrhelp + 1
                  x2cshelp(degrhelp) = icg2
                END IF ! degrHelp
              END IF ! iloc
            END IF ! ABS(fn)
          END IF ! c1, c2                 
        ELSE IF ( ipatch > 0 ) THEN ! Boundary face  
! TO DO 
!       Currently do not add boundary faces for constrained reconstruction
! END TO DO 
        ELSE ! Defensive programming
          CALL errorstop(global,err_reached_default,__line__)
        END IF ! iPatch
      END DO ! ifl
    END DO outerloop

! ******************************************************************************
!   Merge new layer with existing stencil
! ******************************************************************************

    DO isl = 1,degrhelp
      IF ( degr < stencilsizemax ) THEN 
        degr = degr + 1
        x2cs(degr) = x2cshelp(isl)
      END IF ! degr
    END DO ! isl    

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

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_addcelllayer_1d






! ******************************************************************************
!
! Purpose: Extend basic stencil in ONE direction by adding vertex neighbours 
!   of cells based on geometry.
!
! Description: Loop over cells in last layer of stencil and add all cells
!  which share vertices with that layer of cells and are aligned with the
!  given direction.
!
! Input:
!   global              Pointer to global data
!   pGrid               Pointer to grid
!   stencilSizeMax      Maximum allowed size of stencil
!   ixg                 Index of cell for which cells are added
!   degr                Degree of stencil before addition of cell layer
!   x2csBeg             Beginning index for last layer of cells in stencil
!   x2csEnd             Ending index for last layer of cells in stencil
!   x2cs                Stencil before addition of cell layer
!   rc                  Location for which stencil is constructed
!   dir                 Direction in which stencil is to be constructed
!
! Output:
!   degr                Degree of stencil after addition of cell layer
!   x2cs                Stencil after addition of cell layer
!
! Notes: None.
!
! ******************************************************************************
 
  SUBROUTINE rflu_addcelllayer_1d_g(global,pGrid,stencilSizeMax,ixg,degr, &
                                    x2csbeg,x2csend,x2cs,rc,dir)

    USE rflu_modgeometrytools, ONLY: rflu_testvectorcartaxisaligned

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: dir,ixg,stencilsizemax,x2csbeg,x2csend
    INTEGER, INTENT(INOUT) :: degr
    INTEGER, INTENT(INOUT) :: x2cs(stencilsizemax)
    REAL(RFREAL), INTENT(IN) :: rc(xcoord:zcoord)
    TYPE(t_global), POINTER :: global
    TYPE(t_grid), POINTER :: pgrid

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

    LOGICAL :: alignflag
    INTEGER :: degrhelp,icg,icg2,icl,ict,iloc,isl,ivg,ivl,iv2c
    INTEGER :: x2cshelp(stencilsizemax),x2cssort(stencilsizemax)
    REAL(RFREAL) :: dr(xcoord:zcoord)

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

    CALL registerfunction(global,'RFLU_AddCellLayer_1D_G',&
  'RFLU_ModStencilsUtils.F90')

! ******************************************************************************
!   Add layer of cells to existing stencil, store in separate array
! ******************************************************************************

! ==============================================================================      
!   Initialize
! ==============================================================================

    degrhelp = 0

    DO isl = 1,stencilsizemax
      x2cshelp(isl) = 0
      x2cssort(isl) = 0
    END DO ! isl

    DO isl = 1,degr
      x2cssort(isl) = x2cs(isl)
    END DO ! isl

    CALL quicksortinteger(x2cssort(1:degr),degr)

! ==============================================================================
!   Loop over existing members in last layer and find candidates for stencil 
! ==============================================================================

    outerloop: DO isl = x2csbeg,x2csend                       
      icg = x2cs(isl)
      ict = rflu_getglobalcelltype(global,pgrid,icg)
      icl = pgrid%cellGlob2Loc(2,icg)

      SELECT CASE ( ict ) 

! ------------------------------------------------------------------------------
!       Hexahedra
! ------------------------------------------------------------------------------                    

        CASE ( cell_type_hex ) 
          DO ivl = 1,8
            ivg = pgrid%hex2v(ivl,icl)

            DO iv2c = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
              icg2 = pgrid%v2c(iv2c)

              IF ( icg2 /= ixg ) THEN 
                CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)

                IF ( iloc == element_not_found ) THEN 
                  IF ( degrhelp > 0 ) THEN 
                    CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                             icg2,iloc)
                    IF ( iloc == element_not_found ) THEN
                      IF ( degrhelp < stencilsizemax ) THEN
                        dr(xcoord) = pgrid%cofg(xcoord,icg2) - rc(xcoord)
                        dr(ycoord) = pgrid%cofg(ycoord,icg2) - rc(ycoord)
                        dr(zcoord) = pgrid%cofg(zcoord,icg2) - rc(zcoord)
                        
                        alignflag = rflu_testvectorcartaxisaligned(global,dr,dir)
                        
                        IF ( alignflag .EQV. .true. ) THEN                                                 
                          degrhelp = degrhelp + 1
                          x2cshelp(degrhelp) = icg2

                          IF ( degrhelp > 1 ) THEN 
                            CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                          END IF ! degrHelp
                        END IF ! alignFlag
                      ELSE
                        EXIT outerloop
                      END IF ! degrHelp                       
                    END IF ! iloc                          
                  ELSE
                    dr(xcoord) = pgrid%cofg(xcoord,icg2) - rc(xcoord)
                    dr(ycoord) = pgrid%cofg(ycoord,icg2) - rc(ycoord)
                    dr(zcoord) = pgrid%cofg(zcoord,icg2) - rc(zcoord)

                    alignflag = rflu_testvectorcartaxisaligned(global,dr,dir)

                    IF ( alignflag .EQV. .true. ) THEN
                      degrhelp = degrhelp + 1
                      x2cshelp(degrhelp) = icg2
                    END IF ! alignFlag                       
                  END IF ! degrHelp
                END IF ! iloc
              END IF ! icg2

            END DO ! iv2c
          END DO ! ivl

! ------------------------------------------------------------------------------
!       Default
! ------------------------------------------------------------------------------            

        CASE default 
          CALL errorstop(global,err_reached_default,__line__)                
      END SELECT ! ict
    END DO outerloop 

! ******************************************************************************
!   Merge new layer with existing stencil
! ******************************************************************************

    DO isl = 1,degrhelp
      IF ( degr < stencilsizemax ) THEN 
        degr = degr + 1
        x2cs(degr) = x2cshelp(isl)
      END IF ! degr
    END DO ! isl    

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

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_addcelllayer_1d_g



  
  

! ******************************************************************************
!
! Purpose: Extend basic stencil by adding vertex neighbours of cells.
!
! Description: Loop over cells in last layer of stencil and add all cells
!  which share vertices with that layer of cells.
!
! Input:
!   global              Pointer to global data
!   pGrid               Pointer to grid
!   stencilSizeMax      Maximum allowed size of stencil
!   ixg                 Index of cell for which cells are added
!   degr                Degree of stencil before addition of cell layer
!   x2csBeg             Beginning index for last layer of cells in stencil
!   x2csEnd             Ending index for last layer of cells in stencil
!   x2cs                Stencil before addition of cell layer
!
! Output:
!   degr                Degree of stencil after addition of cell layer
!   x2cs                Stencil after addition of cell layer
!
! Notes: None.
!
! ******************************************************************************
 
  SUBROUTINE rflu_addcelllayer(global,pGrid,stencilSizeMax,ixg,degr,x2csBeg, &
                               x2csend,x2cs)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: ixg,stencilsizemax,x2csbeg,x2csend
    INTEGER, INTENT(INOUT) :: degr
    INTEGER, INTENT(INOUT) :: x2cs(stencilsizemax)
    TYPE(t_global), POINTER :: global
    TYPE(t_grid), POINTER :: pgrid

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

    INTEGER :: degrhelp,icg,icg2,icl,ict,iloc,isl,ivg,ivl,iv2c
    INTEGER :: x2cshelp(stencilsizemax),x2cssort(stencilsizemax)

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

    CALL registerfunction(global,'RFLU_AddCellLayer',&
  'RFLU_ModStencilsUtils.F90')

! ******************************************************************************
!   Add layer of cells to existing stencil, store in separate array
! ******************************************************************************

! ==============================================================================      
!   Initialize
! ==============================================================================

    degrhelp = 0

    DO isl = 1,stencilsizemax
      x2cshelp(isl) = 0
      x2cssort(isl) = 0
    END DO ! isl

    DO isl = 1,degr
      x2cssort(isl) = x2cs(isl)
    END DO ! isl

    CALL quicksortinteger(x2cssort(1:degr),degr)

! ==============================================================================
!   Loop over existing members in last layer and find candidates for stencil 
! ==============================================================================

    outerloop: DO isl = x2csbeg,x2csend                       
      icg = x2cs(isl)
      ict = rflu_getglobalcelltype(global,pgrid,icg)
      icl = pgrid%cellGlob2Loc(2,icg)

      SELECT CASE ( ict ) 

! ------------------------------------------------------------------------------
!       Tetrahedra
! ------------------------------------------------------------------------------        

        CASE ( cell_type_tet ) 
          DO ivl = 1,4
            ivg = pgrid%tet2v(ivl,icl)

            DO iv2c = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
              icg2 = pgrid%v2c(iv2c)

              IF ( icg2 /= ixg ) THEN 
                CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)

                IF ( iloc == element_not_found ) THEN 
                  IF ( degrhelp > 0 ) THEN 
                    CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                             icg2,iloc)

                    IF ( iloc == element_not_found ) THEN
                      IF ( degrhelp < stencilsizemax ) THEN  
                        degrhelp = degrhelp + 1
                        x2cshelp(degrhelp) = icg2

                        IF ( degrhelp > 1 ) THEN 
                          CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                        END IF ! degrHelp
                      ELSE
                        EXIT outerloop
                      END IF ! degrHelp                       
                    END IF ! iloc                         
                  ELSE
                    degrhelp = degrhelp + 1
                    x2cshelp(degrhelp) = icg2
                  END IF ! degrHelp
                END IF ! iloc
              END IF ! icg2 

            END DO ! iv2c
          END DO ! ivl

! ------------------------------------------------------------------------------
!       Hexahedra
! ------------------------------------------------------------------------------                    

        CASE ( cell_type_hex ) 
          DO ivl = 1,8
            ivg = pgrid%hex2v(ivl,icl)

            DO iv2c = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
              icg2 = pgrid%v2c(iv2c)

              IF ( icg2 /= ixg ) THEN 
                CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)

                IF ( iloc == element_not_found ) THEN 
                  IF ( degrhelp > 0 ) THEN 
                    CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                             icg2,iloc)
                    IF ( iloc == element_not_found ) THEN
                      IF ( degrhelp < stencilsizemax ) THEN  
                        degrhelp = degrhelp + 1
                        x2cshelp(degrhelp) = icg2

                        IF ( degrhelp > 1 ) THEN 
                          CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                        END IF ! degrHelp
                      ELSE
                        EXIT outerloop
                      END IF ! degrHelp                       
                    END IF ! iloc                          
                  ELSE
                    degrhelp = degrhelp + 1
                    x2cshelp(degrhelp) = icg2
                  END IF ! degrHelp
                END IF ! iloc
              END IF ! icg2

            END DO ! iv2c
          END DO ! ivl

! ------------------------------------------------------------------------------
!       Prisms
! ------------------------------------------------------------------------------                    

        CASE ( cell_type_pri )           
          DO ivl = 1,6
            ivg = pgrid%pri2v(ivl,icl)

            DO iv2c = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
              icg2 = pgrid%v2c(iv2c)

               IF ( icg2 /= ixg ) THEN 
                CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)

                IF ( iloc == element_not_found ) THEN 
                  IF ( degrhelp > 0 ) THEN 
                    CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                             icg2,iloc)

                    IF ( iloc == element_not_found ) THEN
                      IF ( degrhelp < stencilsizemax ) THEN  
                        degrhelp = degrhelp + 1
                        x2cshelp(degrhelp) = icg2

                        IF ( degrhelp > 1 ) THEN 
                          CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                        END IF ! degrHelp
                      ELSE
                        EXIT outerloop
                      END IF ! degrHelp                       
                    END IF ! iloc                          
                  ELSE
                    degrhelp = degrhelp + 1
                    x2cshelp(degrhelp) = icg2
                  END IF ! degrHelp
                END IF ! iloc
              END IF ! icg2

            END DO ! iv2c
          END DO ! ivl  

! ------------------------------------------------------------------------------
!       Pyramids
! ------------------------------------------------------------------------------

        CASE ( cell_type_pyr ) 
          DO ivl = 1,5
            ivg = pgrid%pyr2v(ivl,icl)

            DO iv2c = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
              icg2 = pgrid%v2c(iv2c)

              IF ( icg2 /= ixg ) THEN 
                CALL binarysearchinteger(x2cssort(1:degr),degr,icg2,iloc)

                IF ( iloc == element_not_found ) THEN 
                  IF ( degrhelp > 0 ) THEN 
                    CALL binarysearchinteger(x2cshelp(1:degrhelp),degrhelp, &
                                             icg2,iloc)

                    IF ( iloc == element_not_found ) THEN
                      IF ( degrhelp < stencilsizemax ) THEN  
                        degrhelp = degrhelp + 1
                        x2cshelp(degrhelp) = icg2

                        IF ( degrhelp > 1 ) THEN 
                          CALL quicksortinteger(x2cshelp(1:degrhelp),degrhelp)
                        END IF ! degrHelp
                      ELSE
                        EXIT outerloop
                      END IF ! degrHelp                       
                    END IF ! iloc                           
                  ELSE
                    degrhelp = degrhelp + 1
                    x2cshelp(degrhelp) = icg2
                  END IF ! degrHelp
                END IF ! iloc
              END IF ! icg2

            END DO ! iv2c              
          END DO ! ivl

! ------------------------------------------------------------------------------
!       Default
! ------------------------------------------------------------------------------            

        CASE default 
          CALL errorstop(global,err_reached_default,__line__)                
      END SELECT ! ict
    END DO outerloop 

! ******************************************************************************
!   Merge new layer with existing stencil
! ******************************************************************************

    DO isl = 1,degrhelp
      IF ( degr < stencilsizemax ) THEN 
        degr = degr + 1
        x2cs(degr) = x2cshelp(isl)
      END IF ! degr
    END DO ! isl    

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

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_addcelllayer






! ******************************************************************************
!
! Purpose: Construct initial stencil by adding cells sharing vertices which
!   define face. 
!
! Description: None.
!
! Input:
!   global              Pointer to global data
!   pGrid               Pointer to grid
!   stencilSizeMax      Maximum allowed size of stencil
!   f2v                 List of vertices which define face
!
! Output:
!   degr                Number of cells sharing vertices defining face
!   x2cs                List of cells sharing vertices defining face
!
! Notes:
!   1. It is assumed that there are no members in the stencil when this 
!      routine is called. 
!
! ******************************************************************************

  SUBROUTINE rflu_addfacevertneighbs(global,pGrid,stencilSizeMax,f2v,degr, &
                                     x2cs)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: stencilsizemax
    INTEGER, INTENT(IN) :: f2v(4)
    INTEGER, INTENT(INOUT) :: degr
    INTEGER, INTENT(INOUT) :: x2cs(stencilsizemax)
    TYPE(t_global), POINTER :: global
    TYPE(t_grid), POINTER :: pgrid

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

    INTEGER :: icg,iloc,ivg,ivl,ivl2

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

    CALL registerfunction(global,'RFLU_AddFaceVertNeighbs',&
  'RFLU_ModStencilsUtils.F90')

! ******************************************************************************
!   Add layer of cells to existing stencil, store in separate array
! ******************************************************************************

! ==============================================================================      
!   Initialize
! ==============================================================================

    degr = 0

! ==============================================================================
!   Loop over vertices of face and add cells  
! ==============================================================================

    outerloop: DO ivl = 1,4
      ivg = f2v(ivl)

      IF ( ivg /= vert_none ) THEN 

        DO ivl2 = pgrid%v2cInfo(v2c_beg,ivg),pgrid%v2cInfo(v2c_end,ivg)
          icg = pgrid%v2c(ivl2)

          IF ( degr /= 0 ) THEN 
            CALL binarysearchinteger(x2cs(1:degr),degr,icg,iloc)

            IF ( iloc == element_not_found ) THEN
              IF ( degr < stencilsizemax ) THEN  
                degr = degr + 1
                x2cs(degr) = icg

                IF ( degr > 1 ) THEN 
                  CALL quicksortinteger(x2cs(1:degr),degr)
                END IF ! degr
              ELSE
                EXIT outerloop
              END IF ! degr
            END IF ! iloc
          ELSE 
            degr = degr + 1

            x2cs(degr) = icg
          END IF ! degr
        END DO ! ivl2          

      END IF ! ivg
    END DO outerloop
                                       
! ******************************************************************************
!   End
! ******************************************************************************

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_addfacevertneighbs






! *******************************************************************************
!
! Purpose: Compute stencil size for given order.
!
! Description: None.
!
! Input:
!   global      Global pointer
!   dimens      Dimensionality
!   factor      Safety factor
!   order       Polynomial order
!
! Output: 
!   RFLU_ComputeStencilSize     (Minimum) number of points in stencil
!
! Notes: 
!   1. Include additional support (+1) to allow for interpolation.
!   2. Make larger than necessary to allow for least-squares if order > 1.
!   3. If order = 0, do not include factor. 
!
! ******************************************************************************

  INTEGER FUNCTION rflu_computestencilsize(global,dimens,factor,order)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: dimens,factor,order
    TYPE(t_global), POINTER :: global

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

    CHARACTER(CHRLEN) :: rcsidentstring

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

    CALL registerfunction(global,'RFLU_ComputeStencilSize',&
  'RFLU_ModStencilsUtils.F90')

! ******************************************************************************
!   Compute stencil size
! ******************************************************************************

    IF ( order > 0 ) THEN 
      SELECT CASE ( dimens ) 
        CASE ( 2 ) 
          rflu_computestencilsize = factor*(order + 1)*(order + 2)/2
        CASE ( 3 ) 
          rflu_computestencilsize = factor*(order + 1)*(order + 2)*(order + 3)/6        
        CASE default
          CALL errorstop(global,err_dimens_invalid,__line__)
      END SELECT ! dimens
    ELSE 
      rflu_computestencilsize = 1
    END IF ! order

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

    CALL deregisterfunction(global)

  END FUNCTION rflu_computestencilsize






! ******************************************************************************
!
! Purpose: Compute stencil weights.
!
! Description: Compute stencil weights in iterative fashion. If weights are
!   rejected, requested order is lowered and the weights are recomputed. Once
!   order is equal to zero, the weights are set instead of computed.
!
! Input:
!   global              Global pointer
!   dimens              Dimensionality
!   wtsMode             Mode of computing weights
!   scalMode            Mode of scaling
!   derivDegree         Degree of derivative (0 = interpolation)
!   order               Order of accuracy of approximation
!   nRows               Number of members in stencil
!   dr                  Relative position vectors of members in stencil
!
! Output: 
!   orderNominal        Order of accuracy of approximation
!   wts                 Stencil weights
!   sCount              Number of singular values
!
! Notes:
!   1. It is important to note that the order is the polynomial order of 
!      approximation, that is, order equal to 1 means a linear polynomial
!      which may be first- or second-order accurate...
!   2. Still need to check for singular systems because although initial 
!      stencils ought to be non-singular, they may become singular with grid
!      motion.
!   3. Although this is a routine to compute weights, it is not located in the
!      module RFLU_ComputeWeights, because it requires the function 
!      RFLU_ComputeStencilSize, and this leads to a circular dependency:
!      RFLU_ModWeights depends on RFLU_ModStencils (because of 
!      RFLU_ComputeStencilSize), and RFLU_ModStencils depends on 
!      RFLU_ModWeights (because of RFLU_ComputeStencilWeights)
!
! ******************************************************************************

  SUBROUTINE rflu_computestencilweights(global,dimens,wtsMode,scalMode, &
                                        derivdegree,ordernominal,nrows,dr, &
                                        wts,scount)

    USE modtools, ONLY: compfact

    USE modinterfaces, ONLY: rflu_invertmatrixsvd

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: derivdegree,dimens,nrows,scalmode,wtsmode
    INTEGER, INTENT(IN) :: ordernominal
    INTEGER, INTENT(OUT), OPTIONAL :: scount 
    REAL(RFREAL), INTENT(IN) :: dr(xcoord:zcoord,nrows)
    REAL(RFREAL), INTENT(OUT) :: wts(:,:)
    TYPE(t_global), POINTER :: global

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

    CHARACTER(CHRLEN) :: rcsidentstring
    INTEGER :: errorflag,isl,j,ncols,order,p,q,r
    REAL(RFREAL) :: colscal,term,wtsnorm,wtsnormmax,wtsnormmin
    REAL(RFREAL) :: rowscal(nrows)
    REAL(RFREAL) :: drcopy(xcoord:zcoord,nrows)
    REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: a,ainv

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

    CALL registerfunction(global,'RFLU_ComputeStencilWeights',&
  'RFLU_ModStencilsUtils.F90')

    wtsnormmin = 1.0_rfreal 
    wtsnormmax = 10.0_rfreal 

! ******************************************************************************
!   Check that stencil large enough, reduce order if necessary
! ******************************************************************************

    IF ( nrows < rflu_computestencilsize(global,dimens,1,ordernominal) ) THEN 
      SELECT CASE ( dimens )
        CASE ( 2 ) 
          SELECT CASE ( nrows )
            CASE ( 0:2 ) ! Constant
              order = 0 
            CASE ( 3:5 ) ! Linear
              order = 1
            CASE ( 6:9 ) ! Quadratic
              order = 2
            CASE ( 10:14 ) ! Cubic
              order = 3
            CASE default ! More than cubic
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! nRows           
        CASE ( 3 ) 
          SELECT CASE ( nrows )
            CASE ( 0:3 ) ! Constant
              order = 0 
            CASE ( 4:9 ) ! Linear
              order = 1
            CASE ( 10:19 ) ! Quadratic
              order = 2
            CASE ( 20:34 ) ! Cubic
              order = 3
            CASE default ! More than cubic
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! nRows
        CASE default 
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! dimens
    ELSE 
      order = ordernominal 
    END IF ! nRowsEff

! ******************************************************************************
!   Row scaling
! ******************************************************************************

    IF ( scalmode == compwts_scal_none ) THEN 
      DO isl = 1,nrows
        rowscal(isl) = 1.0_rfreal
      END DO ! isl  
    ELSE IF ( scalmode == compwts_scal_invdist ) THEN 
      DO isl = 1,nrows
        rowscal(isl) = 1.0_rfreal/sqrt(dr(xcoord,isl)*dr(xcoord,isl) & 
                                     + dr(ycoord,isl)*dr(ycoord,isl) & 
                                     + dr(zcoord,isl)*dr(zcoord,isl))
      END DO ! isl
    ELSE 
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! scalMode

! ******************************************************************************
!   Column scaling
! ******************************************************************************

    colscal = maxval(abs(dr(xcoord:zcoord,1:nrows)))

    DO isl = 1,nrows
      drcopy(xcoord,isl) = dr(xcoord,isl)/colscal
      drcopy(ycoord,isl) = dr(ycoord,isl)/colscal
      drcopy(zcoord,isl) = dr(zcoord,isl)/colscal                
    END DO ! isl

! ******************************************************************************
!   Compute weights iteratively
! ******************************************************************************

    DO 

! ==============================================================================
!     Order greater than zero, so compute weights
! ==============================================================================

      IF ( order > 0 ) THEN 
        ncols = rflu_computestencilsize(global,dimens,1,order)        

! ------------------------------------------------------------------------------
!       Allocate temporary memory
! ------------------------------------------------------------------------------

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

        ALLOCATE(ainv(ncols,nrows),stat=errorflag)
        global%error = errorflag   
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_allocate,__line__,'aInv')
        END IF ! global%error      

! ------------------------------------------------------------------------------
!       Build matrix
! ------------------------------------------------------------------------------

        SELECT CASE ( dimens )
          CASE ( 2 ) 
            DO isl = 1,nrows
              a(isl,1) = rowscal(isl)

              j = 2

              DO p = 1,order
                DO q = 0,p
                  term = rowscal(isl)/(compfact(p-q)*compfact(q))                 
                  a(isl,j) = term*drcopy(xcoord,isl)**(p-q) &
                                 *drcopy(ycoord,isl)**q   
                  j = j + 1
                END DO ! q
              END DO ! p                    
            END DO ! isl             
          CASE ( 3 ) 
            DO isl = 1,nrows
              a(isl,1) = rowscal(isl)

              j = 2

              DO p = 1,order
                DO q = 0,p
                  DO r = 0,q
                    term = rowscal(isl)/(compfact(p-q)*compfact(q-r)*compfact(r))                 
                    a(isl,j) = term*drcopy(xcoord,isl)**(p-q) &
                                   *drcopy(ycoord,isl)**(q-r) &
                                   *drcopy(zcoord,isl)**r 
                    j = j + 1
                  END DO ! r
                END DO ! q
              END DO ! p                    
            END DO ! isl
          CASE default
            CALL errorstop(global,err_reached_default,__line__)
        END SELECT ! dimens    

! ------------------------------------------------------------------------------
!       Invert matrix
! ------------------------------------------------------------------------------

        CALL rflu_invertmatrixsvd(global,nrows,ncols,a,ainv,scount)

! ------------------------------------------------------------------------------
!       Extract weights from inverse. NOTE appearance of scalings.
! ------------------------------------------------------------------------------

        IF ( derivdegree == deriv_degree_0 ) THEN 
          DO isl = 1,nrows
            wts(1,isl) = rowscal(isl)*ainv(1,isl) 
          END DO ! isl
        ELSE IF ( derivdegree == deriv_degree_1 ) THEN
          SELECT CASE ( dimens )
            CASE ( 2 ) 
              DO isl = 1,nrows
                wts(xcoord,isl) = rowscal(isl)*ainv(2,isl)/colscal
                wts(ycoord,isl) = rowscal(isl)*ainv(3,isl)/colscal
                wts(zcoord,isl) = 0.0_rfreal
              END DO ! isl                
            CASE ( 3 ) 
              DO isl = 1,nrows
                wts(xcoord,isl) = rowscal(isl)*ainv(2,isl)/colscal
                wts(ycoord,isl) = rowscal(isl)*ainv(3,isl)/colscal
                wts(zcoord,isl) = rowscal(isl)*ainv(4,isl)/colscal 
              END DO ! isl   
            CASE default 
              CALL errorstop(global,err_reached_default,__line__)   
          END SELECT ! dimens
        ELSE
          CALL errorstop(global,err_reached_default,__line__)
        END IF ! derivDegree

! ------------------------------------------------------------------------------
!       Check whether weights make sense
! ------------------------------------------------------------------------------

        wtsnorm = 0.0_rfreal

        DO isl = 1,nrows
          wtsnorm = wtsnorm + abs(ainv(1,isl))
        END DO ! isl

! ------------------------------------------------------------------------------
!       Deallocate temporary memory
! ------------------------------------------------------------------------------

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

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

! ------------------------------------------------------------------------------
!       Decide whether to accept or reject weights
! ------------------------------------------------------------------------------

        IF ( scount == 0 .AND. & 
             wtsnorm > wtsnormmin .AND. wtsnorm < wtsnormmax ) THEN 
          EXIT
        ELSE 
          IF ( wtsmode == compwts_mode_adapt ) THEN 
            order = order - 1                     
          ELSE 
            EXIT
          END IF ! wtsMode
        END IF ! errNorm     

! ==============================================================================
!     Order equal to zero, so set weights 
! ==============================================================================

      ELSE           
        IF ( derivdegree == deriv_degree_0 ) THEN 
          DO isl = 1,nrows
            wts(1,isl) = 1.0_rfreal/REAL(nRows,KIND=RFREAL) 
          END DO ! isl
        ELSE IF ( derivdegree == deriv_degree_1 ) THEN 
          DO isl = 1,nrows
            wts(xcoord,isl) = 0.0_rfreal
            wts(ycoord,isl) = 0.0_rfreal
            wts(zcoord,isl) = 0.0_rfreal
          END DO ! isl      
        ELSE
          CALL errorstop(global,err_reached_default,__line__)
        END IF ! derivDegree         

        EXIT 
      END IF ! order                                
    END DO ! <empty>

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

    CALL deregisterfunction(global)

  END SUBROUTINE rflu_computestencilweights







! ******************************************************************************
!
! Purpose: Sort boundary faces by increasing distance from location at which
!   gradient is reconstructed.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   xyz                 Coordinates at which gradient is reconstructed
!   nBFaceMembs         Number of boundary face stencil members
!   bFaceMembs          Boundary face stencil members
!
! Output:
!   bFaceMembs          Boundary face stencil members
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_sortbfaces(pRegion,xyz,nBFaceMembs,bFaceMembs)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    INTEGER, INTENT(IN) :: nbfacemembs
    INTEGER, INTENT(INOUT) :: bfacemembs(2,nbfacemembs)
    REAL(RFREAL), INTENT(IN) :: xyz(xcoord:zcoord) 
    TYPE(t_region), POINTER :: pregion      

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

    INTEGER :: ifl,ipatch,isl
    INTEGER :: key(nbfacemembs)
    INTEGER :: bfacemembsbak(2,nbfacemembs)      
    REAL(RFREAL) :: dist(nbfacemembs)      
    TYPE(t_global), POINTER :: global 
    TYPE(t_grid), POINTER :: pgrid          
    TYPE(t_patch), POINTER :: ppatch      

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

    global => pregion%global
    pgrid  => pregion%grid

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

    CALL registerfunction(global,'RFLU_SortBFaces',&
  'RFLU_ModStencilsUtils.F90')     

! ******************************************************************************
!   Back up data structure so can copy later on
! ******************************************************************************

    DO isl = 1,nbfacemembs 
      bfacemembsbak(1,isl) = bfacemembs(1,isl)
      bfacemembsbak(2,isl) = bfacemembs(2,isl)
    END DO ! isl   

! ******************************************************************************
!   Compute squared distance and build key for sorting
! ******************************************************************************

    DO isl = 1,nbfacemembs
      ipatch = bfacemembs(1,isl)
      ifl    = bfacemembs(2,isl)

      ppatch => pregion%patches(ipatch)

      key(isl)  = isl
      dist(isl) = (ppatch%fc(xcoord,ifl) - xyz(xcoord))**2 & 
                + (ppatch%fc(ycoord,ifl) - xyz(ycoord))**2 & 
                + (ppatch%fc(zcoord,ifl) - xyz(zcoord))**2
    END DO ! isl    
  
! ******************************************************************************
!   Sort according to increasing distance
! ******************************************************************************

    CALL quicksortrfrealinteger(dist,key,nbfacemembs)  

    DO isl = 1,nbfacemembs
      bfacemembs(1,isl) = bfacemembsbak(1,key(isl))
      bfacemembs(2,isl) = bfacemembsbak(2,key(isl))        
    END DO ! isl        

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

    CALL deregisterfunction(global)      

  END SUBROUTINE rflu_sortbfaces





! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modstencilsutils


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModStencilsUtils.F90,v $
! Revision 1.14  2008/12/06 08:44:24  mtcampbe
! Updated license.
!
! Revision 1.13  2008/11/19 22:17:35  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.12  2007/03/06 18:08:37  haselbac
! Added function to add cell layer in 1d based on geometry (not on hex2f)
!
! Revision 1.11  2006/04/07 15:19:20  haselbac
! Removed tabs
!
! Revision 1.10  2006/04/01 15:56:37  haselbac
! Bug fix: Only add bfaces if have actual faces
!
! Revision 1.9  2006/03/25 21:57:59  haselbac
! Changed bf2cSorted list to contain only actual entries (bcos of sype changes)
!
! Revision 1.8  2006/01/06 22:14:15  haselbac
! Added new routine for adding cell layer in 1d
!
! Revision 1.7  2005/10/09 13:27:18  haselbac
! Temporarily disable addition of v bfaces to stencils
!
! Revision 1.6  2005/10/05 14:10:01  haselbac
! Added routines which used to be in RFLU_ModStencils
!
! Revision 1.5  2005/04/15 15:07:05  haselbac
! Added routines to build and destroy c2c stencil in CSR format
!
! Revision 1.4  2004/12/29 21:09:44  haselbac
! Cosmetics only
!
! Revision 1.3  2004/01/22 16:03:59  haselbac
! Made contents of modules PRIVATE, only procs PUBLIC, to avoid errors on ALC 
! and titan
!
! Revision 1.2  2003/12/05 16:53:11  haselbac
! Added cell-to-cell CSR routines for parallel higher-order scheme
!
! Revision 1.1  2003/12/04 03:29:22  haselbac
! Initial revision
!
! ******************************************************************************