编辑代码

     PROGRAM MAIN
      INTEGER  RR,AA,GG,EP,RP,T
      REAL KK
      CHARACTER*12 FNAME,FNAOT
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A4/KI2,KI3,MI2(20),MI3(20),FI2(20,3),FI3(20,6)
      COMMON/A5/AA(200),KK(6000),NX,NH
      COMMON/A6/FF(60),UV(20,6)
      WRITE(0,101)
 101  FORMAT(///5X,'PLEASE INPUT DATA FILENAME:')
      READ(*,100) FNAME
 100   FORMAT(8A)
      WRITE(0,102)
 102  FORMAT(///5X,'PLEASE INPUT OUTPUT DATA FILENAME:')
      READ(*,103) FNAOT
 103  FORMAT(8A)
      OPEN(5,FILE=FNAME,STATUS='OLD')
      OPEN(6,FILE=FNAOT,STATUS='new')
      CALL INPUT
	CALL FORRR
      CALL FORAA
      CALL FORKK
      CALL DIVKK
      DO 200 N9=1,MX
      CALL FORFF(N9)
      CALL FORDIS
      CALL FORFM
      CALL FORRC
 200  CONTINUE
      STOP
      END

      SUBROUTINE INPUT
      INTEGER  EP,RP
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      READ(5,*)NE,NP,NR,MT,MX
      WRITE(6,50)NE,NP,NR,MT,MX
      WRITE(0,50)NE,NP,NR,MT,MX
50    FORMAT(2X,'NUMBER OF ELEMENT NE=',I2/2X,'NUMBER OF NODE NP=',I2/2X,
     #'NUMBER OF SUPPOT NODE NR=',I2/2X,
     #'NUMBER OF MATERIAL SET MT=',I2/2X,
     #'NUMBER OF WORKING CASE MX=',I2)
      READ(5,*)((AE(I,J),J=1,3),I=1,MT)
      WRITE(6,221)
221   FORMAT(2X,'N0. OF SET',2X,'MODULUS OF ELASTICITY',2X,'CROSS-AREA',2X,'MOMEMT OF INERTIA')
      DO 15 I=1,MT
      WRITE(6,222) I,(AE(I,J),J=1,3)
222   FORMAT(2X,I4,10X,3F16.5)
15    CONTINUE
      READ(5,*)((EP(I,J),J=1,3),I=1,NE)
      WRITE(6,229)
229   FORMAT(10X,'DATA OF ELEMENT')
      DO 17 I=1,NE
      WRITE (6,193)I,(EP(I,J),J=1,3)
193   FORMAT(2X,I4,5X,3I8)
17    CONTINUE
      READ(5,*)((RP(I,J),J=1,4),I=1,NR)
      WRITE(6,224)
224   FORMAT(10X,'DATA OF SUPPORT')
      WRITE(6,96)((RP(I,J),J=1,4),I=1,NR)
      WRITE(0,96)((RP(I,J),J=1,4),I=1,NR)
96    FORMAT(2X,4I6)
      READ(5,*)((XZ(I,J),J=1,2),I=1,NP)
      WRITE(6,226)
226   FORMAT(10X,'NODAL COORDINATE'/5X,'N0. NODE',8X,'X-X'
     #,8X,'Z-Z')
      DO 227 I=1,NP
      WRITE(6,194)I,(XZ(I,J),J=1,2)
194   FORMAT(5X,I4,5X,2F12.6)
227   CONTINUE
      RETURN
      END

      SUBROUTINE FORRR
      INTEGER  RR,AA,GG,EP,RP,T
      REAL KK
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A5/AA(200),KK(6000),NX,NH
      DO 483 I=1,NP
      DO 483 J=1,3
      RR(I,J)=1
483   CONTINUE
      DO 470 I=1,NR
      IZ1=RP(I,1)
      IW1=RP(I,2)
      IW2=RP(I,3)
      IW3=RP(I,4)
      IF(IW1.EQ.0) RR(IZ1,1)=0
      IF(IW2.EQ.0) RR(IZ1,2)=0
      IF(IW3.EQ.0) RR(IZ1,3)=0
470   CONTINUE
      NX=0
      DO 520 I=1,NP
      DO 520 J=1,3
      IF(RR(I,J).EQ.0) GOTO 520
      NX=NX+1
      RR(I,J)=NX
520   CONTINUE
      WRITE(6,540)NX
540   FORMAT(2X,'DEGREE NX=',I4)
      RETURN
      END

      SUBROUTINE FORAA
      INTEGER  RR,AA,GG,EP,RP,T
      REAL KK
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A5/AA(200),KK(6000),NX,NH
      DO 630 I=1,NX
      AA(I)=0
630   CONTINUE
      DO 685 I9=1,NE
      ME(1)=EP(I9,1)
	ME(2)=EP(I9,2)
	T=EP(I9,3)	
      CALL GGT
      DO 685 J3=1,6
      DO 685 J1=1,6
      IF((GG(J3).GT.0).AND.(GG(J1).GT.0))GOTO 671
      GOTO 685
671   CONTINUE
      IGG=GG(J3)
      IG=IGG-GG(J1)
      IF(AA(IGG).LE.IG) AA(IGG)=IG
685   CONTINUE
      AA(1)=1
      DO 711 I=2,NX
      AA(I)=AA(I)+AA(I-1)+1
711   CONTINUE
      NH=AA(NX)
      WRITE(6,9)NH
9     FORMAT(2X,'NUMBER [K] ELEMENT NH=',I4)
      RETURN
      END

      SUBROUTINE FORKK
      INTEGER  RR,AA,GG,EP,RP,T
      REAL LT,KE,KK
	COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A5/AA(200),KK(6000),NX,NH
      DIMENSION LT(3,3),KE(6,6)
      DO 761 I=1,NH
761   KK(I)=0
      DO 855 I9=1,NE
	ME(1)=EP(I9,1)
 	ME(2)=EP(I9,2)
 	T=EP(I9,3)	                    	      
      CALL GGT
      CALL MIC(H,LT,E,A,F,X,Z)
      CALL KET(X,Z,H,E,A,F,KE)
	WRITE(6,11)((KE(I,J),J=1,6),I=1,6)
 11	FORMAT(2X,'KE',6F12.3)
      DO 855 I=1,6
      DO 855 J=1,6
      A1=KE(I,J)
      IF((GG(J).LE.GG(I)).AND.(GG(J).GT.0)) GOTO 841
      GOTO 855
841   IG=GG(I)
      IG1=GG(J)
      IAB=AA(IG)
      IG8=IAB-IG+IG1
      KK(IG8)=KK(IG8)+A1
855   CONTINUE
      CONTINUE
      RETURN
      END

      SUBROUTINE DIVKK
      INTEGER  AA
      REAL KK
      COMMON/A5/AA(200),KK(6000),NX,NH
      DO 825 I=2,NX
      L=AA(I-1)+I-AA(I)+1
      K=I-1
      DO 895 J=L,K
      IF(J-L)900,895,900
900   JP=AA(I)-I+J
      M=AA(J-1)+J-AA(J)+1
      IF(L.GT.M) M=L
      MP=J-1
      DO 850 IP=M,MP
      JJ=AA(I)-I+IP
      JK=AA(J)-J+IP
      KK(JP)=KK(JP)-KK(JJ)*KK(JK)
850   CONTINUE
895   CONTINUE
      DO  906 IP=L,K
      JI=AA(I)-I+IP
      JL=AA(IP)
      KK(JI)=KK(JI)/KK(JL)
      JN=AA(I)
      KK(JN)=KK(JN)-KK(JI)*KK(JI)*KK(JL)
906   CONTINUE
825   CONTINUE
      RETURN
      END

      SUBROUTINE  FORFF(N9)
      INTEGER  RR,AA,GG,EP,RP,T
      REAL LT,KK
	COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A4/KI2,KI3,MI2(20),MI3(20),FI2(20,3),FI3(20,6)
      COMMON/A5/AA(200),KK(6000),NX,NH
      COMMON/A6/FF(60),UV(20,6)
      DIMENSION LT(3,3)
      DO 663 J=1,NX
663   FF(J)=0.0
      WRITE(6,672) N9
672   FORMAT(10X,'WORKING CASE=',I6)
      READ(5,*)KI2
      WRITE(6,675)KI2
675   FORMAT(2X,'NUMBERS OF NODAL LOAD=',I4)
      IF(KI2.EQ.0) GOTO 690
      READ(5,*)(MI2(I),I=1,KI2)
      WRITE(6,677)(MI2(I),I=1,KI2)
677   FORMAT(2X,'NODAL NUMBER OF LOAD=',10I4)
      WRITE(6,676)
676   FORMAT(10X,'VALUE OF NODAL LOAD'/5X,'X-FORCE',8X,
     #'Z-FORCE',8X,'Y-MOMEMT')
      READ(5,*)((FI2(I,J),J=1,3),I=1,KI2)
      WRITE(6,679)((FI2(I,J),J=1,3),I=1,KI2)
679   FORMAT(2X,3F12.4)
      DO 681 I=1,KI2
      M7=MI2(I)
      DO 684 J=1,3
      K=RR(M7,J)
      IF(K.GT.0) FF(K)=FF(K)+FI2(I,J)
684   CONTINUE
681   CONTINUE
690   READ(5,*)KI3
      WRITE(6,692)KI3
692   FORMAT(2X,'NUMBER OF ELEMENTAL LOAD=',I4)
      IF(KI3.EQ.0) GOTO 710
      READ(5,*)(MI3(I),I=1,KI3)
      WRITE(6,693)(MI3(I),I=1,KI3)
693   FORMAT(2X,'ELEMENTAL NUMBER OF LOAD',10I4)
      READ(5,*)((FI3(I,J),J=1,6),I=1,KI3)
      WRITE(6,691)
691   FORMAT(//30X,'VAIUE OF FIXING FORCE'/2X,'ELEMEMT',1X,'I-NX'
     #,4X,'I-QZ',7X,'I-MM',7X,'J-NX',7X,'J-QZ',7X,'J=MM')
      DO 694 I=1,KI3
      J=MI3(I)
      WRITE(6,697)J,(FI3(I,K),K=1,6)
697   FORMAT(1X,I4,6F10.4)
694   CONTINUE
      DO 700 I=1,KI3
      KA1=MI3(I)
      ME(1)=EP(KA1,1)
  	ME(2)=EP(KA1,2)
 	T=EP(KA1,3)
      CALL MIC(H,LT,E,A,F,X,Z)
      A=LT(1,2)
      LT(1,2)=LT(2,1)
      LT(2,1)=A
      M1=ME(1)
      DO 698 K=1,3
       J=RR(M1,K)
	if(j.eq.0) goto 698
      DO 699 L=1,3
      FF(J)=FF(J)-LT(K,L)*FI3(I,L)
699   CONTINUE
698   CONTINUE
      M2=ME(2)
      DO 702 K=1,3
      J=RR(M2,K)
		if(j.eq.0) goto 702
      DO 704 L=1,3
      L1=3+L
      FF(J)=FF(J)-LT(K,L)*FI3(I,L1)
704   CONTINUE
702   CONTINUE
700   CONTINUE
710   WRITE(6,719)
719   FORMAT(//30X,'NODAL EQUIVALENT LOAD')
      WRITE(6,712)(FF(I),I=1,NX)
712   FORMAT(1X,6F11.3)
      RETURN
      END

      SUBROUTINE  FORDIS
      INTEGER  AA,RR,GG,T
      REAL KK
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A5/AA(200),KK(6000),NX,NH
      COMMON/A6/FF(60),UV(20,6)
      DIMENSION DD(3)
      DO 851 I=2,NX
      JJ=AA(I)
      JK=I-1
      JL=AA(JK)+I-JJ+1
      DO 851 J=JL,JK
      JP=AA(I)-I+J
      FF(I)=FF(I)-KK(JP)*FF(J)
851   CONTINUE
      DO 857 I=1,NX
      JJ=AA(I)
      FF(I)=FF(I)/KK(JJ)
857   CONTINUE
      DO 854 J4=2,NX
      I=2+NX-J4
      JJ=AA(I-1)+I-AA(I)+1
      M=AA(I)-I
      JP=I-1
      DO 856 J=JJ,JP
      JM=M+J
      FF(J)=FF(J)-KK(JM)*FF(I)
856   CONTINUE
854   CONTINUE
      WRITE(6,917)
917   FORMAT(//10X,'NODAL DISPLACEMENT'/2X,'ELEMENT',7X,
     #'XX',12X,'ZZ',11X,'ANGLE')
      DO 947 I=1,NP
      DO 945 J=1,3
      IR=RR(I,J)
      IF(IR.EQ.0) DD(J)=0
      IF(IR.NE.0) DD(J)=FF(IR)
945   CONTINUE
      WRITE(6,946)I,(DD(K),K=1,3)
946   FORMAT(2X,I4,3F14.6)
947   CONTINUE
      RETURN
      END

      SUBROUTINE  FORFM
      INTEGER  RR,AA,GG,EP,RP,T
      REAL LT,KE,KK
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A4/KI2,KI3,MI2(20),MI3(20),FI2(20,3),FI3(20,6)
      COMMON/A5/AA(200),KK(6000),NX,NH
      COMMONL/A6/FF(60),UV(20,6)
      DIMENSION LT(3,3),U(6),V(6),FMJ(3),FMK(3),KE(6,6)
      WRITE(6,988)
988   FORMAT(//30X,'MEMBER END FORCE'/14X,'AXIAL FORCE',6X,
     #'SHEARING FORCE',5X,'BENDING MOMEMT')
      DO 683 I=1,NE
      DO 683  L=1,6
683   UV(I,L)=0.0
      DO 845 I9=1,NE
      DO 341 J=1,6
      U(J)=0
341   V(J)=0
      DO 360 J=1,3
      FMJ(J)=0
360   FMK(J)=0
      ME(1)=EP(I9,1)
  	ME(2)=EP(I9,2)
 	T=EP(I9,3)
      CALL MIC(H,LT,E,A,F,X,Z)
      CALL KET(X,Z,H,E,A,F,KE)
      DO 323 L=1,3
      IH1=ME(1)
      IR=RR(IH1,L)
      IF(IR.EQ.0) U(L)=0
      IF(IR.NE.0) U(L)=FF(IR)
323   CONTINUE
      IH2=ME(2)
      DO 343 I=1,3
      IR=RR(IH2,I)
      MIY=3+I
      IF(IR.EQ.0) U(MIY)=0
      IF(IR.NE.0) U(MIY)=FF(IR)
343   CONTINUE
      DO 357 K=1,6
      DO 357 L=1,6
      V(K)=V(K)+KE(K,L)*U(L)
357   CONTINUE
680   DO 370 K=1,3
      DO 370 L=1,3
      UV(I9,K)=UV(I9,K)+LT(K,L)*V(L)
370   FMJ(K)=UV(I9,K)
      DO 730 K=1,3
      J=3+K
      DO 730 L=1,3
      L2=3+L
      UV(I9,J)=UV(I9,J)+LT(K,L)*V(L2)
730   FMK(K)=UV(I9,J)
      IF(KI3.EQ.0) GOTO 820
      DO 860 I8=1,KI3
      KA1=MI3(I8)
      IF(KA1.NE.I9) GOTO 860
      DO 862 I6=1,6
      UV(I9,I6)=UV(I9,I6)+FI3(I8,I6)
862   CONTINUE
      DO 864 I6=1,3
      FMJ(I6)=UV(I9,I6)
864   CONTINUE
      DO 863 K=1,3
      I6=3+K
      FMK(K)=UV(I9,I6)
863   CONTINUE
      GOTO 820
860   CONTINUE
820   WRITE(6,821) I9
821   FORMAT(4X,'NE=',I4)
      WRITE(6,823)(FMJ(K),K=1,3)
      WRITE(6,824)(FMK(K),K=1,3)
824   FORMAT(6X,'FMK=',3F16.4)
823   FORMAT(6X,'FMJ=',3F16.4)
845   CONTINUE
      RETURN
      END

      SUBROUTINE FORRC
      INTEGER  RR,AA,GG,EP,RP,T
      REAL LT,KK
      COMMON/A1/NE,NP,NR,MT,MX
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      COMMON/A4/KI2,KI3,MI2(20),MI3(20),FI2(20,3),FI3(20,6)
      COMMON/A5/AA(200),KK(6000),NX,NH
      COMMON/A6/FF(60),UV(20,6)
      DIMENSION LT(3,3),RC(20,3)
      WRITE(6,901)
901   FORMAT(//10X,'SUPPOT REACTION')
      WRITE(6,904)
904   FORMAT(2X,'NUMBER OF NODE',7X,'X-FORCE',9X,'Z-FORCE',7X,
     #'BENDING MOMEMT')
      DO 902 I=1,NR
      DO 902 J=1,3
902   RC(I,J)=0.0
      DO 995 I=1,NR
      IZ1=RP(I,1)
      DO 930 J=1,NE
      ME(1)=EP(J,1)
  	ME(2)=EP(J,2)
 	T=EP(J,3)
      J1=ME(1)
      J2=ME(2)
      CALL MIC (H,LT,E,A,F,X,Z)
      C1=LT(1,2)
      LT(1,2)=LT(2,1)
      LT(2,1)=C1
      IF(IZ1.NE.J1) GOTO 920
      DO 910 K1=1,3
      DO 910 K2=1,3
      RC(I,K1)=RC(I,K1)+LT(K1,K2)*UV(J,K2)
910   CONTINUE
      GOTO 930
920   IF(IZ1.NE.J2) GOTO 930
      DO 915 K1=1,3
      DO 915 K2=1,3
      L2=3+K2
      RC(I,K1)=RC(I,K1)+LT(K1,K2)*UV(J,L2)
915   CONTINUE
930   CONTINUE
      DO 940 J=1,KI2
      IA3=MI2(J)
      IF(IZ1.NE.IA3) GOTO 940
      DO 908 K1=1,3
      RC(I,K1)=RC(I,K1)-FI2(J,K1)
908   CONTINUE
      GOTO 950
940   CONTINUE
950   WRITE(6,909) IZ1,(RC(I,I1),I1=1,3)
909   FORMAT(5X,I4,5X,3F16.4)
995   CONTINUE
      WRITE(6,992)
992   FORMAT(15X,'*******************************************')
      WRITE(6,993)
993   FORMAT(35X,'END')
      RETURN
      END

      SUBROUTINE KET(X,Z,H,E,A,F,KE)
      DIMENSION KE(6,6)
      REAL KE
      DO 118 K=1,6
      DO 118 J=1,6
118   KE(K,J)=0.0
      H=1/H
      CX=X*H
      CZ=Z*H
      KE(1,1)=E*A*H*CX*CX+12.0*E*F*H*H*H*CZ*CZ
      KE(4,4)=KE(1,1)
      KE(1,2)=(E*A*H-12.0*E*F*H*H*H)*CX*CZ
      KE(4,5)=KE(1,2)
      KE(2,2)=E*A*H*CZ*CZ+12.0*E*F*H*H*H*CX*CX
      KE(5,5)=KE(2,2)
      KE(1,3)=(6.0*E*F*H*H)*CZ
      KE(1,6)=KE(1,3)
      KE(2,3)=-(6.0*E*F*H*H)*CX
      KE(2,6)=KE(2,3)
      KE(3,3)=4.0*E*F*H
      KE(6,6)=KE(3,3)
      KE(2,5)=-KE(2,2)
      KE(1,4)=-KE(1,1)
      KE(2,4)=-KE(1,2)
      KE(1,5)=-KE(1,2)
      KE(3,6)=KE(3,3)/2
      KE(3,4)=-KE(1,3)
      KE(4,6)=-KE(1,3)
      KE(5,6)=-KE(2,3)
      KE(3,5)=-KE(2,3)
      DO 127 I=1,6
      DO 127 J=1,I
127   KE(I,J)=KE(J,I)
      H=1/H
      RETURN
      END

      SUBROUTINE MIC(RL,LT,E,A,RI,X,Z)
      INTEGER T
      REAL LT
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      COMMON/A3/EP(20,3),RP(20,4),XZ(30,2),AE(5,3)
      DIMENSION LT(3,3)
      IH1=ME(1)
      IH2=ME(2)
      X=XZ(IH2,1)-XZ(IH1,1)
      Z=XZ(IH2,2)-XZ(IH1,2)
      RL=SQRT(X*X+Z*Z)
      DO 345 K=1,3
      DO 345 J=1,3
345   LT(K,J)=0.0
      LT(3,3)=1
      LT(1,1)=X/RL
      LT(2,2)=LT(1,1)
      LT(1,2)=Z/RL
      LT(2,1)=-LT(1,2)
      E=AE(T,1)
      A=AE(T,2)
      RI=AE(T,3)
      RETURN
      END

      SUBROUTINE GGT
	 INTEGER T
      COMMON/A2/ME(2),RR(20,3),GG(6),T
      INTEGER RR,GG
      IH1=ME(1)
      IH2=ME(2)
      DO 103 I=1,3
      GG(I)=RR(IH1,I)
103   CONTINUE
      DO 106 I=1,3
      IMM=I+3
       GG(IMM)=RR(IH2,I)
106   CONTINUE
      RETURN
      END