版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

返回列表

黑曼巴1991

铁杆木虫 (正式写手)

应助: 4 (幼儿园)
金币: 7391.8
散金: 335
红花: 2
帖子: 930
在线: 292小时
虫号: 2552663
注册: 2013-07-18
专业: 传动机械学

[求助] fortran编程过程出错了麻烦给解决一下谢谢！！

程序运行过程出错了不知道怎么办麻烦高手给看看。。我用的是FTN95编译器打开.f90的文件的~~ （复件也有程序）
PROGRAM POINTEHL
DIMENSION THETA(15),EALFA(15),EBETA(15)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
COMMON /COMEK/EK,EAL,EBE
DATA N,PAI,Z,EAL,EBE,E1,EDA0,RX,RY,X0,XE,W0,US/65,3.14159265,0.68,1.0,1.0,2.21E11,0.0283,0.05,0.05,-2.5,1.5,39.24,1.5/
DATA THETA/10.,20.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90./
DATA EALFA/6.612,3.778,2.731,2.397,2.136,1.926,1.754,1.611,1.486,1.378,1.284,1.202,1.128,1.061,1.0/
DATA EBETA/0.319,0.408,0.493,0.53,0.567,0.604,0.641,0.678,0.717,0.759,0.802,0.846,0.893,0.944,1.0/
EK=RX/RY
AA=0.5*(1./RX+1./RY)
BB=0.5*ABS(1./RX-1./RY)
CC=ACOS(BB/AA)*180.0/PAI
DO I=1,15
IF(CC.LT.THETA(I))THEN
WRITE(*,*)I
EAL=EALFA(I-1)+(CC-THETA(I))*(EALFA(I)-EALFA(I-1))/(THETA(I)-THETA(I-1))
EBE=EBETA(I-1)+(CC-THETA(I))*(EBETA(I)-EBETA(I-1))/(THETA(I)-THETA(I-1))
GOTO 1
ENDIF
ENDDO
1 EA=EAL*(1.5*W0/AA/E1)**(1./3.0)
EB=EBE*(1.5*W0/AA/E1)**(1./3.0)
PH=1.5*W0/(EA*EB*PAI)
OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
H00=0.0
MM=N-1
U=EDA0*US/(2.*E1*RX)
A1=ALOG(EDA0)+9.67
A2=5.1E-9*PH
A3=0.59/(PH*1.E-9)
B=PAI*PH*RX/E1
W=2.*PAI*PH/(3.*E1)*(B/RX)**2
ALFA=Z*5.1E-9*A1
G=ALFA*E1
AHM=1.0-EXP(-0.68*1.03)
HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)*AHM
ENDA=12.*U*(E1/PH)*(RX/B)**3
UTL=EDA0*US*RX/(B*B*2.E7)
W0=2.0*PAI*EA*EB*PH/3.0
WRITE(*,*)' Wait please'
CALL SUBAK(MM)
CALL MULTI(N,X0,XE,H00)
STOP
END
SUBROUTINE MULTI(N,X0,XE,H00)
DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500)
COMMON /COMEK/EK,EAL,EBE
DATA MK,G00/200,2.0943951/
G0=G00*EAL*EBE
NX=N
NY=N
NN=(N+1)/2
CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
M=0
14 KK=15
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
M=M+1
CALL ERP(N,ER,P,POLD)
WRITE(*,*)'ER=',ER
IF(M.LT.MK.AND.ER.GT.1.E-5)GOTO 14
CALL OUTPUT(N,DX,X,Y,H,P)
RETURN
END
SUBROUTINE ERP(N,ER,P,POLD)
DIMENSION P(N,N),POLD(N,N)
ER=0.0
SUM=0.0
NN=(N+1)/2
DO 10 I=1,N
DO 10 J=1,NN
ER=ER+ABS(P(I,J)-POLD(I,J))
SUM=SUM+P(I,J)
10 CONTINUE
ER=ER/SUM
DO I=1,N
DO J=1,N
POLD(I,J)=P(I,J)
ENDDO
ENDDO
RETURN
END
SUBROUTINE INITI(N,DX,X0,XE,X,Y,P,POLD)
DIMENSION X(N),Y(N),P(N,N),POLD(N,N)
NN=(N+1)/2
DX=(XE-X0)/(N-1.)
Y0=-0.5*(XE-X0)
DO 5 I=1,N
X(I)=X0+(I-1)*DX
Y(I)=Y0+(I-1)*DX
5 CONTINUE
DO 10 I=1,N
D=1.-X(I)*X(I)
DO 10 J=1,NN
C=D-Y(J)*Y(J)
IF(C.LE.0.0)P(I,J)=0.0
10 IF(C.GT.0.0)P(I,J)=SQRT(C)
DO 20 I=1,N
DO 20 J=NN+1,N
JJ=N-J+1
20 P(I,J)=P(I,JJ)
DO I=1,N
DO J=1,N
POLD(I,J)=P(I,J)
ENDDO
ENDDO
RETURN
END
SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
DIMENSION W(150,150),P0(150,150)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0/COMAK/AK(0:65,0:65)
COMMON /COMEK/EK,EAL,EBE
DATA NW,PAI,PAI1/150,3.14159265,0.2026423/
NN=(N+1)/2
CALL VI(NW,N,DX,P,W)
HMIN=1.E3
DO 30 I=1,N
DO 30 J=1,NN
RAD=X(I)*X(I)+EK*Y(J)*Y(J)
W1=0.5*RAD
H0=W1+W(I,J)
IF(H0.LT.HMIN)HMIN=H0
30 H(I,J)=H0
IF(KK.EQ.0)THEN
KG1=0
H01=-HMIN+HM0
DH=0.005*HM0
H02=-HMIN
H00=0.5*(H01+H02)
ENDIF
W1=0.0
DO 32 I=1,N
DO 32 J=1,N
32 W1=W1+P(I,J)
W1=DX*DX*W1/G0
DW=1.-W1
IF(KK.EQ.0)THEN
KK=1
GOTO 50
ENDIF
IF(DW.LT.0.0)THEN
KG1=1
H00=AMIN1(H01,H00+DH)
ENDIF
IF(DW.GT.0.0)THEN
KG2=2
H00=AMAX1(H02,H00-DH)
ENDIF
50 DO 60 I=1,N
DO 60 J=1,NN
H(I,J)=H00+H(I,J)
IF(P(I,J).LT.0.0)P(I,J)=0.0
EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
EDA(I,J)=EDA1
55 RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60 EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA1)
DO 70 J=NN+1,N
JJ=N-J+1
DO 70 I=1,N
H(I,J)=H(I,JJ)
RO(I,J)=RO(I,JJ)
EDA(I,J)=EDA(I,JJ)
70 EPS(I,J)=EPS(I,JJ)
RETURN
END
SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
DIMENSION D(70),A(350),B(210),ID(70)
COMMON /COM1/ENDA,A1,A2,A3,Z,C3/COMAK/AK(0:65,0:65)
DATA KG1,PAI1,C1,C2/0,0.2026423,0.31,0.31/
IF(KG1.NE.0)GOTO 2
KG1=1
AK00=AK(0,0)
AK10=AK(1,0)
AK20=AK(2,0)
BK00=AK00-AK10
BK10=AK10-0.25*(AK00+2.*AK(1,1)+AK(2,0))
BK20=AK20-0.25*(AK10+2.*AK(2,1)+AK(3,0))
2 NN=(N+1)/2
MM=N-1
DX1=1./DX
DX2=DX*DX
DX3=1./DX2
DX4=0.3*DX2
DO 100 K=1,KK
PMAX=0.0
DO 70 J=2,NN
J0=J-1
J1=J+1
JJ=N-J+1
IA=1
8 MM=N-IA
IF(P(MM,J0).GT.1.E-6)GOTO 20
IF(P(MM,J).GT.1.E-6)GOTO 20
IF(P(MM,J1).GT.1.E-6)GOTO 20
IA=IA+1
IF(IA.LT.N)GOTO 8
GOTO 70
20 IF(MM.LT.N-1)MM=MM+1
D2=0.5*(EPS(1,J)+EPS(2,J))
DO 50 I=2,MM
I0=I-1
I1=I+1
II=5*I0
D1=D2
D2=0.5*(EPS(I1,J)+EPS(I,J))
D4=0.5*(EPS(I,J0)+EPS(I,J))
D5=0.5*(EPS(I,J1)+EPS(I,J))
P1=P(I0,JJ)
P2=P(I1,JJ)
P3=P(I,JJ)
P4=P(I,JJ+1)
P5=P(I,JJ-1)
D3=D1+D2+D4+D5
IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
IF(H(I,J).LE.0.0)THEN
ID(I)=2
A(II+1)=0.0
A(II+2)=0.0
A(II+3)=1.0
A(II+4)=0.0
A(II+5)=1.0
A(II-4)=0.0
GOTO 50
ENDIF
IF(D1.GE.DX4)GOTO 30
IF(D2.GE.DX4)GOTO 30
IF(D4.GE.DX4)GOTO 30
IF(D5.GE.DX4)GOTO 30
ID(I)=1
IF(J.EQ.NN)P5=P4
A(II+1)=PAI1*(RO(I0,J)*BK10-RO(I,J)*BK20)
A(II+2)=DX3*(D1+0.25*D3)+PAI1*(RO(I0,J)*BK00-RO(I,J)*BK10)
A(II+3)=-1.25*DX3*D3+PAI1*(RO(I0,J)*BK10-RO(I,J)*BK00)
A(II+4)=DX3*(D2+0.25*D3)+PAI1*(RO(I0,J)*BK20-RO(I,J)*BK10)
A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
GOTO 50
30 ID(I)=0
P4=P(I,J0)
IF(J.EQ.NN)P5=P4
A(II+1)=PAI1*(RO(I0,J)*AK10-RO(I,J)*AK20)
A(II+2)=DX3*D1+PAI1*(RO(I0,J)*AK00-RO(I,J)*AK10)
A(II+3)=-DX3*D3+PAI1*(RO(I0,J)*AK10-RO(I,J)*AK00)
A(II+4)=DX3*D2+PAI1*(RO(I0,J)*AK20-RO(I,J)*AK10)
A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
50 CONTINUE
CALL TRA4(MM,D,A,B)
DO 60 I=2,MM
IF(ID(I).EQ.2)GOTO 60
IF(ID(I).EQ.0)GOTO 52
DD=D(I+1)
IF(I.EQ.MM)DD=0
P(I,J)=P(I,J)+C2*(D(I)-0.25*(D(I-1)+DD))
IF(J0.NE.1)P(I,J0)=P(I,J0)-0.25*C2*D(I)
IF(P(I,J0).LT.0.)P(I,J0)=0.0
IF(J1.GE.NN)GOTO 54
P(I,J1)=P(I,J1)-0.25*C2*D(I)
GOTO 54
52 P(I,J)=P(I,J)+C1*D(I)
54 IF(P(I,J).LT.0.0)P(I,J)=0.0
IF(PMAX.LT.P(I,J))PMAX=P(I,J)
60 CONTINUE
70 CONTINUE
DO 80 J=1,NN
JJ=N+1-J
DO 80 I=1,N
80 P(I,JJ)=P(I,J)
CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
100 CONTINUE
RETURN
END
SUBROUTINE TRA4(N,D,A,B)
DIMENSION D(N),A(5,N),B(3,N)
C=1./A(3,N)
B(1,N)=-A(1,N)*C
B(2,N)=-A(2,N)*C
B(3,N)=A(5,N)*C
DO 10 I=1,N-2
IN=N-I
IN1=IN+1
C=1./(A(3,IN)+A(4,IN)*B(2,IN1))
B(1,IN)=-A(1,IN)*C
B(2,IN)=-(A(2,IN)+A(4,IN)*B(1,IN1))*C
10 B(3,IN)=(A(5,IN)-A(4,IN)*B(3,IN1))*C
D(1)=0.0
D(2)=B(3,2)
DO 20 I=3,N
20 D(I)=B(1,I)*D(I-2)+B(2,I)*D(I-1)+B(3,I)
RETURN
END
SUBROUTINE VI(NW,N,DX,P,V)
DIMENSION P(N,N),V(NW,NW)
COMMON /COMAK/AK(0:65,0:65)
PAI1=0.2026423
DO 40 I=1,N
DO 40 J=1,N
H0=0.0
DO 30 K=1,N
IK=IABS(I-K)
DO 30 L=1,N
JL=IABS(J-L)
30 H0=H0+AK(IK,JL)*P(K,L)
40 V(I,J)=H0*DX*PAI1
RETURN
END
SUBROUTINE SUBAK(MM)
COMMON /COMAK/AK(0:65,0:65)
S(X,Y)=X+SQRT(X**2+Y**2)
DO 10 I=0,MM
XP=I+0.5
XM=I-0.5
DO 10 J=0,I
YP=J+0.5
YM=J-0.5
A1=S(YP,XP)/S(YM,XP)
A2=S(XM,YM)/S(XP,YM)
A3=S(YM,XM)/S(YP,XM)
A4=S(XP,YP)/S(XM,YP)
AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
10 AK(J,I)=AK(I,J)
RETURN
END
SUBROUTINE OUTPUT(N,DX,X,Y,H,P)
DIMENSION X(N),Y(N),H(N,N),P(N,N)
NN=(N+1)/2
A=0.0
WRITE(8,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(8,110)X(I),(H(I,J),J=1,N)
ENDDO
WRITE(10,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(10,110)X(I),(P(I,J),J=1,N)
ENDDO
110 FORMAT(66(E12.6,1X))
RETURN
END

回复此楼

» 本帖附件资源列表

欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
附件 1 : POINTEHL.f90

2013-07-27 10:37:59, 7.56 K

» 猜你喜欢

溴的反应液脱色已经有4人回复
参与限项已经有5人回复
推荐一本书已经有7人回复
有没有人能给点建议已经有5人回复
假如你的研究生提出不合理要求已经有12人回复
全日制（定向）博士已经有5人回复
萌生出自己或许不适合搞科研的想法，现在跑or等等看？已经有4人回复
Materials Today Chemistry审稿周期已经有4人回复
对氯苯硼酸纯化已经有3人回复
所感已经有4人回复

» 本主题相关价值贴推荐，对您同样有帮助:

fortran里一个子程序不运行是怎么回事已经有13人回复
fortran读写问题已经有6人回复
一百个FORTRAN任务怎么能一次提交SSH计算已经有3人回复
急求fortran运行错误原因，在线等已经有7人回复
请教 fortran 运行错误的原因已经有13人回复
FORTRAN新手求助主程序循环问题已经有10人回复
大家帮我看一个fortran的程序，我总是计算不出正确的结果已经有20人回复
请教一个fortran小程序编译出错的问题，谢谢已经有9人回复
如何用fortran编写高斯白噪声程序已经有6人回复
Fortran的格式化输入输出问题已经有14人回复
有关fortran的一次而问题，希望大家能帮帮忙，谢谢已经有4人回复
【求助】如何写FORTRAN程序实现求平均最近邻距离已经有34人回复
【求助】用fortran怎么描述下面这样的情况，我是新手，请大家帮忙已经有8人回复
【求助】初学fortran，已经有15人回复
【求助】请问一下fortran的可视化编程已经有7人回复
【交流】Fortran语言答疑专帖已经有157人回复

1楼 2013-07-27 10:40:50

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询

应助: 413 (硕士)
贵宾: 0.002
金币: 7116.5
散金: 15
红花: 63
帖子: 1639
在线: 798.9小时
虫号: 2469437
注册: 2013-05-14
专业: 计算数学与科学工程计算

【答案】应助回帖

感谢参与，应助指数 +1

怎么错了？错误信息是什么？你期待什么？

赞一下

回复此楼

2楼2013-07-28 08:37:29

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

黑曼巴1991

铁杆木虫 (正式写手)

应助: 4 (幼儿园)
金币: 7391.8
散金: 335
红花: 2
帖子: 930
在线: 292小时
虫号: 2552663
注册: 2013-07-18
专业: 传动机械学

引用回帖:

2楼: Originally posted by pippi6 at 2013-07-28 08:37:29
怎么错了？错误信息是什么？你期待什么？

直接不能运行没有输出结果。。几个DAT文件都是空的。错误信息发到附件里了麻烦看看，期待有结果，谢谢！

D0)FLF7U2802PMM}AX2{6WE.jpg

O05V}WLK6YAUI~~AW6BQU_U.jpg

