* fichier : cube_CJDF3D.dgibi
************************************************************************
* Section : Chimie Combustion
************************************************************************
****************************************************************
****    APPROCHE VF "Cell-Centred Formulation" pour la      ****
****    combustion.                                         ****
****    MODELE RDEM                                         ****
****                                                        ****
****    OPERATEURS 'PRIM', PRET, KONV                       ****
****                                                        ****
****    PROPAGATION D'UNE DEFLAGRATION  DANS UN CUBE        ****
****    We check some symmetry properties.                  ****
****                                                        ****
****    A. BECCANTINI, SFME/LTMF, 22.12.10                  ****
****                                                        ****
****************************************************************
*
* The file is divided into 6 parts
*
* 1) mesh
* 2) initial conditions and gas properties
* 3) parameters used in the computation
* 4) the main part (where the computation is realised)
* 5) the post-treatment
*
 'OPTION' 'ECHO' 1 'DIME' 3
*         'TRAC' 'PSC' ;
         'TRAC' 'X' ;
 GRAPH = VRAI ;
 GRAPH = FAUX ;
* Upwind scheme
* METO = 'VLH' ;
* METO = 'SS' ;
 METO = 'AUSMPUP' ;
* 
******************
**** PRIMCONS ****
******************
 'DEBPROC' PRIMCONS ;
 'ARGUMENT' PGAS*TABLE TN1*'CHPOINT ' TN2*'CHPOINT '
   PN1*'CHPOINT ' PN2*'CHPOINT ' VN1*'CHPOINT ' VN2*'CHPOINT ' ;
