高手帮我指点下这个程序，编译没有错误，执行的时候有错误。初次涉猎子程序，问题应该是出在了主程序和第一个子程序上。第二个子程序确定没问题。急急急！！！！！！！

 dimension v(10,10)
 real dd
 xs=0
 zs=40
 xr=200
 zr=80
 x0=0
   z0=0
 x_e=200
 z_e=200
 nx=6    
 nz=6
 dx=40
 dz=40
 dd=0.5
 k=1
 do i=1,nx
       do j=1,nz
        v(i,j)=2000
       enddo
 enddo
 qq=3.14159
 call subroutine aaa(xs,zs,xr,zr,t)
 print*,t
 end
ccccccccccccccccccccccccccccccccccccccccccccc
      subroutine aaa(xs,zs,xr,zr,t)
 real yy1,yy2
10    call ray1(x0,z0,x_e,z_e,nz,nx,dx,dz,xs,zs,qq,
     *      v,xxx,yyy,ttt,xxx2,yyy2,ttt2)
      yy1=(yyy2-yyy)*(xr-xxx)/(xxx2-xxx)+yyy+dd
  yy2=(yyy2-yyy)*(xr-xxx)/(xxx2-xxx)+yyy-dd
  if(zr.ge.yy2.and.zr.le.yy1)    then
         t=sqrt(((xxx-xr)**2+(yyy-zr)**2)*(ttt2-ttt)**2/
  *     ((xxx2-xxx)**2+(yyy2-yyy)**2))+ttt  
      else
     k=k+1
    if(zr.lt.yy2.or.xxx.le.200.and.yyy.ge.200)    then 
       qq=qq+3.1416/2**k      
            goto 10
    elseif(zr.gt.yy1.or.xxx.le.200.and.yyy.le.0)   then 
            qq=qq-3.1416/2**k 
            goto 10
         endif
      endif
 return
  end
ccccccccccccccccccccccccccccccccccccccccccccc
      subroutine ray1(x0,z0,x_e,z_e,nz,nx,dx,dz,xu,zu,qq,
     *     vel,xxx,yyy,ttt,xxx2,yyy2,ttt2)
      real vel(nx,nz),x_ray(90000),z_ray(90000)
 real t(9000),pp(9000),ppx(9000),ppz(9000)
 kkk=0
 sign=-1
  x=xu
 z=zu
 if(idd.eq.-9) write(*,*) 'sub 1'
 if(x.lt.x0.or.x.gt.z_e) goto 1
 do i=2,nx
    xn1=x0+(i-2)*dx
    xn2=x0+(i-1)*dx
    if(x.ge.xn1.and.x.lt.xn2) then
       vn=vel(i-1,1)+(vel(i,1)-vel(i-1,1))/dx*(x-xn1)
         end if
 end do
 p=1./vn
 px=p*sin(qq)
 pz=p*cos(qq)
 v=vn 
 dt=0.0001
      ix1=int((x-x0)/dx)+1
      izz=int(z/dz)+1
  x1=x0+(ix1-1)*dx
 x2=x1+dx
 z1=z0+(izz-1)*dz
 z2=z1+dz
      tx=(x-x1)/dx
 tz=(z-z1)/dz
 tv1=vel(ix1+1,izz)-vel(ix1,izz)
 tv2=vel(ix1+1,izz+1)-vel(ix1,izz+1)
      vv1=vel(ix1,izz)+tx*tv1
 vv2=vel(ix1,izz+1)+tx*tv2
 v=vv1+tz*(vv2-vv1)
 dvx1=tv1/dx
 dvx2=tv2/dx
 dvx=dvx1+tz*(dvx2-dvx1)
      dvz1=(vel(ix1,izz+1)-vel(ix1,izz))/dz
 dvz2=(vel(ix1+1,izz+1)-vel(ix1+1,izz))/dz
 dvz=dvz1+tx*(dvz2-dvz1)
101 nray_point=1
      x_ray(nray_point)=x
 z_ray(nray_point)=z
 t(nray_point)=0
 pp(nray_point)=p
 ppx(nray_point)=px
 do k=1,1000000
 if(idd.eq.-9) then
 end if
    sign=-1
    x=x+dt*v*v*px
    z=z+dt*v*v*pz
    px=px-dt*p*p*v*dvx
    pz=pz-dt*p*p*v*dvz
    p1=px**2+pz**2
    p1=sqrt(p1)
         ix1=int((x-x0)/dx)+1
         izz=int((z-z0)/dz)+1
    x1=x0+(ix1-1)*dx
    x2=x1+dx
    z1=z0+(izz-1)*dz
    z2=z1+dz
    tx=(x-x1)/dx
    tz=(z-z1)/dz
    tv1=vel(ix1+1,izz)-vel(ix1,izz)
    tv2=vel(ix1+1,izz+1)-vel(ix1,izz+1)
         vv1=vel(ix1,izz)+tx*tv1
    vv2=vel(ix1,izz+1)+tx*tv2
    v=vv1+tz*(vv2-vv1)
    dvx1=tv1/dx
    dvx2=tv2/dx
    dvx=dvx1+tz*(dvx2-dvx1)
    dvz1=(vel(ix1,izz+1)-vel(ix1,izz))/dz
    dvz2=(vel(ix1+1,izz+1)-vel(ix1+1,izz))/dz
    dvz=dvz1+tx*(dvz2-dvz1)
    p=1/v
          a=1.
     px=px*a
     pz=pz*a
     nray_point=nray_point+1
     x_ray(nray_point)=x
     z_ray(nray_point)=z
     t(nray_point)=k*dt
  pp(nray_point)=p
     ppx(nray_point)=px
     ppz(nray_point)=pz
     sign=1
  if(x.le.x0.or.x.ge.x0+(nx-1)*dx) goto 1
     if(z.le.z0.or.z.ge.z0+(nz-1)*dz) goto 1
 end do
1     continue
      xxx=x_ray(nray_point)
      yyy=z_ray(nray_point)
 ttt=t(nray_point)
 ppp=pp(nray_point)
 pppx=ppx(nray_point)
 pppz=ppz(nray_point)
      xxx2=x_ray(nray_point-1)
      yyy2=z_ray(nray_point-1)
 ttt2=t(nray_point-1)
 ppp2=pp(nray_point-1)
 pppx2=ppx(nray_point-1)
 pppz2=ppz(nray_point-1)
      continue
      if(idd.eq.-9) write(*,*) 'sub 4'
      return
      end