c
c *** Find the closest nodes to the source ***
c
      parameter(mnode=300000)
      dimension x(mnode),y(mnode),z(mnode)

c *** Source location ***
c      xs=3.0
c      ys=3.0
c      zs=-5.0
      write(*,*)'Enter Location of the source (km)'
      read(*,*)xs,ys,zs

c *** Open node data ***
      open(10,file='try.2.node',status='old')
        read(10,*)nnode,ntet,ninfo,nflag
        write(*,*)'Number of nodes=',nnode
        if(nnode.ge.mnode)then
          write(*,*)'nnode>mnode !!!'
          stop
        end if
        do i=1,nnode
          read(10,*)n,x(i),y(i),z(i),nf
        end do
      close(10)

c *** Find the domain lengths and the minimum node length ***
      xmin=99999.0
      xmax=0.0      
      ymin=99999.0
      ymax=0.0      
      zmin=99999.0
      zmax=0.0      
      dxmin=99999.0
      dymin=99999.0
      dzmin=99999.0
      do i=1,nnode
        if(x(i).lt.xmin)xmin=x(i)
        if(x(i).gt.xmax)xmax=x(i)
        if(y(i).lt.ymin)ymin=y(i)
        if(y(i).gt.ymax)ymax=y(i)
        if(z(i).lt.zmin)zmin=z(i)
        if(z(i).gt.zmax)zmax=z(i)
        if((abs(x(i)).lt.dxmin).and.(x(i).ne.0.0))dxmin=abs(x(i))
        if((abs(y(i)).lt.dymin).and.(y(i).ne.0.0))dymin=abs(y(i))
        if((abs(z(i)).lt.dzmin).and.(z(i).ne.0.0))dzmin=abs(z(i))
      end do
      write(*,*)'xmin, xmax=',xmin,xmax
      write(*,*)'ymin, ymax=',ymin,ymax
      write(*,*)'zmin, zmax=',zmin,zmax
      write(*,*)'dxmin, dymin, dzmin=',dxmin,dymin,dzmin

c *** Find the closest node number to the source point ***
      dlmin=xmax*ymax
      nmin=0
      do i=1,nnode
        dl=((x(i)-xs)*(x(i)-xs)+(y(i)-ys)*(y(i)-ys)
     *     +(z(i)-zs)*(z(i)-zs))
        if(dl.lt.dlmin)then
          dlmin=dl
          nmin=i
        end if
      end do
      write(*,*)
      write(*,1000)nmin,x(nmin),y(nmin),z(nmin)
 1000 format(' N-source ,x,y,z=',i9,3f9.3)
      write(*,*)

c *** Find minus- and plus-nodes next to the source node along the x-axis ***
      dx=dxmin/30.0
      nminx=nmin
      dlmin=xmax*ymax
      xmins=x(nmin)
 10   xmins=xmins-dx
      do i=1,nnode
        dl=((x(i)-xmins)*(x(i)-xmins)+(y(i)-ys)*(y(i)-ys)
     *     +(z(i)-zs)*(z(i)-zs))
        if(dl.lt.dlmin)then
          dlmin=dl
          nminx=i
        end if
      end do
      if(nminx.eq.nmin)goto 10
      hx=x(nmin)-x(nminx)
      write(*,1100)nminx,x(nminx),y(nminx),z(nminx),hx
 1100 format(' N-x-minus,x,y,z=',i9,3f9.3,',  HX=',f9.3)

      nplsx=nmin
      dlmin=xmax*ymax
      xplus=x(nmin)
 20   xplus=xplus+dx
      do i=1,nnode
        dl=((x(i)-xplus)*(x(i)-xplus)+(y(i)-ys)*(y(i)-ys)
     *     +(z(i)-zs)*(z(i)-zs))
        if(dl.lt.dlmin)then
          dlmin=dl
          nplsx=i
        end if
      end do
      if(nplsx.eq.nmin)goto 20
      hx=x(nplsx)-x(nmin)
      write(*,1110)nplsx,x(nplsx),y(nplsx),z(nplsx),hx
 1110 format(' N-x-plus ,x,y,z=',i9,3f9.3,',  HX=',f9.3)

c *** Find minus- and plus-nodes next to the source node along the y-axis ***
      dy=dymin/30.0
      nminy=nmin
      dlmin=xmax*ymax
      ymins=y(nmin)
 30   ymins=ymins-dy
      do i=1,nnode
        dl=((x(i)-xs)*(x(i)-xs)+(y(i)-ymins)*(y(i)-ymins)
     *     +(z(i)-zs)*(z(i)-zs))
        if(dl.lt.dlmin)then
          dlmin=dl
          nminy=i
        end if
      end do
      if(nminy.eq.nmin)goto 30
      hy=y(nmin)-y(nminy)
      write(*,1200)nminy,x(nminy),y(nminy),z(nminy),hy
 1200 format(' N-y-minus,x,y,z=',i9,3f9.3,',  HY=',f9.3)

      nplsy=nmin
      dlmin=xmax*ymax
      yplus=y(nmin)
 40   yplus=yplus+dy
      do i=1,nnode
        dl=((x(i)-xs)*(x(i)-xs)+(y(i)-yplus)*(y(i)-yplus)
     *     +(z(i)-zs)*(z(i)-zs))
        if(dl.lt.dlmin)then
          dlmin=dl
          nplsy=i
        end if
      end do
      if(nplsy.eq.nmin)goto 40
      hy=y(nplsy)-y(nmin)
      write(*,1210)nplsy,x(nplsy),y(nplsy),z(nplsy),hy
 1210 format(' N-y-plus ,x,y,z=',i9,3f9.3,',  HY=',f9.3)

c *** Find minus- and plus-nodes next to the source node along the z-axis ***
      dz=dzmin/30.0
      nminz=nmin
      dlmin=xmax*ymax
      zmins=z(nmin)
 50   zmins=zmins-dz
      do i=1,nnode
        dl=((x(i)-xs)*(x(i)-xs)+(y(i)-ys)*(y(i)-ys)
     *     +(z(i)-zmins)*(z(i)-zmins))
        if(dl.lt.dlmin)then
          dlmin=dl
          nminz=i
        end if
      end do
      if(nminz.eq.nmin)goto 50
      hz=z(nmin)-z(nminz)
      write(*,1300)nminz,x(nminz),y(nminz),z(nminz),hz
 1300 format(' N-z-minus,x,y,z=',i9,3f9.3,',  HZ=',f9.3)

      nplsz=nmin
      dlmin=xmax*ymax
      zplus=z(nmin)
 60   zplus=zplus+dz
      do i=1,nnode
        dl=((x(i)-xs)*(x(i)-xs)+(y(i)-ys)*(y(i)-ys)
     *     +(z(i)-zplus)*(z(i)-zplus))
        if(dl.lt.dlmin)then
          dlmin=dl
          nplsz=i
        end if
      end do
      if(nplsz.eq.nmin)goto 60
      hz=z(nplsz)-z(nmin)
      write(*,1310)nplsz,x(nplsz),y(nplsz),z(nplsz),hz
 1310 format(' N-z-plus ,x,y,z=',i9,3f9.3,',  HZ=',f9.3)

      stop
      end




