向前辈们请教fortran运行之后错误的原因！

scq86320

小弟刚开始学习fortran编程，最近按例子的提示写了一段有关Stokes流动的程序，在改了很多错误之后开始运行时出现了 ：run-time error M6101:MATH
                                                          -floating-point error: invalid
我把源程序贴出来，请前辈们指点迷津！！
PROGRAM STKS
   CHARACTER SYMBOL(7),GRAPH(70)
     DIMENSION PSI(11,11),PSI1(11,11),ZETA(11,11),ZETA1(11,11),RHS(11,11)
      COMMON H2,M1
     DATA L,L1,U,EPSLON,ERRIN,M,NX/1.0,1.0,1.0,0.001,0.001,11,70/
     DATA SYMBOL/'0',' ','3',' ','5',' ','9'/
       M1=M-1
       H=L/M1
       H2=H*H
       ITER=0
       DO 10 I=1,M
       DO 10 J=1,M
       PSI(I,J)=0.0
    10          ZETA(I,J)=0.0
    20          ITER=ITER+1
                IF (ITER .GT. 30) GOTO 90
       DO 30 J=2,M1
       ZETA(1,J)=(-PSI(2,J-1)+8.0/3.0*PSI(2,J)-PSI(2,J+1)-2.0/3.0*PSI(3,J))/H2
       ZETA(M,J)=(-PSI(M1,J-1)+8.0/3.0*PSI(M1,J)-PSI(M1,J+1)-2.0/3.0*PSI(M1-1,J))/H2
    30       CONTINUE
                DO 40 I=2,M1
    40       ZETA(I,1)= (-PSI(I-1,2)+8.0/3.0*PSI(I,2)-PSI(I+1,2)-2.0/3.0*PSI(I,3))/H2+2.0*U/3.0/H
                DO 50 I=1,M
                      DO 50 J=1,M
                      RHS(I,J)=0.0
       50       ZETA1(I,J)=ZETA(I,J)
                      CALL INITER(ZETA,RHS,M,ERRIN)
                      ERRZ=0.0
                      DO 60 I=1,M
                      DO 60 J=1,M
       60       ERRZ=ERRZ+ABS(ZETA(I,J)-ZETA1(I,J))
                      DO 70 I=1,M
                      DO 70 J=1,M
                      RHS(I,J)=ZETA(I,J)
    70       PSI1(I,J)=PSI(I,J)
                      CALL INITER(PSI,RHS,M,ERRIN)
                      ERRP=0.0
                      DO 80 I=2,M1
                      DO 80 J=2,M1
    80       ERRP=ERRP+ABS(PSI(I,J)-PSI1(I,J))
                      IF (ERRZ .GT. EPSLON .OR. ERRP .GT. EPSLON) GOTO 20
          90       ZETA(1,1)=(ZETA(1,2)+ZETA(2,1))/2
                      ZETA(M,1)=(ZETA(M1,1)+ZETA(M,2))/2
                      ZETA(M,M)=(ZETA(M,M1)+ZETA(M1,M))/2
                      ZETA(1,M)=(ZETA(2,M)+ZETA(1,M1))/2
         100        FORMAT(///20X,'STREAM FUNCTION BASED ON',I3'ITERATION'//33X,'ERRP=',F7.5/)
                     DO 110 K=1,M
                            J=M-K+1
         110         WRITE(6,140)(PSI(I,J),I=1,M)
                       WRITE(6,120)ERRZ
         120        FORMAT(1HI///28X,'VORTICITY DISTRIBUTION'//32X,'ERRZ='F7.5/)
                      DO 130 K=1,M
                             J=M-K+1
         130        WRITE(6,140)(ZETA(I,J),I=I,M)
         140        FORMAT(/1X,11F7.3/)
                      WRITE(6,100)ITER,ERRP
   150            NY=0.6*NX*L1/L
                        DX=L/NX
                        DY=L1/NY
                       DLT=0.015
                     DO 170 JSMBL=1,NY
                                Y=L1-(JSMBL-0.5)*DY
                                J=1+Y/H
                               Y1=Y-(J-1)*H
                               Y2=H-Y1
                     DO 160 ISMBL=1,NX
                                X=(ISMBL-0.5)*DX
                                I=1+X/H
                               X1=X-(I-1)*H
                               X2=H-X1
                               A1=X2*Y2
                               A2=X2*Y1
                               A3=X1*Y1
                               A4=X1*Y2
                               PSIO=(A1*PSI(I,J)+A2*PSI(I,J+1)+A3*PSI(I+1,J+1)+A4*PSI(I+1,J))/H2
                               NR=1+ABS(PSIO)/DLT
                      GRAPH(ISMBL)=SYMBOL(NR)
         160         CONTINUE
         170         WRITE(6,180) (GRAPH(I),I=1,NX)
         180         FORMAT(4X,70A1)
                       WRITE(6,190)
         190         FORMAT(//25X,'FIGURE STREAM FUNCTION PLOT.')
                      STOP
                        END
                     SUBROUTINE INITER(F,Q,M,ERR)
                       DIMENSION F(M,M),Q(M,M)
                         COMMON H2,M1
                                ITERP=0
   200           ERRPZ=0.0
                       ITERP=ITERP+1
                     IF (ITERP.GT.50) GOTO 220
                       DO 210 I=2, M1
                       DO 210 J=2, M1
                             F1=F(I, J)
                        F(I, J)=0.25*(F(I-1,J)+F(I+1,J)+F(I,J-1)+F(I,J+1)+H2*Q(I,J))
         210           ERRPZ=ERRPZ+ABS(F(I,J)-F1)
                     IF (ERRPZ .GT. ERR) GOTO 200
         220        RETURN
                       END

5623728

L,L1,默认是整数变量，要改改。