C FPTAU SOURCE PV 08/09/11 21:15:54 6150 SUBROUTINE FPTAU(IDIM,NP,NBEL,NUM,IPADT,XCOOR, & U,NPTI,NC,IPADS,NPTS,ANUT,AMUEF, & YP,AN,RO,IKR,AMU,IKM,UET,YPS,ITYPEL,TKE,NBINC) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C....................................................................... C Ut(NP,NBEL) module vitesse tangente = |U - Un| C Un(NP,NBEL) module vitesse normale = Un C YP(NP) distance la paroi C UET(*) uet aux noeuds frontieres C AN(IDIM,NP) normale C ANN module de la normale C NUM(NP) connectivite element appuyé C NMS(NS) connectivite relative (frontiere) C NMY(NY) connectivite complementaire C RO(*) densite C AMUEF(NP,NBEL) Viscosite effective aux noeuds/elt C C C1/ ce qu'il faudrait faire c'est attribuer une normale pour chaque point C de l'element appuye -> c'est fait C....................................................................... -INC CCREEL DIMENSION LKKE(1000) real*8 ann(3),Xa(9),Ya(9),Za(9) DIMENSION XCOOR(*),XYZ(3,27) DIMENSION NUM(NP,NBEL),NMS(9),NMY(27) DIMENSION U(NPTI,NC),IPADT(*),VS(27,3),IPADS(*) DIMENSION ANUT(NPTI),AMUEF(NP,NBEL) DIMENSION YP(NP,NBEL),AN(IDIM,NP,NBEL) DIMENSION UT(27),UN(27) DIMENSION RO(*),AMU(*) DIMENSION UET(NPTS),YPS(NPTS),TKE(NPTS,*) DIMENSION NUMFA(27) DATA EPS/1.D-4/,NBIT/100/ c----------------------------------------------------------------------- c------- Caracterisation des elements ---------------------------------- c----------------------------------------------------------------------- DIMENSION NMF9(3,4),N1F9(3,4),N2F9(3,4) DATA NMF9/1,2,3,3,4,5,5,6,7,7,8,1/ DATA N1F9/8,9,4,2,9,6,4,9,8,6,9,2/,N2F9/7,6,5,1,8,7,3,2,1,5,4,3/ c----------------------------------------------------------------------- DIMENSION NMF7(3,3),N1F7(3,3),N2F7(1,3) DATA NMF7/1,2,3,3,4,5,5,6,1/ DATA N1F7/6,7,4,2,7,6,4,7,2/,N2F7/5,1,3/ c----------------------------------------------------------------------- DIMENSION NMF8(4,6),N1F8(4,6),N2F8(4,6) DATA NMF8/1,2,3,4,2,3,7,6,5,6,7,8,1,4,8,5,1,2,6,5,4,3,7,8/ c? DATA N1F8/6,7,4,2,7,6,4,7,2/,N2F8/5,1,3/ c----------------------------------------------------------------------- C CONSTANTES PHYSIQUES DATA CNU,AKER /0.09D0,0.41D0/ c----------------------------------------------------------------------- c----------------------------------------------------------------------- SQCNU=SQRT(CNU) c write(6,*)' DEBUT FPTAU : NBINC=',NBINC c write(6,*)'NES=',nes,' NPG=',npg,' idim=',idim,' np=',np c write(6,*)'NPTS=',NPTS DEUPI=1.D0 cc IF(IAXI.NE.0)DEUPI=2.D0*XPI IERC=0 ERRMAX=0.D0 KT=0 DO 2 K=1,NBEL c write(6,*)' DEBUT BCL sur K ',K NS=0 NY=0 CALL INITI(NMS,NP,0) c CALL INITI(MMS,NP,0) CALL INITI(NMY,NP,0) c write(6,*)' ELEMENT K=',K DO 20 I=1,NP J=NUM(I,K) JI=IPADT(NUM(I,K)) JM=(1-IKM)*(JI-1)+1 c JS=IPADS(NUM(I,K)) c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=', c & IPADS(J) IF(IPADS(J).NE.0)THEN NS=NS+1 NMS(NS)=I c MMS(NS)=I c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=', c & IPADS(J),' NS=',NS,' ',NMS(NS) ELSE NY=NY+1 NMY(NY)=I ENDIF DO 10 N=1,IDIM XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) VS(I,N)=U(JI,N) 10 CONTINUE c UETK(I,K)=VET(JS) c UETK(I,K)=0.D0 c ROS(I) = RO(JR) c AMUS(I)= AMU(JM) 20 CONTINUE c write(6,*)'K=',K,' MMS=',(mms(ijk),ijk=1,nS) c write(6,*)'K=',K,' XYZ=',(XYZ(1,I),I=1,np) c write(6,*)'K=',K,' XYZ=',(XYZ(2,I),I=1,np) c write(6,*)'K=',K,' NP=',NP,' NS=',NS,' NY=',NY c write(6,*)'K=',K,' NUM ',(num(i,k),i=1,np) c write(6,*)'K=',K,' NMSO',(nms(i),i=1,nS) c write(6,*)'K=',K,' NMS ',(num(nms(i),k),i=1,nS) c write(6,*)'K=',K,' NMY ',(num(NMY(I),k),I=1,NY) IF(NY.EQ.1)GO TO 2 IF(NS.EQ.1)GO TO 2 IF(NS.LE.2.AND.IDIM.EQ.3)GO TO 2 KT=KT+1 c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._ On Calcule la normale au SEG2 (NMS) et YP c._._ (La normale est la meme pour les deux pts) IF(IDIM.EQ.2.AND.NS.EQ.2)THEN Xa(1)=XYZ(1,NMS(1)) Ya(1)=XYZ(2,NMS(1)) Xa(2)=XYZ(1,NMS(2)) Ya(2)=XYZ(2,NMS(2)) c a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2)) c a=1. aa=Ya(1)-Ya(2) bb=Xa(1)-Xa(2) AN(1,NMS(1),K)=aa AN(2,NMS(1),K)=(-1.)*bb ann(1)=((aa*aa)+(bb*bb))**0.5D0 AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1) AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1) DO 241 I=1,NY Xq=XYZ(1,NMY(I)) Yq=XYZ(2,NMY(I)) YP(NMY(I),K)=abs((aa*Xq)-(bb*Yq)+(bb*Ya(1))-(aa*Xa(1)))/ann(1) 241 CONTINUE c la normale est la meme pour tous les pts (SEG2) DO 251 I=2,NS AN(1,NMS(I),K)=AN(1,NMS(1),K) AN(2,NMS(I),K)=AN(2,NMS(1),K) 251 CONTINUE DO 252 I=1,NY AN(1,NMY(I),K)=AN(1,NMS(1),K) AN(2,NMY(I),K)=AN(2,NMS(1),K) 252 CONTINUE c._._ FIN On Calcule la normale au SEG2 (NMS) et YP c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._ On Calcule la normale au SEG3 (NMS) et YP c._._ (La normale n'est pas la meme pour tous les pts) ELSEIF(IDIM.EQ.2.AND.NS.EQ.3)THEN C._._ Mais quelle est donc la face? call initi(numfa,np,0) do 2300 I=1,NS numfa(nms(I))=1 2300 continue nfa=numfa(2)+numfa(4)*2+numfa(6)*3+numfa(8)*4 c write(6,*)'NFA=',NFA Xa(1)=XYZ(1,NMF9(1,NFA)) Ya(1)=XYZ(2,NMF9(1,NFA)) Xa(2)=XYZ(1,NMF9(2,NFA)) Ya(2)=XYZ(2,NMF9(2,NFA)) Xa(3)=XYZ(1,NMF9(3,NFA)) Ya(3)=XYZ(2,NMF9(3,NFA)) a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2)) AN(1,NUM(NMF9(1,NFA),K),K)=a AN(2,NUM(NMF9(1,NFA),K),K)=-1.D0 ann(1)=(a*a+1.D0)**0.5D0 AN(1,NUM(NMF9(1,NFA),K),K)=AN(1,NUM(NMF9(1,NFA),K),K)/ann(1) AN(2,NUM(NMF9(1,NFA),K),K)=AN(2,NUM(NMF9(1,NFA),K),K)/ann(1) b=(Ya(1)-Ya(3))/(Xa(1)-Xa(3)) AN(1,NUM(NMF9(2,NFA),K),K)=b AN(2,NUM(NMF9(2,NFA),K),K)=-1.D0 ann(2)=(b*b+1.D0)**0.5D0 AN(1,NUM(NMF9(2,NFA),K),K)=AN(1,NUM(NMF9(2,NFA),K),K)/ann(2) AN(2,NUM(NMF9(2,NFA),K),K)=AN(2,NUM(NMF9(2,NFA),K),K)/ann(2) c=(Ya(3)-Ya(2))/(Xa(3)-Xa(2)) AN(1,NUM(NMF9(3,NFA),K),K)=c AN(2,NUM(NMF9(3,NFA),K),K)=-1.D0 ann(3)=(c*c+1.D0)**0.5D0 AN(1,NUM(NMF9(3,NFA),K),K)=AN(1,NUM(NMF9(3,NFA),K),K)/ann(3) AN(2,NUM(NMF9(3,NFA),K),K)=AN(2,NUM(NMF9(3,NFA),K),K)/ann(3) c write(6,*)'ann=',(ann(kk),kk=1,ns) DO 231 I=1,NS Xq=XYZ(1,N1F9(I,NFA)) Yq=XYZ(2,N1F9(I,NFA)) YP(N1F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I) Xq=XYZ(1,N2F9(I,NFA)) Yq=XYZ(2,N2F9(I,NFA)) YP(N2F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I) 231 CONTINUE c write(6,*)'NMY=',(NMY(I),I=1,NY) c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) c la normale n'est pas la meme pour tous les pts (SEG3) DO 232 I=1,NS AN(1,N1F9(I,NFA),K)=AN(1,NMS(I),K) AN(2,N1F9(I,NFA),K)=AN(2,NMS(I),K) AN(1,N2F9(I,NFA),K)=AN(1,NMS(I),K) AN(2,N2F9(I,NFA),K)=AN(2,NMS(I),K) 232 CONTINUE c._._ On Calcule la normale au SEG3 (NMS) et YP c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._ On Calcule les normales aux CUB8,PRI6,TET4 et PYR5 c._._ (La normale est la meme pour tous les pts) ELSEIF( ITYPEL.EQ.14 & .OR.ITYPEL.EQ.16 & .OR.ITYPEL.EQ.23 & .OR.ITYPEL.EQ.25 )THEN Xa(1)=XYZ(1,NMS(1)) Ya(1)=XYZ(2,NMS(1)) Za(1)=XYZ(3,NMS(1)) Xa(2)=XYZ(1,NMS(2)) Ya(2)=XYZ(2,NMS(2)) Za(2)=XYZ(3,NMS(2)) Xa(3)=XYZ(1,NMS(3)) Ya(3)=XYZ(2,NMS(3)) Za(3)=XYZ(3,NMS(3)) AN(1,NMS(1),K)= & (Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1)) AN(2,NMS(1),K)= & (Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1)) AN(3,NMS(1),K)= & (Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1)) ann(1)=(AN(1,NMS(1),K)*AN(1,NMS(1),K) & +AN(2,NMS(1),K)*AN(2,NMS(1),K) & +AN(3,NMS(1),K)*AN(3,NMS(1),K))**0.5D0 AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1) AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1) AN(3,NMS(1),K)=AN(3,NMS(1),K)/ann(1) DO 431 I=2,NS AN(1,NMS(I),K)=AN(1,NMS(1),K) AN(2,NMS(I),K)=AN(2,NMS(1),K) AN(3,NMS(I),K)=AN(3,NMS(1),K) 431 CONTINUE c write(6,*)'ann=',ann(1) DO 432 I=1,NY Xq=XYZ(1,NMY(I)) Yq=XYZ(2,NMY(I)) Zq=XYZ(3,NMY(I)) YP(NMY(I),K)=ABS( & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K) & +(Zq-Za(1))*AN(3,NMS(1),K)) c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)), c &' (Zq-Za(1))=',(Zq-Za(1)) 432 CONTINUE c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K) c write(6,*)'NMY=',(NMY(I),I=1,NY) c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) c la normale n'est pas la meme pour tous les pts (SEG3) DO 433 I=1,NY AN(1,NMY(I),K)=AN(1,NMS(1),K) AN(2,NMY(I),K)=AN(2,NMS(1),K) AN(3,NMY(I),K)=AN(3,NMS(1),K) 433 CONTINUE c._._ FIN Calcul des normales aux CUB8,PRI6,TET4 et PYR5 c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c._._ On Calcule les normales aux CU27,PR21,TE15 et PY19 c._._ (La normale est la meme pour tous les pts) ELSEIF( ITYPEL.EQ.33 & .OR.ITYPEL.EQ.34 & .OR.ITYPEL.EQ.35 & .OR.ITYPEL.EQ.36 )THEN Xa(1)=XYZ(1,NMS(1)) Ya(1)=XYZ(2,NMS(1)) Za(1)=XYZ(3,NMS(1)) Xa(2)=XYZ(1,NMS(2)) Ya(2)=XYZ(2,NMS(2)) Za(2)=XYZ(3,NMS(2)) Xa(3)=XYZ(1,NMS(3)) Ya(3)=XYZ(2,NMS(3)) Za(3)=XYZ(3,NMS(3)) Xa(4)=XYZ(1,NMS(NS)) Ya(4)=XYZ(2,NMS(NS)) Za(4)=XYZ(3,NMS(NS)) ANx1=(Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1)) ANy1=(Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1)) ANz1=(Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1)) ann(1)=(ANx1*ANx1+ANy1*ANy1+ANz1*ANz1)**0.5D0 ANx2=(Ya(2)-Ya(1))*(Za(4)-Za(1))-(Za(2)-Za(1))*(Ya(4)-Ya(1)) ANy2=(Za(2)-Za(1))*(Xa(4)-Xa(1))-(Xa(2)-Xa(1))*(Za(4)-Za(1)) ANz2=(Xa(2)-Xa(1))*(Za(4)-Za(1))-(Ya(2)-Ya(1))*(Xa(4)-Xa(1)) ann(2)=(ANx2*ANx2+ANy2*ANy2+ANz2*ANz2)**0.5D0 IF(ann(1).GT.ann(2))THEN AN(1,NMS(1),K)=ANx1/ann(1) AN(2,NMS(1),K)=ANy1/ann(1) AN(3,NMS(1),K)=ANz1/ann(1) ELSE AN(1,NMS(1),K)=ANx2/ann(2) AN(2,NMS(1),K)=ANy2/ann(2) AN(3,NMS(1),K)=ANz2/ann(2) ann(1)=ann(2) ENDIF DO 531 I=2,NS AN(1,NMS(I),K)=AN(1,NMS(1),K) AN(2,NMS(I),K)=AN(2,NMS(1),K) AN(3,NMS(I),K)=AN(3,NMS(1),K) 531 CONTINUE c write(6,*)'ann=',ann(1) DO 532 I=1,NY Xq=XYZ(1,NMY(I)) Yq=XYZ(2,NMY(I)) Zq=XYZ(3,NMY(I)) YP(NMY(I),K)=ABS( & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K) & +(Zq-Za(1))*AN(3,NMS(1),K)) c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)), c &' (Zq-Za(1))=',(Zq-Za(1)) 532 CONTINUE c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K) c write(6,*)'NMS=',(NMS(I),I=1,NS) c write(6,*)'NMY=',(NMY(I),I=1,NY) c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) c la normale n'est pas la meme pour tous les pts (SEG3) DO 533 I=1,NY AN(1,NMY(I),K)=AN(1,NMS(1),K) AN(2,NMY(I),K)=AN(2,NMS(1),K) AN(3,NMY(I),K)=AN(3,NMS(1),K) 533 CONTINUE c._._ FIN Calcul des normales aux CU27,PR21,TE15 et PY19 c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. ELSE write(6,*)'Type d elements non prevus dans FPU' write(6,*)'K=',K,' IDIM=',IDIM,' NS=',NS RETURN ENDIF c write(6,*)' FIN CALCUL NORMALES' c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY) DO 244 I=1,NP U0=0.D0 U1=0.D0 call initd(un,np,0.D0) DO 243 N=1,IDIM U0=U0+VS(I,N)*AN(N,I,K) U1=U1+(VS(I,N) - UN(I)*AN(N,I,K))**2.D0 243 CONTINUE UN(I)=U0 UT(I)=U1**0.5D0 + 1.D-10 244 CONTINUE c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY) c write(6,*)'AN=',AN(1,NMS(1),K),AN(2,NMS(1),K) c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) c write(6,*)'UN=',(UN(NMY(I)),I=1,NY) c write(6,*)'UT=',(UT(NMY(I)),I=1,NY) c write(6,*)'VS= ' c write(6,1002)(VS(I,1),I=1,NP) c write(6,1002)(VS(I,2),I=1,NP) cc CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, cc & NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) c._._ On Calcule la matrice masse diagonale maillage appuye cc CALL INITD(DG,NPTA,0.D0) cc DO 221 I=1,NP cc AI=0.D0 cc DO 220 LG=1,NPG cc PDR=PGSQ(LG)*DEUPI*RPG(LG) cc AI=AI+FN(I,LG)*PDR c220 CONTINUE cc JD=IPADA(NUM(I,K)) cc DG(JD)=DG(JD)+AI c221 CONTINUE c._._ On Calcule la matrice masse diagonale maillage appuye UETM=0.D0 YPM =0.D0 DO 225 II=1,NY I=NMY(II) IS=NMS(II) YP(I,K)=YP(I,K)+1.D-20 c write(6,*)'YP ',YP(I,K) JI=IPADT(NUM(I,K)) c JS=IPADS(NUM(IS,K)) JR=(1-IKR)*(JI-1)+1 JM=(1-IKM)*(JI-1)+1 c YPLUS=RO(JR)*YP(I,K)*UETK(I,K)/AMU(JM) YPLUS=RO(JR)*YP(I,K)*UT(I)/AMU(JM) YPLUS=YPLUS+1.D-10 c UET0=UETK(I,K) UET0=1.D-20 c write(6,*)'M=0 YPLUS=',YPLUS,UT(I),AMU(JM),UET0 w=0.5D0 DO 30 M=1,NBIT IF(YPLUS.LE.5.D0)THEN UPLUS=YPLUS UETN=UT(I)/UPLUS YPLUS=w*YPLUS + (1.D0-w)*RO(JR)*YP(I,K)*UETN/AMU(JM) ELSE c Reichardt UPLUS=1.D0/AKER*Log(1.D0+AKER*yplus)+ & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0)) UETN=UT(I)/UPLUS YPLUS=RO(JR)*YP(I,K)*UETN/AMU(JM) ENDIF resid=ABS(UETN-UET0)/UETN UET0=UETN IF(resid.LE.EPS)GO TO 31 30 CONTINUE write(6,*)'Operateur FPU, NON CONVERGENCE en 40 iterations: eps=', & eps,' residu ->',resid 31 CONTINUE c write(6,*)' Convergence en ',m,' iterations' c UETK(I,K)=UETN c VET(JS)=UETN UETM = (UETM*FLOAT(II-1)+UETN)/FLOAT(II) YPM = (YPM *FLOAT(II-1)+YP(I,K))/FLOAT(II) AMUEF(I,K)=AMU(JM) IF(YPLUS.GT.5.D0)THEN AMUEF(I,K)=(RO(JR)*UETN*UETN*YP(I,K)/UT(I))+AMU(JM) ANUT(IPADT(NUM(I,K)))=AMUEF(I,K) ENDIF 225 CONTINUE c write(6,*)' APRES 225 NY=',NY DO 226 I=1,NS JI=IPADT(NUM(I,K)) JM=(1-IKM)*(JI-1)+1 AMUEF(I,K)=AMU(JM) ANUT(IPADT(NUM(NMS(I),K)))=AMUEF(I,K) UET(IPADS(NUM(NMS(I),K)))=UETM YPS(IPADS(NUM(NMS(I),K)))=YPM 226 CONTINUE IF(NBINC.EQ.3)THEN DO 227 I=1,NS JI=IPADT(NUM(I,K)) JM=(1-IKM)*(JI-1)+1 TKE(IPADS(NUM(NMS(I),K)),1)=UETM*UETM/SQCNU TKE(IPADS(NUM(NMS(I),K)),2)=UETM*UETM*UETM/(AKER*YPM) 227 CONTINUE ENDIF c write(6,*)'NS=',NS c write(6,1001)(LKKE(I),I=1,NS) c write(6,1001)(NUM(I,K),I=1,NP) c write(6,*)'UET K=',K,(UETK(I,K),I=1,NP) c write(6,*)'AMUEF K=',K,(AMUEF(I,K),I=1,NP) c._._ On deduit aussi Tau aux sommets c._._ On deduit Nut aux points de Gauss C................................................................. C...... FIN Boucle sur les elements Indice K c write(6,*)'. FIN Boucle sur les elements Indice K ',K 2 CONTINUE c write(6,*)'UET=' c write(6,1001)(UET(i),i=1,npts) c write(6,*)' Nb d element traite KT=',KT c write(6,*)'RETOUR FPTAU' RETURN 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END