Télécharger couplage_TH1D_Th3D.dgibi

Retour à la liste

* fichier :  couplage_TH1D_Th3D.dgibi
************************************************************************
************************************************************************
'OPTION' 'ECHO' 0 ;
* From ~/nlin/sources_dev_new/util_proc :
*BEGINPROCEDUR append
************************************************************************
* NOM         : APPEND
* DESCRIPTION : Rajoute :
*               - un entier à un listentier
*               - un réel   à un listreel
*               - un objet (liste, evolution, matrice ou chpoint)
*                 à un indice de table ('MOT' ou 'ENTIER')
*                    * si l'indice n'existe pas
*                    * 'ET' si l'indice existe 
*                
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 10/09/2004, version initiale
* HISTORIQUE : v1, 10/09/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' APPEND ;
'ARGUMENT' tab/'TABLE' ;
'SI' ('EXISTE' tab) ;
   'ARGUMENT' itab/'MOT' ;
   'SI' ('NON' ('EXISTE' itab)) ;
      'ARGUMENT' itab*'ENTIER' ;
   'FINSI' ;      
   lobj = FAUX ;
   'SI' ('NON' lobj) ;
      'ARGUMENT' lr/'LISTREEL' ;
      'SI' ('EXISTE' lr) ;
         obj = lr ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' le/'LISTENTI' ;
      'SI' ('EXISTE' le) ;
         obj = le ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' lev/'EVOLUTION' ;
      'SI' ('EXISTE' lev) ;
         obj = lev ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' lm/'MAILLAGE' ;
      'SI' ('EXISTE' lm) ;
         obj = lm ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' chpo/'CHPOINT' ;
      'SI' ('EXISTE' chpo) ;
         obj = chpo ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' rig/'RIGIDITE' ;
      'SI' ('EXISTE' rig) ;
         obj = rig ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;               
   'SI' ('NON' lobj) ;
      'ARGUMENT' matk/'MATRIK' ;
      'SI' ('EXISTE' matk) ;
         obj = matk ; lobj = VRAI ;
      'FINSI' ;
   'FINSI' ;
   'SI' ('NON' lobj) ;
      cherr = 'CHAINE'
       'Il faut fournir un objet liste, evolution, matrice ou chpoint.'
        ; 
      'ERREUR' cherr ;
   'FINSI' ;      
   'SI' ('EXISTE' tab itab) ;
      tab . itab = 'ET' (tab . itab) obj ;
   'SINON' ;
      tab . itab = obj ;
   'FINSI' ;
   'RESPRO' tab ;      
'FINSI' ;   
'ARGUMENT' lenti/'LISTENTI' ;
'ARGUMENT' lreel/'LISTREEL' ;
'SI' ('EXISTE' lenti) ;
   'ARGUMENT' enti*'ENTIER' ;
   lenti = 'ET' lenti ('LECT' enti) ;
   'RESPRO' lenti ;
'FINSI' ;   
'SI' ('EXISTE' lreel) ;
   'ARGUMENT' reel*'FLOTTANT' ;
   lreel = 'ET' lreel ('PROG' reel) ;
   'RESPRO' lreel ;
'FINSI' ;   
*
* End of procedure file APPEND
*
'FINPROC' ;
*ENDPROCEDUR append
*BEGINPROCEDUR errrel
************************************************************************
* NOM         : ERRREL
* DESCRIPTION : Calcul d'une erreur relative
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/04/2003, version initiale
* HISTORIQUE : v1, 23/04/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' ERRREL ;
'ARGUMENT' val*'FLOTTANT' ;
'ARGUMENT' valref*'FLOTTANT' ;
*
'SI' ('<' ('ABS' valref) 1.D-10) ;
   echref = 1.D0 ;
'SINON' ;
   echref = valref ;   
'FINSI' ;
*
errabs = 'ABS' ('/' ('-' val valref) echref);
*
'RESPRO' errabs ;
*
* End of procedure file ERRREL
*
'FINPROC' ;
*ENDPROCEDUR errrel
*BEGINPROCEDUR exmomod
************************************************************************
* NOM         : EXMOMOD
* DESCRIPTION : Extraction d'un mot d'un listmots
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/06/2003, version initiale
* HISTORIQUE : v1, 23/06/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXMOMOD ;
'ARGUMENT' lm*'LISTMOTS' i*'ENTIER' ;
j = 'DIME' lm ;
k = '+' (MODULO ('-' i 1) j) 1 ;
lemot = 'EXTRAIRE' lm k ;
* Usage de l'opérateur text pour éviter que lemot
* ne soit interprété comme un opérateur
'RESPRO' 'TEXTE' lemot ;
*
* End of procedure file EXMOMOD
*
'FINPROC' ;
*ENDPROCEDUR exmomod
*BEGINPROCEDUR formar
************************************************************************
* NOM         : FORMAR
* DESCRIPTION : formate un réel de facon courte
*               pratique pour les noms de
*               sauvegarde
*               Exemples :
* 'MESSAGE' ('CHAINE' (formar 2.9e5   1)) ;
* 2.9E5
* 'MESSAGE' ('CHAINE' (formar -2.9e5  1)) ;
* -2.9E5
* 'MESSAGE' ('CHAINE' (formar 2.9e-5  1)) ;
* 2.9E-5
* 'MESSAGE' ('CHAINE' (formar -2.9e-5 1)) ;
* -2.9E-5
* 'MESSAGE' ('CHAINE' (formar 2.9     1)) ;
* 2.9
* 'MESSAGE' ('CHAINE' (formar -2.9    1)) ;
* -2.9
* 'MESSAGE' ('CHAINE' (formar 0       1)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 0       1)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 2.9e5   0)) ;
* 3E5
* 'MESSAGE' ('CHAINE' (formar -2.9e5  0)) ;
* -3E5
* 'MESSAGE' ('CHAINE' (formar 2.9e-5  0)) ;
* 3E-5
* 'MESSAGE' ('CHAINE' (formar -2.9e-5 0)) ;
* -3E-5
* 'MESSAGE' ('CHAINE' (formar 2.9     0)) ;
* 3
* 'MESSAGE' ('CHAINE' (formar -2.9    0)) ;
* -3
* 'MESSAGE' ('CHAINE' (formar 0       0)) ;
* 0
* 'MESSAGE' ('CHAINE' (formar 0       0)) ;
* 0
*               
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/02/2003, version initiale
* HISTORIQUE : v1, 18/02/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' FORMAR ;
'ARGUMENT' fl*'FLOTTANT' ;
'ARGUMENT' vir/'ENTIER  ' ;
'SI' ('NON' ('EXISTE' vir)) ;
   vir = 1 ;
'SINON' ;   
   'SI' ('<' vir 0) ;
      'ERREUR' 'fournir un entier positif' ;
   'FINSI' ;
'FINSI' ;         
'SI' ('<' ('ABS' fl) 10.D-100) ;
   chfl = 'CHAINE' '0' ;
'SINON' ;
*! sans le 1.D-10, ca ne fonctionne pas
*! qd on entre pile poil une puissance de 10
   lfl = LOG10 ('ABS' fl) ;
*   lfl = '+' (LOG10 ('ABS' fl)) 1.D-10  ;
   slfl = 'SIGNE' ('ENTIER' lfl) ;
   'SI' ('EGA' slfl 1) ;
      elfl = 'ENTIER' lfl ;
   'SINON' ;
      elfl = '-' ('ENTIER' lfl) 1 ;
   'FINSI' ;
   man = '/' fl ('**' 10.D0 elfl) ;
*
* Une verrue pour des histoires de précision...
*   
   'SI' ('EGA' man 10.D0 ('**' 10.D0 ('*' vir -1.D0))) ;
      man  = '/' man 10.D0 ;
      elfl = '+' elfl 1 ;
   'FINSI' ;
*   
   sman = 'SIGNE' man ;
   'SI' ('EGA' sman 1) ;   
      fman = 'CHAINE' '(F' ('+' vir 2) '.0' vir ')' ;
   'SINON' ;
      fman = 'CHAINE' '(F' ('+' vir 3) '.0' vir ')' ;
   'FINSI' ;
   'SI' ('NEG' vir 0) ;   
      'SI' ('NEG' elfl 0) ;
         chfl = 'CHAINE' 'FORMAT' fman man 'E' elfl ;
      'SINON' ;
         chfl = 'CHAINE' 'FORMAT' fman man ;
      'FINSI' ;
   'SINON' ;
      man2 = 'ENTIER' ('+' man ('*' 0.5D0 sman)) ;
      'SI' ('NEG' elfl 0) ;
         chfl = 'CHAINE' man2 'E' elfl ;
      'SINON' ;
         chfl = 'CHAINE' man2 ;
      'FINSI' ;
   'FINSI' ;            
'FINSI' ;   
'RESPRO' chfl ;
*
* End of procedure file FORMAR
*
'FINPROC' ;
*ENDPROCEDUR formar
*BEGINPROCEDUR log10
************************************************************************
* NOM         : LOG10
* DESCRIPTION : Log_10
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/02/2003, version initiale
* HISTORIQUE : v1, 18/02/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' LOG10 ;
'ARGUMENT' fl/'FLOTTANT' ;
'ARGUMENT' lr/'LISTREEL' ;
'ARGUMENT' cp/'CHPOINT ' ;
'ARGUMENT' cm/'MCHAML  ' ;
'SI' ('EXISTE' fl) ;
   'RESPRO' ('/' ('LOG' fl) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' lr) ;
   'RESPRO' ('/' ('LOG' lr) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' cp) ;
   'RESPRO' ('/' ('LOG' cp) ('LOG' 10.D0)) ;
'FINSI' ;
'SI' ('EXISTE' cm) ;
   'RESPRO' ('/' ('LOG' cm) ('LOG' 10.D0)) ;
'FINSI' ;
*
* End of procedure file LOG10
*
'FINPROC' ;
*ENDPROCEDUR log10
*BEGINPROCEDUR modulo
************************************************************************
* NOM         : MODULO
* DESCRIPTION : Calcule un entier modulo un autre...
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 15/10/2002, version initiale
* HISTORIQUE : v1, 15/10/2002, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' MODULO ;
'ARGUMENT' i*'ENTIER' j*'ENTIER' ;
'SI' ('EGA' j 0) ;
  'MESSAGE' 'Impossible de faire modulo 0' ;
  'ERREUR' 5 ;
'SINON' ;
  k=i '/' j ;
  mod=i '-' ( k '*'j ) ;
  'RESPRO' mod ;
'FINSI' ;
*
* End of procedure file MODULO
*
'FINPROC' ;
*ENDPROCEDUR modulo
* From . :
*BEGINPROCEDUR calcoeil
************************************************************************
* NOM         : CALCOEIL
* DESCRIPTION : Calcul un oeil absolu convenable pour les tracés
*               en donnant un oeil relatif et un maillage
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 16/11/2004, version initiale
* HISTORIQUE : v1, 16/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' CALCOEIL ;
'ARGUMENT' mt*'MAILLAGE' ;
'ARGUMENT' foeil*'POINT' ;
   foeilx foeily foeilz = 'COORDONNEE' foeil ;
   xmt ymt zmt = 'COORDONNEE' mt ;
   maxmt = 'MAXIMUM' xmt ;
   maymt = 'MAXIMUM' ymt ;
   mazmt = 'MAXIMUM' zmt ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   xoeil = '+' maxmt ('*' dxy foeilx) ;
   yoeil = '+' maymt ('*' dxy foeily) ;
   zoeil = '+' mazmt ('*' dz foeilz) ;
   oeil = xoeil yoeil zoeil ;
'RESPRO' oeil ;
*
* End of procedure file CALCOEIL
*
'FINPROC' ;
*ENDPROCEDUR calcoeil
*BEGINPROCEDUR compute
************************************************************************
* NOM         : COMPUTE
* DESCRIPTION : Do the computations of the stationnary fields 
*               in the reactor
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' COMPUTE ;
*
'ARGUMENT' tabdat*'TABLE' ;
*
* First store all the parameters to make tabres self-contained
*
tabres = 'TABLE' 'INCO' ;
tabres . 'tabdat' = tabdat ;
*
* Second, retrieve parameters from tables
*
tabgeo = tabdat . 'tabgeo' ;
tabphy = tabdat . 'tabphy' ;
tabnum = tabdat . 'tabnum' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
*'SI' logb ;
*   $mtifb   = tabgeo . '$mtifb'   ;
*'FINSI' ;   
mcf     = tabgeo . 'mcf'     ;
mfc     = tabgeo . 'mfc'     ;
*
Dh     = tabphy . 'Dh'    ;
chblock= tabphy . 'chblock' ;
*vhe    = tabphy . 'vhe'   ;
qhe    = tabphy . 'qhe'   ;
prhe   = tabphy . 'prhe'  ;
theek  = tabphy . 'theek' ;
pvol   = tabphy . 'pvol'  ;
pvar   = tabphy . 'pvar'  ;
gamma  = tabphy . 'gamma' ;
irrad  = tabphy . 'irrad' ;
pcb    = tabphy . 'pcb'   ;
*
maxit  = tabnum . 'maxit'  ;
rescvg = tabnum . 'rescvg' ;
solvit = tabnum . 'solvit' ;
omcbst = tabnum . 'omcbst' ;
omhe   = tabnum . 'omhe'   ;
*
* Initialization of the tables for the thermal
* computation in the fuel+structure (rvcbst)
* and in the helium (rvhe)
*
rvcbst = 'EQEX' 'ITMA' 1 'NITER' 1 'OMEGA' omcbst 'FIDT' 1000 ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcbst 'OPER' 'LAPN' 'LCBST' 'INCO' 'TCBS' 
         ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtcb 'OPER' 'FIMP' 'PCB' 'INCO' 'TCBS' 
         ;
rvcbst = 'EQEX' rvcbst
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mtif 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'TCBS' 'TCBS' 
         ;
*
* Iterative solver for the fuel+structure part
*
* 1 : direct solver, more robust, slower
* 3 : iterative solver, less robust, faster
'SI' solvit ;
   rvcbst . 'METHINV' . 'TYPINV' = 3 ;
'SINON' ;   
   rvcbst . 'METHINV' . 'TYPINV' = 1 ;
'FINSI' ;   
*         
rvhe = 'EQEX' 'ITMA' 1 'NITER' 1 'OMEGA' omhe 'FIDT' 1000 ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL' 'CENTREE' 'NOCONS'
         'ZONE' $mthe 'OPER' 'KONV' 'RCPHE' 'UHE' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'LAPN' 'LHE' 'INCO' 'THE'
         ;
rvhe = 'EQEX' rvhe
         'OPTI' 'EF' 'IMPL'
         'ZONE' $mthe 'OPER' 'FIMP' 'PHE' 'INCO' 'THE' 
         ;
rvhe = 'EQEX' rvhe
         'CLIM' 'THE' 'TIMP' ('DOMA' $sthe 'MAILLAGE') theek ;
rvcbst . 'INCO' = tabres ;
rvhe . 'INCO' = tabres ;
tabres . 'TCBS' = 'KCHT'  $mtcbst 'SCAL' 'SOMMET' theek ;
lupu lsic lcb okcb = PRPCB irrad theek ;
tabres . 'LCBST' = 'KCHT' $mtcbst 'SCAL' 'CENTRE' lsic
                    ('KCHT' $mtcb 'SCAL' 'CENTRE' lcb) ;
tabres . 'PCB' = pcb ;
tabres . 'T0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'F0' = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'H'  = 'KCHT'  $mtif 'SCAL' 'CENTRE' 0.D0 ;
tabres . 'THE' = 'KCHT' $mthe 'SCAL' 'SOMMET' theek ;
tabres . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' 0.D0 ;
thep  = 'COPIER' (tabres . 'THE') ;
tcbsp = 'COPIER' (tabres . 'TCBS') ;
*
* Initialization of the physical properties of helium
*
rhe muhe lhe cphe okhe = PRPHE prhe theek ;
*qhe = '*' rhe vhe ;
qheb  = '*' qhe chblock ;
qhebi = '*' mcf ('NOEL' $mthe qheb) ;
* Ce KCHT est indispensable car '*' supprimme les valeurs nulles...
qhebi = 'KCHT' $mtif 'SCAL' 'CENTRE' qhebi ;
uhe  = 'NOMC' 'UZ' ('/' qheb rhe) 'NATURE' 'DISCRET' ;
tabres . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhe cphe) ;
tabres . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhe ;
*tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' (0.D0 0.D0 uhe) ;
*tabres . 'QHE' = qhe ;
tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhe ;
tabres . 'QHE' = qheb ;
***************************************
*
* COMPUTATIONAL LOOP
*
***************************************
* Number of non-linearity iterations
maxiter = vrai ;
*
'REPETER' bouc maxit ;
*   1) Compute a mean temperature in helium
   tmoyhe = 'NOEL' $mthe (tabres . 'THE') ;
      'MESSAGE' 'Maximum mean helium temperature (degrees C) : '
             ('-' ('MAXIMUM' tmoyhe) 273.15D0) ;
*   2) Transfer this temperature on the interfaces
   tmtif  = '*' mcf tmoyhe ;
* Ce KCHT est indispensable car '*' supprimme les valeurs nulles...
   tabres . 'T0' = 'KCHT' $mtif 'SCAL' 'CENTRE' tmtif ;
*   3) Compute the exchange coefficient at the interface
*      and the physical properties of the fuel and the
*      structure with the temperature from the
*      previous iteration
* Mixing temperature
   tm = tabres . 'T0' ;
* Wall temperature
   tp = 'NOEL' $mtif (tabres . 'TCBS') ;
   t0 = '*' ('+' tm tp) 0.5D0 ;
*   t0 = tabres . 'T0' ;
   rho0 mu0 l0 cp0 okhe = PRPHE prhe t0 ;
   pran0 = '*' ('*' mu0 cp0) ('INVERSE' l0) ;
*   reyn0 = '*' ('*' qhe Dh) ('INVERSE' mu0) ;
   reyn0 = '*' ('*' qhebi Dh) ('INVERSE' mu0) ;
*  Colburn correlation
   nuss0 = 0.023 '*' ('**' reyn0 0.8D0)
                 '*' ('**' pran0 ('/' 1.D0 3.D0)) ;
   h0 = '/' ('*' nuss0 l0) Dh ;
* Tests
*   asp = '/' H Dh ;
   minasp  = 60.  ;
   minreyn = 1.D4 ;
   minpran = 0.65 ;
   maxpran = 160. ;
   asp = '/' H ('MAXIMUM' Dh) ;
   mireyn = 'MINIMUM' reyn0 ;
   mareyn = 'MAXIMUM' reyn0 ;
   okreyn = '>' mireyn minreyn ;
   mipran = 'MINIMUM' pran0 ;
   mapran = 'MAXIMUM' pran0 ;
