Flux Function

    function Flux(ql, qr, du, inputs)
      !< Calculate the total flux through face for laminar flow.
      !--------------------------------------
      ! calculate the total flux through face
      !---------------------------------------
      implicit none
      real(wp), dimension(1:n_var), intent(in) :: ql !<left state
      real(wp), dimension(1:n_var), intent(in) :: qr !<right state
      !conservative form of updated neighbour
      real(wp), dimension(1:n_var), intent(in) :: du
      real(wp), dimension(1:7)    , intent(in) :: inputs
      real(wp), dimension(1:n_var)             :: Flux
      real(wp), dimension(1:n_var)             :: U !< conservative variables
      real(wp), dimension(1:n_var)             :: W !< new primitive variables
      real(wp), dimension(1:n_var)             :: P !< primitive variables of right cell

      !for extraction of the inputs
      real(wp) :: area
      real(wp) :: nx
      real(wp) :: ny
      real(wp) :: nz
      real(wp) :: volume
      real(wp) :: mmu
      real(wp) :: tmu


      real(wp)    :: dudx
      real(wp)    :: dudy
      real(wp)    :: dudz
      real(wp)    :: dvdx
      real(wp)    :: dvdy
      real(wp)    :: dvdz
      real(wp)    :: dwdx
      real(wp)    :: dwdy
      real(wp)    :: dwdz
      real(wp)    :: dTdx
      real(wp)    :: dTdy
      real(wp)    :: dTdz
      real(wp)    :: T1, T2
      real(wp)    :: uface
      real(wp)    :: vface
      real(wp)    :: wface
      real(wp)    :: trace
      real(wp)    :: Tauxx
      real(wp)    :: Tauyy
      real(wp)    :: Tauzz
      real(wp)    :: Tauxy
      real(wp)    :: Tauxz
      real(wp)    :: Tauyz
      real(wp)    :: Qx
      real(wp)    :: Qy
      real(wp)    :: Qz
      real(wp)    :: HalfRhoUsquare
      real(wp)    :: RhoHt
      real(wp)    :: K_heat
      real(wp)    :: FaceNormalVelocity
      real(wp)    :: mu

      area   = inputs(1)
      nx     = inputs(2)
      ny     = inputs(3)
      nz     = inputs(4)
      volume = inputs(5)
      mmu    = inputs(6)
      tmu    = inputs(7)


      !save the old stat in P
      P = qr

      ! find conservative variable
      U(1)   =   ql(1)
      U(2)   =   ql(1) * ql(2)
      U(3)   =   ql(1) * ql(3)
      U(4)   =   ql(1) * ql(4)
      U(5)   = ( ql(5) / (gm-1.0) ) + ( 0.5 * ql(1) * sum(ql(2:4)**2) )

      U(1:5) = U(1:5) + du(1:5)
      

      W(1)   =   U(1)
      W(2)   =   U(2) / U(1)
      W(3)   =   U(3) / U(1)
      W(4)   =   U(4) / U(1)
      W(5)   = (gm-1.0) * ( U(5) - ( 0.5 * SUM(U(2:4)**2) / U(1) ) )

      FaceNormalVelocity = (W(2) * nx) + (W(3) * ny) + (W(4) * nz)
      uface = 0.5 * ( W(2) + P(2) )
      vface = 0.5 * ( W(3) + P(3) )
      wface = 0.5 * ( W(4) + P(4) )


      Flux(1) =   W(1) * FaceNormalVelocity
      Flux(2) = ( W(2) * Flux(1) ) + ( W(5) * nx )
      Flux(3) = ( W(3) * Flux(1) ) + ( W(5) * ny )
      Flux(4) = ( W(4) * Flux(1) ) + ( W(5) * nz )

      HalfRhoUsquare = 0.5 * W(1) * ( W(2)*W(2) + W(3)*W(3) + W(4)*W(4) )
      RhoHt          = ( (gm/(gm-1.0)) * W(5) ) + HalfRhoUsquare
      Flux(5)        = RhoHt * FaceNormalVelocity


      ! viscous terms
      mu = mmu + tmu
      T1     =    W(5) / ( W(1) * R_gas )
      T2     =    P(5) / ( P(1) * R_gas )
      dTdx   =  ( T2   - T1   ) * nx * Area / Volume
      dTdy   =  ( T2   - T1   ) * ny * Area / Volume
      dTdz   =  ( T2   - T1   ) * nz * Area / Volume
      dudx   =  ( P(2) - W(2) ) * nx * Area / Volume
      dudy   =  ( P(2) - W(2) ) * ny * Area / Volume
      dudz   =  ( P(2) - W(2) ) * nz * Area / Volume
      dvdx   =  ( P(3) - W(3) ) * nx * Area / Volume
      dvdy   =  ( P(3) - W(3) ) * ny * Area / Volume
      dvdz   =  ( P(3) - W(3) ) * nz * Area / Volume
      dwdx   =  ( P(4) - W(4) ) * nx * Area / Volume
      dwdy   =  ( P(4) - W(4) ) * ny * Area / Volume
      dwdz   =  ( P(4) - W(4) ) * nz * Area / Volume

      trace = dudx + dvdy + dwdz
      Tauxx =  2. * mu * (dudx - trace/3.0)
      Tauyy =  2. * mu * (dvdy - trace/3.0)
      Tauzz =  2. * mu * (dwdz - trace/3.0)
      Tauxy = mu * (dvdx + dudy)
      Tauxz = mu * (dwdx + dudz)
      Tauyz = mu * (dwdy + dvdz)

      K_heat = ( mmu / Pr  + tmu/tpr) * gm * R_gas / ( gm - 1.0 )
      Qx = K_heat*dTdx
      Qy = K_heat*dTdy
      Qz = K_heat*dTdz

      Flux(2) = Flux(2) - ( Tauxx * nx + Tauxy * ny + Tauxz * nz )
      Flux(3) = Flux(3) - ( Tauxy * nx + Tauyy * ny + Tauyz * nz )
      Flux(4) = Flux(4) - ( Tauxz * nx + Tauyz * ny + Tauzz * nz )
      Flux(5) = Flux(5) - ( Tauxx * uface + Tauxy * vface + Tauxz * wface + Qx ) * nx
      Flux(5) = Flux(5) - ( Tauxy * uface + Tauyy * vface + Tauyz * wface + Qy ) * ny
      Flux(5) = Flux(5) - ( Tauxz * uface + Tauyz * vface + Tauzz * wface + Qz ) * nz

      Flux    = Flux * Area

    end function Flux