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

MODULE rflu_modausmflux

  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

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: rflu_ausm_computeflux
  
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
      
  CHARACTER(CHRLEN) :: & 
    RCSIdentString = '$RCSfile: RFLU_ModAUSMFlux.F90,v $ $Revision: 1.11 $' 
              
! *****************************************************************************
! Routines
! *****************************************************************************

  CONTAINS






! ******************************************************************************
!
! Purpose: Wrapper function for AUSM flux functions.
!
! Description: None.
!
! Input: 
!   pRegion     Pointer to region data.
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_ausm_computeflux(pRegion)
                       
  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_AUSM_ComputeFlux',&
  'RFLU_ModAUSMFlux.F90')

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

  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_ausm_computeflux1(pregion)
        CASE ( discr_order_2 )
          CALL rflu_ausm_computeflux2_tcp(pregion)               
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pRegion%mixtInput%spaceOrder    

! ==============================================================================
!   Mixture of thermally and calorically perfect gases
! ==============================================================================

    CASE ( gas_model_mixt_tcperf ) 
      SELECT CASE ( pregion%mixtInput%spaceOrder ) 
        CASE ( discr_order_1 ) 
          CALL rflu_ausm_computeflux1(pregion)
        CASE ( discr_order_2 )
          CALL rflu_ausm_computeflux2_mtcp(pregion)               
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pRegion%mixtInput%spaceOrder

! ==============================================================================
!   Pseudo-gas
! ==============================================================================

    CASE ( gas_model_mixt_pseudo ) 
      SELECT CASE ( pregion%mixtInput%spaceOrder ) 
! TO DO
!        CASE ( DISCR_ORDER_1 ) 
!          CALL RFLU_AUSM_ComputeFlux1_MPSD(pRegion)
! END TO DO 
        CASE ( discr_order_2 )
          CALL rflu_ausm_computeflux2_mpsd(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       

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ausm_computeflux







! ******************************************************************************
!
! Purpose: Compute convective fluxes using AUSM+ scheme.
!
! Description: None.
!
! Input: 
!   nx          x-component of face normal
!   ny          y-component of face normal
!   nz          z-component of face normal
!   nm          Magnitude of face normal
!   fs          Face speed
!   rl          Density of left state
!   ul          x-component of velocity of left state
!   vl          y-component of velocity of left state
!   wl          z-component of velocity of left state   
!   Hl          Total enthalpy of left state
!   al          Speed of sound of left state
!   pl          Pressure of left state
!   rr          Density of right state
!   ur          x-component of velocity of right state
!   vr          y-component of velocity of right state
!   wr          z-component of velocity of right state  
!   pr          Pressure of right state
!   Hr          Total enthalpy of right state
!   ar          Speed of sound of right state
!
! Output: 
!   flx         Fluxes
!   vf          Face velocities
!
! Notes: 
!   1. Liou M.-S., Progress towards an improved CFD method: AUSM+, AIAA Paper
!      95-1701, 1995
!   2. Do not use computation of face speed of sound which leads to exact 
!      capturing of isolated normal shock waves because of robustness problems
!      for unsteady flows and because that formulation is not applicable to 
!      anything but calorically and thermally perfect gases.
!
! ******************************************************************************

SUBROUTINE rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,Hl,al,rr,ur, &
                                  vr,wr,pr,hr,ar,flx,vf)

  IMPLICIT NONE

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

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

  REAL(RFREAL), INTENT(IN) :: al,ar,fs,hl,hr,nm,nx,ny,nz,pl,pr,rl,rr,ul,ur,vl, &
                              vr,wl,wr
  REAL(RFREAL), INTENT(OUT) :: flx(5),vf(3)

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

  REAL(RFREAL) :: af,mf,mfa,mfm,mfp,ml,mla,mlp,mr,mra,mrm,pf,ql,qr,vml,vmr, &
                  wtl,wtr

