Fortran Coder

查看: 11214|回复: 5
打印 上一主题 下一主题

[数值问题] 程序的变量值变化的很奇怪,求助

[复制链接]

10

帖子

4

主题

0

精华

入门

F 币
51 元
贡献
30 点
跳转到指定楼层
楼主
发表于 2018-5-17 15:14:04 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
程序执行过程中并没有给I,J赋值,为什么提示超出范围?
[Fortran] 纯文本查看 复制代码
001PROGRAM PFSAP!ANALYSIS PROGRAM FOR PLANE FRAME
002INTEGER NO,NJ,N,NE,NM,NPJ,NPF,I,J,IE,IP,KK,I1
003REAL K(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100),EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6),FE(6),ADE(6),DE(6),RT(6,6),AFE(6),F(3),D(100)
004INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),M(6),JEAI(100)
005CHARACTER *80 TITLE
0061 OPEN (6,FILE='PFSAP.IN')
007OPEN (8,FILE='PFSAP.OUT')
008READ(6,'(I3)')NO
009IF(NO.EQ.0) STOP
010WRITE(8,"('(NO.='I3,')'/)")NO
011CALL PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)
012DO I=1,N
013   P(I)=0.0
014   DO J=1,N
015      K(I,J)=0
016   END DO
017END DO
018DO IE=1,NE
019   CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
020   CALL MR(R,IE,JE,X,Y)
021   CALL MAKE(KE,R,AKE)
022   CALL CALM(M,IE,JN,JE)
023   CALL MK(K,AKE,M)
024END DO
025DO IP=1,NPF
026   CALL MR(R,JPF(1,IP),JE,X,Y)
027   CALL TRAN(R,RT)
028   CALL PE(FE,IP,JPF,PF,AL)
029   CALL MULV6(RT,FE,AFE)
030   CALL CALM(M,JPF(1,IP),JN,JE)
031   CALL MF(P,AFE,M)
032END DO
033DO I=1,NPJ
034   P(JPJ(I))=P(JPJ(I))+PJ(I)
035END DO
036CALL SLOV(K,P,D,N)
037WRITE(8,*) "/10X,5('* *'),'RESULTS OF CALCULATION',5('* *')"
038WRITE(8,40)
03940 FORMAT(/5X,'NO.N',4X,'X-DISPLACEMENT',2X,'Y-DISPLEMENT',3X,'ANG.ROT.(RAD)')
040DO KK=1,NJ
041   DO I=1,3
042      F(I)=0.0
043      I1=JN(I,KK)
044      IF(I1.GT.0) F(I)=D(I1)
045   END DO
046WRITE(8,70)KK,F(1),F(2),F(3)
04770 FORMAT(I8,2X,3G16.5)
048END DO
049WRITE(8,80)
05080 FORMAT(/'NO.E',5X,'N(1)',8X,'Q(1)',8X,'M(1)',8X,'N(2)',8X,'Q(2)',8X,'M(2)')
051DO IE=1,NE
052   CALL MADE(IE,JN,JE,D,ADE)
053   CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
054   CALL MR(R,IE,JE,X,Y)
055   CALL MULV6(R,ADE,DE)
056   CALL MULV6(KE,DE,FF)
057     DO IP=1,NPF
058        IF(JPF(1,IP).EQ.IE)THEN
059          CALL PE(FE,IP,JPF,PF,AL)
060          DO I=1,6
061             FF(I)=FF(I)-FE(I)
062          END DO
063        END IF
064     END DO
065    WRITE(8,110)IE,(FF(I),I=1,6)
066110    FORMAT(I5,2X,6G12.5)
067END DO
068GOTO 1
069END PROGRAM PFSAP
070SUBROUTINE PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)
071REAL X(100),Y(100),EAI(3,100),PJ(100),PF(2,100)
072INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100)
073CHARACTER*80 TITLE
074READ(6,*)TITLE
075WRITE(8,*)TITLE
076READ(6,*)NJ,N,NE,NM,NPJ,NPF
077WRITE(8,3)NJ,N,NE,NM,NPJ,NPF
0783 FORMAT('NJ='I2,5X,'N='I2,5X,'NE='I2,5X,'NM='I2,5X,'NPJ='I2,5X,'NPF='I2)
079WRITE(8,5)
0805 FORMAT(/3X,'NO. (1)  (2)  (3)',10X,'X',8X,'Y')
081READ(6,10)((JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)
08210 FORMAT(2(3I5,2G16.4))
083WRITE(8,11)(I,(JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)
08411 FORMAT(2X,I2,3I5,4X,2F10.3)
085WRITE(8,30)
08630 FORMAT(/10X,'ELEMENT NO. NODE-1  NODE-2  MATERIALS')
087READ(6,40)(JE(1,I),JE(2,I),JEAI(I),I=1,NE)
08840 FORMAT(5(3I5))
089DO I=1,NE
090   WRITE(8,'(14X,I2,3(7X,I3))')I,JE(1,I),JE(2,I),JEAI(I)
091END DO
092READ(6,*)((EAI(I,J),I=1,3),J=1,NM)
093WRITE(8,60)(J,(EAI(I,J),I=1,3),J=1,NM)
09460 FORMAT(/3X,'NO.MAT',3X,'ELASTIK MOUDULUS',5X,'AREA',5X,'MOMENT OF INERTIA'/(I6,3G16.6))
095IF(NPJ.NE.0)THEN
096  WRITE(8,15)
09715 FORMAT(/20X,'NODEL LOADS')
098  WRITE(8,12)
09912FORMAT(/16X,'NO.DISP.     VALUE')
100  READ(6,70)(JPJ(I),PJ(I),I=1,NPJ)
10170 FORMAT(5(I5,G16.4))
102  DO I=1,NPJ
103   WRITE(8,'(14X,I7,F16.3)')JPJ(I),PJ(I)
104  END DO
105ELSE
106END IF
107IF(NPF.NE.0) THEN
108  WRITE(8,13)
10913 FORMAT(/20X,'NON-NODEL LOADS')
110  WRITE(8,14)
11114 FORMAT(11X,'NO.E NO.LOAD.MODEL',8X,'A',9X,'C')
112  READ(6,100)(JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)
113100 FORMAT(2(2I5,2G16.4))
114  DO I=1,NPF
115     WRITE(8,120)(JPF(J,I),J=1,2),PF(1,I),PF(2,I)
116120 FORMAT(6X,2I8,10X,2F10.3)
117  END DO
118ELSE
119END IF
120RETURN
121END
122SUBROUTINE MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
123INTEGER I,J,II,JJ,MT
124REAL L,A1,A2,A3,A4
125REAL KE(6,6),X(100),Y(100),EAI(3,100),AL(100)
126INTEGER JE(2,100),JEAI(100)
127II=JE(1,IE)
128JJ=JE(2,IE)
129MT=JEAI(IE)
130L=SQRT((X(JJ)-X(II))**2+(Y(JJ)-Y(II))**2)
131AL(IE)=L
132A1=EAI(1,MT)*EAI(2,MT)/L
133A2=EAI(1,MT)*EAI(3,MT)/L**3
134A3=EAI(I,MT)*EAI(3,MT)/L**2
135A4=EAI(1,MT)*EAI(3,MT)/L
136KE(1,1)=A1
137KE(1,4)=-A1
138KE(2,2)=12*A2
139KE(2,3)=6*A3
140KE(2,5)=-12*A2
141KE(2,6)=6*A3
142KE(3,3)=4*A4
143KE(3,5)=-6*A3
144KE(3,6)=2*A4
145KE(4,4)=A1
146KE(5,5)=12*A2
147KE(5,6)=-6*A3
148KE(6,6)=4*A4
149DO I=1,6
150   DO J=I,6
151      KE(J,I)=KE(I,J)
152   END DO
153END DO
154RETURN
155END
156SUBROUTINE MR(R,IE,JE,X,Y)
157INTEGER I,J
158INTEGER JE(2,100)
159REAL L,CX,CY
160REAL R(6,6),X(100),Y(100)
161I=JE(1,IE)
162J=JE(2,IE)
163L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
164CX=(X(J)-X(I))/L
165CY=(Y(J)-Y(I))/L
166DO J=1,6
167   DO I=1,6
168      R(I,J)=0.0
169   END DO
170END DO
171DO I=1,4,3
172   R(I,I)=CX
173   R(I,I+1)=CY
174   R(I+1,I)=-CY
175   R(I+1,I+1)=CX
176   R(I+2,I+2)=1.0
177END DO
178RETURN
179END
180SUBROUTINE MAKE(KE,R,AKE)
181REAL KE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)
182CALL TRAN(R,RT)
183CALL MULV(RT,KE,TMP)
184CALL MULV(TMP,R,AKE)
185RETURN
186END
187SUBROUTINE CALM(M,IE,JN,JE)
188INTEGER I,IE
189INTEGER M(6),JN(3,100),JE(2,100)
190DO I=1,3
191   M(I)=JN(I,JE(1,IE))
192   M(I+3)=JN(I,JE(2,IE))
193END DO
194RETURN
195END
196SUBROUTINE MK(K,AKE,M)
197INTEGER I,J
198REAL K(200,200),AKE(6,6)
199INTEGER M(6)
200DO I=1,6
201   DO J=1,6
202      IF(M(I).NE.0.AND.M(J).NE.0) K(M(I),M(J))=K(M(I),M(J))+AKE(I,J)
203   END DO
204END DO
205RETURN
206END
207SUBROUTINE PE(FE,IP,JPF,PF,AL)
208INTEGER I,IND
209INTEGER JPF(2,100)
210REAL L,A,C
211REAL FE(6),PF(2,100),AL(100)
212A=PF(1,IP)
213C=PF(2,IP)
214L=AL(JPF(1,IP))
215IND=JPF(2,IP)
216DO I=1,6
217   FE(I)=0.0
218END DO
219SELECT CASE(IND)
220  CASE(1)
221   FE(2)=(7.0*A/20.0+3.0*C/20.0)*L
222   FE(3)=(A/20.0+C/30.0)*L**2
223   FE(5)=(3.0*A/20.0+7.0*C/20.0)*L
224   FE(6)=-(A/30.0+C/20.0)*L**2
225   RETURN
226  CASE(2)
227   FE(5)=A*C**3*(2*L-C)/(2*L**3)
228   FE(6)=-A*C**3*(4.0*L-3.0*C)/(12.0*L**2)
229   FE(2)=A*C-FE(5)
230   FE(3)=A*C**2*(6.0*L**2-8.0*C*L+3.0*C**2)/(12.0*L**2)
231   RETURN
232  CASE(3)
233   FE(2)=A*(L-C)**2*(L+2*C)/L**3
234   FE(3)=A*C*(L-C)**2/L**2
235   FE(5)=A-FE(2)
236   FE(6)=-A*C**2*(L-C)/L**2
237   RETURN
238  CASE(4)
239   FE(2)=-6.0*A*C*(L-C)/L**3
240   FE(3)=A*(L-C)*(L-3.0*C)/L**2
241   FE(5)=-FE(2)
242   FE(6)=A*C*(3.0*C-2.0*L)/L**2
243   RETURN
244  CASE(5)
245   FE(1)=A*(1-C/L)
246   FE(4)=A*C/L
247  CASE(6)
248   FE(1)=C*L/2
249   FE(4)=FE(1)
250 END SELECT
251END
252SUBROUTINE MULV6(A,B,C)
253INTEGER I,J
254REAL C(6),A(6,6),B(6)
255DO I=1,6
256   C(I)=0.0
257   DO J=1,6
258      C(I)=C(I)+A(I,J)*B(J)
259   END DO
260END DO
261RETURN
262END
263SUBROUTINE MF(P,AFE,M)
264INTEGER I
265INTEGER M(6)
266REAL P(100),AFE(6)
267DO I=1,6
268   IF(M(I).NE.0)P(M(I))=AFE(I)+P(M(I))
269END DO
270RETURN
271END
272SUBROUTINE SLOV(K,P,D,N)
273INTEGER I,J,KK,N
274REAL K(200,200),P(100),D(100)
275DO I=1,100
276   D(I)=P(I)
277END DO
278DO KK=1,N-1
279   DO I=KK+1,N
280      C=-K(KK,I)/K(KK,KK)
281      DO J=I,N
282         K(I,J)=K(I,J)+C*K(KK,J)
283      END DO
284      D(I)=D(I)+C*D(KK)
285   END DO
286END DO
287D(N)=D(N)/K(N,N)
288DO I=N-1,1,-1
289   DO J=I+1,N
290      D(I)=D(I)-K(I,J)*D(J)
291   END DO
292   D(I)=D(I)/K(I,I)
293END DO
294RETURN
295END
296SUBROUTINE MADE(IE,JN,JE,D,ADE)
297INTEGER I
298REAL ADE(6),D(100)
299INTEGER IE,JN(3,100),JE(2,100)
300DO I=1,6
301   ADE(I)=0
302END DO
303DO I=1,3
304   IF(JN(I,JE(1,IE)).NE.0) ADE(I)=D(JN(I,JE(1,IE)))
305   IF(JN(I,JE(2,IE)).NE.0) ADE(I+3)=D(JN(I,JE(2,IE)))
306END DO
307RETURN
308END
309SUBROUTINE TRAN(R,RT)
310INTEGER I,J
311REAL R(6,6),RT(6,6)
312DO I=1,6
313   DO J=1,6
314      RT(I,J)=R(J,I)
315   END DO
316END DO
317RETURN
318END
319SUBROUTINE MULV(A,B,C)
320INTEGER I,J,K
321REAL A(6,6),B(6,6),C(6,6)
322DO I=1,6
323   DO J=1,6
324      C(I,J)=0.0
325      DO K=1,6
326         C(I,J)=C(I,J)+A(I,K)*B(K,J)
327      END DO
328   END DO
329END DO
330RETURN
331END




分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

27

帖子

0

主题

0

精华

熟手

F 币
212 元
贡献
101 点
沙发
发表于 2018-5-17 21:54:39 | 只看该作者
请提供文件"PFSAP.IN".在没有此文件的情况下,直觉建议您:
1.主程序和子例行程序加上 implicit none;
2.详细检查所有变量、数组的声明和初始化

10

帖子

4

主题

0

精华

入门

F 币
51 元
贡献
30 点
板凳
 楼主| 发表于 2018-5-17 22:16:31 | 只看该作者
吉大渣渣-固体 发表于 2018-5-17 21:54
请提供文件"PFSAP.IN".在没有此文件的情况下,直觉建议您:
1.主程序和子例行程序加上 implicit none;
2.详 ...

2
EXAMPLE----(5-2)
5,7,7,2,3,0
0,1,0,0.0,0.0,2,3,0,150.0,35.0
4,5,0,300.0,70.0,6,7,0,150.0,160.0
0,0,0,300.0,320.0
1,2,1,2,3,1,1,4,1,2,4,2,3,4,2
4,5,2,3,5,2
3.0E6,144.0,0.0
3.0E6,180.0,0.0
1,1610.0,3,3220.0,5,1610.0
0

27

帖子

0

主题

0

精华

熟手

F 币
212 元
贡献
101 点
地板
发表于 2018-5-17 23:04:55 | 只看该作者
chaunceyyou 发表于 2018-5-17 22:16
2
EXAMPLE----(5-2)
5,7,7,2,3,0

坏在您的主程序中CALL PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)这个过程中子例行程序PUTIN中I,I1,IE,IP,J,KK这些变量都没有初始化,以Intel visual Fortran编译器为例,这些变量不初始化直接拿来算编译器默认给的值是-858993460.有这些不确定性的隐患存在,导致后面EAI的值不稳定是正常的。

