Télécharger pplast.eso

Retour à la liste

Numérotation des lignes :

  1. C PPLAST SOURCE CB215821 16/04/21 21:18:03 8920
  2.  
  3. SUBROUTINE PPLAST(XMAT,NMATT,VAR0,VARF,NVARI,SIG0,
  4. & SIGF,DEPST,NSTRS,KERRE)
  5.  
  6. ***************************************
  7. * Routine de traitement en plasticité *
  8. ***************************************
  9.  
  10. * Entrées / Sorties
  11.  
  12. * Entrées
  13.  
  14. INTEGER NMATT, NVARI, NSTRS
  15. * NMATT : Nombre de parametres materiaux
  16. * NVARI : Nombre de variables internes
  17. * NSTRS : Nombre de composantes
  18.  
  19. REAL*8 XMAT(NMATT), VAR0(NVARI), SIG0(NSTRS)
  20. REAL*8 DEPST(NSTRS)
  21. * XMAT : Paramètres matériaux
  22. * VAR0 : Variables internes
  23. * SIG0 : Contraintes reelles
  24. * DEPST : Increment de deformations totales
  25.  
  26. * Sorties
  27.  
  28. REAL*8 SIGF(NSTRS), VARF(NVARI)
  29. * SIGF : Contraintes effective
  30. * VARF : Variables internes en sortie de plasticite
  31.  
  32. * Variables materiau
  33. REAL*8 E,NU
  34. REAL*8 AC, AT, BC, BT, EPSD0, FC, FT
  35. REAL*8 P, AH, BH, CH, GAMMA, ALFA
  36. REAL*8 A, K0
  37. REAL*8 DEUXMU, LAMB, TREPS
  38.  
  39. * Données intermédiaires et autres
  40. INTEGER ITER, I , J, K, KGLOB, KTEST
  41. LOGICAL CONV
  42. REAL*8 TEST,F0
  43. REAL*8 LIM, DEPST2(6)
  44. REAL*8 SIGE0(6), SIGE0P(3), TOL, TOLDYN
  45. REAL*8 VECPC(3,3), NORM, DSIGT(6)
  46. REAL*8 DFDS(6), H, SIG02(6), VAR02(9)
  47. REAL*8 FTB, FH, FM(6), C
  48. REAL*8 TRS, KH, KH0, KH2, TR(3,3)
  49. REAL*8 DMP, RES(8), RESUM(6), RESU(6), KHM, ELAS(6,6)
  50. REAL*8 JAC(7,7)
  51. REAL*8 DFDK, N(6), DHDK, KH1, DHDS(6), DMDK(6)
  52. REAL*8 FHP(6), FMM1(6)
  53. REAL*8 FHM(6)
  54. REAL*8 FMP1(6)
  55. REAL*8 DMDS(6,6), FHPH, FHMH, ETA, ETAP
  56. REAL*8 FHM1
  57. REAL*8 FHP1, SIGP(6), SIGMM(6), DELTA(8)
  58. REAL*8 SIGMP1(6),DSDE(7,7)
  59. REAL*8 M
  60. LOGICAL PLASTI,CONV1,CONV2,TEST5
  61. INTEGER FLAG
  62.  
  63. REAL*8 ELASM(6,6)
  64. REAL*8 ES(7), NUMF, NUM(7), DENOF
  65. REAL*8 DENO(7), MAXI, MAXI2, SIGE00(6)
  66.  
  67.  
  68.  
  69. ***************************************
  70. * Récupération des variables matériau *
  71. ***************************************
  72.  
  73. E = XMAT(1)
  74. NU = XMAT(2)
  75. AC = XMAT(5)
  76. AT = XMAT(6)
  77. BC = XMAT(7)
  78. BT = XMAT(8)
  79. EPSD0 = XMAT(9)
  80. FC = XMAT(10)
  81. FT = XMAT(11)
  82. P = XMAT(12)
  83. AH = XMAT(13)
  84. BH = XMAT(14)
  85. CH = XMAT(15)
  86. GAMMA = XMAT(16)
  87. ALFA = XMAT(17)
  88. A = XMAT(18)
  89. K0 = XMAT(19)
  90.  
  91. DEUXMU=E/(1.D0+NU)
  92. LAMB=NU*DEUXMU/(1.D0-2.D0*NU)
  93.  
  94.  
  95.  
  96. *******************************
  97. * Initialisationdes variables *
  98. *******************************
  99.  
  100. * Calcul de la matrice élastique (une seule fois !)
  101.  
  102. DO I = 1,6
  103. DO J = 1,6
  104. ELAS(I,J) = 0.d0;
  105. ELASM(I,J) = 0.d0;
  106. END DO
  107. END DO
  108.  
  109. ELAS(1,1)=LAMB+DEUXMU
  110. ELAS(2,2)=LAMB+DEUXMU
  111. ELAS(3,3)=LAMB+DEUXMU
  112. ELAS(1,2)=LAMB
  113. ELAS(2,1)=LAMB
  114. ELAS(1,3)=LAMB
  115. ELAS(3,1)=LAMB
  116. ELAS(2,3)=LAMB
  117. ELAS(3,2)=LAMB
  118. ELAS(4,4)=DEUXMU
  119. ELAS(5,5)=DEUXMU
  120. ELAS(6,6)=DEUXMU
  121.  
  122. ELASM(1,1)=1.D0/E
  123. ELASM(2,2)=1.D0/E
  124. ELASM(3,3)=1.D0/E
  125. ELASM(1,2)=(-NU)/E
  126. ELASM(2,1)=(-NU)/E
  127. ELASM(1,3)=(-NU)/E
  128. ELASM(3,1)=-NU/E
  129. ELASM(2,3)=-NU/E
  130. ELASM(3,2)=-NU/E
  131. ELASM(4,4)=(1.D0+NU)/E
  132. ELASM(5,5)=(1.D0+NU)/E
  133. ELASM(6,6)=(1.D0+NU)/E
  134.  
  135.  
  136.  
  137. ******************
  138. * Initialisation *
  139. ******************
  140.  
  141. * Variables tampons
  142.  
  143. DO I=1,NSTRS
  144. DEPST2(I) = DEPST(I)
  145. SIG02(I) = SIG0(I)
  146. END DO
  147.  
  148. DO I=1,NVARI
  149. VAR02(I) = VAR0(I)
  150. END DO
  151.  
  152. KGLOB = 1
  153. KTEST = 1
  154.  
  155. * Test d'iterations globales
  156. TEST5 = .TRUE.
  157.  
  158. ***********************
  159. * Iterations globales *
  160. ***********************
  161.  
  162. DO WHILE (TEST5)
  163.  
  164. * On saisit une petite valeur de KH et on saisit le plus grand à chaque itération...
  165. KH0 = 1.d-3
  166. KH = MAX(VAR0(9), KH0)
  167.  
  168. * Initialisation de la variable KH à l'instant t
  169. KHM = KH
  170.  
  171. * Initialisation de l'incrément du multiplicateur plastique
  172. DMP = 0.d0
  173.  
  174.  
  175. * Variables de test
  176. TEST = 1.d-6
  177. FLAG = 0
  178.  
  179.  
  180.  
  181. ****************************************************
  182. * Calcul de l'incrément de contraintes (élastique) *
  183. ****************************************************
  184.  
  185. * Prediction elastique
  186. TREPS = DEPST(1) + DEPST(2) + DEPST(3)
  187. DO K=1,3
  188. DSIGT(K) = LAMB*TREPS + DEUXMU*DEPST(K)
  189. END DO
  190. DO K=4,NSTRS
  191. DSIGT(K) = DEUXMU*DEPST(K)
  192. END DO
  193.  
  194. * Initialisation des variables plastiques
  195. CONV=.FALSE.
  196.  
  197. DO I=1, NSTRS
  198. DFDS(I) = 0.d0
  199. END DO
  200.  
  201. DO K = 1, NSTRS
  202. SIGE0(K) = VAR0(K+2)
  203. END DO
  204.  
  205. ITER=0
  206.  
  207. * Prediction élastique (NSTRS : nombre de composantes de SIGMA)
  208. DO ISTRS = 1,NSTRS
  209.  
  210. * La prédiction élastique est donnée par la contrainte réelle SIG0
  211. * à laquelle on somme la différence entre contrainte plastique et
  212. * contrainte réelle stockée en variable interne et la prédiction
  213. * élastique
  214.  
  215. * Initialisation pour une itération lorsque l'incrément est divisé :
  216. * - Si c'est la première itération, on reprend la valeur initiale de
  217. * la contrainte plastique pour créer la prédiction
  218. * - Si c'est l'itération suivante, on travaille directement sur
  219. * la contrainte plastique obtenue à l'itération précédente et
  220. * stockée dans SIG0. On ne recalcule pas la contrainte plastique.
  221.  
  222. IF (KTEST.EQ.1) THEN
  223.  
  224. * SigE0 représente la différence entre la valeur de la contrainte dans le cas réel
  225. * et la valeur de la contrainte dans le cas plastique à l'instant t
  226. SIGMP1(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS)
  227. SIGE00(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS)
  228. SIGE0(ISTRS) = SIG0(ISTRS) + SIGE0(ISTRS) + DSIGT(ISTRS)
  229.  
  230. ELSE
  231. SIGMP1(ISTRS) = SIG0(ISTRS)
  232. SIGE00(ISTRS) = SIG0(ISTRS)
  233. SIGE0(ISTRS) = SIG0(ISTRS) + DSIGT(ISTRS)
  234.  
  235. END IF
  236.  
  237. END DO
  238.  
  239. * Apres, SIGE0 represente la solution (contrainte effective)
  240.  
  241. *************************************
  242. * Itération de calcul de plasticité *
  243. *************************************
  244.  
  245. DO WHILE (ITER .LE. 20 .AND. .NOT. CONV)
  246. ITER=ITER+1
  247.  
  248. * Critere sur la prediction permettant de ramener la solution
  249. * a un domaine dans lequel les equations du modele sont definies
  250. * Note : Il a été modifié de telle sorte que l'on applique le critère jusqu'à
  251. * ce que la prédiction soit acceptable (10 itérations max modifiable)
  252. TRS = 0.D0
  253. DO I = 1,3
  254. TRS = TRS + SIGE0(I)
  255. END DO
  256.  
  257. FTB = FT / FC
  258. M = 3.D0 * (1.D0-(FTB)**(2.D0/GAMMA))/(FTB +
  259. & 2.D0*(FTB)**(1.D0/GAMMA))
  260. LIM = 3.D0*FC/M
  261.  
  262. DO WHILE (TRS.GT.LIM .AND. FLAG.LT.10)
  263.  
  264. FLAG = FLAG + 1
  265.  
  266. IF (FLAG .LT. 10) THEN
  267. DO I =1,NSTRS
  268. SIGE0(I) = SIGE00(I)
  269. DSIGT(I) = DSIGT(I)/2.d0
  270. SIGE0(I) = SIGE00(I) + DSIGT(I)
  271. END DO
  272. DMP = 0.d0
  273. KH = max(var0(9),KH0)
  274. ITER = 2
  275.  
  276. TRS = 0.D0
  277.  
  278. DO I = 1,3
  279. TRS = TRS + SIGE0(I)
  280. END DO
  281.  
  282. IF(TRS.GT.LIM) THEN
  283. write(6,*) 'ERREUR'
  284. END IF
  285.  
  286. ELSE
  287.  
  288. DO I =1,6
  289. SIGE0(I) = SIGE00(I)
  290. END DO
  291.  
  292. END IF
  293.  
  294. END DO
  295.  
  296. * Calcul de FM et FH
  297. CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C,
  298. & SIGE0, KH, F0, FM, ETA)
  299.  
  300. CALL FLOWH(SIGE0, KH, AH, BH, CH, FC, FM, FH)
  301.  
  302. * Calcul des résidus + vérification du critère de CV
  303. DO I=1,6
  304. RESU(I) = 0.D0
  305. RESUM(I) = 0.D0
  306. DO J=1,6
  307. RESU(I) = RESU(I) + ELASM(I,J)*SIGMP1(J)
  308. RESUM(I) = RESUM(I) + ELASM(I,J)*SIGE0(J)
  309. END DO
  310. END DO
  311.  
  312. DO I = 1,6
  313. RES(I) = RESUM(I) + DMP*FM(I) - DEPST(I) - RESU(I)
  314. END DO
  315.  
  316. RES(7) = KH - DMP*FH - KHM
  317. RES(8) = F0
  318.  
  319. IF (KH .GE. 1.D0) THEN
  320. RES(7) = 0.D0
  321. KH = 1.d0
  322. END IF
  323.  
  324. * Convergence
  325. * Deux criteres : un absolu, un relatif... Critere relatif conservé
  326. * Les variations de contraintes sont normées avec le max des contraintes
  327. * en fin d'itération et celles de KH et lambda sont normées par les valeurs
  328. * associées en sortie d'itération
  329.  
  330. IF (((ABS(RES(8)).LT.TEST) .OR. (RES(8).LE.0.D0))
  331. & .AND. ITER.EQ.1 .AND. (FLAG.EQ.0)) THEN
  332. CONV = .TRUE.
  333. ELSE
  334. CONV1 = .TRUE.
  335. MAXI = ABS(RES(1))
  336. MAXI2 = ABS(SIGE0(1))
  337.  
  338. DO I=1,8
  339. MAXI = MAX(ABS(RES(I)),MAXI)
  340. END DO
  341. DO I=1,6
  342. MAXI2 = MAX(ABS(SIGE0(I)),MAXI2)
  343. END DO
  344.  
  345. IF (MAXI .LT. TEST) THEN
  346. CONV1 = .TRUE.
  347. ELSE
  348. CONV1 = .FALSE.
  349. END IF
  350.  
  351. CONV2 = .TRUE.
  352. DO I=1,NSTRS
  353. IF ((ABS(DELTA(I)) .LT. MAXI2*TEST)
  354. & .AND. (CONV2)) THEN
  355. CONV2 = .TRUE.
  356. ELSE
  357. CONV2 = .FALSE.
  358. END IF
  359. END DO
  360.  
  361. IF ((ABS(DELTA(7)) .LT. ABS(KH)*TEST) .AND. (CONV2)) THEN
  362. CONV2 = .TRUE.
  363. ELSE
  364. CONV2 = .FALSE.
  365. END IF
  366.  
  367. IF ((ABS(DELTA(8)) .LT. ABS(DMP)*TEST) .AND. (CONV2)) THEN
  368. CONV2 = .TRUE.
  369. ELSE
  370. CONV2 = .FALSE.
  371. END IF
  372.  
  373. IF (CONV2) THEN
  374. CONV = .TRUE.
  375. END IF
  376.  
  377. IF (ITER .EQ. 1) THEN
  378. CONV = .FALSE.
  379. END IF
  380.  
  381. END IF
  382.  
  383. IF (.NOT. CONV) THEN
  384.  
  385. PLASTI = .TRUE.
  386.  
  387.  
  388.  
  389. * Calcul de la dérivée de m : différentiation numérique
  390. * Note : le choix de h joue un role non négligeable
  391. * dans l'efficacité de l'algorithme
  392. CALL ENDOCB(SIGE0,TR,1,2)
  393. CALL JACOB3(TR,3,SIGE0P,VECPC)
  394.  
  395. NORM= 0.D0
  396. DO I = 1,3
  397. NORM=NORM +(SIGE0P(I)**2)
  398. ENDDO
  399. NORM = SQRT(NORM)
  400.  
  401. H= NORM*2.d0**(-12)
  402.  
  403.  
  404.  
  405. * Différenciation numérique
  406. DO J = 1,6
  407.  
  408. DO I = 1,6
  409. SIGP(I) = SIGE0(I)
  410. SIGMM(I)= SIGE0(I)
  411. ENDDO
  412.  
  413. SIGP(J) = SIGP(J) + H
  414. SIGMM(J) = SIGMM(J) - H
  415.  
  416. CALL FLOWM (FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH,
  417. & C, SIGP, KH, F0, FMP1, ETAP)
  418. CALL FLOWM (FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH,
  419. & C, SIGMM, KH, F0, FMM1, ETAP)
  420.  
  421. CALL FLOWH(SIGP, KH, AH, BH, CH, FC, FMP1, FHP1)
  422. CALL FLOWH(SIGMM, KH, AH, BH, CH, FC, FMM1, FHM1)
  423.  
  424.  
  425.  
  426. * Et calcul de DMDS(i,j) = DM(i)/DS(J)
  427.  
  428. DO I = 1,6
  429. DMDS(I,J) = (FMP1(I)-FMM1(I))/(2.D0*H)
  430. DHDS(J) = (FHP1-FHM1)/(2.D0*H)
  431. END DO
  432.  
  433. END DO
  434.  
  435. H=KH*2.d0**(-12)
  436.  
  437. KH1= KH+H
  438.  
  439. IF (KH1.GT.1.D0) THEN
  440. KH1 = KH
  441. END IF
  442.  
  443. CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C,
  444. & SIGE0, KH1, F0, FHP, ETAP)
  445. CALL FLOWH(SIGE0, KH1, AH, BH, CH, FC, FHP, FHPH)
  446.  
  447.  
  448. KH1 = KH -H
  449.  
  450. * Traitement du cas kh = 1 : d2FdK2 = 0
  451.  
  452. IF (KH.GE.1.D0) THEN
  453. KH1=KH
  454. END IF
  455.  
  456. CALL FLOWM(FC, FT, GAMMA, ALFA, P, K0, A, AH, BH, CH, C,
  457. & SIGE0, KH1, F0, FHM, ETAP)
  458.  
  459. CALL FLOWH(SIGE0, KH1, AH, BH, CH, FC, FHM, FHMH)
  460.  
  461.  
  462. DO I=1,6
  463. DMDK(I) = (FHP(I) - FHM(I))/(2.D0*H)
  464. END DO
  465.  
  466. DHDK = (FHPH-FHMH)/(2.D0*H)
  467.  
  468.  
  469.  
  470. * Calcul de N = DFDS
  471.  
  472. DO I = 1,6
  473. N(I) = FM(I)
  474. END DO
  475.  
  476. DFDK = ETA
  477.  
  478. DO I=1,6
  479. DO J=1,6
  480. JAC(I,J)= ELASM(I,J) + DMP*DMDS(I,J)
  481. ENDDO
  482. JAC(I,7) = DMP*DMDK(I)
  483. END DO
  484.  
  485. DO J=1,6
  486. JAC(7,J) = -DMP*DHDS(J)
  487. END DO
  488.  
  489. JAC(7,7) = 1.D0-DMP*DHDK
  490.  
  491. CALL INVALM(JAC,7,7,KERRE,1.d-20)
  492.  
  493. IF (KERRE.EQ.1) THEN
  494. WRITE(6,*) 'Erreur d inversion de matrice'
  495. END IF
  496.  
  497.  
  498.  
  499. ****************************************************************
  500. * Calcul des variations en sortie d'itération (Newton Raphson) *
  501. ****************************************************************
  502.  
  503. DO I=1,7
  504. DENO(I) = 0.d0
  505. NUM(I) = 0.d0
  506. END DO
  507.  
  508. DENOF=0
  509. NUMF=0
  510.  
  511. DO I = 1,8
  512. DELTA(I) = 0.d0
  513. END DO
  514.  
  515. DO I=1,7
  516. DO J=1,6
  517. DENO(I) = DENO(I)+JAC(I,J)*FM(J)
  518. END DO
  519. DENO(I) = DENO(I)-JAC(I,7)*FH
  520. END DO
  521.  
  522. DO I=1,6
  523. DENOF = DENOF+N(I)*DENO(I)
  524. END DO
  525.  
  526. DENOF = DENOF+DFDK*DENO(7)
  527.  
  528. DO I =1,7
  529. DO J = 1,7
  530. NUM(I)= NUM(I)+JAC(I,J)*RES(J)
  531. END DO
  532. END DO
  533.  
  534. DO I =1,6
  535. NUMF = NUMF + N(I)*NUM(I)
  536. END DO
  537.  
  538. NUMF = NUMF + DFDK*NUM(7)
  539. NUMF = RES(8)-NUMF
  540.  
  541. DELTA(8) = NUMF/DENOF
  542.  
  543. DO I =1,6
  544. ES(I) = -RES(I)-DELTA(8)*FM(I)
  545. END DO
  546.  
  547. ES(7) = -RES(7)+DELTA(8)*FH
  548.  
  549. DO I = 1,7
  550. DO J = 1,7
  551. DELTA(I)=DELTA(I)+JAC(I,J)*ES(J)
  552. END DO
  553. END DO
  554.  
  555.  
  556. * B est un vecteur à 8 composantes composées des DELTA (S, kh et LAMBDA)
  557.  
  558. DO I = 1, 6
  559. SIGE0(I) = SIGE0(I) + DELTA(I)
  560. END DO
  561.  
  562. KH2 = KH
  563. KH = KH + DELTA(7)
  564.  
  565.  
  566.  
  567. * Test de KH < 1
  568.  
  569. IF ((KH2.GE. 1.D0).OR.(KH.GE.1.D0)) THEN
  570. KH = 1.D0
  571. END IF
  572.  
  573.  
  574.  
  575. * KH ne diminue pas...
  576.  
  577. IF (KH .LT. KHM) THEN
  578. KH = KHM
  579. END IF
  580.  
  581. DMP = DMP + DELTA(8)
  582.  
  583. END IF
  584.  
  585. END DO
  586.  
  587.  
  588.  
  589. *******************
  590. * Non convergence *
  591. *******************
  592.  
  593. IF ((ITER .EQ. 21) .AND. (.NOT. CONV)
  594. & .AND. (KGLOB.EQ.32)) THEN
  595. WRITE(6,*) '**************************************'
  596. WRITE(6,*) '* Divergence de l algorithme interne *'
  597. WRITE(6,*) '**************************************'
  598. WRITE(6,*) 'DEPST :',DEPST2
  599. WRITE(6,*) 'DELTA :',DELTA
  600. WRITE(6,*) 'MAXSIG,KH,DMP :',MAXI2,KH,DMP
  601. WRITE(6,*) 'VAR0 :',VAR02
  602. WRITE(6,*) 'SIG0 :',SIG02
  603. WRITE(6,*) 'SIGF :',SIGF
  604. * KERRE permet de stopper CAST3M et de reprendre la main
  605. * lorsqu'il est non nul
  606. KERRE = 1
  607. END IF
  608.  
  609.  
  610.  
  611. *****************************************************
  612. * Traitement des données en fin d'itération globale *
  613. *****************************************************
  614.  
  615. DO I=1,NSTRS
  616. SIG0(I) = SIGE0(I)
  617. END DO
  618.  
  619. VAR0(9) = KH
  620.  
  621. * Test de convergence
  622. IF (KGLOB.EQ.KTEST .AND. CONV) THEN
  623. TEST5 = .FALSE.
  624. ELSE IF (CONV) THEN
  625. KTEST = KTEST + 1
  626. END IF
  627.  
  628. IF (.NOT.CONV) THEN
  629. KGLOB = KGLOB*2
  630. KTEST = 1
  631. DO I = 1,NSTRS
  632. SIG0(I) = SIG02(I)
  633. DEPST(I) = DEPST(I)/2.d0
  634. END DO
  635. VAR0(9) = VAR02(9)
  636. END IF
  637.  
  638. * Limite de la division de l'intervalle
  639. * On redivise par 2 pour l'afficher en sortie
  640. IF (KGLOB.GT.32) THEN
  641. TEST5=.FALSE.
  642. KGLOB = KGLOB/2
  643. END IF
  644.  
  645. END DO
  646.  
  647.  
  648.  
  649. *******************************************
  650. * FIN DE LA BOUCLE GLOBALE *
  651. * --> Sorties de la routine de plasticite *
  652. *******************************************
  653.  
  654. * Les résultats finaux sont stockés dans SIG0 et VAR0
  655. * Les valeurs en entrée sont stockés dans SIG02 et VAR02
  656. DO I = 1,NSTRS
  657. SIGF(I) = SIG0(I)
  658. SIG0(I) = SIG02(I)
  659. DEPST(I) = DEPST2(I)
  660. END DO
  661.  
  662. VARF(9) = VAR0(9)
  663. VAR0(9) = VAR02(9)
  664.  
  665. * On affiche s'il y a eu division de l'incrément de déformations
  666. IF (KGLOB.GT.1) THEN
  667. * WRITE(*,*) 'Division de l increment de deformations par',KGLOB
  668. * WRITE(*,*) 'D :',VAR02(2),' et KH :',KH
  669. END IF
  670.  
  671. RETURN
  672.  
  673. END
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  

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