* fichier : comb.dgibi
*************************************************************************
*
* NOM         : comb.dgibi
*
* DESCRIPTION : We   compute   the  AICC  (Adiabatic  Isochoric  Complete
*               Combustion), the  AIBCC,  (Adiabatic  IsoBaric   Complete 
*               Combustion), the  CJ    (Chapman - Jouguet)  and  VN (von 
*               Neumann) states of a ZND detonation  for several mixtures 
*               at NIST-normal conditions involving H2 and air.
*               We check that some physical properties are satisfied.
*               Comparison  of  the  results  obtained  using  the  DETO 
*               operator  and the  ones  given by the CREBCOM  combustion 
*               model (operators PRIM and FLAM).
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               beccantini@cea.fr
* DATE        : 28/09/2006
*************************************************************************
*

   'OPTION' 'ECHO' 1 'DIME' 2
      'TRAC' 'X' 'ELEM' 'QUA4' ;

****************************************************************
*****      1D mesh                                         *****
****************************************************************
*
*
*  A4 ---- A3
*   |      | 
*   |  COM | 
*   |      |
*  A1 ---- A2
*

 A1 = 0.0  0.0 ;
 A2 = 1.0  0.0 ;

 A1A2 = A1 'DROIT' 1 A2 ;
 DOM1 = A1A2 'TRANSLATION' 1 (0.0 1.0) ;

*
**** The table 'DOMAINE'
*

 $DOM1   = 'MODELISER' DOM1 'EULER' ;
 TDOM1   = 'DOMA' $DOM1   'VF' ;
 QDOM1   = TDOM1    . 'QUAF' ;
  

*
***** Initial conditions and gas properties.
*

*
**** 0 from a numerical point of view
*

 zero = 1.0d-9 ;
 
*

 PINI =  'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.013D5
         'NATU' 'DISCRET' ;
 TINI =  'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 298.15
         'NATU' 'DISCRET' ;
 WINI =  'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 2 'UX' 0.0 'UY' 0.0
         'NATU' 'DISCRET' ;
  
* We suppose that the combustion is complete everywhere
 
 CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0 ;

*
***** Gas properties
*

 PGAZ = 'TABLE' ;

* Polynomial degree of specific heats
 
 PGAZ . 'NORD' = 4 ;

* Species explicitly treated in the Euler Equations

 PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;

* Species non explicitly treated

 PGAZ . 'ESPNEULE' = 'N2  ';

* Single gases properties

 PGAZ .  'H2  ' = 'TABLE'  ;
 PGAZ .  'H2O ' = 'TABLE'  ;
 PGAZ .  'N2  ' = 'TABLE'  ;
 PGAZ .  'O2  ' = 'TABLE'  ;

* R (J/Kg/K)

 mH2 = 2. '*' 1.00797E-3 ;
 mo2 = 2. '*' 15.9994E-3 ;
 mH2O = mh2 '+' (0.5 '*' mo2) ;
 mN2 = 2 '*' 14.0067E-3 ;
 RGAS = 8.31441 ;
 
 PGAZ .  'H2  ' . 'R' = RGAS '/' mh2 ;
 PGAZ .  'H2O ' . 'R' = RGAS '/' mh2o ;
 PGAZ .  'N2  ' . 'R' = RGAS '/' mn2 ;
 PGAZ .  'O2  ' . 'R' = RGAS '/' mo2 ;

* Polynomials regressions coefficients

 PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
 PGAZ .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
                              -1.82753232E-08 2.44485692E-12 ;
 PGAZ .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
                              8.78233606E-09 -3.05514485E-13 ;
 PGAZ .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
                              2.33636971E-08 -1.53304905E-12;

* Formation enthalpies at 0K (J/Kg)
*
*     h_i(0K)  = h_i(T0) '-' \int_0^{T0} cp_i(x) dx 
*              = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
*     

 PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
 PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
 PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;

* Passive scalars names

 PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0  ';

*
**** Initial conditions on molar fractions of hydrogen.
*

 LISTXH2 = PROG 0.08 0.134 0.188 0.242 0.29581 0.350 ;
 LISXH2F = PROG  ;
 LISXO2F = PROG  ;
 LISXWAF = PROG  ;
 LISXN2F = PROG  ;
 LISYH2F = PROG  ;
 LISYO2F = PROG  ;
 LISYWAF = PROG  ;
 LISYN2F = PROG  ;
