[求助] 显示数组溢出，实在找不到问题了，请大家帮忙

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



这是留的大作业，研究好久了，感觉不难，但做起来很繁琐，麻烦大家了，另外请问有没有调试程序的好方法，一个一个找太累了！！！