Sources Esope | Cast3M

fbecre
C FBECRE    SOURCE    BECC      09/12/07    21:15:16     6579
      subroutine fbecre(ndim, nesp, nLHS, nordpo,
     &     reacoe, aprop, Runiv,
     &     Tmaxcv,
     &     cnna, ctna, ct1na,
     &     wvec_l,
     &     wvec_r,
     &     wwork, y, acv_l, acv_r, acv_b, acv_w   ,
     &     flu, f_lag, cellt,
     &     logan)
C
c************************************************************************
c
c project           :
c
c name              :  fbecre
C                      reaction shock in WDF, CJDF, CJDT, SDT
c
c description       :  Reactive Euler equations for a mixture of reactive
C                      thermally perfect gases.
c
c                      Shock-shock fds + entropy fix, implemented in a
C                      "Discrete Equation Method" fashion.
c
c                      Voir:
c                      1) Beccantini, report ???
C                      2) Metayer, Massoni, Saurel, "Modelling
C                         evaporations fronts with reactive Riemann
C                          solvers.", JCP 205, 2005
c
c language          :  fortran 77
c
c author            :  a. beccantini den/dm2s/sfme/ltmf
c
c************************************************************************
c
c called sub        :  prithe, racre2, flufun
C
C called by         :  comres
c
c************************************************************************
c
c**** input
C
C     ndim                =  dimension of the domain
c
c     nesp                =  species involved in the Euler equations
C
C     nLHS                =  species involved in the LHS of the global
C                            reaction
c
c     nordp0              =  polynomial order of the cv
c
c     reacoe              =  coefficient involved in the global reaction
c                            (positive in the LHS, negative in the RHS...
c                            the first one should be positive)
c
c     aprop               =  gas properties (cv, molar mass, formation
C                            energy)
c
C     Runiv               =  universal gas constant
c
C     Tmaxcv              =  Tmax for cv expansion
C
C     cnna, ctna, ct1na   = scalar products
C                           (\vec{n},\vec{n}_{\rm alpha})
C                           (\vec{t},\vec{n}_{\rm alpha})
C                           (\vec{t1},\vec{n}_{\rm alpha})
C                           HP. We suppose that \vec{n}_{\rm alpha} is
C                               always from the burnt to the unburnt state
C
c     wvec_l             =  left primitive variables
c
c     wvec_r             =  right primitive variables
C
C                            NB Primitive variables =
C
C                            w(1)             = rho
C                            w(2)             = ux
C                            w(3)             = uy
C                            w(4)             = uz
C                            w(2+ndim)        = p
C                            w(3+ndim)        = csi
C                            w(3+ndim+1)      = yf(1)
C                                ...
C                            w(3+ndim+nLHS)   = yf(nLHS)
C                            w(3+ndim+nLHS+1) = yi(nLHS+1)
C                                ...
C                            w(3+ndim+nesp-1) = yi(nesp-1)
C                            w(3+ndim+nesp)   = dy1t = yi(1) - yf(1)
C                            w(4+ndim+nesp)   = k0
C
C                            For ndim = 1, the total number of variables
C                            is 4 + 1 + nesp = 5 + nesp.
c
c     wwork, y, acv_l,
C     acv_r, acv_b,
C     acv_w               =  temporary work vectors
C                            wwork(4+ndim+nesp), y(1:nesp),
C                            acv_l(0:nordpo),
C                            acv_r(0:nordpo),
C                            acv_b(0:nordpo),
C                            acv_w   (0:nordpo)
C
c**** output
c
c     flu                 =  Eulerian interfacial flux in (n,t1,t2), i.e.
c                            rho*un                 mass
C                            rho*un*un + p          momentum along n
C                            rho*un*ut1             momentum along t1
*                            rho*un*ut2             momentum along t2
C                            rho*un*ht              energy
c                            0.0                    alpha
C                            rho*un*y_{1,f}         mass * y_{1,f}
C                            ...
C                            rho*un*y_{nLHS+1,i}    mass * y_{nLHS+1,i}
c                            ...
C                            rho*un*dy_{1,t}        mass * dy_{1,t}
C                            rho*un*k0              mass * k0
C
C                            NB. flu(1:(4+ndim+nesp))
c
C     f_lag             =  Lagrangian interfacial flux on the reactive
C                            discontinuty in (n,t1,t2), i.e.
c                            rho*(un - D_resh)                 mass
C                            rho*un*un + p - (D_resh rho un)   momentum along n
C                            rho*(un - D_resh)*ut1             momentum along t1
C                            rho*(un - D_resh)*ut2             momentum along t2
C                            rho*un*ht - (D_resh rho et)       energy
c                            - D_resh                          alpha
c                            rho*(un - D_resh)*y_{1,f}         mass * y_{1,f}
C                            ...
C
C                            NB. f_lag(1:(4+ndim+nesp))
C
c
c     cellt               =  stability condition, i.e.
c
c                            dt/dx &lt; cellt (dimension 1/velocity)
c
c     logan               =  si true, anomaly detected
c
c************************************************************************
c
C     07/12/2009      created
C
c************************************************************************
c
c n.b.: all variables are declared
c
c      implicit none
      integer ndim, nesp, nLHS, nordpo, iter
      integer icomp
      real*8 zero
      parameter (zero=1.0d-10)
      real*8 rind_l
      real*8  reacoe(nesp), aprop(0:(nordpo+2),nesp), Runiv, Tmaxcv
     &     , wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp)
     &     , flu(4+ndim+nesp)
     &     , f_lag(4+ndim+nesp),  cellt
      real*8 wwork(4+ndim+nesp)
     &     , acv_l(0:nordpo), acv_r(0:nordpo), acv_b(0:nordpo)
     &     , acv_w(0:nordpo)
     $     , y(1:nesp)
      real*8 csi0, csi1, k0
      real*8 scnna, cnna, ctna, ct1na
      real*8 gam_b, Rgas_b, cv_b, hf_b, rho_b, u_b(1:3), un_b, p_b, T_b
     &     , et_b, csi_b
     &     , gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l
     &     , et_l, csi_l
     &     , gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_r
     &     , et_r, csi_r
     &     , un_r, un_l, q
      real*8 u_l1, u_r1
      real*8 d_l, d_r
      real*8 rho_d, p_d, un_d, T_d, d_d
     &     , rho_c, p_c, un_c, T_c, d_c
     &     , d_b
     &     , rho_a, p_a, un_a, T_a, d_a
      real*8 xst, flurho, fluru, fluret, coefy, omcoey
     &     , trekin, D_resh, k0new
      real*8 coefd, d_o
     &     , k0m
     &     , d_lm, d_rm
     &     , rho_dm, p_dm, un_dm, T_dm, d_dm
     &     , rho_cm, p_cm, un_cm, T_cm, d_cm
     &     , rho_bm, p_bm, un_bm, T_bm, d_bm
     &     , rho_am, p_am, un_am, T_am, d_am
     &     , k0n
     &     , d_ln, d_rn
     &     , rho_dn, p_dn, un_dn, T_dn, d_dn
     &     , rho_cn, p_cn, un_cn, T_cn, d_cn
     &     , rho_bn, p_bn, un_bn, T_bn, d_bn
     &     , rho_an, p_an, un_an, T_an, d_an
     &     , k0o
     &     , d_lo, d_ro
     &     , rho_do, p_do, un_do, T_do, d_do
     &     , rho_co, p_co, un_co, T_co, d_co
     &     , rho_bo, p_bo, un_bo, T_bo, d_bo
     &     , rho_ao, p_ao, un_ao, T_ao, d_ao
