代码运行后显示 The following floating-point exceptions are signalling

sdulv · 发表于 2018-11-14 15:38:43

   CHARACTER*16 FILENAME
   CHARACTER*1 S,S1,S2
   CHARACTER*16 CDATE,CTIME
   CHARACTER(len=8)::ccDATE
   COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,HM0
   DATA MK,PAI,Z/6,3.14159265,0.68/,S1,S2/1HY,1Hy/,
   & N,PH,E1,EDA0,RX,US,X0,XE,C1,C2/65,0.5E9,2.21E11,0.05,
   & 0.02,1.0,-2.5,1.5,0.31,0.31/,
   & CDATE,CTIME/'The date is','The time is'/
   OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
   OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
   OPEN(10,FILE='PRESS.DAT',STATUS='UNKNOWN')
1    FORMAT(20X,A12,I2.2,':',I2.2,':',I4.4)
2    FORMAT(20X,A12,I2.2,':',I2.2,':',I2.2,'.',I2.2)
   WRITE(*,*)'Show the example or not (Y or N)?'
   READ(*,'(A)')S
   IF(S.EQ.S1.OR.S.EQ.S2)GOTO 5
   WRITE(*,*)'N,PH,E1,EDA0,RX,US,X0,XE='
   READ(*,*)N,PH,E1,EDA0,RX,US,X0,XE
   WRITE(*,*)' Change iteration factors or not (Y or N) ?'
   READ(*,'(A)')S
   IF(S.EQ.S1.OR.S.EQ.S2)THEN
   WRITE(*,*)'C1,C2='
   READ(*,*)C1,C2
   ENDIF
5    CALL DATE_AND_TIME(ccDATE)
   READ(ccDATE,"(i4,2i2)")IYR,IMON,IDAY
   WRITE(4,1)CDATE,IMON,IDAY,IYR
   WRITE(*,1)CDATE,IMON,IDAY,IYR
   WRITE(*,20)N,X0,XE,PH,E1,EDA0,RX,US,C1,C2
   WRITE(4,20)N,X0,XE,PH,E1,EDA0,RX,US,C1,C2
20    FORMAT(/////,
   &  15X,'+--------+--------+--------+--------+--------+',/,
   &  15X,'| N | X0 | XE | W(N/M) |  E(Pa) |',/,
   &  15X,'+--------+--------+--------+--------+--------+',/,
   &  15X,'|  ',I4,'  | ',F6.3,' | ',F6.3,' |',E8.3,'|',E8.3,'|',/,
   &  15X,'+--------+--------+--------+--------+--------+',/,
   &  15X,'|Eda(PaS)|  R (M) | Us(M/S)| C1 | C2 |',/,
   &  15X,'+--------+--------+--------+--------+--------+',/,
   &  15X,'| ',F6.3,' | ',F6.3,' | ',F6.3,' | ',F6.3,' | ',F6.3,'|',/,
   &  15X,'+--------+--------+--------+--------+--------+',/)
   H00=0.0
   MM=N-1
   LMIN=ALOG(N-1.)/ALOG(2.)-1.99
   LMAX=LMIN
   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
   AM=W*(2.*U)**(-0.75)
   AL=G*(2.*U)**(0.25)
   HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)
   ENDA=12.*U*(E1/PH)*(RX/B)**3
   WRITE(*,*)'             Wait please  '
   CALL SUBAK(MM)
   CALL MULTI(N,LMIN,LMAX,MK,X0,XE,H00)
   call DATE_AND_TIME(ccDATE)
   READ(ccDATE,"(i4,2i2)")IYR,IMON,IDAY
   WRITE(4,1)CDATE,IMON,IDAY,IYR
   WRITE(*,1)CDATE,IMON,IDAY,IYR
   STOP
   END

sdulv · 发表于 2018-11-14 15:39:24

C**************************************************************
   SUBROUTINE MULTI(N,LMIN,LMAX,MK,X0,XE,H00)
   DIMENSION X(100),Y(100),H(5000),RO(5000),R0(5000),
   & EPS(5000),P(5000),F(5000),R(5000),G(20)
   COMMON /COMK/KL
   DATA NMAX,G0/100,2.0943951/
   NN=(N+1)/2
   KL=LMIN
   CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
   CALL INITI(N,X,Y,P(NS))
