RFLU_ModWeights.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 to compute stencil weights.
!
! Description: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: RFLU_ModWeights.F90,v 1.16 2008/12/06 08:44:25 mtcampbe Exp $
!
! Copyright: (c) 2003-2006 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modweights

  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 modmixture, ONLY: t_mixt_input
  USE modmpi
    
  IMPLICIT NONE
     
  PRIVATE
  PUBLIC :: rflu_createwtsbf2cwrapper, &
            rflu_createwtsc2cwrapper, &
            rflu_createwtsf2cwrapper, &
            rflu_computewtsbf2cwrapper, &
            rflu_computewtsc2cwrapper, &
            rflu_computewtsf2cwrapper, &
            rflu_computewtsx2c_1d, &
            rflu_destroywtsbf2cwrapper, & 
            rflu_destroywtsc2cwrapper, &
            rflu_destroywtsf2cwrapper

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

  CONTAINS
  
  
  



! *******************************************************************************
!
! Purpose: Create weights for boundary face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!   orderInput  Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_createwtsbf2c(pRegion,pPatch,orderInput)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: orderinput  
    TYPE(t_patch), POINTER :: ppatch         
    TYPE(t_region), POINTER :: pregion

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsBF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Creating boundary face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Modify order so that can run with first-order scheme
! ******************************************************************************

    order = max(orderinput,1)

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsbf2c(pregion,ppatch)

! ******************************************************************************
!   Allocate memory
! ******************************************************************************

    SELECT CASE ( order ) 
    
! ==============================================================================
!     Linear approximation
! ==============================================================================    
    
      CASE ( 1 ) 
        IF ( ppatch%bcType /= bc_virtual ) THEN  
          SELECT CASE ( pregion%mixtInput%dimens ) 
            CASE ( 1 ) ! NOTE should never reach here
            CASE ( 2 ) 
              DO ifl = 1,ppatch%nBFaces
                ALLOCATE(ppatch%bf2cs(ifl)%xyzMoms(xyz_mom_11:xyz_mom_33), & 
                         stat=errorflag)  
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_allocate,__line__, & 
                                 'pPatch%bf2cs%xyzMoms')
                END IF ! global%error
              END DO ! ifl            
            CASE ( 3 )            
              DO ifl = 1,ppatch%nBFaces
                ALLOCATE(ppatch%bf2cs(ifl)%xyzMoms(xyz_mom_11:xyz_mom_44), & 
                         stat=errorflag)  
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_allocate,__line__, & 
                                 'pPatch%bf2cs%xyzMoms')
                END IF ! global%error
              END DO ! ifl
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! dimens
        END IF ! pPatch%bcType            

! ==============================================================================
!     Default
! ==============================================================================    

      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order
    
! ******************************************************************************
!   End
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Creating boundary face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsbf2c  
  
  
  
  
  

! *******************************************************************************
!
! Purpose: Wrapper routine for creating boundary-face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_createwtsbf2cwrapper(pRegion,pPatch,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_patch), POINTER :: ppatch
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsBF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensBFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_createwtsbf2c(pregion,ppatch,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensBFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsbf2cwrapper
  
  
  
  
  


! *******************************************************************************
!
! Purpose: Create weights for cell-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_createwtsc2c(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order  
    TYPE(t_region), POINTER :: pregion

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsC2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Creating cell-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsc2c(pregion)

! ******************************************************************************
!   Allocate memory
! ******************************************************************************

    SELECT CASE ( order ) 

! ==============================================================================
!     Linear approximation
! ==============================================================================    

      CASE ( 1 )
        SELECT CASE ( pregion%mixtInput%dimens ) 
          CASE ( 1 ) 
          CASE ( 2 )
            DO icg = 1,pgrid%nCellsTot 
              ALLOCATE(pgrid%c2cs(icg)%xyzMoms(xyz_mom_11:xyz_mom_33), & 
                       stat=errorflag)                   
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__, &
                               'pGrid%c2cs%xyzMoms')
              END IF ! global%error                  
            END DO ! icg           
          CASE ( 3 ) 
            DO icg = 1,pgrid%nCellsTot 
              ALLOCATE(pgrid%c2cs(icg)%xyzMoms(xyz_mom_11:xyz_mom_44), & 
                       stat=errorflag)                   
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__, &
                               'pGrid%c2cs%xyzMoms')
              END IF ! global%error                  
            END DO ! icg
          CASE default 
            CALL errorstop(global,err_reached_default,__line__)
        END SELECT ! pRegion%mixtInput%dimens  

! ==============================================================================
!     Default
! ==============================================================================    

      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Creating cell-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsc2c