* Inside
 LRINSI = PROG ;
 LTINSI = PROG ;
 LPINSI = PROG ;
 LCINSI = PROG ;
 LRGINSI = PROG ;
 LGINSI = PROG ;
 LGAINSI = PROG ;
* VN conditions
 LRVN = PROG ;
 LTVN = PROG ;
 LPVN = PROG ;
 LCVN = PROG ;
 LRGVN = PROG ;
 LGVN = PROG ;
 LGAVN = PROG ;
* AICC
 LRAICC = PROG ;
 LTAICC = PROG ;
 LPAICC = PROG ;
 LCAICC = PROG ;
 LRGAICC = PROG ;
 LGAICC = PROG ;
 LGAAICC = PROG ;
* AIBCC
 LRAIBC = PROG ;
 LTAIBC = PROG ;
 LPAIBC = PROG ;
 LCAIBC = PROG ;
 LRGAIBC = PROG ;
 LGAIBC = PROG ;
 LGAAIBC = PROG ;
* CJ conditions
 LRCJ = PROG ;
 LTCJ = PROG ;
 LPCJ = PROG ;
 LCCJ = PROG ;
 LRGCJ = PROG ;
 LGCJ = PROG ;
 LGACJ = PROG ;
* Released energy per mass unit
 LQ = PROG ;
*
 

 REPE BLCALC (DIME LISTXH2) ;
   
************************************************************************
************************************************************************
*****  LOOP to determine the initial  conditions, AICC, AIBCC, CJ  *****
*****  anv VN states                                               *****
************************************************************************
************************************************************************

*
**** Molar fractions (before combustion)
*

 cell = (EXTR LISTXH2 &BLCALC) ;
 cell1 = (1. - cell) / (1. + (0.79/0.21)) ;
 XH21 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' cell
        'NATU' 'DISCRET' ;
 XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' zero
         'NATURE' 'DISCRET' ;
 XO21 =  'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 
         'O2' cell1  'NATU' 'DISCRET' ;
 XN = XH21 'ET' XH2O1 'ET' XO21 ;
 CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0D0) ;
 XN21 =  CHUN '-' ('PSCAL' XN CHUN  ('MOTS' 'H2' 'O2' 'H2O')
         ('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
 XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
 XN = XN 'ET' XN21 ;

*  
**** We compute the mass fractions
*

 MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' Mh2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' Mo2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' Mh2o
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'N2' Mn2
         'NATURE' 'DISCRET') ;

 MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
    
 YH21 = (XH21 '*' MH2) '/' MTOT ;
 YO21 = (XO21 '*' MO2) '/' MTOT ;
 YN21 = (XN21 '*' MN2) '/' MTOT ;
 YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;

*
**** Burnt gas
*    Here X refers to the mole producted per 1 mole of unburned gas
*    N.B. The number of moles changes after the combustion!
*

* H2MAG = 1 in the region where H2 is the majority species.
* Conversely H2MAG = 0

 H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
    ('EXCO' 'O2'  XN) ;
 O2MAG = CHUN '-' H2MAG ;

 XO22a = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' zero ;
 XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
 
 XO22a = XO22a ;
 XH22a = XH22a ;
 
 XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' zero ;
 XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (xh21 '-' xh22b))) ;

 XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
 XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;

 XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
 XN22 = XN21 ;
*
**** mtot = weight of a mole  of unburned gas
*    mtot2 = weight of a mole of unburned gas after the combustion
*    Of course mtot2 = mtot
*

 MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
   ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;

 mtotcel = 'MAXIMUM' MTOT  ;
 
 
 'SI' ( ('MAXIMUM'  (mtot2 '-' mtot) 'ABS') '>' (zero '*' mtotcel)) ;
    'MESSAGE' 'Problem in the computation of the burnt gas (1).' ;
    'ERREUR' 21 ;
 'FINSI' ;

*
**** xtot2 = mole of burnt gas per mole of unburned gas
*
 
 xtot2 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
   CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');

