c  phase.for - calculates the local phase from dx, dy, and dz grids
c
      character id*56, pgm*8, ifile*50
      dimension x(2048),y(2048),z(2048)
c
      dval=1.e38
      ddval=1.e30
      print *,'enter x-derivative file name:'
      read(5,100) ifile
  100 format(a)
      open(12,file=ifile,form='unformatted',status='old')
      print *,'enter y-derivative file name:'
      read(5,100) ifile
      open(13,file=ifile,form='unformatted',status='old')
      print *,'enter z-derivative file name:'
      read(5,100) ifile
      open(14,file=ifile,form='unformatted',status='old')
      read(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
      read(13) id
      read(14) id
      print *,'enter output file name:'
      read(5,100) ifile
      open(15,file=ifile,form='unformatted',status='unknown')
      pgm='phase'
      id='local phase'
      write(15) id,pgm,nc,nr,nz,xo,dx,yo,dy
      do 20 j=1,nr
      read(12) dlt,(x(i),i=1,nc)
      read(13) dlt,(y(i),i=1,nc)
      read(14) dlt,(z(i),i=1,nc)
      do 10 i=1,nc
      if(x(i).gt.ddval.or.y(i).gt.ddval.or.z(i).gt.ddval) then
              z(i)=dval
              go to 10
      endif
      den=sqrt(y(i)**2 + x(i)**2)
      z(i)=atan2(z(i),den)
10    continue
      write(15) dlt,(z(i),i=1,nc)
20    continue
      stop
      end
