subroutine orbg(ierr)
  use iso_c_binding, only: fp => c_double
  !
  !       bob mccann, 13-sep-84
  !	bob mccann, dec-jan-86	!p-phi, phi(xi)
  !	bob mccann, 23-apr-86	!zdtrob, zdvopt, orbxy
  !
  !       dmc -- hack to shorten timestep near + outside plasma bdy
  !         20 jun 1990
  !
  !       mhr  july 94 -- ripple model coding added
  !
  !       jfs -- replace ichdvp (character) with nchdvp (integer ::) for choice of
  !               method to compute d_vparallel    2 jun 94
  !
  !       jfs -- analytic formulas for dpsi/dxi and d2psi/dxi2  (16 may 95)
  !       rga nov00 -- pulled out spline interpolations to gintrp
  !
  !       --------
  !       comments
  !       --------
  !
  !	this advances the orbit of beam ions one step in (xi,ts)
  !   co-ordinates (ts is theta-star).  these orbit equations agree
  !   with those derived by white and chance, pppl-2089, feb. 1984.
  !   here it is assumed that |b| is a function of xi and theta-star.
  !   the spline coefficients are provided by subroutine magcf.
  !
  !	all 1-d quantities are cubic spline interpolated in the
  !   minor radial direction.  all 2-d quantities are fourier inter-
  !   polated in the angular direction and the fourier coefficients
  !   are cubic spline interpolated in the minor radial direction.
  !   the fourier transforms are done by the imsl routine fftsc.  the
  !   interpolations are performed in the subroutines gintrp,bintrp,
  !   and jintrp.
  !
  !	the (-nsnccw) hack has been removed because only the v_parallel
  !   drift term in the theta and v_parallel equations depends upon the
  !   handedness of the poloidal magnetic field.  this can be accomplished
  !   by assigning a sign to q while using v_parallel .gt. 0 if it is in
  !   the "co" direction (v_phi*j_phi.gt.0). note that psi_t_limiter
  !   is defined to be always positive, and b_toroidal is always positive.
  !   so the (psi_p, theta, phi) coordinate system is either right or
  !   left handed as needed, even though the (xi, theta, phi) system is
  !   always right handed.
  !
  !
  !     -------------------
  !     global declarations
  !     -------------------
  !
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: ion
  !============
  real(fp) :: zgb,zcb,ze,zg1,zpsip1
  real(fp) :: zphi1,zdxdps,zj,zmr1,ztotal1,zeb,zvp,zvcb,zveb,zvvp
  real(fp) :: zvpx,zbtx,zdgdt,z2xpst,z2phxt,z2gxit,zdbdt,zdjdt
  real(fp) :: z2btst,z2bxit,zdcbdt,zdebdt,zdgbdt,zdvpdt,zvcbdt
  real(fp) :: zvebdt,zvvpdt,zdg,zdpsip,zd1ovb,zdbdx,zdbdy,z2b2x
  real(fp) :: z2b2y,z2bxy,zb,zg2,zpsip2,zphi2,zdphi,zdb,zzdvp,zdvopt
  real(fp) :: zdummy,zvlim,zmr2,ztotal2
  !
  integer, intent(out) :: ierr  ! status code, 0=ok
  !
  !     ---------------
  !     input variables
  !     ---------------
  !
  !     ctr denotes variables passed via transp common - dmc apr 85
  !
  !tr	integer ::	nbi		!# ions to orbit advance one step
  !tr	real	biv(mmaxion)	!total velocity of each ion(cm/sec)(conserved)
  !tr	real	bia(mmaxion)	!atomic number of beam ions
  !tr	real	biz(mmaxion)	!charge state of beam ions
  !tr	real	bismin		!min. arclength/step (xi units)
  !tr	real	bismax		!max. arelength/step (xi units)
  !tr	real	bieps		!max. fraction of radius of
  !					curvature/step
  !tr	real	bimag		!dtnew.le.dtold*bimag
  !tr	real	bidtnv	!norm vel.*time step for 1 kev passing ion
  !
  !	--------------
  !	from trans.com
  !	--------------
  !
  !	integer ::	nximax		!dim of # radial zones
  !	integer ::	momdim		!dim of # moments
  !	integer ::	nxi		!# radial zones for splines
  !	integer ::	nmomp1		!(1 + #moments used) .le. momdim
  !	real	tflux		!toroidal flux at limiter (tesla-m**2)
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !tr	real	bimodb(mmaxion)!|b| at each ion's guiding center(tesla)
  !
  !    the following are input as well as output
  !
  !tr	real	bidt(mmaxion)	!input=previous dt, output=current dt
  !tr	real	biksid(mmaxion)!v_par/v of each ion
  !tr	real	bithst(mmaxion)!theta-star of each ion
  !tr	real	bivpar(mmaxion)!parallel velocity of each ion (cm/sec)
  !tr	real	bivprp(mmaxion)!perpendicular velocity of each ion
  !tr	real	bixi(mmaxion)	!xi of each ion = sqrt(psi_t/psi_t_lim)
  !
  !tr	real	bi1ovb(mmaxion)	!1/b (1/gauss)
  !tr	real	bibdxi(mmaxion)	!d(1/b)/d(xi)
  !tr	real	bibdts(mmaxion)	!d(1/b)/d(theta-star)
  !tr	real	bib2xi(mmaxion)	!d2(1/b)/d2(xi)
  !tr	real	bib2ts(mmaxion)	!d2(1/b)/d2(theta-star)
  !tr	real	bibtsx(mmaxion)	!d2(1/b)/d(theta-star)d(xi)
  !
  !tr     real
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  parameter	(ion = 1) !vec	needed to force compiler opt.
  !
  real(fp) ::		zminxi		!minimum xi (avoid div. by zero)
  logical		ilfrst		!cancel second order terms(debug)
  logical		ildiff		!use differential differences
  !
  !
  !       ---------------
  !       data statements
  !       ---------------
  !
  !
  !
  data	zminxi /1.0e-3_fp/		!min xi: avoid div. by 0
  data	ilfrst /.false./	!true=> use only first order terms
  data	ildiff /.false./	!true=> use differential differences
  !             which method to use to compute d_vparallel (preset or namelist):
  !			nchdvp=1 => use conservation of angular momentum
  !			nchdvp=2 => use conservation of energy
  !			nchdvp=3 => use v_parallel equation
  !
  !=========================================
  !
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  ierr=0
  !
  bixi(ion) = max(bixi(ion),zminxi) !small jiggles don't count
  if(abs(bivpar(ion)).gt.abs(biv(ion)))  then
     write(nonlin,*) '  vpll > v in orbg: ',bivpar(ion),biv(ion)
     ierr=1
     return
  end if
  !
  ze = zel*biz(ion)
  bim(ion)  = zmp*bia(ion)
  bike(ion) = 0.5_fp*bim(ion)*biv(ion)**2
  bimu(ion) = 0.5_fp*bim(ion)*bi1ovb(ion)*bivprp(ion)**2 !0.5*m*(v_perp)**2/b
  !
  ! load up biq, big, bipsi, biphi and derivatives
  !
  call gintrp(1)

  ! save some data for later differencing
  zg1    = big(ion)
  zpsip1 = bipsi(ion)
  zphi1  = biphi(ion)

  zdxdps = 1.0_fp/bi1psi1x(ion)   ! d(xi)/d(psi)
  !
  !     ----------------------------
  !     1/b and its first and second
  !     derivatives are computed in
  !     bintrp
  !     ----------------------------
  !
  !     ----------------------------
  !     j(theta,xi) (jacobian) and
  !     derivates are computed in
  !     jintrp
  !     ----------------------------
  !
  !     -------------------------
  !     compute roscoe's jacobian
  !     j = (i+gq)/b**2
  !     -------------------------
  !
  !     **************************
  !     this affects d(j)/dt below
  !     **************************
  !
  !     zj = (zi + big(ion)*biq(ion))*bi1ovb(ion)**2	!roscoe's jacobian
  !
  zj = bijac(ion)           !arbitrary jacobian
  !
  !
  !     ----------------------------------------------------
  !     2.4	compute conserved canonical angular momentum
  !     ----------------------------------------------------
  !
  !     roscoe's conserved momentum:
  !
  bipphi(ion) = bim(ion)*bivpar(ion)*big(ion)*bi1ovb(ion) - &
       (ze/zc)*bipsi(ion)
  !
  !     debug stuf:
  !
  !        momentum:
  !
  zmr1 = bipphi(ion)
  !
  !        total energy:
  !
  ztotal1 = 0.5_fp*bim(1)*bivpar(1)**2 + bimu(1)/bi1ovb(1) + &
       ze*biphi(ion)
  !
  !
  !     ------------------------------------------
  !     3.0     compute the first time derivatives
  !     ------------------------------------------
  !
  !     this section needs the following constants:
  !
  !               zc  = speed of light (cgs)
  !		bimu(ion) = magnetic moment of ion larmour orbit
  !		ze  = charge on ion (cgs)
  !		bim(ion)  = mass of ion (cgs)
  !		zdxdps = d(xi)/d(psi_p)		!all but v_par drift
  !		bibdxi(ion) = d(1/b)/d(xi)	!for grad-b drift
  !		bibdts(ion) = d(1/b)/d(theta)	!for grad-b drift
  !               bi1g1x(ion) = d(g)/d(xi)	!new curl-b drift
  !		bi1phi1x(ion) = d(phi)/d(xi)	!for e x b drift
  !
  !     !curl-b drift cf
  !
  zcb = zdxdps*(zc*bi1ovb(ion)**2/(ze*zj))* &
       bim(ion)*bivpar(ion)**2
  !
  !     !exb drift cf
  !
  zeb = zdxdps*big(ion)*zc*bi1ovb(ion)**2/zj
  !
  !     !grad-b drift cf
  !
  zgb = zdxdps*(big(ion)*zc/(ze*zj))* &
       (bimu(ion) + bim(ion)*bi1ovb(ion)*bivpar(ion)**2)
  !
  !
  !     !v_parallel drift term
  !
  zvp = bivpar(ion)*bi1ovb(ion)/zj
  !
  !     drift velocity if vpll/v=0.25, for timestep control
  !
  !ap      bivpd(ion) = biv(ion)*0.25*bi1ovb(ion)/zj
  !
  !     !velocity curl-b drift cf
  !
  zvcb = zdxdps*bimu(ion)*zc*(bivpar(ion)/(ze*zj))*bi1g1x(ion)
  !
  !     !velocity exb drift cf
  !
  zveb = zdxdps*big(ion)*zc*bi1ovb(ion)*(bivpar(ion)/zj) &
       *bi1phi1x(ion)
  !
  !     !velocity v_parallel drift cf
  !
  zvvp = bimu(ion)/(bim(ion)*zj*bi1ovb(ion))
  !
  !     --------
  !     d(xi)/dt
  !     --------
  !
  bixidt(ion) = zgb*bibdts(ion) !d(xi)/dt
  !
  !     ----------------
  !     d(theta-star)/dt
  !     ----------------
  !
  bithdt(ion) =   - zcb*bi1g1x(ion)& !new curl-b drift &
       + zeb*bi1phi1x(ion)&           !(e_radial x b) drift &
       - zgb*bibdxi(ion)&	         !new grad-b drift form &
       + zvp                         !v_parallel motion
  !
  !
  zvpx=zvp*bixi(ion)
  zbtx=bithdt(ion)*bixi(ion)
  !
  !     -----------
  !     d(v-par)/dt
  !     -----------
  !
  bivpdt(ion) = ( - zvcb &    !new curl-b drift &
       - zveb &	 	!e x b drift &
       + zvvp)*bibdts(ion)	!v-parallel drift
  !
  !
  !     -----------------------------------------------------
  !     4.0     compute the moving derivatives for next order
  !     -----------------------------------------------------
  !
  !     (d/dt) = d(xi)/dt*(d/d(xi)) + d(th)/dt*(d/d(th))
  !
  !     where	(d/dt) = total time derivative
  !     and	(d/dt) = partial time derivative
  !                      = 0.0 for now
  !
  !
  !     -------
  !     d(g)/dt
  !     -------
  !
  zdgdt = bi1g1x(ion)*bixidt(ion)
  !
  !     -------
  !     d(q)/dt
  !     -------
  !
  !     zdqdt = zdqdxi*bixidt(ion)	!only used in old z2xpst
  !
  !     --------------------
  !     (d/dt)d(xi)/d(psi-p)
  !     --------------------
  !
  !     z2xpst = zdxdps*
  !     1	  (zdqdt/biq(ion) - bixidt(ion)/bixi(ion))
  !
  z2xpst = (-bi2psi2x(ion)/bi1psi1x(ion)**2)*bixidt(ion)
  !
  !     --------------------
  !     (d/dt)(d(phi)/d(xi))
  !     --------------------
  !
  z2phxt = bixidt(ion)*bi2phi2x(ion)
  !
  !     ------------------
  !     (d/dt)(d(g)/d(xi))
  !     ------------------
  !
  z2gxit = bixidt(ion)*bi2g2x(ion)
  !
  !     ---------
  !     d(1/b)/dt
  !     ---------
  !
  !     here i substitute to cancel the
  !     grad(b) drift terms analytically (rob's trick)
  !
  !
  zdbdt = (-zcb*bi1g1x(ion)+zeb*bi1phi1x(ion)+zvp)*bibdts(ion)
  !
  !     ----------------------------
  !     d(j)/dt (arbitrary jacobian)
  !     ----------------------------
  !
  zdjdt = bixidt(ion)*bijdxi(ion) + bithdt(ion)*bijdts(ion)
  !
  !     -----------------------
  !     (d/dt)(d(1/b)/d(theta))
  !     -----------------------
  !
  z2btst = bixidt(ion)*bibtsx(ion) + bithdt(ion)*bib2ts(ion)
  !
  !     --------------------
  !     (d/dt)(d(1/b)/d(xi))
  !     --------------------
  !
  z2bxit = bixidt(ion)*bib2xi(ion) + bithdt(ion)*bibtsx(ion)
  !
  !     -------------------
  !     (d/dt)(drift terms)
  !     -------------------
  !
  !     n.b. bivpar can goto zero!	!new form 9 dec.86
  !
  !     !(d/dt)(curl-b drift cf)
  !
  zdcbdt = zcb*(z2xpst/zdxdps + 2.0_fp*zdbdt/bi1ovb(ion) - zdjdt/zj) &
       + zdxdps*(zc*bi1ovb(ion)**2/(ze*zj))*bim(ion)* &
       2.0_fp*bivpar(ion)*bivpdt(ion)
  !
  !     !(d/dt)(exb drift cf)
  !
  zdebdt = zeb*(z2xpst/zdxdps + zdgdt/big(ion) &
       + 2.0_fp*zdbdt/bi1ovb(ion) &
       - zdjdt/zj)
  !
  !     !(d/dt)(grad-b drift cf)
  !
  zdgbdt = zgb*(z2xpst/zdxdps + zdgdt/big(ion) - zdjdt/zj) + &
       zdxdps*(big(ion)*zc/(ze*zj))*bim(ion)* &
       (zdbdt*bivpar(ion)**2 &
       + bi1ovb(ion)*2.0_fp*bivpar(ion)*bivpdt(ion))
  !
  !     !(d/dt)(v_parallel drift term)
  !
  zdvpdt = zvp*(zdbdt/bi1ovb(ion) - zdjdt/zj) + &
       bi1ovb(ion)*bivpdt(ion)/zj
  !
  !     !(d/dt)(velocity curl-b drift cf)
  !     n.b. bi1g1x(ion) can goto zero also
  !
  zvcbdt = zvcb*(z2xpst/zdxdps - zdjdt/zj) + &
       (zdxdps*bimu(ion)*zc/(ze*zj))* &
       (bivpdt(ion)*bi1g1x(ion) + bivpar(ion)*z2gxit)
  !
  !     !(d/dt)(velocity exb drift cf)
  !     n.b. bi1phi1x(ion) can goto zero also
  !
  zvebdt = zveb*(z2xpst/zdxdps + zdgdt/big(ion) &
       + zdbdt/bi1ovb(ion) - zdjdt/zj) + &
       (zdxdps*big(ion)*zc*bi1ovb(ion)/zj)* &
       (bivpdt(ion)*bi1phi1x(ion) + bivpar(ion)*z2phxt)
  !
  !     !(d/dt)(velocity v_parallel drift cf)
  !
  zvvpdt = zvvp*(-zdjdt/zj - zdbdt/bi1ovb(ion))
  !
  !
  !     -------------------------------------------
  !     5.0     compute the second time derivatives
  !     -------------------------------------------
  !
  !     ----------
  !     d2(xi)/d2t
  !     ----------
  !
  bixi2t(ion) = zdgbdt*bibdts(ion) + zgb*z2btst
  !
  !     ------------------
  !     d2(theta-star)/d2t
  !     ------------------
  !
  bith2t(ion) = &
       - (zdcbdt*bi1g1x(ion) + zcb*z2gxit) &  !curl-b drift &
       + (zdebdt*bi1phi1x(ion) + zeb*z2phxt)& !exb drift &
       - (zdgbdt*bibdxi(ion) + zgb*z2bxit)  & !grad-b drift &
       + zdvpdt                              !v_parallel motion
  !
  !     ------------
  !     d2(v-par)/dt
  !     ------------
  !
  bivp2t(ion) = &
       (- zvcbdt &		 !curl-b drift &
       - zvebdt &		 !exb drift &
       + zvvpdt)*bibdts(ion)&  !v-parallel drift &
       + ( - zvcb &		 !curl-b drift &
       - zveb &		 !exb drift &
       + zvvp)*z2btst	 !v-parallel drift
  !
  if(ilfrst) then		!use only the first order terms(debug only)
     bixi2t(1) = 0.0_fp
     bith2t(1) = 0.0_fp
     bivp2t(1) = 0.0_fp
  end if
  !
  !     ***note that (d/dt)(vpar*(p/pt)(jb)) is missing***
  !
  !
  !     ------------------------------------------------
  !     6.0     integrate the 2-nd order orbit equations
  !     ------------------------------------------------
  !
  !
  !     for now use rob's time step.  so far, everything i've tried ends
  !     up being much more complicated and gives only mixed results.
  !     rtm...26 june 1985.
  !
  !     dmc -- timestep calculation moved to integration routines
  !     orbxit and orbxy -- 20 jul 1994
  !
  if(bixi(ion).ge.dxitk)	then    !use xi-theta
     !
     call orbxit(ierr)              !compute d(xi) and d(theta_star)
     if(ierr.ne.0) return
     !
     if(ildiff)	then            !use differential differences
        !
        zdg = bidxi(ion)*(bi1g1x(ion)+0.5_fp*bidxi(ion)*bi2g2x(ion))
        zdpsip = bidxi(ion)*(bi1psi1x(ion) &
             + 0.5_fp*bidxi(ion)*bi2psi2x(ion))
        bidphi(ion) = bidxi(ion)* &
             (bi1phi1x(ion) + 0.5_fp*bidxi(ion)*bi2phi2x(ion))
        !
        zd1ovb = bidxi(ion)* &
             (bibdxi(ion) + 0.5_fp*bidxi(ion)*bib2xi(ion)) + &
             bidts(ion)* &
             (bibdts(ion) + 0.5_fp*bidts(ion)*bib2ts(ion)) + &
             bidxi(ion)*bidts(ion)*bibtsx(ion)
        !
     end if
     !
  else                      !use x-y (twiddle)
     !
     bicos(ion) = cos(bithst(ion))
     bisin(ion) = sin(bithst(ion))
     !
     call orbxy(ierr)               !compute d(x) and d(y)
     if(ierr.ne.0) return
     !
     if(ildiff) then        !use differential differences
        !
        bidxi(ion) = bicos(ion)*bidx(ion) + bisin(ion)*bidy(ion)
        !
        zdg = bidxi(ion)*(bi1g1x(ion)+0.5_fp*bidxi(ion)*bi2g2x(ion))
        zdpsip = bidxi(ion)*(bi1psi1x(ion) &
             + 0.5_fp*bidxi(ion)*bi2psi2x(ion))
        bidphi(ion) = bidxi(ion)* &
             (bi1phi1x(ion) + 0.5_fp*bidxi(ion)*bi2phi2x(ion))
        !
        zdbdx = bibdxi(ion)*bicos(ion) &
             - bibdts(ion)*bisin(ion)/bixi(ion)
        !
        zdbdy = bibdxi(ion)*bisin(ion) &
             + bibdts(ion)*bicos(ion)/bixi(ion)
        !
        z2b2x = &
             bibdxi(ion)*bisin(ion)**2/bixi(ion) &
             + bibdts(ion)*2.0_fp*bisin(ion)*bicos(ion)/bixi(ion)**2 &
             + bib2xi(ion)*bicos(ion)**2 &
             + bib2ts(ion)*bisin(ion)**2/bixi(ion)**2 &
             - bibtsx(ion)*2.0_fp*bisin(ion)*bicos(ion)/bixi(ion)
        !
        z2b2y = &
             bibdxi(ion)*bicos(ion)**2/bixi(ion) &
             - bibdts(ion)*bisin(ion)*bicos(ion)/bixi(ion)**2 &
             + bib2xi(ion)*bisin(ion)**2 &
             + bib2ts(ion)*bicos(ion)**2/bixi(ion)**2 &
             + bibtsx(ion)*2.0_fp*bisin(ion)*bicos(ion)/bixi(ion)
        !
        z2bxy = &
             - bibdxi(ion)*bisin(ion)*bicos(ion)/bixi(ion) &
             + bibdts(ion)*(bisin(ion)**2 - bicos(ion)**2) &
             /bixi(ion)**2 &
             + bib2xi(ion)*bisin(ion)*bicos(ion) &
             - bib2ts(ion)*bisin(ion)*bicos(ion)/bixi(ion)**2 &
             + bibtsx(ion)*(bicos(ion)**2 - bisin(ion)**2) &
             /bixi(ion)
        !
        zd1ovb = bidx(ion)*(zdbdx + 0.5_fp*bidx(ion)*z2b2x) + &
             bidy(ion)*(zdbdy + 0.5_fp*bidy(ion)*z2b2y) + &
             bidx(ion)*bidy(ion)*z2bxy
        !
     end if
     !
  end if
  !
  !
  !     ----------------------------------------------
  !     now recompute vpar from the conserved k.e.
  !     (energy changes in the compression subroutine)
  !     ----------------------------------------------
  !
  zb = bi1ovb(ion)          !store current value
  !
  !
  !     -----------------------------------
  !     recompute the necessary 1-d splines
  !     -----------------------------------
  !
  call gintrp(2)

  !
  !     -----------------------------------
  !     recompute the necessary 2-d splines
  !     -----------------------------------
  !
  call bintrp                       !that gets 1/b and its derivatives
  !				!and puts them into transp common
  call jintrp                       !that gets j and its derivatives
  !				!and puts them into transp common
  if(bijac(ion).le.0.0_fp) then
     write(nonlin,*) ' ?? orbg -- jintrp error, jacobian <= 0'
     ierr=1
     return
  end if
  !
  zg2    = big(ion)
  zpsip2 = bipsi(ion)
  zphi2  = biphi(ion)
  !
  zdphi = zphi2 - zphi1     !d(phi) in statvolts
  zdb = 1._fp/bi1ovb(1) - 1._fp/zb  !d(b) in gauss
  !
  if(ildiff)	then    !use differential differences
     !
     bidb(ion) = - zd1ovb/(zb*bi1ovb(ion))
     !
  else                      !use finite differences
     !
     bidb(ion)   = zdb
     bidphi(ion) = zdphi
     zdg         = zg2 - zg1
     zdpsip      = zpsip2 - zpsip1
     zd1ovb      = bi1ovb(1) - zb
     !
  end if
  !
  !     compute the optimum d(v_parallel) which conserves ke
  !
  if(nchdvp.eq.1) then      !use p_phi conservation
     !
     zzdvp = (bipphi(ion) + (ze/zc)*zpsip2)/ &
          (bim(ion)*zg2*bi1ovb(ion)) &
          - bivpar(ion)
     !
  elseif(nchdvp.eq.2) then  !use energy conservation
     !
     zzdvp = zdvopt(zdummy)
     !
  elseif(nchdvp.eq.3) then  !use v_parallel equation
     !
     zzdvp = bidvpa(ion)
     !
  else
     write(nonlin,68001) nchdvp
68001 format(/ &
          ' ?orbg -- invalid choice for d_vparallel calc.:',i6/)
     ierr=1
     return
  end if
  !
  bivpar(ion) = bivpar(ion) + zzdvp
  !
  bivprp(ion) = sqrt(2.0_fp*bimu(ion)/(bim(ion)*bi1ovb(ion)))
  !
  !     limit bivpar to .le. biv in magnitude
  !
  if(zphi1.eq.zphi2) then   !conserve biv exactly (phi constant)
     !
     zvlim = abs(biv(ion))
     bivpar(ion) = max(-zvlim,min(zvlim,bivpar(ion))) !prevent trouble
     !
  else                      !conserve energy including phi
     !
     biv(ion) = sqrt(bivpar(ion)*(1.0_fp*bivpar(ion))& !forced &
          + bivprp(ion)*(1.0_fp*bivprp(ion))) !dbl.p.
     !
  end if
  !
  !
  !     ---------------------------
  !     now compute biksid
  !     ---------------------------
  !
  biksid(ion) = bivpar(ion)/biv(ion)
  !
  if(biksid(ion)**2.gt.1.0_fp)  then
     write(nonlin,*) ' ?orbg: vpll/v > 1: ',biksid(ion)
     ierr=1
     return
  end if
  !
  !
  !     ------------------------
  !     7.0     cleanup & return
  !       ------------------------
  !     debug stuf:
  !
  !     recompute canonical angular momentum
  !
  zmr2 = bim(ion)*bivpar(ion)*big(ion)*bi1ovb(ion) - &
       (ze/zc)*zpsip2
  bidpph(ion) = zmr2 - zmr1 !change in p_phi (gm*cm^2/sec))
  !xx   zdmr = bidpph(ion)/&
  !xx         max(abs(zmr1),abs(zmr2),epslon)
  !!zdmr is relative change in canonical angular momentum
  !
  !     total energy:
  !
  ztotal2 = 0.5_fp*bim(1)*bivpar(1)**2 + bimu(1)/bi1ovb(1) + &
       ze*zphi2
  bite(ion) = ztotal2
  bidte(ion) = ztotal2 - ztotal1 !change in total energy (ergs)
  !xx   zdtotal_rel = bidte(ion)/ztotal1
  !***
  !
  !
  return
end subroutine orbg



