real(fp) function zdtrob(ierr)
  !
  !
  !  jfs  21 sep 94 -- massive clean-out of no-longer-used code:  keep only
  !                           code for ndtorb=1 and =4
  !
  !  jfs  23 aug 94 -- multipy by time-step magnification factor
  !                           (last thing before return)
  !
  !       bob mccann, 11-apr-86
  !	changes coming soon...rtm 24 feb.86
  !	rtm: min/max time step limits added may 1990 (sects. 2.0 & 3.0)
  !       pam stubberfield 15-jul-94 theta indexing for gxtabl change for uda
  !
  !	there are two ways to pick dsmax:
  !	vary it while plotting some orbits:
  !	  largest passing orbit (xi=1.28,e=30kev,xksid=1.0,goose=1.0)
  !	  large trapped orbit   (xi=1.28,e=30kev,xksid=???,goose=1.0)
  !	  tiny trapped orbit    (xi=1.10,e=30kev,xksid=???,goose=1.0)
  !	plot j and 1/b and look at the characteristic size of the
  !	structure of the highest derivatives used on orbg:
  !	  the second derivatives of 1/b and
  !	  the first derivatives of j
  !
  !	clearly, it would be nice to have a way of computing the
  !	characteristic sizes of the structure in the derivatives of
  !	j and 1/b.
  !
  !	it would also be nice to have a way of estimating
  !	the size of a trapped orbit (dsmax must be small enough to
  !	allow a particle to transcribe it's orbit).
  !
  !       --------
  !       comments
  !       --------
  !
  !
  !	this function is based upon the empirical time step rob initially
  !    used in orb50 (which got transformed into orbc).  it is somewhat more
  !    complicated due to the observation that orbits near the plasma edge,
  !    for q>~3 at the edge, needed smaller time steps to close the orbits.
  !    existence of this function allows easier modification of the timestep
  !    by more complicated means if needed without the confusion of trying to
  !    "do it consistently" in many different places.
  !	several different time step picking methods were tried, and they all
  !    failed to be a practical as this one.  trying to use conservation of
  !    either energy of momentum failed due to the polynomial nature of the
  !    error function---they have severe wiggles for time steps near the size
  !    of interest.  using pseudo-curvature of the drift orbit and the assumption
  !    that integration error was proportional to the curvature failed because
  !    it wasted time giving very clear definition of the turning points while
  !    taking too large a step at the outside midplane crossings.  the degear
  !    method may work, but that was scrapped due to the large time it would take
  !    to try it out.
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  use nbi_com
  !
  use physconst_mod, only:
  implicit none
  integer,intent(out) :: ierr  ! status code, 0=ok
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !	bim(ion)	    !ion mass - grams
  !	bimu(ion)	    !magnetic moment of guiding center - cgs
  !	biq(ion)	    !qmhd at the current guiding center location
  !
  !	odsmax		    !maximum differential step size (def=0.1)
  !	odsmin		    !minimum differential step size (def=0.001)
  !	ndtorb		    !=1 => original zdtrob with global q
  !			    !=4 => original zdtrob with local q
  !	odshrink	    !amount by which to optionally shrink dt
  !	nlshrink            !external logical :: control over dt shrinking
  !	ondrel              !relative cancellation of descriminant for shrinking
  !	odpowtip            !power law for localizing banana tip shrinking
  !	odshrtip	    !amount by which to shrink dt at banana tip
  !	rmjion		    !r at ion's last position
  !	bixi(ion)	    !xi
  !	dxi(j,2)	    !dxi
  !	0.5*(omegag(j)+omegag(j+1))    !rotation frequency of j-th zone
  !	splg(4,jh)	    !spline cf's of g (r*bphi) j-th zone (1,1+1/dxi)
  !	bixidt(ion)	    !(d/dt)xi
  !	bithdt(ion)	    !(d/dt)th
  !	bixi2t(ion)	    !(d2/dt2)xi
  !	bith2t(ion)	    !(d2/dt2)th
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !	zdtrob	    !(1/2)*(rob's optimum time step for orbc)
  !		    !limited so that ds.le.dsmax
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !
  integer, parameter :: ion = 1
  !
  real(fp) :: zomeg
  !
  real(fp) :: zemin,zvmin,zvprp2,zv2,zv1
  integer :: j8,jx,ith,ingc
  real(fp) :: zqfudg,zdtold,zh,zth,zqfac
  real(fp) :: ztrans,zadj
  !
  !       ---------------
  !       data statements
  !       ---------------
  !=========================================
  !
  ierr = 0
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  zemin=ti(ngc,2)
  ingc=min(ledge,ngc)
  zomeg=half*(omegag(ingc)+omegag(ingc+1))
  zvmin=0.2_fp*zomeg*raxis
  !
  !  non-thermal ions with v=0 usually only occur in very fast
  !  (pathological ::?) rotators.
  !
  zvprp2 = two*bimu(ion)/(bim(ion)*bi1ovb(ion))
  zv2 = biv(ion)**2
  zv1 = biv(ion)
  if(zv1.lt.zvmin) then
     zv1=zvmin
  end if
  if(zv1*zv1*v2toev.lt.zemin) then
     zv1=zv1*sqrt(zemin/(zv1*zv1*v2toev))
  end if
  !
  !       -----------------------------------------------------
  !       1.0     compute the optimum time step using surface q
  !       -----------------------------------------------------
  !
  !
  !  set up to interpolate doub's gxtabl (1/(r*grad(xi)))/<1/(r*grad(xi))>
  !
  jx = int(bixi(ion)/dxi(lcentr,2)) + 1
  zh = bixi(ion) - real(jx-1)*dxi(lcentr,2)
  !
  !   theta index
  !
  zth=bithst(ion)
  if(abs(zth).gt.(1.01_fp*xpi)) then
     write(nonlin,*) ' ?zdtrob: theta out of range [-pi,pi]: ',zth
     ierr=1
     return
  end if
  !
  if(nlsym2b) then
     ith = 1 + int( nthzons *abs(zth)/xpi )
  else
     ith = 1 + int((zth-thbdy0)*nthzons/twopi)
  end if
  ith = min(nthzons,ith)
  !
  zqfac = gxtabl(ith,jx)
  zqfudg = max(0.2_fp,min(8.0_fp, abs(biq(ion)*zqfac)))
  !!don't let "q" get huge
  !
  zdtold = half*zqfudg*(bidtnv/zv1)* &
       (one + 3.0_fp*zvprp2/zv2)
  !
  if(zdtold.le.zero) then
     write(nonlin,80101) zdtold, zqfudg, bidtnv, biv(ion), zvprp2
80101 format(' zdtrob error: zdtrob negative!',/, &
           ' zdtrob=',1pe15.5,/, &
           ' zqfudg=',e15.5,/, &
           ' bidtnv=',e15.5,/, &
           ' biv(ion)=',e15.5,/, &
           ' zvprp2=',/,e15.5)
     ierr=1
     zdtrob = bidtnv/abs(biv(ion))	!dummy default
  else
     zdtrob = zdtold
  endif
  !  dmc -- timestep reduction in edge region
  !         unity at xi<=0.9  --- 0.25 at xi>=1.0
  !
  if(bixi(ion).gt.0.9_fp) then
     if(bixi(ion).ge.one) then
        zdtrob=0.25_fp*zdtrob
     else
        ztrans=(bixi(ion)-0.9_fp)/0.1_fp
        zadj=0.25_fp*ztrans+one*(one-ztrans)
        zdtrob=zadj*zdtrob
     end if
  end if
  !
  return
end function zdtrob
! 11sep2004 fgtok