Y4CBXDGY`CI6513HB{85[WL.jpg

赞一下

回复此楼

3楼2013-07-28 09:43:49

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询

应助: 413 (硕士)
贵宾: 0.002
金币: 7116.5
散金: 15
红花: 63
帖子: 1639
在线: 798.9小时
虫号: 2469437
注册: 2013-05-14
专业: 计算数学与科学工程计算

【答案】应助回帖

引用回帖:

3楼: Originally posted by 黑曼巴1991 at 2013-07-28 09:43:49
直接不能运行没有输出结果。。几个DAT文件都是空的。错误信息发到附件里了麻烦看看，期待有结果，谢谢！

D0)FLF7U2802PMM}AX2{6WE.jpg

O05V}WLK6YAUI~~AW6BQU_U.jpg

Y4CBXDGY`CI6513HB{85[WL.jpg
...

把几个do loop修改了一下，把标号去掉了。可以运算了。没出现你说的问题

PROGRAM POINTEHL
  DIMENSION THETA(15),EALFA(15),EBETA(15)
  COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
  COMMON /COMEK/EK,EAL,EBE
  DATA N,PAI,Z,EAL,EBE,E1,EDA0,RX,RY,X0,XE,W0,US/65,3.14159265,0.68,1.0,1.0,2.21E11,0.0283,0.05,0.05,-2.5,1.5,39.24,1.5/
  DATA THETA/10.,20.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90./
  DATA EALFA/6.612,3.778,2.731,2.397,2.136,1.926,1.754,1.611,1.486,1.378,1.284,1.202,1.128,1.061,1.0/
  DATA EBETA/0.319,0.408,0.493,0.53,0.567,0.604,0.641,0.678,0.717,0.759,0.802,0.846,0.893,0.944,1.0/
  EK=RX/RY
  AA=0.5*(1./RX+1./RY)
  BB=0.5*ABS(1./RX-1./RY)
  CC=ACOS(BB/AA)*180.0/PAI
  DO I=1,15
   IF(CC.LT.THETA(I))THEN
      WRITE(*,*)I
      EAL=EALFA(I-1)+(CC-THETA(I))*(EALFA(I)-EALFA(I-1))/(THETA(I)-THETA(I-1))
      EBE=EBETA(I-1)+(CC-THETA(I))*(EBETA(I)-EBETA(I-1))/(THETA(I)-THETA(I-1))
      GOTO 1
   ENDIF
  ENDDO
1 EA=EAL*(1.5*W0/AA/E1)**(1./3.0)
  EB=EBE*(1.5*W0/AA/E1)**(1./3.0)
  PH=1.5*W0/(EA*EB*PAI)
  OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
  OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
  OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
  WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
  WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
  H00=0.0
  MM=N-1
  U=EDA0*US/(2.*E1*RX)
  A1=ALOG(EDA0)+9.67
  A2=5.1E-9*PH
  A3=0.59/(PH*1.E-9)
  B=PAI*PH*RX/E1
  W=2.*PAI*PH/(3.*E1)*(B/RX)**2
  ALFA=Z*5.1E-9*A1
  G=ALFA*E1
  AHM=1.0-EXP(-0.68*1.03)
  HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)*AHM
  ENDA=12.*U*(E1/PH)*(RX/B)**3
  UTL=EDA0*US*RX/(B*B*2.E7)
  W0=2.0*PAI*EA*EB*PH/3.0
  WRITE(*,*)'             Wait please'
  CALL SUBAK(MM)
  CALL MULTI(N,X0,XE,H00)
  STOP
END PROGRAM POINTEHL

!flag
SUBROUTINE MULTI(N,X0,XE,H00)
  DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500)
  COMMON /COMEK/EK,EAL,EBE
  DATA MK,G00/200,2.0943951/
  G0=G00*EAL*EBE
  NX=N
  NY=N
  NN=(N+1)/2
  CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
  CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
  M=0
14 KK=15
  CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
  M=M+1
  CALL ERP(N,ER,P,POLD)
  WRITE(*,*)'ER=',ER
  IF(M.LT.MK.AND.ER.GT.1.E-5)GOTO 14
  CALL OUTPUT(N,DX,X,Y,H,P)
  RETURN
END SUBROUTINE MULTI

!flag
SUBROUTINE ERP(N,ER,P,POLD)
  DIMENSION P(N,N),POLD(N,N)
  ER=0.0
  SUM=0.0
  NN=(N+1)/2
  DO  I=1,N
   DO  J=1,NN
      ER=ER+ABS(P(I,J)-POLD(I,J))
      SUM=SUM+P(I,J)
   end DO
  end DO
  ER=ER/SUM
  DO I=1,N
   DO J=1,N
      POLD(I,J)=P(I,J)
   ENDDO
  ENDDO
  RETURN
END SUBROUTINE ERP

!flag
SUBROUTINE INITI(N,DX,X0,XE,X,Y,P,POLD)
  DIMENSION X(N),Y(N),P(N,N),POLD(N,N)
  NN=(N+1)/2
  DX=(XE-X0)/(N-1.)
  Y0=-0.5*(XE-X0)
  DO  I=1,N
   X(I)=X0+(I-1)*DX
   Y(I)=Y0+(I-1)*DX

  end DO

  DO I=1,N
   D=1.-X(I)*X(I)
   DO  J=1,NN
      C=D-Y(J)*Y(J)
      IF(C.LE.0.0)P(I,J)=0.0
      IF(C.GT.0.0)P(I,J)=SQRT(C)
   end DO
  end DO

  DO  I=1,N
   DO  J=NN+1,N
      JJ=N-J+1

      P(I,J)=P(I,JJ)
   ENDDO
  ENDDO

  DO I=1,N
   DO J=1,N
      POLD(I,J)=P(I,J)
   ENDDO
  ENDDO
  RETURN
