subroutine chkphip
  use nbi_com
  use nbi_com, only : zpi=>xpi
  use nbi_dimensions
  !
  !  debug routine -- look at toroidal precession of trapped orbits
  !
  !  user interactively selects a surface, and for fixed energy can
  !  plot precession as a function of bounce angle, or, for fixed bounce
  !  angle plot precession as a function of energy.
  !
  !  an approximate expression for precession can be compared with results
  !  of explicit orbit integration.
  !
  !------------------------
  !  local arrays for plots
  !
  !============
  implicit none
  integer :: imaxpl,igbd,inth,ine,luntrm,lnout,isurf,idcod,iert
  integer :: isb,iab,izb,inmax,idir,iaflag,ie,it,iaoff
  !============
  real(fp) :: zfacdt,ze,zth,zthmin,zthmax,zemin,zemax,zxi,ze1
  real(fp) :: ze2,zth1,zth2,zbth,zthx,zpp,zrat,zfac,z_el,zv,zrlara
  real(fp) :: zrho,zf,zthl
  !============
  parameter (imaxpl=1000)
  real(fp) :: zepl(imaxpl),zthpl(imaxpl),zrlpl(imaxpl),zrhpl(imaxpl)
  real(fp) :: zphip(imaxpl),zphipr(imaxpl)
  real(fp) :: zx(imaxpl),zy(imaxpl),zy2(imaxpl)
  real(fp) :: zbtime(imaxpl)
  integer :: iamin(imaxpl)
  !
  character(len=30) :: zlabx,zlaby,zunsx,zunsy,zzzlbl
  !
  character(len=60) :: zrunid,zchkid
  !
  logical :: lxout,lyes,lchar
  !
  logical :: idbgorb,iloge,ilvphi,ilbtime
  !
  !------------------------
  !  defaults
  !
  zzzlbl=' '
  !
  ilbtime=.false.
  !
  zfacdt=1.0_fp
  !
  igbd=0
  !
  ilvphi=nlvphi
  !
  idbgorb=.false.
  !
  ze=1000.0_fp  ! default orbit energy, ev
  zth=1.57_fp  ! default bounce angle, approx pi/2.
  !
  inth=100
  zthmin=-zpi
  zthmax=zpi
  !
  ine=100
  iloge=.true.
  zemin=10.0_fp
  zemax=10.0e6_fp
  !
  !  select surface or quit
  !
  luntrm=lnout(0)
10 continue
  if(lxout(0)) write(luntrm,1001) nphipc
1001 format(' select numeric surface btw 1 and ',i3,' or "q":')
  call uread('chkphip:  select xi surface:$')
  if(lchar('q')) then
     nlvphi=ilvphi
     return
  end if
  isurf=idcod(iert)
  if(iert.ne.0) isurf=0
  if((isurf.le.0).or.(isurf.gt.nphipc)) then
     call uermsi('surface index out of range',isurf,0)
     goto 10
  end if
  !
  !  select species or back up
  !
12 continue
  if(lxout(0)) write(luntrm,1002) nsbeam
1002 format(' select fast ion species btw 1 and ',i1,' or "q":')
  isb=0
  call ureadi('chkphip:  species select:$',isb)
  if(isb.eq.0) goto 10
  if((isb.lt.0).or.(isb.gt.nsbeam)) then
     call uermsi('species index out of range',isb,0)
     goto 12
  end if
  !
  iab=abeams(isb)+0.5_fp
  izb=xzbeams(isb)+0.5_fp
  !
  call uread('chkphip:  enter "y" to disable exb effects:$')
  if(lyes(0)) then
     nlvphi=.false.
  end if
  !  species selection done
  !
  lsbeam=isb
  call inispc
  !
  !  xi value at surface
  !
  zxi=real(isurf,fp)/real(nphipc,fp)
  write(zrunid,'(a,1x,'' xi='',1pe11.4)') runid,zxi
  !
  !  select plot type
  !
20 continue
  if(idbgorb) then
     inmax=20
  else
     inmax=imaxpl
  end if
  if(lxout(0)) write(luntrm,1003) idbgorb,zfacdt,igbd,ilbtime,zzzlbl
