recursive subroutine em_field_custom_inj(ele, param, s_rel,t_rel, orb, local_ref,field,calcd, err_flag, &
                                  calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)
use bmad
use em_field_mod
use em_field_interface, dummy =>em_field_custom_inj
use magfield
use parameters_bmad
   implicit none
   type(ele_struct) :: ele
   type(lat_param_struct) param
   type(coord_struct) orb
   type(em_field_struct) field
   type (ele_struct), pointer :: lord
   type (ele_struct), pointer :: slave
   type (ele_pointer_struct), allocatable, optional :: used_eles(:)
   type (grid_field_struct), pointer ::  g_field

   real(rp), intent(in) :: s_rel, t_rel
   real(rp) xx(3), BB(3),zz, b1(3,3),b2(3,3)
!   real(rp) xoffset/0.089/, yoffset/0.0/, zoffset
   real(rp) xoffset/0.0/, yoffset/0.0/, zoffset
   real(rp) zoff/-9999999./
   real(rp), save :: zref_hole, zref_free, zref
   real(rp), save :: length
   real(rp), optional :: rf_time
   real(rp) :: s_body,  x, y, z, x_save, s_norm
   real(rp) cos_ang, sin_ang

   logical  out_of_range,local_ref
   logical first/.true./, first_out/.true./
   logical grid/.false./
   logical local_ref_frame/.false./, grid_out_bounds
   logical, optional :: calcd, err_flag, calc_potential, grid_allow_s_out_of_bounds, use_overlap
   logical err

   character(32) :: r_name = 'em_field_custom'
   character*120 file_inj, file_inf
   character*16 last/' '/
   integer n/0/
   integer i
   integer ig0, ig1, i0


 ! then check to see if this element has a grid map
  if ((associated(ele%grid_field)))then
     grid = .true.
    elseif(ele%field_calc == refer_to_lords$)then
         do i=1,ele%n_lord
          lord =>pointer_to_lord(ele,i)
          if(associated(lord%grid_field))then
           grid = .true.
           exit
          endif
         end do
   endif

If(grid)then
    ele%field_calc = fieldmap$ 
    grid_out_bounds = .true.
    call em_field_calc(ele, param, s_rel, orb, local_ref, field,  calcd, err_flag, &
                              calc_potential, use_overlap, .true., rf_time, used_eles)

! check to see if out of bounds

     g_field => ele%grid_field(1)
     call to_fieldmap_coords (ele, orb, s_rel, g_field%ele_anchor_pt, g_field%r0, g_field%curved_ref_frame, x, y, z, &
                                           cos_ang, sin_ang, err)

    ig0 = lbound(g_field%ptr%pt, 3)
    ig1 = ubound(g_field%ptr%pt, 3)

     s_norm = s_rel / g_field%dr(3)  ! Note that to_fieldmap_coords has already been called.
     i0 = floor(s_norm)          ! index of lower 1 data point
!     rel_s0 = s_norm - i0        ! Relative distance from lower x1 grid point

!      if (i0 < ig0 .or. i0 >= ig1) then
!        write(109, '(a,a,a,3es12.4,a,3es12.4)')' out of bounds:',ele%name,' x,y,z = ',x,y,z, ' field%B =',field%B
!      endif



!    grid_allow_s_out_of_bounds = .false.
    ele%field_calc = custom$
    if(all(field%B(:) == 0)) then
!       field%B(2) =  -1.451
!     print '(a,3es12.4)',' em_field_custom_inj: field%B = ',field%B
    endif
   return
endif

 if(first)then
   do i=1,2
   print '(i10,a,a)',i,' field_file(i)%name = ',field_file(i)%name
   print '(i10,a,i10)',i,' field_file(i)%type = ',field_file(i)%type
   print '(i10,a,es12.4)',i,' field_file(i)%grid_spacing = ',field_file(i)%grid_spacing
   end do
   first=.false.
   zref=0.
 endif  ! first time 

 ! meters to cm
   xx(1) = (orb%vec(1)+xoffset)*100. 
   xx(2) = (orb%vec(3)+yoffset)*100.

  if(s_rel < 1.e-10 .and. zoff < -99)then
     zoff = ele%s-ele%value(l$)
     print *, ' zoff = ', zoff
   endif
   zoffset = ele%s - ele%value(l$) - zoff 
   xx(3) = (zoffset + s_rel)*100.                
   out_of_range = .false.

!   if(ele%name /= 'SKIP' .and. (s_rel < zref .or. s_rel > zref+length))then
   if(ele%name /= 'SKIP')then

     call get_g2_fields(xx, BB,b1,b2,out_of_range)
     field%B(1:3) = -BB(1:3) * 1.e-4  ! gauss to tesla
   else
    field%B=0.

  endif

   if(out_of_range .and. ele%alias == 'trajectory')then
      if(first_out)then
         print *,' write out of range to unit 19'
         write(19, '(/,a)')'Out of range'
         first_out = .false.
       write(19, '(2(a,3(es12.4)),1x,a,3i10)')' xx ', xx, ' field%B', field%B, ele%name,  int((xx(1:3))*2+1)
       endif
   endif
   return
   end subroutine