*
* ETHER = int_0^T cv(T') dT'           T <  TMAX
*       = int_0^TMAX cv(T') dT' '+'
*         cv(TMAX)                     T >= TMAX
 ESP1 = 'EXTRAIRE' (PGAS . 'SPECIES') 1 ;
* DY1 = y_i - y_f for the species 1 
 DY1 = (('EXTRAIRE' (PGAS . 'MASSFRA') 1) '-'
   ('EXTRAIRE' (PGAS . 'MASSFRA') 2)) ;
 COEF1 = ('EXTRAIRE' (PGAS . 'CHEMCOEF') 1) '*'
    (PGAS . ESP1 . 'W') ;
 YFINPH1 = 1.0 ;
 YFINPH2 = 1.0 ;
 'SI' (COEF1 > 0) ;
    YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
    YPH1 = YPH2 '+' DY1 ;
 'SINON' ;
    YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
    YPH2 = YPH1 '-' DY1 ;
 'FINSI' ;
 YFINPH1 = YFINPH1 '-' YPH1 ;
 YFINPH2 = YFINPH2 '-' YPH2 ;
 PRYPH1 = 'PROG' YPH1 ;
 PRYPH2 = 'PROG' YPH2 ;
 'REPETER' BLESP (('DIME' (PGAS . 'SPECIES')) '-' 2) ;
    ESP = 'EXTRAIRE' (PGAS . 'SPECIES') (&BLESP '+' 1)  ;
    COEF = ('EXTRAIRE' (PGAS . 'CHEMCOEF') (&BLESP '+' 1)) 
       '*' (PGAS . ESP . 'W') ;
    DY = (DY1 * (COEF '/' COEF1)) ;
    'SI' (COEF > 0) ;
       YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
       YPH1 = YPH2 '+' DY ;
    'SINON' ;
       YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
       YPH2 = YPH1 '-' DY ;
    'FINSI' ;
    PRYPH1 = PRYPH1 'ET' ('PROG' YPH1) ;
    PRYPH2 = PRYPH2 'ET' ('PROG' YPH2) ;
    YFINPH1 = YFINPH1 '-' YPH1 ;
    YFINPH2 = YFINPH2 '-' YPH2 ;
 'FIN' BLESP ;
 PRYPH1 = PRYPH1 'ET' ('PROG' YFINPH1) ; 
 PRYPH2 = PRYPH2 'ET' ('PROG' YFINPH2) ;
* 'LISTE' PRYPH1 ;
* 'LISTE' PRYPH2 ; 
*
 TMAX = (PGAS . 'TMAX') ;
* TCAL1 = MIN TN1, TMAX 
 TCAL1 = 0.5D0 '*' ((TMAX '+' TN1) '-' ('ABS' (TN1 '-' TMAX))) ;
 DTN1 = TN1 '-' TCAL1 ;
* TCAL1 = MIN TN1, TMAX 
 TCAL2 = 0.5D0 '*' ((TMAX '+' TN2) '-' ('ABS' (TN2 '-' TMAX))) ;
 DTN2 = TN2 '-' TCAL2 ;
*
* Internal energy (J/kg in SI)
*
 ETHER1 = 0.0 ;
 CV1 = 0.0 ;
 ETHER2 = 0.0 ;
 CV2 = 0.0 ;
 FUNTN1 = 1.0 ;
 FUNTN2 = 1.0 ;
 'REPETER' BLPO ((PGAS . 'NORD') '+' 1) ;
    'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
       ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
       YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
       YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
       AA = 'EXTRAIRE' (PGAS . ESP . 'A') &BLPO ;
       DCV1 = (AA * YCEL1 * FUNTN1) ;
       DCV2 = (AA * YCEL2 * FUNTN2) ;
       CV1 = CV1 '+' DCV1 ;
       CV2 = CV2 '+' DCV2 ;
       ETHER1 = ETHER1 '+' (DCV1 * TCAL1 '/' (&BLPO)) ;
       ETHER2 = ETHER2 '+' (DCV2 * TCAL2 '/' (&BLPO)) ;
    'FIN' BLESP ;
    FUNTN1 = FUNTN1 '*' TCAL1 ;
    FUNTN2 = FUNTN2 '*' TCAL2 ;
 'FIN' BLPO ;
 ETHER1 = ETHER1 '+' (CV1 '*' DTN1) ;
 ETHER2 = ETHER2 '+' (CV2 '*' DTN2) ;
*
* Formation energy/enthalpy (J/kg in SI) and gas constant (J/kg/K)
*
 EFORM1 = 0.0 ;
 EFORM2 = 0.0 ;
 RGAS1 = 0.0 ;
 RGAS2 = 0.0 ;
 'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
    ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
    YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
    YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
    EFORM1 = EFORM1 '+' (YCEL1 * (PGAS . ESP . 'H0K')) ;
    EFORM2 = EFORM2 '+' (YCEL2 * (PGAS . ESP . 'H0K')) ;
    RGAS1 = RGAS1 '+' (YCEL1 *  (PGAS . 'RUNIV') '/' 
      (PGAS . ESP . 'W')) ;
    RGAS2 = RGAS2 '+' (YCEL2 *  (PGAS . 'RUNIV') '/' 
      (PGAS . ESP . 'W')) ;
 'FIN' BLESP ;
*
* Computation of the conservative variables
*
 RN1 = PN1 '/' (RGAS1 '*' TN1) ;
 RN2 = PN2 '/' (RGAS2 '*' TN2) ;
 GN1 = RN1 * VN1 ;
 GN2 = RN2 * VN2 ;
 LVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
 ECIN1 = 0.5D0 '*' ('PSCAL' VN1 VN1 LVEL LVEL) ;
 ECIN2 = 0.5D0 '*' ('PSCAL' VN2 VN2 LVEL LVEL) ;
 RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ; 
 RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ;
*
 'RESPRO' RN1 RN2 GN1 GN2 RETN1 RETN2 ;
 'FINPROC' ;

* 
****************************************************************
*****      Cube                                            *****
****************************************************************
*
*
*   A4 ------ A3
*   |         | 
*   |         | 
*   | ZONE    | 
*   A1 ------ A2
*   |         
*   |   L1    | 
*  >|<------->|
*
*  AR = activation region close to A1
*  ZONE1
*  ZONE2
*

 RAF = 10 ;
 
 L1 = 1. ;
 L2 = 1. ;
 L3 = 1. ;

 DX = L1 '/' RAF ;
 N1 = RAF ;
 N2 = RAF ;
 N3 = RAF ;

 A1 = 0.0 0.0 0.0 ;
 A2 = L1 0.0 0.0 ;
 A3 = L1  L1 0.0 ;
 A4 = 0.0  L2 0.0 ;
* 
 'OPTION' 'ELEM' 'SEG2' ;
 A1A2 = A1 'DROIT' N1 A2 ;
 A2A3 = A2 'DROIT' N2 A3 ;
 A3A4 = A3 'DROIT' N1 A4 ;
 A4A1 = A4 'DROIT' N2 A1 ;
*
*
**** 2D mesh
*
 'OPTION' 'ELEM' 'QUA4' ;
 DOMP = 'DALLER' A1A2 A2A3 A3A4 A4A1 ;
* 
 'OPTION' 'ELEM' 'CUB8' ;
 DOMTOT = DOMP 'VOLUME' 'TRANSLATION' N3 (0.0 0.0 L3) ;
 DOM1 = DOMTOT 'ELEM' 'APPU' 'LARG' A1 ;
 DOM2 = 'DIFF' DOMTOT DOM1 ;
*
* DOMTOT = total region
*
 DOMLIM = 'ENVE' DOMTOT ;
*
**** The tables 'DOMAINE'
*
 $DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
 $DOMLIM = 'MODELISER' DOMLIM 'EULER' ;
 $DOM1   = 'MODELISER' DOM1 'EULER' ;
 $DOM2   = 'MODELISER' DOM2 'EULER' ;
* 
 TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
 TDOMLIM = 'DOMA' $DOMLIM 'VF' ;
 TDOM1   = 'DOMA' $DOM1   'VF' ;
 TDOM2   = 'DOMA' $DOM2   'VF' ;
* 
 QDOMTOT = TDOMTOT  . 'QUAF' ;
 QDOMLIM = TDOMLIM  . 'QUAF' ;
 QDOM1   = TDOM1    . 'QUAF' ;
 QDOM2   = TDOM2    . 'QUAF' ;
*
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMLIM ;
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM1 ;
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM2 ;
*
**** MOD1 = model (created just to display the CHAMELEMs)
*
 MOD1 =  'MODELISER'  ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
*
**** Lines for graphics
*
 DOMX = 'POIN' ('COORDONNEE' 2 DOMTOT) 'MINI' ;
 DOMX = 'POIN' ('COORDONNEE' 3 DOMX) 'MINI' ;
 DOMX = DOMTOT 'ELEM' 'APPU' 'LARG' DOMX ;
 DOMY = 'POIN' ('COORDONNEE' 3 DOMTOT) 'MINI' ;
 DOMY = 'POIN' ('COORDONNEE' 1 DOMY) 'MINI' ;
 DOMY = DOMTOT 'ELEM' 'APPU' 'LARG' DOMY ;
 DOMZ = 'POIN' ('COORDONNEE' 1 DOMTOT) 'MINI' ;
 DOMZ = 'POIN' ('COORDONNEE' 2 DOMZ) 'MINI' ;
 DOMZ = DOMTOT 'ELEM' 'APPU' 'LARG' DOMZ ;
 $DOMX   = 'MODELISER' DOMX 'EULER' ;
 $DOMY   = 'MODELISER' DOMY 'EULER' ;
 $DOMZ   = 'MODELISER' DOMZ 'EULER' ;
 TDOMX   = 'DOMA' $DOMX   'VF' ;
 TDOMY   = 'DOMA' $DOMY   'VF' ;
 TDOMZ   = 'DOMA' $DOMZ   'VF' ;
 QDOMX   = TDOMX    . 'QUAF' ;
 QDOMY   = TDOMY    . 'QUAF' ;
 QDOMZ   = TDOMZ    . 'QUAF' ;
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMX ;
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMY ;
 'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMZ ;
 
 LIGEX = (TDOMX . 'CENTRE') ;
*  P1 = 'POIN' (COOR 1 LIGEX) 'MINIMUM' ;
*  'ORDONNER' LIGEX P1 ;
 'ORDONNER' LIGEX ;
 LIGEX = 'QUELCONQUE' 'SEG2' LIGEX 'COULEUR' 'VERT';
*
 LIGEY = (TDOMY . 'CENTRE') ;
*  P1 = 'POIN' (COOR 2 LIGEY) 'MINIMUM' ;
*  'ORDONNER' LIGEY P1 ;
 'ORDONNER' LIGEY ;
 LIGEY = 'QUELCONQUE' 'SEG2' LIGEY 'COULEUR' 'ROSE';
*
 LIGEZ = (TDOMZ . 'CENTRE') ;
*  P1 = 'POIN' (COOR 3 LIGEZ) 'MINIMUM' ;
*  'ORDONNER' LIGEZ P1 ;
 'ORDONNER' LIGEZ ;
 LIGEZ = 'QUELCONQUE' 'SEG2' LIGEZ 'COULEUR' 'ROUG';

 'SI' GRAPH ;
    'TRACER' (QDOMTOT 'ET' LIGEX 'ET' LIGEY 'ET' LIGEZ) ;
 'FINSI' ;

* 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;

****************************************************************
****************************************************************
*****      Initial conditions and gas properties.          *****
****************************************************************
****************************************************************
*
*************************************************
**** The table for the properties of the gas ****
*************************************************
*
 PGAS = 'TABLE' ;
*
**** Order of the polynomial order for cv = cv(T) 
*    For T > TMAX, cv(T) = cv(Tmax)
* 
 PGAS . 'TMAX' = 6000.0 ;
 PGAS . 'NORD' = 4 ;
*
**** Species involved in the mixture (before or after
*    the chemical reaction)
*
 PGAS . 'SPECIES' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' 'N2  ' ;
*
*
**** Coefficient of the chemical reaction.
*    Note that for the first species this coefficient should be positive
*    Normal, we take it equal to 1.
*
*    H2 '+' 0.5 O2 ---> H2O
*
 PGAS . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 0.0 ;
* 
**** Mass fraction of the first species before and after the combustion
*    Final mass fractions of the species with positive coefficients.
*    Final mass fractions of the species with non-positive coefficient.
*    The mass fraction of the last species is not given.
*
 PGAS . 'MASSFRA' = 'PROG' 0.285219E-01  0.964039E-11  0.765104E-10
                           0.127442E-10 ;
*
**** Coef with the gas properties
*
 PGAS .  'H2  ' = 'TABLE'  ;
 PGAS .  'H2O ' = 'TABLE'  ;
 PGAS .  'N2  ' = 'TABLE'  ;
 PGAS .  'O2  ' = 'TABLE'  ;
*
**** Runiv (J/mole/K)
*
 PGAS . 'RUNIV' = 8.31441 ;
*
**** W (kg/mole). Gas constant (J/kg/K = Runiv/W)
*
 PGAS .  'H2  ' . 'W' =  2. * 1.00797E-3 ;
 PGAS .  'O2  ' . 'W' =  2. * 15.9994E-3 ;
 PGAS .  'H2O ' . 'W' = (PGAS .  'H2  ' . 'W' ) '+'
    (0.5 * (PGAS .  'O2  ' . 'W' )) ;
 PGAS .  'N2  ' . 'W' = 2 * 14.0067E-3 ;
*
**** Polynomial coefficients 
*
 PGAS .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
 PGAS .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
                              -1.82753232E-08 2.44485692E-12 ;
 PGAS .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
                              8.78233606E-09 -3.05514485E-13 ;
 PGAS .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
                              2.33636971E-08 -1.53304905E-12;
*
**** Formation enthalpies (energies) at 0K (J/Kg)
*
 PGAS .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAS .  'H2O ' . 'H0K' = -1.395D7 ;
 PGAS .  'N2  ' . 'H0K' = -2.953D5 ;
 PGAS .  'O2  ' . 'H0K' = -2.634D5 ;
*
*************************************************
**** The initial conditions *********************
*************************************************
*
 eps = 1.0D-64 ;
 K0 = 200. ;
* 
 T1   = 293. ;
 alph11 = 1 '-' 0.0001 ;
 alph12 =  0.0001 ;
 alph21 = 1.0 '-' alph11 ;
 alph22 = 1.0 - alph12 ;
 T2   = 2800. ;
 un1  = 100.;
 un2  = 100.;
 ut1  = 100.;
 ut2  = 100.;
 uv1  = 100.;
 uv2  = 100.;
 pre1 = 1.013D5 ;
 pre2 = 2.013D5 ;
*
 CHVN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 3 'UX' un1
   'UY' ut1 'UZ' uv1) ;
 CHVN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 3 'UX' un2
   'UY' ut2 'UZ' uv2) ;