END SUBROUTINE INITI

SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
  DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
  DIMENSION W(150,150),P0(150,150)
  COMMON /COM1/ENDA,A1,A2,A3,Z,HM0/COMAK/AK(0:65,0:65)
  COMMON /COMEK/EK,EAL,EBE
  DATA NW,PAI,PAI1/150,3.14159265,0.2026423/
  NN=(N+1)/2
  CALL VI(NW,N,DX,P,W)
  HMIN=1.E3
  DO  I=1,N
   DO  J=1,NN
      RAD=X(I)*X(I)+EK*Y(J)*Y(J)
      W1=0.5*RAD
      H0=W1+W(I,J)
      IF(H0.LT.HMIN)HMIN=H0
30    H(I,J)=H0
   ENDDO
  ENDDO

  IF(KK.EQ.0)THEN
   KG1=0
   H01=-HMIN+HM0
   DH=0.005*HM0
   H02=-HMIN
   H00=0.5*(H01+H02)
  ENDIF
  W1=0.0
  DO  I=1,N
   DO  J=1,N
32    W1=W1+P(I,J)
   ENDDO
  ENDDO
  W1=DX*DX*W1/G0
  DW=1.-W1
  IF(KK.EQ.0)THEN
   KK=1
   GOTO 50
  ENDIF
  IF(DW.LT.0.0)THEN
   KG1=1
   H00=AMIN1(H01,H00+DH)
  ENDIF
  IF(DW.GT.0.0)THEN
   KG2=2
   H00=AMAX1(H02,H00-DH)
  ENDIF

50 continue
  DO  I=1,N
   DO  J=1,NN
      H(I,J)=H00+H(I,J)
      IF(P(I,J).LT.0.0)P(I,J)=0.0
      EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
      EDA(I,J)=EDA1
55    RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60    EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA1)
   end DO
  end DO

  DO  J=NN+1,N
   JJ=N-J+1
   DO  I=1,N
      H(I,J)=H(I,JJ)
      RO(I,J)=RO(I,JJ)
      EDA(I,J)=EDA(I,JJ)
70    EPS(I,J)=EPS(I,JJ)
   end DO

  end DO

  RETURN
END SUBROUTINE HREE

SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
  DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N)
  DIMENSION D(70),A(350),B(210),ID(70)
  COMMON /COM1/ENDA,A1,A2,A3,Z,C3/COMAK/AK(0:65,0:65)
  DATA KG1,PAI1,C1,C2/0,0.2026423,0.31,0.31/
  IF(KG1.NE.0)GOTO 2
  KG1=1
  AK00=AK(0,0)
  AK10=AK(1,0)
  AK20=AK(2,0)
  BK00=AK00-AK10
  BK10=AK10-0.25*(AK00+2.*AK(1,1)+AK(2,0))
  BK20=AK20-0.25*(AK10+2.*AK(2,1)+AK(3,0))