! ******************************************************************************
! Start, compute face state
! ******************************************************************************

  ql = ul*nx + vl*ny + wl*nz - fs
  qr = ur*nx + vr*ny + wr*nz - fs
    
  af = 0.5_rfreal*(al+ar) ! NOTE not using original formulation, see note

  ml  = ql/af
  mla = abs(ml)

  mr  = qr/af
  mra = abs(mr)    

  IF ( mla <= 1.0_rfreal ) THEN 
    mlp = 0.25_rfreal*(ml+1.0_rfreal)*(ml+1.0_rfreal) & 
        + 0.125_rfreal*(ml*ml-1.0_rfreal)*(ml*ml-1.0_rfreal)
    wtl = 0.25_rfreal*(ml+1.0_rfreal)*(ml+1.0_rfreal)*(2.0_rfreal-ml) & 
        + 0.1875_rfreal*ml*(ml*ml-1.0_rfreal)*(ml*ml-1.0_rfreal)
  ELSE
    mlp = 0.5_rfreal*(ml+mla)
    wtl = 0.5_rfreal*(1.0_rfreal+ml/mla)
  END IF ! mla

  IF ( mra <= 1.0_rfreal ) THEN 
    mrm = -0.25_rfreal*(mr-1.0_rfreal)*(mr-1.0_rfreal) & 
          -0.125_rfreal*(mr*mr-1.0_rfreal)*(mr*mr-1.0_rfreal)
    wtr = 0.25_rfreal*(mr-1.0_rfreal)*(mr-1.0_rfreal)*(2.0_rfreal+mr) & 
        - 0.1875_rfreal*mr*(mr*mr-1.0_rfreal)*(mr*mr-1.0_rfreal)
  ELSE
    mrm = 0.5_rfreal*(mr-mra)
    wtr = 0.5_rfreal*(1.0_rfreal-mr/mra)
  END IF ! mla

  mf  = mlp + mrm
  mfa = abs(mf)
  mfp = 0.5_rfreal*(mf+mfa)
  mfm = 0.5_rfreal*(mf-mfa)

  pf = wtl*pl + wtr*pr 

! ******************************************************************************
! Compute fluxes
! ******************************************************************************

  vf(1) = mfp*ul + mfm*ur
  vf(2) = mfp*vl + mfm*vr
  vf(3) = mfp*wl + mfm*wr    

  flx(1) = (af*(mfp*rl    + mfm*rr   )        )*nm
  flx(2) = (af*(mfp*rl*ul + mfm*rr*ur) + pf*nx)*nm
  flx(3) = (af*(mfp*rl*vl + mfm*rr*vr) + pf*ny)*nm
  flx(4) = (af*(mfp*rl*wl + mfm*rr*wr) + pf*nz)*nm
  flx(5) = (af*(mfp*rl*hl + mfm*rr*hr) + pf*fs)*nm

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

END SUBROUTINE rflu_ausm_fluxfunction









! ******************************************************************************
!
! Purpose: Compute convective fluxes using first-order accurate AUSM+ scheme.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: 
!   1. This routine can also be used for a mixture of thermally and calorically 
!      perfect gases because, the face states being identical to the cell 
!      states, the mixture properties are constant.
!
! ******************************************************************************

SUBROUTINE rflu_ausm_computeflux1(pRegion)

  USE modinterfaces, ONLY: mixtperf_ho_cptuvw, & 
                           rflu_centralfirstpatch

  IMPLICIT NONE
  
! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indcp,indgs,indmf,indmol,indsd,ipatch
  REAL(RFREAL) :: al,ar,cpl,cpr,fs,hl,hr,irl,irr,nm,nx,ny,nz,pl,pr,rl,rr, &
                  tl,tr,ul,ur,vl,vr,wl,wr
  REAL(RFREAL) :: flx(5),vf(3)
  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_AUSM_ComputeFlux1',&
  'RFLU_ModAUSMFlux.F90')

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

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

  pgrid => pregion%grid

  indgs  = pregion%grid%indGs

  indcp  = pregion%mixtInput%indCp
  indmf  = pregion%mixtInput%indMfMixt  
  indmol = pregion%mixtInput%indMol
  indsd  = pregion%mixtInput%indSd

  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  pgv  => pregion%mixt%gv
  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd  

! ******************************************************************************
! 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
! ------------------------------------------------------------------------------

    cpl = pgv(gv_mixt_cp,indcp*c1)
    
    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)
    al  = pdv(dv_mixt_soun,c1)

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

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

    cpr = pgv(gv_mixt_cp,indcp*c2)

    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)    
    ar  = pdv(dv_mixt_soun,c2)
    
    hr  = mixtperf_ho_cptuvw(cpr,tr,ur,vr,wr)

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

    CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al,rr,ur,vr, &
                                wr,pr,hr,ar,flx,vf)

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

    pmf(indmf*ifg) = flx(1)