*
 CHTN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T1) ;
 CHTN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T2) ;
*
 CHPN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre1) ;
 CHPN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre2) ;
* CHPN1 = CHPN1 '+' ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' pre1) ;
* 'ERREUR' 21 ;
*
 CHAL1  = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
   alph11) '+'
   ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
   alph12) ;
 CHAL2  = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
   alph21) '+'
   ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
   alph22) ;
*
 CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS
   PGAS CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;
* 
 RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS CHAL1 CHAL2
   (CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1)
   (CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2)
    CHTN1 CHTN2 EPS ;
 TN1M = COPIER TN1 ;
 TN2M = COPIER TN2 ;
*
 ERRO = 'MAXIMUM' (PRE1 '-' ('REDU' PN1 (TDOM1 . 'CENTRE'))) 'ABS' ;
 ERRO = ERRO '/' (pre1) ;
 'SI' (ERRO > 1.0D-6) ;
     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
     'ERREUR' 21 ;
 'FINSI' ; 
 ERRO = 'MAXIMUM' (PRE2 '-' ('REDU' PN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
 ERRO = ERRO '/' (pre2) ;
 'SI' (ERRO > 1.0D-6) ;
     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
     'ERREUR' 21 ;
 'FINSI' ;
 ERRO = 'MAXIMUM' (T1 '-' ('REDU' TN1 (TDOM1 . 'CENTRE'))) 'ABS' ;
 ERRO = ERRO '/' (T1) ;
 'SI' (ERRO > 1.0D-6) ;
     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
     'ERREUR' 21 ;
 'FINSI' ;
 ERRO = 'MAXIMUM' (T2 '-' ('REDU' TN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
 ERRO = ERRO '/' (T2) ;
 'SI' (ERRO > 1.0D-6) ;
     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
     'ERREUR' 21 ;
 'FINSI' ;
*
**** Parameter for the time loop 
*
*    Names of the residuum components
*
 LISTINCO = ('MOTS' 'ALF1' 'RN1' 'RNX1' 'RNY1' 'RNZ1' 'RET1'
          'ALF2' 'RN2' 'RNX2' 'RNY2' 'RNZ2' 'RET2') ;
*
**** BC
*
 KONLIM = 'DIFF' DOMTOT DOMTOT ;
 CHPVID CHMVID = 'KOPS' 'MATRIK' ;
 RESLIM = 'COPIER' CHPVID ;
 
* 
****************************************************************
****************************************************************
*****           Parameters for the computations           ******
****************************************************************
****************************************************************
* 
 
* Iterations
* Final time
* Safety Factor fot the time step
* Second order reconstruction?
* We compute the gradients during the correction or not ?

 NITER = 10000 ;
 TFINAL = 0.4D-3 ;
 SAFFAC = 0.25 ;
 LOGSO = VRAI ;
 LOGGRA = VRAI ;


* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER'  METO  NITER  TFINAL  SAFFAC LOGSO ;

****************************************************************
****************************************************************
*****           The computation                           ******
****************************************************************
****************************************************************

 AL1   = CHAL1          ;
 AL2   = CHAL2          ;
 ARN1  = CHAL1 * CHRN1  ;
 ARN2  = CHAL2 * CHRN2  ;
 AGN1  = CHAL1 * CHGN1  ;
 AGN2  = CHAL2 * CHGN2  ;
 AREN1 = CHAL1 * CHRET1 ;
 AREN2 = CHAL2 * CHRET2 ;
*
 AL10   =  'COPIER' AL1   ;
 AL20   =  'COPIER' AL2   ;
 ARN10  =  'COPIER' ARN1  ;
 ARN20  =  'COPIER' ARN2  ;
 AGN10  =  'COPIER' AGN1  ;
 AGN20  =  'COPIER' AGN2  ;
 AREN10 =  'COPIER' AREN1 ;
 AREN20 =  'COPIER' AREN2 ;
 TN10   =  'COPIER' TN1   ;
 TN20   =  'COPIER' TN2   ;
 
*
**** Geometrical coefficient to compute grad(alpha)/|grad(alpha)|
*

 EPSGR = 1.0D-6 ;
 CHPL1 = CHPVID ;
 CHPL2 = 'MANUEL' 'CHPO' (TDOMLIM . 'CENTRE') 3
         'P1DX' 0.0 'P1DY'  0.0 'P1DZ'  0.0 ;
 GRALP1 GRAL = 'PENT' $DOMTOT 'FACE'
     'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
     AL1  CHPL1 CHPL2 ;
* 'SI' VRAI ;
*    V1 = 'VECTEUR'
*       ('NOMC' ('MOTS' 'P1DX' 'P1DY' 'P1DZ') GRALP1
*       ('MOTS' 'UX' 'UY' 'UZ')) ;
*    'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
 GRALP1 = GRALP1 '+' EPSGR ;
 GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1
    ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
    ('MOTS' 'P1DX' 'P1DY' 'P1DZ')) '**' 0.5) ;
* 'SI' VRAI ;
*    V1 = 'VECTEUR'
*       ('NOMC' ('MOTS' 'P1DX' 'P1DY' 'P1DZ') GRALP1
*       ('MOTS' 'UX' 'UY' 'UZ')) ;
*    'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
 
*
**** Geometrical coefficient to compute gradients
*
 GRADR0 ALR0 COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
          ('MOTS' 'SCAL') AL10 ;
 GRADR0 = GRADR0 * 0.0 ;
 ALR0   = ALR0   * 0.0 ;
          
 GRADV0 ALV0 COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
          ('MOTS' 'UX' 'UY' 'UZ') AGN10 ;
 GRADV0 = GRADV0 * 0.0 ;
 ALV0   = ALV0   * 0.0 ;

 VINF1 = 'MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' 10. ;
 VINF2 = 'MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' 100. ;
*
**** Temporal loop
*

 TPS = 0.0 ;
 
 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
 'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
 'MESSAGE' ;

 'TEMPS' 'ZERO' ;
 'REPETER' BLITER NITER ;

*
**** Primitive variables
*
    
    RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
       ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;

    TN1M = 'COPIER' TN1 ;
    TN2M = 'COPIER' TN2 ;

**** Gradient of alpha

    GRALP1 = 'PENT' $DOMTOT 'FACE'
    'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
        AL1  CHPL1 CHPL2 'GRADGEO' GRAL  ; 
*    'SI' VRAI ;
*       V1 = 'VECTEUR'
*          ('NOMC' ('MOTS' 'P1DX' 'P1DY' 'P1DZ') GRALP1
*          ('MOTS' 'UX' 'UY' 'UZ')) ;
*       'TRACER' DOMTOT V1 DOMLIM ;
*    'FINSI' ;
    GRALP1 = GRALP1 '+' EPSGR ;
    GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1
       ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
       ('MOTS' 'P1DX' 'P1DY' 'P1DZ')) '**' 0.5) ;
    
    'SI' LOGSO ;
    
       GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') AL1  'GRADGEO' COEFSCAL ;
          
       GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN1  'GRADGEO' COEFSCAL ;

       GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN1  'GRADGEO' COEFSCAL ;

       GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY' 'UZ') VN1   'GRADGEO'  COEFVECT ;
          
       GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') AL2  'GRADGEO' COEFSCAL ;
          
       GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN2  'GRADGEO' COEFSCAL ;

       GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN2  'GRADGEO' COEFSCAL ;

       GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY' 'UZ') VN2   'GRADGEO'  COEFVECT ;

       CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
         'PRET' 'DEM'  $DOMTOT
          AL1  GRADA1  ALA1
          AL2  GRADA2  ALA2
