平面彈性力學有限元源程序(FORTRAN)（doc X頁）

1、平面彈性力學有限元源程序（FORTRAN）說明：1 基本控制參數信息：NG,NE,MC,NX,NB,EO,VO,DENSITY ,T(共計5個整形數，4個實型數)  NG：結構的結點總數；  NE：結構的單元總數；  MC：平面問題的類型，MC=0，為平面應力，MC=1，為平面應變；  NX：荷載工況數；  NB：支承位移數；  EO：材料彈性模量(Pa)；  VO：材料泊松比；  DENSITY ：容重(N/m3) 

2、0;T ：材料厚度(m)；2 打印輸出控制參數：NWA,NEW,NWK,NWP(4個整形數)  等于1時，輸出，否則不輸出。3 單元結點信息：(K,(IJK(I,K),I=1,3),K=1,NE) (每行4個整形數，共計NE行)  K：單元號；  IJK(1,K)：K單元I結點編號；  IJK(2,K)：K單元J結點編號；  IJK(3,K)：K單元K結點編號；4  結點坐標信息：(K,XY(1,K),XY(2,K),K=1,NG)(每行3個整形數，共計NG行) &

3、#160;K：結點號  XY(1,K)：K結點X坐標；  XY(2,K)：K結點Y坐標；5 支承信息：(K,MB(1,K),MB(2,K),ZB(K),K=1,NB)(每行3個整形數，1個實型數，共計NB行)  K：支承位移序號；  MB(1,K)：第K個支承位移所在的結點號；  MB(2,K)：第K個支承位移的坐標方向；  ZB(K)：  第K個支承位移的數值；6 按NX荷載工況數輸入荷載信息：每一荷載工況如下 : (1) NF,NP

4、,NM(3個整型數)     NF：集中荷載個數；     NP：分布荷載個數；     NM：計自重單元數；  (2) 若NF0，則輸入下面數據      K,MF(1,K),MF(2,K),ZF(K)(每行3個整形數，1個實型數，共計NF行)       K：集中荷載序號；       MF(1,K)：第K個集中荷載作用的結點號；&#

5、160;      MF(2,K)：第K個集中荷載的坐標方向；       ZF(K)：  第K個集中荷載的數值；(3) 若NP0，則輸入下面數據      K,MP(1,K),MP(2,K),ZP(K)(每行3個整形數，1個實型數，共計NP行)       K：分布荷載序號；       MP(1,K)：第K個分布荷載作用的結點號；      

6、MP(2,K)：第K個分布荷載的坐標方向；       ZP(K)：  第K個分布荷載的數值；  (4) 若NM0，則輸入下面數據      若NMNE，則表示計所有單元的自重，不需輸入計自重的單元號；      若NM<NE，則需要輸入計自重的單元號；$DEBUG         PROGRAM PLANE         IMP

7、LICIT REAL*8(A-H,O-Z),INTEGER(I-N)         ALLOCATABLE:IJK(:,:),XY(:,:),BCA(:,:),SK(:,:),STR(:,:),MB(:,:),ZB(:),B(:)         ALLOCATABLE:DELD(:,:,:),TOD(:,:),DELST(:,:,:),TOST(:,:),DELSUP(:,:),TOTSUP(:)             

8、   DIMENSION EK(6,6)        CHARACTER PN*40,FN*12        WRITE(*,'(A)') ' 本程序為計算平面問題的有限元程序'        WRITE(*,'(A)') ' 特點：(1)采用三結點三角形單元；'        W

9、RITE(*,'(A)') '  (2)采用等帶寬存貯技術；'        WRITE(*,'(A)') '    (3)采用高斯消元法解線性方程組。'        WRITE(*,'(/A)') ' 輸入計算問題名(PN):'        READ(*,'(A)') P

