Télécharger la plaquette Cast3M


fr en

Télécharger transgen.procedur

Retour à la liste

Numérotation des lignes : 

  1. * TRANSGEN  PROCEDUR  JC220346  12/09/12    21:15:09     7501           
  2. *TRANSGEN   PROCEDUR  GBM    04/01/03
  3. 'DEBP' TRANSGEN TRANSI*'TABLE';
  4.  
  5.  
  6. *   SUPPRIMER TOUS LES 'H', stocker composante de CHRG 'ET'
  7. *   REMPLACER LES 'H' par cela. A faire.
  8.  
  9.  
  10.  
  11. *---------------------------------------------------------------------
  12. * Résolution de l'équation de darcy pour un problème d'écoulement ou
  13. * de transport par une méthode d'éléments
  14. * finis mixtes hybrides ou VF. Les inconnues du problème sont
  15. *  - en EFMH, la concentration (H), la trace de
  16. *    concentration (TH) et le débit diffusif (FLUX).
  17. *  - en VF, la concentration (H)
  18. *---------------------------------------------------------------------
  19. *
  20. *------------------------------
  21. * Phrase d'appel (en GIBIANE) :
  22. *------------------------------
  23. *
  24. *  TRANSGEN TABLE ;
  25. *
  26. *----------------------------------
  27. * Opérandes (à mettre dans TABLE) :
  28. *----------------------------------
  29. *
  30. *  ___________________________________________________________________
  31. *  |                                                                 |
  32. *  | Indice      Contenu                                             |
  33. *  |                                                                 |
  34. *  -------------------------------------------------------------------
  35. *  |                                                                 |
  36. *  |------------------------------------------------                 |
  37. *  |Données physiques, géométriques et materielles :                 |
  38. *  |------------------------------------------------                 |
  39. *  |                                                                 |
  40. *  |'MODELE'      Objet modèle (MMODEL créé par MODE)                |
  41. *  |                                                                 |
  42. *  |'CARACTERISTIQUES' Données physiques et materielles :            |
  43. *  |              diffusivité effective (CHAMELEM créé par MATE)     |
  44. *  |                                                                 |
  45. *  |'POROSITE'    Valeur de la porosité (Type Champoint, Comp        |
  46. *  |              'SCAL', ou FLOTTANT) - Défaut 1.                   |
  47. *  |                                                                 |
  48. *  |'DECROISSANCE' Valeur du terme de décroissance (Type FLOTTANT)   |
  49. *  |              Tel que dC/dt = - Lambda * C     - Défaut 0.       |
  50. *  |                                                                 |
  51. *  |'COEF_RETARD' Coefficient de retard linéaire dans le cas simple, |
  52. *  |              ou  Pente à l'origine de la fonction F(C) dans le  |
  53. *  |                  cas d'isotherme non linéaire de Langmuir       |
  54. *  |              ou  Coefficient K de l'isotherme de Freundlich     |
  55. *  |              (Type CHPO Centre 'SCAL', ou FLOTTANT)             |
  56. *  |                                                                 |
  57. *  |'LANGMUIR'    Quantité maximale 'Fsat' adsorbée sur le solide    |
  58. *  |              rapportée à l'unité de volume du fluide et exprimée|
  59. *  |              dans la meme unité que le soluté.                  |
  60. *  |              (Type CHPO Centre 'SCAL', ou FLOTTANT).            |
  61. *  |              F = (R-1) C / [1 + ((R-1) C / Fsat)]               |
  62. *  |              Si cet indice et le suivant sont absents,          |
  63. *  |              l'équilibre d'adsorption est linéaire. Cet indice a|
  64. *  |              priorité sur l'indice FREUNDLICH.                  |
  65. *  |                                                                 |
  66. *  |'FREUNDLICH'  Exposant de la loi de Freundlich F = K (C ^ 1/n)   |
  67. *  |              (Type FLOTTANT).                                   |
  68. *  |              Dans ce cas (et si l'indice LANGMUIR n'existe pas),|
  69. *  |              l'indice 'COEF_RETARD' contient le coefficient     |
  70. *  |              K ramené à une unité de volume de fluide.          |
  71. *  |              - Non disponible pour l'instant -                  |
  72. *  |                                                                 |
  73. *  |'LIMITE_SOLUBILITE' Limite de solubilité (Type chpoin, Comp 'H') |
  74. *  |                                                                 |
  75. *  |'COEF_DISSOLUTION' Coef. de dissolution (Type CHPO Centre, Comp  |
  76. *  |              'SCAL'). Tel que dC/dt = Coef * (Csat - C) - Par   |
  77. *  |              défaut, la dissolution est instantanée             |
  78. *  |                                                                 |
  79. *  |'CONVECTION'  Débit intégré de la vitesse convective à travers   |
  80. *  |              chaque face des éléments (Type CHPO Face, comp.    |
  81. *  |              'FLUX')                                            |
  82. *  |                                                                 |
  83. *  |----------------------                                           |
  84. *  |Conditions initiales :                                           |
  85. *  |----------------------                                           |
  86. *  |                                                                 |
  87. *  |'TEMPS'       TABLE contenant à l'indice 0 la valeur du temps    |
  88. *  |              initial (FLOTTANT)                                 |
  89. *  |                                                                 |
  90. *  |                                                                 |
  91. *  |'CONCENTRATION' TABLE contenant à l'indice 0 la concentration    |
  92. *  |              (quantité d'élément par unité de volume d'eau)     |
  93. *  |              (Type CHPO Centre, Comp 'H')                       |
  94. *  |                                                                 |
  95. *  |'FLUX'        TABLE contenant à l'indice 0 le flux diffusif      |
  96. *  |              initial intégré sur chaque face (Type CHPO Face,   |
  97. *  |              comp. 'FLUX')                                      |
  98. *  |                                                                 |
  99. *  |'PRECIPITE'   TABLE contenant à l'indice 0 la quantité initiale  |
  100. *  |              de précipité par unité de volume de milieu solide  |
  101. *  |              (Type CHPO Centre, Comp 'H')                       |
  102. *  |                                                                 |
  103. *  |'DISSOLUTION' TABLE contenant à l'indice 0 la quantité initiale  |
  104. *  |              pour estimer la dissolution au premier pas de temps|
  105. *  |              (Type CHPO, Comp 'H'), voir plus loin.             |
  106. *  |                                                                 |
  107. *  |--------------------------------------                           |
  108. *  |Conditions aux limites / chargements :                           |
  109. *  |--------------------------------------                           |
  110. *  |                                                                 |
  111. *  |'BLOCAGE'     Contient les matrices de blocage (RIGIDITE)        |
  112. *  |                                                                 |
  113. *  |'TRACE_IMPOSE' Valeurs des traces imposées (charge ou concentra- |
  114. *  |              -tion) (CHARGEMENT 'TH' - Obligatoire si BLOCAGE)  |
  115. *  |                                                                 |
  116. *  |'FLUX_IMPOSE' Valeurs des flux imposés intégrés par face         |
  117. *  |              (Type CHARGEMENT de CHPO Face, comp. 'FLUX'-       |
  118. *  |              défaut 0.)                                         |
  119. *  |                                                                 |
  120. *  |'FLUXTOT_IMP' Valeurs des flux totaux imposés intégrés par face  |
  121. *  |              (Type CHARGEMENT de CHPO Face, comp. 'FLUX'-       |
  122. *  |              défaut 0.). Il s'agit du flux diffusif + convectif |
  123. *  |                                                                 |
  124. *  |'MIXTES'      Table : indice C contient les valeurs des flux     |
  125. *  |                               mixtes imposés intégrés par face  |
  126. *  |              (Type CHARGEMENT de CHPO Face, comp. 'FLUX'-       |
  127. *  |              défaut 0.) il est égal à A * flux diffusif +       |
  128. *  |                                       B * Concentration         |
  129. *  |              Indice A et B sont les réels A et B                |
  130. *  |                                                                 |
  131. *  |'SOURCE'      Valeurs du terme source par maille et par unité de |
  132. *  |              temps (ex : puits, filiation)                      |
  133. *  |              Les valeurs à l'indice i sont les valeurs entre    |
  134. *  |              les temps i-1 et i.                                |
  135. *  |              (CHARGEMENT de CHPO Centre, comp 'SOUR'- défaut 0.)|
  136. *  |                                                                 |
  137. *  |                                                                 |
  138. *  |'DISSOLUTION_IMPOSEE' Valeurs des dissolutions imposées par unité|
  139. *  |              de temps et par maille. (Type CHARGEMENT de CHPO,  |
  140. *  |              Comp 'H'). Les valeurs à l'indice i sont les       |
  141. *  |              valeurs moyennes de dissolution par unité de temps |
  142. *  |              entre les temps i-1 et i.                          |
  143. *  |              Priorité de la dissolution imposée sur les         |
  144. *  |              cinétiques.                                        |
  145. *  |                                                                 |
  146. *  |--------------------                                             |
  147. *  |Données numériques :                                             |
  148. *  |--------------------                                             |
  149. *  |                                                                 |
  150. *  |                                                                 |
  151. *  |'TEMPS_CALCULES' Valeur des temps calculés (LISTREEL)            |
  152. *  |              Contient obligatoirement le temps final.           |
  153. *  |                                                                 |
  154. *  |'TEMPS_SAUVES' Valeur des temps sauvegardés (LISTREEL - défaut : |
  155. *  |              on sauve tous les pas de temps)                    |
  156. *  |                                                                 |
  157. *  |                                                                 |
  158. *  |                                                                 |
  159. *  |'THETA_DIFF'  Coefficient de relaxation compris entre 0. et 1.   |
  160. *  |              (theta-méthode diffusion) ('FLOTTANT' - défaut 1.) |
  161. *  |                                                                 |
  162. *  |'THETA_CONV' Idem pour la convection                             |
  163. *  |              ('FLOTTANT', Défaut = THETA_DIFF)                  |
  164. *  |'THETA_DEC'   Idem mais pour la décroissance                     |
  165. *  |              ('FLOTTANT' - défaut 1/2)                          |
  166. *  |                                                                 |
  167. *  |'THETA_DISS'   Idem mais pour la dissolution                     |
  168. *  |              ('FLOTTANT' - défaut 1.)                           |
  169. *  |                                                                 |
  170. *  |'PENALISATION' Coefficient de pénalisation pour la prise en      |
  171. *  |              compte de la limite de solubilité. La présence de  |
  172. *  |              cet indice ou du suivant indique quel schéma a été |
  173. *  |              choisi.                                            |
  174. *  |              (Type 'FLOTTANT') - Valeur conseillée 1.D7         |
  175. *  |                                                                 |
  176. *  |'EPSI_LIM'    Précision relative d'arrêt pour le shéma limite de |
  177. *  |              solubilité prédicteur-correcteur itératif          |
  178. *  |              (Type FLOTTANT) - Valeur conseillée 5.D-3          |
  179. *  |                                                                 |
  180. *  |'ITMAX_LIM'   Nombre maxi d'itérations correspondant aux modules |
  181. *  |              de dissolution avant d'abandonner                  |
  182. *  |              (Type 'ENTIER') - Défaut 50                        |
  183. *  |                                                                 |
  184. *  |'EPSI_RET'    Précision relative d'arrêt pour la résolution      |
  185. *  |              itérative (Picard) de l'adsorption non linéaire    |
  186. *  |              (Type FLOTTANT) - Défaut 1.D-4                     |
  187. *  |                                                                 |
  188. *  |'EPSI_COR'    Petit saut de concentration pour calculer le coef. |
  189. *  |              de retard par la méthode de la corde lorsque le    |
  190. *  |              retard est non-linéaire.                           |
  191. *  |              (Type FLOTTANT) - Défaut 1.D-4                     |
  192. *  |                                                                 |
  193. *  |'ITMAX_RET'   Nombre maxi d'itérations correspondant au retard   |
  194. *  |              non linéaire avant d'abandonner.                   |
  195. *  |              (Type 'ENTIER') - Défaut 20                        |
  196. *  |_________________________________________________________________|
  197. *
  198. *
  199. *
  200. *---------------------------------
  201. * Résultats (stockés dans TABLE) :
  202. *---------------------------------
  203. *
  204. *  ___________________________________________________________________
  205. *  |                                                                 |
  206. *  | Indice      Contenu                                             |
  207. *  |                                                                 |
  208. *  -------------------------------------------------------------------
  209. *  |                                                                 |
  210. *  |                                                                 |
  211. *  |'TEMPS'       TABLE contenant les temps sauvegardés (FLOTTANT)   |
  212. *  |                                                                 |
  213. *  |'CONCENTRATION' TABLE contenant les concentrations               |
  214. *  |              (Type CHPO Centre, Comp 'H')                       |
  215. *  |                                                                 |
  216. *  |'FLUX'        TABLE contenant les débits diffusifs intégrés      |
  217. *  |              par face (Type CHPO Face, comp. 'FLUX')            |
  218. *  |                                                                 |
  219. *  |'PRECIPITE'   TABLE contenant la quantité de précipité par maille|
  220. *  |              (Type CHPO Centre, Comp 'H')                       |
  221. *  |                                                                 |
  222. *  |'DISSOLUTION' TABLE contenant la quantité de précipité dissoute  |
  223. *  |              entre deux pas de temps par unité de volume et par |
  224. *  |              unité de temps. La valeur stockée à l'indice i,    |
  225. *  |              est valable entre les temps i-1 et i               |
  226. *  |              (Type CHPO, Comp 'H').                             |
  227. *  |              ATTENTION, les valeurs de cette table résultat     |
  228. *  |              n'ont aucun sens lorsque les temps sauvegardés ne  |
  229. *  |              sont pas les memes que les temps calculés. Toute   |
  230. *  |              tentative d'exploitation donnera alors des         |
  231. *  |              résultats incohérents (erreurs de bilan)           |
  232. *  |                                                                 |
  233. *  |'RETARD'      Si cet indice a été préalablement défini comme une |
  234. *  |              TABLE, alors il contient les valeurs du coefficient|
  235. *  |              de retard (Type 'CHPO' centre, Comp 'SCAL'). Sinon,|
  236. *  |              les valeurs du coefficient de retard ne sont pas   |
  237. *  |              sauvegardées.                                      |
  238. *  |_________________________________________________________________|
  239. *
  240. *
  241. *  ___________________________________________________________________
  242. *  |                                                                 |
  243. *  | Les tables résultats sont indicées par des entiers variant de 0 |
  244. *  | à N .                                                           |
  245. *  | A l'indice 0 on stocke les valeurs initiales, aux indices       |
  246. *  | suivants les champs correspondant au temps de sortie TEMPS.I .  |
  247. *  | Les champs servant en cas de reprise sont ceux correpondant au  |
  248. *  | dernier indice.                                                 |
  249. *  |_________________________________________________________________|
  250. *
  251. *
  252. *                                                   TYPE DE RESOLUTION
  253.  
  254. NOMINC  = 'MOT' 'CONCENTRATION'                                      ;
  255. NOMTETA = 'MOT' 'THETA_DIFF'                                         ;
  256. NOMDDT  = 'MOT' 'POROSITE'                                           ;
  257.  
  258. *
  259. *--------------------------------------------------------------------*
  260. * RECUPERATION DES DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES    *
  261. *--------------------------------------------------------------------*
  262. *
  263. *                                                               MODELE
  264. 'SI' ( 'EXISTE' TRANSI 'MODELE' )                                    ;
  265.      MODDARCY = TRANSI . 'MODELE'                                    ;
  266. * METTRE LE PRECOND à 1 pour ne pas recalculer les préconditionnements
  267. * sur le maillage
  268.      ('DOMA' MODDARCY TABLE) . 'PRECONDI' = 1                        ;
  269. 'SINON'                                                              ;
  270.      'ERREUR' 'Il manque le modele.'                                 ;
  271.      'QUITTER' TRANSGEN                                              ;
  272. 'FINSI'                                                              ;
  273.  
  274. *                                                               VOLUME
  275. VOLU1 = 'DOMA' MODDARCY 'VOLUME'                                     ;
  276. *                                                          ORIENTATION
  277. MCHYB = 'DOMA' MODDARCY 'ORIENTAT'                                   ;
  278. *
  279. MCENT = 'DOMA' MODDARCY 'CENTRE'                                     ;
  280. *
  281. MFACE = 'DOMA' MODDARCY 'FACE'                                       ;
  282. *
  283. MAILL = 'DOMA' MODDARCY 'MAILLAGE'                                   ;
  284. *                                                            FILIATION
  285. 'SI' ('NON' ('EXISTE' TRANSI 'PERE'))                                ;
  286.    TRANSI . 'PERE' = FAUX                                            ;
  287. 'FINSI'                                                              ;
  288. *                                                     CARACTERISTIQUES
  289. 'SI' ( 'EXISTE' TRANSI 'CARACTERISTIQUES' )                          ;
  290.    'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHPOINT')       ;
  291.         MAT1 = TRANSI . 'CARACTERISTIQUES'                           ;
  292.    'FINSI'                                                           ;
  293. * Si la diffusion est un chargement, le traitement à lieu dans le
  294. * corps du programme
  295. 'SINON'                                                              ;
  296.      'ERREUR'
  297.        'Indice CARACTERISTIQUES absent de la table de données.'      ;
  298.      'QUITTER' TRANSGEN                                              ;
  299. 'FINSI'                                                              ;
  300. *                                           EMMAGASINEMENT OU POROSITE
  301. 'SI' ( 'EXISTE' TRANSI NOMDDT )                                      ;
  302.      POROS = TRANSI . NOMDDT                                         ;
  303. 'SINON'                                                              ;
  304.      POROS = 1.D0                                                    ;
  305. 'FINSI'                                                              ;
  306. 'SI' ('EGA' ('TYPE' POROS) 'FLOTTANT')                               ;
  307.      POROS = 'MANU' 'CHPO' MCENT 'SCAL' POROS                        ;
  308.      POROS = 'CHAN' 'ATTRIBUT' POROS  'NATURE' 'DISCRET'             ;
  309. * VERRUE
  310. *     POROS = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0 NATURE DISCRET          ;
  311. *
  312. 'FINSI'                                                              ;
  313. 'SI' ('NEG' ('TYPE' POROS) 'CHPOINT ')                               ;
  314.      'ERREUR' ('CHAINE' 'L indice ' NOMDDT ' doit etre de type '
  315.                'FLOTTANT ou CHAMPOINT')                              ;
  316. 'SINON'                                                              ;
  317.    POROS = 'NOMC' 'SCAL' POROS                                       ;
  318. 'FINSI'                                                              ;
  319.  
  320.  
  321. *'SI' ((('MINI' POROS) < 0.D0) 'OU' (('MAXI' POROS) > 1.D0))     ;
  322. *     'ERREUR' 'La porosité doit etre comprise entre 0 et 1'     ;
  323. *'FINSI'                                                         ;
  324.  
  325. *                                                               RETARD
  326. 'SI' ( 'EXISTE' TRANSI 'COEF_RETARD' )                               ;
  327.      RETAR1 = TRANSI . 'COEF_RETARD'                                 ;
  328.      'SI' ('EGA' ('TYPE' RETAR1) 'CHPOINT ')                         ;
  329.           RETAR1 = NOMC 'SCAL' RETAR1                                ;
  330.           MINRET = ('MINI' RETAR1)                                   ;
  331.      'SINON'                                                         ;
  332.           MINRET = RETAR1                                            ;
  333.           RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' MINRET                 ;
  334.           RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1  'NATURE' 'DISCRET'      ;
  335.      'FINSI'                                                         ;
  336.      'SI' (MINRET < 1.D0)                                            ;
  337.           'ERREUR' 'Le coefficient de retard doit etre supérieur à 1';
  338.      'FINSI'                                                         ;
  339. 'SINON'                                                              ;
  340. *     RETAR1 = kcht modhyb scal centre 1.D0                           ;
  341.      RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0                        ;
  342.      RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1  'NATURE' 'DISCRET'           ;
  343. 'FINSI'                                                              ;
  344. RETAR1M1 = RETAR1 - 1.                                               ;
  345. SAUVRET = 'EXISTE' TRANSI 'RETARD'                                   ;
  346. *                                                  RETARD NON LINEAIRE
  347. RNONLINL = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  348.            ('EXISTE' TRANSI 'LANGMUIR')                              ;
  349. RNONLINF = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  350.            ('EXISTE' TRANSI 'FREUNDLICH')  'ET' ('NON' RNONLINL)     ;
  351. RNONLIN  = RNONLINL 'OU' RNONLINF                                    ;
  352. * Cas Langmuir :
  353. 'SI' RNONLINL                                                        ;
  354.      FSAT = TRANSI . 'LANGMUIR'                                      ;
  355.      'SI' ('EGA' ('TYPE' FSAT) 'CHPOINT ')                           ;
  356.           'SI' (('MINI' FSAT) < 0.D0)                                ;
  357.             'ERREUR'
  358.               'La masse adsorbee a saturation doit etre positive'    ;
  359.             'QUITTER' TRANSGEN                                       ;
  360.           'FINSI'                                                    ;
  361.      'SINON'                                                         ;
  362.           'SI' (FSAT < 0.D0)                                         ;
  363.             'ERREUR'
  364.               'La masse adsorbee a saturation doit etre positive'    ;
  365.             'QUITTER' TRANSGEN                                       ;
  366.           'FINSI'                                                    ;
  367.      'FINSI'                                                         ;
  368.      RM1SURF  = RETAR1M1 / FSAT                                      ;
  369. 'FINSI'                                                              ;
  370. * Cas Freundlich :
  371. 'SI' RNONLINF                                                        ;
  372.      UNSURN = TRANSI . 'FREUNDLICH'                                  ;
  373.      'SI' ((UNSURN &lt;EG  0.D0) 'OU' (UNSURN > 1.D0))                  ;
  374.           'ERREUR'
  375.             'L exposant de Freundlich doit etre compris entre 0'
  376.             ' (exclu) et 1'                                          ;
  377.           'QUITTER' TRANSGEN                                         ;
  378.      'FINSI'                                                         ;
  379. 'FINSI'                                                              ;
  380. *                                                         DECROISSANCE
  381. DECRO = ( 'EXISTE' TRANSI 'DECROISSANCE' )                           ;
  382. 'SI' DECRO                                                           ;
  383.      LAMBD0 = TRANSI . 'DECROISSANCE'                                ;
  384.      'SI' (LAMBD0 < 0.D0)                                            ;
  385.          'ERREUR' 'Le coefficient de décroissance doit etre positif.';
  386.      'FINSI'                                                         ;
  387. 'SINON'                                                              ;
  388.      LAMBD0 = 0.                                                     ;
  389. 'FINSI'                                                              ;
  390. *                                                  DISSOLUTION_IMPOSEE
  391. SOLUA = ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE' )                     ;
  392. 'SI' SOLUA                                                           ;
  393.      DISIMP = TRANSI . 'DISSOLUTION_IMPOSEE'                         ;
  394. 'FINSI'                                                              ;
  395. *                                                    LIMITE_SOLUBILITE
  396. * (Priorité de la dissolution imposée sur les autres processus)
  397. SOLUL = (('EXISTE' TRANSI 'LIMITE_SOLUBILITE' ) 'ET' ('NON' SOLUA))  ;
  398. 'SI' SOLUL                                                           ;
  399.      'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT')    ;
  400.         LIMSOL = 'NOMC' 'H' TRANSI . 'LIMITE_SOLUBILITE'             ;
  401.      'FINSI'                                                         ;
  402.      'SI' (('MINI' LIMSOL) < 0.D0)                                   ;
  403.           mess 'minimum de limsol' (MINI LIMSOL);
  404.           'ERREUR' 'La limite de solubilité doit etre positive.'     ;
  405.      'FINSI'                                                         ;
  406. 'FINSI'                                                              ;
  407. *                                                     COEF_DISSOLUTION
  408. SOLUP = SOLUL 'ET' ('EXISTE' TRANSI 'COEF_DISSOLUTION' )             ;
  409. SOLUI = SOLUL 'ET' ('NON' SOLUP)                                     ;
  410. 'SI' SOLUP                                                           ;
  411.      codis = TRANSI . 'COEF_DISSOLUTION' ;
  412.      'SI' ('EGA' ('TYPE' CODIS) 'CHPOINT ')                          ;
  413.      CODIS = 'NOMC' 'SCAL' TRANSI . 'COEF_DISSOLUTION' NATU DIFFUS   ;
  414.           'SI' (('MINI' CODIS) < 0.D0)                               ;
  415.           'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  416.           'FINSI'                                                    ;
  417.      'SINON'                                                         ;
  418.           CODIS = TRANSI . 'COEF_DISSOLUTION'                        ;
  419.           CODIS = MANU CHPO MCENT SCAL CODIS                         ;
  420.           'SI' ((MINI CODIS)  < 0.D0)                                ;
  421.           'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  422.           'FINSI'                                                    ;
  423.      'FINSI'                                                         ;
  424. 'FINSI'                                                              ;
  425.  
  426. 'SI' (SOLUP)                                                         ;
  427. *  Là où la précipitation est nulle (codis =0) on choisit
  428. *  une limite de solubilité élevée à des fins d'optimisation
  429. *  algorithmique
  430.    'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT')      ;
  431.       dum = 'MASQUE' CODIS INFERIEUR 1.D-30                          ;
  432.       LIMSOL = ((1.D0 '-' dum) * LIMSOL) '+' ('NOMC' 'H' (1.D3*dum)) ;
  433.    'FINSI'                                                           ;
  434. 'FINSI'                                                              ;
  435. SOLUB = SOLUL 'OU' SOLUA                                             ;
  436. *
  437. * ------------------------------------------------------------------ *
  438. *  VERIFICATION DE LA BONNE STRUCTURE DES TABLES RESULTAT D'ORIGINE  *
  439. * ------------------------------------------------------------------ *
  440. *
  441. *                                                                TEMPS
  442. 'SI' ( 'NON' ('EXISTE' TRANSI 'TEMPS' ) )                            ;
  443.      'ERREUR'
  444.        'Indice TEMPS absent de la table de données.'                 ;
  445.      'QUITTER' TRANSGEN                                              ;
  446. 'FINSI'                                                              ;
  447. *                                                             INCONNUE
  448. 'SI' ( 'NON' ('EXISTE' TRANSI NOMINC ) )                             ;
  449.      'ERREUR' ('CHAINE'
  450.        'Indice ' NOMINC ' absent de la table de données.')           ;
  451.      'QUITTER' TRANSGEN                                              ;
  452. 'FINSI'                                                              ;
  453. *                                                             FLUXDIFF
  454. 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXDIFF' ) )                         ;
  455.      'ERREUR'
  456.        'Indice FLUX diffusif absent de la table de données.'         ;
  457.      'QUITTER' TRANSGEN                                              ;
  458. 'FINSI'                                                              ;
  459. *                                                             FLUXCONV
  460. 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXCONV' ) )                         ;
  461.      'ERREUR'
  462.        'Indice FLUX convectif absent de la table de données.'        ;
  463.      'QUITTER' TRANSGEN                                              ;
  464. 'FINSI'                                                              ;
  465. *                                               DISSOLUTION, PRECIPITE
  466. 'SI' SOLUB                                                           ;
  467. *    GBM mettre a jour dissolution .
  468. *     'SI' ( 'NON' ('EXISTE' TRANSI 'DISSOLUTION' ) )                 ;
  469. *          'ERREUR'
  470. *            'Indice DISSOLUTION absent de la table de données.'      ;
  471. *          'QUITTER' TRANSGEN                                         ;
  472. *     'FINSI'                                                         ;
  473.      'SI' ( 'NON' ('EXISTE' TRANSI 'PRECIPITE' ) )                   ;
  474.           'ERREUR'
  475.             'Indice PRECIPITE absent de la table de données.'        ;
  476.           'QUITTER' TRANSGEN                                         ;
  477.      'FINSI'                                                         ;
  478. 'FINSI'                                                              ;
  479. *                                            TEST DES TAILLES DE TABLE
  480. IND1 = 'INDEX' ( TRANSI . 'TEMPS' )                                  ;
  481. IND3 = 'INDEX' ( TRANSI . NOMINC )                                   ;
  482. IND4 = 'INDEX' ( TRANSI . 'FLUXDIFF' )                               ;
  483. LIN1 = 'DIME' IND1                                                   ;
  484. LIN3 = 'DIME' IND3                                                   ;
  485. LIN4 = 'DIME' IND4                                                   ;
  486. 'SI' ( LIN1 'NEG' LIN3 )                                             ;
  487.      'ERREUR' ('CHAINE'
  488.           'Longueur des tables TEMPS et ' NOMINC ' différente.')     ;
  489.      'QUITTER' TRANSGEN                                              ;
  490. 'FINSI'                                                              ;
  491. 'SI' ( LIN1 'NEG' LIN4 )                                             ;
  492.      'ERREUR'
  493.        'Longueur des tables TEMPS et FLUXDIFF différente.'           ;
  494.      'QUITTER' TRANSGEN                                              ;
  495. 'FINSI'                                                              ;
  496. 'SI' SOLUB                                                           ;
  497. *    GBM mettre a jour dissolution .
  498. *     IND5 = 'INDEX' ( TRANSI . 'DISSOLUTION' )                       ;
  499.      IND6 = 'INDEX' ( TRANSI . 'PRECIPITE' )                         ;
  500. *     LIN5 = 'DIME' IND5                                              ;
  501.      LIN6 = 'DIME' IND6                                              ;
  502. *     'SI' ( LIN1 'NEG' LIN5 )                                        ;
  503. *          'ERREUR'
  504. *            'Longueur des tables TEMPS et DISSOLUTION différente.'   ;
  505. *          'QUITTER' TRANSGEN                                         ;
  506. *     'FINSI'                                                         ;
  507.      'SI' ( LIN1 'NEG' LIN6 )                                        ;
  508.           'ERREUR'
  509.              'Longueur des tables TEMPS et PRECIPITE différente.'    ;
  510.           'QUITTER' TRANSGEN                                         ;
  511.      'FINSI'                                                         ;
  512. 'FINSI'                                                              ;
  513. *                                            TEST DES INDICES DE TABLE
  514. IPO1 = 0                                                             ;
  515. 'REPETER' BOU1 LIN1                                                  ;
  516.      IPO1  = IPO1 + 1                                                ;
  517.      LAST1 = IND1 . IPO1                                             ;
  518.      LAST3 = IND3 . IPO1                                             ;
  519.      LAST4 = IND4 . IPO1                                             ;
  520.      'SI' SOLUB                                                      ;
  521. *    GBM mettre a jour dissolution .
  522. *          LAST5 = IND5 . IPO1                                        ;
  523.           LAST6 = IND6 . IPO1                                        ;
  524.      'FINSI'                                                         ;
  525.      'SI' ( 'NEG' LAST1 LAST3 )                                      ;
  526.           'ERREUR'
  527.             'Indices des tables TEMPS et CONCENTRATION incohérents.' ;
  528.           'QUITTER' TRANSGEN                                         ;
  529.      'FINSI'                                                         ;
  530.      'SI' ( 'NEG' LAST1 LAST4 )                                      ;
  531.           'ERREUR'
  532.             'Indices des tables TEMPS et FLUX incohérents.'          ;
  533.           'QUITTER' TRANSGEN                                         ;
  534.      'FINSI'                                                         ;
  535.      'SI' SOLUB                                                      ;
  536. *    GBM mettre a jour dissolution .
  537. *        'SI' ( 'NEG' LAST1 LAST5 )                                   ;
  538. *             'ERREUR'
  539. *               'Indices des tables TEMPS et DISSOLUTION incohérents.';
  540. *             'QUITTER' TRANSGEN                                      ;
  541. *        'FINSI'                                                      ;
  542.         'SI' ( 'NEG' LAST1 LAST6 )                                   ;
  543.              'ERREUR'
  544.                'Indices des tables TEMPS et PRECIPITE incohérents.'  ;
  545.              'QUITTER' TRANSGEN                                      ;
  546.         'FINSI'                                                      ;
  547.      'FINSI'                                                         ;
  548. 'FIN' BOU1                                                           ;
  549. *
  550. * ------------------------------------------------------------------ *
  551. *  RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS SAUVE     *
  552. * ------------------------------------------------------------------ *
  553. *
  554. TPSINI = TRANSI . 'TEMPS' . LAST1                                    ;
  555.  
  556. *Momentanée pour la trace
  557.  
  558. * On extrait les noms de composante des espèces
  559. NOMESP = 'EXTRAIRE' ('EXTRAIRE' TRANSI . NOMINC . LAST1 'COMP') 1    ;
  560.  
  561. * nomespl est un nom d'espece locale à la procédure car
  562. * les EFMH veulent voir une composante 'H'
  563. NOMESPL = 'MOT' 'H' ;
  564.  
  565. CHRG   = 'NOMC' 'H' TRANSI . NOMINC . LAST1                          ;
  566. *VERRUE
  567. FLU0   = 'NOMC' 'H' TRANSI . 'FLUXDIFF' . LAST1                   ;
  568. FLUCO0 = 'NOMC' 'H' TRANSI . 'FLUXCONV' . LAST1                   ;
  569.  
  570.  
  571. *
  572. *--------------------------------------------------------------------*
  573. * RECUPERATION DES THETA SCHEMAS DIFFUSION-CONVECTION                *
  574. *--------------------------------------------------------------------*
  575. *                                                      THETA DIFFUSION
  576.   'SI' ( 'EXISTE' TRANSI 'THETA_DIFF' )                              ;
  577.        TETA = TRANSI . 'THETA_DIFF'                                  ;
  578.   'SINON'                                                            ;
  579.        TETA = 1.D0                                                   ;
  580.   'FINSI'                                                            ;
  581. *                                                     THETA CONVECTION
  582. 'SI' ( 'EXISTE' TRANSI 'THETA_CONVECTION' )                          ;
  583.      TETAC = TRANSI . 'THETA_CONVECTION'                             ;
  584. 'SINON'                                                              ;
  585.      TETAC = TETA                                                    ;
  586. 'FINSI'                                                              ;
  587. *
  588. *--------------------------------------------------------------------*
  589. *  RECUPERATION DES CONDITIONS AUX LIMITES ET DES CHARGEMENTS        *
  590. *--------------------------------------------------------------------*
  591. *
  592. *                                  CL de type Dirichlet :      BLOCAGE
  593. *                                                         TRACE_IMPOSE
  594.  
  595. * On débranche les tests sur les longueurs de flux non nuls
  596. * en CL pour les cas ou ils ne sont pas intégrés en temps
  597. startflu = 1.D30;
  598. startflt = 1.D30;
  599. startmix = 1.D30;
  600. 'SI' ( 'EXISTE' TRANSI 'TRACE_IMPOSE')                               ;
  601.    CHDIRI = TRANSI . 'TRACE_IMPOSE'                                  ;
  602. 'FINSI'                                                              ;
  603. *                                     CL de type Neumann : FLUX_IMPOSE
  604. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' )                               ;
  605.      FLUIMP = TRANSI . 'FLUX_IMPOSE'                                 ;
  606.      'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  607.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  608.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  609. *       On écrit les flux sous forme intégrale sauf en explicite et
  610. *       kranck-Nickholson pour EFMH
  611.         FLUIMP startflu = CHAMINT FLUIMP                             ;
  612.      'FINSI'                                                         ;
  613. 'FINSI'                                                              ;
  614. *                                   CL de type flux total: FLUXTOT_IMP
  615. 'SI' ( 'EXISTE' TRANSI 'FLUXTOT_IMP' )                               ;
  616.      FLTOTIMP = TRANSI . 'FLUXTOT_IMP'                               ;
  617.      'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  618.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  619.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  620. *       On écrit les flux sous forme intégrale sauf en explicite et
  621. *       kranck-Nickholson pour EFMH
  622.         FLTOTIMP startflt = CHAMINT FLTOTIMP                         ;
  623.      'FINSI'                                                         ;
  624. 'FINSI'                                                              ;
  625. *                                            CL de type mixte : MIXTES
  626. 'SI' ( 'EXISTE' TRANSI 'MIXTES' )                                    ;
  627.      FLUMIX = TRANSI . 'MIXTES'                                      ;
  628.      'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  629.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  630.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  631. *       On écrit les flux sous forme intégrale sauf en explicite et
  632. *       kranck-Nickholson pour EFMH
  633.         FLUMIX startmix = CHAMINT FLUMIX                             ;
  634.      'FINSI'                                                         ;
  635. 'FINSI'                                                              ;
  636.  
  637. 'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
  638.      TERSOU = TRANSI . 'SOURCE'                                      ;
  639. * on crée une evolution intégrée temporellement pour assurer la
  640. * conservation lors de liste de pas de temps de chargement differente
  641. * des temps de discrétisation
  642.      TERSOU startsou = CHAMINT TERSOU                                ;
  643. 'FINSI'                                                              ;
  644.  
  645.  
  646. *
  647. *--------------------------------------------------------------------*
  648. * RECUPERATION DES données NUMERIQUES                                *
  649. *--------------------------------------------------------------------*
  650. *
  651. *                                             PARAMETRES DE SOLUBILITE
  652. 'SI' (SOLUI 'OU' SOLUP)                                              ;
  653.      'SI' ( 'EXISTE' TRANSI 'ITMAX_LIM')                             ;
  654.           ITMAXI = TRANSI . 'ITMAX_LIM'                              ;
  655.      'SINON'                                                         ;
  656.           ITMAXI = 50                                                ;
  657.      'FINSI'                                                         ;
  658. 'FINSI'                                                              ;
  659. 'SI' SOLUI                                                           ;
  660.      'SI' ( 'EXISTE' TRANSI 'PENALISATION')                          ;
  661.           PENAL = TRANSI . 'PENALISATION'                            ;
  662.           MET0  = 'PENALISATION'                                     ;
  663.      'SINON'                                                         ;
  664.           'SI' ( 'EXISTE' TRANSI 'EPSI_LIM')                         ;
  665.                EPS1 = TRANSI . 'EPSI_LIM'                            ;
  666.                MET0 = 'PF DISSOLUTION'                               ;
  667.           'SINON'                                                    ;
  668.                'MESS' 'Il manque le paramètre '                      ;
  669.                'ERREUR' 'du schéma numérique de limite de solubilité';
  670.           'FINSI'                                                    ;
  671.      'FINSI'                                                         ;
  672. 'FINSI'                                                              ;
  673. *                                    PARAMETRES DE RETARD NON LINEAIRE
  674. 'SI' RNONLIN                                                         ;
  675.      'SI' ( 'EXISTE' TRANSI 'EPSI_RET' )                             ;
  676.           EPSRNL = TRANSI . 'EPSI_RET'                               ;
  677.      'SINON'                                                         ;
  678.           EPSRNL = 1.D-4                                             ;
  679.      'FINSI'                                                         ;
  680.      'SI' ( 'EXISTE' TRANSI 'ITMAX_RET' )                            ;
  681.           ITMAXRNL = TRANSI . 'ITMAX_RET'                            ;
  682.      'SINON'                                                         ;
  683.           ITMAXRNL = 20                                              ;
  684.      'FINSI'                                                         ;
  685. *    Petit saut de concentration minimal pour le calcul de dérivée
  686.      'SI' ( 'EXISTE' TRANSI 'EPSI_COR' )                             ;
  687.           EPSCORD = TRANSI 'EPSI_COR'                                ;
  688.      'SINON'                                                         ;
  689.           'SI' RNONLINL                                              ;
  690.                EPSCORD = ('MANU' 'CHPO' MCENT 1 'SCAL' 1.D-4)
  691.                        / RETAR1 * FSAT                               ;
  692.           'FINSI'                                                    ;
  693.           'SI' RNONLINF                                              ;
  694.                EPSCORD = 1.D-4                                       ;
  695.           'FINSI'                                                    ;
  696.      'FINSI'                                                         ;
  697. *    Pour les isothermes linéaires la boucle de prise en compte de
  698. *    retard non linéaire n'est parcourue qu'une fois :
  699.      NPicard = ITMAXRNL                                              ;
  700. 'SINON'                                                              ;
  701.      NPicard = 1                                                     ;
  702. 'FINSI'                                                              ;
  703.  
  704. *                                                   THETA_DECROISSANCE
  705. 'SI' ( 'EXISTE' TRANSI 'THETA_DEC' )                                 ;
  706.      BETA = TRANSI . 'THETA_DEC'                                     ;
  707. 'SINON'                                                              ;
  708.      BETA = 0.5                                                      ;
  709. 'FINSI'                                                              ;
  710. *                                                    THETA_DISSOLUTION
  711. 'SI' ( 'EXISTE' TRANSI 'THETA_DIS' )                                 ;
  712.      GAMMA = TRANSI . 'THETA_DIS'                                    ;
  713. 'SINON'                                                              ;
  714.      GAMMA = 1.                                                      ;
  715. 'FINSI'                                                              ;
  716. *                                                       TEMPS_CALCULES
  717. 'SI' ( 'EXISTE' TRANSI 'TEMPS_CALCULES' )                            ;
  718.      TPCAL  = 'ORDO' (TRANSI . 'TEMPS_CALCULES')                     ;
  719.      TRANSI . 'TEMPS_CALCULES' = TPCAL                               ;
  720.      DCAL   = 'DIME' TPCAL                                           ;
  721.      TPSFIN = 'EXTR' TPCAL DCAL                                      ;
  722.      TPS0   = 'EXTR' TPCAL 1                                         ;
  723.      'SI' ('NEG' DCAL 1)                                             ;
  724.           TPSANT = 'EXTR' TPCAL (DCAL - 1)                           ;
  725.           TPS1   = 'EXTR' TPCAL 2                                    ;
  726.           DT10   = TPS1 - TPS0                                       ;
  727.           DT21   = TPSFIN - TPSANT                                   ;
  728.           PROPS  = 'PROG' DT10 DT21                                  ;
  729.           DTEPS  = 'MINIMUM' PROPS                                   ;
  730.      'SINON'                                                         ;
  731.           DTEPS  = TPSFIN - TPS0                                     ;
  732.      'FINSI'                                                         ;
  733. *
  734.      EPS0  = DTEPS * 1.D-6                                           ;
  735. *    ICAL est le premier indice des temps à calculer supérieur à
  736. *    TpsIni + eps
  737.      IOK1 = 0                                                        ;
  738.      'REPETER' BOU2 DCAL                                             ;
  739.           ICAL = &BOU2                                               ;
  740.           TEMS = 'EXTR' TPCAL ICAL                                   ;
  741.           'SI' ( TPSINI '<' (TEMS - EPS0) )                          ;
  742.               IOK1 = 1                                               ;
  743.               'QUITTER' BOU2                                         ;
  744.           'FINSI'                                                    ;
  745.      'FIN' BOU2                                                      ;
  746.      'SI' ( IOK1 'EGA' 0)                                            ;
  747.           'MESS' 'Listreel des temps de calcul inférieur à tini.'    ;
  748.           'MESS' 'On ne fait donc rien.                         '    ;
  749.           'QUITTER' TRANSGEN                                         ;
  750.      'FINSI'                                                         ;
  751.      DELTAT0 = TEMS - TPSINI                                         ;
  752. 'SINON'                                                              ;
  753.      'ERREUR'
  754.        'Indice TEMPS_CALCULES absent de la table de données.'        ;
  755.      'QUITTER' TRANSGEN                                              ;
  756. 'FINSI'                                                              ;
  757. *                                                         TEMPS_SAUVES
  758. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' )                              ;
  759.      TPSOR = 'ORDO' (TRANSI . 'TEMPS_SAUVES')                        ;
  760.      TRANSI . 'TEMPS_SAUVES' = TPSOR                                 ;
  761.      DSOR   = 'DIME' TPSOR                                           ;
  762.      ISOR = 0                                                        ;
  763.      IOK1 = 0                                                        ;
  764.      'REPETER' BOU3 DSOR                                             ;
  765.           ISOR = ISOR + 1                                            ;
  766.           TEMS = 'EXTR' TPSOR ISOR                                   ;
  767.           'SI' ( ( TPSINI '<' (TEMS - EPS0) ) 'ET'
  768.                  ( TEMS   '&lt;EG' (TPSFIN +EPS0)) )                    ;
  769.                IOK1 = 1                                              ;
  770.                'QUITTER' BOU3                                        ;
  771.           'FINSI'                                                    ;
  772.      'FIN' BOU3                                                      ;
  773.      'SI' ( IOK1 'EGA' 0)                                            ;
  774.           'ERREUR' 'Listreel des temps de sauvegarde hors tmin,tmax.';
  775.           'QUITTER' TRANSGEN                                         ;
  776.      'FINSI'                                                         ;
  777. 'FINSI'                                                              ;
  778. *
  779. *--------------------------------------------------------------------*
  780. * RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE                       *
  781. *--------------------------------------------------------------------*
  782. *
  783. *
  784. 'SI' (EXISTE TRANSI 'AFFICH') ;
  785. affich = TRANSI . 'AFFICH' ;
  786. SINON;
  787. affich = VRAI;
  788. 'FINSI'                       ;
  789.  
  790. 'SI' (affich) ;
  791. 'MESS' ' '                                                           ;
  792. 'MESS' 'MODELISATION Darcy en transitoire.'
  793.               (TRANSI . 'TYPDISCRETISATION')                         ;
  794. 'MESS' '--------------------------------------------'                ;
  795. 'MESS' ' '                                                           ;
  796. 'MESS' 'Donnees présentes en entrée : '                              ;
  797. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' )                                   ;
  798.      'MESS' 'Calcul avec conditions aux limites sur Th (Dirichlet)'  ;
  799. 'FINSI'                                                              ;
  800. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' )                               ;
  801.      'MESS' 'Calcul avec conditions aux limites de flux (Neumann)'   ;
  802. 'FINSI'                                                              ;
  803. 'SI' ( 'EXISTE' TRANSI 'FLUXTOT_IMP' )                               ;
  804.      'MESS' 'Calcul avec conditions aux limites de flux total)'      ;
  805. 'FINSI'                                                              ;
  806. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
  807.      'MESS' 'Ce problème comporte un terme convectif'                ;
  808. 'FINSI'                                                              ;
  809. 'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
  810.      'MESS' 'Ce problème comporte un terme source'                   ;
  811. 'FINSI'                                                              ;
  812. 'SI' DECRO                                                           ;
  813.      'MESS' 'Ce problème comporte un terme de décroissance'          ;
  814. 'FINSI'                                                              ;
  815. 'SI' RNONLIN                                                         ;
  816.      'MESS' 'Ce problème comporte la prise en compte d isothermes '
  817.             'd adsorption non-linéaires'                             ;
  818.      'SI' RNONLINL                                                   ;
  819.           'MESS' '  de type Langmuir (méthode de Picard)'            ;
  820.      'FINSI'                                                         ;
  821.      'SI' RNONLINF                                                   ;
  822.           'MESS' '  de type Freundlich (méthode de Picard)'          ;
  823.      'FINSI'                                                         ;
  824. 'FINSI'                                                              ;
  825. 'SI' SOLUL                                                           ;
  826.      'MESS' 'Ce problème comporte une limite de solubilité'          ;
  827.      'SI' SOLUP                                                      ;
  828.           'MESS' '  avec dissolution progressive d ordre 1'          ;
  829.      'FINSI'                                                         ;
  830.      'SI' (EGA MET0 'PENALISATION')                                  ;
  831.           'MESS' '  avec coefficient de pénalisation = ' PENAL       ;
  832.      'FINSI'                                                         ;
  833.      'SI' (EGA MET0 'PF DISSOLUTION')                                ;
  834.           'MESS' '  avec critère de convergence = ' EPS1             ;
  835.      'FINSI'                                                         ;
  836. 'FINSI'                                                              ;
  837. 'SI' SOLUA                                                           ;
  838.      'MESS' 'Ce probleme comporte une dissolution imposée'           ;
  839. 'FINSI'                                                              ;
  840. 'MESS' ' '                                                           ;
  841. 'MESS' 'Valeur des paramètres des schémas numériques :'              ;
  842. 'MESS' '  (0:Schéma explicite, 0.5:Crank-Nicholson, 1:Implicite)'    ;
  843. 'MESS' '  diffusion    : ' TETA                                      ;
  844. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
  845.      'MESS' '  convection   : ' TETAC                                ;
  846. 'FINSI'                                                              ;
  847. 'SI' DECRO                                                           ;
  848.      'MESS' '  décroissance : ' BETA                                 ;
  849. 'FINSI'                                                              ;
  850. 'SI' SOLUP                                                           ;
  851.      'MESS' '  dissolution d ordre 1 : ' GAMMA                       ;
  852. 'FINSI'                                                              ;
  853. 'MESS' ' '                                                           ;
  854. 'MESS' 'Valeur du temps initial    : ' TPSINI                        ;
  855. 'MESS' 'Valeur du temps final      : ' TPSFIN                        ;
  856. 'FINSI' ;
  857. *
  858. *
  859. *--------------------------------------------------------------------*
  860. * BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS               *
  861. *--------------------------------------------------------------------*
  862. *
  863. *==============
  864. * Préparation préliminaire de la boucle sur les pas de temps :
  865. *==============
  866. *
  867. NBIT    = 0                                                          ;
  868. NBITRNL = 0                                                          ;
  869. PRECED  = TPSINI                                                     ;
  870. IPAS    = ICAL                                                       ;
  871. DELOLD  = 0.D0                                                       ;
  872. TPS1    = 'EXTR' TPCAL ICAL                                          ;
  873. DELTAT  = TPS1 - TPSINI                                              ;
  874. * initialisation de l'indice de stockage des intégrales de concentrat
  875. LASTINT = LAST1                                                      ;
  876. 'MESS' 'Incrément de temps initial : ' DELTAT                        ;
  877. 'MESS' ' '                                                           ;
  878.  
  879. *
  880. * initialisation du terme source intégral.
  881. *
  882.  
  883. 'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
  884.    TERSC2M1 = 'NOMC'  (NOMESPL)  ('TIRE' TERSOU TPSINI)              ;
  885. 'FINSI'                                                              ;
  886.  
  887. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE')                                 ;
  888.    FLUIMPM1 = 'TIRE' FLUIMP TPSINI                                   ;
  889. 'FINSI'                                                              ;
  890.  
  891. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP')                                 ;
  892.    FLUTMPM1 = 'TIRE' FLTOTIMP TPSINI                                 ;
  893. 'FINSI'                                                              ;
  894.  
  895. 'SI' ('EXISTE' TRANSI 'MIXTES')                                      ;
  896.    FLUMMPM1 = 'TIRE' FLUMIX TPSINI                                   ;
  897. 'FINSI'                                                              ;
  898.  
  899. *
  900. * Initialisation du materiau - diffusivité
  901. *
  902.  
  903. 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHARGEME')         ;
  904.    MAT1 = 'TIRE' TRANSI . 'CARACTERISTIQUES' TPSINI                  ;
  905. 'FINSI'                                                              ;
  906.  
  907.  
  908. *
  909. * Initialisation de la quantité dissoute
  910. *
  911. 'SI' SOLUB                                                           ;
  912.      PRCI = 'NOMC' 'H' TRANSI . 'PRECIPITE' . LAST1                  ;
  913.      ANT0 = 'EXCO' 'H' ('MASQ' PRCI 'SUPERIEUR' 1.D-14) 'SCAL'       ;
  914. 'FINSI'                                                              ;
  915. *
  916. * Calcul de la fraction adsorbée initiale :
  917. *
  918. CC0  = 'NOMC' CHRG 'SCAL'                                            ;
  919. 'SI' RNONLIN                                                         ;
  920.      'SI' RNONLINL                                                   ;
  921. *         Langmuir :
  922.           FF0 = RETAR1M1 * CC0
  923.                                 /(1.D0 + ( RM1SURF * CC0 ))          ;
  924.      'FINSI'                                                         ;
  925.      'SI' RNONLINF                                                   ;
  926. *         Freundlich :
  927. *         La bidouille utilisant les masques sert à obtenir le champ point
  928. *         des signes de CC0.
  929.           a = 'MASQUE' CC0 'SUPERIEUR' 0.                            ;
  930.           b = 'MASQUE' CC0 'INFERIEUR' 0.                            ;
  931.           FF0 = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b)           ;
  932.           'OUBLI' a                                                  ;
  933.           'OUBLI' b                                                  ;
  934.      'FINSI'                                                         ;
  935. 'SINON'                                                              ;
  936. *    Isotherme linéaire :
  937.      FF0 =  RETAR1M1 * CC0                                           ;
  938. 'FINSI'                                                              ;
  939. FF  = FF0                                                            ;
  940.  
  941. *
  942. * Calcul du coefficient de retard initial :
  943. * Le facteur de retard s'exprime comme
  944. *     R(C) = 1 + ( F'(Ct+dt + eps) - F'(Ct) ) / ( Ct+dt + eps - Ct) )
  945. * C'est la méthode de la corde employée ici pour ses vertus
  946. * stabilisantes. Le petit 'eps' (EpsCord) intervient pour évacuer les
  947. * questions de conditions initiales ou de régime permanent.
  948. 'SI' RNONLIN                                                         ;
  949.      CC1 = CC0 + EpsCord                                             ;
  950.      'SI' RNONLINL                                                   ;
  951. *         Langmuir
  952.           FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 ))          ;
  953.      'SINON'                                                         ;
  954. *         Freundlich
  955.           a = 'MASQUE' CC1 'SUPERIEUR' 0.                            ;
  956.           b = 'MASQUE' CC1 'INFERIEUR' 0.                            ;
  957.           FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b)           ;
  958.           'DETRUIT' a                                                ;
  959.           'DETRUIT' b                                                ;
  960.      'FINSI'                                                         ;
  961.      RETARC = 1.D0 + ( (FF1 - FF0) / EpsCord )                       ;
  962.      'OUBLI' CC1                                                     ;
  963.      'OUBLI' FF1                                                     ;
  964. 'SINON'                                                              ;
  965.      RETARC = RETAR1                                                 ;
  966. 'FINSI'                                                              ;
  967. RETAR0 = RETARC                                                      ;
  968. 'SI' SAUVRET                                                         ;
  969.      TRANSI . 'RETARD' . 0 = RETARC                                  ;
  970. 'FINSI'                                                              ;
  971. *
  972. * Fonctions de volume :
  973. PORETSU = RETARC * POROS                                             ;
  974. DECROI0  = 'EXCO' ( (-1.) * LAMBD0 * PORETSU ) 'SCAL' 'SCAL'         ;
  975. DECROI0  = VOLU1 * DECROI0                                           ;
  976. BDECROI0 = BETA * DECROI0                                            ;
  977. BMDCROI0 = (1.D0 - BETA) * DECROI0                                   ;
  978. POROM1V = (1.D0 - POROS) * VOLU1                                     ;
  979. POROV   = POROS * VOLU1                                              ;
  980. *
  981. * Indicateur les recalculs de matrices à effectuer
  982. *
  983. TABMODI = TABLE                                                      ;
  984. TABMODI . 'POROSITE' = VRAI                                          ;
  985. TABMODI . 'CONVECTI' = VRAI                                          ;
  986. TABMODI . 'DELTAT'   = VRAI                                          ;
  987. TABMODI . 'COEF_LIN' = VRAI                                          ;
  988. TABMODI . 'DIFFUSIV' = VRAI                                          ;
  989. *
  990. * Initialisation de la table de calcul de trangeol
  991. * pour le calcul d'un itéré de l'équation de transport
  992. *
  993. CHCLIM = 'TABLE'                                                     ;
  994. GEOL1 = 'TABLE'                                                      ;
  995. * l'option abandon indique qu'en dessous de 10**-13 pour la concentr
  996. * et pour des sources ou flux n'impliquant pas des variation de C
  997. * superieurs à 10**-16 on sort 0 sans calcul
  998. 'SI' ('EXISTE' TRANSI 'SEUILCALC')                                   ;
  999. GEOL1 . 'ABANDON'           = VRAI                                   ;
  1000. GEOL1 . 'SEUILCALC'         = TRANSI . 'SEUILCALC'                   ;
  1001. 'SINON'                                                              ;
  1002. GEOL1 . 'ABANDON'           = FAUX                                   ;
  1003. GEOL1 . 'SEUILCALC'         = 1.D-30                                 ;
  1004. 'FINSI'                                                              ;
  1005.  
  1006. 'SI' ('EXISTE' TRANSI 'NUM_PECLET' )                                 ;
  1007. GEOL1 . 'NUM_PECLET'        = TRANSI . 'NUM_PECLET'                  ;
  1008. 'SINON'                                                              ;
  1009. GEOL1 . 'NUM_PECLET'        = 2.D0                                   ;
  1010. 'FINSI'                                                              ;
  1011. GEOL1 . 'CONCENTRATION'     = CHRG                                   ;
  1012. GEOL1 . 'LUMP'              = TRANSI . 'LUMP'                        ;
  1013. GEOL1 . 'TYPDISCRETISATION' = TRANSI . 'TYPDISCRETISATION'           ;
  1014. GEOL1 . 'THETA_DIFFUSION'   = TETA                                   ;
  1015. GEOL1 . 'THETA_CONVECTION'  = TETAC                                  ;
  1016. GEOL1 . 'DECENTREMENT'      = TRANSI . 'DECENTR'                     ;
  1017. GEOL1 . 'DIFFUSIVITE'       = MAT1                                   ;
  1018. *GEOL1 . 'SOLVEUR'           = TRANSI . 'METHINV' . TYPINV            ;
  1019. *GEOL1 . 'PRECONDITIONNEUR'  = TRANSI . 'METHINV' . PRECOND           ;
  1020.  
  1021.  
  1022. GEOL1 . 'MODIFICATI'        = TABMODI                                ;
  1023. GEOL1 . 'POROSITE'          = PORETSU                                ;
  1024. GEOL1 . 'SOURCE'      = 'NOMC' NOMESPL
  1025.      ('KCHT' (TRANSI.'MODELE') SCAL 'CENTRE' 0.D0)                   ;
  1026.  
  1027.  
  1028.  
  1029.  
  1030.  
  1031. *=================================
  1032. 'REPETE' BOUTPS (DCAL - ICAL + 1);    COMM 'Boucle sur le temps'     ;
  1033. *=================================
  1034. *MESS 'TEMPS1' ; TEMPS;
  1035.  
  1036. *
  1037. IPAS = ICAL + &BOUTPS - 1                                            ;
  1038. *
  1039. *-- Initialisation en-tête de boucle sur le temps
  1040. *
  1041. TPS                = 'EXTR' TPCAL IPAS                               ;
  1042. MESS 'TEMPS T =' TPS                                                 ;
  1043.  
  1044. *list ((PLACE) * 4. / 1000000.)  ;
  1045.  
  1046.  
  1047. DELTAT             = TPS - PRECED                                    ;
  1048. GEOL1 . 'DELTAT'   = DELTAT                                          ;
  1049. EPSDT              = DELTAT + DELOLD / 2.D0 * 1.D-6                  ;
  1050. TABMODI . 'DELTAT' = TABMODI . 'DELTAT'
  1051.                      'OU' ( DELOLD 'NEG' DELTAT EPSDT ) ;
  1052.  
  1053. *                                                           CONVECTION
  1054. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' )                                ;
  1055.    GEOL1 . 'CONVECTION' = TRANSI . 'CONVECTION'                      ;
  1056. 'FINSI'                                                              ;
  1057.  
  1058. *                                                      VITESSE CENTRE
  1059. 'SI' ( 'EXISTE' TRANSI 'VITELEM' )                                   ;
  1060.    GEOL1 . 'VITELEM' = TRANSI . 'VITELEM'                            ;
  1061. 'FINSI'                                                              ;
  1062.  
  1063. *                                                               ALPHAL
  1064. 'SI' ( 'EXISTE' TRANSI 'ALPHAL' )                                    ;
  1065.    GEOL1 . 'ALPHAL' = TRANSI . 'ALPHAL'                              ;
  1066. 'FINSI'                                                              ;
  1067.  
  1068. 'SI' ( 'EXISTE' TRANSI 'ALPHAT')                                     ;
  1069.    GEOL1 . 'ALPHAT' = TRANSI . 'ALPHAT'                              ;
  1070. 'FINSI'                                                              ;
  1071.  
  1072. *
  1073. *- Calcul de la contribution des termes sources
  1074. *
  1075. TERSCE = 'MANU' 'CHPO' MCENT 1 'SOUR' 0.                             ;
  1076. TERSCE = 'NOMC' NOMESPL TERSCE                                       ;
  1077.  
  1078. tpsm1 =  'EXTR' TPCAL (IPAS - 1)   ;
  1079. * Terme source propre :
  1080. 'SI' ( 'EXISTE' TRANSI 'SOURCE' )                                    ;
  1081.     'SI' (tpsm1 &lt;EG startsou)                                        ;
  1082. *       la source varie encore et est non nulle
  1083.         zozo = 'TIRE' TERSOU TPS                                     ;
  1084.     'SINON'                                                          ;
  1085. *       la source est nulle on garde l'intégrale précédente
  1086.         zozo = 'COPIER' TERSC2M1                                     ;
  1087.     'FINSI'                                                          ;
  1088.     TERSC2 =  'NOMC' NOMESPL  zozo                                   ;
  1089.     'DETRUIT' zozo                                                   ;
  1090.     zozo = TERSC2 '-' TERSC2M1                                       ;
  1091.     TERSCV = zozo '/' DELTAT                                         ;
  1092.     'DETRUIT' zozo                                                   ;
  1093.     zozo = TERSCV + TERSCE                                           ;
  1094.     'DETRUIT' TERSCE                                                 ;
  1095.     TERSCE = 'COPIER' zozo                                           ;
  1096.     'DETRUIT' zozo                                                   ;
  1097.     'DETRUIT' TERSCV                                                 ;
  1098. 'FINSI'                                                              ;
  1099.  
  1100. *
  1101. * Terme source de décroissance explicite :
  1102. * du soluté et de l'adsorbat
  1103. 'SI' DECRO                                                           ;
  1104. * modif GBM suppression de FF, revient à compter R * R
  1105. *     TERSC3 = ( ('NOMC' CHRG 'SCAL') + FF ) * DECROI0               ;
  1106.      zozo   = 'NOMC' CHRG 'SCAL'                                     ;
  1107.      zuzu = zozo * DECROI0                                           ;
  1108.      'DETRUIT' zozo                                                  ;
  1109.      TERSC3 =  'NOMC' NOMESPL zuzu                                   ;
  1110.      'DETRUIT' zuzu                                                  ;
  1111.      zozo = TERSC3 + TERSCE                                          ;
  1112.      'DETRUIT' TERSCE                                                ;
  1113.      TERSCE = 'COPIER' zozo                                          ;
  1114.      'DETRUIT' zozo                                                  ;
  1115. 'SINON'                                                              ;
  1116. *    pas de filiation explicite
  1117.      TERSC3 = 0.D0 * TERSCE                                          ;
  1118. 'FINSI'                                                              ;
  1119.  
  1120.  
  1121. * Chaine de filiation source venant du pere
  1122. 'SI' ('EXISTE' TRANSI 'FILIATION')                                   ;
  1123. *   TERSC4 =  VOLU1 * (TRANSI . 'FILIATION' . (IPAS '-' 1)
  1124. *         '-' TRANSI . 'FILIATION' . (IPAS '-' 2)) '/' DELTAT         ;
  1125.    zozo   = TRANSI . 'FILIATION' . (IPAS '-' 1)
  1126.          '-' TRANSI . 'FILIATION' . (IPAS '-' 2)                     ;
  1127. *   GBM attention à cette destruction
  1128.    'DETRUIT' (TRANSI . 'FILIATION' . (IPAS '-' 2))                   ;
  1129.    zuzu = zozo '/' DELTAT                                            ;
  1130.    'DETRUIT' zozo                                                    ;
  1131.    zozo = VOLU1 * zuzu                                               ;
  1132.    'DETRUIT' zuzu                                                    ;
  1133.    TERSC4 = 'NOMC'  NOMESPL zozo                                     ;
  1134.    'DETRUIT' zozo                                                    ;
  1135.    zozo =  TERSC4 '+' TERSCE                                         ;
  1136.    TERSCE = 'COPIER' zozo                                            ;
  1137.    'DETRUIT' zozo                                                    ;
  1138.    'DETRUIT' TERSC4                                                  ;
  1139. 'FINSI'                                                              ;
  1140.  
  1141. *
  1142. * Gestion d'un CSAT variable (t) - chargement
  1143. *
  1144. 'SI' ('EXISTE' TRANSI 'LIMITE_SOLUBILITE')                           ;
  1145. 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHARGEME')        ;
  1146. * GBM ON PEUT OPTIMISER EN TESTANT SI LIMSOL CHANGE REELEMENT EN ECART
  1147. * RELATIF PAR RAPPORT AU TEMPS PRECEDENT
  1148.    TABMODI . 'LIMSOL' = VRAI                                         ;
  1149.    LIMSOL = 'TIRE' TRANSI . 'LIMITE_SOLUBILITE' TPS                  ;
  1150.    'SI' (SOLUP)                                                      ;
  1151. *  Là où la précipitation est nulle (codis =0) on choisit
  1152. *  une limite de solubilité élevée à des fins d'optimisation
  1153. *  algorithmique
  1154.       dum = 'MASQUE' CODIS INFERIEUR 1.D-30                          ;
  1155.       LIMSOL = ((1.D0 '-' dum) * LIMSOL) '+' ('NOMC' 'H' (1.D3*dum)) ;
  1156.    'FINSI'                                                           ;
  1157. 'SINON'                                                              ;
  1158.    TABMODI . 'LIMSOL' = FAUX                                         ;
  1159. 'FINSI'                                                              ;
  1160. 'FINSI'                                                              ;
  1161.  
  1162. *
  1163. * Gestion de la diffusivité variable en temps - chargement. Donnée
  1164. * explicite, ie au temps n lors du traitement vers n+1
  1165. *
  1166.  
  1167. 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHARGEME')         ;
  1168. *  On peut optimiser en testant la variation réelle de la diffusion
  1169. *  elle pourrait en effet être constante par palliers
  1170.    TABMODI . 'DIFFUSIV' = VRAI                                       ;
  1171.    MAT1 = 'TIRE' TRANSI . 'CARACTERISTIQUES' TPS                     ;
  1172. * Estelle : affectation de la diffusion dans geol1
  1173.    GEOL1 . 'DIFFUSIVITE'       = MAT1                                ;
  1174. 'SINON'                                                              ;
  1175.    TABMODI . 'DIFFUSIV' = FAUX                                       ;
  1176. 'FINSI'                                                              ;
  1177.  
  1178. *
  1179. *- Incorporation des CLs
  1180. *
  1181. 'SI' ('EXISTE' TRANSI 'TRACE_IMPOSE')                                ;
  1182.    CHARIMPO = 'TIRE' CHDIRI TPS                                      ;
  1183.    CHCLIM . 'DIRICHLET' = 'NOMC' NOMESPL CHARIMPO                    ;
  1184.    'DETRUIT' CHARIMPO                                                ;
  1185. 'FINSI'                                                              ;
  1186.  
  1187. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE')                                 ;
  1188.    'SI' (tpsm1 '&lt;EG' startflu)                                       ;
  1189.        FLUIMPO = 'TIRE' FLUIMP TPS                                   ;
  1190.    'SINON'                                                           ;
  1191.        FLUIMPO = 'COPIER' FLUIMPM1                                   ;
  1192.    'FINSI'                                                           ;
  1193.    FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET'            ;
  1194.    'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1195.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1196.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  1197. *     On écrit les flux sous forme intégrale sauf en explicite et
  1198. *     kranck-Nickholson pour EFMH. On les différentie ici
  1199.       zozo = FLUIMPO '-' FLUIMPM1                                    ;
  1200.       FLUMP = zozo '/' DELTAT                                        ;
  1201.       'DETRUIT' zozo                                                 ;
  1202.    'SINON'                                                           ;
  1203.       FLUMP = 'COPIER' FLUIMPO                                       ;
  1204.    'FINSI'                                                           ;
  1205.    CHCLIM . 'NEUMANN' = 'NOMC' NOMESPL FLUMP                         ;
  1206.    'DETRUIT' FLUMP                                                   ;
  1207. 'FINSI'                                                              ;
  1208.  
  1209.  
  1210. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP')                                 ;
  1211.    'SI' (tpsm1 '&lt;EG' startflt)                                       ;
  1212.        FLUTMPO = 'TIRE' FLTOTIMP TPS                                 ;
  1213.    'SINON'                                                           ;
  1214.        FLUTMPO = 'COPIER' FLUTMPM1                                   ;
  1215.    'FINSI'                                                           ;
  1216.    FLUTMPO = 'CHAN' 'ATTRIBUT' FLUTMPO 'NATURE' 'DISCRET'            ;
  1217.    'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1218.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1219.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  1220. *     On écrit les flux sous forme intégrale sauf en explicite et
  1221. *     kranck-Nickholson pour EFMH. On les différentie ici
  1222.       zozo = FLUTMPO '-' FLUTMPM1                                    ;
  1223.       FLUTMP = zozo '/' DELTAT                                       ;
  1224.       'DETRUIT' zozo                                                 ;
  1225.    'SINON'                                                           ;
  1226.       FLUTMP = 'COPIER'  FLUTMPO                                     ;
  1227.    'FINSI'                                                           ;
  1228.    CHCLIM . 'FLUTOTAL' = 'NOMC' NOMESPL  FLUTMP                      ;
  1229.    'DETRUIT' FLUTMP                                                  ;
  1230. 'FINSI'                                                              ;
  1231.  
  1232.  
  1233. 'SI' ('EXISTE' TRANSI 'MIXTES')                                      ;
  1234.    'SI' (tpsm1 '&lt;EG' startmix)                                       ;
  1235.        FLUMMPO = 'TIRE' FLUMIX TPS                                   ;
  1236.    'SINON'                                                           ;
  1237.        FLUMMPO = 'COPIER' FLUMMPM1                                   ;
  1238.    'FINSI'                                                           ;
  1239.    FLUMMPO = 'CHAN' 'ATTRIBUT' FLUMMPO 'NATURE' 'DISCRET'            ;
  1240.    'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1241.            'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1242.                  'ET' ('EGA' (teta * tetac) 1.D0)))                  ;
  1243. *     On écrit les flux sous forme intégrale sauf en explicite et
  1244. *     kranck-Nickholson pour EFMH. On les différentie ici
  1245.       zozo = FLUMMPO '-' FLUMMPM1                                    ;
  1246.       FLUMMP = zozo '/' DELTAT                                       ;
  1247.       'DETRUIT' zozo                                                 ;
  1248.    'SINON'                                                           ;
  1249.       FLUMMP = 'COPIER'  FLUMMPO                                     ;
  1250.    'FINSI'                                                           ;
  1251.    CHCLIM . 'FLUMIXTE' = TABLE                                       ;
  1252.    CHCLIM . 'FLUMIXTE' . 'VAL' = 'NOMC' NOMESPL  FLUMMP              ;
  1253.    CHCLIM . 'FLUMIXTE' . 'COEFA' = TRANSI . 'MIXCOFA'                ;
  1254.    CHCLIM . 'FLUMIXTE' . 'COEFB' = TRANSI . 'MIXCOFB'                ;
  1255.    'DETRUIT' FLUMMP                                                  ;
  1256. 'FINSI'                                                              ;
  1257.  
  1258.  GEOL1 . 'CLIMITES'          = CHCLIM                                ;
  1259. * initialisation de la table de préconditionnement et eventuellement
  1260. * des traces de concentration en EFMH au premier temps
  1261. * A virer dès modification arguments transgeol
  1262.  
  1263. 'SI' ('EGA' &BOUTPS 1)                                               ;
  1264. *  VERRUE d'initialisation à virer
  1265. *  il s'agit d'une avancée de DT = 0.
  1266.    GEOL1 . 'DELTAT' = 1.D-15                                         ;
  1267.    GEOL1 . 'METHINV'           = TRANSI . 'METHINV'                  ;
  1268. * Dans ce cas particulier d'un dt petit, precond diagonal, bcgs
  1269.    GEOL1 . 'SOLVEUR' = 3                                             ;
  1270.    GEOL1 . 'PRECONDITIONNEUR' = 1                                    ;
  1271.    GEOL1 GEOL2 = TRANGEOL transi.'MODELE' GEOL1                      ;
  1272. * On remet les choix utilisateur pour la suite
  1273.    GEOL2 . 'METHINV' . 'TYPINV' = TRANSI . 'METHINV' . 'TYPINV'      ;
  1274.    GEOL2 . 'METHINV' . 'PRECOND' = TRANSI . 'METHINV' . 'PRECOND'    ;
  1275.    SI (EGA (TRANSI . 'METHINV' . 'PRECOND') 8)                       ;
  1276.       GEOL2 . 'METHINV' . 'ILUTDTOL' = 1.D-2                         ;
  1277.    FINSI                                                             ;
  1278.  
  1279. *  remise à jour des paramètres d'entrée de EFMHTgen
  1280.    GEOL1 . 'DELTAT'   = DELTAT                                       ;
  1281.  
  1282.    TABMODI . 'POROSITE' = FAUX                                       ;
  1283.    TABMODI . 'CONVECTI' = FAUX                                       ;
  1284.    TABMODI . 'COEF_LIN' = FAUX                                       ;
  1285. * TABMODI . DIFFUSIV est écrasé à FAUX meme si chargement car
  1286. * un premier appel à trangeol a déjà été fait pour initialiser les
  1287. * préconditionnements, seul le pas de temps change ici. C'est vrai
  1288. * uniquement pour le premier pas de temps.
  1289.    TABMODI . 'DIFFUSIV' = FAUX                                       ;
  1290.    TABMODI . 'DELTAT'   = VRAI                                       ;
  1291.  
  1292. 'FINSI'                                                              ;
  1293.  
  1294. *
  1295. *
  1296. *|-------------------------------------------------------------------|
  1297. *| Boucle de Retard Non Linéaire                                     |
  1298. *|-------------------------------------------------------------------|
  1299. *
  1300. * Boucle de Picard pour la prise en compte
  1301. * d'isothermes d'adsorption non linéaires : isothermes de Langmuir et
  1302. * de Freundlich. Ce point nécessite d'effectuer une boucle sur
  1303. * le facteur de retard qui devient fonction de la concentration
  1304. * en solution.
  1305. *
  1306. *------------------------
  1307. 'REPETE' BOURNL NPicard ;     COMM 'Boucle sur le coef. de retard'   ;
  1308. *------------------------
  1309. *
  1310. IRNL = &BOURNL                                                       ;
  1311. *
  1312. * Quatre grands blocs selon que l'on traite :
  1313. * -------------------------------------------
  1314. *  -1-  sans limite de solubilité ou avec dissolution arbitraire
  1315. *  -2-  avec limite de solubilité et dissolution d'ordre 1
  1316. *  -3-  avec limite suivant méthode de pénalisation,
  1317. *  -4-  avec limite suivant schéma prédicteur/correcteur itératif
  1318.  
  1319.  
  1320. *|----------------------------------------------------|
  1321. *| Pas de limite de solubilité ou Dissolution imposée |
  1322. *|----------------------------------------------------|
  1323. 'SI' (('NON' SOLUB) 'OU' SOLUA)                                      ;
  1324. *
  1325. * On recalcule la matrice hybride si DeltaT ou le coefdt ont changé :
  1326. 'SI' (TABMODI.'DELTAT' 'OU' TABMODI . 'POROSITE')                    ;
  1327.  
  1328. * GBM FAIRE MEILLEUR TEST SI DELTAT BOUGE SEULEMENT TABMODI.PORO
  1329. * NE CHANGE PAS
  1330.  
  1331. *    Coefficients correspondant à la décroissance expl. et implicite :
  1332.      NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
  1333.      DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
  1334. *    Calcul du coefficient devant (Ct+dt - Ct)/dt
  1335.      COFDC = PORETSU * DEN0                                          ;
  1336. *    On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1337.      GEOL1 . 'POROSITE' = COFDC                                      ;
  1338. *    On indique que le coef devant DC/DT a changé
  1339.      TABMODI . 'POROSITE' = VRAI                                     ;
  1340. 'FINSI'                                                              ;
  1341. *
  1342. 'SI' SOLUA                                                           ;
  1343. *    Evaluation des termes sources de dissolution :
  1344. *    Dissolution arbitraire (ordre 0 par exemple), limitée au total
  1345. *    précipité :
  1346.      DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1                        ;
  1347.      DISSOLP = 'TIRE' DISIMP TPS                                     ;
  1348.      TERSC4 = 0.5 * (DISSOLT + DISSOLP - ('ABS' (DISSOLT - DISSOLP)));
  1349.      TERSC4  = 'EXCO' 'H' TERSC4 'H'                                 ;
  1350.      TERSCT  = TERSCE + TERSC4                                       ;
  1351. 'SINON'                                                              ;
  1352.      TERSCT  = TERSCE                                                ;
  1353. 'FINSI'                                                              ;
  1354.  
  1355. * On charge le terme source dans GEOL1
  1356.  
  1357. GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSCT                             ;
  1358.  
  1359. * options de calcul
  1360. GEOL1 . 'MODIFICATI' = TABMODI                                       ;
  1361.  
  1362. *- Résolution - APPEL TRANGEOL
  1363. * OPTIMISER STOCKAGE DOUBLE DES MATRICES PRECOND ?????
  1364. * faire transgeol mat1 mat2 et stocké tracini tracfin, conini concfin...
  1365. * on pourrait alors écraser les matrices sytématiquement et faire un
  1366. * advance comme philippe : stocke cfin dans cini
  1367. * transgeol doit pouvoir accepter matrice vide precond en option
  1368. GEOLPF1  GEOLPF2 = TRANGEOL transi.'MODELE' GEOL1 GEOL2              ;
  1369. CHA2 = GEOLPF1 . 'CONCENTRATION'                                     ;
  1370. FLU2 = GEOLPF1 . 'FLUXDIFF'                                          ;
  1371. FLUCO2 = GEOLPF1 . 'FLUXCONV'                                        ;
  1372.  
  1373. * Remise à jour des coef de recalcul des matrices - dans le point
  1374. * fixe, ces propriétés ne changent pas
  1375. TABMODI . 'POROSITE' = FAUX                                          ;
  1376. TABMODI . 'CONVECTI' = FAUX                                          ;
  1377. TABMODI . 'DELTAT' = FAUX                                            ;
  1378. TABMODI . 'COEF_LIN' = FAUX                                          ;
  1379. TABMODI . 'DIFFUSIV' = FAUX                                          ;
  1380.  
  1381. 'SI' SOLUA                                                           ;
  1382.      DIS1 = 'EXCO' 'H' (TERSC4 * DELTAT) 'H'                         ;
  1383.      PRE2 = ((NUM0 * PRCI) - (DIS1 / VOLU1)) / DEN0                  ;
  1384. 'FINSI'                                                              ;
  1385.  
  1386.  
  1387. * 'Fin module dissolution arbitraire'
  1388. 'FINSI'                                                              ;
  1389.  
  1390. *list ((PLACE) * 4. / 1000000.)  ;
  1391.  
  1392. *|---------------------------|
  1393. *|   Dissolution d'ordre 1   |
  1394. *|---------------------------|
  1395. 'SI' SOLUP;
  1396.  
  1397. * On recalcule la matrice hybride si DeltaT ou le coefdt ont changé ou
  1398. * la limite de solubilité a changé:
  1399. 'SI' (TABMODI.'DELTAT' 'OU' TABMODI . 'POROSITE'
  1400.           'OU' TABMODI . 'LIMSOL')                                   ;
  1401. *    présence de précipitation dissolution
  1402.      zozo = 'MASQ' PRCI 'SUPERIEUR' 1.D-14                           ;
  1403.      ANTP = 'EXCO' 'H' zozo  'SCAL'                                  ;
  1404.      'DETRUIT' zozo                                                  ;
  1405.  
  1406.      zozo = (1.D0 '+' 1.D-10) * LIMSOL                               ;
  1407.      zuzu = 'MASQ' CHRG 'SUPERIEUR' zozo                             ;
  1408.      'DETRUIT' zozo                                                  ;
  1409.      ANTC = 'EXCO' 'H' zuzu 'SCAL'                                   ;
  1410.      'DETRUIT' zuzu                                                  ;
  1411.  
  1412.      zozo =  ANTP '+' ANTC                                           ;
  1413.      'DETRUIT' ANTC                                                  ;
  1414.      ANT0 =  'MASQ' zozo 'SUPERIEUR' 0.1                             ;
  1415.      'DETRUIT' zozo                                                  ;
  1416.  
  1417. *    Coefficients correspondant à la décroissance expl. et implicite :
  1418.      NUM0  = (BETA - 1.) * LAMBD0 * DELTAT + 1.                      ;
  1419.      DEN0  = LAMBD0 * BETA * DELTAT + 1.                             ;
  1420.  
  1421. *    Terme implicite de dissolution, schéma Euler
  1422.      zozo = DELTAT * CODIS                                           ;
  1423.      DEN1 = 'KOPS' zozo * ANT0                                       ;
  1424.      'DETRUIT' zozo                                                  ;
  1425.  
  1426. *    Calcul du coefficient devant (Ct+dt - Ct)/dt
  1427. *    solub en poros*codis Id et non Retard*poros*codis
  1428.      zozo = 'KOPS' PORETSU * DEN0                                    ;
  1429.      zeze = 'KOPS' POROS * DEN1                                      ;
  1430.      COFDC = 'KOPS' zeze '+' zozo                                    ;
  1431.      'DETRUIT' zozo                                                  ;
  1432.      'DETRUIT' zeze                                                  ;
  1433.  
  1434. *    On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1435.      GEOL1 . 'POROSITE' = COFDC                                      ;
  1436. *    On indique que le coef devant DC/DT a changé
  1437.      TABMODI . 'POROSITE' = VRAI                                     ;
  1438.  
  1439. 'FINSI'                                                              ;
  1440.  
  1441. *  ON charge dans une variable point fixe le precipité
  1442. *  initial
  1443. PRCIPF = 'COPIER' PRCI                                               ;
  1444. *  On initialise le terme correspondant à une dissolution
  1445. *  complete du précipité le derniere itération avant sa disparition
  1446. *  Pour eviter des précipités négatifs
  1447. SRCPRCI = 0.D0 * PRCI                                                ;
  1448.  
  1449. *---------------
  1450. 'REPETER' BOU8 ;     COMM 'Boucle dissolution ordre 1'               ;
  1451. *---------------
  1452.  
  1453.  
  1454. * terme source modifié en raison du coef de solubilité :
  1455. * on calcul deltat codis (cn+1 - cn) /deltat, astuce informatique
  1456. * pour créer le terme deprécipitation. Il faut retirer le - codis cn
  1457. * au second membre : - codis * posos * volume (car source integrale)
  1458. *  TERSC5 = TERSCE '-' ( CODIS * ANT0 * POROS * VOLU1 * CHRG)        ;
  1459. * On a également la source liée au terme codis poros limsol
  1460. *   TERSC5 = TERSC5 '+' ( CODIS * ANT0 * POROS * VOLU1 * LIMSOL)      ;
  1461.  
  1462.    zozo = 'KOPS' CODIS * ANT0    ;
  1463.    zeze = 'KOPS' zozo * POROV    ;
  1464.    'DETRUIT' zozo         ;
  1465.    zuzu = LIMSOL '-' CHRG ;
  1466.    zozo = zeze * zuzu     ;
  1467.    'DETRUIT' zeze         ;
  1468.    'DETRUIT' zuzu         ;
  1469.    TERSC5 = TERSCE '+' zozo ;
  1470.    'DETRUIT' zozo         ;
  1471.  
  1472. * On a également la source liée a la dissolution
  1473. * instantanée du précipité là 'OU' il est devenu négatif
  1474. * artifice numérique permettant d'avoir toujours un
  1475. * précipité négatif. Bilan de matiere verifie, erreur en
  1476. * O(DT) su le reste de la physique
  1477. *   TERSC5 = TERSC5 '+' (VOLU1 * SRCPRCI)                             ;
  1478.  
  1479.    zozo = VOLU1 * SRCPRCI                                            ;
  1480.    zuzu = TERSC5 '+' zozo                                            ;
  1481.    'DETRUIT' zozo                                                    ;
  1482.    'DETRUIT' TERSC5                                                  ;
  1483.    TERSC5 = 'COPIER' zuzu                                            ;
  1484.    'DETRUIT' zuzu                                                    ;
  1485.  
  1486. * le fait que codis puisse etre gigantesque est sans influence
  1487. * sur la précision car le précipité est calulé par des bilans de flux
  1488. * et n'utilisera plus ces valeurs.
  1489.  
  1490.  
  1491. * On charge le terme source dans GEOL1.
  1492.  
  1493.    GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSC5                          ;
  1494.  
  1495. * options de calcul
  1496.    GEOL1 . 'MODIFICATI' = TABMODI                                    ;
  1497.  
  1498. *- Résolution - APPEL TRANGEOL
  1499. * OPTIMISER STOCKAGE DOUBLE DES MATRICES PRECOND ?????
  1500. * faire transgeol mat1 mat2 et stocké tracini tracfin, conini concfin...
  1501. * on pourrait alors écraser les matrices sytématiquement et faire un
  1502. * advance comme philippe : stocke cfin dans cini
  1503. * transgeol doit pouvoir accepter matrice vide precond en option
  1504. *list ((PLACE) * 4. / 1000000.)  ;
  1505.    GEOLPF1  GEOLPF2 = TRANGEOL transi.'MODELE' GEOL1 GEOL2           ;
  1506. *list ((PLACE) * 4. / 1000000.)  ;
  1507.    CHA2 = GEOLPF1 . 'CONCENTRATION'                                  ;
  1508.    FLU2 = GEOLPF1 . 'FLUXDIFF'                                       ;
  1509.    FLUCO2 = GEOLPF1 . 'FLUXCONV'                                     ;
  1510.  
  1511. *  remise à jour de XINIT pour le point fixe. Pourra etre retirer
  1512. *  quand transgeol stockera les état initiaux- finaux et intermédiaire
  1513. *  pour un point fixe. Permet uniquement d'accélérer la convergence
  1514. *  des méthodes itératives. Sans incidence sur les résultats physiques
  1515.    GEOL2 . 'METHINV' . 'XINIT' = GEOLPF2 . 'METHINV' . 'XINIT'       ;
  1516.    'SI' ('EXISTE' GEOL2 TRACE_CONC)                                  ;
  1517.  
  1518. *     On conserva absolument la trace en EFMH car c'est la vraie
  1519. *     inconnue et on doit garder celle du début du point fixe
  1520. *     (sinon on incrémente en temps dans le point fixe).
  1521. **      'MESSAGE' 'yoyououou';
  1522.       dum = GEOL2 . 'TRACE_CONC'                                     ;
  1523.       GEOL2 = GEOLPF2                                                ;
  1524.       GEOL2 . 'TRACE_CONC' = dum                                     ;
  1525.    'FINSI'                                                           ;
  1526.  
  1527. * Remise à jour des coef de recalcul des matrices
  1528.    TABMODI . 'POROSITE' = FAUX                                       ;
  1529.    TABMODI . 'CONVECTI' = FAUX                                       ;
  1530.    TABMODI . 'DELTAT' = FAUX                                         ;
  1531.    TABMODI . 'COEF_LIN' = FAUX                                       ;
  1532.    TABMODI . 'DIFFUSIV' = FAUX                                       ;
  1533.  
  1534. * SI ANT0 nul partout alors pas de prec dissolution donc le precipité
  1535. * est forcément nul partout. On ne fait pas de calcul.
  1536.  
  1537.   'SI' (('MAXIMUM' ANT0) 'EGA' 0)                                    ;
  1538. *    on met à 0 le précipité, on utilise le champ de concentration
  1539. *    pour avoir la géométrie du support.
  1540.      PRE2 = 0.D0 * CHRG                                              ;
  1541.   'SINON'                                                            ;
  1542.  
  1543. * Calcul du précipité par bilan des variations des quantités et des
  1544. * flux diffusifs et convectifs sur la maille. On n'utilise surtout
  1545. * pas la loi de dissolution codis * (C - limsol) car si codis
  1546. * est grand (pénalisation) le calcul est très imprécis.
  1547. * le précipité est en mole/volume de milieu solide.
  1548. * On décompose le calcul en plusieurs termes :
  1549. * Faire ATTENTION rien n'est intuitif. Si des schémas sont modifiés
  1550. * en amont, que les termes sources sont modifiés, il faudra
  1551. * s'assurer que cela ne modifie pas les lignes suivantes.
  1552.  
  1553. * la variation de concentration dans le milieu (et du sorba)
  1554. * poros * retar * (delta COncentration) * volum / deltat
  1555.  
  1556. *list ((PLACE) * 4. / 1000000.)  ;
  1557.    PRR1 = (-1.D0 '/' DELTAT) * VOLU1      ;
  1558.    PRR2 = 'KOPS' PRR1 * PORETSU                  ;
  1559.    'DETRUIT'  PRR1;
  1560.    zozo = 'KOPS' CHA2 - CHRG ;
  1561.    PRR1 = 'KOPS' PRR2 * zozo     ;
  1562.    'DETRUIT'  zozo;
  1563.    'DETRUIT'  PRR2;
  1564.  
  1565.  
  1566. * la variation de concentration disparue par filiation. Il faut
  1567. * écrire les équations de bilan pour se convaincre que tout n'est pas
  1568. * compté plusieurs fois de trop: - lambda * poros * retard * volume *
  1569. * Cn+1  multiplié par beta (valeur du theta schéma pour décroissance).
  1570. * plkus partie explicite : - lambda * poros * retard * volume *
  1571. * Cn * (1 - beta). Le signe '-' est dans DECROI0
  1572.  
  1573.    zozo = 'KOPS' BDECROI0 * CHA2 ;
  1574.    PRR2 = 'KOPS' PRR1 '+' zozo ;
  1575.    'DETRUIT'  zozo;
  1576.    'DETRUIT'  PRR1;
  1577.    zozo = 'KOPS' BMDCROI0 * CHRG ;
  1578.    PRR1 = 'KOPS' PRR2 '+' zozo;
  1579.    'DETRUIT'  zozo;
  1580.    'DETRUIT'  PRR2;
  1581.  
  1582.  
  1583. * Le terme source de C qui comprend d'ailleurs, la filiation venant
  1584. * du pere (concentration + précipité cumulés), et les termes
  1585. * source de l'utilisateur issu de l'indice SOURCE de la table
  1586. * d'entrée. La multiplication par le volume de l'élément
  1587. * est déjà effectuée. On a retirer TERSC3 (source de filiation)
  1588. * déjà comptée ligne précédente. Encore une fois faire très
  1589. * attention lors de modif de ces lignes, on peut s'y perdre.
  1590. * A la limite tout coder en Fortran
  1591.  
  1592.    zaza = 'KOPS' TERSCE - TERSC3 ;
  1593.    PRR2 = 'KOPS' PRR1 + zaza     ;
  1594.    'DETRUIT'  zaza;
  1595.    'DETRUIT'  PRR1 ;
  1596.  
  1597. * le variation par filiation du précipité - partie explicite
  1598. * du theta schéma de décroissance radioactive : - (1 - poros)
  1599. * fois  LAMBDA * (1 - BETA) * volume maille * precipité (instant n)
  1600.  
  1601.    PRR1 = 'KOPS' POROM1V * PRCIPF ;
  1602.    zozo = (LAMBD0 * (1.D0 '-' BETA)) * PRR1 ;
  1603.    'DETRUIT'  PRR1 ;
  1604.    PRR1 = 'KOPS' PRR2 - zozo ;
  1605.    'DETRUIT'  zozo;
  1606.    'DETRUIT'  PRR2 ;
  1607.  
  1608. * Les flux diffusifs et convectifs intégrés au travers des faces de
  1609. * l'élément. flux SORTANT donc diminuant le précipité
  1610. *   fluttot = (TETA * FLU2) + (1.D0 '-' TETA) * FLU0) '+'
  1611. *            ((TETAC * FLUCO2) '+' ((1.D0 '-' TETAC) * FLUCO0))       ;
  1612.  
  1613.    zozo = TETA '*' FLU2 ;
  1614.    zizi = ((1.D0 '-' TETA) * FLU0) ;
  1615.    zaza = 'KOPS' zizi + zozo ;
  1616.    'DETRUIT'  zizi;
  1617.    'DETRUIT'  zozo;
  1618.  
  1619.    zozo = TETAC * FLUCO2 ;
  1620.    zizi = ((1.D0 '-' TETAC) * FLUCO0) ;
  1621.    zeze = 'KOPS' zizi + zozo;
  1622.    'DETRUIT'  zizi;
  1623.    'DETRUIT'  zozo;
  1624.  
  1625.    fluttot1 = 'KOPS' zeze + zaza;
  1626.    'DETRUIT'  zeze;
  1627.    'DETRUIT'  zaza;
  1628.  
  1629. *  DIVU veut une composante FLUX
  1630.  
  1631.    fluttot = 'NOMC' 'FLUX' fluttot1  ;
  1632.    'DETRUIT' fluttot1;
  1633.  
  1634.    divflt = 'DIVU' MODDARCY fluttot MCHYB ;
  1635.    'DETRUIT' fluttot;
  1636.  
  1637.    divfll = 'NOMC' 'SCAL' divflt;
  1638.    'DETRUIT' divflt;
  1639.  
  1640.    PRR2 =  'KOPS' PRR1 '-' divfll ;
  1641.    'DETRUIT' divfll;
  1642.    'DETRUIT' PRR1;
  1643.  
  1644.  
  1645.  
  1646. * La quantité de précipité à l'instant précédent :
  1647. * (1 '-' POROS) PRCI * VOLU1 '/' DELTAT
  1648. * AFIN d'éviter les problemes lorsque DELTAT = 0
  1649. * on multiplie tout par DELTAT, donc ne pas le faire 'PLUS' loin
  1650. *  PRE2 = (deltat * PRE2) '+' ((1.D0 '-' POROS) * VOLU1 * PRCIPF)    ;
  1651.  
  1652.    PRR1 = deltat * PRR2 ;
  1653.    'DETRUIT' PRR2 ;
  1654.    zozo = 'KOPS' POROM1V * PRCIPF ;
  1655.    PRR2 = 'KOPS' PRR1 '+' zozo ;
  1656.    'DETRUIT' zozo;
  1657.    'DETRUIT' PRR1;
  1658.  
  1659.  
  1660. * On résout alors (vol * (1 - poros) /deltat) + (1 - poros) vol lambda
  1661. * beta     le tout fois precipité instant N+1 = PRE2 calculé au dessus
  1662. *
  1663. *   coff = ((DELTAT * LAMBD0 * BETA) '+' 1.D0) * (1.D0 '-' POROS)
  1664. *       * VOLU1                                                       ;
  1665. *   PRE2 = PRE2 '/' coff                                              ;
  1666.  
  1667.  
  1668.     coff = ((DELTAT * LAMBD0 * BETA) '+' 1.D0) * POROM1V;
  1669.     PRR1 = 'KOPS' PRR2 '/' coff;
  1670.     'DETRUIT' coff;
  1671.     'DETRUIT' PRR2;
  1672.  
  1673.  
  1674. *     ATTENTION SI IL Y A DU PR2CIPITE, POROSITE = 1 INTERDIT CAR
  1675. *     DIT QU'il NY A PAS DE SOLIDE
  1676.  
  1677. *  On multiplie par la fonction caractéristique de présence du précipité
  1678. *  inutile en pratique si les schémas sont conservatifs sur l'élément
  1679. *  c'est mis pour ne pas propager d'erreurs machines.
  1680.  
  1681.    PRE2 = 'KOPS' ANT0 * PRR1                ;
  1682.    'DETRUIT' PRR1 ;
  1683.  
  1684.  
  1685.    zozo = 'NOMC' NOMESPL PRE2 ;
  1686.    PRE2 = 'COPIER' zozo;
  1687.    'DETRUIT' zozo ;
  1688.  
  1689. *   Fin du calcul du précipité par bilan
  1690.    'FINSI'                                                           ;
  1691.  
  1692. *    mess 'maximum precipite' (maxi pre2);
  1693. *    mess 'maximum concentration' (maxi cha2);
  1694.  
  1695. *    Critère de sortie et mise à jour de la nouvelle distribution :
  1696. *  On sort la distribution de précipité n'a pas changé (ie le précpité
  1697. * est toujours positif aux memes endoits) et que la concentration ne
  1698. * dépasse pas limsol dans la partie ou il n'y a pas de précipité
  1699.  
  1700. *    LTI2 = 'CHAINE' 'Front 1D-h temps '    ;
  1701. *   'TITR' LTI2                                 ;
  1702. *   drmil = 'QUELCONQUE' 'SEG2' ('EXTRAIRE' CHA2 maillage) ;
  1703. *   drmil = 'INVERSE' drmil                     ;
  1704. *   AV2  = 'EVOL' 'ROUG' 'CHPO' CHA2 'H'  drmil ;
  1705. *   'DESS' AV2 'MIMA'                           ;
  1706. *   AV2  = 'EVOL' 'ROUG' 'CHPO' PRE2 'H'  drmil ;
  1707. *   'DESS' AV2 'MIMA'  'NCLK'                   ;
  1708.  
  1709.  
  1710. *  ANT2 caracterise la présence de dissolution quant vaut 1
  1711.    zozo = 'MASQ' PRE2 'SUPERIEUR' 1.D-14                             ;
  1712.    ANT2  = 'EXCO' 'H' zozo 'SCAL'                                    ;
  1713.    'DETRUIT' zozo                                                    ;
  1714. *  ANT3 caracterise la précipitation quand vaut 1
  1715. *  On autorise une seuil supérieur legerement car on relache du
  1716. *  precipité en solution dans le point fixe
  1717.    zozo = 'MASQ' CHA2 'SUPERIEUR'  ((1.D00 '+' 1.D-7) * LIMSOL)      ;
  1718.    ANT3  = 'EXCO' 'H' zozo 'SCAL'                                    ;
  1719.    'DETRUIT' zozo                                                    ;
  1720.  
  1721.  
  1722. *   'MESSAGE' 'maxi cha2' ('MAXIMUM' cha2);
  1723. *   'MESSAGE' 'maxi ant2' ('MAXIMUM' ant2);
  1724. *   'MESSAGE' 'maxi ant3' ('MAXIMUM' ant3);
  1725.  
  1726. *     'SI' (TPS 'EGA' 8.0D5 1.) ;
  1727. *     DEPLA (DOMA moddarcy 'MAILLAGE')  AFFI (1./15.) (0. 0.) (1. 0.) ;
  1728. *     toto = kcha MODDARCY pre2 'CHAM';
  1729. *     trac MODDARCY toto;
  1730. *     toto =  kcha MODDARCY ANT2 'CHAM';
  1731. **     trac MODDARCY toto;
  1732. *    tutu =  kcha MODDARCY ANT3 'CHAM';
  1733. *     trac MODDARCY (toto '-' tutu);
  1734. *     DEPLA (DOMA moddarcy 'MAILLAGE')  AFFI (15.) (0. 0.) (1. 0.) ;
  1735. *     'FINSI'  ;
  1736.  
  1737.  
  1738. *  le nouveau ANT2 caracterise la précipitation dissolution
  1739.  
  1740.    zozo = ANT2 '+' ANT3 ;
  1741.    'DETRUIT' ANT3   ;
  1742.    ANT4 =  'MASQ' zozo 'SUPERIEUR' 0.1                    ;
  1743.    'DETRUIT' zozo    ;
  1744.  
  1745.    zozo = ANT4 '-' ANT0 ;
  1746.    CRIT  = 'ABS' zozo                                       ;
  1747.    'DETRUIT' zozo ;
  1748.  
  1749.    'DETRUIT' ANT0 ;
  1750.    ANT0  = 'COPIER' ANT4                                             ;
  1751.    'DETRUIT' ANT4                                                    ;
  1752.  
  1753.    MAXCRIT = 'MAXIMUM' CRIT                                          ;
  1754.    'DETRUIT' CRIT                                                    ;
  1755.  
  1756.      'SI' (MAXCRIT < 0.1)                                            ;
  1757. *         On a convergé dans le point fixe.
  1758. *         On stocke la nouvelle distribution du précipité, ANT0
  1759. *         est déjà à jour (prec et dissolution)
  1760.           'DETRUIT'  ANTP                                            ;
  1761.           ANTP = 'COPIER' ANT2                                       ;
  1762.           'DETRUIT'  ANT2                                            ;
  1763.           'DETRUIT'  PRCIPF                                          ;
  1764.           'DETRUIT'  SRCPRCI                                         ;
  1765.           'QUITTER' BOU8                                             ;
  1766.      'SINON'                                                         ;
  1767.  
  1768. *       On résout par point fixe
  1769.  
  1770. *       La où precipitation apparait, on ne fait rien de
  1771. *       special. C'est pris en compte dans le nouveau ant0
  1772.  
  1773. *       La 'OU' le precipité a diparu, le fait qu'il devienne
  1774. *       négatif est génant. ON dissout en fait le précipité
  1775. *       entierement au début du nouvel itéré de point fixe.
  1776. *       Dans la mesure où il est censé disparaître, on espere
  1777. *       que la concentration finale sera en dessous de la
  1778. *       saturation. Sinon plusieurs itérés auront lieu avant
  1779. *       convergence. Par contre on ne veut SURTOUT pas faire
  1780. *       précipiter cet exces temporaire de concentration. DOnc
  1781. *       on met valeur de ant0 indiquant qu'il n'a 'PLUS'
  1782. *       de précipité ('ET' pas de précipitation 'NON' 'PLUS')
  1783.  
  1784. *       1 là 'OU' le précipité vient de diparaitre par
  1785. *       rapport à l'état initial du point fixe !!!!
  1786.         ANTDMM = ANTP '-' ANT2                                       ;
  1787. *       on teste les valeurs positives correspondant à une disparition
  1788. *       effective du precipité (le cas écheant c'est une apparition)
  1789.         ANTDUM = 'MASQUE' ANTDMM SUPERIEUR 0.1                       ;
  1790.         'DETRUIT' ANTDMM                                             ;
  1791. *       0 la 'OU' le precipité vient de disparaitre, 1 ailleurs
  1792.         ANTDM1 = 1.D0 '-' ANTDUM                                     ;
  1793. *       On s'assure d'interdire toute précipitation là 'OU' le
  1794. *       précipité à disparu 'ET' l'on relache son état initial
  1795. *       sous forme d'exces de concentration dans le point fixe.
  1796.         zozo = ANTDM1 * ANT0                                         ;
  1797.         'DETRUIT' ANT0                                               ;
  1798.         ANT0 = 'COPIER' zozo                                         ;
  1799.         'DETRUIT' zozo                                               ;
  1800. *       relâchement de l'état initial du précipité par volume
  1801. *       quantité (1 '-' poros) prec * volum '/' DELTAT
  1802. *       on retire la patie explicite de décroissance du précipité
  1803. *       SRCPRCI = coff
  1804. *                 * ((1.D0 '-' POROS)) *  ANTDUM * PRCI  '/' DELTAT   ;
  1805.         coff = (1.D0 '-' (DELTAT * LAMBD0 * (1.D0 '-' BETA)))        ;
  1806.         zozo = 1.D0 '-' POROS                                        ;
  1807.         zizi = zozo * ANTDUM                                         ;
  1808.         'DETRUIT' zozo                                               ;
  1809.         zozo = zizi * PRCI                                           ;
  1810.         'DETRUIT' zizi                                               ;
  1811.         SRCPRCI = (coff '/' DELTAT) * zozo                           ;
  1812.         'DETRUIT' zozo                                               ;
  1813. *       On retire cela du precipité initial
  1814.         'DETRUIT' PRCIPF                                             ;
  1815.         PRCIPF = ANTDM1 * PRCI                                       ;
  1816.  
  1817.         'DETRUIT'  ANTDUM                                            ;
  1818.         'DETRUIT'  ANTDM1                                            ;
  1819.         'DETRUIT'  ANT2                                              ;
  1820.  
  1821. *      Coefficient correspondant à la décroissance implicite :
  1822.        DEN0  = LAMBD0 * BETA * DELTAT + 1.                           ;
  1823. *      Terme implicite de dissolution, schéma Euler
  1824.        DEN1 = CODIS * ANT0 * DELTAT                                  ;
  1825. *      Calcul du coefficient devant (Ct+dt - Ct)/dt
  1826. *      solub en poros*codis Id et non Retard*poros*codis
  1827.        zozo = DEN0 * PORETSU                                         ;
  1828.        zeze = DEN1 * POROS                                           ;
  1829.        'DETRUIT' COFDC                                               ;
  1830.        COFDC = zozo '+' zeze                                         ;
  1831.        'DETRUIT' zozo                                                ;
  1832.        'DETRUIT' zeze                                                ;
  1833. *      On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1834.        GEOL1 . 'POROSITE' = COFDC                                    ;
  1835. *      On indque que le coef devant DC/DT a changé
  1836.        TABMODI . 'POROSITE' = VRAI                                   ;
  1837.        NBIT = NBIT + 1                                               ;
  1838. *
  1839.      'FINSI'                                                         ;
  1840.      'SI' (&BOU8 EGA ITMAXI)                                         ;
  1841.           'MESS' 'Sortie sans convergence'                           ;
  1842.           'QUITTER' BOU8                                             ;
  1843.      'FINSI'                                                         ;
  1844.  
  1845.  
  1846.  
  1847. *-----------
  1848. 'FIN' BOU8 ;
  1849. *-----------
  1850.  
  1851.  
  1852.  
  1853. NB0 = &BOU8 - 2                                                      ;
  1854. 'SI' (NB0 >EG 1)                                                     ;
  1855.     'MESS' '    ' NB0  ' itérations supp de solubilité'              ;
  1856. 'FINSI'                                                              ;
  1857.  
  1858. * 'Fin module Dissolution d ordre 1'
  1859. 'FINSI'                                                              ;
  1860.  
  1861. *|-------------------------------------|
  1862. *|                Retard Non Linéaire  |
  1863. *|-------------------------------------|
  1864. *- Calcul du nouveau coefficient de retard et de l'adsorbât
  1865. 'SI' RNONLIN                                                         ;
  1866. *    1.  F(C)
  1867.      CC0  = 'NOMC' CHRG 'SCAL'                                       ;
  1868.      'SI' RNONLINL                                                   ;
  1869. *         Langmuir
  1870.           FF = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 ))           ;
  1871.      'SINON'                                                         ;
  1872. *         Freundlich
  1873. *         La bidouille utilisant les masques sert à obtenir le champ
  1874. *         des signes de CC0.
  1875.           a = 'MASQUE' CC0 'SUPERIEUR' 0.                            ;
  1876.           b = 'MASQUE' CC0 'INFERIEUR' 0.                            ;
  1877.           FF = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b)            ;
  1878.           'DETRUIT' a                                                ;
  1879.           'DETRUIT' b                                                ;
  1880.      'FINSI'                                                         ;
  1881.      FF0 = FF                                                        ;
  1882. *    2.  F(C + dC)
  1883.      CC1  = ('NOMC' CHA2 'SCAL') + EpsCord                           ;
  1884.      'SI' RNONLINL                                                   ;
  1885. *         Langmuir
  1886.           FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 ))          ;
  1887.      'SINON'                                                         ;
  1888. *         Freundlich
  1889.           a = 'MASQUE' CC1 'SUPERIEUR' 0.                            ;
  1890.           b = 'MASQUE' CC1 'INFERIEUR' 0.                            ;
  1891.           FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b)           ;
  1892.           'DETRUIT' a                                                ;
  1893.           'DETRUIT' b                                                ;
  1894.      'FINSI'                                                         ;
  1895. *    3.  R = 1 + dF/dC
  1896.      RETAR2 = 1.D0 + ( (FF1 - FF0) / (CC1 - CC0) )                   ;
  1897.      'DETRUIT' CC1                                                   ;
  1898.      'DETRUIT' FF1                                                   ;
  1899. *
  1900. *-   Critère de sortie de la boucle de retard non linéaire :
  1901.      CRIT = ( 'ABS' (RETAR2 - RETARC) ) / RETARC                     ;
  1902.      CRIT = 'MAXIMUM' CRIT                                           ;
  1903.      'SI' (CRIT < EPSRNL)                                            ;
  1904.           'QUITTER' BOURNL                                           ;
  1905.      'SINON'                                                         ;
  1906. *         Le coef. de retard ne change que si le critère est mauvais.
  1907. *         Réévaluation des variables affectées devant DC/DT:
  1908.           RETARC   = RETAR2                                          ;
  1909.           PORETSU  = RETARC  * POROS                                 ;
  1910. *         coef devant dc/dt modifié
  1911.           TABMODI . 'POROSITE' = VRAI                                ;
  1912.  
  1913. *         Fonctions de volume réactualisées
  1914.           PORETSU = RETARC * POROS                                   ;
  1915.           DECROI0  = 'EXCO' ((-1.) * LAMBD0 * PORETSU ) 'SCAL' 'SCAL';
  1916.           DECROI0  = VOLU1 * DECROI0                                 ;
  1917.  
  1918.      'FINSI'                                                         ;
  1919. 'SINON'                                                              ;
  1920. *    NE sert a rien GBM
  1921. *    Retard linéaire F = (R - 1) C
  1922. *     zozo = 'NOMC' CHRG 'SCAL'                                       ;
  1923. *     FF  =  RETAR1M1 * zozo                                          ;
  1924. *     'DETRUIT' zozo                                                  ;
  1925. 'FINSI'                                                              ;
  1926.  
  1927. * Si on atteint le nombre maximum d'itérations - en particulier si le
  1928. * retard est linéaire - on quitte la boucle :
  1929. 'SI' (IRNL 'EGA' ITMAXRNL)                                           ;
  1931.      '     sortie boucle de retard non linéaire sans convergence'    ;
  1932.      'QUITTER' BOURNL                                                ;
  1933. 'FINSI'                                                              ;
  1934.  
  1935. *-------------
  1936. 'FIN' BOURNL ;
  1937. *-------------
  1938.  
  1939. 'SI' (IRNL > 1 )                                                     ;
  1940.      'MESS' '    ' IRNL ' itérations de retard non linéaire.'
  1941.             ' Critère = ' crit                                       ;
  1942. 'FINSI'                                                              ;
  1943. NBITRNL = NBITRNL + (IRNL - 1)                                       ;
  1944.  
  1945.  
  1946. *|-----------------------------------|
  1947. *|   Fin des modules de résolution   |
  1948. *|-----------------------------------|
  1949. *
  1950. *- Calcul de la dissolution par unité de volume et par unité de temps :
  1951. *
  1952. 'SI' SOLUA                                                           ;
  1953.      DIS2 = DIS1 / VOLU1 / DELTAT                                    ;
  1954. 'FINSI'                                                              ;
  1955. *
  1956. *- Archivage des resultats
  1957. *
  1958. *  si le prochain temps à sauver tombe entre les deux
  1959. *  incréments de temps, on obtient les valeurs à sauver
  1960. *  par interpolation linéaire puis on les stocke.
  1961. *
  1962.  
  1963. *list ((PLACE) * 4. / 1000000.)  ;
  1964.  
  1965. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' )                              ;
  1966.  
  1967. *    Sauvegarde de tous les temps intermédiaires, s'il y en a :
  1968.      'REPETER' BOU9                                                  ;
  1969.           TEMS = 'EXTR' TPSOR ISOR                                   ;
  1970.  
  1971.           'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '&lt;EG' TPS ) )       ;
  1972. *              on interpole les variables linéairement entre
  1973. *              les pas de temps
  1974.                LAST1   = LAST1 + 1                                   ;
  1975.                DTEM    = ( TPS - TEMS ) / DELTAT                     ;
  1976.                UNMO    = 1.D0 - DTEM                                 ;
  1977.                CHBIS   = 'COLI' CHRG DTEM CHA2 UNMO                  ;
  1978.                FLUBIS  = 'COLI' FLU0 DTEM FLU2 UNMO                  ;
  1979.                FLUBISCO= 'COLI' FLUCO0 DTEM FLUCO2 UNMO              ;
  1980.                'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ')               ;
  1981.                     RETBIS  = 'COLI' RETAR0 DTEM RETARC UNMO         ;
  1982.                'SINON'                                               ;
  1983.                     RETBIS  = (DTEM * RETAR0) + (UNMO * RETARC)      ;
  1984.                'FINSI'                                               ;
  1985.                'SI' SOLUB                                            ;
  1986.                     PRBIS   = 'COLI' PRCI DTEM PRE2 UNMO             ;
  1987.                'FINSI'                                               ;
  1988.                TRANSI . 'TEMPS'            . LAST1 = TEMS            ;
  1989.                TRANSI . NOMINC             . LAST1 = 'NOMC' NOMESP
  1990.                                                             CHBIS    ;
  1991.                'DETRUIT' CHBIS                                       ;
  1992.                TRANSI . 'FLUXDIFF'        . LAST1 = 'NOMC' NOMESP
  1993.                                                             FLUBIS   ;
  1994.                'DETRUIT' FLUBIS                                      ;
  1995.                TRANSI . 'FLUXCONV'           . LAST1 = 'NOMC' NOMESP
  1996.                                                             FLUBISCO ;
  1997.                'DETRUIT' FLUBISCO                                    ;
  1998.                'SI' SOLUB                                            ;
  1999.                     TRANSI . 'PRECIPITE'   . LAST1 = 'NOMC' NOMESP
  2000.                                                             PRBIS    ;
  2001.                'DETRUIT' PRBIS                                       ;
  2002.                'FINSI'                                               ;
  2003.                'SI' SAUVRET                                          ;
  2004.                     TRANSI . 'RETARD'      . LAST1 = 'NOMC' NOMESP
  2005.                                                             RETBIS   ;
  2006.                'DETRUIT' RETBIS                                      ;
  2007.                'FINSI'                                               ;
  2008.                ISOR = ISOR + 1                                       ;
  2009.           'SINON'                                                    ;
  2010.                'QUITTER' BOU9                                        ;
  2011.           'FINSI'                                                    ;
  2012.           'SI' ( ISOR '>' DSOR )                                     ;
  2013.                ISOR = DSOR                                           ;
  2014.                'QUITTER' BOU9                                        ;
  2015.           'FINSI'                                                    ;
  2016.      'FIN' BOU9                                                      ;
  2017. 'SINON'                                                              ;
  2018. *    sinon, on sauvegarde tous les temps :
  2019.      LAST1 = LAST1 + 1                                               ;
  2020.      TRANSI . 'TEMPS'            . LAST1 = TPS                       ;
  2021.      TRANSI . NOMINC             . LAST1 = 'NOMC' NOMESP CHA2        ;
  2022.      TRANSI . 'FLUXDIFF'         . LAST1 = 'NOMC' NOMESP FLU2        ;
  2023.      TRANSI . 'FLUXCONV'         . LAST1 = 'NOMC' NOMESP FLUCO2      ;
  2024.      'SI' SOLUB                                                      ;
  2025.           TRANSI . 'PRECIPITE'   . LAST1 = 'NOMC' NOMESP PRE2        ;
  2026.      'FINSI'                                                         ;
  2027.      'SI' SAUVRET                                                    ;
  2028.           TRANSI . 'RETARD'      . LAST1 = 'NOMC' NOMESP RETARC      ;
  2029.      'FINSI'                                                         ;
  2030. 'FINSI'                                                              ;
  2031.  
  2032.  
  2033. *On sauve l'intégrale de decroissance*concentration et
  2034. *                        decroissance*précipité au cours du temps
  2035. * stockée pour l'instant à tous les temps. Peut s'améliorer.
  2036.  
  2037. * indice pour le stockage des intégrales de concentration
  2038. LASTINT = LASTINT '+' 1;
  2039.  
  2040. BM1 = 1.D0 '-' BETA                                                  ;
  2041. * si l'espece fait l'objet d'une filiation on stocke les intégrales
  2042. * du précipité et de la concentration.
  2043. SI ((EXISTE TRANSI DECROISSANCE) 'ET' (TRANSI . 'PERE'))             ;
  2044. *TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1)
  2045. *                   '+' (DELTAT * PORETSU * TRANSI . 'DECROISSANCE' *
  2046. *                          ((BM1 * ('NOMC' NOMESP CHRG))
  2047. *                       '+' ((BETA) * ('NOMC' NOMESP CHA2))))         ;
  2048.    zozo = (DELTAT * TRANSI . 'DECROISSANCE') *  PORETSU              ;
  2049.    zaza = BM1 * CHRG                                                 ;
  2050.    zuzu = BETA *  CHA2                                               ;
  2051.    zeze = 'KOPS' zaza '+' zuzu                                       ;
  2052.    'DETRUIT' zaza                                                    ;
  2053.    'DETRUIT' zuzu                                                    ;
  2054.    zaza = 'KOPS' zozo * zeze                                         ;
  2055.    'DETRUIT' zozo                                                    ;
  2056.    'DETRUIT' zeze                                                    ;
  2057.    zozo = 'NOMC' NOMESP zaza                                         ;
  2058.    'DETRUIT' zaza                                                    ;
  2059.    TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1)
  2060.                           '+' zozo                                   ;
  2061.    'DETRUIT' zozo                                                    ;
  2062.  
  2063. 'SI' SOLUB                                                           ;
  2064.  
  2065. *   TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . LASTINT
  2066. *          '+' (DELTAT * (1.D0 '-' POROS) * TRANSI . 'DECROISSANCE' *
  2067. *                            ((BM1 * ('NOMC' NOMESP PRCI))
  2068. *                       '+' ((BETA) * ('NOMC' NOMESP PRE2))))         ;
  2069.    zizi = (1.D0 '-' POROS)                                           ;
  2070.    zozo = (DELTAT * TRANSI . 'DECROISSANCE') *  zizi                 ;
  2071.    'DETRUIT' zizi                                                    ;
  2072.    zaza = BM1 * PRCI                                                 ;
  2073.    zuzu = BETA *  PRE2                                               ;
  2074.    zeze = 'KOPS' zaza '+' zuzu                                       ;
  2075.    'DETRUIT' zaza                                                    ;
  2076.    'DETRUIT' zuzu                                                    ;
  2077.    zaza = 'KOPS' zozo * zeze                                         ;
  2078.    'DETRUIT' zozo                                                    ;
  2079.    'DETRUIT' zeze                                                    ;
  2080.    zozo = 'NOMC' NOMESP zaza                                         ;
  2081.    'DETRUIT' zaza                                                    ;
  2082.    zaza = TRANSI . 'INTFILI' . LASTINT '+' zozo                      ;
  2083.    'DETRUIT' TRANSI . 'INTFILI' . LASTINT                            ;
  2084.    'DETRUIT' zozo                                                    ;
  2085.    TRANSI . 'INTFILI' . LASTINT = zaza                               ;
  2086.    'OUBLIER' zaza                                                    ;
  2087.  
  2088. 'FINSI'                                                              ;
  2089. FINSI                                                                ;
  2090.  
  2091.  
  2092. * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait.
  2093. * Le dernier pas de temps ne peut pas faire l'objet d'interpolation
  2094. * linéaire si le temps sauvegardé ne tombe pas pile dessus.
  2095. 'SI' (IPAS 'EGA' DCAL)                                               ;
  2096.      TPSDER = TRANSI. 'TEMPS' . LAST1                                ;
  2097.      'SI' ( TPSDER 'NEG' TPSFIN EPSDT )                              ;
  2098.           LAST1   = LAST1 + 1                                        ;
  2099.           TRANSI . 'TEMPS'            . LAST1 = TPS                  ;
  2100.           TRANSI . NOMINC             . LAST1 = 'NOMC' NOMESP CHA2   ;
  2101.           TRANSI . 'FLUXDIFF'         . LAST1 = 'NOMC' NOMESP FLU2   ;
  2102.           TRANSI . 'FLUXCONV'         . LAST1 = 'NOMC' NOMESP FLUCO2 ;
  2103.           'SI' SOLUB                                                 ;
  2104.                TRANSI . 'PRECIPITE'   . LAST1 = 'NOMC' NOMESP PRE2   ;
  2105.           'FINSI'                                                    ;
  2106.           'SI' SAUVRET                                               ;
  2107.                TRANSI . 'RETARD'      . LAST1 = 'NOMC' NOMESP RETARC ;
  2108.           'FINSI'                                                    ;
  2109.      'FINSI'                                                         ;
  2110. 'FINSI'                                                              ;
  2111. *
  2112. *- Initialisations pour le pas suivant
  2113. *
  2114. PRECED = TPS                                                         ;
  2115. 'DETRUIT' CHRG                                                       ;
  2116. CHRG   = CHA2                                                        ;
  2117. 'DETRUIT' FLU0                                                       ;
  2118. FLU0   = FLU2                                                        ;
  2119. 'DETRUIT' FLUCO0                                                     ;
  2120. FLUCO0 = FLUCO2                                                      ;
  2121. 'OUBLIER' GEOL1                                                      ;
  2122. 'OUBLIER' GEOL2                                                      ;
  2123. GEOL1  = 'TABLE' GEOLPF1                                             ;
  2124. GEOL2  = 'TABLE' GEOLPF2                                             ;
  2125. 'DETRUIT' RETAR0                                                     ;
  2126. RETAR0 = RETARC                                                      ;
  2127. DELOLD = DELTAT                                                      ;
  2128. 'SI' SOLUB                                                           ;
  2129.      'DETRUIT' PRCI                                                  ;
  2130.      PRCI = PRE2                                                     ;
  2131. 'FINSI'                                                              ;
  2132. 'SI' ('EXISTE' TRANSI 'SOURCE')                                      ;
  2133.    'DETRUIT' TERSC2M1                                                ;
  2134.    TERSC2M1 = TERSC2                                                 ;
  2135. 'FINSI'                                                              ;
  2136. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE')                                 ;
  2137.    'DETRUIT' FLUIMPM1                                                ;
  2138.    FLUIMPM1 = FLUIMPO                                                ;
  2139. 'FINSI'                                                              ;
  2140. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP')                                 ;
  2141.    'DETRUIT' FLUTMPM1                                                ;
  2142.    FLUTMPM1 = FLUTMPO                                                ;
  2143. 'FINSI'                                                              ;
  2144. 'SI' ('EXISTE' TRANSI 'MIXTES')                                    ;
  2145.    'DETRUIT' FLUMMPM1                                                ;
  2146.    FLUMMPM1 = FLUMMPO                                                ;
  2147. 'FINSI'                                                              ;
  2148.  
  2149. *MESS 'TEMPS2' ; TEMPS;
  2150. *'MENAGE'                                                             ;
  2151.  
  2152. * destruction d'objets temporaires
  2153.  
  2154. 'DETRUIT' TERSC3 ;
  2155.  
  2156. *list ((PLACE) * 4. / 1000000.)  ;
  2157.  
  2158. *temps impr;
  2159.  
  2160. *
  2161. *=============
  2162. 'FIN' BOUTPS ;    COMM 'Boucle sur le temps'                         ;
  2163. *=============
  2164. *
  2165. *- Indication du nombre de passages :
  2166. *
  2167. 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION'))                        ;
  2168.   'MESS' 'Nombre total de matrices calculées : ' NBIT                ;
  2169. 'FINSI'                                                              ;
  2170. 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION'))                      ;
  2171.   'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  2172.           NBIT ;
  2173. 'FINSI'                                                              ;
  2174. 'SI' SOLUP                                                           ;
  2175.   'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  2176.           NBIT;
  2177. 'FINSI'                                                              ;
  2178. 'SI' RNONLIN                                                         ;
  2179.   'MESS' 'Nombre total d itérations supplémentaires de '
  2180.          'retard non linéaire : ' NBITRNL;
  2181. 'FINSI'                                                              ;
  2182.  
  2183. 'FINP'                                                               ;
  2184.  

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