*   okpr   = 'ET' ('>' mipran 0.7) ('<' mapran 160.) ;
   okipr =  '>' mipran minpran ;
   okapr =  '<' mapran maxpran ; 
   okpr   = 'ET' okipr okapr ;
   okasp  = '>' asp minasp ;
   oknuss = okreyn 'ET' okpr 'ET' okasp ;
   'SI' ('NON' oknuss) ;
      cherr = 'CHAINE' '!! Correlation de Colburn non valide' ;
      'MESSAGE' cherr ;
      'SI' ('NON' okreyn) ;
         cherr = 'CHAINE' ' min Re=' (formar mireyn 2)
                          ' (should be '  (formar minreyn 2) ')' ;
         'MESSAGE' cherr ;
         cherr = 'CHAINE' '  (Note : max Re=' (formar mareyn 2) ')' ;
        'MESSAGE' cherr ;
      'FINSI' ;          
      'SI' ('NON' okipr) ;         
         cherr = 'CHAINE' ' min Pr=' (formar mipran 2)
                          ' (should be '  (formar minpran 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
      'SI' ('NON' okapr) ;                                   
         cherr = 'CHAINE' ' max Pr=' (formar mapran 2)
                          ' (should be '  (formar maxpran 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
      'SI' ('NON' okasp) ;                                   
         cherr = 'CHAINE' ' min H/Dh=' (formar asp 2)
                          ' (should be '  (formar minasp 2) ')' ;
         'MESSAGE' cherr ;
      'FINSI' ;
   'FINSI' ;      
*   tabres . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' ('*' h0 chblock) ;
   tabres . 'H' =  'KCHT' $mtif 'SCAL' 'CENTRE' h0 ;
   tabres . 'nuss0' = 'KCHT' $mtif 'SCAL' 'CENTRE' nuss0 ;
   tabres . 'reyn0' = 'KCHT' $mtif 'SCAL' 'CENTRE' reyn0 ;
   tabres . 'pran0' = 'KCHT' $mtif 'SCAL' 'CENTRE' pran0 ;
   tabres . 'l0' = 'KCHT' $mtif 'SCAL' 'CENTRE' l0 ;
   tabres . 'Dh' = 'KCHT' $mtif 'SCAL' 'CENTRE' Dh ;
* Physical properties of fuel and structure
   tcbck = 'NOEL' $mtcbst (tabres . 'TCBS') ;
   lupu lsic lcb okcb = PRPCB irrad tcbck ;
   tabres . 'LCBST' = 'KCHT' $mtcbst 'SCAL' 'CENTRE' lsic
                       ('KCHT' $mtcb 'SCAL' 'CENTRE' lcb) ;
*   4) Thermal computation in fuel + structure
   EXEC rvcbst ;
*   5) Flux computation on the surface of the fuel
   rvmdia1 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'TCBS' 'F0' ;
   rvmdia1 . 'INCO' = tabres ;                   
   smdia1 mmdia1 = 'MDIA' (rvmdia1 . '1MDIA') ;
   rvmdia2 = 'EQEX' 'OPTI' 'VF' 'IMPL'
             'ZONE' $mtif 'OPER' 'MDIA' 'H' 'INCO' 'T0'  'F0' ;
   rvmdia2 . 'INCO' = tabres ;                   
   smdia2 mmdia2 = 'MDIA' (rvmdia2 . '1MDIA') ;
   f0 = 'NOMC' 'SCAL' ('-' ('KOPS' mmdia1 '*'
                            ('NOMC' 'TCBS' (tabres . 'TCBS')))
                           ('KOPS' mmdia2 '*'
                            ('NOMC' 'T0'  (tabres . 'T0'))))
                ;
      'MESSAGE' 'Integrated flux on the surface of the fuel (W) : '
             ('MAXIMUM' ('RESULT' f0)) ;
*   6) Transfer it to a volume power in the helium
   tabres . 'F0' = f0 ;
   phe = '*' mfc f0 ;
   tabres . 'PHE' = 'KCHT' $mthe 'SCAL' 'CENTRE' phe ;
*
*   7) Compute the physical properties of helium with temperature and
*      pressure from the previous iteration
*
   theck = 'NOEL' $mthe (tabres . 'THE') ;
   thesk = tabres . 'THE' ;
   rhec muhec lhec cphec okhe = PRPHE prhe theck ;
   rhes muhes lhes cphes okhe = PRPHE prhe thesk ;
*   uhes = 'NOMC' 'UZ' ('*' qhe ('INVERSE' rhes)) ;
   uhes = 'NOMC' 'UZ' ('*' qheb ('INVERSE' rhes)) ;
   tabres . 'RCPHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' ('*' rhec cphec) ;
   tabres . 'LHE' =  'KCHT' $mthe 'SCAL' 'CENTRE' lhec ;
   tabres . 'UHE' =  'KCHT' $mthe 'VECT' 'SOMMET' uhes ;
*   8) Thermal computation in the helium
   EXEC rvhe ;
*
* Convergence check
*
   then  = 'COPIER' (tabres . 'THE') ;
   tcbsn = 'COPIER' (tabres . 'TCBS') ;
   errhe  = '/' ('MAXIMUM'  ('-' then thep) 'ABS')
               ('MAXIMUM'  thep 'ABS') ;
   errcbs = '/' ('MAXIMUM'  ('-' tcbsn tcbsp) 'ABS')
               ('MAXIMUM'  tcbsp 'ABS') ;
   'MESSAGE' ('CHAINE' 'Iteration ' &bouc) ;
   'MESSAGE' ('CHAINE' '   errhe  = ' errhe) ;
   'MESSAGE' ('CHAINE' '   errcbs = ' errcbs) ;
   thep = then ;
   tcbsp = tcbsn ;
   testcvg = 'ET' (errhe '<' rescvg)  (errcbs '<' rescvg) ;
*
* Plots 
*
*   explores tabres ;
   'SI' testcvg ;
      maxiter = faux ;
      'QUITTER' bouc ;
    'FINSI' ;      
'FIN' bouc ;
'SI' ('EGA' maxiter vrai) ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'Warning !!!!!!!!!!!!') ;
   'MESSAGE' ('CHAINE' 'The computation has not converged') ;
   'ERREUR' 'Check compute.procedur !' ;
'FINSI' ;
'RESPRO' tabres ;
*
* End of procedure file COMPUTE
*
'FINPROC' ;
*ENDPROCEDUR compute
*BEGINPROCEDUR cpowdis
************************************************************************
* NOM         : CPOWDIS
* DESCRIPTION : Compute cos (z) power distribution
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' CPOWDIS ;
'ARGUMENT' tabgeo*'TABLE' ;
'ARGUMENT' pvol*'FLOTTANT' ;
'ARGUMENT' gamma*'FLOTTANT' ;
'ARGUMENT' pvar*'LOGIQUE' ;
*
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
*
vcbst = 'MAXIMUM' ('RESULT' ('DOMA' $mtcbst 'VOLUME')) ;
vcb   = 'MAXIMUM' ('RESULT' ('DOMA' $mtcb 'VOLUME'))   ;
vhe   = 'MAXIMUM' ('RESULT' ('DOMA' $mthe 'VOLUME'))   ;
vcoeur = '+' vcbst vhe ;
*
pmoy = pvol '*' ('/' vcoeur vcb) ;
'SI' ('NON' pvar) ;
   pcb = 'KCHT' $mtcb 'SCAL' 'CENTRE' pmoy ;
'SINON' ;   
   raddeg = '/' 180.D0 PI ;
   zcb = 'COORDONNEE' 3 ('DOMA' $mtcb 'CENTRE') ;
   zmoy = '/' H 2.D0 ;
   zcb = '+' zcb ('*' zmoy -1.D0) ;
   pmax = '*' pmoy ('/' gamma ('SIN' ('*' raddeg gamma))) ;
   'MESSAGE' ('CHAINE' 'gamma=' gamma) ;
   'MESSAGE' ('CHAINE' 'pmax=' pmax) ;
   fccos = '/' (2.D0 '*' gamma '*' raddeg) H ;
   pcb1 = '*' ('COS' ('*' zcb fccos)) pmax ;
   chvol = 'DOMA' $mtcb 'VOLUME' ;
   puis1 = 'MAXIMUM' ('RESULT' ('*' chvol pmoy)) ;
   puis2 = 'MAXIMUM' ('RESULT' ('*' chvol pcb1)) ;
  'SI' ('NON' ('EGA' puis1 puis2 ('*' 0.1 ('ABS' puis2)))) ;
      'MESSAGE' ('CHAINE' 'puis1 = ' puis1 ) ;
      'MESSAGE' ('CHAINE' 'puis2 = ' puis2 ) ;
      'ERREUR' 'Les deux puissances devraient etre a peu pres egales' ;
  'FINSI' ;         
* Correction
   facor = '/' puis1 puis2 ;
   pcb2 = '*' pcb1 facor ;
   pcb = 'KCHT' $mtcb 'SCAL' 'CENTRE' pcb2 ;
'FINSI' ;
'RESPRO' pcb ;
*
* End of procedure file CPOWDIS
*
'FINPROC' ;
*ENDPROCEDUR cpowdis
*BEGINPROCEDUR dessevol
************************************************************************
* NOM         : DESSEVOL
* DESCRIPTION : Dessine des évolutions : choisit automatiquement
*               les options, marqueurs, couleurs...
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 16/11/2004, version initiale
* HISTORIQUE : v1, 16/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' DESSEVOL ;
'ARGUMENT' evtot*'EVOLUTION' ;
'ARGUMENT' tabt*'TABLE' ;
'ARGUMENT' tit*'MOT' ;
'ARGUMENT' tix*'MOT' ;
'ARGUMENT' tiy*'MOT' ;
'ARGUMENT' lnclk*'LOGIQUE' ;
'ARGUMENT' nb/'LOGIQUE' ;
*
'SI' ('NON' ('EXISTE' nb)) ;
   nb = FAUX ;
'FINSI' ;   
*
nt  = 'DIME' tabt ;
nev = 'DIME' evtot ;
*
* Attention, dans evtot, il y a une évolution avec des noms de points ?
*
*'SI' ('NEG' nev nt) ;
*   cherr = 'CHAINE' 'Evolution and title table : not same dim.' ;
*   'ERREUR' cherr ;
*'FINSI' ;
*
tev = 'TABLE' ;
tev . 'TITRE' = tabt ;
*
toto = 'TABLE' ;
*
lcoul = 'MOTS' 'TURQ' 'VERT' 'JAUN' 'ROSE' 'ROUG' 'BLEU' ;
lmarq = 'MOTS' 'TRIB' 'TRIA' 'LOSA' 'CARR' 'ETOI' 'PLUS' 'CROI' ;
ltirr = 'MOTS' 'TIRR' 'TIRC' 'TIRL' 'TIRM' ;
*
'SI' nb ;
   ev2 = evtot ;
'SINON' ;
   icou = 0 ;
   'REPETER' iev nev ;
      ii = &iev ;
      evi = 'EXTRAIRE' evtot 'COUR' ii ;
      'SI' ('NEG' ('TYPE' ('EXTRAIRE' evi 'ORDO')) 'LISTMOTS') ;
         icou = '+' icou 1 ;
      'FINSI' ;         
*      ii2 = '/' ('+' ii 1) 2 ;
*      ci  = EXMOMOD lcoul ii2 ;
*      ci  = EXMOMOD lcoul ii ;
      ci  = EXMOMOD lcoul icou ;
      APPEND toto 'EVOLUTION' ('COULEUR' evi ci) ;
   'FIN' iev ;
   ev2 = toto . 'EVOLUTION' ;
'FINSI' ;
*      
'REPETER' iev nev ;
   ii = &iev ;
   mi = EXMOMOD lmarq ii ;
   ti = EXMOMOD ltirr ii ;
   'SI' nb ;
      tev . ii = 'CHAINE' 'MARQ ' mi ' ' ti ;
   'SINON' ;      
      tev . ii = 'CHAINE' 'MARQ ' mi ;
   'FINSI' ;      
'FIN' iev ;   
*
'SI' lnclk ;
   'DESSIN' ev2 'TITR' tit 'TITX' tix 'TITY' tiy 'LEGE' 'MIMA' tev
   'NCLK' ;
'SINON' ;   
   'DESSIN' ev2 'TITR' tit 'TITX' tix 'TITY' tiy 'LEGE' 'MIMA' tev ;
'FINSI' ;   
*
* End of procedure file DESSEVOL
*
'FINPROC' ;
*ENDPROCEDUR dessevol
*BEGINPROCEDUR explodat
************************************************************************
* NOM         : EXPLODAT
* DESCRIPTION : Explore the data interactively or not.
*
* - first argument (mandatory) : the tables with data
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                  an integer  
*              an integer corresponds to a particular view
*   View 1 : the base mesh
*   View 2 : the mesh
*   View 3 : the blocked channels
*   View 4 : the power distribution
*   View 5 : fuel and structure thermal conductivity (T)
*   View 6 : Helium density
*   View 7 : Helium thermal conductivity
*   View 8 : Helium heat capacity
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXPLODAT ;
'ARGUMENT' tabdat*'TABLE';
'ARGUMENT' interact/'LOGIQUE' ;
'SI' ('NON' ('EXISTE' interact)) ;
   interact = VRAI ;
'FINSI' ;   
*
'SI' ('NON' interact) ;
   'ARGUMENT' liview/'LISTENTI' ;
   'SI' ('NON' ('EXISTE' liview)) ;
      'ARGUMENT' iview*'ENTIER' ;
      liview = 'LECT' iview ;
   'FINSI' ;      
   niview = 'DIME' liview ;
   'SI' ('<' niview 1) ;
      'QUITTER' EXPLODAT ;
   'FINSI' ;       
   iiview = 1 ;
'FINSI' ;
*
* Retrieve parameters from tables
*
tabgeo = tabdat . 'tabgeo' ;
tabphy = tabdat . 'tabphy' ;
tabin  = tabgeo . 'tabin' ;
logst  = 'EXISTE' tabin 'STRUCTURE' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
cfc     = tabgeo . 'cfc'     ;
amtcb   = tabgeo . 'amtcb' ;
amthe   = tabgeo . 'amthe'   ;
* QUAF meshes
_mtcbst = 'DOMA' $mtcbst 'QUAF' ;
_mthe   = 'DOMA' $mthe   'QUAF' ;
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mtcb   = 'DOMA' $mtcb   'MAILLAGE' ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
'SI' logst ;
   mtst   = 'DIFF' mtcbst mtcb ;
'FINSI' ;   
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
mtif   = 'DOMA' $mtif   'MAILLAGE' ;
*
Dh     = tabphy . 'Dh'    ;
chblock= tabphy . 'chblock' ;
*vhe    = tabphy . 'vhe'   ;
qhe    = tabphy . 'qhe'   ;
prhe   = tabphy . 'prhe'  ;
theek  = tabphy . 'theek' ;
pvol   = tabphy . 'pvol'  ;
pvar   = tabphy . 'pvar'  ;
gamma  = tabphy . 'gamma' ;
irrad  = tabphy . 'irrad' ;
pcb    = tabphy . 'pcb'   ;
* compute bottom He speed
rhohe muhe lhe cphe dum = PRPHE prhe theek ;
vhe = '/' qhe rhohe ;
*
* Precomputations
*
* One wants the height (z) not to be more than four times x or y dimension
*   targetz = 4.D0 ;
   targetz = 1.5D0 ;
*   targetz = 0.5D0 ;
   ftour = 0.1D0 ;
   _mtot = 'ET' _mtcbst _mthe ;
   xmt ymt zmt = 'COORDONNEE' _mtot ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   dzdx = '/' dz dx ;
   chzscale = '>' dzdx targetz ;
   'SI' chzscale ;
      zfactor = '/' dzdx targetz ;
      factech =  ('-' ('/' 1.D0 zfactor) 1.D0) ;
      dzmt = 'NOMC' 'UZ' ('*' zmt factech) 'NATURE' 'DISCRET' ;
      orig = 'FORME' ;
      defrm = 'FORME' dzmt ;
   'FINSI' ;      
*
   'SI' interact ;                             
      mail = ('COULEUR' mtcb 'ROUG')
        'ET' ('COULEUR' mthe 'BLEU') ;
      'SI' (logst) ;
         mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
      'FINSI' ;                 
      npomail = 'NBNO' mail ;
      nelmail = 'NBEL' mail ;
      titmail = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                                ' NBELEM=' nelmail ;
      TRACMESH mail titmail VRAI ;
   'FINSI' ;      
*
* Precomp of physical properties
*
*  Temperature evolutions (in °C !!!)
   Tmin = 500.D0 ;
   Tmax = 1600.D0 ;
   npasT = 10.D0 ;
   pasT = '/' ('-' Tmax Tmin) npasT ;
   prT = 'PROG' Tmin 'PAS' pasT Tmax ;
   prlupu = 'PROG' ;
   prlsic = 'PROG' ;
   prlcb  = 'PROG' ;
   nbc = 'DIME' prT ;
   'REPETER' ibc nbc ;
      T = 'EXTRAIRE' pRT &ibc ;
      T = '+' T 273.15D0 ;
      lupu lsic lcb dum = PRPCB FAUX  T ;
      prlupu  = 'ET' prlupu  ('PROG' lupu) ;
      prlsic  = 'ET' prlsic  ('PROG' lsic) ;
      prlcb   = 'ET' prlcb   ('PROG' lcb) ; 
   'FIN' ibc ;   
   evlupu = 'EVOL' 'MANU' 'T' prT 'Lambda' prlupu ;
   evlsic = 'EVOL' 'MANU' 'T' prT 'Lambda' prlsic ;
   evlcb =  'EVOL' 'MANU' 'T' prT 'Lambda' prlcb ;
* Helium prop
*  Temperature evolutions (in °C !!!)
Tmin = 450.D0 ;
Tmax = 1150.D0 ;
npasT = 10.D0 ;
pasT = '/' ('-' Tmax Tmin) npasT ;
prT    = 'PROG' Tmin 'PAS' pasT Tmax ;
prrho  = 'PROG' ;
prlamb = 'PROG' ;
prcp   = 'PROG' ;
nbc = 'DIME' prT ;
'REPETER' ibc nbc ;
   T = 'EXTRAIRE' pRT &ibc ;
   T = '+' T 273.15D0 ;
   rho mu lamb cp dum = PRPHE prhe T ;
   prrho = 'ET' prrho ('PROG' rho) ;
   prlamb = 'ET' prlamb ('PROG' lamb) ;
   prcp = 'ET' prcp ('PROG' cp) ;
'FIN' ibc ;
evrho  = 'EVOL' 'MANU' 'T' prT 'Rho' prrho ;
evlamb = 'EVOL' 'MANU' 'T' prT 'Lambda' prlamb ;
evcp   = 'EVOL' 'MANU' 'T' prT 'Cp' prcp ;
*
* Arete precomp
*   
*amtcb = 'ARETE' mtcb ;
*
   tcha = 'TABLE' ;
   tcha . 1 = 'BMesh'  ;
   tcha . 2 = 'Mesh'   ;
   tcha . 3 = 'BChan'  ;
   tcha . 4 = 'FuelT2' ;
   tcha . 5 = 'LFuel'  ;
   tcha . 6 = 'RhoHe'  ;
   tcha . 7 = 'LamHe'  ;
   tcha . 8 = 'CpHe'   ;
   ntcha = 'DIME' tcha ;
   'SI' interact ;
      itcha = 1 ;
      noclic = VRAI ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
      noclic = FAUX ;
   'FINSI' ;            
