代码运行显示the following floating-point exceptions are signalling.

sdulv · 发表于 2018-8-1 09:42:56

为什么我这个代码显示浮点异常呢？应该怎么改才行呢？

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

        PROGRAM POINTEHL
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
    DATA PAI,Z/3.14159265,0.68/
        DATA N,PH,E1,EDA0,RX,US,X0,XE/33,0.8E9,2.21E11,0.05,0.02,1.0,-2.5,1.5/
        OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
        OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
        OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
        MM=N-1
        A1=ALOG(EDA0)+9.67
        A2=5.1E-9*PH
        A3=0.59/(PH*1.E-9)
        U=EDA0*US/(2.*E1*RX)
        B=PAI*PH*RX/E1
        W0=2.*PAI*PH/(3.*E1)*(B/RX)**2
        ALFA=Z*5.1E-9*A1
        G=ALFA*E1
        HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W0**(-0.073)
        ENDA=12.*U*(E1/PH)*(RX/B)**3
        WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
        WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
        WRITE(*,*)'               Wait please'
        CALL SUBAK(MM)
        CALL EHL(N,X0,XE)
        STOP
        END
        SUBROUTINE EHL(N,X0,XE)
        DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500),V(4500)
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
        DATA MK,G0/1,2.0943951/
        CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
        KK=0
        CALL HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
14        KK=15
        CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
        MK=MK+1
        CALL ERP(N,ER,P,POLD)
        WRITE(*,*)'ER=',ER
        IF(ER.GT.1.E-5.AND.DH.GT.1.E-6)THEN
        IF(MK.GE.20)THEN
        MK=1
        DH=0.5*DH
        ENDIF
        GOTO 14
        ENDIF
        IF(DH.LE.1.E-6)WRITE(*,*)'Pressures are not convergent!!!'
        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
        DO 10 I=1,N
        DO 10 J=1,N
        ER=ER+ABS(P(I,J)-POLD(I,J))
        POLD(I,J)=P(I,J)
        SUM=SUM+P(I,J)
10        CONTINUE
        ER=ER/SUM
        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 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)
        POLD(I,J)=P(I,J)
        ENDDO
        ENDDO
        RETURN
        END
        SUBROUTINE HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
        DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COMAK/AK(0:65,0:65)
        DATA PAI,PAI1/3.14159265,0.2026423/
        NN=(N+1)/2
        CALL VI(N,DX,P,V)
        HMIN=1.E3
        DO 30 I=1,N
        DO 30 J=1,NN
        RAD=X(I)*X(I)+Y(J)*Y(J)
        W1=0.5*RAD
        H0=W1+V(I,J)
        IF(H0.LT.HMIN)HMIN=H0
30        H(I,J)=H0
        IF(KK.EQ.0)THEN
        KK=1
        DH=0.01*HM0
        H00=-HMIN+HM0
        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(DW.LT.0.0)H00=H00+DH
        IF(DW.GT.0.0)H00=H00-DH
        DO 60 I=1,N
        DO 60 J=1,NN
        H(I,J)=H00+H(I,J)
        EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
        EDA(I,J)=EDA1
        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,V)
        DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
        COMMON /COMAK/AK(0:65,0:65)
        DATA KG1,PAI/0,3.14159265/
        IF(KG1.NE.0)GOTO 2
        KG1=1
        AK00=AK(0,0)
        AK10=AK(1,0)
2        NN=(N+1)/2
        DO 100 K=1,KK
        DO 70 J=2,NN
        J0=J-1
        J1=J+1
        D2=0.5*(EPS(1,J)+EPS(2,J))
        DO 70 I=2,N-1
        I0=I-1
        I1=I+1
        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))
        D8=2.0*RO(I,J)*AK00/PAI**2
        D9=2.0*RO(I0,J)*AK10/PAI**2
        D10=D1+D2+D4+D5+D8*DX-D9*DX
        D11=D1*P(I0,J)+D2*P(I1,J)+D4*P(I,J0)+D5*P(I,J1)
        D12=(RO(I,J)*H(I,J)-D8*P(I,J)-RO(I0,J)*H(I0,J)+D9*P(I,J))*DX
        P(I,J)=(D11-D12)/D10
        IF(P(I,J).LT.0.0)P(I,J)=0.0
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,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
100        CONTINUE
        RETURN
        END
        SUBROUTINE VI(N,DX,P,V)
        DIMENSION P(N,N),V(N,N)
        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)
        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

