      subroutine mag_b(x0,y0,z0,bave,zrep,flag)

      implicit none

      include 'common.inc'

c
      if(flag.eq.1) open(11,file='x.dat',status='old')
      if(flag.eq.2) open(11,file='y.dat',status='old')
      if(flag.eq.3) open(11,file='z.dat',status='old')
      i=1
 20   read(11,*,end=30) x(i), y(i), z(i), b(i)
      i=i+1
      goto 20
 30   tmp = 1
c
      bave=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)) bp(1)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bp(2)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bp(3)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bp(4)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) bp(5)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bp(6)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bp(7)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bp(8)=b(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)**2
      enddo
c
      rattot=0.
      bmrat=0.
      do j=1,8
         rattot=rattot+rat(j)
         bmrat=bmrat+rat(j)*bp(j)
      enddo
      bave=bmrat/rattot

      close(11)

 100  return
      end