*          AL1  GRADA1  ALR0
*          AL2  GRADA2  ALR0
          RN1  GRADR1  ALR1
          RN2  GRADR2  ALR2 
          VN1  GRADV1  ALV1
          VN2  GRADV2  ALV2
          PN1  GRADP1  ALP1
          PN2  GRADP2  ALP2 ;
                
    'SINON' ;
                              
       CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
         'PRET' 'DEM'  $DOMTOT
          AL1  GRADR0  ALR0
          AL2  GRADR0  ALR0
          RN1  GRADR0  ALR0
          RN2  GRADR0  ALR0
          VN1  GRADV0  ALV0
          VN2  GRADV0  ALV0
          PN1  GRADR0  ALR0
          PN2  GRADR0  ALR0 ;

    'FINSI' ;

    SI ('EGA' METO 'AUSMPUP') ;
       RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
          $DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
          CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM VINF1 VINF2 ;
    'SINON' ;
       RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
          $DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
          CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
    'FINSI' ;
*
    RESIDU = RESIDU '+' RESLIM ;
*    
*    'REPETER' BL1 ('DIME' LISTINCO) ;  
*       mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
*       valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
*       tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
*       evaa = 'EVOL' 'CHPO' RESIDU mot1  LIGEVO ;
*       'DESSIN' evaa 'TITRE' tit1 ;
*    'FIN' BL1 ;
 
