RFLU_ModNSCBC.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 routines for NCSBC boundary treatment
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModNSCBC.F90,v 1.7 2008/12/06 08:44:22 mtcampbe Exp $
!
! Copyright: (c) 2006 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modnscbc

  USE modparameters
  USE moddatatypes  
  USE modglobal, ONLY: t_global
  USE modbndpatch, ONLY: t_patch  
  USE moddatastruct, ONLY: t_region
  USE modmixture, ONLY: t_mixt_input
  USE moderror
  USE modmpi

  USE rflu_modboundconvertcv, ONLY: rflu_bxv_convertcvcons2prim, &
                                    rflu_bxv_convertcvprim2cons
  USE rflu_modconvertcv, ONLY: rflu_convertcvcons2prim, &
                               rflu_convertcvprim2cons
  USE rflu_moddifferentiationbfaces
  USE rflu_modentropyfixes  
  USE rflu_modrindstates, ONLY: rflu_setrindstatefarfieldperf
  
  USE modinterfaces, ONLY: bcondfarfieldperf, &
                           bcondinflowperf, &
                           bcondoutflowperf, &
                           mixtperf_c_dgp, & 
                           mixtperf_c_ghovm2, &
                           mixtperf_d_prt, &                                                      
                           mixtperf_eo_dgpuvw, &
                           mixtperf_g_cpr, &
                           mixtperf_ho_cptuvw, &
                           mixtperf_p_deogvm2, &
                           mixtperf_p_drt, &
                           mixtperf_r_cpg, &                           
                           mixtperf_r_m, &
                           mixtperf_t_dpr, & 
                           rflu_decideneedbgradface, &
                           rflu_decideprint, &
                           rflu_printlocinfo

  IMPLICIT NONE

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

! ==============================================================================
! Private data
! ==============================================================================

  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: RFLU_ModNSCBC.F90,v $ $Revision: 1.7 $'

! ==============================================================================
! Public functions
! ==============================================================================

  PUBLIC :: rflu_nscbc_compfirstpatchflux, &
            rflu_nscbc_comppatchflux, &
            rflu_nscbc_compsecondpatchflux, &
            rflu_nscbc_comprhs, &
            rflu_nscbc_decidehavenscbc, &
            rflu_nscbc_initff, &            
            rflu_nscbc_initiftotang, &
            rflu_nscbc_initifveltemp, &
            rflu_nscbc_initij, &             
            rflu_nscbc_initnswheat, &
            rflu_nscbc_initnswtemp, &    
            rflu_nscbc_initof, &            
            rflu_nscbc_initsw                    

! ==============================================================================
! Private functions
! ==============================================================================

  PRIVATE :: rflu_nscbc_comprhsff, &
             rflu_nscbc_comprhsiftotang, &
             rflu_nscbc_comprhsifveltemp, &             
             rflu_nscbc_comprhsij, & 
             rflu_nscbc_comprhsnswheat, &
             rflu_nscbc_comprhsnswtemp, &             
             rflu_nscbc_comprhsof, &             
             rflu_nscbc_comprhssw

! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  
  
 
 
  
  
  
! ******************************************************************************
!
! Purpose: Compute patch flux with first-order accuracy.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_compfirstpatchflux(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  LOGICAL :: decideprintflag
  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER, PARAMETER :: max_inoutflow_locs = 10
  INTEGER :: c1,ifl,indgs,nlocs,loc(max_inoutflow_locs,min_val:max_val)
  REAL(RFREAL) :: ah,cpgas,de,dissfact,dp,dq,dr,du,dv1,dv2,dv5,dw,efc,epsentr, &
                  fs,ggas,hh,hl,hr,icmassref,icpref,irl,irr,l1,l2,l5,nm,nx,ny, &
                  nz,pl,pr,pref,qh,ql,qr,rgas,rh,rl,rr,rref,sh,tl,tr,t1,t2,t3, &
                  t5,uh,ul,ur,vh,vl,vr,vref,wh,wl,wr,wt,xc,yc,zc
  REAL(RFREAL) :: flxconv(5),flxdiss(5),flx(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: mfmixt
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,rhs,vals,sd
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompFirstPatchFlux',&
  'RFLU_ModNSCBC.F90')

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

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

  cv   => pregion%mixt%cv
  dv   => pregion%mixt%dv
  gv   => pregion%mixt%gv
  rhs  => pregion%mixt%rhs
  sd   => pregion%mixt%sd

  nlocs = 0

  decideprintflag = rflu_decideprint(global)

  pref = global%refPressure
  rref = global%refDensity
  vref = global%refVelocity

  icpref    = 2.0_rfreal/(rref*vref*vref)
  icmassref = 1.0_rfreal/(rref*vref)

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

  cpgas = global%refCp
  ggas  = global%refGamma
  rgas  = mixtperf_r_cpg(cpgas,ggas)

! ******************************************************************************
! Main Body
! ******************************************************************************

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)
    nm = ppatch%fn(xyzmag,ifl)

    xc = ppatch%fc(xcoord,ifl)
    yc = ppatch%fc(ycoord,ifl)
    zc = ppatch%fc(zcoord,ifl)

    fs = ppatch%gs(indgs*ifl)

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

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

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

    ul  = cv(cv_mixt_xmom,c1)*irl
    vl  = cv(cv_mixt_ymom,c1)*irl
    wl  = cv(cv_mixt_zmom,c1)*irl
    pl  = dv(dv_mixt_pres,c1)

    tl = mixtperf_t_dpr(rl,pl,rgas)
    hl = mixtperf_ho_cptuvw(cpgas,tl,ul,vl,wl)

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

