Sources Esope | Cast3M

racre2
C RACRE2    SOURCE    BECC      09/12/07    21:16:00     6579
      subroutine racre2(nordpo, Tmaxcv,
     &     Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
     &     q, k0,
     &     r_l, p_l, u_l, T_l, gam_l,
     &     r_r, p_r, u_r, T_r, gam_r,
     &     d_l, d_r,
     &     r_d, p_d, u_d, T_d, d_d,
     &     r_c, p_c, u_c, T_c, d_c,
     &     r_b, p_b, u_b, T_b, d_b,
     &     r_a, p_a, u_a, T_a, d_a,
     &     logan, lognc)
C
c************************************************************************
c
c project           :
c
c name              :  racre2
c
c description       :  computation of the shock-shock field-by-field
c                      decomposition of the reactive Riemann problem for
c                      thermally perfect gases.
C
C                      We have to distinguish between 6 different states
C
C                      l    = left state;
C                             burnt state (csi = csimax)
C                                         |
C                                         | Left GNL wave
C                                         |
C                      d    = left star state;
C                             connected to l by a GNL wave
C                             (csi = csimax)
C                                         |
C                                         | Contact discontinuity
C                                         |
C                      c    = right star star star state;
C                             connected to rss state by a rarefaction
C                             wave
C                             (csi = csimax)
C                                         |
C                                         | Rarefaction wave (weak defl.)
C                                         |
C                      b    = right star star state;
C                             connected to rs by a reactive wave
C                             (csi = csimax)
C                                         |
C                                         | Reactive shock (weak defl.)
C                                         |
C                      a    = right star state;
C                             connected to r by a GNL wave
C                             (csi = csimin)
C                                         |
C                                         | Right GNL wave
C                                         |
C                      r    = right state;
C                             unburnt state (csi = csimin)
c
c language          :  fortran 77
c
c author            :  a. beccantini den/dm2s/sfme/ltmf
c
c************************************************************************
c
c called by         :  fbecre
c
c called program    :  ...
c
c************************************************************************
c
c
c**** input
c
c             nordpo    = polynomials degree of cv(t)
c             tmaxcv    = threshold temperature for cv(t)
c             rgas_l    = gas constant (J/kg/K)in l (and in ls)
c             acv_l     = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i)
c             rgas_rss  = gas constant in rss
c             acv_z   = cv_rss = \sum_{i=0,nordpo} (acv_z{i} T^i)
c             rgas_r    = gas constant r (and in rs)
c             acv_r     = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i)
c
c             q         = released heat per mass unit from rs to rss
c             k0        = fundamental flame speed
c
C             r_l       = density in l
C             p_l       = pressure in
c             u_l       = velocity in l
C             T_l       = temperature in l (redundant with r_l and
c                         p_l)
c             gam_l     = specific heat ratio in l (redundant
c                         with r_l and p_l)
c
C             r_r       = density in r
C             p_r       = pressure in r
c             u_r       = velocity in r
C             T_r       = temperature in r (redundant with r_r and
c                         p_r)
c             gam_r     = specific heat ratio in r (redundant
c                         with r_r and p_r)
c
c**** output
C
C             d_r, d_l  = speed of states l and r
c
C             r_d       = density in d
C             p_d       = pressure in d
c             u_d       = velocity in d
c             T_d       = temperature in d
c             d_d       = shock speed of the left shock l&lt;-d
c
C             r_c       = density in c
C             p_c       = pressure in c
c             u_c       = velocity in c
c             T_c       = temperature in c
c             d_c       = shock speed of the (non-entropy respecting)
C                         right shock c->b
c
C             r_b       = density in b
C             p_b       = pressure in b
c             u_b       = velocity in b
c             T_b       = temperature in b
c             d_b       = shock speed of the reactive shock
C                         right shock b->a
c
C             r_a       = density in a
C             p_a       = pressure in a
c             u_a       = velocity in a
c             T_a       = temperature in a
c             d_a       = shock speed of the
C                         right shock a->l
c
c             logan     = an anomaly is detected
c             lognc     = the approach has some problems
C
c
c************************************************************************
c
c     07/12/2009      Created.
c
c      implicit none
      integer nordpo
      integer nmax
      parameter (nmax = 20)
      integer iter, indreg
      real*8 eps
      parameter(eps = 1.0D-6)
      real*8  Tmaxcv, Rgas_l, Rgas_b, Rgas_r
     &     , acv_l(0:nordpo), acv_b(0:nordpo), acv_r(0:nordpo)
     &     , q, k0
     &     , r_l, p_l, u_l, T_l, gam_l
     &     , r_r, p_r, u_r, T_r, gam_r
     &     , d_l, d_r
     &     , r_d, p_d, u_d, T_d, d_d
     &     , r_c, p_c, u_c, T_c, d_c
     &     , r_b, p_b, u_b, T_b, d_b
     &     , r_a, p_a, u_a, T_a, d_a
      real*8 k0cjdt, k0sdt
     &     , d_acj, d_bcj
     &     , r_acj, T_acj, p_acj, u_acj
     &     , r_bcj, T_bcj, p_bcj, u_bcj
     &     , P_asdt, r_asdt, u_asdt, d_sdt
     &     , T_csdt, P_csdt, r_csdt, u_csdt
     &     , T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
     &     , dt1, dt2
