Télécharger darcytra.procedur

Retour à la liste

* DARCYTRA  PROCEDUR  JC220346  14/02/19    21:15:02     7941           
'DEBP' DARCYTRA TRANSI*'TABLE';
 
*---------------------------------------------------------------------
* Résolution de l'équation de darcy pour un problème d'écoulement ou
* de transport par une méthode d'éléments
* finis mixtes hybrides. Les inconnues du problème sont
*  - pour l'écoulement, la charge (H) et trace de charge (TH) et le
*    débit hydraulique (FLUX)
*  - pour le transport, la concentration (H), la trace de
*    concentration (TH) et le débit diffusif (FLUX).
*---------------------------------------------------------------------
*
*------------------------------
* Phrase d'appel (en GIBIANE) :
*------------------------------
*
*  DARCYTRA TABLE ;
*
*----------------------------------
* Opérandes (à mettre dans TABLE) :
*----------------------------------
*
*  ___________________________________________________________________
*  |                                                                 |
*  | Indice      Contenu                                             |
*  |                                                                 |
*  -------------------------------------------------------------------
*  |                                                                 |
*  |------------------------------------------------                 |
*  |Données physiques, géométriques et materielles :                 |
*  |------------------------------------------------                 |
*  |                                                                 |
*  |   ------ Indices communs à l'écoulement et au transport ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'SOUSTYPE'   'DARCY' (type MOT)                                  |
*  |                                                                 |
*  |'MODELE'      Objet modèle (MMODEL créé par MODE)                |
*  |                                                                 |
*  |                                                                 |
*  |   ------ 1ère possibilité : Résolution de l'écoulement ------   |
*  |   -----------------------------------------------------------   |
*  |                                                                 |
*  |'CARACTERISTIQUES' Données physiques et materielles :            |
*  |              conductivité hydraulique (CHAMELEM créé par MATE)  |
*  |                                                                 |
*  |'EMMAGASINEMENT' Valeur du coefficient d'emmagasinement          |
*  |              (Type CHPO Centre, Comp 'CK', ou FLOTTANT)         |
*  |              - Défaut 1.                                        |
*  |                                                                 |
*  |   ------ 2ème possibilité : Résolution du transport     ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'CARACTERISTIQUES' Données physiques et materielles :            |
*  |              diffusivité effective (CHAMELEM créé par MATE)     |
*  |                                                                 |
*  |'POROSITE'    Valeur de la porosité (Type CHAMELEM, Comp         |
*  |              'CK', ou FLOTTANT) - Défaut 1.                     |
*  |                                                                 |
*  |'DECROISSANCE' Valeur du terme de décroissance (Type FLOTTANT)   |
*  |              Tel que dC/dt = - Lambda * C     - Défaut 0.       |
*  |                                                                 |
*  |'COEF_RETARD' Coefficient de retard linéaire dans le cas simple, |
*  |              ou  Pente à l'origine de la fonction F(C) dans le  |
*  |                  cas d'isotherme non linéaire de Langmuir       |
*  |              ou  Coefficient K de l'isotherme de Freundlich     |
*  |              (Type CHPO Centre 'SCAL', ou FLOTTANT)             |
*  |                                                                 |
*  |'LANGMUIR'    Quantité maximale 'Fsat' adsorbée sur le solide    |
*  |              rapportée à l'unité de volume du fluide et exprimée|
*  |              dans la meme unité que le soluté.                  |
*  |              (Type CHPO Centre 'SCAL', ou FLOTTANT).            |
*  |              F = (R-1) C / [1 + ((R-1) C / Fsat)]               |
*  |              Si cet indice et le suivant sont absents,          |
*  |              l'équilibre d'adsorption est linéaire. Cet indice a|
*  |              priorité sur l'indice FREUNDLICH.                  |
*  |                                                                 |
*  |'FREUNDLICH'  Exposant de la loi de Freundlich F = K (C ^ 1/n)   |
*  |              (Type FLOTTANT).                                   |
*  |              Dans ce cas (et si l'indice LANGMUIR n'existe pas),|
*  |              l'indice 'COEF_RETARD' contient le coefficient     |
*  |              K ramené à une unité de volume de fluide.          |
*  |              - Non disponible pour l'instant -                  |
*  |                                                                 |
*  |'LIMITE_SOLUBILITE' Limite de solubilité (Type MCHAML, Comp 'H') |
*  |                                                                 |
*  |'COEF_DISSOLUTION' Coef. de dissolution (Type CHPO Centre, Comp  |
*  |              'SCAL'). Tel que dC/dt = Coef * (Csat - C) - Par   |
*  |              défaut, la dissolution est instantanée             |
*  |                                                                 |
*  |'CONVECTION'  Débit intégré de la vitesse convective à travers   |
*  |              chaque face des éléments (Type CHPO Face, comp.    |
*  |              'FLUX')                                            |
*  |                                                                 |
*  |----------------------                                           |
*  |Conditions initiales :                                           |
*  |----------------------                                           |
*  |                                                                 |
*  |   ------ Indices communs à l'écoulement et au transport ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'TEMPS'       TABLE contenant à l'indice 0 la valeur du temps    |
*  |              initial (FLOTTANT)                                 |
*  |                                                                 |
*  |                                                                 |
*  |   ------ 1ère possibilité : Résolution de l'écoulement ------   |
*  |   -----------------------------------------------------------   |
*  |                                                                 |
*  |'CHARGE'      TABLE contenant à l'indice 0 la charge hydraulique |
*  |              (quantité d'élément par unité de volume d'eau)     |
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'TRACE_CHARGE' TABLE contenant à l'indice 0 la trace de          |
*  |              charge initiale (CHPO, 'TH')                       |
*  |                                                                 |
*  |'FLUX'        TABLE contenant à l'indice 0 le flux hydraulique   |
*  |              initial intégré sur chaque face (Type CHPO Face,   |
*  |              comp. 'FLUX')                                      |
*  |                                                                 |
*  |   ------ 2ème possibilité : Résolution du transport     ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'CONCENTRATION' TABLE contenant à l'indice 0 la concentration    |
*  |              (quantité d'élément par unité de volume d'eau)     |
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'TRACE_CONC'  TABLE contenant à l'indice 0 la trace de           |
*  |              concentration initiale (Type CHPO Face, Comp 'TH') |
*  |                                                                 |
*  |'FLUX'        TABLE contenant à l'indice 0 le flux diffusif      |
*  |              initial intégré sur chaque face (Type CHPO Face,   |
*  |              comp. 'FLUX')                                      |
*  |                                                                 |
*  |'PRECIPITE'   TABLE contenant à l'indice 0 la quantité initiale  |
*  |              de précipité par unité de volume de milieu         |
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'DISSOLUTION' TABLE contenant à l'indice 0 la quantité initiale  |
*  |              pour estimer la dissolution au premier pas de temps|
*  |              (Type CHPO, Comp 'H'), voir plus loin.             |
*  |                                                                 |
*  |--------------------------------------                           |
*  |Conditions aux limites / chargements :                           |
*  |--------------------------------------                           |
*  |                                                                 |
*  |   ------ Indices communs à l'écoulement et au transport ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'BLOCAGE'     Contient les matrices de blocage (RIGIDITE)        |
*  |                                                                 |
*  |'TRACE_IMPOSE' Valeurs des traces imposées (charge ou concentra- |
*  |              -tion) (CHARGEMENT 'TH' - Obligatoire si BLOCAGE)  |
*  |                                                                 |
*  |'FLUX_IMPOSE' Valeurs des flux imposés intégrés par face         |
*  |              (Type CHARGEMENT de CHPO Face, comp. 'FLUX'-       |
*  |              défaut 0.)                                         |
*  |                                                                 |
*  |'SOURCE'      Valeurs du terme source par maille et par unité de |
*  |              temps (ex : puits, filiation)                      |
*  |              Les valeurs à l'indice i sont les valeurs entre    |
*  |              les temps i-1 et i.                                |
*  |              (CHARGEMENT de CHPO Centre, comp 'SOUR'- défaut 0.)|
*  |                                                                 |
*  |   ------ 2ème possibilité : Résolution du transport     ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'DISSOLUTION_IMPOSEE' Valeurs des dissolutions imposées par unité|
*  |              de temps et par maille. (Type CHARGEMENT de CHPO,  |
*  |              Comp 'H'). Les valeurs à l'indice i sont les       |
*  |              valeurs moyennes de dissolution par unité de temps |
*  |              entre les temps i-1 et i.                          |
*  |              Priorité de la dissolution imposée sur les         |
*  |              cinétiques.                                        |
*  |                                                                 |
*  |--------------------                                             |
*  |Données numériques :                                             |
*  |--------------------                                             |
*  |                                                                 |
*  |   ------ Indices communs à l'écoulement et au transport ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'TEMPS_CALCULES' Valeur des temps calculés (LISTREEL)            |
*  |              Contient obligatoirement le temps final.           |
*  |                                                                 |
*  |'TEMPS_SAUVES' Valeur des temps sauvegardés (LISTREEL - défaut : |
*  |              on sauve tous les pas de temps)                    |
*  |                                                                 |
*  |   ------ 1ère possibilité : Résolution de l'écoulement ------   |
*  |   -----------------------------------------------------------   |
*  |                                                                 |
*  |'THETA'       Coefficient de relaxation compris entre 0. et 1.   |
*  |              (theta-méthode diffusion) (FLOTTANT - défaut 1.)   |
*  |              Possibilité de non-convergence lorsque theta<1/2   |
*  |              Valeurs de theta généralement utilisées :          |
*  |                   Schéma de Euler explicite : 0.                |
*  |                   Schéma de Crank-Nicholson : 1/2               |
*  |                   Schéma de Galerkin        : 2/3               |
*  |                   Schéma de Euler implicite : 1.                |
*  |                                                                 |
*  |   ------ 2ème possibilité : Résolution du transport     ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'THETA_DIFF'  Coefficient de relaxation compris entre 0. et 1.   |
*  |              (theta-méthode diffusion) ('FLOTTANT' - défaut 1.) |
*  |                                                                 |
*  |'THETA_CONVECTION' Idem pour la convection                       |
*  |              ('FLOTTANT', Défaut = THETA_DIFF)                  |
*  |                                                                 |
*  |'THETA_DEC'   Idem mais pour la décroissance                     |
*  |              ('FLOTTANT' - défaut 1/2)                          |
*  |                                                                 |
*  |'THETA_DISS'   Idem mais pour la dissolution                     |
*  |              ('FLOTTANT' - défaut 1.)                           |
*  |                                                                 |
*  |'PENALISATION' Coefficient de pénalisation pour la prise en      |
*  |              compte de la limite de solubilité. La présence de  |
*  |              cet indice ou du suivant indique quel schéma a été |
*  |              choisi.                                            |
*  |              (Type 'FLOTTANT') - Valeur conseillée 1.D7         |
*  |                                                                 |
*  |'EPSI_LIM'    Précision relative d'arrêt pour le shéma limite de |
*  |              solubilité prédicteur-correcteur itératif          |
*  |              (Type FLOTTANT) - Valeur conseillée 5.D-3          |
*  |                                                                 |
*  |'ITMAX_LIM'   Nombre maxi d'itérations correspondant aux modules |
*  |              de dissolution avant d'abandonner                  |
*  |              (Type 'ENTIER') - Défaut 50                        |
*  |                                                                 |
*  |'EPSI_RET'    Précision relative d'arrêt pour la résolution      |
*  |              itérative (Picard) de l'adsorption non linéaire    |
*  |              (Type FLOTTANT) - Défaut 1.D-4                     |
*  |                                                                 |
*  |'EPSI_COR'    Petit saut de concentration pour calculer le coef. |
*  |              de retard par la méthode de la corde lorsque le    |
*  |              retard est non-linéaire.                           |
*  |              (Type FLOTTANT) - Défaut 1.D-4                     |
*  |                                                                 |
*  |'ITMAX_RET'   Nombre maxi d'itérations correspondant au retard   |
*  |              non linéaire avant d'abandonner.                   |
*  |              (Type 'ENTIER') - Défaut 20                        |
*  |_________________________________________________________________|
*
*
*
*---------------------------------
* Résultats (stockés dans TABLE) :
*---------------------------------
*
*  ___________________________________________________________________
*  |                                                                 |
*  | Indice      Contenu                                             |
*  |                                                                 |
*  -------------------------------------------------------------------
*  |                                                                 |
*  |   ------ Indices communs à l'écoulement et au transport ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'TEMPS'       TABLE contenant les temps sauvegardés (FLOTTANT)   |
*  |                                                                 |
*  |   ------ 1ère possibilité : Résolution de l'écoulement  ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'CHARGE'      TABLE contenant les charges                        |
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'TRACE_CHARGE' TABLE contenant les traces de charge              |
*  |              (Type CHPO Face, Comp 'TH')                        |
*  |                                                                 |
*  |'FLUX'        TABLE contenant les débits hydrauliques intégrés   |
*  |              par face (Type CHPO Face, comp. 'FLUX')            |
*  |                                                                 |
*  |   ------ 2ème possibilité : Résolution du transport     ------  |
*  |   ------------------------------------------------------------  |
*  |                                                                 |
*  |'CONCENTRATION' TABLE contenant les concentrations               |
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'TRACE_CONC'  TABLE contenant les traces de concentration        |
*  |              (Type CHPO Face, Comp 'TH')                        |
*  |                                                                 |
*  |'FLUX'        TABLE contenant les débits diffusifs intégrés      |
*  |              par face (Type CHPO Face, comp. 'FLUX')            |
*  |                                                                 |
*  |'PRECIPITE'   TABLE contenant la quantité de précipité par maille|
*  |              (Type CHPO Centre, Comp 'H')                       |
*  |                                                                 |
*  |'DISSOLUTION' TABLE contenant la quantité de précipité dissoute  |
*  |              entre deux pas de temps par unité de volume et par |
*  |              unité de temps. La valeur stockée à l'indice i,    |
*  |              est valable entre les temps i-1 et i               |
*  |              (Type CHPO, Comp 'H').                             |
*  |              ATTENTION, les valeurs de cette table résultat     |
*  |              n'ont aucun sens lorsque les temps sauvegardés ne  |
*  |              sont pas les memes que les temps calculés. Toute   |
*  |              tentative d'exploitation donnera alors des         |
*  |              résultats incohérents (erreurs de bilan)           |
*  |                                                                 |
*  |'RETARD'      Si cet indice a été préalablement défini comme une |
*  |              TABLE, alors il contient les valeurs du coefficient|
*  |              de retard (Type 'CHPO' centre, Comp 'SCAL'). Sinon,|
*  |              les valeurs du coefficient de retard ne sont pas   |
*  |              sauvegardées.                                      |
*  |_________________________________________________________________|
*
*
*  ___________________________________________________________________
*  |                                                                 |
*  | Les tables résultats sont indicées par des entiers variant de 0 |
*  | à N .                                                           |
*  | A l'indice 0 on stocke les valeurs initiales, aux indices       |
*  | suivants les champs correspondant au temps de sortie TEMPS.I .  |
*  | Les champs servant en cas de reprise sont ceux correpondant au  |
*  | dernier indice.                                                 |
*  |_________________________________________________________________|
*
*
*                                                             SOUSTYPE
'SI' ( 'NON' ('EXISTE' TRANSI 'SOUSTYPE' ) )                         ;
     'ERREUR'
     'Indice SOUSTYPE absent de la table de données.'                ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
