* fichier :  shock3d.dgibi
************************************************************************
* NAME        : main3d.dgibi
* DESCRIPTION : file to solve 3D tests:
*               2. shock tube 
*
* LANGUAGE    : GIBIANE-CAST3M
* AUTHOR      : José R. Garcia Cascales 
*               Universidad Politecnica de Cartagena,
*               jr.garcia@upct.es
************************************************************************
* VERSION    : v1, 04/04/2002, version initiale
* HISTORIQUE : vfinal, 22/11/2005
************************************************************************
*
* MESH DEFINITION
*
************************************************************************  
'OPTION'  'DIME' 3 'ELEM' 'CUB8' 'ISOV' 'SULI'
    'ECHO' 1 'TRAC' 'X';

 GRAPH  = VRAI ;
 GRAPH  = FAUX ;
*
** number of elements in the z axis (after rotation) = NZ
*
 NZ = 100 ;

 NX = 2;
 NY = 2;

 H  = 10.D0  ;
 DZ = H/NZ ;

 L  =  NX*DZ ;
 P  =  NY*DZ ;
 xR =  0.5D0*L  ;
 xL = -1.0D0*xR ;
 yR =  0.5D0*P  ; 
 yL = -1.0D0*yR ;

 P1 = XR YR 0. ;
 P2 = XL YR 0. ;
 P3 = XL YL 0. ;
 P4 = XR YL 0. ;

 P1P2 = 'DROIT' NX P1 P2 ;
 P2P3 = 'DROIT' NY P2 P3 ;
 P3P4 = 'DROIT' NX P3 P4 ;
 P4P1 = 'DROIT' NY P4 P1 ;

 BASEJR = 'DALLER' P1P2 P2P3 P3P4 P4P1 ;

 POSH = H/2. ;
 NEGH = (-1.)*POSH ;
 MIDDIV = NZ/2 ;
 
 V1 = 0. 0. POSH ;
 V2 = 0. 0. NEGH ;

 'OPTI' 'ELEM' 'CUB8';
 DOM1 = BASEJR 'VOLUME' MIDDIV 'TRAN' V1 ;
 DOM2 = BASEJR 'VOLUME' MIDDIV 'TRAN' V2 ;
 DOMTOT = DOM1 'ET' DOM2 ;

 MDOMTOT = 'MODELISER' DOMTOT 'EULER';
 MDOM1   = 'MODELISER' DOM1 'EULER';
 MDOM2   = 'MODELISER' DOM2 'EULER';

 $DOMTOT = 'DOMA' MDOMTOT 'VF' ;
 $DOM1 = 'DOMA' MDOM1 'VF' ;
 $DOM2 = 'DOMA' MDOM2 'VF' ;

 TDOMTOT = $DOMTOT. 'QUAF';
 TDOM1   = $DOM1. 'QUAF';
 TDOM2   = $DOM2. 'QUAF'; 
 'ELIMINATION' (TDOMTOT 'ET' TDOM1 'ET' TDOM2) 1D-6;
*
******* Creation de la ligne Utilisée pour le Post-Traitement
*       reliant les points centres
*
  XINIT = 0.5D0*DZ ;
  YINIT = 0.5D0*DZ ;
  ZINIT = NEGH '+' (0.5D0*DZ)  ;
  XFIN  = XINIT ;
  YFIN  = YINIT ;
  ZFIN  = POSH '-' (0.5D0*DZ)  ;
  PINI  = XINIT YINIT ZINIT;
  PFIN  =  XFIN  YFIN  ZFIN;
  DIV   = NZ '-' 1 ;
  COURB = PINI 'DROIT' DIV PFIN;
  COURB = COURB 'COULEUR' 'VERT';
  'ELIMINATION'  (1.0D-6) COURB $DOMTOT.CENTRE;

  'SI' GRAPH ;
    'TRACER' (DOMTOT 'ET' COURB)  'TITRE' 'Domaine total' ;
  'FINSI' ;
************************************************************************
*
* INITIAL CONDITIONS
*
************************************************************************
  CP = 2.d0 ;
  CVM = 0.d0 ;
*
*
 VCH1 = ('MANU' 'CHPO' ($DOM1.'CENTRE') 1 'SCAL' 0.1d0 'NATURE'
         'DISCRET') 'ET'
        ('MANU' 'CHPO' ($DOM2.'CENTRE') 1 'SCAL' 0.25d0 'NATURE'
         'DISCRET') ;
 VCH2 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 3 'UVX' 0.d0 'UVY' 0.d0
         'UVZ' 0.d0 'NATURE' 'DISCRET') ;
 VCH3 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 3 'ULX' 0.d0 'ULY' 0.d0
         'ULZ' 0.d0 'NATURE' 'DISCRET') ;
 VCH4 = ('MANU' 'CHPO' ($DOM1.'CENTRE') 1 'SCAL' 10.d6 'NATURE'
         'DISCRET') 'ET'
        ('MANU' 'CHPO' ($DOM2.'CENTRE') 1 'SCAL' 20.d6 'NATURE'
         'DISCRET') ;
 VCH5 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 308.d0 'NATURE'
         'DISCRET') ;
 VCH6 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 308.d0 'NATURE'
         'DISCRET');
 GAMV = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL'  1.4d0 'NATURE'
 'DISCRET');
 CPV  = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 1008.d0 'NATURE'
 'DISCRET');
 GAML = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL'  2.8d0 'NATURE'
 'DISCRET');
 CPL  = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 4186.d0 'NATURE'
 'DISCRET');
 UNO  = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 1.d0 'NATURE'
 'DISCRET');
 PIL  = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 1 'SCAL' 8.5d8 'NATURE'
 'DISCRET');
 
* vch7 (rhov) = gammav*p/((gammav - 1.D0)*Cpv*T)
* vch8 (rhol) = gammal*(p + pil)/((gammal - 1.D0)*Cpl*T)
* IMPROVE, VCH7 & VCH8 ARE NOT ASSOCIATED TO ANY MESH!!!!!!!!!

 NUM1 = 'PSCAL' GAMV  VCH4 ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
 NUM2 = 0.4d0 * ('PSCAL' CPV  VCH5 ('MOTS' 'SCAL') ('MOTS' 'SCAL')) ;
 VCH7 = NUM1 '/' NUM2 ;
 VCH7 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' VCH7 ;
 NUM3 = 2.8d0 * (VCH4 '+' 8.5d8) ;
 NUM4 = 1.8d0 * ('PSCAL' CPL VCH6 ('MOTS' 'SCAL') ('MOTS' 'SCAL')) ; 
 VCH8 = NUM3 '/' NUM4 ;
 VCH8 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' VCH8 ;
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
 'DEBPROC' CONS ;
  'ARGUMENT'VCH1*'CHPOINT' VCH2*'CHPOINT' VCH3*'CHPOINT' VCH4*'CHPOINT'
            VCH5*'CHPOINT' VCH6*'CHPOINT' VCH7*'CHPOINT' VCH8*'CHPOINT'
            ;
            
    WCH1 = 'PSCAL' VCH1 VCH7 ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
    WCH1 = 'EXCO' 'SCAL' WCH1 'SCAL' 'NATURE' 'DISCRET' ;
    AL   = VCH1 '-' 1.0 ;
    AL   = AL * (-1.0) ;
    WCH2 = 'PSCAL' AL VCH8 ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
    WCH2 = 'EXCO' 'SCAL' WCH2 'SCAL' 'NATURE' 'DISCRET' ;    
    WCH3 = WCH1 '*' VCH2 ('MOTS' 'SCAL' 'SCAL' 'SCAL') ('MOTS' 'UVX'
    'UVY' 'UVZ')  ('MOTS' 'UVX' 'UVY' 'UVZ') ;
    WCH3 = 'EXCO' ('MOTS' 'UVX' 'UVY' 'UVZ') WCH3 ('MOTS' 'UVX'
    'UVY''UVZ') 'NATURE' 'DISCRET' ;    
    WCH4 = WCH2 '*' VCH3 ('MOTS' 'SCAL' 'SCAL' 'SCAL') ('MOTS' 'ULX'
    'ULY' 'ULZ') ('MOTS' 'ULX' 'ULY' 'ULZ') ;
    WCH4 = 'EXCO' ('MOTS' 'ULX' 'ULY' 'ULZ') WCH4 ('MOTS' 'ULX' 'ULY'
     'ULZ')'NATURE' 'DISCRET' ;    
    UVX2 = 'PSCAL' VCH2 VCH2 ('MOTS' 'UVX') ('MOTS' 'UVX') ;
    UVY2 = 'PSCAL' VCH2 VCH2 ('MOTS' 'UVY') ('MOTS' 'UVY') ;
    UVZ2 = 'PSCAL' VCH2 VCH2 ('MOTS' 'UVZ') ('MOTS' 'UVZ') ;
    ECV  = 0.5d0 * (UVX2 '+' UVY2 '+' UVY2) ;