10、N        CALL FNAME(PN,'.DAT',FN)        WRITE(*,'(2A)') '  輸入數據文件名為：',FN        OPEN(5,FILE=FN,STATUS='OLD')        CALL FNAME(PN,'.O

11、UT',FN)        WRITE(*,'(/2A)') ' 結果輸出數據文件名為： ',FN        OPEN(6,FILE=FN,STATUS='UNKNOWN')        CALL FNAME(PN,'.OU1',FN)        WRITE(*,'(

12、/2A)') ' 參數輸出數據文件名為： ',FN        OPEN(7,FILE=FN,STATUS='UNKNOWN')                READ(5,*) NG,NE,MC,NX,NB,EO,VO,DENSITY,T                WRITE(6,120) NG,NE,MC,NX,NB 

13、               WRITE(6,130) EO,VO,DENSITY,T                READ(5,*) NWA,NWE,NWK,NWP                NT=2*NG         ALLOCATE (IJK(3,NE),XY(2,NG),BCA(7,NE),STR(3

14、,NE),MB(2,NB),ZB(NB),B(NT)ALLOCATE(DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB)                CALL CLEAR(2,NG,TOD)                CALL CLEAR(3,NE,TOST)         

15、60;     CALL CLEAR1(NB,TOTSUP)                IF (NG.EQ.0) THEN                 STOP 000                ENDIF              &#

16、160; CALL INPUT(NG,NE,NB,IJK,XY,MB,ZB)                CALL CALND(NG,NE,IJK,ND)                ALLOCATE (SK(NT,ND)                IF(MC.EQ.0) THEN         

17、      E=EO                V=VO                ELSE                E=EO/(1-VO*VO)                V=VO/(1-VO)    

18、0;           ENDIF                IP=0                NX1=NX                A1=E/(1-V*V)/4.0                A2=0.5*

19、(1-V)                CALL ABC(NG,NE,IJK,XY,BCA,NWA)                CALL CLEAR(NT,ND,SK)                DO 100 K=1,NE                CALL KE(

20、K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)                CALL SUMK(K,NE,NT,ND,IJK,EK,SK)100                CONTINUE                CALL CHECK(NT,ND,SK)110          &#

21、160;     CONTINUE                IP=IP+1                WRITE(6,'(/1X,A,I4)') "荷載工況=",IP                READ(5,*) NF,NP,NM        &#

22、160;       DO I=1,NT                B(I)=0.0                ENDDO                IF(NF.GT.0) THEN                CALL PF(NF,N

23、T,B)                ENDIF                IF(NP.GT.0) THEN                CALL PP(NP,NG,NT,T,XY,B)                ENDIF      &

24、#160;         IF(NM.GT.0) THEN                CALL PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)                ENDIF                DO I=1,NB       

25、60;       I1=2*(MB(1,I)-1)+2-MB(2,I)                DELSUP(I,IP)=-B(I1)                ENDDO                IF(NF.NE.0.OR.NP.NE.0.OR.NM.NE.0) THEN     

26、;           CALL DBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B)                CALL GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP)                CALL STRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)         

27、60;     CALL TRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST)            CALL SUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DEL SUP,        TOTSUP) CALL OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD ,DELST,TOST,DELSUP, TOTSUP)       

28、        ELSE                WRITE(*,'(2x,A,I4,A)') "荷載工況=",IP,"沒有荷載!"                WRITE(6,'(2x,A,I4,A)') "荷載工況=",IP,"沒有荷載!"    

29、0;           ENDIF                NX1=NX1-1                IF (NX1.GT.0) GOTO 110120                FORMAT(1X, '結點總數=', I3, 1X, '單元總數=', I3

30、, 1X, '問題類型=', &      & I1, 1X, '荷載工況數=', I2, 1X, '支承位移數=', I2)130                FORMAT(/1X, '彈性模量=', E10.3, 1X, '泊松比=', F5.3, 1X, '密度=', E10.3, &      & 1

31、X, '厚度=', F6.3)                END                                                SUBROUTINE GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP) 

32、              IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION SK(NT,ND),A(NT,NT),B(NT)                CALL CLEAR(NT,NT,A)             