衷心给您的敬告:
1.尽管您的程序以Release模式可以得过且过运行出结果并正常退出,但是不初始化就拿来算这个隐患希望您心里有数。
2.不写implicit none,不重视初始化的程序都是在耍流氓。

评分

参与人数 1F 币 +5 贡献 +5 收起 理由
fcode + 5 + 5 赞一个!

查看全部评分

10

帖子

4

主题

0

精华

入门

F 币
51 元
贡献
30 点
5#
 楼主| 发表于 2018-5-18 09:50:13 | 只看该作者
吉大渣渣-固体 发表于 2018-5-17 23:04
坏在您的主程序中CALL PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)这个过程中子例行程 ...

非常谢谢你,我刚接触fortran没多久,书上的程序用的版本比较老,再继续时一定注重规范

27

帖子

0

主题

0

精华

熟手

F 币
212 元
贡献
101 点
6#
发表于 2018-5-19 19:21:52 | 只看该作者
chaunceyyou 发表于 2018-5-18 09:50
非常谢谢你,我刚接触fortran没多久,书上的程序用的版本比较老,再继续时一定注重规范 ...

客气 祝学习愉快
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2025-5-4 12:59

Powered by Discuz! X3.4

© 2013-2025 Comsenz Inc.

快速回复 返回顶部 返回列表