2 NN=(N+1)/2
  MM=N-1
  DX1=1./DX
  DX2=DX*DX
  DX3=1./DX2
  DX4=0.3*DX2
  DO  K=1,KK
   PMAX=0.0
   DO  J=2,NN
      J0=J-1
      J1=J+1
      JJ=N-J+1
      IA=1
8    MM=N-IA
      IF(P(MM,J0).GT.1.E-6)GOTO 20
      IF(P(MM,J).GT.1.E-6)GOTO 20
      IF(P(MM,J1).GT.1.E-6)GOTO 20
      IA=IA+1
      IF(IA.LT.N)GOTO 8
      cycle
20    IF(MM.LT.N-1)MM=MM+1
      D2=0.5*(EPS(1,J)+EPS(2,J))
      DO  I=2,MM
         I0=I-1
         I1=I+1
         II=5*I0
         D1=D2
         D2=0.5*(EPS(I1,J)+EPS(I,J))
         D4=0.5*(EPS(I,J0)+EPS(I,J))
         D5=0.5*(EPS(I,J1)+EPS(I,J))
         P1=P(I0,JJ)
         P2=P(I1,JJ)
         P3=P(I,JJ)
         P4=P(I,JJ+1)
         P5=P(I,JJ-1)
         D3=D1+D2+D4+D5
         IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
         IF(H(I,J).LE.0.0)THEN
            ID(I)=2
            A(II+1)=0.0
            A(II+2)=0.0
            A(II+3)=1.0
            A(II+4)=0.0
            A(II+5)=1.0
            A(II-4)=0.0
            GOTO 50
         ENDIF
         IF(D1.GE.DX4)GOTO 30
         IF(D2.GE.DX4)GOTO 30
         IF(D4.GE.DX4)GOTO 30
         IF(D5.GE.DX4)GOTO 30
         ID(I)=1
         IF(J.EQ.NN)P5=P4
         A(II+1)=PAI1*(RO(I0,J)*BK10-RO(I,J)*BK20)
         A(II+2)=DX3*(D1+0.25*D3)+PAI1*(RO(I0,J)*BK00-RO(I,J)*BK10)
         A(II+3)=-1.25*DX3*D3+PAI1*(RO(I0,J)*BK10-RO(I,J)*BK00)
         A(II+4)=DX3*(D2+0.25*D3)+PAI1*(RO(I0,J)*BK20-RO(I,J)*BK10)
         A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
         GOTO 50
