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

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

为什么我这个代码显示浮点异常呢?应该怎么改才行呢?
      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

我试了一下,可以了,谢谢你了。:-loveliness:
如果你方便的话,能帮我看一下,下面另一个程序么?
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
也许是你想要的。


我试了一下,可以了,谢谢你了。:-loveliness:
如果你方便的话,能帮我看一下,下面另一个程序么?
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

小姐姐你好厉害啊:-P
谢谢你了:-loveliness:

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

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

Kyra是對的, 佩服
页: [1]
查看完整版本: 代码运行显示the following floating-point exceptions are signalling.