C FUNREA    SOURCE    BECC      09/12/07    21:15:32     6579           
      subroutine funrea(nordp1, tmaxcv, rg_rss, acv_z, 
     &     rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues, 
     &     rho_z, p_rss, u_rss, Dsw, et_rss, logan)
c************************************************************************
c
c projet            :  
c
c name              :  funrea
c
c description       :  computation of the weak state behind a reactive 
c                      shock
c
c language          :  fortran 77 
c      
c author            :  a. beccantini den/dm2s/sfme/ltmf  
c
c************************************************************************
c
c callen by         :  fcolre
c
c called program    :  none
c
c************************************************************************
c
c
c**** input  
c
c             nordp1    = polynomials degree + 1
c             tmaxcv    = threshold temperature for cv(t)
c             rg_rss    = r of the gas in rss
c             acv_z     = \sum_{i=1,nesp} y_i*acv_z{i,j} 
C             rho_rs    = density in rs
C             p_rs      = pressure in rs
c             u_rs      = velocity in rs
c             et_rs     = (sensible) internal energy in rs
c             q         = realeased heat per mass unit
c             k0        = fundamental flame speed
c             t_gues    = guess temperature in rss around which we 
c                         linearised the internal energy (enthalpy)
c**** output
c             rho_z     = approximate density in rss (behind the shock)
c             p_rss     = approximate pressure in rss
c             u_rss     = approximate speed in rss
C             Dsw       = approximate shock wave speed... it could 
C                         differs from u_rs + k0 if (k0lim < k0)
c             cv_rss    = approximate cv in rss
c             et_rss    = approximate (sensible) internal energy in rs
c             logan     = an anomaly is detected
c
c************************************************************************
c
c     created the 07/12/2009
C
C      implicit none
      integer nordp1
      integer i1
      real*8 tmaxcv, rg_rss, acv_z(1:nordp1)
      real*8 rho_rs, p_rs, u_rs,  et_rs, q, k0, t_gues
      real*8 rho_z, p_rss, u_rss, cv_rss
      real*8 et_rss, funt, funcv, t2, ga_rss
      real*8 acoef, bcoef, ccoef, acoef1, bcoef1, ccoef1, k0lim2, k0lim
     &     ,Dsw, ratio
      logical logan
CC
CC**** We write the input
CC
C      write(*,*) 'We enter funrea'
C      write(*,'(I2,E12.4)') nordp1, tmaxcv
C      write(*,'(6E12.4)') rg_rss, (acv_z(icon), icon=1,nordp1)
C      write(*,'(7E12.4)') rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues
C
C**** Computation of ga_rss, ether_ss
C     
C     t2 = min (t_gues, Tmaxcv)
      t2 = 0.5d0 * ((t_gues + tmaxcv) - abs (t_gues - tmaxcv))
      cv_rss = acv_z(1)
      et_rss = t2 * cv_rss
      funt = 1.0d0
      do i1 = 2, nordp1, 1
         funt = funt * t2
         funcv = acv_z(i1) * funt
         cv_rss = cv_rss + funcv
         et_rss = et_rss + (funcv * t2/ i1)
      enddo
c     t2 = max ((t_gues - tmaxcv),0)
      t2 = 0.5d0 * ((t_gues - tmaxcv) + abs (t_gues - tmaxcv))
      et_rss = et_rss + (cv_rss * t2)
      ga_rss = (cv_rss + Rg_rss) / cv_rss
C
C**** Computation of rho_z
C
      k0lim2 = k0 * k0
      bcoef = ga_rss / (ga_rss + 1.0d0)
      bcoef = bcoef * (1.0d0 + (p_rs / (rho_rs * k0lim2)))
C      
      ccoef = (p_rs / rho_rs) + et_rs + q 
      ccoef = ccoef - (et_rss - (cv_rss * t_gues))
      ccoef = (ccoef * 2.0D0 / k0lim2) + 1.0d0
      ccoef = ccoef * (ga_rss - 1.0d0) / (ga_rss + 1.0d0)
      acoef = (bcoef * bcoef) - ccoef
      if (acoef .lt. 0.0d0) then
C
C******* k0 is too big.
C        Computation of K0lim2 = k0lim * k0lim corresponding to
C        CJ deflagration
C
         acoef1 = (ga_rss * (p_rs / rho_rs)) / (ga_rss + 1.0d0)
         acoef1 = acoef1 * acoef1
         bcoef1 = ga_rss / (ga_rss + 1.0d0)
         bcoef1 = bcoef1 * bcoef1 * (p_rs / rho_rs)
         ccoef1 = (p_rs / rho_rs) + et_rs + q 
         ccoef1 = ccoef1 - (et_rss - (cv_rss * t_gues))
         ccoef1 = ccoef1 * (ga_rss - 1.0d0) / (ga_rss + 1.0d0) 
         bcoef1 = ccoef1 - bcoef1
         ccoef1 = 1.0d0 / ((ga_rss + 1.0d0) ** 2)
         k0lim2 = (bcoef1 * bcoef1) - (acoef1 * ccoef1)
         if ( (k0lim2 .lt. 0.0d0) ) then
C
C           Discrimant zero... quite bizarre for exothermic reaction...
C           it can happen but for numerical reasons... 
C           in that case, the best is to force it to zero!!!
C
            logan = .true.
            write(*,'(/A8)') 'funrea.f'
            write(*,*) 'Negative discrimant for k0lim'
            write(*,*) 'delta ', k0lim2
            write(*,*) 'ANOMALY DETECTED'
            goto 9999
         endif
         k0lim2 = k0lim2 ** 0.5d0
C        Before this line, k0lim2 is just the square root of the discrimant
         k0lim2 = (bcoef1 + k0lim2) / acoef1
         k0lim2 = 1.0d0 / k0lim2
         bcoef1 = ga_rss / (ga_rss + 1.0d0)
         bcoef1 = bcoef1 * (1.0d0 + (p_rs / (rho_rs * k0lim2)))
         ratio = bcoef1
      else
         acoef1 = acoef ** 0.5d0
         ratio = bcoef - acoef1
      endif
C
      rho_z = rho_rs / ratio
      if ( (rho_z .lt. 0.0d0) ) then
C
C        Negative density
C
         logan = .true.
         write(*,'(/A8)') 'funrea.f'
C         write(*,*)  'rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues'
C         write(*,'(7E16.6)')  rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues
C         write(*,*) 'rho_z, p_rss, u_rss, Dsw, et_rss'
C         write(*,'(5E16.6)') rho_z, p_rss, u_rss, Dsw, et_rss
         write(*,*) 'Negative density'
         write(*,*) 'rho_z', rho_z
         write(*,*) 'ANOMALY DETECTED'
         goto 9999
      endif
      k0lim = k0lim2 ** 0.5d0
      p_rss = p_rs + (rho_rs * k0lim2 * (1.0d0 - ratio))
      Dsw = k0lim + u_rs
      u_rss = u_rs + ((p_rss - p_rs) / (k0lim * rho_rs))
CC
C      write(*,*) 'rho_z, p_rss, u_rss, Dsw, et_rss'
C      write(*,'(5E12.4)') rho_z, p_rss, u_rss, Dsw, et_rss
C      write(*,*) 'We exit funrea'
C
 9999 continue
      return
      end