C
      real*8 e_b, cv_b
      real*8 T_ag, T_bg, T_ag1, T_bg1, T_cg, T_dg
      real*8 Tref, pref, Dtmin, coefder
      real*8 P_ag, r_ag, u_ag, d_ag
      real*8 r_bg, p_bg, u_bg, d_bg
      real*8 r_cg, p_cg, u_cg, d_cg
      real*8 T_agp, P_agp, r_agp, u_agp, d_agp
      real*8 T_agpp, P_agpp, r_agpp, u_agpp, d_agpp
      real*8 r_bgp, p_bgp, u_bgp, T_bgp, d_bgp
      real*8 r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp
      real*8 dtudtb
      real*8 dpduc, dpdud
      real*8 P_dg, r_dg, u_dg, d_dg
      real*8 T_dgp, P_dgp, r_dgp, u_dgp, d_dgp
      real*8 a, b
      logical lognc, logan, logdef, logde1, logdeb
      real*8 erro
      parameter(logdeb = .false.)
CC
CC**** We write the inputs
CC
C      write(*,*) 'racre2.f'
C      write(*,'(I2,E12.4)') nordpo, tmaxcv
C      write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo)
C      write(*,'(6E12.4)') rgas_b, (acv_b(icon), icon=0,nordpo)
C      write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo)
C      write(*,'(2E12.4)') q, k0
C      write(*,'(6E12.4)') r_l, p_l, u_l, T_l, gam_l
C      write(*,'(6E12.4)') r_r, p_r, u_r, T_r, gam_r
C      stop
C
C**** Reference scales for a compressible solver.
C     For a low Mach number solver, the cut off speed is better than
C     the sound speed
C
      Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r))
      Dtmin = 1.0D-3 * Tref
      coefder = 1.001D0
      pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r))
C
C     Initialisation of temperatures in a, b, c, d
C
      T_a = T_r
      T_d = T_l
C
C     T_b and T_c computed using AIbCC
C     i.e.
C     h_b = h_r + q
c     =>
      call prith1(nordpo, acv_b, Tmaxcv, T_l, e_b, cv_b)
      T_b = gam_r / (gam_r - 1.0D0) * Rgas_r * T_r
      T_b = T_b + q
      T_b = T_b / (cv_b + Rgas_b)
      T_c = T_b
C
C
C     TEST
      if (logdeb) then
         write(*,*) ' '
         write(*,*) 'racre2.f'
         write(*,*) 'Guess values'
         write(*,*) 'T_l, T_d, T_c, T_b, T_a, T_r '
         write(*,*) T_l,  T_d, T_c, T_b, T_a, T_r
      endif
C
      do iter = 1, nmax
C
C        Guess temperatures at this iteration
C
         T_ag = T_a
         T_bg = T_b
         T_cg = T_c
         T_dg = T_d
C
C        states g1 to compute the properties (gamma)
C
         T_ag1 = T_a
         T_bg1 = T_b
C
         if (logdeb) then
            write(*,*) ' '
            write(*,*) 'iter = ', iter
            write(*,*) 'T_ag, T_bg, T_cg, T_dg ='
            write(*,*) T_ag, T_bg, T_cg, T_dg
            write(*,*) 'T_ag1, T_bg1 ='
            write(*,*) T_ag1, T_bg1
         endif