*
**** X = real molar fractions
*
 
 xh22 = xh22 '/' xtot2 ;
 xo22 = xo22 '/' xtot2 ;
 xh2o2 = xh2o2 '/' xtot2 ;
 xn22 = xn22 '/' xtot2 ;

 LISXH2F = LISXH2F 'ET' (PROG (MAXI XH22)) ;
 LISXO2F = LISXO2F 'ET' (PROG (MAXI XO22)) ;
 LISXWAF = LISXWAF 'ET' (PROG (MAXI XH2O2)) ;
 LISXN2F = LISXN2F 'ET' (PROG (MAXI XN22)) ;


 xtot3 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
   CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');

 'SI' ( ('MAXIMUM' (xtot3 '-' 1.0D0) 'ABS' ) '>'
    zero) ;
    'MESSAGE' 'Problem in the computation of the burnt gas (2).' ;
    'ERREUR' 21 ;
 'FINSI' ;

*
**** mtot2 = weight of 1 mole of burned gas
*

 MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
   ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;

*
**** Mass fractions of burned gas after the complete combustion
*
    
 YH22 = (XH22 '*' MH2) '/' MTOT2 ;
 YO22 = (XO22 '*' MO2) '/' MTOT2 ;
 YN22 = (XN22 '*' MN2) '/' MTOT2 ;
 YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;

 LISYH2F = LISYH2F 'ET' (PROG (MAXI YH22)) ;
 LISYO2F = LISYO2F 'ET' (PROG (MAXI YO22)) ;
 LISYWAF = LISYWAF 'ET' (PROG (MAXI YH2O2)) ;
 LISYN2F = LISYN2F 'ET' (PROG (MAXI YN22)) ;


*
**** The conservative variables
*