! *******************************************************************************
!
! Purpose: Wrapper routine for creating cell-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_createwtsc2cwrapper(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsC2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensCells )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_createwtsc2c(pregion,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensCells

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsc2cwrapper








! *******************************************************************************
!
! Purpose: Create weights for face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   orderInput  Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_createwtsf2c(pRegion,orderInput)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: orderinput  
    TYPE(t_region), POINTER :: pregion

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Creating face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Modify order so that can run with first-order scheme
! ******************************************************************************

    order = max(orderinput,1)

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsf2c(pregion)

! ******************************************************************************
!   Allocate memory
! ******************************************************************************

    SELECT CASE ( order ) 
    
! ==============================================================================
!     Linear approximation
! ==============================================================================    
    
      CASE ( 1 ) 
        SELECT CASE ( pregion%mixtInput%dimens ) 
          CASE ( 1 ) ! NOTE should never reach here
          CASE ( 2 ) 
            DO ifg = 1,pgrid%nFaces
              ALLOCATE(pgrid%f2cs(ifg)%xyzMoms(xyz_mom_11:xyz_mom_33), & 
                       stat=errorflag)                   
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__, &
                               'pGrid%f2cs%xyzMoms')
              END IF ! global%error        
            END DO ! ifg          
          CASE ( 3 )           
            DO ifg = 1,pgrid%nFaces
              ALLOCATE(pgrid%f2cs(ifg)%xyzMoms(xyz_mom_11:xyz_mom_44), & 
                       stat=errorflag)                   
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__, &
                               'pGrid%f2cs%xyzMoms')
              END IF ! global%error        
            END DO ! ifg
          CASE default 
            CALL errorstop(global,err_reached_default,__line__)
        END SELECT ! dimens
    
! ==============================================================================
!     Default
! ==============================================================================    

      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order
    
! ******************************************************************************
!   End
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Creating face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsf2c







! *******************************************************************************
!
! Purpose: Wrapper routine for creating face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_createwtsf2cwrapper(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_CreateWtsF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_createwtsf2c(pregion,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_createwtsf2cwrapper






! *******************************************************************************
!
! Purpose: Compute weights for boundary face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!   orderInput  Desired order 
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_computewtsbf2c(pRegion,pPatch,orderInput)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: orderinput  
    TYPE(t_patch), POINTER :: ppatch    
    TYPE(t_region), POINTER :: pregion

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

    INTEGER :: errorflag,nmembs,icg,ifl,isl,order
    REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: dr    
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid        

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsBF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Computing boundary face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Modify order so that can run with first-order scheme
! ******************************************************************************

    order = max(orderinput,1)
    
! ******************************************************************************
!   Compute weights
! ******************************************************************************

    SELECT CASE ( order ) 
    
! ==============================================================================
!     Linear approximation
! ==============================================================================    
    
      CASE ( 1 )         
        IF ( ppatch%bcType /= bc_virtual ) THEN 
          SELECT CASE ( pregion%mixtInput%dimens ) 
            CASE ( 1 ) 
            CASE ( 2 ) 
              DO ifl = 1,ppatch%nBFaces
                nmembs = ppatch%bf2cs(ifl)%nCellMembs

                ALLOCATE(dr(xcoord:ycoord,nmembs),stat=errorflag)
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_allocate,__line__,'dr')
                END IF ! global%error

                DO isl = 1,nmembs
                  icg = ppatch%bf2cs(ifl)%cellMembs(isl)

                  dr(xcoord,isl) = pgrid%cofg(xcoord,icg) & 
                                 - ppatch%fc(xcoord,ifl)
                  dr(ycoord,isl) = pgrid%cofg(ycoord,icg) & 
                                 - ppatch%fc(ycoord,ifl)
                END DO ! isl        

                CALL rflu_computestencilmoments2d1(global,nmembs,dr, & 
                                                   ppatch%bf2cs(ifl)%xyzMoms)

                DEALLOCATE(dr,stat=errorflag)
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_deallocate,__line__,'dr')
                END IF ! global%error              
              END DO ! ifl            
            CASE ( 3 ) 
              DO ifl = 1,ppatch%nBFaces
                nmembs = ppatch%bf2cs(ifl)%nCellMembs

                ALLOCATE(dr(xcoord:zcoord,nmembs),stat=errorflag)
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_allocate,__line__,'dr')
                END IF ! global%error

                DO isl = 1,nmembs
                  icg = ppatch%bf2cs(ifl)%cellMembs(isl)

                  dr(xcoord,isl) = pgrid%cofg(xcoord,icg) & 
                                 - ppatch%fc(xcoord,ifl)
                  dr(ycoord,isl) = pgrid%cofg(ycoord,icg) & 
                                 - ppatch%fc(ycoord,ifl)
                  dr(zcoord,isl) = pgrid%cofg(zcoord,icg) & 
                                 - ppatch%fc(zcoord,ifl)                    
                END DO ! isl        

                CALL rflu_computestencilmoments3d1(global,nmembs,dr, & 
                                                   ppatch%bf2cs(ifl)%xyzMoms)

                DEALLOCATE(dr,stat=errorflag)
                global%error = errorflag
                IF ( global%error /= err_none ) THEN 
                  CALL errorstop(global,err_deallocate,__line__,'dr')
                END IF ! global%error              
              END DO ! ifl
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! pRegion%mixtInput%dimens
        END IF ! pPatch%bcType 