*   Problem with RNY1, RNY2
                              
    DT_CON = SAFFAC '*' DELTAT ;     
*
**** The time step linked to tps
*

    DTTPS = (TFINAL '-' TPS) * (1. '+' 1.0D-6) ;
    
*
**** Total time step
*

    DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;

*
**** Increment of the variables (predictor)
*
*    RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
*       ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 CHTN1 CHTN2 EPS ;

    RESIDU = (0.5 * DTMIN) '*' RESIDU ;

    DALP1  = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
    DRN1   = 'EXCO' 'RN1'  RESIDU 'SCAL' ;
    DGN1   = 'EXCO' ('MOTS' 'RNX1' 'RNY1' 'RNZ1') RESIDU
       ('MOTS' 'UX' 'UY' 'UZ') ;
    DRET1  = 'EXCO' 'RET1' RESIDU 'SCAL' ;
    DALP2  = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
    DRN2   = 'EXCO' 'RN2'  RESIDU 'SCAL' ;
    DGN2   = 'EXCO' ('MOTS' 'RNX2' 'RNY2' 'RNZ2') RESIDU
       ('MOTS' 'UX' 'UY' 'UZ') ;
    DRET2  = 'EXCO' 'RET2' RESIDU 'SCAL' ;

    AL1B   = AL1   '+'  DALP1 ;
    ARN1B  = ARN1  '+'  DRN1  ;
    AGN1B  = AGN1  '+'  DGN1  ;
    AREN1B = AREN1 '+'  DRET1 ;
    AL2B   = AL2   '+'  DALP2 ;
    ARN2B  = ARN2  '+'  DRN2  ;
    AGN2B  = AGN2  '+'  DGN2  ;
    AREN2B = AREN2 '+'  DRET2 ;

*    
**** Corrector
*

    RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1B AL2B
       ARN1B ARN2B AGN1B AGN2B AREN1B AREN2B TN1M TN2M EPS ;

**** Gradient of alpha

    GRALP1 = 'PENT' $DOMTOT 'FACE'
    'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
        AL1B  CHPL1 CHPL2 'GRADGEO' GRAL  ; 
*    'SI' VRAI ;
*       V1 = 'VECTEUR'
*          ('NOMC' ('MOTS' 'P1DX' 'P1DY' 'P1DZ') GRALP1
*          ('MOTS' 'UX' 'UY' 'UZ')) ;
*       'TRACER' DOMTOT V1 DOMLIM ;
*    'FINSI' ;
    GRALP1 = GRALP1 '+' EPSGR ;
    GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1
      ('MOTS' 'P1DX' 'P1DY' 'P1DZ')
       ('MOTS' 'P1DX' 'P1DY' 'P1DZ')) '**' 0.5) ;
       
    'SI' LOGSO ;

       SI LOGGRA ;
       GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') AL1B  'GRADGEO' COEFSCAL ;
          
       GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN1  'GRADGEO' COEFSCAL ;

       GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN1  'GRADGEO' COEFSCAL ;

       GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY' 'UZ') VN1   'GRADGEO'  COEFVECT ;
          
       GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') AL2B  'GRADGEO' COEFSCAL ;
          
       GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN2  'GRADGEO' COEFSCAL ;

       GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN2  'GRADGEO' COEFSCAL ;

       GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY' 'UZ') VN2   'GRADGEO'  COEFVECT ;
       'FINSI' ;

       CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
         'PRET' 'DEM'  $DOMTOT
          AL1B  GRADA1  ALA1
          AL2B  GRADA2  ALA2