C
      logical logan
      logical lognc
      logical logdeb, logk0n
C     logdeb = debugging ?
C     logk0n = .true.  : instead of k0, we take k0 * abs(cnna)
C              .false. : instead of k0, we take k0 * abs(cnna)
      parameter (logdeb = .false., logk0n = .true.)
C
C     HP $\vec{n}_{\rm \alpha}$ is always from the burnt to the unburnt
C     region
C
      scnna = sign(1.0D0, cnna)
C
C     In the reactive case, we have to distingish between 6 different
C     states
C
C     l    = burnt state (csi = csimax)
C                        | Left GNL wave
C     ls   = connected to l by a GNL wave (csi = csimax)
C          = d
C                        | Contact discontinuty
C     rsss = connected to rss by a rarefaction wave (true
C            in CJDF and CJDT regime, = rss in WDF and SDT)
C            (csi = csimax)
C          = c
C                        | Rarefaction wave
C     rss  = connected to rs by a reactive wave (csi = csimax)
C          = b
C                        | Reactive shock (WDF, CJDF, CJDT, SDT)
C     rs   = connected to r by a GNL wave (csi = csimin)
C            in SDT, state behind the non-reactive shock travelling
C            at the same speed than the reactive shock
C          = a
C                        | Right GNL wave
C     r    = unburnt state (csi = csimin)
C
C
C     In the non-reactive case, we have to distingish between 4 different
C     states:
C
C     l    = burnt/unburnt state
C                        | Left GNL wave
C     ls   = connected to l by a GNL wave
C          = b
C                        | Contact discontinuty
C     rs   = connected to r by a GNL wave
C          = a
C                        | Right GNL wave
C     r    = burnt/unburnt state
C
      csi0 = wvec_l(ndim+3)
      csi1 = wvec_r(ndim+3)