'SI' ( 'NEG' ( TRANSI.'SOUSTYPE' ) 'DARCY_TRANSITOIRE' )             ;
     'ERREUR' 'Le sous-type de la table d entree est incorrect'      ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                   TYPE DE RESOLUTION
* Le nom des variables suit la logique du transport.
* Par la suite, si la résolution est hydro, la porosité tiendra lieu
* de coef. d'emmagasinement, la concentration de charge (idem pour les
* traces), et theta-diff de theta.
* Pour chaque variable spécifique, on vérifie que sa présence
* correspond bien au type de résolution envisagé.
HYDRO  = 'EXISTE' TRANSI 'CHARGE'                                    ;
TRANSP ='NON' HYDRO                                                  ;
'SI' HYDRO                                                           ;
     NOMINC  = 'MOT' 'CHARGE'                                        ;
     NOMTINC = 'MOT' 'TRACE_CHARGE'                                  ;
     NOMTETA = 'MOT' 'THETA'                                         ;
     NOMDDT  = 'MOT' 'EMMAGASINEMENT'                                ;
     'SI' ('EXISTE' TRANSI 'POROSITE')                               ;
          'ERREUR' ('CHAINE' 'L indice POROSITE ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'DECROISSANCE')                           ;
          'ERREUR' ('CHAINE' 'L indice DECROISSANCE ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'COEF_RETARD')                            ;
          'ERREUR' ('CHAINE' 'L indice COEF_RETARD ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'LANGMUIR')                               ;
          'ERREUR' ('CHAINE' 'L indice LANGMUIR ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'FREUNDLICH')                             ;
          'ERREUR' ('CHAINE' 'L indice FREUNDLICH ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'LIMITE_SOLUBILITE')                      ;
          'ERREUR' ('CHAINE' 'L indice LIMITE_SOLUBILITE ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'COEF_DISSOLUTION')                       ;
          'ERREUR' ('CHAINE' 'L indice COEF_DISSOLUTION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'CONVECTION')                             ;
          'ERREUR' ('CHAINE' 'L indice CONVECTION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'CONCENTRATION')                          ;
          'ERREUR' ('CHAINE' 'L indice CONCENTRATION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'TRACE_CONC')                             ;
          'ERREUR' ('CHAINE' 'L indice TRACE_CONC ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'PRECIPITE')                              ;
          'ERREUR' ('CHAINE' 'L indice PRECIPITE ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'DISSOLUTION')                            ;
          'ERREUR' ('CHAINE' 'L indice DISSOLUTION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE')                    ;
        'ERREUR' ('CHAINE' 'L indice DISSOLUTION_IMPOSEE ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'THETA_DIFF')                             ;
        'ERREUR' ('CHAINE' 'L indice THETA_DIFF ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'THETA_CONVECTION')                       ;
        'ERREUR' ('CHAINE' 'L indice THETA_CONVECTION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'THETA_DEC')                              ;
        'ERREUR' ('CHAINE' 'L indice THETA_DEC ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'THETA_DISS')                             ;
        'ERREUR' ('CHAINE' 'L indice THETA_DISS ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'PENALISATION')                           ;
        'ERREUR' ('CHAINE' 'L indice PENALISATION ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'EPSI_LIM')                               ;
        'ERREUR' ('CHAINE' 'L indice EPSI_LIM ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'EPSI_RET')                               ;
        'ERREUR' ('CHAINE' 'L indice EPSI_RET ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'EPSI_COR')                               ;
        'ERREUR' ('CHAINE' 'L indice EPSI_COR ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'ITMAX_RET')                              ;
        'ERREUR' ('CHAINE' 'L indice ITMAX_RET ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'RETARD')                                 ;
        'ERREUR' ('CHAINE' 'L indice RETARD ne correspond '
                    'pas à la résolution hydraulique')               ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
'SI' TRANSP                                                          ;
     NOMINC  = 'MOT' 'CONCENTRATION'                                 ;
     NOMTINC = 'MOT' 'TRACE_CONC'                                    ;
     NOMTETA = 'MOT' 'THETA_DIFF'                                    ;
     NOMDDT  = 'MOT' 'POROSITE'                                      ;
     'SI' ('EXISTE' TRANSI 'THETA')                                  ;
          'ERREUR' ('CHAINE' 'L indice THETA ne correspond '
                    'pas à la résolution du transport')              ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'TRACE_CHARGE')                           ;
          'ERREUR' ('CHAINE' 'L indice TRACE_CHARGE ne correspond '
                    'pas à la résolution du transport')              ;
     'FINSI'                                                         ;
     'SI' ('EXISTE' TRANSI 'EMMAGASINEMENT')                         ;
          'ERREUR' ('CHAINE' 'L indice EMMAGASINEMENT ne correspond '
                    'pas à la résolution du transport')              ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*
*--------------------------------------------------------------------*
* RECUPERATION DES DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES    *
*--------------------------------------------------------------------*
*
*                                                               MODELE
'SI' ( 'EXISTE' TRANSI 'MODELE' )                                    ;
     MODHYB = TRANSI . 'MODELE'                                      ;
     HYTOT = DOMA MODHYB 'TABLE' ;
'SINON'                                                              ;
     'ERREUR' 'Il manque le modele.'                                 ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                               VOLUME
VOLU1 = 'DOMA' MODHYB 'VOLUME'                                      ;
*                                                          ORIENTATION
MCHYB = 'DOMA' MODHYB 'ORIENTAT'                                    ;
*                                                     CARACTERISTIQUES
'SI' ( 'EXISTE' TRANSI 'CARACTERISTIQUES' )                          ;
     MAT1 = TRANSI . 'CARACTERISTIQUES'                              ;
'SINON'                                                              ;
     'ERREUR'
       'Indice CARACTERISTIQUES absent de la table de données.'      ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                           EMMAGASINEMENT OU POROSITE
'SI' ( 'EXISTE' TRANSI NOMDDT )                                      ;
     POROS = TRANSI . NOMDDT                                         ;
'SINON'                                                              ;
     POROS = 1.D0                                                    ;
'FINSI'                                                              ;
'SI' ('EGA' ('TYPE' POROS) 'FLOTTANT')                               ;
     POROS = 'MANU' 'CHML' HYTOT.'MAILLAGE' 'CK' POROS         ;
'FINSI'                                                              ;
'SI' ('NEG' ('TYPE' POROS) 'MCHAML  ')                               ;
     'ERREUR' ('CHAINE' 'L indice ' NOMDDT ' doit etre de type '
               'FLOTTANT ou CHAMELEM')                               ;
'FINSI'                                                              ;
'SI' HYDRO                                                           ;
     'SI' (('MINI' POROS) < 0.D0 )                                   ;
          'ERREUR' 'Le coef. d emmagasinement doit etre positif'     ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
'SI' TRANSP                                                          ;
     'SI' ((('MINI' POROS) < 0.D0) 'OU' (('MAXI' POROS) > 1.D0))     ;
          'ERREUR' 'La porosité doit etre comprise entre 0 et 1'     ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
POROS1 = 'KCHA' MODHYB POROS 'CHPO'                                   ;
*                                                           CONVECTION
'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
     CONV1 = TRANSI . 'CONVECTION'                                   ;
'FINSI'                                                              ;
*                                                               RETARD
'SI' ( 'EXISTE' TRANSI 'COEF_RETARD' )                               ;
     RETAR1 = TRANSI . 'COEF_RETARD'                                 ;
     'SI' ('EGA' ('TYPE' RETAR1) 'CHPOINT ')                         ;
          MINRET = ('MINI' RETAR1)                                   ;
     'SINON'                                                         ;
          MINRET = RETAR1                                            ;
     'FINSI'                                                         ;
     'SI' (MINRET < 1.D0)                                            ;
          'ERREUR' 'Le coefficient de retard doit etre supérieur à 1';
     'FINSI'                                                         ;
'SINON'                                                              ;
     RETAR1 = 1.                                                     ;
'FINSI'                                                              ;
RETAR1M1 = RETAR1 - 1.                                               ;
SAUVRET = 'EXISTE' TRANSI 'RETARD'                                   ;
*                                                  RETARD NON LINEAIRE
RNONLINL = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
           ('EXISTE' TRANSI 'LANGMUIR')                              ;
RNONLINF = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
           ('EXISTE' TRANSI 'FREUNDLICH')  'ET' ('NON' RNONLINL)     ;
RNONLIN  = RNONLINL 'OU' RNONLINF                                    ;
* Cas Langmuir :
'SI' RNONLINL                                                        ;
     FSAT = TRANSI . 'LANGMUIR'                                      ;
     'SI' ('EGA' ('TYPE' FSAT) 'CHPOINT ')                           ;
          'SI' (('MINI' FSAT) < 0.D0)                                ;
            'ERREUR'
              'La masse adsorbee a saturation doit etre positive'    ;
            'QUITTER' DARCYTRA                                       ;
          'FINSI'                                                    ;
     'SINON'                                                         ;
          'SI' (FSAT < 0.D0)                                         ;
            'ERREUR'
              'La masse adsorbee a saturation doit etre positive'    ;
            'QUITTER' DARCYTRA                                       ;
          'FINSI'                                                    ;
     'FINSI'                                                         ;
     RM1SURF  = RETAR1M1 / FSAT                                      ;
'FINSI'                                                              ;
* Cas Freundlich :
'SI' RNONLINF                                                        ;
     UNSURN = TRANSI . 'FREUNDLICH'                                  ;
     'SI' ((UNSURN &lt;EG  0.D0) 'OU' (UNSURN > 1.D0))                  ;
          'ERREUR'
            'L exposant de Freundlich doit etre compris entre 0'
            ' (exclu) et 1'                                          ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                                         DECROISSANCE
DECRO = ( 'EXISTE' TRANSI 'DECROISSANCE' )                           ;
'SI' DECRO                                                           ;
     LAMBD0 = TRANSI . 'DECROISSANCE'                                ;
     'SI' (LAMBD0 < 0.D0)                                            ;
          'ERREUR' 'Le coefficient de décroissance doit etre positif.';
     'FINSI'                                                         ;
'SINON'                                                              ;
     LAMBD0 = 0.                                                     ;
'FINSI'                                                              ;
*                                                  DISSOLUTION_IMPOSEE
SOLUA = ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE' )                     ;
'SI' SOLUA                                                           ;
     DISIMP = TRANSI . 'DISSOLUTION_IMPOSEE'                         ;
'FINSI'                                                              ;
*                                                    LIMITE_SOLUBILITE
* (Priorité de la dissolution imposée sur les autres processus)
SOLUL = (('EXISTE' TRANSI 'LIMITE_SOLUBILITE' ) 'ET' ('NON' SOLUA))  ;
'SI' SOLUL                                                           ;
     LIMSOL = TRANSI . 'LIMITE_SOLUBILITE'                           ;
     'SI' (('MINI' LIMSOL) < 0.D0)                                   ;
          'ERREUR' 'La limite de solubilité doit etre positive.'     ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                                     COEF_DISSOLUTION
SOLUP = SOLUL 'ET' ('EXISTE' TRANSI 'COEF_DISSOLUTION' )             ;
SOLUI = SOLUL 'ET' ('NON' SOLUP)                                     ;
'SI' SOLUP                                                           ;
     CODIS = TRANSI . 'COEF_DISSOLUTION'                             ;
     'SI' ('EGA' ('TYPE' CODIS) 'CHPOINT ')                          ;
          'SI' (('MINI' CODIS) < 0.D0)                               ;
            'ERREUR' 'Le coefficient de dissolution doit etre positif.';
          'FINSI'                                                    ;
     'SINON'                                                         ;
          'SI' (CODIS < 0.D0)                                        ;
            'ERREUR' 'Le coefficient de dissolution doit etre positif.';
          'FINSI'                                                    ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
SOLUB = SOLUL 'OU' SOLUA                                             ;
*
* ------------------------------------------------------------------ *
*  VERIFICATION DE LA BONNE STRUCTURE DES TABLES RESULTAT D'ORIGINE  *
* ------------------------------------------------------------------ *
*
*                                                                TEMPS
'SI' ( 'NON' ('EXISTE' TRANSI 'TEMPS' ) )                            ;
     'ERREUR'
       'Indice TEMPS absent de la table de données.'                 ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                       TRACE_INCONNUE
'SI' ( 'NON' ('EXISTE' TRANSI NOMTINC ) )                            ;
     'ERREUR' ('CHAINE'
       'Indice ' NOMTINC ' absent de la table de données.')          ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                             INCONNUE
'SI' ( 'NON' ('EXISTE' TRANSI NOMINC ) )                             ;
     'ERREUR' ('CHAINE'
       'Indice ' NOMINC ' absent de la table de données.')           ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                                 FLUX
'SI' ( 'NON' ('EXISTE' TRANSI 'FLUX' ) )                             ;
     'ERREUR'
       'Indice FLUX absent de la table de données.'                  ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                               DISSOLUTION, PRECIPITE
'SI' SOLUB                                                           ;
     'SI' ( 'NON' ('EXISTE' TRANSI 'DISSOLUTION' ) )                 ;
          'ERREUR'
            'Indice DISSOLUTION absent de la table de données.'      ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     'SI' ( 'NON' ('EXISTE' TRANSI 'PRECIPITE' ) )                   ;
          'ERREUR'
            'Indice PRECIPITE absent de la table de données.'        ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                            TEST DES TAILLES DE TABLE
IND1 = 'INDEX' ( TRANSI . 'TEMPS' )                                  ;
IND2 = 'INDEX' ( TRANSI . NOMTINC )                                  ;
IND3 = 'INDEX' ( TRANSI . NOMINC )                                   ;
IND4 = 'INDEX' ( TRANSI . 'FLUX' )                                   ;
LIN1 = 'DIME' IND1                                                   ;
LIN2 = 'DIME' IND2                                                   ;
LIN3 = 'DIME' IND3                                                   ;
LIN4 = 'DIME' IND4                                                   ;
'SI' ( LIN1 'NEG' LIN2 )                                             ;
     'ERREUR' ('CHAINE'
          'Longueur des tables TEMPS et ' NOMTINC ' différente.')    ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
'SI' ( LIN1 'NEG' LIN3 )                                             ;
     'ERREUR' ('CHAINE'
          'Longueur des tables TEMPS et ' NOMINC ' différente.')     ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
'SI' ( LIN1 'NEG' LIN4 )                                             ;
     'ERREUR'
       'Longueur des tables TEMPS et FLUX différente.'               ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
'SI' SOLUB                                                           ;
     IND5 = 'INDEX' ( TRANSI . 'DISSOLUTION' )                       ;
     IND6 = 'INDEX' ( TRANSI . 'PRECIPITE' )                         ;
     LIN5 = 'DIME' IND5                                              ;
     LIN6 = 'DIME' IND6                                              ;
     'SI' ( LIN1 'NEG' LIN5 )                                        ;
          'ERREUR'
            'Longueur des tables TEMPS et DISSOLUTION différente.'   ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     'SI' ( LIN1 'NEG' LIN6 )                                        ;
          'ERREUR'
             'Longueur des tables TEMPS et PRECIPITE différente.'    ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                            TEST DES INDICES DE TABLE
IPO1 = 0                                                             ;
'REPETER' BOU1 LIN1                                                  ;
     IPO1  = IPO1 + 1                                                ;
     LAST1 = IND1 . IPO1                                             ;
     LAST2 = IND2 . IPO1                                             ;
     LAST3 = IND3 . IPO1                                             ;
     LAST4 = IND4 . IPO1                                             ;
     'SI' SOLUB                                                      ;
          LAST5 = IND5 . IPO1                                        ;
          LAST6 = IND6 . IPO1                                        ;
     'FINSI'                                                         ;
     'SI' ( 'NEG' LAST1 LAST2 )                                      ;
          'ERREUR'
            'Indices des tables TEMPS et TRACE_CONC incohérents.'    ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     'SI' ( 'NEG' LAST1 LAST3 )                                      ;
          'ERREUR'
            'Indices des tables TEMPS et CONCENTRATION incohérents.' ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     'SI' ( 'NEG' LAST1 LAST4 )                                      ;
          'ERREUR'
            'Indices des tables TEMPS et FLUX incohérents.'          ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     'SI' SOLUB                                                      ;
        'SI' ( 'NEG' LAST1 LAST5 )                                   ;
             'ERREUR'
               'Indices des tables TEMPS et DISSOLUTION incohérents.';
             'QUITTER' DARCYTRA                                      ;
        'FINSI'                                                      ;
        'SI' ( 'NEG' LAST1 LAST6 )                                   ;
             'ERREUR'
               'Indices des tables TEMPS et PRECIPITE incohérents.'  ;
             'QUITTER' DARCYTRA                                      ;
        'FINSI'                                                      ;
     'FINSI'                                                         ;
'FIN' BOU1                                                           ;
*
* ------------------------------------------------------------------ *
*  RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS SAUVE     *
* ------------------------------------------------------------------ *
*
TPSINI = TRANSI . 'TEMPS' . LAST1                                    ;
TRP    = TRANSI . NOMTINC . LAST1                                    ;
CHRG   = TRANSI . NOMINC . LAST1                                     ;
FLU0   = TRANSI . 'FLUX' . LAST1                                     ;
* On enlève les multiplicateurs de Lagrange éventuels
LIST1 = 'EXTR' TRP 'COMP'                                            ;
LLIS1 = 'DIME' LIST1                                                 ;
'SI' ( LLIS1 'NEG' 1 )                                               ;
     TRP = 'EXCO' TRP 'TH' 'TH'                                      ;
'FINSI'                                                              ;
*
*--------------------------------------------------------------------*
*  RECUPERATION DES CONDITIONS AUX LIMITES ET DES CHARGEMENTS        *
*--------------------------------------------------------------------*
*
*                                  CL de type Dirichlet :      BLOCAGE
*                                                         TRACE_IMPOSE
'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                                   ;
     MATBLOC = TRANSI . 'BLOCAGE'                                    ;
     'SI' ( 'EXISTE' TRANSI 'TRACE_IMPOSE')                          ;
          CHAIMP = TRANSI . 'TRACE_IMPOSE'                           ;
     'SINON'                                                         ;
          'ERREUR'
            'Valeurs de traces de concentration imposées absentes'   ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                     CL de type Neumann : FLUX_IMPOSE
'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' )                               ;
     FLUIMP = TRANSI . 'FLUX_IMPOSE'                                 ;
'FINSI'                                                              ;
'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
     TERSOU = TRANSI . 'SOURCE'                                      ;
'FINSI'                                                              ;
*
*--------------------------------------------------------------------*
* RECUPERATION DES DONNEES NUMERIQUES                                *
*--------------------------------------------------------------------*
*                                                      THETA DIFFUSION
  'SI' ( 'EXISTE' TRANSI NOMTETA )                                   ;
       TETA = TRANSI . NOMTETA                                       ;
  'SINON'                                                            ;
       TETA = 1.D0                                                   ;
  'FINSI'                                                            ;
*                                                     THETA CONVECTION
'SI' ( 'EXISTE' TRANSI 'THETA_CONVECTION' )                          ;
     TETAC = TRANSI . 'THETA_CONVECTION'                             ;
'SINON'                                                              ;
     TETAC = TETA                                                    ;
'FINSI'                                                              ;
*                                             PARAMETRES DE SOLUBILITE
'SI' (SOLUI 'OU' SOLUP)                                              ;
     'SI' ( 'EXISTE' TRANSI 'ITMAX_LIM')                             ;
          ITMAXI = TRANSI . 'ITMAX_LIM'                              ;
     'SINON'                                                         ;
          ITMAXI = 50                                                ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
'SI' SOLUI                                                           ;
     'SI' ( 'EXISTE' TRANSI 'PENALISATION')                          ;
          PENAL = TRANSI . 'PENALISATION'                            ;
          MET0  = 'PENALISATION'                                     ;
     'SINON'                                                         ;
          'SI' ( 'EXISTE' TRANSI 'EPSI_LIM')                         ;
               EPS1 = TRANSI . 'EPSI_LIM'                            ;
               MET0 = 'PF DISSOLUTION'                               ;
          'SINON'                                                    ;
               'MESS' 'Il manque le paramètre '                      ;
               'ERREUR' 'du schéma numérique de limite de solubilité';
          'FINSI'                                                    ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*                                    PARAMETRES DE RETARD NON LINEAIRE
'SI' RNONLIN                                                         ;
     'SI' ( 'EXISTE' TRANSI 'EPSI_RET' )                             ;
          EPSRNL = TRANSI . 'EPSI_RET'                               ;
     'SINON'                                                         ;
          EPSRNL = 1.D-4                                             ;
     'FINSI'                                                         ;
     'SI' ( 'EXISTE' TRANSI 'ITMAX_RET' )                            ;
          ITMAXRNL = TRANSI . 'ITMAX_RET'                            ;
     'SINON'                                                         ;
          ITMAXRNL = 20                                              ;
     'FINSI'                                                         ;
*    Petit saut de concentration minimal pour le calcul de dérivée
     'SI' ( 'EXISTE' TRANSI 'EPSI_COR' )                             ;
          EPSCORD = TRANSI 'EPSI_COR'                                ;
     'SINON'                                                         ;
          'SI' RNONLINL                                              ;
               EPSCORD = ('MANU' 'CHPO' HYTOT.CENTRE 1 'SCAL' 1.D-4)
                       / RETAR1 * FSAT                               ;
          'FINSI'                                                    ;
          'SI' RNONLINF                                              ;
               EPSCORD = 1.D-4                                       ;
          'FINSI'                                                    ;
     'FINSI'                                                         ;
*    Pour les isothermes linéaires la boucle de prise en compte de
*    retard non linéaire n'est parcourue qu'une fois :
     NPicard = ITMAXRNL                                              ;
'SINON'                                                              ;
     NPicard = 1                                                     ;
'FINSI'                                                              ;
*                                                   THETA_DECROISSANCE
'SI' ( 'EXISTE' TRANSI 'THETA_DEC' )                                 ;
     BETA = TRANSI . 'THETA_DEC'                                     ;
'SINON'                                                              ;
     BETA = 0.5                                                      ;
'FINSI'                                                              ;
*                                                    THETA_DISSOLUTION
'SI' ( 'EXISTE' TRANSI 'THETA_DIS' )                                 ;
     GAMMA = TRANSI . 'THETA_DIS'                                    ;
'SINON'                                                              ;
     GAMMA = 1.                                                      ;
'FINSI'                                                              ;
*                                                       TEMPS_CALCULES
'SI' ( 'EXISTE' TRANSI 'TEMPS_CALCULES' )                            ;
     TPCAL  = ORDO (TRANSI.'TEMPS_CALCULES')                         ;
     TRANSI . 'TEMPS_CALCULES' = TPCAL                               ;
     DCAL   = 'DIME' TPCAL                                           ;
     TPSFIN = 'EXTR' TPCAL DCAL                                      ;
     TPS0   = 'EXTR' TPCAL 1                                         ;
     'SI' ('NEG' DCAL 1)                                             ;
          TPSANT = 'EXTR' TPCAL (DCAL - 1)                           ;
          TPS1   = 'EXTR' TPCAL 2                                    ;
          DT10   = TPS1 - TPS0                                       ;
          DT21   = TPSFIN - TPSANT                                   ;
          PROPS  = 'PROG' DT10 DT21                                  ;
          DTEPS  = 'MINIMUM' PROPS                                   ;
     'SINON'                                                         ;
          DTEPS  = TPSFIN - TPS0                                     ;
     'FINSI'                                                         ;
     EPS0  = DTEPS * 1.D-6                                           ;
*
*    ICAL est le premier indice des temps à calculer supérieur à
*    TpsIni + eps
     IOK1 = 0                                                        ;
     'REPETER' BOU2 DCAL                                             ;
          ICAL = &BOU2                                               ;
          TEMS = 'EXTR' TPCAL ICAL                                   ;
          'SI' ( TPSINI '<' (TEMS - EPS0) )                          ;
              IOK1 = 1                                               ;
              'QUITTER' BOU2                                         ;
          'FINSI'                                                    ;
     'FIN' BOU2                                                      ;
     'SI' ( IOK1 'EGA' 0)                                            ;
          'MESS' 'Listreel des temps de calcul inférieur à tini.'    ;
          'MESS' 'On ne fait donc rien.                         '    ;
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
     DELTAT0 = TEMS - TPSINI                                         ;
'SINON'                                                              ;
     'ERREUR'
       'Indice TEMPS_CALCULES absent de la table de données.'        ;
     'QUITTER' DARCYTRA                                              ;
'FINSI'                                                              ;
*                                                         TEMPS_SAUVES
'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' )                              ;
     TPSOR = ORDO (TRANSI . 'TEMPS_SAUVES')                          ;
     TRANSI . 'TEMPS_SAUVES' = TPSOR                                 ;
     DSOR   = 'DIME' TPSOR                                           ;
     ISOR = 0                                                        ;
     IOK1 = 0                                                        ;
     'REPETER' BOU3 DSOR                                             ;
          ISOR = ISOR + 1                                            ;
          TEMS = 'EXTR' TPSOR ISOR                                   ;
          'SI' ( ( TPSINI '<' (TEMS - EPS0) ) 'ET'
                 ( TEMS   '&lt;EG' (TPSFIN +EPS0)) )                    ;
               IOK1 = 1                                              ;
               'QUITTER' BOU3                                        ;
          'FINSI'                                                    ;
     'FIN' BOU3                                                      ;
     'SI' ( IOK1 'EGA' 0)                                            ;
          'ERREUR' 'Listreel des temps de sauvegarde hors tmin,tmax.';
          'QUITTER' DARCYTRA                                         ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*
*--------------------------------------------------------------------*
* RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE                       *
*--------------------------------------------------------------------*
*
'MESS' ' '                                                           ;
'MESS' 'MODELISATION Darcy EFMH en transitoire.'                     ;
'MESS' '---------------------------------------'                     ;
'MESS' ' '                                                           ;
'MESS' 'Donnees présentes en entrée : '                              ;
'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                                   ;
     'MESS' 'Calcul avec conditions aux limites sur Th (Dirichlet)'  ;
'FINSI'                                                              ;
'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' )                               ;
     'MESS' 'Calcul avec conditions aux limites de flux (Neumann)'   ;
'FINSI'                                                              ;
'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
     'MESS' 'Ce problème comporte un terme convectif'                ;
'FINSI'                                                              ;
'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
     'MESS' 'Ce problème comporte un terme source'                   ;
'FINSI'                                                              ;
'SI' DECRO                                                           ;
     'MESS' 'Ce problème comporte un terme de décroissance'          ;
'FINSI'                                                              ;
'SI' RNONLIN                                                         ;
     'MESS' 'Ce problème comporte la prise en compte d isothermes '
            'd adsorption non-linéaires'                             ;
     'SI' RNONLINL                                                   ;
          'MESS' '  de type Langmuir (méthode de Picard)'            ;
     'FINSI'                                                         ;
     'SI' RNONLINF                                                   ;
          'MESS' '  de type Freundlich (méthode de Picard)'          ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
'SI' SOLUL                                                           ;
     'MESS' 'Ce problème comporte une limite de solubilité'          ;
     'SI' SOLUP                                                      ;
          'MESS' '  avec dissolution progressive d ordre 1'          ;
     'FINSI'                                                         ;
     'SI' (EGA MET0 'PENALISATION')                                  ;
          'MESS' '  avec coefficient de pénalisation = ' PENAL       ;
     'FINSI'                                                         ;
     'SI' (EGA MET0 'PF DISSOLUTION')                                ;
          'MESS' '  avec critère de convergence = ' EPS1             ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
'SI' SOLUA                                                           ;
     'MESS' 'Ce probleme comporte une dissolution imposée'           ;
'FINSI'                                                              ;
'MESS' ' '                                                           ;
'MESS' 'Valeur des paramètres des schémas numériques :'              ;
'MESS' '  (0:Schéma explicite, 0.5:Crank-Nicholson, 1:Implicite)'    ;
'MESS' '  diffusion    : ' TETA                                      ;
'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
     'MESS' '  convection   : ' TETAC                                ;
'FINSI'                                                              ;
'SI' DECRO                                                           ;
     'MESS' '  décroissance : ' BETA                                 ;
'FINSI'                                                              ;
'SI' SOLUP                                                           ;
     'MESS' '  dissolution d ordre 1 : ' GAMMA                       ;
'FINSI'                                                              ;
'MESS' ' '                                                           ;
'MESS' 'Valeur du temps initial    : ' TPSINI                        ;
'MESS' 'Valeur du temps final      : ' TPSFIN                        ;
*
*
*--------------------------------------------------------------------*
* BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS               *
*--------------------------------------------------------------------*
*
*==============
* Préparation préliminaire de la boucle sur les pas de temps :
*==============
*
NBIT    = 0                                                          ;
NBITRNL = 0                                                          ;
MASHYB  = 'MHYB' MODHYB MAT1                                   ;
PRECED  = TPSINI                                                     ;
IPAS    = ICAL                                                       ;
DELOLD  = 0.D0                                                       ;
TPS1    = 'EXTR' TPCAL ICAL                                          ;
DELTAT  = TPS1 - TPSINI                                              ;
'MESS' 'Incrément de temps initial : ' DELTAT                        ;
'MESS' ' '                                                           ;
*
* Initialisation de la quantité dissoute
*
'SI' SOLUB                                                           ;
*    DIS et DIS1 sont des valeurs de dissolution valables pour
*    l'ensemble d'une maille, et intégré sur tout le pas de temps
     DIS  = (TRANSI . 'DISSOLUTION' . LAST1) * VOLU1 * DELTAT0       ;
     PRCI = TRANSI . 'PRECIPITE' . LAST1                             ;
     ANT0 = 'EXCO' 'H' ('MASQ' PRCI 'SUPERIEUR' 1.D-14) 'SCAL'       ;
'FINSI'                                                              ;
*
* Calcul de la fraction adsorbée initiale :
*
CC0  = 'NOMC' CHRG 'SCAL'                                            ;
'SI' RNONLIN                                                         ;
     'SI' RNONLINL                                                   ;
*         Langmuir :
          FF0 = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 ))          ;
     'FINSI'                                                         ;
     'SI' RNONLINF                                                   ;
*         Freundlich :
*         La bidouille utilisant les masques sert à obtenir le champ point
*         des signes de CC0.
          a = 'MASQUE' CC0 'SUPERIEUR' 0.                            ;
          b = 'MASQUE' CC0 'INFERIEUR' 0.                            ;
          FF0 = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b)           ;
          'DETRUIT' a                                                ;
          'DETRUIT' b                                                ;
     'FINSI'                                                         ;
'SINON'                                                              ;
*    Isotherme linéaire :
     FF0 = RETAR1M1 * CC0                                            ;
'FINSI'                                                              ;
FF  = FF0                                                            ;
 
*
* Calcul du coefficient de retard initial :
* Le facteur de retard s'exprime comme
*     R(C) = 1 + ( F'(Ct+dt + eps) - F'(Ct) ) / ( Ct+dt + eps - Ct) )
* C'est la méthode de la corde employée ici pour ses vertus
* stabilisantes. Le petit 'eps' (EpsCord) intervient pour évacuer les
* questions de conditions initiales ou de régime permanent.
'SI' RNONLIN                                                         ;
     CC1 = CC0 + EpsCord                                             ;
     'SI' RNONLINL                                                   ;
*         Langmuir
          FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 ))          ;
     'SINON'                                                         ;
*         Freundlich
          a = 'MASQUE' CC1 'SUPERIEUR' 0.                            ;
          b = 'MASQUE' CC1 'INFERIEUR' 0.                            ;
          FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b)           ;
          'DETRUIT' a                                                ;
          'DETRUIT' b                                                ;
     'FINSI'                                                         ;
     RETARC = 1.D0 + ( (FF1 - FF0) / EpsCord )                       ;
     'DETRUIT' CC1                                                   ;
     'DETRUIT' FF1                                                   ;
'SINON'                                                              ;
     RETARC = RETAR1                                                 ;
'FINSI'                                                              ;
RETAR0 = RETARC                                                      ;
'SI' SAUVRET                                                         ;
     TRANSI . 'RETARD' . 0 = RETARC                                  ;
'FINSI'                                                              ;
*
* Fonctions de volume :
PORSUF1  = POROS1 * VOLU1                                            ;
PORETSU1 = PORSUF1 * RETARC                                          ;
DECROI0  = 'EXCO' ( (-1.) * LAMBD0 * PORSUF1 ) 'CK' 'SCAL'           ;
*
* Indicateur stipulant si l'on doit recalculer la matrice hybride :
* Au départ, vrai car il faut la calculer au moins la première fois,
* puis faux par la suite jusqu-à notification contraire.
Recalhyb = VRAI                                                      ;
*
* Table utilisée par les opérateurs MATP, SMTP, HYBP et HDEB
* TAB.'SURF' est le coefficient devant dC/dt (défini dans chaque module).
* On y intègre la contribution de la décroissance implicite
* ainsi que celle du terme de pénalisation
TAB = 'TABLE'                                                        ;
TAB.'SOUSTYPE' = 'DARCY_TRANSITOIRE'                                 ;
TAB.'THETA'    = TETA                                                ;
TAB.'THETA_CONVECTION' = TETAC                                       ;
*
*
*=================================
'REPETE' BOUTPS (DCAL - ICAL + 1);    COMM 'Boucle sur le temps'     ;
*=================================
*
IPAS = ICAL + &BOUTPS - 1                                            ;
*
*-- Initialisation en-tête de boucle sur le temps
*
TPS          = 'EXTR' TPCAL IPAS                                     ;
MESS 'TEMPS T=' TPS                                                  ;
DELTAT       = TPS - PRECED                                          ;
EPSDT        = DELTAT + DELOLD / 2.D0 * 1.D-6                        ;
TAB.'PAS'    = DELTAT                                                ;
TAB.'TRACE'  = TRP                                                   ;
TAB.'CHARGE' = CHRG                                                  ;
Recalhyb     = Recalhyb 'OU' ( DELOLD 'NEG' DELTAT EPSDT )           ;
 
*
*- Calcul de la contribution des termes sources et de la convection
*
TERSCE = 'MANU' 'CHPO' HYTOT.CENTRE 1 'SOUR' 0.                      ;
*
* terme source de convection :
'SI' ((TETAC < 1.D-4) 'ET' ( 'EXISTE' TRANSI 'CONVECTION' ))         ;
    CONC1  = 'NOMC' 'SCAL' TRP                                       ;
    CONC2  = (-1.) * CONC1 * CONV1                                   ;
    TERSC1 = 'DIVU' MODHYB  CONC2 MCHYB                         ;
    TERSC1 = 'NOMC' 'SOUR' TERSC1                                    ;
    TERSCE = TERSC1 + TERSCE                                         ;
'FINSI'                                                              ;
*
* Terme source propre :
'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
    TERSC2 = 'EXCO' ('TIRE' TERSOU TPS) 'SOUR' 'SOUR'                ;
    TERSCE = TERSC2 + TERSCE                                         ;
'FINSI'                                                              ;
*
* Terme source de décroissance explicite :
* du soluté et de l'adsorbat
'SI' DECRO                                                           ;
     TERSC3 = ( ('NOMC' CHRG 'SCAL') + FF ) * DECROI0                ;
     TERSC3 = 'NOMC' 'SOUR' TERSC3                                   ;
     TERSCE = TERSC3 + TERSCE                                        ;
'FINSI'                                                              ;
*
*- Incorporation des CLs
*
SMTR0 = 'MANU' 'CHPO' HYTOT.FACE 1 'FLUX' 0. 'NATURE' 'DISCRET'      ;
'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                                   ;
     CHARIMPO = 'TIRE' CHAIMP TPS                                    ;
     CHARIMPO = 'CHAN' 'ATTRIBUT' CHARIMPO 'NATURE' 'DISCRET'        ;
     SMTR0    = SMTR0 'ET' CHARIMPO                                  ;
'FINSI'                                                              ;
*
'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' )                               ;
     FLUIMPO = 'TIRE' FLUIMP TPS                                     ;
     FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET'          ;
     SMTR0   = SMTR0 'ET' FLUIMPO                                    ;
'FINSI'                                                              ;
 
*
*
*|-------------------------------------------------------------------|
*| Boucle de Retard Non Linéaire                                     |
*|-------------------------------------------------------------------|
*
* Boucle de Picard pour la prise en compte
* d'isothermes d'adsorption non linéaires : isothermes de Langmuir et
* de Freundlich. Ce point nécessite d'effectuer une boucle sur
* le facteur de retard qui devient fonction de la concentration
* en solution.
*
*------------------------
'REPETE' BOURNL NPicard ;     COMM 'Boucle sur le coef. de retard'   ;
*------------------------
*
IRNL = &BOURNL                                                       ;
*
* Quatre grands blocs selon que l'on traite :
* -------------------------------------------
*  -1-  sans limite de solubilité ou avec dissolution arbitraire
*  -2-  avec limite de solubilité et dissolution d'ordre 1
*  -3-  avec limite suivant méthode de pénalisation,
*  -4-  avec limite suivant schéma prédicteur/correcteur itératif
 
 
*|----------------------------------------------------|
*| Pas de limite de solubilité ou Dissolution imposée |
*|----------------------------------------------------|
'SI' (('NON' SOLUB) 'OU' SOLUA)                                      ;
*
* On recalcule la matrice hybride si DeltaT a changé :
'SI' Recalhyb                                                        ;
*    Coefficients correspondant à la décroissance expl. et implicite :
     NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
     DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
*    Calcul du coefficient devant (Ct+dt - Ct)/dt
     COFDC = PORETSU1 * DEN0                                         ;
     TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM'                       ;
*    Matrice du problème
     HNDTR = 'MATP' MODHYB  MASHYB TAB                          ;
     'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                              ;
           HNDT1 = HNDTR 'ET' MATBLOC                                ;
     'SINON'                                                         ;
           HNDT1 = HNDTR                                             ;
     'FINSI'                                                         ;
*    On ajoute la rigidité associée à la convection implicite
*    Elle ne dépend ni du terme source, ni de la masse hybride,
*    ni de la valeur de la trace de concentration. On doit néanmoins
*    fournir des valeurs de ces grandeurs :
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCE CONV1   ;
          HNDT1 = HNDT1 'ET' RIGC                                    ;
     'FINSI'                                                         ;
     Recalhyb = FAUX                                                 ;
'FINSI'                                                              ;
*
'SI' SOLUA                                                           ;
*    Evaluation des termes sources de dissolution :
*    Dissolution arbitraire (ordre 0 par exemple), limitée au total précipité :
     DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1                        ;
     DISSOLP = 'TIRE' DISIMP TPS                                     ;
     TERSC4 = 0.5 * (DISSOLT + DISSOLP - ('ABS' (DISSOLT - DISSOLP)));
     TERSC4  = 'EXCO' 'H' TERSC4 'SOUR'                              ;
     TERSCT  = TERSCE + TERSC4                                       ;
'SINON'                                                              ;
     TERSCT  = TERSCE                                                ;
'FINSI'                                                              ;
 
* Rigidité et flux convectifs :
* SMTR0 contient la contribution en termes de flux des conditions aux
* limites
* Dans tous les cas, SMTR1 contient les termes sources et la convection,
* mais TERSCT ne contient la convection que si TETAC > 10-4
'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))         ;
*    TERSCT contient tous les T.S. sauf la convection,
     RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT CONV1        ;
'SINON'                                                              ;
*    TERSCT contient tous les T.S. convection comprise
     SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT                   ;
'FINSI'                                                              ;
SMTR  = SMTR0 + SMTR1                                                ;
 
*- Résolution en trace de concentration :
TP1   = 'RESO' HNDT1 SMTR                                            ;
TP2   = 'EXCO' TP1 'TH' 'TH'                                         ;
 
*- Terme source de convection :
* TERSCVT contient tous les T.S. y compris la convection
'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))         ;
*    On l'ajoute au terme source total maintenant seulement car on
*    a besoin de C(t+dt).
     CONC1   = ((1.D0 - TETAC) * TRP) + (TETAC * TP2)                ;
     CONC2   = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1               ;
     CONC3   = 'DIVU' MODHYB  CONC2 MCHYB                       ;
     TERSCV  = 'NOMC' 'SOUR' CONC3                                   ;
     TERSCVT = TERSCT + TERSCV                                       ;
'SINON'                                                              ;
*    Le terme source de convection est déjà dans TERSCT
     TERSCVT = TERSCT                                                ;
'FINSI'                                                              ;
 
* Calcul des nouveaux concentration, flux et précipité :
CHA2 = 'HYBP' MODHYB  MASHYB TP2 TAB TERSCVT                    ;
FLU2 = 'HDEB' MODHYB  MASHYB CHA2 TP2                           ;
'SI' SOLUA                                                           ;
     DIS1 = 'EXCO' 'SOUR' (TERSC4 * DELTAT) 'H'                      ;
     PRE2 = ((NUM0 * PRCI) - (DIS1 / VOLU1)) / DEN0                  ;
'FINSI'                                                              ;
 
* 'Fin module dissolution arbitraire'
'FINSI'                                                              ;
 
*|---------------------------|
*|   Dissolution d'ordre 1   |
*|---------------------------|
'SI' SOLUP;
 
* On recalcule la matrice hybride si DeltaT a changé :
'SI' Recalhyb                                                        ;
*    Coefficients correspondant à la décroissance expl. et implicite :
     NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
     DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
*    Terme implicite de dissolution :
     'SI' (GAMMA > 1.D-4)                                            ;
          DEN1 = ANT0 * GAMMA * CODIS * DELTAT                       ;
     'SINON'                                                         ;
          DEN1 = 0.                                                  ;
     'FINSI'                                                         ;
*    Calcul du coefficient devant (Ct+dt - Ct)/dt
     COFDC = PORETSU1 * (DEN0 + DEN1)                                ;
     TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM'                       ;
*    Matrice du problème
     HNDTR = 'MATP' MODHYB  MASHYB TAB                          ;
     'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                              ;
           HNDT1 = HNDTR 'ET' MATBLOC                                ;
     'SINON'                                                         ;
           HNDT1 = HNDTR                                             ;
     'FINSI'                                                         ;
*    On ajoute la rigidité associée à la convection implicite
*    Elle ne dépend ni du terme source, ni de la masse hybride,
*    ni de la valeur de la trace de concentration. On doit néanmoins
*    fournir des valeurs de ces grandeurs :
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCE CONV1   ;
          HNDT1 = HNDT1 'ET' RIGC                                    ;
     'FINSI'                                                         ;
     NBIT = NBIT + 1                                                 ;
     Recalhyb = FAUX                                                 ;
'FINSI'                                                              ;
*
* Evaluation des termes sources de dissolution (ici partie explicite) :
* terme de dissolution totale puis terme de dissolution d'ordre 1
* ramenés à une unité de temps :
DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1                             ;
DIFCHG  = LIMSOL - CHRG                                              ;
DISSOLP = ('EXCO' 'CK' PORETSU1 'SCAL') * CODIS * (LIMSOL - CHRG)    ;
 
*---------------
'REPETER' BOU8 ;     COMM 'Boucle dissolution ordre 1'               ;
*---------------
*    On prend celui des deux termes correspondant à la présence de
*    précipité de l'itération d'avant.
     TERSC4  = 'EXCO' 'H' (((1.D0 - ANT0) * DISSOLT)
                           + (ANT0 * DISSOLP)) 'SOUR'                ;
     TERSCT  = TERSCE + TERSC4                                       ;
 
*    Rigidité et flux convectifs :
*    SMTR0 contient la contribution en termes de flux des conditions aux
*    limites
*    Dans tous les cas, SMTR1 contient les termes sources et la
*    convection, mais TERSCT ne contient la convection que si TETAC>10-4
    'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))     ;
*         TERSCT contient tous les T.S. sauf la convection,
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT CONV1   ;
     'SINON'                                                         ;
*         TERSCT contient tous les T.S. convection comprise
          SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT              ;
     'FINSI'                                                         ;
     SMTR  = SMTR0 + SMTR1                                           ;
 
*-   Résolution en trace de concentration :
     TP1   = 'RESO' HNDT1 SMTR                                       ;
     TP2   = 'EXCO' TP1 'TH' 'TH'                                    ;
 
*-   Terme source de convection :
*    TERSCVT contient tous les T.S. y compris la convection
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
*         On l'ajoute au terme source total maintenant seulement car on
*         a besoin de C(t+dt).
          CONC1   = ((1.D0 - TETAC) * TRP) + (TETAC * TP2)           ;
          CONC2   = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1          ;
          CONC3   = 'DIVU' MODHYB  CONC2 MCHYB                  ;
          TERSCV  = 'NOMC' 'SOUR' CONC3                              ;
          TERSCVT = TERSCT + TERSCV                                  ;
     'SINON'                                                         ;
*         Le terme source de convection est déjà dans TERSCT
          TERSCVT = TERSCT                                           ;
     'FINSI'                                                         ;
 
*    Calcul des nouveaux concentration et précipité :
     CHA2 = 'HYBP' MODHYB  MASHYB TP2 TAB TERSCVT               ;
     CHA3 = (1.D0 - GAMMA) * CHRG + (GAMMA * CHA2)                   ;
     DIS3 = ('EXCO' 'CK' PORETSU1 'SCAL') / VOLU1
            * CODIS * (LIMSOL - CHA3) * DELTAT;
     PRE2 = ANT0 / DEN0 * ( (NUM0 * PRCI) - DIS3 )                   ;
*    On calcule le flux ici pour etre sur qu'il correspond aux données
*    en cours de la résolution. En effet, en cas de sortie sans
*    convergence, si on avait calculé le flux après les itérations,
*    MCHYB et TAB.SURF auraient tenu compte de la nouvelle distribution
*    ANT2 alors que les calculs ont été effectués avec la courante ANT0.
     FLU2 = 'HDEB' MODHYB  MASHYB CHA2 TP2                      ;
 
*    Critère de sortie et mise à jour de la nouvelle distribution :
     DIFOK = (NUM0 * PRCI) - DIS3                                    ;
     ANT2  = 'EXCO' 'H' ('MASQ' DIFOK 'SUPERIEUR' 0.) 'SCAL'         ;
     CRIT  = 'ABS' (ANT2 - ANT0)                                     ;
     ANT0  = ANT2                                                    ;
     'SI' ( ('MAXI' CRIT) < 1.D-14)                                  ;
          'QUITTER' BOU8                                             ;
     'SINON'                                                         ;
*         on reconstruit la matrice hybride
          'SI' (GAMMA > 1.D-4)                                       ;
               DEN1 = ANT0 * GAMMA * CODIS * DELTAT                  ;
          'SINON'                                                    ;
               DEN1 = 0.                                             ;
          'FINSI'                                                    ;
          COFDC = PORETSU1 * (DEN0 + DEN1)                           ;
          TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM'                 ;
          HNDTR = 'MATP' MODHYB  MASHYB TAB                     ;
          'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                         ;
                HNDT1 = HNDTR 'ET' MATBLOC                           ;
          'SINON'                                                    ;
                HNDT1 = HNDTR                                        ;
          'FINSI'                                                    ;
*         N'oublions pas la rigidité associée à la convection implicite
*         calculée avec le nouveau terme de pénalisation
          'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
               HNDT1 = HNDT1 'ET' RIGC                               ;
          'FINSI'                                                    ;
           NBIT = NBIT + 1                                           ;
     'FINSI'                                                         ;
     'SI' (&BOU8 EGA ITMAXI)                                         ;
          'MESS' 'Sortie sans convergence'                           ;
          'QUITTER' BOU8                                             ;
     'FINSI'                                                         ;
 
*-----------
'FIN' BOU8 ;
*-----------
 
* Calcul de la dissolution :
DIS1 = ((NUM0 * PRCI) - (PRE2 * DEN0)) * VOLU1                       ;
 
NB0 = &BOU8 - 1                                                      ;
'SI' (NB0 > 1);
    'MESS' '    ' NB0  ' itérations de solubilité'                   ;
'FINSI';
 
* 'Fin module Dissolution d ordre 1'
'FINSI'                                                              ;
 
*|-------------------------------------------------------------|
*|   Dissolution instantanée - Point fixe Présence précipité   |
*|-------------------------------------------------------------|
'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION'))                        ;
 
* On recalcule la matrice hybride si DeltaT a changé :
'SI' Recalhyb                                                        ;
*    Coefficients correspondant à la décroissance expl. et implicite :
     NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
     DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
*    Calcul du coefficient devant (Ct+dt - Ct)/dt
     PEN0  = ANT0 * PENAL * DELTAT                                   ;
     COFDC = PORETSU1 + (('EXCO' 'SCAL' PEN0 'CK') * VOLU1) * DEN0   ;
     TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM'                       ;
*    Matrice du problème
     HNDTR = 'MATP' MODHYB  MASHYB TAB                          ;
     'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                              ;
          HNDT1 = HNDTR 'ET' MATBLOC                                 ;
     'SINON'                                                         ;
          HNDT1 = HNDTR                                              ;
     'FINSI'                                                         ;
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         On ajoute la rigidité associée à la convection implicite
*         Elle ne dépend ni du terme source, ni de la masse hybride,
*         ni de la valeur de la trace de concentration. On doit néanmoins
*         fournir des valeurs de ces grandeurs :
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCE CONV1   ;
          HNDT1 = HNDT1 'ET' RIGC                                    ;
     'FINSI'                                                         ;
     NBIT = NBIT + 1                                                 ;
     Recalhyb = FAUX                                                 ;
'FINSI'                                                              ;
*
*- Résolution en trace de concentration et calcul des inconnues secondaires
*  Méthode du point fixe sur la présence de précipité
*  par pénalisation (dissolution instantanée)
 
* Décroissance et dissolution du précipité avant les itérations
* La filiation du précipité est en principe comprise dans
* celle du soluté (à la charge de l'utilisateur)
* On prendra au début la répartition de précipité du pas précédent
DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1                             ;
TERSC4  = 'EXCO' 'H' DISSOLT 'SOUR'                                  ;
TERSCE  = TERSCE + TERSC4                                            ;
TERSCT  = TERSCE                                                     ;
 
*--------------
'REPETE' BOU6 ;     'COMM Boucle limite de sol., Pénalisation'       ;
*--------------
*    Terme de pénalisation :
     TERSC5 = PEN0 * (DEN0 / DELTAT) * (LIMSOL - CHRG) * VOLU1       ;
     TERSC5 = 'EXCO' 'H' TERSC5 'SOUR'                               ;
     TERSCT = TERSCE + TERSC5                                        ;
 
*    Calcul du second membre :
*    SMTR0 contient la contribution en termes de flux des conditions aux
*    limites. Dans tous les cas, SMTR1 contient le second membre du
*    à la diffusion, les termes sources et le second membre du à
*    la convection, mais TERSCT ne contient la convection que si
*    TETAC est supérieur à 10-4
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         TERSCT contient tous les T.S. sauf la convection, on
*         regénére la rigidité au passage mais on ne s'en sert pas.
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT CONV1   ;
     'SINON'                                                         ;
*         TERSCT contient tous les T.S. convection comprise
          SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT              ;
     'FINSI'                                                         ;
     SMTR = SMTR0 + SMTR1                                            ;
 
*-   Résolution :
     TP1  = 'RESO' HNDT1 SMTR                                        ;
     TP2  = 'EXCO' TP1 'TH' 'TH'                                     ;
 
*-   Terme source de convection pour calculer la nouvelle concentration
*    TERSCVT contient tous les T.S. y compris la convection
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         On l'ajoute au terme source total maintenant seulement car on
*         a besoin de C(t+dt).
          CONC1   = ((1.D0 - TETAC) * TRP) + (TETAC * TP2)           ;
          CONC2   = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1          ;
          CONC3   = 'DIVU' MODHYB CONC2 MCHYB                  ;
          TERSCV  = 'NOMC' 'SOUR' CONC3                              ;
          TERSCVT = TERSCT + TERSCV                                  ;
     'SINON'                                                         ;
*         Le terme source de convection est déjà dans TERSCT
          TERSCVT = TERSCT                                           ;
     'FINSI'                                                         ;
     CHA2   = 'HYBP' MODHYB  MASHYB TP2 TAB TERSCVT             ;
     DIFCHG = CHA2 - LIMSOL                                          ;
     PRE2   = PEN0 * DIFCHG                                          ;
*    On calcule le flux ici pour etre sur qu'il correspond aux données
*    en cours de la résolution. En effet, en cas de sortie sans
*    convergence, si on avait calculé le flux après les itérations,
*    MCHYB et TAB.SURF auraient tenu compte de la nouvelle distribution
*    ANT2 alors que les calculs ont été effectués avec la courante ANT0.
     FLU2 = 'HDEB' MODHYB  MASHYB CHA2 TP2                      ;
 
*    Critère de sortie de la boucle d'itération :
     MASK1 = 'MASQ' DIFCHG 'SUPERIEUR' 0.                            ;
     MASK2 = 'MASQ' PRE2 'SUPERIEUR' 0.                              ;
     ANT2  = 'EXCO' 'H'
             ('MASQ' (MASK1 + MASK2) 'EGSUPE' .5) 'SCAL'             ;
     CRIT  = 'ABS' (ANT2 - ANT0)                                     ;
     ANT0  = ANT2                                                    ;
     'SI' ( ('MAXI' CRIT) < 1.D-14)                                  ;
          'QUITTER' BOU6                                             ;
     'SINON'                                                         ;
*         Calcul du coefficient devant (Ct+dt - Ct)/dt
          PEN0  = ANT0 * PENAL * DELTAT                              ;
          COFDC = PORETSU1 + (('EXCO' 'SCAL' PEN0 'CK') * VOLU1) * DEN0;
          TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM'                ;
*         Matrice du problème
          HNDTR = 'MATP' MODHYB  MASHYB TAB                     ;
          'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                         ;
                HNDT1 = HNDTR 'ET' MATBLOC                           ;
          'SINON'                                                    ;
                HNDT1 = HNDTR                                        ;
          'FINSI'                                                    ;
*         N'oublions pas la rigidité associée à la convection implicite
*         calculée avec le nouveau terme de pénalisation
          'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
               HNDT1 = HNDT1 'ET' RIGC                               ;
          'FINSI'                                                    ;
           NBIT = NBIT + 1                                           ;
     'FINSI'                                                         ;
     'SI' (&BOU6 EGA ITMAXI)                                         ;
          'MESS' 'Sortie sans convergence'                           ;
          'QUITTER' BOU6                                             ;
     'FINSI'                                                         ;
*-----------
'FIN' BOU6 ;
*-----------
 
* Calcul de la dissolution :
DIS1 = ((NUM0 * PRCI) - (PRE2 * DEN0)) * VOLU1                       ;
 
NB0 = &BOU6 - 1                                                      ;
'SI' (NB0 > 1);
    'MESS' '    ' NB0 ' itérations de solubilité'                    ;
'FINSI';
 
'MENA'                                                               ;
* Fin module PENALISATION
'FINSI'                                                              ;
 
*|------------------------------------------------------|
*|   Dissolution instantanée - Point fixe Dissolution   |
*|------------------------------------------------------|
'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION'))                      ;
 
* On recalcule la matrice hybride si DeltaT a changé :
'SI' Recalhyb                                                        ;
*    Coefficients correspondant à la décroissance expl. et implicite :
     NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
     DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
*    Calcul du coefficient devant (Ct+dt - Ct)/dt
     TAB.'SURF' = ('KCHA' MODHYB PORETSU1 'CHAM') * DEN0           ;
*    Matrice du problème
     HNDTR = 'MATP' MODHYB  MASHYB TAB                          ;
     'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                              ;
           HNDT1 = HNDTR 'ET' MATBLOC                                ;
     'SINON'                                                         ;
           HNDT1 = HNDTR                                             ;
     'FINSI'                                                         ;
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         On ajoute la rigidité associée à la convection implicite
*         Elle ne dépend ni du terme source, ni de la masse hybride,
*         ni de la valeur de la trace de concentration. On doit néanmoins
*         fournir des valeurs de ces grandeurs :
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCE CONV1   ;
          HNDT1 = HNDT1 'ET' RIGC                                    ;
     'FINSI'                                                         ;
     NBIT = NBIT + 1                                                 ;
     Recalhyb = FAUX                                                 ;
'FINSI'                                                              ;
*
*- Résolution en trace de concentration et calcul des inconnues secondaires
*  Méthode du point fixe sur la dissolution (dissolution instantanée)
 
* Décroissance et dissolution du précipité avant les itérations
* La filiation du précipité est en principe comprise dans
* celle du soluté (à la charge de l'utilisateur)
PRE1 = PRCI * (NUM0 / DEN0)                                          ;
* On prend pour commencer le T.S. de solubilité du pas précédent,
* limité au maximum précipité :
DIS0 = PRCI * VOLU1 * NUM0                                           ;
DIS1 = 0.5 * (DIS0 + DIS - ('ABS' (DIS0 - DIS)));
 
*--------------
'REPETE' BOU7 ;     'COMM Boucle limite de sol., Picard dissolution' ;
*--------------
     TERSC4 = DIS1 / DELTAT                                          ;
     TERSC4 = 'EXCO' TERSC4 'H' 'SOUR'                               ;
     TERSCT = TERSCE + TERSC4                                        ;
     PRE2   = PRE1 - (DIS1 / VOLU1 / DEN0)                           ;
 
*    Calcul du second membre :
*    SMTR0 contient la contribution en termes de flux des conditions aux
*    limites. Dans tous les cas, SMTR1 contient le second membre du
*    à la diffusion, les termes sources et le second membre du à
*    la convection, mais TERSCT ne contient la convection que si
*    TETAC est supérieur à 10-4
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         TERSCT contient tous les T.S. sauf la convection, on
*         regénére la rigidité au passage mais on ne s'en sert pas.
          RIGC SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT CONV1   ;
     'SINON'                                                         ;
*         TERSCT contient tous les T.S. convection comprise
          SMTR1 = 'SMTP' MODHYB  MASHYB TAB TERSCT              ;
     'FINSI'                                                         ;
     SMTR   = SMTR0 + SMTR1                                          ;
 
*-   Résolution
     TP1  = 'RESO' HNDT1 SMTR                                        ;
     TP2  = 'EXCO' TP1 'TH' 'TH'                                     ;
 
*-   Terme source de convection pour calculer la nouvelle concentration
*    TERSCVT contient tous les T.S. y compris la convection
     'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4))    ;
*         On l'ajoute au terme source total maintenant seulement car on
*         a besoin de C(t+dt).
          CONC1   = ((1.D0 - TETAC) * TRP) + (TETAC * TP2)           ;
          CONC2   = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1          ;
          CONC3   = 'DIVU' MODHYB  CONC2 MCHYB                  ;
          TERSCV  = 'NOMC' 'SOUR' CONC3                              ;
          TERSCVT = TERSCT + TERSCV                                  ;
     'SINON'                                                         ;
*         Le terme source de convection est déjà dans TERSCT
          TERSCVT = TERSCT                                           ;
     'FINSI'                                                         ;
     CHA2 = 'HYBP' MODHYB  MASHYB TP2 TAB TERSCVT               ;
     DIS3 = DEN0 * (LIMSOL - CHA2) * ('EXCO' PORETSU1 'CK' 'SCAL')   ;
*    On limite la dissolution au maximum précipité :
     DIS0 = PRE2 * VOLU1 * DEN0                                      ;
     DIS4 = 0.5 * (DIS0 + DIS3 - ('ABS' (DIS0 - DIS3)));
 
*    Critère de sortie de la boucle d'itération :
     CRIT = (ABS DIS4) -
            (DEN0 * EPS1 * ('EXCO' PORETSU1 'CK' 'SCAL') * LIMSOL)   ;
     'SI' ( (MAXI CRIT) < 0.)                                        ;
          'QUITTER' BOU7                                             ;
     'FINSI'                                                         ;
     'SI' (&BOU7 EGA ITMAXI)                                         ;
          'MESS' 'Sortie sans convergence'                           ;
          'QUITTER' BOU7                                             ;
     'FINSI'                                                         ;
 
     DIS1 = DIS1 + DIS4                                              ;
 
*-----------
'FIN' BOU7 ;
*-----------
 
NB0 = &BOU7 - 1                                                      ;
'SI' (NB0 > 1);
    'MESS' '    ' NB0 ' itérations de solubilité'                    ;
'FINSI';
NBIT = NBIT + NB0                                                    ;
 
* Calcul des nouveaux flux :
FLU2 = 'HDEB' MODHYB MASHYB CHA2 TP2                           ;
 
'MENA'                                                               ;
* Fin module POINT FIXE DISSOLUTION
'FINSI'                                                              ;
 
*|-------------------------------------|
*|  Fin boucle de Retard Non Linéaire  |
*|-------------------------------------|
*- Calcul du nouveau coefficient de retard et de l'adsorbât
'SI' RNONLIN                                                         ;
*    1.  F(C)
     CC0  = 'NOMC' CHRG 'SCAL'                                       ;
     'SI' RNONLINL                                                   ;
*         Langmuir
          FF = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 ))           ;
     'SINON'                                                         ;
*         Freundlich
*         La bidouille utilisant les masques sert à obtenir le champ point
*         des signes de CC0.
          a = 'MASQUE' CC0 'SUPERIEUR' 0.                            ;
          b = 'MASQUE' CC0 'INFERIEUR' 0.                            ;
          FF = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b)            ;
          'DETRUIT' a                                                ;
          'DETRUIT' b                                                ;
     'FINSI'                                                         ;
     FF0 = FF                                                        ;
*    2.  F(C + dC)
     CC1  = ('NOMC' CHA2 'SCAL') + EpsCord                           ;
     'SI' RNONLINL                                                   ;
*         Langmuir
          FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 ))          ;
     'SINON'                                                         ;
*         Freundlich
          a = 'MASQUE' CC1 'SUPERIEUR' 0.                            ;
          b = 'MASQUE' CC1 'INFERIEUR' 0.                            ;
          FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b)           ;
          'DETRUIT' a                                                ;
          'DETRUIT' b                                                ;
     'FINSI'                                                         ;
*    3.  R = 1 + dF/dC
     RETAR2 = 1.D0 + ( (FF1 - FF0) / (CC1 - CC0) )                   ;
     'DETRUIT' CC1                                                   ;
     'DETRUIT' FF1                                                   ;
*
*-   Critère de sortie de la boucle de retard non linéaire :
     CRIT = ( 'ABS' (RETAR2 - RETARC) ) / RETARC                     ;
     CRIT = 'MAXIMUM' CRIT                                           ;
     'SI' (CRIT < EPSRNL)                                            ;
          'QUITTER' BOURNL                                           ;
     'SINON'                                                         ;
*         Le coef. de retard ne change que si le critère est mauvais.
*         Réévaluation des variables affectées :
          RETARC   = RETAR2                                          ;
          PORETSU1 = PORSUF1 * RETARC                                ;
          Recalhyb = VRAI                                            ;
     'FINSI'                                                         ;
'SINON'                                                              ;
*    Retard linéaire F = (R - 1) C
     CC0 = 'NOMC' CHRG 'SCAL'                                        ;
     FF  = RETAR1M1 * CC0                                            ;
'FINSI'                                                              ;
 
* Si on atteint le nombre maximum d'itérations - en particulier si le
* retard est linéaire - on quitte la boucle :
'SI' (IRNL 'EGA' ITMAXRNL)                                           ;
     MESS
     '     sortie boucle de retard non linéaire sans convergence'    ;
     'QUITTER' BOURNL                                                ;
'FINSI'                                                              ;
 
*-------------
'FIN' BOURNL ;
*-------------
 
'SI' (IRNL > 1 )                                                     ;
     'MESS' '    ' IRNL ' itérations de retard non linéaire.'
            ' Critère = ' crit                                       ;
'FINSI'                                                              ;
NBITRNL = NBITRNL + (IRNL - 1)                                       ;
 
 
*|-----------------------------------|
*|   Fin des modules de résolution   |
*|-----------------------------------|
*
*- Calcul de la dissolution par unité de volume et par unité de temps :
*
'SI' SOLUB                                                           ;
     DIS2 = DIS1 / VOLU1 / DELTAT                                    ;
'FINSI'                                                              ;
*
*- Archivage des resultats
*
*  si le prochain temps à sauver tombe entre les deux
*  incréments de temps, on obtient les valeurs à sauver
*  par interpolation linéaire puis on les stocke.
*
'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' )                              ;
*    Sauvegarde de tous les temps intermédiaires, s'il y en a :
     'REPETER' BOU9                                                  ;
          TEMS = 'EXTR' TPSOR ISOR                                   ;
          'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '&lt;EG' TPS ) )       ;
               LAST1   = LAST1 + 1                                   ;
               DTEM    = ( TPS - TEMS ) / DELTAT                     ;
               UNMO    = 1. - DTEM                                   ;
               TPBIS   = 'COLI' TRP  DTEM TP2  UNMO                  ;
               CHBIS   = 'COLI' CHRG DTEM CHA2 UNMO                  ;
               FLUBIS  = 'COLI' FLU0 DTEM FLU2 UNMO                  ;
               'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ')               ;
                    RETBIS  = 'COLI' RETAR0 DTEM RETARC UNMO         ;
               'SINON'                                               ;
                    RETBIS  = (DTEM * RETAR0) + (UNMO * RETARC)      ;
               'FINSI'                                               ;
               'SI' SOLUB                                            ;
                    PRBIS   = 'COLI' PRCI DTEM PRE2 UNMO             ;
                    DISBIS  = 'COLI' (DIS / VOLU1 / DELTAT) DTEM
                                     DIS2 UNMO                       ;
               'FINSI'                                               ;
               TRANSI . 'TEMPS'            . LAST1 = TEMS            ;
               TRANSI . NOMTINC            . LAST1 = TPBIS           ;
               TRANSI . NOMINC             . LAST1 = CHBIS           ;
               TRANSI . 'FLUX'             . LAST1 = FLUBIS          ;
               'SI' SOLUB                                            ;
                    TRANSI . 'PRECIPITE'   . LAST1 = PRBIS           ;
                    TRANSI . 'DISSOLUTION' . LAST1 = DISBIS          ;
               'FINSI'                                               ;
               'SI' SAUVRET                                          ;
                    TRANSI . 'RETARD'      . LAST1 = RETBIS          ;
               'FINSI'                                               ;
               ISOR = ISOR + 1                                       ;
          'SINON'                                                    ;
               'QUITTER' BOU9                                        ;
          'FINSI'                                                    ;
          'SI' ( ISOR '>' DSOR )                                     ;
               ISOR = DSOR                                           ;
               'QUITTER' BOU9                                        ;
          'FINSI'                                                    ;
     'FIN' BOU9                                                      ;