12    CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),0)
14    DO 100 L=LMIN,LMAX
   KL=L
   G(KL)=G0
20    KK=2
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
   KK=1
   CALL CHAN(N,R0,R(NS))
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R0,1)
   CALL GCAL(N,G1,P(NS))
   G(KL-1)=G(KL)-DX*DX*G1
   N2=N
   NS0=NS
   KL=KL-1
   CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
   CALL TRANS(N,N2,N2*N2,H,RO,EPS,P(NS),P(NS0),R(NS),R0)
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),2)
   CALL GCAL(N,G1,P(NS))
   G(KL)=G(KL)+DX*DX*G1
   G0=G(KL)
   CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),1)
   IF(KL.NE.1)GOTO 20
   KK=15
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
40    NS0=NS
   KL=KL+1
   N1=N
   CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
   G0=G(KL)
   CALL TRAP(N1,N,P(NS0),P(NS))
   CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),0)
   KK=1
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
   IF(KL.LT.L)GOTO 40
   CALL CHAN(N,R0,R(NS))
   CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R0,1)
100 CONTINUE
   M=M+1
   IF(M.LT.MK)GOTO 14
120 CALL OUPT(N,X,Y,H,P(NS))
   RETURN
   END
C**************************************************************
   SUBROUTINE OUPT(N,X,Y,H,P)
   DIMENSION X(N),Y(N),H(N,N),P(N,N)
   NN=(N+1)/2
   A=0.0
   WRITE(8,*)N,NN
   WRITE(8,110)A,(Y(I),I=1,N)
   DO I=1,N
   WRITE(8,110)X(I),(H(I,J),J=1,N)
   END DO
   WRITE(10,*)N,NN
   WRITE(10,110)A,(Y(I),I=1,N)
   DO I=1,N
   WRITE(10,110)X(I),(P(I,J),J=1,N)
   END DO
110 FORMAT(34(1X,E12.6))
   RETURN
   END
C**************************************************************
   SUBROUTINE TRAP(N1,N2,P1,P2)
   DIMENSION P1(N1,N1),P2(N2,N2)
   DO 10 I=1,N1-1
   II=2*I
   DO 10 J=1,N1-1
   JJ=2*J
10    P2(II,JJ)=P2(II,JJ)+0.25*(P1(I,J)+P1(I+1,J)+
   & P1(I+1,J+1)+P1(I,J+1)-P2(II-1,JJ-1)-P2(II-1,JJ+1)-
   & P2(II+1,JJ+1)-P2(II+1,JJ-1))
   DO 20 I=1,N1-1
   II=2*I
   DO 20 J=1,N1
   JJ=2*J-1
   P2(II,JJ)=P2(II,JJ)+0.5*(P1(I,J)+P1(I+1,J)-
   & P2(II-1,JJ)-P2(II+1,JJ))
20    P2(JJ,II)=P2(JJ,II)+0.5*(P1(J,I)+P1(J,I+1)-
   & P2(JJ,II-1)-P2(JJ,II+1))
   DO 30 I=1,N1
   II=2*I-1
   DO 30 J=1,N1
   JJ=2*J-1
30    P2(II,JJ)=P1(I,J)
   RETURN
   END
C**************************************************************
   SUBROUTINE INITI(N,X,Y,P)
   DIMENSION X(N),Y(N),P(N,N)
   NN=(N+1)/2
   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)
   RETURN
   END
C**************************************************************
   SUBROUTINE CHAN(N,R1,R2)
   DIMENSION R1(N,N),R2(N,N)
   DO 10 I=1,N
   DO 10 J=1,N
10    R1(I,J)=R2(I,J)
   RETURN
   END
C**************************************************************
   SUBROUTINE GCAL(N,G,P)
   DIMENSION P(N,N)
   G=0.
   DO 10 I=1,N
   DO 10 J=1,N
10    G=G+P(I,J)
   RETURN
   END
