subroutine ptclsave(zxi,zth,zvion,zksid,zrmjion,zion,ztflux,zmu,zpphi,zmo)
  !
  !  compute conserved quantities or reference quantities for use on
  !  next timestep; see ptclrgen.f90
  !  dmc 16 nov 2001
  use nbi_com
 
  implicit none
  real(fp), intent(in) :: zxi,zth         ! flux coordinate position of ptcl
  real(fp), intent(in) :: zvion,zksid     ! v,vpll/v (cm/sec & dimensionless)
 
  real(fp), intent(out) :: zrmjion,zion   ! (r,z) position (cm)
  real(fp), intent(out) :: ztflux         ! flux surface (tor. flux, wb.)
  real(fp), intent(out) :: zmu            ! mu (cgs)
  real(fp), intent(out) :: zpphi          ! canonical momentum (cgs)
  real(fp), intent(out) :: zmo            ! mechanical momentum (cgs)
  real(fp) ::  zf1d,zpsint
 
  !--------------------------------------
 
  real(fp) zjdum(2,2)
 
  real(fp) zvpll,zvperp
  real(fp) zbi    ! 1/b, 1/gauss
 
  !--------------------------------------
 
  integer :: jx       ! radial zone index
 
  real(fp)    :: zx0      ! xi at radial zone for spline interpolation
  real(fp)    :: zh       ! offset into radial zone
 
  real(fp)    :: zkonst   ! toroidal flux / pi with normalization
 
  real(fp)    :: zg       ! (r*bt) (from spline)
  real(fp)    :: zpsi     ! poloidal flux (from spline integration)
 
  real(fp)    :: zb       ! z of beam ion
  real(fp) zdum1,zdum2,zdum3,zdum4,zdum5
  !--------------------------------------
  ! (r,z) position
 
  if(nfast_map.eq.1) then
     call nbmomry(zxi,zth,0,zrmjion,zion,zjdum)  ! use bilinear maps
  else
     call nbmomry_xs(zxi,zth,0,zrmjion,zion,zjdum)  !use xplasma spline maps
  endif
  ztflux=tflux*zxi*zxi
 
  ! (1/b)
 
  call nbsplb(zxi,zth,0,zbi,zdum1,zdum2,zdum3,zdum4,zdum5)
 
  ! mu
 
  zvpll=zksid*zvion
  zvperp=sqrt(max(0.0_fp,zvion*zvion-zvpll*zvpll))
  zmu = 0.5_fp*zmp*abeam*zvperp*zvperp*zbi
 
  ! get g (r*bt) & psi for pphi
 
  jx  = int(zxi/dxi(lcentr,2)) + 1
  zx0 = real((jx-1),fp)*dxi(lcentr,2)
  zh  = zxi - zx0
 
  zkonst = 1.e8_fp*tflux/xpi
 
  zg = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  zpsi = zpsint(plflxi(jx),zkonst,zx0,zh,spliot(1,jx),  &
       spliot(2,jx),spliot(3,jx),spliot(4,jx))
 
  ! momentum
 
  zb = xzbeams(lsbeam)
  zmo = abeam*zmp*zvpll*zg*zbi
  zpphi = zmo - zb*zel*zpsi/zc
 
  return
end subroutine ptclsave
! 11sep2004 fgtok -s r8_precision.sub nubeam_test_r8.sub "r8con.csh conversion"
! 11sep2004 fgtok