*          AL1B  GRADA1  ALR0
*          AL2B  GRADA2  ALR0
          RN1  GRADR1  ALR1
          RN2  GRADR2  ALR2 
          VN1  GRADV1  ALV1
          VN2  GRADV2  ALV2
          PN1  GRADP1  ALP1
          PN2  GRADP2  ALP2 ;
                
    'SINON' ;
                              
       CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
         'PRET' 'DEM'  $DOMTOT
          AL1B  GRADR0  ALR0
          AL2B  GRADR0  ALR0
          RN1  GRADR0  ALR0
          RN2  GRADR0  ALR0
          VN1  GRADV0  ALV0
          VN2  GRADV0  ALV0
          PN1  GRADR0  ALR0
          PN2  GRADR0  ALR0 ;

    'FINSI' ;

    SI ('EGA' METO 'AUSMPUP') ;
       RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
          $DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
          CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM VINF1 VINF2 ;
    'SINON' ;
       RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
          $DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
          CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
    'FINSI' ;

    RESIDU = RESIDU '+' RESLIM ;
*    
*    'REPETER' BL1 ('DIME' LISTINCO) ;  
*       mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
*       valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
*       tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
*       evaa = 'EVOL' 'CHPO' RESIDU mot1  LIGEVO ;
*       'DESSIN' evaa 'TITRE' tit1 ;
*    'FIN' BL1 ;
 
*    'OPTION' DONN 5 ;
*   Problem with RNY1, RNY2
                              

    RESIDU = DTMIN '*' RESIDU ;

    DALP1  = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
    DRN1   = 'EXCO' 'RN1'  RESIDU 'SCAL' ;
    DGN1   = 'EXCO' ('MOTS' 'RNX1' 'RNY1' 'RNZ1') RESIDU
       ('MOTS' 'UX' 'UY' 'UZ') ;
    DRET1  = 'EXCO' 'RET1' RESIDU 'SCAL' ;
    DALP2  = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
    DRN2   = 'EXCO' 'RN2'  RESIDU 'SCAL' ;
    DGN2   = 'EXCO' ('MOTS' 'RNX2' 'RNY2' 'RNZ2')
       RESIDU ('MOTS' 'UX' 'UY' 'UZ') ;
    DRET2  = 'EXCO' 'RET2' RESIDU 'SCAL' ;
    
    TPS   = TPS   '+'  DTMIN ;
    AL1   = AL1   '+'  DALP1 ;
    ARN1  = ARN1  '+'  DRN1  ;
    AGN1  = AGN1  '+'  DGN1  ;
    AREN1 = AREN1 '+'  DRET1 ;
    AL2   = AL2   '+'  DALP2 ;
    ARN2  = ARN2  '+'  DRN2  ;
    AGN2  = AGN2  '+'  DGN2  ;
    AREN2 = AREN2 '+'  DRET2 ;

    'SI' (((&BLITER '/' 20) '*' 20) 'EGA' &BLITER) ;
        'MESSAGE' ('CHAINE' 'ITER =' &BLITER '  TPS =' TPS) ;
    'FINSI' ;

    'SI' (TPS > TFINAL) ;
       'QUITTER' BLITER ;
    'FINSI' ;

 'FIN' BLITER ;
* 'TEMPS' 'IMPR' ;
 'TEMPS' ;

* 'OPTION' 'SAUVER' ('CHAINE' 'result.sauv_' RAF 'tps_' TPS ) ;
* 'SAUVER'  ;

****************************************************************
****************************************************************
*****           The post-treatment                        ******
****************************************************************
****************************************************************

* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;

* 
**** The mesh
*

 'SI' GRAPH ;
    'TRACER' ('DOMA' $DOMTOT 'MAILLAGE')
         'TITR' ('CHAINE' 'Maillage: ' ('NBEL' DOMTOT) ' elts');
 'FINSI' ;

