显示数组溢出，实在找不到问题了，请大家帮忙
program DVM!
implicit none !-->all the variables MUST be declared
!
integer,parameter :: n = 10
!Declaration of variable
real(8) :: U,a,dt
real(8) :: xs,ys
integer :: p,q
integer :: i
!p for which vortex it is, q for how long it exist
real(8),dimension(n-1):: gama
!the circulation of each vortex
real(8),dimension(n,n) :: xx,yy,xxi,yyi
!the location of the vortices and their corresponding images, the first index refers to which vortex it is, and the second one refers to how long it have existed
print*,'input the Velocity of the free stream, the Radius of the cylinder'
read*,U,a
print*,'input the Unit Time Step'
read*,dt
xs = 0
ys = a
!uu = U(1-((xx(p,q)**2 - yy(p,q)**2)*a**2)/((xx(p,q)**2 + yy(p,q)**2)**2))
!vv = -2*U*xx(p,q)*yy(p,q)*a**2/(xx(p,q)**2+yy(p,q)**2)**2
!formula (14) (15) in the lecture note
do p = 1,n-1
!set q = 1 to get the new pop-up vortex---------------------
i = 1
q = 1
!define the value of x of every pop-up vortex---------------
xx(p,q) = 0
!define the value of y of every pop-up vortex--------------
if (MOD(p,2)/=0) then
yy(p,q) = 1.1*a
else
yy(p,q) = -1.1*a
end if
!calculate gamma of the pop-up vortex----------------------
q = p
call getgama()
q = p
!calculate the new location of the vortices at present-----
do i = 1,p
call getlocation()
q = q - 1
end do
end do
print*,gama
contains
subroutine image(xxx,yyy,xxxi,yyyi)
implicit none
real(8) :: xxx,yyy,xxxi,yyyi
xxxi = a**2*(xxx/(xxx**2 + yyy**2))
yyyi = a**2*(yyy/(xxx**2 + yyy**2))
return
end subroutine image
subroutine getgama()
implicit none
real(8), dimension(p) :: simp
!to simplify the calculation, use simp to represent the iterative part of formula (18)
real(8), parameter :: Pi = 3.1415927
if(p==1) then
simp(1) = -2*Pi*U*(1-((xs**2 - ys**2)*(a**2))/((xs**2 + ys**2)**2))
call image(xx(1,1),yy(1,1),xxi(1,1),yyi(1,1))
gama(1) = simp(1)/(-(ys-yy(1,1))/((xs-xx(1,1))**2+(ys-yy(1,1))**2)+(ys-yyi(1,1))/((xs-xxi(1,1))**2+(ys-yyi(1,1))**2))
else
if (MOD(p,2)/=0) then
do while (i<p)
call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
simp(i+1) = simp(i)+gama(i)*(ys-yy(i,q))/((xs-xx(i,q))**2+(ys-yy(i,q))**2)-gama(i)*(ys-yyi(i,q))/((xs-xxi(i,q))**2+(ys-yyi(i,q))**2)
q = q-1
i = i+1
end do
gama(p) = simp(p)/(-(ys-yy(p,1))/((xs-xx(p,1))**2+(ys-yy(p,1))**2)+(ys-yyi(p,1))/((xs-xxi(p,1))**2+(ys-yyi(p,1))**2))
else
do while (i<p)
call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
simp(i+1) = simp(i)+gama(i)*(-ys-yy(i,q))/((xs-xx(i,q))**2+(-ys-yy(i,q))**2)-gama(i)*(-ys-yyi(i,q))/((xs-xxi(i,q))**2+(-ys-yyi(i,q))**2)
q = q-1
i = i+1
end do
gama(p) = simp(p)/(-(-ys-yy(p,1))/((xs-xx(p,1))**2+(-ys-yy(p,1))**2)+(-ys-yyi(p,1))/((xs-xxi(p,1))**2+(-ys-yyi(p,1))**2))
end if
end if
return
end subroutine getgama
subroutine getlocation()
implicit none
!p tells how many vortices do we have now
!we should move (p-1) vortices
!q tells which stage is each vortex at
real, parameter :: Pi = 3.1415927
real(8) :: Ustream,Vstream
real(8) :: simpxm = 0,simpxp = 0,simpym = 0,simpyp = 0
integer :: j,m
m = q
do j = 1,p
if(j/=i) then
q = p + 1 - j
simpxm = simpxm + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yy(j,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
end if
end do
do j = 1,p
q = p + 1 - j
simpxp = simpxp + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yyi(j,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
end do
do j = 1,p
if(j/=i) then
q = p + 1 - j
simpym = simpym + dt*( Gama(j) / (2*Pi) ) * (xx(j,q) - xx(i,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
end if
end do
do j = 1,p
q = p + 1 - j
simpyp = simpyp + dt*( Gama(j) / (2*Pi) ) * (xxi(j,q) - xx(i,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
end do
Ustream = U * (1-(xx(i,m)**2-yy(i,m)**2) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2)
Vstream = -2 * U * xx(i,m) * yy(i,m) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2
xx(i,m+1) = xx(i,m) + Ustream * dt - simpxm + simpxp
yy(i,m+1) = yy(i,m) + Vstream * dt - simpym + simpyp
return
end subroutine getlocation
end program DVM
这是留的大作业，研究好久了，感觉不难，但做起来很繁琐，麻烦大家了，另外请问有没有调试程序的好方法，一个一个找太累了！！！
因为不知道你的算法，所以只是猜想。
主程序里 43 到 49 行，q 在循环中变小。每次循环 q = q - 1
但是 getlocation 函数里，又改变了 q，比如 q = p + 1 - j，这样很危险，不是吗？ 我觉得问题也许是在subroutine getgama中。一般我的经验是按步骤测试，首先去掉sub，运行主程序，然后一步一步加入call。加入getgamma后，print*, gama 没有正确显示数值(出现NaN)。
纯属个人意见，仅供参考。
页:
[1]