! ==============================================================================
!   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)
    
! ==============================================================================
!   Accumulate into substantial derivative
! ==============================================================================

    psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + vf(1)*flx(1)
    psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vf(2)*flx(1)
    psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + vf(3)*flx(1)
    
    psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - vf(1)*flx(1)
    psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vf(2)*flx(1)
    psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - vf(3)*flx(1)   
  END DO ! ifg

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

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

    CALL rflu_centralfirstpatch(pregion,ppatch)
  END DO ! iPatch

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ausm_computeflux1







! ******************************************************************************
!
! Purpose: Compute convective fluxes using second-order accurate AUSM+ scheme 
!   for mixture of thermally and calorically perfect gaseous species and 
!   particulate species.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_ausm_computeflux2_mpsd(pRegion)

#ifdef SPEC
  USE modspecies, ONLY: t_spec_type
#endif

  USE modinterfaces, ONLY: mixtperf_c_grt, &
                           mixtperf_g_cpr, &
                           mixtperf_ho_cptuvw, & 
                           mixtperf_m_r, & 
                           mixtperf_r_m, & 
                           mixtperf_t_dpr, & 
                           rflu_centralfirstpatch, &
                           rflu_centralsecondpatch

  IMPLICIT NONE
  
! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indcp,indgs,indmf,indmol,indsd,ipatch
  REAL(RFREAL) :: al,ar,cpl,cpr,dx,dy,dz,fs,gcl,gcr,gl,gr,hl,hr,irl,irr,mml, &
                  mmr,nm,nx,ny,nz,pl,pr,rl,rr,tl,tr,ul,ur,vl,vr,wl,wr,xc,yc,zc
  REAL(RFREAL) :: flx(5),vf(3)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs,psd
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

#ifdef SPEC
  INTEGER :: ispec
  REAL(RFREAL) :: gcg,immg,mmg,phip,yg,y1,y2
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvspec 
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgcspec 
  TYPE(t_spec_type), POINTER :: pspectype
#endif  

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_AUSM_ComputeFlux2_MPSD',&
  'RFLU_ModAUSMFlux.F90')

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

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

  pgrid => pregion%grid

  indgs  = pgrid%indGs
  
  indcp  = pregion%mixtInput%indCp
  indmf  = pregion%mixtInput%indMfMixt  
  indmol = pregion%mixtInput%indMol 
  indsd  = pregion%mixtInput%indSd
  
  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  pgv  => pregion%mixt%gv
  
  pgc  => pregion%mixt%gradCell  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd  

#ifdef SPEC
  pcvspec => pregion%spec%cv  
  pgcspec => pregion%spec%gradCell    
#endif

  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN 
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! pRegion%spec%cvState  

  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_pseudo ) THEN
    CALL errorstop(global,err_gasmodel_invalid,__line__)  
  END IF ! pRegion%mixtInput%gasModel

  IF ( (indcp /= 1) .OR. (indmol /= 1) ) THEN 
    CALL errorstop(global,err_reached_default,__line__)  
  END IF ! indCp

! ******************************************************************************
! 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
    
    dx  = xc - pgrid%cofg(xcoord,c1)
    dy  = yc - pgrid%cofg(ycoord,c1)
    dz  = zc - pgrid%cofg(zcoord,c1)

#ifdef SPEC
    immg = 0.0_rfreal
    cpl  = 0.0_rfreal
    yg   = 0.0_rfreal
    phip = 0.0_rfreal

    DO ispec = 1,pregion%specInput%nSpecies
      pspectype => pregion%specInput%specType(ispec)

      y1 = irl*pcvspec(ispec,c1) + pgcspec(xcoord,ispec,c1)*dx & 
                                 + pgcspec(ycoord,ispec,c1)*dy & 
                                 + pgcspec(zcoord,ispec,c1)*dz

      IF ( pspectype%discreteFlag .EQV. .true. ) THEN 
        cpl  = cpl  +    y1*pspectype%pMaterial%spht
        phip = phip + rl*y1/pspectype%pMaterial%dens
      ELSE 
        immg = immg + y1/pspectype%pMaterial%molw
        cpl  = cpl  + y1*pspectype%pMaterial%spht
        yg   = yg   + y1                                                              
      END IF ! pSpecType%discreteFlag
    END DO ! iSpec

    mmg = yg/immg
    gcg = mixtperf_r_m(mmg)
    gcl = yg/(1.0_rfreal-phip)*gcg
    mml = mixtperf_m_r(gcl)