*
'REPETER' bouc2 ;
'SI' interact ;
   'MESSAGE' 'Noclic allows not to click between each view' ;   
   'MESSAGE'
'However, one will deactivate it in order to zoom or save' ;
   noclic = 'CHOI' 'Choose desired options' noclic ;
'FINSI' ;   
*
'REPETER' bouc1;
*
'SI' ('EGA' itcha 1) ;
   bmesh = ('COULEUR' (tabin . 'FUEL') 'ROUG') 'ET' 
           ('COULEUR' (@STBL (tabin . 'CHANNEL')) 'BLEU') ;
   'SI' logst ;
      bmesh = 'ET' bmesh (tabin . 'STRUCTURE') ;
   'FINSI' ;
   oeil = (0. 0. 1000.) ;
   npomail = 'NBNO' bmesh ;
   nelmail = 'NBEL' bmesh ;
   tit = 'CHAINE' 'Base Mesh ' 'NBPO=' npomail
                               ' NBELEM=' nelmail ;
   'SI' noclic ;
      TRAC oeil bmesh 'TITRE' tit 'NCLK' ;
   'SINON' ;      
      TRAC oeil bmesh 'TITRE' tit  ;
   'FINSI' ;      
'FINSI' ;
'SI' ('EGA' itcha 2) ;
   mail = ('COULEUR' mtcb 'ROUG')
     'ET' ('COULEUR' mthe 'BLEU') ;
   'SI' (logst) ;
      mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
   'FINSI' ;                 
   npomail = 'NBNO' mail ;
   nelmail = 'NBEL' mail ;
   tit = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                          ' NBELEM=' nelmail ;
   TRACMESH mail tit noclic ;
'FINSI' ;
'SI' ('EGA' itcha 3) ;
*   tit = 'CHAINE' 'Blocked channels'  ;
*   tracchml $mtif chblock tit noclic ;
   tit = 'CHAINE' 'Blocked channels'  ;
   tracchpo chblock mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 4) ;
   tit = 'CHAINE' 'Power Distribution (W.m-3)' ;
   tracchml $mtcb pcb amtcb tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 5) ;
   evt = evlupu 'ET' evlsic 'ET' evlcb ;
   tt = 'TABLE' ;
   'SI' irrad ;
      tt . 1 = 'CHAINE' '(U,Pu)C irr.' ;
      tt . 2 = 'CHAINE' 'SiC irr.'     ;
      tt . 3 = 'CHAINE' 'CerCer irr.'  ;
   'SINON' ;      
      tt . 1 = 'CHAINE' '(U,Pu)C' ;
      tt . 2 = 'CHAINE' 'SiC' ;
      tt . 3 = 'CHAINE' 'CerCer' ;
   'FINSI' ;      
   tit = 'CHAINE' 'Lambda(Temperature)' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Lambda (W.m-1.K-1)' ;
   DESSEVOL evt tt tit titx tity noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 6) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Density(Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Rho (kg.m-3)' ;
   DESSEVOL evrho tt tit titx tity noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 7) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Thermal Conductivity (Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Lambda (W.m-1.K-1)' ;
   DESSEVOL evlamb tt tit titx tity noclic ;
'FINSI' ;   
*
'SI' ('EGA' itcha 8) ;
   tt = 'TABLE' ;
   tt . 1 = 'CHAINE' 'Helium' ;
   tit = 'CHAINE' 'Specific Heat (Temperature)'
         ' pres= ' (formar prhe 2) ' bar' ;
   titx = 'CHAINE' 'T (°C)' ;
   tity = 'CHAINE' 'Cp (J.kg-1.K-1)' ;
   DESSEVOL evcp tt tit titx tity noclic ;
'FINSI' ;   
*
'SI' interact ;
   cha = tit ;
   ret = 'MENU' cha 'Next' 'Previous' 'Options' ;
   'SI' ('EGA' ret 'Options') ; 'QUITTER' bouc1  ; 'FINSI';
   'SI' ('EGA' ret 'Quitter') ; 'QUITTER' bouc2; 'FINSI';
   'SI' ('EGA' ret 'Next') ; itcha = '+' itcha 1 ; 'FINSI' ; 
   'SI' ('EGA' ret 'Previous') ; itcha = '-' itcha 1 ; 'FINSI' ; 
   'SI' (itcha > ntcha); itcha = 1      ; 'FINSI';
   'SI' (itcha < 1); itcha = ntcha      ; 'FINSI';
'SINON' ;
   iiview = '+' iiview 1 ;
   'SI' ('>' iiview niview) ;
      'QUITTER' bouc2 ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
   'FINSI' ;
'FINSI' ;            
*
'FIN' bouc1;
'FIN' bouc2 ;
*
'SI' chzscale ;
   'FORME' orig ;
'FINSI' ;      
*
* End of procedure file EXPLODAT
*
'FINPROC' ;
*ENDPROCEDUR explodat
*BEGINPROCEDUR explores
************************************************************************
* NOM         : EXPLORES
* DESCRIPTION : Explore the Results interactively or not.
*
* EXPLORES does graphical output of results
* - first argument (mandatory) : the table with the fields
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                 an integer  
*              an integer corresponds to a particular view
*   View 1 : temperature (fuel + He)
*   View 2 : temperature (fuel)
*   View 3 : temperature (fuel : 2 viewpoints)
*   View 4 : exchange coefficient
*   View 5 : temperature (helium)
*   View 6 : vertical speed (helium)
*   View 7 : temperature profiles per channel (helium)
*   View 8 : vertical speed profiles per channel (helium)
*   View 9  : Nusselt number
*   View 10 : Reynolds number
*   View 11 : Prandtl number
*   View 12 : Lambda Helium
*   View 13 : Hydraulic diameter
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' EXPLORES ;
'ARGUMENT' tabres*'TABLE' ;
'ARGUMENT' interact/'LOGIQUE' ;
'SI' ('NON' ('EXISTE' interact)) ;
   interact = VRAI ;
'FINSI' ;   
*
'SI' ('NON' interact) ;
   'ARGUMENT' liview/'LISTENTI' ;
   'SI' ('NON' ('EXISTE' liview)) ;
      'ARGUMENT' iview*'ENTIER' ;
      liview = 'LECT' iview ;
   'FINSI' ;      
   niview = 'DIME' liview ;
   'SI' ('<' niview 1) ;
      'QUITTER' EXPLORES ;
   'FINSI' ;       
   iiview = 1 ;
'FINSI' ;
*
* Retrieve parameters from tables
*
tabdat = tabres . 'tabdat' ;
tabgeo = tabdat . 'tabgeo' ;
tabin  = tabgeo . 'tabin' ;
logst  = 'EXISTE' tabin 'STRUCTURE' ;
*
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
cfc     = tabgeo . 'cfc'     ;
amtcbst = tabgeo . 'amtcbst' ;
amthe   = tabgeo . 'amthe'   ;
* QUAF meshes
_mtcbst = 'DOMA' $mtcbst 'QUAF' ;
_mthe   = 'DOMA' $mthe   'QUAF' ;
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mtcb   = 'DOMA' $mtcb   'MAILLAGE' ;
'SI' logst ;
   mtst   = 'DIFF' mtcbst mtcb ;
'FINSI' ;   
mtst   = 'DIFF' mtcbst mtcb ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
mtif   = 'DOMA' $mtif   'MAILLAGE' ;
*
* Precomputations
*
* One wants the height (z) not to be more than four times x or y dimension
*   targetz = 4.D0 ;
   targetz = 1.5D0 ;
   ftour = 0.1D0 ;
   _mtot = 'ET' _mtcbst _mthe ;
   xmt ymt zmt = 'COORDONNEE' _mtot ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dxy = 'MAXIMUM' ('PROG' dx dy) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   dzdx = '/' dz dx ;
   chzscale = '>' dzdx targetz ;
   'SI' chzscale ;
      zfactor = '/' dzdx targetz ;
      factech =  ('-' ('/' 1.D0 zfactor) 1.D0) ;
      dzmt = 'NOMC' 'UZ' ('*' zmt factech) 'NATURE' 'DISCRET' ;
      orig = 'FORME' ;
      defrm = 'FORME' dzmt ;
   'FINSI' ;      
*
* manipulation maillage et champ de temperature
   tcbsc = ('-' (tabres . 'TCBS') 273.15D0) ;
   thec  = ('-' (tabres . 'THE')  273.15D0) ;
   cech  = tabres . 'H' ;
   nuss0 = tabres . 'nuss0' ;
   reyn0 = tabres . 'reyn0' ;
   pran0 = tabres . 'pran0' ;
   l0 = tabres . 'l0' ;
   Dh = tabres . 'Dh' ;
   uz = 'EXCO' 'UZ' (tabres . 'UHE') ;
* tournage 
   mtot = 'ET' mtcbst mthe ;
   xmt ymt zmt = 'COORDONNEE' mtot ;
   mixmt = 'MINIMUM' xmt ;
   miymt = 'MINIMUM' ymt ;
   mizmt = 'MINIMUM' zmt ;
   dx = '-' ('MAXIMUM' xmt) ('MINIMUM' xmt) ;
   dy = '-' ('MAXIMUM' ymt) ('MINIMUM' ymt) ;
   dz = '-' ('MAXIMUM' zmt) ('MINIMUM' zmt) ;
   p1 = ('-' mixmt ('*' dx ftour))
        ('-' miymt ('*' dy ftour))
        mizmt ;
   p2 = ('-' mixmt ('*' dx ftour))
        ('-' miymt ('*' dy ftour))
        ('+' mizmt 1.D0) ;
*   amtcbst  = 'ARETE' mtcbst  ;
*   amthe    = 'ARETE' mthe    ;
   tcbsc2 thec2 mtcbst2 mthe2 amtcbst2 amthe2 =
    'TOURNER' tcbsc thec mtcbst mthe amtcbst amthe 180. p1 p2 ;
   tcbsct  = 'ET' tcbsc tcbsc2 ;
   thect   = 'ET' thec  thec2 ;
   mtcbstt = 'ET' mtcbst mtcbst2 ;
   mthet   = 'ET' mthe   mthe2 ;
   amtcbstt = 'ET' amtcbst amtcbst2 ;
   amthet   = 'ET' amthe   amthe2 ;
*
   'SI' interact ;                             
      mail = ('COULEUR' mtcb 'ROUG')
        'ET' ('COULEUR' mthe 'BLEU') ;
      'SI' (logst) ;
         mail = 'ET' mail ('COULEUR' mtst 'BLANC') ;
      'FINSI' ;                 
      npomail = 'NBNO' mail ;
      nelmail = 'NBEL' mail ;
      titmail = 'CHAINE' 'Mesh ' 'NBPO=' npomail
                                ' NBELEM=' nelmail ;
      TRACMESH mail titmail VRAI ;
   'FINSI' ;      
*
* Evol precalc
*
   'SI' chzscale ;
      'FORME' orig ;
   'FINSI' ;      
   nchan = 'DIME' (tabin . 'CHANNEL') ;
   tt = 'TABLE' ;
   tev = 'TABLE' ;
   'REPETER' ichan nchan ;
      lighe = tabgeo . 'tlchan' . &ichan ;
      evthec = 'EVOL' 'CHPO' thec lighe ;
      APPEND tev 'TEMPERATURE' evthec ;
      evuhe = 'EVOL' 'CHPO' uz lighe ;
      APPEND tev 'VITESSE' evuhe ;
      tt . &ichan = 'CHAINE' 'Canal ' &ichan ;
   'FIN' ichan ;
   'SI' chzscale ;
      'FORME' defrm ;
   'FINSI' ;      
*
   tcha = 'TABLE' ;
   tcha . 1 = 'T' ;
   tcha . 2 = 'FuelT' ;
   tcha . 3 = 'FuelT2' ;
   tcha . 4 = 'Excoef' ;
   tcha . 5 = 'HeT' ;
   tcha . 6 = 'HeUz' ;
   tcha . 7 = 'HeevT' ;
   tcha . 8 = 'HeevUz' ;
   tcha . 9 = 'Nuss' ;
   tcha . 10 = 'Reyn' ;
   tcha . 11 = 'Prand' ;
   tcha . 12 = 'LamHe' ;
   tcha . 13 = 'Dh' ;
   ntcha = 'DIME' tcha ;
   'SI' interact ;
      itcha = 1 ;
      noclic = VRAI ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
      noclic = FAUX ;
   'FINSI' ;            
*
'REPETER' bouc2 ;
'SI' interact ;
   'MESSAGE' 'Noclic allows not to click between each view' ;   
   'MESSAGE'
'However, one will deactivate it in order to zoom or save' ;
   noclic = 'CHOI' 'Choose desired options' noclic ;
'FINSI' ;   
*
'REPETER' bouc1;
*
'SI' ('EGA' itcha 1) ;
   tit = 'CHAINE' 'Temp. (°C)' ;
   tracchpo  (tcbsct 'ET' thect)
             (mtcbstt 'ET' mthet)
             (amtcbstt 'ET' amthet)
             tit noclic ;
'FINSI' ;
'SI' ('EGA' itcha 2) ;
   tit = 'CHAINE' 'Fuel Temp. (°C)' ;
   tracchpo tcbsc mtcbst amtcbst tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 3) ;
   tit = 'CHAINE' 'Fuel Temp. (°C)' ;
   tracchpo tcbsct mtcbstt amtcbstt tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 4) ;
   tit = 'CHAINE' 'Exchange coeff. (W.m-2.K-1)'  ;
   tracchml $mtif cech ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 5) ;
   tit = 'CHAINE' 'Temp. He (°C)' ;
   tracchpo thect mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 6) ;
   tit = 'CHAINE' 'Uz He (m.s-1)' ;
   tracchpo uz  mthe amthe tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 7) ;
    tit = 'CHAINE' 'Helium Temperature' ;
    titxthe = 'CHAINE' 's (m)' ;
    titythe = 'CHAINE' 'T (°C)' ;
    DESSEVOL (tev . 'TEMPERATURE') tt tit titxthe titythe noclic ;
'FINSI' ;
'SI' ('EGA' itcha 8) ;
    tit = 'CHAINE' 'Helium Vertical Speed' ;
    titxuhe = 'CHAINE' 's (m)' ;
    tityuhe = 'CHAINE' 'uz (m.s-1)' ;
    DESSEVOL (tev . 'VITESSE') tt tit titxuhe tityuhe noclic ;
'FINSI' ;       
'SI' ('EGA' itcha 9) ;
   tit = 'CHAINE' 'Nusselt number'  ;
   tracchml $mtif nuss0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 10) ;
   tit = 'CHAINE' 'Reynolds number'  ;
   tracchml $mtif reyn0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 11) ;
   tit = 'CHAINE' 'Prandtl number'  ;
   tracchml $mtif pran0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 12) ;
   tit = 'CHAINE' 'Lambda (He)'  ;
   tracchml $mtif l0 ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
'SI' ('EGA' itcha 13) ;
   tit = 'CHAINE' 'Hydraulic diameter (m)'  ;
   tracchml $mtif Dh ('CONTOUR' mtif) tit noclic ;
'FINSI' ;   
*
'SI' interact ;
   cha = tit ;
   ret = 'MENU' cha 'Next' 'Previous' 'Options' ;
   'SI' ('EGA' ret 'Options') ; 'QUITTER' bouc1  ; 'FINSI';
   'SI' ('EGA' ret 'Quitter') ; 'QUITTER' bouc2; 'FINSI';
   'SI' ('EGA' ret 'Next') ; itcha = '+' itcha 1 ; 'FINSI' ; 
   'SI' ('EGA' ret 'Previous') ; itcha = '-' itcha 1 ; 'FINSI' ; 
   'SI' (itcha > ntcha); itcha = 1      ; 'FINSI';
   'SI' (itcha < 1); itcha = ntcha      ; 'FINSI';
'SINON' ;
   iiview = '+' iiview 1 ;
   'SI' ('>' iiview niview) ;
      'QUITTER' bouc2 ;
   'SINON' ;
      itcha = 'EXTRAIRE' liview iiview ;
   'FINSI' ;
'FINSI' ;            
*
'FIN' bouc1;
'FIN' bouc2 ;
*
'SI' chzscale ;
   'FORME' orig ;
'FINSI' ;      
*
* End of procedure file EXPLORES
*
'FINPROC' ;
*ENDPROCEDUR explores
*BEGINPROCEDUR genmema
************************************************************************
* NOM         : GENMEMA
* DESCRIPTION : Generate meshes and coupling matrix given base surface
*               mesh
*
* Input for GENMEMA procedure
* tabin . 'INTERFACE'  . i : i-th interface       (line mesh)
* tabin . 'CHANNEL'    . i : i-th helium channel  (surface mesh)
* tabin . 'FUEL'           : fuel mesh (surface mesh)
* tabin . 'STRUCTURE'      : structure mesh (surface mesh) (optional)
* H : height (m) and nh : number of elements in the height
* (regular meshing in z direction)
*
* Temporary table for GENMEMA procedure
* tabout . 'INTERFACE'  . i : i-th interface       (surface mesh)
* tabout . 'CHANNEL'    . i : i-th helium channel  (volume mesh)
* tabout . 'CCONINCH'   . i : i-th coupling connectivity intreface <-> channel
* tabout . 'CMATCHIN'   . i : i-th coupling matrix channel -> interface
* tabout . 'CMATINCH'   . i : i-th coupling matrix interface -> channel
* tabin . 'FUEL'           : fuel mesh (volume mesh)
* tabin . 'STRUCTURE'      : structure mesh (volume mesh) (optional)
*
* Output table for GENMEMA procedure
*
* output tabgeo should contain all the necessary geometric information
* needed further
* tabgeo . 'tabin'   = tabin
* tabgeo . 'H'       = H   : height
* tabgeo . 'nh'      = nh  : number of elements in the height
* tabgeo . '$mtcbst' = NS model objet (fuel+structure)
* tabgeo . '$mtcb'   = NS model objet (fuel only)
* tabgeo . '$mthe'   = NS model objet (helium channels)
* tabgeo . '$sthe'   = NS model objet (bottom of the helium channels)
* tabgeo . '$mtif'   = NS model objet (interfaces)
* tabgeo . '$mtifb'  = NS model objet (blocked interfaces)
* tabgeo . 'mcf'     = coupling matrix helium -> fuel
* tabgeo . 'mfc'     = coupling matrix fuel -> helium
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/11/2004, version initiale
* HISTORIQUE : v1, 09/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GENMEMA ;
'ARGUMENT' tabin*'TABLE' ;
'ARGUMENT' H*'FLOTTANT' ;
'ARGUMENT' nh*'ENTIER' ;
*
tabout = 'TABLE' ;
tabout . 'INTERFACE' = 'TABLE' ;
tabout . 'CHANNEL'   = 'TABLE' ;
tabout . 'CCONINCH'  = 'TABLE' ;
tabout . 'CMATCHIN'  = 'TABLE' ;
tabout . 'CMATINCH'  = 'TABLE' ;
*
hnh    = '/' H ('FLOTTANT' nh) ;
vechnh = (0. 0. hnh) ;
*
* Treatment of everything except 'OTHER'
*
ninterf = 'DIME' (tabin . 'INTERFACE') ;
nchanne = 'DIME' (tabin . 'CHANNEL') ;
'SI' ('NEG' ninterf nchanne) ;
   cherr = 'CHAINE' 'NEG (DIME (tabin . INTERFACE)) '
                    '(DIME (tabin . CHANNEL))' ;
   'ERREUR' cherr ;                    
