[求助] 有限元调试程序

program main
 common NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION X(100),Y(100),NCO(200),PROP(100),IB(60),TK(200,200),&
 &         AL(200),FORC(100),REAC(200),ELST(4,4),V(20)
 CHARACTER(len=20) FILE1,FILE2

 !初始化程序参数
 NRMX=200
 NCMX=20
 NDF=2
 NNE=2
 NDFEL=NDF*NNE

 IN=5
 IO=6

 !打开输入输出文件
 READ(*,*)FILE1,FILE2
 OPEN(UNIT=IO,FILE=FILE1,FORM='FORMATTED',STATUS='UNKNOWN')
 OPEN(UNIT=IN,FILE=FILE2,FORM='FORMATTED',STATUS='OLD')

 !数据输入
 CALL INPUT(X,Y,NCO,PROP,AL,IB,REAC)

 !组装总体刚度矩阵
 CALL ASSEM(X,Y,NCO,PROP,TK,ELST,AL)

 !引入边界条件
 CALL BOUND(TK,AL,REAC,IB)

 !总体方程求解
 CALL SLBSI(TK,AL,V,N,MS,NRMX,NCMX)

 !计算单元内力与结点合力
 CALL FORCE(NCO,PROP,FORC,REAC,X,Y,AL)

 !输出
 CALL OUTPT(AL,FORC,REAC)

 STOP
 END
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 SUBROUTINE INPUT(X,Y,NCO,PROP,AL,IB,REAC)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION X(1),Y(1),PROP(1),AL(1),IB(1),REAC(1),W(3),IC(2),NCO(1)

 WRITE(IO,20)
 20  FORMAT('',70('*'))
 !读入基本常数
 READ(IN,*)NN,NE,NLN,NBN,E
 WRITE(IO,21)NN,NE,NLN,NBN,E
 21  FORMAT(//'INTERNAL DATA'//'NUMBER OF NODES:',I5/&
            & 'NUMBER OF ELEMENTS:',I5/'NUMBER OF LOADED NODES:',&
                    & I5/'NUMBER OF SUPPORT NODES:',I5/'MODULUS OF ELASTICITY:',&
                    & F15.0//'NODAL COORDINATES'/7X,'NODE',6X,'X',9X,'Y' )
 !输入结点x和y坐标值
 READ(IN,*)(I,X(I),Y(I),I=1,NN)
 WRITE(IO,2)(I,X(I),Y(I),I=1,NN)
 2 FORMAT(I10,2F10.2)
 !输入结点号和截面几何性质
 WRITE(IO,22)
 22 FORMAT(/'ELEMENT CONNECTIVITY AND PROPERTIES'/4X,'ELEMENT',3X,&
 & 'START NODE END NODE',5X,'AREA')

 DO I=1,NE
 READ(IN,*)NUM,IC(1),IC(2),PROP(NUM)
 WRITE(IO,34)NUM,IC(1),IC(2),PROP(NUM)
 N1=NNE*(NUM-1)
 NCO(N1+1)=IC(1)
 NCO(N1+2)=IC(2)
 END DO
 34 FORMAT(3I10,F15.5)

 N=NN*NDF
 DO I=1,N
 REAC(I)=0
 AL(I)=0
 END DO
 WRITE(IO,23)
 23 FORMAT(/'NODAL LOADS'/7X,'NODE',5X,'PX',8X,'PY')
 DO I=1,NLN
 READ(IN,*)NUM,(W(K),K=1,NDF)
 WRITE(IO,2)NUM,(W(K),K=1,NDF)
   DO K=1,NDF
   L=NDF*(NUM-1)+K
   AL(L)=W(K)
   END DO
 END DO

 !READ BOUNDARY DATA.STORE UNKNOWN STATUS INDICATORS
 WRITE(IO,24)
 24 FORMAT(/'BOUNDARY CONDITION DATA'/23X,'STATUS',14X,'PRESCRIBED VALUES'/15X,'(0:PRESCRIBED,1:FREE)'/7X,'NODE',8X,'U',9X,&
            & 'V',16X,'U',9X,'V')
 DO I=1,NBN
    READ(IN,*)NUM,(IC(K),K=1,NDF),(W(K),K=1,NDF)
    WRITE(IO,9) NUM,(IC(K),K=1,NDF),(W(K),K=1,NDF)
    L1=(NDF+1)*(I-1)+1
    L2=NDF*(NUM-1)
    IB(L1)=NUM
    DO K=1,NDF
    N1=L1+K
    N2=L2+K
    IB(N1)=IC(K)
    REAC(N2)=W(K)
    END DO
 END DO
 9 FORMAT(3I10,10X,2F10.4)
 RETURN 
 END
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 SUBROUTINE ASSEM(X,Y,NCO,PROP,TK,ELST,AL)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION X(1),Y(1),NCO(1),TK(200,20),ELST(NDFEL,NDFEL),PROP(1),AL(1)
 !COMPUTE HALF BAND WIDTH AND STORE IN MS
 N1=NNE-1
 MS=0
 DO I=1,NE
    L1=NNE*(I-1)
    DO J=1,N1
       L2=L1+J
       J1=J+1
           DO K=J1,NNE
           L3=L1+K
           L=ABS(NCO(L2)-NCO(L3))
           IF((MS-L).LT.0)THEN
           MS=L
           ENDIF
           ENDDO
    ENDDO
 ENDDO
 MS=NDF*(MS+1)
 !CLEAR THE TOTAL STIFFENSS MATRIX
 DO I=1,N
    DO J=1,MS
    TK(I,J)=0
    END DO
 END DO

 DO NEL=1,NE
    CALL STIFF(NEL,X,Y,NCO,PROP,ELST,AL)
    CALL ELASS(NEL,NCO,TK,ELST)
 END DO
 RETURN
 END
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 SUBROUTINE STIFF(NEL,X,Y,NCO,PROP,ELST,AL)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION X(1),Y(1),NCO(1),PROP(1),ELST(4,4),AL(1)
 L=NNE*(NEL-1)
 N1=NCO(L+1)
 N2=NCO(L+2)

 D=SQRT((X(N2)-X(N1))**2+(Y(N2)-Y(N1))**2)
 CO=(X(N2)-X(N1))/D
 SI=(Y(N2)-Y(N1))/D

 COFE=E*PROP(NEL)/D
 ELST(1,1)=COFE*CO*CO
 ELST(1,2)=COFE*CO*SI
 ELST(2,2)=COFE*SI*SI
 DO I=1,2
    DO J=1,2
    K1=I+NDF
    K2=J+NDF
    ELST(K1,K2)=ELST(I,J)
    ELST(I,K2)=-ELST(I,J)
    ENDDO
 ENDDO
 ELST(2,3)=-ELST(1,2)
 RETURN
 END

 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ELASS(NEL,NCO,TM,ELMAT)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION NCO(1),TM(200,20),ELMAT(NDFEL,NDFEL)
L1=NNE*(NEL-1)
DO I=1,NNE
   L2=L1+I
   N1=NCO(L2)
   I1=NDF*(I-1)
   J1=NDF*(N1-1)
   DO J=1,NNE
      L2=L1+J
      N2=NCO(L2)
      I2=NDF*(J-1)
      J2=NDF*(N2-1)
      DO K=1,NDF
          KI=1
          IF((N1-N2).EQ.0)THEN
          KI=K
          ENDIF
          IF((N1-N2).LE.0)THEN
          KR=J1+K
          IC=J2-KR+1
          K1=I1+K
          ELSE
          KR=J2+K
          IC=J1-KR+1
          K2=I2+K
          ENDIF
          DO L=KI,NDF
             KC=IC+L
                 IF((N1-N2).LE.0)THEN
                 K2=I2+L
                 ELSE
                 K1=I1+L
                 ENDIF
                 TM(KR,KC)=TM(KR,KC)+ELMAT(K1,K2)
          ENDDO
          ENDDO
  ENDDO
ENDDO
RETURN
END  
 SUBROUTINE BOUND(TK,AL,REAC,IB)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION AL(1),IB(1),REAC(1),TK(200,20)
 DO L=1,NBN
    L1=(NDF+1)*(L-1)+1
    NO=IB(L1)
    K1=NDF*(NO-1)
    DO I=1,NDF
       L2=L1+I
           IF(IB(L2).EQ.0)THEN
           KR=K1+I
           DO J=2,MS
           KV=KR+J-1
           IF((N-KV).GE.0)THEN
           AL(KV)=AL(KV)-TK(KR,J)*REAC(KR)
           TK(KR,J)=0
           ENDIF
           KV=KR-J+1
           IF(KV.GT.0)THEN
           AL(KV)=AL(KV)-TK(KV,J)*REAC(KR)
           TK(KV,J)=0
           ENDIF
    ENDDO
    TK(KR,1)=1
    AL(KR)=REAC(KR)
    ENDIF
    ENDDO
 ENDDO
 RETURN
 END


 SUBROUTINE SLBSI(A,B,D,N,MS,NX,MX)
 DIMENSION A(NX,MX),B(NX),D(MX)
 N1=N-1
 DO K=1,N1
    C=A(K,1)
    K1=K+1
    IF((ABS(C)-0.00000001).LE.0)THEN
    WRITE(6,2)K
 2  FORMAT('*****SINGULARITY IN ROW',I5)
    STOP
    ELSE
    N1=K1+MS-2
    L=MIN(N1,N)
    DO J=2,MS
    D(J)=A(K,J)
    ENDDO
    DO J=K1,L
    A(K,K2)=A(K,K2)/C
    ENDDO
    B(K)=B(K)/C

    DO I=K1,L
    K2=I-K1+2
    C=D(K2)
    DO J=I,L
    K2=J-I+1
    K3=J-K+1
    A(I,K2)=A(I,K2)-C*A(K,K3)
    ENDDO
    B(I)=B(I)-C*B(K)
    ENDDO
 ENDIF
 ENDDO

 IF((ABS(C)-0.0000001).LE.0)THEN
 WRITE(6,7)K
 7     FORMAT('*****SINGULARITY IN ROW',I5)
 STOP
 ELSE
 B(N)=B(N)/A(N,1)
 DO I=1,N1
    K=N-1
    K1=K+1
    N1=K1+MS-2
    L=MIN(N1,N)
    DO J=K1,L
    K2=J-K+1
    B(K)=B(K)-A(K,K2)*B(J)
    ENDDO
 ENDDO
 ENDIF
 RETURN
 END


 SUBROUTINE FORCE(NCO,PROP,FORC,REAC,X,Y,AL)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION NCO(1),PROP(1),FORC(1),REAC(1),X(1),Y(1),AL(1)
 DO  I=1,N
    REAC(I)=0
 ENDDO
 DO NEL=1,NE
    L=NNE*(NEL-1)
    N1=NCO(L+1)
    N2=NCO(L+2)
    K1=NDF*(N1-1)
    K2=NDF*(N2-1)
    D=SQRT((X(N2)-X(N1))**2+(Y(N2)-Y(N1))**2)
    CO=(X(N2)-X(N1))/D
    SI=(Y(N2)-Y(N1))/D
    COEF=E*PROP(NEL)/D
    FORC(NEL)=COEF*((AL(K2+1)-AL(K1+1))*CO+(AL(K2+2)-AL(K1+2))*SI)
    REAC(K1+1)=REAC(K1+1)-FORC(NEL)*CO
    REAC(K1+2)=REAC(K1+2)-FORC(NEL)*SI
    REAC(K2+1)=REAC(K2+1)+FORC(NEL)*CO
    REAC(K2+2)=REAC(K2+2)+FORC(NEL)*SI
 ENDDO
 RETURN
 END


 SUBROUTINE OUTPT(AL,FORC,REAC)
 COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
 DIMENSION AL(1),REAC(1),FORC(1)
 WRITE(IO,1)
 1     FORMAT(//1X,70('*')//'RESULTS'//'NODAL DISPLACEMENTS'/7X,'NODE' &
      & ,11X,'U',14X,'V')
         DO I=1,NN
          K1=NDF*(I-1)+1
          K2=K1+NDF-1
          WRITE(IO,2)I,(AL(J),J=K1,K2)
          ENDDO
 2     FORMAT(I10,6F15.4)
 WRITE(IO,3)
 3     FORMAT(/'NODAL REACTIONS'/7X,'NODE',10X,'PX',13X,'PY')
      DO I=1,NN
          K1=NDF*(I-1)+1
          K2=K1+NDF-1
          WRITE(IO,2)I,(REAC(J),J=K1,K2)
          ENDDO
          WRITE(IO,4)
 4     FORMAT(/'MEMBER FORCES'/6X,'MEMBER AXIAL FORCE')
 DO I=1,NE
 WRITE(IO,2)I,FORC(I)
 ENDDO
 WRITE(IO,5)
 5     FORMAT(//1X,70('*'))
 RETURN
 END