'SINON'                                                              ;
*    sinon, on sauvegarde tous les temps :
     LAST1 = LAST1 + 1                                               ;
     TRANSI . 'TEMPS'            . LAST1 = TPS                       ;
     TRANSI . NOMTINC            . LAST1 = TP2                       ;
     TRANSI . NOMINC             . LAST1 = CHA2                      ;
     TRANSI . 'FLUX'             . LAST1 = FLU2                      ;
     'SI' SOLUB                                                      ;
          TRANSI . 'PRECIPITE'   . LAST1 = PRE2                      ;
          TRANSI . 'DISSOLUTION' . LAST1 = DIS2                      ;
     'FINSI'                                                         ;
     'SI' SAUVRET                                                    ;
          TRANSI . 'RETARD'      . LAST1 = RETARC                    ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
* Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait.
* Le dernier pas de temps ne peut pas faire l'objet d'interpolation
* linéaire si le temps sauvegardé ne tombe pas pile dessus.
'SI' (IPAS 'EGA' DCAL)                                               ;
     TPSDER = TRANSI. 'TEMPS' . LAST1                                ;
     'SI' ( TPSDER 'NEG' TPSFIN EPSDT )                              ;
          LAST1   = LAST1 + 1                                        ;
          TRANSI . 'TEMPS'            . LAST1 = TPS                  ;
          TRANSI . NOMTINC            . LAST1 = TP2                  ;
          TRANSI . NOMINC             . LAST1 = CHA2                 ;
          TRANSI . 'FLUX'             . LAST1 = FLU2                 ;
          'SI' SOLUB                                                 ;
               TRANSI . 'PRECIPITE'   . LAST1 = PRE2                 ;
               TRANSI . 'DISSOLUTION' . LAST1 = DIS2                 ;
          'FINSI'                                                    ;
          'SI' SAUVRET                                               ;
               TRANSI . 'RETARD'      . LAST1 = RETARC               ;
          'FINSI'                                                    ;
     'FINSI'                                                         ;
'FINSI'                                                              ;
*
*- Initialisations pour le pas suivant
*
PRECED = TPS                                                         ;
TRP    = TP2                                                         ;
CHRG   = CHA2                                                        ;
FLU0   = FLU2                                                        ;
RETAR0 = RETARC                                                      ;
DELOLD = DELTAT                                                      ;
'SI' SOLUB                                                           ;
     DIS  = DIS1                                                     ;
     PRCI = PRE2                                                     ;
'FINSI'                                                              ;
 
'MENAGE'                                                             ;
*
*=============
'FIN' BOUTPS ;    COMM 'Boucle sur le temps'                         ;
*=============
*
*- Indication du nombre de passages :
*
'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION'))                        ;
  'MESS' 'Nombre total de matrices calculées : ' NBIT                ;
'FINSI'                                                              ;
'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION'))                      ;
  'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
          NBIT ;
'FINSI'                                                              ;
'SI' SOLUP                                                           ;
  'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
          NBIT;
'FINSI'                                                              ;
'SI' RNONLIN                                                         ;
  'MESS' 'Nombre total d itérations supplémentaires de '
         'retard non linéaire : ' NBITRNL;
'FINSI'                                                              ;
 
'FINP'                                                               ;