'FINSI' ;
*
* Loop on interfaces
*
*   Préconditionnement
*
   mott = 'MOTS' 'T' ;
   pr1  = 'PROG' 0.D0 0.D0 ;
   pr3  = 'PROG' 0.D0 1.D0 ;
*
'REPETER' iinterf ninterf ;
   lpreums = vrai ;
   linterf = tabin  . 'INTERFACE' . &iinterf ;
   scan    = tabin  . 'CHANNEL'   . &iinterf ;
* Loop on height   
   'REPETER' ih nh ;
      vectran = (0. 0. ('*' hnh ('FLOTTANT' ('-' &ih 1)))) ;
      tli = 'PLUS' linterf vectran ;
      tsc = 'PLUS' scan    vectran ;
      sli = 'TRANSLATION' tli 1 vechnh ;
      vsc  = 'VOLUME' tsc 'TRAN' 1 vechnh ;
      cvsc = 'BARYCENTRE' vsc ;
      vvsc = 'MESURE' vsc ;
      ivvsc = '/' 1.D0 vvsc ;
      pr2i  = 'PROG' ivvsc 0.D0 ;
*      'LISTE' cvsc ; 'LISTE' cvsc2 ;
      APPEND (tabout . 'CHANNEL')  &iinterf vsc ;
* Loop on interface elements
      nsli = 'NBEL' sli ;
      'REPETER' isli nsli ;
         esli  = 'ELEM' sli &isli ;
         cesli = 'BARYCENTRE' esli ;
         cco = 'MANUEL' 'SEG2' cesli cvsc ;
         mci = 'MANUEL' 'RIGIDITE' cco mott 'QUEL' pr3 pr1 ;
*                ('PROG' 0.D0 1.D0) ('PROG' 0.D0 0.D0) ;
         mic = 'MANUEL' 'RIGIDITE' cco mott 'QUEL' pr1 pr2i ;
         'SI' lpreums ;
            tabout . 'INTERFACE' . &iinterf = esli ;
            tabout . 'CCONINCH'  . &iinterf = cco ;
            tabout . 'CMATCHIN'  . &iinterf = mci ;
            tabout . 'CMATINCH'  . &iinterf = mic ;
            lpreums = faux ;
         'SINON' ;
            tabout . 'INTERFACE' . &iinterf = 'ET'
               (tabout . 'INTERFACE' . &iinterf) esli ;
            tabout . 'CCONINCH'  . &iinterf = 'ET'
               (tabout . 'CCONINCH'  . &iinterf) cco ;
            tabout . 'CMATCHIN'  . &iinterf = 'ET' 
               (tabout . 'CMATCHIN'  . &iinterf) mci ;
            tabout . 'CMATINCH'  . &iinterf = 'ET'
               (tabout . 'CMATINCH' . &iinterf) mic ;
         'FINSI' ;
      'FIN' isli ;
   'FIN' ih ;
'FIN' iinterf ;                            
*
* Treatment of 'FUEL' and 'STRUCTURE'
*
sother = tabin . 'FUEL' ;
vother = 'VOLUME' sother 'TRAN' nh (0. 0. H) ;
tabout . 'FUEL' = vother ;
tabout . 'FUEL2' = 'VOLUME' sother 'TRAN' 1 (0. 0. H) ;
'SI' ('EXISTE' tabin 'STRUCTURE') ;
   sother = tabin . 'STRUCTURE' ;
   vother = 'VOLUME' sother 'TRAN' nh (0. 0. H) ;
   tabout . 'STRUCTURE' = vother ;
   tabout . 'STRUCT2'   = 'VOLUME' sother 'TRAN' 1 (0. 0. H) ;
'FINSI' ;   
*
* Some eliminations
*
tinter = @STBL (tabout . 'INTERFACE') ;
tother = tabout . 'FUEL' ;
'SI' ('EXISTE' tabout 'STRUCTURE') ;
   tother = 'ET' tother (tabout . 'STRUCTURE') ;
'FINSI' ;   
'ELIMINATION' ('ET' tinter tother) 1.D-6 ;
*
* Second part : tabgeo
*
tabgeo = 'TABLE' ;
*
* Build the global meshes
*
*  Structure & Combustible
mtcb   = tabout . 'FUEL' ;
amtcb  = 'ARETE' (tabout . 'FUEL2') ;
'SI' ('EXISTE' tabout 'STRUCTURE') ;
   mtcbst  = 'ET' mtcb (tabout . 'STRUCTURE') ;
   amtcbst = 'ET' amtcb ('ARETE' (tabout . 'STRUCT2')) ;
'SINON' ;
   mtcbst  = mtcb ;
   amtcbst = amtcb ;
'FINSI' ;   
*  Channels
*mthe = @STBL (tabout . 'CHANNEL') ;
tmthe = tabout . 'CHANNEL' ;
sthe =  @STBL (tabin  . 'CHANNEL') ;
amthe = 'TRANSLATION' ('CONTOUR' sthe) 1 (0. 0. H) ;
*  Line in channels
nchan = 'DIME' (tabin . 'CHANNEL') ;
tlchan = 'TABLE' ;
tpechan = 'TABLE' ;
tpschan = 'TABLE' ;
'REPETER' ichan nchan ;
   selem1 = 'ELEM' (tabin . 'CHANNEL' . &ichan) 1 ;
   po1    = 'POIN' ('CONTOUR' selem1) 1 ;
   po2    =  'PLUS' po1 (0. 0. H) ;
   lighe  = 'DROIT' nh po1 po2 ;
   tlchan . &ichan  = lighe ;
   tpechan . &ichan = po1 ;
   tpschan . &ichan = po2 ;
'FIN' ichan ;   
*  Interfaces
tmtif = tabout . 'INTERFACE' ;
mtif = @STBL (tabout . 'INTERFACE') ;
* Blocked interfaces
lblock = tabin . 'BLOCKED' ;
dblock = 'DIME' lblock ;
logb =  ('>' dblock 0) ;
'SI' logb ;
   tchan = 'TABLE' ;
   'REPETER' iblock dblock ;
      iib = 'EXTRAIRE' lblock &iblock ;
      APPEND tchan 'mtifb' (tabout . 'INTERFACE' . iib) ;
   'FIN' iblock ;
   mtifb = tchan . 'mtifb' ;      
'FINSI' ;
*  Coupling
cfc = @STBL (tabout . 'CCONINCH') ;
mcf = @STBL (tabout . 'CMATCHIN') ;
mfc = @STBL (tabout . 'CMATINCH') ;
mcf = 'CHANGER' mcf 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
mfc = 'CHANGER' mfc 'INCO' ('MOTS' 'T') ('MOTS' 'SCAL')
                           ('MOTS' 'Q') ('MOTS' 'SCAL') ;
*  Generate the QUAF meshes
*   Channels
nhe = 'DIME' tmthe ;
_tmthe = 'TABLE' ;
'REPETER' ihe nhe ;
   _tmthe . &ihe = QUAFME (tmthe . &ihe) ;
'FIN' ihe ;
_mthe = @STBL _tmthe ;
*   Interface
nif = 'DIME' tmtif ;
_tmtif = 'TABLE' ;
'REPETER' iif nif ;
   _tmtif . &iif = QUAFME (tmtif . &iif) ;
'FIN' iif ;
_mtif = @STBL _tmtif ;
*   Other
 _mtcbst _mtcb _sthe =
   QUAFME mtcbst mtcb sthe ;
mcbt = _mtcbst 'ET' _mtcb 'ET' _mtif 'ET' cfc
       'ET' amtcb 'ET' amtcbst ;
'SI' logb ;
   _mtifb = QUAFME mtifb ;
   mcbt = 'ET' mcbt _mtifb ;
'FINSI' ;   
*  Elimination of points with same coordinates
'ELIMINATION' mcbt 1.D-6 ;
ttlchan = @STBL tlchan ;
'ELIMINATION' (_mthe 'ET' _sthe 'ET' cfc 'ET' amthe 'ET' ttlchan)
              1.D-6 ;
* Create the Navier-Stokes model objects
*  Channels
$tmthe = 'TABLE' ;
'REPETER' ihe nhe ;
   $tmthe . &ihe = MODNAV 'LINE' (_tmthe . &ihe) ;
'FIN' ihe ;
*  Interface
$tmtif = 'TABLE' ;
'REPETER' iif nif ;
   $tmtif . &iif = MODNAV 'LINE' (_tmtif . &iif) ;
'FIN' iif ;
$mtcbst $mtcb $mthe $mtif $sthe =
   MODNAV 'LINE' _mtcbst _mtcb _mthe _mtif _sthe ;
* Create the hydraulic diameter table 
tdhif = 'TABLE' ;
'REPETER' iif nif ;
   tdhif . &iif = 'KCHT' ($tmtif . &iif) 'SCAL' 'CENTRE'
                          (tabin . 'DH' .&iif) ;
'FIN' iif ;
dh = @STBL tdhif ;
*
* Store results
*
tabgeo . 'tabin'   = tabin   ;
tabgeo . 'H'       = H       ;
tabgeo . 'nh'      = nh      ;
tabgeo . '$mtcbst' = $mtcbst ;
tabgeo . '$mtcb'   = $mtcb   ;
tabgeo . '$tmthe'  = $tmthe   ;
tabgeo . '$mthe'   = $mthe   ;
tabgeo . '$sthe'   = $sthe   ;
tabgeo . '$mtif'   = $mtif   ;
tabgeo . 'amtcb'   = amtcb   ;
tabgeo . 'amtcbst' = amtcbst ;
tabgeo . 'amthe'   = amthe   ;
tabgeo . 'mcf'     = mcf     ;
tabgeo . 'mfc'     = mfc     ;
tabgeo . 'cfc'     = cfc     ;
'SI' logb ;
   $mtifb = MODNAV 'LINE' _mtifb ;
   tabgeo . '$mtifb' = $mtifb ;
