Sources Esope | Cast3M

funre1
C FUNRE1    SOURCE    BECC      09/12/07    21:15:31     6579
      subroutine funre1(nordpo, tmaxcv, rgas_a, acv_a,
     &     rgas_b, acv_b, q, k0, T_ag, T_bg,
     &     r_a, p_a, u_a, T_a,
     &     r_b, p_b, u_b, T_b, Dsw, logan)
c************************************************************************
c
c projet            :
c
c name              :  funre1
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         :  racre2
c
c called program    :  none
c
c************************************************************************
c
c
c**** input
c
c             nordpo    = polynomials degree
c             tmaxcv    = threshold temperature for cv(t)
c             rgas_a    = r of the gas in r and rs
c             acv_a     = \sum_{i=1,nesp} y_i*acv_a{i,j}
c             rgas_b    = r of the gas in rss
c             acv_b     = \sum_{i=1,nesp} y_i*acv_b{i,j}
c             q         = realeased heat per mass unit
c             k0        = fundamental flame speed
C             T_ag,T_bg = guess temperatures for a and b
C                         (to compute properties)
C             r_a       = density in rs
C             p_a       = pressure in rs
c             u_a       = velocity in rs
C             T_a       = temperature in rs
c**** output
c             r_b       = approximate density in rss (behind the shock)
c             p_b       = approximate pressure in rss
c             u_b       = approximate speed in rss
C             T_b       = approximate temperature in rss
C                         (NB if T_bg = T_b and T_ag = T_a, rss thus
C                          computed is exact)
C             Dsw       = approximate shock wave speed... it could
C                         differs from u_a + k0 if (k0lim &lt; k0)
C                         Nevertheless, if logsto = .true., if k0lim
C                         &lt; k0 we stop!
c             logan     = an anomaly is detected
c
c************************************************************************
c
c     created the 07/12/2009
c
c      implicit none
      integer nordpo
      real*8 tmaxcv, rgas_a, acv_a(1:(nordpo+1))
     &     ,rgas_b, acv_b(1:(nordpo+1))
      real*8 q, k0
      real*8 T_ag, T_bg
      real*8 r_a, p_a, u_a,  T_a
      real*8 r_b, p_b, u_b, T_b, Dsw
      real*8 gam_a, e_a, cv_a, gam_b, e_b, cv_b
      real*8 ttq, k02ad, h_a
      real*8 a, b, c, delta, delta1, k0cj2, k0new
      real*8 csi
      logical logan
CC
CC**** We write the input
CC
C      write(*,*) 'We enter funre1'
C
C**** Computation of gam_a and gam_b
C
      call prith1(nordpo, acv_a, Tmaxcv, T_ag, e_a, cv_a)
      gam_a = (cv_a + Rgas_a) / cv_a
      call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b)
      gam_b = (cv_b + Rgas_b) / cv_b
      ttq = q + ((cv_b * T_bg) - e_b) -
     &     ((cv_a * T_ag) - e_a)
C     h_a and h_b computed as we deal with cp gases
      h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_a * T_a
C
      r_a = P_a / (Rgas_a * T_a)
C
C**** Computation of r_b
C
      k02ad =  k0 * k0 / (Rgas_a * T_a)
      if (k02ad .lt. 1.0D0)then
C      if ( .true. ) then
C      if ( .false. ) then
C         write(*,*) 'funre1.f. Parte 1'
         c = k02ad
         b = k02ad + 1.0D0
         b = b * gam_b / (gam_b + 1.0D0)
         a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
         a = a + k02ad
         a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
         delta = (b * b) - (a * c)
         if (delta .lt. 0.0D0) then
