RFLU_ModRoeFlux.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: Collection of Roe flux routines.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModRoeFlux.F90,v 1.10 2008/12/06 08:44:24 mtcampbe Exp $
!
! Copyright: (c) 2005-2006 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modroeflux

  USE moddatatypes
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moderror
  USE moddatastruct, ONLY: t_region
  USE modbndpatch, ONLY: t_patch
  USE modgrid, ONLY: t_grid
  USE modmpi

  USE rflu_modentropyfixes

  USE rflu_modnscbc

  IMPLICIT NONE

  PRIVATE 
  PUBLIC :: rflu_roe_computeflux, &
            rflu_roe_computefluxc1, & 
            rflu_roe_computefluxc2_tcp, & 
            rflu_roe_computefluxd1_tcp, & 
            rflu_roe_computefluxd2_tcp, &         
            rflu_roe_computeflux1_gl, & 
            rflu_roe_computeflux1_tcp, &             
            rflu_roe_computeflux2_gl, &
            rflu_roe_computeflux2_tcp
  
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
      
  CHARACTER(CHRLEN) :: & 
    RCSIdentString = '$RCSfile: RFLU_ModRoeFlux.F90,v $ $Revision: 1.10 $' 
              
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS








! ******************************************************************************
!
! Purpose: Wrapper function for Roe flux functions.
!
! Description: None.
!
! Input: 
!   pRegion     Pointer to region data
!   fluxPart    Part of flux (central or dissipation)
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computeflux(pRegion,fluxPart)
                       
  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

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

! ==============================================================================
! Locals
! ==============================================================================
  
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFlux',&
  'RFLU_ModRoeFlux.F90')

! ******************************************************************************
! Call flux functions
! ******************************************************************************

  SELECT CASE ( global%flowType ) 

! ==============================================================================
!   Unsteady flow
! ==============================================================================

    CASE ( flow_unsteady )   
      SELECT CASE ( pregion%mixtInput%gasModel ) 

! ------------------------------------------------------------------------------    
!       Thermally and calorically perfect gas         
! ------------------------------------------------------------------------------    

        CASE ( gas_model_tcperf )
          SELECT CASE ( pregion%mixtInput%spaceOrder ) 
            CASE ( discr_order_1 ) 
              CALL rflu_roe_computeflux1_tcp(pregion)
            CASE ( discr_order_2 )
              CALL rflu_roe_computeflux2_tcp(pregion)   
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! pRegion%mixtInput%spaceOrder          

! ------------------------------------------------------------------------------    
!       Mixture of gas, vapor, and liquid     
! ------------------------------------------------------------------------------    

        CASE ( gas_model_mixt_gasliq )
          SELECT CASE ( pregion%mixtInput%spaceOrder ) 
            CASE ( discr_order_1 ) 
              CALL rflu_roe_computeflux1_gl(pregion)
            CASE ( discr_order_2 )
              CALL rflu_roe_computeflux2_gl(pregion)               
            CASE default
              CALL errorstop(global,err_reached_default,__line__)
          END SELECT ! pRegion%mixtInput%spaceOrder           

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

        CASE default
          CALL errorstop(global,err_reached_default,__line__)             
      END SELECT ! pRegion%mixtInput%gasModel  

! ==============================================================================
!   Steady flow. NOTE for explicit solver, distinguish between central and
!   dissipative parts of flux.
! ==============================================================================

    CASE ( flow_steady )
      SELECT CASE ( global%solverType ) 

! ------------------------------------------------------------------------------    
!       Explicit solver
! ------------------------------------------------------------------------------    

        CASE ( solv_explicit ) 
          SELECT CASE ( fluxpart ) 

! --------- Central part of flux -----------------------------------------------

            CASE ( flux_part_central ) 
              SELECT CASE ( pregion%mixtInput%gasModel ) 
                CASE ( gas_model_tcperf )
                  SELECT CASE ( pregion%mixtInput%spaceOrder ) 
                    CASE ( discr_order_1 ) 
                      CALL rflu_roe_computefluxc1(pregion)
                    CASE ( discr_order_2 )
                      CALL rflu_roe_computefluxc2_tcp(pregion)              
                    CASE default
                      CALL errorstop(global,err_reached_default,__line__)
                  END SELECT ! pRegion%mixtInput%spaceOrder          
                CASE default
                  CALL errorstop(global,err_reached_default,__line__)             
              END SELECT ! pRegion%mixtInput%gasModel              

! --------- Dissipative part of flux -------------------------------------------

            CASE ( flux_part_dissip ) 
              SELECT CASE ( pregion%mixtInput%gasModel ) 
                CASE ( gas_model_tcperf )
                  SELECT CASE ( pregion%mixtInput%spaceOrder ) 
                    CASE ( discr_order_1 ) 
                      CALL rflu_roe_computefluxd1_tcp(pregion)
                    CASE ( discr_order_2 )
                      CALL rflu_roe_computefluxd2_tcp(pregion)           
                    CASE default
                      CALL errorstop(global,err_reached_default,__line__)
                  END SELECT ! pRegion%mixtInput%spaceOrder          
                CASE default
                  CALL errorstop(global,err_reached_default,__line__)             
              END SELECT ! pRegion%mixtInput%gasModel            

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

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

! ------------------------------------------------------------------------------    
!       Implicit solver. NOTE compute entire flux.             
! ------------------------------------------------------------------------------    

        CASE ( solv_implicit_nk ) 
          SELECT CASE ( pregion%mixtInput%gasModel ) 
            CASE ( gas_model_tcperf )
              SELECT CASE ( pregion%mixtInput%spaceOrder ) 
                CASE ( discr_order_1 ) 
                  CALL rflu_roe_computeflux1_tcp(pregion)
                CASE ( discr_order_2 )
                  CALL rflu_roe_computeflux2_tcp(pregion)               
                CASE default
                  CALL errorstop(global,err_reached_default,__line__)
              END SELECT ! pRegion%mixtInput%spaceOrder          
            CASE default
              CALL errorstop(global,err_reached_default,__line__)             
          END SELECT ! pRegion%mixtInput%gasModel              

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

        CASE default
          CALL errorstop(global,err_reached_default,__line__)   
      END SELECT ! global%solverType                               
    CASE default
      CALL errorstop(global,err_reached_default,__line__) 
  END SELECT ! global%flowType      

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computeflux








! ******************************************************************************
!
! Purpose: Compute central convective fluxes using first-order accurate
!   accurate approximation by average of variables.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computefluxc1(pRegion)

  USE modinterfaces, ONLY: rflu_centralfirstpatch

  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,ipatch
  REAL(RFREAL) :: el,er,fs,irl,irr,nm,nx,ny,nz,ql,qr,rl,rr,ul,ur,vl,vr,wl,wr, &
                  pl,pr
  REAL(RFREAL) :: flx(5)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,prhs
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFluxC1',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFluxC1")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs = pgrid%indGs

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  
  prhs => pregion%mixt%rhs

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl

    ul  = pcv(cv_mixt_xmom,c1)*irl
    vl  = pcv(cv_mixt_ymom,c1)*irl
    wl  = pcv(cv_mixt_zmom,c1)*irl
    el  = pcv(cv_mixt_ener,c1)*irl
    pl  = pdv(dv_mixt_pres,c1)

    ql  = ul*nx + vl*ny + wl*nz - fs

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr

    ur  = pcv(cv_mixt_xmom,c2)*irr
    vr  = pcv(cv_mixt_ymom,c2)*irr
    wr  = pcv(cv_mixt_zmom,c2)*irr
    er  = pcv(cv_mixt_ener,c2)*irr
    pr  = pdv(dv_mixt_pres,c2)

    qr  = ur*nx + vr*ny + wr*nz - fs

! ==============================================================================
!   Compute fluxes 
! ==============================================================================

    flx(1) = 0.5_rfreal*(ql* rl                  + qr* rr                 )*nm
    flx(2) = 0.5_rfreal*(ql* rl*ul       + pl*nx + qr* rr*ur       + pr*nx)*nm
    flx(3) = 0.5_rfreal*(ql* rl*vl       + pl*ny + qr* rr*vr       + pr*ny)*nm
    flx(4) = 0.5_rfreal*(ql* rl*wl       + pl*nz + qr* rr*wr       + pr*nz)*nm
    flx(5) = 0.5_rfreal*(ql*(rl*el + pl) + pl*fs + qr*(rr*er + pr) + pr*fs)*nm

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      CALL rflu_nscbc_compfirstpatchflux(pregion,ppatch)
    ELSE
      CALL rflu_centralfirstpatch(pregion,ppatch)
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFluxC1")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computefluxc1










! ******************************************************************************
!
! Purpose: Compute central convective fluxes using second-order accurate
!   accurate approximation by average of variables for thermally and 
!   calorically perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes:
!   1. Assumes ideal gas and constant Cp and R! This assumption is necessary
!      because need to compute face states from reconstructed variables.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computefluxc2_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_eo_dgpuvw, & 
                           rflu_centralfirstpatch, &
                           rflu_centralsecondpatch

  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,ipatch
  REAL(RFREAL) :: dx,dy,dz,el,er,fs,g,irl,irr,nm,nx,ny,nz,ql,qr,rl,rr,ul,ur, &
                  vl,vr,wl,wr,pl,pr,xc,yc,zc
  REAL(RFREAL) :: flx(5)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFluxC2_TCP',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFluxC2_TCP")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs = pgrid%indGs

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  
  pgc  => pregion%mixt%gradCell    
  prhs => pregion%mixt%rhs  