! ==============================================================================
!     Default
! ==============================================================================    
        
      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Computing boundary face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsbf2c








! *******************************************************************************
!
! Purpose: Wrapper routine for computing boundary-face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch 
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_computewtsbf2cwrapper(pRegion,pPatch,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_patch), POINTER :: ppatch
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsBF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensBFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_computewtsbf2c(pregion,ppatch,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensBFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsbf2cwrapper







! *******************************************************************************
!
! Purpose: Compute weights for cell-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_computewtsc2c(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order  
    TYPE(t_region), POINTER :: pregion

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

    INTEGER :: errorflag,nmembs,icg,icg2,isl
    REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: dr    
    TYPE(t_grid), POINTER :: pgrid        
    TYPE(t_global), POINTER :: global  

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsC2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Computing cell-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

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

    SELECT CASE ( order ) 

! ==============================================================================
!     Linear approximation
! ==============================================================================    

      CASE ( 1 ) 
        SELECT CASE ( pregion%mixtInput%dimens ) 
          CASE ( 1 ) ! NOTE should never reach here
          CASE ( 2 ) 
            DO icg = 1,pgrid%nCellsTot
              nmembs = pgrid%c2cs(icg)%nCellMembs
            
              ALLOCATE(dr(xcoord:ycoord,nmembs),stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__,'dr')
              END IF ! global%error

              DO isl = 1,nmembs
                icg2 = pgrid%c2cs(icg)%cellMembs(isl)

                dr(xcoord,isl) = pgrid%cofg(xcoord,icg2)-pgrid%cofg(xcoord,icg)
                dr(ycoord,isl) = pgrid%cofg(ycoord,icg2)-pgrid%cofg(ycoord,icg)
              END DO ! isl        

              CALL rflu_computestencilmoments2d1(global,nmembs,dr, & 
                                                 pgrid%c2cs(icg)%xyzMoms)

              DEALLOCATE(dr,stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_deallocate,__line__,'dr')
              END IF ! global%error        
            END DO ! icg          
          CASE ( 3 )           
            DO icg = 1,pgrid%nCellsTot
              nmembs = pgrid%c2cs(icg)%nCellMembs            
            
              ALLOCATE(dr(xcoord:zcoord,nmembs),stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__,'dr')
              END IF ! global%error

              DO isl = 1,nmembs
                icg2 = pgrid%c2cs(icg)%cellMembs(isl)

                dr(xcoord,isl) = pgrid%cofg(xcoord,icg2)-pgrid%cofg(xcoord,icg)
                dr(ycoord,isl) = pgrid%cofg(ycoord,icg2)-pgrid%cofg(ycoord,icg)
                dr(zcoord,isl) = pgrid%cofg(zcoord,icg2)-pgrid%cofg(zcoord,icg)                    
              END DO ! isl        

              CALL rflu_computestencilmoments3d1(global,nmembs,dr, & 
                                                 pgrid%c2cs(icg)%xyzMoms)

              DEALLOCATE(dr,stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_deallocate,__line__,'dr')
              END IF ! global%error        
            END DO ! icg
          CASE default
            CALL errorstop(global,err_reached_default,__line__)
        END SELECT ! pRegion%mixtInput%dimens

! ==============================================================================
!     Default
! ==============================================================================    

      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Computing cell-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsc2c








! *******************************************************************************
!
! Purpose: Wrapper routine for computing cell-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_computewtsc2cwrapper(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsC2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensCells )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_computewtsc2c(pregion,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensCells

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsc2cwrapper







! *******************************************************************************
!
! Purpose: Compute weights for face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   orderInput  Desired order 
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_computewtsf2c(pRegion,orderInput)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: orderinput  
    TYPE(t_region), POINTER :: pregion

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

    INTEGER :: errorflag,nmembs,icg,ifg,isl,order
    REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: dr    
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid     

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Computing face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Modify order so that can run with first-order scheme
! ******************************************************************************

    order = max(orderinput,1)
    
! ******************************************************************************
!   Compute weights
! ******************************************************************************

    SELECT CASE ( order ) 
    
! ==============================================================================
!     Linear approximation
! ==============================================================================    
    
      CASE ( 1 ) 
      
! ------------------------------------------------------------------------------
!       Interior faces
! ------------------------------------------------------------------------------      
      
        SELECT CASE ( pregion%mixtInput%dimens ) 
          CASE ( 1 ) ! NOTE should never reach here
          CASE ( 2 )
            DO ifg = 1,pgrid%nFaces
              nmembs = pgrid%f2cs(ifg)%nCellMembs
              
              ALLOCATE(dr(xcoord:ycoord,nmembs),stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__,'dr')
              END IF ! global%error

              DO isl = 1,nmembs
                icg = pgrid%f2cs(ifg)%cellMembs(isl)

                dr(xcoord,isl) = pgrid%cofg(xcoord,icg) - pgrid%fc(xcoord,ifg)
                dr(ycoord,isl) = pgrid%cofg(ycoord,icg) - pgrid%fc(ycoord,ifg)
              END DO ! isl        

              CALL rflu_computestencilmoments2d1(global,nmembs,dr, & 
                                                 pgrid%f2cs(ifg)%xyzMoms)

              DEALLOCATE(dr,stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_deallocate,__line__,'dr')
              END IF ! global%error        
            END DO ! ifg           
          CASE ( 3 ) 
            DO ifg = 1,pgrid%nFaces
              nmembs = pgrid%f2cs(ifg)%nCellMembs            
            
              ALLOCATE(dr(xcoord:zcoord,nmembs),stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_allocate,__line__,'dr')
              END IF ! global%error

              DO isl = 1,nmembs
                icg = pgrid%f2cs(ifg)%cellMembs(isl)

                dr(xcoord,isl) = pgrid%cofg(xcoord,icg)-pgrid%fc(xcoord,ifg)
                dr(ycoord,isl) = pgrid%cofg(ycoord,icg)-pgrid%fc(ycoord,ifg)
                dr(zcoord,isl) = pgrid%cofg(zcoord,icg)-pgrid%fc(zcoord,ifg)                    
              END DO ! isl        

              CALL rflu_computestencilmoments3d1(global,nmembs, & 
                                                 dr,pgrid%f2cs(ifg)%xyzMoms)

              DEALLOCATE(dr,stat=errorflag)
              global%error = errorflag
              IF ( global%error /= err_none ) THEN 
                CALL errorstop(global,err_deallocate,__line__,'dr')
              END IF ! global%error        
            END DO ! ifg
          CASE default
            CALL errorstop(global,err_reached_default,__line__)
        END SELECT ! pRegion%mixtInput%dimens

! ==============================================================================
!     Default
! ==============================================================================    
        
      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! order

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Computing face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsf2c







! *******************************************************************************
!
! Purpose: Wrapper routine for computing face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   order       Desired order
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_computewtsf2cwrapper(pRegion,order)

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: order
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_ComputeWtsF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_computewtsf2c(pregion,order)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsf2cwrapper







! *******************************************************************************
!
! Purpose: Compute weights for 1D x-to-cell stencil.
!
! Description: Compute Lagrangian polynomial stencil weights from Fornberg 
!   algorithm.
!
! Input:
!   global      Pointer to global data
!   m           Order of derivative for which weights are sought
!   nMembs      Number of stencil members
!   x           Locations of stencil members
!   z           Location at which weights are sough
!
! Output:
!   w           Stencil weights (corresponding to c array in Fornberg)
!
! Notes: 
!   1. See Fornberg, SIAM Rev., Vol.40, No.3, pp.685-691, Sep. 1998.
!
! ******************************************************************************
  
  SUBROUTINE rflu_computewtsx2c_1d(global,m,nMembs,x,z,w) 

    IMPLICIT NONE

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

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

    INTEGER, INTENT(IN) :: m,nmembs
    REAL(RFREAL), INTENT(IN) :: z
    REAL(RFREAL), INTENT(IN) :: x(0:nmembs-1)
    REAL(RFREAL), INTENT(OUT) :: w(0:nmembs-1)
    TYPE(t_global), POINTER :: global  

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

    INTEGER :: errorflag,i,j,k,mn,n
    REAL(RFREAL) :: c1,c2,c3,c4,c5 ! Same notation as Fornberg
    REAL(RFREAL) :: c(0:nmembs-1,0:m)    

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

    CALL registerfunction(global,'RFLU_ComputeWtsX2C_1D',&
  'RFLU_ModWeights.F90')

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

    n = nmembs-1

! TO DO 
!   Add check for order (given nMembs)
! END TO DO  

! ******************************************************************************
!   Initialize
! ******************************************************************************
    
    c1 = 1.0_rfreal
    c4 = x(0) - z

    DO k = 0,m
      DO j = 0,n
        c(j,k) = 0.0_rfreal
      END DO ! j
    END DO ! k

    c(0,0) = 1.0_rfreal

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

    DO i = 1,n
      mn = min(i,m)
      c2 = 1.0_rfreal
      c5 = c4
      c4 = x(i) - z
      
      DO j = 0,i-1
        c3 = x(i) - x(j)
        c2 = c2*c3
        
        IF ( j == (i-1) ) THEN
          DO k = mn,1,-1
            c(i,k) = c1*(k*c(i-1,k-1) - c5*c(i-1,k))/c2
          END DO ! k
          
          c(i,0) = -c1*c5*c(i-1,0)/c2
        END IF ! j
        
        DO k = mn,1,-1
          c(j,k) = (c4*c(j,k) - k*c(j,k-1))/c3
        END DO ! k
        
        c(j,0) = c4*c(j,0)/c3
      END DO ! j
      
      c1 = c2
    END DO ! i

    DO j = 0,nmembs-1
      w(j) = c(j,m)
    END DO ! j
    
! ******************************************************************************
!   End
! ******************************************************************************

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_computewtsx2c_1d






! *******************************************************************************
!
! Purpose: Compute weights from coordinate differences.
!
! Description: None.
!
! Input:
!   global      Pointer to global data
!   nMembs      Number of members in stencil
!   dr          Coordinate differences
!
! Output: 
!   xyzMoms     Coordinate moments
!
! Notes:
!   1. Use inverse-distance weighting.
!
! ******************************************************************************

  SUBROUTINE rflu_computestencilmoments2d1(global,nMembs,dr,xyzMoms)
  
    IMPLICIT NONE
    
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================
    
    INTEGER, INTENT(IN) :: nmembs
    REAL(RFREAL), INTENT(INOUT) :: dr(xcoord:ycoord,nmembs)
    REAL(RFREAL), INTENT(INOUT) :: xyzmoms(xyz_mom_11:xyz_mom_33)
    TYPE(t_global), POINTER :: global
  
! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: isl  
    REAL(RFREAL) :: dx,dy,ir11,ir22,ir33,r11,r12,r13,r22,r23, &
                    r33,wt
  
! ******************************************************************************
!   Start
! ******************************************************************************

    CALL registerfunction(global,'RFLU_ComputeStencilMoments2D1',&
  'RFLU_ModWeights.F90')

! ******************************************************************************
!   Compute moments
! ******************************************************************************
                
! ==============================================================================
!   Initialize weights
! ==============================================================================    
    
    xyzmoms(xyz_mom_11) = 0.0_rfreal

    xyzmoms(xyz_mom_12) = 0.0_rfreal
    xyzmoms(xyz_mom_22) = 0.0_rfreal

    xyzmoms(xyz_mom_13) = 0.0_rfreal
    xyzmoms(xyz_mom_23) = 0.0_rfreal
    xyzmoms(xyz_mom_33) = 0.0_rfreal

! ==============================================================================
!   Compute weights
! ==============================================================================

    DO isl = 1,nmembs          
      dx = dr(xcoord,isl) 
      dy = dr(ycoord,isl)                      

      wt = 1.0_rfreal/sqrt(dx*dx + dy*dy)

      dx = wt*dx 
      dy = wt*dy

      xyzmoms(xyz_mom_11) = xyzmoms(xyz_mom_11) + dx*dx

      xyzmoms(xyz_mom_12) = xyzmoms(xyz_mom_12) + dx*dy
      xyzmoms(xyz_mom_22) = xyzmoms(xyz_mom_22) + dy*dy

      xyzmoms(xyz_mom_13) = xyzmoms(xyz_mom_13) + wt*dx
      xyzmoms(xyz_mom_23) = xyzmoms(xyz_mom_23) + wt*dy
      xyzmoms(xyz_mom_33) = xyzmoms(xyz_mom_33) + wt*wt  
    END DO ! isl

    r11  = sqrt(xyzmoms(xyz_mom_11)) 
    ir11 = 1.0_rfreal/r11

    r12  = ir11*xyzmoms(xyz_mom_12)
    r22  = sqrt(xyzmoms(xyz_mom_22) - r12*r12)
    ir22 = 1.0_rfreal/r22

    r13  = ir11*xyzmoms(xyz_mom_13)
    r23  = ir22*(xyzmoms(xyz_mom_23) -  r12*r13           ) 
    r33  = sqrt(xyzmoms(xyz_mom_33)  - (r13*r13 + r23*r23))
    ir33 = 1.0_rfreal/r33

! ==============================================================================
!   Store weights
! ==============================================================================

    xyzmoms(xyz_mom_11) = r11        

    xyzmoms(xyz_mom_12) = r12
    xyzmoms(xyz_mom_22) = r22

    xyzmoms(xyz_mom_13) = r13
    xyzmoms(xyz_mom_23) = r23
    xyzmoms(xyz_mom_33) = r33 
    
! ******************************************************************************
!   End
! ******************************************************************************

    CALL deregisterfunction(global)    
  
  END SUBROUTINE rflu_computestencilmoments2d1








! *******************************************************************************
!
! Purpose: Compute weights from coordinate differences.
!
! Description: None.
!
! Input:
!   global      Pointer to global data
!   nMembs      Number of members in stencil
!   dr          Coordinate differences
!
! Output: 
!   xyzMoms     Coordinate moments
!
! Notes: 
!   1. Use inverse-distance weighting.
!
! ******************************************************************************

  SUBROUTINE rflu_computestencilmoments3d1(global,nMembs,dr,xyzMoms)
  
    IMPLICIT NONE
    
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================
    
    INTEGER, INTENT(IN) :: nmembs
    REAL(RFREAL), INTENT(INOUT) :: dr(xcoord:zcoord,nmembs)
    REAL(RFREAL), INTENT(INOUT) :: xyzmoms(xyz_mom_11:xyz_mom_44)
    TYPE(t_global), POINTER :: global
  
! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: isl  
    REAL(RFREAL) :: dx,dy,dz,ir11,ir22,ir33,r11,r12,r13,r14,r22,r23,r24, &
                    r33,r34,r44,wt
  
! ******************************************************************************
!   Start
! ******************************************************************************

    CALL registerfunction(global,'RFLU_ComputeStencilMoments3D1',&
  'RFLU_ModWeights.F90')

! ******************************************************************************
!   Compute moments
! ******************************************************************************
                
! ==============================================================================
!   Initialize weights
! ==============================================================================    
    
    xyzmoms(xyz_mom_11) = 0.0_rfreal

    xyzmoms(xyz_mom_12) = 0.0_rfreal
    xyzmoms(xyz_mom_22) = 0.0_rfreal

    xyzmoms(xyz_mom_13) = 0.0_rfreal
    xyzmoms(xyz_mom_23) = 0.0_rfreal
    xyzmoms(xyz_mom_33) = 0.0_rfreal

    xyzmoms(xyz_mom_14) = 0.0_rfreal
    xyzmoms(xyz_mom_24) = 0.0_rfreal                              
    xyzmoms(xyz_mom_34) = 0.0_rfreal 
    xyzmoms(xyz_mom_44) = 0.0_rfreal                                             

! ==============================================================================
!   Compute weights
! ==============================================================================

    DO isl = 1,nmembs          
      dx = dr(xcoord,isl) 
      dy = dr(ycoord,isl) 
      dz = dr(zcoord,isl)                     

      wt = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)

      dx = wt*dx 
      dy = wt*dy
      dz = wt*dz

      xyzmoms(xyz_mom_11) = xyzmoms(xyz_mom_11) + dx*dx

      xyzmoms(xyz_mom_12) = xyzmoms(xyz_mom_12) + dx*dy
      xyzmoms(xyz_mom_22) = xyzmoms(xyz_mom_22) + dy*dy

      xyzmoms(xyz_mom_13) = xyzmoms(xyz_mom_13) + dx*dz
      xyzmoms(xyz_mom_23) = xyzmoms(xyz_mom_23) + dy*dz
      xyzmoms(xyz_mom_33) = xyzmoms(xyz_mom_33) + dz*dz 

      xyzmoms(xyz_mom_14) = xyzmoms(xyz_mom_14) + wt*dx
      xyzmoms(xyz_mom_24) = xyzmoms(xyz_mom_24) + wt*dy
      xyzmoms(xyz_mom_34) = xyzmoms(xyz_mom_34) + wt*dz
      xyzmoms(xyz_mom_44) = xyzmoms(xyz_mom_44) + wt*wt  
    END DO ! isl

    r11  = sqrt(xyzmoms(xyz_mom_11)) 
    ir11 = 1.0_rfreal/r11

    r12  = ir11*xyzmoms(xyz_mom_12)
    r22  = sqrt(xyzmoms(xyz_mom_22) - r12*r12)
    ir22 = 1.0_rfreal/r22

    r13  = ir11*xyzmoms(xyz_mom_13)
    r23  = ir22*(xyzmoms(xyz_mom_23) -  r12*r13           ) 
    r33  = sqrt(xyzmoms(xyz_mom_33)  - (r13*r13 + r23*r23))
    ir33 = 1.0_rfreal/r33

    r14  = ir11*xyzmoms(xyz_mom_14)
    r24  = ir22*(xyzmoms(xyz_mom_24) -  r12*r14                     )
    r34  = ir33*(xyzmoms(xyz_mom_34) - (r13*r14 + r23*r24          ))
    r44  = sqrt(xyzmoms(xyz_mom_44)  - (r14*r14 + r24*r24 + r34*r34))

! ==============================================================================
!   Store weights
! ==============================================================================

    xyzmoms(xyz_mom_11) = r11        

    xyzmoms(xyz_mom_12) = r12
    xyzmoms(xyz_mom_22) = r22

    xyzmoms(xyz_mom_13) = r13
    xyzmoms(xyz_mom_23) = r23
    xyzmoms(xyz_mom_33) = r33 

    xyzmoms(xyz_mom_14) = r14
    xyzmoms(xyz_mom_24) = r24
    xyzmoms(xyz_mom_34) = r34
    xyzmoms(xyz_mom_44) = r44
    
! ******************************************************************************
!   End
! ******************************************************************************

    CALL deregisterfunction(global)    
  
  END SUBROUTINE rflu_computestencilmoments3d1





! *******************************************************************************
!
! Purpose: Destroy weights for boundary face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pPatch              Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_destroywtsbf2c(pRegion,pPatch)

    IMPLICIT NONE
    
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    TYPE(t_region), POINTER :: pregion
    TYPE(t_patch), POINTER :: ppatch      

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_DestroyBWtsF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Destroying boundary face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

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

    IF ( ppatch%bcType /= bc_virtual ) THEN 
      DO ifl = 1,ppatch%nBFaces
        DEALLOCATE(ppatch%bf2cs(ifl)%xyzMoms,stat=errorflag)  
        global%error = errorflag
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_deallocate,__line__,'pPatch%bf2cs%xyzMoms')
        END IF ! global%error
      END DO ! ifl
    END IF ! pPatch%bcType

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsbf2c(pregion,ppatch)   

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Destroying boundary face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsbf2c






! *******************************************************************************
!
! Purpose: Wrapper routine for destroying boundary-face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_destroywtsbf2cwrapper(pRegion,pPatch)

    IMPLICIT NONE

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

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

    TYPE(t_patch), POINTER :: ppatch
    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_DestroyWtsBF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensBFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_destroywtsbf2c(pregion,ppatch)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensBFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsbf2cwrapper






! *******************************************************************************
!
! Purpose: Destroy weights for cell-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_destroywtsc2c(pRegion)

    IMPLICIT NONE
    
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    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,'RFLU_DestroyWtsC2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Destroying cell-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Allocate memory
! ******************************************************************************

    DO icg = 1,pgrid%nCellsTot 
      DEALLOCATE(pgrid%c2cs(icg)%xyzMoms,stat=errorflag)                   
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pGrid%c2cs%xyzMoms')
      END IF ! global%error        
    END DO ! icg

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsc2c(pregion)   

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Destroying cell-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsc2c









! *******************************************************************************
!
! Purpose: Wrapper routine for destroying cell-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_destroywtsc2cwrapper(pRegion)

    IMPLICIT NONE

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

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

    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_DestroyWtsC2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensCells )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_destroywtsc2c(pregion)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensCells

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsc2cwrapper