* ev  = Cpv*v(5)/gammav         
*    IEV  = ('PSCAL' CPV VCH5 ('MOTS' 'SCAL') ('MOTS' 'SCAL'))'/' 1.4d0 ;
    IEV  = (1008.d0 *  VCH5)'/' 1.4d0 ;
    TIEV = IEV '+' ECV ;
    WCH5 = WCH1 * TIEV ;
    WCH5 = 'EXCO' 'SCAL' WCH5 'SCAL' 'NATURE' 'DISCRET' ;    
    ULX2 = 'PSCAL' VCH3 VCH3 ('MOTS' 'ULX') ('MOTS' 'ULX') ;
    ULY2 = 'PSCAL' VCH3 VCH3 ('MOTS' 'ULY') ('MOTS' 'ULY') ;
    ULZ2 = 'PSCAL' VCH3 VCH3 ('MOTS' 'ULZ') ('MOTS' 'ULZ') ;
    ECL  = 0.5d0 * (ULX2 '+' ULY2 '+' ULZ2) ;
* el  = Cpl*v(6)/gammal + pil/rl
*    NUM9 = ('PSCAL' CPL VCH5 ('MOTS' 'SCAL') ('MOTS' 'SCAL'))'/' 2.8d0 ;
    NUM9 = (4186.d0 * VCH6)'/' 2.8d0 ;
    NUM0 = 8.5d8 '*'(VCH8 '**' -1.) ;
    NUM0 = 'NOMC' 'SCAL' NUM0 ;
    IEL  = NUM9 '+' NUM0 ;
    TIEL = IEL '+' ECL ;
    WCH6 = WCH2 * TIEL ;
    WCH6 = 'EXCO' 'SCAL' WCH6 'SCAL' 'NATURE' 'DISCRET' ;
     
 'RESPRO' WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 ;
 
 'FINPROC' ;
***************************************************************
 WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 = CONS
 VCH1 VCH2 VCH3 VCH4 VCH5 VCH6 VCH7 VCH8 ;
*
* We create old alpha, to satisfy the prim operator
* necessities
*  
 OALP = VCH1 ;
 OTV  = VCH5 ;
 OTL  = VCH6 ;
 
 PINT RL RV TL TV P UL UV ALP = 'PRIM' 'TWOFLUID'
 WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 OALP OTV OTL CP CVM;

 ERRO = 'MAXIMUM' (VCH4 '-' P) 'ABS' ;
 
 'SI' (ERRO > 1.0D-6) ;
    'MESSAGE' 'Problem in the stic file!!!'
    'ERREUR' 5 ;
 'FINSI' ;
*
**** Plot of IC
*
 'SI' GRAPH ;

    MOD1  =  'MODELISER'  ($DOMTOT . 'MAILLAGE' ) 'THERMIQUE';
    
    CHM_V1   =  'KCHA' $DOMTOT 'CHAM' VCH1 ;
    CHM_V2   =  'KCHA' $DOMTOT 'CHAM' VCH2 ;
    CHM_V3   =  'KCHA' $DOMTOT 'CHAM' VCH3 ;
    CHM_V4   =  'KCHA' $DOMTOT 'CHAM' VCH4 ;
    CHM_V5   =  'KCHA' $DOMTOT 'CHAM' VCH5 ;
    CHM_V6   =  'KCHA' $DOMTOT 'CHAM' VCH6 ;
    CHM_V7   =  'KCHA' $DOMTOT 'CHAM' VCH7 ;
    CHM_V8   =  'KCHA' $DOMTOT 'CHAM' VCH8 ;

    'TRACER'  CHM_V1  MOD1 
       'TITR'  ('CHAINE'  'alpha at t='  0.0);
    'TRACER'  CHM_V2  MOD1 
       'TITR'  ('CHAINE'  'uv at t='  0.0);
    'TRACER'  CHM_V3 MOD1
       'TITR'  ('CHAINE'  'ul at t='  0.0);
    'TRACER'  CHM_V4 MOD1
       'TITR'  ('CHAINE'  'p at t='  0.0);
    'TRACER'  CHM_V5 MOD1
       'TITR'  ('CHAINE'  'Tv at t='  0.0);
    'TRACER'  CHM_V6 MOD1
       'TITR'  ('CHAINE'  'Tl at t='  0.0);
    'TRACER'  CHM_V7 MOD1
       'TITR'  ('CHAINE'  'rv at t='  0.0);
    'TRACER'  CHM_V8 MOD1
       'TITR'  ('CHAINE'  'rl at t='  0.0);       

       
 'FINSI' ;
*************************************************************************
*
* MAIN PROGRAM
*
*************************************************************************
* Different upwind scheme available
* AUSM + 
*  METO = 'AUSMP1' ;
* Preconditioned AUSM +
  METO = 'AUSMP2';
* AUSMDV
*  METO = 'AUSMDV1';
* Preconditioned AUSMDV
*  METO = 'AUSMDV2';
* 
* Iterations
*
 NITER = 10000000 ;
*
* Gravity, Final time, CFL-like parameter
*
  g = 0.d0 ;
  TFINAL = 0.005d0 ;
  MCFL = 0.1D0 ; 
*
* Order in space, order in time and a useful counter
*
 ORDESP = 2;
 ORDTEM = 2;    
 COUNT  = 1;
*
* TIME IS RESETED TO 0 seconds
*
  TPS = 0. ;
*
* Names of the residuum components
*
 LISTINCO = 'MOTS' 'RVN' 'RLN' 'RUVX' 'RUVY' 'RUVZ'
                               'RULX' 'RULY' 'RULZ'
                               'RETVN' 'RETLN' ;

 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
 'MESSAGE' ;
*
* Some preliminary calculation for the evaluation of gradients
* 1st. for the gradient of void fraction
* 2nd. for the gradients in 2nd order
* 
 DALP LIMCH GRG = 'PENT' MDOMTOT 'CENTRE'
    'EULESCAL' 'NOLIMITE' ('MOTS' 'SCAL') ALP ;

 'SI' (ORDESP 'EGA' 2) ;
 
    GRADAL CACA COEFSCAL =  'PENT' MDOMTOT 'CENTRE' 'EULESCAL'
                            'NOLIMITE' ('MOTS' 'SCAL') ALP ;
    GRADUV CACA COEFVECT =  'PENT' MDOMTOT 'CENTRE' 'EULEVECT'
                            'NOLIMITE' ('MOTS' 'UX' 'UY' 'UZ')
                            ('EXCO' ('MOTS' 'UVX' 'UVY' 'UVZ')
                              UV    ('MOTS' 'UX' 'UY' 'UZ'))  ;    
 'FINSI' ;

*
**** Temporal loop
*
 'REPETER' BL1 NITER ;