C**************************************************************
   SUBROUTINE KNDX(K,N,NS,NMAX,DX,X0,XE,X,Y)
   DIMENSION INDEX(6),NNDEX(6),X(NMAX),Y(NMAX)
   DATA INDEX,NNDEX/1,82,371,1460,5685,16616,
   & 9,17,33,65,129,257/
   N=NNDEX(K)
   NS=INDEX(K)
   DX=(XE-X0)/(N-1.)
   Y0=-0.5*(XE-X0)
   DO 10 I=1,N
   X(I)=X0+(I-1)*DX
10    Y(I)=Y0+(I-1)*DX
   RETURN
   END
C**************************************************************
   SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P,F,KG)
   DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),F(N,N)
   DIMENSION W(150,150),P0(150,150)

   COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,HM0/COMK/KL
   & /COMAK/AK(0:100,0:100)
   DATA KK,NW,PAI1/0,150,0.2026423/
   NN=(N+1)/2
   CALL VI(N,DX,P,W,NW)
   HMIN=1.E3
   DO 30 I=1,N
   W1=0.5*X(I)*X(I)
   DO 30 J=1,NN
   H0=W1+0.5*Y(J)*Y(J)+W(I,J)
   IF(KG.EQ.1)GOTO 20
   IF(H0+F(I,J).LT.HMIN)HMIN=H0+F(I,J)
   H(I,J)=H0
   GOTO 30
20    F(I,J)=H(I,J)-H00-H0
   F(I,N-J+1)=F(I,J)
30    CONTINUE
   IF(KK.EQ.0)H00=HM0-HMIN
   KK=1
   IF(KG.EQ.1)RETURN
   H0=H00+HMIN
   IF(H0.LE.0.0)GOTO 40
   IF(KL.NE.1)GOTO 50
   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
   HM0=H0*W1
40    H00=-HMIN+HM0
50    DO 60 I=1,N
   DO 60 J=1,NN
   IF(P(I,J).LT.0.0)P(I,J)=0.0
   EDA=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
   H(I,J)=H00+H(I,J)+F(I,J)
   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*EDA)
   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)
70    EPS(I,J)=EPS(I,JJ)
   RETURN
   END
C**************************************************************
   SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P,F,R,KG)
   DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),
   & EPS(N,N),F(N,N),R(N,N)

   DIMENSION D(200),A(1000),B(600),ID(200)
   COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,C3
   & /COMAK/AK(0:100,0:100)
   DATA KG1,PAI1/0,0.2026423/
   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
   DO 70 J=2,NN
   J0=J-1
   J1=J+1
   JJ=N-J+1
   IF(KG.NE.0)GOTO 20
   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(KG.NE.0)GOTO 32
   IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
   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))+R(I,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))+R(I,J)
   GOTO 50
32    IF(KG.EQ.2)GOTO 40
   R(I,J)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+
   & DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
   GOTO 50
40    R(I,J)=DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)-
   & DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
50    CONTINUE
   IF(KG.NE.0)GOTO 70
   CALL TRA4(MM,D,A,B)
   DO 54 I=2,MM
   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
70    CONTINUE
   IF(KG.NE.0)GOTO 100
   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,P,F,0)
100 CONTINUE
   RETURN
   END

sdulv · 发表于 2018-11-14 15:39:54

C**************************************************************
   SUBROUTINE TRANS(N1,N2,NMAX,H,RO,EPS,P1,P2,R1,R2)
   DIMENSION H(NMAX),RO(NMAX),EPS(NMAX),
   & P1(N1,N1),P2(N2,N2),R1(N1,N1),R2(N2,N2)
   DO 10 I=1,N1
   II=2*I-1
   NI1=(I-1)*N1
   NI2=2*(I-1)*N2
   DO 10 J=1,N1
   JJ=2*J-1
   NJ1=NI1+J
   NJ2=NI2+JJ
   H(NJ1)=H(NJ2)
   RO(NJ1)=RO(NJ2)
   EPS(NJ1)=EPS(NJ2)
   P1(I,J)=P2(II,JJ)
10    R1(I,J)=R2(II,JJ)
   RETURN
   END
C**************************************************************
   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
C**************************************************************
SUBROUTINE VI(N,DX,P,V,NW)
DIMENSION P(N,N),V(NW,NW)
COMMON /COMAK/AK(0:100,0:100)
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
C**************************************************************
   SUBROUTINE SUBAK(MM)
   COMMON /COMAK/AK(0:100,0:100)
   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
C*********************************************************