! *******************************************************************************
!
! Purpose: Destroy weights for face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_destroywtsf2c(pRegion)

    IMPLICIT NONE
    
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

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

    TYPE(t_region), POINTER :: pregion

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_DestroyWtsF2C',&
  'RFLU_ModWeights.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Destroying face-to-cell weights...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

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

    DO ifg = 1,pgrid%nFaces
      DEALLOCATE(pgrid%f2cs(ifg)%xyzMoms,stat=errorflag)                   
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pGrid%f2cs%xyzMoms')
      END IF ! global%error        
    END DO ! ifg

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    CALL rflu_nullifywtsf2c(pregion)   

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,1X,A)') solver_name, &
                               'Destroying face-to-cell weights done.'
    END IF ! global%verbLevel  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsf2c






! *******************************************************************************
!
! Purpose: Wrapper routine for destroying face-to-cell weights.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_destroywtsf2cwrapper(pRegion)

    IMPLICIT NONE

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

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

    TYPE(t_region), POINTER :: pregion

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

    TYPE(t_global), POINTER :: global  
    TYPE(t_mixt_input), POINTER :: pmixtinput   

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_DestroyWtsF2CWrapper',&
  'RFLU_ModWeights.F90')

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

    pmixtinput => pregion%mixtInput 

! ******************************************************************************
!   Call routines to compute weights
! ******************************************************************************

    SELECT CASE ( pmixtinput%stencilDimensFaces )
      CASE ( 1 )               
      CASE ( 2,3 ) 
        CALL rflu_destroywtsf2c(pregion)
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pMixtInput%stencilDimensFaces

