* fichier :  crebe12.dgibi
************************************************************************
* Section : Fluides Euler
************************************************************************
*--------------------------------------------
* GIBIANE file for the HDR test E12.3.2     *
*               in 2D                       *
* We make 100 iterations on the coarse mesh *
* and compare the overpressure in the room  *
* R1904 with earlier computed values        *
*  CREBCOM combustion model is used         * 
*--------------------------------------------
*--------------------------------------------
*   Here we construct the mesh for the      *
*           HDR test E12.3.2                *
*--------------------------------------------

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

*-------------------------------------------
*  Given a line and a mesh, this procedure
*  computes a new line which containes the
*  points of the mesh closest to the points
*  of the given line
*-------------------------------------------
 'DEBPROC' LIGPROC ;
******* inputs
    'ARGUMENT' LIGG*'MAILLAGE' MAI1*'MAILLAGE' ;
    NPOI = 'NBNO' LIGG ;
    MAI2 = 'CHANGER' MAI1 'POI1' ;
* First element of the line
    P1 = LIGG 'POIN' 1 ;
    P1P = MAI2 'POIN' 'PROC' P1 ;
    P2 = LIGG 'POIN' 2 ;
    P2P = MAI2 'POIN' 'PROC' P2 ;
    'REPETER' BL1 (NPOI '-' 1) ;
       'SI' ('NEG' P1P P2P) ;
          NPOI = (NPOI '-' &BL1) '+' 1 ;
          'QUITTER' BL1 ;
       'FINSI' ;
       P2 = LIGG 'POIN' (&BL1 '+' 2) ;
       P2P = LIGG 'POIN' 'PROC' P2 ;
    'FIN' BL1 ;
    LIGRES = 'MANUEL' 'SEG2' P1P P2P ;
* The other element of the line
    'REPETER' BL1 (NPOI '-' 2) ;
       P1P = P2P ;
       P2 = LIGG 'POIN' (&BL1 '+' 2) ;
       P2P = MAI2 'POIN' 'PROC' P2 ;
      'SI' ('NEG' P1P P2P) ;
        LIGRES = LIGRES 'ET' ('MANUEL' 'SEG2' P1P P2P) ;
      'FINSI' ;
    'FIN' BL1 ;
    'RESPRO' LIGRES ;
 'FINPROC' ;
*--------------------------------------------------------
*---refinement of the mesh (10 or 16)   
 RAF = 10 ;
*-- Length of the thirst and second room
 L1 = 7.32D0 ;
*-- Height of the first window
 H1 = 0.5D0 ;
*-- Height of the first and second room
 H2 = 4.0D0 ;    
*---------------------------------
* Number of cells in each region
*---------------------------------   
 N1 = RAF ;
*---- number of points in window1
** It is BETTER to have NY1 even!
 NY1 = 'ENTIER' (RAF '/' 4) ;
*--------------------------------- 
 DX = L1 '/' N1 ;
*---- cell size in window 
 DX1 = H1 '/' NY1 ;
*---- number of points in the outer window
 NY2 = 'ENTIER' ((H2 '-' H1) '/' DX) ;
***---- wall thickness
 DW = H1 ; 
*-----------------------------------------
*  First and Second Room (R1904 and R1905)
*-----------------------------------------
 LB = L1 '-' (7.5 '*' H1) ;
 A1 = 0.0 0.0 ;
 A12 = LB 0.0 ;
 A13 = LB (H1 '/' 2) ;
 A14 = 0.0 (H1 '/' 2) ;
*------------------------------------
*  "Lower-left" domain
*-------------------------------------
 NY1M = NY1 '/' 2 ;
 A1A12 = A1 'DROIT'  A12 'DINI' DX1 'DFIN' (DW '/' NY1) ;
 A12A13 = A12 'DROIT'  NY1M A13 ;
 A13A14 = A13 'DROIT'  A14 'DINI' (DW '/' NY1) 'DFIN' DX1 ; 
 A14A1 = A14 'DROIT'  NY1M A1 ;

 DOM1L = 'DALLER' A1A12 A12A13 A13A14 A14A1 'PLAN' ;
*--------------------------------------
* "Lower-right" domain
*--------------------------------------
 LB1 = L1 '-' (2.5 '*' H1) ;
 AI1 = LB1 0.0 ;
 A2 = L1 0.0 ;
 AI3 = L1 (H1 '/' 2) ;
 AI4 = LB1 (H1 '/' 2) ;
*------------------------------------
 DENF = (DW '/' NY1) ;
 AI1A2 = AI1 'DROIT' A2 'DINI' DENF 'DFIN' DENF ;
 A2AI3 = A2 'DROIT'  NY1M AI3 ;
 AI3AI4 = AI3 'DROIT' AI4 'DINI' DENF 'DFIN' DENF ;
 AI4AI1 = AI4 'DROIT'  NY1M AI1 ;

 DOM1R = 'DALLER' AI1A2 A2AI3 AI3AI4 AI4AI1 'PLAN' ;

 DOM1  = DOM1L 'ET' DOM1R ;
 'ELIMINATION'  DOM1 1.0D-6 ;
*----------------------------------------------------
*  Line for postreatment
*----------------------------------------------------
*** NNN is the number of points in the disc
 NNN = NY1 '*' 4 ;
 xcel = LB '+' (H1 '/' 2) ;
 ycel = H1 '/' 2 ;
 ACEL = xcel ycel ;
 LIGB = A13 'DROIT' NY1M ACEL  ;
*------------------
  xmm = LB '+' (2.5 '*' H1) ;
 'REPETER' BL1 NNN ;
   ACEL0 = ACEL ;
   xcel = xcel '+' DX1 ;
   dop = 18.0625 '*' H1 '*' H1 ;
   dop = dop '-' ((xcel-xmm) '*' (xcel-xmm)) ;
   dop = dop '**' 0.5 ;
   ycel = dop '-' (3.25 '*' H1) ;
   ACEL = xcel ycel ;
   LIGB = LIGB 'ET' (ACEL0 'DROIT' 1 ACEL) ;
 'FIN'  BL1 ;
 LIGB = LIGB 'ET' (ACEL 'DROIT' AI4 'DINI' DX1 'DFIN' DX1 ) ;
*--------------------------------------------
 LIGBG = ('INVERSE' A13A14) 'ET' LIGB 'ET'
     ('INVERSE' AI3AI4) ;
*-------------------------------------
*  "HIGH" domain
*-------------------------------------
 A3 = L1 H1 ;
 A4 = 0.0 H1 ;
 A5 = L1 H2 ;
 A6 = 0.0 H2 ;
 AI5 = LB H2 ;
*---------------------------------
 AI3A5 = AI3 'DROIT'  A5 'DINI' DX1 'DFIN' DX1 ;
 A5AI5 = A5 'DROIT' AI5 'DINI' DX1 'DFIN' DX1 ;
 AI5A6 = AI5 'DROIT'  A6 'DINI' DX1 'DFIN' DX1 ;
 LIGH = A5AI5 'ET' AI5A6 ;
 A6A14 = A6 'DROIT'  A14 'DINI' DX1 'DFIN' DX1 ;
*---------------------------
 DOM2 = 'DALLER' LIGBG AI3A5 LIGH A6A14 'PLAN' ;
*---------------------------------------------
*  The Channel (junction between rooms)
*---------------------------------------------
 A21 = A2 'PLUS' (DW 0.0) ;
 A31 = A3 'PLUS' (DW 0.0) ;
 A2A21 = A2 'DROIT' NY1 A21 ;
 A21A31 = A21 'DROIT' NY1 A31 ;
 A31A3 = A31 'DROIT'  NY1 A3 ;
 A3A2 = A3 'DROIT'  NY1 A2 ;
*-------------------------
 DOMC = 'DALLER' A2A21 A21A31 A31A3 A3A2 'PLAN' ;
*--------------------------------------------
* Third room (R1801)
*--------------------------------------------
*---- Length of the third room
 L3 = 5.0D0 ;
*--- length before the window 
 LB = 1.8D0 ;
*--- width of window
 LW = 0.668D0 ;
*--- length after the window
 LA = L3 '-' LB '-' LW ;  
*----------------------------------------
*  Low Domain
*----------------------------------------
 H3 = (9.5 '-' H2) '-' H1  ;
 MH3 = (-1.0D0) '*' H3 ;
 A70 = (L1 '+' DW) MH3 ;
 A71 = (L1 '+' LB '+' DW) MH3 ;
 A72 = (L1 '+' LB '+' DW) 0.0 ;