1003 format(' plot type option:'/ &
       ' d -- debug switch:  plot each orbit (currently ',l1,')'/ &
       ' e -- plot vs. energy, bounce angle fixed.'/ &
       ' f -- set orbit timestep reduction factor, currently: ',1pe11.4/ &
       ' g -- grad(b) set phi drift multiplier (currently ',i2,')'/ &
       ' b -- plot vs. bounce angle, energy fixed.'/ &
       ' t -- toggle bounce time plot switch:  now ',l1/ &
       ' l -- set (30 char) label to appear before runid on plots.'/ &
       '      currently set to:  ',a/ &
       ' q -- quit, or reset surface and/or species index'/ &
       ' i -- test interpolations.')
  call uread('chkphip:  enter plot type option:$')
  if(lchar('q')) then
     goto 10
  else if(lchar('i')) then
     call chkphipi(zzzlbl)
  else if(lchar('t')) then
     ilbtime=.not.ilbtime
  else if(lchar('d')) then
     idbgorb=.not.idbgorb
  else if(lchar('l')) then
     call ureadl(' enter label string:$', zzzlbl)
  else if(lchar('f')) then
     call ureadr8(' orbit timestep adjust factor, btw 0.2 & 1:$', &
          zfacdt)
     zfacdt=max(0.2_fp,min(1.0_fp,zfacdt))
  else if(lchar('g')) then
     call uread(' enter "+", "-" or "0":$')
     if(lchar('+')) then
        igbd=1
     else if(lchar('-')) then
        igbd=-1
     else if(lchar('0')) then
        igbd=0
     else
        call uermsg('invalid input',' ',0)
     end if
  else if(lchar('e')) then
     !
     call ureadr8( &
          'chkphip: enter fixed bounce angle, btw -pi and pi:$',zth)
     zth=max(-zpi,min(zpi,zth))
     write(zchkid,2201) iab,izb,zth
2201 format('chkphip test plot:  a=',i2,' z=',i2,' th(bounce)=', &
          1pe10.3)
     if(lxout(0)) write(luntrm,'('' nmax='',i5)') inmax
     call ureadi('chkphip:  no. of e pts, btw 2 and nmax:$',ine)
     ine=max(2,min(inmax,ine))
     !
22   continue
     ze1=zemin
     call ureadr8('chkphip:  min energy, .gt. 1 ev:$',ze1)
     ze2=zemax
     call ureadr8('chkphip:  max energy, .lt. 1.0e8 ev:$',ze2)
     ze1=max(1.0_fp,min(1.0e8_fp,ze1))
     ze2=max(1.0_fp,min(1.0e8_fp,ze2))
     if(ze1.eq.ze2) then
        call uermsr8('invalid energy range, e1=e2',ze1,0)
        goto 22
     end if
     zemin=min(ze1,ze2)
     zemax=max(ze1,ze2)
     !
     call uread('chkphip:  enter "y" for log spacing of e pts:$')
     if(lyes(0)) then
        iloge=.true.
     else
        iloge=.false.
     end if
     goto 100
     !
  else if(lchar('b')) then
     !
     call ureadr8( &
          'chkphip: enter fixed energy, ev:$',ze)
     ze=max(1.0_fp,min(1.0e8_fp,ze))
     write(zchkid,2204) iab,izb,ze
2204 format('chkphip test plot:  a=',i2,', z=',i2,' e=', &
          1pe10.3,' ev.')
     if(lxout(0)) write(luntrm,'('' nmax='',i5)') inmax
     call ureadi('chkphip:  no. of th pts, btw 2 and nmax:$',inth)
     inth=max(2,min(inmax,inth))
     !
24   continue
     zth1=zthmin
     call ureadr8('chkphip:  min bounce angle (th), .gt. -pi:$', &
          zth1)
     zth2=zthmax
     call ureadr8('chkphip:  max bounce angle (th), .lt. -pi:$', &
          zth2)
     zth1=max(-zpi,min(zpi,zth1))
     zth2=max(-zpi,min(zpi,zth2))
     if(zth1.eq.zth2) then
        call uermsg('invalid angle range, th1=th2',zth1,0)
        goto 24
     end if
     zthmin=min(zth1,zth2)
     zthmax=max(zth1,zth2)
     goto 200
     !
  else
     call uermsg('?chkphip:  unrecognized plot type option.',' ',0)
  end if
  goto 20
  !
100 continue
  !
  !  first eval roscoe's energy independent p(psi) integral
  !
  bixi(1)=zxi
  bithst(1)=zth
  call bintrp
  zbth=bimodb(1)
  call bphip_phip(isurf,zth,zbth,zthx,zpp,idir,iaflag)
  !
  zrat=zemax/zemin
  do ie=1,ine
     zfac=real(ie-1,fp)/real(ine-1,fp)
     if(iloge) then
        z_el=zemin*zrat**zfac
     else
        z_el=zemin+(zemax-zemin)*zfac
     end if
     zepl(ie)=z_el
     call chkphip_std(z_el,zv,zrlara,zrho)
     zrlpl(ie)=zrlara
     zrhpl(ie)=zrho
     !
     call chkphip_orb(isurf,idbgorb,igbd,zth,zv,zphip(ie),zfacdt, &
          zbtime(ie),zrho)
     zphipr(ie)=zphip(ie)/zrhpl(ie)
     zy2(ie)=zpp  ! roscoe's
  end do
  !