! ******************************************************************************
!   End
! ******************************************************************************
    
    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_destroywtsf2cwrapper






! *******************************************************************************
!
! Purpose: Nullify weights for boundary face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_nullifywtsbf2c(pRegion,pPatch)

    IMPLICIT NONE

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

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

    TYPE(t_patch), POINTER :: ppatch   
    TYPE(t_region), POINTER :: pregion

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

    INTEGER :: ifl
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid      

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_NullifyWtsBF2C',&
  'RFLU_ModWeights.F90')

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    IF ( ppatch%bcType /= bc_virtual ) THEN 
      DO ifl = 1,ppatch%nBFaces
        nullify(ppatch%bf2cs(ifl)%xyzMoms)  
      END DO ! ifl
    END IF ! pPatch%bcType

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

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_nullifywtsbf2c  


  
  
  
  
! *******************************************************************************
!
! Purpose: Nullify weights for cell-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_nullifywtsc2c(pRegion)

    IMPLICIT NONE

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

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

    TYPE(t_region), POINTER :: pregion

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

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

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_NullifyWtsC2C',&
  'RFLU_ModWeights.F90')

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    DO icg = 1,pgrid%nCellsTot 
      nullify(pgrid%c2cs(icg)%xyzMoms)                   
    END DO ! icg

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

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_nullifywtsc2c  
  




