subroutine orbxy(ierr)
  !
  !
  !       bob mccann, 10-apr-86
  !
  !
  !       --------
  !       comments
  !       --------
  !
  !
  !	this routine advances the guiding center position in pseudo (x,y)
  !    coordinates for the sake of numerical accuracy.  the v_parallel
  !    equation is included for sake of having all differential advances
  !    happen consistently (note that bidt is halved here).
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  use nbi_com
  !
  implicit none
  integer, intent(out) :: ierr  ! status code, 0=ok
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !	bidt	    !both input and output
  !	bicos
  !	dxizc
  !	ledge
  !	bisin
  !	xbmbnd
  !
  !	bith2t(ion)
  !	bivp2t(ion)
  !	bixi2t(ion)
  !	bithdt(ion)
  !	bivpdt(ion)
  !	bixidt(ion)
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !	bithst
  !	bixi	    !both input and output
  !
  !	bidvpa(ion)
  !
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !       ---------------
  !       data statements
  !       ---------------
  !
  integer, parameter :: ion=1
  !
  real(fp), parameter :: zminxi = 1.0e-5_fp  !minimum xi allowed
  !
  real(fp) :: zximin,zximax,zx,zy,zdummy,zdtrob,zfac,zmax,zmin
  !
  !       -------------------
  !       statement functions
  !       -------------------
  !
  !
  !***
  !
  !
  !=========================================
  !
  ierr = 0
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  !
  !
  zximin = zminxi				!min xi allowed
  zximax = xbmbnd + 0.5_fp*dxizc(ledge)	!max xi allowed
  !
  !
  !       -----------------------------------------------------------
  !       1.0     convert the calculation to pseudo (x,y) coordinates
  !       -----------------------------------------------------------
  !
  !	note that this section computes the integrals in (x,y)
  !	coordinates.  this is a temporary kluge to enhance the
  !	numerical accuracy near xi=0.0.  this version still has
  !	a numerical singularity there.
  !
  !		zx = bixi(ion)*cos(zth)
  !		zy = bixi(ion)*sin(zth)
  !	thus,
  !		d(x)/dt = cos(ts)*bixidt(ion) - sin(ts)*bixi(ion)*bithdt(ion)
  !		d(y)/dt = sin(ts)*bixidt(ion) + cos(ts)*bixi(ion)*bithdt(ion)
  !
  !	and similarly for the second derivative
  !
  !
  !   compute the initial (x,y) position
  !
  zx = bixi(ion)*bicos(ion)
  zy = bixi(ion)*bisin(ion)
  !
  !   compute the first derivatives
  !
  bidxdt(ion) = bicos(ion)*bixidt(ion) &
       - bisin(ion)*bixi(ion)*bithdt(ion)
  bidydt(ion) = bisin(ion)*bixidt(ion) &
       + bicos(ion)*bixi(ion)*bithdt(ion)
  !
  !   compute the second derivatives
  !
  bi2x2t(ion) = bicos(ion)*bixi2t(ion) &
       - bisin(ion)*bixi(ion)*bith2t(ion) &
       - 2.0_fp*bisin(ion)*bixidt(ion)*bithdt(ion) &
       - bicos(ion)*bixi(ion)*bithdt(ion)**2
  !
  bi2y2t(ion) = bisin(ion)*bixi2t(ion) &
       + bicos(ion)*bixi(ion)*bith2t(ion) &
       + 2.0_fp*bicos(ion)*bixidt(ion)*bithdt(ion) &
       - bisin(ion)*bixi(ion)*bithdt(ion)**2
  !
  !  compute the timestep
  !
  bidt(ion) = zdtrob(ierr)
  if(ierr.ne.0) return
  !
  bidt(ion) = fdtnxy*bidt(ion)      !half the usual time step
  !						!for (x,y) integration
  !  shorten bidt further to avoid dominance of 2nd order terms, which seem
  !  to be subject to fluctuations near the axis
  !
  bidt(ion)=min(bidt(ion), &
       max(abs(bidxdt(ion)),abs(bidydt(ion)))/ &
       max(epslon,abs(bi2x2t(ion)),abs(bi2y2t(ion))) )
  !
  !       --------------------------------------
  !       2.0     compute the new (x,y) position
  !       --------------------------------------
  !
20000 continue
  !
  bidx(ion) = bidt(ion)*(bidxdt(ion) + 0.5_fp*bidt(ion)*bi2x2t(ion))
  bidy(ion) = bidt(ion)*(bidydt(ion) + 0.5_fp*bidt(ion)*bi2y2t(ion))
  bidvpa(ion)= bidt(ion)*(bivpdt(ion) +0.5_fp*bidt(ion)*bivp2t(ion))
  !
  zx = zx + bidx(ion)
  zy = zy + bidy(ion)
  !
  !
  !       -------------------------------------------
  !       3.0     convert back to (xi,th) coordinates
  !       -------------------------------------------
  !
  zmax = max(abs(zx),abs(zy))
  zmin = min(abs(zx),abs(zy))
  !
  bixi(ion) = zmax*sqrt(1.0_fp+ (zmin/zmax)**2)	!avoid underflow
  !
  !    now limit bixi to useful values
  !
  bixi(ion) = max(zximin,min(zximax,bixi(ion)))
  !
  !  automatically in range (-pi,pi)
  !
  bithst(ion) = atan2(zy,zx)

  return
end subroutine orbxy