*----------------------------
 DW1 = LW '/' NY1 ;
 A70A71 = A70 'DROIT'  A71 'DINI' (DW '/' NY1) 'DFIN' DW1 ;
 A71A72 = A71 'DROIT' A72 'DINI' DX1 'DFIN' DX1 ;
 A72A21 = A72 'DROIT' A21 'DINI' DW1 'DFIN' (DW '/' NY1) ;
 A21A70 = A21 'DROIT' A70 'DINI' DX1 'DFIN' DX1 ;
*------------------------------
 DOM3A = 'DALLER' A70A71 A71A72 A72A21 A21A70 'PLAN' ;
*------------------------------
 A7W = A71 'PLUS' (LW 0.0) ;
 A7U = A72 'PLUS' (LW 0.0) ;
*--------------------------
 A71A7W = A71 'DROIT' NY1 A7W ;
 A7WA7U = A7W 'DROIT' A7U 'DINI' DX1 'DFIN' DX1 ;
 A7UA72 = A7U 'DROIT' NY1 A72 ;
 A72A71 = 'INVERSE' (A71A72) ;
*-------------------------
 DOM3B = 'DALLER' A71A7W A7WA7U A7UA72 A72A71 'PLAN' ;
*-------------------------
 A81 = A7W 'PLUS' (LA 0.0) ;
 A82 = A7U 'PLUS' (LA 0.0) ;
*-------------------------
 A7WA81 = A7W 'DROIT' A81 'DINI' DW1 'DFIN' DX1 ;
 A81A82 = A81 'DROIT' A82 'DINI' DX1 'DFIN' DX1 ;
 A82A7U = A82 'DROIT' A7U 'DINI' DX1 'DFIN' DW1 ;
 A7UA7W = 'INVERSE' (A7WA7U) ;
*-------------------------
 DOM3C = 'DALLER' A7WA81 A81A82 A82A7U A7UA7W 'PLAN' ;
*-------------------------
 DOM3 = DOM3A 'ET' DOM3B 'ET' DOM3C ;
 'ELIMINATION'  DOM3 1.0D-6 ;
*-----------------------------------------
*  Middle domain
*-----------------------------------------
 A73 = A72 'PLUS' (0.0 H1) ;
***--- first layer 
 A21A72 = A21 'DROIT' A72 'DINI' (DW '/' NY1) 'DFIN' DW1 ;
 A72A73 = A72 'DROIT' NY1 A73 ;
 A73A31 =  A73 'DROIT' A31 'DINI' DW1 'DFIN' (DW '/' NY1) ;
 A31A21 = A31 'DROIT' NY1 A21 ; 
*-----------------------------
 DOM4A = 'DALLER' A21A72 A72A73 A73A31 A31A21 'PLAN' ;
*-------------------------------------------
 A7M = A7U 'PLUS' (0.0 H1) ;
*---------------
 A72A7U = A72 'DROIT' NY1 A7U ;
 A7UA7M = A7U 'DROIT' NY1 A7M ;
 A7MA73 = A7M 'DROIT' NY1 A73 ;
 A73A72 = 'INVERSE' (A72A73) ;
*---------------------
 DOM4B = 'DALLER' A72A7U A7UA7M A7MA73 A73A72 'PLAN' ;
*---------------------------
 A83 = A82 'PLUS' (0.0 H1) ;
*--------------------------- 
 A7UA82 = A7U 'DROIT' A82 'DINI' DW1 'DFIN' DX1 ;
 A82A83 = A82 'DROIT' NY1 A83 ;
 A83A7M = A83 'DROIT' A7M 'DINI' DX1 'DFIN' DW1 ; 
 A7MA7U = 'INVERSE' (A7UA7M) ;
*---------------------------
 DOM4C = 'DALLER' A7UA82 A82A83 A83A7M A7MA7U 'PLAN' ;
*---------------------------
 DOM4 = DOM4A 'ET' DOM4B 'ET' DOM4C ;
 'ELIMINATION'  DOM4 1.0D-6 ;
*-----------------------------------------
*  High domain
*-----------------------------------------
 A32 = A31 'PLUS' (0.0 H2) ;
 A74 = A73 'PLUS' (0.0 H2) ; 
*---------------------------
 A31A73 = 'INVERSE' (A73A31) ;
 A73A74 = A73 'DROIT' A74 'DINI' DX1 'DFIN' DW1 ;
 A74A32 = A74 'DROIT' A32 'DINI' DW1 'DFIN' (DW '/' NY1) ;
 A32A31 = A32 'DROIT' A31 'DINI' DW1 'DFIN' DX1 ;
*---------------------------
 DOM5A = 'DALLER' A31A73 A73A74 A74A32 A32A31 'PLAN' ;
*---------------------------
 A7H = A7M 'PLUS' (0.0 H2) ;
*---------------------------
 A73A7M = 'INVERSE' (A7MA73) ;
 A7MA7H = A7M 'DROIT' A7H 'DINI' DX1 'DFIN' DW1 ;
 A7HA74 = A7H 'DROIT' NY1 A74 ;
 A74A73 = 'INVERSE' (A73A74) ;
*---------------------------
 DOM5B = 'DALLER' A73A7M A7MA7H A7HA74 A74A73 'PLAN' ;
*---------------------------
 A84 = A83 'PLUS' (0.0 H2) ;
*---------------------------
 A7MA83 = A7M 'DROIT' A83 'DINI' DW1 'DFIN' DX1 ;
 A83A84 = A83 'DROIT' A84 'DINI' DX1 'DFIN' DW1 ;
 A84A7H = A84 'DROIT' A7H 'DINI' DX1 'DFIN' DW1 ;
 A7HA7M = 'INVERSE' (A7MA7H) ; 
*-----------------------------
 DOM5C = 'DALLER' A7MA83 A83A84 A84A7H A7HA7M 'PLAN' ;    
*------------------------------------- 
 DOM5 = DOM5A 'ET' DOM5B 'ET' DOM5C ;
 'ELIMINATION'  DOM5 1.0D-6 ;
****----------*****----------------------******------- 
*           Upper vent
*---**********-----**********************------*******
 A75 = A74 'PLUS' (0.0 LW) ;
 A7F = A7H 'PLUS' (0.0 LW) ;
*---------------------------
 A74A7H = A74 'DROIT' NY1 A7H ;
 A7HA7F = A7H 'DROIT' NY1 A7F ;
 A7FA75 = A7F 'DROIT' NY1 A75 ;
 A75A74 = A75 'DROIT' NY1 A74 ;
*---------------------------
 DOMW = 'DALLER' A74A7H A7HA7F A7FA75 A75A74 'PLAN' ;

 FRONTH = A7FA75 ;
*----------------------------------*
*      Ignition Zone               *
*----------------------------------*
**  height at which we ignite
 HIGN = H2 '/' 2.0D0 ;
*-------------------------
 PP1 = 0.0 HIGN ;
 PP2 = PP1 'PLUS' (0.0 (DX1)) ;
 PP3 = PP2 'PLUS' (0.0 (DX1)) ;
 gg = (-1.0D0) '*' DX '*' 0.6 ;
 PP5 = PP1 'PLUS' (0.0 gg) ;
 PP51 = PP5 'PLUS' (0.0 gg) ;
 PO1 = A6A14 'POIN' 'PROCHE' PP1 ;
 PO2 = A6A14 'POIN' 'PROCHE' PP2 ;
 PO31 = A6A14 'POIN' 'PROCHE' PP3 ;
 PO5 = A6A14 'POIN' 'PROCHE' PP5 ;
 PO51 = A6A14 'POIN' 'PROCHE' PP51 ;
**--------------------------------
 PPP = DX1 0.0 ;
 PA1 = A1A12 'POIN' 'PROCHE' PPP ;
 PPP = (2.0 '*' DX) 0.0 ;
 PA2 = A1A12 'POIN' 'PROCHE' PPP ;
 XXX1 = 'COORDONNEE' 1 PA1 ;
 XXX2 = 'COORDONNEE' 1 PA2 ;
 YYY2 = 'COORDONNEE' 2 PO2 ;
 YYY3 = 'COORDONNEE' 2 PO31 ;
 PO4 = XXX1 YYY2 ;
 PO7 = XXX1 YYY3 ;
*-----------------------------
*  Just eight cells
*-----------------------------
 PO1PO2 = PO1 'DROIT' 1 PO2 ;
 PO2PO31 = PO2 'DROIT' 1 PO31 ;
 PO4PO7 = PO4 'DROIT' 1 PO7 ;
 PO5PO51 = PO5 'DROIT' 1 PO51 ;