C
C******* Right  state
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC Test for funcj state CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         acv computed for a mixture of h2-air at stoichiometric
C         condition (burnt vs unburn)
C         if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
C         write(*,*) 'q =', q
C         p_r = 1.0D5
C         t_r = 290.0D0
C         u_r = 20.0D0
C         q = 3376364.59D0
C         k0 = 1000.0D0
C         T_ag1 = 1726.278895493752D0
C         T_ag = T_ag1
C         T_bg1 = 3519.399539515807D0
C         T_bg = T_bg1
C         T_cg = T_bg1 * 1.01D0
CC
C        a) With these values and k0 = 1000 -> CJDT states
C           r_acj, T_acj, p_acj, u_acj, d_acj
C           4.84002641182335  1665.81525651246
C           3205701.36057718  1734.45181775933  2108.72324760979
C           r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
C           1.53143449865192  3355.53701100291  1740943.41687983
C           925.856068333306  2108.72324760979
C           k0cjdt = 374.271429850455
C
C        b) With these values and k0 = 100 -> CJDF states
C           r_acj, T_acj, p_acj, u_acj, d_acj
C           1.89601907217341  420.523548481071  317016.244413484
C           388.470810477263  699.102996003551
C           r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
C           0.193932466089758 2292.31572201026  150608.384363281
C           -489.198947662085 488.470810477263
C           k0cjdt = 374.271429850455
C
C        c) With these values and k0 = 50 -> CJDF states
C           r_acj, T_acj, p_acj, u_acj, d_acj
C           0.428757456680544  207.721167691767  35411.2237965742
C           -255.985416642147  85.8302953725479
C           r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
C           2.283207457276537E-002  2114.28097729292 16354.3117888704
C           -1144.92212016881  -205.985416642147
C           k0cjdt =    374.271429850455
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC Test for funsdt state CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C         acv computed for a mixture of h2-air at stoichiometric
C         condition (burnt vs unburn)
C         if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
C         write(*,*) 'q =', q
C         p_r = 1.0D5
C         t_r = 290.0D0
C         u_r = 20.0D0
CC         q = 3376364.59D0
C         k0 = 1000.0D0
C         T_ag1 = 1726.278895493752D0
C         T_ag = 1800.0D0
C         T_bg1 = 3519.399539515807D0
C         T_bg = T_bg1
C         T_cg = T_bg1 * 1.2D0
C
C        With this value, for funcj we find the same values as for
C        case a.
C        Moreover, for funsdt, we find
C        T_r, P_r, u_r
C        290.000000000000     100000.000000000        20.0000000000
C        T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt
C        1800.00000000000     3547944.85591171        4.95742167645132
C        1831.11181823155        2215.13654202499        384.024723793441
C        T_csdt, P_csdt, r_csdt, u_csdt
C        3625.49600977864 2494329.73747622 2.03077700146451
C        1277.67640304670
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC Test for funre1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC        acv computed for a mixture of h2-air at stoichiometric
CC        condition (burnt vs unburn)
C         if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
C         write(*,*) 'q =', q
C         p_r = 1.0D5
C         t_r = 290.0D0
C         u_r = 20.0D0
CCC         q = 3376364.59D0
C         k0 = 100.0D0
C         T_ag1 = 1726.278895493752D0
CC        T_ag = 420.7D0
C         T_ag = 600.7D0
C         T_bg1 = 3519.399539515807D0
C         T_bg = T_bg1
C         T_cg = T_bg1 * 1.2D0
CC
CC        With these values
CC
CC        WDF
CC        T_r, P_r, u_r
CC        290.000000000000  100000.000000000  20.0000000000000
CC        T_ag, P_ag, r_ag, u_ag, d_ag
CC        600.700000000000  688354.959139395  2.88208190644972
CC        708.663754615697  1005.09629332136
CC        T_bg, P_bg, r_bg, u_bg, d_bg
CC        2669.18760214319  589828.308238861  0.652262358864745
CC        366.804452560316  808.663754615697
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC Test for fundt  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
CC        acv computed for a mixture of h2-air at stoichiometric
CC        condition (burnt vs unburn)
C         if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
C         write(*,*) 'q =', q
C         p_r = 1.0D5
C         t_r = 290.0D0
C         u_r = 20.0D0
CC         q = 3376364.59D0
C         k0 = 1000.0D0
C         T_ag1 = 1726.278895493752D0
C         T_ag = 600.7D0
C         T_bg1 = 3519.399539515807D0
C         T_bg = T_bg1
C         T_cg = 4000.0D0
CC
CC        With these values
CC
CC        SDT
CC        T_r, P_r, u_r
CC        290.000000000000  100000.000000000  20.0000000000000
CC        T_ag, P_ag, r_ag, u_ag
CC        2133.61520897134  4396903.39501966  5.18301338812827
CC        2051.13020586146
CC        T_cg, P_cg, r_cg, u_cg, d_cg
CC        4000.00000000000  3553808.96027808  2.62246614378061
CC        1652.60261159855  2459.29495487122
CC        T_agp, P_agp, r_agp, u_agp
CC        2137.48514727336  4406760.43821584  5.18522780966239
CC        2053.54581926517
CC        T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
CC        4004.00000000000  3565189.78187338  2.62823616977283
CC        1656.17695829123  2461.98642456893
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC END OF TEST CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C        Computation of CJDF/CJDT states
C
         call funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
     &        k0, q,
     &        T_ag1, T_bg1,
     &        p_r, T_r, u_r,
     &        logdef,
     &        r_acj, T_acj, p_acj, u_acj, d_acj,
     &        r_bcj, T_bcj, p_bcj, u_bcj, d_bcj,
     &        k0cjdt,
     &        logan)