* Gas constant

 RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
        (('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
        (('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
        (('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
   ;
   
* Density

 RNINI = PINI '/' (RINI '*' TINI) ;
 RN = 'COPIER' RNINI ;

* Internal energy

 eini = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
 Tfois = Tini ;
 
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
 
    acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
            &BL1))  ;
            
    eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
    
    Tfois = Tfois * tini ;
    
 'FIN' BL1 ;

 RETNINI = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
   ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
 RETN = RETNINI ;
 GN = RN '*' WINI ;
 YN = YH21 '+' YO21 '+' YH2O1 ;
 RYN = RN '*' YN ;
*
**** K0/DX = chemical source term
*
*    i.e. dcsi/dt  = K0 '/' DX . {criterium function}
*         where
*         {criterium function} = 1.0 or 0.0
*
*         d(rho etot)/dt = -rho (dcsi/dt) (h_f '-' h_i)
*
*         h_f reference enthalpy (0K) of the burned gas
*         h_i reference enthalpy (0K) of the unburned gas
*
**** Mass fraction of H2 before and after the combustion
*
*    CSIMAX = 1.0 -> complete combustion
*

 YININ = yh21 ;
 YFINN = YININ '+' ((yh22 '-' yh21) '*' csimax) ;
 K0N = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 500 ;

*
**** Passive scalars  
*

 RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
               
* Test of prim
 VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
 PCEL = 'MAXIMUM' PINI ;
 TCEL = 'MAXIMUM' TINI ;
 'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
    'MESSAGE' 'Problem in the computation of the energy  (1).' 
    'ERREUR' 21 ;
 'FINSI' ;
 'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
    'MESSAGE' 'Problem in the computation of the energy  (2).' 
    'ERREUR' 21 ;
 'FINSI' ;

 PINSI1 = MAXI P ;
 TINSI1 = MAXI T ;
 RINSI1 = MAXI RN ; 
 GAMINSI = MAXI GAMN ;
 GAMAINSI = 1 + (MAXI (P / RETN)) ;  
 RGAS = MAXI (P / (RN * T)) ; 

 LRINSI = LRINSI 'ET' (PROG RINSI1) ;
 LTINSI = LTINSI 'ET' (PROG TINSI1) ;
 LPINSI = LPINSI 'ET' (PROG PINSI1) ;
 LRGINSI = LRGINSI 'ET' (PROG RGAS) ;
 LGINSI = LGINSI 'ET' (PROG GAMINSI) ;
 LGAINSI = LGAINSI 'ET' (PROG GAMAINSI) ;
 LCINSI = LCINSI 'ET' (PROG ((GAMINSI * PINSI1 / RINSI1) ** 0.5)) ;

*
* VN state
*
* VN state via DETO operator (not really complete combustion...)
*
 CHPO1 = 'MANUEL' 'CHPO'  ('DOMA' $DOM1 'CENTRE') 6
         'H2' (MAXI  XH21) 
         'O2' (MAXI XO21) 
         'N2' (1. - ((MAXI  XH21) + (MAXI XO21)))
         'H2O' 0.0
         'P' (MAXI PINI) 'T' (MAXI TINI) ;
         
 CHPCJ CHPVN CHPAICC = 'DETO' CHPO1 ;

 PVN =  (EXCO CHPVN 'PZND') ;
 RVN =  (EXCO CHPVN 'RZND') ;
 TVN =  (EXCO CHPVN 'TZND') ;
 DVN =  (EXCO CHPCJ 'VCJ') ;
* VVN via RH conditions
 VVN = DVN '*' (RNINI / RVN) ;
 VVN = DVN - VVN ;
 VVN = (NOMC VVN 'UX') + 
   ('MANUEL' 'CHPO'  ('DOMA' $DOM1 'CENTRE') 1 'UY' 0.0) ;

* Internal energy

 EVN = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
 Tfois = TVN ;
 
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
 
    acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
            &BL1))  ;
            
    eVN = eVN '+' (acel '*' (Tfois '/' &BL1)) ;
    
    Tfois = Tfois * tVN ;
    
 'FIN' BL1 ;

 RETNVN = RVN '*' (eVN '+' (0.5 '*' ('PSCAL' VVN VVN
   ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;

 V1VN P1VN T1VN YNVN SNVN GAMVN = 'PRIM' 'PERFTEMP' PGAZ RVN (RVN * VVN)
      RETNVN (RVN * YN)  RSN ;

 ERRO = MAXI ((TVN - T1VN) / T1VN) ABS ;

 'SI' (ERRO > 1.0d-3) ;
     MESS 'Problem in VN' ;
     ERREUR 21 ;
 FINS ;
*
* We check that the RH on the energy is satisfied
*
 VVN = DVN '*' (RNINI / RVN) ;
 HVN = (EVN + (PVN / RVN)) + (0.5 * VVN * VVN) ;
 HINI = (EINI + (PINI / RNINI)) + (0.5 * DVN * DVN) ;
 ERRO = MAXI ((HVN - HINI) / HVN) ABS ;
 'SI' (ERRO > 1.0d-2) ;
    MESS 'Problem in VN (2)' ;
    ERREUR 21 ;
 FINS ;

 PVN1 = MAXI PVN ;
 TVN1 = MAXI TVN ;
 RVN1 = MAXI RVN ; 
 GAMVN = MAXI GAMVN ;
 GAMAVN = 1 + (MAXI (PVN / (RVN * EVN))) ;  
 RGAS = MAXI (PVN / (RVN * TVN)) ; 

 LRVN = LRVN 'ET' (PROG RVN1) ;
 LTVN = LTVN 'ET' (PROG TVN1) ;
 LPVN = LPVN 'ET' (PROG PVN1) ;
 LRGVN = LRGVN 'ET' (PROG RGAS) ;
 LGVN = LGVN 'ET' (PROG GAMVN) ;
 LGAVN = LGAVN 'ET' (PROG GAMAVN) ;
 LCVN = LCVN 'ET' (PROG ((GAMVN * PVN1 / RVN1) ** 0.5)) ;

*
***** Combustion via the operator FLAM
*
        
 DELTAXN = ('DOMA' $DOM1 'VOLUME') '**' (1. '/' 2.)  ;
 YN = RYN '/' RN ;
 
*
**** We want the complete combustion in DOM1.
*    We can use the operator 'FLAM', with epsilon=0 and
*    DELTAT >> DELTATC
*    where DELTATC is the characteristic time step linked to
*    the source term.
*

 SN = RSN '/' RN ;
 K0N = 'EXCO' 'K0' SN ;
 YININ = 'EXCO' 'H2IN' SN 'H2' ;
 YFINN = 'EXCO' 'H2FI' SN 'H2' ;

 DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
 LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
 LCOEF = 'PROG' 1.0 0.5 -1.0 ;

 DELTARE DELTARY = 'FLAM' 'CREBCOM2'  $DOM1 PGAZ LMOT1 LCOEF
   RN YN  YININ  YFINN  K0N  DELTAXN 
   0.0 (1D4 '*' DELTATC)  0.3 ;

 Q = ((MAXI (DELTARE / RN))) ;
 Q1 = (('NOMC' 'SCAL' (yh21 '-' yh22)) * (PGAZ .'H2  ' . 'H0K')) + 
      (('NOMC' 'SCAL' (yo21 '-' yo22)) * (PGAZ .'O2  ' . 'H0K')) + 
      (('NOMC' 'SCAL' (yn21 '-' yn22)) * (PGAZ .'N2  ' . 'H0K')) + 
      (('NOMC' 'SCAL' (yh2o1 '-' yh2o2)) * (PGAZ .'H2O ' . 'H0K')) ;

 ERRO = MAXI ((Q - Q1) / Q1) ABS ;
 'SI' (ERRO > 1.0d-2) ;
    MESS 'Problem in Q' ;
    ERREUR 21 ;
 FINS ;

 LQ = LQ 'ET' (PROG Q) ;

 RYN = RYN '+' DELTARY ;
 YN = RYN '/' RN ;
 RETN = RETN '+' DELTARE ;
*
*
 
*
 VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
 PAICC1 = MAXI PN ;
 TAICC1 = MAXI TN ;
 RAICC1 = MAXI RN ; 
 GAMAICC = MAXI GAMMAN ;
 GAMAAICC = 1 + (MAXI (PN / RETN)) ;  
 RGAS = MAXI (PN / (RN * TN)) ; 
*
* AICC, CJ via DETO operator (not really complete combustion...)
*
 CHPO1 = 'MANUEL' 'CHPO'  ('DOMA' $DOM1 'CENTRE') 6
         'H2' (MAXI  XH21) 
         'O2' (MAXI XO21) 
         'N2' (1. - ((MAXI  XH21) + (MAXI XO21)))
         'H2O' 0.0
         'P' (MAXI PINI) 'T' (MAXI TINI) ;
         
 CHPCJ CHPZND CHPAICC = 'DETO' CHPO1 ;
 PAICC2 = MAXI (EXCO CHPAICC 'PAIC') ;
 RAICC2 = MAXI (EXCO CHPAICC 'RAIC') ;
 TAICC2 = MAXI (EXCO CHPAICC 'TAIC') ;
*
 LERRO = ABS (PROG ((PAICC1 - PAICC2) / PAICC1) 
                   ((RAICC1 - RAICC2) / RAICC1) 
                   ((TAICC1 - TAICC2) / TAICC1)) ;
 ERRO = MAXI LERRO ;
 'SI' (ERRO > 1.0d-2) ;
    MESS 'Problem in AICC' ;
    ERREUR 21 ;
 FINS ;
 LRAICC = LRAICC 'ET' (PROG RAICC1) ;
 LTAICC = LTAICC 'ET' (PROG TAICC1) ;
 LPAICC = LPAICC 'ET' (PROG PAICC1) ;
 LRGAICC = LRGAICC 'ET' (PROG RGAS) ;
 LGAICC = LGAICC 'ET' (PROG GAMAICC) ;
 LGAAICC = LGAAICC 'ET' (PROG GAMAAICC) ;
 LCAICC = LCAICC 'ET' (PROG ((GAMAICC * PAICC1 / RAICC1) ** 0.5)) ;
*
* CJ
*
 PCJ =  (EXCO CHPCJ 'PCJ') ;
 RCJ =  (EXCO CHPCJ 'RCJ') ;
 TCJ =  (EXCO CHPCJ 'TCJ') ;
 DCJ =  (EXCO CHPCJ 'VCJ') ;
 
 Tfois = TCJ ;
 eCJ = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;

 
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
 
    acel = (('NOMC' 'SCAL' yh22) '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yo22) '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yh2o2) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yn22) '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
            &BL1))  ;
            
    eCJ = eCJ '+' (acel '*' (Tfois '/' &BL1)) ;
    
    Tfois = Tfois * TCJ ;
    
 'FIN' BL1 ;

 RETNCJ = RCJ '*' eCJ  ;

 RWCJ = RCJ * VN * 0 ;

 V1CJ P1CJ T1CJ YNCJ SNCJ GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RCJ RWCJ 
     RETNCJ (RCJ * YN)  (RCJ * SN) ;
 CCJ = (GAMMAN * (P1CJ / RCJ)) '**' 0.5 ;
 VCJ = DCJ * (RNINI / RCJ) ;
 erro = MAXI ((VCJ / CCJ) - 1.0) ABS ; 

* We check weather the speed in the shock reference is equal to the speed of
* sound.
 'SI' (erro > 2.0d-2) ;
    MESS 'Problem in CJ' ;
    ERREUR 21 ;
 FINS ;
 GAMACJ = 1 + (MAXI (PCJ / RETNCJ)) ;  

 LRCJ = LRCJ 'ET' (PROG (maxi RCJ)) ;
 LTCJ = LTCJ 'ET' (PROG (maxi TCJ)) ;
 LPCJ = LPCJ 'ET' (PROG (maxi PCJ)) ;
 LRGCJ = LRGCJ 'ET' (PROG RGAS) ;
 LGCJ = LGCJ 'ET' (PROG (MAXI GAMMAN)) ;
 LGACJ = LGACJ 'ET' (PROG GAMACJ) ;
 LCCJ = LCCJ 'ET' (PROG (maxi CCJ)) ;
 
*
* AIBCC
*
* Gas constant
*
 YO2 = 'MAXIMUM' ('EXCO' 'O2' YN) ;
 YH2 = 'MAXIMUM' ('EXCO' 'H2' YN) ;
 YH2O = 'MAXIMUM' ('EXCO' 'H2O' YN) ;
 YN2 = (-1 '*' ((YO2 + YH2 + YH2O) '-' 1)) ;
 R =  ((PGAZ .  'H2O ' . 'R') '*' YH2O) '+' 
      ((PGAZ .  'O2  ' . 'R') '*' YO2) '+' 
      ((PGAZ .  'H2  ' . 'R') '*' YH2) '+' 
      ((PGAZ .  'N2  ' . 'R') '*' YN2) ;
 T0 = 500 ;
 T1 = 3500 ;
 hfin = 'MAXIMUM' ((RETN '/' RN) '+' (PINI '/' RN)) ;
* hfin = EINI '+' q '+'  (PINI '/' RN) ;
*
* Regula falsi
*
 'REPETER' BLITER 100 ;
    T = (T0 '+' T1) '/' 2. ;
    eini =  0.0 ;
    Tfois = T ;
    'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
        acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
           &BL1)) '+'
           (YN2 '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
           &BL1))  '+'
           (YO2 '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
           &BL1))  '+'
           (YH2 '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
           &BL1))  ;
        eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
        Tfois = Tfois * T ;
    'FIN' BL1 ;
     hini = eini '+' (R * T) ;
     'MESS' ('CHAINE' 'h, T ' hini ' ' T) ;
     'SI' (< ('ABS' (hini '-' hfin)) 1.0D-4) ;
         'MESSAGE' OK ;
         'QUITTER' BLITER ;
     'SINON' ;
         'SI' (hini > hfin) ;
             T1 = T ;
         'SINON' ;
             T0 = T ;
         'FINSI' ;
     'FINSI' ;
 'FIN' BLITER ;
*
 TAIBC1=T;
 PAIBC1=(MAXI PINI) ;
 RNAIBCC = PINI '/' (R '*' T) ;
 RAIBC1=(MAXI RNAIBCC) ;
 eini = 0.0 ;
 Tfois = T ;
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
    acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
           &BL1)) '+'
           (YN2 '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
           &BL1)) '+'
           (YO2 '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
           &BL1)) '+'
           (YH2 '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
           &BL1))  ;
    eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
    Tfois = Tfois * T ;
 'FIN' BL1 ;
 RETAIBCC = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' eini) '*'
           RNAIBCC  ;
 RYNAIBCC = RNAIBCC '*' YN ;
 VN2 P2 T2 Y2 S2 GAMN = 'PRIM' 'PERFTEMP' PGAZ RNAIBCC GN RETAIBCC
    RYNAIBCC  RSN ;
 TAIBC2= MAXI T2 ;
 PAIBC2= MAXI P2 ;
 LERRO = ABS (PROG ((PAIBC1 - PAIBC2) / PAIBC1) 
                   ((TAIBC1 - TAIBC2) / TAIBC1)) ;
 ERRO = MAXI LERRO ;
 'SI' (ERRO > 1.0d-2) ;
    MESS 'Problem in AIBCC' ;
    ERREUR 21 ;
 FINS ;
 LRAIBC = LRAIBC 'ET' (PROG RAIBC1) ;
 LTAIBC = LTAIBC 'ET' (PROG TAIBC1) ;
 LPAIBC = LPAIBC 'ET' (PROG PAIBC1) ;
 GAMAIBC = MAXI GAMN ;
 GAMAAIBC = 1 + (MAXI (P2 / RETAIBCC)) ;  
 RGAS = MAXI (P2 / (RNAIBCC * T2)) ; 
 LRGAIBC = LRGAIBC 'ET' (PROG RGAS) ;
 LGAIBC = LGAIBC 'ET' (PROG GAMAIBC) ;
 LGAAIBC = LGAAIBC 'ET' (PROG GAMAAIBC) ;
 LCAIBC = LCAIBC 'ET' (PROG ((GAMAIBC * PAIBC1 / RAIBC1) ** 0.5)) ;
*
*************************************************************************
*************************************************************************
*********************** END OF THE LOOP *********************************
*************************************************************************
*************************************************************************
*
 FIN BLCALC ;

*
**** Some LaTex tables.
*
 SI FAUX ;
 OPTI ECHO 0 ;
 REPE BLCALC (DIME LISTXH2) ;
*
    R = EXTR LRINSI &BLCALC ;
    P = EXTR LPINSI &BLCALC ;
    T = EXTR LTINSI &BLCALC ;
    C = EXTR LCINSI &BLCALC ;
    GAM = EXTR LGINSI &BLCALC ;
    GAMA = EXTR LGAINSI &BLCALC ;
    RGAS = EXTR LRGINSI &BLCALC ;
    AA = CHAI &blcalc ' &  Inside & ' 'FORMAT' '(E9.3)'  R  ' & ' 
      'FORMAT' '(E9.3)'  P ' & ' 'FORMAT' '(E9.3)'  T ' & ' 
      'FORMAT' '(E9.3)'  C ' & ' 'FORMAT' '(E9.3)'  RGAS  ' & ' 
      'FORMAT' '(E9.3)'  GAM ' & ' 'FORMAT' '(E9.3)'  GAMA ' \\';
    MESS AA ;
*
    R = EXTR LRAICC &BLCALC ;
    P = EXTR LPAICC &BLCALC ;
    T = EXTR LTAICC &BLCALC ;
    C = EXTR LCAICC &BLCALC ;
    GAM = EXTR LGAICC &BLCALC ;
    GAMA = EXTR LGAAICC &BLCALC ;
    RGAS = EXTR LRGAICC &BLCALC ;
    AA = CHAI &blcalc ' &  AICC & ' 'FORMAT' '(E9.3)'  R  ' & ' 
      'FORMAT' '(E9.3)'  P ' & ' 'FORMAT' '(E9.3)'  T ' & ' 
      'FORMAT' '(E9.3)'  C ' & ' 'FORMAT' '(E9.3)'  RGAS  ' & ' 
      'FORMAT' '(E9.3)'  GAM ' & ' 'FORMAT' '(E9.3)'  GAMA ' \\';
    MESS AA ;