33、  DO I=1,NT                 DO J=1,ND                  IF(I+J-1).LE.NT) THEN                  A(I,I+J-1)=SK(I,J)         

34、;         ENDIF                 ENDDO                ENDDO                IF(NWK.EQ.1.AND.NX1.EQ.NX) THEN           

35、; WRITE(7,'(A)') "結構總剛(列出右上三角元素)"                DO I=1,NT                WRITE(7,100) (A(I,J),J=1,NT)                ENDDO          

36、0;     ENDIF                IF (NWP.EQ.1) THEN                WRITE(7,'(1X,A,I3,A)') "第",NX-NX1+1,"荷載工況的荷載列陣"                DO I=1,NT  &

37、#160;             WRITE(7,'(1x,E11.4)') B(I)                ENDDO                ENDIF                DO 50 K=1,NT-1       

38、;          DO 60 I=K+1,NT                  C1=A(K,I)/A(K,K)                  DO 70 J=I,NT                  A(I,J)=A(I,J)

39、-C1*A(K,J)70                  CONTINUE                  B(I)=B(I)-C1*B(K)60                  CONTINUE50             

40、;   CONTINUE                B(NT)=B(NT)/A(NT,NT)                DO 80 I=NT-1,1,-1                 DO 90 J=I+1,NT              

41、0;  B(I)=B(I)-A(I,J)*B(J)90                 CONTINUE                 B(I)=B(I)/A(I,I)80                CONTINUE100          &

42、#160;     FORMAT(1x, 4000E10.3)                END                SUBROUTINE CHECK(NT,ND,SK)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)      

43、0;         DIMENSION SK(NT,ND)                M=0                DO I=1,NT                 IF(SK(I,1).LE.0.0) THEN          &

44、#160;      M=M+1                 ENDIF                ENDDO                IF (M.GT.0) THEN                WRI

45、TE(*,*) "總剛矩陣的對角元素為負值!"                STOP 2222                ENDIF                END                SUBROUTINE STRESS(NE,NT,V,

46、A1,A2,IJK,BCA,B,STR)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION D1(6),D2(3),B(NT),S(3,6),IJK(3,NE),BCA(7,NE),STR(3,NE)                CALL CLEAR(

47、3,NE,STR)                DO 10 K=1,NE                 DO I=1,3                 II=IJK(I,K)                 D1(

48、2*(I-1)+1)=B(2*(II-1)+1)                 D1(2*(I-1)+2)=B(2*(II-1)+2)                 ENDDO                 CALL CLEAR(3,6,S)       

49、;          C1=2*A1/BCA(7,K)                 DO 20 I=1,3                  S(1,2*(I-1)+1)=C1*BCA(I,K)                  S(1

50、,2*(I-1)+2)=C1*V*BCA(I+3,K)                  S(2,2*(I-1)+1)=C1*V*BCA(I,K)                  S(2,2*(I-1)+2)=C1*BCA(I+3,K)                  S(3,2*(I-1)+1

51、)=C1*A2*BCA(I+3,K)                  S(3,2*(I-1)+2)=C1*A2*BCA(I,K)20                 CONTINUE                 CALL MATMUL(3,6,1,S,D1,D2)    &

52、#160;            DO I=1,3                 STR(I,K)=D2(I)                 ENDDO10                CONTINUE     

53、          END                SUBROUTINE OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)            &

54、#160;   DIMENSION MB(2,NB),DELD(2,NG,NX),TOD(2,NG)                DIMENSION DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB)                CHARACTER VE*6                WRI

55、TE(6,20) IP                WRITE(6,30)                DO I=1,NG                WRITE(6,40) I,DELD(1,I,IP),TOD(1,I),DELD(2,I,IP),TOD(2,I)            &#

56、160;   ENDDO                WRITE(6,50)                DO J=1,NE                WRITE(6,60) J,(DELST(L,J,IP),TOST(L,J),L=1,3)              &