sdulv · 发表于 2018-8-1 09:45:53

我的是win7系统，用的是simply fortran

kyra · 发表于 2018-8-1 11:08:33

ITER 函数里加一句
real , save :: AK00 , AK10
也许是你想要的。

顺便说，如果这代码是你自己书写的，那么你应该考虑换一本比较新的书，学习一下比较新的语法了。

sdulv · 发表于 2018-8-1 11:58:00

我试了一下，可以了，谢谢你了。

如果你方便的话，能帮我看一下，下面另一个程序么？
ITER 函数里的这一句real , save :: AK00 , AK10是按照你说的方法加的，但是它运行后还是出现浮点异常

PROGRAM ELLIPEHL
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
COMMON /COMEK/EK,BX,BY,PTR
REAL*8 RX,RY,KA,KB
DATA N,PAI,Z,E1,EDA0,RX,RY,X0,XE,W0,US/65,3.14159265,0.68,2.21E11,0.02,0.03,0.06,-2.5,1.5,1000.,3./
OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
EK=RX/RY
AA=0.5*(1./RX+1./RY)
BB=0.5*ABS(1./RX-1./RY)
CALL HERTZELLIPTIC(RX,RY,KA,KB)
EA=KA*(1.5*W0/AA/E1)**(1./3.0)
EB=KB*(1.5*W0/AA/E1)**(1./3.0)
PH=1.5*W0/(EA*EB*PAI)
WRITE(*,*)"N,X0,XE,W0,PH,E1,EDA0,RX,US=",N,X0,XE,W0,PH,E1,EDA0,RX,US
MM=N-1
U=EDA0*US/(E1*RX)
A1=ALOG(EDA0)+9.67
A2=5.1E-9*PH
A3=0.59/(PH*1.E-9)
BX=EB
BY=EA
IF(RX.GT.RY) THEN
BX=EA
BY=EB
ENDIF
W=W0/(E1*RX**2)
PTR=3*W*(RX/BY)*(RX/BX)**2/(PAI**2)
ALFA=Z*5.1E-9*A1
G=ALFA*E1
AHM=1.0-EXP(-0.68*1.03)
HM0=3.63*(RX/BX)**2*G**0.49*U**0.68*W**(-0.073)*AHM
ENDA=12.*U*(E1/PH)*(RX/BX)**3
WRITE(*,*)'             Wait please'
CALL SUBAK(MM)
CALL MULTI(N,X0,XE)
STOP
END
SUBROUTINE MULTI(N,X0,XE)
DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500),V(4500)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
COMMON /COMEK/EK,BX,BY,PTR
DATA MK,G00/1,2.0943951/
G0=G00*BY/BX
NX=N
NY=N
NN=(N+1)/2
CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
CALL HREE(N,DX,G0,X,Y,H,RO,EPS,EDA,P,V)
14 KK=15
CALL ITER(N,KK,DX,G0,X,Y,H,RO,EPS,EDA,P,V)
MK=MK+1
CALL ERP(N,ER,P,POLD)
WRITE(*,*)'ER=',ER
IF(ER.GT.1.E-5.AND.DH.GT.1.E-7)THEN
IF(MK.GE.10)THEN
WRITE(*,*)'ER,DH=',ER,DH
MK=1
DH=0.5*DH
ENDIF
GOTO 14
ENDIF
IF(DH.LE.1.E-7)WRITE(*,*)'Pressures are not convergent!!!'
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)
COMMON /COMEK/EK,BX,BY,PTR
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-(BX/BY)**2*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,G0,X,Y,H,RO,EPS,EDA,P,V)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COMAK/AK(0:65,0:65)
COMMON /COMEK/EK,BX,BY,PTR
DATA KK,PAI,PAI1/0,3.14159265,0.2026423/
NN=(N+1)/2
CALL VI(N,DX,P,V)
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+V(I,J)
IF(H0.LT.HMIN)HMIN=H0
30 H(I,J)=H0
W1=0.0
DO 40 I=1,N
DO 40 J=1,N
40 W1=W1+P(I,J)
W1=DX*DX*W1/G0
DW=1.-W1
IF(KK.NE.0)GOTO 50
KK=1
DH=0.1*HM0
H00=-HMIN+HM0
50 IF(DW.LT.0.0)H00=H00+DH
IF(DW.GT.0.0)H00=H00-DH
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
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,G0,X,Y,H,RO,EPS,EDA,P,V)
real , save :: AK00 , AK10
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
COMMON /COMAK/AK(0:65,0:65)
COMMON /COMEK/EK,BX,BY,PTR
DATA KG1,PAI/0,3.14159265/
IF(KG1.NE.0)GOTO 2
KG1=1
AK00=AK(0,0)
AK10=AK(1,0)
2 NN=(N+1)/2
DO 100 K=1,KK
DO 70 J=2,NN
J0=J-1
J1=J+1
D2=0.5*(EPS(1,J)+EPS(2,J))
DO 70 I=2,N-1
I0=I-1
I1=I+1
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))
D8=PTR*RO(I,J)*AK00
D9=PTR*RO(I0,J)*AK10
D10=D1+D2+D4+D5+D8*DX-D9*DX
D11=D1*P(I0,J)+D2*P(I1,J)+D4*P(I,J0)+D5*P(I,J1)
D12=(RO(I,J)*H(I,J)-D8*P(I,J)-RO(I0,J)*H(I0,J)+D9*P(I,J))*DX
P(I,J)=(D11-D12)/D10
IF(P(I,J).LT.0.0)P(I,J)=0.0
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,G0,X,Y,H,RO,EPS,EDA,P,V)
100 CONTINUE
RETURN
END
SUBROUTINE VI(N,DX,P,V)
DIMENSION P(N,N),V(N,N)
COMMON /COMEK/EK,BX,BY,PTR
COMMON /COMAK/AK(0:65,0:65)
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*PTR
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
SUBROUTINE HERTZELLIPTIC(RX,RY,KA,KB)
IMPLICIT NONE
REAL*8, EXTERNAL :: EE,KE
REAL*8 RX,RY,BPA,BMA,CTH,THT,PAI,E1,KA,KB
DATA PAI/3.1415926/
BPA=0.5*(1./RX+1./RY)
BMA=0.5*ABS(1./RX-1./RY)
CTH=BMA/BPA
THT=ACOS(CTH)*180.0/PAI
CALL CACUE(CTH,E1)
KA=(2.*EE(E1)/(PAI*(1-E1**2)))**(1/3.)
KB=KA*(1.-E1**2)**(1/2.)
END SUBROUTINE
SUBROUTINE CACUE(CTH,E1)
IMPLICIT NONE
REAL*8, EXTERNAL :: EE,KE
REAL*8, EXTERNAL :: FAB
INTEGER FLG,I
REAL*8 PAI,CTH,E1,E11,E12,DX,A,B,A1,A2,A3,A4,A5,T1,T2,T3,T4,T5,ER0
DATA PAI,DX,FLG,I,T1,T5,ER0/3.1415926,0.0001,1,1,1.E-30,1.,1.E-12/
IF(CTH.LT.1.E-6)THEN
E1=0.
RETURN
ENDIF
IF(CTH.GT.0.9999999999)THEN
E1=1.
RETURN
ENDIF
A1=FAB(T1,CTH)
A5=FAB(T5,CTH)
DO WHILE(FLG.EQ.1)
T3=T1+I*DX
A3=FAB(T3,CTH)
I=I+1
IF((A1*A3.LT.0.).AND.(A3*A5.LT.0.)) THEN
FLG=0
END IF
END DO
DO WHILE((T3-T1).GT.ER0)
T2=(T1+T3)/2.
A2=FAB(T2,CTH)
IF(A2.GT.0.) T1=T2
IF(A2.LT.0.) T3=T2
IF(A2.EQ.0.)THEN
E11=T2
EXIT
END IF
END DO
E11=T2
DO WHILE((T5-T3).GT.ER0)
T4=(T3+T5)/2.
A4=FAB(T4,CTH)
IF(A4.GT.0.) T5=T4
IF(A4.LT.0.) T3=T4
IF(A4.EQ.0.)THEN
E12=T2
EXIT
END IF
END DO
E12=T4
E1=E11
IF(E11.LT.E12) E1=E12
RETURN
END SUBROUTINE
REAL*8 FUNCTION FAB(E1,CTH)
IMPLICIT NONE
REAL*8 E1,CTH,T1,T2
REAL*8, EXTERNAL :: EE,KE
T1=EE(E1)
T2=KE(E1)
FAB=2*(1-E1**2)*(T1-T2)+(1.-CTH)*E1**2*T1
END FUNCTION
REAL*8 FUNCTION KE(E1)
IMPLICIT NONE
INTEGER N,I,FLG
REAL*8 E1,PAI,H,T,T1,T2,S1,S2,P,Q
PAI=3.1415926
IF(E1.EQ.1) THEN
KE=1.E10
RETURN
ENDIF
IF(E1.LT.1.E-20) THEN
KE=PAI/2.
RETURN
ENDIF
N=1
H=PAI/2.
Q=SQRT(1.-E1*E1*SIN(H)*SIN(H))
IF(Q.LT.1.E-35) Q=1.E35
Q=1./Q
T1=.5*H*(1+Q)
S1=T1
FLG=1
DO WHILE(FLG.EQ.1)
P=0.
DO I=0,N-1
T=(I+0.5)*H
Q=SQRT(1.-E1*E1*SIN(T)*SIN(T))
IF(Q.LT.1.E-35) Q=1.E35
Q=1./Q
P=P+Q
END DO
T2=(T1+H*P)/2.
S2=(4.*T2-T1)/3.
IF(ABS(S2-S1).LT.ABS(S2)*1.E-7) FLG=0
T1=T2
S1=S2
N=N+N
H=.5*H
END DO
KE=S2
RETURN
END FUNCTION
REAL*8 FUNCTION EE(E1)
IMPLICIT NONE
INTEGER N,I,FLG
REAL*8 E1,PAI,H,T,T1,T2,S1,S2,P,Q
PAI=3.1415926
N=1
H=PAI/2.
IF(E1.EQ.1) THEN
EE=1.
RETURN
ENDIF
IF(E1.LT.1.E-20) THEN
EE=PAI/2.
RETURN
ENDIF
Q=SQRT(1.-E1*E1*SIN(H)*SIN(H))
T1=.5*H*(1+Q)
S1=T1
FLG=1
DO WHILE(FLG==1)
P=0.
DO I=0,N-1
T=(I+0.5)*H
Q=SQRT(1.-E1*E1*SIN(T)*SIN(T))
P=P+Q
END DO
T2=(T1+H*P)/2.
S2=(4.*T2-T1)/3.
IF(ABS(S2-S1).LT.ABS(S2)*1.E-7) FLG=0
T1=T2
S1=S2
N=N+N
H=.5*H
END DO
EE=S2
RETURN
END FUNCTION