C
C     csi0 = 0 or 1
C     csi1 = 0 or 1
C
      if (abs(csi0 - csi1) .lt. 1.0D-1)then
C      if (.false.) then
C
C************************************************************************
C******* NON REACTIVE CASE **********************************************
C************************************************************************
C
         call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $        Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
     $        u_l, p_l, T_l,y, et_l, csi_l, logan)
         if (logan) then
            write(*,'(/A8)') 'fbecre.f'
            write(*,*) 'ANOMALY DETECTED 01'
            goto 9999
         endif
C
         call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $        Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
     $        u_r, p_r, T_r,y, et_r, csi_r, logan)
         if (logan) then
            write(*,'(/A8)') 'fcolre.f'
            write(*,*) 'ANOMALY DETECTED 02'
            goto 9999
         endif
C
         un_l = u_l(1)
         un_r = u_r(1)
C
         call racnre(nordpo, Tmaxcv,
     &        Rgas_l, acv_l, Rgas_r, acv_r,
     &        rho_l, p_l, un_l, T_l, gam_l,
     &        rho_r, p_r, un_r, T_r, gam_r,
     &        d_l, d_r,
     &        rho_b, p_b, un_b, T_b, d_b,
     &        rho_a, p_a, un_a, T_a, d_a,
     &        logan, lognc)
C
         cellt = max (abs(d_l), abs(d_r))
C
C******* States
C
C         l
C        ***
C          |*                              r
C          | *   b                 *********
C          |  *******             *|
C          |  |     *   a        * |
C          |  |     *************  |
C          |  |     |           |  |
C          |  |     |           |  |
C        d_l d_b  un_b=un_a   d_a d_r
C
C
         xst = 0.0D0
         call flufu2(nordpo,
     &        Tmaxcv, Rgas_r, acv_r,
     &        Rgas_l, acv_l, acv_w,
     &        xst,
     &        d_l, d_b, d_a, d_r,
     &        rho_l, T_l, un_l, hf_l,
     &        rho_b, T_b, un_b,
     &        rho_a, T_a, un_a,
     &        rho_r, T_r, un_r, hf_r,
     &        flurho, fluru, fluret)
C
         flu(1) = flurho
         flu(2) = fluru
         flu(2+ndim) = fluret
         flu(3+ndim) = 0.0d0
C        The Larrouturou theorem does hold
C        Namely un_a (= un_b) has the same sign as flurho
C        write (*,*) un_d, un_c, flurho
C        un_ls > 0 => coefy = 1
         coefy = 0.5d0 * (1.0d0 + sign(1.0d0, flurho))
         omcoey = 1.0d0 - coefy
C        Tangential speeds (ndim .ge. 2)
         do icomp = 3, (ndim + 1), 1
            flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
     &           (omcoey * wvec_r(icomp)))
         enddo
         do icomp = (ndim + 4), (ndim + nesp + 4), 1
            flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
     &           (omcoey * wvec_r(icomp)))
         enddo
C        Kinetic energy flux
         trekin = 0.0D0
         do icomp = 2, ndim, 1
            trekin = trekin + (((coefy * wvec_l(icomp + 1)) +
     &           (omcoey * wvec_r(icomp + 1))) ** 2)
         enddo
         trekin = 0.5D0 * trekin
         flu(2+ndim) = flu(2+ndim) + (flurho * trekin)
C
C******* f_lag is here not defined !!!
C
      else