C            if ( .true. ) then
C            write(*,*) 'funre1.f. Parte 1.2'
C           We are in CJDF regime
            a = 1.0D0 / (gam_b * gam_b)
            b = 1.0D0 - a
            b = b * (h_a + ttq)
            b = b - (Rgas_a * T_a)
            c = (Rgas_a * T_a)**2
            delta1 = (b * b) - (a * c)
            if (delta1 .lt. 0.0D0) then
               write(*,*) 'subroutine funre1.f'
               write(*,*) 'We cannot compute k0cj'
               write(*,*) 'ANOMALY DETECTED'
               logan = .true.
               goto 9999
            endif
            k0cj2 = (b - sqrt(delta1)) / a
            if (k0cj2 .lt. 0.0D0) then
               write(*,*) 'subroutine funre1.f'
               write(*,*) 'k0cj^2 &lt; 0'
               write(*,*) 'ANOMALY DETECTED'
               logan = .true.
               goto 9999
            endif
            k02ad = k0cj2 / (Rgas_a * T_a)
            b = k02ad + 1.0D0
            b = b * gam_b / (gam_b + 1.0D0)
            a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
            a = a + k02ad
            a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
            csi = b / a
            k0new = k0cj2 ** 0.5D0
         else
            k0new = k0
            csi = b + sqrt(delta)
            csi = csi / a
         endif
         r_b = r_a * csi
      else
C         write(*,*) 'funre1.f. Parte 2'
         c = 2.0D0 * (h_a + ttq)
         c = c / (k0 * k0)
         c = c + 1.0D0
         c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
         b = Rgas_a * T_a / (k0 * k0)
         b = b + 1.0D0
         b = b * gam_b / (gam_b + 1.0D0)
         delta = (b * b) - c
         if (delta .lt. 0.0D0) then
C         if ( .true. ) then
C            write(*,*) 'funre1.f. Parte 2.2'
            a = 1.0D0 / (gam_b * gam_b)
            b = 1.0D0 - a
            b = b * (h_a + ttq)
            b = b - (Rgas_a * T_a)
            c = (Rgas_a * T_a)
            c = c * c
            delta1 = (b * b) - (a * c)
            if (delta1 .lt. 0.0D0) then
               write(*,*) 'subroutine funre1.f'
               write(*,*) 'We cannot compute k0cj'
               write(*,*) 'ANOMALY DETECTED'
               logan = .true.
               goto 9999
            endif
            k0cj2 = (b - sqrt(delta1)) / a
            if (k0cj2 .lt. 0.0D0) then
               write(*,*) 'subroutine funre1.f'
               write(*,*) 'k0cj^2 &lt; 0'
               write(*,*) 'ANOMALY DETECTED'
               logan = .true.
               goto 9999
            endif
            c = 2.0D0 * (h_a + ttq)
            c = c / k0cj2
            c = c + 1.0D0
            c = c *  (gam_b - 1.0D0) / (gam_b + 1.0D0)
            b = Rgas_a * T_a / k0cj2
            b = b + 1.0D0
            b = b * gam_b / (gam_b + 1.0D0)
            delta1 = b * b - c
C            if (abs (delta1) > 1.0D-6) then
            if (abs (delta1) .gt. 1.0D-6) then
               write(*,*) 'subroutine funre1.f'
               write(*,*) 'CJDF'
               write(*,*) 'delta != 0'
               write(*,*) 'ANOMALY DETECTED'
               logan = .true.
               goto 9999
            endif
            csi = b
            k0new = k0cj2 ** 0.5D0
         else
            csi = b - sqrt (delta)
            k0new = k0
         endif
         r_b = r_a / csi
      endif
      csi = r_a / r_b
      Dsw = k0new + u_a
      u_b = Dsw - (k0new * csi)
      p_b = p_a + (r_a * k0new * k0new * (1.0d0 - csi))
      if (p_b .lt. 0.0D0) then
         write(*,*) 'subroutine funre1.f'
         write(*,*) 'p_b &lt; 0'
         write(*,*) 'ANOMALY DETECTED'
         logan = .true.
         goto 9999
      endif
      T_b = P_b / (Rgas_b * r_b)
C      write(*,*) 'r_b, p_b, u_b, T_b, Dsw, k0new'
C      write(*,'(6E12.4)') r_b, p_b, u_b, T_b, Dsw
C      write(*,*) r_b, p_b, u_b, T_b, Dsw, k0new
C      stop
CC      write(*,*) 'We exit funre1'
CC
 9999 continue
      return
      end