C
C             Now I know the CJ state.
C             If logdef, CJ is the CJDF state.
C             Otherwise CJ is the CJDT state.
C
         if (logan) then
            write(*,*) 'subroutine racre2.f'
            write(*,*) 'an error has been detected'
            goto 9999
         endif
         if (logdeb) then
            write(*,*) 'CJ state'
            write(*,*) 'Deflagration regime = ', logdef
            write(*,*) 'r_acj, T_acj, p_acj, u_acj, d_acj'
            write(*,*) r_acj, T_acj, p_acj, u_acj, d_acj
            write(*,*) 'r_bcj, T_bcj, p_bcj, u_bcj, d_bcj'
            write(*,*) r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
         endif
CC
CC        Test
CC
C         write(*,*) 'deflagration ', logdef
C         write(*,*) ' r_acj, T_acj, p_acj, u_acj, d_acj '
C         write(*,*) r_acj, T_acj, p_acj, u_acj, d_acj
C         write(*,*) ' r_bcj, T_bcj, p_bcj, u_bcj, d_bcj '
C         write(*,*) r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
C         write(*,*) 'k0cjdt = ', k0cjdt
CC         stop
C
         if (T_cg .gt. T_bcj) then
C
C********** SDT/WDF
C
            T_bg = max(T_bg, T_bcj)
            T_ag = max(T_ag, T_acj)
C
            if (logdeb ) then
               write(*,*) 'WDF/SDT regime'
               write(*,*) 'We change T_ag, T_bg'
               write(*,*) T_ag, T_bg
            endif
C
            logde1 = .true.
C
            if( .not. logdef )then
C              We are here if k0 > k0cjdt
C              In this case T_ag >= T_a,cjdt
C              In this case T_bg >= T_b,cjdt
C              NB. We can have k0 > k0cjdt and being in WDF regime
C
C************* Computation of the SDT state (necessary to select the
C              regime) as function of T_ag
C
               call funsdt(nordpo, Tmaxcv, acv_r, Rgas_r,
     &              acv_b, Rgas_b, q, T_ag1, T_bg1,
     &              P_r, T_r, u_r,
     &              T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt,
     &              T_csdt, P_csdt, r_csdt, u_csdt,
     &              logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
CC
CC              Test. We write the values thus determined
CC
C               write(*,*) 'T_r, P_r, u_r'
C               write(*,*) T_r, P_r, u_r
C               write(*,*) 'T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt'
C               write(*,*) T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt
C               write(*,*) 'T_csdt, P_csdt, r_csdt, u_csdt'
C               write(*,*) T_csdt, P_csdt, r_csdt, u_csdt
CC               stop
               if (k0 .lt. k0sdt) then
                  logde1 = .true.
               else
                  logde1 = .false.
               endif
            endif
C
            if (logde1) then
               indreg = 1
C
C************* WDF
C
C              Computation of the state ag
C
               call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
     &              P_r, T_r, u_r, .true.,
     &              T_ag, P_ag, r_ag, u_ag, d_ag)
