subroutine read_field(field_file, map,xmin)

   use magfield
   use magfield_interface

   implicit none
   type(magfield_struct), allocatable :: position(:), map(:,:,:)
   type(magfield_struct) x0, xmax, xmin, dxyz(-1:1)
   integer iu, istat/0/, nlines, ind(3),i,j,k,l, ntot(3), ix,jx
   integer k1, k2
   integer delta
   character(*)  field_file
   character*140 string
   real(rp)  b1(3), b2(3)

   iu = lunget()
   print '(a19,a)',' Open field file =', field_file
   open(unit=iu, file = field_file, status = 'old', action='read')
   nlines = 0
    xmin%x=1000.
   xmax%x=-1000.
   do while(.true.)
      string(1:140) = '                                                                                                                                                                       '
      read(iu, '(a140)', IOSTAT = istat ) string
      if(istat < 0) exit
     if(index(string,'#') /= 0)cycle
     if(string(1:10) == '          ')cycle
     if(index(field_file,'inj')/=0)read(string, *)x0%x(1:3), x0%B(1:3), x0%H
     if(index(field_file,'inf')/=0)read(string, *)x0%x(1:3), x0%B(1:3)
     if(index(field_file,'712')/=0)read(string, *)x0%x(1:3), x0%B(1:3)
     xmin%x(1:3) = min(xmin%x(1:3),x0%x(1:3))
     xmax%x(1:3) = max(xmax%x(1:3),x0%x(1:3))

     nlines = nlines + 1
   end do

!   print *, ' nlines = ', nlines
   print *,' xmax ',xmax%x, 'xmin', xmin%x
   if(index(field_file,'inj')/=0 .or. index(field_file,'inf') /=0)ntot(1:3) = (xmax%x(1:3)-xmin%x(1:3))*2 + 1
   if(index(field_file,'712')/=0 )ntot(1:3) = (xmax%x(1:3)-xmin%x(1:3))*4 + 1
!   print *, ntot
   allocate(map(ntot(1),ntot(2),ntot(3)))
   rewind(unit=iu)

   do while(.true.)
      string(1:140) = '                                                                                                                                                                       '
     read(iu,'(a140)',IOSTAT = istat) string
     if(istat < 0)exit
     if(index(string,'#') /= 0)cycle
     if(string(1:10) == '          ')cycle
     if(index(field_file,'inj')/=0)read(string, *)x0%x(1:3), x0%B(1:3), x0%H
     if(index(field_file,'inj')==0)read(string, *)x0%x(1:3), x0%B(1:3)

    ind(1:3) = (x0%x(1:3) - xmin%x(1:3))*2 + 1
    map(ind(1),ind(2),ind(3))%B(1:3) = x0%B(1:3)
!    if(index(field_file,'inj')/=0) &         !temporary to get plot of field with fringe or inflector only
!     map(ind(1),ind(2),ind(3))%B(1:3) = 0.   !20140727 to zero out either inflector field or fringe

    map(ind(1),ind(2),ind(3))%x(1:3) = x0%x(1:3)
  end do
! correct the discontinuity at z=259.5. Set the field at this point to be average of 259 and 260
!try skipping averaging now that there are two elements
if(1<0)then
  k1 = 259.5*2+1
  k2 = 430.5*2+1
  do i = 1,ntot(1)
   do j = 1,ntot(2)
!    print *,'map(i,j,k)%x(3)',map(i,j,k1)%x(3), map(i,j,k2)%x(3)
    map(i,j,k1)%B(1:3) = (map(i,j,k1-1)%B(1:3) + map(i,j,k1+1)%B(1:3))*0.5
    map(i,j,k2)%B(1:3) = (map(i,j,k2-1)%B(1:3) + map(i,j,k2+1)%B(1:3))*0.5
   delta = 10
   do k = -delta,delta
    map(i,j,k1+k)%B(1) = (map(i,j,k1-delta)%B(1) + map(i,j,k1+delta)%B(1))*0.5
    map(i,j,k2+k)%B(1) = (map(i,j,k2-delta)%B(1) + map(i,j,k2+delta)%B(1))*0.5
   end do
   end do
 end do
endif
return
end