**------------------------------
 'SI' (RAF 'EGA' 10) ;
    NZD = 1 ;
    KEF = -1.0D0 ;
 'SINON'  ;
    KEF = -2.0D0 ;
    NZD = 2 ;
 'FINSI'  ;      
 DIG1 = 'TRANSLATION' PO1PO2 NZD ((KEF '*' XXX1) 0.0) ;
 DIG5 = 'TRANSLATION' PO2PO31 NZD ((KEF '*' XXX1) 0.0) ;
 DIG6 = 'TRANSLATION' PO4PO7 1 ((XXX2-XXX1) 0.0) ;
 DIG7 = 'TRANSLATION' PO5PO51 1 (XXX1 0.0) ;
 'SI' (RAF '>' 10) ;
   YO31 = 'COORDONNEE' 2 PO31 ;
   YO2 = 'COORDONNEE' 2 PO1 ;
   DIFY = YO31 '-' YO2 ;
   VECUP = (0.0 DIFY) ;
   DIG6 = DIG1 'PLUS' VECUP ;
   DIG7 = DIG5 'PLUS' VECUP ;
 'FINSI'  ; 
*-----------------------------
 'SI' (RAF '>' 10) ; 
   DOMI = DIG1 'ET' DIG5 'ET' DIG6 'ET' DIG7  ;
 'SINON' ;
   DOMI = DIG1 'ET' DIG5 ;
 'FINSI' ; 
     
 'ELIMINATION'  DOMI 1.0D-6;
* 'TRACER' DOMI ;
*----------------------------------------------------
*                  Total Domain  
*---------------------------------------------------- 
 DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 'ET' DOM4 'ET' DOM5 'ET'
          DOMC 'ET' DOMW 'ET' DOMI ;
 'ELIMINATION' DOMTOT 1D-6 ;
*-------------------------------------
* Total Domain without ignition region
*------------------------------------- 
 DREST = 'DIFF' DOMTOT DOMI ;
*----------------------------------
* The fird room (R1801)
*---------------------------------- 
 DTHIRD = DOM3 'ET' DOM4 'ET' DOM5 'ET' DOMC 'ET' DOMW ;
 'ELIMINATION'  DTHIRD 1.0D-6 ;
*----------------------------------------
* First and second room (R1904 and R1905)
*----------------------------------------
 DSEC = DOM1 'ET' DOM2 ; 
*-------------------------------
* Lignes de postreatment
*-------------------------------
 LIGA = LIGBG 'PLUS' (0.0 (DX1 '/' 2.0)) ;
 A3A31 = 'INVERSE'  A31A3 ;
 A31A73 = 'INVERSE' A73A31 ;
 A73A7M = 'INVERSE' A7MA73 ;
 A7MA83 = 'INVERSE' A83A7M ;
 LIGT = A3A31 'ET' A31A73 'ET' A73A7M 'ET' A7MA83 ;
 'ELIMINATION'  LIGT 1.0D-6 ;
 dis1 = (DX1 '/' 2) '-' (DW '/' 2.0) ;  
 LIGB2 = LIGT 'PLUS' (0.0 dis1) ;
 LIGC = LIGA 'ET' LIGB2 ;
 'ELIMINATION' LIGC 1.0D-6; 
 LIG1 = LIGC ;

* 'TRACER' (DOMTOT 'ET' LIG1) ;
*-------------------------------
* Points for capturing variables
*-------------------------------
*--------- inside the first two rooms
 dish = (H2 '/' 2) '-' (DX1 '/' 2) ;
 disl = (H1 '/' 2) '+' (DX1 '/' 2) ;
 PCEL1 = (L1 '/' 6) dish ;
 PCEL2 = (L1 '/' 6) disl ;
 PCEL3 = (L1 '/' 2) dish ;
 PCEL4 = (L1 '/' 2) disl ;
 PCEL5 = (L1 '-' (L1 '/' 6)) dish ;
 PCEL6 = (L1 '-' (L1 '/' 6)) disl ;
*--------  inside the third room
 dis1 = L1 '+' (2.0 '*' DW) ;
 PCEL7 = dis1 (2.0 '*' DW) ;
 dis2 = (LB '+' (LW '/' 2) '+' (DX1 '/' 2.0)) '+' L1 '+' DW ;
 PCEL8 = dis2 (DW '-' 5.0) ;
 PCEL9 = dis2 (4.5 '-' DW) ;
  