! *******************************************************************************
!
! Purpose: Nullify weights for face-to-cell stencil.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
  
  SUBROUTINE rflu_nullifywtsf2c(pRegion)

    IMPLICIT NONE

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

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

    TYPE(t_region), POINTER :: pregion

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

    INTEGER :: ifg
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid          

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

    global => pregion%global

    CALL registerfunction(global,'RFLU_NullifyWtsF2C',&
  'RFLU_ModWeights.F90')

! ******************************************************************************
!   Set grid pointer
! ******************************************************************************

    pgrid => pregion%grid

! ******************************************************************************
!   Nullify memory
! ******************************************************************************

    DO ifg = 1,pgrid%nFaces 
      nullify(pgrid%f2cs(ifg)%xyzMoms)                   
    END DO ! ifg

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

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_nullifywtsf2c  

  
  
  



! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modweights


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModWeights.F90,v $
! Revision 1.16  2008/12/06 08:44:25  mtcampbe
! Updated license.
!
! Revision 1.15  2008/11/19 22:17:36  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.14  2007/02/27 13:08:27  haselbac
! Enabled 1d computations
!
! Revision 1.13  2006/04/07 16:03:33  haselbac
! Changed computation of bf2c wts to be done patch-wise
!
! Revision 1.12  2006/04/07 15:19:21  haselbac
! Removed tabs
!
! Revision 1.11  2006/04/07 14:52:40  haselbac
! Adapted to changes in stencilDimens params, bug fixes in 1D routine
!
! Revision 1.10  2006/03/09 20:50:56  haselbac
! Bug fix: Removed icg from INTEGER, INTENT(IN)
!
! Revision 1.9  2006/03/09 14:09:22  haselbac
! Wrapperified module bcos of 1D routines
!
! Revision 1.8  2006/01/06 22:14:45  haselbac
! Renamed routines bcos added 1D routines
!
! Revision 1.7  2005/10/05 14:16:51  haselbac
! Added routines for bfaces, clean-up
!
! Revision 1.6  2005/03/09 15:06:03  haselbac
! Added 2d option
!
! Revision 1.5  2004/07/06 15:14:48  haselbac
! Cosmetics only
!                                                
! Revision 1.4  2004/03/18 03:32:39  haselbac                                  
! Bug fix for i4 terms, joined init and computation                            
!
! 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  2004/01/13 16:22:32  haselbac                                  
! Changed to inverse-distance weighting                                        
!
! Revision 1.1  2003/12/04 03:28:43  haselbac                                  
! Initial revision                                                             
!
! ******************************************************************************