*
**** Initial conditions
*

  MOD1 =  'MODELISER'  ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
 'SI' (VRAI  ET GRAPH) ;

    RN10 RN20 VN10 VN20 PN10 PN20 TN10 TN20 = 'PRIM' 'DEM' PGAS
      AL10 AL20 ARN10 ARN20 AGN10 AGN20 AREN10 AREN20 TN1M TN2M EPS ;

    CHM_AL10  =  'KCHA'  $DOMTOT  'CHAM'  AL10  ;
    CHM_RN10  =  'KCHA'  $DOMTOT  'CHAM'  RN10  ;
    CHM_VN10  =  'KCHA'  $DOMTOT  'CHAM'  VN10  ;
    CHM_PN10  =  'KCHA'  $DOMTOT  'CHAM'  PN10  ;
    CHM_TN10  =  'KCHA'  $DOMTOT  'CHAM'  TN10  ;
    CHM_AL20  =  'KCHA'  $DOMTOT  'CHAM'  AL20  ;
    CHM_RN20  =  'KCHA'  $DOMTOT  'CHAM'  RN20  ;
    CHM_VN20  =  'KCHA'  $DOMTOT  'CHAM'  VN20  ;
    CHM_PN20  =  'KCHA'  $DOMTOT  'CHAM'  PN20  ;
    CHM_TN20  =  'KCHA'  $DOMTOT  'CHAM'  TN20  ;
    
    'TRAC' CHM_AL10  MOD1 
       'TITR' ('CHAINE' 'alp_1 at t=' 0.0)   ;
    'TRAC' CHM_RN10  MOD1 
       'TITR' ('CHAINE' 'rho_1 at t=' 0.0)   ;
    'TRAC' CHM_VN10  MOD1 
       'TITR' ('CHAINE' 'v_1 at t= '  0.0)   ;
    'TRAC' CHM_PN10  MOD1 
       'TITR' ('CHAINE' 'p_1 at t= '  0.0)   ;
    'TRAC' CHM_TN10  MOD1 
       'TITR' ('CHAINE' 'T_1 at t= '  0.0)   ;
    'TRAC' CHM_AL20  MOD1 
       'TITR' ('CHAINE' 'alp_2 at t=' 0.0)   ;
    'TRAC' CHM_RN20  MOD1 
       'TITR' ('CHAINE' 'rho_2 at t=' 0.0)   ;
    'TRAC' CHM_VN20  MOD1 
       'TITR' ('CHAINE' 'v_2 at t= '  0.0)   ;
    'TRAC' CHM_PN20  MOD1 
       'TITR' ('CHAINE' 'p_2 at t= '  0.0)   ;
    'TRAC' CHM_TN20  MOD1 
       'TITR' ('CHAINE' 'T_2 at t= '  0.0)   ;
    
 'FINSI' ;

*
**** Results
*

 RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
      ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;
  
 'SI' (VRAI ET GRAPH) ;
 
    CHM_AL1  =  'KCHA'  $DOMTOT  'CHAM'  AL1  ;
    CHM_RN1  =  'KCHA'  $DOMTOT  'CHAM'  RN1  ;
    CHM_VN1  =  'KCHA'  $DOMTOT  'CHAM'  VN1  ;
    CHM_PN1  =  'KCHA'  $DOMTOT  'CHAM'  PN1  ;
    CHM_TN1  =  'KCHA'  $DOMTOT  'CHAM'  TN1  ;
    CHM_AL2  =  'KCHA'  $DOMTOT  'CHAM'  AL2  ;
    CHM_RN2  =  'KCHA'  $DOMTOT  'CHAM'  RN2  ;
    CHM_VN2  =  'KCHA'  $DOMTOT  'CHAM'  VN2  ;
    CHM_PN2  =  'KCHA'  $DOMTOT  'CHAM'  PN2  ;
    CHM_TN2  =  'KCHA'  $DOMTOT  'CHAM'  TN2  ;
    
    'TRAC' CHM_AL1  MOD1 
       'TITR' ('CHAINE' 'alp_1 at t=' TPS)   ;
    'TRAC' CHM_RN1  MOD1 
       'TITR' ('CHAINE' 'rho_1 at t=' TPS)   ;
    'TRAC' CHM_VN1  MOD1 
       'TITR' ('CHAINE' 'v_1 at t= ' TPS)   ;
    'TRAC' CHM_PN1  MOD1 
       'TITR' ('CHAINE' 'p_1 at t= ' TPS)   ;
    'TRAC' CHM_TN1  MOD1 
       'TITR' ('CHAINE' 'T_1 at t= ' TPS)   ;
    'TRAC' CHM_AL2  MOD1 
       'TITR' ('CHAINE' 'alp_2 at t=' TPS)  ;
    'TRAC' CHM_RN2  MOD1 
       'TITR' ('CHAINE' 'rho_2 at t=' TPS)  ;
    'TRAC' CHM_VN2  MOD1 
       'TITR' ('CHAINE' 'v_2 at t= ' TPS)   ;
    'TRAC' CHM_PN2  MOD1 
       'TITR' ('CHAINE' 'p_2 at t= ' TPS)   ;
    'TRAC' CHM_TN2  MOD1 
       'TITR' ('CHAINE' 'T_2 at t= ' TPS)   ;
    
 'FINSI' ;

*
*** Evolution Objects
*
* LIGEX
 evxal1 = ('EVOL' 'CHPO' al1 LIGEX) 'COULEUR' 'ROUG' ;
 evxal2 = ('EVOL' 'CHPO' al2 LIGEX) 'COULEUR' 'BLEU' ;
 evxone = (evxal1 '+' evxal2) 'COULEUR' 'VERT' ;

 evxrn1 = ('EVOL' 'CHPO' rn1 LIGEX) 'COULEUR' 'ROUG' ;
 evxrn2 = ('EVOL' 'CHPO' rn2 LIGEX) 'COULEUR' 'BLEU' ;

 evxpn1 = ('EVOL' 'CHPO' pn1 LIGEX) 'COULEUR' 'ROUG' ;
 evxpn2 = ('EVOL' 'CHPO' pn2 LIGEX) 'COULEUR' 'BLEU' ;
 
 evxux1 = ('EVOL' 'CHPO' vn1 LIGEX 'UX') 'COULEUR' 'ROUG' ;
 evxux2 = ('EVOL' 'CHPO' vn2 LIGEX 'UX') 'COULEUR' 'BLEU' ;
 
 evxuy1 = ('EVOL' 'CHPO' vn1 LIGEX 'UY') 'COULEUR' 'ROUG' ;
 evxuy2 = ('EVOL' 'CHPO' vn2 LIGEX 'UY') 'COULEUR' 'BLEU' ;

 evxuz1 = ('EVOL' 'CHPO' vn1 LIGEX 'UZ') 'COULEUR' 'ROUG' ;
 evxuz2 = ('EVOL' 'CHPO' vn2 LIGEX 'UZ') 'COULEUR' 'BLEU' ;

 evxtn1 = ('EVOL' 'CHPO' tn1 LIGEX) 'COULEUR' 'ROUG' ;
 evxtn2 = ('EVOL' 'CHPO' tn2 LIGEX) 'COULEUR' 'BLEU' ;