C
C
C              Computation of the state bg (we change the value of T_bg)
C
               call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
     &              r_ag, p_ag, u_ag, T_ag,
     &              r_bg, p_bg, u_bg, T_bg, d_bg, logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C
               if (logdeb) then
                  write(*,*) 'WDF'
                  write(*,*) 'Parameter T_ag = ', T_ag
                  write(*,*)
     &                 'New values for T_bg = T_cg = ', T_bg
               endif
C
CC
CC              Test. We write the values thus determined
CC
C               write(*,*) 'WDF'
C               write(*,*) 'T_r, P_r, u_r'
C               write(*,*) T_r, P_r, u_r
C               write(*,*) 'T_ag, P_ag, r_ag, u_ag, d_ag'
C               write(*,*) T_ag, P_ag, r_ag, u_ag, d_ag
C               write(*,*) 'T_bg, P_bg, r_bg, u_bg, d_bg'
C               write(*,*) T_bg, P_bg, r_bg, u_bg, d_bg
C               stop
C
               if (d_bg .gt. d_ag) then
                  if (logdeb) then
                     write(*,*) 'WDF'
                     write(*,*) 'd_bg,  d_ag = ', d_bg, d_ag
                     write(*,*) 'subroutine racre2.f'
                     write(*,*) 'an error has been detected'
                  endif
C                  logan = .true.
C                  goto 9999
               endif
C
C              cg = bg
C
               r_cg = r_bg
               p_cg = p_bg
               u_cg = u_bg
               T_cg = T_bg
               d_cg = d_bg
C
               T_agp = max (T_ag + Dtmin, (T_ag * coefder))
               T_agpp = 2.0D0 * T_agp - T_ag
C
C              Computation of states agp and agpp
C
               call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
     &              P_r, T_r, u_r, .true.,
     &              T_agp, P_agp, r_agp, u_agp, d_agp)
               call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
     &              P_r, T_r, u_r, .true.,
     &              T_agpp, P_agpp, r_agpp, u_agpp, d_agpp)
C
               call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
     &              r_agp, p_agp, u_agp, T_agp,
     &              r_bgp, p_bgp, u_bgp, T_bgp, d_bgp, logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C
               call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
     &              r_agpp, p_agpp, u_agpp, T_agpp,
     &              r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp, logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C
               r_cgp = r_bgp
               p_cgp = p_bgp
               u_cgp = u_bgp
               T_cgp = T_bgp
CC
C               write(*,*) 'WDF'
C               write(*,*) 'T_r, P_r, u_r'
C               write(*,*) T_r, P_r, u_r
C               write(*,*) 'T_agp, P_agp, r_agp, u_agp, d_agp'
C               write(*,*) T_agp, P_agp, r_agp, u_agp, d_agp
C               write(*,*) 'T_bgp, P_bgp, r_bgp, u_bgp, d_bgp'
C               write(*,*) T_bgp, P_bgp, r_bgp, u_bgp, d_bgp
C               write(*,*) 'T_agpp, P_agpp, r_agpp, u_agpp, d_agpp'
C               write(*,*) T_agpp, P_agpp, r_agpp, u_agpp, d_agpp
C               write(*,*) 'T_bgpp, P_bgpp, r_bgpp, u_bgpp, d_bgpp'
C               write(*,*) T_bgpp, P_bgpp, r_bgpp, u_bgpp, d_bgpp
C
C              Second order formula (derivative evaluated in gp)
C
               dTudTb = (T_agpp - T_ag) / (T_bgpp - T_bg)
C
C               write(*,*) 'dTudTb = ', dTudTb
C               stop
            else
               indreg = 4
C
C************* SDT
C              Free parameter T_cg
C              Here we change the value of T_ag and T_bg
C
               call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
     &              acv_b, Rgas_b, q, T_ag1,
     &              P_r, T_r, u_r,
     &              T_cg, P_cg, r_cg, u_cg, d_cg, k0sdt,
     &              T_ag, P_ag, r_ag, u_ag,
     &              logan)
