subroutine mag_bz(x0,y0,z0,bzave)

      implicit none

      integer i, j, k, n, repno
      parameter(n=1500)
      double precision x0, y0, z0
      double precision x(n), y(n), z(n), bz(n), tmp
      double precision dr(8),rat(8),bzp(8), rattot, bzmrat, bzave
      integer xrep, yrep, zrep
      parameter(xrep=3)
      parameter(yrep=1)
      parameter(zrep=60)
c
      open(11,file='z.dat',status='old')
      i=1
 20   read(11,*,end=30) x(i), y(i), z(i), bz(i)
      i=i+1
      goto 20
 30   tmp = 1
c
      bzave=0
      repno=0
      do k=1,i-1
         if(  (x0.ge.x(k).and.x0.le.x(k+xrep) .and.
     &        abs(x(k)-x(k+xrep)).lt.1. )     .and.
     &        (y0.ge.y(k).and.y0.le.y(k+yrep) .and.
     &        abs(y(k)-y(k+yrep)).lt.3. )     .and.
     &        (z0.ge.z(k).and.z0.le.z(k+zrep) .and.
     &        abs(z(k)-z(k+zrep)).lt.2. )   ) then
            repno=k
         endif
      enddo
c
      do k=1,i-1
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) bzp(1)=bz(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bzp(2)=bz(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bzp(3)=bz(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bzp(4)=bz(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) bzp(5)=bz(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bzp(6)=bz(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bzp(7)=bz(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bzp(8)=bz(k)
      enddo
c
      dr(1)=sqrt((x0-x(repno))**2      + (y0-y(repno))**2      +
     &     (z0-z(repno))**2)
      dr(2)=sqrt((x0-x(repno))**2      + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno))**2)
      dr(3)=sqrt((x0-x(repno))**2      + (y0-y(repno))**2      +
     &     (z0-z(repno+zrep))**2)
      dr(4)=sqrt((x0-x(repno))**2      + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno+zrep))**2)
      dr(5)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno))**2      +
     &     (z0-z(repno))**2)
      dr(6)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno))**2)
      dr(7)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno))**2      +
     &     (z0-z(repno+zrep))**2)
      dr(8)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno+zrep))**2)
c
      do j=1,8
         if(dr(j).eq.0) dr(j)=1e-10
         rat(j)=1./dr(j)
      enddo
c
      rattot=0.
      bzmrat=0.
      do j=1,8
         rattot=rattot+rat(j)
         bzmrat=bzmrat+rat(j)*bzp(j)
      enddo
      bzave=bzmrat/rattot
      
      close(11)

 100  return
      end