! ------------------------------------------------------------------------------
!   Right state : get right values from boundary arrays cv and dv
! ------------------------------------------------------------------------------

    rr   = ppatch%mixt%cv(cv_mixt_dens,ifl)
    irr  = 1.0_rfreal/rr

    ur   = ppatch%mixt%cv(cv_mixt_xmom,ifl)*irr
    vr   = ppatch%mixt%cv(cv_mixt_ymom,ifl)*irr
    wr   = ppatch%mixt%cv(cv_mixt_zmom,ifl)*irr
    pr   = ppatch%mixt%dv(dv_mixt_pres,ifl)

    tr = mixtperf_t_dpr(rr,pr,rgas)
    hr = mixtperf_ho_cptuvw(cpgas,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(ggas,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

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

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

    flx(1) = flxconv(1) - flxdiss(1)
    flx(2) = flxconv(2) - flxdiss(2)
    flx(3) = flxconv(3) - flxdiss(3)
    flx(4) = flxconv(4) - flxdiss(4)
    flx(5) = flxconv(5) - flxdiss(5)

    ppatch%cp(ifl)          = icpref*(0.5_rfreal*(pl + pr) - pref)
    ppatch%cmass(ifl)       = icmassref*flx(1)/nm
    ppatch%cmom(xcoord,ifl) = icpref*0.5_rfreal*(ul+ur)*flx(1)/nm
    ppatch%cmom(ycoord,ifl) = icpref*0.5_rfreal*(vl+vr)*flx(1)/nm
    ppatch%cmom(zcoord,ifl) = icpref*0.5_rfreal*(wl+wr)*flx(1)/nm

    IF ( pregion%irkStep == 1 ) THEN
      global%massIn  = global%massIn  - min(flx(1),0.0_rfreal)
      global%massOut = global%massOut + max(flx(1),0.0_rfreal)
    END IF ! pRegion

    IF ( (global%checkLevel == check_high) .AND. &
         (global%verbLevel >= verbose_high) .AND. &
         (global%myProcid == masterproc) .AND. &
         (decideprintflag .EQV. .true.) ) THEN
      IF ( ppatch%bcType == bc_outflow ) THEN
        IF ( flx(1) < 0.0_rfreal ) THEN
          nlocs = nlocs + 1

          IF ( nlocs == 1 ) THEN
            global%warnCounter = global%warnCounter + 1
  
            WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, &
                  '*** WARNING *** Inflow detected at outflow boundary!'
            WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
                                           pregion%iRegionGlobal
            IF ( global%flowType == flow_unsteady ) THEN
              WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name, &
                                                  'Current time:', &
                                                  global%currentTime
            ELSE
              WRITE(stdout,'(A,3X,A,1X,I6.6)') solver_name, &
                                               'Current iteration number:', &
                                               global%currentIter
            END IF ! global%flowType

            WRITE(stdout,'(A,3X,A,1X,I1)') solver_name, &
                                           'Runge-Kutta stage:', &
                                           pregion%irkStep
          END IF ! nLocs

          IF ( nlocs <= max_inoutflow_locs ) THEN
            loc(nlocs,min_val:max_val) = c1
          END IF ! nLocs
        END IF ! flx(1)
      END IF ! pPatch%bcType
    END IF ! global%checkLevel
  END DO ! ifl

! ------------------------------------------------------------------------------
! Write info on inflow at outflow boundary
! ------------------------------------------------------------------------------

  IF ( (global%checkLevel == check_high) .AND. &
       (global%verbLevel >= verbose_high) .AND. &
       (global%myProcid == masterproc) .AND. &
       (decideprintflag .EQV. .true.) .AND. &
       (nlocs > 0) ) THEN
    IF ( ppatch%bcType == bc_outflow ) THEN
      IF ( nlocs > max_inoutflow_locs ) THEN
        WRITE(stdout,'(A,3X,A,1X,I3,1X,A,1X,I9,1X,A)') solver_name, &
               'Only wrote the first',max_inoutflow_locs,'of',nlocs, &
               'outflow faces with inflow.'
        CALL rflu_printlocinfo(pregion,loc,max_inoutflow_locs, &
                             locinfo_mode_silent,output_mode_anybody)
      ELSE
        CALL rflu_printlocinfo(pregion,loc(1:nlocs,min_val:max_val),nlocs, &
                               locinfo_mode_silent,output_mode_anybody)
      END IF ! nLocs
    END IF ! pPatch%bcType
  END IF ! nLocs

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_compfirstpatchflux
  









 
! ******************************************************************************
!
! Purpose: Compute patch flux.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comppatchflux(pRegion,pPatch)

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

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

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

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

  LOGICAL :: decideprintflag
  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER, PARAMETER :: max_inoutflow_locs = 10
  INTEGER :: c1,ifl,indgs,nlocs,loc(max_inoutflow_locs,min_val:max_val)
  REAL(RFREAL) :: ah,cpgas,de,dissfact,dp,dq,dr,du,dv1,dv2,dv5,dw,efc,epsentr, &
                  fs,ggas,hh,hl,hr,icmassref,icpref,irl,irr,l1,l2,l5,nm,nx,ny, &
                  nz,pl,pr,pref,qh,ql,qr,rgas,rh,rl,rr,rref,sh,tl,tr,t1,t2,t3, &
                  t5,uh,ul,ur,vh,vl,vr,vref,wh,wl,wr,wt,xc,yc,zc
  REAL(RFREAL) :: flxconv(5),flxdiss(5),flx(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: mfmixt
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,rhs,vals,sd
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompPatchFlux',&
  'RFLU_ModNSCBC.F90')

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

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

  cv   => pregion%mixt%cv
  dv   => pregion%mixt%dv
  gv   => pregion%mixt%gv
  rhs  => pregion%mixt%rhs
  sd   => pregion%mixt%sd

  nlocs = 0

  decideprintflag = rflu_decideprint(global)

  pref = global%refPressure
  rref = global%refDensity
  vref = global%refVelocity

  icpref    = 2.0_rfreal/(rref*vref*vref)
  icmassref = 1.0_rfreal/(rref*vref)

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

  cpgas = global%refCp
  ggas  = global%refGamma
  rgas  = mixtperf_r_cpg(cpgas,ggas)

! ******************************************************************************
! Main Body
! ******************************************************************************

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)
    nm = ppatch%fn(xyzmag,ifl)

    xc = ppatch%fc(xcoord,ifl)
    yc = ppatch%fc(ycoord,ifl)
    zc = ppatch%fc(zcoord,ifl)

    fs = ppatch%gs(indgs*ifl)

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

! ------------------------------------------------------------------------------
!  get right values from boundary arrays cv and dv
! ------------------------------------------------------------------------------

    rr   = ppatch%mixt%cv(cv_mixt_dens,ifl)
    irr  = 1.0_rfreal/rr

    ur   = ppatch%mixt%cv(cv_mixt_xmom,ifl)*irr
    vr   = ppatch%mixt%cv(cv_mixt_ymom,ifl)*irr
    wr   = ppatch%mixt%cv(cv_mixt_zmom,ifl)*irr
    pr   = ppatch%mixt%dv(dv_mixt_pres,ifl)

    tr = mixtperf_t_dpr(rr,pr,rgas)
    hr = mixtperf_ho_cptuvw(cpgas,tr,ur,vr,wr)

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

! ------------------------------------------------------------------------------
!  get Left state also from boundary arrays cv and dv
! ------------------------------------------------------------------------------

    rl  = rr 

    ul  = ur
    vl  = vr
    wl  = wr
    pl  = pr

    tl = tr
    hl = hr

    ql = qr 

! ==============================================================================
!   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(ggas,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

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

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

    flx(1) = flxconv(1) - flxdiss(1)
    flx(2) = flxconv(2) - flxdiss(2)
    flx(3) = flxconv(3) - flxdiss(3)
    flx(4) = flxconv(4) - flxdiss(4)
    flx(5) = flxconv(5) - flxdiss(5)

    ppatch%cp(ifl)          = icpref*(0.5_rfreal*(pl + pr) - pref)
    ppatch%cmass(ifl)       = icmassref*flx(1)/nm
    ppatch%cmom(xcoord,ifl) = icpref*0.5_rfreal*(ul+ur)*flx(1)/nm
    ppatch%cmom(ycoord,ifl) = icpref*0.5_rfreal*(vl+vr)*flx(1)/nm
    ppatch%cmom(zcoord,ifl) = icpref*0.5_rfreal*(wl+wr)*flx(1)/nm

    IF ( pregion%irkStep == 1 ) THEN
      global%massIn  = global%massIn  - min(flx(1),0.0_rfreal)
      global%massOut = global%massOut + max(flx(1),0.0_rfreal)
    END IF ! pRegion

    IF ( (global%checkLevel == check_high) .AND. &
         (global%verbLevel >= verbose_high) .AND. &
         (global%myProcid == masterproc) .AND. &
         (decideprintflag .EQV. .true.) ) THEN
      IF ( ppatch%bcType == bc_outflow ) THEN
        IF ( flx(1) < 0.0_rfreal ) THEN
          nlocs = nlocs + 1

          IF ( nlocs == 1 ) THEN
            global%warnCounter = global%warnCounter + 1

            WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, &
                  '*** WARNING *** Inflow detected at outflow boundary!'
            WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
                                             pregion%iRegionGlobal

            IF ( global%flowType == flow_unsteady ) THEN
              WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name, &
                                                'Current time:', &
                                                  global%currentTime
            ELSE
              WRITE(stdout,'(A,3X,A,1X,I6.6)') solver_name, &
                                             'Current iteration number:', &
                                               global%currentIter
            END IF ! global%flowType

            WRITE(stdout,'(A,3X,A,1X,I1)') solver_name, &
                                           'Runge-Kutta stage:', &
                                           pregion%irkStep
          END IF ! nLocs

          IF ( nlocs <= max_inoutflow_locs ) THEN
            loc(nlocs,min_val:max_val) = c1
          END IF ! nLocs
        END IF ! flx(1)
      END IF ! pPatch%bcType
    END IF ! global%checkLevel
  END DO ! ifl

! ------------------------------------------------------------------------------
! Write info on inflow at outflow boundary
! ------------------------------------------------------------------------------

  IF ( (global%checkLevel == check_high) .AND. &
       (global%verbLevel >= verbose_high) .AND. &
       (global%myProcid == masterproc) .AND. &
       (decideprintflag .EQV. .true.) .AND. &
       (nlocs > 0) ) THEN
    IF ( ppatch%bcType == bc_outflow ) THEN
      IF ( nlocs > max_inoutflow_locs ) THEN
        WRITE(stdout,'(A,3X,A,1X,I3,1X,A,1X,I9,1X,A)') solver_name, &
               'Only wrote the first',max_inoutflow_locs,'of',nlocs, &
               'outflow faces with inflow.'
        CALL rflu_printlocinfo(pregion,loc,max_inoutflow_locs, &
                               locinfo_mode_silent,output_mode_anybody)
      ELSE
        CALL rflu_printlocinfo(pregion,loc(1:nlocs,min_val:max_val),nlocs, &
                               locinfo_mode_silent,output_mode_anybody)
      END IF ! nLocs
    END IF ! pPatch%bcType
  END IF ! nLocs

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comppatchflux
  







! ******************************************************************************
!
! Purpose: Compute patch fluxes with second-order accuracy.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_compsecondpatchflux(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  LOGICAL :: decideprintflag
  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER, PARAMETER :: max_inoutflow_locs = 10
  INTEGER :: c1,ifl,indgs,nlocs,loc(max_inoutflow_locs,min_val:max_val)
  REAL(RFREAL) :: ah,cpgas,de,dissfact,dp,dq,dr,dx,dy,dz,du,dv1,dv2,dv5,dw,efc, &
                  epsentr,fs,ggas,hh,hl,hr,icmassref,icpref,irl,irr,l1,l2,l5,  &
                  nm,nx,ny,nz,pl,pr,pref,qh,ql,qr,rgas,rh,rl,rr,rref,sh,tl,tr, &
                  t1,t2,t3,t5,uh,ul,ur,vh,vl,vr,vref,wh,wl,wr,wt,xc,yc,zc
  REAL(RFREAL) :: flxconv(5),flxdiss(5),flx(5)
  REAL(RFREAL), DIMENSION(:), POINTER :: mfmixt
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,rhs,vals,sd
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompSecondPatchFlux',&
  'RFLU_ModNSCBC.F90')

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

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

  cv   => pregion%mixt%cv
  dv   => pregion%mixt%dv
  gv   => pregion%mixt%gv
  grad => pregion%mixt%gradCell
  rhs  => pregion%mixt%rhs
  sd   => pregion%mixt%sd

  nlocs = 0

  decideprintflag = rflu_decideprint(global)

  pref = global%refPressure
  rref = global%refDensity
  vref = global%refVelocity

  icpref    = 2.0_rfreal/(rref*vref*vref)
  icmassref = 1.0_rfreal/(rref*vref)

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

  cpgas = global%refCp
  ggas  = global%refGamma
  rgas  = mixtperf_r_cpg(cpgas,ggas)

! ******************************************************************************
! Main Body
! ******************************************************************************

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)
    nm = ppatch%fn(xyzmag,ifl)

    xc = ppatch%fc(xcoord,ifl)
    yc = ppatch%fc(ycoord,ifl)
    zc = ppatch%fc(zcoord,ifl)

    fs = ppatch%gs(indgs*ifl)

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

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

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

    ul  = cv(cv_mixt_xmom,c1)*irl
    vl  = cv(cv_mixt_ymom,c1)*irl
    wl  = cv(cv_mixt_zmom,c1)*irl
    pl  = dv(dv_mixt_pres,c1)

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

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

    tl = mixtperf_t_dpr(rl,pl,rgas)
    hl = mixtperf_ho_cptuvw(cpgas,tl,ul,vl,wl)

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

! ------------------------------------------------------------------------------
!   Get right values from boundary arrays cv and dv
! ------------------------------------------------------------------------------

    rr   = ppatch%mixt%cv(cv_mixt_dens,ifl)
    irr  = 1.0_rfreal/rr

    ur   = ppatch%mixt%cv(cv_mixt_xmom,ifl)*irr
    vr   = ppatch%mixt%cv(cv_mixt_ymom,ifl)*irr
    wr   = ppatch%mixt%cv(cv_mixt_zmom,ifl)*irr
    pr   = ppatch%mixt%dv(dv_mixt_pres,ifl)

    tr = mixtperf_t_dpr(rr,pr,rgas)
    hr = mixtperf_ho_cptuvw(cpgas,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(ggas,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

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

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

    flx(1) = flxconv(1) - flxdiss(1)
    flx(2) = flxconv(2) - flxdiss(2)
    flx(3) = flxconv(3) - flxdiss(3)
    flx(4) = flxconv(4) - flxdiss(4)
    flx(5) = flxconv(5) - flxdiss(5)

    ppatch%cp(ifl)          = icpref*(0.5_rfreal*(pl + pr) - pref)
    ppatch%cmass(ifl)       = icmassref*flx(1)/nm
    ppatch%cmom(xcoord,ifl) = icpref*0.5_rfreal*(ul+ur)*flx(1)/nm
    ppatch%cmom(ycoord,ifl) = icpref*0.5_rfreal*(vl+vr)*flx(1)/nm
    ppatch%cmom(zcoord,ifl) = icpref*0.5_rfreal*(wl+wr)*flx(1)/nm

    IF ( pregion%irkStep == 1 ) THEN
      global%massIn  = global%massIn  - min(flx(1),0.0_rfreal)
      global%massOut = global%massOut + max(flx(1),0.0_rfreal)
    END IF ! pRegion

    IF ( (global%checkLevel == check_high) .AND. &
         (global%verbLevel >= verbose_high) .AND. &
         (global%myProcid == masterproc) .AND. &
         (decideprintflag .EQV. .true.) ) THEN
      IF ( ppatch%bcType == bc_outflow ) THEN
        IF ( flx(1) < 0.0_rfreal ) THEN
          nlocs = nlocs + 1

          IF ( nlocs == 1 ) THEN
            global%warnCounter = global%warnCounter + 1

            WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, &
                  '*** WARNING *** Reverse flow detected at global patch', &
                  ppatch%iPatchGlobal
            WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
                                             pregion%iRegionGlobal
            IF ( global%flowType == flow_unsteady ) THEN
              WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name, &
                                                  'Current time:', &
                                                  global%currentTime
            ELSE
              WRITE(stdout,'(A,3X,A,1X,I6.6)') solver_name, &
                                               'Current iteration number:', &
                                               global%currentIter
            END IF ! global%flowType

            WRITE(stdout,'(A,3X,A,1X,I1)') solver_name, &
                                           'Runge-Kutta stage:', &
                                           pregion%irkStep
          END IF ! nLocs

          IF ( nlocs <= max_inoutflow_locs ) THEN
            loc(nlocs,min_val:max_val) = c1
          END IF ! nLocs
        END IF ! flx(1)
      END IF ! pPatch%bcType
    END IF ! global%checkLevel
  END DO ! ifl

! ------------------------------------------------------------------------------
! Write info on inflow at outflow boundary
! ------------------------------------------------------------------------------

  IF ( (global%checkLevel == check_high) .AND. &
       (global%verbLevel >= verbose_high) .AND. &
       (global%myProcid == masterproc) .AND. &
       (decideprintflag .EQV. .true.) .AND. &
       (nlocs > 0) ) THEN
    IF ( ppatch%bcType == bc_outflow ) THEN
      IF ( nlocs > max_inoutflow_locs ) THEN
        WRITE(stdout,'(A,3X,A,1X,I3,1X,A,1X,I9,1X,A)') solver_name, &
               'Only wrote the first',max_inoutflow_locs,'of',nlocs, &
               'outflow faces with inflow.'
        CALL rflu_printlocinfo(pregion,loc,max_inoutflow_locs, &
                               locinfo_mode_silent,output_mode_anybody)
      ELSE
        CALL rflu_printlocinfo(pregion,loc(1:nlocs,min_val:max_val),nlocs, &
                               locinfo_mode_silent,output_mode_anybody)
      END IF ! nLocs
    END IF ! pPatch%bcType
  END IF ! nLocs

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_compsecondpatchflux




 

  

! ******************************************************************************
!
! Purpose: Wrapper for computing RHS.
!
! Description: None.
!
! Input:
!   pRegion        Region pointer
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhs(pRegion)

  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: errorflag,ipatch
  INTEGER :: varinfo(cv_mixt_xvel:cv_mixt_temp)
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhs',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Compute gradients at boundary faces
! ******************************************************************************

  CALL rflu_convertcvcons2prim(pregion,cv_mixt_state_duvwp)

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

    IF ( rflu_decideneedbgradface(pregion,ppatch) .EQV. .true. ) THEN
      CALL rflu_bxv_convertcvcons2prim(pregion,ppatch,cv_mixt_state_duvwp)

      CALL rflu_computegradbfaceswrapper(pregion,ppatch,cv_mixt_dens,cv_mixt_pres, &
                                         grbf_mixt_dens,grbf_mixt_pres, &
                                         pregion%mixt%cv,ppatch%mixt%gradFace)
  
      CALL rflu_bxv_convertcvprim2cons(pregion,ppatch,cv_mixt_state_cons)
    END IF ! RFLU_DecideNeedBGradFace
  END DO ! iPatch

  CALL rflu_convertcvprim2cons(pregion,cv_mixt_state_cons)

! ******************************************************************************
! Loop over patches and call Rhs computation routine based on type of BC
! ******************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      SELECT CASE( ppatch%bcType )
        CASE (bc_slipwall)
          CALL rflu_nscbc_comprhssw(pregion,ppatch)
        CASE (bc_noslipwall_hflux)
          CALL rflu_nscbc_comprhsnswheat(pregion,ppatch)
        CASE (bc_noslipwall_temp)
          CALL rflu_nscbc_comprhsnswtemp(pregion,ppatch)
        CASE (bc_inflow_totang)
          CALL rflu_nscbc_comprhsiftotang(pregion,ppatch)
        CASE (bc_inflow_veltemp)
          CALL rflu_nscbc_comprhsifveltemp(pregion,ppatch)
        CASE (bc_outflow)
          CALL rflu_nscbc_comprhsof(pregion,ppatch)
        CASE (bc_farfield)
          CALL rflu_nscbc_comprhsff(pregion,ppatch)
        CASE (bc_injection)
          CALL rflu_nscbc_comprhsij(pregion,ppatch)
        CASE default
          CALL errorstop(global,err_reached_default,__line__)
      END SELECT ! pPatch%bcType
    END IF ! pPatch%bcKind
  END DO ! iPatch

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhs  
  
    
  
  
  
  
  
   
  
  
  
! ******************************************************************************
!
! Purpose: Compute RHS for farfield.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsff(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  LOGICAL :: farfieldinflow,farfieldsupersonic
  INTEGER :: c1,distrib,errorflag,ifl,indcp,indmol
  REAL(RFREAL) :: a,aoa,aos,cp,d1,d2,d3,d4,d5,dotproduct,dndx,dndy,dndz,dpdn, &
                  dpds,dpdt,dpdx,dpdy,dpdz,drdn,drds,drdt,drdx,drdy,drdz,dsdx, &
                  dsdy,dsdz,dtdx,dtdy,dtdz,dudn,duds,dudt,dudx,dudy,dudz,dundn, &
                  dunds,dundt,dusdn,dusds,dusdt,dutdn,dutds,dutdt,dvdn,dvds, &
                  dvdt,dvdx,dvdy,dvdz,dwdn,dwds,dwdt,dwdx,dwdy,dwdz,dxdn,dxds, &
                  dxdt,dydn,dyds,dydt,dzdn,dzds,dzdt,g,l1,l2,l3,l4,l5,lambda1, &
                  lambda2,lambda3,lambda4,lambda5,mm,mf,mfnormal,nx,ny,nz,p,r, &
                  rgas,rhs1,rhs2,rhs3,rhs4,rhs5,ru,rv,rw,sx,sy,sz,tx,ty,tz,u, &
                  un,us,ut,v,w
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcv,pdv,prhs,pgv,vals
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradface
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsFF',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body - loop over all boundary faces
! ******************************************************************************

  vals => ppatch%mixt%vals

  distrib   = ppatch%mixt%distrib

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  pcv  => ppatch%mixt%cv
  pdv  => ppatch%mixt%dv
  pgv  => pregion%mixt%gv ! NOTE get gv from volume
  prhs => ppatch%mixt%rhs

  pgradface => ppatch%mixt%gradFace

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    cp   = pgv(gv_mixt_cp ,indcp *c1)
    mm   = pgv(gv_mixt_mol,indmol*c1)
    rgas = mixtperf_r_m(mm)
    g    = mixtperf_g_cpr(cp,rgas)

    r  = pcv(cv_mixt_dens,ifl)
    ru = pcv(cv_mixt_xmom,ifl)
    rv = pcv(cv_mixt_ymom,ifl)
    rw = pcv(cv_mixt_zmom,ifl)
    p  = pdv(dv_mixt_pres,ifl)
    a  = pdv(dv_mixt_soun,ifl)

    u = ru/r
    v = rv/r
    w = rw/r

    drdx = pgradface(xcoord,grbf_mixt_dens,ifl) 
    drdy = pgradface(ycoord,grbf_mixt_dens,ifl) 
    drdz = pgradface(zcoord,grbf_mixt_dens,ifl) 

    dudx = pgradface(xcoord,grbf_mixt_xvel,ifl) 
    dudy = pgradface(ycoord,grbf_mixt_xvel,ifl) 
    dudz = pgradface(zcoord,grbf_mixt_xvel,ifl) 

    dvdx = pgradface(xcoord,grbf_mixt_yvel,ifl) 
    dvdy = pgradface(ycoord,grbf_mixt_yvel,ifl) 
    dvdz = pgradface(zcoord,grbf_mixt_yvel,ifl) 

    dwdx = pgradface(xcoord,grbf_mixt_zvel,ifl) 
    dwdy = pgradface(ycoord,grbf_mixt_zvel,ifl) 
    dwdz = pgradface(zcoord,grbf_mixt_zvel,ifl) 

    dpdx = pgradface(xcoord,grbf_mixt_pres,ifl) 
    dpdy = pgradface(ycoord,grbf_mixt_pres,ifl) 
    dpdz = pgradface(zcoord,grbf_mixt_pres,ifl) 

! ==============================================================================
!   Computation of normal and tangential vectors
! ==============================================================================

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)    

    SELECT CASE ( pregion%mixtInput%dimens )
      CASE ( 2 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE ( 3 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens    

    dxdn = nx
    dydn = ny
    dzdn = nz

    dxds = sx
    dyds = sy
    dzds = sz

    dxdt = tx
    dydt = ty
    dzdt = tz

    dndx = dxdn
    dndy = dydn
    dndz = dzdn

    dsdx = dxds
    dsdy = dyds
    dsdz = dzds

    dtdx = dxdt
    dtdy = dydt
    dtdz = dzdt

! ==============================================================================
!   Determine if current face is inflow or outflow  AND if supersonic in normal 
!   face direction
! ==============================================================================

    aoa = vals(bcdat_farf_attack,distrib*ifl)
    aos = vals(bcdat_farf_slip  ,distrib*ifl)

    dotproduct = nx*(cos(aoa)*cos(aos)) + ny*(sin(aoa)*cos(aos)) + nz*sin(aos)

    IF ( dotproduct < 0.0_rfreal ) THEN
      farfieldinflow = .true.
    ELSE
      farfieldinflow = .false.
    END IF ! dotProduct

    mf = vals(bcdat_farf_mach  ,distrib*ifl)
    mfnormal = mf*abs(dotproduct)

    IF ( mfnormal > 1.0_rfreal ) THEN
      farfieldsupersonic = .true.
    ELSE
      farfieldsupersonic = .false.
    END IF ! dotProduct

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    drdn = drdx*dxdn + drdy*dydn + drdz*dzdn
    drds = drdx*dxds + drdy*dyds + drdz*dzds
    drdt = drdx*dxdt + drdy*dydt + drdz*dzdt

    dudn = dudx*dxdn + dudy*dydn + dudz*dzdn
    duds = dudx*dxds + dudy*dyds + dudz*dzds
    dudt = dudx*dxdt + dudy*dydt + dudz*dzdt

    dvdn = dvdx*dxdn + dvdy*dydn + dvdz*dzdn
    dvds = dvdx*dxds + dvdy*dyds + dvdz*dzds
    dvdt = dvdx*dxdt + dvdy*dydt + dvdz*dzdt

    dwdn = dwdx*dxdn + dwdy*dydn + dwdz*dzdn
    dwds = dwdx*dxds + dwdy*dyds + dwdz*dzds
    dwdt = dwdx*dxdt + dwdy*dydt + dwdz*dzdt

    dpdn = dpdx*dxdn + dpdy*dydn + dpdz*dzdn
    dpds = dpdx*dxds + dpdy*dyds + dpdz*dzds
    dpdt = dpdx*dxdt + dpdy*dydt + dpdz*dzdt

! ==============================================================================
!   Computation of normal and tangential vector components
! ==============================================================================
 
    un = u*dndx + v*dndy + w*dndz
    us = u*dsdx + v*dsdy + w*dsdz 
    ut = u*dtdx + v*dtdy + w*dtdz 

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    dundn = dudn*dndx + dvdn*dndy + dwdn*dndz
    dunds = duds*dndx + dvds*dndy + dwds*dndz
    dundt = dudt*dndx + dvdt*dndy + dwdt*dndz

    dusdn = dudn*dsdx + dvdn*dsdy + dwdn*dsdz
    dusds = duds*dsdx + dvds*dsdy + dwds*dsdz
    dusdt = dudt*dsdx + dvdt*dsdy + dwdt*dsdz

    dutdn = dudn*dtdx + dvdn*dtdy + dwdn*dtdz
    dutds = duds*dtdx + dvds*dtdy + dwds*dtdz
    dutdt = dudt*dtdx + dvdt*dtdy + dwdt*dtdz

! ==============================================================================
!   Computation of eigenvalues 
! ==============================================================================

    lambda1 = un - a
    lambda2 = un
    lambda3 = un
    lambda4 = un
    lambda5 = un + a

! ==============================================================================
!   Computations of L
! ==============================================================================

    IF ( farfieldinflow .EQV. .true. ) THEN 
      IF ( farfieldsupersonic .EQV. .true. ) THEN
        l1 = 0.0_rfreal
        l2 = 0.0_rfreal
        l3 = 0.0_rfreal
        l4 = 0.0_rfreal
        l5 = 0.0_rfreal
      ELSE
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = 0.0_rfreal
        l2 = 0.0_rfreal 
        l3 = 0.0_rfreal 
        l4 = 0.0_rfreal
      END IF ! farFieldSupersonic 
    ELSE
      IF ( farfieldsupersonic .EQV. .true. ) THEN
        l2 = lambda2*(a*a*drdn - dpdn)
        l3 = lambda3*dusdn
        l4 = lambda4*dutdn
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = lambda1*(dpdn - r*a*dundn)
      ELSE
        l2 = lambda2*(a*a*drdn - dpdn)
        l3 = lambda3*dusdn
        l4 = lambda4*dutdn
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = 0.0_rfreal 
      END IF ! farFieldSupersonic 
    END IF ! farFieldInflow

! ==============================================================================
!   Computations of d
! ==============================================================================

    d1 = (l2 + 0.5_rfreal*(l1+l5))/(a*a)
    d2 = 0.5_rfreal*(l1+l5)
    d3 = (l5-l1)/(2.0_rfreal*r*a)
    d4 = l3
    d5 = l4 

! ==============================================================================
!   Computations of Rhs
! ==============================================================================

    rhs1 = d1 + r*dusds + us*drds + r*dutdt + ut*drdt
    rhs2 = un*d1 + r*d3 + r*un*dusds + r*us*dunds + un*us*drds &
                        + r*un*dutdt + r*ut*dundt + un*ut*drdt
    rhs3 = us*d1 + r*d4 + dpds &
         + 2.0_rfreal*r*us*dusds + us*us*drds &
         + r*us*dutdt + r*ut*dusdt + us*ut*drdt
    rhs4 = ut*d1 + r*d5 + dpdt &
         + r*ut*dusds + r*us*dutds + ut*us*drds &
         + 2.0_rfreal*r*ut*dutdt + ut*ut*drdt
    rhs5 = 0.5_rfreal*(un*un+us*us+ut*ut)*d1 + d2/(g-1.0_rfreal) &
         + r*un*d3 + r*us*d4 + r*ut*d5 &
         + us*(  0.5_rfreal*(un*un+us*us+ut*ut)*drds &
              + r*(un*dunds+us*dusds+ut*dutds) &
              + g*dpds/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dusds &
         + ut*(  0.5_rfreal*(un*un+us*us+ut*ut)*drdt &
              + r*(un*dundt+us*dusdt+ut*dutdt) &
              + g*dpdt/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dutdt

    prhs(cv_mixt_dens,ifl) = rhs1 
    prhs(cv_mixt_xmom,ifl) = rhs2*dxdn + rhs3*dxds + rhs4*dxdt 
    prhs(cv_mixt_ymom,ifl) = rhs2*dydn + rhs3*dyds + rhs4*dydt 
    prhs(cv_mixt_zmom,ifl) = rhs2*dzdn + rhs3*dzds + rhs4*dzdt
    prhs(cv_mixt_ener,ifl) = rhs5 
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsff


  
  
  
  

! ******************************************************************************
!
! Purpose: Compute RHS for inflow with imposed total quantities.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsiftotang(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsIFTotAng',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body
! ******************************************************************************

  ppatch%mixt%rhs = 0.0_rfreal

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsiftotang 


  
  
  
  
  
  
  
! ******************************************************************************
!
! Purpose: Compute RHS for inflow with imposed velocity and temperature.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsifveltemp(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  INTEGER :: c1,distrib,errorflag,ifl,indcp,indmol
  REAL(RFREAL) :: a,aoa,aos,bc_dudt,bc_dvdt,bc_dwdt,bc_dtdt,bc_dundt,bc_dusdt, &
                  bc_dutdt,cp,d1,d2,d3,d4,d5,dotproduct,dndx,dndy,dndz,dpdn, &
                  dpds,dpdt,dpdx,dpdy,dpdz,drdn,drds,drdt,drdx,drdy,drdz,dsdx, &
                  dsdy,dsdz,dtdx,dtdy,dtdz,dudn,duds,dudt,dudx,dudy,dudz,dundn, &
                  dunds,dundt,dusdn,dusds,dusdt,dutdn,dutds,dutdt,dvdn,dvds, &
                  dvdt,dvdx,dvdy,dvdz,dwdn,dwds,dwdt,dwdx,dwdy,dwdz,dxdn,dxds, &
                  dxdt,dydn,dyds,dydt,dzdn,dzds,dzdt,g,l1,l2,l3,l4,l5,lambda1, &
                  lambda2,lambda3,lambda4,lambda5,mm,mf,mfnormal,nx,ny,nz,p,r, &
                  rgas,rhs1,rhs2,rhs3,rhs4,rhs5,ru,rv,rw,sx,sy,sz,t,tx,ty,tz, &
                  u,un,us,ut,v,w
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcv,pdv,prhs,pgv
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradface
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsIFVelTemp',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body - loop over all boundary faces
! ******************************************************************************

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  pcv  => ppatch%mixt%cv
  pdv  => ppatch%mixt%dv
  pgv  => pregion%mixt%gv ! NOTE get gv from volume
  prhs => ppatch%mixt%rhs

  pgradface => ppatch%mixt%gradFace

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    cp   = pgv(gv_mixt_cp ,indcp *c1)
    mm   = pgv(gv_mixt_mol,indmol*c1)
    rgas = mixtperf_r_m(mm)
    g    = mixtperf_g_cpr(cp,rgas)

    r  = pcv(cv_mixt_dens,ifl)
    ru = pcv(cv_mixt_xmom,ifl)
    rv = pcv(cv_mixt_ymom,ifl)
    rw = pcv(cv_mixt_zmom,ifl)
    p  = pdv(dv_mixt_pres,ifl)
    t  = pdv(dv_mixt_temp,ifl)
    a  = pdv(dv_mixt_soun,ifl)

    u = ru/r
    v = rv/r
    w = rw/r

    drdx = pgradface(xcoord,grbf_mixt_dens,ifl) 
    drdy = pgradface(ycoord,grbf_mixt_dens,ifl) 
    drdz = pgradface(zcoord,grbf_mixt_dens,ifl) 

    dudx = pgradface(xcoord,grbf_mixt_xvel,ifl) 
    dudy = pgradface(ycoord,grbf_mixt_xvel,ifl) 
    dudz = pgradface(zcoord,grbf_mixt_xvel,ifl) 

    dvdx = pgradface(xcoord,grbf_mixt_yvel,ifl) 
    dvdy = pgradface(ycoord,grbf_mixt_yvel,ifl) 
    dvdz = pgradface(zcoord,grbf_mixt_yvel,ifl) 

    dwdx = pgradface(xcoord,grbf_mixt_zvel,ifl) 
    dwdy = pgradface(ycoord,grbf_mixt_zvel,ifl) 
    dwdz = pgradface(zcoord,grbf_mixt_zvel,ifl) 

    dpdx = pgradface(xcoord,grbf_mixt_pres,ifl) 
    dpdy = pgradface(ycoord,grbf_mixt_pres,ifl) 
    dpdz = pgradface(zcoord,grbf_mixt_pres,ifl) 

! ==============================================================================
!   Computation of boundary condition at inflow
! ==============================================================================

    bc_dudt = 0.0_rfreal
    bc_dvdt = 0.0_rfreal
    bc_dwdt = 0.0_rfreal
    bc_dtdt = 0.0_rfreal

! ==============================================================================
!   Computation of normal and tangential vectors
! ==============================================================================

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)    

    SELECT CASE ( pregion%mixtInput%dimens )
      CASE ( 2 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE ( 3 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens    

    dxdn = nx
    dydn = ny
    dzdn = nz

    dxds = sx
    dyds = sy
    dzds = sz

    dxdt = tx
    dydt = ty
    dzdt = tz

    dndx = dxdn
    dndy = dydn
    dndz = dzdn

    dsdx = dxds
    dsdy = dyds
    dsdz = dzds

    dtdx = dxdt
    dtdy = dydt
    dtdz = dzdt

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    drdn = drdx*dxdn + drdy*dydn + drdz*dzdn
    drds = drdx*dxds + drdy*dyds + drdz*dzds
    drdt = drdx*dxdt + drdy*dydt + drdz*dzdt

    dudn = dudx*dxdn + dudy*dydn + dudz*dzdn
    duds = dudx*dxds + dudy*dyds + dudz*dzds
    dudt = dudx*dxdt + dudy*dydt + dudz*dzdt

    dvdn = dvdx*dxdn + dvdy*dydn + dvdz*dzdn
    dvds = dvdx*dxds + dvdy*dyds + dvdz*dzds
    dvdt = dvdx*dxdt + dvdy*dydt + dvdz*dzdt

    dwdn = dwdx*dxdn + dwdy*dydn + dwdz*dzdn
    dwds = dwdx*dxds + dwdy*dyds + dwdz*dzds
    dwdt = dwdx*dxdt + dwdy*dydt + dwdz*dzdt

    dpdn = dpdx*dxdn + dpdy*dydn + dpdz*dzdn
    dpds = dpdx*dxds + dpdy*dyds + dpdz*dzds
    dpdt = dpdx*dxdt + dpdy*dydt + dpdz*dzdt

! ==============================================================================
!   Computation of normal and tangential vector components
! ==============================================================================
 
    un = u*dndx + v*dndy + w*dndz
    us = u*dsdx + v*dsdy + w*dsdz 
    ut = u*dtdx + v*dtdy + w*dtdz 

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    dundn = dudn*dndx + dvdn*dndy + dwdn*dndz
    dunds = duds*dndx + dvds*dndy + dwds*dndz
    dundt = dudt*dndx + dvdt*dndy + dwdt*dndz

    dusdn = dudn*dsdx + dvdn*dsdy + dwdn*dsdz
    dusds = duds*dsdx + dvds*dsdy + dwds*dsdz
    dusdt = dudt*dsdx + dvdt*dsdy + dwdt*dsdz

    dutdn = dudn*dtdx + dvdn*dtdy + dwdn*dtdz
    dutds = duds*dtdx + dvds*dtdy + dwds*dtdz
    dutdt = dudt*dtdx + dvdt*dtdy + dwdt*dtdz

! ==============================================================================
!   Computation of boundary condition at inflow in normal and tangential direc.
! ==============================================================================

    bc_dundt = bc_dudt*dndx + bc_dvdt*dndy + bc_dwdt*dndz
    bc_dusdt = bc_dudt*dsdx + bc_dvdt*dsdy + bc_dwdt*dsdz
    bc_dutdt = bc_dudt*dtdx + bc_dvdt*dtdy + bc_dwdt*dtdz 

! ==============================================================================
!   Computation of eigenvalues 
! ==============================================================================

    lambda1 = un - a
    lambda2 = un
    lambda3 = un
    lambda4 = un
    lambda5 = un + a

! ==============================================================================
!   Computations of L
! ==============================================================================

    IF ( ppatch%mixt%switches(bcswi_inflow_type) /= bcopt_supersonic ) THEN
      IF ( ppatch%reflect == bc_reflecting ) THEN
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = l5 + 2.0_rfreal*r*a*bc_dundt
        l2 = 0.5_rfreal*(g-1.0_rfreal)*(l1+l5) + (r*a*a/t)*bc_dtdt
        l3 = -bc_dusdt 
        l4 = -bc_dutdt
      ELSE
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = 0.0_rfreal
        l2 = 0.0_rfreal 
        l3 = 0.0_rfreal 
        l4 = 0.0_rfreal
      END IF ! pPatch%reflect
    ELSE
      l1 = 0.0_rfreal
      l2 = 0.0_rfreal
      l3 = 0.0_rfreal
      l4 = 0.0_rfreal
      l5 = 0.0_rfreal
    END IF ! pPatch%mixt%switches 

! ==============================================================================
!   Computations of d
! ==============================================================================

    d1 = (l2 + 0.5_rfreal*(l1+l5))/(a*a)
    d2 = 0.5_rfreal*(l1+l5)
    d3 = (l5-l1)/(2.0_rfreal*r*a)
    d4 = l3
    d5 = l4 

! ==============================================================================
!   Computations of Rhs
! ==============================================================================

    rhs1 = d1 + r*dusds + us*drds + r*dutdt + ut*drdt
    rhs2 = un*d1 + r*d3 + r*un*dusds + r*us*dunds + un*us*drds &
                        + r*un*dutdt + r*ut*dundt + un*ut*drdt
    rhs3 = us*d1 + r*d4 + dpds &
         + 2.0_rfreal*r*us*dusds + us*us*drds &
         + r*us*dutdt + r*ut*dusdt + us*ut*drdt
    rhs4 = ut*d1 + r*d5 + dpdt &
         + r*ut*dusds + r*us*dutds + ut*us*drds &
         + 2.0_rfreal*r*ut*dutdt + ut*ut*drdt
    rhs5 = 0.5_rfreal*(un*un+us*us+ut*ut)*d1 + d2/(g-1.0_rfreal) &
         + r*un*d3 + r*us*d4 + r*ut*d5 &
         + us*(  0.5_rfreal*(un*un+us*us+ut*ut)*drds &
              + r*(un*dunds+us*dusds+ut*dutds) &
              + g*dpds/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dusds &
         + ut*(  0.5_rfreal*(un*un+us*us+ut*ut)*drdt &
              + r*(un*dundt+us*dusdt+ut*dutdt) &
              + g*dpdt/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dutdt

    prhs(cv_mixt_dens,ifl) = rhs1 
    prhs(cv_mixt_xmom,ifl) = rhs2*dxdn + rhs3*dxds + rhs4*dxdt 
    prhs(cv_mixt_ymom,ifl) = rhs2*dydn + rhs3*dyds + rhs4*dydt 
    prhs(cv_mixt_zmom,ifl) = rhs2*dzdn + rhs3*dzds + rhs4*dzdt
    prhs(cv_mixt_ener,ifl) = rhs5 

    IF ( ppatch%mixt%switches(bcswi_inflow_type) /= bcopt_supersonic ) THEN
      IF ( ppatch%reflect == bc_reflecting ) THEN
        prhs(cv_mixt_xmom,ifl) = 0.0_rfreal 
        prhs(cv_mixt_ymom,ifl) = 0.0_rfreal 
        prhs(cv_mixt_zmom,ifl) = 0.0_rfreal 
        prhs(cv_mixt_ener,ifl) = 0.0_rfreal 
      END IF ! pPatch%reflect
    END IF ! pPatch%mixt%switches
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsifveltemp 
   
  
  
  
  
  
  
  
! ******************************************************************************
!
! Purpose: Compute RHS for injection boundary.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsij(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsIJ',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body
! ******************************************************************************

  ppatch%mixt%rhs = 0.0_rfreal

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsij


  
  
  


! ******************************************************************************
!
! Purpose: Compute RHS for no-slip wall with imposed heat flux.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsnswheat(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsNSWHeat',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body
! ******************************************************************************

  ppatch%mixt%rhs = 0.0_rfreal

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsnswheat 


  
  
  
  
  
! ******************************************************************************
!
! Purpose: Compute RHS for no-slip wall with imposed temperature.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsnswtemp(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsNSWTemp',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body
! ******************************************************************************

  ppatch%mixt%rhs = 0.0_rfreal

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsnswtemp   
  
  
  
  
  
  

! ******************************************************************************
!
! Purpose: Compute RHS for outflow.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhsof(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  INTEGER :: bcopttype,c1,distrib,errorflag,ifl,indcp,indmol
  REAL(RFREAL) :: a,aoa,aos,cp,d1,d2,d3,d4,d5,dotproduct,dndx,dndy,dndz,dpdn, &
                  dpds,dpdt,dpdx,dpdy,dpdz,drdn,drds,drdt,drdx,drdy,drdz,dsdx, &
                  dsdy,dsdz,dtdx,dtdy,dtdz,dudn,duds,dudt,dudx,dudy,dudz,dundn, &
                  dunds,dundt,dusdn,dusds,dusdt,dutdn,dutds,dutdt,dvdn,dvds, &
                  dvdt,dvdx,dvdy,dvdz,dwdn,dwds,dwdt,dwdx,dwdy,dwdz,dxdn,dxds, &
                  dxdt,dydn,dyds,dydt,dzdn,dzds,dzdt,g,l1,l2,l3,l4,l5,lambda1, &
                  lambda2,lambda3,lambda4,lambda5,mm,mf,mfnormal,nscbck,nx,ny, &
                  nz,p,pinf,r,rgas,rhs1,rhs2,rhs3,rhs4,rhs5,ru,rv,rw,sx,sy,sz, &
                  t,tx,ty,tz,u,un,us,ut,v,w
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcv,pdv,prhs,pgv,vals
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradface
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsOF',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body - loop over all boundary faces
! ******************************************************************************

  distrib   = ppatch%mixt%distrib
  bcopttype = ppatch%mixt%switches(bcswi_outflow_type)

  IF ( bcopttype /= bcopt_supersonic ) THEN
    vals => ppatch%mixt%vals
  END IF ! bcOptType

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  pcv  => ppatch%mixt%cv
  pdv  => ppatch%mixt%dv
  pgv  => pregion%mixt%gv ! NOTE get gv from volume
  prhs => ppatch%mixt%rhs

  pgradface => ppatch%mixt%gradFace

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    cp   = pgv(gv_mixt_cp ,indcp *c1)
    mm   = pgv(gv_mixt_mol,indmol*c1)
    rgas = mixtperf_r_m(mm)
    g    = mixtperf_g_cpr(cp,rgas)

    r  = pcv(cv_mixt_dens,ifl)
    ru = pcv(cv_mixt_xmom,ifl)
    rv = pcv(cv_mixt_ymom,ifl)
    rw = pcv(cv_mixt_zmom,ifl)
    p  = pdv(dv_mixt_pres,ifl)
    a  = pdv(dv_mixt_soun,ifl)

    u = ru/r
    v = rv/r
    w = rw/r

    drdx = pgradface(xcoord,grbf_mixt_dens,ifl) 
    drdy = pgradface(ycoord,grbf_mixt_dens,ifl) 
    drdz = pgradface(zcoord,grbf_mixt_dens,ifl) 

    dudx = pgradface(xcoord,grbf_mixt_xvel,ifl) 
    dudy = pgradface(ycoord,grbf_mixt_xvel,ifl) 
    dudz = pgradface(zcoord,grbf_mixt_xvel,ifl) 

    dvdx = pgradface(xcoord,grbf_mixt_yvel,ifl) 
    dvdy = pgradface(ycoord,grbf_mixt_yvel,ifl) 
    dvdz = pgradface(zcoord,grbf_mixt_yvel,ifl) 

    dwdx = pgradface(xcoord,grbf_mixt_zvel,ifl) 
    dwdy = pgradface(ycoord,grbf_mixt_zvel,ifl) 
    dwdz = pgradface(zcoord,grbf_mixt_zvel,ifl) 

    dpdx = pgradface(xcoord,grbf_mixt_pres,ifl) 
    dpdy = pgradface(ycoord,grbf_mixt_pres,ifl) 
    dpdz = pgradface(zcoord,grbf_mixt_pres,ifl) 

! ==============================================================================
!   Computation of normal and tangential vectors
! ==============================================================================

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)    

    SELECT CASE ( pregion%mixtInput%dimens )
      CASE ( 2 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE ( 3 )
        sx = -ny
        sy =  nx
        sz =  nz    

        tx = 0.0_rfreal
        ty = 0.0_rfreal
        tz = 1.0_rfreal    
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens    

    dxdn = nx
    dydn = ny
    dzdn = nz

    dxds = sx
    dyds = sy
    dzds = sz

    dxdt = tx
    dydt = ty
    dzdt = tz

    dndx = dxdn
    dndy = dydn
    dndz = dzdn

    dsdx = dxds
    dsdy = dyds
    dsdz = dzds

    dtdx = dxdt
    dtdy = dydt
    dtdz = dzdt

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    drdn = drdx*dxdn + drdy*dydn + drdz*dzdn
    drds = drdx*dxds + drdy*dyds + drdz*dzds
    drdt = drdx*dxdt + drdy*dydt + drdz*dzdt

    dudn = dudx*dxdn + dudy*dydn + dudz*dzdn
    duds = dudx*dxds + dudy*dyds + dudz*dzds
    dudt = dudx*dxdt + dudy*dydt + dudz*dzdt

    dvdn = dvdx*dxdn + dvdy*dydn + dvdz*dzdn
    dvds = dvdx*dxds + dvdy*dyds + dvdz*dzds
    dvdt = dvdx*dxdt + dvdy*dydt + dvdz*dzdt

    dwdn = dwdx*dxdn + dwdy*dydn + dwdz*dzdn
    dwds = dwdx*dxds + dwdy*dyds + dwdz*dzds
    dwdt = dwdx*dxdt + dwdy*dydt + dwdz*dzdt

    dpdn = dpdx*dxdn + dpdy*dydn + dpdz*dzdn
    dpds = dpdx*dxds + dpdy*dyds + dpdz*dzds
    dpdt = dpdx*dxdt + dpdy*dydt + dpdz*dzdt

! ==============================================================================
!   Computation of normal and tangential vector components
! ==============================================================================
 
    un = u*dndx + v*dndy + w*dndz
    us = u*dsdx + v*dsdy + w*dsdz 
    ut = u*dtdx + v*dtdy + w*dtdz 

! ==============================================================================
!   Computation of normal and tangential derivatives
! ==============================================================================

    dundn = dudn*dndx + dvdn*dndy + dwdn*dndz
    dunds = duds*dndx + dvds*dndy + dwds*dndz
    dundt = dudt*dndx + dvdt*dndy + dwdt*dndz

    dusdn = dudn*dsdx + dvdn*dsdy + dwdn*dsdz
    dusds = duds*dsdx + dvds*dsdy + dwds*dsdz
    dusdt = dudt*dsdx + dvdt*dsdy + dwdt*dsdz

    dutdn = dudn*dtdx + dvdn*dtdy + dwdn*dtdz
    dutds = duds*dtdx + dvds*dtdy + dwds*dtdz
    dutdt = dudt*dtdx + dvdt*dtdy + dwdt*dtdz

! ==============================================================================
!   Computation of eigenvalues 
! ==============================================================================

    nscbck = ppatch%nscbcK

    IF ( bcopttype /= bcopt_supersonic ) THEN
      pinf = vals(bcdat_outflow_press,distrib*ifl)
    ELSE
      pinf = p
    END IF ! bcOptType

    lambda1 = un - a
    lambda2 = un
    lambda3 = un
    lambda4 = un
    lambda5 = un + a

! ==============================================================================
!   Computations of L
! ==============================================================================

    IF ( ppatch%mixt%switches(bcswi_outflow_type) == bcopt_supersonic ) THEN
      l2 = lambda2*(a*a*drdn - dpdn)
      l3 = lambda3*dusdn
      l4 = lambda4*dutdn
      l5 = lambda5*(dpdn + r*a*dundn)
      l1 = lambda1*(dpdn - r*a*dundn)
    ELSE
      IF ( ppatch%reflect == bc_reflecting ) THEN
        l2 = lambda2*(a*a*drdn - dpdn)
        l3 = lambda3*dusdn
        l4 = lambda4*dutdn
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = -l5
      ELSE ! non reflecting BC
        l2 = lambda2*(a*a*drdn - dpdn)
        l3 = lambda3*dusdn
        l4 = lambda4*dutdn
        l5 = lambda5*(dpdn + r*a*dundn)
        l1 = nscbck*(p-pinf)
      END IF ! pPatch%reflect
    END IF ! pPatch%mixt%switches

! ==============================================================================
!   Computations of d
! ==============================================================================

    d1 = (l2 + 0.5_rfreal*(l1+l5))/(a*a)
    d2 = 0.5_rfreal*(l1+l5)
    d3 = (l5-l1)/(2.0_rfreal*r*a)
    d4 = l3
    d5 = l4 

! ==============================================================================
!   Computations of Rhs
! ==============================================================================

    rhs1 = d1 + r*dusds + us*drds + r*dutdt + ut*drdt
    rhs2 = un*d1 + r*d3 + r*un*dusds + r*us*dunds + un*us*drds &
                        + r*un*dutdt + r*ut*dundt + un*ut*drdt
    rhs3 = us*d1 + r*d4 + dpds &
         + 2.0_rfreal*r*us*dusds + us*us*drds &
         + r*us*dutdt + r*ut*dusdt + us*ut*drdt
    rhs4 = ut*d1 + r*d5 + dpdt &
         + r*ut*dusds + r*us*dutds + ut*us*drds &
         + 2.0_rfreal*r*ut*dutdt + ut*ut*drdt
    rhs5 = 0.5_rfreal*(un*un+us*us+ut*ut)*d1 + d2/(g-1.0_rfreal) &
         + r*un*d3 + r*us*d4 + r*ut*d5 &
         + us*(  0.5_rfreal*(un*un+us*us+ut*ut)*drds &
              + r*(un*dunds+us*dusds+ut*dutds) &
              + g*dpds/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dusds &
         + ut*(  0.5_rfreal*(un*un+us*us+ut*ut)*drdt &
              + r*(un*dundt+us*dusdt+ut*dutdt) &
              + g*dpdt/(g-1.0_rfreal) ) &
         + (0.5_rfreal*r*(un*un+us*us+ut*ut) &
              + g*p/(g-1.0_rfreal))*dutdt

    prhs(cv_mixt_dens,ifl) = rhs1 
    prhs(cv_mixt_xmom,ifl) = rhs2*dxdn + rhs3*dxds + rhs4*dxdt 
    prhs(cv_mixt_ymom,ifl) = rhs2*dydn + rhs3*dyds + rhs4*dydt 
    prhs(cv_mixt_zmom,ifl) = rhs2*dzdn + rhs3*dzds + rhs4*dzdt

    IF ( ppatch%mixt%switches(bcswi_outflow_type) == bcopt_supersonic ) THEN
      prhs(cv_mixt_ener,ifl) = rhs5 
    ELSE IF ( ppatch%mixt%switches(bcswi_outflow_type) == bcopt_subsonic ) THEN
      IF ( ppatch%reflect == bc_reflecting ) THEN
        prhs(cv_mixt_ener,ifl) = 0.0_rfreal
      ELSE ! non reflecting BC
        prhs(cv_mixt_ener,ifl) = rhs5 
      END IF ! pPatch%reflect
    END IF ! pPatch%mixt%switches
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhsof
 

  
  

  
  
  



! ******************************************************************************
!
! Purpose: Compute RHS for slip wall.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_comprhssw(pRegion,pPatch) 

  IMPLICIT NONE

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

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

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

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

  INTEGER :: bcopttype,c1,distrib,errorflag,ifl,indcp,indmol
  REAL(RFREAL) :: a,aoa,aos,cp,d1,d2,d3,d4,d5,dotproduct,dndx,dndy,dndz,dpdn, &
                  dpds,dpdt,dpdx,dpdy,dpdz,drdn,drds,drdt,drdx,drdy,drdz,dsdx, &
                  dsdy,dsdz,dtdx,dtdy,dtdz,dudn,duds,dudt,dudx,dudy,dudz,dundn, &
                  dunds,dundt,dusdn,dusds,dusdt,dutdn,dutds,dutdt,dvdn,dvds, &
                  dvdt,dvdx,dvdy,dvdz,dwdn,dwds,dwdt,dwdx,dwdy,dwdz,dxdn,dxds, &
                  dxdt,dydn,dyds,dydt,dzdn,dzds,dzdt,g,l1,l2,l3,l4,l5,lambda1, &
                  lambda2,lambda3,lambda4,lambda5,mm,mf,mfnormal,nscbck,nx,ny, &
                  nz,p,pinf,r,rgas,rhs1,rhs2,rhs3,rhs4,rhs5,ru,rv,rw,sx,sy,sz, &
                  t,tx,ty,tz,u,un,us,ut,v,w
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcv,pdv,prhs,pgv,vals
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradface
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_CompRhsSW',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body - loop over all boundary faces
! ******************************************************************************

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  pcv  => ppatch%mixt%cv
  pdv  => ppatch%mixt%dv
  pgv  => pregion%mixt%gv ! NOTE get gv from volume
  prhs => ppatch%mixt%rhs

  pgradface => ppatch%mixt%gradFace

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    cp   = pgv(gv_mixt_cp ,indcp *c1)
    mm   = pgv(gv_mixt_mol,indmol*c1)
    rgas = mixtperf_r_m(mm)
    g    = mixtperf_g_cpr(cp,rgas)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)

    r  = pcv(cv_mixt_dens,ifl)
    ru = pcv(cv_mixt_xmom,ifl)
    rv = pcv(cv_mixt_ymom,ifl)
    rw = pcv(cv_mixt_zmom,ifl)
    p  = pdv(dv_mixt_pres,ifl)
    a  = pdv(dv_mixt_soun,ifl)

    u = ru/r
    v = rv/r
    w = rw/r

    drdx = pgradface(xcoord,grbf_mixt_dens,ifl)
    drdy = pgradface(ycoord,grbf_mixt_dens,ifl)
    drdz = pgradface(zcoord,grbf_mixt_dens,ifl)

    dudx = pgradface(xcoord,grbf_mixt_xvel,ifl)
    dudy = pgradface(ycoord,grbf_mixt_xvel,ifl)
    dudz = pgradface(zcoord,grbf_mixt_xvel,ifl)

    dvdx = pgradface(xcoord,grbf_mixt_yvel,ifl)
    dvdy = pgradface(ycoord,grbf_mixt_yvel,ifl)
    dvdz = pgradface(zcoord,grbf_mixt_yvel,ifl)

    dwdx = pgradface(xcoord,grbf_mixt_zvel,ifl)
    dwdy = pgradface(ycoord,grbf_mixt_zvel,ifl)
    dwdz = pgradface(zcoord,grbf_mixt_zvel,ifl)

    dpdx = pgradface(xcoord,grbf_mixt_pres,ifl)
    dpdy = pgradface(ycoord,grbf_mixt_pres,ifl)
    dpdz = pgradface(zcoord,grbf_mixt_pres,ifl)

! ASSUMPTION : only horizontal slip walls are there
! slip wall : normal component of velocity is zero
! TEMPORARY : make sure this is the case
    v = 0.0_rfreal 

    lambda1 = v - a
    lambda2 = v
    lambda3 = v
    lambda4 = v
    lambda5 = v + a

! ==============================================================================
!   Computations of L
! ==============================================================================

    l2 = lambda2*(a*a*drdy - dpdy)
    l3 = lambda3*dudy
    l4 = lambda4*dwdy

! reflecting slip wall 
!    IF ( ny > 0.0_RFREAL ) THEN ! upper wall
      l5 = lambda5*(dpdy + r*a*dvdy)
      l1 = l5
!    ELSE ! lower wall
      l1 = lambda1*(dpdy - r*a*dvdy)
      l5 = l1
!    ENF IF

! ==============================================================================
!   Computations of d
! ==============================================================================

    d1 = (l2 + 0.5_rfreal*(l1+l5))/(a*a)
    d2 = 0.5_rfreal*(l1+l5)
    d3 = l3
    d4 = (l5-l1)/(2.0_rfreal*r*a)
    d5 = l4

! ==============================================================================
!   Computations of Rhs
! ==============================================================================

    prhs(cv_mixt_dens,ifl) = r*dudx + u*drdx + d1 + r*dwdz + w*drdz
    prhs(cv_mixt_xmom,ifl) = u*d1 + r*d3 + dpdx &
                           + 2.0_rfreal*r*u*dudx + u*u*drdx &
                           + r*u*dwdz + r*w*dudz + u*w*drdz

! No need to solve y-momentum equation
! y momentum = 0, always
    prhs(cv_mixt_ymom,ifl) = 0.0_rfreal 

    prhs(cv_mixt_zmom,ifl) = w*d1 + r*d5 + dpdz &
                           + r*w*dudx + r*u*dwdx + w*u*drdx &
                           + 2.0_rfreal*r*w*dwdz + w*w*drdz

    prhs(cv_mixt_ener,ifl) = 0.5_rfreal*(u*u+v*v+w*w)*d1 + d2/(g-1.0_rfreal) &
                     + r*u*d3 + r*v*d4 + r*w*d5 &
                     + u*(  0.5_rfreal*(u*u+v*v+w*w)*drdx &
                          + r*(u*dudx+v*dvdx+w*dwdx) &
                          + g*dpdx/(g-1.0_rfreal) ) &
                     + (0.5_rfreal*r*(u*u+v*v+w*w)+g*p/(g-1.0_rfreal))*dudx &
                     + w*(  0.5_rfreal*(u*u+v*v+w*w)*drdz &
                          + r*(u*dudz+v*dvdz+w*dwdz) &
                          + g*dpdz/(g-1.0_rfreal) ) &
                     + (0.5_rfreal*r*(u*u+v*v+w*w)+g*p/(g-1.0_rfreal))*dwdz

  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_comprhssw 
  
  
  
  
  
  
 


! ******************************************************************************
!
! Purpose: Determine whether need boundary face gradients and related functions.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region data
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

LOGICAL FUNCTION rflu_nscbc_decidehavenscbc(pRegion)

  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: ipatch
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_DecideHaveNSCBC',&
  'RFLU_ModNSCBC.F90')

! *****************************************************************************
! Initialize
! *****************************************************************************

  rflu_nscbc_decidehavenscbc = .false.

! *****************************************************************************
! Determine whether need boundary face gradients
! *****************************************************************************

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

    IF ( ppatch%bcKind == bc_kind_nscbc ) THEN
      rflu_nscbc_decidehavenscbc = .true.
  
      EXIT patchloop
    END IF ! pPatch%bcKind
  END DO patchloop

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

  CALL deregisterfunction(global)

END FUNCTION rflu_nscbc_decidehavenscbc







! ******************************************************************************
!
! Purpose: Initialize farfield boundary.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initff(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  LOGICAL :: corrflag
  INTEGER :: errorflag,ifl,c1,bcopttype,distrib,gasmodel,indcp,indgs,indmol
  REAL(RFREAL) :: aoa,aos,corr,cp,irl,eo,g,gc,liftcoef,mf,mw,nm,nx,ny,nz,pf, &
                  pl,pr,rel,rer,rl,rr,rul,rur,rvl,rvr,rwl,rwr,tf,ul,vl,wl,xc, &
                  yc,zc
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,vals
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitFF',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body 
! ******************************************************************************

  corr = ppatch%mixt%switches(bcswi_farf_corr)

  vals => ppatch%mixt%vals
      
  IF ( corr == bcopt_corr_yes ) THEN 
    corrflag = .true. 
  ELSE 
    corrflag = .false.
    liftcoef = 0.0_rfreal
  END IF ! corr

  indcp    = pregion%mixtInput%indCp
  indmol   = pregion%mixtInput%indMol
  gasmodel = pregion%mixtInput%gasModel

  distrib  =  ppatch%mixt%distrib

  ppatch%mixt%cvState = cv_mixt_state_cons

  IF (pregion%mixt%cvState /= cv_mixt_state_cons ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__, &
                   'region mixt cvState is incompatible...')
  END IF ! pRegion%mixt%cvState

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)
    nm = ppatch%fn(xyzmag,ifl)

    xc = ppatch%fc(xcoord,ifl)
    yc = ppatch%fc(ycoord,ifl)
    zc = ppatch%fc(zcoord,ifl)

    rl  = pregion%mixt%cv(cv_mixt_dens,c1)
    irl = 1.0_rfreal/rl
    rul = pregion%mixt%cv(cv_mixt_xmom,c1)*irl
    rvl = pregion%mixt%cv(cv_mixt_ymom,c1)*irl
    rwl = pregion%mixt%cv(cv_mixt_zmom,c1)*irl
    rel = pregion%mixt%cv(cv_mixt_ener,c1)*irl

    pl  = pregion%mixt%dv(dv_mixt_pres,c1)

    mf  = vals(bcdat_farf_mach  ,distrib*ifl)
    aoa = vals(bcdat_farf_attack,distrib*ifl)
    aos = vals(bcdat_farf_slip  ,distrib*ifl)
    pf  = vals(bcdat_farf_press ,distrib*ifl)
    tf  = vals(bcdat_farf_temp  ,distrib*ifl)

    IF ( gasmodel == gas_model_tcperf ) THEN

      mw = pregion%mixt%gv(gv_mixt_mol,indmol*c1)
      cp = pregion%mixt%gv(gv_mixt_cp ,indcp *c1)
      gc = mixtperf_r_m(mw)
      g  = mixtperf_g_cpr(cp,gc)

      rel = rl*mixtperf_eo_dgpuvw(rl,g,pl,ul,vl,wl)
      
      CALL rflu_setrindstatefarfieldperf(global,cp,mw,nx,ny,nz,mf,pf,tf, &
                                         aoa,aos,corrflag,liftcoef,xc,yc, &
                                         zc,rl,rul,rvl,rwl,rel,rr,rur, & 
                                         rvr,rwr,rer,pr)     
    ELSE
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! gasModel

    ppatch%mixt%cv(cv_mixt_dens,ifl) = rr 
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = rur
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = rvr
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = rwr
    ppatch%mixt%cv(cv_mixt_ener,ifl) = rer
    
! TEMPORARY : overwriting above computation of farfield BC data by cell data
!             this is done to avoid waves at farfield just in begining of run
    ppatch%mixt%cv(cv_mixt_dens,ifl) = pregion%mixt%cv(cv_mixt_dens,c1)
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = pregion%mixt%cv(cv_mixt_xmom,c1)
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = pregion%mixt%cv(cv_mixt_ymom,c1)
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = pregion%mixt%cv(cv_mixt_zmom,c1)
    ppatch%mixt%cv(cv_mixt_ener,ifl) = pregion%mixt%cv(cv_mixt_ener,c1) 
! END TEMPORARY
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initff







! ******************************************************************************
!
! Purpose: Initialize inflow boundary with imposed total quantities.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initiftotang(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  INTEGER :: errorflag,ifl,c1
  INTEGER :: bcoptfixed,bcopttype,distrib,gasmodel,indcp,indgs,indmol
  REAL(RFREAL) :: cp,betah,betav,g,gc,mach,mw,nx,ny,nz,nm,pr,ptot,rer,rl,rr, &
                  rul,rvl,rwl,rur,rvr,rwr,ttot
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,vals
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitIFTotAng',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body 
! ******************************************************************************

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  distrib    = ppatch%mixt%distrib
  bcopttype  = ppatch%mixt%switches(bcswi_inflow_type)
  bcoptfixed = ppatch%mixt%switches(bcswi_inflow_fixed)

  vals => ppatch%mixt%vals

  ppatch%mixt%cvState = cv_mixt_state_cons

  IF (pregion%mixt%cvState /= cv_mixt_state_cons ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__, &
                   'region mixt cvState is incompatible...')
  END IF ! pRegion%mixt%cvState

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    mw = pregion%mixt%gv(gv_mixt_mol,indmol*c1)
    cp = pregion%mixt%gv(gv_mixt_cp ,indcp *c1)

    gc = mixtperf_r_m(mw)
    g  = mixtperf_g_cpr(cp,gc)

    nx = ppatch%fn(xcoord,ifl)
    ny = ppatch%fn(ycoord,ifl)
    nz = ppatch%fn(zcoord,ifl)
    nm = ppatch%fn(xyzmag,ifl)

    rl  = pregion%mixt%cv(cv_mixt_dens,c1)
    rul = pregion%mixt%cv(cv_mixt_xmom,c1)
    rvl = pregion%mixt%cv(cv_mixt_ymom,c1)
    rwl = pregion%mixt%cv(cv_mixt_zmom,c1)

    ptot  = vals(bcdat_inflow_ptot, distrib*ifl)
    ttot  = vals(bcdat_inflow_ttot, distrib*ifl)
    betah = vals(bcdat_inflow_betah,distrib*ifl)
    betav = vals(bcdat_inflow_betav,distrib*ifl)

    IF ( bcopttype == bcopt_supersonic ) THEN
      mach = vals(bcdat_inflow_mach,distrib*ifl)
    ELSE
      mach = 0.0_rfreal
    END IF ! bcOptType

    CALL bcondinflowperf(bcopttype,bcoptfixed,ptot,ttot,betah,betav, &
                         mach,nx,ny,nz,cp,mw,rl,rul,rvl,rwl,rr,rur, &
                         rvr,rwr,rer,pr)

    ppatch%mixt%cv(cv_mixt_dens,ifl) = rr 
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = rur
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = rvr
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = rwr
    ppatch%mixt%cv(cv_mixt_ener,ifl) = rer
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initiftotang





! ******************************************************************************
!
! Purpose: Initialize inflow with imposed velocity and temperature.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initifveltemp(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  INTEGER :: errorflag,ifl,c1,bcopttype,distrib,gasmodel,indcp,indgs,indmol
  REAL(RFREAL) :: cp,eo,g,gc,mw,pr,rl,rr,rur,rvr,rwr,tr,ur,vr,wr
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,vals
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitIFVelTemp',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body 
! ******************************************************************************

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  distrib   = ppatch%mixt%distrib
  bcopttype = ppatch%mixt%switches(bcswi_inflow_type)

  vals => ppatch%mixt%vals

  ppatch%mixt%cvState = cv_mixt_state_cons

  IF (pregion%mixt%cvState /= cv_mixt_state_cons ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__, &
                   'region mixt cvState is incompatible...')
  END IF ! pRegion%mixt%cvState

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    mw = pregion%mixt%gv(gv_mixt_mol,indmol*c1)
    cp = pregion%mixt%gv(gv_mixt_cp ,indcp *c1)

    gc = mixtperf_r_m(mw)
    g  = mixtperf_g_cpr(cp,gc)

    rl  = pregion%mixt%cv(cv_mixt_dens,c1)

    ur = vals(bcdat_inflow_u,distrib*ifl)
    vr = vals(bcdat_inflow_v,distrib*ifl)
    wr = vals(bcdat_inflow_w,distrib*ifl)
    tr = vals(bcdat_inflow_t,distrib*ifl)

    IF ( bcopttype == bcopt_supersonic ) THEN
      pr = vals(bcdat_inflow_p,distrib*ifl)
      rr = mixtperf_d_prt(pr,gc,tr)
    ELSE
      rr = rl
      pr = mixtperf_p_drt(rr,gc,tr)
    END IF ! bcOptType

    rur = rr*ur
    rvr = rr*vr
    rwr = rr*wr

    ppatch%mixt%cv(cv_mixt_dens,ifl) = rr 
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = rur
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = rvr
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = rwr

    eo = mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)

    ppatch%mixt%cv(cv_mixt_ener,ifl) = rr*eo
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initifveltemp







! ******************************************************************************
!
! Purpose: Initialize injection boundary.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initij(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitIJ',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body 
! ******************************************************************************

  ppatch%mixt%cv = 0.0_rfreal

  ppatch%mixt%cvState = cv_mixt_state_cons

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initij







! ******************************************************************************
!
! Purpose: Initialize no-slip wall with imposed heat flux.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initnswheat(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitNSWHeat',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body 
! ******************************************************************************

  ppatch%mixt%cv = 0.0_rfreal

  ppatch%mixt%cvState = cv_mixt_state_cons

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initnswheat






! ******************************************************************************
!
! Purpose: Initialize no-slip wall with imposed temperature.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initnswtemp(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitNSWTemp',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main body 
! ******************************************************************************

  ppatch%mixt%cv = 0.0_rfreal

  ppatch%mixt%cvState = cv_mixt_state_cons

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initnswtemp






! ******************************************************************************
!
! Purpose: Initialize outflow.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initof(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  INTEGER :: errorflag,ifl,c1,bcopttype,distrib,gasmodel,indcp,indgs,indmol
  REAL(RFREAL) :: cp,eo,g,gc,mw,pl,pr,rl,rul,rvl,rwl,rel,ul,vl,wl
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,vals
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitOF',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body 
! ******************************************************************************

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  distrib =  ppatch%mixt%distrib
  bcopttype = ppatch%mixt%switches(bcswi_outflow_type)

  IF ( bcopttype /= bcopt_supersonic ) THEN
    vals => ppatch%mixt%vals
  END IF ! bcOptType

  ppatch%mixt%cvState = cv_mixt_state_cons

  IF (pregion%mixt%cvState /= cv_mixt_state_cons ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__, &
                   'region mixt cvState is incompatible...')
  END IF ! pRegion%mixt%cvState

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    mw = pregion%mixt%gv(gv_mixt_mol,indmol*c1)
    cp = pregion%mixt%gv(gv_mixt_cp ,indcp *c1)

    gc = mixtperf_r_m(mw)
    g  = mixtperf_g_cpr(cp,gc)

    rl  = pregion%mixt%cv(cv_mixt_dens,c1)
    rul = pregion%mixt%cv(cv_mixt_xmom,c1)
    rvl = pregion%mixt%cv(cv_mixt_ymom,c1)
    rwl = pregion%mixt%cv(cv_mixt_zmom,c1)
    rel = pregion%mixt%cv(cv_mixt_ener,c1)

    ul = rul/rl
    vl = rvl/rl
    wl = rwl/rl

    pl  = pregion%mixt%dv(dv_mixt_pres,c1)

    IF ( bcopttype /= bcopt_supersonic ) THEN
      pr = vals(bcdat_outflow_press,distrib*ifl)
    ELSE
      pr = pl
    END IF ! bcOptType

    ppatch%mixt%cv(cv_mixt_dens,ifl) = rl 
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = rul
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = rvl
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = rwl

    eo = mixtperf_eo_dgpuvw(rl,g,pr,ul,vl,wl)

    ppatch%mixt%cv(cv_mixt_ener,ifl) = rl*eo
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initof







! ******************************************************************************
!
! Purpose: Initialize slip wall.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   pPatch      Pointer to patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nscbc_initsw(pRegion,pPatch)

  IMPLICIT NONE

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

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

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

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

  INTEGER :: errorflag,ifl,c1,bcopttype,distrib,gasmodel,indcp,indgs,indmol
  REAL(RFREAL) :: cp,eo,g,gc,mw,rl,rr,rul,rvl,rwl,rel,pl,pr,ul,ur,vl,vr,wl,wr
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,vals
  TYPE(t_global), POINTER :: global

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_NSCBC_InitSW',&
  'RFLU_ModNSCBC.F90')

! ******************************************************************************
! Main Body 
! ******************************************************************************

! ASSUMPTION : slip wall is horizontal
! => v=0

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol

  distrib =  ppatch%mixt%distrib

  ppatch%mixt%cvState = cv_mixt_state_cons

  IF (pregion%mixt%cvState /= cv_mixt_state_cons ) THEN
    CALL errorstop(global,err_cv_state_invalid,__line__, &
                   'region mixt cvState is incompatible...')
  END IF ! pRegion%mixt%cvState

  DO ifl = 1,ppatch%nBFaces
    c1 = ppatch%bf2c(ifl)

    mw = pregion%mixt%gv(gv_mixt_mol,indmol*c1)
    cp = pregion%mixt%gv(gv_mixt_cp ,indcp *c1)

    gc = mixtperf_r_m(mw)
    g  = mixtperf_g_cpr(cp,gc)

    rl  = pregion%mixt%cv(cv_mixt_dens,c1)
    rul = pregion%mixt%cv(cv_mixt_xmom,c1)
    rvl = pregion%mixt%cv(cv_mixt_ymom,c1)
    rwl = pregion%mixt%cv(cv_mixt_zmom,c1)
    rel = pregion%mixt%cv(cv_mixt_ener,c1)

    ul = rul/rl
    vl = rvl/rl
    wl = rwl/rl

    pl  = pregion%mixt%dv(dv_mixt_pres,c1)

    rr = rl
    ur = ul
    vr = 0.0_rfreal
    wr = wl
    pr = pl
! TEMPORARY : pressure at wall should be computed using rind states
!             but for now, copying from cell
 
    ppatch%mixt%cv(cv_mixt_dens,ifl) = rr
    ppatch%mixt%cv(cv_mixt_xmom,ifl) = rr*ur
    ppatch%mixt%cv(cv_mixt_ymom,ifl) = rr*vr
    ppatch%mixt%cv(cv_mixt_zmom,ifl) = rr*wr

    eo = mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)

    ppatch%mixt%cv(cv_mixt_ener,ifl) = rr*eo
  END DO ! ifl

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nscbc_initsw







END MODULE rflu_modnscbc

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModNSCBC.F90,v $
! Revision 1.7  2008/12/06 08:44:22  mtcampbe
! Updated license.
!
! Revision 1.6  2008/11/19 22:17:33  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.5  2007/02/27 13:05:23  haselbac
! Cosmetics only
!
! Revision 1.4  2006/10/20 21:18:45  mparmar
! Added mass/momentum coeffs in flux routines and cosmetic clean-up
!
! Revision 1.3  2006/08/21 16:11:19  haselbac
! Clean-up: comments, order of routines, missing headers, indentation, etc
!
! Revision 1.2  2006/08/19 19:44:08  haselbac
! Significant clean-up and cosmetic changes
!
! Revision 1.1  2006/08/19 15:37:46  mparmar
! Initial revision
!
! ******************************************************************************