C
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C
               if (logdeb) then
                  write(*,*) 'SDT'
                  write(*,*) 'Parameter T_cg = ', T_cg
                  write(*,*) 'New values for T_ag = ', T_ag
                  write(*,*) 'T_ag, P_ag, r_ag, u_ag'
                  write(*,*) T_ag, P_ag, r_ag, u_ag
                  write(*,*) 'T_cg, P_cg, r_cg, u_cg'
                  write(*,*) T_cg, P_cg, r_cg, u_cg
                  write(*,*) 'd_cg, k0sdt'
                  write(*,*) d_cg, k0sdt
               endif
C
               T_bg = T_cg
               P_bg = P_cg
               r_bg = r_cg
               u_bg = u_cg
               d_bg = d_cg
               d_ag = d_cg
C
               T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
               call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
     &              acv_b, Rgas_b, q, T_ag1,
     &              P_r, T_r, u_r,
     &              T_cgp, P_cgp, r_cgp, u_cgp, d_cgp, k0sdt,
     &              T_agp, P_agp, r_agp, u_agp,
     &              logan)
C
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C
               T_bgp = T_cgp
               P_bgp = P_cgp
               r_bgp = r_cgp
               u_bgp = u_cgp
               d_bgp = d_cgp
               d_agp = d_cgp
CC
CC              Test. We write the values thus determined
CC
C               write(*,*) 'SDT'
C               write(*,*) 'T_r, P_r, u_r'
C               write(*,*) T_r, P_r, u_r
C               write(*,*) 'T_ag, P_ag, r_ag, u_ag'
C               write(*,*) T_ag, P_ag, r_ag, u_ag
C               write(*,*) 'T_cg, P_cg, r_cg, u_cg, d_cg'
C               write(*,*) T_cg, P_cg, r_cg, u_cg, d_cg
C               write(*,*) 'T_agp, P_agp, r_agp, u_agp'
C               write(*,*) T_agp, P_agp, r_agp, u_agp
C               write(*,*) 'T_cgp, P_cgp, r_cgp, u_cgp, d_cgp'
C               write(*,*) T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
C               stop
C
            endif
         else
C
C********** CJDT/CJDF
C
            r_ag = r_acj
            T_ag = T_acj
            P_ag = P_acj
            u_ag = u_acj
            d_ag = d_acj
C
            r_bg = r_bcj
            T_bg = T_bcj
            P_bg = P_bcj
            u_bg = u_bcj
            d_bg = d_bcj
C
            if (logdef) then
C
C              CJDF
C
               indreg = 2
C
               if (logdeb) then
                  write(*,*) 'CJDF'
                  write(*,*) 'Parameter T_cg = ', T_cg
                  write(*,*) 'New values for T_ag = ', T_ag
               endif
C              T_cg &lt;= T_bg
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_bg, T_bg, u_bg, .true.,
     &              T_cg, P_cg, r_cg, u_cg, d_cg)
C
C              We take T_cgp > T_cg since a lower value gives problem
C              as T_cg goes to zero.
C              Note that as T_c changes, T_a and T_b do not
C
               T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_bg, T_bg, u_bg, .true.,
     &              T_cgp, P_cgp, r_cgp, u_cgp, d_cgp)
C
C              Even in this case, I need dTudTb and states
C              agp, agpp, bgp, bgpp
C
               T_agp = max((T_ag + Dtmin), (T_ag*coefder))
               T_agpp = (2.0D0 * T_agp) - T_ag
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_r, T_r, u_r, .true.,
     &              T_agp, P_agp, r_agp, u_agp, d_agp)
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_r, T_r, u_r, .true.,
     &              T_agpp, P_agpp, r_agpp, u_agpp, d_agpp)
C
               call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
     &              r_agp, p_agp, u_agp, T_agp,
     &              r_bgp, p_bgp, u_bgp, T_bgp, d_bgp, logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
               call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
     &              r_agpp, p_agpp, u_agpp, T_agpp,
     &              r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp, logan)
               if (logan) then
                  write(*,*) 'subroutine racre2.f'
                  write(*,*) 'an error has been detected'
                  goto 9999
               endif
