一个小程序哪位帮忙修改一下

草木

是数论网格法计算三重积分的。计算网格部分没问题，是从书上抄来的，是基本程序语句有问题。当H=R=W=1时，积分真值是0.011858263……
EXTERNAL SLF
         (M,I,N,N0,A,C,D,X,F,GW)
DIMENSION A(I),C(I),D(I),X(I),F(M),GW(M)
DOUBLE PRECISION F,GW,U
COMMON/DIM/M,I,N,N0
COMMON/DIM/H,R,W
M=1
I=3
N=100063
N0=N
A(1)=1
A(2)=28036
A(3)=22431
WRITE(*,*) 'LEASE INPUT H,R,W VALUE:';
READ(5,*) H,R,W
    WRITE(*,*) 'LEASE INPUT 积分下限 C(1),C(2),C(3):';
READ(5,*) C(1),C(2),C(3)
WRITE(*,*) 'LEASE INPUT 积分上限 D(1),D(2),D(3):';
READ(5,*) D(1),D(2),D(3)
      DO 10 J=1,M
10   GW(J)=0.0
CONTINUE
      DO 20 K=1,N
G=1.0
CONTINUE
DO 40 L=1,I
U0=FLOAT(K)*A(L)/FLOAT(N0)
U=U0-IFIX(U0)
V=(C(L)-D(L))*U
W=(U-1.0)*V
G=6.0*W*G
40X(L)=(W+W-V)*U+C(L)
CALL FCT(X,F)
      CONTINUE
DO 20 J=1,M
20       GW(J)=G*F(J)/FLOAT(N)+GW(J)
RETURN
WRITE(*,*) ';GW= ';, GW(J)
      FORMAT(1X,';GW=';,D16.9)
END
SUBROUTINE FCT(X,F)
DOUBLE PRECISION X,F
F=((X(1)+X(2)+X(3))*H*R*W)**3
F=1.0/F
RETURN
END

草木

草木

积分区间漏了，X1, X2, X3的积分区间都是[1，2]