'FINSI' ;
tabgeo . 'Dh' = dh ;
tabgeo . 'tlchan' = tlchan ;
tabgeo . 'tpechan' = tpechan ;
tabgeo . 'tpschan' = tpschan ;
'RESPRO' tabgeo ;
*
* End of procedure file GENMEMA
*
'FINPROC' ;
*ENDPROCEDUR genmema
*BEGINPROCEDUR geomhex1
************************************************************************
* NOM         : GEOMHEX1
* DESCRIPTION : 30° symmetry cell of an hexagonal geometry
*               in the x-y plane
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, ??/??/2004, version initiale
* HISTORIQUE : v1, ??/??/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMHEX1 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' pas*'FLOTTANT' ;
*
* Geometric quantities deduced from the one given above
*
q = '/' pas 2.D0 ;
R = '/' Dh 2.D0 ;
*
*  Geometric discretization parameters 
*  (number of elements)
*  nb (base) : number of elements for the line in (x-y) plane
*              with smallest length
*
nb = 7 ;
s30 = 'SIN' 30.D0 ;
c30 = 'COS' 30.D0 ;
t30 = '/' s30 c30 ;
* length of each line
lcer = ('*' ('/' PI 6.D0) R) ;
lpdia = '-' q R ;
lgdia = '-' ('/' q c30) R ;
lhex  = ('*' q t30) ;
* scaling based on shortest line
lbase = 'MINIMUM' ('PROG' lcer lpdia lgdia lhex) ;
* factors in order to have an almost regular mesh
fcer  = '/' lcer lbase ;
fpdia = '/' lpdia lbase ;
fgdia = '/' lgdia lbase ;
fhex  = '/' lhex lbase ;
* mesh
rbase = 'FLOTTANT' nb ;
ncer  = '*' rbase fcer ;
nhex  = '*' rbase fhex ;
npdia  = '*' rbase fpdia ;
ngdia  = '*' rbase fgdia ;
*
ncer = 'ENTIER' ('+' ncer 0.99D0) ;
nhex = 'ENTIER' ('+' nhex 0.99D0) ;
npdia = 'ENTIER' ('+' npdia 0.99D0) ;
ngdia = 'ENTIER' ('+' ngdia 0.99D0) ;
* Points definition for fuel
p0 = 0.D0 0.D0 0.D0 ;
p1 = R    0.D0 0.D0 ;
p2 = q    0.D0 0.D0 ;
p3 = q    ('*' q t30) 0.D0 ;
p4 = ('*' R c30) ('*' R s30) 0.D0 ;
* Lines definition for fuel
l1 = 'DROIT' npdia p1 p2 ;
l2 = 'DROIT' nhex  p2 p3 ;
l3 = 'DROIT' ngdia p3 p4 ;
l4 = 'CERCLE' ncer p4 p0 p1 ;
* Surfaces definition for fuel
ctn  = l1 'ET' l2 'ET' l3 'ET' l4 ;
stcb = 'SURFACE' ctn 'PLANE' ;
* Points and surface definitions for helium (an isocele triangle with
* surface equal to the corresponding circular sector
Rp = '*' R ('**' ('/' PI 3.D0) 0.5D0) ;
p5 = Rp   0.D0 0.D0 ;
p6 = ('*' Rp c30) ('*' Rp s30) 0.D0 ;
sthe = 'MANUEL' 'TRI3' p0 p5 p6 ;
*
* Output results in a table suitable for GENMEMA
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
*
tabin . 'INTERFACE' . 1 = l4 ;
tabin . 'CHANNEL'   . 1 = sthe ;
tabin . 'FUEL'          = stcb ;
'RESPRO' tabin ;
*
* End of procedure file GEOMHEX1
*
'FINPROC' ;
*ENDPROCEDUR geomhex1
*BEGINPROCEDUR geomhex2
************************************************************************
* NOM         : GEOMHEX2
* DESCRIPTION : Same as GEOMHEX1 but larger symmetry cell
*               in the x-y plane
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 18/11/2004, version initiale
* HISTORIQUE : v1, 18/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMHEX2 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' pas*'FLOTTANT' ;
*
* Geometric quantities deduced from the one given above
*
q = '/' pas 2.D0 ;
R = '/' Dh 2.D0 ;
*
*  Geometric discretization parameters 
*  (number of elements)
*  nb (base) : number of elements for the line in (x-y) plane
*              with smallest length
*
nb = 7 ;
s30 = 'SIN' 30.D0 ;
c30 = 'COS' 30.D0 ;
t30 = '/' s30 c30 ;
* length of each line
lcera = ('*' ('/' PI 6.D0) R) ;
lcerb = ('*' ('/' PI 2.D0) R) ;
lpdia = '-' pas Dh ;
lgdia = '-' ('/' pas c30) R ;
lhex  = '-' ('*' pas t30) R ;
* scaling based on shortest line
lbase = 'MINIMUM' ('PROG' lcera lcerb lpdia lgdia lhex) ;
* factors in order to have an almost regular mesh
fcera = '/' lcera lbase ;
fcerb = '/' lcerb lbase ;
fpdia = '/' lpdia lbase ;
fgdia = '/' lgdia lbase ;
fhex  = '/' lhex lbase ;
* mesh
rbase = 'FLOTTANT' nb ;
ncera  = '*' rbase fcera ;
ncerb  = '*' rbase fcerb ;
nhex  = '*' rbase fhex ;
npdia  = '*' rbase fpdia ;
ngdia  = '*' rbase fgdia ;
*
ncera = 'ENTIER' ('+' ncera 0.99D0) ;
ncerb = 'ENTIER' ('+' ncerb 0.99D0) ;
nhex = 'ENTIER' ('+' nhex 0.99D0) ;
npdia = 'ENTIER' ('+' npdia 0.99D0) ;
ngdia = 'ENTIER' ('+' ngdia 0.99D0) ;
* Points definition for fuel
p0 = 0.D0 0.D0 0.D0 ;
p1 = R    0.D0 0.D0 ;
p2 = ('-' pas R) 0.D0 0.D0 ;
p3 = pas  0.D0 0.D0 ;
p4 = pas  R 0.D0    ;
p5 = pas  ('*' pas t30) 0.D0 ;
p6 = ('*' R c30) ('*' R s30) 0.D0 ;
* Lines definition for fuel
l1 = 'DROIT' npdia p1 p2 ;
l2 = 'CERCLE' ncerb p2 p3 p4 ;
l3 = 'DROIT' nhex  p4 p5 ;
l4 = 'DROIT' ngdia p5 p6 ;
l5 = 'CERCLE' ncera p6 p0 p1 ;
* Surfaces definition for fuel
ctn  = l1 'ET' l2 'ET' l3 'ET' l4 'ET' l5 ;
stcb = 'SURFACE' ctn 'PLANE' ;
* Points and surface definitions for helium (two isocele triangles with
* surface equal to the corresponding circular sector
Rpa = '*' R ('**' ('/' PI 3.D0) 0.5D0) ;
p7 = Rpa   0.D0  0.D0 ;
p8 = ('*' Rpa c30) ('*' Rpa s30) 0.D0 ;
sthea = 'MANUEL' 'TRI3' p0 p7 p8 ;
Rpb = '*' R ('**' ('/' PI 2.D0) 0.5D0) ;
p9  = ('-' pas Rpb) 0.D0 0.D0 ;
p10 = pas Rpb 0.D0 ;
stheb = 'MANUEL' 'TRI3' p9 p3 p10 ;
*
* Output results in a table suitable for GENMEMA
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
*
tabin . 'INTERFACE' . 1 = l5 ;
tabin . 'CHANNEL'   . 1 = sthea ;
tabin . 'INTERFACE' . 2 = l2 ;
tabin . 'CHANNEL'   . 2 = stheb ;
tabin . 'FUEL'          = stcb ;
'RESPRO' tabin ;
*
* End of procedure file GEOMHEX2
*
'FINPROC' ;
*ENDPROCEDUR geomhex2
*BEGINPROCEDUR geompl21
************************************************************************
* NOM         : GEOMPL21
* DESCRIPTION : Type 1 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL21 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   'SI' ('<' &iplate nplate) ;
      APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   'FINSI' ;      
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
   'FINSI' ;      
'FIN' iplate ;    
* Part1 : right structure
vt = ef 0. 0. ;
lcct = 'PLUS' ('ET' lcor lcb) bp ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' lcct nf vt) ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL21
*
'FINPROC' ;
*ENDPROCEDUR geompl21
*BEGINPROCEDUR geompl22
************************************************************************
* NOM         : GEOMPL22
* DESCRIPTION : Type 2 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL22 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = ec2 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p2t nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL22
*
'FINPROC' ;
*ENDPROCEDUR geompl22
*BEGINPROCEDUR geompl23
************************************************************************
* NOM         : GEOMPL23
* DESCRIPTION : Type 3 symmetry cell for GEOMPLA2
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPL23 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
'ARGUMENT' nc2*'ENTIER' ;
'ARGUMENT' nf*'ENTIER' ;
'ARGUMENT' np*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc = '*' nc2 2 ;
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = ef 0. 0. ;
lcorg = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' iplate nplate ;
   p2 = 'PLUS' bp (ep 0. 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
   'SI' ('<' &iplate nplate) ;   
      p2 = 'PLUS' bp (ec 0. 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'FINSI' ;
'FIN' iplate ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
bp = p2 ;
p3 = 'PLUS' bp (ef 0. 0.) ;
lcord = 'DROIT' nf p2 p3 ;
lcor = 'ET' lcorg lcord ;
*
* Part1 : bottom structure
*
bp = 0. 0. 0. ;
vt = 0. ef 0. ;
APPEND tabin 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = 0. ep 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = 0. ec 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t    = 'PLUS' p1   bp ;
      p1tt   = 'PLUS' p1t  vt ;
      p2t    = 'PLUS' p2   bp ;
      p2tt   = 'PLUS' p2t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p1tt nc) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p2t p2tt nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = 0. ec2 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t    = 'PLUS' p1   bp ;
      p1tt   = 'PLUS' p1t  vt ;
      p2t    = 'PLUS' p2   bp ;
      p2tt   = 'PLUS' p2t  vt ;
      APPEND tabin 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p1t p1tt nc2) ;
      APPEND (tabin . 'INTERFACE') iitf ('DROIT' p2t p2tt nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPL23
*
'FINPROC' ;
*ENDPROCEDUR geompl23
*BEGINPROCEDUR geompla1
************************************************************************
* NOM         : GEOMPLA1
* DESCRIPTION : Mesh of a plate geometry : small symmetry cell
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/11/2004, version initiale
* HISTORIQUE : v1, 22/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLA1 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
L4 = '/' L2 2.D0 ;
*
* Discretisation
* 
nplate2 = '/' nplate 2 ;
nc2 = 1 ;
nc = '*' nc2 2 ;
nf = 2 ;
np = 2 ;
*
* Tables
*
tabin1 = 'TABLE' ;
tabin1 . 'INTERFACE' = 'TABLE' ;
tabin1 . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* Part1 : Lines
*
p0 = 0. 0. 0. ; 
p1 = 0. ('*' ef -1.D0) 0. ;
lcor = 'DROIT' nf p0 p1 ;
tablig = 'TABLE' ;
bp = p1 ;
'REPETER' ipl nplate2 ;
   p2 = 'PLUS' bp (0. ('*' ep -1.D0) 0.) ;
   APPEND tablig 'LIG' ('DROIT' bp p2 np) ;
   bp = p2 ;
* The last one is an half channel
   'SI' ('NEG' &ipl nplate2) ;   
      p2 = 'PLUS' bp (0. ('*' ec -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc) ;
      bp = p2 ;
   'SINON' ;      
      p2 = 'PLUS' bp (0. ('*' ec2 -1.D0) 0.) ;
      APPEND tablig 'LIG' ('DROIT' bp p2 nc2) ;
      bp = p2 ;
   'FINSI' ;
'FIN' ipl ;
lcb = tablig . 'LIG' ;
p2  = 'PLUS' p0 bp ;
lch = 'DROIT' p1 p2 1 ;
*
* Part1 : left structure
*
bp = 0. 0. 0. ;
vt = ef 0. 0. ;
APPEND tabin1 'STRUCTURE' ('TRANSLATION' (lcor 'ET' lcb) nf vt) ;
bp = 'PLUS' bp vt ;
*
* Loop on plates
*
'REPETER' iplate nplate2 ;
*  Structure part
   vt = ep 0. 0. ;
   lcort = 'PLUS' lcor bp ;
   APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort np vt) ;
* Fuel part
   lcbt = 'PLUS' lcb bp ;
   APPEND tabin1 'FUEL'  ('TRANSLATION' lcbt np vt) ;
   lcbtt = 'PLUS' lcbt vt ;
   'SI' ('>' &iplate 1) ;
      APPEND (tabin1 . 'INTERFACE') ('-' iitf 1) lcbt ;
   'FINSI' ;      
   APPEND (tabin1 . 'INTERFACE') iitf lcbtt ;
   bp = 'PLUS' bp vt ;
   'SI' ('<' &iplate nplate2) ;
* Structure part 
      vt = ec 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort nc vt) ;
      APPEND (tabin1 . 'INTERFACE') iitf ('DROIT' p1t p2t nc) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin1 . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
      bp = 'PLUS' bp vt ;
* The last one is a half channel      
   'SINON' ;      
* Structure part 
      vt = ec2 0. 0. ;
      lcort = 'PLUS' lcor bp ;
      p1t   = 'PLUS' p1   bp ;
      p2t   = 'PLUS' p1t  vt ;
      APPEND tabin1 'STRUCTURE' ('TRANSLATION' lcort nc2 vt) ;
      APPEND (tabin1 . 'INTERFACE') iitf ('DROIT' p1t p2t nc2) ;
* Channel part
      lcht = 'PLUS' lch bp ;
      APPEND (tabin1 . 'CHANNEL') iitf ('TRANSLATION'  lcht 1 vt) ;
      iitf = '+' iitf 1 ;
   'FINSI' ;      
'FIN' iplate ;
nitf = ('-' iitf 1) ;
*
* Perform rotations
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
pr = 0. 0. 1. ;
'REPETER' bclrot 4 ;
   arot = '*' 90. ('-' &bclrot 1) ;
   append tabin 'STRUCTURE'
      ('TOURNER' (tabin1 . 'STRUCTURE') arot p0 pr) ;
   append tabin 'FUEL'
      ('TOURNER' (tabin1 . 'FUEL') arot p0 pr) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         ('TOURNER' (tabin1 . 'INTERFACE' . &bitf) arot p0 pr) ;
      append (tabin . 'CHANNEL') iitf
         ('TOURNER' (tabin1 . 'CHANNEL' . &bitf) arot p0 pr) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' bclrot ;            
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPLA1
*
'FINPROC' ;
*ENDPROCEDUR geompla1
*BEGINPROCEDUR geompla2
************************************************************************
* NOM         : GEOMPLA2
* DESCRIPTION : Mesh of a plate geometry
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/11/2004, version initiale
* HISTORIQUE : v1, 22/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLA2 ;
'ARGUMENT' Dh*'FLOTTANT' ;
'ARGUMENT' ef*'FLOTTANT' ;
'ARGUMENT' ep*'FLOTTANT' ;
'ARGUMENT' nchan*'ENTIER' ;
*
'SI' ('NEG' (MODULO nchan 2) 1) ;
   cherr = 'CHAINE' 'nchan should be odd' ;
'FINSI' ;   
*
* Geometric quantities deduced from the one given above
*
ec  = '/' Dh 2.D0 ;
ec2 = '/' ec 2 ;
nplate = '+' nchan 1 ;
L2 = '+' ('*' nplate ep) ('*' nchan ec) '+' ('*' ef 2) ;
mL2 = '*' L2 -1.D0 ;
L4 = '/' L2 2.D0 ;
*
* Discretization
*
nc2 = 1 ;
nf  = 2 ;
np  = 2 ;
*
* Type 1 cells 
*
tabin1 = GEOMPL21 Dh ef ep nchan nc2 nf np ;
*
* Type 2 cells
*
tabin2 = GEOMPL22 Dh ef ep nchan nc2 nf np ;
*
* Type 3 cells
*
tabin3 = GEOMPL23 Dh ef ep nchan nc2 nf np ;
*
* Perform rotations and translations to generate total mesh
*
tabin = 'TABLE' ;
tabin . 'INTERFACE' = 'TABLE' ;
tabin . 'CHANNEL'   = 'TABLE' ;
iitf = 1 ;
*
* part1 : Type 1 cells
*
nitf = 'DIME' (tabin1 . 'INTERFACE') ;
append tabin 'STRUCTURE' (tabin1 . 'STRUCTURE') ;
append tabin 'FUEL'      (tabin1 . 'FUEL') ;
'REPETER' bitf nitf ;
   append (tabin . 'INTERFACE') iitf
          (tabin1 . 'INTERFACE' . &bitf) ;
   append (tabin . 'CHANNEL') iitf
          (tabin1 . 'CHANNEL' . &bitf) ;
   iitf = '+' iitf 1 ;
'FIN' bitf ;
*
* part2 : two Type2 cells
*
'REPETER' btype2 2 ;
   'SI' ('EGA' &btype2 1) ;
      arot = 90.D0 ;
      vtra = (L2 0.D0 0.D0) ;
      lsym = VRAI ;
   'FINSI' ;    
   'SI' ('EGA' &btype2 2) ;
      arot = 180.D0 ;
      vtra = (L2 0.D0 0.D0) ;
      lsym = VRAI ;
   'FINSI' ;    
   nitf = 'DIME' (tabin2 . 'INTERFACE') ;
   append tabin 'STRUCTURE'
      (SYTOPL (tabin2 . 'STRUCTURE') arot vtra lsym) ;
   append tabin 'FUEL'
      (SYTOPL (tabin2 . 'FUEL') arot vtra lsym) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         (SYTOPL (tabin2 . 'INTERFACE' . &bitf) arot vtra lsym) ;
      append (tabin . 'CHANNEL') iitf
         (SYTOPL (tabin2 . 'CHANNEL' . &bitf) arot vtra lsym) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' btype2 ;   
*
* part3 : Type 1 cells
*
nitf = 'DIME' (tabin3 . 'INTERFACE') ;
append tabin 'STRUCTURE' (tabin3 . 'STRUCTURE') ;
append tabin 'FUEL'      (tabin3 . 'FUEL') ;
'REPETER' bitf nitf ;
   append (tabin . 'INTERFACE') iitf
          (tabin3 . 'INTERFACE' . &bitf) ;
   append (tabin . 'CHANNEL') iitf
          (tabin3 . 'CHANNEL' . &bitf) ;
   iitf = '+' iitf 1 ;
'FIN' bitf ;
*
* part4 : two Type2 cells
*
'REPETER' btype2 2 ;
   'SI' ('EGA' &btype2 1) ;
      arot = 180.D0 ;
      vtra = (0.D0 0.D0 0.D0) ;
      lsym = FAUX ;
   'FINSI' ;    
   'SI' ('EGA' &btype2 2) ;
      arot = 270.D0 ;
      vtra = (0.D0 0.D0 0.D0) ;
      lsym = FAUX ;
   'FINSI' ;    
   nitf = 'DIME' (tabin2 . 'INTERFACE') ;
   append tabin 'STRUCTURE'
      (SYTOPL (tabin2 . 'STRUCTURE') arot vtra lsym) ;
   append tabin 'FUEL'
      (SYTOPL (tabin2 . 'FUEL') arot vtra lsym) ;
   'REPETER' bitf nitf ;
      append (tabin . 'INTERFACE') iitf
         (SYTOPL (tabin2 . 'INTERFACE' . &bitf) arot vtra lsym) ;
      append (tabin . 'CHANNEL') iitf
         (SYTOPL (tabin2 . 'CHANNEL' . &bitf) arot vtra lsym) ;
      iitf = '+' iitf 1 ;
   'FIN' bitf ;
'FIN' btype2 ;   
*
'RESPRO' tabin ;
*
* End of procedure file GEOMPLA2
*
'FINPROC' ;
*ENDPROCEDUR geompla2
*BEGINPROCEDUR geomplh1
************************************************************************
* NOM         : GEOMPLH1
* DESCRIPTION : Géométrie hexagonale à plaques
*               Technical report DM2S/SFME/LTMF/RT/05-050/A
*               typgeo = 1 small losange cell
*               typgeo = 2 normal hexagonal cell
*               typgeo = 3 large two hexagonal cells
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 21/03/2005, version initiale
* HISTORIQUE : v1, 21/03/2005, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' GEOMPLH1 ;
'ARGUMENT' ept*'FLOTTANT' ;
'ARGUMENT' eps*'FLOTTANT' ;
'ARGUMENT' epp*'FLOTTANT' ;
'ARGUMENT' epg*'FLOTTANT' ;
'ARGUMENT' npp*'ENTIER' ;
'ARGUMENT' typgeo*'ENTIER' ;
*
* On commence par projeter toutes les quantités sur une génératrice
*
ff = '/' 2.D0 ('**' 3.D0 0.5D0) ;
l   = '*' ('/' ept 2.D0) ff ;
eps = '*' eps ff ;
epp = '*' epp ff ;
epg = '*' epg ff ;
eps2 = '/' eps 2.D0 ;
*
res = '-' l (eps2 '+' ('*' epp npp) '+' ('*' epg ('-' npp 1))
                  '+' eps) ;
epg2 = '/' res 2.D0 ;
*
* Discrétisation 
*
ns2 = 1 ;
*ns2 = 2 ;
ns  = '*' ns2 2 ;
ng  = 2 ;
ng2 = 'ENTIER' ('+' ('*' ng ('/' epg2 epg)) 0.5D0) ;
*'MESSAGE' ('CHAINE' 'ng2 = ' ng2) ;
*'MESSAGE' ('CHAINE' 'epg = ' epg) ;
*'MESSAGE' ('CHAINE' 'epg2 = ' epg2) ;
*'MESSAGE' ('CHAINE' 'epg '-' epg2  = ' ('-' epg epg2)) ;
*ng2 = 1 ;
*ng2 = 2 ;
np  = 3 ;
*ns2 = 3 ;
*ns  = 3 ;
*ng  = 3 ;
*ng2 = 3 ;
*np  = 3 ;
*
* Build the horizontal generator
*
itg = 0 ;
tabgen = 'TABLE' ;
bp  = 0.D0 0.D0 0.D0 ;
*  left structure
bpp = 'PLUS'  bp (eps2 0. 0.) ;
ss  = 'DROIT' bp ns2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'ST' ;
tabgen . itg . 'SEG'  = ss ;
*  left gas channel
bpp = 'PLUS' bp (epg2 0. 0.) ;
ss  = 'DROIT' bp ng2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
tabgen . itg . 'SEG'  = ss ;
*  plates and channel
'REPETER' ipp npp ;
   bpp = 'PLUS' bp (epp 0. 0.) ;
   ss  = 'DROIT' bp np bpp ;
   bp  = bpp ;
   itg = '+' itg 1 ;
   tabgen . itg = 'TABLE' ;
   tabgen . itg . 'TYPE' = 'CHAINE' 'FU' ;
   tabgen . itg . 'SEG'  = ss ;
   'SI' ('NEG' &ipp npp) ;
      bpp = 'PLUS' bp (epg 0. 0.) ;
      ss  = 'DROIT' bp ng bpp ;
      bp  = bpp ;
      itg = '+' itg 1 ;
      tabgen . itg = 'TABLE' ;
      tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
      tabgen . itg . 'SEG'  = ss ;
   'FINSI' ;
'FIN' ipp ;   
*
*  right gas channel
bpp = 'PLUS' bp (epg2 0. 0.) ;
ss  = 'DROIT' bp ng2 bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'CH' ;
tabgen . itg . 'SEG'  = ss ;
*  right structure
bpp = 'PLUS'  bp (eps 0. 0.) ;
ss  = 'DROIT' bp ns bpp ;
bp  = bpp ;
itg = '+' itg 1 ;
tabgen . itg = 'TABLE' ;
tabgen . itg . 'TYPE' = 'CHAINE' 'ST' ;
tabgen . itg . 'SEG'  = ss ;
*
* Build the vertical generators
*
prot1 = 0.D0 0.D0 0.D0 ;
prot2 = 0.D0 0.D0 1.D0 ;
arot = 120.D0 ;
tabtmp = 'TABLE' ;
'REPETER' iitg ('DIME' tabgen) ;
   titg = tabgen . &iitg ;
   'SI' ('EGA' (titg . 'TYPE') 'ST') ;
      APPEND tabtmp 'GENSTR' (titg . 'SEG') ;
   'SINON' ;      
      APPEND tabtmp 'GENFU' (titg . 'SEG') ;
   'FINSI' ;
'FIN' iitg ;
genstr  = tabtmp . 'GENSTR' ;   
genfu = tabtmp . 'GENFU' ;
*
genstr  = 'TOURNER' genstr  arot prot1 prot2 ;
genfu   = 'TOURNER' genfu   arot prot1 prot2 ;
gench   = 'DROIT' ('POIN' genfu 'INITIAL') 1 ('POIN' genfu 'FINAL') ;
*gent = 'ET' (tabtmp . 'GENSTR') (tabtmp . 'GENFU') ;
*'LISTE' genstr ;
*'LISTE' genfu ;
*'LISTE' gent ;
*strgen = 'COULEUR' (genstr 'GENERATRICE' gent) 'BLAN' ;
*fugen  = 'COULEUR' (genfu  'GENERATRICE' gent) 'ROUG' ;
*'OPTION' 'DIME' 2 ;
*'TRACER' ('ET' genstr ('COULEUR' gent 'ROUG')) ;
*'TRACER' fugen ;
*'TRACER' (fugen 'ET' strgen) ;
*
* Construire une cellule élémentaire dans tabin1
*
tabin1 = 'TABLE' ;
tabin1 . 'INTERFACE' = 'TABLE' ;
tabin1 . 'CHANNEL'   = 'TABLE' ;
tabin1 . 'DH'        = 'TABLE' ;
ich = 0 ;
bp = 0. 0. 0. ;
'REPETER' iitg ('DIME' tabgen) ;
   titg = tabgen . &iitg ;
   'SI' ('EGA' (titg . 'TYPE') 'ST') ;
      li = titg . 'SEG' ;
      su = ('PLUS' ('ET' genstr genfu) bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
   'SI' ('EGA' (titg . 'TYPE') 'CH') ;
      li  = titg . 'SEG' ;
      ct  = 'CONTOUR' (('PLUS' genfu bp) 'GENERATRICE' li) ;
      dh  = '/' ('*' ('MESURE' ct 'SURF') 4.D0) ('MESURE' ct) ;
*      'MESSAGE' ('CHAINE' 'Dh = ' dh) ;
      li2 = 'DROIT' ('POIN' li 'INITIAL') 1 ('POIN' li 'FINAL') ;
      su  = ('PLUS' gench bp) 'GENERATRICE' li2 ;
      ich = '+' ich 1 ;
      APPEND (tabin1 . 'INTERFACE') ich ct ;
      APPEND (tabin1 . 'CHANNEL')   ich su ;
      tabin1 . 'DH' . ich = dh ;
      su  = ('PLUS' genstr bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
   'SI' ('EGA' (titg . 'TYPE') 'FU') ;
      li  = titg . 'SEG' ;
      su  = ('PLUS' genfu bp) 'GENERATRICE' li ;
      APPEND tabin1 'FUEL' su ;
      su  = ('PLUS' genstr bp) 'GENERATRICE' li ;
      APPEND tabin1 'STRUCTURE' su ;
      bp = 'POIN' li 'FINAL' ;
   'FINSI' ;
'FIN' iitg ;
*
* Construire un bloc de trois cellules 
*
tabin2 = 'TABLE' ;
tabin2 . 'INTERFACE' = 'TABLE' ;
tabin2 . 'CHANNEL'   = 'TABLE' ;
tabin2 . 'DH'        = 'TABLE' ;
iitf = 0 ;
p0 = 0. 0. 0. ;
pr = 0. 0. 1. ;
'REPETER' bclrot 3 ;
   arot = '*' 120.D0 ('-' &bclrot 1) ;
   append tabin2 'STRUCTURE'
      ('TOURNER' (tabin1 . 'STRUCTURE') arot p0 pr) ;
   append tabin2 'FUEL'
      ('TOURNER' (tabin1 . 'FUEL') arot p0 pr) ;
   'REPETER' bich ich ;
      iitf = '+' iitf 1 ;
      append (tabin2 . 'INTERFACE') iitf
         ('TOURNER' (tabin1 . 'INTERFACE' . &bich) arot p0 pr) ;
      append (tabin2 . 'CHANNEL') iitf
         ('TOURNER' (tabin1 . 'CHANNEL' . &bich) arot p0 pr) ;
      tabin2 . 'DH' . iitf = tabin1 . 'DH' . &bich ;         
   'FIN' bich ;
'FIN' bclrot ;
*
* Construire un bloc de 2 hexagones
* 
tabin3 = 'TABLE' ;
tabin3 . 'INTERFACE' = 'TABLE' ;
tabin3 . 'CHANNEL'   = 'TABLE' ;
tabin3 . 'DH'        = 'TABLE' ;
nitf = iitf ;
iitf = 0 ;
lrot = '/' ('*' l 1.5D0) ('COS' 30.) ;
'REPETER' ibloc 2 ;
   'SI' ('EGA' &ibloc 1) ;
      ltra = 0. 0. 0. ;
   'SINON' ;
      arot = '+' 90.D0 ('*' 60.D0 ('-' &ibloc 2)) ;
      ltra = ('*' lrot ('COS' arot))  ('*' lrot ('SIN' arot)) 0. ;
   'FINSI' ;
   append tabin3 'STRUCTURE'
      ('PLUS' (tabin2 . 'STRUCTURE') ltra) ;
   append tabin3 'FUEL'
      ('PLUS' (tabin2 . 'FUEL') ltra) ;
   'REPETER' bitf nitf ;
      iitf = '+' iitf 1 ;
      append (tabin3 . 'INTERFACE') iitf
         ('PLUS' (tabin2 . 'INTERFACE' . &bitf) ltra) ;
      append (tabin3 . 'CHANNEL') iitf
         ('PLUS' (tabin2 . 'CHANNEL' . &bitf) ltra) ;
      tabin3 . 'DH' . iitf = tabin2 . 'DH' . &bitf ;         
   'FIN' bitf ;
'FIN' ibloc ;   
*
* Compute hydraulic diameter
*
'SI' ('EGA' typgeo 1) ;
   tabin = tabin1 ;
'FINSI' ;   
'SI' ('EGA' typgeo 2) ;
   tabin = tabin2 ;
'FINSI' ;   
'SI' ('EGA' typgeo 3) ;
   tabin = tabin3 ;
'FINSI' ;   
*'LISTE' (tabin . 'DH') ;
*fu = 'COULEUR' (tabin . 'FUEL') 'ROUG' ;
*st = 'COULEUR' (tabin . 'STRUCTURE') 'BLAN' ;
*ch = 'COULEUR' (@STBL (tabin . 'CHANNEL')) 'BLEU' ;
*it = 'COULEUR' (@STBL (tabin . 'INTERFACE')) 'JAUN' ;
*'OPTION' 'DIME' 2 ;
*'TRACER' (fu 'ET' st 'ET' ch) ;
*'TRACER' it ;
*'FIN' ;
*
*tabin = 'TABLE' ;
'RESPRO' tabin epg2 ;
*
* End of procedure file GEOMPLH1
*
'FINPROC' ;
*ENDPROCEDUR geomplh1
*BEGINPROCEDUR modnav
************************************************************************
* NOM         : MODNAV
* DESCRIPTION :
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 04/06/2003, version initiale
* HISTORIQUE : v1, 04/06/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' MODNAV ;
'ARGUMENT' dd*'MOT' ;
'REPETER' bcl ;
   'ARGUMENT' mquaf/'MAILLAGE' ;
   'SI' ('EXISTE' mquaf) ;
      modn = 'MODELISER' mquaf 'NAVIER_STOKES' dd ;
      'RESPRO' modn ;                  
   'SINON' ;
      'QUITTER' bcl ;
   'FINSI' ;
'FIN' bcl ;
*
* End of procedure file MODNAV
*
'FINPROC' ;
*ENDPROCEDUR modnav
*BEGINPROCEDUR prpcb
************************************************************************
* NOM         : PRPCB
*-----------------------------------------------------------------------
*                Thermal conductivity of a CerCer fuel
*-----------------------------------------------------------------------
* Input : T ou TS : temperature (CHPOINT or FLOTTANT)
*                   (degree K)
*         Validity domaine : 500 to 1300 degree C
*         IRRAD : Irradiated fuel ?
* Output : LUPU : thermal conductivity of (U,Pu) same type as T
*          LSIC : thermal conductivity of SiC same type as T
*          LCB  : thermal conductivity of CerCer fuel same type as T
*                 (W.m-1.K-1)
*-----------------------------------------------------------------------
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : 
**********************************************************************
* VERSION    : v1.1, 06/10/2004 Test domaine de validité
* VERSION    : v1, 20/10/2003, version initiale
* HISTORIQUE : v1, 20/10/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
*
*
'DEBPROC' PRPCB ;
'ARGUMENT' irrad*'LOGIQUE' ;
'ARGUMENT' t/'CHPOINT' ;
'SI' ('EXISTE' t) ;
   typarg = 'MOT' 'CHPOINT' ;
'SINON' ;
   'ARGUMENT' t*'FLOTTANT' ;
   typarg = 'MOT' 'FLOTTANT' ;
'FINSI';
*
* Validity domain (in °C)
*
vdinfc = 499.D0 ;
vdsupc = 1301.D0 ;
*         PORO : Porosity (controled size)
*         BETA : Coefficient to take into account the porosity
*                in the Maxwell-Eucken law
*         FVOL : Volume fraction of (U,Pu)C inclusions in the SiC
*                matrix
*                 If irrad=true, we choose a
*         combrat (at%) combustion ratio of 10 at% for (U,PU)C
*         tirradc (°C)  and an irradiation temperature of
*                       1125°C for the SiC .
poro = 0.15D0 ;
beta = 2.D0 ;
fvol = 0.7D0 ;
combrat = 10.D0 ;
tirradc = 1125.D0 ;
*
* Coefficients for the correlations
*
* (U,Pu)C
c1 = 12.76D0 ;
c2 = 8.71D-3 ;
c3 = -1.88D-6 ;
* SiC
a = 42.58D0 ;
b = -1.5564D4 ;
c = 1.2977D7 ;
d = -1.8458D9 ;
* SiC irradié
e = 3.91112D-2 ;
f = 2.24732D-3 ;
*
'SI' ('EGA' typarg 'FLOTTANT') ;
   t1 = 1.D0 ;
'SINON' ;
   t = 'NOMC' 'SCAL' t ;
   tmail = 'EXTRAIRE'  t 'MAIL' ;
   t1 = 'MANUEL' 'CHPO' tmail 1 'SCAL' 1.D0 'NATURE' 'DISCRET' ;
'FINSI' ;
*
tc = '-' t ('*' t1 273.15D0) ;
*
* Test if temperature is in the validity range 
*
'SI' ('EGA' typarg 'FLOTTANT') ;
   mintc = tc ; maxtc = tc ;
'SINON' ;
   mintc = 'MINIMUM' tc ;
   maxtc = 'MAXIMUM' tc ;
'FINSI' ;
testval = 'ET' ('>' mintc vdinfc) ('<' maxtc vdsupc) ;
*
tc2 = '**' tc 2 ;
* (U,Pu)C non irradiated and 100% dense
lupu1 = '+' ('+' ('*' tc c2) ('*' tc2 c3))
               c1 ;
* (U,Pu)C non irradiated, porous
cpor = '/' ('-' 1.d0 poro) ('+' 1.d0 ('*' beta poro)) ;
lupu2 = '*' lupu1 cpor ;
* (U,Pu)C irradiated porous
'SI' irrad ;
   cirrad = ('-' 1.d0 ('*' 0.02d0 combrat)) ;
   lupu = '*' lupu2 cirrad ;
'SINON' ;
   lupu = lupu2 ;
'FINSI' ;
* SiC non irradiated and 100% dense
it = '/' t1 t ;
it2 = '**' it 2 ;
it3 = '*' it2 it ;
lsic1 = '+' ('+' ('*' it b)
                 ('+' ('*' it2 c) ('*' it3 d)))
            a ;
* SiC irradiated               
'SI' irrad ;
   cirrad = '*' ('EXP' ('*' tirradc f))
                e ;
   lsic = '*' lsic1 cirrad ;
'SINON' ;
   lsic = lsic1 ;                                  
'FINSI' ;
* CerCer
lmat  = lsic ;
lincl = lupu ;
lmat2 = '*' lmat 2.d0 ;
lmi = '-' lmat lincl ;
numer = '+' lincl
            ('+' lmat2
                 ('*' lmi ('*' -2.d0 fvol))) ;
denom = '+' lincl
            ('+' lmat2
                 ('*' lmi fvol)) ;
lcb = '*' lmat ('/' numer denom) ;
'RESPRO' lupu lsic lcb testval ;
*
* End of procedure file PRPCB
*
'FINPROC' ;
*ENDPROCEDUR prpcb
*BEGINPROCEDUR prphe
************************************************************************
* NOM         : PRPHE
*-----------------------------------------------------------------------
*               Thermal properties of Helium 4
*-----------------------------------------------------------------------
* Input : ptot : pressure (bars)
*         T    : temperature (K) (CHPOINT or FLOTTANT)
*
* Output : RHOHE : density (kg.m-3)
*          MUHE  : viscosity (kg.m-1.s-1)
*          LHE   : thermal conductivity (W.m-1.K-1)
*          CPHE  : specific heat (at constant pressure) (J.kg-1.K-1)
*-----------------------------------------------------------------------
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : 
**********************************************************************
* VERSION    : v1.1, 06/10/2004, test + simplification
* VERSION    : v1, 26/02/2003, version initiale
* HISTORIQUE : v1, 26/02/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' PRPHE ;
'ARGUMENT' ptot*'FLOTTANT' ;
'ARGUMENT' t/'CHPOINT' ;
'SI' ('EXISTE' t) ;
   typarg = 'MOT' 'CHPOINT' ;
'SINON' ;
   'ARGUMENT' t*'FLOTTANT' ;
   typarg = 'MOT' 'FLOTTANT' ;
'FINSI';
*
pref = 101325. ;
ppa = '*' ptot pref ;
*
* Validity domain
*
vdpinf = 1. ;
vdpsup = 100. ;
vdtinf = 273. ;
vdtsup = 1500. ;
* Perfect gas constant for helium
rhe = 2077.22 ;
cphe = 5193. ;
C1 = 9.489433E-04 ;
C2 = 9.528079E-04 ;
C3 = 3.420680E-02 ;
C4 = 2.739470E-03 ;
C5 = 9.409120E-04 ;
*
'SI' ('EGA' typarg 'FLOTTANT') ;
   t1 = 1.D0 ;
'SINON' ;
   t = 'NOMC' 'SCAL' t ;
   tmail = 'EXTRAIRE'  t 'MAIL' ;
   t1 = 'MANUEL' 'CHPO' tmail 1 'SCAL' 1.D0 'NATURE' 'DISCRET' ;
'FINSI' ;
* Test validity range 
*
minp = ptot ;
maxp = ptot ;
'SI' ('EGA' typarg 'FLOTTANT') ;
   mint = t ; maxt = t ;
'SINON' ;
   mint = 'MINIMUM' t ;
   maxt = 'MAXIMUM' t ;
'FINSI' ;
testval = 'ET' ('ET' ('>' mint vdtinf) ('<' maxt vdtsup))
               ('ET' ('>' minp vdpinf) ('<' maxp vdpsup))
 ;
*
b =     ('/' ('*' t1 c4) ('+' t1 ('*' t c5)))
    '+' ('/' ('*' t1 c2) ('-' t1 ('*' t c3)))
    '+' c1 ;
rhohe = '*' ('/' t1 ('+' ('*' t rhe) ('*' b ppa))) ppa ;
muhe  = '*' ('**' t 0.687) 3.953D-7 ;
lhe   = '*' ('**' t 0.701) 2.774D-3 ;
'RESPRO' rhohe muhe lhe cphe testval ;
*
* End of procedure file PRPHE
*
'FINPROC' ;
*ENDPROCEDUR prphe
*BEGINPROCEDUR quafme
************************************************************************
* NOM         : QUAFME
* DESCRIPTION :
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 04/06/2003, version initiale
* HISTORIQUE : v1, 04/06/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' QUAFME ;
'REPETER' bcl ;
   'ARGUMENT' mquad/'MAILLAGE' ;
   'SI' ('EXISTE' mquad) ;
      mquaf = 'CHANGER' mquad 'QUAF' ;
      'RESPRO' mquaf ;                  
   'SINON' ;
      'QUITTER' bcl ;
   'FINSI' ;
'FIN' bcl ;
*
* End of procedure file QUAFME
*
'FINPROC' ;
*ENDPROCEDUR quafme
*BEGINPROCEDUR summdat
************************************************************************
* NOM         : SUMMDAT
* DESCRIPTION : 
* SUMMDAT prints a summary of the datas.
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' SUMMDAT ;
'ARGUMENT' tabdat*'TABLE' ;
*
* Retrieve parameters from tables
*
tabgeo = tabdat . 'tabgeo' ;
tabin  = tabgeo . 'tabin'  ;
tabphy = tabdat . 'tabphy' ;
tabnum = tabdat . 'tabnum' ;
*
ikas = tabdat . 'ikas' ; 
*
lblock  = tabin . 'BLOCKED' ;
*
H       = tabgeo . 'H'       ;
$mtcbst = tabgeo . '$mtcbst' ;
$mthe   = tabgeo . '$mthe'   ;
* QUAF meshes
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
*
tabphy = tabdat . 'tabphy' ;
* Interface
Dh = tabphy . 'Dh' ;
* Helium
* vhe   = tabphy . 'vhe'   ;
 qhe   = tabphy . 'qhe'   ;
 prhe  = tabphy . 'prhe'  ;
 theek = tabphy . 'theek' ;
 theec = '-' theek 273.15D0 ;
* compute bottom He speed
 rhohe muhe lhe cphe dum = PRPHE prhe theek ;
 vhe = '/' qhe rhohe ;
* Fuel 
 pvol  = tabphy . 'pvol'  ;
 pvar  = tabphy . 'pvar'  ;
 gamma = tabphy . 'gamma' ;
 irrad = tabphy . 'irrad' ;
 maxit  = tabnum . 'maxit'  ;
 rescvg = tabnum . 'rescvg' ;
 solvit = tabnum . 'solvit' ;
*
'MESSAGE' '**************' ;
'MESSAGE' ('CHAINE' 'Case ' ikas) ;
'SI' ('ET' ('>EG' ikas 1)  ('&lt;EG' ikas 3)) ;
   chcore = 'CHAINE' 'Hexagonal core' ;
'FINSI' ;   
'SI' ('ET' ('>EG' ikas 4)  ('&lt;EG' ikas 6)) ;
   chcore = 'CHAINE' 'Plate core' ;
'FINSI' ;   
'SI' ('ET' ('>EG' ikas 7)  ('&lt;EG' ikas 29)) ;
   chcore = 'CHAINE' 'Hexagonal plate core' ;
'FINSI' ;   
'SI' ('>EG' ikas 30) ;
   chcore = 'CHAINE' 'Modified hexagonal plate core' ;
'FINSI' ;   
dlblock = 'DIME' lblock ;
'SI' ('EGA' dlblock 0) ;
   cblock = 'CHAINE' 'with no blocked channels' ;
'SINON' ;   
   cblock = 'CHAINE' 'with ' dlblock ' blocked channels :' ;
'FINSI' ;
'MESSAGE' ('CHAINE' chcore ' ' cblock) ;
'SI' ('NEG' dlblock 0) ;
   cnum = 'CHAINE' ' '  ;
   'REPETER' ilblock dlblock ;
      cnum = 'CHAINE' cnum ('EXTRAIRE' lblock &ilblock) ' ' ;
   'FIN' ilblock ;
   'MESSAGE' cnum ;
   pcflow= 'ENTIER' ('*' fblock 100.) ;
   'SI' ('>' pcflow 0) ;
      'MESSAGE' ('CHAINE' 'Residual flow rate = ' pcflow ' %') ;
   'SINON' ;      
      'MESSAGE' ('CHAINE' 'Residual flow rate = ' (formar pcflow 2)) ;
   'FINSI' ;
'FINSI' ;
*
* Geometry
*
'MESSAGE' '-----------' ;
'MESSAGE' 'Geometry :' ;
'MESSAGE' ('CHAINE' 'Height = ' (formar H 2) ' m') ;
mail = mtcbst 'ET' mthe ;
npomail = 'NBNO' mail ;
nelmail = 'NBEL' mail ;
'MESSAGE' ('CHAINE' 'Mesh : nb. points=' npomail
                    ' ; nb. elem=' nelmail) ;
*
* Physics
*
'MESSAGE' '-----------' ;
'MESSAGE' 'Physics :' ;
dhmin = 'MINIMUM' Dh ;
dhmax = 'MAXIMUM' Dh ;
'MESSAGE' ('CHAINE' 'Min hydraulic diameter = ' (formar dhmin 2) ' m') ;
'MESSAGE' ('CHAINE' 'Max hydraulic diameter = ' (formar dhmax 2) ' m') ;
'MESSAGE' ('CHAINE' 'He pressure            = '
                    (formar prhe 2) ' bar') ;
'MESSAGE' ('CHAINE' 'He flow rate          = ' (formar qhe  2)
                    ' kg.m-2.s-1') ;
'MESSAGE' ('CHAINE' 'Bottom He speed       = ' (formar vhe  2)
                    ' m.s-1') ;
'MESSAGE' ('CHAINE' 'Bottom He temp.       = ' (formar theec 2)
                    ' degrees C') ;
'MESSAGE' ('CHAINE' 'Volumetric core power = ' (formar pvol  2)
                                               ' W.m-3') ;
'SI' pvar ;
   chpvar = 'CHAINE'
        'Z-cosine power profile with peaking factor gamma = '
        (formar gamma 2) ;
'SINON' ;        
   chpvar = 'Z-constant power profile ' ;
'FINSI' ;   
'MESSAGE' chpvar ;
'SI' irrad ;
   chirr = 'The fuel is irradiated ' ;
'SINON' ;        
   chirr = 'The fuel is new ' ;
'FINSI' ;   
'MESSAGE' chirr ;
*
* Numerics
*
'MESSAGE' '-----------' ;
'MESSAGE' 'Numerics :' ;
'MESSAGE' ('CHAINE' 'Max. nb of non-linear iterations : ' maxit) ;
'MESSAGE' ('CHAINE' 'Max. temperature residual        : '
                    (formar rescvg 2)) ;
'SI' solvit ;
   chsol = 'CHAINE' 'Iterative linear solver' ;
'SINON' ;
   chsol = 'CHAINE' 'Direct linear solver' ;
'FINSI' ;
'MESSAGE' chsol ;   
'MESSAGE' ('CHAINE' '3D Thermics relaxation factor    : '
                    (formar omcbst 2)) ;
'MESSAGE' ('CHAINE' '1D Thermohydraulics rel. factor  : '
                    (formar omhe 2)) ;
'MESSAGE' '**************' ;
*
* End of procedure file SUMMDAT
*
'FINPROC' ;
*ENDPROCEDUR summdat
*BEGINPROCEDUR summres
************************************************************************
* NOM         : SUMMRES
* DESCRIPTION : 
* SUMMRES prints a summary of the results of the computation
*         it also does some validation tests (hence the return value).
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 17/11/2004, version initiale
* HISTORIQUE : v1, 17/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' SUMMRES ;
'ARGUMENT' tabres*'TABLE' ;
*
* Retrieve parameters from tables
*
tabdat = tabres . 'tabdat' ;
ikas   = tabdat . 'ikas' ;
tabgeo = tabdat . 'tabgeo' ;
tabin  = tabgeo . 'tabin' ;
H       = tabgeo . 'H'       ;
nh      = tabgeo . 'nh'      ;
$mtcbst = tabgeo . '$mtcbst' ;
$mtcb   = tabgeo . '$mtcb'   ;
$mthe   = tabgeo . '$mthe'   ;
$sthe   = tabgeo . '$sthe'   ;
$mtif   = tabgeo . '$mtif'   ;
* QUAF meshes
_mtcbst = 'DOMA' $mtcbst 'QUAF' ;
_mthe   = 'DOMA' $mthe   'QUAF' ;
mtcbst = 'DOMA' $mtcbst 'MAILLAGE' ;
mthe   = 'DOMA' $mthe   'MAILLAGE' ;
mtif   = 'DOMA' $mtif   'MAILLAGE' ;
*
tabphy = tabdat . 'tabphy' ;
prhe = tabphy . 'prhe' ;
*
tcpus = tabres . 'tcpus' ;
*
testerr = FAUX ;
'MESSAGE' '**************' ;
'MESSAGE' ('CHAINE' 'Case ' ikas) ;
mail = mtcbst 'ET' mthe ;
npomail = 'NBNO' mail ;
nelmail = 'NBEL' mail ;
'MESSAGE' ('CHAINE' 'nb. points=' npomail ' ; nb. elem=' nelmail) ;
'MESSAGE' ('CHAINE' 'CPU Time (s)=' tcpus) ;
*
* Verification of energy balance
*
'MESSAGE' '-----------' ;
'MESSAGE' 'Energy balance :' ;
*  Power in the fuel
pcb = tabres . 'PCB' ;
chvolcb = 'DOMA' $mtcb 'VOLUME' ;
pdegcb = 'MAXIMUM' ('RESULT' ('*' chvolcb pcb)) ;
*  Integrated flux on the channel boundaries
f0 = tabres . 'F0' ;
pdegif = 'MAXIMUM' ('RESULT' f0) ;
errpif = '*' (ERRREL pdegif pdegcb) 100.D0 ;
*  Power gained by helium..
the = tabres . 'THE' ;
qhe = tabres . 'QHE' ;
*    .. in each channels
nchan = 'DIME' (tabin . 'CHANNEL') ;
pgcan = 0.D0 ;
'REPETER' ichan nchan ;
   po1    = tabgeo . 'tpechan' . &ichan ;
   po2    = tabgeo . 'tpschan' . &ichan ;
   tent = 'EXTRAIRE' the 'SCAL' po1 ;
   tsor = 'EXTRAIRE' the 'SCAL' po2 ;
   qent = 'EXTRAIRE' qhe 'SCAL' po1 ;
   qsor = 'EXTRAIRE' qhe 'SCAL' po2 ;
   'SI' ('NEG' qent qsor) ;
      cherr = 'CHAINE' 'Débit pas constant dans le canal ' &ichan
                       ' ???' ;
      'ERREUR' cherr ;
   'FINSI' ;                             
* cp is supposed to be constant; this is ok.
* However we suppose h = rho cp T, which is dubious
* we should have a procedure to compute h by integration
   selem1 = 'ELEM' (tabin . 'CHANNEL' . &ichan) 1 ;
   she  = 'MESURE' selem1 ;
   rhe muhe lhe cphe dum = PRPHE prhe tent ;
*   qhe = tabres . 'QHE' ;
*   pdeghe = qhe '*' she '*' cphe '*' ('-' tsor tent) ;
   pdeghe = qent '*' she '*' cphe '*' ('-' tsor tent) ;
   pgcan = '+' pgcan pdeghe ;
'FIN' ichan ;
pdeghe = pgcan ;
errphe = '*' (ERRREL pdeghe pdegcb) 100.D0 ;
'MESSAGE' ('CHAINE' 'Fuel power                 = ' (formar pdegcb 2)
          ' W') ;
'MESSAGE' ('CHAINE' 'Flux on channel boundaries = ' (formar pdegif 2)
          ' W') ;
'MESSAGE' ('CHAINE' '                    errpif = ' (formar errpif 2)
                                          ' %') ;
'MESSAGE' ('CHAINE' 'Power gained by Helium     = ' (formar pdeghe 2)
          'W') ;
'MESSAGE' ('CHAINE' '                    errphe = ' (formar errphe 2)
                                          ' %') ;
*  % error (i.e. limit to 0.1% = 1.D-4 error)
testene1 = > errpif 1.D0 ;
testene2 = > errphe 1.D0 ;
'SI' ('OU' testene1 testene2) ;
   'MESSAGE' '!!Pb. energy conservation' ;
'FINSI' ;
testerr = testerr 'OU' testene1 'OU' testene2 ;
'SI' testerr ;
   cherr = 'CHAINE' 'Error on the energy conservation !' ;
   'ERREUR' cherr ;
'FINSI' ;
* MinMax Helium Exit Temperature Speed
'MESSAGE' '-----------' ;
'MESSAGE' 'Helium Results :' ;
the = tabres . 'THE' ;
uhe = 'EXCO' 'UZ' (tabres . 'UHE') ;
miu = 'MINIMUM' uhe ;
mau = 'MAXIMUM' uhe ;
mit = 'MINIMUM' the ;
mat = 'MAXIMUM' the ;
'MESSAGE' ('CHAINE' 'Min. He speed  = ' (formar miu 2) ' m.s-1') ;
'MESSAGE' ('CHAINE' 'Max. He speed  = ' (formar mau 2) ' m.s-1') ;
'MESSAGE' ('CHAINE' 'Min. He temp.  = ' (formar ('-' mit 273.15D0) 2)
              ' degrees C') ;
'MESSAGE' ('CHAINE' 'Max. He temp.  = ' (formar ('-' mat 273.15D0) 2)
              ' degrees C') ;
* MinMax Fuel temperature + locations
'MESSAGE' '-----------' ;
'MESSAGE' 'Fuel Results :' ;
tcb = tabres . 'TCBS' ;
nmaxloc = 4 ;
tcbmi = 'MINIMUM' tcb ;
tcbmi = '*' tcbmi 1.0001D0 ;
mtcbmi = tcb 'POIN' 'INFERIEUR' tcbmi ;
nmtcbmi = 'NBNO' mtcbmi ;
'MESSAGE' ('CHAINE' 'Min fuel temp. = ' (formar ('-' tcbmi 273.15D0) 2)
                    ' degres C') ;
'MESSAGE' ('CHAINE' 'Location(m)'/1 '  x'/15  '  y'/30   '  z'/45) ;
'REPETER' imaxloc ('MINIMUM' ('LECT' nmaxloc nmtcbmi)) ;
   poi = 'POIN' mtcbmi &imaxloc ;
   xp yp zp = 'COORDONNEE' poi ;
   cxp = formar xp 2 ;
   cyp = formar yp 2 ;
   czp = formar zp 2 ;
   'MESSAGE' ('CHAINE' cxp*20  cyp*35 czp*50) ;
'FIN' imaxloc ;
'SI' ('>' nmtcbmi nmaxloc) ;
   'MESSAGE' ('CHAINE' 'etc...'/20) ;
'FINSI' ;   
*
nmaxloc = 4 ;
tcbma = 'MAXIMUM' tcb ;
tcbma = '*' tcbma 0.9999D0 ;
mtcbma = tcb 'POIN' 'SUPERIEUR' tcbma ;
nmtcbma = 'NBNO' mtcbma ;
'MESSAGE' ('CHAINE' 'Max fuel temp. = ' (formar ('-' tcbma 273.15D0) 2)
                    ' degrees C') ;
'MESSAGE' ('CHAINE' 'Location(m)'/1 '  x'/15  '  y'/30   '  z'/45) ;
'REPETER' imaxloc ('MINIMUM' ('LECT' nmaxloc nmtcbma)) ;
   poi = 'POIN' mtcbma &imaxloc ;
   xp yp zp = 'COORDONNEE' poi ;
   cxp = formar xp 2 ;
   cyp = formar yp 2 ;
   czp = formar zp 2 ;
   'MESSAGE' ('CHAINE' cxp*20  cyp*35 czp*50) ;
'FIN' imaxloc ;
'SI' ('>' nmtcbma nmaxloc) ;
   'MESSAGE' ('CHAINE' 'etc...'/20) ;
'FINSI' ;   
'MESSAGE' '**************' ;
*
* End of procedure file SUMMRES
*
'FINPROC' ;
*ENDPROCEDUR summres
*BEGINPROCEDUR sytopl
************************************************************************
* NOM         : SYTOPL
* DESCRIPTION : eventual symmetry + rotation + translation
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 23/11/2004, version initiale
* HISTORIQUE : v1, 23/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' SYTOPL ;
'ARGUMENT' mail*'MAILLAGE' ;
'ARGUMENT' arot*'FLOTTANT' ;
'ARGUMENT' vtra*'POINT' ;
'ARGUMENT' lsym*'LOGIQUE' ;
p0 = 0. 0. 0. ;
p1 = 0. 1. 0. ;
pr = 0. 0. 1. ;
'SI' lsym ;
   mail = 'SYMETRIE' mail 'DROIT' p0 p1 ;
'FINSI' ;
mail3 = 'TOURNER' mail arot p0 pr ;
mail2 = 'PLUS' mail3 vtra ;
'RESPRO' mail2 ;
*
* End of procedure file SYTOPL
*
'FINPROC' ;
*ENDPROCEDUR sytopl
*BEGINPROCEDUR tracchml
************************************************************************
* NOM         : TRACCHML
* DESCRIPTION : Trace un CHPOINT défini sur les centres des éléments
*               (valeur constante par élément) avec titre optionnel
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 12/03/2003, version initiale
* HISTORIQUE : v1, 12/03/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' TRACCHML ;
'ARGUMENT' chml*'CHPOINT ' ;
'ARGUMENT' $mt*'MMODEL  ' ;
'ARGUMENT' tit*'MOT     ' ;
'ARGUMENT' mt2/'MAILLAGE' ;
'ARGUMENT' foeil/'POINT   ' ;
'ARGUMENT' lnclk/'LOGIQUE' ;
*
rescal = 'KCHA' $mt chml 'CHAM' ;
mt = 'DOMA' $mt 'MAILLAGE' ;
'SI' ('NON' ('EXISTE' foeil)) ;
   foeil = (2.D0 -3.D0 0.7D0) ;
'FINSI' ;
oeil = CALCOEIL mt foeil ;
'SI' ('NON' ('EXISTE' lnclk)) ;
   lnclk = FAUX ;
'FINSI' ;   
modbid = 'MODELISER' ('DOMA' $mt 'MAILLAGE') 'THERMIQUE' ;
'SI' ('EXISTE' mt2) ;
   'SI' lnclk ;
      'TRACER' oeil 'CACH' rescal modbid mt2 'TITR' tit nbisov 'NCLK' ;
   'SINON' ;
      'TRACER' oeil 'CACH' rescal modbid mt2 'TITR' tit nbisov ;
   'FINSI' ;      
'SINON' ;
   'SI' lnclk ;
      'TRACER' oeil 'CACH' rescal modbid mt 'TITR' tit nbisov 'NCLK' ;
   'SINON' ;      
      'TRACER' oeil 'CACH' rescal modbid mt 'TITR' tit nbisov ;
   'FINSI' ;      
'FINSI' ;
*
* End of procedure file TRACCHML
*
'FINPROC' ;
*ENDPROCEDUR tracchml
*BEGINPROCEDUR tracchpo
************************************************************************
* NOM         : TRACCHPO
* DESCRIPTION : Tracé d'un chpoint avec titre optionnel.
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 14/10/2002, version initiale
* HISTORIQUE : v1, 14/10/2002, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' TRACCHPO ;
'ARGUMENT' pn*'CHPOINT ' ;
'ARGUMENT' mt*'MAILLAGE' ;
'ARGUMENT' mt2/'MAILLAGE' ;
'ARGUMENT' tit*'MOT     ' ;
'ARGUMENT' foeil/'POINT   ' ;
'ARGUMENT' lnclk/'LOGIQUE' ;
rescal = pn ;
'SI' ('NON' ('EXISTE' foeil)) ;
   foeil = (2.D0 -3.D0 0.7D0) ;
'FINSI' ;
oeil = CALCOEIL mt foeil ;
'SI' ('NON' ('EXISTE' lnclk)) ;
   lnclk = FAUX ;
'FINSI' ;   
'SI' ('EXISTE' mt2) ;
   'SI' lnclk ;
      'TRACER' oeil 'CACH' rescal mt mt2 'TITR' tit nbisov 'NCLK' ;
   'SINON' ;      
      'TRACER' oeil 'CACH' rescal mt mt2 'TITR' tit nbisov ;
   'FINSI' ;      
'SINON' ;
   'SI' lnclk ;
      'TRACER' oeil 'CACH' rescal mt mt 'TITR' tit nbisov 'NCLK' ;
   'SINON' ;      
      'TRACER' oeil 'CACH' rescal mt mt 'TITR' tit nbisov ;
   'FINSI' ;      
'FINSI' ;
*
* End of procedure file TRACCHPO
*
'FINPROC' ;
*ENDPROCEDUR tracchpo
*BEGINPROCEDUR tracmesh
************************************************************************
* NOM         : TRACMESH
* DESCRIPTION : Tracé d'un maillage
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 16/11/2004, version initiale
* HISTORIQUE : v1, 16/11/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' TRACMESH ;
'ARGUMENT' mt*'MAILLAGE' ;
'ARGUMENT' tit*'MOT' ;
'ARGUMENT' foeil/'POINT   ' ;
'ARGUMENT' lnclk/'LOGIQUE' ;
*
'SI' ('NON' ('EXISTE' foeil)) ;
   foeil = (2.D0 -3.D0 0.7D0) ;
'FINSI' ;
oeil = CALCOEIL mt foeil ;
'SI' ('NON' ('EXISTE' lnclk)) ;
   lnclk = FAUX ;
'FINSI' ;   
'SI' lnclk ;
   'TRACER' oeil mt 'CACH' 'TITRE' tit 'NCLK' ;
'SINON' ;   
   'TRACER' oeil mt 'CACH' 'TITRE' tit ;
'FINSI' ;   
*
* End of procedure file TRACMESH
*
'FINPROC' ;
*
* 
*
*ENDPROCEDUR tracmesh
'OPTION' 'ECHO' 1 ;
************************************************************************
* NOM         : COUPLAGE_TH1D_TH3D
* DESCRIPTION : 
* Idem COUPLAGE_TH1D_TH3D_x, (x=1 or) but better :
* we try to write an easy to use dgibi file with well
* specified inputs and outputs.
*
* See technical reports DM2S/SFME/LTMF/RT/03-012/A
*                       DM2S/SFME/LTMF/RT/03-012/B
*                       DM2S/SFME/LTMF/RT/04-024/A
*                       DM2S/SFME/LTMF/RT/05-050/A
* Coupled modelization 3D thermic-1D thermohydraulics
* of a gas-cooled reactor.
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 09/10/2004, version initiale
* HISTORIQUE : v1, 09/10/2004, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
*
* Some global options and settings
*
* 
* complet  = VRAI : runs all the cases
* interact = VRAI : interactive use
* graph    = VRAI : if one wants graphic output (TRAC, DESS)
complet = FAUX ;
interact= FAUX ;
graph   = FAUX ;
*
'OPTION' 'DIME' 3 'ELEM' 'CUB8' ;
'OPTION' 'ISOV' 'SULI' ;
* nbisov : number of isovalues 
nbisov = 25 ;
'SI' ('NON' interact) ;
  'OPTION' 'TRAC' 'PSC' ;
   'OPTION' 'ECHO' 0 ;
'SINON' ;
  'OPTION' 'TRAC' 'X' ;
*  'OPTION' 'TRAC' 'OPEN' ;
'FINSI' ;  
***************************************
*
* Geometry-related settings
*
***************************************
*
* Geometry data for hexagonal configuration
* (quantities to be specified)
*  pas : dimension of the unit cell (m)
*  Dh     : Hydraulique diameter (m)
*  H      : Core height (m)   
* Geometry data for plate configuration
* (quantities to be specified)
*  Dh     : Hydraulique diameter (m)
*  ef     : width of the structure part (m)
*  ep     : width of a fuel plate (m) ;
*  nchan  : number of channels
*  H      : Core height (m)
*
* Geometry data for hexagonal plate configuration
*
* ept : Entreplat (m)
* eps : Epaisseur de structure (m)
* epp : Epaisseur de plaque (m)
* epg : Epaisseur lame de gaz (m)
* epge : Epaisseur lame de gaz extrême (m)
* npp : Nombre de plaques (m)
*
*
*  Geometric discretization parameters 
*  (number of elements)
*  nh (height) : number of elements in the height
*
***************************************
*
* Physics-related settings
*
***************************************
*
* Helium
*
* prhe : Helium pressure (bar)
* qhe  : Helium flow rate (bottom of the core) (kg.m-2.s-1)
* theec, theek : Cold Helium temperature in degree C and in K
* fblock : relative residual flow rate in blocked channels (between 0.
*          and 1.)
*
* Fuel
*
* pvol : Volumic power
* It is the power per unit of total core (fuel+helium+structure)
* and not the power per unit of fuel (W.m-3)
* gamma : Gamma factor (power peaking factor)
*         1.3 < gamma < 2.0
* pvar     = VRAI : sinusoidal power profile in the fuel
*          = FAUX : constant power profile in the fuel
* irrad    = VRAI : for used nuclear fuel (new otherwise)
*
***************************************
*
* Solver-related settings
*
***************************************
*
* Solver parameters
*
* maxit  : maximum number of non-linear iterations
* rescvg : convergence occur for a temperature residual < rescvg
* solvit = FAUX : direct solver, more robust, slower
*        = VRAI : iterative solver, less robust, faster
* omcbst : relaxation factor for 3D thermics
* omhe   : relaxation factor for 1D thermo-hydraulics
*
*
*************************************
*
* Definition of the test cases :
****
*   ikas : 1 to 3  Technical report DM2S/SFME/LTMF/RT/03-012/B
*   ikas = 1 : Hexagonal core, nominal case, small  symmetry cell
*   ikas = 2 : Hexagonal core, nominal case, larger symmetry cell
*   ikas = 3 : Hexagonal core, accidental case, larger symmetry cell
****
*   ikas : 4 to 6  Technical report DM2S/SFME/LTMF/RT/04-024/A
*   ikas = 4 : Square core, nominal case, small  symmetry cell
*   ikas = 5 : Square core, nominal case, larger symmetry cell
*   ikas = 6 : Square core, accidental case, larger symmetry cell
****
*   ikas : 7 to    Technical report DM2S/SFME/LTMF/RT/05-050/A
*   ikas = 7 : Hexagonal core, nominal case, small losange cell
*   ikas = 8 : Hexagonal core, channel 2 blocked, small losange cell
*      chapter 2
*   ikas = 9  : Hexagonal core, nominal case, normal hexa cell
*   ikas = 10 : Hexagonal core, channel 1 blocked, normal hexa cell
*   ikas = 11 : Hexagonal core, channel 2 blocked, normal hexa cell
*   ikas = 12 : Hexagonal core, channel 5 blocked, normal hexa cell
*   ikas = 13 : Hexagonal core, channel 4 and 6 blocked, normal hexa cell
*   ikas = 14 : Hexagonal core, channel 5 and 6 blocked, normal hexa
*               cell (doesn't work)
*   ikas = 15 : Hexagonal core, channel 5 and 6 95% blocked, normal hexa
*               cell
*   ikas = 16 : Hexagonal core, channel 2 blocked, large 2 hexa cells
*    chapter 4
*   ikas = 17 : Hexagonal core, channel 2 blocked, normal hexa cell
*               Reduced Pvol and Qhe : * 75 %
*   ikas = 18 : Hexagonal core, channel 2 blocked, normal hexa cell
*               Reduced Pvol and Qhe : * 40 %
*   ikas = 19 : Hexagonal core, channel 2 blocked, normal hexa cell
*               Qhe Increase : + 50 %
*   ikas = 20 : Hexagonal core, channel 2 blocked, normal hexa cell
*               ((Pvol Qhe)*75%) + Qhe Increase : + 25 %
*    chapter 3
*   ikas = 30 : Modified hexagonal core, nominal case, normal hexa cell
*   ikas = 31 : Modified hexagonal core, channel 1 blocked, normal hexa cell 
*   ikas = 32 : Modified hexagonal core, channel 2 blocked, normal hexa cell
*   ikas = 33 : Modified hexagonal core, channel 5 blocked, normal hexa cell
*   ikas = 34 : Modified hexagonal core, channel 4 and 6 blocked,
*               normal hexa cell
*   ikas = 35 : Modified hexagonal core, channel 5 and 6 blocked, normal hexa
*               cell (doesn't work)
*   ikas = 36 : Modified hexagonal core, channel 5 and 6 95% blocked,
*               normal hexa cell
 
 
*
* Tabtot stores tabres + maxi temperature
*
tabtot = 'TABLE' ;
* Non-regression testing on the maximum temperature in fuel and
* structure in listem 
* Obtained on 01/12/05
listem = 'PROG'
  1154. 1154. 1595. 1101. 1101. 1865. 1080. 1910. 1110. 1700.
  1840. 1700. 1810. 0.    4470. 1820. 1690. 1500. 1580. 1550.
  0.    0.    0.    0.    0.    0.    0.    0.    0.    1070.
  1760. 1760. 1760. 1870. 0.    4680. ;
'SI' complet ;
   lkas   = 'LECT' 1 2 3 4 5 6 7 8 9 10 11 12 13 15
                   16 17 18 19 20 30 31 32 33 34 36 ;
'SINON' ;   
   lkas   = 'LECT' 1 2 3 4 ;
'FINSI' ;
listemc = 'EXTRAIRE' listem lkas ;
nkas = 'DIME' lkas ;
'REPETER' bikas nkas ;
ikas = 'EXTRAIRE' lkas &bikas ;
tabtot . ikas = 'TABLE' ;
'TEMPS' 'ZERO' ;
petitb = 0.D0 ;
omcbst = 0.9 ;
omhe   = 0.5 ;
*
'SI' ('ET' ('>EG' ikas 1)  ('&lt;EG' ikas 3)) ;
   pas  = 21.1D-3 ;
   Dh   = 16.5D-3 ;
   H    = 2.D0 ;
   nh = 7 ;
   prhe = 90.D0 ;
   qhe  = 210.D0 ;
   theec = 480.D0 ;
   theek = '+' theec 273.15D0 ;
   pvol = 100.D6  ;
   gamma = 1.5D0 ;
   pvar    = VRAI ;
   irrad   = VRAI ;
   maxit  = 40 ;
   rescvg = 1.D-4 ;
   solvit = FAUX ;
   'SI' ('EGA' ikas 1) ;   
      tabin = GEOMHEX1 Dh pas ;
      tabin . 'BLOCKED'   = 'LECT' ;
   'FINSI' ;      
   'SI' ('EGA' ikas 2) ;   
      tabin = GEOMHEX2 Dh pas ;
      tabin . 'BLOCKED'   = 'LECT' ;
   'FINSI' ;      
   'SI' ('EGA' ikas 3) ;   
      tabin = GEOMHEX2 Dh pas ;
      tabin . 'BLOCKED'   = 'LECT' 1 ;
   'FINSI' ;
   nchan = 'DIME' (tabin . 'CHANNEL') ;
   tabin . 'DH' = 'TABLE' ;
   'REPETER' ichan nchan ;
      tabin . 'DH' . &ichan = Dh ;
   'FIN' ichan ;      
   fblock = petitb ;
'FINSI' ;   
'SI' ('ET' ('>EG' ikas 4)  ('&lt;EG' ikas 6)) ;
   Dh = 2.D-2 ;
   ef = 3.04D-3 ;
   ep = 1.04D-2 ;
   nchan = 5 ;
   H  = 4.2D0 ;
   nh = 7 ;
   prhe = 90.D0 ;
   qhe  = 278.D0 ;
   theec = 480.D0 ;
   theek = '+' theec 273.15D0 ;
   pvol = 40.D6  ;
   gamma = 1.3D0 ;
   pvar    = VRAI ;
   irrad   = VRAI ;
   maxit  = 40 ;
   rescvg = 1.D-4 ;
   solvit = VRAI ;
   'SI' ('EGA' ikas 4) ;   
      tabin = GEOMPLA1 Dh ef ep nchan ;
      tabin . 'BLOCKED'   = 'LECT' ;
   'FINSI' ;      
   'SI' ('EGA' ikas 5) ;   
      tabin = GEOMPLA2 Dh ef ep nchan ;
      tabin . 'BLOCKED'   = 'LECT' ;
   'FINSI' ;      
   'SI' ('EGA' ikas 6) ;
      tabin = GEOMPLA2 Dh ef ep nchan ;
      tabin . 'BLOCKED'   = 'LECT' 1 ;
   'FINSI' ;
   nchan = 'DIME' (tabin . 'CHANNEL') ;
   tabin . 'DH' = 'TABLE' ;
   'REPETER' ichan nchan ;
      tabin . 'DH' . &ichan = Dh ;
   'FIN' ichan ;      
   fblock = petitb ;
'FINSI' ;
'SI' ('>EG' ikas 7) ;
   H = 1.55D0 ;
   nh = 7 ;
   prhe = 90.D0 ;
   qhe  = 210.D0 ;
   theec = 480.D0 ;
   theek = '+' theec 273.15D0 ;
   pvol = 100.D6  ;
   gamma = 1.5D0 ;
   pvar    = VRAI ;
   irrad   = VRAI ;
   maxit  = 40 ;
   rescvg = 1.D-4 ;
   solvit = VRAI ;
   ept = 211.3D-3 ;
   eps = 1.5D-3 ;
   epp = 7.D-3  ;
   epg = 4.6D-3 ;
   npp  = 9 ;
   ff = '/' 2.D0 ('**' 3.D0 0.5D0) ;
   'SI' ('>EG' ikas 30) ;
      ept = '+' ept ('*' ('-' 4.6D-3 ('/' 2.08D-3 ff)) 4.) ;
   'FINSI' ;      
   'SI' ('EGA' ikas 7) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 1 ;
      tabin . 'BLOCKED' = 'LECT' ;
      fblock = 0.D0 ;
   'FINSI' ;      
   'SI' ('EGA' ikas 8) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 1 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 9) ('EGA' ikas 30)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 10) ('EGA' ikas 31)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 1 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 11) ('EGA' ikas 32)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 12) ('EGA' ikas 33)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 5 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 13) ('EGA' ikas 34)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 4 6 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 14) ('EGA' ikas 35)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 5 6 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('OU' ('EGA' ikas 15) ('EGA' ikas 36)) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 5 6 ;
      fblock = 0.05D0 ;
   'FINSI' ;
   'SI' ('EGA' ikas 16) ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 3 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('EGA' ikas 17) ;
      fq = 1.D0 ;
      fp = 0.75D0 ;
      pvol = '*' pvol fp ;
      qhe  = '*' ('*' qhe fp) fq ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('EGA' ikas 18) ;
      fq = 1.D0 ;
      fp = 0.4D0 ;
      pvol = '*' pvol fp ;
      qhe  = '*' ('*' qhe fp) fq ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('EGA' ikas 19) ;
      fq = 1.5D0 ;
      fp = 1.0D0 ;
      pvol = '*' pvol fp ;
      qhe  = '*' ('*' qhe fp) fq ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
   'SI' ('EGA' ikas 20) ;
      fq = 1.25D0 ;
      fp = 0.75D0 ;
      pvol = '*' pvol fp ;
      qhe  = '*' ('*' qhe fp) fq ;
      tabin epge = GEOMPLH1 ept eps epp epg npp 2 ;
      tabin . 'BLOCKED' = 'LECT' 2 ;
      fblock = 0.D0 ;
   'FINSI' ;