57、#160; ENDDO                WRITE(6,70)                DO J=1,NB                IF (MB(2,J).EQ.1) THEN                VE='x方向' 

58、              ELSE                VE='Y方向'                ENDIF                WRITE(6,80) MB(1,J),VE,DELSUP(J,IP),TOTSUP(J)   

59、             ENDDO20                FORMAT(/1X, '荷載工況=', I3, '的計算結果')30                FORMAT(/, 1X, '結點號', 5X, 'X位移增量', 5X, 'X累計位移', 5x, & 

60、    & 'Y位移增量', 5X, 'Y累計位移')40                FORMAT(1X, I3, 6X, F10.4, 4x, F10.4, 4X, F10.4, 4x, F10.4)50                FORMAT(/, 1X, '單元號', 6x, 'X應力增量', 6x, 'X累計應力

61、', 6x, &      & 'Y應力增量', 6x, 'Y累計應力', 6x, '剪應力增量', 6x, &      & '累計剪應力')60                FORMAT(2X, I3, 6X, E10.3, 5X, E10.3, 5X, E10.3, 5X, E10.3, 6X, E10.3, &

62、60;     & 6X, E10.3)70                FORMAT(/, 1X, '支座結點號', 6x, '支反力方向', 6x, '支反力增量', &      & 6x, '支反力累計量')80                FORMAT(5X, I3,

63、 12X, A, 6x, E10.3, 8X, E10.3)            END                                                          &#

64、160;                SUBROUTINE MATMUL(M,N,L,A,B,C)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)        DIMENSION A(M,N),B(N,L),C(M,L)        DO 100 I=1,

65、M         DO 100 J=1,L          C(I,J)=0.0          DO 100 K=1,N100     C(I,J)=C(I,J)+A(I,K)*B(K,J)        END              &

66、#160; SUBROUTINE CLEAR(M,N,A)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION A(M,N)                DO I=1,M           

67、    DO J=1,N                A(I,J)=0.0                ENDDO                ENDDO                END       

68、60;       SUBROUTINE ABC(NG,NE,IJK,XY,BCA,NWA)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION X(2,5),XY(2,NG),IJK(3,NE),B(7),BCA(7,NE)           

69、;     DO 10 I=1,NE                 DO 20 K=1,3                 DO 20 J=1,2                 X(J,K)=XY(J,IJK(K,I)20       

70、          CONTINUE                 DO 30 J=1,2                 X(J,4)=X(J,1)                 X(J,5)=X(J,2)30   

71、;              CONTINUE                 DO 40 K=1,3                 B(K)=X(2,K+1)-X(2,K+2)                 B(K

72、+3)=X(1,K+2)-X(1,K+1)40                 CONTINUE                 B(7)=(B(1)*B(5)-B(4)*B(2)*0.5                 IF (B(7).LE.0.0) THEN    

73、0;            WRITE(*,*) "單元號=",I                 WRITE(*,*) "單元的結點I,J AND K編號：",(IJK(K,I),K=1,3)                 STOP 1111      &

74、#160;          ELSE                 DO J=1,7                 BCA(J,I)=B(J)                 ENDDO     

75、60;          ENDIF10                CONTINUE                IF(NWA.EQ.1) THEN                WRITE(7,100)           

76、60;   DO I=1,NE                WRITE(7,110) I,(BCA(J,I),J=1,7)                ENDDO                ENDIF100                 F

77、ORMAT(1X, '單元號', 3x, 'bi', 7x, 'bj', 7x, 'bk', 7x, 'ci', 7x, &      & 'cj', 7x, 'ck', 7x, '面積A')110                 FORMAT(2X, I3, 3X, F6.3, 3X, F6.3, 3X, F6.

78、3, 3X, F6.3, 3X, F6.3, &      & 3X, F6.3, 3X, F6.3)             END                 SUBROUTINE SUMK(K,NE,NT,ND,IJK,EK,SK)             IMPLICIT

79、 REAL*8(A-H,O-Z),INTEGER(I-N)                 DIMENSION        IJ(3),EK(6,6),IJK(3,NE),SK(NT,ND)                 DO I=1,3               

80、0; IJ(I)=IJK(I,K)                 ENDDO                 DO 10 I=1,3                  DO 20 J=1,3             

81、    IF (IJ(I).GT.IJ(J) GOTO 20                  M=2*(IJ(I)-1)+1                  N=2*(IJ(J)-1)+1-(2*(IJ(I)-1)+1)+1                  MO=

82、2*(I-1)+1                  NO=2*(J-1)+1                  SK(M,N)=SK(M,N)+EK(MO,NO)                  SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)     

83、;             SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1)                  IF (IJ(I).EQ.IJ(J) GOTO 20                  SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO)20     

84、;             CONTINUE10             CONTINUE                 END                SUBROUTINE KE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE) 