* 'TRACER' (DOMTOT 'ET' PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
*    'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9) ;
*
* 'TRACER' (('CONTOUR' DREST)
*     'ET' PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
*    'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9) ;
*----------------------------------
* Creating the models
*----------------------------------
 MDNS = 'EULER' ;
*---------------------------------
 MDOMI = 'MODELISER' DOMI MDNS ;
 MDREST = 'MODELISER' DREST MDNS ;
 MFRONTH = 'MODELISER' FRONTH MDNS ;
 MDOMTOT = 'MODELISER' DOMTOT MDNS ;
 MDTHIRD = 'MODELISER'  DTHIRD MDNS ;
 MDSEC = 'MODELISER'  DSEC MDNS ; 
*----------------------------------
 $DOMI  = 'DOMA' MDOMI 'VF' ;
 $DREST = 'DOMA'  MDREST 'VF' ;
 $FRONTH = 'DOMA' MFRONTH 'VF' ;
 $DOMTOT  = 'DOMA' MDOMTOT 'VF' ;
 $DTHIRD = 'DOMA' MDTHIRD 'VF' ;
 $DSEC = 'DOMA' MDSEC 'VF' ;
*----------------------------------
 QDOMI  = $DOMI . 'QUAF' ;
 QDREST = $DREST . 'QUAF' ;
 QFRONTH = $FRONTH . 'QUAF' ;
 QDOMTOT  = $DOMTOT . 'QUAF' ;
 QDTHIRD = $DTHIRD . 'QUAF' ;
 QDSEC = $DSEC . 'QUAF' ;
*----------------------------------
 'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDOMI ;
 'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDREST ;
 'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QFRONTH ;
 'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDTHIRD ;
 'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDSEC ;
*-------------------------------
 LIN1 = LIGPROC LIG1 ('DOMA' MDOMTOT 'CENTRE') ;

* 'TRACER' (DOMTOT 'ET' LIN1) ;
*----------------------------------
*  Ligne at the "coridor"
*----------------------------------
  VECL = ((-0.5 '*' DX1) (0.5 '*' DX1)) ;
  PL2 = A21 'PLUS'  VECL ;
  PL2 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PL2 ;
  VECL = ((-0.5 '*' DX1) (-0.5 '*' DX1)) ;
  PL3 = A31 'PLUS'  VECL ;
  PL3 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PL3 ;
  LIG2 = PL2 'DROIT' (NY1 '-' 1)  PL3 ;
  'SI' (RAF 'EGA' 10) ;
    LIN2 = LIG2 ;
  'SINON' ;  
    LIN2 = LIGPROC LIG2 ('DOMA' MDOMTOT 'CENTRE') ;
  'FINSI' ;
*  'TRACER'  (DOMTOT 'ET' LIG2) ;
*  'TRACER'  (DOMTOT 'ET' LIN2) ;
*--------------------------------
  MDLIN2 = 'MODELISER' LIN2 MDNS ;
  $DLIN2  = 'DOMA' MDLIN2 'VF' ;
  QDLIN2  = $DLIN2 . 'QUAF' ;
  'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDLIN2 ;
*----------------------------------
* "Centralising the capteurs"
*---------------------------------- 
 PCE1 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL1 ;
 PCE2 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL2 ;
 PCE3 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL3 ;
 PCE4 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL4 ;
 PCE5 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL5 ;
 PCE6 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL6 ;
 PCE7 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL7 ;
 PCE8 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL8 ;
 PCE9 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL9 ;

* 'TRACER' (DOMTOT 'ET' PCE1 'ET' PCE2 'ET' PCE3 'ET' PCE4
*    'ET' PCE5 'ET' PCE6 'ET' PCE7 'ET' PCE8 'ET' PCE9) ;
* 'TRACER' (DOMTOT 'ET' ('DOMA' MDOMI 'CENTRE')) ;

  numbel = 'NBEL' DOMTOT ;
*-------------------------------------------
*  Cells sizes: min, max and averaged
*-------------------------------------------
  VOLC = 'DOMA' MDOMTOT 'VOLUME' ;
  SIZE = VOLC '**' 0.5 ;
*------- minimal  
  SMIN = 'MINIMUM'  SIZE ;
*------- maximal  
  SMAX = 'MAXIMUM'  SIZE ;
*------- averaged
  CHUN = 'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE')
             1 'SCAL' 1.0 'NATU' 'DISCRET' ;
  SIZEAV = 'XTY' SIZE CHUN ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
  SIZEAV = SIZEAV '/' numbel ;
*-------------------------------
* Saving the results
*-------------------------------
 
*'OPTION' 'SAUVER' ('CHAINE' '/serg/HDR12.3.2_dir/'
* 'mail' RAF '.sauv') ;
* 'SAUV' ;
*-----------------------------------------------
*  Initial conditions '+' gas properties
*-----------------------------------------------
*  'OPTION'  'ECHO' 1 'DIME' 2 'ELEM' 'QUA4'
*           'TRAC' 'X' ;
*
*  RAF = 10 ;
*
*  'OPTION'  'REST' ('CHAINE' '/serg/HDR12.3.2_dir/mail' RAF '.sauv') ;
*  'RESTITUER' ;         

 zero = 1.0D-9 ;
 GRAPH = FAUX ;
 
***************************************************************
***** PROCEDURE POUR CALCULER LES 'CP' ET 'CV'            *****
***************************************************************
'DEBPROC' CPCV TN*'FLOTTANT' PGAZ*'TABLE' ;
*
 CPP = 'TABLE' ;
 CVV = 'TABLE' ;

 TT = TN ;
 T2 = TT '*' TT ;
 T3 = T2 '*' TT ;
 T4 = T3 '*' TT ;
 
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE')) ;
    NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
    A0 = ('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) ;
    A1 = ('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '/' 2.0 ;
    A2 = ('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '/' 3.0 ;
    A3 = ('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '/' 4.0 ;
    A4 = ('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '/' 5.0 ;
    CVV . NOMCOM  = A0 '+' (A1*TT) '+' (A2 '*' T2) '+'
        (A3 '*' T3) '+' (A4 '*' T4) ;
    CPP . NOMCOM = (CVV . NOMCOM) '+' (PGAZ . NOMCOM . 'R') ;     
'FIN' BL1 ;

 NOMCOM = PGAZ . 'ESPNEULE' ;
 A0 = ('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) ;
 A1 = ('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '/' 2.0 ;
 A2 = ('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '/' 3.0 ;
 A3 = ('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '/' 4.0 ;
 A4 = ('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '/' 5.0 ;
 CVV . NOMCOM  = A0 '+' (A1*TT) '+' (A2 '*' T2) '+'
        (A3 '*' T3) '+' (A4 '*' T4) ;
 CPP . NOMCOM = (CVV . NOMCOM) '+' (PGAZ . NOMCOM . 'R') ;   
 

'FINPROC' CPP CVV ;
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
'DEBPROC' CONS RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT'
               YN*'CHPOINT' $DM*'MMODEL' PGAZ*'TABLE' ;
*

A0 =  'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A1 =  'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A2 =  'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A3 =  'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A4 =  'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
RTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
YTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;

'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE')) ;
    NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
    YCEL = 'EXCO'  NOMCOM YN ;
    RYCEL = 'NOMC' (YCEL * RN) NOMCOM 'NATU' 'DISCRET' ;
    YTOT = YTOT '+' YCEL ;
    RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') * YCEL) ;
    A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
    A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
    A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
    A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
    A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
    'SI' ('EGA' &BL1 1) ;
        RYN = 'KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL ;
    'SINON' ;
        RYN = RYN 'ET'
        ('KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL) ;
    'FINSI' ;
'FIN' BL1 ;

YCEL = ('KCHT' $DM 'SCAL' 'CENTRE' 1.0) '-' YTOT ;
NOMCOM = PGAZ . 'ESPNEULE' ;
A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL)  ;
RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') '*' YCEL);

TN = PN '/' (RN '*' RTOT) ;
T2 = TN '*' TN ;
T3 = T2 '*' TN ;
T4 = T3 '*' TN ;
T5 = T4 '*' TN ;


ETHER =  (A0 * TN) '+' ((A1 '/' 2.0) * T2) '+'
         ((A2 '/' 3.0) * T3) '+' ((A3 '/' 4.0) * T4) '+'
         ((A4 '/' 5.0) * T5) ; 
         
GN = RN '*' VN ('MOTS' 'SCAL' 'SCAL')  ('MOTS' 'UX  ' 'UY  ')
         ('MOTS' 'UX  '  'UY  ') ;

GN = 'KCHT' $DM 'VECT' 'CENTRE' GN ;         
         
ECIN = 0.5D0 '*' ('PSCAL' GN VN ('MOTS' 'UX  ' 'UY  ')
                 ('MOTS' 'UX  ' 'UY  ')) ;        

REN = (RN '*' ETHER) '+' ECIN ;

REN = 'KCHT' $DM 'SCAL' 'CENTRE' REN ;

'FINPROC' GN REN RYN ;
************************************************************
************************************************************
*                                                          * 
*****                Gas properties                    *****
*                                                          *
************************************************************
************************************************************
 PGAZ = 'TABLE' ;
 PGAZB = 'TABLE' ;
* Polynomial degree of specific heats
 PGAZ . 'NORD' = 4 ;
 PGAZB . 'NORD' = 0 ;
* Species explicitly treated in the Euler Equations
 PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;
 PGAZB . 'ESPEULE' =  PGAZ . 'ESPEULE' ;  
* Species non explicitly treated
 PGAZ . 'ESPNEULE' = 'N2  ';
 PGAZB . 'ESPNEULE' = PGAZ . 'ESPNEULE' ;
* Gas species properties
 PGAZ .  'H2  ' = 'TABLE'  ;
 PGAZ .  'H2O ' = 'TABLE'  ;
 PGAZ .  'N2  ' = 'TABLE'  ;
 PGAZ .  'O2  ' = 'TABLE'  ;
*----------------
 PGAZB .  'H2  ' = 'TABLE'  ;
 PGAZB .  'H2O ' = 'TABLE'  ;
 PGAZB .  'N2  ' = 'TABLE'  ;
 PGAZB .  'O2  ' = 'TABLE'  ; 
* R (J/Kg/K) and molar masses
 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 ;
*-------------------------------------
 PGAZB .  'H2  ' . 'R' = RGAS '/' mh2 ;
 PGAZB .  'H2O ' . 'R' = RGAS '/' mh2o ;
 PGAZB .  'N2  ' . 'R' = RGAS '/' mn2 ;
 PGAZB .  '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 ;
*-------------------------------------
 PGAZB .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAZB .  'H2O ' . 'H0K' = -1.395D7 ;
 PGAZB .  'N2  ' . 'H0K' = -2.953D5 ;
 PGAZB .  'O2  ' . 'H0K' = -2.634D5 ; 
*-----------------------------------
*** Names of the passive scalars
*-----------------------------------
 PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ' ;
*-----------------------------------------------
 PGAZB . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ' ;
****************************************************
****************************************************
*             Initial conditions                   *
****************************************************
****************************************************
*** inside the domain
 tg    = 337.0D0 ;
 pg    = 100000.0D0 ;
 xh2o = 0.26D0 ;
 xh2  = 0.1D0 ;
 xo2  = 0.1344D0 ;
 xn2  = 1.0D0  - (xh2  + xo2 + xh2O) ; 
 uxg   = 0.0     ;
 uyg   = 0.0     ;
*-------------------------
 td   =  337.0D0 ;
 uxd  = 0.0D0 ;
 uyd  = 0.0D0 ;
 pd   = 100000.0D0 ;
*---- values of K_0 
 K0IN = 0.41D0 ;
 K0IG = 0.41D0 ;
 K03 = 3.7D0 ;
*--------------------------------------------
* After burning: AIBC (same pressure)
*--------------------------------------------
 tc = 1086.38864D0 ;
 xh2c = zero ;
 xo2c = 0.0888421D0 ;
 xh2oc = 0.3789469D0 ;
 xn2c = 1.0D0 '-' (xh2c '+' xo2c '+' xh2oc) ;
 pc = pg ;
*---------------------------------------------
 CP1 CV1 = CPCV tg PGAZ ;
*'-'*'-'*'-'*'-'*'-'*- 
 PGAZB . 'CP' = 'TABLE' ;
 PGAZB . 'CP' . 'H2  '  = CP1 . 'H2' ;
 PGAZB . 'CP' . 'H2O '  = CP1 . 'H2O' ;
 PGAZB . 'CP' . 'O2  '  = CP1 . 'O2' ;
 PGAZB . 'CP' . 'N2  '  = CP1 . 'N2' ;
*----------------------------------------
 PGAZB . 'CV' = 'TABLE' ;
 PGAZB . 'CV' . 'H2  '  = CV1 . 'H2' ;
 PGAZB . 'CV' . 'H2O '  = CV1 . 'H2O' ;
 PGAZB . 'CV' . 'O2  '  = CV1 . 'O2' ;
 PGAZB . 'CV' . 'N2  '  = CV1 . 'N2' ;
*---------------------------------------------
* Creating CHAMPOINTS
*---------------------------------------------
*** Pressure
 MAILC = 'DOMA' MDOMTOT 'CENTRE' ;
 MAILCI = 'DOMA'  MDOMI 'CENTRE' ;
 MAILCR = 'DOMA'  MDREST 'CENTRE' ;
 MAILC2 = 'DOMA'  MDSEC 'CENTRE' ;
 MAILC3 = 'DOMA'  MDTHIRD 'CENTRE' ; 
*********************************  
 PINI = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' pc 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCR 1 'SCAL' pg 'NATU' 'DISCRET') ;
*** Temperature
 TINI = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' tc
         'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCR 1 'SCAL' tg
         'NATU' 'DISCRET') ;
**** Velocity
 WINI = 'MANUEL' 'CHPO' MAILC 2 'UX' uxg 'UY' uyg ;
*-------------------------------------------------------
* Combustion is complete, i.e. \ksi = 1 after combustion 
*-------------------------------------------------------
 CSIMAX = 'MANUEL' 'CHPO' MAILC 1 'SCAL' 1.0D0 ;
**** Molar fractions before combustion
 XH21 = ('MANUEL' 'CHPO' MAILCI  1 'H2' xh2c
         'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCR  1 'H2' xh2
         'NATU' 'DISCRET') ; 
 XH2O1 = ('MANUEL' 'CHPO' MAILCI 1 'H2O' xh2oc
         'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCR 1 'H2O' xh2o
         'NATU' 'DISCRET') ;
 XO21 =  ('MANUEL' 'CHPO' MAILCI 1 'O2' xo2c
         'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCR 1 'O2' xo2
         'NATU' 'DISCRET') ;
 XN21 =  ('MANUEL' 'CHPO' MAILCI 1 'N2' xn2c
         'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCR 1 'N2' xn2
         'NATU' 'DISCRET') ;

 XN = XH21 'ET' XH2O1 'ET' XO21 'ET' XN21 ;        
**** Mass Fractions
 MMOL = ('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'H2' mH2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'O2' mO2
         'NATURE' 'DISCRET') 'ET'   
        ('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'H2O' mH2O
         'NATURE' 'DISCRET') 'ET'    
        ('MANUEL' 'CHPO' ('DOMA' MDOMTOT '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 ;
  aaa = 'MAXIMUM'  (YH21) ;
*----------------------------------------------
*  Burned gas
* First we find how much of the burned gas
* produces using 1 mole of unburned gas
*----------------------------------------------
 CHUN = 'MANUEL' 'CHPO' MAILC 1 'SCAL' 1.0D0 'NATU' 'DISCRET' ; 
 H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
          ('EXCO' 'O2' XN) ;
 O2MAG = CHUN '-' H2MAG ;          
 
 XO22a = 'MANUEL' 'CHPO' MAILC 1 'O2' zero ;            
 XH22a = XH21 '-' ('NOMC' 'H2' (2.0 '*' (XO21 '-' XO22a))) ;

 XO22a = XO22a ;
 XH22a = XH22a ;

 XH22b = 'MANUEL' 'CHPO' MAILC 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 ;
*--------------------------------------------------- 
* Verifying that mtot2 = mtot
*---------------------------------------------------
 MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
         ('MOTS' 'H2' 'O2' 'H2O' 'N2')
         ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;

 'SI' (('MAXIMUM' (mtot2 '-' mtot) 'ABS') '>' zero) ;
     'MESSAGE' 'problem1' ;          
 'FINSI' ;              
***** xtot2 = mole of burned gas per mole of unburned gas
 xtot2 = 'PSCAL' (XH22 '+' XO22 '+'  XH2O2 '+' XN22)
         CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
         ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL') ;
***** X now is the real molar fractions 
 XH22 = XH22 '/' xtot2 ;
 XO22 = XO22 '/' xtot2 ; 
 XH2O2 = XH2O2 '/' xtot2 ;
 XN22 = XN22 '/' xtot2 ;

 xtot3 = 'PSCAL' (XH22 '+' XO22 '+' XH2O2 '+' XN22)
          CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
         ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL') ;

 'SI' (('MAXIMUM' (xtot3 '-' 1) 'ABS') '>' zero) ;
     'MESSAGE' 'problem2' ;          
 '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 complete combustion
 YH22 = (XH22 '*' MH2) '/' MTOT2 ;
 YO22 = (XO22 '*' MO2) '/' MTOT2 ;
 YN22 = (XN22 '*' MN2) '/' MTOT2 ;
 YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
*---------------------------------------------
*  Conservative variables
*---------------------------------------------
  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')) ;

  RN = PINI '/' (RINI '*' TINI) ;       
  YN = YH21 '+' YO21 '+' YH2O1 ;
  
  GN RETN RYN = CONS RN WINI PINI YN MDOMTOT PGAZ ;
*-------------------------------------- 
*           Passive scalars
*--------------------------------------
 YININ = 'MANUEL' 'CHPO' MAILC 1 'H2' aaa 'NATU' 'DISCRET' ; 
 YFINN = 'MANUEL' 'CHPO' MAILC 1 'H2' zero 'NATU' 'DISCRET' ; 
 K0N = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' K0IG
         'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILC2 1 'SCAL' K0IN 
         'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILC3 1 'SCAL' K03 
         'NATU' 'DISCRET') ; 
*--------------------------------------- 
 RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ; 
*-------------------
***** Verification 
*-------------------
 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') '>' (zero '*' 10.0)) ;
   'MESSAGE' 'problem3' ;
   'LISTE' ('MAXIMUM' ((P '-' PINI) '/' PCEL) 'ABS') ; 
 'FINSI'  ;
 'SI'(('MAXIMUM' ((T '-' TINI) '/' TCEL) 'ABS') '>' (zero '*' 10.0)) ;
   'MESSAGE' 'problem4' ;
   'LISTE' ('MAXIMUM' ((T '-' TINI) '/' TCEL) 'ABS') ; 
 'FINSI'  ;
 'SI' (('MAXIMUM'  VN 'ABS') '>' zero) ;
   'MESSAGE' 'problem5' ;
   'LISTE' ('MAXIMUM'  VN 'ABS') ; 
 'FINSI'  ;
 'SI'(('MAXIMUM' (Y '-' YN) 'ABS') '>' zero) ;
   'MESSAGE' 'problem6' ;
   'LISTE' ('MAXIMUM' (Y '-' YN) 'ABS') ;  
 'FINSI'  ;
*--------------------------------------
* Graphics
*--------------------------------------               
  'SI' GRAPH ;
    CHM_TN   =  'KCHA' MDOMTOT 'CHAM' TINI ;
    CHM_RN   =  'KCHA' MDOMTOT 'CHAM' RN ;
    CHM_PN   =  'KCHA' MDOMTOT 'CHAM' PINI ;
    CHM_VN   =  'KCHA' MDOMTOT 'CHAM' WINI ;
    CHM_H2   =  'KCHA' MDOMTOT 'CHAM' YN ;
    CHM_YIN  =  'KCHA' MDOMTOT 'CHAM' YININ ;
    CHM_YFI  =  'KCHA' MDOMTOT 'CHAM' YFINN ;
    CHM_K0   =  'KCHA' MDOMTOT 'CHAM' K0N ;
    'TRACER'  CHM_TN  MDOMTOT  
        'TITR'  ('CHAINE'  'Temperature at t='  0.0) ;      
    'TRACER'  CHM_RN  MDOMTOT  
        'TITR'  ('CHAINE'  'RN at t='  0.0) ;
    'TRACER'  CHM_PN  MDOMTOT 
        'TITR'  ('CHAINE' 'PN at t='  0.0) ;
    'TRACER'  CHM_VN MDOMTOT
        'TITR'  ('CHAINE' 'VN at t='  0.0) ;
    'TRACER'  CHM_H2 MDOMTOT
        'TITR'  ('CHAINE' 'H2 at t='  0.0) ;
    'TRACER'  CHM_YIN MDOMTOT
        'TITR'  ('CHAINE' 'YIN at t='  0.0) ;
    'TRACER'  CHM_YFI MDOMTOT
        'TITR'  ('CHAINE' 'YFI at t='  0.0) ;
    'TRACER'  CHM_K0 MDOMTOT
        'TITR'  ('CHAINE' 'K0 at t='  0.0) ;                    
 'FINSI' ; 


*-----------------------------------------------
* Parameters for the computation
*-----------------------------------------------
*** Upwind scheme
 METO = 'SS' ;
**** Iterations 
 NITER = 100 ;
**** Final time (no restriction)   
 TFINAL = 10000.1 ;
**** Safety factor for the time step
 SAFFAC = 0.7 ;
*** Second order ?
 LOGSO = VRAI ;
 LOGSO = FAUX ;
*** EPS_CC '-' epsilon in CREBCOM criterion
 EPS_CC = 0.5 ;
**** "explosive" initial conditions?
 LOGEXP = VRAI ;
 LOGEXP = FAUX ;
*** Are there any Boundary conditions
 LOGBC = VRAI ;
*-------------------------------
*  Time intervals for the output
*-------------------------------
 LISTT = 'PROG' 4.5 ;
*** Mesh size
 DELTAXN = ('DOMA'  MDOMTOT 'VOLUME') '**' 0.5 ;
 DELX = 'MAXIMUM' DELTAXN ;
*-------------------------------------
 TABLIM = 'TABLE' ;
 pwin = (('MAXIMUM'  RN) '*' 9.8) '*' 4.668 ;
 TABLIM . 'CHPOUT' =
        'MANUEL' 'CHPO' ('DOMA' MFRONTH 'CENTRE') 1 'PN' (pg '-' pwin) ;
 TABLIM . 'H' = 50.0D0 ;
 TABLIM . 'TNUL' =
        'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'SCAL' tg ;
 TABLIM . 'FREQ' = 1 ;
*----- gravity ----------
 TABLIM . 'GRAV' = 'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE')
       'UX' 0.0D0 'UY' (-9.8D0) ; 
*----------------------------------------------
* Table for evolution of the variables with time
* at the fixed points
*---------------------------------------------- 
 TABC = 'TABLE' 'CAPTEUR' ;
 TABC . 'LPOINTS' = PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
  'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 ;

 TABC . 'LTPS'    = 'PROG';
 TABC . 'NOMC'    = 'MOTS' 'C1' 'C2' 'C3' 'C4' 'C5' 'C6'
                     'C7' 'C8' 'C9' ;

 TABC . 'COMP'    = 'MOTS' 'PN' 'TN' ;
 TABC . 'FREQ'    = 10 ;
 MSLIN2 = 'DOMA'  MDLIN2 'SOMMET' ;
 bnb = 'NBEL' MSLIN2 ;
 VOLGOR = DX1 '*' DX1 '*' bnb ;
*-----------------------------------------------
* In MDOM1 we 'burn' the gas;
* In order to do so we shall apply the
* operator 'FLAM' in MDOM1 with eps=0
* and DELTAT >> DELTATC
*----------------------------------------------
 YN = RYN '/' RN ;
 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' ;
**** Stocheom. coefficients
 LCOEF = 'PROG' 1.0 0.5 -1.0 ; 
*-------------------------------------------
* This to activate in case of constant volume
*    combustion
*-------------------------------------------
 'SI' LOGEXP ; 
   DELTARE DELTARY = 'FLAM' 'CREBCOM2' MDOMI PGAZ LMOT1 LCOEF 
    ('REDU' RN MAILCI) ('REDU' YN MAILCI)
    ('REDU' YININ MAILCI) ('REDU' YFINN MAILCI)
    ('REDU' K0N MAILCI) ('REDU' DELTAXN MAILCI)
      0.0 (1D10 '*' DELTATC)  zero ;

   RYN = RYN '+' DELTARY ;
   YN = RYN '/' RN ;
   RETN = RETN '+' DELTARE ;
 'FINSI'  ;  
*-----------------------------------------------
*  Names of the components    
*----------------------------------------------- 
 LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2'
                   'H2O' 'H2IN' 'H2FI' 'K0' ;
****** names of the primitive variables
 LISTINCP = 'MOTS' 'RN' 'UX' 'UY' 'PN' 'H2' 'O2'
                   'H2O' 'H2IN' 'H2FI' 'K0' ;
*--------------------------------------------                                      
 LMOT111 = 'MOTS' 'SCAL' ;
 GRADR CCC COEFSCAL = 'PENT' MDOMTOT 'CENTRE'
 'EULESCAL' 'NOLIMITE' LMOT111  RN ;

 LMOT2 = 'MOTS' 'UX' 'UY' ;
 GRADV CCC COEFVECT = 'PENT' MDOMTOT 'CENTRE'
 'EULEVECT' 'NOLIMITE' LMOT2 GN ;
*--------------------------------------------
*   Volumes of the rooms ;
*--------------------------------------------
 VV2 = 'DOMA' MDSEC 'VOLUME' ;
 CHUN = ('MANUEL' 'CHPO' MAILC2 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
 VV2S = ('XTY' CHUN VV2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
*-------------------------------
 VV3 = 'DOMA' MDTHIRD 'VOLUME' ;
 CHUN = ('MANUEL' 'CHPO' MAILC3 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
 VV3S = ('XTY' CHUN VV3 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
*-------------------------------------------
 LISTPS2 = 'PROG' ;
 LISTPS3 = 'PROG' ;
 LISTTS2 = 'PROG' ;
 LISTTS2 = 'PROG' ;
*---- 'LISTE' for the discharge at the coridor
 LISCOR1 = 'PROG' ;
*---- 'LISTE' for the averaged velocity at the coridor
 LISCOR2 = 'PROG' ;
*---- 'LISTE' for the discharge at the vent
 LISVEN1 = 'PROG' ;
*---- 'LISTE' for the averaged velocity at the vent
 LISVEN2 = 'PROG' ;  
*--------------------------------------------
 TPS = 0.0 ;
*------------------------------------------
*        Chemical time step
*------------------------------------------
 DT_CHEM = SAFFAC '*' DELTATC ; 
*-------------------------------------------
**** Temporal loop
*-------------------------------------------
 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
 'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
 'MESSAGE' ;
*-----------------------------------
 'REPETER' BLG ('DIME' LISTT) ;
*----------------------------------
  TFINAL = 'EXTRAIRE' &BLG LISTT ; 
*---------------------------------- 
  VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;

  TNM1 = 'COPIER' TN ;

 'TEMPS' 'ZERO' ; 
 'REPETER' BL1 NITER ;
*
**** Primitive variables
*
 FREQC = TABC.'FREQ' ;
 'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
******* The "coridor" **********************
   RL2 = 'REDU' RN MSLIN2 ;
   GL2 = 'REDU' GN MSLIN2 ;
   GLX2 = 'EXCO' 'UX' GL2 ;
   CHUN = ('MANUEL' 'CHPO' MSLIN2 1
             'SCAL' (DX1) 'NATU' 'DISCRET') ;
*------- rashod discharge ----------------
   GLXAV = ('XTY' CHUN GLX2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
   GLXAV = GLXAV '/' H1 ;
   LISCOR1 = LISCOR1 'ET' ('PROG' GLXAV) ;
*---- redifining CHUN -------------------
   CHUN = ('MANUEL' 'CHPO' MSLIN2 1
             'SCAL' (DX1 '*' DX1) 'NATU' 'DISCRET') ;           
*------- averaged density ----------------             
   RLAV = ('XTY' CHUN RL2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;              
   RLAV = RLAV '/' VOLGOR ;
*------- averaged momentum ---------------
   GLXAV = ('XTY' CHUN GLX2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;              
   GLXAV = GLXAV '/' VOLGOR ;
*-----------------------------------------   
*   UXL = GLXAV '/' RLAV ;
   CHVIT = GLX2 '*' ('INVERSE' RL2) ;
   UXL = 'MAXIMUM' CHVIT ;
   LISCOR2 = LISCOR2 'ET' ('PROG' UXL) ;
******* The second room **********************   
*----------------------------------------------------
*   Computation of averaged pressure and temperature
*----------------------------------------------------
  RN2 = 'REDU' RN MAILC2 ;
  GN2 = 'REDU' GN MAILC2 ;
  RETN2 = 'REDU' RETN MAILC2 ;
  RYN2 = 'REDU' RYN MAILC2 ;
  RSN2 = 'REDU' RSN MAILC2 ;
*--- we average the conservative data
  RN2S = ('XTY' RN2 VV2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  RN2S = RN2S '/' VV2S ;
*----------------------
  RETN2S = ('XTY' RETN2 VV2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  RETN2S = RETN2S '/' VV2S ;
*-----------------------
  GNX2 = 'EXCO' 'UX' GN2 'UX' ;   
  GNY2 = 'EXCO' 'UY' GN2 'UY' ;
*-------------------  
  GNX2S = ('XTY' GNX2 VV2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  GNX2S = GNX2S '/' VV2S ;
*-------------------
  GNY2S = ('XTY' GNY2 VV2 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  GNY2S = GNY2S '/' VV2S ;
*--------------------
  RYNH2 = 'EXCO' 'H2' RYN2 'H2' ;   
  RYNO2 = 'EXCO' 'O2' RYN2 'O2' ;
  RYNH2O = 'EXCO' 'H2O' RYN2 'H2O' ;  
*-------------------  
  RYH2S = ('XTY' RYNH2 VV2 ('MOTS' 'H2')
       ('MOTS' 'SCAL')) ;
  RYH2S = RYH2S '/' VV2S ;
*-------------------
  RYO2S = ('XTY' RYNO2 VV2 ('MOTS' 'O2')
       ('MOTS' 'SCAL')) ;
  RYO2S = RYO2S '/' VV2S ;
*-------------------
  RYHO2S = ('XTY' RYNH2O VV2 ('MOTS' 'H2O')
       ('MOTS' 'SCAL')) ;
  RYHO2S = RYHO2S '/' VV2S ;  
*--------------------------------------
*--------------------------------------
  R2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RN2S 'NATU' 'DISCRET') ;
  G2S = ('MANUEL' 'CHPO' MAILCI 1 'UX' GNX2S 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'UY' GNY2S 'NATU' 'DISCRET') ;
  RET2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RETN2S 'NATU' 'DISCRET') ;      
  RY2S = ('MANUEL' 'CHPO' MAILCI 1 'H2' RYH2S 'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCI 1 'O2' RYO2S 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'H2O' RYHO2S 'NATU' 'DISCRET') ;
        
  RS2S = ('MANUEL' 'CHPO' MAILCI 1 'H2IN' (aaa) 'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCI 1 'H2FI' zero 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'K0' K0IN 'NATU' 'DISCRET') ;
*--------------------------------------------
  VS PS TS YS SS GAMS = 'PRIM' 'PERFTEMP' PGAZ R2S G2S RET2S RY2S RS2S ;
  PN2S = 'MAXIMUM'  PS ;
  TN2S = 'MAXIMUM'  TS ;
*-------------------------------
   LISTPS2 = LISTPS2 'ET' ('PROG' PN2S) ;
   LISTTS2 = LISTTS2 'ET' ('PROG' TN2S) ;  
*-------------------------------------------------
******** IN THE BIGGEST ROOM  ********************
  RN2 = 'REDU' RN MAILC3 ;
  GN2 = 'REDU' GN MAILC3 ;
 RETN2 = 'REDU' RETN MAILC3 ;
 RYN2 = 'REDU' RYN MAILC3 ;
 RSN2 = 'REDU' RSN MAILC3 ;
*--- we average the conservative data
  RN2S = ('XTY' RN2 VV3 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  RN2S = RN2S '/' VV3S ;
*----------------------
  RETN2S = ('XTY' RETN2 VV3 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  RETN2S = RETN2S '/' VV3S ;
*-----------------------
  GNX2 = 'EXCO' 'UX' GN2 'UX' ;   
  GNY2 = 'EXCO' 'UY' GN2 'UY' ;
*-------------------  
  GNX2S = ('XTY' GNX2 VV3 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  GNX2S = GNX2S '/' VV3S ;
*-------------------
  GNY2S = ('XTY' GNY2 VV3 ('MOTS' 'SCAL')
       ('MOTS' 'SCAL')) ;
  GNY2S = GNY2S '/' VV3S ;
*--------------------
  RYNH2 = 'EXCO' 'H2' RYN2 'H2' ;   
  RYNO2 = 'EXCO' 'O2' RYN2 'O2' ;
  RYNH2O = 'EXCO' 'H2O' RYN2 'H2O' ;  
*-------------------  
  RYH2S = ('XTY' RYNH2 VV3 ('MOTS' 'H2')
       ('MOTS' 'SCAL')) ;
  RYH2S = RYH2S '/' VV3S ;
*-------------------
  RYO2S = ('XTY' RYNO2 VV3 ('MOTS' 'O2')
       ('MOTS' 'SCAL')) ;
  RYO2S = RYO2S '/' VV3S ;
*-------------------
  RYHO2S = ('XTY' RYNH2O VV3 ('MOTS' 'H2O')
       ('MOTS' 'SCAL')) ;
  RYHO2S = RYHO2S '/' VV3S ;  
*--------------------------------------
*--------------------------------------
  R2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RN2S 'NATU' 'DISCRET') ;
  G2S = ('MANUEL' 'CHPO' MAILCI 1 'UX' GNX2S 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'UY' GNY2S 'NATU' 'DISCRET') ;
  RET2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RETN2S 'NATU' 'DISCRET') ;      
  RY2S = ('MANUEL' 'CHPO' MAILCI 1 'H2' RYH2S 'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCI 1 'O2' RYO2S 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'H2O' RYHO2S 'NATU' 'DISCRET') ;
        
  RS2S = ('MANUEL' 'CHPO' MAILCI 1 'H2IN' (aaa) 'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' MAILCI 1 'H2FI' zero 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' MAILCI 1 'K0' K0IN 'NATU' 'DISCRET') ;
*--------------------------------------------
  VS PS TS YS SS GAMS = 'PRIM' 'PERFTEMP' PGAZ R2S G2S RET2S RY2S RS2S ;
  PN2S = 'MAXIMUM'  PS ;
  TN2S = 'MAXIMUM'  TS ;          
*-------------------------------------------------
    LISTPS3 = LISTPS3 'ET' ('PROG' PN2S) ;
    LISTTS3 = LISTTS3 'ET' ('PROG' TN2S) ;
 'FINSI' ;
*--------------------------------------------------
*-------------------------------------------
    VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN 
      TNM1  ;
*--------------------------------------------
*   Boundary conditions
*--------------------------------------------
*** find new CP and CV *******
*    CP1 CV1 = CPCV tg PGAZ ;
*------------------------------------------
  'SI' LOGBC ;
   'SI' (&BLG 'EGA' 1) ;    
     SHIN = 'EXCO' SN 'H2IN' ;
     SHFI = 'EXCO' SN 'H2FI' ;
     SHK0 = 'EXCO' SN 'K0' ;
     RCHLIM RCHRES = 'KONV' 'VF' 'PERFMULT' 'CLIM'  'RESI'
        MDOMTOT MFRONTH PGAZB LISTINCO LISTINCP
        RN VN PN YN SHIN SHFI SHK0
        (TABLIM . 'CHPOUT') 'OUTP' ;
   'FINSI' ;     
*------------------------------------
    MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
    MAIL111 = 'EXTRAIRE' RCHRES 'MAILLAGE' ;
*------------------------------------
    NBB = 'NBEL'  MAIL111 ;
    TTT = 'REDU'  TN MAIL111 ;
    UNIF = 'MANUEL' CHPO ('DOMA' MDOMTOT 'CENTRE')
       1 'SCAL' 1.0 ;
    UNI1 = 'REDU' UNIF MAIL111 ;
    NUNI = UNI1 '*' (1.0 '/' NBB) ;
**** computing average temperature at the boundary    
    TAV = 'XTY' TTT NUNI ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
*** find new CP and CV *******
    CPN CVN = CPCV TAV PGAZ ;
*** Updating CP and CV for the boundary conditions    
     PGAZB . 'CP' . 'H2  '  = CPN . 'H2' ;
     PGAZB . 'CP' . 'H2O '  = CPN . 'H2O' ;
     PGAZB . 'CP' . 'O2  '  = CPN . 'O2' ;
     PGAZB . 'CP' . 'N2  '  = CPN . 'N2' ;
*----------------------------------------
     PGAZB . 'CV' . 'H2  '  = CVN . 'H2' ;
     PGAZB . 'CV' . 'H2O '  = CVN . 'H2O' ;
     PGAZB . 'CV' . 'O2  '  = CVN . 'O2' ;
     PGAZB . 'CV' . 'N2  '  = CVN . 'N2' ;
*  'FINSI' ;   
*--------------------------------
    SHIN = 'EXCO' SN 'H2IN' ;
    SHFI = 'EXCO' SN 'H2FI' ;
    SHK0 = 'EXCO' SN 'K0' ;
    RCHLIM RCHRES = 'KONV' 'VF' 'PERFMULT' 'CLIM'  'RESI'
        MDOMTOT MFRONTH PGAZB LISTINCO LISTINCP
        RN VN PN YN SHIN SHFI SHK0
        (TABLIM . 'CHPOUT') 'OUTP' ;
    MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
    MAIL111 = 'EXTRAIRE' RCHRES 'MAILLAGE' ;
*--------------------------------------------
*  We compute the debit et vitesse
*--------------------------------------------
  'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
   VOLUP = 'REDU' ('DOMA' MDOMTOT 'VOLUME')
       MAIL111 ;
   RLIM1 = 'EXCO' 'RN' RCHRES ;
   DEBUP = 'XTY' VOLUP RLIM1 ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
   DEBUP = (DEBUP '/' LW) '*' (-1.0) ;
   LISVEN1 = LISVEN1 'ET' ('PROG' DEBUP) ;
*--------------------------------------------
   RNUP = 'REDU' RN MAIL111 ;
   GNUP = 'REDU' GN MAIL111 ;
   GXUP = 'EXCO' 'UY' GNUP ;
   CUN = ('MANUEL' 'CHPO' MAIL111 1
             'SCAL' 1.0D0 'NATU' 'DISCRET') ;
*------------------------
   GXUPS = 'XTY' GXUP VOLUP ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
   RNUPS = 'XTY' RNUP VOLUP ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
   VITUP = GXUPS '/' RNUPS ;
   LISVEN2 = LISVEN2 'ET' ('PROG' VITUP) ;
  'FINSI' ; 
*--------------------------------------------    
 'FINSI'  ; 
*------------------------------------    
*--------------------------------
    TNM1 = 'COPIER' TN ;
    'SI' LOGSO ;    
       GRADR ALR = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN  'GRADGEO' COEFSCAL ;

       GRADP ALP = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN  'GRADGEO' COEFSCAL ;

       GRADV ALV = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY') VN   'GRADGEO'  COEFVECT ;
          
       GRADY ALY = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'H2' 'O2' 'H2O') YN  'GRADGEO' COEFSCAL ;

       GRADS ALS = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'H2IN' 'H2FI' 'K0') SN  'GRADGEO' COEFSCAL ;      

       ROF VITF PF YF SF = 'PRET' 'PERFTEMP'  2 1
                            MDOMTOT  PGAZ
                            RN    GRADR  ALR
                            VN    GRADV  ALV
                            PN    GRADP  ALP
                            YN    GRADY  ALY 
                            SN    GRADS  ALS ;
    'SINON' ;
                              
       ROF VITF PF YF SF  = 'PRET' 'PERFTEMP'  1 1 
                              MDOMTOT  PGAZ
                              RN   
                              VN   
                              PN   
                              YN    
                              SN ; 
    'FINSI' ;
*-----------------------------------------------    
  'SI' LOGBC ;
    RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
         MDOMTOT PGAZ LISTINCO  ROF VITF PF YF SF MAILLIM ;
    RESIDU  = RESIDU 'ET' RCHRES ;
  'SINON' ;
     RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
         MDOMTOT PGAZ LISTINCO  ROF VITF PF YF SF ;
  'FINSI' ;       
*-----------------------------------------------
    DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
    DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;    
*
**** Total time step
*
    DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM ) ;
*
**** Increment of the variables (convection)
*
    RESIDU = DTMIN '*' RESIDU ;
    
    DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
    DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
    DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
    DRYN  = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
            ('MOTS' 'H2' 'O2' 'H2O') ;
    DRSN  = 'EXCO' ('MOTS' 'H2IN' 'H2FI' 'K0') RESIDU
             ('MOTS' 'H2IN' 'H2FI' 'K0') ;        
    
    
    TPS  = TPS '+' DTMIN ;

    RN   = RN '+' DRN ;
    GN   = GN '+' DGN ;
    RETN = RETN '+' DRETN ;
    RYN = RYN '+' DRYN ;
    RSN = RSN '+' DRSN ;
    YN = RYN '/' RN ;
    SN = RSN '/' RN ;
*
**** Increment of the variables (chemical source)
*
    K0N = 'EXCO'  'K0' SN 'SCAL' ;
    YININ = 'EXCO' 'H2IN' SN 'H2' ;
    YFINN = 'EXCO' 'H2FI' SN 'H2' ; 

    DRETN DRYN = 'FLAM' 'CREBCOM2' MDOMTOT PGAZ LMOT1 LCOEF RN YN
     YININ YFINN K0N DELTAXN EPS_CC DTMIN (zero '*' 1.0D3) ;

    RYN = RYN '+' DRYN ;    
    RETN = RETN '+' DRETN ;
*---------------------------------------------------
*  Rate of loosing energy due to boundary conditions 
*---------------------------------------------------
*    'SI' (TPS '>' 3.0) ;   
    VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
    HH = TABLIM . 'H' ;
    PERT = HH '*' (TN '-' (TABLIM . 'TNUL')) ;
    RETN = RETN '-' (DTMIN '*' PERT) ;
*    'FINSI' ;
*--------------------------------------------------
*   Gravity forces
*--------------------------------------------------
    LIMOT1 = 'MOTS' 'SCAL' 'SCAL' ;
    LIMOT2 = 'MOTS' 'UX' 'UY' ;
    LIMOT3 = 'MOTS' 'UX' 'UY' ;
    ROG = RN '*' (TABLIM . 'GRAV') LIMOT1 LIMOT2 LIMOT3 ;
    GN = GN '+' (DTMIN '*' ROG) ;
    ROGU = 'PSCAL'  ROG VN LIMOT3 LIMOT2 ;
    RETN = RETN '+' (DTMIN '*' ROGU) ;
*---------------------------------------------------
*---------------------------------------------
*  Capturing variables at the capteurs
*---------------------------------------------
 'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
   CENTOT = 'DOMA' MDOMTOT 'CENTRE';
   PTS = TABC . 'LPOINTS'   ;
***** M1 is the number of capteurs   
   M1  = 'NBEL' PTS ;

   LTPS = TABC . 'LTPS' ;
   NTPS = 'DIME' LTPS     ;
**** N1 is the number of variables captured
   L1 = TABC . 'COMP' ;
   N1 = 'DIME' L1 ;
* On etudie uniquement la pression
* peut evoluer par la suite !!!
   LOG1 = 'EXISTE' TABC 'NOMC' ;
*   Initialisation des tables
  'SI' (NTPS EGA 0 );
   'REPETER' BOU1 M1;
**** extract the names of the capteurs   
      PT1  = PTS 'POIN' &BOU1;
      'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
         CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
      'SINON' ;
         CAPI = 'CHAINE' 'CAPT' &BOU1;
      'FINSI';
      PT2 = CENTOT 'POIN' 'PROCHE' PT1;
      'SI' (&BOU1 'EGA' 1);
         TABC.'RPOINTS' = PT2 ;
      'SINON';
         TABC.'RPOINTS' = TABC .'RPOINTS' 'ET' PT2;
      'FINSI';
      TABC . CAPI = 'TABLE';
        'REPETER' BOU2 N1;
           MOI = EXTR L1 &BOU2 ;
           TABC . CAPI . MOI = 'PROG';
        'FIN' BOU2;
   'FIN' BOU1;
 'FINSI';
*--------------------------------------------
  'REPETER' BOU1 M1 ;
     PT1  = TABC.'RPOINTS' 'POIN' &BOU1;
    'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
       CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
    'SINON' ;
       CAPI = 'CHAINE' 'CAPT' &BOU1;
    'FINSI';

   'REPETER' BOU2 N1;
     MOI = EXTR L1 &BOU2 ;
*** captage de la pression
    'SI' ('EGA' MOI 'PN');
        X1 = 'EXTR' PN 'SCAL' PT1 ;
*        'LISTE' X1 ;
        TABC . CAPI . MOI = (TABC . CAPI . MOI) 'ET' ('PROG' X1);
    'FINSI';
*** captage de temperature
    'SI' ('EGA' MOI 'TN');
        X1 = 'EXTR' TN 'SCAL' PT1 ;
        TABC . CAPI . MOI = TABC . CAPI . MOI 'ET' ('PROG' X1);
    'FINSI';
  'FIN' BOU2;

 'FIN' BOU1;
 'FINSI' ;
*****---- creating listreel of time steps 
  'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
      TABC . 'LTPS' = (TABC . 'LTPS') 'ET' ('PROG' TPS);
  'FINSI' ;    
*---------------------------------------------    
    'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
        'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
    'FINSI' ;
*-------------------------------------------------------
   'SI' (TPS '>EG' TFINAL) ;
       'QUITTER' BL1 ;
    'FINSI' ;

    
 'FIN' BL1 ;
*------------------------------------
*-----------------------------------   
 'FIN'  BLG ;
 'TEMPS' ;
*************************************************
*                  Posttreatment                *
*************************************************
    LISTP0 = 'PROG' 10 '*' 100000.0D0 ;
    OVPRES = (TABC . 'C1' . 'PN')  '-' LISTP0 ;
    EVERP1 = 'EVOL' 'MANU' 'Niter' (TABC . 'LTPS' ) 'PC1'
         OVPRES ;
*------------------------------------------------
    LPGOOD = 'PROG' 1.3086  21.753  43.070  41.450  30.604  26.187
               32.615  43.088   48.922   47.665 ;

    DIFFER = OVPRES '-' LPGOOD ;

    'SI' (('MAXIMUM' DIFFER) <EG 1.D-3) ;
      'MESSAGE' ;
      'MESSAGE' 'OVERPRESSURE IS OK' ;
      'MESSAGE' ;
    'SINON' ;
      'MESSAGE' ;
      'MESSAGE' 'THERE IS A PROBLEM' ;
      'MESSAGE' ;  
    'FINSI' ;               
 
  'FIN' ;

 
 
 


 

 

