Fortran Coder

查看: 21288|回复: 9

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

[复制链接]

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
发表于 2018-8-1 09:42:56 | 显示全部楼层 |阅读模式
为什么我这个代码显示浮点异常呢?应该怎么改才行呢?
[Fortran] 纯文本查看 复制代码
        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

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
 楼主| 发表于 2018-8-1 09:45:53 | 显示全部楼层
我的是win7系统,用的是simply fortran

235

帖子

0

主题

0

精华

版主

World Analyser

F 币
630 元
贡献
464 点

新人勋章美女勋章元老勋章热心勋章规矩勋章管理勋章

QQ
发表于 2018-8-1 11:08:33 | 显示全部楼层
ITER 函数里加一句
real , save :: AK00 , AK10
也许是你想要的。

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

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
 楼主| 发表于 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

130

帖子

10

主题

0

精华

大师

F 币
617 元
贡献
372 点

贡献勋章管理勋章帅哥勋章元老勋章星光勋章规矩勋章

发表于 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者幾多?

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
 楼主| 发表于 2018-8-1 16:05:32 | 显示全部楼层
chiangtp 发表于 2018-8-1 13:05
你對"DATA"的認知可能有誤

SUBROUTINE ABC(...)

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

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
 楼主| 发表于 2018-8-1 16:07:53 | 显示全部楼层
kyra 发表于 2018-8-1 11:08
ITER 函数里加一句
real , save :: AK00 , AK10
也许是你想要的。

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

               

235

帖子

0

主题

0

精华

版主

World Analyser

F 币
630 元
贡献
464 点

新人勋章美女勋章元老勋章热心勋章规矩勋章管理勋章

QQ
发表于 2018-8-1 17:16:36 | 显示全部楼层
HERR函数里
real , save :: H00

19

帖子

4

主题

0

精华

入门

F 币
115 元
贡献
57 点
 楼主| 发表于 2018-8-1 17:37:37 | 显示全部楼层
kyra 发表于 2018-8-1 17:16
HERR函数里
real , save :: H00

小姐姐你好厉害啊
谢谢你了

130

帖子

10

主题

0

精华

大师

F 币
617 元
贡献
372 点

贡献勋章管理勋章帅哥勋章元老勋章星光勋章规矩勋章

发表于 2018-8-1 18:08:40 | 显示全部楼层
sdulv 发表于 2018-8-1 17:37
小姐姐你好厉害啊
谢谢你了

Kyra是對的, 佩服
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-3-19 11:34

Powered by Tencent X3.4

© 2013-2024 Tencent

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