*
    R = EXTR LRAIBC &BLCALC ;
    P = EXTR LPAIBC &BLCALC ;
    T = EXTR LTAIBC &BLCALC ;
    C = EXTR LCAIBC &BLCALC ;
    GAM = EXTR LGAIBC &BLCALC ;
    GAMA = EXTR LGAAIBC &BLCALC ;
    RGAS = EXTR LRGAIBC &BLCALC ;
    AA = CHAI &blcalc ' &  AIBCC & ' 'FORMAT' '(E9.3)'  R  ' & ' 
      'FORMAT' '(E9.3)'  P ' & ' 'FORMAT' '(E9.3)'  T ' & ' 
      'FORMAT' '(E9.3)'  C ' & ' 'FORMAT' '(E9.3)'  RGAS  ' & ' 
      'FORMAT' '(E9.3)'  GAM ' & ' 'FORMAT' '(E9.3)'  GAMA ' \\';
    MESS AA ;
*
    R = EXTR LRCJ &BLCALC ;
    P = EXTR LPCJ &BLCALC ;
    T = EXTR LTCJ &BLCALC ;
    C = EXTR LCCJ &BLCALC ;
    GAM = EXTR LGCJ &BLCALC ;
    GAMA = EXTR LGACJ &BLCALC ;
    RGAS = EXTR LRGCJ &BLCALC ;
    AA = CHAI &blcalc ' &  CJ & ' 'FORMAT' '(E9.3)'  R  ' & ' 
      'FORMAT' '(E9.3)'  P ' & ' 'FORMAT' '(E9.3)'  T ' & ' 
      'FORMAT' '(E9.3)'  C ' & ' 'FORMAT' '(E9.3)'  RGAS  ' & ' 
      'FORMAT' '(E9.3)'  GAM ' & ' 'FORMAT' '(E9.3)'  GAMA ' \\';
    MESS AA ;
