小弟这有个程序，不知道哪位大侠能帮忙改正！

program main
implicit none
  real,parameter :: width=120.0,heighth=120.0                           !模型大小
  double precision r                                                   !种子数
  real :: totarea=0.0                                                   !多边形的总面积
  real :: area=0.0                                                       !单个多边形的面积
  real,parameter :: e=0.6                                              !面积比 
  integer :: number=1                                                   !记录多边形的编号
  real :: xx(9,10000),yy(9,10000)
  integer Edges(10000)
  common xx,yy,Edges

  open(unit=10,file='point.txt')
  write(*,*) "Please input the number of seed!"
  read(*,*) r

  do while(totarea<e*heighth*width)
    call point(r,width,heighth,area,number)
    totarea=totarea+area
  end do
  close(10)

  stop
end program

!产生多边形的顶点
subroutine point(r,width,heighth,area,number)
implicit none
  double precision r
  real width,heighth
  real area
  real,external :: NRND1                                                !产生随机数的子程序
  real theta,beta                                                       !定义N边形的角度和随机旋转的角度
  real,parameter :: pi=3.1415926
  integer i                                                               !循环用
  integer n                                                                  !多边形边长
  real,parameter :: Nmax=8.0,Nmin=3.0                                   !岩石块的最大,最小边数
  real,parameter :: Bmax=20.0,Bmin=2.0                                   !岩石块的最大,最小半径
  real x0,y0                                                           !岩石块的圆心坐标
  real,allocatable :: x(:),y(:)                                           !暂时存放生成的多边形顶点坐标
  real num                                                               !存放生成的随机数
  real angle                                                           !顶点的角度
  real lenth                                                           !多边形半径的长
  integer number                                                       !记录多边形的编号

  if(allocated(x))deallocate(x)
  if(allocated(y))deallocate(y)
  num=NRND1(r)
  n=nint(num*(Nmax-Nmin)+Nmin)                                           !生成多边形的边的数目
  theta=2*pi/n
  allocate(x(n+1),y(n+1))
  num=NRND1(r)
  x0=num*width                                                           !多边形圆心的x坐标值
  num=NRND1(r)
  y0=num*heighth                                                       !多边形圆心的y坐标值
  beta=NRND1(r)                                                           !多边形第一条边偏转的角度
  
  do i=1,n
    angle=beta+theta*(i-1)
    num=NRND1(r)
    lenth=num*(Bmax-Bmin)+Bmin                                           !多边形的半径
    x(i)=x0+lenth*cos(angle)
    y(i)=y0+lenth*sin(angle)
  end do
  x(n+1)=x(1)
  y(n+1)=y(1)

  call calculate(x0,y0,x,y,n,r,width,heighth,area,number)
  call checkin(x,y,n,width,heighth,r,area,number)
  deallocate(x,y)
  
return
end subroutine

!产生随机数的子程序        
real function NRND1(R)
double precision S,U,V,R
S=65536.0
U=2053.0
V=13849.0
M=R/S
R=R-M*S
R=U*R+V
M=R/S
R=R-M*S
NRND1=R/S
return
end    function

!判断多边形顶点是否在模型内
subroutine checkin(x,y,n,width,heighth,r,area,number)
implicit none
  integer n,i 
  real x(n+1),y(n+1)
  real width,heighth
  double precision r
  real area
  integer number

  do i=1,n
    if(x(i)<=0 .or. x(i)>=width)then                     
      area=0.0
      call point(r,width,heighth,area,number)
    end if

    if(y(i)<=0 .or. y(i)>=heighth)then                     
      area=0.0
      call point(r,width,heighth,area,number)
    end if
  end do
  call overlap(r,x,y,number,n,width,heighth,area)

return
end subroutine

!验证是否重叠的子程序
subroutine overlap(r,x,y,number,n,width,heighth,area)
implicit none
   double precision r
   integer n,number
   real x(n+1),y(n+1)
   real xx(9,10000),yy(9,10000)
   integer Edges(10000)
   common xx,yy,Edges
   real t(9)
   integer i,j,k,m,num
   real p,q
   real xxx(8)
   real width,heighth,area

   if(number>1)then
        do i=1,n
       do j=1,number-1
            m=Edges(j)
         num=0
         p=0.0
         do k=1,m
           t(k)=yy(k,j)-y(i)
         end do
         t(m+1)=t(1)
         
         do k=1,m
           if(t(k)*t(k+1)<0)then
                num=num+1
             xxx(num)=(xx(k+1,j)*y(i)-xx(k,j)*y(i)-xx(k+1,j)*yy(k,j)+xx(k,j)*yy(k+1,j))/(yy(k+1,j)-yy(k,j))
           end if
         end do

         if(num>0)then
           do k=1,num
             if(xxx(k)>x(i))then
               p=p+1.0
             end if
           end do
         end if
         
         if(mod(p,2.0) .neqv. 0)then
           area=0.0
           call point(r,width,heighth,area,number)
         end if
       end do
     end do
     
     do i=1,number-1
       m=Edges(i)
       do j=1,m
         num=0
         q=0.0
         do k=1,n
           t(k)=y(k)-yy(j,i)
         end do
         t(n+1)=t(1)

         do k=1,n
           if(t(k)*t(k+1)<0)then
                num=num+1
             xxx(num)=(x(k+1)*yy(j,i)-x(k)*yy(j,i)-x(k+1)*y(k)+x(k)*y(k+1))/(y(k+1)-y(k))
           end if
         end do

         if(num>0)then
           do k=1,num
             if(xxx(k)>xx(j,i))then
                  q=q+1.0
             end if
           end do
         end if
         
         if(mod(q,2.0) .neqv. 0)then
           area=0.0
           call point(r,width,heighth,area,number)
         end if
       end do
     end do

   end if

   Edges(number)=n
   write(10,*) "Number=",number
   do i=1,n
     xx(i,number)=x(i)
     yy(i,number)=y(i)
     write(10,*) i,x(i),y(i)
   end do
   xx(n+1,number)=xx(1,number)
   yy(n+1,number)=yy(1,number)
   write(10,*)
   number=number+1

return
end subroutine

!计算多边形面积
subroutine calculate(x0,y0,x,y,n,r,width,heighth,area,number)
implicit none
  integer n,i
  real x0,y0
  real s(n)                                                    !暂时存放多边形各小三角形的面积
  real area                                                    !多边形面积
  real x(n+1),y(n+1)
  double precision r
  real width,heighth
  integer number

  area=0.0
  do i=1,n
    s(i)=x0*y(i)+x(i)*y(i+1)+x(i+1)*y0-x(i+1)*y(i)-x0*y(i+1)-x(i)*y0
    area=area+abs(s(i))
  end do
  
  if(area>30 .and. area<1500)then
    return
  else
    area=0.0
    call point(r,width,heighth,area,number)
  end if
  
return
end subroutine

用debug检查发现，程序运行几次之后就出错了