Question

Fortran语言中call语句在do循环中不跟着循环是怎么回事？求大神赐教！！

program w
real z,o,v
do a1=0.0,1.0,0.1

z=3*a1+0.1
call q
v=answer
o=sqrt((z-v)**2)
if(o<0.000000001) then
print*,a1,z,answer
stop
else
print*,'没有',z,a1,v,o
endif
enddo
end program

module GSLD !高斯—勒让德积分高斯点及权重的求解模块

implicit none
integer,parameter::n=5 !设置求解高斯点的个数
integer,parameter::nn=16 !设置kind的数值
contains

real(kind=nn) function f(x) !定义四阶精度的函数f（x）

implicit none
integer::i
real(kind=nn)::a(n),x !a(n)代表n阶勒让德多项式
a(1)=x !1阶勒让德多项式
a(2)=1.5_16*(x**2)-0.5_16 !2阶勒让德多项式
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
f=a(n) !生成的n阶勒让德多项式
end function

real(kind=nn) function f1(x)
implicit none
integer::i
real(kind=nn)::a(n),x
a(1)=x
a(2)=1.5_16*x**2-0.5_16
do i=3,n-1
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
f1=a(n-1) !生成的（n-1）阶勒让德多项式
endfunction

real(kind=nn) function g(x)
implicit none
integer::i
real(kind=nn)::a(n),x
a(1)=x
a(2)=1.5_16*x**2-0.5_16
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
g=n*a(n-1)/(1-x**2)-n*x*a(n)/(1-x**2)
endfunction

real(kind=nn) function bis(a,b) !二分法求解函数的解
implicit none
real(kind=nn)::a,b,c
do while(.true.)
c=(a+b)/2.0_16
if(f(c)*f(a)<0)then
b=c
else
a=c
endif
if((b-a)<1E-34)exit !这个精度设置的。。（有点狠了）
enddo
bis=c !bis即是利用二分法求得的解
end function

subroutine fn0(fn,ak)
implicit none
real(kind=nn)::m,fn(n),ak(n)
integer::i,j
j=0 !赋值控制循环变量的初值
m=-1.000001 !设置计算域[-1，1] 的下限，即代替-1
do i=1,200000 !这个循环次数应该是由步长0.00001决 定,计算方法：200000=2/0.000001
if(f(m)*f(m+0.00001)<0)then !从下限处开始往上逐步累加，
j=j+1 !记录这是第几个解
fn(j)=bis(m,m+0.00001)
ak(j)=2.0_16/(n*f1(fn(j))*g(fn(j))) !高斯点的权重
write(*,*) '高斯点序号',j
endif
m=m+0.00001 !执行完一次判断m向前推进一步
enddo
endsubroutine

real(kind=nn) function y(x) !被积函数
implicit none
real(kind=nn)::x,a2

!y=1/((sin(x*2*10**(9))/(x**2*4)-(cos(x*2*10**(9))/(x*2*10**(9)))))
y=x+a2+0.1

end function
end module

subroutine q
use GSLD
implicit none
real(kind=nn)::fn(n),ak(n),x,a,b,answer,a2
integer::i
call fn0(fn,ak) !调用求高斯零点和权重的子函数
a=0.0 !积分上限
b=1.0 !积分下限
answer=0.0_16 !求积分结果赋初始值
do i=1,n
answer=answer+ak(i)*y((a+b)/2.0_16+(b-a)/2.0_16*fn(i))
enddo
answer=answer*(b-a)/2.0_16
print*,answer
end

Answer 1

GRKT1(Y,W,F,D)是个子程序，Y,W,F,D分别为子程序的参数，call为调用它