85、60;             IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION  IJ(3),IJK(3,NE),BCA(7,NE),EK(6,6)                CALL CLEAR(6,6,EK)         

86、;       DO I=1,3                IJ(I)=IJK(I,K)                ENDDO                DO 10 I=1,3                 DO 2

87、0 J=I,3                 IO=2*(I-1)                 JO=2*(J-1)                 EK(IO+1,JO+1)=BCA(I,K)*BCA(J,K)+A2*BCA(I+3,K)*BCA(J+3,K)    &#

88、160;            EK(IO+1,JO+2)=V*BCA(I,K)*BCA(J+3,K)+A2*BCA(I+3,K)*BCA(J,K)                 EK(IO+2,JO+2)=BCA(I+3,K)*BCA(J+3,K)+A2*BCA(I,K)*BCA(J,K)                 I

89、F (I.EQ.J) GOTO 20                 EK(IO+2,JO+1)=V*BCA(I+3,K)*BCA(J,K)+A2*BCA(I,K)*BCA(J+3,K)20                 CONTINUE10                CONTINUE   

90、            DO I=2,6                 DO J=1,I-1                 EK(I,J)=EK(J,I)                 ENDDO     

91、;           ENDDO                DO I=1,6                 DO J=1,6                 EK(I,J)=EK(I,J)*A1*T/BCA(7,K)      &#

92、160;          ENDDO                ENDDO                IF(NWE.EQ.1) THEN             WRITE(7,40)        K, "號單元的單剛"

93、0;               DO I=1,6                 WRITE(7,50) (EK(I,J),J=1,6)                ENDDO                ENDIF40      &

94、#160;         FORMAT(1X, I3, A)50                FORMAT(1X, 6(E10.4, 4X)                END                SUBROUTINE PF(NF,NT,B)       

95、60;       IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION MF(2,NF),ZF(NF),B(NT)                DO I=1,NF                READ(5,*) K,MF(1,I),MF(2,I),ZF

96、(I)                 ENDDO                WRITE(6,'(/A,I3)') "集中荷載個數NF=",NF                WRITE(6,20)                DO I=

97、1,NF                WRITE(6,30) MF(1,I),MF(2,I),ZF(I)                ENDDO                DO 10 I=1,NF                 II=2*(MF(1,

98、I)-1)-MF(2,I)+2                 B(II)=B(II)+ZF(I)10                CONTINUE20                FORMAT(1X, '結點號', 4x, '方向', 4x, '荷載值')30   

99、             FORMAT(2X, I3, 5X, I3, 4X, E10.4)                END                        SUBROUTINE FNAME(PN,FN2,FN)        CHARAC

100、TER PN*40,FN2*4,FN*12!       去掉PN中前面的空格        DO 10 I=1,40         IF(PN(I:I).EQ.' ') GOTO 10         IP=I         GOTO 2010      C

101、ONTINUE20      CONTINUE        FN(1:8)=PN(IP:IP+7)!       去掉FN中后面的空格        DO 30 I=8,1,-1         IF(FN(I:I).EQ.' ') GOTO 30        

