C FLUFU1    SOURCE    BECC      09/12/07    21:15:21     6579
      subroutine flufu1(nordpo,
     &     Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
     &     Rgas_l, acv_l, acv_z,
     &     z,
     &     d_l, d_d, d_c, d_b, d_a, d_r,
     &     rho_l, T_l, u_l, hf_l,
     &     rho_d, T_d, u_d,
     &     rho_c, T_c, u_c,
     &     rho_b, T_b, u_b, hf_b,
     &     rho_a, T_a, u_a,
     &     rho_r, T_r, u_r, hf_r,
     &     flurho, fluru, fluret)
C
c************************************************************************
c
c project           :
c
c name              :  flufu1
c
c description       :  Reactive Euler equations for a mixture of reactive
C                      thermally perfect gases.
c
c                      Given z = x/t, we compute the flux
C                      F(x/t) - ((x/t) * U)
C                      We compute density, temperature and velocity by
C                      averaging the intermediate states using the
C                      weve speeds.
C                      We compute the other variables using these
C                      quantities and the EOS
c
c language          :  fortran 77
c
c author            :  a. beccantini den/dm2s/sfme/ltmf
c
c************************************************************************
c
C called by         :  fbecre
C
c
c************************************************************************
c
c**** input
c
c     nordpo          =  polynomial order of the cv
c
C     Tmaxcv          =  Tmax for cv expansion
C
C     z               =  x / t
C
C     Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
C     Rgas_l, acv_l   = quantities to compute cv and cp
C
C     acv_z           = work vector
C
C     d_l, d_d, un_d, un_c, d_c, d_b, d_a, d_r =
C                     = wave speeds
C
C     rho_, u_, T_, hf_
C                     = density, velocity, temperature,
C                       formation enthalpy
C**** output
C
c     flurho, flurhou, fluret = fluxes for rho, rho u, rho e_t
C
c
c************************************************************************
c
C     07/12/2009      created
C
c************************************************************************
c
c n.b.: all variables are declared
c
c      implicit none
      integer nordpo, iord
      real*8 Tmaxcv, Rgas_r, acv_r(0:nordpo)
     &     , Rgas_b, acv_b(0:nordpo)
     &     , Rgas_l, acv_l(0:nordpo)
     &     , z, acv_z(0:nordpo)
     &     , d_l, d_d, d_c, d_b, d_a, d_r
     &     , rho_l, u_l, T_l, hf_l
     &     , rho_d, T_d, u_d
     &     , rho_c, T_c, u_c
     &     , rho_b, T_b, u_b, hf_b
     &     , rho_a, T_a, u_a
     &     , rho_r, T_r, u_r, hf_r
     &     , flurho, fluru, fluret
      real*8 u_cd, fp_l, fp_d, fp_c, fp_b, fp_a, fp_r
     &     , Rgas_z, rho_z, T_z, u_z, p_z, e_z, cv_z, hf_z
     &     , ecin_z
      logical logdeb
      parameter(logdeb = .false.)
CC
CC     Test
CC
C      z = 1.13D0
C      d_l = -1.0D0
C      d_d = 0.10D0
C      u_c = 0.2D0
C      u_d = 0.2D0
C      d_c = 1.09D0
C      d_b = 1.10D0
C      d_a = 1.11D0
C      d_r = 1.12D0
C
C     Speed of the contact discontinuty u_cd = u_c = u_d
C
      u_cd = 0.5D0 * (u_c + u_d)
C
C**** Weight functions
C
      fp_l = 0.0d0
      fp_d = 0.0D0
      fp_c = 0.0d0
      fp_b = 0.0d0
      fp_a = 0.0d0
      fp_r = 0.0d0
C
C     This can be simplified from a computational point of view
C
      if (z .lt. d_l) then
         fp_l = 1.0d0
      elseif(z .lt. d_d) then
C        Rarefaction if d_l < d_d
         fp_l = (d_d - z) / (d_d - d_l)
         fp_d = 1.0d0 - fp_l
      elseif(z .lt. u_cd)then
C        Contact discontinuty, left
         fp_d = 1.0d0
      elseif(z .lt. d_c)then
C        Contact discontinuty, right
         fp_c = 1.0d0
      elseif(z .lt. d_b)then
C        Rarefaction if d_b < d_c
         fp_c = (d_b - z) / (d_b - d_c)
         fp_b = 1.0d0 - fp_c
      elseif(z .lt. d_a)then
C        RS, right
         fp_a = 1.0D0
      elseif(z .lt. d_r)then
C        Rgnl rarefaction if d_r < d_a
         fp_a = (d_r - z) / (d_r - d_a)
         fp_r = 1.0d0 - fp_a
      else
C        Rgnl, right
         fp_r = 1.0d0
      endif
CC
C
C**** Shock-shock + (e.f) intermediate solution
C
      Rgas_z = ((fp_l + fp_d) * Rgas_l) +
     &     ((fp_c + fp_b) * Rgas_b)  +
     &     ((fp_a + fp_r) * Rgas_r)
      rho_z = (fp_l * rho_l) + (fp_d * rho_d) +
     &     (fp_c * rho_c) + (fp_b * rho_b) +
     &     (fp_a * rho_a) + (fp_r * rho_r)
      T_z = (fp_l * T_l) + (fp_d * T_d) +
     &     (fp_c * T_c) + (fp_b * T_b) +
     &     (fp_a * T_a) + (fp_r * T_r)
      p_z = Rgas_z * rho_z * T_z
      u_z = (fp_l * u_l) + (fp_d * u_d) +
     &     (fp_c * u_c) + (fp_b * u_b) +
     &     (fp_a * u_a) + (fp_r * u_r)
      if (logdeb) then
         write(*,*) 'd_l, d_d, d_c, d_b, d_a, d_r'
         write(*,'(6E12.4)') d_l, d_d, d_c, d_b, d_a, d_r
         write(*,*) 'fp_l, fp_d, fp_c, fp_b, fp_a, fp_r'
         write(*,'(6E12.4)') fp_l, fp_d, fp_c, fp_b, fp_a, fp_r
         write(*,*) 'State in x/t'
         write(*,*) 'Rgas_z, rho_z, T_z, p_z, u_z'
         write(*,*) Rgas_z, rho_z, T_z, p_z, u_z
      endif
C
C**** Computation of thermal energy
C
      do iord = 0, nordpo , 1
         acv_z(iord) = ((fp_l + fp_d) * acv_l(iord)) +
     &        ((fp_c + fp_b) * acv_b(iord)) +
     &        ((fp_a + fp_r) * acv_r(iord))
      enddo
C
      call prith1(nordpo, acv_z, Tmaxcv, T_z, e_z, cv_z)
C
      hf_z = ((fp_l + fp_d) * hf_l) + ((fp_c + fp_b) * hf_b)
     &     + ((fp_a + fp_r) * hf_r)
      ecin_z = 0.5d0 * u_z * u_z
C
C**** Interfacial flux
C     According to NKONGA, Comput Methods Appl. Mech Engnr 190, 2000
C     \dep{U}{x} + \dep{F(U)}{x} = 0
C     z = speed of the surface
C     f_z = F - z U
C
      flurho = rho_z * (u_z - z)
C      write(*,*) 'flurho', flurho
      fluru = (flurho * u_z) + p_z
      fluret = (flurho * (e_z + hf_z +
     &     ecin_z)) + (p_z * u_z)
C      write(*,*) (et_z + ef_z +
C     &     ecin_z), flurho
C      write(*,*) p_z, u_z
      return
      end