105 continue
  if(lxout(0)) write(luntrm,1050)
1050 format(' enter "p" for phi(p) plot, "r" for phi(p)/rho plot.')
  call uread('chkphip:  enter "p" or "r":$')
  if(lchar('p')) then
     call chkphip_p(zphip,zy,ine,'phi(p)','radians',zlaby,zunsy)
     do ie=1,ine
        zy2(ie)=zrhpl(ie)*zy2(ie)
     end do
  else if(lchar('r')) then
     call chkphip_p(zphipr,zy,ine, &
          'phi(p)/rho','radians',zlaby,zunsy)
  else if(lchar('q')) then
     goto 105
  else
     call uermsg(' ? enter "p" or "r"',' ',0)
     goto 105
  end if
  !
106 continue
  if(lxout(0)) write(luntrm,1060)
1060 format( &
       ' enter "e", "r", "l" for plot vs. energy, rho, larmour radius.')
  call uread('chkphip:  enter "e", "r", or "l":$')
  if(lchar('e')) then
     call chkphip_p(zepl,zx,ine,'energy','ev',zlabx,zunsx)
  else if(lchar('r')) then
     call chkphip_p(zrhpl,zx,ine,'rho',' ',zlabx,zunsx)
  else if(lchar('l')) then
     call chkphip_p(zrlpl,zx,ine, &
          'axial larmour radius','cm',zlabx,zunsx)
  else if(lchar('q')) then
     goto 106
  else
     call uermsg(' ? enter "e" or "r" or "l"',' ',0)
     goto 106
  end if
  !
  call r8_grafx2(zx,zy,zy2,ine,zunsx,zunsy,zzzlbl//zrunid,zchkid, &
       zlaby//' vs. '//zlabx//' (roscoe value dashed)')
  !
  if(ilbtime) then
     call r8_grafx1(zx,zbtime,ine,zunsx,'seconds',zzzlbl//zrunid, &
          zchkid,'bounce time (orbg) vs. '//zlabx)
  end if
  !
  goto 20
  !
200 continue
  call chkphip_std(ze,zv,zrlara,zrho)
  do it=1,inth
     zf=real(it-1,fp)/real(inth-1,fp)
     zthl=zthmin+zf*(zthmax-zthmin)
     zthpl(it)=zthl
     !
     !  roscoe integral
     bixi(1)=zxi
     bithst(1)=zthl
     call bintrp
     zbth=bimodb(1)
     call bphip_phip(isurf,zthl,zbth,zthx,zpp,idir,iaflag)
     zy2(it)=zpp*zrho
     iamin(it)=iaflag
     if(iaflag.ne.0) then
        write(6,8801) zthl
8801    format(' %chkphip:  deeply trapped at theta=',1pe11.4)
     end if
     !
     !  orbit integrator
     call chkphip_orb(isurf,idbgorb,igbd,zthl,zv,zphip(it),zfacdt, &
          zbtime(it),zrho)
     !
  end do
  !
  !  patch the deeply trapped regions
  !
  iaoff=ndatspl(isurf)-1
  call bphip_minp(zthpl,zy2,iamin,inth,iaoff)
  !
  call r8_grafx2(zthpl,zphip,zy2,inth,'radians','radians', &
       zzzlbl//zrunid, &
       zchkid,'phi(p) (precession) vs. bounce angle -- roscoe dashed')
  !
  if(ilbtime) then
     call r8_grafx1(zthpl,zbtime,inth,'radians','seconds', &
          zzzlbl//zrunid, &
          zchkid,'bounce time vs. bounce angle (orbg)')
  end if
  !
  goto 20
  !
end subroutine chkphip
!-------------------------------
!
subroutine chkphip_std(ze,zv,zrlara,zrho)
  use iso_c_binding, only: fp => c_double
  use nbi_com
  use nbi_dimensions
  !
  !  from the energy evaluate the speed, axial vperp=v larmour radius
  !  (zrlara) and the normalized larmour radious rho=(zrlara)/(mjr radius)
  !
  !============
  implicit none
  real(fp) :: zv,zrlara,zrho,ze,zv2,zrb
  !============
  zv2=ze/v2toev
  zv=sqrt(zv2)
  !
  zrb=bzxr                          ! vacuum r*b on axis
  !
  zrho=zv*abeam/(9.578e7_fp*xzbeami*zrb)  !  rlar/raxis
  zrlara=raxis*zrho                    !  rlar
  !
  return
end subroutine chkphip_std
!--------------------------------
!
subroutine chkphip_p(zain,zaout,inum,zli1,zli2,zlo1,zlo2)
  use iso_c_binding, only: fp => c_double
  !
  !  copy data and labels
  !
  !============
  implicit none
  integer :: inum,ii
  !============
  real(fp) :: zain(inum),zaout(inum)
  character(len=*) zli1,zli2,zlo1,zlo2
  !
  !---------
  !
  zlo1=zli1
  zlo2=zli2
  !
  do ii=1,inum
     zaout(ii)=zain(ii)
  end do
  !
  return
end subroutine chkphip_p
!
!------------------------------------------------------------------------
!
subroutine chkphip_orb(isurf,idbg,igradbd,zth,zv,zphip,zfacdt, &
     zbtime,zrho)
  use iso_c_binding, only: fp => c_double
  use nbi_com
  use nbi_dimensions
  !
  !  this subroutine follows one full trapped orbit, given as input
  !  the orbit parameters at the bounce point (v=vperp, vpll=0)
  !
  !  if idbg=t, plot the orbit trajectory.
  !
  !  isurf identifies the starting flux surface
  !  zxi = float(isurf)/float(nphipc)
  !
  !  (zxi,zth) -- mag. flux starting coordinates.
  !  (zv) -- ion speed, vpll/v = 0 assumed initially
  !
  !  igradbd=-1,0,1 -- optional factor for correction for grad(b)
  !                    drift in toroidal direction
  !
  ! output:
  !
  !  (zphip) = toroidal precession arc per half bounce
  !          = 1/2 the precession for the whole bounce followed here.
  !            (radians)
  !
  !  zbtime  = bounce time (secs)
  !
  !============
  implicit none
  integer :: igradbd,intot,imaxpl,i,itrm,jx,ilinx,ilinth,ier1,ierr
  integer :: istepo,ijac,ibsign,ibdel,iprev,ibounce,idir,ibs,ii
  !============
  real(fp) :: zfacdt,zbtime,zrho,zgbthdot,zgb1,zgb2,zgb,zcb,zpi,zf1d
  real(fp) :: zzh,zf0,zf1,zf2,zf3,zdf1d,zpsint,zzprev,zzkonst,zzx0
  real(fp) :: zkonst,zdbdxi,zx0,zh,zg,zi,zdidxi,zrb,zpsi0,zdthcon
  real(fp) :: zgbfac,zthasum,zthxsum,zbmax,zthstart,zphipgb,zphipg2
  real(fp) :: zintg,zbmodp,zintgp,zxip,zthp,zpsi,zvpll,zgbfacp
  real(fp) :: zdphigb,zd,zgx2,zr2i,zgxgt,zdphig1,zdphig2,zgbtmp
  real(fp) :: zcbtmp,zdphib,zphdr,zcorr,zthcur,zphib,zthtr,zxiav
  real(fp) :: zxiavo
  !============
  common/zorbgtst/ zgbthdot,zgb1,zgb2,zgb,zcb         ! temp for debugging
  !
  integer :: isurf
  logical :: idbg
  real(fp) :: zxi,zth
  real(fp) :: zv
  real(fp) :: zphip
  !
  parameter (intot=400)
  real(fp) :: ztha(intot),zterm(intot,5)
  real(fp) :: ziterm(5)
  common/zbphipg2/ ziterm,ztha,zterm
  !
  !-----------------------
  !  local
  !
  parameter (imaxpl=1000)
  !
  real(fp) :: zr(imaxpl),zy(imaxpl)
  real(fp) :: ztime(imaxpl),zepl(imaxpl),zmupl(imaxpl),zvvpl(imaxpl)
  real(fp) :: zvplla(imaxpl),zvpllb(imaxpl),zvpllc(imaxpl)
  real(fp) :: zxipl(imaxpl),zthpl(imaxpl)
  !
  real(fp) :: zigpl(intot),zigth(imaxpl),zthipl(imaxpl)
  !
  integer :: inumpl
  !
  logical :: idbgl
  !
  logical :: lyes,lchar
  !
  real(fp) :: zsum
  !
  integer :: ideriv(6)
  data ideriv/0,0,1,0,0,0/
  data zpi/3.141593_fp/
  !
  !  statement fcn
  !
  zf1d(zzh,zf0,zf1,zf2,zf3) &
       = ((zf3*zzh + zf2)*zzh + zf1)*zzh + zf0
  !
  !
  zdf1d(zzh,zf1,zf2,zf3) = (3.0_fp*zf3*zzh + 2.0_fp*zf2)*zzh + zf1
  !
  !  poloidal flux integral formula, from iota(bar) spline.
  !  based on d(psi)/d(xi) = iota(bar)*xi*(tor.flux.at.bdy)/xpi
  !
  zpsint(zzprev,zzkonst,zzx0,zzh,zf0,zf1,zf2,zf3) = zzprev + &
       zzkonst * zzh * (zzx0*zf0 + zzh*((zf0 + zzx0*zf1)/2.0_fp+ &
       zzh*((zf1+zzx0*zf2)/3.0_fp+ zzh*((zf2 + zzx0*zf3)/4.0_fp+ &
       zzh*zf3/5.0_fp))))
  !
  !-----------------------
  !  initialize
  !
  idbgl=.false.
  if(idbg) then
     write(6,'('' th='',1pe11.4)') zth
     call uread('make debug plots? (y/n):$')
     if(lyes(0)) then
        idbgl=.true.
        do i=1,intot
           zigpl(i)=0.0_fp
           do itrm=1,5
              zigpl(i)=zigpl(i)+zterm(i,itrm)*zrho
           end do
        end do
     end if
  end if
  !
  zxi=real(isurf,fp)/real(nphipc,fp)
  !
  zkonst=1.e8_fp*tflux/zpi
  !
  xion=zxi
  th=zth
  !
  !  evaluate b, 1/b from bifourier spline
  !  also evaluates derivatives for first step in orbit eqn.
  !
  bixi(1)=xion
  bithst(1)=th
  call bintrp
  call jintrp
  zdbdxi=bibdxi(1)
  !
  jx = int(bixi(1)/dxi(lcentr,2)) + 1
  zx0 = real(jx-1,fp)*dxi(lcentr,2)
  zh = bixi(1) - zx0
  zg = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  zi = zf1d(zh,splimhd(1,jx), &
       splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  zdidxi = zdf1d(zh,splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  zrb=1.0e-4_fp*zg   ! tesla*cm
  !
  zpsi0=zpsint(plflxi(jx),zkonst,zx0,zh,spliot(1,jx), &
       spliot(2,jx),spliot(3,jx),spliot(4,jx))
  !
  ilinx=1
  ilinth=1
  call bcspeval(xion,th,ideriv,zdthcon, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspthets),nzp1,ier1)
  !
  zgbfac=zi/zg*zdthcon
  !
  inumpl=1
  istepo=1
  ijac=1
  call nbmomry(xion,th,ijac,rmjion,yion,xbjaco)
  zr(1)=rmjion
  zy(1)=yion
  zxipl(1)=xion
  zthpl(1)=th
  !
  zthasum=0.0_fp
  zthxsum=0.0_fp
  !
  vion=zv
  vpll=0.0_fp
  vperp=vion
  xksid=0.0_fp
  !
  ztime(1)=0.0_fp
  zepl(1)=v2toev*vion*vion
  zmupl(1)=zepl(1)         ! v2toev*vperp*vperp*bmax/b ~ mu
  zvvpl(1)=zepl(1)         ! (vion+vpll(pphi error))*vion*v2toev
  zbmax=bimodb(1)
  !
  zvplla(1)=vpll   ! will be from orbg
  zvpllb(1)=vpll   ! will be from pphi conservation
  zvpllc(1)=vpll   ! vpll=vion*sqrt(1-b/bmax)*sign(vpll)
  !                      ! from e & mu conservation
  !
  !-------
  !  integrate the orbit
  !
  zbtime=0.0_fp
  !
  zthstart=th
  !
  zphip=0.0_fp
  zphipgb=0.0_fp
  zphipg2=0.0_fp
  !
  zsum=0.0_fp
  ibsign=0
  ibdel=0
  iprev=0
  !
  ibounce=0
  !
  zintg=vpll*zrb/(bimodb(1)*rmjion*rmjion)
  !
100 continue
  !
  zbmodp=bimodb(1)
  zintgp=zintg
  !
  zxip=xion
  zthp=th
  !
  call orbpre
  !
  bidtnv=zfacdt*bidtnv              ! option to shorten timestep
  !
  call orbg(ierr)
  !
  call orbout
  !
  jx = int(bixi(1)/dxi(lcentr,2)) + 1
  zx0 = real(jx-1,fp)*dxi(lcentr,2)
  zh = bixi(1) - zx0
  zg = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  zi = zf1d(zh,splimhd(1,jx), &
       splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  zdidxi = zdf1d(zh,splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  zrb=1.0e-4_fp*zg   ! tesla*cm
  !
  zpsi=zpsint(plflxi(jx),zkonst,zx0,zh,spliot(1,jx), &
       spliot(2,jx),spliot(3,jx),spliot(4,jx))
  !
  zvpll=(zpsi-zpsi0)*9578.0_fp*xzbeami*bimodb(1)/(abeam*zrb)
  !
  !x      write(6,1999) vpll,zvpll
  !x 1999 format(' chkphip_orb:  vpll(orbg)=',1pe11.4,' vpll(psi)=',1pe11.4)
  !
  call bcspeval(xion,th,ideriv,zdthcon, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspthets),nzp1,ier1)
  !
  zgbfacp=zgbfac
  zgbfac=zi/zg*zdthcon
  !
  zdphigb=zgbthdot*0.5_fp*(zgbfac+zgbfacp)
  !
  call nbmomry(xion,th,ijac,rmjion,yion,xbjaco)
  !
  zd=xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1)
  zgx2=(xbjaco(1,2)**2+xbjaco(2,2)**2)/(zd*zd)
  !
  zr2i=1.0_fp/rmjion**2
  !
  zgxgt=(-xbjaco(2,1)*xbjaco(2,2)-xbjaco(1,1)*xbjaco(1,2))/(zd*zd)
  !
  zdphig1=zgb1*zr2i*zgx2
  zdphig2=zgb2*zr2i*zgxgt
  zdphigb=zdphig1+zdphig2
  !
  zphipgb=zphipgb+zdphigb*bidt(1)
  !
  zgbtmp=zi*zgb*zdbdxi/zg
  zcbtmp=zdidxi*zcb
  zphipg2=zphipg2+(zgbtmp+zcbtmp)*bidt(1)
  zdbdxi=bibdxi(1)

  !
  zintg=vpll*zrb/(bimodb(1)*rmjion*rmjion)
  !
  if(zthp.lt.th) then
     idir=1
     if((th-zthp).gt.zpi) idir=-1
  else
     idir=-1
     if((zthp-th).gt.zpi) idir=1
  end if
  !
  !  field line length along step
  !
  zthasum=zthasum+abs(th-zthp)
  zthxsum=zthxsum+0.5_fp*(xion+zxip)*abs(th-zthp)
  !
  call bphip_phib(isurf,zthp,th,zdphib,idir,0.0_fp)
  !
  zphdr=0.5_fp*(zintg+zintgp)*bidt(1)
  zcorr=sign(1.0_fp,zphdr)*zdphib
  zsum=zsum+zphdr-zcorr
  !
  if(bimodb(1).gt.zbmodp) then
     ibs=1
  else
     ibs=-1
  end if
  !
  if(ibs.ne.ibsign) then
     ibdel=ibdel+1
     ibsign=ibs
  end if
  !
  if(vpll.gt.0) then
     ii=1
  else
     ii=-1
  end if
  !
  if(ii.ne.iprev) then
     if(iprev.ne.0) then
        ibounce=ibounce+1
        zthcur=th
        zphib=0.0_fp
        if(iprev.lt.0) then
           if(idbgl) write(6,'('' ctr sum='',1pd11.4)') zsum
           zphip=zphip+igradbd*zphipg2+0.5_fp*(zsum+zphib)
        else
           if(idbgl) write(6,'(''  co sum='',1pd11.4)') zsum
           zphip=zphip+igradbd*zphipg2+0.5_fp*(zsum-zphib)
        end if
        zthstart=zthcur
        zsum=0.0_fp
        zphipgb=0.0_fp
        zphipg2=0.0_fp
     end if
     iprev=ii
  end if
  !
  istepo=istepo+1
  if(istepo.gt.10*inumpl) then
     write(nonlin,'('' ?chkphip:  too many orbit steps.'')')
     zphip=0.0_fp
     goto 200
  end if
  inumpl=min(imaxpl,(inumpl+1))
  !
  zr(inumpl)=rmjion
  zy(inumpl)=yion
  !
  zxipl(inumpl)=xion
  zthpl(inumpl)=th
  !
  zepl(inumpl)=v2toev*vion*vion
  zmupl(inumpl)=v2toev*vperp*vperp*zbmax/bimodb(1)
  zvvpl(inumpl)=(vion+zvpll-vpll)*vion*v2toev
  !
  zvplla(inumpl)=vpll
  zvpllb(inumpl)=zvpll
  zvpllc(inumpl)=vion*sqrt(max(0.0_fp,(1.0_fp-bimodb(1)/zbmax)))* &
       sign(1.0_fp,vpll)
  !
  ztime(inumpl)=ztime(inumpl-1)+bidt(1)
  zbtime=zbtime+bidt(1)
  !
  zthtr=0.5_fp*(th+zthp)
  zxiav=0.5_fp*(xion+zxip)
  !
  zthipl(inumpl)=zthtr
  zigth(inumpl)=(zphdr-zcorr)/(zdthcon*max(1.e-6_fp,abs(th-zthp)))
  !
  if(ibounce.eq.2) goto 200
  if(ibdel.gt.8) then
     write(6,'('' %chkphip_orb:  orbit did not trap.'')')
     zphip=0.0_fp
     goto 200
  end if
  if(xion.gt.1.29_fp) then
     write(6,'('' %chkphip_prb:  orbit escaped!'')')
     zphip=0.0_fp
     goto 200
  end if
  !
  !  return to orbit
  !
  goto 100
  !
  !  orbit done ................
  !
200 continue
  !
  zxiavo=zthxsum/max(epslon,zthasum)
  !
  if(idbgl) then
     !
     call r8_grafx1(zr,zy,inumpl,'cm','cm',runid, &
          'chkphip_orb debug plot','orbit trajectory')
     !
     call r8_grafx1(zthpl,zxipl,inumpl,'th','xi',runid, &
          'chkphip_orb debug plot','orbit trajectory')
     !
     call r8_grafx3(ztime,zvplla,zvpllb,zvpllc,inumpl, &
          'seconds','cm/sec',runid,'chkphip_orb debug plot', &
          'vpll:  orbg:solid, pphi-conserv, mu&e-conserv long dash')
     !
     call r8_grafx3(ztime,zepl,zmupl,zvvpl,inumpl, &
          'seconds','cm/sec',runid, &
          'chkphip_orb debug plot (k.e. solid, mu*bmax short dash)', &
          ' kin. energy, mu*bmax, m/2*vion*(vion-(vpll error))')
     !
     call r8_grafxf(ztha,zigpl,intot,zthipl,zigth,inumpl, &
          'th','dphip/dth',runid,'chkphip_orb debug plot', &
          ' integrand comparison (roscoe = solid)')
  end if
  !
  return
end subroutine chkphip_orb
!---------------------------------------------------------------------
! chkphipi -- check phi(p), phi(b) interpolation
!
subroutine chkphipi(zzzlbl)
  use iso_c_binding, only: fp => c_double
  use nbi_com
  use nbi_dimensions
  !
  !============
  implicit none
  integer :: imaxx,imaxt,inx,ix,imaxth,inumth,iengy,ith,idcod,iert
  integer :: ia
  !============
  real(fp) :: zpi,zthlimlo,zthlimhi,zrho,zdel,zev,zv2,zv,zth1,r8dcod
  real(fp) :: zth2,zphib,zpp,zthb
  !============
  character(len=*) zzzlbl
  !
  !-----------------------------------
  !  local:
  !
  parameter(imaxx=400)
  parameter(imaxt=16000)
  real(fp) :: zxi(imaxx)
  real(fp) :: zth(imaxx)
  real(fp) :: zwk(imaxx)
  !
  real(fp) :: zfb(imaxt)
  real(fp) :: zfp(imaxt)
  real(fp) :: zft(imaxt)
  !
  character(len=40) :: zlbphip,zlbphib,zunphip,zunphib,zlbtheb,zuntheb
  character(len=40) :: zlblp,zunsp
  !
  data zlbphip/'precession phi(p)/rho'/
  data zlbphib/'orbit toroidal arclength phi(b)'/
  data zlbtheb/'orbit poloidal arclength'/
  data zunphip/'radians/cm'/
  data zunphib/'radians'/
  data zuntheb/'radians'/
  !
  logical :: lxout,lyes,lchar
  !
  logical :: ixfirst
  !
  data zpi/3.141593_fp/
  !---------------------
  inx=nphipc
  do ix=1,inx
     zxi(ix)=real(ix,fp)/real(inx,fp)
  end do
  !
  imaxth=min(imaxx,(imaxt/inx))
  inumth=imaxth
  zthlimlo=-zpi
  zthlimhi=zpi
  !
  iengy=0
  zrho=1.0_fp
  zlblp=zlbphip
  zunsp=zunphip
  !
10 continue
  zdel=zthlimhi-zthlimlo
  do ith=1,inumth
     zth(ith)=zthlimlo+zdel*real(ith-1,fp)/real(inumth-1,fp)
  end do
  if(lxout(0)) then
     write(nonlin,1001) inumth,zthlimlo,zthlimhi,iengy
1001 format(/ &
          ' plotting choices:'/ &
          '  a -- plots vs. x,theta(bounce)'/ &
          '  b -- plots vs. theta,x'/ &
          '  c -- reset theta bounds and/or #pts, currently:  ',i4/ &
          '       theta bounds:  ',1pe11.4,' to ',1pe11.4/ &
          '  d -- normalize phi(p) plot to fixed d+ ion energy (ev):  ',i8/ &
          '       (do not normalize if energy is set to zero).'/ &
          '  q -- quit')
  end if
  !
  call uread('chkphipi:  enter opcode:$')
  if(lchar('q')) then
     return
  else if(lchar('a')) then
     ixfirst=.true.
  else if(lchar('b')) then
     ixfirst=.false.
  else if(lchar('d')) then
     write(nonlin,1002) iengy
1002 format(/' current d+ energy = ',i5,' ev'/ &
          ' enter an energy btw 0 and 10,000,000 ev, integer ::.'/ &
          ' 0 for phi(p)/rho plots, .gt. 0 for phi(p) plot'/ &
          ' for the selected energy.')
     call ureadi('enter integer :: d+ energy, ev:$',iengy)
     iengy=max(0,min(10000000,iengy))
     if(iengy.eq.0) then
        zlblp=zlbphip
        zunsp=zunphip
     else
        write(zlblp, &
             '('' precesion phi(p), d+ '',i8,'' ev'')') iengy
        zunsp=zunphib
     end if
     if(iengy.gt.0) then
        zev=real(iengy,fp)
        zv2=zev*(1.384e6_fp)**2/2.0_fp  ! d: a=2.0, z=1.0
        zv=sqrt(zv2)
        zrho=zv*2.0_fp/(9.578e7_fp*1.0_fp*bzxr)  ! rho is rlar(axis)/rmajor
     else
        zrho=1.0_fp
     end if
     goto 10
  else if(lchar('c')) then
     if(lxout(0)) write(nonlin,'('' limit: '',i4)') imaxth
     call uread('no. of theta pts (btw 10 and limit):$')
     inumth=idcod(iert)
     if(iert.ne.0) goto 10
     inumth=max(10,min(imaxx,inumth))
     !
20   continue
     call uread('enter first theta limit (btw -pi and pi):$')
     zth1=r8dcod(iert)
     if(iert.ne.0) goto 20
     zth1=max(-zpi,min(zpi,zth1))
     call uread('enter second theta limit (btw -pi and pi):$')
     zth2=r8dcod(iert)
     if(iert.ne.0) goto 20
     zth2=max(-zpi,min(zpi,zth2))
     !
     if(zth1.eq.zth2) then
        write(nonlin,'('' zth1,zth2 = '',2(1x,1pe11.4))') zth1,zth2
        call uermsg(' ? theta limits must be distinct.',' ',0)
        goto 20
     end if
     zthlimhi=max(zth1,zth2)
     zthlimlo=min(zth1,zth2)
     goto 10
  else
     call uermsg('illegal opcode',' ',0)
     goto 10
  end if
  !
  !  ok get the data for phi(b)
  !
  if(.not.ixfirst) then
     do ix=1,inx
        do ith=1,inumth
           ia=(ix-1)*inumth+ith
           call bphifun(ix,zth(ith),zphib,zpp,zthb)
           zfb(ia)=zphib
           zfp(ia)=zpp*zrho
           zft(ia)=zthb
        end do
     end do
     call r8_grf3f1(zth,zft,zxi,zwk,inumth,inx,inumth, &
          'bounce angle',zlbtheb,'xi', &
          'radians',zuntheb,'  ', &
          'chkphipi debug test plot',zzzlbl//runid,0)
     call r8_grf3f1(zth,zfb,zxi,zwk,inumth,inx,inumth, &
          'bounce angle',zlbphib,'xi', &
          'radians',zunphib,'  ', &
          'chkphipi debug test plot',zzzlbl//runid,0)
     call r8_grf3f1(zth,zfp,zxi,zwk,inumth,inx,inumth, &
          'bounce angle',zlblp,'xi', &
          'radians',zunsp,'  ', &
          'chkphipi debug test plot',zzzlbl//runid,0)
  else
     do ith=1,inumth
        do ix=1,inx
           ia=(ith-1)*inx+ix
           call bphifun(ix,zth(ith),zphib,zpp,zthb)
           zfb(ia)=zphib
           zfp(ia)=zpp*zrho
           zft(ia)=zthb
        end do
     end do
     call r8_grf3f1(zxi,zft,zth,zwk,inx,inumth,inx, &
          'xi',zlbtheb,'bounce angle', &
          '  ',zuntheb,'radians', &
          'chkphipi debug test plot',zzzlbl//runid,0)
     call r8_grf3f1(zxi,zfb,zth,zwk,inx,inumth,inx, &
          'xi',zlbphib,'bounce angle', &
          '  ',zunphib,'radians', &
          'chkphipi debug test plot',zzzlbl//runid,0)
     call r8_grf3f1(zxi,zfp,zth,zwk,inx,inumth,inx, &
          'xi',zlblp,'bounce angle', &
          '  ',zunsp,'radians', &
          'chkphipi debug test plot',zzzlbl//runid,0)
  end if
  !
  goto 10
end subroutine chkphipi
! 11sep2004 fgtok