*
**** Evaluation of variables at each side of the interface (PRET operator)
*
   'SI' ((ORDESP 'EGA' 1) 'OU' (&BL1 'EGA' 1)) ;

       RLF  RVF TLF TVF PF ULF UVF ALPF = 'PRET'
       'TWOFLUID' 1 1 $DOMTOT
       ALP UV UL P TV TL RV RL ;

   'SINON' ;
             
       GRADAL LIMALP = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'SCAL') ALP 'GRADGEO' COEFSCAL ;
       GRADUV LIMUV  = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'UVX' 'UVY' 'UVZ') UV  'GRADGEO' COEFVECT ;
       GRADUL LIMUL  = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'ULX' 'ULY' 'ULZ') UL  'GRADGEO' COEFVECT ;
       GRADP  LIMP   = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'SCAL') P   'GRADGEO' COEFSCAL ;
       GRADTV LIMTV  = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'SCAL') TV  'GRADGEO' COEFSCAL ;
       GRADTL LIMTL  = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
       ('MOTS' 'SCAL') TL  'GRADGEO' COEFSCAL ;
          
       RLF  RVF TLF TVF PF ULF UVF ALPF = 'PRET'
       'TWOFLUID' ORDESP ORDTEM $DOMTOT
       ALP GRADAL LIMALP
       UV  GRADUV LIMUV
       UL  GRADUL LIMUL
       P   GRADP  LIMP
       TV  GRADTV LIMTV
       TL  GRADTL LIMTL
       RV  RL DELTAT ;
       
    'FINSI' ;
*
**** Evaluation of Residual and Dt
*    
    RESIDU DELTAT =  'KONV' 'VF' 'TWOFLUID' 'RESI' METO
    $DOMTOT LISTINCO  ALPF UVF ULF PF TVF TLF RVF RLF ;

    DT = MCFL '*' DELTAT ;

    'SI' (COUNT 'EGA' 1 ) ;
      COUNT = 0 ;
      OALP = ALP ;
      OLDDT = DT ;
    'FINSI' ;
*  
**** Increment of the variables (convection)
*
    RESIDU = DT '*' RESIDU ;

    DRVN   = 'EXCO' 'RVN' RESIDU 'SCAL' ;
    DRLN   = 'EXCO' 'RLN' RESIDU 'SCAL' ;
    DRUVN  = 'EXCO' ('MOTS' 'RUVX' 'RUVY' 'RUVZ') RESIDU ('MOTS' 'UVX'
    'UVY' 'UVZ') ;
    DRULN  = 'EXCO' ('MOTS' 'RULX' 'RULY' 'RULZ') RESIDU ('MOTS' 'ULX'
    'ULY' 'ULZ') ;
    DRETVN = 'EXCO' 'RETVN' RESIDU 'SCAL' ;
    DRETLN = 'EXCO' 'RETLN' RESIDU 'SCAL' ;
*
* Analizamos la conservacion de la masa
*
    MASAV = 'XTY' WCH1 UNO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
    MASAL = 'XTY' WCH2 UNO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
    MOMV  = 'XTY' WCH3 UNO ('MOTS' 'UVY')  ('MOTS' 'SCAL') ; 
    MOML  = 'XTY' WCH4 UNO ('MOTS' 'ULY')  ('MOTS' 'SCAL') ;
    ENERV = 'XTY' WCH5 UNO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
    ENERL = 'XTY' WCH6 UNO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ; 
*
**** Previous update of conserved variables
*
    TPS  = TPS  '+' DT ;
    WCH1 = WCH1 '+' DRVN ;
    WCH2 = WCH2 '+' DRLN ;
    WCH3 = WCH3 '+' DRUVN ;
    WCH4 = WCH4 '+' DRULN ;
    WCH5 = WCH5 '+' DRETVN ;
    WCH6 = WCH6 '+' DRETLN ;
************************************************************************
**** Adittion of the source terms, final update of conserved variables
************************************************************************
    DALP1 LIMCH = 'PENT' MDOMTOT 'CENTRE'
    'EULESCAL' 'LIMITEUR' ('MOTS' 'SCAL') ALP 'GRADGEO' GRG  ;
 
    DADP = DT * PINT * DALP1 ;