C
C************************************************************************
C******* REACTIVE CASE **************************************************
C************************************************************************
C
C        We define the left state such that csi_l > csi_r
C
         if (csi0 .gt. csi1) then
            k0 = wvec_r(4+ndim+nesp)
            rind_l = 1.0D0
C
C           The rss (b) state is the burnt state relatively to the
C           right and right star (a) state. It has the same value of csi as
C           the left state.
C           Note that this state has not necessary the same mass
C           fractions as the left state.
C           We compute the physical properties of the rss (b) state
C           namely
C           gam_b, Rgas_b, cv_b, acv_b, hf_b.
C           The other physical values are of course meaningless!
C
            do icomp = 1, (nesp + 4 + ndim), 1
               wwork(icomp) = wvec_r(icomp)
            enddo
            wwork(3+ndim) = wvec_l(3+ndim)
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wwork,gam_b, Rgas_b, cv_b, acv_b, hf_b,
     $           rho_b, u_b, p_b, T_b ,y, et_b, csi_b, logan)
            un_b = u_b(1)
            if (logan) then
               write(*,'(/A8)') 'fbecre.f'
               write(*,*) 'ANOMALY DETECTED 1'
               goto 9999
            endif
C
C           The left state is considered as the burnt state.
C           We compute the physical properties of the left state, namely
C           gam_l, Rgas_l, cv_l, acv_l, hf_l.
C           We also compute the physical variable on the left, namely
C           rho_l, u_l, p_l, T_L, et_l and also csi_r
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
     $           u_l, p_l, T_l,y, et_l, csi_l, logan)
            if (logan) then
               write(*,'(/A8)') 'fbecre.f'
               write(*,*) 'ANOMALY DETECTED 2'
               goto 9999
            endif
C
C           We compute the normal speed in the frame
C           (\vec{n}_{\rm alpha}, \vec{t}_{\rm alpha},
C            \vec{t1}_{\rm alpha})
C
            u_l1 = u_l(1) * cnna * scnna
            if (ndim .ge. 2) then
               u_l1 = u_l1 + (u_l(2) * ctna * scnna)
            endif
            if (ndim .eq. 3) then
               u_l1 = u_l1 + (u_l(3) * ct1na * scnna)
            endif
C
C           The right state is considered as the unburnt state.
C           We compute the physical properties of the right state, namely
C           gam_r, Rgas_r, cv_r, acv_r, hf_r.
C           We also compute the physical variable on the left, namely
C           rho_r, u_r, p_r, T_r, et_r, csi_r
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
     $           u_r, p_r, T_r,y, et_r, csi_r, logan)
            if (logan) then
               write(*,'(/A8)') 'fcolre.f'
               write(*,*) 'ANOMALY DETECTED 3'
               goto 9999
            endif
            u_r1 = u_r(1) * cnna * scnna
            if (ndim .ge. 2) then
               u_r1 = u_r1 + (u_r(2) * ctna * scnna)
            endif
            if (ndim .eq. 3) then
               u_r1 = u_r1 + (u_r(3) * ct1na * scnna)
            endif
         else
            rind_l = -1.0D0
            k0 = wvec_l(4+ndim+nesp)
C
            do icomp = 1, (nesp + 4 + ndim), 1
               wwork(icomp)=wvec_l(icomp)
            enddo
            wwork(3+ndim)=wvec_r(3+ndim)
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wwork,gam_b, Rgas_b, cv_b, acv_b, hf_b,
     $           rho_b, u_b, p_b, T_b, y, et_b, csi_b, logan)
            un_b = u_b(1)
            if (logan) then
               write(*,'(/A8)') 'fcolre.f'
               write(*,*) 'ANOMALY DETECTED 4'
               goto 9999
            endif
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wvec_r,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
     $           u_l, p_l, T_l,y, et_l, csi_l, logan)
            if (logan) then
               write(*,'(/A8)') 'fcolre.f'
               write(*,*) 'ANOMALY DETECTED 5'
               goto 9999
            endif
            u_l1 = u_l(1) * cnna * scnna
            if (ndim .ge. 2) then
               u_l1 = u_l1 + (u_l(2) * ctna * scnna)
            endif
            if (ndim .eq. 3) then
               u_l1 = u_l1 + (u_l(3) * ct1na * scnna)
            endif