! ******************************************************************************
! Define constants
! ******************************************************************************

  g = global%refGamma

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    xc = pgrid%fc(xcoord,ifg)
    yc = pgrid%fc(ycoord,ifg)
    zc = pgrid%fc(zcoord,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl

    ul  = pcv(cv_mixt_xmom,c1)*irl
    vl  = pcv(cv_mixt_ymom,c1)*irl
    wl  = pcv(cv_mixt_zmom,c1)*irl
    pl  = pdv(dv_mixt_pres,c1)

    dx  = xc - pgrid%cofg(xcoord,c1)
    dy  = yc - pgrid%cofg(ycoord,c1)
    dz  = zc - pgrid%cofg(zcoord,c1)

    rl = rl + pgc(xcoord,grc_mixt_dens,c1)*dx &
            + pgc(ycoord,grc_mixt_dens,c1)*dy &
            + pgc(zcoord,grc_mixt_dens,c1)*dz
    ul = ul + pgc(xcoord,grc_mixt_xvel,c1)*dx &
            + pgc(ycoord,grc_mixt_xvel,c1)*dy &
            + pgc(zcoord,grc_mixt_xvel,c1)*dz
    vl = vl + pgc(xcoord,grc_mixt_yvel,c1)*dx &
            + pgc(ycoord,grc_mixt_yvel,c1)*dy &
            + pgc(zcoord,grc_mixt_yvel,c1)*dz
    wl = wl + pgc(xcoord,grc_mixt_zvel,c1)*dx &
            + pgc(ycoord,grc_mixt_zvel,c1)*dy &
            + pgc(zcoord,grc_mixt_zvel,c1)*dz
    pl = pl + pgc(xcoord,grc_mixt_pres,c1)*dx &
            + pgc(ycoord,grc_mixt_pres,c1)*dy &
            + pgc(zcoord,grc_mixt_pres,c1)*dz

    el = mixtperf_eo_dgpuvw(rl,g,pl,ul,vl,wl)
    ql = ul*nx + vl*ny + wl*nz - fs

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr

    ur  = pcv(cv_mixt_xmom,c2)*irr
    vr  = pcv(cv_mixt_ymom,c2)*irr
    wr  = pcv(cv_mixt_zmom,c2)*irr
    pr  = pdv(dv_mixt_pres,c2)

    dx  = xc - pgrid%cofg(xcoord,c2)
    dy  = yc - pgrid%cofg(ycoord,c2)
    dz  = zc - pgrid%cofg(zcoord,c2)

    rr = rr + pgc(xcoord,grc_mixt_dens,c2)*dx &
            + pgc(ycoord,grc_mixt_dens,c2)*dy &
            + pgc(zcoord,grc_mixt_dens,c2)*dz
    ur = ur + pgc(xcoord,grc_mixt_xvel,c2)*dx &
            + pgc(ycoord,grc_mixt_xvel,c2)*dy &
            + pgc(zcoord,grc_mixt_xvel,c2)*dz
    vr = vr + pgc(xcoord,grc_mixt_yvel,c2)*dx &
            + pgc(ycoord,grc_mixt_yvel,c2)*dy &
            + pgc(zcoord,grc_mixt_yvel,c2)*dz
    wr = wr + pgc(xcoord,grc_mixt_zvel,c2)*dx &
            + pgc(ycoord,grc_mixt_zvel,c2)*dy &
            + pgc(zcoord,grc_mixt_zvel,c2)*dz
    pr = pr + pgc(xcoord,grc_mixt_pres,c2)*dx &
            + pgc(ycoord,grc_mixt_pres,c2)*dy &
            + pgc(zcoord,grc_mixt_pres,c2)*dz

    er = mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)
    qr = ur*nx + vr*ny + wr*nz - fs