#endif                    

    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)
        
    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,gcl)
    al = mixtperf_c_grt(gl,gcl,tl)
    hl = mixtperf_ho_cptuvw(cpl,tl,ul,vl,wl)
    
! ------------------------------------------------------------------------------
!   Right state
! ------------------------------------------------------------------------------

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    
    dx  = xc - pgrid%cofg(xcoord,c2)
    dy  = yc - pgrid%cofg(ycoord,c2)
    dz  = zc - pgrid%cofg(zcoord,c2)

#ifdef SPEC
    immg = 0.0_rfreal
    cpr  = 0.0_rfreal
    yg   = 0.0_rfreal
    phip = 0.0_rfreal

    DO ispec = 1,pregion%specInput%nSpecies
      pspectype => pregion%specInput%specType(ispec)

      y2 = irr*pcvspec(ispec,c2) + pgcspec(xcoord,ispec,c2)*dx & 
                                 + pgcspec(ycoord,ispec,c2)*dy & 
                                 + pgcspec(zcoord,ispec,c2)*dz

      IF ( pspectype%discreteFlag .EQV. .true. ) THEN 
        cpr  = cpr  +    y2*pspectype%pMaterial%spht
        phip = phip + rr*y2/pspectype%pMaterial%dens
      ELSE 
        immg = immg + y2/pspectype%pMaterial%molw
        cpr  = cpr  + y2*pspectype%pMaterial%spht
        yg   = yg   + y2                                                              
      END IF ! pSpecType%discreteFlag
    END DO ! iSpec

    mmg = yg/immg
    gcg = mixtperf_r_m(mmg)
    gcr = yg/(1.0_rfreal-phip)*gcg
    mmr = mixtperf_m_r(gcr)
#endif  
    
    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)
        
    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,gcr)
    ar = mixtperf_c_grt(gr,gcr,tr)
    hr = mixtperf_ho_cptuvw(cpr,tr,ur,vr,wr)
    
! ==============================================================================
!   Compute fluxes 
! ==============================================================================

    gl = cpl/(cpl-gcl)
    gr = cpr/(cpr-gcr)
    
    IF ( abs(gr-gl) < 0.05_rfreal ) THEN 
      CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al, & 
                                  rr,ur,vr,wr,pr,hr,ar,flx,vf)

      pmf(indmf*ifg) = flx(1)

      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)

      psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + vf(1)*flx(1)
      psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vf(2)*flx(1)
      psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + vf(3)*flx(1)

      psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - vf(1)*flx(1)
      psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vf(2)*flx(1)
      psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - vf(3)*flx(1)
    ELSE 
      CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al, & 
                                  rr,ur,vr,wr,pr,hr,ar,flx,vf)

      pmf(indmf*ifg) = 0.5_rfreal*flx(1)

      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)

      psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + vf(1)*flx(1)
      psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vf(2)*flx(1)
      psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + vf(3)*flx(1)
      
      CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al, & 
                                  rr,ur,vr,wr,pr,hr,ar,flx,vf)

      pmf(indmf*ifg) = pmf(indmf*ifg) + 0.5_rfreal*flx(1)

      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)            
      
      psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - vf(1)*flx(1)
      psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vf(2)*flx(1)
      psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - vf(3)*flx(1)      
    END IF ! ABS(gr-gl)   
  END DO ! ifg

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

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

    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 DO ! iPatch

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ausm_computeflux2_mpsd









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

SUBROUTINE rflu_ausm_computeflux2_mtcp(pRegion)

#ifdef SPEC
  USE modspecies, ONLY: t_spec_type
#endif

  USE modinterfaces, ONLY: mixtperf_c_grt, &
                           mixtperf_g_cpr, &
                           mixtperf_ho_cptuvw, & 
                           mixtperf_r_m, & 
                           mixtperf_t_dpr, & 
                           rflu_centralfirstpatch, & 
                           rflu_centralsecondpatch

  IMPLICIT NONE
  
! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,ifg,indcp,indgs,indmf,indmol,indsd,ipatch
  REAL(RFREAL) :: al,ar,cpl,cpr,dx,dy,dz,fs,gcl,gcr,gl,gr,hl,hr,irl,irr, &
                  nm,nx,ny,nz,pl,pr,rl,rr,tl,tr,ul,ur,vl,vr,wl,wr,xc,yc,zc
  REAL(RFREAL) :: flx(5),vf(3)
  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs,psd
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

#ifdef SPEC
  INTEGER :: ispec
  REAL(RFREAL) :: mml,mmr,y1,y2
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvspec 
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgcspec 
  TYPE(t_spec_type), POINTER :: pspectype
#endif  

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_AUSM_ComputeFlux2_MTCP',&
  'RFLU_ModAUSMFlux.F90')

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

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

  pgrid => pregion%grid

  indgs  = pgrid%indGs
  
  indcp  = pregion%mixtInput%indCp
  indmf  = pregion%mixtInput%indMfMixt  
  indmol = pregion%mixtInput%indMol 
  indsd  = pregion%mixtInput%indSd
  
  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  pgv  => pregion%mixt%gv
  
  pgc  => pregion%mixt%gradCell  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd  

#ifdef SPEC
  pcvspec => pregion%spec%cv  
  pgcspec => pregion%spec%gradCell    
#endif

  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN 
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! pRegion%spec%cvState  

  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_tcperf ) THEN
    CALL errorstop(global,err_gasmodel_invalid,__line__)  
  END IF ! pRegion%mixtInput%gasModel

  IF ( (indcp /= 1) .OR. (indmol /= 1) ) THEN 
    CALL errorstop(global,err_reached_default,__line__)  
  END IF ! indCp

! ******************************************************************************
! 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
    
    dx  = xc - pgrid%cofg(xcoord,c1)
    dy  = yc - pgrid%cofg(ycoord,c1)
    dz  = zc - pgrid%cofg(zcoord,c1)

#ifdef SPEC
    mml = 0.0_rfreal
    cpl = 0.0_rfreal

    DO ispec = 1,pregion%specInput%nSpecies
      pspectype => pregion%specInput%specType(ispec)

      y1 = irl*pcvspec(ispec,c1) + pgcspec(xcoord,ispec,c1)*dx & 
                                 + pgcspec(ycoord,ispec,c1)*dy & 
                                 + pgcspec(zcoord,ispec,c1)*dz

      mml = mml + y1/pspectype%pMaterial%molw
      cpl = cpl + y1*pspectype%pMaterial%spht
    END DO ! iSpec
    
    mml = 1.0_rfreal/mml
    gcl = mixtperf_r_m(mml)
    gl  = mixtperf_g_cpr(cpl,gcl)    
#endif                   

    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)
        
    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,gcl)
    al = mixtperf_c_grt(gl,gcl,tl)
    hl = mixtperf_ho_cptuvw(cpl,tl,ul,vl,wl)

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

    rr  = pcv(cv_mixt_dens,c2)
    irr = 1.0_rfreal/rr
    
    dx  = xc - pgrid%cofg(xcoord,c2)
    dy  = yc - pgrid%cofg(ycoord,c2)
    dz  = zc - pgrid%cofg(zcoord,c2)

#ifdef SPEC
    mmr = 0.0_rfreal
    cpr = 0.0_rfreal


    DO ispec = 1,pregion%specInput%nSpecies
      pspectype => pregion%specInput%specType(ispec)

      y2 = irr*pcvspec(ispec,c2) + pgcspec(xcoord,ispec,c2)*dx & 
                                 + pgcspec(ycoord,ispec,c2)*dy & 
                                 + pgcspec(zcoord,ispec,c2)*dz


      mmr = mmr + y2/pspectype%pMaterial%molw
      cpr = cpr + y2*pspectype%pMaterial%spht
    END DO ! iSpec
    
    mmr = 1.0_rfreal/mmr 
    gcr = mixtperf_r_m(mmr)
    gr  = mixtperf_g_cpr(cpr,gcr)           
#endif
              
    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)
        
    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,gcr)
    ar = mixtperf_c_grt(gr,gcr,tr)
    hr = mixtperf_ho_cptuvw(cpr,tr,ur,vr,wr)

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

    CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al,rr,ur,vr, &
                                wr,pr,hr,ar,flx,vf)

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

    pmf(indmf*ifg) = flx(1)

! ==============================================================================
!   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)
    
