Télécharger grad1x.eso

Retour à la liste

Numérotation des lignes :

grad1x
  1. C GRAD1X SOURCE CB215821 19/08/20 21:18:13 10287
  2. C
  3. subroutine GRAD1X (IMODEL,IVADEP,LRE,IVAGRA,NGRA,
  4. & MINT1,MINT2,IIPDPG,IOK)
  5. c
  6. c CALCUL DU GRADIENT DANS LE CAS D'ELEMENTS XFEM (MFR=63)
  7. c
  8. c largement inspiré de epsix.eso + bgrmas.eso
  9. C
  10. C
  11. C*********************************************************
  12. C PARTIE DECLARATIVE
  13. C*********************************************************
  14.  
  15. C
  16. IMPLICIT REAL*8(A-H,O-Z)
  17. C
  18. PARAMETER (NDDLMAX=30,NBNIMAX=10)
  19. C
  20. PARAMETER (NBENRMAX=5)
  21. DIMENSION MLRE(NBENRMAX+1)
  22. C
  23.  
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC SMCOORD
  27. -INC SMMODEL
  28. -INC SMCHAML
  29. -INC SMINTE
  30. -INC SMELEME
  31. -INC SMLREEL
  32. C
  33. POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  34. C
  35. SEGMENT WRK1
  36. REAL*8 XE(3,NBBB)
  37. REAL*8 XDDL(LRE),GRADI(NGRA)
  38. ENDSEGMENT
  39. c
  40. SEGMENT WRK2
  41. REAL*8 SHPWRK(6,NBNO),BGR(NGRA,LRE)
  42. REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  43. ENDSEGMENT
  44. C
  45. SEGMENT MPTVAL
  46. INTEGER IPOS(NS),NSOF(NS)
  47. INTEGER IVAL(NCOSOU)
  48. CHARACTER*16 TYVAL(NCOSOU)
  49. ENDSEGMENT
  50. C
  51. SEGMENT MRACC
  52. INTEGER TLREEL(NBENRMA2,NBI)
  53. ENDSEGMENT
  54. C
  55. C write (*,*) '############################'
  56. C write (*,*) '##### DEBUT DE GRAD1X #####'
  57. C write (*,*) '############################'
  58. C
  59. C
  60. C*********************************************************
  61. c
  62. C RECUPERATION + ACTIVATIONS + VALEURS UTILES
  63. c
  64. C*********************************************************
  65. c
  66. IOK = 0
  67. C++++ Recup de la geometrie deja activée par grad1 ++++++
  68. MELEME= IMAMOD
  69. C nbre de noeuds par element, nbre d elements
  70. NBNN1 = NUM(/1)
  71. NBEL1 = NUM(/2)
  72.  
  73.  
  74. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  75. MELE = NEFMOD
  76. NGAU1 = MINT1.POIGAU(/1)
  77. c C sous decoupage et points de Gauss de l element geometrique de base
  78. c IF(MINT2.NE.0) SEGACT,MINT2
  79. c NGAU2 = MINT2.POIGAU(/1)
  80.  
  81.  
  82. C++++ Recup des infos d enrichissement +++++++++++++++++++
  83. c recup du MCHEX1 d'enrichissement
  84. MCHAM1 = 0
  85. NOBMO1 = IVAMOD(/1)
  86. if (NOBMO1.ne.0) then
  87. do iobmo1=1,NOBMO1
  88. if ((TYMODE(iobmo1)).eq.'MCHAML') then
  89. MCHEX1 = IVAMOD(iobmo1)
  90. segact,MCHEX1
  91. if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  92. MCHAM1 = MCHEX1.ICHAML(1)
  93. segact,MCHAM1
  94. goto 1000
  95. endif
  96. endif
  97. enddo
  98. write(*,*) 'Le modele est vide (absence d enrichissement)'
  99. else
  100. write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
  101. endif
  102.  
  103. 1000 continue
  104. c niveau d enrichissement(s) du modele (ddls std u exclus)
  105. c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc...
  106. if (MCHAM1.ne.0) then
  107. NBENR1= MCHAM1.IELVAL(/1)
  108. else
  109. NBENR1 = 0
  110. endif
  111. C write(*,*) 'niveau d enrichissement(s) du modele',NBENR1
  112. c
  113. C
  114. C
  115. C*********************************************************
  116. C INITIALISATIONS...
  117. C*********************************************************
  118. c
  119. c preparation des tables avec:
  120.  
  121. if (NBENR1.ne.0) then
  122. do ienr=1,NBENR1
  123. c -le nombre d'inconnues de chaque sous-zone
  124. c determinee depuis le nombre de fonction de forme
  125. c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  126. nbniJ = 2 + ((ienr-1)*4)
  127. MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  128. enddo
  129. endif
  130. C Tables + longues car 1er indice correspond au fontion de forme std
  131. MLRE(1) = IDIM*NBNN1*1
  132.  
  133. if (NBENR1.lt.(NBENRMAX+1)) then
  134. do ienr=(NBENR1+1),(NBENRMAX)
  135. MLRE(1+ienr) = 0
  136. enddo
  137. endif
  138. c
  139. c ...DU SEGMENT WRK1
  140. NBENRMA2 = NBENRMAX
  141. NBBB = NBNN1
  142. SEGINI,WRK1
  143.  
  144. c ...DU SEGMENT WRK2
  145. c NBNO = NBNI
  146. NBNO = LRE/IDIM
  147. SEGINI,WRK2
  148. C
  149. c ...DU SEGMENT MRACC
  150. NBENRMA2 = NBENRMAX
  151. NBI = NBNN1
  152. segini,MRACC
  153. C
  154. C du nombre d erreur sur le noms de composantes
  155. NBERR1=0
  156.  
  157.  
  158.  
  159. C*********************************************************
  160. C
  161. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS
  162. C
  163. DO 2000 J=1,NBEL1
  164. C write(*,*) '========= EF',J,' /',NBEL1,'========='
  165. C
  166. C
  167. C*********************************************************
  168. C POUR CHAQUE ELEMENT, ...
  169. C*********************************************************
  170. C
  171. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  172. CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  173. C
  174. c
  175. C++++ NBENRJ = niveau d enrichissement de l element ++++
  176. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  177. NBENRJ=0
  178. if(NBENR1.ne.0) then
  179. do IENR=1,NBENR1
  180. MELVA1 = MCHAM1.IELVAL(IENR)
  181. segact,MELVA1
  182. N2PTEL=MELVA1.IELCHE(/1)
  183. N2EL =MELVA1.IELCHE(/2)
  184. do I=1,NBNN1
  185. IF (N2PTEL.EQ.1 .AND. N2EL.EQ.1)THEN
  186. C Cas du champ constant sur tout le MAILLAGE support
  187. mlree1 = MELVA1.IELCHE(1,1)
  188. ELSEIF(N2PTEL.EQ.1)THEN
  189. C Cas du champ constant par element
  190. mlree1 = MELVA1.IELCHE(1,J)
  191. ELSE
  192. C Cas du champ variable sur les noeuds et les elements
  193. mlree1 = MELVA1.IELCHE(I,J)
  194. ENDIF
  195.  
  196. C on en profite pour remplir MRACC table de raccourcis pour cet element
  197. TLREEL(IENR,I) = mlree1
  198. if (mlree1.ne.0) then
  199. NBENRJ = max(NBENRJ,IENR)
  200. C et on active la listreel
  201. segact,mlree1
  202. endif
  203. enddo
  204. enddo
  205. endif
  206. if (NBENRMA2.gt.NBENR1) then
  207. do IENR2=(NBENR1+1),NBENRMA2
  208. do I=1,NBNN1
  209. TLREEL(IENR2,I) = 0
  210. enddo
  211. enddo
  212. endif
  213. C
  214. c
  215. c++++ quelques indices pour dimensionner au plus juste
  216. c nbre total de ddl de l'élément considéré
  217. c NLIGRD = MLRE(1+NBENRJ)
  218. c NLIGRD = MLRE(1+NBENRJ)
  219. NDDL = MLRE(1+NBENRJ)
  220. c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  221. NBNI = (MLRE(1+NBENRJ)) / IDIM
  222. * nbre d inconnue / noeud
  223. NCOMP = NDDL/NBNN1
  224. C write(6,*) 'EF',J,' NBENRJ=',NBENRJ,
  225. C &'donc',NDDL,' ddls et ',NBNI,' fonctions de forme'
  226. c
  227. c
  228. C
  229. C++++ ON CALCULE XDDL ++++
  230. MPTVAL = IVADEP
  231. NCOSOU = IVAL(/1)
  232. C
  233. INITOT = 0
  234. C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
  235. DO IENR=0,NBENRJ
  236. *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
  237. if (IENR .le. 1) then
  238. NBNIENR = 1
  239. else
  240. NBNIENR = 4
  241. endif
  242. C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
  243. DO INI=1,NBNIENR
  244. INITOT = INITOT + 1
  245. C------>> BOUCLE SUR LA DIMENSION
  246. DO 2220 KDIM=1,IDIM
  247. ICOMP = (INITOT-1)*IDIM + KDIM
  248.  
  249. c --cas ou on n'a pas trouvé assez de composantes--
  250. if(ICOMP.GT.NCOSOU) GOTO 2221
  251.  
  252. MELVAL = IVAL(ICOMP)
  253. c --cas ou on a pas trouvé cette composante dans cette zone du
  254. c chpoint solution devenu mchaml --
  255. if(MELVAL.eq.0)then
  256. c Avait on besoin de cette composante?
  257. c oui, si c'est une composante obligatoire
  258. if(IENR.eq.0) goto 2991
  259. c oui, si l'un des noeuds est enrichi
  260. do I=1,NBNN1
  261. if(TLREEL(IENR,I).ne.0) goto 2991
  262. enddo
  263. c non, si c'est facultatif et qu'on n'est pas enrichi -> on saute
  264. goto 2220
  265. c ->AVERTISSEMENT puis on saute
  266. 2991 NBERR1=NBERR1+1
  267. if(IIMPI.lt.1) goto 2220
  268. c write(IOIMP,*) 'PB OPERATEUR GRAD :'
  269. write(IOIMP,991) ICOMP,IENR,INI,KDIM
  270. 991 format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3,
  271. $ 'ieme COMPOSANTE (enrichissement',I3,
  272. $ ', fonction',I3,', direction ',I3,
  273. $ ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :')
  274. write(IOIMP,*)' noeuds :',(NUM(iou,J),iou=1,NBNN1)
  275. goto 2220
  276. endif
  277.  
  278. C---------->> BOUCLE SUR LES NOEUDS
  279. N1PTEL = VELCHE(/1)
  280. N1EL = VELCHE(/2)
  281. DO I=1,NBNN1
  282. IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
  283. IF (N1PTEL.EQ.1 .AND. N1EL.EQ.1)THEN
  284. C Cas du champ constant sur tout le MAILLAGE support
  285. XDDL(IQ) = VELCHE(1,1)
  286. ELSEIF(N1PTEL.EQ.1) THEN
  287. C Cas du champ constant par element
  288. XDDL(IQ) = VELCHE(1,J)
  289. ELSE
  290. C Cas du champ variable sur les noeuds et les elements
  291. XDDL(IQ) = VELCHE(I,J)
  292. ENDIF
  293. ENDDO
  294.  
  295. 2220 CONTINUE
  296. ENDDO
  297. ENDDO
  298.  
  299. c --cas normal (toutes les composantes souhaitees etaient presentes)--
  300. GOTO 2223
  301.  
  302. c --cas ou on n'a pas trouvé assez de composantes--
  303. 2221 CONTINUE
  304. if (IIMPI.ge.1) then
  305. WRITE(IOIMP,2222) J,NCOSOU,NCOMP
  306. 2222 FORMAT('OPERATEUR GRAD : ELEMENT ',I6,
  307. & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES',
  308. & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES')
  309. endif
  310. NDDL=NCOSOU*NBNN1
  311. NBENRJ=IENR
  312.  
  313. 2223 CONTINUE
  314. c --cas ou on a une ou des erreurs--
  315. IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN
  316. write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ',
  317. & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  318. ENDIF
  319.  
  320. c
  321. c
  322. c rem: il serait probablement interessant au niveau du tems cpu
  323. c d'utiliser moins de pts de Gauss lorsque l element est élastique
  324. c On pourrait par ex. utiliser MINT2 = infell(12) qui contient
  325. c le segment d'integration de l'EF std (QUA4 par ex.)
  326. * if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  327. * MINTE = MINT2
  328. * NBPGAU= NGAU2
  329. * else
  330. MINTE = MINT1
  331. NBPGAU= NGAU1
  332. * endif
  333. c
  334. C
  335. C*********************************************************
  336. C
  337. C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS
  338. C
  339. DO 2500 KGAU=1,NBPGAU
  340. C
  341. C*********************************************************
  342. C Initialisation à 0
  343. C*********************************************************
  344.  
  345. c ZERO ne serait-il pas facultatif?
  346. CALL ZERO(SHPWRK,6,NBNO)
  347. XGAU = 0.D0
  348. YGAU = 0.D0
  349. ZGAU = 0.D0
  350. C
  351. i6zz = 3
  352. IF(IDIM.EQ.3) i6zz = 4
  353. C
  354. do ienr7=1,NBENRJ
  355. do inod7=1,NBNN1
  356. do i6=1,i6zz
  357. LV7WRK(ienr7,1,i6,inod7) = 0.D0
  358. LV7WRK(ienr7,2,i6,inod7) = 0.D0
  359. enddo
  360. enddo
  361. enddo
  362.  
  363.  
  364. c*********************************************************
  365. c Calcul des fonction de forme std dans repere local
  366. c*********************************************************
  367.  
  368. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  369. c (et donc sur les Ni std)
  370. DO 2510 I=1,NBNN1
  371.  
  372. C++++ Calcul des Ni std
  373. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  374. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  375. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  376. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  377. IF(IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  378.  
  379. C++++ Calcul des coordonnees dans le repere global
  380. XGAU = XGAU + (SHPWRK(1,I)*XE(1,I))
  381. YGAU = YGAU + (SHPWRK(1,I)*XE(2,I))
  382. IF(IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I))
  383.  
  384. 2510 CONTINUE
  385. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  386. c if(J.eq.77) write(6,*) 'xgau_',KGAU,' =',XGAU,YGAU,ZGAU
  387.  
  388.  
  389.  
  390. c*********************************************************
  391. c Passage des fonctions de forme std dans repere global
  392. c*********************************************************
  393.  
  394. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  395. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  396. c if(J.eq.1.and.KGAU.eq.1)
  397. c & write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  398.  
  399. c*********************************************************
  400. c Si on est pas du tout enrichi on peut sauter une grosse partie
  401. c*********************************************************
  402. if(NBENRJ.eq.0) goto 2999
  403.  
  404. c*********************************************************
  405. c Calcul des level set + leurs derivees dans repere global
  406. c*********************************************************
  407.  
  408. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  409. do 2520 ienr=1,NBENRJ
  410.  
  411. c MELVA1=MCHAM1.IELVAL(IENR)
  412. c segact,MELVA1
  413.  
  414. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  415. DO 2521 I=1,NBNN1
  416.  
  417. C++++ Le I eme noeud est-il ienr-enrichi?
  418. mlree1= TLREEL(IENR,I)
  419. c if(J.eq.77) write(6,*) 'ienr,I',ienr,I,' mlree1=',mlree1
  420. if(mlree1.eq.0) goto 2521
  421.  
  422.  
  423. c Calcul du repere local de fissure(=PSI,PHI)
  424. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  425. do 2522 inode=1,NBNN1
  426. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  427. if (ienr.ne.1) then
  428. c valeur de PSI au inode^ieme noeud
  429. xpsi1 = mlree1.PROG(inode)
  430. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  431. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  432. & + (SHPWRK(1,inode)*xpsi1)
  433. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  434. & + (SHPWRK(2,inode)*xpsi1)
  435. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  436. & + (SHPWRK(3,inode)*xpsi1)
  437. IF(IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  438. & + (SHPWRK(4,inode)*xpsi1)
  439. c valeur de PHI au inode^ieme noeud
  440. xphi1 = mlree1.PROG(NBNN1+inode)
  441. else
  442. xphi1 = mlree1.PROG(inode)
  443. endif
  444. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  445. & + (SHPWRK(1,inode)*xphi1)
  446. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  447. & + (SHPWRK(2,inode)*xphi1)
  448. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  449. & + (SHPWRK(3,inode)*xphi1)
  450. IF(IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  451. & + (SHPWRK(4,inode)*xphi1)
  452. 2522 continue
  453.  
  454. 2521 continue
  455. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  456.  
  457.  
  458. 2520 CONTINUE
  459. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  460.  
  461. c on a construit
  462. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  463.  
  464.  
  465. c*********************************************************
  466. c Ajout des fonctions de forme d enrichissement
  467. c + leurs derivees dans repere global
  468. c*********************************************************
  469. CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  470.  
  471.  
  472.  
  473. c retour a la partie commune aux EF enrichi et non enrichi
  474. 2999 continue
  475.  
  476. C*********************************************************
  477. C CALCUL DE BGR
  478. C*********************************************************
  479. c ZERO ne serait-il pas facultatif?
  480. call ZERO(BGR,9,NDDL)
  481. KB=1
  482. C boucle sur tous les Ni
  483. DO 3001 II=1,NBNI
  484.  
  485. c partie commune 2D 3D
  486. BGR(1,KB) = SHPWRK(2,II)
  487. BGR(2,KB) = SHPWRK(3,II)
  488. BGR(4,KB+1) = SHPWRK(2,II)
  489. BGR(5,KB+1) = SHPWRK(3,II)
  490.  
  491. c cas 3D
  492. IF (IDIM.EQ.3) THEN
  493. BGR(3,KB) = SHPWRK(4,II)
  494. BGR(6,KB+1) = SHPWRK(4,II)
  495. BGR(7,KB+2) = SHPWRK(2,II)
  496. BGR(8,KB+2) = SHPWRK(3,II)
  497. BGR(9,KB+2) = SHPWRK(4,II)
  498. ENDIF
  499.  
  500. c cas 2D PLAN GENE
  501. IF (IFOUR.EQ.-3) THEN
  502. IREF=(IIPDPG-1)*(IDIM+1)
  503. BGR(9,KB) = 1.D0
  504. BGR(9,KB+1) = XCOOR(IREF+1)-XGAU
  505. BGR(9,KB+2) = YGAU-XCOOR(IREF+2)
  506. ENDIF
  507.  
  508. KB = KB + IDIM
  509.  
  510. 3001 CONTINUE
  511. C
  512. c
  513.  
  514. C*********************************************************
  515. C CALCUL DE BGR.q
  516. C*********************************************************
  517.  
  518. C APPEL A LA PROCEDURE DE CALCUL
  519. CALL BGRDEP(BGR,NGRA,XDDL,NDDL,GRADI)
  520. c
  521.  
  522. C*********************************************************
  523. C ECRITURE DES DIFFERENTES COMPOSANTES DU GRADIENT
  524. C*********************************************************
  525. MPTVAL=IVAGRA
  526. DO i=1,NGRA
  527. MELVAL=IVAL(i)
  528. IGMN=MIN(KGAU,VELCHE(/1))
  529. IBMN=MIN(J,VELCHE(/2))
  530. VELCHE(IGMN,IBMN)=GRADI(i)
  531. ENDDO
  532. c
  533. c
  534. 2500 CONTINUE
  535. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  536. C
  537. cc
  538. c
  539. c
  540. 2000 CONTINUE
  541. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  542. c
  543. c
  544.  
  545. C*********************************************************
  546. C SUPPRESSION ET DESACTIVATION DE SEGMENTS
  547. C*********************************************************
  548. C
  549. SEGSUP,WRK1,WRK2
  550. SEGSUP,MRACC
  551.  
  552. c on met mint2 a zero pour eviter pb dans grad1
  553. MINT2 = 0
  554. c
  555. c pas de retour en erreur pour linstant
  556. IOK=1
  557. c write(6,*) 'iok=',IOK
  558. C
  559. END
  560.  
  561.  
  562.  

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