Télécharger la plaquette Cast3M


fr en

Télécharger satutils.procedur

Retour à la liste

Numérotation des lignes : 

  1. * SATUTILS  PROCEDUR  FANDEUR   13/03/25    21:15:00     7446           
  2. ************************************************************************
  3. * NOM         : SATUTILS
  4. * DESCRIPTION :
  5. *
  6. *
  7. *
  8. * LANGAGE     : GIBIANE-CAST3M
  9. * AUTEUR      : Gilles BERNARD - MICHEL (CEA/DEN/DM2S/SFME/LTMF)
  10. *               mél : gbm@semt2.smts.cea.fr
  11. **********************************************************************
  12. * VERSION    : v1, 25/12/2006, version initiale
  13. * HISTORIQUE : v1, 25/12/2006, création
  14. * HISTORIQUE :
  15. * HISTORIQUE :
  16. ************************************************************************
  17. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  18. * en cas de modification de ce sous-programme afin de faciliter
  19. * la maintenance !
  20. ************************************************************************
  21. *
  22. 'DEBPROC' SATUTILS                                                    ;
  23. 'ARGUMENT' roro*'MOT'                                                 ;
  24.  
  25. *************************************************************************
  26.  
  27. 'SI' ('EGA' roro LICODINI)                                            ;
  28.  
  29. 'ARGUMENT'  SATUR*'TABLE'                                             ;
  30.  
  31. * GBM INITIALISE CORRECT TRACE ET FLUX
  32.  
  33. * ----------------------------------------------------------------- *
  34. *  RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS CALCULE  *
  35. * ----------------------------------------------------------------- *
  36. *
  37. *                                                                TEMPS
  38. 'SI' ( 'NON' ('EXISTE' SATUR 'TEMPS' ) )                              ;
  39.      SATUR.'TEMPS' = 'TABLE'                                          ;
  40.      SATUR.'TEMPS'. 0 = 0.                                            ;
  41. 'FINSI'                                                               ;
  42. *                                                               CHARGE
  43. 'SI' ( 'NON' ('EXISTE' SATUR 'CHARGE' ) )                             ;
  44.      'ERREUR'
  45.      'Indice CHARGE absent de la table de données.'                   ;
  46.      'QUITTER' DARCYSAT                                               ;
  47. 'FINSI'                                                               ;
  48. *                                                                 FLUX
  49. * La présence de l'indice FLUX est obligatoire.
  50. *                                          Tests des tailles de tables
  51. IND1 = 'INDEX'  ( SATUR . 'TEMPS' )                                   ;
  52. IND3 = 'INDEX'  ( SATUR . 'CHARGE' )                                  ;
  53. LIN1 = 'DIME' IND1                                                    ;
  54. LIN3 = 'DIME' IND3                                                    ;
  55. 'SI' ( LIN1 'NEG' LIN3 )                                              ;
  56.      'ERREUR'
  57.      'Longueur des tables TEMPS et CHARGE différente.'                ;
  58.      'QUITTER' DARCYSAT                                               ;
  59. 'FINSI'                                                               ;
  60. *                                          Tests des indices de tables
  61. IPO1 = 0                                                              ;
  62. 'REPETER' BOU1 LIN1                                                   ;
  63.    IPO1  = IPO1 + 1                                                   ;
  64.    LAST1 = IND1 . IPO1                                                ;
  65.    LAST3 = IND3 . IPO1                                                ;
  66.    'SI' ( 'NEG' LAST1 LAST3 )                                         ;
  67.         'ERREUR'
  68.         'Indices des tables TEMPS et CHARGE incohérents.'             ;
  69.         'QUITTER' DARCYSAT                                            ;
  70.    'FINSI'                                                            ;
  71. 'FIN' BOU1                                                            ;
  72. *                                 Récupération des inconnues initiales
  73. TPSINI = SATUR . 'TEMPS' . LAST1                                      ;
  74. CHRG   = SATUR . 'CHARGE' . LAST1                                     ;
  75. FLU0 = SATUR . 'FLUX' . LAST1                                         ;
  76. 'SI' (EXISTE SATUR  'TRACE_CHARGE')                                   ;
  77. *GBM METTRE A JOUR INITIALISATION - TRES IMPORTANT
  78.    SI (EXISTE SATUR . 'TRACE_CHARGE' LAST1)                           ;
  79.       TRH   = SATUR . 'TRACE_CHARGE' . LAST1                          ;
  80.    SINON                                                              ;
  81.       TRH    = 0.D0 * SATUR . 'FLUX' . LAST1                          ;
  82.       SATUR . 'TRACE_CHARGE' . LAST1 = TRH                            ;
  83.    FINSI                                                              ;
  84. 'SINON'                                                               ;
  85.    TRH    = 0.D0 * SATUR . 'FLUX' . LAST1                             ;
  86.    SATUR . 'TRACE_CHARGE' = table                                     ;
  87.    SATUR . 'TRACE_CHARGE' . LAST1 = TRH                               ;
  88. 'FINSI'                                                               ;
  89.  
  90. 'RESPRO'  TPSINI CHRG TRH FLU0 LAST1                                  ;
  91. 'FINSI'                                                               ;
  92.  
  93. ***********************************************************************
  94.  
  95. 'SI' ('EGA' roro LIMODELE)                                            ;
  96.  
  97. 'ARGUMENT' SATUR*'TABLE'                                              ;
  98.  
  99. *
  100. 'SI' ( 'EXISTE' SATUR 'MODELE' )                                      ;
  101.    MDMH = 'EXTR' (SATUR . 'MODELE') 'FORMULATION' 'DARCY'             ;
  102. 'SINON'                                                               ;
  103.    'ERREUR' 'Le modèle doit contenir une formulation DARCY'           ;
  104.    'QUITTER' DARCYSAT                                                 ;
  105. 'FINSI'                                                               ;
  106. MCHYB    = 'DOMA' MDMH 'ORIENTAT'                                     ;
  107. VHYB     = 'DOMA' MDMH 'VOLUME'                                       ;
  108. MAILSOM  = 'DOMA' MDMH 'SOMMET'                                       ;
  109. MAILCENT = 'DOMA' MDMH 'CENTRE'                                       ;
  110. ABMC     = 'ABS'  MCHYB                                               ;
  111.  
  112. 'RESPRO'  MDMH MCHYB VHYB MAILSOM MAILCENT ABMC                       ;
  113. 'FINSI'                                                               ;
  114.  
  115. ***********************************************************************
  116.  
  117. 'SI' ('EGA' roro LIPHYSIK)                                            ;
  118.  
  119. 'ARGUMENT'  SATUR*'TABLE'                                             ;
  120.  
  121. *                                                  COEF_EMMAGASINEMENT
  122. 'SI' ( 'EXISTE' SATUR 'COEF_EMMAGASINEMENT')                          ;
  123.    COFEMMAG = SATUR.'COEF_EMMAGASINEMENT'                             ;
  124.    'SI' ( ('NEG' ('TYPE' COFEMMAG) 'CHPOINT ') 'ET'
  125.           ('NEG' ('TYPE' COFEMMAG) 'FLOTTANT') )                      ;
  126.        'ERREUR' ('CHAINE' 'Le coefficient d emmagasinement doit etre'
  127.                           ' de type CHPOINT ou FLOTTANT' )            ;
  128.    'FINSI' ;
  129.    'SI' ('EGA' ('TYPE' COFEMMAG) 'FLOTTANT')                          ;
  130.       COFEMMAG = 'MANU' 'CHPO' MAILCENT 1 'SCAL' COFEMMAG
  131.                         'NATURE' 'DISCRET'                            ;
  132.    'FINSI'                                                            ;
  133. 'SINON'                                                               ;
  134.    COFEMMAG = 'MANU' 'CHPO' MAILCENT 1 'SCAL' 0.D0 'NATURE' 'DISCRET' ;
  135. 'FINSI'                                                               ;
  136.  
  137. *                                                            PESANTEUR
  138. * Valeur de l'accélération de la pesanteur
  139. L_GRAV = 'EXISTE' SATUR 'FORCE_GRAVITE'                               ;
  140. 'SI' L_GRAV                                                           ;
  141.    RHOWG = SATUR . 'FORCE_GRAVITE'                                    ;
  142. 'SINON'                                                               ;
  143.    RHOWG = 0.D0                                                       ;
  144. 'FINSI'                                                               ;
  145.  
  146. 'SI' ( 'EXISTE' SATUR 'AXE_G')                                        ;
  147.    GAXE = SATUR.'AXE_G'                                               ;
  148. *  normalisation de GAXE
  149.    NORGAXE = 'NORM' GAXE                                              ;
  150.    'SI' (NORGAXE 'NEG' 0.)                                            ;
  151.       XAXE = ('COOR' 1 GAXE) / NORGAXE                                ;
  152.       YAXE = ('COOR' 2 GAXE) / NORGAXE                                ;
  153.       'SI' (NDIME 'EGA' 3)                                            ;
  154.          ZAXE = ('COOR' 3 GAXE) / NORGAXE                             ;
  155.       'FINSI'                                                         ;
  156.    'SINON'                                                            ;
  157. *     GBM PENSE DOIT ETRE FAUX - A VOIR
  158.       XAXE YAXE ZAXE = 0. 0. 0.                                       ;
  159.    'FINSI'                                                            ;
  160. 'SINON'                                                               ;
  161.    'SI' (NDIME 'EGA' 2)                                               ;
  162.       XAXE YAXE = 0. 1.                                               ;
  163.    'SINON'                                                            ;
  164.       XAXE YAXE ZAXE = 0. 0. 1.                                       ;
  165.    'FINSI'                                                            ;
  166. 'FINSI'                                                               ;
  167.  
  168. 'SI' (NDIME 'EGA' 3)                                                  ;
  169.          DAXE = (XAXE YAXE ZAXE)                                      ;
  170. 'SINON'                                                               ;
  171.          DAXE = (XAXE YAXE)                                           ;
  172. 'FINSI'                                                               ;
  173.  
  174. *                                                    CONVERSION_CHARGE
  175. * Facteur de convertion de l'unité d'expression de la charge vers des
  176. * Pascals. Par défaut, on suppose que les pressions sont en Pascal.
  177. 'SI' ('EXISTE' SATUR 'CONVERSION_CHARGE')                             ;
  178.    CONVH = SATUR.'CONVERSION_CHARGE'                                  ;
  179. 'SINON'                                                               ;
  180.    CONVH = 1.                                                         ;
  181. 'FINSI'                                                               ;
  182. 'RESPRO'  COFEMMAG L_GRAV RHOWG DAXE CONVH                            ;
  183.  
  184. 'FINSI'                                                               ;
  185.  
  186. **************************************************************************
  187.  
  188. * Procédure d'affichage
  189. 'SI' ('EGA' roro AFFICH)                                              ;
  190. 'ARGU'  IPAL*'ENTIER' ITRANSI*'ENTIER' LAST1*'ENTIER'
  191.           TPS*'FLOTTANT' TPSANC*'FLOTTANT' DELTAT*'FLOTTANT'
  192.           DTI*'FLOTTANT' debug*'LOGIQUE'                              ;
  193.  
  194.     'SI' debug                                                        ;
  195. *      Affichage de l'en-tête de messages, complété à chaque itération
  196.        'MESS' ' '                                                     ;
  197.        'SI' (('EGA' IPAL 1) 'ET' ('EGA' ITRANSI 1)) ;
  198.          'MESS' 'PAS :' (LAST1 + 1)  'temps (s) : ' TPS '=' TPSANC
  199.                 '+'  (DELTAT)                                         ;
  200.        'SINON'                                                        ;
  201. * rendre cet affichage optionnel
  202.          'MESS' 'PAS :' (LAST1 + 1)  'CFL : ' (DELTAT/DTI)
  203.                 'temps (s) : ' TPS '=' TPSANC '+' DELTAT              ;
  204.        'FINSI'                                                        ;
  205.  
  206.        'MESS' ('CHAINE' ' ---------------------------------------------'
  207.                    '------------------------------------------------') ;
  208.        'MESS' ('CHAINE' '|Picard| max(Résidu) |    min tp    |  '
  209.              '  max tp    |   max capa  |   min Perm  | faces désat.|');
  210.        'MESS' ('CHAINE' ' ---------------------------------------------'
  211.                    '------------------------------------------------') ;
  212.     'SINON'                                                           ;
  213. *      Affichage de temps en cours
  214.        mess '  DARCYSAT : PAS ' (LAST1 + 1) ', temps ' TPS '(s)'      ;
  215.     'FINSI'                                                           ;
  216. 'FINSI'                                                               ;
  217.  
  218.  
  219. *************************************************************************
  220.  
  221. 'SI' ('EGA' roro KALSAT)                                              ;
  222. 'ARGUMENT'  MDMH*'MMODEL' LOIS*'TABLE' NIVEAU*'MOT'
  223.          PNS*'CHPOINT';
  224.  
  225. * GBM  REGARDER SI PAS MEILLEURE BOUCLE QUAND TOUTES MEMES LOIS MAIS UN
  226. *      COEF PAR ZONE
  227. mailcent = 'DOMA' MDMH CENTRE                                         ;
  228. *- calcul teneur en eau, saturation et premier terme capacité capillaire
  229. *  Attention, la teneur calculée utilise une porosité constante au cours
  230. *  du pas de temps.
  231. SAT   = 0.                                                            ;
  232. TEN   = 0.                                                            ;
  233. CAPA  = 0.                                                            ;
  234. NBDS  = 0                                                             ;
  235. PORO1 = 'MANU' 'CHPO' mailcent 1 'SCAL' 0. 'NATURE' 'DISCRET'         ;
  236. 'REPETER' BB ('DIME' LOIS.'INDEX')                                    ;
  237.    NOMT  = 'CHAINE' LOIS.'INDEX'.&BB                                  ;
  238.    LOI1  = LOIS . NOMT                                                ;
  239.    MOD1  = LOI1 . 'MODELE'                                            ;
  240.    NOMHT = 'MOT' ('TEXTE' ('CHAINE'  LOI1.'NOM_PROCEDURE')) ;
  241.  
  242. *  réduction de la pression au sous-domaine (local)
  243.    PNSL  = 'KCHT' MOD1 'SCAL' NIVEAU PNS                              ;
  244. *  calcul teneur et saturation
  245. *
  246.    SATL TENL CAPL= ('TEXTE' NOMHT) LOI1 PNSL                        ;
  247. *
  248. *  récupération porosité
  249.    PORL  = 'KCHT' LOI1.'MODELE' 'SCAL' NIVEAU LOI1.'PORO'             ;
  250.    PORL  = RECENTRE PORL MOD1 NIVEAU                                  ;
  251.  
  252. *- Calcul capacité capillaire = d teneur / d h
  253. *  méthode de la corde
  254.  
  255. **  avec la séparation en pression minimale :
  256. *   PNSL2 =  PNSL '-' LOI1.'PREC'                                      ;
  257. *   SATL2 TENL2 = ('TEXTE' nmproc) LOI1 PNSL2                          ;
  258. **  Dérivation
  259. *   CAP1  = TENL2 '-' TENL                                             ;
  260. *   CAP2  = 'ABS' (CAP1 '/' LOI1.'PREC')                               ;
  261. *  La capacité est nulle à saturation :
  262.    MP1   = 'MASQUE' PNSL 'INFERIEUR' 0.                                ;
  263.    MP    = 'KCHT' LOI1.'MODELE' 'SCAL' NIVEAU MP1                      ;
  264.    CAPL = CAPL '*' MP                                                  ;
  265.    CAPL = RECENTRE CAPL MOD1 NIVEAU                                    ;
  266. * ON EST ICI
  267.  
  268. *  nombre de faces/centres désaturés :
  269.    nbds  = nbds + ('ENTIER' ('MAXIMUM' ('RESULT' MP))) ;
  270.  
  271. *  concaténation
  272. * GBM a changé CENTRE en NIVEAU sur SAT !!!!!!!!
  273.    SAT   = 'KCHT' MDMH 'SCAL' NIVEAU SAT SATL                         ;
  274.    TEN   = 'KCHT' MDMH 'SCAL' NIVEAU TEN TENL                         ;
  275.    CAPA  = 'KCHT' MDMH 'SCAL' 'CENTRE' CAPA CAPL                      ;
  276.    PORO1 = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' PORO1 PORL       ;
  277.  
  278.    'DETRUIT' PNSL ;
  279.    'DETRUIT' SATL ;
  280.    'DETRUIT' TENL ;
  281.    'DETRUIT' PORL ;
  282. *   'DETRUIT' PNSL2 ;
  283. *   'DETRUIT' SATL2 ;
  284. *   'DETRUIT' TENL2 ;
  285. *   'DETRUIT' CAP1 ;
  286. *   'DETRUIT' CAP2 ;
  287.    'DETRUIT' MP1 ;
  288.    'DETRUIT' MP ;
  289.    'DETRUIT' CAPL ;
  290. 'FIN' BB                                                              ;
  291. 'RESPRO'  nbds SAT TEN CAPA PORO1                                     ;
  292. 'FINSI'                                                               ;
  293.  
  294. ***********************************************************************
  295.  
  296. 'SI' ('EGA' roro KALSOT)                                              ;
  297. 'ARGU' MDMH*'MMODEL' LOIS*'TABLE' NIVEAU*'MOT'
  298.          PNS*'CHPOINT'                                                ;
  299.  
  300.  
  301. * GBM  REGARDER SI PAS MEILLEURE BOUCLE QUAND TOUTES MEMES LOIS MAIS UN
  302. *      COEF PAR ZONE
  303. mailcent = 'DOMA' MDMH CENTRE                                         ;
  304. *- calcul teneur en eau, saturation et premier terme capacité capillaire
  305. *  Attention, la teneur calculée utilise une porosité constante au cours
  306. *  du pas de temps.
  307. *     GBM AU LIEU DE LA BOUCLE ON POURRAIT PLUTOT CONCATENER AU
  308. *     niveau des lois les coef et appliquer directement au champ
  309. *     c'est la proc NOM_HT qui ferait ce qu'il faut, tout le reste
  310. *     est faisable en creant des champs adéquats de coef.
  311. *--   Calcul teneur en eau, saturation et capacité capillaire
  312.       nbds = 0                                                        ;
  313.       TEN  = 0.                                                       ;
  314.       SAT  = 0.                                                       ;
  315.       CAPA = 0.                                                       ;
  316.       'REPETER' BB ('DIME' LOIS.'INDEX')                              ;
  317.          NOMT  = 'CHAINE' LOIS.'INDEX'.&BB                            ;
  318.          LOI1  = LOIS . NOMT                                          ;
  319.          MOD1  = LOI1 . 'MODELE'                                      ;
  320.          NOMHT = 'MOT' ('TEXTE' ('CHAINE' LOI1.'NOM_PROCEDURE'))      ;
  321.  
  322. *        pression et teneur en eau actuelles réduites au sous-domaine :
  323. *        PNS pression en eau
  324.          PNSL  = 'KCHT' MOD1 'SCAL' NIVEAU PNS                        ;
  325.          SATL TENL CAPL = ('TEXTE' NOMHT) LOI1 PNSL                   ;
  326. *        concaténation teneur et saturation
  327.          TEN   = 'KCHT' MDMH 'SCAL' NIVEAU TEN TENL                   ;
  328.          SAT   = 'KCHT' MDMH 'SCAL' NIVEAU SAT SATL                   ;
  329.  
  330. **        Détermination des zones saturées :
  331. **        GBM - attention test fouareux sur les inférieurs - pas précis
  332.          MP1   = 'MASQUE' PNSL 'INFERIEUR' -1.D-15                    ;
  333.          MP    = 'KCHT' MOD1 'SCAL' NIVEAU MP1                        ;
  334.  
  335. **        La capacité est nulle à saturation :
  336. **        GBM la remise à 0 est elle nécessaire ????????
  337. **            le calcul de Cap2 devrait aboutir naturellement à cela
  338. *            sinon C que le schéma n'est pas stable
  339.          CAPL = CAPL '*' MP                                           ;
  340.  
  341. *        GBM application du "decentrement"
  342.          CAPL = RECENTRE CAPL LOI1 . 'MODELE'  NIVEAU            ;
  343.  
  344. *        nombre de faces/centres désaturés :
  345. *        fouareux car peut etre sur les traces ou sommet
  346. *        il faut faire sur pression et non trace de pression
  347. *        GBM GBM GBM
  348.          nbds  = nbds + ('ENTIER' ('MAXIMUM' ('RESULT' MP)))          ;
  349.  
  350. *        concaténation capacité
  351.          CAPA = 'KCHT' MDMH 'SCAL' 'CENTRE' CAPA CAPL                 ;
  352.  
  353.  
  354.          'DETRUIT' PNSL                                               ;
  355.          'DETRUIT' SATL                                               ;
  356.          'DETRUIT' TENL                                               ;
  357.          'DETRUIT' CAPL                                               ;
  358.          'DETRUIT' MP                                                 ;
  359.          'DETRUIT' MP1                                                ;
  360.       'FIN' BB                                                        ;
  361. 'RESPRO' nbds SAT TEN CAPA                                            ;
  362. 'FINSI'                                                               ;
  363.  
  364. ************************************************************************
  365.  
  366. 'SI' ('EGA' roro LILOIPER)                                            ;
  367.  
  368. 'ARGUMENT'  SATUR*'TABLE'                                             ;
  369. 'SI' ('NON' ('EXISTE' SATUR 'LOI_PERMEABILITE'))                      ;
  370.   'ERREUR' 'La table LOI_PERMEABILITE est absente'                    ;
  371.   'QUITTER' DARCYSAT                                                  ;
  372. 'FINSI'                                                               ;
  373. LOIP    = 'TABLE' SATUR.'LOI_PERMEABILITE'                            ;
  374. 'SI' ('NON' ('EXISTE' LOIP 'SOUSTYPE'))                               ;
  375.    'ERREUR' 'La table LOI_PERMEABILITE n a pas de soustype'           ;
  376.     mess '(soustype : PUISSANCE, EXPONENTIELLE, LOGARITHMIQUE,'       ;
  377.     mess  ' MUALEM, BURDINE, MUALEM_BURDINE,BROOKS_COREY,'            ;
  378.     mess  ' MULTIZONE ou PERSONNELLE)'                                ;
  379.    'QUITTER' DARCYSAT                                                 ;
  380. 'FINSI'                                                               ;
  381.  
  382. * Gestion des sous-zones éventuelles
  383. 'SI' ('EGA' LOIP.'SOUSTYPE' 'MULTIZONE')                              ;
  384. *  Il y a plusieurs lois suivant les zones.
  385. *  On crée l'index qui permettra de parcourir les zones
  386. *  et on en retranche l'indice contenant le mot 'SOUSTYPE'
  387. *  De cette façon, il ne reste dans l'index que les sous-zones,
  388. *  et pas d'indice parasite.
  389.    'SI' ('NON' ('EXISTE' LOIP 'INDEX'))                               ;
  390.       LOIP.'INDEX' = 'INDEX' LOIP                                     ;
  391.    'FINSI'                                                            ;
  392.    dd = 'DIME' LOIP.'INDEX'                                           ;
  393.    'REPETER' bcl dd                                                   ;
  394.       i = &bcl                                                        ;
  395.       'SI' ('EGA' LOIP.'INDEX'.i 'SOUSTYPE')                          ;
  396. *        on écrase cet indice par les suivants :
  397.          'SI' (i < dd)                                                ;
  398.             'REPETER' bcl2 (dd - i)                                   ;
  399.                j = i + &bcl2                                          ;
  400.                LOIP.'INDEX'.(j-1) =  LOIP.'INDEX'.j                   ;
  401.             'FIN' bcl2                                                ;
  402.          'FINSI'                                                      ;
  403. *        en oubliant le dernier
  404.          'OUBLIER' (LOIP.'INDEX') dd                                  ;
  405.          'QUITTER' bcl                                                ;
  406.       'FINSI'                                                         ;
  407.    'FIN' bcl                                                          ;
  408.  
  409.    'REPETER' BB ('DIME' LOIP.'INDEX')                                 ;
  410. *     On parcourt toutes les sous-lois
  411.       NOMT = 'MOT' LOIP.'INDEX'.&BB                                   ;
  412.       LOI1 = LOIP . NOMT                                              ;
  413.  
  414. *     nom de la zone
  415.       LOI1.'NOMZONE' = 'CHAINE' NOMT                                  ;
  416. *     sous-type
  417.       'SI'  ('NON' ('EXISTE' LOI1 'SOUSTYPE'))                        ;
  418.          'ERREUR' ('CHAINE' 'La soustable ' NOMT
  419.                             ' de LOI_PERMEABILITE na pas de soustype');
  420.          mess '(soustype : PUISSANCE, EXPONENTIELLE, MUALEM, BURDINE,';
  421.          mess '       MUALEM_BURDINE,BROOKS_COREY,'                   ;
  422.          mess '       LOGARITHMIQUE, ou PERSONNELLE)'                 ;
  423.          'QUITTER' DARCYSAT                                           ;
  424.       'FINSI'                                                         ;
  425.       'SI' (('NEG' LOI1.'SOUSTYPE' 'PUISSANCE') 'ET'
  426.             ('NEG' LOI1.'SOUSTYPE' 'EXPONENTIELLE') 'ET'
  427.             ('NEG' LOI1.'SOUSTYPE' 'LOGARITMIQUE') 'ET'
  428.             ('NEG' LOI1.'SOUSTYPE' 'MUALEM') 'ET'
  429.             ('NEG' LOI1.'SOUSTYPE' 'BURDINE') 'ET'
  430.             ('NEG' LOI1.'SOUSTYPE' 'MUALEM_BURDINE') 'ET'
  431.             ('NEG' LOI1.'SOUSTYPE' 'BROOKS_COREY') 'ET'
  432.             ('NEG' LOI1.'SOUSTYPE' 'PERSONNELLE'))                    ;
  433.          'ERREUR' 'SOUSTYPE de loi inconnu'                           ;
  434.          'QUITTER' DARCYSAT                                           ;
  435.       'FINSI'                                                         ;
  436. *     nom procédure
  437.       'SI' ('EGA' LOI1.'SOUSTYPE' 'PERSONNELLE')                      ;
  438.          'SI' ('NON' ('EXISTE' LOI1 'NOM_PROCEDURE'))                 ;
  439.             'ERREUR' 'Il manque le nom de la procédure personnelle de'
  440.                      ' perméabilité pour la zone ' NOMT               ;
  441.          'FINSI'                                                      ;
  442.       'SINON'                                                         ;
  443. *        nom de la procédure standard
  444.          LOI1.'NOM_PROCEDURE' = 'MOT' 'KR_PRO'                        ;
  445.       'FINSI'                                                         ;
  446. *     modèle
  447.       'SI' ('NON' ('EXISTE' LOI1 'MODELE'))                           ;
  448.          'ERREUR' ('CHAINE' 'Il manque le MODELE '
  449.                             'dans LOI_PERMEABILITE.' NOMT)            ;
  450.          'QUITTER' DARCYSAT                                           ;
  451.       'FINSI'                                                         ;
  452.    'FIN' BB                                                           ;
  453. 'SINON'                                                               ;
  454. *  Il y a une seule loi pour tout le domaine
  455. *  on imite la structure d'une sous-loi pour généraliser
  456.  
  457. *  nom de la zone
  458.    NOMT  = 'MOT' 'DOMAINE_ENTIER'                                     ;
  459. *  on recopie coefficients et paramètres, et on crée un index
  460.    idtmp = 'INDEX' LOIP                                               ;
  461.    LOIP.NOMT = 'TABLE'                                                ;
  462.    'REPETER' bcl ('DIME' idtmp)                                       ;
  463.       id = 'MOT' idtmp.&bcl                                           ;
  464.       'SI' ('NEG' id NOMT)                                            ;
  465.          LOIP.NOMT.id = LOIP.id                                       ;
  466.       'FINSI'                                                         ;
  467.    'FIN' bcl                                                          ;
  468.    LOIP.'INDEX' = 'TABLE'                                             ;
  469.    LOIP.'INDEX'. 1 = 'MOT' NOMT                                       ;
  470.    NZONE = 1                                                          ;
  471.    LOIP.NOMT.'NOMZONE' = 'MOT' NOMT                                   ;
  472. *  sous-type
  473.    'SI' (('NEG' LOIp . 'SOUSTYPE' 'PUISSANCE') 'ET'
  474.             ('NEG' LOIp . 'SOUSTYPE' 'EXPONENTIELLE') 'ET'
  475.             ('NEG' LOIp . 'SOUSTYPE' 'LOGARITMIQUE') 'ET'
  476.             ('NEG' LOIp . 'SOUSTYPE' 'MUALEM') 'ET'
  477.             ('NEG' LOIp . 'SOUSTYPE' 'BURDINE') 'ET'
  478.             ('NEG' LOIp . 'SOUSTYPE' 'MUALEM_BURDINE') 'ET'
  479.             ('NEG' LOIp . 'SOUSTYPE' 'BROOKS_COREY') 'ET'
  480.             ('NEG' LOIp . 'SOUSTYPE' 'PERSONNELLE'))                  ;
  481.       'ERREUR' 'SOUSTYPE de loi inconnu'                              ;
  482.       'QUITTER' DARCYSAT                                              ;
  483.    'FINSI'                                                            ;
  484. *  nom procédure
  485.    'SI' ('EGA' LOIP.'SOUSTYPE' 'PERSONNELLE')                         ;
  486.       'SI' ('NON' ('EXISTE' LOIP 'NOM_PROCEDURE'))                    ;
  487.          'ERREUR' 'Il manque le nom de la procédure personnelle de'
  488.                   ' perméabilité pour la zone ' NOMT                  ;
  489.       'FINSI'                                                         ;
  490.    'SINON'                                                            ;
  491.       LOIP.NOMT.'NOM_PROCEDURE' =  'MOT' 'KR_PRO'                     ;
  492.    'FINSI'                                                            ;
  493. *  modèle
  494.    LOIP.NOMT.'MODELE'  = MDMH                                         ;
  495.  
  496. 'FINSI'                                                               ;
  497. 'RESPRO'  LOIP                                                        ;
  498.  
  499. 'FINSI'                                                               ;
  500.  
  501. ***********************************************************************
  502.  
  503. 'SI' ('EGA' roro 'LILOISAT')                                          ;
  504. 'ARGU' SATUR*'TABLE'                                                  ;
  505.  
  506. 'SI' ('NON' ('EXISTE' SATUR 'LOI_SATURATION'))                        ;
  507.   'ERREUR' 'La table LOI_SATURATION est absente'                      ;
  508.   'QUITTER' DARCYSAT                                                  ;
  509. 'FINSI'                                                               ;
  510.  
  511. LOIS    = 'TABLE' SATUR.'LOI_SATURATION'                              ;
  512.  
  513. 'SI' ('NON' ('EXISTE' LOIS 'SOUSTYPE'))                               ;
  514.    'ERREUR' 'La table LOI_SATURATION n a pas de soustype'             ;
  515.    'MESS' 'Soustype MOT = VAN_GENUCHTEN, EXPONENTIELLE, LOGARITHMIQUE';
  516.    'MESS' '               MULTIZONE ou PERSONNELLE'                   ;
  517.    'QUITTER' DARCYSAT                                                 ;
  518. 'FINSI'                                                               ;
  519.  
  520. LOISAT  = 'CHAINE' LOIS.'SOUSTYPE'                                    ;
  521.  
  522. * Gestion des sous-zones éventuelles
  523. 'SI' ('EGA' LOISAT 'MULTIZONE')                                       ;
  524. *  Il y a plusieurs lois suivant les zones.
  525. *  On crée l'index qui permettra de parcourir les zones
  526. *  et on en retranche l'indice contenant le mot 'SOUSTYPE'
  527. *  De cette façon, il ne reste dans l'index que les sous-zones,
  528. *  et pas d'indice parasite.
  529.    'SI' ('NON' ('EXISTE' LOIS 'INDEX'))                               ;
  530.       LOIS.'INDEX' = 'INDEX' LOIS                                     ;
  531.    'FINSI'                                                            ;
  532.    dd = 'DIME' LOIS.'INDEX'                                           ;
  533.    'REPETER' bcl dd                                                   ;
  534.       i = &bcl                                                        ;
  535.       'SI' ('EGA' LOIS.'INDEX'.i 'SOUSTYPE')                          ;
  536. *        on écrase cet indice par les suivants :
  537.          'SI' (i < dd)                                                ;
  538.             'REPETER' bcl2 (dd - i)                                   ;
  539.                j = i + &bcl2                                          ;
  540.                LOIS.'INDEX'.(j-1) =  LOIS.'INDEX'.j                   ;
  541.             'FIN' bcl2                                                ;
  542.          'FINSI'                                                      ;
  543. *        en oubliant le dernier
  544.          'OUBLIER' (LOIS.'INDEX') dd                                  ;
  545.          'QUITTER' bcl                                                ;
  546.       'FINSI'                                                         ;
  547.    'FIN' bcl                                                          ;
  548.  
  549.    'REPETER' BB ('DIME' LOIS.'INDEX') ;
  550. *     On parcourt toutes les sous-lois
  551. * CS      NOMT = 'MOT' LOIS.'INDEX'.&BB                                   ;
  552.       NOMT = LOIS.'INDEX'.&BB                                         ;
  553.       LOI1 = LOIS . NOMT                                              ;
  554.  
  555. *     nom de la zone
  556.       LOI1.'NOMZONE' = 'CHAINE' NOMT                                  ;
  557. *     sous-type
  558.       'SI'  ('NON' ('EXISTE' LOI1 'SOUSTYPE'))                        ;
  559.          'ERREUR' ('CHAINE' 'La soustable ' NOMT
  560.                             ' de LOI_SATURATION na pas de soustype')  ;
  561.          mess '(soustype : VAN_GENUCHTEN, EXPONENTIELLE,'             ;
  562.          mess '       LOGARITHMIQUE, ou PERSONNELLE)'                 ;
  563.          'QUITTER' DARCYSAT                                           ;
  564.       'FINSI'                                                         ;
  565.       LOISAT1 = 'CHAINE' LOI1.'SOUSTYPE'                              ;
  566.       'SI' (('NEG' LOISAT1 'VAN_GENUCHTEN') 'ET'
  567.             ('NEG' LOISAT1 'EXPONENTIELLE') 'ET'
  568.             ('NEG' LOISAT1 'LOGARITHMIQUE') 'ET'
  569.             ('NEG' LOISAT1 'PERSONNELLE'))                            ;
  570.          'MESS' 'soustype actuel 1' LOISAT1                           ;
  571.          'ERREUR' 'SOUSTYPE de loi inconnu'                           ;
  572.          'QUITTER' DARCYSAT                                           ;
  573.       'FINSI'                                                         ;
  574. *     nom procédure
  575.       'SI' ('EGA' LOISAT1 'PERSONNELLE')                              ;
  576.          'SI' ('NON' ('EXISTE' LOI1 'NOM_PROCEDURE'))                 ;
  577.             'ERREUR' 'Il manque le nom de la procédure personnelle de'
  578.                      ' saturation pour la zone ' NOMT                 ;
  579.  
  580.           'FINSI'                                                     ;
  581.       'SINON'                                                         ;
  582. *        nom de la procédure standard
  583.         LOI1.'NOM_PROCEDURE'  = 'MOT' 'HT_PRO'                        ;
  584.       'FINSI'                                                         ;
  585. *     modèle
  586.       'SI' ('NON' ('EXISTE' LOI1 'MODELE'))                           ;
  587.          'ERREUR' ('CHAINE' 'Il manque le MODELE'
  588.                             'dans LOI_SATURATION.' NOMT)              ;
  589.          'QUITTER' DARCYSAT                                           ;
  590.       'FINSI'                                                         ;
  591.    'FIN' BB                                                           ;
  592.  
  593. 'SINON'                                                               ;
  594. *  Il y a une seule loi pour tout le domaine
  595. *  on imite la structure d'une sous-loi pour généraliser
  596.  
  597. *  nom de la zone
  598.    NOMT  = 'MOT' 'DOMAINE_ENTIER'                                     ;
  599. *  on recopie coefficients et paramètres, et on crée un index
  600.    idtmp = 'INDEX' LOIS                                               ;
  601.    LOIS.NOMT = 'TABLE'                                                ;
  602.    'REPETER' bcl ('DIME' idtmp)                                       ;
  603.       id = 'MOT' idtmp.&bcl                                           ;
  604.       'SI' ('NEG' id NOMT)                                            ;
  605.          LOIS.NOMT.id = LOIS.id                                       ;
  606.       'FINSI'                                                         ;
  607.    'FIN' bcl                                                          ;
  608.    LOIS.'INDEX' = 'TABLE'                                             ;
  609.    LOIS.'INDEX'. 1 = 'MOT' NOMT                                       ;
  610.    NZONE = 1                                                          ;
  611.    LOIS.NOMT.'NOMZONE' = 'MOT' NOMT                                   ;
  612. *  sous-type
  613.    'SI' (('NEG' LOISAT 'VAN_GENUCHTEN') 'ET'
  614.          ('NEG' LOISAT 'EXPONENTIELLE') 'ET'
  615.          ('NEG' LOISAT 'LOGARITHMIQUE') 'ET'
  616.          ('NEG' LOISAT 'PERSONNELLE'))                                ;
  617.       'MESS' 'soustype actuel' LOISAT                                 ;
  618.       'ERREUR' 'SOUSTYPE de loi inconnu'                              ;
  619.       'QUITTER' DARCYSAT                                              ;
  620.    'FINSI'                                                            ;
  621. *  nom procédure
  622.    'SI' ('EGA' LOISAT 'PERSONNELLE')                                  ;
  623.       'SI' ('NON' ('EXISTE' LOIS 'NOM_PROCEDURE'))                    ;
  624.          'ERREUR' 'Il manque le nom de la procédure personnelle de'
  625.                   ' saturation pour la zone ' NOMT                    ;
  626. *      'SINON'         ;
  627.      LOIS.NOMT.'NOM_PROCEDURE'= LOIS.'NOM_PROCEDURE'                  ;
  628.       'FINSI'                                                         ;
  629.    'SINON'                                                            ;
  630. *     nom de la procédure standard
  631.      LOIS.NOMT.'NOM_PROCEDURE' = 'MOT' 'HT_PRO'                       ;
  632.    'FINSI'                                                            ;
  633. *  modèle
  634.    LOIS.NOMT.'MODELE'  = MDMH                                         ;
  635.  
  636. 'FINSI'                                                               ;
  637.  
  638. 'RESPRO' LOIS                                                         ;
  639.  
  640. 'FINSI'                                                               ;
  641.  
  642. ***********************************************************************
  643.  
  644. ***********************************************************************
  645. *--------------------------------------------------------------------*
  646. * RECUPERATION DES DONNEES NUMERIQUES                                *
  647. *--------------------------------------------------------------------*
  648.  
  649. 'SI' ('EGA' roro LIPARANU)                                            ;
  650.  
  651. 'ARGU'  SATUR*'TABLE'                                                 ;
  652. *
  653. *                                                      HOMOGENEISATION
  654. 'SI' ('NON' ('EXISTE' SATUR 'HOMOGENEISATION'))                       ;
  655.    NIVEAU = 'TEXTE' (CHAINE 'CENTRE')                                 ;
  656.    NIVEAU = CHAINE 'CENTRE'                                           ;
  657. 'SINON'                                                               ;
  658.    'SI' ('EGA' SATUR.'HOMOGENEISATION' 'CENTRE')                      ;
  659.       NIVEAU = 'TEXTE' ('CHAINE' 'CENTRE')                            ;
  660.       NIVEAU = 'CHAINE' 'CENTRE'                                      ;
  661.    'SINON'                                                            ;
  662.       'SI' ('EGA' SATUR.'HOMOGENEISATION' 'DECENTRE')                 ;
  663.          NIVEAU = 'TEXTE' ('CHAINE' 'MSOMMET')                        ;
  664.          NIVEAU = 'CHAINE' 'MSOMMET'                                  ;
  665.       'SINON'                                                         ;
  666.            'SI' ('EGA' SATUR.'HOMOGENEISATION' 'TRACE')               ;
  667.               NIVEAU =  'TEXTE' ('CHAINE' 'FACE')                     ;
  668. *              NIVEAU = 'CHAINE' 'FACE'                                ;
  669.            'SINON'                                                    ;
  670.          'ERREUR' 'la méthode d homogénisation n est pas reconnue'    ;
  671.          'QUITTER' DARCYSAT                                           ;
  672.             'FINSI'                                                   ;
  673.       'FINSI'                                                         ;
  674.    'FINSI'                                                            ;
  675. 'FINSI'                                                               ;
  676. *                                                                THETA
  677. 'SI' ( 'EXISTE' SATUR 'THETA' )                                       ;
  678.      TETA = SATUR . 'THETA'                                           ;
  679. 'SINON'                                                               ;
  680.      TETA = 1.D0                                                      ;
  681. 'FINSI'                                                               ;
  682.  
  683. *                                                                 NPAS
  684. 'SI' ( 'EXISTE' SATUR 'NPAS' )                                        ;
  685.      NPAS0 = SATUR . 'NPAS'                                           ;
  686. 'SINON'                                                               ;
  687. *     'SI' ( SATUR.'HYDRO_MECA' ) ;
  688. **        cas du couplage avec la méca, Darcysat sert de façon incrémentale
  689. *         NPAS0 = 1 ;
  690. *     'SINON' ;
  691. *        boucle infinie (arrêtée par TEMPS_FINAL)
  692.          NPAS0 = 0                                                    ;
  693. *     'FINSI' ;
  694. 'FINSI'                                                               ;
  695.  
  696. *                                                            SOUS RELAX
  697.  
  698. 'SI' ('NON' ('EXISTE' SATUR 'SOUS_RELAXATION'))                       ;
  699.    SATUR.'SOUS_RELAXATION' = 1.                                       ;
  700. 'FINSI'                                                               ;
  701.  
  702. *                                                                NITER
  703. 'SI' ( 'EXISTE' SATUR 'NITER' )                                       ;
  704.      NITER0 = SATUR . 'NITER'                                         ;
  705. 'SINON'                                                               ;
  706.      NITER0 = 10                                                      ;
  707. 'FINSI'                                                               ;
  708.  
  709.  
  710. *                                                           RESIDU_MAX
  711. 'SI' ( 'EXISTE' SATUR 'RESIDU_MAX' )                                  ;
  712.    ERR0 = SATUR . 'RESIDU_MAX'                                        ;
  713. 'SINON'                                                               ;
  714.    ERR0 = 1.D-4                                                       ;
  715. 'FINSI'                                                               ;
  716.  
  717. SATUR . 'CFL' = 1.D0                                                  ;
  718. CFL0 = 1.D0                                                           ;
  719.  
  720. *                                                               ITMAXI
  721. * nb max d'itérations non linéaires
  722. 'SI' ( 'EXISTE' SATUR 'ITMAX')                                        ;
  723.     itmaxi  = SATUR . 'ITMAX'                                         ;
  724. 'SINON'                                                               ;
  725.     itmaxi  = 40                                                      ;
  726. 'FINSI'                                                               ;
  727.  
  728. *                                                                   XI
  729. * coefficient de pénalisation à utiliser éventuellement
  730. * La pénalisation a priorité sur la réduction de pas de temps.
  731. OKPENAL = 'EXISTE' SATUR 'XI'                                         ;
  732. 'SI' OKPENAL                                                          ;
  733.    cofpenal = SATUR . 'XI'                                            ;
  734.    NPENALDT = 2                                                       ;
  735. *  seront inutilisés mais plus pratique pour les sorties
  736.    COFDIV = 1.D0                                                      ;
  737.    NMAXDT = 2                                                         ;
  738. 'SINON'                                                               ;
  739. *                                                               COFDIV
  740. *  astuce par défaut : on divise le pas de temps,
  741. *  d'un facteur COFDIV et un nombre NMAXDT maximum de fois.
  742.    'SI' ('EXISTE' SATUR 'COFDIV')                                     ;
  743.       COFDIV = SATUR.'COFDIV'                                         ;
  744.    'SINON'                                                            ;
  745.       COFDIV = 0.5                                                    ;
  746.    'FINSI'                                                            ;
  747. *                                                               NMAXDT
  748.    'SI' ('EXISTE' SATUR 'NMAXDT')                                     ;
  749.       NMAXDT = SATUR.'NMAXDT'                                         ;
  750.    'SINON'                                                            ;
  751.       NMAXDT = 6                                                      ;
  752.    'FINSI'                                                            ;
  753.    NPENALDT = NMAXDT                                                  ;
  754. * pas de penalisation
  755.    cofpenal = 0.D0                                                    ;
  756. 'FINSI'                                                               ;
  757.  
  758. 'SI' ('NON' ('EXISTE' SATUR 'TYPDISCRETISATION' ))                    ;
  759. *  On choisit VF par défaut
  760.    SATUR . 'TYPDISCRETISATION' = 'VF'                                 ;
  761. 'FINSI'                                                               ;
  762.  
  763. 'SI' ( ('EGA' SATUR . 'TYPDISCRETISATION'  'VF') 'ET'
  764.        ('EGA' SATUR . 'HOMOGENEISATION'  'TRACE') )                   ;
  765. 'ERREUR' 'méthode d homogénisation imcompatible avec les VF'       ;
  766. 'QUITTER' DARCYSAT                                                    ;
  767. 'FINSI'                                                               ;
  768.  
  769. 'RESPRO' NIVEAU TETA NPAS0 NITER0 ERR0 CFL0 ITMAXI OKPENAL
  770.           cofpenal npenaldt cofdiv nmaxdt                     ;
  771. *'RESPRO' NIVEAU;
  772.  
  773. 'FINSI'                                                               ;
  774.  
  775.  
  776.  
  777.  
  778. ***********************************************************************
  779.  
  780.  
  781.  
  782. 'SI' ('EGA' roro LITABTPS)                                            ;
  783.  
  784. 'ARGUMENT'   SATUR*'TABLE' TPSI*'FLOTTANT' debug*'LOGIQUE'
  785.           NPAS0*'ENTIER';
  786.  
  787. tmin   = TPSI ;
  788. DTAUTO = FAUX ;
  789. 'SI' ('EXISTE' SATUR 'TEMPS_CALCULES') ;
  790. *  Si les temps calculés sont spécifiés, on en déduit le temps final
  791. *  et le pas de temps initial
  792.    LTCALCUL = 'ORDONNER' (SATUR.'TEMPS_CALCULES') ;
  793.    DCAL     = 'DIME' LTCALCUL ;
  794.  
  795. *  détermination du 1er indice de temps à calculer > TPSINI
  796. *  c'est le prochain temps à calculer TPS = LTCALCUL(ICAL)
  797.    'REPETER' bouc1 DCAL ;
  798.       ICAL = &bouc1 ;
  799.       tps   = 'EXTR' LTCALCUL ICAL ;
  800.       'SI' (tps > TPSINI)  ;
  801.          'QUITTER' bouc1 ;
  802.       'FINSI' ;
  803.    'FIN' bouc1 ;
  804.  
  805.    tmin   = 'MINI' LTCALCUL ;
  806.    NPAS0  = DCAL - ICAL + 1 ;
  807.    DELTAT = TPS - TPSINI ;
  808.  
  809.    'SI' ('EXISTE' SATUR 'TEMPS_FINAL');
  810.       'SI' (SATUR.'TEMPS_FINAL' < ('MAXIMUM' LTCALCUL)) ;
  811.          TFINAL = SATUR.'TEMPS_FINAL' ;
  812.          'SI' debug ;
  813.            mess 'Le calcul sera conduit jusqu a TEMPS_FINAL, malgré que'
  814.                 ' la liste TEMPS_CALCULES aille au-delà' ;
  815.          'FINSI' ;
  816.       'SINON' ;
  817.          TFINAL = 'MAXIMUM' LTCALCUL ;
  818.          'SI' debug ;
  819.             mess 'indice TEMPS_FINAL ignoré car supérieur ou égal au '
  820.                  'dernier temps spécifié par TEMPS_CALCULES';
  821.          'FINSI' ;
  822.       'FINSI' ;
  823.    'SINON' ;
  824.       TFINAL = 'MAXIMUM' LTCALCUL ;
  825.    'FINSI' ;
  826.    'SI' ('EXISTE' SATUR 'DT_INITIAL');
  827.       'SI' debug ;
  828.         mess 'indice DT_INITIAL ignoré, car on respecte TEMPS_CALCULES';
  829.       'FINSI' ;
  830.    'FINSI' ;
  831.  
  832. 'SINON' ;
  833. *   Sinon, on a besoin du temps final pour savoir où s'arrêter
  834.    'SI' ('EXISTE' SATUR 'TEMPS_FINAL');
  835.       TFINAL = SATUR . 'TEMPS_FINAL' ;
  836.    'SINON'                          ;
  837.       'ERREUR' 'Renseigner TEMPS_FINAL si TEMPS_CALCULES laissé libre ';
  838.       'QUITTER' DARCYSAT             ;
  839.    'FINSI'                          ;
  840.  
  841.    'SI' ('EXISTE' SATUR 'DT_INITIAL');
  842.       DELTAT = SATUR .'DT_INITIAL';
  843.    'SINON' ;
  844. *     pas de temps calculé automatiquement plus tard
  845. *     DELTAT inutil mais initialisé pour sortie plus simple
  846.       DELTAT = 0.D0;
  847.       DTAUTO = VRAI ;
  848.    'FINSI' ;
  849.  
  850. * parametre inutils mais pour sortie unifiee. Leur valeur est
  851. * toutefois non pipot au cas où on s'en servirait.
  852.    DCAL = 0.D0;
  853.    LTCALCUL = 'PROG' ;
  854.    ICAL = 0;
  855. 'FINSI' ;
  856.  
  857. 'SI' ( 'EXISTE' SATUR 'TEMPS_SAUVES' ) ;
  858.    TPSOR = 'ORDONNER' (SATUR . 'TEMPS_SAUVES') ;
  859.    DSOR  = 'DIME' TPSOR ;
  860.    ISOR  = 0 ;
  861.    IOK1  = FAUX ;
  862.    'REPETER' BOU3 DSOR ;
  863.       ISOR = ISOR + 1 ;
  864.       TEMS = 'EXTR' TPSOR ISOR ;
  865.       'SI' ( ( tmin '<' (TEMS - 1.e-6) ) 'ET'
  866.              ( TEMS   '&lt;EG' (TFINAL + 1.e-6)) ) ;
  867.          IOK1 = VRAI ;
  868.          'QUITTER' BOU3 ;
  869.       'FINSI' ;
  870.    'FIN' BOU3 ;
  871.    'SI' ('NON' IOK1) ;
  872.       'ERREUR'
  873.         'Aucun temps de sauvegarde imposé dans l intervalle à calculer';
  874.       'QUITTER' DARCYTRA ;
  875.    'FINSI' ;
  876. 'SINON' ;
  877.    ISOR = 0;
  878.    DSOR = 0;
  879. *  on met le temps final dans liste tps sortie - GBM tester
  880.    TPSOR = 'PROG' TFINAL ;
  881. 'FINSI' ;
  882. 'RESPRO'  DCAL LTCALCUL DTAUTO ICAL tmin NPAS0 DELTAT TFINAL isor tpsor
  883. dsor ;
  884.  
  885. 'FINSI'                                                               ;
  886.  
  887.  
  888.  
  889.  
  890. ***********************************************************************
  891.  
  892.  
  893.  
  894. 'SI' ('EGA' roro REKAPITU)                                            ;
  895.  
  896. 'ARGUMENT'  SATUR*'TABLE' L_GRAV*'LOGIQUE'
  897.          rhowg*'FLOTTANT' DAXE*'POINT' convh*'FLOTTANT'
  898.          cfl0*'FLOTTANT' teta*'FLOTTANT' LAST1*'ENTIER'
  899.           debug*'LOGIQUE'                                             ;
  900.  
  901.  
  902.  
  903. 'SI' debug                                                            ;
  904.    'SAUTER' 1 'LIGNE'                                                 ;
  905.    'MESS' 'Modélisation Darcy Eléments Finis Mixtes Hybrides'
  906.           ' en transitoire en milieu saturé et non-saturé.'           ;
  907.    'MESS' '-------------------------------------------------'
  908.           '-----------------------------------------------'           ;
  909.    'SAUTER' 1 'LIGNE'                                                 ;
  910.    'MESS' '    --------------------------------------------- '        ;
  911.    'MESS' '   |              d h          --->              |'        ;
  912.    'MESS' '   |  ( Ca + Ce ) ---  - div(K grad H) = Q       |'        ;
  913.    'MESS' '   |              d t                            |'        ;
  914.    'MESS' '   |                                             |'        ;
  915.    'MESS' '    --------------------------------------------- '        ;
  916.    'MESS'                                                             ;
  917.    'MESS' 'avec :'                                                    ;
  918.    'MESS' '       d, dérivée ronde '                                  ;
  919.    'MESS' '       Ca, capacité capillaire [1/m]'                      ;
  920.    'MESS' '       Ce, coefficient d emmagasinement [1/m]'             ;
  921.    'MESS' '       h, succion capillaire (Pgaz-Peau)/(rho*g) [m]'      ;
  922.    'MESS' '       H, charge h+z [m]'                                  ;
  923.    'MESS' '       Q, terme source [s-1]'                              ;
  924.    'MESS' ' '                                                         ;
  925.    'MESS' 'Données présentes en entrée : '                            ;
  926.  
  927.    'SI' ( 'EXISTE' SATUR 'SOURCE' )                                   ;
  928.         'MESS' '  Ce problème comporte un terme source'               ;
  929.    'FINSI'                                                            ;
  930.  
  931.    'SI' L_GRAV                                                        ;
  932.       'MESS' '  Ce problème comporte un effet gravitationnel '
  933.              'de valeur ' (@ARR RHOWG 2) ' N'                         ;
  934.          XAXE = 'COORDONNEE' 1 DAXE                                   ;
  935.          YAXE = 'COORDONNEE' 2 DAXE                                   ;
  936.          'SI' (NDIME 'EGA' 2)                                         ;
  937.             'MESS' '    Direction de l axe de pesanteur :'            ;
  938.             'MESS' ('CHAIN' '                    -> '
  939.                              '                  ->')                  ;
  940.             'MESS' '    ' XAXE 'X + ' YAXE ' Y'                       ;
  941.          'SINON'                                                      ;
  942.             ZAXE = 'COORDONNEE' 3 DAXE                                ;
  943.             'MESS' '    Direction de l axe de pesanteur :'            ;
  944.             'MESS' ( 'CHAIN' '                    -> '
  945.                         '                  ->                   ->' ) ;
  946.             'MESS' '    ' XAXE 'X + ' YAXE ' Y + ' ZAXE ' Z'          ;
  947.          'FINSI'                                                      ;
  948.    'SINON'                                                            ;
  949.       'MESS' '  Ce problème est SANS effet gravitationnel'            ;
  950.    'FINSI'                                                            ;
  951.  
  952.    'SI' ('EGA' CONVH 1.)                                              ;
  953.       'MESS' '  Les pressions sont exprimées en Pascal'               ;
  954.    'SINON'                                                            ;
  955.       'MESS' '  Les pressions sont exprimées en Pascal/'
  956.              (@ARR CONVH 0)                                           ;
  957.    'FINSI'                                                            ;
  958.  
  959.    'SI' ( 'EXISTE' SATUR 'COEF_EMMAGASINEMENT' )                      ;
  960.       'MESS' '  Ce problème comporte un coefficient d emmagasinement' ;
  961.    'SINON'                                                            ;
  962.       'MESS' '  Ce problème est SANS coefficient d emmagasinement'    ;
  963.    'FINSI'                                                            ;
  964.  
  965.  
  966.  
  967.    'SI' ( 'EXISTE' SATUR 'XI' )                                       ;
  968.       'MESS' '  Paramètre de pénalisation : ' SATUR.'XI'              ;
  969.    'FINSI'                                                            ;
  970.    'MESS' '  Relaxation : ' SATUR.'SOUS_RELAXATION'                   ;
  971.    'MESS' '  CFL initial : ' cfl0                                     ;
  972.    'MESS' '  Méthode d homogénisation de K et C : '                   ;
  973.    'MESS'   SATUR.'HOMOGENEISATION'                                   ;
  974.    'MESS' ' '                                                         ;
  975.    'MESS' '  Valeur de THETA, parametre du schéma diffusif : ' TETA   ;
  976.    'MESS' '    (0:Schéma explicite, 0.5:Crank-Nicholson, 1:Implicite)';
  977. 'FINSI'                                                               ;
  978.  
  979.  
  980. 'MESSAGE' 'hihihihiihih ' roro;
  981. 'SI' (LAST1 > 0)                                                      ;
  982.    'MESSAGE' roro                                                     ;
  983.    'MESS' 'Reprise de calcul au pas : ' LAST1                         ;
  984. 'FINSI'                                                               ;
  985. 'MESS' ' '                                                            ;
  986.  
  987. 'FINSI'                                                               ;
  988.  
  989.  
  990.  
  991.  
  992.  
  993.  
  994. ***********************************************************************
  995.  
  996.  
  997.  
  998.  
  999.  
  1000. 'SI' ('EGA' roro SAVRESU)                                             ;
  1001.  
  1002. 'ARGUMENT'  SATUR*'TABLE' TPSOR/'LISTREEL'
  1003.                   ISOR*'ENTIER'
  1004.                   TPSANC*'FLOTTANT' TPS*'FLOTTANT' ISAUV*'ENTIER'
  1005.                   DELTAT*'FLOTTANT' TRHANC*'CHPOINT' TRH*'CHPOINT'
  1006.                   CHRGANC*'CHPOINT' CHRG*'CHPOINT' TENANC*'CHPOINT'
  1007.                   TEN*'CHPOINT' PNSANC*'CHPOINT' PNS*'CHPOINT'
  1008.                   SATANC*'CHPOINT' SAT*'CHPOINT' FLUANC*'CHPOINT'
  1009.                   FLU0*'CHPOINT' NIVEAU*'MOT' DSOR*'ENTIER'
  1010.                   nnind*'ENTIER' MDMH*'MMODEL'                       ;
  1011.  
  1012.  
  1013.   'SI' ( 'EXISTE' SATUR 'TEMPS_SAUVES' )                             ;
  1014. *    Sauvegarde de tous les temps intermédiaires, s'il y en a :
  1015.      'REPETER' BOU9                                                  ;
  1016.         TEMS = 'EXTR' TPSOR ISOR                                     ;
  1017.  
  1018.         'SI' (isauv < 0)                                             ;
  1019. * sauvegarde initiale
  1020.             isauv = isauv + 1                                        ;
  1021.  
  1022.             DTEM  = ( TPS - TEMS ) / DELTAT                          ;
  1023.             UNMO  = 1. - DTEM                                        ;
  1024. *          'SI' ('NON' (SATUR.'HYDRO_MECA')) ;
  1025. *            Lorsque DARCYSAT est appellé par PASAPAS, c'est PASAPAS qui
  1026. *            sauvegarde le temps calculé, sinon, c'est DARCYSAT.
  1027.             SATUR.'TEMPS'    .isauv = TEMS                           ;
  1028. *          'FINSI' ;
  1029.             SATUR.'TRACE_CHARGE'.isauv = 'COLI' TRHANC  DTEM TRH  UNMO ;
  1030.             SATUR.'CHARGE'      .isauv = 'COLI' CHRGANC DTEM CHRG UNMO ;
  1031.             SATUR.'FLUX'     .isauv = 'COLI' FLUANC  DTEM FLU0 UNMO    ;
  1032.             SATUR. 'TH2O'       .isauv = RECENTRE
  1033.              ('COLI' TENANC  DTEM TEN  UNMO) MDMH NIVEAU               ;
  1034. *            SATUR. 'PRESSION'   .isauv = RECENTRE
  1035. *                          ('COLI' PNSANC  DTEM PNS  UNMO) MDMH NIVEAU ;
  1036.             SATUR. 'PRESSION'   .isauv =
  1037.                           'COLI' PNSANC  DTEM PNS  UNMO                ;
  1038.  
  1039.             SATUR. 'SATURATION' .isauv = RECENTRE
  1040.             ('COLI' SATANC  DTEM SAT  UNMO) MDMH         NIVEAU        ;
  1041.             ISOR = ISOR + 1                                            ;
  1042.  
  1043.         'FINSI'                                                        ;
  1044.         'SI' ((TPSANC '<' TEMS) 'ET' (TEMS '&lt;EG' TPS))                 ;
  1045.  
  1046.             isauv = isauv + 1                                          ;
  1047.  
  1048.             DTEM  = ( TPS - TEMS ) / DELTAT                            ;
  1049.             UNMO  = 1. - DTEM                                          ;
  1050. *          'SI' ('NON' (SATUR.'HYDRO_MECA')) ;
  1051. *            Lorsque DARCYSAT est appellé par PASAPAS, c'est PASAPAS qui
  1052. *            sauvegarde le temps calculé, sinon, c'est DARCYSAT.
  1053.             SATUR.'TEMPS'    .isauv = TEMS                             ;
  1054. *          'FINSI' ;
  1055.             SATUR.'TRACE_CHARGE'.isauv = 'COLI' TRHANC  DTEM TRH  UNMO ;
  1056.             SATUR.'CHARGE'      .isauv = 'COLI' CHRGANC DTEM CHRG UNMO ;
  1057.             SATUR.'FLUX'     .isauv = 'COLI' FLUANC  DTEM FLU0 UNMO    ;
  1058.             SATUR. 'TH2O'       .isauv = RECENTRE
  1059.             ('COLI' TENANC  DTEM TEN  UNMO) MDMH  NIVEAU               ;
  1060. *            SATUR. 'PRESSION'   .isauv = RECENTRE
  1061. *                          ('COLI' PNSANC  DTEM PNS  UNMO) MDMH NIVEAU ;
  1062.             SATUR. 'PRESSION'   .isauv =
  1063.                           'COLI' PNSANC  DTEM PNS  UNMO                ;
  1064.  
  1065.             SATUR. 'SATURATION' .isauv = RECENTRE
  1066.             ('COLI' SATANC  DTEM SAT  UNMO) MDMH         NIVEAU        ;
  1067.             ISOR = ISOR + 1                                            ;
  1068.        'SINON'                                                         ;
  1069.             'QUITTER' BOU9                                             ;
  1070.        'FINSI'                                                         ;
  1071.        'SI' ( ISOR '>' DSOR )                                          ;
  1072.             ISOR = DSOR                                                ;
  1073.             'QUITTER' BOU9                                             ;
  1074.        'FINSI'                                                         ;
  1075.     'FIN' BOU9                                                         ;
  1076. *   on ne sauvegarde le nb d'itérations de Picard que si les temps
  1077. *   de sauvegarde correspondent à des temps calculés.
  1078.     SATUR.'ITERPICARD' = SATUR.'ITERPICARD' 'ET' ('LECT' &linear)      ;
  1079.   'SINON'                                                              ;
  1080.  
  1081.  
  1082.  
  1083.      isauv = isauv + 1                                                 ;
  1084.      SATUR. 'TEMPS'     . isauv = TPS                                  ;
  1085.      SATUR. 'TRACE_CHARGE' . isauv = TRH                               ;
  1086.      SATUR. 'CHARGE'       . isauv = CHRG                              ;
  1087.         SATUR. 'FLUX'      . isauv = FLU0                              ;
  1088.      SATUR. 'TH2O'         . isauv = RECENTRE TEN MDMH  NIVEAU         ;
  1089. * correction pas de recentrage de la pression dejà au centre
  1090.      SATUR. 'PRESSION'     . isauv =  PNS          ;
  1091.      SATUR. 'SATURATION'   . isauv = RECENTRE SAT MDMH  NIVEAU         ;
  1092.   'FINSI' ;
  1093.  
  1094. * on met de côté les temps effectivement calculés
  1095.   SATUR.'TEMPS_EFFECTUES' = SATUR.'TEMPS_EFFECTUES' 'ET' ('PROG' TPS)  ;
  1096.  
  1097. 'RESPRO'  ISOR ISAUV                                                   ;
  1098.  
  1099. 'FINSI'                                                                ;
  1100.  
  1101.  
  1102.  
  1103. ***********************************************************************
  1104.  
  1105.  
  1106.  
  1107.  
  1108. 'SI' ('EGA' roro TESARRET)                                            ;
  1109.  
  1110. 'ARGUMENT'  ipen*'ENTIER' npen*'ENTIER' ila*'ENTIER'
  1111.                    ilin*'ENTIER' delt*'FLOTTANT' maxx*'FLOTTANT'
  1112.                    okpen*'LOGIQUE'  cofpe*'FLOTTANT'
  1113.                    CFLL*'FLOTTANT' COFD*'FLOTTANT' debg*'LOGIQUE'
  1114.                    SATUR*'TABLE' tpss*'FLOTTANT'                      ;
  1115.  
  1116. * ipen indice de boucle penalisation
  1117. * npen nbmax de boucle penalisation
  1118. * ila+1 ILAST + 1 indice boucle tps
  1119. * ilin indice de boucle LINEAR point fixe
  1120. * delt pas de tps
  1121. * maxx maximum residu de equation a resoudre
  1122. * okpen vrai si penalisation
  1123. * penn vrai si on effectue concretement une penalisation
  1124. *      au cours du pas traité
  1125. * cofpe coef de penalisation devant d/dt
  1126. * cfll cfl de Benet
  1127. * COFD coef de diminution des pas de temps
  1128. * debg vrai si messages affichés
  1129. * SATUR table du calcul
  1130. * tpss temps en cours
  1131.  
  1132. 'SI' (ipen < npen)                                                    ;
  1133.    'SI' debg                                                          ;
  1134.       'MESS' ' non convergence au temps ' tpss 's'                    ;
  1135.       'MESS' ' pas ' (ila+1) ', it. ' ilin
  1136.              ', résidu :' maxx ' delta t : ' delt                     ;
  1137.    'FINSI'                                                            ;
  1138.  
  1139.    'SI' OKPEN                                                         ;
  1140. *     pénalisation
  1141.       PENN = VRAI                                                     ;
  1142.       'SI' debg                                                       ;
  1143.          'MESS' ' on relance le calcul avec le coef. de '
  1144.                 ' pénalisation XI :' cofpe                            ;
  1145.       'FINSI'                                                         ;
  1146.    'SINON'                                                            ;
  1147.  
  1148. *     diminution du pas de temps
  1149.       CFLL      = CFLL * COFD                                         ;
  1150.       delt    = DELT * COFD                                           ;
  1151.       'SI' debg                                                       ;
  1152.          'MESS' ' on relance le calcul en contractant le pas de '
  1153.                  'temps d un facteur ' (@ARR COFD 1) ' ('
  1154.                   (@ARR ipen 0) 'eme essai)'                          ;
  1155.       'FINSI'                                                         ;
  1156.    'FINSI'                                                            ;
  1157.  
  1158.  
  1159. 'SINON'                                                               ;
  1160. *  si vraiment on ne veut pas converger, on le dit et on continue
  1161.    'SI' debg                                                          ;
  1162.       txt = 'CHAINE' 'Non convergence, meme avec '                    ;
  1163.       'SI' OKPEN                                                      ;
  1164.          txt = 'CHAINE' txt 'pénalisation'                            ;
  1165.       'SINON'                                                         ;
  1166.          txt = 'CHAINE' txt 'division du pas de temps par '
  1167.                              (COFD ** (-1. * (npen - 1)))             ;
  1168.        'FINSI'                                                        ;
  1169.        mess txt                                                       ;
  1170.        mess 'On passe au temps suivant, instationnaire faux.'         ;
  1171.    'FINSI'                                                            ;
  1172. 'FINSI'                                                               ;
  1173.  
  1174. *  on sauvegarde les temps où on a usé d'artifice (pénal ou division dt)
  1175. SATUR.'TEMPS_XI' = SATUR.'TEMPS_XI' 'ET' ('PROG' TPSS)                ;
  1176.  
  1177. 'RESPRO'  delt CFLL                                                   ;
  1178.  
  1179. 'FINSI'                                                               ;
  1180.  
  1181.  
  1182.  
  1183.  
  1184.  
  1185. ***********************************************************************
  1186.  
  1187.  
  1188.  
  1189.  
  1190.  
  1191. 'SI' ('EGA' roro TRONKP)                                            ;
  1192.  
  1193. 'ARGUMENT'  MDMH*'MMODEL' CHRG*'CHPOINT' TRH*'CHPOINT'
  1194.           L_GRAV*'LOGIQUE' NIVEAU*'MOT' RHOWG*'FLOTTANT'
  1195.           CONVH*'FLOTTANT'                                           ;
  1196.  
  1197. P1  = 'NOMC' CHRG 'SCAL'                                             ;
  1198. TR1  = 'NOMC' TRH 'SCAL'                                             ;
  1199.  
  1200. 'SI' ('EGA' NIVEAU FACE)                                             ;
  1201. TP1 = TR1                                                            ;
  1202. 'SINON'                                                              ;
  1203. 'SI' ('EGA' NIVEAU CENTRE)                                           ;
  1204.    TP1 = P1                                                          ;
  1205. SINON                                                                ;
  1206.    'SI' ('EXISTE' SATUR 'TRACE_IMPOSE')                              ;
  1207.    CHARIMPO = 'TIRE' SATUR . 'TRACE_IMPOSE' TPS                      ;
  1208.    CHARIMPO = 'NOMC' 'H' CHARIMPO                                    ;
  1209.    rchpo1 rchpo2 rchelem1 = PENT MDMH  'CENTRE' 'EULESCAL' 'NOLIMITE'
  1210.                                ('MOTS' 'H') CHRG
  1211.                                'CLIM' charimpo;
  1212. *   TP1 = 'ELNO' MDMH CHRG 'VOLUMF' rchpo1 rchpo2                    ;
  1213.    TP1 = 'ELNO' 'VF' (doma MDMH TABLE)
  1214.            (NOMC SCAL CHRG) rchpo1 rchpo2                            ;
  1215.    'SINON'                                                           ;
  1216.    rchpo1 rchpo2 rchelem1 = PENT MDMH 'CENTRE' 'EULESCAL' 'NOLIMITE'
  1217.                                ('MOTS' 'H') CHRG                     ;
  1218. *   TP1 = 'ELNO' MDMH CHRG 'VOLUMF' rchpo1 rchpo2                     ;
  1219.    TP1 = 'ELNO' 'VF' (doma MDMH TABLE) ( NOMC SCAL CHRG)
  1220.                                              rchpo1 rchpo2           ;
  1221.    'FINSI'                                                           ;
  1222. 'FINSI'                                                              ;
  1223. 'FINSI'                                                              ;
  1224.  
  1225.  
  1226. 'SI' L_GRAV                                                          ;
  1227.    P1  = P1  '-' ( RHOWG / CONVH * ZCC )                             ;
  1228.    'SI' (('EGA' NIVEAU 'MSOMMET') 'OU' ('EGA' NIVEAU 'FACE'))        ;
  1229.       TP1 =  TP1 '-' ( RHOWG / CONVH * ZFF )                         ;
  1230.     'FINSI'                                                          ;
  1231. 'FINSI'                                                              ;
  1232. P1  = 'KCHT' MDMH 'SCAL' 'CENTRE' P1                                 ;
  1233. 'SI' ('EGA' NIVEAU 'MSOMMET')                                        ;
  1234.    TP1 = 'KCHT' MDMH 'SCAL' 'SOMMET'  TP1                            ;
  1235.    PNS = 'KOPS' TP1 '>|' -1.D-15                                     ;
  1236. 'SINON'                                                              ;
  1237. 'SI' ('EGA' NIVEAU 'FACE')                                           ;
  1238.    TP1 = 'KCHT' MDMH 'SCAL' 'FACE'  TP1                              ;
  1239.    PNS = 'KOPS' TP1 '>|' -1.D-15                                     ;
  1240.  'SINON'                                                             ;
  1241.    PNS = 'KOPS' P1 '>|' -1.D-15                                      ;
  1242.  'FINSI'                                                             ;
  1243. 'FINSI'                                                              ;
  1244. * GBM P1 charge avec unité differente - PNS idem tronquée
  1245. * P1 ne devrait en théorie plus servir après nettoyage de la proc.
  1246. 'RESPRO'  PNS                                                         ;
  1247.  
  1248. 'FINSI'                                                               ;
  1249.  
  1250.  
  1251.  
  1252.  
  1253.  
  1254. ***********************************************************************
  1255.  
  1256.  
  1257.  
  1258.  
  1259.  
  1260. 'SI' ('EGA' roro YNTGFLUX)                                            ;
  1261.  
  1262. 'ARGUMENT' SATUR*'TABLE'                                              ;
  1263.  
  1264. startflu = 1.D30                                                      ;
  1265. startsou = 1.D30                                                      ;
  1266. startmix = 1.D30                                                      ;
  1267. *                                     integrale des cond de flux
  1268. *                                     CL de type Neumann : FLUX_IMPOSE
  1269. 'SI' ( 'EXISTE' SATUR 'FLUX_IMPOSE' )                                 ;
  1270.      FLUIMP = SATUR . 'FLUX_IMPOSE'                                   ;
  1271.      'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  1272.            'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  1273.                  'ET' ('EGA' teta 1.D0)))                             ;
  1274. *       On écrit les flux sous forme intégrale sauf en explicite et
  1275. *       kranck-Nickholson pour EFMH
  1276.         FLUIMP startflu = CHAMINT FLUIMP                              ;
  1277.      'FINSI'                                                          ;
  1278. 'SINON'                                                               ;
  1279.    FLUIMP = 0.D0                                                      ;
  1280. 'FINSI'                                                               ;
  1281. *                                            CL de type mixte : MIXTES
  1282. *                                 GBM IL Y AURA MODIF CAR HYDRAU PURE
  1283. 'SI' ( 'EXISTE' SATUR 'FLUMIXTE' )                                    ;
  1284.      FLUMIX = SATUR . 'FLUMIXTE' . 'VAL'                              ;
  1285.      'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  1286.            'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  1287.                  'ET' ('EGA' teta 1.D0)))                             ;
  1288. *       On écrit les flux sous forme intégrale sauf en explicite et
  1289. *       kranck-Nickholson pour EFMH
  1290.         FLUMIX startmix = CHAMINT FLUMIX                              ;
  1291.      'FINSI'                                                          ;
  1292. 'SINON'                                                               ;
  1293.    FLUMIX = 0.D0                                                      ;
  1294. 'FINSI'                                                               ;
  1295.  
  1296. 'SI' ( 'EXISTE' SATUR 'SOURCE' )                                      ;
  1297.      TERSOU = SATUR . 'SOURCE'                                        ;
  1298. * on crée une evolution intégrée temporellement pour assurer la
  1299. * conservation lors de liste de pas de temps de chargement differente
  1300. * des temps de discrétisation
  1301.      TERSOU startsou = CHAMINT TERSOU                                 ;
  1302. 'FINSI'                                                               ;
  1303.  
  1304. 'RESPRO'  startflu startmix startsou FLUIMP FLUMIX TERSOU             ;
  1305.  
  1306. 'FINSI'                                                               ;
  1307.  
  1308.  
  1309.  
  1310.  
  1311. ***********************************************************************
  1312.  
  1313.  
  1314.  
  1315. 'SI' ('EGA' roro YNYTDT)                                              ;
  1316.  
  1317.  
  1318. 'ARGUMENT'  MMDD*'MMODEL' FLU0*'CHPOINT' ABBM*'MCHAML'
  1319.           VVOL*'CHPOINT' CFFL*'FLOTTANT' deb*'LOGIQUE'
  1320.           daut*'LOGIQUE' DELTAT*'FLOTTANT' TERSCS/'CHPOINT'           ;
  1321.  
  1322.  
  1323. 'SI' (deb 'OU' daut)                                                  ;
  1324.  
  1325.    SI ('NEG' ('TYPE' TERSCS) 'CHPOINT')                               ;
  1326.       TERSCS = 0.D0 * FLU0                                            ;
  1327.    'FINSI'                                                            ;
  1328. * GBM MODIFIE MDMH EN MMDD
  1329.    IDT    = ('DIVU' MMDD ('NOMC' 'FLUX' FLU0) ABBM)                   ;
  1330.    IDT    = IDT '+' ('NOMC' 'SCAL' TERSCS)                            ;
  1331. *  temps qu'il faut pour remplir une maille avec ce débit moyen :
  1332.    DTI    = ((('MAXIMUM' ('ABS' ('RESULT' IDT)))
  1333.              '/' ('MAXIMUM' ('RESULT' VVOL)))  + 1.D-15) ** -1.       ;
  1334. *  on s'aligne sur la cinétique la plus rapide :
  1335.    DTI    = ('ABS'  DTI)                                              ;
  1336. * TEmps pour remplir 1/300 eme du domaine.
  1337. * 300 peut etre le nombre de pas viser.
  1338. * inclure la quantité d'eau présente dans le bilan de remplissage
  1339.    ndim = 'VALEUR' 'DIME'                                             ;
  1340. *  On considere minimum 50 mailles dans une direction plus
  1341. *  le nombre moyen de mailles par direction
  1342.    nnpas = 50.D0 '+' (('NBEL' ('DOMA' mmdd centre)) '**' (1.D0 / ndim));
  1343.    DTI = DTI '/' nnpas                                                ;
  1344. 'SINON'                                                               ;
  1345. * GBM mettre deltat en arg. C juste pour affichage ce truc.
  1346.    DTI = DELTAT / CFFL                                                ;
  1347. 'FINSI'                                                               ;
  1348. *  pas de temps automatique :
  1349. 'SI' (daut)                                                           ;
  1350. *  GBM
  1351. *   DELTAT = 'MINIMUM' ('PROG' DTI (DELTAT * 1.2D0))                   ;
  1352. *   DELTAT = DELTAT * CFFL                                             ;
  1353.     DELTAT = 'MINIMUM' ('PROG' DTI (CFFL * DELTAT * 1.1D0)
  1354.                          (DELTAT * 1.1D0))           ;
  1355. *   DELTAT = DELTAT * CFFL                                             ;
  1356.  
  1357.    'SI' ('EGA' DELTAT 0.D0 1.D-15)                                    ;
  1358. *     GBM - on prend un très petit pas de temps pour initialiser les
  1359. *           flux - TESTER 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1360.       DELTAT = 1.D-15                                                 ;
  1361.    'FINSI'                                                            ;
  1362. 'FINSI'                                                               ;
  1363. 'SI' daut                                                             ;
  1364.    'MESSAGE'  'Pas de temps calculé automatiquement.'                 ;
  1365. 'FINSI'                                                               ;
  1366.  
  1367.  
  1368. 'RESPRO'  DELTAT DTI                                                  ;
  1369.  
  1370. 'FINSI'                                                               ;
  1371.  
  1372.  
  1373. ***********************************************************************
  1374.  
  1375. 'SI' ('EGA' roro YNYTYAL)                                            ;
  1376.  
  1377.  
  1378. 'ARGUMENT'   SATUR*'TABLE' TETA*'FLOTTANT'                            ;
  1379.  
  1380. SATUR. 'PRESSION'   = 'TABLE'                                         ;
  1381. SATUR. 'TH2O'       = 'TABLE'                                         ;
  1382. SATUR. 'SATURATION' = 'TABLE'                                         ;
  1383. * Liste des temps où il aura fallu user de pénalisation :
  1384. SATUR.'TEMPS_XI'      = 'PROG'                                        ;
  1385. * sauvegarde du nombre d'itérations pour chaque pas de temps :
  1386. SATUR.'ITERPICARD'    = 'LECT' 0                                      ;
  1387. * Liste des temps calculés a posteriori
  1388. SATUR.'TEMPS_EFFECTUES' = 'PROG'                                      ;
  1389.  
  1390. * Initialisation des tables la proc TRANGEOL
  1391. * Indicateur les recalculs de matrices à effectuer
  1392. * A virer quand procedure 1 pas incluse
  1393. TABMODI = TABLE                                                      ;
  1394. TABMODI . 'POROSITE' = VRAI                                          ;
  1395. TABMODI . 'CONVECTI' = VRAI                                          ;
  1396. TABMODI . 'DELTAT'   = VRAI                                          ;
  1397. TABMODI . 'COEF_LIN' = VRAI                                          ;
  1398. TABMODI . 'DIFFUSIV' = VRAI                                          ;
  1399. *
  1400. CHCLIM = 'TABLE'                                                     ;
  1401. GEOL1 = 'TABLE'                                                      ;
  1402. * l'option abandon indique qu'en dessous de 10**-13 pour la concentr
  1403. * et pour des sources ou flux n'impliquant pas des variation de C
  1404. * superieurs à 10**-16 on sort 0 sans calcul
  1405. 'SI' ('EXISTE' SATUR 'SEUILCALC')                                    ;
  1406. GEOL1 . 'ABANDON'           = VRAI                                   ;
  1407. GEOL1 . 'SEUILCALC'         = SATUR . 'SEUILCALC'                    ;
  1408. 'SINON'                                                              ;
  1409. GEOL1 . 'ABANDON'           = FAUX                                   ;
  1410. GEOL1 . 'SEUILCALC'         = 1.D-30                                 ;
  1411. 'FINSI'                                                              ;
  1412.  
  1413. GEOL1 . 'NUM_PECLET'        = 2.D0                                   ;
  1414. 'SI' ('EXISTE' SATUR LUMP)                                           ;
  1415. GEOL1 . 'LUMP'              = SATUR . 'LUMP'                         ;
  1416. 'SINON'                                                              ;
  1417. GEOL1 . 'LUMP'              = FAUX                                   ;
  1418. 'FINSI'                                                              ;
  1419.  
  1420. 'SI' ('EXISTE' SATUR 'TYPDISCRETISATION' )                           ;
  1421. GEOL1 . 'TYPDISCRETISATION' = SATUR . 'TYPDISCRETISATION'            ;
  1422. 'SINON'                                                              ;
  1423. GEOL1 . 'TYPDISCRETISATION' = 'VF'                                   ;
  1424. 'FINSI'                                                              ;
  1425.  
  1426. GEOL1 . 'THETA_DIFFUSION'   = TETA                                   ;
  1427. * GBM inutile car pas de convection
  1428. GEOL1 . 'THETA_CONVECTION'  = 1.D0                                   ;
  1429. GEOL1 . 'DECENTREMENT'      = FAUX                                   ;
  1430.  
  1431. GEOL1 . 'MODIFICATI'        = TABMODI                                ;
  1432.  
  1433.  
  1434. 'RESPRO'  CHCLIM GEOL1 TABMODI                                       ;
  1435.  
  1436.  
  1437. 'FINSI'                                                               ;
  1438.  
  1439. ***********************************************************************
  1440.  
  1441. 'FINPROC'                                                             ;
  1442.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales