[請益] 幫忙Debug, P幣1000

看板Fortran作者wop (nothing)時間10年前 (2014/05/02 22:20)推噓1(1推 0噓 9→)

留言10則, 4人參與討論串1/1

以下的程式碼執行後會出現以下錯誤訊息取對數時的定義域有問題, 煩請各位先進指點迷津小弟感恩不盡 !!! 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=1,MAX-1 DO 102 I=1,MAX-1 DO 102 J=1,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(S)**2.)) XJK=(XNS**2.)*S*V1*(AA**0.5)*EXP(BB) CC=0.6666667*THETA/DLOG(S) XKSTAR=CC**3 DPSTAR=4.*SIGMA*V1/(XKB*T*DLOG(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 -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.113.223.47 ※ 文章網址: http://www.ptt.cc/bbs/Fortran/M.1399040419.A.28F.html