phi(2,1,1)=0

call calculation (xp(1,1,ix,iy),x(1,ix,iy,1),phi(1,ix,iy),phi(2,ix,iy)) 
[color=red]此时 ix，iy 均为 1
所以，phi(2,ix,iy) = phi(2,1,1) = 0[/color]

subroutine calculation (x,xc,phi1,phi2) [color=red]此时，phi2 =0[/color]

f(x)=1.06366481d10+7.94800172d9+1.72992961d9*log(abs(x))-1.72992961d9*log(abs(1-x))-1.182988286d10*(1-2*x)+2.75714286d9*(6*x-6*x**2-1)+3.083731429d9*(24*x**2-10*x-16*x**3+1)-1d11*(xc-phi1*x)/phi2 + 9.229804981d10

[color=red]phi2 出现在分母上，分母为 0，非法了。[/color]

谢谢楼上的回复。

我是fortran的初学者，现在用别人的程序。但我现在需要用newton-raphson解一个非线性方程，例如：lnx+x^3+A=0, 这里我要解的未知数x和常量A都是一个数组，像x是这样一个数组：(1,1,100,100), A是（1,100,100）. 我只能从网上找到关于解x不是数组的code,如下：
subroutine calculation (x,A)
integer ::k
real(8) :: x,x0,A
real(8):: tol

f(x)=lnx+x**3+A
df(x)=1/x+3*x**2
tol=1e-6
x0=0.2d0
k=0
10 k=k+1
x=x0-f(x0)/df(x0)
if(abs(x-x0).lt.tol)goto 20
x0= x
goto 10
20 write(*,*) x
end do
end do

end subroutine calculation

但是对于数组来说应该怎么改呢？我想把这个作为一个subroutine。非常感谢各位的帮组！

如果能解决，可以有尝的，呵呵，谢谢！

试试slatec包

http://people.sc.fsu.edu/~jburkardt/f_src/slatec/slatec.html

里边有例程SOS (FNC, NEQ, X, RTOLX, ATOLX, TOLF, IFLAG, RW, LRW, &
IW, LIW)和例程SNSQE (FCN, JAC, IOPT, N, X, FVEC, TOL, NPRINT, INFO, &
WA, LWA)可以试试

SOS好像用的牛顿法解非线性方程组，NSQE用的Powell方法解方程组。

双精度版本，可以试试DSOS和DNSQE，源代码slatec.f90里NSQE有调用的例子，你自己看看。应该可以

谢谢楼上的回复！
不过这个对我来说太难了啊。有没有什么简单的啊，谢谢！