[数值问题] 程序运行出现NAN

program ex
                implicit none
        integer E,P,I,J,K,PP,NE,NP,NDT,L,MDT,ICOUNT,IPRINT,IJK(20,3),MP(3)
        REAL MU,DT,DF,TXX,TYY,AREA,MPP,D,TIME,DFMAX,OLDVAL,SUMGST,X(22),Y(22),H0(22),H(22),BB(22),B(3),C(3),G(22,22),ST(22,22)
        NE=20
        NP=18
        NDT=20
        DT=0.5
        TXX=10.0
        TYY=10.0
        MU=0.005
        OPEN(1,FILE="XZB.txt")
        READ(1,*)(X(I),I=1,NP)
        CLOSE(1)
        OPEN(1,FILE="YZB.txt")
        READ(1,*)(Y(I),I=1,NP)
        CLOSE(1)
        OPEN(1,FILE="IJKDY.txt")
        READ(1,*)((IJK(E,P),P=1,3),E=1,NE)
        CLOSE(1)
        WRITE(*,5)(X(I),I=1,NP)
   5    FORMAT(1X,"X(NP)=",/,10(9F5.1,/))
        WRITE(*,6)(Y(I),I=1,NP)
   6    FORMAT(1X,"Y(NP)=",/,10(9F5.1,/))
        WRITE(*,10)((IJK(E,P),P=1,3),E=1,NE)
   10   FORMAT(1X,"IJK(NE,3)=",/,100(4(3I4,2X),/))

                DO 20 P=1,NP
                H0(P)=16.
                20   H(P)=16.
                DO 25 P=16,18
                H0(P)=11.
   25   H(P)=11.
        DO 30 P=1,NP
                DO 30 PP=1,NP
                ST(P,PP)=0.0
   30   G(P,PP)=0.0
        DO 100 E=1,NE
                I=IJK(E,1)
                J=IJK(E,2)
                K=IJK(E,3)
                B(1)=Y(J)-Y(K)
                B(2)=Y(K)-Y(I)
                B(3)=Y(I)-Y(J)
                C(1)=X(K)-X(J)
                C(2)=X(I)-X(J)
                C(3)=X(J)-X(I)
                AREA=ABS(B(2)*C(3)-B(3)*C(2))/2.
                MP(1)=I
                MP(2)=J
                MP(3)=K
                DO 50 P=1,3
                L=MP(P)
                G(L,I)=G(L,I)+(TXX*B(P)*B(1)+TYY*C(P)*C(1))/(4.*AREA)
                G(L,J)=G(L,J)+(TXX*B(P)*B(2)+TYY*C(P)*C(2))/(4.*AREA)
                G(L,K)=G(L,K)+(TXX*B(P)*B(3)+TYY*C(P)*C(3))/(4.*AREA)
                DO 40 PP=1,3
                MPP=MP(PP)
                IF(L.EQ.MPP)GO TO 36
                D=12.0
                GO TO 38
   36   D=6.0
   38   ST(L,MPP)=ST(L,MPP)+MU*AREA/D
   40   CONTINUE
   50   CONTINUE
   100  CONTINUE
        WRITE(*,110)
   110  FORMAT(/,5X,"TIME",31X,"HEAD")
        IPRINT=2
                ICOUNT=1
                TIME=DT
                DO 500 MDT=1,NDT
                DO 150 L=1,NP
                BB(L)=0.0
                DO 150 P=1,NP
   150  BB(L)=BB(L)+ST(L,P)*H0(P)/DT

   200  DFMAX=0.0
        DO 400 L=1,NP
                IF(L.EQ.1.OR.L.EQ.2.OR.L.EQ.3.OR.L.EQ.16.OR.L.EQ.17.OR.L.EQ.18) GOTO 400
                OLDVAL=H(L)
                SUMGST=0.0
                DO 300 P=1,NP
                IF(P.EQ.L) GO TO 300
                SUMGST=SUMGST+(G(L,P)+ST(L,P)/DT)*H(P)
   300  CONTINUE
        H(L)=(BB(L)-SUMGST)/(G(L,L)+ST(L,L)/DT)
                DF=ABS(OLDVAL-H(L))
                IF(DF.GT.DFMAX) DFMAX=DF
   400  CONTINUE
        IF(DFMAX.GT.0.01) GO TO 490
                DO 450 L=1,NP
   450  H0(L)=H(L)
        IF(ICOUNT.NE.IPRINT) GO TO 490
                WRITE(*,460)TIME,(H(I),I=1,NP)
   460  FORMAT(1X,F8.2,1X,9F7.2,/,10X,9F7.2)
        ICOUNT=0
   490  TIME=TIME+DT
        ICOUNT=ICOUNT+1
   500  CONTINUE
        STOP
                END