30       ID(I)=0
         P4=P(I,J0)
         IF(J.EQ.NN)P5=P4
         A(II+1)=PAI1*(RO(I0,J)*AK10-RO(I,J)*AK20)
         A(II+2)=DX3*D1+PAI1*(RO(I0,J)*AK00-RO(I,J)*AK10)
         A(II+3)=-DX3*D3+PAI1*(RO(I0,J)*AK10-RO(I,J)*AK00)
         A(II+4)=DX3*D2+PAI1*(RO(I0,J)*AK20-RO(I,J)*AK10)
         A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))
50       CONTINUE
      end DO
      CALL TRA4(MM,D,A,B)
      DO  I=2,MM
         IF(ID(I).EQ.2)GOTO 60
         IF(ID(I).EQ.0)GOTO 52
         DD=D(I+1)
         IF(I.EQ.MM)DD=0
         P(I,J)=P(I,J)+C2*(D(I)-0.25*(D(I-1)+DD))
         IF(J0.NE.1)P(I,J0)=P(I,J0)-0.25*C2*D(I)
         IF(P(I,J0).LT.0.)P(I,J0)=0.0
         IF(J1.GE.NN)GOTO 54
         P(I,J1)=P(I,J1)-0.25*C2*D(I)
         GOTO 54
52       P(I,J)=P(I,J)+C1*D(I)
54       IF(P(I,J).LT.0.0)P(I,J)=0.0
         IF(PMAX.LT.P(I,J))PMAX=P(I,J)