! ==============================================================================
!   Compute fluxes 
! ==============================================================================

    flx(1) = 0.5_rfreal*(ql* rl                  + qr* rr                 )*nm
    flx(2) = 0.5_rfreal*(ql* rl*ul       + pl*nx + qr* rr*ur       + pr*nx)*nm
    flx(3) = 0.5_rfreal*(ql* rl*vl       + pl*ny + qr* rr*vr       + pr*ny)*nm
    flx(4) = 0.5_rfreal*(ql* rl*wl       + pl*nz + qr* rr*wr       + pr*nz)*nm
    flx(5) = 0.5_rfreal*(ql*(rl*el + pl) + pl*fs + qr*(rr*er + pr) + pr*fs)*nm

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_nscbc_compfirstpatchflux(pregion,ppatch)
        CASE ( 2 ) 
          CALL rflu_nscbc_compsecondpatchflux(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    ELSE
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_centralfirstpatch(pregion,ppatch)      
        CASE ( 2 ) 
          CALL rflu_centralsecondpatch(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFluxC2_TCP")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computefluxc2_tcp







! ******************************************************************************
!
! Purpose: Compute dissipative fluxes using first-order accurate Roe scheme
!   for thermally and calorically perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computefluxd1_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_c_ghovm2, &
                           mixtperf_ho_cptuvw

  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs
  REAL(RFREAL) :: ah,betrk,cp,dissfact,dp,dq,dr,du,de,dv1,dv2,dv5,dw,efc, &
                  epsentr,fs,g,hh,hl,hr,irl,irr,l1,l2,l5,nm,nx,ny,nz,pl, &
                  pr,qh,ql,qr,rl,rh,rr,sh,term,tl,tr,t1,t2,t3,t5,ul,uh,ur, &
                  vl,vh,vr,wl,wh,wr,wt
  REAL(RFREAL) :: flx(5)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdiss,pdv
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid    

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFluxD1_TCP',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFluxD1_TCP")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs = pgrid%indGs

  pcv   => pregion%mixt%cv
  pdv   => pregion%mixt%dv
  
  pdiss => pregion%mixt%diss

! ******************************************************************************
! Define constants
! ******************************************************************************

  cp = global%refCp
  g  = global%refGamma

  epsentr  = pregion%mixtInput%epsentr
  betrk    = pregion%mixtInput%betrk(pregion%irkStep)
  dissfact = pregion%mixtInput%dissFact
  term     = 0.5_rfreal*betrk*dissfact

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl
    ul  = pcv(cv_mixt_xmom,c1)*irl
    vl  = pcv(cv_mixt_ymom,c1)*irl
    wl  = pcv(cv_mixt_zmom,c1)*irl
    pl  = pdv(dv_mixt_pres,c1)
    tl  = pdv(dv_mixt_temp,c1)

    hl = mixtperf_ho_cptuvw(cp,tl,ul,vl,wl)

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    ur  = pcv(cv_mixt_xmom,c2)*irr
    vr  = pcv(cv_mixt_ymom,c2)*irr
    wr  = pcv(cv_mixt_zmom,c2)*irr
    pr  = pdv(dv_mixt_pres,c2)
    tr  = pdv(dv_mixt_temp,c2)

    hr = mixtperf_ho_cptuvw(cp,tr,ur,vr,wr)

! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)
    wt = rl/(rl + rh)

    uh = wt*ul + (1.0_rfreal-wt)*ur
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    hh = wt*hl + (1.0_rfreal-wt)*hr

    qh = uh*nx + vh*ny + wh*nz - fs
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    ah = mixtperf_c_ghovm2(g,hh,2.0_rfreal*sh) ! NOTE factor of 2

! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh - ah)
    l2 = abs(qh     )
    l5 = abs(qh + ah)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr = rr - rl
    du = ur - ul
    de = vr - vl ! use de to denote dv because dv is already used
    dw = wr - wl
    dp = pr - pl
    dq = du*nx + de*ny + dw*nz ! Note: no gs contribution

    dv1 = (dp - rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv5 = (dp + rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv2 = dr - dp/(ah*ah)

    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l2*rh
    t5 = l5*dv5

! ==============================================================================
!   Compute fluxes
! ==============================================================================

    flx(1) = term*(t1            + t2                                   + t5           )*nm
    flx(2) = term*(t1*(uh-nx*ah) + t2*uh + t3*(du-nx*dq)                + t5*(uh+nx*ah))*nm
    flx(3) = term*(t1*(vh-ny*ah) + t2*vh + t3*(de-ny*dq)                + t5*(vh+ny*ah))*nm
    flx(4) = term*(t1*(wh-nz*ah) + t2*wh + t3*(dw-nz*dq)                + t5*(wh+nz*ah))*nm
    flx(5) = term*(t1*(hh-qh*ah) + t2*sh + t3*(uh*du+vh*de+wh*dw-qh*dq) + t5*(hh+qh*ah))*nm

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    pdiss(cv_mixt_dens,c1) = pdiss(cv_mixt_dens,c1) + flx(1)
    pdiss(cv_mixt_xmom,c1) = pdiss(cv_mixt_xmom,c1) + flx(2)
    pdiss(cv_mixt_ymom,c1) = pdiss(cv_mixt_ymom,c1) + flx(3)
    pdiss(cv_mixt_zmom,c1) = pdiss(cv_mixt_zmom,c1) + flx(4)
    pdiss(cv_mixt_ener,c1) = pdiss(cv_mixt_ener,c1) + flx(5)

    pdiss(cv_mixt_dens,c2) = pdiss(cv_mixt_dens,c2) - flx(1)
    pdiss(cv_mixt_xmom,c2) = pdiss(cv_mixt_xmom,c2) - flx(2)
    pdiss(cv_mixt_ymom,c2) = pdiss(cv_mixt_ymom,c2) - flx(3)
    pdiss(cv_mixt_zmom,c2) = pdiss(cv_mixt_zmom,c2) - flx(4)
    pdiss(cv_mixt_ener,c2) = pdiss(cv_mixt_ener,c2) - flx(5)
  END DO ! ifg

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFluxD1_TCP")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computefluxd1_tcp








! ******************************************************************************
!
! Purpose: Compute dissipative fluxes using second-order accurate Roe scheme
!   for thermally and calorically perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes:
!   1. Assumes ideal gas and constant Cp and R!
!
! ******************************************************************************

SUBROUTINE rflu_roe_computefluxd2_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_c_ghovm2, & 
                           mixtperf_ho_cptuvw, &
                           mixtperf_r_cpg, &
                           mixtperf_t_dpr

  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs
  REAL(RFREAL) :: ah,betrk,cp,de,dissfact,dp,dq,dr,du,dv1,dv2,dv5,dw,dx, &
                  dy,dz,efc,epsentr,fs,g,gc,hh,hl,hr,irl,irr,l1,l2,l5,nm,nx, &
                  ny,nz,pl,pr,qh,ql,qr,rl,rh,rr,sh,term,tl,tr,t1,t2,t3, &
                  t5,ul,uh,ur,vl,vh,vr,wl,wh,wr,wt,xc,yc,zc
  REAL(RFREAL) :: flx(5)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdiss,pdv
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid    

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFluxD2_TCP',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFluxD2_TCP")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs = pgrid%indGs

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  
  pgc   => pregion%mixt%gradCell  
  pdiss => pregion%mixt%diss
  
! ******************************************************************************
! Define constants
! ******************************************************************************

  cp = global%refCp
  g  = global%refGamma
  gc = mixtperf_r_cpg(cp,g)

  epsentr  = pregion%mixtInput%epsentr
  betrk    = pregion%mixtInput%betrk(pregion%irkStep)
  dissfact = pregion%mixtInput%dissFact
  term     = 0.5_rfreal*betrk*dissfact

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    xc = pgrid%fc(xcoord,ifg)
    yc = pgrid%fc(ycoord,ifg)
    zc = pgrid%fc(zcoord,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl

    ul  = pcv(cv_mixt_xmom,c1)*irl
    vl  = pcv(cv_mixt_ymom,c1)*irl
    wl  = pcv(cv_mixt_zmom,c1)*irl
    pl  = pdv(dv_mixt_pres,c1)

    dx  = xc - pgrid%cofg(xcoord,c1)
    dy  = yc - pgrid%cofg(ycoord,c1)
    dz  = zc - pgrid%cofg(zcoord,c1)

    rl = rl + pgc(xcoord,grc_mixt_dens,c1)*dx &
            + pgc(ycoord,grc_mixt_dens,c1)*dy &
            + pgc(zcoord,grc_mixt_dens,c1)*dz
    ul = ul + pgc(xcoord,grc_mixt_xvel,c1)*dx &
            + pgc(ycoord,grc_mixt_xvel,c1)*dy &
            + pgc(zcoord,grc_mixt_xvel,c1)*dz
    vl = vl + pgc(xcoord,grc_mixt_yvel,c1)*dx &
            + pgc(ycoord,grc_mixt_yvel,c1)*dy &
            + pgc(zcoord,grc_mixt_yvel,c1)*dz
    wl = wl + pgc(xcoord,grc_mixt_zvel,c1)*dx &
            + pgc(ycoord,grc_mixt_zvel,c1)*dy &
            + pgc(zcoord,grc_mixt_zvel,c1)*dz
    pl = pl + pgc(xcoord,grc_mixt_pres,c1)*dx &
            + pgc(ycoord,grc_mixt_pres,c1)*dy &
            + pgc(zcoord,grc_mixt_pres,c1)*dz

    tl = mixtperf_t_dpr(rl,pl,gc)
    hl = mixtperf_ho_cptuvw(cp,tl,ul,vl,wl)

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    ur  = pcv(cv_mixt_xmom,c2)*irr
    vr  = pcv(cv_mixt_ymom,c2)*irr
    wr  = pcv(cv_mixt_zmom,c2)*irr
    pr  = pdv(dv_mixt_pres,c2)

    dx  = xc - pgrid%cofg(xcoord,c2)
    dy  = yc - pgrid%cofg(ycoord,c2)
    dz  = zc - pgrid%cofg(zcoord,c2)

    rr = rr + pgc(xcoord,grc_mixt_dens,c2)*dx &
            + pgc(ycoord,grc_mixt_dens,c2)*dy &
            + pgc(zcoord,grc_mixt_dens,c2)*dz
    ur = ur + pgc(xcoord,grc_mixt_xvel,c2)*dx &
            + pgc(ycoord,grc_mixt_xvel,c2)*dy &
            + pgc(zcoord,grc_mixt_xvel,c2)*dz
    vr = vr + pgc(xcoord,grc_mixt_yvel,c2)*dx &
            + pgc(ycoord,grc_mixt_yvel,c2)*dy &
            + pgc(zcoord,grc_mixt_yvel,c2)*dz
    wr = wr + pgc(xcoord,grc_mixt_zvel,c2)*dx &
            + pgc(ycoord,grc_mixt_zvel,c2)*dy &
            + pgc(zcoord,grc_mixt_zvel,c2)*dz
    pr = pr + pgc(xcoord,grc_mixt_pres,c2)*dx &
            + pgc(ycoord,grc_mixt_pres,c2)*dy &
            + pgc(zcoord,grc_mixt_pres,c2)*dz

    tr = mixtperf_t_dpr(rr,pr,gc)
    hr = mixtperf_ho_cptuvw(cp,tr,ur,vr,wr)

! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)
    wt = rl/(rl + rh)

    uh = wt*ul + (1.0_rfreal-wt)*ur
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    hh = wt*hl + (1.0_rfreal-wt)*hr

    qh = uh*nx + vh*ny + wh*nz - fs
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    ah = mixtperf_c_ghovm2(g,hh,2.0_rfreal*sh) ! NOTE factor of 2

! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh - ah)
    l2 = abs(qh     )
    l5 = abs(qh + ah)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr = rr - rl
    du = ur - ul
    de = vr - vl ! use de to denote dv because dv is already used
    dw = wr - wl
    dp = pr - pl
    dq = du*nx + de*ny + dw*nz ! Note: no gs contribution

    dv1 = (dp - rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv5 = (dp + rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv2 = dr - dp/(ah*ah)

    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l2*rh
    t5 = l5*dv5

! ==============================================================================
!   Compute fluxes
! ==============================================================================

    flx(1) = term*(t1            + t2                                   + t5           )*nm
    flx(2) = term*(t1*(uh-nx*ah) + t2*uh + t3*(du-nx*dq)                + t5*(uh+nx*ah))*nm
    flx(3) = term*(t1*(vh-ny*ah) + t2*vh + t3*(de-ny*dq)                + t5*(vh+ny*ah))*nm
    flx(4) = term*(t1*(wh-nz*ah) + t2*wh + t3*(dw-nz*dq)                + t5*(wh+nz*ah))*nm
    flx(5) = term*(t1*(hh-qh*ah) + t2*sh + t3*(uh*du+vh*de+wh*dw-qh*dq) + t5*(hh+qh*ah))*nm

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    pdiss(cv_mixt_dens,c1) = pdiss(cv_mixt_dens,c1) + flx(1)
    pdiss(cv_mixt_xmom,c1) = pdiss(cv_mixt_xmom,c1) + flx(2)
    pdiss(cv_mixt_ymom,c1) = pdiss(cv_mixt_ymom,c1) + flx(3)
    pdiss(cv_mixt_zmom,c1) = pdiss(cv_mixt_zmom,c1) + flx(4)
    pdiss(cv_mixt_ener,c1) = pdiss(cv_mixt_ener,c1) + flx(5)

    pdiss(cv_mixt_dens,c2) = pdiss(cv_mixt_dens,c2) - flx(1)
    pdiss(cv_mixt_xmom,c2) = pdiss(cv_mixt_xmom,c2) - flx(2)
    pdiss(cv_mixt_ymom,c2) = pdiss(cv_mixt_ymom,c2) - flx(3)
    pdiss(cv_mixt_zmom,c2) = pdiss(cv_mixt_zmom,c2) - flx(4)
    pdiss(cv_mixt_ener,c2) = pdiss(cv_mixt_ener,c2) - flx(5)
  END DO ! ifg

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFluxD2_TCP")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computefluxd2_tcp










! ******************************************************************************
!
! Purpose: Compute convective fluxes using first-order accurate Roe scheme for
!   multiphase flow.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!  
! ******************************************************************************

SUBROUTINE rflu_roe_computeflux1_gl(pRegion)

  USE modinterfaces, ONLY: mixtgasliq_c, &
                           mixtliq_c2_bp, &
                           mixtliq_d_dobpppobttto, &
                           mixtperf_c2_grt, &
                           mixtperf_cv_cpr, & 
                           mixtperf_d_prt, &
                           mixtperf_r_cpg, &
                           mixtperf_r_m, &
                           rflu_centralfirstpatch_gl
  
  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,indmf,ipatch
  REAL(RFREAL) :: bgh2,blh2,bp,bt,bvh2,cmh,cml,cmr,cvg,cvl,cvv,cgh2,clh2, &
                  cvh2,cvm,dissfact,de,dp,dq,dr,dtemp,du,dv1,dv2,dv3,dv4, &
                  dv5,dvfg,dvfv,dw,efc,egh,el,elh,epsentr,er,evh,factor1, &
                  factor2,factor4,factor5,fdv1,fdv2,fdv3,fs,fxn1,fxn2, &
                  fxn3,fxn4,fxn5,fxn6,invc2p,invc2pb,irl,irr,l1,l2,l3,l4, &  
                  l5,nm,nx,ny,nz,ph,pl,po,pr,qh,ql,qr,ra11,ra12x,ra12y, &
                  ra12z,ra13,ra14,ra15,ra21,ra22x,ra22y,ra22z,ra23,ra24, &
                  ra25,ra31,ra32x,ra32y,ra32z,ra33,ra34,ra35,ra41,ra42x, &
                  ra42y,ra42z,ra43,ra44,ra45,ra51,ra52x,ra52y,ra52z,ra53, &
                  ra54,ra55,rgh,rgl,rgr,rh,rl,rlh,ro,rr,rvh,rvl,rvr, &
                  rygl,rygr,ryvl,ryvr,rg,rv,sh,t1,t2,t3,t4,t5,t6,th,tl, &
                  to,tr,uh,ul,ur,vfgh,vfgl,vfgr,vflh,vfvh,vfvl,vfvr,vh, & 
                  vl,vr,wh,wl,wr,wt              
  REAL(RFREAL) :: flxconv(5),flxdiss(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: fn,pcvmixt,pcvspec,pdvmixt,prhs
  TYPE(t_global), POINTER :: global
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_grid), POINTER :: pgrid   
 
! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFlux1_GL',&
  'RFLU_ModRoeFlux.F90')

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************
 
  pgrid => pregion%grid 
 
  indgs    = pregion%grid%indGs
  indmf    = pregion%mixtInput%indMfMixt
  epsentr  = pregion%mixtInput%epsentr
  dissfact = pregion%mixtInput%dissFact

  pcvmixt => pregion%mixt%cv
  pdvmixt => pregion%mixt%dv
  pcvspec => pregion%spec%cv
  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs

! ******************************************************************************
! Define constants
! ******************************************************************************

  ro  = global%refDensityLiq
  po  = global%refPressLiq
  to  = global%refTempLiq
  bp  = global%refBetaPLiq
  bt  = global%refBetaTLiq
  cvl = global%refCvLiq

  rg  = mixtperf_r_m(pregion%specInput%specType(1)%pMaterial%molw)
  cvg = mixtperf_cv_cpr(pregion%specInput%specType(1)%pMaterial%spht,rg)

  rv  = mixtperf_r_m(pregion%specInput%specType(2)%pMaterial%molw)
  cvv = mixtperf_cv_cpr(pregion%specInput%specType(2)%pMaterial%spht,rv)

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------
    
    rl  = pcvmixt(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl
    
    ul  = pcvmixt(cv_mixt_xmom,c1)*irl
    vl  = pcvmixt(cv_mixt_ymom,c1)*irl
    wl  = pcvmixt(cv_mixt_zmom,c1)*irl
    el  = pcvmixt(cv_mixt_ener,c1)*irl
    pl  = pdvmixt(dv_mixt_pres,c1)
    tl  = pdvmixt(dv_mixt_temp,c1)
    cml = pdvmixt(dv_mixt_soun,c1)
    
    rygl = pcvspec(1,c1)
    ryvl = pcvspec(2,c1)    
    rvl  = mixtperf_d_prt(pl,rv,tl)
    rgl  = mixtperf_d_prt(pl,rg,tl)    
    vfvl = ryvl/rvl
    vfgl = rygl/rgl
 
    ql = ul*nx + vl*ny + wl*nz - fs

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------
   
    rr  = pcvmixt(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    
    ur  = pcvmixt(cv_mixt_xmom,c2)*irr
    vr  = pcvmixt(cv_mixt_ymom,c2)*irr
    wr  = pcvmixt(cv_mixt_zmom,c2)*irr
    er  = pcvmixt(cv_mixt_ener,c2)*irr
    pr  = pdvmixt(dv_mixt_pres,c2)
    tr  = pdvmixt(dv_mixt_temp,c2)
    cmr = pdvmixt(dv_mixt_soun,c2)
 
    rygr = pcvspec(1,c2)
    ryvr = pcvspec(2,c2)   
    rvr  = mixtperf_d_prt(pr,rv,tr)
    rgr  = mixtperf_d_prt(pr,rg,tr)            
    vfvr = ryvr/rvr
    vfgr = rygr/rgr
 
    qr = ur*nx + vr*ny + wr*nz - fs
   
! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)                     
    wt = rl/(rl + rh)                     
    uh = wt*ul + (1.0_rfreal-wt)*ur     
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    ph = wt*pl + (1.0_rfreal-wt)*pr      
    th = wt*tl + (1.0_rfreal-wt)*tr       

    vfgh = wt*vfgl + (1.0_rfreal-wt)*vfgr
    vfvh = wt*vfvl + (1.0_rfreal-wt)*vfvr
    vflh = 1.0_rfreal - vfvh - vfgh

    qh = uh*nx + vh*ny + wh*nz - fs 
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    elh = cvl*th + sh
    evh = cvv*th + sh
    egh = cvg*th + sh

    clh2 = mixtliq_c2_bp(bp) 
    cvh2 = mixtperf_c2_grt(1.0_rfreal,rv,th)
    cgh2 = mixtperf_c2_grt(1.0_rfreal,rg,th)

    rlh = mixtliq_d_dobpppobttto(ro,bp,bt,ph,po,th,to)
    rvh = mixtperf_d_prt(ph,rv,th)
    rgh = mixtperf_d_prt(ph,rg,th)

    blh2 = -bt/bp
    bvh2 = rvh*rv
    bgh2 = rgh*rg

    cvm = (rlh*vflh*cvl + rvh*vfvh*cvv + rgh*vfgh*cvg)/rh
 
    cmh = mixtgasliq_c(cvm,rh,ph,rlh,rvh,rgh,vflh,vfvh,vfgh,clh2,cvh2,cgh2, &
                       blh2,bvh2,bgh2)
    
! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh      )
    l2 = abs(qh      )
    l3 = abs(qh      )
    l4 = abs(qh - cmh)
    l5 = abs(qh + cmh)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5 
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l3  = entropyfixhartenhyman(l3,efc)
    l4  = entropyfixhartenhyman(l4,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr    = rr - rl     
    du    = ur - ul
    de    = vr - vl
    dw    = wr - wl  
    dq    = du*nx + de*ny + dw*nz ! Note: no gs contribution
    dp    = pr - pl 
    dtemp = tr - tl
    dvfv  = vfvr - vfvl
    dvfg  = vfgr - vfgl

    fdv1 = ph/(rh*cmh*cmh*rh*cvm*th)
    fdv2 = (vfgh*(rgh*cgh2 - rh*cmh*cmh  &
         + (bgh2*ph)/(rh*cvm)))/(rgh*cgh2*rh*cmh*cmh)
    fdv3 = (vfvh*(rvh*cvh2 - rh*cmh*cmh  &
         + (bvh2*ph)/(rh*cvm)))/(rvh*cvh2*rh*cmh*cmh)

    dv1  = (dtemp/th) - (dp*fdv1)
    dv2  = dvfg - dp*fdv2
    dv3  = dvfv - dp*fdv3
    dv4  = -dq/cmh + dp/(rh*cmh*cmh) 
    dv5  =  dq/cmh + dp/(rh*cmh*cmh) 

    invc2p  = vflh/clh2 + vfgh/cgh2 + vfvh/cvh2 
    invc2pb = (blh2*vflh)/clh2 + (bgh2*vfgh)/cgh2 & 
            + (bvh2*vfvh)/cvh2

    ra11  = -th*invc2pb
    ra12x = -th*invc2pb*uh
    ra12y = -th*invc2pb*vh
    ra12z = -th*invc2pb*wh
    ra13  = rlh*vflh*cvl*th + bt*cvl*th*th*vflh - th*invc2pb*sh
    ra14  = -(bgh2*vfgh*th)/cgh2
    ra15  = -(bvh2*vfvh*th)/cvh2

    ra21  =  rgh - rlh
    ra22x = (rgh - rlh)*uh
    ra22y = (rgh - rlh)*vh
    ra22z = (rgh - rlh)*wh
    ra23  =  rgh*egh - rlh*elh
    ra24  =  rgh
    ra25  =  0.0_rfreal    

    ra31  =  rvh - rlh
    ra32x = (rvh - rlh)*uh
    ra32y = (rvh - rlh)*vh
    ra32z = (rvh - rlh)*wh
    ra33  =  rvh*evh - rlh*elh
    ra34  =  0.0_rfreal
    ra35  =  rvh    
    
    fxn1  = ((rgh - rlh)*vfgh*(rgh*cgh2 - rh*cmh*cmh &
          + (bgh2*ph)/(rh*cvm)))/(2.0_rfreal*rgh*cgh2)
    fxn2  = ((rvh - rlh)*vfvh*(rvh*cvh2 - rh*cmh*cmh &
          + (bvh2*ph)/(rh*cvm)))/(2.0_rfreal*rvh*cvh2)
    factor1 = 0.5_rfreal*rh*cmh*cmh*invc2p & 
            - (ph*invc2pb)/(2.0_rfreal*rh*cvm) + fxn1 + fxn2

    fxn3  = 0.5_rfreal*rh*cmh*cmh*((elh*vflh)/clh2 &
          + (egh*vfgh)/cgh2 + (evh*vfvh)/cvh2)
    fxn4  = (ph*(rlh*vflh*cvl + bt*cvl*th*vflh &
          -   (sh*invc2pb)))/(2.0_rfreal*rh*cvm)
    fxn5  = ((rgh*egh - rlh*elh)*vfgh*(rgh*cgh2 - rh*cmh*cmh &
          + (bgh2*ph)/(rh*cvm)))/(2.0_rfreal*rgh*cgh2)
    fxn6  = ((rvh*evh - rlh*elh)*vfvh*(rvh*cvh2 - rh*cmh*cmh &
          + (bvh2*ph)/(rh*cvm)))/(2.0_rfreal*rvh*cvh2)
    factor2 = fxn3 + fxn4 + fxn5 + fxn6

    factor4 = 0.5_rfreal*rgh*vfgh
    factor5 = 0.5_rfreal*rvh*vfvh

    ra41  = factor1
    ra42x = factor1*uh - 0.5_rfreal*rh*cmh*nx 
    ra42y = factor1*vh - 0.5_rfreal*rh*cmh*ny
    ra42z = factor1*wh - 0.5_rfreal*rh*cmh*nz
    ra43  = factor2    - 0.5_rfreal*rh*cmh*qh
    ra44  = factor4
    ra45  = factor5

    ra51  = factor1
    ra52x = factor1*uh + 0.5_rfreal*rh*cmh*nx 
    ra52y = factor1*vh + 0.5_rfreal*rh*cmh*ny
    ra52z = factor1*wh + 0.5_rfreal*rh*cmh*nz
    ra53  = factor2    + 0.5_rfreal*rh*cmh*qh
    ra54  = factor4
    ra55  = factor5   
 
    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l3*dv3
    t4 = l4*dv4
    t5 = l5*dv5
    t6 = l1*rh 

! ==============================================================================
!   Compute fluxes
! ==============================================================================

! ------------------------------------------------------------------------------
!   Central part
! ------------------------------------------------------------------------------
  
    flxconv(1) = 0.5_rfreal*(ql*rl            + qr*rr           )*nm
    flxconv(2) = 0.5_rfreal*(ql*rl*ul + pl*nx + qr*rr*ur + pr*nx)*nm
    flxconv(3) = 0.5_rfreal*(ql*rl*vl + pl*ny + qr*rr*vr + pr*ny)*nm
    flxconv(4) = 0.5_rfreal*(ql*rl*wl + pl*nz + qr*rr*wr + pr*nz)*nm
    flxconv(5) = 0.5_rfreal*(ql*(rl*el + pl) + pl*fs  + qr*(rr*er + pr) & 
                           + pr*fs)*nm

! ------------------------------------------------------------------------------
!   Dissipative part
! ------------------------------------------------------------------------------

    flxdiss(1) = 0.5_rfreal*dissfact*(ra11*t1  + ra21*t2  + ra31*t3 &
                                    + ra41*t4  + ra51*t5)*nm 
    flxdiss(2) = 0.5_rfreal*dissfact*(ra12x*t1 + ra22x*t2 + ra32x*t3 &
                                    + ra42x*t4 + ra52x*t5 + t6*(du-nx*dq))*nm
    flxdiss(3) = 0.5_rfreal*dissfact*(ra12y*t1 + ra22y*t2 + ra32y*t3 &
                                    + ra42y*t4 + ra52y*t5 + t6*(de-ny*dq))*nm
    flxdiss(4) = 0.5_rfreal*dissfact*(ra12z*t1 + ra22z*t2 + ra32z*t3 &
                                    + ra42z*t4 + ra52z*t5 + t6*(dw-nz*dq))*nm
    flxdiss(5) = 0.5_rfreal*dissfact*(ra13*t1  + ra23*t2  + ra33*t3 &
                                    + ra43*t4  + ra53*t5 &
                                    + t6*(uh*du+vh*de+wh*dw-qh*dq))*nm

! ==============================================================================
!   Store mass flux
! ==============================================================================

    pmf(indmf*ifg) = flxconv(1) - flxdiss(1)

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + (flxconv(5) - flxdiss(5))

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - (flxconv(5) - flxdiss(5))
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    CALL rflu_centralfirstpatch_gl(pregion,ppatch)
  END DO ! iPatch

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computeflux1_gl 








! ******************************************************************************
!
! Purpose: Compute convective fluxes using first-order accurate Roe scheme
!  for calorically and thermally perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computeflux1_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_c_ghovm2, &
                           mixtperf_ho_cptuvw, &
                           rflu_centralfirstpatch

  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
  REAL(RFREAL) :: ah,cp,dissfact,dp,dq,dr,du,de,dv1,dv2,dv5,dw,efc, &
                  epsentr,fs,g,hh,hl,hr,irl,irr,l1,l2,l5,nm,nx,ny,nz,pl, &
                  pr,qh,ql,qr,rl,rh,rr,sh,term,tl,tr,t1,t2,t3,t5,ul,uh,ur,vl, &
                  vh,vr,wl,wh,wr,wt
  REAL(RFREAL) :: flxconv(5),flxdiss(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs,psd
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid    
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFlux1_TCP',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFlux1_TCP")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs    = pregion%grid%indGs
  indmf    = pregion%mixtInput%indMfMixt
  indsd    = pregion%mixtInput%indSd
  epsentr  = pregion%mixtInput%epsentr
  dissfact = pregion%mixtInput%dissFact

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv

  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd   

! ******************************************************************************
! Define constants
! ******************************************************************************

  cp = global%refCp
  g  = global%refGamma

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl

    ul = pcv(cv_mixt_xmom,c1)*irl
    vl = pcv(cv_mixt_ymom,c1)*irl
    wl = pcv(cv_mixt_zmom,c1)*irl
    pl = pdv(dv_mixt_pres,c1)
    tl = pdv(dv_mixt_temp,c1)

    hl = mixtperf_ho_cptuvw(cp,tl,ul,vl,wl)

    ql = ul*nx + vl*ny + wl*nz - fs

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr

    ur = pcv(cv_mixt_xmom,c2)*irr
    vr = pcv(cv_mixt_ymom,c2)*irr
    wr = pcv(cv_mixt_zmom,c2)*irr
    pr = pdv(dv_mixt_pres,c2)
    tr = pdv(dv_mixt_temp,c2)

    hr = mixtperf_ho_cptuvw(cp,tr,ur,vr,wr)

    qr = ur*nx + vr*ny + wr*nz - fs

! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)
    wt = rl/(rl + rh)

    uh = wt*ul + (1.0_rfreal-wt)*ur
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    hh = wt*hl + (1.0_rfreal-wt)*hr

    qh = uh*nx + vh*ny + wh*nz - fs
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    ah = mixtperf_c_ghovm2(g,hh,2.0_rfreal*sh) ! NOTE factor of 2

! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh - ah)
    l2 = abs(qh     )
    l5 = abs(qh + ah)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr = rr - rl
    du = ur - ul
    de = vr - vl ! use de to denote dv because dv is already used
    dw = wr - wl
    dp = pr - pl
    dq = du*nx + de*ny + dw*nz ! Note: no gs contribution

    dv1 = (dp - rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv5 = (dp + rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv2 = dr - dp/(ah*ah)

    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l2*rh
    t5 = l5*dv5

! ==============================================================================
!   Compute fluxes
! ==============================================================================

! ------------------------------------------------------------------------------
!   Central part
! ------------------------------------------------------------------------------

    flxconv(1) = 0.5_rfreal*(ql*rl            + qr*rr           )*nm
    flxconv(2) = 0.5_rfreal*(ql*rl*ul + pl*nx + qr*rr*ur + pr*nx)*nm
    flxconv(3) = 0.5_rfreal*(ql*rl*vl + pl*ny + qr*rr*vr + pr*ny)*nm
    flxconv(4) = 0.5_rfreal*(ql*rl*wl + pl*nz + qr*rr*wr + pr*nz)*nm
    flxconv(5) = 0.5_rfreal*(ql*rl*hl + pl*fs + qr*rr*hr + pr*fs)*nm

! ------------------------------------------------------------------------------
!   Dissipative part
! ------------------------------------------------------------------------------

    flxdiss(1) = 0.5_rfreal*dissfact*(t1            + t2                                   + t5           )*nm
    flxdiss(2) = 0.5_rfreal*dissfact*(t1*(uh-nx*ah) + t2*uh + t3*(du-nx*dq)                + t5*(uh+nx*ah))*nm
    flxdiss(3) = 0.5_rfreal*dissfact*(t1*(vh-ny*ah) + t2*vh + t3*(de-ny*dq)                + t5*(vh+ny*ah))*nm
    flxdiss(4) = 0.5_rfreal*dissfact*(t1*(wh-nz*ah) + t2*wh + t3*(dw-nz*dq)                + t5*(wh+nz*ah))*nm
    flxdiss(5) = 0.5_rfreal*dissfact*(t1*(hh-qh*ah) + t2*sh + t3*(uh*du+vh*de+wh*dw-qh*dq) + t5*(hh+qh*ah))*nm

! ==============================================================================
!   Store mass flux
! ==============================================================================

    pmf(indmf*ifg) = flxconv(1) - flxdiss(1)

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + (flxconv(5) - flxdiss(5))

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - (flxconv(5) - flxdiss(5))

! ==============================================================================
!   Accumulate into substantial derivative
! ==============================================================================

    psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + uh*pmf(indmf*ifg)
    psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vh*pmf(indmf*ifg)
    psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + wh*pmf(indmf*ifg)

    psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - uh*pmf(indmf*ifg)
    psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vh*pmf(indmf*ifg)
    psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - wh*pmf(indmf*ifg)
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      CALL rflu_nscbc_compfirstpatchflux(pregion,ppatch)
    ELSE
      CALL rflu_centralfirstpatch(pregion,ppatch)
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFlux1_TCP")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computeflux1_tcp









! ******************************************************************************
!
! Purpose: Compute convective fluxes using second-order accurate Roe scheme for
!   multiphase flow.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computeflux2_gl(pRegion)

  USE modinterfaces, ONLY: mixtgasliq_c, &
                           mixtliq_c2_bp, &
                           mixtliq_d_dobpppobttto, &
                           mixtperf_c2_grt, &
                           mixtperf_cv_cpr, & 
                           mixtperf_d_prt, &
                           mixtperf_r_cpg, &
                           mixtperf_r_m, &
                           rflu_centralfirstpatch_gl, &
                           rflu_centralsecondpatch_gl

  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

 TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,indmf,ipatch
  REAL(RFREAL) :: bgh2,bgl2,bgr2,blh2,bll2,blr2,bp,bt,bvh2,bvl2,bvr2,cmh, &
                  cml,cmr,cvg,cvl,cvv,cgh2,cgl2,cgr2,clh2,cll2,clr2,cvh2, &
                  cvl2,cvr2,cvm,dissfact,de,dp,dq,dr,dtemp,du,dv1,dv2,dv3, &
                  dv4,dv5,dvfg,dvfv,dw,dx,dy,dz,efc,egh,el,elh,epsentr,er, &
                  evh,factor1,factor2,factor4,factor5,fdv1,fdv2,fdv3,fs, &
                  fxn1,fxn2,fxn3,fxn4,fxn5,fxn6,invc2p,invc2pb,irl,irr,l1, &
                  l2,l3,l4,l5,nm,nx,ny,nz,ph,pl,po,pr,qh,ql,qr,ra11,ra12x, &
                  ra12y,ra12z,ra13,ra14,ra15,ra21,ra22x,ra22y,ra22z,ra23, &
                  ra24,ra25,ra31,ra32x,ra32y,ra32z,ra33,ra34,ra35,ra41, &
                  ra42x,ra42y,ra42z,ra43,ra44,ra45,ra51,ra52x,ra52y,ra52z, &
                  ra53,ra54,ra55,rel,rer,rgh,rgl,rgr,rh,rl,rlh,rll,rlr,ro, &
                  rr,rvh,rvl,rvr,rygl,rygr,ryvl,ryvr,rg,rv,sh,t1,t2,t3,t4, &
                  t5,t6,th,tl,to,tr,uh,ul,ur,vfgh,vfgl,vfgr,vflh,vfll,vflr, &
                  vfvh,vfvl,vfvr,vh,vl,vr,vl2,vr2,ygl,ygr,yll,ylr,yvl,yvr, & 
                  wh,wl,wr,wt,xc,yc,zc  
  REAL(RFREAL) :: flxconv(5),flxdiss(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvmixt,pcvspec,pdvmixt,prhs
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradmixt,pgradspec
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid  
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFlux2_GL',&
  'RFLU_ModRoeFlux.F90')

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs    = pregion%grid%indGs
  indmf    = pregion%mixtInput%indMfMixt
  epsentr  = pregion%mixtInput%epsentr
  dissfact = pregion%mixtInput%dissFact

  pcvmixt => pregion%mixt%cv
  pdvmixt => pregion%mixt%dv
  pcvspec => pregion%spec%cv

  pgradmixt => pregion%mixt%gradCell
  pgradspec => pregion%spec%gradCell

  pmf   => pregion%mixt%mfMixt
  prhs  => pregion%mixt%rhs

! ******************************************************************************
! Define constants
! ******************************************************************************

  ro  = global%refDensityLiq
  po  = global%refPressLiq
  to  = global%refTempLiq
  bp  = global%refBetaPLiq
  bt  = global%refBetaTLiq
  cvl = global%refCvLiq

  rg  = mixtperf_r_m(pregion%specInput%specType(1)%pMaterial%molw)
  cvg = mixtperf_cv_cpr(pregion%specInput%specType(1)%pMaterial%spht,rg)

  rv  = mixtperf_r_m(pregion%specInput%specType(2)%pMaterial%molw)
  cvv = mixtperf_cv_cpr(pregion%specInput%specType(2)%pMaterial%spht,rv)

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    xc = pgrid%fc(xcoord,ifg)
    yc = pgrid%fc(ycoord,ifg)
    zc = pgrid%fc(zcoord,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! -----------------------------------------------------------------------------
!   Left state
! -----------------------------------------------------------------------------    
  
    rl  = pcvmixt(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl
    
    ul = pcvmixt(cv_mixt_xmom,c1)*irl
    vl = pcvmixt(cv_mixt_ymom,c1)*irl
    wl = pcvmixt(cv_mixt_zmom,c1)*irl
    pl = pdvmixt(dv_mixt_pres,c1)
    tl = pdvmixt(dv_mixt_temp,c1)
    
    ygl = pcvspec(1,c1)*irl
    yvl = pcvspec(2,c1)*irl
         
    dx = xc - pgrid%cofg(xcoord,c1)
    dy = yc - pgrid%cofg(ycoord,c1)
    dz = zc - pgrid%cofg(zcoord,c1)    

    pl = pl + pgradmixt(xcoord,grc_mixt_pres,c1)*dx & 
            + pgradmixt(ycoord,grc_mixt_pres,c1)*dy & 
            + pgradmixt(zcoord,grc_mixt_pres,c1)*dz 
    ul = ul + pgradmixt(xcoord,grc_mixt_xvel,c1)*dx & 
            + pgradmixt(ycoord,grc_mixt_xvel,c1)*dy & 
            + pgradmixt(zcoord,grc_mixt_xvel,c1)*dz     
    vl = vl + pgradmixt(xcoord,grc_mixt_yvel,c1)*dx & 
            + pgradmixt(ycoord,grc_mixt_yvel,c1)*dy & 
            + pgradmixt(zcoord,grc_mixt_yvel,c1)*dz
    wl = wl + pgradmixt(xcoord,grc_mixt_zvel,c1)*dx & 
            + pgradmixt(ycoord,grc_mixt_zvel,c1)*dy & 
            + pgradmixt(zcoord,grc_mixt_zvel,c1)*dz
    tl = tl + pgradmixt(xcoord,grc_mixt_temp,c1)*dx & 
            + pgradmixt(ycoord,grc_mixt_temp,c1)*dy & 
            + pgradmixt(zcoord,grc_mixt_temp,c1)*dz 
                    
    ygl = ygl + pgradspec(xcoord,1,c1)*dx &
              + pgradspec(ycoord,1,c1)*dy &
              + pgradspec(zcoord,1,c1)*dz
    yvl = yvl + pgradspec(xcoord,2,c1)*dx &
              + pgradspec(ycoord,2,c1)*dy &
              + pgradspec(zcoord,2,c1)*dz

    yll = 1.0_rfreal - ygl - yvl

    vl2 = ul*ul + vl*vl + wl*wl
    ql  = ul*nx + vl*ny + wl*nz - fs

    rll = mixtliq_d_dobpppobttto(ro,bp,bt,pl,po,tl,to)
    rvl = mixtperf_d_prt(pl,rv,tl)
    rgl = mixtperf_d_prt(pl,rg,tl) 
    rl  = 1.0_rfreal/(yll/rll + yvl/rvl + ygl/rgl)

    vfgl = rl*ygl/rgl
    vfvl = rl*yvl/rvl
    vfll = 1.0_rfreal - vfgl - vfvl

    rel = (rll*vfll*cvl + rgl*vfgl*cvg + rvl*vfvl*cvv)*tl + 0.5_rfreal*rl*vl2   
    el = rel/rl
      
    cll2 = mixtliq_c2_bp(bp) 
    cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
    cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)

    bll2 = -bt/bp
    bvl2 = rvl*rv
    bgl2 = rgl*rg
    
    cvm = (rll*vfll*cvl + rvl*vfvl*cvv + rgl*vfgl*cvg)/rl
 
    cml = mixtgasliq_c(cvm,rl,pl,rll,rvl,rgl,vfll,vfvl,vfgl,cll2,cvl2,cgl2, &
                       bll2,bvl2,bgl2)
    
!----------------------------------------------------------------------------
!   Right state
! -----------------------------------------------------------------------------    
    
    rr  = pcvmixt(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    
    ur = pcvmixt(cv_mixt_xmom,c2)*irr
    vr = pcvmixt(cv_mixt_ymom,c2)*irr
    wr = pcvmixt(cv_mixt_zmom,c2)*irr
    pr = pdvmixt(dv_mixt_pres,c2)
    tr = pdvmixt(dv_mixt_temp,c2)
    
    ygr = pcvspec(1,c2)*irr
    yvr = pcvspec(2,c2)*irr 
       
    dx = xc - pgrid%cofg(xcoord,c2)
    dy = yc - pgrid%cofg(ycoord,c2)
    dz = zc - pgrid%cofg(zcoord,c2)    

    pr = pr + pgradmixt(xcoord,grc_mixt_pres,c2)*dx & 
            + pgradmixt(ycoord,grc_mixt_pres,c2)*dy & 
            + pgradmixt(zcoord,grc_mixt_pres,c2)*dz 
    ur = ur + pgradmixt(xcoord,grc_mixt_xvel,c2)*dx & 
            + pgradmixt(ycoord,grc_mixt_xvel,c2)*dy & 
            + pgradmixt(zcoord,grc_mixt_xvel,c2)*dz     
    vr = vr + pgradmixt(xcoord,grc_mixt_yvel,c2)*dx & 
            + pgradmixt(ycoord,grc_mixt_yvel,c2)*dy & 
            + pgradmixt(zcoord,grc_mixt_yvel,c2)*dz
    wr = wr + pgradmixt(xcoord,grc_mixt_zvel,c2)*dx & 
            + pgradmixt(ycoord,grc_mixt_zvel,c2)*dy & 
            + pgradmixt(zcoord,grc_mixt_zvel,c2)*dz
    tr = tr + pgradmixt(xcoord,grc_mixt_temp,c2)*dx & 
            + pgradmixt(ycoord,grc_mixt_temp,c2)*dy & 
            + pgradmixt(zcoord,grc_mixt_temp,c2)*dz 
                    
    ygr = ygr + pgradspec(xcoord,1,c2)*dx &
              + pgradspec(ycoord,1,c2)*dy &
              + pgradspec(zcoord,1,c2)*dz
    yvr = yvr + pgradspec(xcoord,2,c2)*dx &
              + pgradspec(ycoord,2,c2)*dy &
              + pgradspec(zcoord,2,c2)*dz

    ylr = 1.0_rfreal - ygr - yvr
 
    vr2 = ur*ur + vr*vr + wr*wr
    qr  = ur*nx + vr*ny + wr*nz - fs
   
    rlr = mixtliq_d_dobpppobttto(ro,bp,bt,pr,po,tr,to)
    rvr = mixtperf_d_prt(pr,rv,tr)
    rgr = mixtperf_d_prt(pr,rg,tr)
    rr  = 1.0_rfreal/(ylr/rlr + yvr/rvr + ygr/rgr)  
   
    vfgr = rr*ygr/rgr
    vfvr = rr*yvr/rvr
    vflr = 1.0_rfreal - vfgr - vfvr

    rer = (rlr*vflr*cvl + rgr*vfgr*cvg + rvr*vfvr*cvv)*tr + 0.5_rfreal*rr*vr2
    er  = rer/rr
    
    clr2 = mixtliq_c2_bp(bp) 
    cvr2 = mixtperf_c2_grt(1.0_rfreal,rv,tr)
    cgr2 = mixtperf_c2_grt(1.0_rfreal,rg,tr)

    blr2 = -bt/bp
    bvr2 = rvr*rv
    bgr2 = rgr*rg
    
    cvm = (rlr*vflr*cvl + rvr*vfvr*cvv + rgr*vfgr*cvg)/rr
 
    cmr = mixtgasliq_c(cvm,rr,pr,rlr,rvr,rgr,vflr,vfvr,vfgr,clr2,cvr2,cgr2, &
                       blr2,bvr2,bgr2)

! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)                     
    wt = rl/(rl + rh)                     
    uh = wt*ul + (1.0_rfreal-wt)*ur     
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    ph = wt*pl + (1.0_rfreal-wt)*pr      
    th = wt*tl + (1.0_rfreal-wt)*tr       

    vfgh = wt*vfgl + (1.0_rfreal-wt)*vfgr
    vfvh = wt*vfvl + (1.0_rfreal-wt)*vfvr
    vflh = 1.0_rfreal - vfvh - vfgh

    qh = uh*nx + vh*ny + wh*nz - fs 
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    elh = cvl*th + sh
    evh = cvv*th + sh
    egh = cvg*th + sh

    clh2 = mixtliq_c2_bp(bp) 
    cvh2 = mixtperf_c2_grt(1.0_rfreal,rv,th)
    cgh2 = mixtperf_c2_grt(1.0_rfreal,rg,th)

    rlh = mixtliq_d_dobpppobttto(ro,bp,bt,ph,po,th,to)
    rvh = mixtperf_d_prt(ph,rv,th)
    rgh = mixtperf_d_prt(ph,rg,th)

    blh2 = -bt/bp
    bvh2 = rvh*rv
    bgh2 = rgh*rg

    cvm = (rlh*vflh*cvl + rvh*vfvh*cvv + rgh*vfgh*cvg)/rh
 
    cmh = mixtgasliq_c(cvm,rh,ph,rlh,rvh,rgh,vflh,vfvh,vfgh,clh2,cvh2,cgh2, &
                       blh2,bvh2,bgh2) 

! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh      )
    l2 = abs(qh      )
    l3 = abs(qh      )
    l4 = abs(qh - cmh)
    l5 = abs(qh + cmh)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5 
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l3  = entropyfixhartenhyman(l3,efc)
    l4  = entropyfixhartenhyman(l4,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr    = rr - rl     
    du    = ur - ul
    de    = vr - vl
    dw    = wr - wl  
    dq    = du*nx + de*ny + dw*nz   ! Note: no gs contribution
    dp    = pr - pl 
    dtemp = tr - tl
    dvfv  = vfvr - vfvl
    dvfg  = vfgr - vfgl

    fdv1 = ph/(rh*cmh*cmh*rh*cvm*th)
    fdv2 = (vfgh*(rgh*cgh2 - rh*cmh*cmh &
         + (bgh2*ph)/(rh*cvm)))/(rgh*cgh2*rh*cmh*cmh)
    fdv3 = (vfvh*(rvh*cvh2 - rh*cmh*cmh &
         + (bvh2*ph)/(rh*cvm)))/(rvh*cvh2*rh*cmh*cmh)

    dv1  = (dtemp/th) - (dp*fdv1)
    dv2  = dvfg - dp*fdv2
    dv3  = dvfv - dp*fdv3
    dv4  = -dq/cmh + dp/(rh*cmh*cmh) 
    dv5  =  dq/cmh + dp/(rh*cmh*cmh) 

    invc2p  = vflh/clh2 + vfgh/cgh2 + vfvh/cvh2 
    invc2pb = (blh2*vflh)/clh2 + (bgh2*vfgh)/cgh2 &
            + (bvh2*vfvh)/cvh2

    ra11  = -th*invc2pb
    ra12x = -th*invc2pb*uh
    ra12y = -th*invc2pb*vh
    ra12z = -th*invc2pb*wh
    ra13  = rlh*vflh*cvl*th + bt*cvl*th*th*vflh - th*invc2pb*sh
    ra14  = -(bgh2*vfgh*th)/cgh2
    ra15  = -(bvh2*vfvh*th)/cvh2

    ra21  =  rgh - rlh
    ra22x = (rgh - rlh)*uh
    ra22y = (rgh - rlh)*vh
    ra22z = (rgh - rlh)*wh
    ra23  =  rgh*egh - rlh*elh
    ra24  =  rgh
    ra25  =  0.0_rfreal    

    ra31  =  rvh - rlh
    ra32x = (rvh - rlh)*uh
    ra32y = (rvh - rlh)*vh
    ra32z = (rvh - rlh)*wh
    ra33  =  rvh*evh - rlh*elh
    ra34  =  0.0_rfreal
    ra35  =  rvh    
    
    fxn1 = ((rgh - rlh)*vfgh*(rgh*cgh2 - rh*cmh*cmh &
         + (bgh2*ph)/(rh*cvm)))/(2.0_rfreal*rgh*cgh2)
    fxn2 = ((rvh - rlh)*vfvh*(rvh*cvh2 - rh*cmh*cmh &
         + (bvh2*ph)/(rh*cvm)))/(2.0_rfreal*rvh*cvh2)
    factor1 = 0.5_rfreal*rh*cmh*cmh*invc2p & 
            - (ph*invc2pb)/(2.0_rfreal*rh*cvm) + fxn1 + fxn2

    fxn3 = 0.5_rfreal*rh*cmh*cmh*((elh*vflh)/clh2 &
         + (egh*vfgh)/cgh2 + (evh*vfvh)/cvh2)
    fxn4 = (ph*(rlh*vflh*cvl + bt*cvl*th*vflh &
         - (sh*invc2pb)))/(2.0_rfreal*rh*cvm)
    fxn5 = ((rgh*egh - rlh*elh)*vfgh*(rgh*cgh2 - rh*cmh*cmh &
         + (bgh2*ph)/(rh*cvm)))/(2.0_rfreal*rgh*cgh2)
    fxn6 = ((rvh*evh - rlh*elh)*vfvh*(rvh*cvh2 - rh*cmh*cmh &
         +  (bvh2*ph)/(rh*cvm)))/(2.0_rfreal*rvh*cvh2)
    factor2 = fxn3 + fxn4 + fxn5 + fxn6

    factor4 = 0.5_rfreal*rgh*vfgh
    factor5 = 0.5_rfreal*rvh*vfvh

    ra41  = factor1
    ra42x = factor1*uh - 0.5_rfreal*rh*cmh*nx 
    ra42y = factor1*vh - 0.5_rfreal*rh*cmh*ny
    ra42z = factor1*wh - 0.5_rfreal*rh*cmh*nz
    ra43  = factor2    - 0.5_rfreal*rh*cmh*qh
    ra44  = factor4
    ra45  = factor5

    ra51  = factor1
    ra52x = factor1*uh + 0.5_rfreal*rh*cmh*nx 
    ra52y = factor1*vh + 0.5_rfreal*rh*cmh*ny
    ra52z = factor1*wh + 0.5_rfreal*rh*cmh*nz
    ra53  = factor2    + 0.5_rfreal*rh*cmh*qh
    ra54  = factor4
    ra55  = factor5   
 
    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l3*dv3
    t4 = l4*dv4
    t5 = l5*dv5
    t6 = l1*rh 

! ==============================================================================
!   Compute fluxes
! ==============================================================================

! ------------------------------------------------------------------------------
!   Central part
! ------------------------------------------------------------------------------
  
    flxconv(1) = 0.5_rfreal*(ql*rl            + qr*rr           )*nm
    flxconv(2) = 0.5_rfreal*(ql*rl*ul + pl*nx + qr*rr*ur + pr*nx)*nm
    flxconv(3) = 0.5_rfreal*(ql*rl*vl + pl*ny + qr*rr*vr + pr*ny)*nm
    flxconv(4) = 0.5_rfreal*(ql*rl*wl + pl*nz + qr*rr*wr + pr*nz)*nm
    flxconv(5) = 0.5_rfreal*(ql*(rl*el + pl) + pl*fs  + qr*(rr*er + pr) &
                           + pr*fs)*nm

! ------------------------------------------------------------------------------
!   Dissipative part
! ------------------------------------------------------------------------------

    flxdiss(1) = 0.5_rfreal*dissfact*(ra11*t1  + ra21*t2  + ra31*t3 &
                                    + ra41*t4  + ra51*t5)*nm 
    flxdiss(2) = 0.5_rfreal*dissfact*(ra12x*t1 + ra22x*t2 + ra32x*t3 &
                                    + ra42x*t4 + ra52x*t5 + t6*(du-nx*dq))*nm
    flxdiss(3) = 0.5_rfreal*dissfact*(ra12y*t1 + ra22y*t2 + ra32y*t3 &
                                    + ra42y*t4 + ra52y*t5 + t6*(de-ny*dq))*nm
    flxdiss(4) = 0.5_rfreal*dissfact*(ra12z*t1 + ra22z*t2 + ra32z*t3 &
                                    + ra42z*t4 + ra52z*t5 + t6*(dw-nz*dq))*nm
    flxdiss(5) = 0.5_rfreal*dissfact*(ra13*t1  + ra23*t2  + ra33*t3 &
                                    + ra43*t4  + ra53*t5 &
                                    + t6*(uh*du+vh*de+wh*dw-qh*dq))*nm

! ==============================================================================
!   Store mass flux
! ==============================================================================

    pmf(indmf*ifg) = flxconv(1) - flxdiss(1)

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + (flxconv(5) - flxdiss(5))

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - (flxconv(5) - flxdiss(5))
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_nscbc_compfirstpatchflux(pregion,ppatch)
        CASE ( 2 ) 
          CALL rflu_nscbc_compsecondpatchflux(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    ELSE
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_centralfirstpatch_gl(pregion,ppatch)      
        CASE ( 2 ) 
          CALL rflu_centralsecondpatch_gl(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computeflux2_gl 










! ******************************************************************************
!
! Purpose: Compute convective fluxes using second-order accurate Roe scheme
!  for calorically and thermally perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_roe_computeflux2_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_c_ghovm2, &
                           mixtperf_ho_cptuvw, &
                           mixtperf_r_cpg, &
                           mixtperf_t_dpr, &
                           rflu_centralfirstpatch, & 
                           rflu_centralsecondpatch

  IMPLICIT NONE

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
  REAL(RFREAL) :: ah,cp,de,dissfact,dp,dq,dr,du,dv1,dv2,dv5,dw,dx,dy,dz, &
                  efc,epsentr,fs,g,gc,hh,hl,hr,irl,irr,l1,l2,l5,nm,nx,ny,nz, &
                  pl,pr,qh,ql,qr,rl,rh,rr,sh,term,tl,tr,t1,t2,t3,t5,ul, &
                  uh,ur,vl,vh,vr,wl,wh,wr,wt,xc,yc,zc
  REAL(RFREAL) :: flxconv(5),flxdiss(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,prhs,psd
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ROE_ComputeFlux2_TCP',&
  'RFLU_ModRoeFlux.F90')

#ifdef ROCPROF 
  CALL fprofiler_begins("RFLU::ROE_ComputeFlux2_TCP")
#endif

! ******************************************************************************
! Set dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs    = pregion%grid%indGs
  indmf    = pregion%mixtInput%indMfMixt
  indsd    = pregion%mixtInput%indSd
  epsentr  = pregion%mixtInput%epsentr
  dissfact = pregion%mixtInput%dissFact

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  
  pgc  => pregion%mixt%gradCell  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd   
    
! ******************************************************************************
! Define constants
! ******************************************************************************

  cp = global%refCp
  g  = global%refGamma
  gc = mixtperf_r_cpg(cp,g)

! ******************************************************************************
! Compute fluxes through interior faces
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================
!   Get face geometry and grid speed
! ==============================================================================

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)

    xc = pgrid%fc(xcoord,ifg)
    yc = pgrid%fc(ycoord,ifg)
    zc = pgrid%fc(zcoord,ifg)

    fs = pgrid%gs(indgs*ifg)

! ==============================================================================
!   Compute left and right states
! ==============================================================================

! ------------------------------------------------------------------------------
!   Left state
! ------------------------------------------------------------------------------

    rl  = pcv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl

    ul = pcv(cv_mixt_xmom,c1)*irl
    vl = pcv(cv_mixt_ymom,c1)*irl
    wl = pcv(cv_mixt_zmom,c1)*irl
    pl = pdv(dv_mixt_pres,c1)

    dx = xc - pgrid%cofg(xcoord,c1)
    dy = yc - pgrid%cofg(ycoord,c1)
    dz = zc - pgrid%cofg(zcoord,c1)

    rl = rl + pgc(xcoord,grc_mixt_dens,c1)*dx &
            + pgc(ycoord,grc_mixt_dens,c1)*dy &
            + pgc(zcoord,grc_mixt_dens,c1)*dz
    ul = ul + pgc(xcoord,grc_mixt_xvel,c1)*dx &
            + pgc(ycoord,grc_mixt_xvel,c1)*dy &
            + pgc(zcoord,grc_mixt_xvel,c1)*dz
    vl = vl + pgc(xcoord,grc_mixt_yvel,c1)*dx &
            + pgc(ycoord,grc_mixt_yvel,c1)*dy &
            + pgc(zcoord,grc_mixt_yvel,c1)*dz
    wl = wl + pgc(xcoord,grc_mixt_zvel,c1)*dx &
            + pgc(ycoord,grc_mixt_zvel,c1)*dy &
            + pgc(zcoord,grc_mixt_zvel,c1)*dz
    pl = pl + pgc(xcoord,grc_mixt_pres,c1)*dx &
            + pgc(ycoord,grc_mixt_pres,c1)*dy &
            + pgc(zcoord,grc_mixt_pres,c1)*dz

    tl = mixtperf_t_dpr(rl,pl,gc)
    hl = mixtperf_ho_cptuvw(cp,tl,ul,vl,wl)

    ql = ul*nx + vl*ny + wl*nz - fs

! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr

    ur = pcv(cv_mixt_xmom,c2)*irr
    vr = pcv(cv_mixt_ymom,c2)*irr
    wr = pcv(cv_mixt_zmom,c2)*irr
    pr = pdv(dv_mixt_pres,c2)

    dx = xc - pgrid%cofg(xcoord,c2)
    dy = yc - pgrid%cofg(ycoord,c2)
    dz = zc - pgrid%cofg(zcoord,c2)

    rr = rr + pgc(xcoord,grc_mixt_dens,c2)*dx &
            + pgc(ycoord,grc_mixt_dens,c2)*dy &
            + pgc(zcoord,grc_mixt_dens,c2)*dz
    ur = ur + pgc(xcoord,grc_mixt_xvel,c2)*dx &
            + pgc(ycoord,grc_mixt_xvel,c2)*dy &
            + pgc(zcoord,grc_mixt_xvel,c2)*dz
    vr = vr + pgc(xcoord,grc_mixt_yvel,c2)*dx &
            + pgc(ycoord,grc_mixt_yvel,c2)*dy &
            + pgc(zcoord,grc_mixt_yvel,c2)*dz
    wr = wr + pgc(xcoord,grc_mixt_zvel,c2)*dx &
            + pgc(ycoord,grc_mixt_zvel,c2)*dy &
            + pgc(zcoord,grc_mixt_zvel,c2)*dz
    pr = pr + pgc(xcoord,grc_mixt_pres,c2)*dx &
            + pgc(ycoord,grc_mixt_pres,c2)*dy &
            + pgc(zcoord,grc_mixt_pres,c2)*dz

    tr = mixtperf_t_dpr(rr,pr,gc)
    hr = mixtperf_ho_cptuvw(cp,tr,ur,vr,wr)

    qr = ur*nx + vr*ny + wr*nz - fs

! ==============================================================================
!   Compute Roe-averaged face variables (h for hat)
! ==============================================================================

    rh = sqrt(rl*rr)
    wt = rl/(rl + rh)

    uh = wt*ul + (1.0_rfreal-wt)*ur
    vh = wt*vl + (1.0_rfreal-wt)*vr
    wh = wt*wl + (1.0_rfreal-wt)*wr
    hh = wt*hl + (1.0_rfreal-wt)*hr

    qh = uh*nx + vh*ny + wh*nz - fs
    sh = 0.5_rfreal*(uh*uh + vh*vh + wh*wh)

    ah = mixtperf_c_ghovm2(g,hh,2.0_rfreal*sh) ! NOTE factor of 2

! ==============================================================================
!   Compute eigenvalues
! ==============================================================================

    l1 = abs(qh - ah)
    l2 = abs(qh     )
    l5 = abs(qh + ah)

! ==============================================================================
!   Compute entropy fix
! ==============================================================================

    efc = epsentr*l5
    l1  = entropyfixhartenhyman(l1,efc)
    l2  = entropyfixhartenhyman(l2,efc)
    l5  = entropyfixhartenhyman(l5,efc)

! ==============================================================================
!   Compute wavestrengths
! ==============================================================================

    dr = rr - rl
    du = ur - ul
    de = vr - vl ! use de to denote dv because dv is already used
    dw = wr - wl
    dp = pr - pl
    dq = du*nx + de*ny + dw*nz ! Note: no gs contribution

    dv1 = (dp - rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv5 = (dp + rh*ah*dq)/(2.0_rfreal*ah*ah)
    dv2 = dr - dp/(ah*ah)

    t1 = l1*dv1
    t2 = l2*dv2
    t3 = l2*rh
    t5 = l5*dv5

! =============================================================================
!   Compute fluxes
! =============================================================================

! ------------------------------------------------------------------------------
!   Central part
! ------------------------------------------------------------------------------

    flxconv(1) = 0.5_rfreal*(ql*rl            + qr*rr           )*nm
    flxconv(2) = 0.5_rfreal*(ql*rl*ul + pl*nx + qr*rr*ur + pr*nx)*nm
    flxconv(3) = 0.5_rfreal*(ql*rl*vl + pl*ny + qr*rr*vr + pr*ny)*nm
    flxconv(4) = 0.5_rfreal*(ql*rl*wl + pl*nz + qr*rr*wr + pr*nz)*nm
    flxconv(5) = 0.5_rfreal*(ql*rl*hl + pl*fs + qr*rr*hr + pr*fs)*nm

! ------------------------------------------------------------------------------
!   Dissipative part
! ------------------------------------------------------------------------------

    flxdiss(1) = 0.5_rfreal*dissfact*(t1            + t2                                   + t5           )*nm
    flxdiss(2) = 0.5_rfreal*dissfact*(t1*(uh-nx*ah) + t2*uh + t3*(du-nx*dq)                + t5*(uh+nx*ah))*nm
    flxdiss(3) = 0.5_rfreal*dissfact*(t1*(vh-ny*ah) + t2*vh + t3*(de-ny*dq)                + t5*(vh+ny*ah))*nm
    flxdiss(4) = 0.5_rfreal*dissfact*(t1*(wh-nz*ah) + t2*wh + t3*(dw-nz*dq)                + t5*(wh+nz*ah))*nm
    flxdiss(5) = 0.5_rfreal*dissfact*(t1*(hh-qh*ah) + t2*sh + t3*(uh*du+vh*de+wh*dw-qh*dq) + t5*(hh+qh*ah))*nm

! =============================================================================
!   Store mass flux
! =============================================================================

    pmf(indmf*ifg) = flxconv(1) - flxdiss(1)

! ==============================================================================
!   Accumulate into residual
! ==============================================================================

    prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + (flxconv(5) - flxdiss(5))

    prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - (flxconv(1) - flxdiss(1))
    prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - (flxconv(2) - flxdiss(2))
    prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - (flxconv(3) - flxdiss(3))
    prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - (flxconv(4) - flxdiss(4))
    prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - (flxconv(5) - flxdiss(5))

! ==============================================================================
!   Accumulate into substantial derivative
! ==============================================================================

    psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + uh*pmf(indmf*ifg)
    psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vh*pmf(indmf*ifg)
    psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + wh*pmf(indmf*ifg)
   
    psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - uh*pmf(indmf*ifg)
    psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vh*pmf(indmf*ifg)
    psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - wh*pmf(indmf*ifg)
  END DO ! ifg

! ******************************************************************************
! Compute fluxes through boundary faces
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_nscbc_compfirstpatchflux(pregion,ppatch)
        CASE ( 2 ) 
          CALL rflu_nscbc_compsecondpatchflux(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    ELSE
      SELECT CASE ( ppatch%spaceOrder )
        CASE ( 1 ) 
          CALL rflu_centralfirstpatch(pregion,ppatch)      
        CASE ( 2 ) 
          CALL rflu_centralsecondpatch(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%spaceOrder
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

#ifdef ROCPROF 
  CALL fprofiler_ends("RFLU::ROE_ComputeFlux2_TCP")
#endif

  CALL deregisterfunction(global)

END SUBROUTINE rflu_roe_computeflux2_tcp








! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modroeflux

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModRoeFlux.F90,v $
! Revision 1.10  2008/12/06 08:44:24  mtcampbe
! Updated license.
!
! Revision 1.9  2008/11/19 22:17:35  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.8  2006/08/19 15:39:18  mparmar
! Added computation of boundary flux using boundary variables
!
! Revision 1.7  2006/05/03 17:55:42  haselbac
! Bug fix: Missing pGrid def in RFLU_ROE_ComputeFlux2_TCP
!
! Revision 1.6  2006/05/01 21:00:47  haselbac
! Rewrite for consistency and cleanliness
!
! Revision 1.5  2006/04/15 16:58:42  haselbac
! Added capability of running 1st order boundary fluxes with 2nd order volume fluxes
!
! Revision 1.4  2006/04/07 15:19:20  haselbac
! Removed tabs
!
! Revision 1.3  2006/03/26 20:22:09  haselbac
! Added support for GL model
!
! Revision 1.2  2005/07/07 22:45:01  haselbac
! Added profiling calls
!
! Revision 1.1  2005/05/16 20:36:29  haselbac
! Initial revision
!
! ******************************************************************************