C              Second order formula (derivative evaluated in gp)
               dTudTb = (T_agpp - T_ag) / (T_bgpp - T_bg)
            else
C
C              CJDT
C
               indreg = 3
C
               if (logdeb) then
                  write(*,*) 'CJDT'
                  write(*,*) 'Parameter T_cg = ', T_cg
                  write(*,*) 'New values for T_ag = ', T_ag
               endif
C              T_cg &lt;= T_bg
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_bg, T_bg, u_bg, .true.,
     &              T_cg, P_cg, r_cg, u_cg, d_cg)
C
C              We take T_cgp > T_cg.
C              Note that as T_c changes, T_a and T_b do not
C
               T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
               call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &              P_bg, T_bg, u_bg, .true.,
     &              T_cgp, P_cgp, r_cgp, u_cgp, d_cgp)
C
            endif
         endif
C
         dpduc = (P_cgp - P_cg) / (u_cgp - u_cg)
C         write(*,*) dpduc
C
C
C******* Left state
C
         call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
     &        P_l, T_l, u_l, .false.,
     &        T_dg, P_dg, r_dg, u_dg, d_dg)
         T_dgp =  max((T_dg + Dtmin), (T_dg * coefder))
         call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
     &        P_l, T_l, u_l, .false.,
     &        T_dgp, P_dgp, r_dgp, u_dgp, d_dgp)
         if (dpduc .lt. 0.0d0) then
            write(*,*) ' dpduc &lt; 0 , ', dpduc
            logan = .true.
            write(*,*) 'subroutine racre2.f'
            write(*,*) 'an error has been detected'
            goto 9999
         endif
C
         dpdud = (P_dgp - P_dg) / (u_dgp - u_dg)
         if (dpdud .gt. 0.0d0) then
            write(*,*) 'dpdud > 0 , ', dpdud
            logan = .true.
            write(*,*) 'subroutine racre2.f'
            write(*,*) 'an error has been detected'
            goto 9999
         endif
C
C******* Intersection
C
         a = dpduc - dpdud
         b = ((dpduc * u_cg) - (dpdud * u_dg)) + (p_dg - p_cg)
         u_d =  b / a
         p_d = p_dg + (dpdud * (u_d - u_dg))
         p_c = p_cg + (dpduc * (u_d - u_cg))
         if ( abs (p_d - p_c) .gt. (eps * pref)) then
            write(*,*) 'pc - pd too big, ', (p_c - p_d)
            logan = .true.
            write(*,*) 'subroutine racre2.f'
            write(*,*) 'an error has been detected'
            goto 9999
         endif
C        p_c, p_d should be positive
         p_c = max (p_c, (eps * p_cg))
         p_d = max (p_d, (eps * p_dg))
C
C******* Computation of temperature T_d
C
         T_d = T_dg + ((p_d - p_dg) * (T_dgp - T_dg) / (p_dgp - p_dg))
         T_d = max (T_d, (eps * T_dg))
C
         T_c = T_cg + ((p_c - p_cg) * (T_cgp - T_cg) / (p_cgp - p_cg))
         T_c = max (T_c, (eps * T_cg))
C
         if ((indreg .eq. 3) .or. (indreg .eq. 4)) then
            if( indreg .eq. 3) then
C              CJDT
               T_a = T_ag
               T_b = T_bg
            else
               T_b = T_c
               T_a = T_ag
            endif
         else
C
C           indreg = 1 (WDF) or  (indreg = 2) CJDF
C
            if (T_c .le. T_bcj) then
               T_a = T_acj
               T_b = T_bcj
            else
               T_b = T_c
               b = dTudTb
               b = b / (2.0D0 * (T_bgp - T_bcj))
               a = -2.0D0 * b * (T_bcj - T_bg)
               dt1 = ((a + (b * (T_b - T_bg))) *
     &              (T_b - T_bg))
               dt2 = dTudTb * (T_b - T_bg)
               if (abs(dt1) .lt. abs(dt2)) then
                  T_a = T_ag + dt1
               else
                  T_a = T_ag + dt2
               endif
               T_a = max (T_a, (eps * T_a))
            endif
         endif