C
            call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
     $           Tmaxcv,wvec_l,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
     $           u_r, p_r, T_r,y, et_r, csi_r, logan)
            if (logan) then
               write(*,'(/A8)') 'fcolre.f'
               write(*,*) 'ANOMALY DETECTED 6'
               goto 9999
            endif
            u_r1 = u_r(1) * cnna * scnna
            if (ndim .ge. 2) then
               u_r1 = u_r1 + (u_r(2) * ctna * scnna)
            endif
            if (ndim .eq. 3) then
               u_r1 = u_r1 + (u_r(3) * ct1na * scnna)
            endif
         endif
         un_r = rind_l * u_r(1)
         un_l = rind_l * u_l(1)
         u_r1 = rind_l * u_r1
         u_l1 = rind_l * u_l1
C
C**** Test
C
         if (logdeb) then
            write(*,'(A23)') 'rho_r,  p_r,  t_r,  u_r, u_r1'
            write(*,'(4E16.6)') rho_r, p_r, t_r, un_r, u_r1
            write(*,'(A23)') 'rho_l,  p_l,  t_l,  u_l, u_l1'
            write(*,'(4E16.6)') rho_l, p_l, t_l, un_l, u_l1
C           stop
         endif
C
C
C        The reaction should be exothermic
C
         q = hf_r - hf_b
         if (q .lt. (-1.0D0 * zero * max( abs(hf_b), abs(hf_r) )) ) then
            write(*,'(/A8)') 'fcolre.f'
            write(*,*) 'Endothermic reaction ???                '
            write(*,*) 'q =', q
            logan = .true.
            goto 9999
         endif
         q = 0.5d0 * (q + abs(q))
C
C        Computation of the speed of the wave
C        It can deal with vacuum
C
         call racre2(nordpo, Tmaxcv,
     &        Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &        q, k0,
     &        rho_l, p_l, u_l1, T_l, gam_l,
     &        rho_r, p_r, u_r1, T_r, gam_r,
     &        d_l, d_r,
     &        rho_d, p_d, un_d, T_d, d_d,
     &        rho_c, p_c, un_c, T_c, d_c,
     &        rho_b, p_b, un_b, T_b, d_b,
     &        rho_a, p_a, un_a, T_a, d_a,
     &        logan, lognc)
C
         if (ndim .eq. 1)then
C
C***************************
C********** 1D *************
C***************************
C
            cellt = max (abs(d_l), abs(d_r))
            D_resh = d_b
         elseif (logk0n) then
C
C******************************
C********** K0 n \cdot n_a ****
C******************************
C
C
C           Intermediate state
C           It can deal with vacuum
C
            k0new = max ((abs ((d_b - un_a) * cnna)), (k0*1.0D-3))
C            k0new = k0
            call racre2(nordpo, Tmaxcv,
     &           Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &           q, k0new,
     &           rho_l, p_l, un_l, T_l, gam_l,
     &           rho_r, p_r, un_r, T_r, gam_r,
     &           d_l, d_r,
     &           rho_d, p_d, un_d, T_d, d_d,
     &           rho_c, p_c, un_c, T_c, d_c,
     &           rho_b, p_b, un_b, T_b, d_b,
     &           rho_a, p_a, un_a, T_a, d_a,
     &           logan, lognc)
C
            D_resh = d_b
C
            if (logan) then
               write(*,'(/A8)') 'fbecre.f'
               write(*,*) 'ANOMALY DETECTED 7'
               goto 9999
            endif
C
            cellt = max (abs(d_l), abs(d_r))
C
         else
C
C******************************
C********** D n \cdot n_a *****
C******************************
C
C
C            coefd = rind_l * cnna
C            write(*,*) 'd_b, k0 = ', d_b, k0
            coefd = max( abs (cnna), 1.0D-3)
            d_o = d_b * coefd
C
C********** Linear interpolation
C
            k0m = 1.0D-3 * k0
            call racre2(nordpo, Tmaxcv,
     &           Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &           q, k0m,
     &           rho_l, p_l, un_l, T_l, gam_l,
     &           rho_r, p_r, un_r, T_r, gam_r,
     &           d_lm, d_rm,
     &           rho_dm, p_dm, un_dm, T_dm, d_dm,
     &           rho_cm, p_cm, un_cm, T_cm, d_cm,
     &           rho_bm, p_bm, un_bm, T_bm, d_bm,
     &           rho_am, p_am, un_am, T_am, d_am,
     &           logan, lognc)