* LIGEY
 evyal1 = ('EVOL' 'CHPO' al1 LIGEY) 'COULEUR' 'ROUG' ;
 evyal2 = ('EVOL' 'CHPO' al2 LIGEY) 'COULEUR' 'BLEU' ;
 evyone = (evyal1 '+' evyal2) 'COULEUR' 'VERT' ;

 evyrn1 = ('EVOL' 'CHPO' rn1 LIGEY) 'COULEUR' 'ROUG' ;
 evyrn2 = ('EVOL' 'CHPO' rn2 LIGEY) 'COULEUR' 'BLEU' ;

 evypn1 = ('EVOL' 'CHPO' pn1 LIGEY) 'COULEUR' 'ROUG' ;
 evypn2 = ('EVOL' 'CHPO' pn2 LIGEY) 'COULEUR' 'BLEU' ;
 
 evyux1 = ('EVOL' 'CHPO' vn1 LIGEY 'UX') 'COULEUR' 'ROUG' ;
 evyux2 = ('EVOL' 'CHPO' vn2 LIGEY 'UX') 'COULEUR' 'BLEU' ;
 
 evyuy1 = ('EVOL' 'CHPO' vn1 LIGEY 'UY') 'COULEUR' 'ROUG' ;
 evyuy2 = ('EVOL' 'CHPO' vn2 LIGEY 'UY') 'COULEUR' 'BLEU' ;

 evyuz1 = ('EVOL' 'CHPO' vn1 LIGEY 'UZ') 'COULEUR' 'ROUG' ;
 evyuz2 = ('EVOL' 'CHPO' vn2 LIGEY 'UZ') 'COULEUR' 'BLEU' ;

 evytn1 = ('EVOL' 'CHPO' tn1 LIGEY) 'COULEUR' 'ROUG' ;
 evytn2 = ('EVOL' 'CHPO' tn2 LIGEY) 'COULEUR' 'BLEU' ;
* LIGEZ
 evzal1 = ('EVOL' 'CHPO' al1 LIGEZ) 'COULEUR' 'ROUG' ;
 evzal2 = ('EVOL' 'CHPO' al2 LIGEZ) 'COULEUR' 'BLEU' ;
 evzone = (evzal1 '+' evzal2) 'COULEUR' 'VERT' ;

 evzrn1 = ('EVOL' 'CHPO' rn1 LIGEZ) 'COULEUR' 'ROUG' ;
 evzrn2 = ('EVOL' 'CHPO' rn2 LIGEZ) 'COULEUR' 'BLEU' ;

 evzpn1 = ('EVOL' 'CHPO' pn1 LIGEZ) 'COULEUR' 'ROUG' ;
 evzpn2 = ('EVOL' 'CHPO' pn2 LIGEZ) 'COULEUR' 'BLEU' ;
 
 evzux1 = ('EVOL' 'CHPO' vn1 LIGEZ 'UX') 'COULEUR' 'ROUG' ;
 evzux2 = ('EVOL' 'CHPO' vn2 LIGEZ 'UX') 'COULEUR' 'BLEU' ;
 
 evzuy1 = ('EVOL' 'CHPO' vn1 LIGEZ 'UY') 'COULEUR' 'ROUG' ;
 evzuy2 = ('EVOL' 'CHPO' vn2 LIGEZ 'UY') 'COULEUR' 'BLEU' ;

 evzuz1 = ('EVOL' 'CHPO' vn1 LIGEZ 'UZ') 'COULEUR' 'ROUG' ;
 evzuz2 = ('EVOL' 'CHPO' vn2 LIGEZ 'UZ') 'COULEUR' 'BLEU' ;

 evztn1 = ('EVOL' 'CHPO' tn1 LIGEZ) 'COULEUR' 'ROUG' ;
 evztn2 = ('EVOL' 'CHPO' tn2 LIGEZ) 'COULEUR' 'BLEU' ;
*
 erro1 = (evxtn1 '-' evytn1) ;
 erro2 = (evxtn1 '-' evztn1) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 1, problem with T' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
 erro1 = (evxtn2 '-' evytn2) ;
 erro2 = (evxtn2 '-' evztn2) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 2, problem with T' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
 erro1 = (evxpn1 '-' evypn1) ;
 erro2 = (evxpn1 '-' evzpn1) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 1, problem with P' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
 erro1 = (evxpn2 '-' evypn2) ;
 erro2 = (evxpn2 '-' evzpn2) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 2, problem with P' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
*
 erro1 = (evxux1 '-' evyuy1) ;
 erro2 = (evxux1 '-' evzuz1) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 1, problem with un' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
 erro1 = (evxux2 '-' evyuy2) ;
 erro2 = (evxux2 '-' evzuz2) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 2, problem with un' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
*
 erro1 = (evxuy1 '-' evyuz1) ;
 erro2 = (evxuy1 '-' evzux1) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 1, problem with un' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
 erro1 = (evxuy2 '-' evyuz2) ;
 erro2 = (evxuy2 '-' evzux2) ;
 erro1 = 'MAXIMUM' ('ABS' ('EXTR' erro1 'ORDO')) ;
 erro2 = 'MAXIMUM' ('ABS' ('EXTR' erro2 'ORDO')) ;
 erro = 'MAXIMUM' ('PROG' erro1 erro2) ;
 SI (ERRO > 1.0D-8) ;
    'Phase 2, problem with un' ;
    'ERREUR' 21 ;
 'FINSI' ;
* 
 'FIN' ;