C FLUFUN    SOURCE    BECC      09/12/07    21:15:24     6579           
      subroutine flufun(nordpo,
     &     Tmaxcv, acv_z,
     &     z,
     &     vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r,
     &     Rgas_l, acv_l, rho_l, u_l, T_l, hf_l, 
     &     rho_ls, T_ls,
     &     Rgas_t, acv_t, rho_t, u_t, T_t, hf_t,
     &     rho_rs, u_rs, T_rs,
     &     Rgas_r, acv_r, rho_r, u_r, T_r, hf_r,
     &     flurho, fluru, fluret)
C
c************************************************************************
c
c project           :  
c
c name              :  flufun
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 state
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         :  fcolte
C
c
c************************************************************************
c
c**** input
c
c     nordp0          =  polynomial order of the cv
c
C     Tmaxcv          =  Tmax for cv expansion
C
C     z               =  x / t
C
C     vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r =
C                     =  wave speeds
C
C     Rgas_l, acv_l, rho_l, u_l, T_l, hf_l =
C                     = some quantities in l
C    
C     rho_ls, T_ls    = some quantities in ls
C
C     Rgas_t, acv_t, rho_t, u_t, T_t, hf_t =
C                     = some quantities in rss
C
C     rho_rs, u_rs, T_rs =
C                     = some quantities in rss
C
C     Rgas_r, acv_r, rho_r, u_r, T_r, hf_r =
C                     = some quantities in r
C
C     acv_z = work vector
C                     
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, acv_z(0:nordpo), z
     &     , vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r
     &     , Rgas_l, acv_l(0:nordpo), rho_l, u_l, T_l, hf_l
     &     , rho_ls, T_ls
     &     , rho_t, u_t, T_t, hf_t
     &     , Rgas_t, acv_t(0:nordpo), rho_rs, u_rs, T_rs
     &     , Rgas_r, acv_r(0:nordpo), rho_r, u_r, T_r, hf_r
     &     , flurho, fluru, fluret
      real*8 fp_l, fp_ls, fp_t, fp_rs, fp_r
      real*8 Rgas_z, rho_z, T_z, p_z, u_z, et_z
     &     , ef_z, ecin_z
      real*8 T1, funt, cvtot, dcv
      logical logdeb
      parameter (logdeb = .false.)
C
C**** Weight functions
C
      fp_l = 0.0d0
      fp_ls = 0.0d0 
      fp_t = 0.0d0
      fp_rs = 0.0d0
      fp_r = 0.0d0
C
C     This can be simplified from a computational point of view
C
      if (z .lt. vit_l) then
         fp_l = 1.0d0 
      else
         if(z .lt. vit_ls) then
            fp_l = (vit_ls - z) / (vit_ls - vit_l)
            fp_ls = 1.0d0 - fp_l
         else
            if(z .lt. u_ls)then
               fp_ls = 1.0d0
            else
               if(z .lt. D_resh)then
                  fp_t = 1.0d0
               else
                  if(z .lt. vit_rs)then
                     fp_rs = 1.0d0
                  else
                     if(z .lt. vit_r)then
                        fp_r = (z - vit_rs) / (vit_r - vit_rs)
                        fp_rs = 1.0d0 - fp_r
                     else
                        fp_r = 1.0d0
                     endif
                  endif
               endif
            endif
         endif
      endif
C
C**** Shock-shock + (e.f) intermediate solution
C
      Rgas_z = ((fp_l + fp_ls) * Rgas_l) + 
     &     (fp_t * Rgas_t) + 
     &     ((fp_rs + fp_r) * Rgas_r)
      rho_z = (fp_l * rho_l) + (fp_ls * rho_ls) + 
     &     (fp_t * rho_t) + (fp_rs * rho_rs) + 
     &     (fp_r * rho_r)
      T_z = (fp_l * T_l) + (fp_ls * T_ls) + 
     &     (fp_t * T_t) + (fp_rs * T_rs) + 
     &     (fp_r * T_r)
      p_z = Rgas_z * rho_z * T_z 
      u_z = (fp_l * u_l) + (fp_ls * u_ls) + 
     &     (fp_t * u_t) + (fp_rs * u_rs) + 
     &     (fp_r * u_r)
      if (logdeb) then 
         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_ls) * acv_l(iord)) +
     &        (fp_t * acv_t(iord)) +
     &        ((fp_rs + fp_r) * acv_r(iord))
      enddo
C     T1 = min (T, Tmaxcv)
      T1 = 0.5d0 * ((T_z + Tmaxcv) - abs (T_z - Tmaxcv))
C         
      funt = 1.0D0
      cvtot = acv_z(0)
      et_z = cvtot * T1
      do iord = 1, nordpo, 1
         funt = funt * T1
         dcv = acv_z(iord) * funt
         cvtot = cvtot + dcv
         et_z = et_z + (dcv * T1 / (iord + 1))
      enddo
C
C     if ( T > Tmaxcv ) => T1 = (T - Tmaxcv)
C     else  T1 = 0
C     i.e.
C     T1 = max ((T - Tmaxcv),0)
C
      T1 = 0.5d0 * ((T_z - Tmaxcv) + abs (T_z - Tmaxcv))
      et_z = et_z + (cvtot * T1)
C
      ef_z = ((fp_l + fp_ls) * hf_l) + (fp_t * hf_t) 
     &     + ((fp_r + fp_rs) * hf_r)
      ecin_z = 0.5d0 * u_z * u_z
C      write(*,*) 'et, ecin, ef', et_z, ecin_z, ef_z
C      write(*,*) 'rho_z, u_z, p_z', rho_z, u_z, p_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 * (et_z + ef_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