C
            if (logan) then
               write(*,'(/A8)') 'fbecre.f'
               write(*,*) 'ANOMALY DETECTED 8'
               goto 9999
            endif
            k0m = abs (d_bm - un_am)
C            write(*,*) 'd_bm, k0m = ', d_bm, k0m
C
            k0n = 1.0D3 * k0
            call racre2(nordpo, Tmaxcv,
     &           Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &           q, k0n,
     &           rho_l, p_l, un_l, T_l, gam_l,
     &           rho_r, p_r, un_r, T_r, gam_r,
     &           d_ln, d_rn,
     &           rho_dn, p_dn, un_dn, T_dn, d_dn,
     &           rho_cn, p_cn, un_cn, T_cn, d_cn,
     &           rho_bn, p_bn, un_bn, T_bn, d_bn,
     &           rho_an, p_an, un_an, T_an, d_an,
     &           logan, lognc)
C
            if (logan) then
               write(*,'(/A8)') 'fbecre.f'
               write(*,*) 'ANOMALY DETECTED 9'
               goto 9999
            endif
            k0n = abs (d_bn - un_an)
C            write(*,*) 'd_bn, k0n = ', d_bn, k0n
C
            if (d_o .le. d_bm) then
               d_l = d_lm
               d_r = d_rm
C
               rho_d = rho_dm
               p_d = p_dm
               un_d = un_dm
               T_d = T_dm
               d_d = d_dm
C
               rho_c = rho_cm
               p_c  = p_cm
               un_c = un_cm
               T_c = T_cm
               d_c = d_cm
C
               rho_b = rho_bm
               p_b = p_bm
               un_b = un_bm
               T_b = T_bm
               d_b = d_bm
C
               rho_a = rho_am
               p_a = p_am
               un_a = un_am
               T_a = T_am
               d_a = d_am
C
C               write(*,*) 'Sono piccolino come il Masakino !'
C               write(*,*) 'coefd , k0 =', coefd, k0m
C               write(*,*) 'do =', d_o
C               write(*,*) 'dn =', d_bn
C               write(*,*) 'dm =', d_bm
C               write(*,*) 'cnna =', cnna
C
C               k0o = 1.0D-3 * k0
C               do iter = 1, 10
C                  call racre2(nordpo, Tmaxcv,
C     &                 Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
C     &                 q, k0o,
C     &                 rho_l, p_l, un_l, T_l, gam_l,
C     &                 rho_r, p_r, un_r, T_r, gam_r,
C     &                 d_lo, d_ro,
C     &                 rho_do, p_do, un_do, T_do, d_do,
C     &                 rho_co, p_co, un_co, T_co, d_co,
C     &                 rho_bo, p_bo, un_bo, T_bo, d_bo,
C     &                 rho_ao, p_ao, un_ao, T_ao, d_ao,
C     &                 logan, lognc)
C
C                  if (logan) then
C                     write(*,'(/A8)') 'fbecre.f'
C                     write(*,*) 'ANOMALY DETECTED 10'
C                     goto 9999
C                  endif
C
C                  k0o = d_bo - un_ao
C                  write(*,*) k0o
C
C                  write(*,*) 'd_b0, k0o = ', d_bo, k0o
C                  k0o = k0o + 0.1D0 * (k0n - k0m)
C               enddo
C
C               stop
C
            elseif (d_o .ge. d_bn) then
               d_l = d_ln
               d_r = d_rn
C
               rho_d = rho_dn
               p_d = p_dn
               un_d = un_dn
               T_d = T_dn
               d_d = d_dn
C
               rho_c = rho_cn
               p_c  = p_cn
               un_c = un_cn
               T_c = T_cn
               d_c = d_cn
C
               rho_b = rho_bn
               p_b = p_bn
               un_b = un_bn
               T_b = T_bn
               d_b = d_bn
C
               rho_a = rho_an
               p_a = p_an
               un_a = un_an
               T_a = T_an
               d_a = d_an