102、0;IL=I         GOTO 4030      CONTINUE40      CONTINUE!       生成文件名FN=PN+FN2        FN(IL+1:IL+4)=FN2(1:4)        END         

103、;       SUBROUTINE PP(NP,NG,NT,T,XY,B)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION MP(2,NP),ZP(NP),XY(2,NG),B(NT)                DO

104、 I=1,NP                READ(5,*) K,MP(1,I),MP(2,I),ZP(I)                ENDDO                WRITE(6,'(/A,I3)') "分布荷載數NP=",NP          &#

105、160;     WRITE(6,20)                DO I=1,NP                WRITE(6,30) I,MP(1,I),MP(2,I),ZP(I)                ENDDO              

106、0; DO 10 I=1,NP                 II=2*(MP(1,I)-1)                 JJ=2*(MP(2,I)-1)                 PX=0.5*ZP(I)*T*(XY(2,MP(1,I)-XY(2,MP(2,I)  

107、0;              PY=0.5*ZP(I)*T*(XY(1,MP(2,I)-XY(1,MP(1,I)                 B(II+1)=B(II+1)+PX                 B(II+2)=B(II+2)+PY         

108、;        B(JJ+1)=B(JJ+1)+PX                 B(JJ+2)=B(JJ+2)+PY10                CONTINUE20                FORMAT(1X, '序號', 4x, '首結點號&

109、#39;, 4x, '末結點號', 7x, '荷載值')30                FORMAT(1X, I3, 8x, I3, 9X, I3, 7X, F10.3)                END                           

110、    SUBROUTINE PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)                DIMENSION IJK(3,NE),BCA(7,NE),MG(NM),B(NT)             

111、60; IF (NM.GE.NE) THEN                NM=NE                DO I=1,NM                MG(I)=I                ENDDO       

112、        ELSE                READ(5,*) (MG(I),I=1,NM)                ENDIF                WRITE(6,'(/A,I3)') "計自重單元數NM=",NM     

113、           WRITE(6,'(A)') "計自重的單元號為："                K1=0                K2=0                DO WHILE(NM-K2).GE.10)      &

114、#160;         WRITE(6,100) (MG(K2+J),J=1,10)                K1=K1+1                K2=10*K1                ENDDO             

115、   WRITE(6,100) (MG(K2+J),J=1,NM-K2)                DO 10 I=1,NM                 II=2*(IJK(1,MG(I)-1)                 JJ=2*(IJK(2,MG(I)-1)   

116、60;            KK=2*(IJK(3,MG(I)-1)                 PE1=-DENSITY*T*BCA(7,MG(I)/3.0                 B(II+2)=B(II+2)+PE1             

117、;    B(JJ+2)=B(JJ+2)+PE1                 B(KK+2)=B(KK+2)+PE110                CONTINUE100                FORMAT(10(1X,I3)         

118、       END                                SUBROUTINE INPUT(NG,NE,NB,IJK,XY,MB,ZB)                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)       

119、60;       DIMENSION IJK(3,NE),XY(2,NG),MB(2,NB),ZB(NB)                DO I=1,NE                READ(5,*) K,(IJK(J,I),J=1,3)                 ENDDO     

120、           DO I=1,NG                READ(5,*) K,XY(1,I),XY(2,I)                ENDDO                DO I=1,NB            

121、0;    READ(5,*) K,MB(1,I),MB(2,I),ZB(I)                ENDDO                WRITE(6,15)                DO K=1,NE              

122、0; WRITE(6,20) K,(IJK(J,K),J=1,3)                 ENDDO                WRITE(6,30)                DO K=1,NG                WRITE(6,35) K,XY(

123、1,K),XY(2,K)                ENDDO                WRITE(6,40)                DO K=1,NB                 WRITE(6,45) K,MB(1,K),MB(2,K),ZB(K)                ENDDO15                FORMAT(/1X, '單元號', 4X, 'I結點號',