*   
   'SAUTER' 2 'LIGN' ;
   'MESSAGE' ('CHAINE' '***************')  ;
   'MESSAGE' ('CHAINE' 'Epai. gaz extreme : ' (formar ('/' epge ff) 2)
                       ' m') ;
   'MESSAGE' ('CHAINE' 'Entreplat         : ' (formar ept 5)
                       ' m') ;
   'MESSAGE' ('CHAINE' '***************')  ;
   'SAUTER' 2 'LIGN' ;
'FINSI' ;
* Input for GENMEMA procedure
* tabin . 'INTERFACE'  . i : i-th interface       (line mesh)
* tabin . 'CHANNEL'    . i : i-th helium channel  (surface mesh)
* tabin . 'FUEL'           : fuel mesh (surface mesh)
* tabin . 'STRUCTURE'      : structure mesh (surface mesh) (optional)
* tabin . 'BLOCKED'        : LISTENTI : list of the numbers
*                            of the blocked channels (can be empty)
* H : height (m) and nh : number of elements in the height
* (regular meshing in z direction)
*
* output tabgeo should contain all the necessary geometric information
* needed further
* tabgeo . 'tabin'   = tabin
* tabgeo . 'H'       = H   : height
* tabgeo . 'nh'      = nh  : number of elements in the height
* tabgeo . '$mtcbst' = NS model objet (fuel+structure)
* tabgeo . '$mtcb'   = NS model objet (fuel only)
* tabgeo . '$mthe'   = NS model objet (helium channels)
* tabgeo . '$sthe'   = NS model objet (bottom of the helium channels)
* tabgeo . '$mtif'   = NS model objet (interfaces)
* tabgeo . '$mtifb'  = NS model objet (blocked interfaces) (if any)
* tabgeo . 'mcf'     = coupling matrix helium -> fuel
* tabgeo . 'mfc'     = coupling matrix fuel -> helium
* GENMEMA is slow
tabgeo = GENMEMA tabin H nh ;
*
* Compute power distribution
*
pcb = CPOWDIS tabgeo pvol gamma pvar ;
*
* chblock
* Volumic field = 1 in unblocked channels
*               = fblock (between 0. and 1.) in blocked channels
*
$tmthe  = tabgeo . '$tmthe' ;
nhe     = 'DIME' $tmthe ;
lblock  = tabin . 'BLOCKED' ;
tablock = 'TABLE' ;
'REPETER' ihe nhe ;
   $mthei = $tmthe . &ihe ;
   'SI' ('DANS' ('LECT' &ihe) lblock) ;
      APPEND tablock 'chblock' ('KCHT' $mthei 'SCAL' 'SOMMET' fblock) ;
   'SINON' ;      
      APPEND tablock 'chblock' ('KCHT' $mthei 'SCAL' 'SOMMET' 1.D0) ;
   'FINSI' ;
