[数值问题] run-time error M6021:MATH - log:DOMAIN error

      PROGRAM NGDE
       IMPLICIT DOUBLE PRECISION(A-H,O-Z)

       PARAMETER (MAX=42)
       PARAMETER (MAX1=42)

         DIMENSION V(MAX1+1),XN(MAX1),X(MAX1,MAX1,MAX1),COLF(MAX1,MAX1), &
        XN1S(MAX1),ZETA(MAX1),DP(MAX1),XM(MAX1)

       PI=3.1415926
       XKB=1.38D-23
       R=8.314
       XNA=6.02225D23

       T=1773.D0
       COOLRATE=1000.D0
       P=1.
       XMW=26.95D-3
       RHO=2700.D0
       A=948.D0
       B=0.202
       C=13.07
       D=36373.D0
       AMU=1.966D-6
       BMU=147.47

       Q=10.**(12./(MAX-2))
       V1=XMW/(RHO*XNA)

       V(1)=V1

       DO 100 I=1,MAX-1
           XN(I)=0.
           XN1S(I)=0.
 100     CONTINUE

       DO 101 I=2,MAX
           V(I)=V(1)*(Q**(I-1))


           DP(I)=(6.*V(I)/PI)**0.333333
           XM(I)=RHO*V(I)
 101     CONTINUE

       DP(1)=(6.*V(1)/PI)**0.333333
         XM(1)=RHO*V(1)

         DO 102 K=2,MAX-1
         DO 102 I=2,MAX-1
         DO 102 J=2,MAX-1

           IF (V(K) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K+1)) THEN
                    X(I,J,K)=(V(K+1)-V(I)-V(J))/(V(K+1)-V(K))

         ELSE IF (V(K-1) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K)) THEN
                    X(I,J,K)=(V(I)+V(J)-V(K-1))/(V(K)-V(K-1))

                 ELSE
                    X(I,J,K)=0.D0
           END IF
 102     CONTINUE

       TIME=0.D0
       STEP=1.D-4
       S=1.001
       TEMP1=6.*V1/PI
       S1=PI*(TEMP1**0.66666666667)
       XMASS=RHO*V1
       PS=EXP(C-D/T)*101325.*P
       XNS=PS/(XKB*T)
       XN(1)=S*XNS


 200   FORMAT(  'Time',5X,'Temperature',5X,'XN(1)',5X,'XJK',5X,'S',5X,'Diameter',5X,'Particle volume',5X,'Total number')
       WRITE(*,200)

       ICOUNTER=1
       ICOUNTER2=1

         DO WHILE (T .GT. 1600)
         SIGMA=(A-B*T)*1.D-3
         PS=EXP(C-D/T)*101325.*P
         XNS=PS/(XKB*T)
         S=XN(1)/XNS
         THETA=(S1*SIGMA)/(XKB*T)
         AA=(2.*SIGMA)/(PI*XMASS)
         BB=THETA-(4.*(THETA**3.))/(27.*(DLOG(DABS(S))**2.))
         XJK=(XNS**2.)*S*V1*(AA**0.5)*EXP(BB)
         CC=0.6666667*THETA/DLOG(DABS(S))
         XKSTAR=CC**3
         DPSTAR=4.*SIGMA*V1/(XKB*T*DLOG(DABS(S)))
         VSTAR=PI*(DPSTAR**3.)/6.

         XNTOT=0.D0
           DO I=2,MAX-1
             XNTOT=XNTOT+XN(I)
           END DO

         VTOT=0.D0
         VTOTAL=0.D0
           DO I=2,MAX-1
             VTOT=VTOT+XN(I)*V(I)

           END DO
             VTOTAL=VTOT+XN(1)*V(1)

         VAV=VTOT/XNTOT
         DPAV=(6.*VAV/PI)**0.3333333

           IF (XNTOT .GT. 100.) STEP=1.D-5
         IF (TIME .GT. 0.14) STEP=5.D-8

 201     FORMAT(8(1X,E14.8))
 202     FORMAT(8(1X,E14.8))
 203     FORMAT(/'time=',E14.8)
 204     FORMAT('N',I2,2(2X,E14.8))
         OPEN(11,FILE='output1.DAT',STATUS='UNKNOWN',ACCESS='APPEND')
         OPEN(12,FILE='output2.DAT',STATUS='UNKNOWN',ACCESS='APPEND')


         ICOUNTER=ICOUNTER+1
         ICOUNTER2=ICOUNTER2+1

           IF (ICOUNTER .EQ. 400) THEN
                 WRITE(*,201)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
           WRITE(11,202)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
           ICOUNTER=1
         END IF

           IF (ICOUNTER2 .EQ. 40000) THEN
           WRITE(*,203)TIME
           WRITE(12,203)TIME

            DO I=1,MAX-1
              WRITE(*,204)I,DP(I),XN(I)
              WRITE(12,204)I,DP(I),XN(I)
            END DO
           ICOUNTER2=1
         END IF

         XKMIN=1.D-9
         XMU=AMU*(T**1.5)/(BMU+T)
         XLAMBDA=(XMU/(P*101325.))*SQRT(PI*R*T/(2.*0.04))
         DO 103 I=1,MAX-1
           DO 103 J=1,MAX-1
           XKN1=(2.*XLAMBDA)/DP(I)
           XKN2=(2.*XLAMBDA)/DP(J)
           D1=(XKB*T)/(3.*PI*XMU*DP(I))*((5.+4.*XKN1+6.*(XKN1**2)+18.*(XKN1**3))/(5.-XKN1+(8.+PI)*(XKN1**2)))
           D2=(XKB*T)/(3.*PI*XMU*DP(J))*((5.+4.*XKN2+6.*(XKN2**2)+18.*(XKN2**3))/(5.-XKN2+(8.+PI)*(XKN2**2)))
           C1=SQRT((8.*XKB*T)/(PI*XM(I)))
           C2=SQRT((8.*XKB*T)/(PI*XM(J)))
           XL1=(8.*D1)/(PI*C1)
           XL2=(8.*D2)/(PI*C2)
           G1=((DP(I)+XL1)**3)-((DP(I)**2+XL1**2)**1.5)/(3.*DP(I)*XL1)-DP(I)
           G2=((DP(J)+XL2)**3)-((DP(J)**2+XL2**2)**1.5)/(3.*DP(J)*XL2)-DP(J)
           COLF(I,J)=2.*PI*(D1+D2)*(DP(I)+DP(J))/((DP(I)+DP(J))/(DP(I)+DP(J)+2.*SQRT(G1**2+G2**2))+(8.*(D1+D2))/(SQRT(C1**2+C2**2)*(DP(I)+DP(J))))


           IF (COLF(I,J) .LT. XKMIN) THEN
             XKMIN=COLF(I,J)
           END IF

 103       CONTINUE


         DO 104 K=2,MAX-1
           IF (VSTAR .LT. V1) THEN
                ZETA(2)=VSTAR/V(2)
               ELSE IF (V(K-1) .LE. VSTAR .AND. VSTAR .LT. V(K)) THEN
                ZETA(K)=VSTAR/V(K)
             ELSE
                ZETA(K)=0.D0
              END IF
 104       CONTINUE


         DO 105 I=1,MAX-1
           DP(I)=(6.*V(I)/PI)**0.333333
           XN1S(I)=XNS*EXP(4.*SIGMA*XMW/(R*T*RHO*DP(I)))
 105       CONTINUE

         T1=0.D0
         T2=0.D0
         T3=0.D0
         T4=0.D0

         DO 106 K=2,MAX-1
           SUM1=0.D0
           SUM2=0.D0

             IF (STEP .NE. 1.D-4) THEN
             IF (XN(1) .GT. XN1S(K-1)) THEN

               IF (K .EQ. 2) THEN
                 ADDTERM=0.D0
                 ELSE
                 ADDTERM=(V(1)/(V(K)-V(K-1)))*COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
                 T1=T1+COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
               END IF
               END IF

             IF (XN(1) .LT. XN1S(K+1)) THEN
                 ADDTERM=-(V(1)/(V(K+1)-V(K)))*COLF(1,K+1)*(XN(1)-XN1S(K+1))*XN(K+1)
                 T2=T2+COLF(1,K+1)*(-XN(1)+XN1S(K+1))*XN(K+1)
               END IF

             IF (XN(1) .LT. XN1S(K)) THEN
                 SUBTERM=-(V(1)/(V(K)-V(K-1)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
                 T3=T3+COLF(1,K)*(-XN(1)+XN1S(K))*XN(K)
               END IF

             IF (XN(1) .GT. XN1S(K)) THEN
                 SUBTERM=(V(1)/(V(K+1)-V(K)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
                 T4=T4+COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
               END IF
           END IF
           DO I=2,MAX-1
               SUM2=SUM2+COLF(K,I)*XN(I)
              DO J=2,K
                SUM1=SUM1+X(I,J,K)*COLF(I,J)*XN(I)*XN(J)
              END DO
           END DO
           XN(K)=XN(K)+STEP*(0.5*SUM1-XN(K)*SUM2+XJK*ZETA(K)+ADDTERM-SUBTERM)

 106       CONTINUE
         XN(1)=XN(1)-XJK*XKSTAR*STEP-STEP*(T1+T4-T3-T2)
         T=T-STEP*COOLRATE
         TIME=TIME+STEP


         CLOSE(11)
           CLOSE(12)

       END DO

         STOP
         END