!----------------------------------------------------------------------
!  nbimov
!
subroutine nbimov(zdr,zdy,ierr)
  !
  !
  !  move beam ion by small increment
  !   use difference in fast inverse map to get change in mag. coords
  !
  !  change in ion (xi,theta) resulting...
  !  if far from axis use jacobian
  !
  use xstraln_calls

  use nbi_com

  !============
  implicit none
  integer :: ier,ierp,ierl,ierr
  !============
  real(fp) :: zdy,zdr,zrp,zyp,zphiarg,zxi,zth,zxip,zthp
  !============
  integer, parameter :: ivec = 2  ! 2 points from xplasma
  integer :: iregion(ivec)        ! output region code
  integer :: ixer                 ! error flag from xplasma

  real(fp) :: zzr(ivec), zzy(ivec), zzphi(ivec) ! convert from these coordinates
  real(fp) :: zzrho(ivec), zzchi(ivec)          ! convert to   these coordinates
  real(fp) :: ztol                              ! tolerance for conversion
  real(fp), parameter :: zconv = .01_fp  ! cm to meter conversion
  !
  !----------------------------------------------
  ierr = 0

  zrp=rmjion
  zyp=yion
  !
  rmjion = rmjion + zdr
  yion = yion + zdy
  !
  !  use the fast inverse map
  !  to reduce error from this approximation use the difference and add
  !  to the orig. magnetic coordinates
  !
  zphiarg=0.0_fp
  !test      call nbimap(rmjion,yion,zxi,zth,ier)
  call xs_fxget(0.01_fp,zphiarg,rmjion,yion,zxi,zth,ier)
  if(ier.gt.1) then
     rmjion = zrp
     yion = zyp
     ierr = 999 !repeat
     return
  end if
  !test      call nbimap(zrp,zyp,zxip,zthp,ierp)
  call xs_fxget(0.01_fp,zphiarg,zrp,zyp,zxip,zthp,ierp)
  ierl=max(ier,ierp)
  !
  if(ierl.gt.1) then
     rmjion = zrp
     yion = zyp
     ierr = 999
     return
     !        write(nonlin,9001) xion,th

9001 format( &
          ' %nbimov - inverse map error, mag coords not changed:'/ &
          '   xi=',1pe11.4,' th=',1pe11.4)
     !
  else
     th=th+(zth-zthp)
     !  keep the range in order
     th = th + twopi + xpi
     th = th - twopi*(int(th/twopi))
     th = th - xpi
     !
     if(abs(th).gt.(xpi+0.0001_fp)) then
        call nbi_message( &
             ' ?nbimov: th range error even after correction.')
        ierr=1
     end if
     !
     xion=xion+(zxi-zxip)
  end if
  !
  !  if xion .lt. 0.0 ...
  if(xion.lt.0.0_fp) then
     xion=-xion
     if(th.le.0.0_fp) then
        th=th+xpi
     else
        th=th-xpi
     end if
  end if
  !
  return
end subroutine nbimov