C
C
C******* TEST
C
         erro = abs (T_a - T_ag) + abs (T_b - T_bg) +
     &        abs(T_c - T_cg) + abs (T_d - T_dg) +
     &        abs (T_a - T_ag1) + abs (T_b - T_bg1)
         erro = erro / Tref
         if (logdeb) then
            write(*,*) 'Intersection'
            write(*,*) 'T_ag1, T_bg1 = ', T_ag1, T_bg1
            write(*,*) 'dpduc, dpdud = ', dpduc, dpdud
            write(*,*) 'p_c, p_d = ', p_c, p_d
            write(*,*) 'T_d, T_c = ', T_d, T_c
            write(*,*) 'T_a, T_b = ', T_a, T_b
            WRITE(*,*) 'erro =', erro
         endif
         if (erro .Le. 1.0D-14) then
            lognc = .false.
            goto 9998
         endif
      enddo
      if (logdeb ) then
         write(*,*) 'Warning. Convergence non-reached'
         write(*,*) 'erro = ', erro
      endif
      lognc = .true.
 9998 continue
      if (indreg .eq.  1) then
         call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
     &        P_r, T_r, u_r, .true.,
     &        T_a, P_a, r_a, u_a, d_a)
         call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
     &              Rgas_b, acv_b, q, k0, T_a, T_b,
     &              r_a, p_a, u_a, T_a,
     &              r_b, p_b, u_b, T_b, d_b, logan)
         d_b = min(d_a,d_b)
         r_c = r_b
         p_c = p_b
         u_c = u_b
         T_c = T_b
         d_c = d_b
      endif
      if (indreg .eq. 4) then
         call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
     &        acv_b, Rgas_b, q, T_a,
     &        P_r, T_r, u_r,
     &        T_c, P_c, r_c, u_c, d_c, k0sdt,
     &        T_a, P_a, r_a, u_a,
     &        logan)
         d_a = d_c
         r_b = r_c
         p_b = p_c
         u_b = u_c
         T_b = T_c
         d_b = d_c
         if (logdeb) then
            write(*,*) 'SDT'
            write(*,*) 'T_a, P_a, r_a, u_a'
            write(*,*) T_a, P_a, r_a, u_a
            write(*,*) 'T_c, P_c, r_c, u_c'
            write(*,*) T_c, P_c, r_c, u_c
            write(*,*) 'd_c, k0sdt'
            write(*,*) d_c, k0sdt
         endif
      endif
      if ((indreg .eq. 2) .or. (indreg .eq. 3)) then
         call funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
     &        k0, q,
     &        T_a, T_b,
     &        p_r, T_r, u_r,
     &        logdef,
     &        r_a, T_a, p_a, u_a, d_a,
     &        r_b, T_b, p_b, u_b, d_b,
     &        k0cjdt,
     &        logan)
         if (logan) goto 9999
         call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
     &        P_b, T_b, u_b, .true.,
     &        T_c, P_c, r_c, u_c, d_c)
         d_b = min(d_a, d_b)
         d_c = min(d_b, d_c)
      endif
      call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
     &     P_l, T_l, u_l, .false.,
     &     T_d, P_d, r_d, u_d, d_d)
      d_d = min(d_c, d_d)
      d_l = gam_l * P_l / r_l
      d_l = u_l - sqrt(d_l)
      d_l = min(d_l, d_d)
      d_r = gam_r * P_r / r_r
      d_r = u_r + sqrt(d_r)
      d_r = max(d_r, d_a)
C
      if (logdeb) then
         write(*,*) 'Final result'
         write(*,*) 'r'
         write(*,*) r_r, p_r, u_r, T_r, d_r
         write(*,*) 'a'
         write(*,*) r_a, p_a, u_a, T_a, d_a
         write(*,*) 'b'
         write(*,*) r_b, p_b, u_b, T_b, d_b
         write(*,*) 'c'
         write(*,*) r_c, p_c, u_c, T_c, d_c
         write(*,*) 'd'
         write(*,*) r_d, p_d, u_d, T_d, d_d
         write(*,*) 'l'
         write(*,*) r_l, p_l, u_l, T_l, d_l
C         write(*,*) 'acv'
C         do iter = 0, nordpo, 1
C            write(*,*) acv_b(iter), acv_l(iter)
C         enddo
      endif
 9999 continue
      return
      end