C
C               write(*,*) 'Sono piu grande del Masakino !'
C               write(*,*) 'coefd , k0 =', coefd, k0n
C               write(*,*) 'do =', d_o
C               write(*,*) 'dn =', d_bn
C               write(*,*) 'dm =', d_bm
C               write(*,*) 'cnna =', cnna
CC
C               k0o = 1.0D-3 * k0
C               do iter = 1, 11
C                  call racre2(nordpo, Tmaxcv,
C     &                 Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
C     &                 q, k0o,
C     &                 rho_l, p_l, un_l, T_l, gam_l,
C     &                 rho_r, p_r, un_r, T_r, gam_r,
C     &                 d_lo, d_ro,
C     &                 rho_do, p_do, un_do, T_do, d_do,
C     &                 rho_co, p_co, un_co, T_co, d_co,
C     &                 rho_bo, p_bo, un_bo, T_bo, d_bo,
C     &                 rho_ao, p_ao, un_ao, T_ao, d_ao,
C     &                 logan, lognc)
C
C                  if (logan) then
C                     write(*,'(/A8)') 'fbecre.f'
C                     write(*,*) 'ANOMALY DETECTED 10'
C                     goto 9999
C                  endif
C
C                  k0o = d_bo - un_ao
C                  write(*,*) k0o
C
C                  write(*,*) 'd_b0, k0o = ', d_bo, k0o
C                  k0o = k0o + 0.1D0 * (k0n - k0m)
C               enddo
CC
C               stop
C
            else
               k0o = ((d_o - d_bm) * (k0n - k0m) / (d_bn - d_bm)) + k0m
               call racre2(nordpo, Tmaxcv,
     &              Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &              q, k0o,
     &              rho_l, p_l, un_l, T_l, gam_l,
     &              rho_r, p_r, un_r, T_r, gam_r,
     &              d_lo, d_ro,
     &              rho_do, p_do, un_do, T_do, d_do,
     &              rho_co, p_co, un_co, T_co, d_co,
     &              rho_bo, p_bo, un_bo, T_bo, d_bo,
     &              rho_ao, p_ao, un_ao, T_ao, d_ao,
     &              logan, lognc)
C
               if (logan) then
                  write(*,'(/A8)') 'fbecre.f'
                  write(*,*) 'ANOMALY DETECTED 10'
                  goto 9999
               endif
C
               k0o = abs (d_bo - un_ao)
C
C               write(*,*) 'd_bo, k0o = ', d_bo, k0o
C
               do iter = 1, 10
                  if (d_bo .le. d_o) then
                     k0m = k0o
                     d_bm = d_bm
                  else
                     k0n = k0o
                     d_bn = d_bo
                  endif
C
                  k0o = ((d_o - d_bm) * (k0n - k0m) / (d_bn - d_bm)) +
     $                 k0m
                  call racre2(nordpo, Tmaxcv,
     &                 Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &                 q, k0o,
     &                 rho_l, p_l, un_l, T_l, gam_l,
     &                 rho_r, p_r, un_r, T_r, gam_r,
     &                 d_lo, d_ro,
     &                 rho_do, p_do, un_do, T_do, d_do,
     &                 rho_co, p_co, un_co, T_co, d_co,
     &                 rho_bo, p_bo, un_bo, T_bo, d_bo,
     &                 rho_ao, p_ao, un_ao, T_ao, d_ao,
     &                 logan, lognc)
C
                  if (logan) then
                     write(*,'(/A8)') 'fbecre.f'
                     write(*,*) 'ANOMALY DETECTED 10'
                     goto 9999
                  endif
C
                  k0o = d_bo - un_ao
C                  write(*,*) k0o
C
C                  write(*,*) 'd_b0, k0o = ', d_bo, k0o
C
                  d_l = d_lo
                  d_r = d_ro
C
                  rho_d = rho_do
                  p_d = p_do
                  un_d = un_do
                  T_d = T_do
                  d_d = d_do
C
                  rho_c = rho_co
                  p_c  = p_co
                  un_c = un_co
                  T_c = T_co
                  d_c = d_co
C
                  rho_b = rho_bo
                  p_b = p_bo
                  un_b = un_bo
                  T_b = T_bo
                  d_b = d_bo
C
                  rho_a = rho_ao
                  p_a = p_ao
                  un_a = un_ao
                  T_a = T_ao
                  d_a = d_ao
C
               enddo
C
            endif
            D_resh = d_b
            cellt = max (abs(d_l), abs(d_r))
C
         endif