! ==============================================================================
!   Accumulate into substantial derivative
! ==============================================================================

    psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + vf(1)*flx(1)
    psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vf(2)*flx(1)
    psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + vf(3)*flx(1)
    
    psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - vf(1)*flx(1)
    psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vf(2)*flx(1)
    psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - vf(3)*flx(1)   
  END DO ! ifg

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

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

    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 DO ! iPatch

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ausm_computeflux2_mtcp








! ******************************************************************************
!
! Purpose: Compute convective fluxes using second-order accurate AUSM+ scheme
!   for thermally and calorically perfect gas.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: 
!  1. Only applicable to thermally and calorically perfect gas because would
!     need to recompute mixture properties at faces depending on face states.
!
! ******************************************************************************

SUBROUTINE rflu_ausm_computeflux2_tcp(pRegion)

  USE modinterfaces, ONLY: mixtperf_c_grt, &
                           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) :: al,ar,cp,dx,dy,dz,fs,gc,g,hl,hr,irl,irr,nm,nx,ny,nz,pl,pr, &
                  rl,rr,tl,tr,ul,ur,vl,vr,wl,wr,xc,yc,zc
  REAL(RFREAL) :: flx(5),vf(3)
  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_AUSM_ComputeFlux2_TCP',&
  'RFLU_ModAUSMFlux.F90')

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

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

  pgrid => pregion%grid

  indgs = pgrid%indGs
  
  indmf = pregion%mixtInput%indMfMixt  
  indsd = pregion%mixtInput%indSd
  
  pcv  => pregion%mixt%cv
  pdv  => pregion%mixt%dv
  
  pgc  => pregion%mixt%gradCell  
  pmf  => pregion%mixt%mfMixt  
  prhs => pregion%mixt%rhs
  psd  => pregion%mixt%sd  

  IF ( pregion%mixtInput%gasModel /= gas_model_tcperf ) THEN
    CALL errorstop(global,err_gasmodel_invalid,__line__)  
  END IF ! pRegion%mixtInput%gasModel

  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)
    al = mixtperf_c_grt(g,gc,tl)
    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)
    ar = mixtperf_c_grt(g,gc,tr)
    hr = mixtperf_ho_cptuvw(cp,tr,ur,vr,wr)
    
! ==============================================================================
!   Compute fluxes 
! ==============================================================================

    CALL rflu_ausm_fluxfunction(nx,ny,nz,nm,fs,rl,ul,vl,wl,pl,hl,al,rr,ur,vr, & 
                                wr,pr,hr,ar,flx,vf)

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

    pmf(indmf*ifg) = flx(1)

! ==============================================================================
!   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)
    
! ==============================================================================
!   Accumulate into substantial derivative
! ==============================================================================

    psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + vf(1)*flx(1)
    psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vf(2)*flx(1)
    psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + vf(3)*flx(1)
    
    psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - vf(1)*flx(1)
    psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vf(2)*flx(1)
    psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - vf(3)*flx(1)   
  END DO ! ifg

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

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

    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 DO ! iPatch

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ausm_computeflux2_tcp






  
! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modausmflux


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModAUSMFlux.F90,v $
! Revision 1.11  2008/12/06 08:44:20  mtcampbe
! Updated license.
!
! Revision 1.10  2008/11/19 22:17:31  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.9  2006/05/01 22:20:28  haselbac
! Removed debug statements
!
! Revision 1.8  2006/05/01 21:00:47  haselbac
! Rewrite for consistency and cleanliness
!
! Revision 1.7  2006/04/15 16:58:42  haselbac
! Added capability of running 1st order boundary fluxes with 2nd order volume 
! fluxes
!
! Revision 1.6  2006/04/07 15:19:18  haselbac
! Removed tabs
!
! Revision 1.5  2005/11/15 17:25:07  haselbac
! Bug fix: Missing declarations lead to compilation failure without SPEC=1
!
! Revision 1.4  2005/11/14 16:58:05  haselbac
! Added flux for pseudo-gas, reordered routines
!
! Revision 1.3  2005/11/10 02:26:06  haselbac
! Complete rewrite to allow for computations with variable properties
!
! Revision 1.2  2005/07/19 20:06:29  haselbac
! Modified af to prevent crashes for Skews problem
!
! Revision 1.1  2005/07/14 21:39:02  haselbac
! Initial revision
!
! ******************************************************************************