real(fp) function zdvopt(zdummy)
  !
  !
  !       bob mccann, 11-apr-86
  !
  !	rtm, march-87: fixed round-off trouble case zround=0
  !	rtm, 30-jan-87: fixed double root at zero when zzb=zzc=0.0
  !	rtm, dec-86: fixed nonphysical case zdescr<0.0
  !
  !
  !       --------
  !       comments
  !       --------
  !
  !
  !	this function was created in order to isolate the orbit routines
  !    from the considerable numerical subtlety involved in picking the "best"
  !    value for dv based upon conservation of energy.  note that even with
  !    the best numerical methods i know, it is still possible that energy
  !    will not be conserved at the banana tips.  that was too much to ask.
  !    it turns out that doing that based upon the sign of db caused orbits
  !    to "get stuck" on a banana tip which then drifted out of the plasma
  !    on a mod(b) surface; the
  !    integration error in xi and theta put the orbit on the wrong side of
  !    the mod(b) surface, giving db the wrong sign; the computed dv put the
  !    orbit back into the same condition as the previous step; etc., etc.
  !    always picking the root nearest the one computed from the dv equation
  !    seems to avoid this difficulty.
  !
  !    		note that double precision is not needed on the cray!
  !
  !	the use of double precision is intended to result in getting the
  !    "correct" answer to full machine precision on the vax.
  !
  !       please note that this
  !    still is no more accurate than the accuracy in d(b) and d(phi)
  !    which have fewer and fewer significant digits as the step size
  !    deminishes, whereas the dv equation answer still has full
  !    precision since it uses the derivatives of b and phi directly.
  !    however, the dv equation answer does not conserve energy exactly
  !    for any stepsize, as there is an error which is proportional to
  !    dt**3.
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  use iso_c_binding, only: fp => c_double
  use nbi_com
  use physconst_mod, only: zone=>one,ztwo=>two
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !	bim(ion)	    !ion mass - cgs
  !	bimu(ion)	    !magnetic moment of the guiding center
  !	bivpar(ion)	    !v_parallel
  !	bidvpa(ion)	    !d(v_parallel) computed from v_parallel equation
  !	bidb(ion)	    !(b_new - b_old) in gauss
  !	bidphi(ion)	    !(phi_new - phi_old) in statvolts
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !	zdvopt	    !dv nearest bidvpa(ion) which conserves ke
  !
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !
  !============
  implicit none
  integer :: ion
  !============
  real(fp) :: zdummy,ze,zround,zzdvp1,zzdvp2,zd1,zd2
  !============
  parameter   (ion=1)
  !
  real(fp) ::	zza,zzb,zzc,zdescr,zzden
  real(fp) ::	zfour
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  !
  ze = zel*biz(ion)
  !
  !
  !       --------------------------------------------------------
  !       1.0     compute the optimum d(v_parallel) to conserve ke
  !       --------------------------------------------------------
  !
  !
  !
  !	rtm, july 1985
  !
  !    here i make use of energy conservation:
  !    de = e*d(phi) + mu*d(b) + m*vp*d(vp) + 0.5*m*d(vp)**2 == 0.0,
  !    first mult. by 2/m, then solve for d(vp):
  !
  !    n.b.  the solutions will be as accurate as bidb(ion) = b(new)-b(old)
  !	   for large time steps, the taylor expansion doesn't work.
  !	   for very small steps, the direct difference fails.
  !	   it is concievable that a weighted combination would be
  !	   better, but the direct difference seems good enough for now.
  !
  zza = zone				!*d(vp)**2
  zzb = ztwo*bivpar(ion)			!*d(vp)
  zzc = ztwo*(bimu(ion)*bidb(ion) + &
       ze*bidphi(ion))/bim(ion)	  	!constant term
  !
  !		***dec-86; rtm***
  !
  !    treat troublesome cases:
  !
  zfour=ztwo*ztwo
  zdescr = zzb**2 - zfour*zza*zzc		!avoid sqrt(negative)
  !
  zround = abs(bidb(ion)) + abs(bidphi(ion)) !avoid inaccuracy
  !
  if(zdescr.lt.zero)	then	!nonphysical location
     !
     zdvopt = bidvpa(ion)			!use (d/dt)v_parallel
     return			!done
     !
  elseif(zround.eq.0.0_fp.and. zzb.ne.0.0_fp)	then	!probable round-off trouble
     !
     zdvopt = bidvpa(ion)			!use (d/dt)v_parallel
     return 			!done
     !
  end if
  !
  !   prepare to solve the quadratic equation
  !
  zzden  = abs(zzb) + sqrt(zdescr)	!avoid div. by zero
  !
  if(zzden.eq.zero)	then	!zzb=zzc=0 => double root at 0
     !
     zdvopt = 0.0_fp
     return		!done
     !
  else				!solve the quadratic equation
     !
     !    note that there are two solutions:
     !
     !    the following is newton's method of expressing the roots of
     !    the quadratic equation.  this is the best numerically conditioned
     !    way of computing the roots that i know.  rtm, july 1985.
     !
     zzdvp1 = -ztwo*zzc*sign(zone,zzb)/zzden
     !
     zzdvp2 = -zzden*sign(zone,zzb)/(ztwo*zza)
     !
     !    choose between them:
     !
     !	pick the one which is closest to bidvpa(ion) computed in the v-parallel
     !    equation.  n.b. picking based upon sign of db and v-parallel leads to
     !    orbits getting stuck to mob(b) surfaces at banana tips!
     !
     zd1 = abs(zzdvp1-bidvpa(ion))
     zd2 = abs(zzdvp2-bidvpa(ion))
     !
     if(zd1.le.zd2)	then
        zdvopt = zzdvp1		!first root is closer
     else
        zdvopt = zzdvp2		!second root is closer
     end if
     !
  end if
  !
  !
  !       ------------------------
  !       7.0     cleanup & return
  !       ------------------------
  !
  !
  !
 return
end function zdvopt