C
C******* States
C
C         l
C        ***                       b
C          |*                      *                 r
C          | *   d                **         *********
C          |  *******            * *        *|
C          |  |     *   c       *  *       * |
C          |  |     ************   *  a   *  |
C          |  |     |          |   *******   |
C          |  |     |          |   |     |   |
C        d_l d_d  un_d=un_c   d_c d_b   d_a d_r
C
C
C
C******* Eulerian flux
C
         xst = 0.0D0
         call flufu1(nordpo,
     &        Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
     &        Rgas_l, acv_l, acv_w,
     &        xst,
     &        d_l, d_d, d_c, d_b, d_a, d_r,
     &        rho_l, T_l, un_l, hf_l,
     &        rho_d, T_d, un_d,
     &        rho_c, T_c, un_c,
     &        rho_b, T_b, un_b, hf_b,
     &        rho_a, T_a, un_a,
     &        rho_r, T_r, un_r, hf_r,
     &        flurho, fluru, fluret)
C
C        In order to obtain the flux in the true left-right frame, I have
C        to  multiply the speed * rind_l
         flurho = flurho * rind_l
         fluret = fluret * rind_l
         flu(1) = flurho
         flu(2) = fluru
         flu(2+ndim) = fluret
         flu(3+ndim) = 0.0d0
C        Note that the Larrouturou theorem does not hold
C        Namely un_d (= un_c) has not the same sign as flurho
C        write (*,*) un_d, un_c, flurho
C        un_ls > 0 => coefy = 1
         coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (rind_l * un_d)))
         omcoey = 1.0d0 - coefy
C        Tangential speeds (ndim .ge. 2)
         do icomp = 3, (ndim + 1), 1
            flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
     &           (omcoey * wvec_r(icomp)))
         enddo
         do icomp = (ndim + 4), (ndim + nesp + 4), 1
            flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
     &           (omcoey * wvec_r(icomp)))
         enddo
C        Kinetic energy flux
         trekin = 0.0D0
         do icomp = 2, ndim, 1
            trekin = trekin + (((coefy * wvec_l(icomp + 1)) +
     &           (omcoey * wvec_r(icomp + 1))) ** 2)
         enddo
         trekin = 0.5D0 * trekin
         flu(2+ndim) = flu(2+ndim) + (flurho * trekin)
C
C******* Lagrangian flux over the reactive shock
C
         xst = D_resh
         call flufu1(nordpo,
     &        Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
     &        Rgas_l, acv_l, acv_w,
     &        xst,
     &        d_l, d_d, d_c, d_b, d_a, d_r,
     &        rho_l, T_l, un_l, hf_l,
     &        rho_d, T_d, un_d,
     &        rho_c, T_c, un_c,
     &        rho_b, T_b, un_b, hf_b,
     &        rho_a, T_a, un_a,
     &        rho_r, T_r, un_r, hf_r,
     &        flurho, fluru, fluret)
C
         flurho = flurho * rind_l
         fluret = fluret * rind_l
         f_lag(1) = flurho
         f_lag(2) = fluru
         f_lag(2+ndim) = fluret
         f_lag(3+ndim) = -1.0D0 * rind_l * D_resh
C
C        Since we suppose that left state is the burnt one,
C        D_resh > un_d, i.e. over (x/t) = D_resh we have the
C        right passive scalar !
C        rind_l = 1   =>  burnt state in wvec_l => passive scalar in r
C        rind_l = -1  =>  burnt state in wvec_r => passive scalar in l
C        Note that the result does not change as q -> 0
C        (virtual combustion wave).
C
         coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (1.0D0 * rind_l)))
         omcoey = 1.0d0 - coefy
C
         do icomp = 3, (ndim + 1), 1
            f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
     &           (omcoey * wvec_l(icomp)))
         enddo
         do icomp = (ndim + 4), (ndim + nesp + 4), 1
            f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
     &           (omcoey * wvec_l(icomp)))
         enddo
C        Kinetic energy flux
         trekin = 0.0D0
         do icomp = 2, ndim, 1
            trekin = trekin + (((coefy * wvec_r(icomp + 1)) +
     &           (omcoey * wvec_l(icomp + 1))) ** 2)
         enddo
         trekin = 0.5D0 * trekin
C        write(*,*) trekin
         f_lag(2+ndim) = f_lag(2+ndim) + (flurho * trekin)
      endif
C
 9999 continue
      return
      end