*
* Source term corresponding to the vapour momentum equation
*    
    S31 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 2 'UVX' 0.d0  'UVY' 0.d0
    'NATURE' 'DISCRET'); 
    S32 =  DT * ALP * RV * g ;
    S321 = 'EXCO' 'SCAL' S32 'UVZ' 'NATURE' 'DISCRET' ;
    S3 = S31 'ET' S321 ;
    DADP = 'EXCO' ('MOTS' 'P1DX' 'P1DY' 'P1DZ') DADP ('MOTS' 'UVX' 'UVY'
    'UVZ') 'NATURE' 'DISCRET' ;
    SO3 =  (S3 '+' DADP) ;
    SO3 = 'NOMC' ('MOTS' 'UVX' 'UVY' 'UVZ') SO3  ('MOTS' 'UVX' 'UVY'
     'UVZ') 'NATURE' 'DISCRET' ;
    WCH3 = WCH3 '+' SO3 ;

* Source term corresponding to the liquid momentum equation

    S41 = ('MANU' 'CHPO' ($DOMTOT.'CENTRE') 2 'ULX' 0.d0 'ULY' 0.d0
    'NATURE' 'DISCRET'); 
    S42 =  DT *(ALP '-' 1.d0) * RL * g ;
    S42 = S42 * (-1.) ;
    S421 = 'EXCO' 'SCAL' S42 'ULZ' 'NATURE' 'DISCRET' ;
    S4 = S41 'ET' S421 ;
    DADP = 'EXCO' ('MOTS' 'UVX' 'UVY' 'UVZ') DADP ('MOTS' 'ULX' 'ULY'
    'ULZ') 'NATURE' 'DISCRET' ;
    SO4 = (S4 '-' DADP) ;
    SO4 = 'NOMC' ('MOTS' 'ULX' 'ULY' 'ULZ') SO4 ('MOTS' 'ULX' 'ULY'
    'ULZ') 'NATURE' 'DISCRET' ;
    WCH4 = WCH4 '+' SO4 ;

* Source term corresponding to the vapour total energy equation               

    VY = 'EXCO' 'UVZ' UV ;
    S5 = S32 * VY ;  
    DADT = DT * PINT * (ALP '-'  OALP) '/' OLDDT ;
    SO5 = (S5 '-' DADT) ;
    SO5 = 'EXCO' 'SCAL' SO5 'SCAL' 'NATURE' 'DISCRET' ;
    WCH5 = WCH5 '+' SO5  ;
    
* Source term corresponding to the liquid total equation

    LY = 'EXCO' 'ULZ' UL ;
    S6 = S42 * LY ;
    SO6 = (S6 '+' DADT) ;
    SO6 = 'EXCO' 'SCAL' SO6 'SCAL' 'NATURE' 'DISCRET' ;              
    WCH6 = WCH6 '+' SO6 ;

    OALP  = ALP ;
    OLDDT = DT  ;
    OTV   = TV  ;
    OTL   = TL  ;
************************************************************************
**** Updating primitive variables
************************************************************************
    PINT RL RV TL TV P UL UV ALP = 'PRIM' 'TWOFLUID'
    WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 OALP OTV OTL CP CVM;
* ANADIR CVM TRAS COMPILAR    
************************************************************************
**** Graphical outputs
************************************************************************
   'SI' (((&BL1 '/' 1000) '*' 1000) 'EGA' &BL1) ;
       'MESSAGE' &BL1  TPS ;
   'FINSI' ;
    
   'SI' (TPS > TFINAL) ;
   
      EVALP = 'EVOL' 'CHPO' ALP COURB ;
      EVUVX = 'EVOL' 'CHPO' UV 'UVX' COURB ;
      EVUVY = 'EVOL' 'CHPO' UV 'UVY' COURB ;
      EVUVZ = 'EVOL' 'CHPO' UV 'UVZ' COURB ;
      EVULX = 'EVOL' 'CHPO' UL 'ULX' COURB ;
      EVULY = 'EVOL' 'CHPO' UL 'ULY' COURB ;
      EVULZ = 'EVOL' 'CHPO' UL 'ULZ' COURB ;
      EVP   = 'EVOL' 'CHPO' P  COURB ;
      EVTV  = 'EVOL' 'CHPO' TV COURB ;
      EVTL  = 'EVOL' 'CHPO' TL COURB ;
      EVRV  = 'EVOL' 'CHPO' RV COURB ;
      EVRL  = 'EVOL' 'CHPO' RL COURB ;

      'MESSAGE' &BL1  TPS ;

      'SI' GRAPH ;
      'DESSIN' EVALP 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'alpha' ;
      'DESSIN' EVUVX 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'uvx (m/s)' ;
      'DESSIN' EVUVY 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'uvy (m/s)' ;
      'DESSIN' EVUVZ 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'uvz (m/s)' ;
      'DESSIN' EVULX 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'ulx (m/s)' ;
      'DESSIN' EVULY 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'uly (m/s)' ;
      'DESSIN' EVULZ 'TITR' 'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'ulz (m/s)' ;
      'DESSIN' EVP  'TITR'  'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'p (Pa)' ;     
      'DESSIN' EVTV 'TITR'  'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'Tv (K)' ;
      'DESSIN' EVTL 'TITR'  'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'Tl (K)' ;
      'DESSIN' EVRV 'TITR'  'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'rv (kg/m3)' ;
      'DESSIN' EVRL 'TITR'  'Shock tube test' 'TITX' 'x (m)' 'TITY'
      'rl (kg/m3)' ;
      'FINSI' ;

      'QUITTER' BL1 ;
      
   'FINSI' ;


   'FIN' BL1 ;
************************************************************************
**** End of loop
************************************************************************

*
 'SI' GRAPH ;
    MOD1  =  'MODELISER'  ($DOMTOT . 'MAILLAGE' ) 'THERMIQUE';
    
    CHM_V1   =  'KCHA' $DOMTOT 'CHAM' alp ;
    CHM_V2   =  'KCHA' $DOMTOT 'CHAM' uv  ;
    CHM_V3   =  'KCHA' $DOMTOT 'CHAM' ul  ;
    CHM_V4   =  'KCHA' $DOMTOT 'CHAM' p   ;
    CHM_V5   =  'KCHA' $DOMTOT 'CHAM' Tv  ;
    CHM_V6   =  'KCHA' $DOMTOT 'CHAM' TL  ;
    CHM_V7   =  'KCHA' $DOMTOT 'CHAM' RV  ;
    CHM_V8   =  'KCHA' $DOMTOT 'CHAM' RL  ;

    'TRACER'  CHM_V1  MOD1 
       'TITR'  ('CHAINE'  'alpha at t='  TPS);
    'TRACER'  CHM_V2  MOD1
           'TITR'  ('CHAINE'  'uv at t=' TPS);
    'TRACER'  CHM_V3 MOD1
       'TITR'  ('CHAINE'  'ul at t='  TPS);
    'TRACER'  CHM_V4 MOD1
       'TITR'  ('CHAINE'  'p at t='  TPS);
    'TRACER'  CHM_V5 MOD1
       'TITR'  ('CHAINE'  'Tv at t='  TPS);
    'TRACER'  CHM_V6 MOD1
       'TITR'  ('CHAINE'  'Tl at t='  TPS);
    'TRACER'  CHM_V7 MOD1
       'TITR'  ('CHAINE'  'rv at t='  TPS);
    'TRACER'  CHM_V8 MOD1
       'TITR'  ('CHAINE'  'rl at t='  TPS);       
       
 'FINSI' ;
************************************************************************
*** saving the results in a file
************************************************************************

* 'OPTION' 'SAUVER' ('CHAINE' 'shock3d.sauv') ;
* 'SAUVER' ;

 'FIN' ;