*
    R = EXTR LRVN &BLCALC ;
    P = EXTR LPVN &BLCALC ;
    T = EXTR LTVN &BLCALC ;
    C = EXTR LCVN &BLCALC ;
    GAM = EXTR LGVN &BLCALC ;
    GAMA = EXTR LGAVN &BLCALC ;
    RGAS = EXTR LRGVN &BLCALC ;
    AA = CHAI &blcalc ' &  VN & ' 'FORMAT' '(E9.3)'  R  ' & ' 
      'FORMAT' '(E9.3)'  P ' & ' 'FORMAT' '(E9.3)'  T ' & ' 
      'FORMAT' '(E9.3)'  C ' & ' 'FORMAT' '(E9.3)'  RGAS  ' & ' 
      'FORMAT' '(E9.3)'  GAM ' & ' 'FORMAT' '(E9.3)'  GAMA ' \\';
    MESS AA ;
 FIN BLCALC ;
 OPTI ECHO 1 ;
 OPTI ECHO 0 ;
 REPE BLCALC (DIME LISTXH2) ;
*
    XH2 = EXTR LISXH2F &BLCALC ;
    XO2 = EXTR LISXO2F &BLCALC ;
    XWA = EXTR LISXWAF &BLCALC ;
    XN2 = EXTR LISXN2F &BLCALC ;
    YH2 = EXTR LISYH2F &BLCALC ;
    YO2 = EXTR LISYO2F &BLCALC ;
    YWA = EXTR LISYWAF &BLCALC ;
    YN2 = EXTR LISYN2F &BLCALC ;
    Q   = EXTR LQ &BLCALC ;
*    AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)'  XH2  ' & ' 
*      'FORMAT' '(E9.3)'  XO2 ' & ' 'FORMAT' '(E9.3)'  XWA ' & ' 
*      'FORMAT' '(E9.3)'  XN2 ' & ' 'FORMAT' '(E9.3)'  YH2  ' & ' 
*      'FORMAT' '(E9.3)'  YO2 ' & ' 'FORMAT' '(E9.3)'  YWA ' & ' 
*      'FORMAT' '(E9.3)'  YN2  ' \\';
    AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)'  XH2  ' & ' 
      'FORMAT' '(E9.3)'  XO2 ' & ' 'FORMAT' '(E9.3)'  XWA ' & ' 
      'FORMAT' '(E9.3)'  YH2 ' & ' 
      'FORMAT' '(E9.3)'  YO2 ' & ' 'FORMAT' '(E9.3)'  YWA ' & ' 
      'FORMAT' '(E9.3)'  Q   ' \\';
    MESS AA ;
    MESS '\hline' ;
 FIN BLCALC ;
 OPTI ECHO 1 ;
 FINSI ;
 FIN ;