'FIN' ihe ;
chblock = 'KCHT' (tabgeo . '$mthe') 'SCAL' 'SOMMET'
                                    (tablock . 'chblock') ;
*
* Put the physics in a table 
*
tabphy = 'TABLE' ;
* Interface
tabphy . 'Dh'      = tabgeo . 'Dh' ;
tabphy . 'chblock' = chblock    ;
* Helium
tabphy . 'qhe'   = qhe    ;
tabphy . 'prhe'  = prhe   ;
tabphy . 'theek' = theek  ;
* Fuel 
tabphy . 'pvol'  = pvol   ;
tabphy . 'pvar'  = pvar   ;
tabphy . 'gamma' = gamma  ;
tabphy . 'irrad' = irrad  ;
tabphy . 'pcb'   = pcb    ;
*
* Put the numerics in a table
*
tabnum = 'TABLE' ;
tabnum . 'maxit'  = maxit  ;
tabnum . 'rescvg' = rescvg ;
tabnum . 'solvit' = solvit ;
tabnum . 'omcbst' = omcbst ;
tabnum . 'omhe'   = omhe ;
*
* Put the tables in a table !
*
tabdat = 'TABLE' ;
tabdat . 'ikas' = ikas ;
tabdat . 'tabgeo' = tabgeo ;
tabdat . 'tabphy' = tabphy ;
tabdat . 'tabnum' = tabnum ;
*
* Text summary of data
*
'SAUTER' 2 'LIGN' ;
SUMMDAT tabdat ;
'SAUTER' 2 'LIGN' ;
*
* EXPLODAT does graphical output of data
* - first argument (mandatory) : the tables with data
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                  an integer  
*              an integer corresponds to a particular view
*   View 1 : the base mesh
*   View 2 : the mesh
*   View 3 : the blocked channels
*   View 4 : the power distribution
*   View 5 : fuel and structure thermal conductivity (T)
*   View 6 : Helium density
*   View 7 : Helium thermal conductivity
*   View 8 : Helium heat capacity
*
* Graphical output of data
*
'SI' graph ;
   'SI' interact ;
     EXPLODAT tabdat ;
   'SINON' ;   