chiangtp · 发表于 2018-8-1 13:05:13

sdulv 发表于 2018-8-1 11:58
我试了一下，可以了，谢谢你了。
如果你方便的话，能帮我看一下，下面另一个程序么？
ITER 函 ...

你對"DATA"的認知可能有誤

SUBROUTINE ABC(...)
DATA MK, G00 /100, 2.0/

副程式"ABC"被呼叫多次, 如果每次進入"ABC", MK與G00的值都要是100與2.0
請將
DATA MK, G00 /100, 2.0/
改為
MK=100
G00=2.0

Fortran人識DATA與SAVE者幾多?

sdulv · 发表于 2018-8-1 16:05:32

chiangtp 发表于 2018-8-1 13:05
你對"DATA"的認知可能有誤

SUBROUTINE ABC(...)

我试了一下你的方法，直接在编译的时候就显示错误了。

sdulv · 发表于 2018-8-1 16:07:53

kyra 发表于 2018-8-1 11:08
ITER 函数里加一句
real , save :: AK00 , AK10
也许是你想要的。

我试了一下，可以了，谢谢你了。

如果你方便的话，能帮我看一下，下面另一个程序么？
ITER 函数里real , save :: AK00 , AK10是按照你说的方法加的，但是它运行后还是出现浮点异常

kyra · 发表于 2018-8-1 17:16:36

HERR函数里
real , save :: H00

sdulv · 发表于 2018-8-1 17:37:37

kyra 发表于 2018-8-1 17:16
HERR函数里
real , save :: H00

小姐姐你好厉害啊

谢谢你了

chiangtp · 发表于 2018-8-1 18:08:40

sdulv 发表于 2018-8-1 17:37
小姐姐你好厉害啊
谢谢你了

Kyra是對的, 佩服

[数值问题] 代码运行显示the following floating-point exceptions are signalling.

新人勋章

美女勋章

元老勋章

热心勋章

规矩勋章

管理勋章

贡献勋章

帅哥勋章

星光勋章