求助……，关于一个平面杆系结构分析程序
PROGRAM PBSAP! C PLANE BAR STRUCTURE ANALYSIS PROGRAM
REAL KS
DIMENSION ID(100,3),EF(100,7),PJ(50,3),PE(50,4)
DIMENSION KD(100),KS(5000),PS(100),PF(100,6)
CHARACTER NAME*12
WRITE(*,10)
10 FORMAT(1X,22HINPUT DATA FILE NAME=>)
READ(*,*) NAME
OPEN(1,FILE=NAME,STATUS='OLD')
OPEN(2,FILE='PBO.DAT')
WRITE(*,1000)
WRITE(2,1000)
CALL INDATA(NJ,NE,NP,NF,ID,EF,PJ,PE)
CALL PTDATA(NJ,NE,NP,NF,ID,EF,PJ,PE)
CALL SETJID(NJ,NN,ID)
CALL SETKKD(NE,NN,LL,ID,EF,KD)
CALL TOTLKS(NE,ID,EF,KD,KS)
CALL TOTLPS(NP,NF,ID,EF,PJ,PE,PS)
CALL GAUSS3(NN,KD,KS,PS)
CALL INTEPF(NE,NF,ID,EF,PE,PS,PF)
CALL PTDPPF(NJ,NE,NN,LL,ID,PS,PF)
WRITE(*,2000)
WRITE(2,2000)
STOP
1000 FORMAT(//5X,30(1H.)/6X,1H.,57X,1H./6X,1H.,24X,9HP B S A P,24X,1H./6X,1H.,57X,1H./6X,1H.,10X,35HP L A N EB A RS T R U C T U R E,10X,1H./6X,1H.,57X,1H./6X,1H.,12X,30HA N A L Y S I SP R O G R A M,14X,1H./6X,1H.,57X,1H./5X,30(1H.)///)
2000 FORMAT(///22X,25H* * * * * E N D * * * * */)
END
SUBROUTINE INDATA(NJ,NE,NP,NF,ID,EF,PJ,PE)
DIMENSION ID(100,3),EF(100,7),PJ(50,3),PE(50,4)
READ(1,*) NJ,NE,NP,NF
DO 10N=1,NJ
READ(1,*) (ID(N,I),I=1,3)
10 CONTINUE
DO 20N=1,NE
READ(1,*) (EF(N,I),I=1,7)
20 CONTINUE
IF(NP.EQ.0)GOTO 40
DO 30N=1,NP
READ(1,*) (PJ(N,I),I=1,3)
30 CONTINUE
40 IF(NF.EQ.0)GOTO 60
DO 50N=1,NF
READ(1,*) (PE(N,I),I=1,4)
50 CONTINUE
60 RETURN
END
SUBROUTINE PTDATA(NJ,NE,NP,NF,ID,EF,PJ,PE)
DIMENSION ID(100,3),EF(100,7),PJ(50,3),PE(50,4)
WRITE(*,1000) NJ,NE,NP,NF
WRITE(2,1000) NJ,NE,NP,NF
WRITE(*,1100)
WRITE(2,1100)
DO 10N=1,NJ
WRITE(*,1200)N,(ID(N,I),I=1,3)
WRITE(2,1200)N,(ID(N,I),I=1,3)
10 CONTINUE
WRITE(*,1300)
WRITE(2,1300)
DO 20 N=1,NE
WRITE(*,1400) N,(INT(EF(N,I)),I=1,2),(EF(N,I),I=3,7)
WRITE(2,1400) N,(INT(EF(N,I)),I=1,2),(EF(N,I),I=3,7)
20 CONTINUE
IF(NP.EQ.0)GOTO 40
WRITE(*,1500)
WRITE(2,1500)
DO 30N=1,NP
WRITE(*,1600)(INT(PJ(N,I)),I=1,2),PJ(N,3)
WRITE(2,1600)(INT(PJ(N,I)),I=1,2),PJ(N,3)
30 CONTINUE
40 IF(NF.EQ.0) GOTO 60
WRITE(*,1700)
WRITE(2,1700)
DO 50N=1,NF
WRITE(*,1800)(INT(PE(N,I)),I=1,2),(PE(N,I),I=3,4)
WRITE(2,1800)(INT(PE(N,I)),I=1,2),(PE(N,I),I=3,4)
50 CONTINUE
60 RETURN
1000 FORMAT(/15X,40H* * * * * I N P U T D A T A * * * * *//5X,11HNODE NUMBER,16X,3HNE=,I3/5X,14HELEMENT NUMBER,13X,3HNE=,I3/5X,16HNODE LOAD NUMBER,11X,3HNP=,I3/5X,19HELEMENT LOAD NUMBER,8X,3HNF=,I3/)
1100 FORMAT(/5X,5(1H-),16HNODE INFORMATION,5(1H-)//3X,7HNODE NO,4X,5HX-DIR,3X,5HY-DIR,3X,5HQ-DIR/)
1200 FORMAT(5X,I2,8X,3(I2,6X))
1300 FORMAT(//5X,5(1H-),19HELEMENT INFORMATION,5(1H-)//1X,2HNO,3X,1HI,3X,1HJ,4X,6HLENGTH,7X,4HAREA,7X,6HMOMENT,6X,7HMODULUS,6X,5HANGLE/)
1400 FORMAT(1X,3(I2,2X),5(E10.4,2X))
1500 FORMAT(//5X,5(1H-),21HNODE LOAD INFORMATION,5(1H-)//3X,7HNODE NO,2X,14HLOAD - DIRECTION,4X,10HLOAD VALUS/)
1600 FORMAT(5X,I2,11X,I2,11X,F8.3)
1700 FORMAT(//5X,5(1H-),24HELEMENT LOAD INFORMATION,5(1H-)//1X,10HELEMENT NO,3X,9HLOAD TYPE,5X,10HLOAD VALUS,7X,13HLOAD DISTANCE/)
1800 FORMAT(5X,I2,10X,I2,10X,2(F8.3,11X))
END
SUBROUTINE SETJID(NJ,NN,ID)
DIMENSION ID(100,3)
NN=0
DO 50N=1,NJ
DO 40I=1,3
IF(ID(N,I)) 10,20,30
10 ID(N,I)=0
GOTO 40
20 NN=NN+1
ID(N,I)=NN
GOTO 40
30 NI=ID(N,I)
ID(N,I)=ID(NI,I)
40 CONTINUE
50 CONTINUE
RETURN
END
SUBROUTINE SETKKD(NE,NN,LL,ID,EF,KD)
DIMENSION ID(100,3),EF(100,7),KD(100),IJ(6)
KD(1)=1
DO 50N=2,NN
MA=0
DO 40M=1,NE
IE=EF(M,1)
JE=EF(M,2)
DO 10I=1,3
IJ(I)=ID(IE,I)
IJ(I+3)=ID(JE,I)
10 CONTINUE
NA=0
DO 20I=1,6
IF(N.EQ.IJ(I)) NA=1
20 CONTINUE
IF(NA.EQ.0) GOTO 40
DO 30I=1,6
IF(IJ(I).EQ.0) GOTO 30
NA=N-IJ(I)
IF(NA.GT.MA) MA=NA
30 CONTINUE
40 CONTINUE
KD(N)=KD(N-1)+MA+1
50 CONTINUE
LL=KD(NN)
RETURN
END
SUBROUTINE TOTLKS(NE,ID,EF,KD,KS)
REAL KS,KE
DIMENSION ID(100,3),EF(100,7)
DIMENSION KD(100),KS(5000),KE(6,6),IJ(6)
DO 40N=1,NE
CALL TOTLKE(N,EF,KE)
IE=EF(N,1)
JE=EF(N,2)
DO 10I=1,3
IJ(I)=ID(IE,I)
IJ(I+3)=ID(JE,I)
10 CONTINUE
DO 30I=1,6
IF(IJ(I).EQ.0) GOTO 30
DO 20J=1,6
IF(IJ(J).LT.IJ(I))GOTO 20
MW=KD(IJ(J))-(IJ(J)-IJ(I))
KS(MW)=KS(MW)+KE(I,J)
20 CONTINUE
30 CONTINUE
40 CONTINUE
RETURN
END
SUBROUTINE TOTLKE(N,EF,KE)
REAL KE
DIMENSION EF(100,7),KE(6,6),AE(4),BE(6)
AE(1)=EF(N,4)*EF(N,6)/EF(N,3)
AE(2)=12*EF(N,5)*EF(N,6)/EF(N,3)/EF(N,3)/EF(N,3)
AE(3)=6*EF(N,5)*EF(N,6)/EF(N,3)/EF(N,3)
AE(4)=4*EF(N,5)*EF(N,6)/EF(N,3)
EH=EF(N,7)*3.14159/180.
CX=COS(EH)
CY=SIN(EH)
BE(1)=AE(1)*CX*CX+AE(2)*CY*CY
BE(2)=(AE(1)-AE(2))*CX*CY
BE(3)=AE(1)*CY*CY+AE(2)*CX*CX
BE(4)=-AE(3)*CY
BE(5)=AE(3)*CX
BE(6)=AE(4)
KE(1,1)=BE(1)
KE(1,2)=BE(2)
KE(1,3)=BE(4)
KE(1,4)=-BE(1)
KE(1,5)=-BE(2)
KE(1,6)=BE(4)
KE(2,2)=BE(3)
KE(2,3)=BE(5)
KE(2,4)=-BE(2)
KE(2,5)=-BE(3)
KE(2,6)=BE(5)
KE(3,3)=BE(6)
KE(3,4)=-BE(4)
KE(3,5)=-BE(5)
KE(3,6)=BE(6)/2
KE(4,4)=BE(1)
KE(4,5)=BE(2)
KE(4,6)=-BE(4)
KE(5,5)=BE(3)
KE(5,6)=-BE(5)
KE(6,6)=BE(6)
DO 10I=1,6
DO 10J=I,6
KE(J,I)=KE(I,J)
10 CONTINUE
RETURN
END
SUBROUTINE TOTLPS(NP,NF,ID,EF,PJ,PE,PS)
DIMENSION ID(100,3),EF(100,7),PJ(50,3),PE(50,4)
DIMENSION PS(100),IJ(6),P(6)
IF(NP.EQ.0)GOTO 20
DO 10N=1,NP
JJ=PJ(N,1)
JD=PJ(N,2)
JQ=ID(JJ,JD)
PS(JQ)=PS(JQ)+PJ(N,3)
10 CONTINUE
20 IF(NF.EQ.0)GOTO 50
DO 40N=1,NF
ME=PE(N,1)
JT=PE(N,2)
EP=PE(N,3)
EA=PE(N,4)
EL=EF(ME,3)
EH=EF(ME,7)
CALL SETEEP(JT,EP,EA,EL,P)
CALL TRANST(EH,P)
IE=EF(ME,1)
JE=EF(ME,2)
DO 30I=1,3
IJ(I)=ID(IE,I)
IJ(I+3)=ID(JE,I)
30 CONTINUE
DO 40I=1,6
IF(IJ(I).EQ.0)GOTO 40
PS(IJ(I))=PS(IJ(I))-P(I)
40 CONTINUE
50 RETURN
END
SUBROUTINE SETEEP(JT,EP,EA,EL,P)
DIMENSION P(6)
C1=EA/EL
C2=1-C1
GOTO (1,2,3),JT
1 P(1)=0
P(2)=EP*(1.+2.*C1)*C2*C2
P(3)=EP*C1*C2*C2*EL
P(4)=0
P(5)=EP*(1.+2.*C2)*C1*C1
P(6)=-EP*C2*C1*C1*EL
GOTO 10
2 P(1)=0
P(2)=EP*EL/2
P(3)=EP*EL*EL/12.
P(4)=0
P(5)=EP*EL/2
P(6)=-EP*EL*EL/12.
GOTO 10
3 P(1)=0
P(2)=-6*EP*C1*C2/EL
P(3)=-EP*C2*(3.*C1-1.)
P(4)=0
P(5)=6*EP*C1*C2/EL
P(6)=-EP*C1*(3.*C2-1.)
10 RETURN
END
SUBROUTINE TRANST(EH,P)
DIMENSION P(6),F(6)
EH=EH*3.14159/180
CX=COS(EH)
CY=SIN(EH)
DO 10N=1,6
F(N)=P(N)
10 CONTINUE
P(1)=F(1)*CX-F(2)*CY
P(2)=F(1)*CY+F(2)*CX
P(4)=F(4)*CX-F(5)*CY
P(5)=F(4)*CY+F(5)*CX
RETURN
END
SUBROUTINE GAUSS3(NN,KD,KS,PS)
REAL KS
DIMENSION KD(100),KS(5000),PS(100)
DO 20N=1,NN-1
DO 20I=N+1,NN
IF((KD(I)-(I-N)).LE.DK(I-1)) GOTO 20
C=KS(KD(I)-(I-N))/KS(KD(N))
DO 10J=I,NN
IF((KD(J)-(J-N)).LE.KD(J-1)) GOTO 10
KS(KD(J)-(J-I))=KS(KD(J)-(J-I))-KS(KD(J)-(J-N))*C
10 CONTINUE
PS(I)=PS(I)-C*PS(N)
20 CONTINUE
PS(NN)=PS(NN)/KS(KD(NN))
DO 40I=NN-1,1,-1
DO 30J=I+1,NN
IF((KD(J)-(J-I)).LE.KD(J-I))GOTO 30
PS(I)=PS(I)-KS(KD(J)-(J-I))*PS(J)
30 CONTINUE
P(I)=PS(I)/KS(KD(I))
40 CONTINUE
RETURN
END
SUBROUTINE INTEPF(NE,NF,ID,EF,PE,PS,PF)
DIMENSION ID(100,3),EF(100,7),PE(50,4),PS(100)
DIMENSION PF(100,6),AD(4),IJ(6),P(6)
DO 50N=1,NE
IE=EF(N,1)
JE=EF(N,2)
DO 10I=1,3
IJ(I)=ID(IE,I)
IJ(I+3)=ID(JE,I)
10 CONTINUE
DO 20I=1,6
P(I)=0
IF(IJ(I).EQ.0)GOTO 20
P(I)=PS(IJ(I))
20 CONTINUE
EH=-EF(N,7)
CALL TRANST(EH,P)
AE(1)=EF(N,4)*EF(N,6)/EF(N,3)
AE(2)=12*EF(N,5)*EF(N,6)/EF(N,3)/EF(N,3)/EF(N,3)
AE(3)=6*EF(N,5)*EF(N,6)/EF(N,3)/EF(N,3)
AE(4)=4*EF(N,5)*EF(N,6)/EF(N,3)
PF(N,1)=AE(1)*(P(1)-P(4))
PF(N,2)=AE(2)*(P(2)-P(5))+AE(3)*(P(3)+P(6))
PF(N,3)=AE(3)*(P(2)-P(5))+AE(4)*(P(3)+P(6)/2)
PF(N,4)=-AE(1)*(P(1)-P(4))
PF(N,5)=-AE(2)*(P(2)-P(5))-AE(3)*(P(3)+P(6))
PF(N,6)=AE(3)*(P(2)-P(5))+AE(4)*(P(3)/2+P(6))
IF(NF.EQ.0) GOTO 45
DO 40I=1,NF
ME=PE(I,1)
IF(ME.NE.N) GOTO 40
JT=PE(I,2)
EP=PE(I,3)
EA=PE(I,4)
EL=EF(N,3)
CALL SETEEP(JT,EP,EA,EL,P)
DO 30J=1,6
PF(N,J)=PF(N,J)+P(J)
30 CONTINUE
40 CONTINUE
45 PF(N,1)=-PF(N,1)
PF(N,3)=-PF(N,3)
PF(N,5)=-PF(N,5)
PF(N,6)=-PF(N,6)
50 CONTINUE
RETURN
END
SUBROUTINE PTDPPF(NJ,NE,NN,LL,ID,PS,PF)
DIMENSION ID(100,3),PS(100),PF(100,6),P(3)
WRITE(*,1000) NN,LL
WRITE(2,1000) NN,LL
DO 10N=1,NJ
WRITE(*,1100) N,(ID(N,I),I=1,3)
WRITE(2,1100) N,(ID(N,I),I=1,3)
10 CONTINUE
WRITE(*,1200)
WRITE(2,1200)
DO 30N=1,NJ
DO 20I=1,3
P(I)=0
IF(ID(N,I).EQ.0) GOTO 20
P(I)=PS(IID(N,I))
20 CONTINUE
WRITE(*,1300)N,(P(I),I=1,3)
WRITE(2,1300)N,(P(I),I=1,3)
30 CONTINUE
WRITE(*,1400)
WRITE(2,1400)
DO 40 N=1,NE
WRITE(*,1500)N,(PF(N,I),I=1,6)
WRITE(2,1500)N,(PF(N,I),I=1,6)
40 CONTINUE
RETURN
1000 FORMAT(//14X,41H* * * * * O U T P U TD A T A * * * * * //5X,17HDEGREE OF FREEDOM,10X,3HNN=,I3/5X,22HTOTAL STIFFNESS NUMBER,5X,3HLL=,I3//5X,5(1H-),29HNODE LOCATION VECATION VECTOR,5(1H-)//3X,7HNODE NO,4X,5HX-DIR,3X,5HY-DIR,3X,5HQ-DIR/)
1100 FORMAT(5X,I2,8X,3(I2,6X))
1200 FORMAT(//5X,5(1H-),17HNODE DISPLACEMENT,5(1H-)//3X,7HNODE NO,8X,6HX-DISP,12X,6HY-DISP,11X,8HROTATION/)
1300 FORMAT(5X,I2,7X,3(E14.8,4X))
1400 FORMAT(//5X,5(1H-),22HELEMENT INTERNAL FORCE,5(1H-)//1X,2HNO,5X,3HI-N,8X,3HI-Q,8X,3HI-M,8X,3HJ-N,8X,3HJ-Q,8X,3HJ-M/)
1500 FORMAT(1X,I2,1X,6(E10.4,1X))
END 刚才是全部的源程序，好像在SUBROUTINE INTEP里面有错误.
帮帮忙 :'( :'( :'(
清帮帮忙，拜托拉。 :lol
各位，我搞定拉，不用劳烦大家拉。
请帮我一下，你是怎么搞定的？求助！
朋友，你是怎么搞定的，我正好也在学习这一个程序，要不请将你的正确代码发我一份吧，coldmoonflower@163.com，先谢谢啊！！！
页:
[1]