Download satutils.procedur

Back to the list

   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   '<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 '<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 : 
1443 :  
1444 :  

© Cast3M 2003 - All rights reserved.
Disclaimer