60       CONTINUE
      end DO
70    CONTINUE
   end do
   DO J=1,NN
      JJ=N+1-J
      DO I=1,N
80       P(I,JJ)=P(I,J)
      end do
   end do
   CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,EDA,P)
  end do
  RETURN
END SUBROUTINE ITER

SUBROUTINE TRA4(N,D,A,B)
  DIMENSION D(N),A(5,N),B(3,N)
  C=1./A(3,N)
  B(1,N)=-A(1,N)*C
  B(2,N)=-A(2,N)*C
  B(3,N)=A(5,N)*C
  DO  I=1,N-2
   IN=N-I
   IN1=IN+1
   C=1./(A(3,IN)+A(4,IN)*B(2,IN1))
   B(1,IN)=-A(1,IN)*C
   B(2,IN)=-(A(2,IN)+A(4,IN)*B(1,IN1))*C
10 B(3,IN)=(A(5,IN)-A(4,IN)*B(3,IN1))*C
  end DO
  D(1)=0.0
  D(2)=B(3,2)
  DO  I=3,N
20 D(I)=B(1,I)*D(I-2)+B(2,I)*D(I-1)+B(3,I)
  end DO
  RETURN
END SUBROUTINE TRA4

SUBROUTINE VI(NW,N,DX,P,V)
  DIMENSION P(N,N),V(NW,NW)
  COMMON /COMAK/AK(0:65,0:65)
  PAI1=0.2026423
  DO  I=1,N
   DO  J=1,N
      H0=0.0
      DO  K=1,N
         IK=IABS(I-K)
         DO  L=1,N
            JL=IABS(J-L)
30          H0=H0+AK(IK,JL)*P(K,L)
         end DO
      end DO
40    V(I,J)=H0*DX*PAI1
   end DO
  end DO
  RETURN
END SUBROUTINE VI

SUBROUTINE SUBAK(MM)
  COMMON /COMAK/AK(0:65,0:65)
  S(X,Y)=X+SQRT(X**2+Y**2)
  DO  I=0,MM
   XP=I+0.5
   XM=I-0.5
   DO J=0,I
      YP=J+0.5
      YM=J-0.5
      A1=S(YP,XP)/S(YM,XP)
      A2=S(XM,YM)/S(XP,YM)
      A3=S(YM,XM)/S(YP,XM)
      A4=S(XP,YP)/S(XM,YP)
      AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
10    AK(J,I)=AK(I,J)
   end DO
  end DO
  RETURN
END SUBROUTINE SUBAK

SUBROUTINE OUTPUT(N,DX,X,Y,H,P)
  DIMENSION X(N),Y(N),H(N,N),P(N,N)
  NN=(N+1)/2
  A=0.0
  WRITE(8,110)A,(Y(I),I=1,N)
  DO I=1,N
   WRITE(8,110)X(I),(H(I,J),J=1,N)
  ENDDO
  WRITE(10,110)A,(Y(I),I=1,N)
  DO I=1,N
   WRITE(10,110)X(I),(P(I,J),J=1,N)
  ENDDO
110 FORMAT(66(E12.6,1X))
  RETURN
END SUBROUTINE OUTPUT

赞一下

回复此楼

4楼2013-07-28 15:08:23

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

黑曼巴1991

铁杆木虫 (正式写手)

应助: 4 (幼儿园)
金币: 7391.8
散金: 335
红花: 2
帖子: 930
在线: 292小时
虫号: 2552663
注册: 2013-07-18
专业: 传动机械学

引用回帖:

4楼: Originally posted by pippi6 at 2013-07-28 15:08:23
把几个do loop修改了一下，把标号去掉了。可以运算了。没出现你说的问题

PROGRAM POINTEHL
  DIMENSION THETA(15),EALFA(15),EBETA(15)
  COMMON /COM1/ENDA,A1,A2,A3,Z,HM0
  COMMON /COMEK/EK,EAL,EBE
...

我用你的试了一下还是一样的问题啊。。

奇怪了是不是我用的FTN95的原因？麻烦你给我发一下修改后的.90文件？、万分感谢！！！