Télécharger la plaquette Cast3M


fr en

Télécharger kfm12.eso

Retour à la liste

Numérotation des lignes : 

kfm12

  1. C KFM12     SOURCE    CB215821  20/11/25    13:31:14     10792          
  2.       SUBROUTINE KFM12(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
  3.      &     ICHPSU,ICHPDI,ICHPVO,INORM,
  4.      &     MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL,
  5.      &     MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
  6.      &     IDU,IPROBL)
  7. C************************************************************************
  8. C
  9. C PROJET            :  CASTEM 2000
  10. C
  11. C NOM               :  KFM12
  12. C
  13. C DESCRIPTION       :  Voir aussi KFM1
  14. C
  15. C                      Cas deux dimensions, gaz "calorically perfect"
  16. C                      Méthode de relaxation de type SGS
  17. C
  18. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  19. C
  20. C AUTEUR            :  T. KLOCZKO DEN/DM2S/SFME/LTMF
  21. C
  22. C************************************************************************
  23. C ENTREES
  24. C
  25. C
  26. C  1) Pointeurs des CHPOINTs
  27. C
  28. C     IRES    : CHPOINT 'CENTRE' contenant le residu
  29. C
  30. C     IRN     : CHPOINT 'CENTRE' contenant la masse volumique
  31. C
  32. C     IGN     : CHPOINT 'CENTRE' contenant la q.d.m.
  33. C
  34. C     IRETN   : CHPOINT 'CENTRE' contenant l'energie totale
  35. C
  36. C     IGAMN   : CHPOINT 'CENTRE' contenant le gamma
  37. C
  38. C     IVCO    : CHPOINT 'CENTRE' contenant le cut-off
  39. C
  40. C     ICHLIM  : CHPOINT conditions aux bords
  41. C
  42. C     ICOEV   : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho)
  43. C               coeff pour calculer le ray. spect. de la matrice
  44. C               visc.
  45. C
  46. C  3) Pointeurs de CHPOINTs de la table DOMAINE
  47. C
  48. C     ICHPSU  : CHPOINT "FACE" contenant la surface des faces
  49. C
  50. C     ICHPDI  : CHPOINT "CENTRE" contenant le diametre minimum
  51. C               de chaque element
  52. C
  53. C     ICHPVO  : CHPOINT "CENTRE" contenant le volume
  54. C               de chaque element
  55. C
  56. C     INORM   : CHPOINT "CENTRE" contenant le normales aux faces
  57. C
  58. C
  59. C  4) Pointeurs de MELEME de la table DOMAINE
  60. C
  61. C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
  62. C
  63. C     MELEMF  : MELEME 'FACE' du SPG des FACES
  64. C
  65. C     MELEFL  : MELEME 'FACEL' du connectivité FACES -> ELEM
  66. C
  67. C     MELLIM  : MELEME SPG conditions aux bords
  68. C
  69. C  5)
  70. C
  71. C     DCFL = le double de la CFL
  72. C
  73. C     DTPS = pas de temps physique
  74. C
  75. C     NJAC = nombre d'iterations dans le jacobien
  76. C
  77. C
  78. C SORTIES
  79. C
  80. C     IDU  : resultat
  81. C
  82. C     IPROBL : si > 0, il y a un probleme
  83. C
  84. C************************************************************************
  85. C
  86. C HISTORIQUE (Anomalies et modifications éventuelles)
  87. C
  88. C HISTORIQUE : 2004  : creation
  89. C
  90. C
  91. C************************************************************************
  92. C
  93. C
  94. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  95. C                                         GAMMA \in (1,3)
  96. C                                         Y \in (0,1)
  97. C       Si non il faut le faire!!!
  98. C
  99. C************************************************************************
  100. C
  101. C
  102. C**** Variables de COOPTIO
  103. C
  104. C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  105. C     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
  106. C     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
  107. C     &        ,IECHO,   IIMPI,  IOSPI
  108. C     &        ,IDIM
  109. C     &        ,MCOORD
  110. C     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  111. C     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  112. C     &        ,NORINC,NORVAL,NORIND,NORVAD
  113. C     &        ,NUCROU, IPSAUV, IFICLE, IPREFI, IREFOR, ISAFOR
  114. C
  115.       IMPLICIT INTEGER(I-N)
  116.       INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM
  117.      &     ,MELEMF,IDU,NFAC, NFACE
  118.      &     ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCEG, NLCED
  119.      &     ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP
  120.      &     ,I1,I2, IBLJAC, IPROBL, IVCO, MELTFA
  121.      &     ,NSOUPO,IFACT,NELT, ISOUP, IELT, IFAC, IFACG, JG
  122.      &     ,ISOUP1, ISOUP2, ICOEF,ICOEV
  123.      &     , ID, IELT1, IELT2, IFAC1, IFAC2
  124.       REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, RECG, RECD, PG, VOLG
  125.      &     , DIAMG, ROD, RUXD, RUYD, UND, RETD, GAMD, PD
  126.      &     , CNX, CNY, DCFL, DTPS
  127.      &     , SURF, SIGMAF
  128.      &     , UNSDT, UNSDTF, UXD, UYD
  129.      &     , FR, FRUX, FRUY, FRET, COEF
  130.      &     , CTX, CTY, RHTD, UTD, UTG, UXG, UYG, QG, CG, UCOG
  131.      &     , ROF, PF, UNF,  GAMF, UCOF, CF, PHI2, COEF1
  132.      &     , QN
  133. C
  134.       CHARACTER*8 TYPE
  135.       LOGICAL LOGINV, LOGJAC
  136. C
  137. C**** LES INCLUDES
  138. C
  139.  
  140. -INC PPARAM
  141. -INC CCOPTIO
  142. -INC SMCHPOI
  143.       POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
  144.      &       , MPODU.MPOVAL, MPOVOL.MPOVAL
  145.      &       , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
  146.      &       , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL
  147.      &       , MPOCO1.MPOVAL, MELEFL.MELEME, MELEMC.MELEME
  148.      &       , MPCOEV.MPOVAL
  149. -INC SMCOORD
  150. -INC SMELEME
  151. -INC SMLMOTS
  152. -INC SMLENTI
  153.       POINTEUR MLELIM.MLENTI, MLESUP.MLENTI
  154. C
  155. C**** POINTEUR vecteurs
  156. C
  157.       SEGMENT VEC
  158.       REAL*8 VALVEC(I1)
  159.       ENDSEGMENT
  160.       POINTEUR SIGMA0.VEC, AN0.VEC, DN0.VEC, D1.VEC
  161.      &     , USDTI.VEC, PHI2V.VEC, SIGMAV.VEC
  162. C
  163. C**** POINTEUR matriciels
  164. C
  165.       SEGMENT MAT
  166.       REAL*8 VAL(I1,I2)
  167.       ENDSEGMENT
  168.       POINTEUR CONS0.MAT, CONS1.MAT, BN0.MAT, BN1.MAT
  169.      &     , CN0.MAT, CN1.MAT, DU.MAT
  170.      &     , RES.MAT
  171.      &     , Q1.MAT, Q2.MAT
  172. C
  173. C
  174. C**** IPROBL: si different de 0, un probleme a été detecté
  175. C
  176.       IPROBL=0
  177. C
  178. C**** Initialisation des MLENTI des conditions aux limites
  179. C
  180. C
  181.       CALL KRIPAD(MELLIM,MLELIM)
  182. C     SEGINI MLELIM
  183. C
  184. C**** Initialisation des MELEMEs
  185. C
  186. C     'CENTRE', 'FACEL'
  187. C
  188.       SEGACT MELEFL
  189.       NFAC = MELEFL.NUM(/2)
  190. C
  191. C**** KRIPAD pour la correspondance global/local de centre
  192. C
  193.       CALL KRIPAD(MELEMC,MLENT1)
  194. C
  195. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  196. C
  197. C     Si i est le numero global d'un noeud de ICEN,
  198. C     MLENT1.LECT(i) contient sa position, i.e.
  199. C
  200. C     I              = numero global du noeud centre
  201. C     MLENT1.LECT(i) = numero local du noeud centre
  202. C
  203. C     MLENT1 déjà activé, i.e.
  204. C
  205. C     SEGINI MLENT1
  206. C
  207. C
  208. C**** KRIPAD pour la correspondance global/local de 'FACE'
  209. C
  210.       CALL KRIPAD(MELEMF,MLENT2)
  211. C     SEGINI MLENT2
  212. C
  213. C
  214. C**** CHPOINTs de la table DOMAINE
  215. C
  216.       CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  217.       CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  218.       CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
  219.       CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
  220. C
  221. C**** LICHT active les MPOVALs en *MOD
  222. C
  223. C     i.e.
  224. C
  225. C     SEGACT MPOVSU*MOD
  226. C     SEGACT MPOVOL*MOD
  227. C     SEGACT MPOVDI*MOD
  228. C     SEGACT MPNORM*MOD
  229. C
  230.       CALL LICHT(IDU,MPODU,TYPE,IGEOMC)
  231. C     SEGACT MPODU*MOD
  232. C
  233. C     USDTI initialisé a zero; USDTI = 1 / DT
  234. C
  235.       NCEN=MPOVOL.VPOCHA(/1)
  236.       I1=NCEN
  237.       SEGINI USDTI
  238. C
  239.       CALL LICHT(IRES,MPORES,TYPE,IGEOMC)
  240.       CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
  241.       CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
  242.       CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
  243.       CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
  244.       CALL LICHT(IVCO,MPOCO,TYPE,IGEOMC)
  245.       CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
  246.       CALL LICHT(ICOEV,MPCOEV,TYPE,IGEOMC)
  247. C
  248. C     SEGACT MPORES*MOD
  249. C     SEGACT MPRN*MOD
  250. C     SEGACT MPGN*MOD
  251. C     SEGACT MPRETN*MOD
  252. C     SEGACT MPGAMN*MOD
  253. C     SEGACT MPOCO
  254. C     SEGACT MPLIM*MOD
  255. C     SEGACT MPCOEV*MOD
  256. C
  257.       SEGINI, MPOCO1=MPOCO
  258. C
  259. C**** MPOCO1 is the "corrected cut-off" speed.
  260. C     It is bigger than the given cut-off and the local
  261. C     speed.
  262. C
  263.       IF(IGEOMF .NE. MELEMF)THEN
  264.          WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  265.          WRITE(IOIMP,*) 'Probleme de SPG'
  266. C        21 2
  267. C        Données incompatibles
  268.          CALL ERREUR(21)
  269.          GOTO 9999
  270.       ENDIF
  271.       IF(IGEOMC .NE. MELEMC)THEN
  272.          WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  273.          WRITE(IOIMP,*) 'Probleme de SPG'
  274. C        21 2
  275. C        Données incompatibles
  276.          CALL ERREUR(21)
  277.          GOTO 9999
  278.       ENDIF
  279. C
  280. C**** On initialise GAMMA=I+(1-phi2)/phi2*Q1.(Q2)^t
  281. C
  282.       I1=4
  283.       I2=NCEN
  284.       SEGINI Q1
  285.       SEGINI Q2
  286. C
  287. C**** On initialise le segment de travail D1
  288. C
  289.       I1=4
  290.       SEGINI D1
  291. C
  292. C**** On initialise AN0, DN0, PHI2V
  293. C
  294.       I1=NCEN
  295.       SEGINI AN0
  296.       SEGINI DN0
  297.       SEGINI PHI2V
  298. C
  299. C**** On initialise BN0
  300. C
  301.       I1=4
  302.       I2=NCEN
  303.       SEGINI BN0
  304. C
  305. C**** On initialise BN1
  306. C
  307.       I1=4
  308.       I2=NCEN
  309.       SEGINI BN1
  310. C
  311. C**** On initialise CN0
  312. C
  313.       I1=4
  314.       I2=NCEN
  315.       SEGINI CN0
  316. C
  317. C**** On initialise CN1
  318. C
  319.       I1=4
  320.       I2=NCEN
  321.       SEGINI CN1
  322. C
  323. C**** On initialise SIGMA0, SIGMAV
  324. C
  325.       I1=NFAC
  326.       SEGINI SIGMA0
  327.       SEGINI SIGMAV
  328. C
  329. C**** RES
  330. C
  331.       I1=4
  332.       I2=NCEN
  333.       SEGINI RES
  334.       DO ICEN=1,NCEN,1
  335.          DO ICOMP=1,4,1
  336.             RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP)
  337.          ENDDO
  338.       ENDDO
  339.       SEGDES MPORES
  340. C
  341. C***** DU (increment des variables conservatives)
  342. C
  343.       I1=4
  344.       I2=NCEN
  345.       SEGINI DU
  346. C
  347. C**** Les variables conservatives
  348. C
  349. C     CONS1.VAL contient les variables conservatives + gamma
  350. C     i.e. rho(i) = CONS.VAL(1,i)
  351. C          rhoux(i) = CONS.VAL(2,i)
  352. C          rhouy(i) = CONS.VAL(3,i)
  353. C          rhoet(i) = CONS.VAL(4,i)
  354. C          gam(i)   = CONS.VAL(5,i)
  355. C
  356.       I1=5
  357.       I2=NCEN
  358.       SEGINI CONS0
  359.       DO ICEN=1,NCEN,1
  360.          CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
  361.          CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
  362.          CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
  363.          CONS0.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1)
  364.          CONS0.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1)
  365.       ENDDO
  366.       SEGDES MPRN
  367.       SEGDES MPRETN
  368.       SEGDES MPGN
  369.       SEGDES MPGAMN
  370. C
  371.       SEGINI, CONS1=CONS0
  372. C
  373. C**** Storage of l-1 values of F at the l Jacobi iteration
  374. C
  375.       IPT1=MELTFA
  376.       SEGACT IPT1
  377.       NSOUPO= IPT1.LISOUS(/1)
  378.       NSOUPO=MAX(NSOUPO,1)
  379.       JG=NSOUPO
  380.       SEGINI MLESUP
  381.       IFACT=0
  382.       IF(NSOUPO .EQ. 1)THEN
  383.          MLESUP.LECT(1)=IPT1
  384.          NFACE=IPT1.NUM(/1)
  385.          NELT=IPT1.NUM(/2)
  386.          IFACT=IFACT+(NFACE*NELT)
  387.       ELSE
  388.          DO ISOUP = 1, NSOUPO, 1
  389.             IPT2=IPT1.LISOUS(ISOUP)
  390.             MLESUP.LECT(ISOUP)=IPT2
  391.             SEGACT IPT2
  392.             NFACE=IPT2.NUM(/1)
  393.             NELT=IPT2.NUM(/2)
  394.             IFACT=IFACT+(NFACE*NELT)
  395.          ENDDO
  396.       ENDIF
  397. C
  398. C***********************************************************************
  399. C**** FIRST JACOBI ITERATION *******************************************
  400. C***********************************************************************
  401. C
  402. C
  403. C**** LOOP on ELTFA to compute:
  404. C
  405. C     AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j
  406. C     SIGMA0_j
  407. C     The cut-off in each cell
  408. C     The reciprocal of the dual time step USDTI
  409. C
  410.       ICEN=0
  411.       IFACG=0
  412.       DO ISOUP=1,NSOUPO,1
  413.          IPT2=MLESUP.LECT(ISOUP)
  414.          NFACE=IPT2.NUM(/1)
  415.          NELT=IPT2.NUM(/2)
  416.          DO IELT=1,NELT,1
  417.             ICEN=ICEN+1
  418.             DO IFAC=1,NFACE,1
  419. C              IFACG increase with IFAC,IELT,ISOUP
  420.                IFACG = IFACG+1
  421.                NGCF  = IPT2.NUM(IFAC,IELT)
  422.                NLCF  = MLENT2.LECT(NGCF)
  423. C
  424. C************* NLCF  = numero local du centre de facel
  425. C              NGCF  = numero global du centre de facel
  426. C              NGCEG = numero global du centre ELT "gauche"
  427. C              NLCEG = numero local du centre ELT "gauche"
  428. C              NGCED = numero global du centre ELT "droite"
  429. C              NLCED = numero local du centre ELT "droite"
  430. C
  431.                NGCEG = MELEFL.NUM(1,NLCF)
  432.                NGCED = MELEFL.NUM(3,NLCF)
  433.                NLCEG = MLENT1.LECT(NGCEG)
  434.                NLCED = MLENT1.LECT(NGCED)
  435. C
  436.                CNX = MPNORM.VPOCHA(NLCF,1)
  437.                CNY = MPNORM.VPOCHA(NLCF,2)
  438.                CTX = -1.0D0 * CNY
  439.                CTY = CNX
  440. C
  441. C************* Left or right?
  442. C
  443.                IF(NLCEG .EQ. ICEN)THEN
  444.                   LOGINV=.FALSE.
  445.                ELSEIF(NLCED .EQ. ICEN)THEN
  446.                   LOGINV=.TRUE.
  447.                ELSE
  448.                   WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
  449.                   WRITE(IOIMP,*) 'Probleme de SPG'
  450. C                 21 2
  451. C                 Données incompatibles
  452.                   CALL ERREUR(21)
  453.                   GOTO 9999
  454.                ENDIF
  455. C************* If right, inversion
  456.                IF(LOGINV) THEN
  457.                   NLCED=NLCEG
  458.                   NLCEG=ICEN
  459.                   CNX = -1.0D0*CNX
  460.                   CNY = -1.0D0*CNY
  461.                   CTX = -1.0D0*CTX
  462.                   CTY = -1.0D0*CTY
  463.                ENDIF
  464. C
  465.                SURF = MPOVSU.VPOCHA(NLCF,1)
  466. C
  467. C************* Récupération des Etats "gauche" et "droite"
  468. C
  469. C
  470.                ROG  = CONS0.VAL(1,NLCEG)
  471.                RUXG = CONS0.VAL(2,NLCEG)
  472.                RUYG = CONS0.VAL(3,NLCEG)
  473.                UXG  = RUXG / ROG
  474.                UYG  = RUYG / ROG
  475.                UNG  = (UXG * CNX) + (UYG * CNY)
  476.                UTG  = (UXG * CTX) + (UYG * CTY)
  477.                RETG = CONS0.VAL(4,NLCEG)
  478.                GAMG = CONS0.VAL(5,NLCEG)
  479.                RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
  480.                RECG = RECG / ROG
  481.                PG   = (GAMG - 1.0D0) * (RETG - RECG)
  482.                VOLG = MPOVOL.VPOCHA(NLCEG,1)
  483.                DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  484. C
  485.                ROD  = CONS0.VAL(1,NLCED)
  486.                RUXD = CONS0.VAL(2,NLCED)
  487.                RUYD = CONS0.VAL(3,NLCED)
  488.                UXD  = RUXD / ROD
  489.                UYD  = RUYD / ROD
  490.                RETD = CONS0.VAL(4,NLCED)
  491.                GAMD = CONS0.VAL(5,NLCED)
  492.                RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
  493.                RECD = RECD / ROD
  494.                PD   = (GAMD - 1.0D0) * (RETD - RECD)
  495.                UND  = (UXD * CNX) + (UYD * CNY)
  496.                UTD  = (UXD * CTX) + (UYD * CTY)
  497. C
  498.                IF(NLCEG .EQ. NLCED)THEN
  499. C           Murs ou condition limite
  500.                   NLLIM=MLELIM.LECT(NGCF)
  501.                   IF(NLLIM .EQ.0)THEN
  502. C              Mur
  503.                      UND  = -1.0D0 * UNG
  504.                      UTD  = UTG
  505.                      UXD  = UND * CNX + UTD * CTX
  506.                      UYD  = UND * CNY + UTD * CTY
  507.                      RUXD = ROD*UXD
  508.                      RUYD = ROD*UYD
  509.                   ELSE
  510.                      ROD  = MPLIM.VPOCHA(NLLIM,1)
  511.                      UXD  = MPLIM.VPOCHA(NLLIM,2)
  512.                      UYD  = MPLIM.VPOCHA(NLLIM,3)
  513.                      RUXD = ROD*UXD
  514.                      RUYD = ROD*UYD
  515.                      PD   = MPLIM.VPOCHA(NLLIM,4)
  516.                      GAMD = GAMG
  517.                      UND  = (UXD * CNX) + (UYD * CNY)
  518.                      UTD  = (UXD * CTX) + (UYD * CTY)
  519.                      RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  520.      &                    +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
  521.                   ENDIF
  522.                ENDIF
  523. C
  524. C************* We check that density and pressure are positive
  525. C
  526.                IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  527.      &              (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  528.                   WRITE(IOIMP,*) 'FMM'
  529.                   WRITE(IOIMP,*) 'TEST'
  530. C                 22 0
  531. C**************** Opération malvenue. Résultat douteux
  532. C
  533.                   CALL ERREUR(22)
  534.                   IPROBL = 1
  535.                   GOTO 9998
  536.                ENDIF
  537. C
  538. C************* We compute the interfacial state
  539. C
  540.                ROF  = 0.5D0*(ROG+ROD)
  541.                PF   = 0.5D0*(PG+PD)
  542.                UNF  = 0.5D0*(UNG+UND)
  543.                GAMF = 0.5D0*(GAMG+GAMD)
  544. C
  545. C************* We compute SIGMA and the corrected cut-off speed
  546. C
  547.                UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1))
  548.                UCOF=MAX(UCOF,(((UNG*UNG)+(UTG*UTG))**0.5D0))
  549.                UCOF=MAX(UCOF,(((UND*UND)+(UTD*UTD))**0.5D0))
  550.                IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1))
  551.      &            MPOCO1.VPOCHA(NLCEG,1)= UCOF
  552. C
  553.                QN=ABS(UNF)
  554.                CF=(GAMF*PF/ROF)**0.5D0
  555.                IF(UCOF .GE. CF)THEN
  556.                   PHI2 = 1.0D0
  557.                ELSEIF(QN .GT. CF)THEN
  558.                   PHI2 = 1.0D0
  559.                ELSEIF(QN .LT. UCOF)THEN
  560.                   PHI2 = UCOF / CF
  561.                ELSE
  562.                   PHI2 = QN / CF
  563.                ENDIF
  564.                PHI2=PHI2*PHI2
  565. C
  566.                SIGMAF=(1.0D0-PHI2)*QN
  567.                SIGMAF=SIGMAF*SIGMAF
  568.                SIGMAF=SIGMAF+(4.0D0*PHI2*CF*CF)
  569.                SIGMAF=SIGMAF**0.5D0
  570.                SIGMAF=SIGMAF+((1.0D0+PHI2)*QN)
  571.                SIGMAF=0.5D0*SIGMAF
  572. C
  573.                SIGMA0.VALVEC(NLCF)=SIGMAF
  574. C
  575. C************* We compute AN0
  576. C
  577.                AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG)
  578.      &                          +(0.5D0*SIGMAF*SURF/VOLG)
  579. C
  580. C************* We compute the dual time step
  581. C
  582.                UNSDT=USDTI.VALVEC(NLCEG)
  583.                UNSDTF = SIGMAF / DIAMG
  584.                IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF
  585. C
  586. C************* We compute SIGMAV
  587. C
  588. C              COEF=FG
  589.                COEF=(MCOORD.XCOOR(((NGCEG-1)*3)+1)-
  590.      &                 MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
  591.                COEF=COEF+
  592.      &              ((MCOORD.XCOOR(((NGCEG-1)*3)+2)-
  593.      &              MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
  594. C              COEF1=FD
  595.                COEF1=(MCOORD.XCOOR(((NGCED-1)*3)+1)-
  596.      &              MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
  597.                COEF1=COEF1+
  598.      &              ((MCOORD.XCOOR(((NGCED-1)*3)+2)-
  599.      &              MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
  600. C              COEF=|FG|+|FD|
  601.                COEF=ABS(COEF)+ABS(COEF1)
  602.                SIGMAF=0.5D0*(MPCOEV.VPOCHA(NLCEG,1)+
  603.      &              MPCOEV.VPOCHA(NLCED,1))/COEF
  604.                SIGMAV.VALVEC(NLCF)=SIGMAF
  605. C
  606. C************* We compute the contribution of SIGMAV to
  607. C              the diagonal part of the matrix to invert
  608. C
  609. C              In the given report, DNO is a_i^0, i.e.
  610. C              DN0 = AN0 (previously defined) +
  611. C                    + (\frac{1}{{\Omega_i}}
  612. C                       \sum_j \sigma_{V,i,j}^{0} S_j)
  613. C                    + 1 / \Delta t^n
  614. C
  615.                COEF=SURF/VOLG
  616.                DN0.VALVEC(NLCEG)=DN0.VALVEC(NLCEG)+COEF*SIGMAF
  617. C
  618. C
  619. C************* We compute the flux
  620. C              FLU.VALVEC(1:4,IFAC) flux in IFAC face
  621. C              1 -> mass
  622. C              2 -> momentum (x)
  623. C              3 -> momentum (y)
  624. C              4 -> energy
  625. C
  626.                RHTD = RETD+PD
  627.                FR   = (ROD*UND)
  628.                FRUX = (ROD*UND*UXD)+(PD*CNX)
  629.                FRUY = (ROD*UND*UYD)+(PD*CNY)
  630.                FRET = (UND*RHTD)
  631. C
  632.                COEF = SURF/(2.0D0*VOLG)
  633. C
  634.                BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  635.      &                          + FR   * COEF
  636.                BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  637.      &                          + FRUX * COEF
  638.                BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  639.      &                          + FRUY * COEF
  640.                BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  641.      &                          + FRET * COEF
  642. C
  643. C************* Viscous contribution
  644. C
  645.                COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
  646. C
  647.                BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  648.      &                          - (COEF * ROD)
  649.                BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  650.      &                          - (COEF * RUXD)
  651.                BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  652.      &                          - (COEF *  RUYD)
  653.                BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  654.      &                          - (COEF *  RETD)
  655. C
  656. C************* On calcule CN0
  657. C
  658. C              CN0 in the referred report is given by
  659. C              \sum_j
  660. C              \left\{
  661. C              \frac{S_j}{2 \Omega_i} ~
  662. C              \left(
  663. C              \sigma_{i,j}^{0}
  664. C              \mathbf{U}_{j,R}^{n+1,m}
  665. C              \right)
  666. C              \right\}
  667. C
  668. C
  669.                COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  670. C
  671.                CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG)
  672.      &                          + (COEF * ROD)
  673.                CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG)
  674.      &                          + (COEF * RUXD)
  675.                CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG)
  676.      &                          + (COEF *  RUYD)
  677.                CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG)
  678.      &                          + (COEF *  RETD)
  679.  
  680.             ENDDO
  681. C
  682. C********** End of loop on the faces of the element
  683. C
  684. C           N.B.: ICEN = NLCEG
  685. C
  686. C
  687. C           Computation of the low-Mach preconditioning Matrix
  688. C
  689.             CG=(GAMG*PG/ROG)**0.5D0
  690.             UCOG=MPOCO1.VPOCHA(ICEN,1)
  691.             QG=2.0D0*RECG/ROG
  692.             QG=QG**0.5D0
  693. C
  694.             IF(UCOG .GE. CG)THEN
  695.                PHI2 = 1.0D0
  696.             ELSEIF(QG .GT. CG)THEN
  697.                PHI2 = 1.0D0
  698.             ELSEIF(QG .LT. UCOG)THEN
  699.                PHI2 = UCOG / CG
  700.             ELSE
  701.                PHI2 = QG / CG
  702.             ENDIF
  703.             PHI2=PHI2*PHI2
  704.             PHI2V.VALVEC(ICEN)=PHI2
  705. C
  706.             Q2.VAL(1,ICEN)=(0.5D0*QG*QG)
  707.             Q2.VAL(2,ICEN)=-1.0D0*RUXG/ROG
  708.             Q2.VAL(3,ICEN)=-1.0D0*RUYG/ROG
  709.             Q2.VAL(4,ICEN)=1.0D0
  710. C
  711.             COEF1=(GAMG-1.0D0)/(CG*CG)
  712.             COEF=(1.0D0/COEF1)+(0.5D0*QG*QG)
  713. C
  714.             Q1.VAL(1,ICEN)=COEF1
  715.             Q1.VAL(2,ICEN)=COEF1*RUXG/ROG
  716.             Q1.VAL(3,ICEN)=COEF1*RUYG/ROG
  717.             Q1.VAL(4,ICEN)=COEF1*COEF
  718. C
  719. C
  720. C           Computation of the diagonal coefficients
  721. C
  722.             AN0.VALVEC(ICEN)=AN0.VALVEC(ICEN)
  723. C     &                      +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
  724. C
  725.             DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS)
  726.      &                      +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
  727. C
  728.         ENDDO
  729.       ENDDO
  730. C
  731. C***********************************************************************
  732. C**** END FIRST JACOBI ITERATION ***************************************
  733. C***********************************************************************
  734. C
  735. C***********************************************************************
  736. C**** JACOBI LOOP ******************************************************
  737. C***********************************************************************
  738. C
  739. C
  740.       LOGJAC=.TRUE.
  741. C
  742.       DO IBLJAC=1,NJAC,1
  743. C
  744. C
  745. C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso))
  746. C            alors relaxation SGS avec uniquement des balayages "forward"
  747. C
  748. C Si ICOEF=3 (i.e. si on a "LJACOB")
  749. C            alors relaxation SGS avec uniquement des balayages "backward"
  750. C
  751. C Sinon      (i.e. si on a "LJACOFB")
  752. C            alors relaxation SGS avec des balayages "forward" et "backward"
  753. C
  754.          IF(ICOEF .EQ. 2)THEN
  755.             LOGJAC=.TRUE.
  756.          ELSEIF(ICOEF .EQ. 3)THEN
  757.             LOGJAC=.FALSE.
  758.          ELSE
  759.             LOGINV=LOGJAC
  760.             LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC
  761.             IF(LOGJAC)THEN
  762.                LOGJAC= .NOT. LOGINV
  763.             ELSE
  764.                LOGJAC= LOGINV
  765.             ENDIF
  766.          ENDIF
  767. C
  768.          IF(LOGJAC)THEN
  769.             ICEN=0
  770.             IFACG=0
  771.             ID=1
  772.             ISOUP1=1
  773.             ISOUP2=NSOUPO
  774.          ELSE
  775.             ICEN=NCEN+1
  776.             IFACG=IFACT+1
  777.             ID=-1
  778.             ISOUP1=NSOUPO
  779.             ISOUP2=1
  780.          ENDIF
  781. C
  782.          DO ISOUP=ISOUP1,ISOUP2,ID
  783.             IPT2=MLESUP.LECT(ISOUP)
  784.             NFACE=IPT2.NUM(/1)
  785.             NELT=IPT2.NUM(/2)
  786. C
  787.             IF(LOGJAC)THEN
  788.                IELT1=1
  789.                IELT2=NELT
  790.                IFAC1=1
  791.                IFAC2=NFACE
  792.             ELSE
  793.                IELT2=1
  794.                IELT1=NELT
  795.                IFAC2=1
  796.                IFAC1=NFACE
  797.             ENDIF
  798. C
  799.             DO IELT=IELT1,IELT2,ID
  800.                ICEN=ICEN+ID
  801.                DO IFAC=IFAC1,IFAC2,ID
  802. C                 IFACG increase with IFAC,IELT,ISOUP
  803.                   IFACG = IFACG+ID
  804.                   NGCF  = IPT2.NUM(IFAC,IELT)
  805.                   NLCF  = MLENT2.LECT(NGCF)
  806.                   NGCEG = MELEFL.NUM(1,NLCF)
  807.                   NGCED = MELEFL.NUM(3,NLCF)
  808.                   NLCEG = MLENT1.LECT(NGCEG)
  809.                   NLCED = MLENT1.LECT(NGCED)
  810. C
  811.                   CNX = MPNORM.VPOCHA(NLCF,1)
  812.                   CNY = MPNORM.VPOCHA(NLCF,2)
  813.                   CTX = -1.0D0 * CNY
  814.                   CTY = CNX
  815. C
  816. C**************** Left or right?
  817. C
  818.                   IF(NLCEG .EQ. ICEN)THEN
  819.                      LOGINV=.FALSE.
  820.                   ELSEIF(NLCED .EQ. ICEN)THEN
  821.                      LOGINV=.TRUE.
  822.                   ELSE
  823.                      WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  824.                      WRITE(IOIMP,*) 'Probleme de SPG'
  825. C                    21 2
  826. C                    Données incompatibles
  827.                      CALL ERREUR(21)
  828.                      GOTO 9999
  829.                   ENDIF
  830. C**************** If right, inversion
  831.                   IF(LOGINV) THEN
  832.                      NLCED=NLCEG
  833.                      NLCEG=ICEN
  834.                      CNX = -1.0D0*CNX
  835.                      CNY = -1.0D0*CNY
  836.                      CTX = -1.0D0*CTX
  837.                      CTY = -1.0D0*CTY
  838.                   ENDIF
  839. C
  840.                   SURF   = MPOVSU.VPOCHA(NLCF,1)
  841.                   SIGMAF = SIGMA0.VALVEC(NLCF)
  842. C
  843. C**************** Recuperation des Etats "gauche" et "droite"
  844. C
  845.                   ROG  = CONS1.VAL(1,NLCEG)
  846.                   RUXG = CONS1.VAL(2,NLCEG)
  847.                   RUYG = CONS1.VAL(3,NLCEG)
  848.                   ROG  = CONS1.VAL(1,NLCEG)
  849.                   RUXG = CONS1.VAL(2,NLCEG)
  850.                   RUYG = CONS1.VAL(3,NLCEG)
  851.                   UXG  = RUXG / ROG
  852.                   UYG  = RUYG / ROG
  853.                   UNG  = (UXG * CNX) + (UYG * CNY)
  854.                   UTG  = (UXG * CTX) + (UYG * CTY)
  855.                   RETG = CONS1.VAL(4,NLCEG)
  856.                   GAMG = CONS1.VAL(5,NLCEG)
  857.                   RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
  858.                   RECG = RECG / ROG
  859.                   PG   = (GAMG - 1.0D0) * (RETG - RECG)
  860.                   VOLG = MPOVOL.VPOCHA(NLCEG,1)
  861. C
  862.                   ROD  = CONS1.VAL(1,NLCED)
  863.                   RUXD = CONS1.VAL(2,NLCED)
  864.                   RUYD = CONS1.VAL(3,NLCED)
  865.                   UXD  = RUXD / ROD
  866.                   UYD  = RUYD / ROD
  867.                   RETD = CONS1.VAL(4,NLCED)
  868.                   GAMD = CONS1.VAL(5,NLCED)
  869.                   RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
  870.                   RECD = RECD / ROD
  871.                   PD   = (GAMD - 1.0D0) * (RETD - RECD)
  872.                   UND  = (UXD * CNX) + (UYD * CNY)
  873.                   UTD  = (UXD * CTX) + (UYD * CTY)
  874. C
  875.                   IF(NLCEG .EQ. NLCED)THEN
  876. C                    Murs ou condition limite
  877.                      NLLIM=MLELIM.LECT(NGCF)
  878.                      IF(NLLIM .EQ.0)THEN
  879. C                       Mur
  880.                         UND  = -1.0D0 * UNG
  881.                         UTD  = UTG
  882.                         UXD  = UND * CNX + UTD * CTX
  883.                         UYD  = UND * CNY + UTD * CTY
  884.                         RUXD = ROD*UXD
  885.                         RUYD = ROD*UYD
  886.                      ELSE
  887.                         ROD  = MPLIM.VPOCHA(NLLIM,1)
  888.                         UXD  = MPLIM.VPOCHA(NLLIM,2)
  889.                         UYD  = MPLIM.VPOCHA(NLLIM,3)
  890.                         RUXD = ROD*UXD
  891.                         RUYD = ROD*UYD
  892.                         PD   = MPLIM.VPOCHA(NLLIM,4)
  893.                         GAMD = GAMG
  894.                         UND  = (UXD * CNX) + (UYD * CNY)
  895.                         UTD  = (UXD * CTX) + (UYD * CTY)
  896.                         RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  897.      &                       +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
  898.                      ENDIF
  899.                   ENDIF
  900. C
  901. C**************** We check that density and pressure are positive
  902. C
  903.                   IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  904.      &                 (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  905.                      WRITE(IOIMP,*) 'FMM'
  906.                      WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC
  907.                      WRITE(IOIMP,*) 'ITERATION ICEN',ICEN
  908.                      WRITE(IOIMP,*) ROG,PG
  909.                      WRITE(IOIMP,*) ROD,PD
  910. C
  911. C                    22 0
  912. C******************* Opération malvenue. Résultat douteux
  913. C
  914.                      CALL ERREUR(22)
  915.                      IPROBL = 1
  916.                      GOTO 9998
  917.                   ENDIF
  918. C
  919. C
  920. C**************** We compute the flux
  921. C                 FLU.VALVEC(1:4,IFAC) flux in IFAC face
  922. C                 1 -> mass
  923. C                 2 -> momentum (x)
  924. C                 3 -> momentum (y)
  925. C                 4 -> energy
  926. C
  927.                   RHTD = RETD+PD
  928.                   FR   = (ROD*UND)
  929.                   FRUX = (ROD*UND*UXD)+(PD*CNX)
  930.                   FRUY = (ROD*UND*UYD)+(PD*CNY)
  931.                   FRET = (UND*RHTD)
  932. C
  933.                   COEF = SURF/(2.0D0*VOLG)
  934. C
  935.                   BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  936.      &                             + FR   * COEF
  937.                   BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  938.      &                             + FRUX * COEF
  939.                   BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  940.      &                             + FRUY * COEF
  941.                   BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  942.      &                             + FRET * COEF
  943. C
  944. C**************** Viscous contribution
  945. C
  946.                   COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
  947. C
  948.                   BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  949.      &                             - (COEF * ROD)
  950.                   BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  951.      &                             - (COEF * RUXD)
  952.                   BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  953.      &                             - (COEF *  RUYD)
  954.                   BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  955.      &                             - (COEF *  RETD)
  956. C
  957. C**************** On calcule DIS1
  958. C
  959. C                 DIS1 in the referred report is given by
  960. C                 \sum_j
  961. C                 \left\{
  962. C                 \frac{S_j}{2 \Omega_i} ~
  963. C                 \left(
  964. C                 \sigma_{i,j}^{0}
  965. C                 \mathbf{U}_{j,R}^{n+1,m,l}
  966. C                 \right)
  967. C                 \right\}
  968. C
  969. C
  970.                   COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  971. C
  972.                   CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG)
  973.      &                             + (COEF * ROD)
  974.                   CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG)
  975.      &                             + (COEF * RUXD)
  976.                   CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG)
  977.      &                             + (COEF *  RUYD)
  978.                   CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG)
  979.      &                             + (COEF *  RETD)
  980.  
  981.                ENDDO
  982. C
  983. C************* End of the loop on the faces of the element
  984. C
  985.                DO ICOMP=1,4,1
  986.                   D1.VALVEC(ICOMP) = 0.0D0
  987.                ENDDO
  988. C
  989.                COEF=0.0D0
  990.                DO I1=1,4,1
  991.                   COEF=COEF+(Q2.VAL(I1,ICEN)
  992.      &                *(CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN)))
  993.                ENDDO
  994.                COEF=COEF*(1.0D0-PHI2V.VALVEC(ICEN))
  995.      &                     /PHI2V.VALVEC(ICEN)
  996. C
  997.                DO I1=1,4,1
  998.                   D1.VALVEC(I1)=(CN1.VAL(I1,ICEN)
  999.      &                         - CN0.VAL(I1,ICEN))
  1000.      &                         +(COEF*Q1.VAL(I1,ICEN))
  1001. c                  write(6,*)D1.VALVEC(I1)
  1002.                ENDDO
  1003. C
  1004. C************* At this moment D1 = \Gamma \cdot (\Delta CN1)
  1005. C
  1006.                DO ICOMP=1,4,1
  1007. C
  1008. C************* First increment
  1009. C
  1010.                   D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP)
  1011.      &                 +RES.VAL(ICOMP,ICEN)
  1012.      &                 +(BN0.VAL(ICOMP,ICEN)-BN1.VAL(ICOMP,ICEN)))
  1013.      &                 /(AN0.VALVEC(ICEN)+DN0.VALVEC(ICEN))
  1014.                ENDDO
  1015. C
  1016. C************* At this moment, in the referred report,
  1017. C              D1 is (\delta \mathbf{U}')^{l+1}_i in section A3
  1018. C
  1019. C************* Second increment
  1020. C
  1021.                COEF=0.0D0
  1022.                DO I1=1,4,1
  1023.                   COEF=COEF+(Q2.VAL(I1,ICEN)*D1.VALVEC(I1))
  1024.                ENDDO
  1025. C
  1026. C************* COEF is the scalar product between Q2 and the first
  1027. C              increment.
  1028. C              Let us multiply it by
  1029. C              b_i^0 / (a_i^0 + b_i^0)
  1030. C
  1031.                COEF=COEF*AN0.VALVEC(ICEN)
  1032.      &                  *(PHI2V.VALVEC(ICEN)-1.0D0)
  1033.                COEF=COEF/(AN0.VALVEC(ICEN)+(DN0.VALVEC(ICEN)
  1034.      &                  *PHI2V.VALVEC(ICEN)))
  1035. C
  1036.                DO I1=1,4,1
  1037.                   D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN)
  1038.                ENDDO
  1039. C
  1040. C************* D1 = \delta \mathbf{U}
  1041. C
  1042.                DO ICOMP=1,4,1
  1043.                   DU.VAL(ICOMP,ICEN)    = D1.VALVEC(ICOMP)
  1044. C
  1045.                   CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN)
  1046.      &                                   + D1.VALVEC(ICOMP)
  1047. C
  1048.                   BN1.VAL(ICOMP,ICEN)   = 0.0D0
  1049.                   CN1.VAL(ICOMP,ICEN)   = 0.0D0
  1050.                ENDDO
  1051. C
  1052. C********** End of the loop on the elements
  1053.             ENDDO
  1054. C
  1055. C******* End of the loop on the Domains
  1056.          ENDDO
  1057. C
  1058. C**** End of the Jacobi loop
  1059.       ENDDO
  1060. C
  1061.  9998 CONTINUE
  1062.       DO ICEN=1,NCEN,1
  1063.          DO ICOMP=1,4,1
  1064.             MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
  1065.          ENDDO
  1066.       ENDDO
  1067. C
  1068.       SEGSUP DU
  1069.       SEGSUP RES
  1070.       SEGSUP CONS0
  1071.       SEGSUP CONS1
  1072.       SEGSUP Q1
  1073.       SEGSUP Q2
  1074.       SEGSUP SIGMA0
  1075.       SEGSUP AN0
  1076.       SEGSUP AN0
  1077.       SEGSUP BN0
  1078.       SEGSUP BN1
  1079.       SEGSUP CN0
  1080.       SEGSUP CN1
  1081.       SEGSUP DN0
  1082.       SEGSUP PHI2V
  1083.       SEGSUP USDTI
  1084. C
  1085.       SEGDES MPOVSU
  1086.       SEGDES MPOVDI
  1087.       SEGDES MPOVOL
  1088.       SEGDES MPNORM
  1089.       SEGDES MPOCO
  1090.       SEGDES MPCOEV
  1091. C
  1092.       SEGSUP MPOCO1
  1093. C
  1094.       SEGDES MPODU
  1095. C
  1096.       SEGDES MELEFL
  1097.       SEGSUP MLENT1
  1098.       SEGSUP MLENT2
  1099. C
  1100.       SEGSUP MLELIM
  1101.       IF(MPLIM .GT. 0) SEGDES MPLIM
  1102. C
  1103.       DO ISOUP=1,NSOUPO,1
  1104.          IPT2=MLESUP.LECT(ISOUP)
  1105.          SEGDES IPT2
  1106.       ENDDO
  1107.       IPT2 = MELTFA
  1108.       IF(ISOUP .GT. 1) SEGDES IPT2
  1109.       SEGSUP MLESUP
  1110. C
  1111.  9999 CONTINUE
  1112.       RETURN
  1113.       END
  1114. C
  1115.  
  1116.  
  1117.  
  1118.  
  1119.  

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