*      EXPLODAT tabdat FAUX ('LECT' 1 2 3 4) ;
*      EXPLODAT tabdat FAUX ('LECT' 1) ;
   'FINSI' ;
'FINSI' ;
***************************************
*
* DO THE COMPUTATIONS
*
***************************************
tabres = COMPUTE tabdat ;
 
TABTPS = TEMP 'NOEC';
tcpu = TABTPS.'TEMPS_CPU'.'INITIAL' ;
tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
'MESSAGE' ('CHAINE' 'tcpus=' tcpus) ;
tabres . 'tcpus' = tcpus ;
***************************************
*
* POST-TREATEMENT
*
***************************************
*
* Text summary of results :
* SUMMRES prints a summary of the results of the computation
*         it also does some tests (energy conservation...)
*
'SAUTER' 2 'LIGN' ;
SUMMRES tabres ;
'SAUTER' 2 'LIGN' ;
*
* EXPLORES does graphical output of results
* - first argument (mandatory) : the table with the results
* - second argument (optional) : VRAI (default) interactive use
*                                     (mouse clicks)
*                                FAUX batch use
* - third argument (mandatory if 2ndarg = FAUX) : a list of integer or
*                                                 an integer  
*              an integer corresponds to a particular view
*   View 1 : temperature (fuel + He)
*   View 2 : temperature (fuel)
*   View 3 : temperature (fuel : 2 viewpoints)
*   View 4 : exchange coefficient
*   View 5 : temperature (helium)
*   View 6 : vertical speed (helium)
*   View 7 : temperature profiles per channel (helium)
*   View 8 : vertical speed profiles per channel (helium)
*   View 9  : Nusselt number
*   View 10 : Reynolds number
*   View 11 : Prandtl number
*   View 12 : Lambda Helium
*   View 13 : Hydraulic diameter
'SI' graph ;
   'SI' interact ;
      EXPLORES tabres ;
   'SINON' ;   
*      EXPLORES tabres FAUX ('LECT' 2 4 5 7) ;
      EXPLORES tabres FAUX ('LECT' 3 4 5 6 7 8) ;
   'FINSI' ;
'FINSI' ;
*interact = VRAI ;
*explores tabres ;
*
* Save the results (Maximum temperature & cpu time)
*
tabtot . ikas . 'maxtem' = 'MAXIMUM' (tabres . 'TCBS') ;
tabtot . ikas . 'tcpus'  = tabres . 'tcpus' ;
*tabtot . ikas . 'tabdat' = tabdat ;
*tabtot . ikas . 'tabres' = tabres ;
'FIN' bikas ;
*'SI' interact ;
*   'OPTION' 'ECHO' 1 ;
*   'OPTION' 'DONN' 5 ;
*'FINSI' ;   
*'TEMPS' 'IMPR' ;
*
* Non-regression testing on the maximum temperature in fuel and structure.
*
'MESSAGE' '**********************' ;
'MESSAGE' 'Non regression testing' ;
testerr = FAUX ;
* 1 % tolerance
toler   = 0.01 ;
'REPETER' bikas nkas ;
   ikas = 'EXTRAIRE' lkas &bikas ;
   maxtem   = tabtot . ikas . 'maxtem' ;
   tcpus    = tabtot . ikas . 'tcpus' ;
   maxexpc  = 'EXTRAIRE' listemc &bikas ;
   maxexp = '+' maxexpc 273.15D0 ;
   testerr = 'OU' testerr ('NEG' maxtem maxexp ('*' maxexp toler)) ; 
   'MESSAGE' ('CHAINE' 'Ikas            = ' ikas) ;
   'MESSAGE' ('CHAINE' 'CPU Time (s)    =' tcpus) ;
   maxtemc = '-' maxtem 273.15D0 ;
   'MESSAGE' ('CHAINE' 'Maxtem          = ' (formar maxtemc 3)
              ' degrees C') ;
   'MESSAGE' ('CHAINE' 'Expected Maxtem = ' (formar maxexpc 3)
              ' degrees C') ;
'FIN' bikas ;
'MESSAGE' '**********************' ;
*!
*! Error checking
*!
'SI' testerr ;
   'MESSAGE' 'There were some errors' ;
   'SI' interact ;
      'OPTION' 'ECHO' 1 ;
      'OPTION' 'DONN' 5 ;
   'FINSI' ;         
   'ERREUR' 5 ;
'SINON' ;   
   'MESSAGE' 'Everything went ok' ;
'FINSI' ;
'SI' interact ;
   'OPTION' 'ECHO' 1 ;
   'OPTION' 'DONN' 5 ;
'FINSI' ;   
*
* End of dgibi file COUPLAGE_TH1D_TH3D
*
'FIN' ;