subroutine chkdt
  use nbi_com
  use nbi_com, only: zpi=>xpi
  use nbi_dimensions
  !
  !   nov00 rga -- removed orbs option
  !
  !       increased zemax to value for dt fusion product protons
  !        also fixed bad label in dialog for zth2                (jfs 7 jul 94)
  !
  !	this routine checks the time step control of orbit errors
  !   by allowing the user to plot the timestep (and orbit errors)
  !   as a function of position for a given vion and xksid or as a
  !   function of vion and xksid for a fixed location. in addition,
  !   the user can plot the errors vs. dt**3 for a fixed position,
  !   vion, and xksid.
  !
  !	the best way to use this is to set break/return inispc
  !   after doing a debug restart of transp.
  !
  !	this routine is patterned after chkspl.fdb and works by
  !   asking the user what is wanted and then making repeated calls to
  !   orbg. (i presume no one needs to check orbc or orbs since they are
  !   so much simpler).
  !
  !	when used it is necessary to do so from the interactive debugger
  !   at a point in transp execution after the magnetics have been defined
  !   by a call to magcf. in fact, the use of abeam requires a previous
  !   call of inispc which is called from depall and orball which are called
  !   from fimain.
  !
  !   rtm, june 1990
  !
  !-----------------------------
  implicit none
  LOGICAL, EXTERNAL :: LCHAR,LXOUT,LYES,LYESC                !
  !-----------------------------
  !
  !-------------------------
  !	local declarations
  !-------------------------
  !
  !============
  integer :: ion,ixmax,iymax,iscalmax,ix,iy,lnout,jx,i,j,imid
  integer :: isteps,iexp,nlshrink
  integer :: iloc,ierr,index,ndtorb
  !============
  real(fp) :: zzdes,zdescr,zdpphi,zpphiabs,zrelpphi,zemin,zemax
  real(fp) :: zminxi,zmaxxi,zminks,zmaxks,zminth,zmaxth
  real(fp) :: zf1d,zdf1d,zd2f1d,zk1,zk2,ze1,ze2,zxifix,zthfix,zh
  real(fp) :: zg,zdgdxi,ziotb,zdpsdx,zdphdx,zr,zzemin,zzemax,zvmin
  real(fp) :: zvmax,zksmin,zksmax,zm,zomegaci,zb,zvpar,zmu,zdeterm
  real(fp) :: zcosb,zvrot,zvpphi,zxip,zvparp,zpphimax,zdescrim
  real(fp) :: odshrink,zdrdt,zdydt,zdrdt2,zdydt2,zvd,zad,zzefix
  real(fp) :: zvfix,zksfix,zxi1,zxi2,zximin,zximax,zth1,zth2,zthmin
  real(fp) :: zthmax,zdtsav1,zdtsav2,odsmax,odpowtip,odshrtip,odthmx
  real(fp) :: odpsmx,odvpmx,odxymx
  !============
  parameter (ion=1) !old vector ion number
  parameter (ixmax = 100)    !max xi values (same as iymax!)
  parameter (iymax = 100)    !max theta values (same as ixmax!)
  parameter (iscalmax = 65)  !number of points on scaling plots
  !
  real(fp) ::  zxi(ixmax)
  real(fp) ::  zv(ixmax)
  real(fp) ::  zke(ixmax)
  real(fp) ::  zavdt(ixmax)
  real(fp) ::  zavdp(ixmax)
  real(fp) ::  zavde(ixmax)
  real(fp) ::  zth(iymax)
  real(fp) ::  zks(iymax)
  real(fp) ::  ztemp(ixmax)
  !
  real(fp) ::  zf1(ixmax*iymax)	  !plotting arrays (2-d plots)
  real(fp) ::  zf2(ixmax*iymax)
  real(fp) ::  zf3(ixmax*iymax)
  real(fp) ::  zf4(ixmax*iymax)
  real(fp) ::  zf5(ixmax*iymax)
  real(fp) ::  zf6(ixmax*iymax)
  !
  real(fp) ::  zf7(ixmax*iymax)
  real(fp) ::  zf8(ixmax*iymax)
  real(fp) ::  zf9(ixmax*iymax)
  real(fp) ::  zf11(ixmax*iymax)
  real(fp) ::  zf12(ixmax*iymax)
  real(fp) ::  zf13(ixmax*iymax)
  !
  real(fp) :: zfdtcod(ixmax*iymax)
  !
  real(fp) ::  zdes(iscalmax)
  real(fp) ::  zdps(iscalmax)
  real(fp) ::  zdts(iscalmax)  !dt
  real(fp) ::  zdts3(iscalmax) !dt**3
  real(fp) ::  zdxi(iscalmax)
  real(fp) ::  zdth(iscalmax)
  real(fp) ::  zds(iscalmax)
  real(fp) ::  zds3(iscalmax)
  real(fp) ::  zfactor(iscalmax)
  !
  character(len=10) ::  iorblab
  character(len=20) ::  ifix1, ifix2, ifix3, ifix4
  !
  !
  real(fp) :: zzha,zf0a,zf1a,zf2a,zf3a  ! for stmt functions
  !
  !----------------------
  !	data statements
  !----------------------
  !
  data  ix/ixmax/	    !default # x values
  data  iy/iymax/	    !default # y values
  data  zemin/100.0_fp/  !min k.e. for an ion (ev)
  data  zemax/3.02e6_fp/  !max k.e. for an ion (ev)
  data  zminxi/1.0e-5_fp/  !this had better match those in orbg,orbs
  data  zmaxxi/1.3_fp/   !the grid isn't defined this far out
  data  zminks/-1.0_fp/  !min xksid
  data  zmaxks/ 1.0_fp/  !max xksid
  !
  !
  !============================================
  !
  !
  !	------------------
  !	0.1 initialization
  !	------------------
  !
  !	verify that inispc has been called
  !
  ndtorb=4  ! orbit timestep control algorithm selection switch;
            ! default value should be used.
  if(abeam.eq.0.0_fp) then
     write(nonlin, 1)
1    format(' chkdt error: abeam = 0',/, &
          ' ***it is likely that inispc has not finished',/, &
          '    initializing things at this point. it is',/, &
          '    best to use ctrl-y to get back to debug',/, &
          '    and to set break/return inispc before using',/, &
          '    chkdt again. otherwise you may try to set',/, &
          '    abeam and continue, but this is questionable.')
     !
     call ureadr8('chkdt: please enter abeam or ctrl-y:$$',abeam)
     !
  else
  end if
  !
  !	set defaults for energy, xksid and other limits
  !
  zk1 = -1.0_fp
  zk2 = 1.0_fp
  ze1 = 100.0_fp
  ze2 = 100000.0_fp
  zminth = -zpi
  zmaxth =  zpi
  !	zmaxxi = xiblo(nznblo+lcentr)	!last bloated zone
  !
  do
     call neupag !clear the screen
     !
     if(lxout(0)) write(lnout(0),100)
100  format(/,' chkdt: pick an option',/, &
          ' 1 - setup 2-d grid (default=100 x 100)',/, &
          ' 2 - examine things for fixed (xi,theta)',/, &
          ' 3 - examine things for fixed (v,xksid)',/, &
          ' 4 - examine error scaling given (xi,th,v,xksid)',/, &
!1 disabled          ' 5 - change zdtrob parameters (ndtorb, etc.)',/, &
          ' q - quit')
     !
     call uread('chkdt: enter 1-4, or q (quit):$')
     !
     if    (lchar('1')) then
        goto 10000
     elseif(lchar('2')) then
        goto 20000
     elseif(lchar('3')) then
        goto 30000
     elseif(lchar('4')) then
        goto 40000
!1     elseif(lchar('5')) then
!1        goto 50000
     elseif(lchar('q')) then
        goto 70000
     else
        call uermsg('chkdt: invalid option',' ',0)
        cycle
     end if
     !
     !
     !	----------------------
     !	1.0 inquire about grid
     !	----------------------
     !
     !
10000 continue
     !
     if(lxout(0)) write(lnout(0),10100) ixmax, iymax
10100 format(/,' chkdt: enter the number of (xi,theta) or', &
           ' (v,vpar/v) points',/, &
           ' for evaluating orbit steps',/, &
           ' ixmax = ',i5,' iymax = ',i5)
     !
10200 continue
     call ureadi('chkdt: enter # xi (or v) points (.le.ixmax):$',ix)
     if(ix.lt.1 .or. ix.gt.ixmax) then
        call uermsi( &
             '?chkdt: # xi (v) points must be btw 1 and ixmax',ix,0)
        goto 10200
     end if
     !
10300 continue
     call ureadi('chkdt: enter # th (or vp/v) points (.le.iymax):$', &
          iy)
     if(iy.lt.1 .or. iy.gt.iymax) then
        call uermsi( &
             '?chkdt: # th (vp/v) points must be btw 1 and ixmax',ix,0)
        goto 10300
     end if
     !
     cycle
     !
     !
     !	-----------------------------------------------
     !	2.0 compute & display things for fixed location
     !	-----------------------------------------------
     !
     !
20000 continue
     !
     call ureadr8('chkdt: enter fixed xi (1.0e-5, 1.3):$$',zxifix)
     zxifix = max(zminxi, min(zmaxxi, zxifix))
     call ureadr8('chkdt: enter fixed th (-pi, pi):$$',zthfix)
     zthfix = max(zminth, min(zmaxth, zthfix))
     !
     jx = int(zxifix/dxi(lcentr,2)) + 1
     zh = zxifix - real(jx-1,fp)*dxi(lcentr,2)
     zg     = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
     zdgdxi = zdf1d(zh,splg(2,jx),splg(3,jx),splg(4,jx))
     ziotb  = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx), &
          spliot(4,jx))
     zdpsdx = ziotb*zxifix*tflux*1.0e8_fp/zpi		!d(psi-p)/d(xi)
     zdphdx = zdf1d(zh,splphi(2,jx),splphi(3,jx), &
          splphi(4,jx))			!d(phi)/d(xi)
     !
     !  compute r,y, jacobean
     !
     call nbmomry(zxifix,zthfix,1,rmjion,yion,xbjaco)
     zr=rmjion
     !
     call ureadr8('chkdt: enter min ion energy (ev):$$',ze1)
     call ureadr8('chkdt: enter max ion energy (ev):$$',ze2)
     ze1 = max(zemin,min(zemax,ze1))
     ze2 = max(zemin,min(zemax,ze2))
     zzemin = min(ze1,ze2)
     zzemax = max(ze1,ze2)
     zzemin = max(zemin,zzemin)
     zzemax = min(zemax,zzemax)
     !
     zvmin = sqrt(2.0_fp*zzemin/(1.04414e-12_fp*abeam))
     zvmax = sqrt(2.0_fp*zzemax/(1.04414e-12_fp*abeam))
     !
     call ureadr8('chkdt: enter min xksid (-1,+1):$$',zk1)
     call ureadr8('chkdt: enter max xksid (-1,+1):$$',zk2)
     zksmin = min(zk1,zk2)
     zksmax = max(zk1,zk2)
     zksmin = max(zminks,min(zmaxks,zksmin))
     zksmax = max(zminks,min(zmaxks,zksmax))
     !
     !
     !	    ----------------------------
     !	    2.1	loop over velocity space
     !	    ----------------------------
     !
     !
21000 continue
     !
     zm = zmp*abeam
     do i = 1, ix
        zv(i) = zvmin + real(i-1,fp)*(zvmax-zvmin)/real(ix-1,fp)
        zke(i) = 0.5_fp*abeam*1.04414e-12_fp*zv(i)**2 !ev
        
     end do
     !
     do j = 1, iy
        zks(j) = zksmin + real(j-1,fp)*(zksmax-zksmin)/real(iy-1,fp)
        
     end do
     !
     do i = 1, ix
        zavdt(i) = 0.0_fp  !average time step
        zavdp(i) = 0.0_fp  !average |dpphi|
        zavde(i) = 0.0_fp  !average |de|
        zomegaci = zel*zb/(zm*zc)
        if(zdgdxi.eq.0.0_fp) then
           zf7(i) = 0.0_fp
        else
           zf7(i) = zomegaci*zdpsdx/(zv(i)*zdgdxi)  !critical xi_d
        end if
        do j = 1, iy
           !
           !   reset for each call to orbg
           !
           bixi(ion) = zxifix
           bithst(ion) = zthfix
           call bintrp
           zb = 1.0_fp/bi1ovb(ion)
           call jintrp
           !
           !   set up for orbg call
           !
           biv(ion)    = zv(i)
           biksid(ion) = zks(j)
           bivpar(ion) = zv(i)*zks(j)
           zvpar = bivpar(ion)
           bivprp(ion) = zv(i)*sqrt(1.0_fp-zks(j)**2)
           zmu = 0.5_fp*zm*bivprp(ion)**2*bi1ovb(ion)
           !
           !	--------------------------------------------------- tbt  6/26/90
           
           xion	= bixi(ion)
           th	= bithst(ion)
           xksid	= biksid(ion)
           vpll	= bivpar(ion)
           vperp = bivprp(ion)
           vion	= biv(ion)
           bimu(ion) = zmu
           bidtnv= dtnv
           
           !	----------------------------------------------------------------
           
           iloc = (i-1)*iy + j
           !
           zdeterm = (zm*zg*zel*bivpar(ion)*bi1ovb(ion)**2)*zdphdx + &
                (bimu(ion)*zm*bivpar(ion)*bi1ovb(ion))*zdgdxi  - &
                (bimu(ion)*zel/zc)*zdpsdx
           zf5(iloc) = zdeterm
           !
           zcosb = zg*bi1ovb(ion)/zr
           index = min(ledge,jx+lcentr-1)
           zvrot = zr*0.5_fp*(omegag(index)+omegag(index+1))
           zvpphi = bivpar(ion)*zcosb
           zf6(iloc) = zvpphi - zvrot
           !
           !---------------------------------
           zxip=bixi(ion)
           zvparp=bivpar(ion)
           !
           call orbg(ierr)
           !
           !---------------------------------
           zpphimax = zm*biv(ion)*zr !max angular momentum for a neutral
           !
           zf1(iloc) = bidte(ion)/bite(ion)
           zf2(iloc) = bidpph(ion)/zpphimax  !relative error
           zf3(iloc) = bidt(ion)
           !
           zavde(i) = zavde(i) + abs(zf1(iloc))
           zavdp(i) = zavdp(i) + abs(zf2(iloc))
           zavdt(i) = zavdt(i) + zf3(iloc)
           !
           zdescrim = zvpar**2 - (2.0_fp/zm)* &
                (zmu*bidb(ion) + zel*bidphi(ion))
           zf4(iloc) = max(0.0_fp,-zdescrim)  !plot minus the descriminant
           !
           !   now plot some things from zdtrob
           !
           !             change in mechanical momentum / m*v*r (approx)
           zf8(iloc) = zm*zr*(bivpar(ion)-zvparp)/zpphimax
           !
           zf9(iloc) = zzdes   !descriminant for computing zdt1, zdt2
           zf9(iloc) = max(0.0_fp,-zf9(iloc))
           zf11(iloc) = zdescr !descriminant predicted from derivatives
           zf11(iloc) = max(0.0_fp,-zf11(iloc))
           zf12(iloc) = zdpphi/zpphimax !predicted change in p-phi
           !
           zf13(iloc) = bixi(ion)-zxip  ! change in flux surface locn.
           !
           zfdtcod(iloc) = odshrink
           !
           zdrdt=xbjaco(1,1)*bixidt(ion)+xbjaco(1,2)*bithdt(ion)
           zdydt=xbjaco(2,1)*bixidt(ion)+xbjaco(2,2)*bithdt(ion)
           zdrdt2=xbjaco(1,1)*bixi2t(ion)+xbjaco(1,2)*bith2t(ion)
           zdydt2=xbjaco(2,1)*bixi2t(ion)+xbjaco(2,2)*bith2t(ion)
           !
           zvd=sqrt(zdrdt**2+zdydt**2)
           zad=sqrt(zdrdt2**2+zdydt2**2)
           !
           zf4(iloc)=zvd
           zf5(iloc)=zad
           !
           zf6(iloc)=zvd/max(1.0e-20_fp,zad)
           !
           
        end do
     end do
     !
     do i = 1, ix
        zavdt(i) = zavdt(i)/iy
        zavdp(i) = zavdp(i)/iy
        zavde(i) = zavde(i)/iy
        
     end do
     !
     !
     !	    ----------------
     !	    2.2	plot results
     !	    ----------------
     !
     !
22000 continue
     !
     write(unit=iorblab,fmt=22100)ndtorb
22100 format(' ndtorb=',i2)
     write(unit=ifix1,fmt=22110)zxifix
22110 format(' xi=',1pe12.4,' ')
     write(unit=ifix2,fmt=22120)zthfix
22120 format(' th=',1pe12.4,' ')
     !
     if(lxout(0)) write(lnout(0),22200)
22200 format(/,' chkdt: pick a function to plot',/, &
           ' 1 - de/t.e. vs. vp/v and ke',/, &
           ' 2 - dpphi/mvr vs. vp/v and ke',/, &
           ' 3 - |vdrift| vs. vp/v and ke'/ &
           ' 4 - |adrift| vs. vp/v and ke'/ &
           ' 5 - |vd|/|ad| vs. vp/v and ke'/ &
           ' 6 - dt vs. vp/v and ke',/, &
           ' 7 - <dt> vs. ke',/, &
           ' 8 - <|de|/te> vs. ke',/, &
           ' 9 - max(0,-zzdes) from zdtrob',/, &
           ' a - max(0,-zdescr) from zdtrob',/, &
           ' b - zdpphi/mvr from zdtrob',/, &
           ' c - <|dpphi|/mvr> vs. ke',/, &
           ' d - non-rotating critical xksid',/, &
           ' p - d(m*r*vp)/mvr vs. vp/v and ke'/ &
           ' x - d(xi) vs. vp/v and ke'/ &
           ' t - min(dt) indicator code vs. vp/v and ke'/ &
           ' q - quit')
     !
     !x     3		 ' 3 - max(0,-descr) vs. vp/v and ke',/,
     !x     4		 ' 4 - determinant vs. vp/v and ke',/,
     !x     5		 ' 5 - (vphi-vrot) vs. vp/v and ke',/,
     call uread('chkdt: enter 1-9, a-d or q (quit):$')
     !
     if(lchar('1')) then
        goto 22210
     elseif(lchar('2')) then
        goto 22220
     elseif(lchar('3')) then
        goto 22230
     elseif(lchar('4')) then
        goto 22240
     elseif(lchar('5')) then
        goto 22250
     elseif(lchar('6')) then
        goto 22260
     elseif(lchar('7')) then
        goto 22270
     elseif(lchar('8')) then
        goto 22280
     elseif(lchar('9')) then
        goto 22290
     elseif(lchar('a')) then
        goto 22211
     elseif(lchar('b')) then
        goto 22212
     elseif(lchar('c')) then
        goto 22213
     elseif(lchar('d')) then
        goto 22214
     else if(lchar('p')) then
        goto 22221
     else if(lchar('x')) then
        goto 22222
     else if(lchar('t')) then
        goto 22223
     elseif(lchar('q')) then
        cycle
     else
        call uermsg('chkdt: invalid option',' ',0)
        goto 22000
     end if
     !
     goto 22000
     !
22210 continue
     call r8_grf3f1(zks,zf1,zke,ztemp,iy,ix,iy, &
          'xksid','de','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'de/te vs. vp/v and ke',0)
     !
     goto 22000
     !
22220 continue
     call r8_grf3f1(zks,zf2,zke,ztemp,iy,ix,iy, &
          'xksid','dp-phi','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dp-phi/mvr vs. vp/v and ke',0)
     !
     goto 22000
     !
22230 continue
     call r8_grf3f1(zks,zf4,zke,ztemp,iy,ix,iy, &
          'xksid','|vdrift|','k.e.ion',' ','cm/sec',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'vdrift vs. vp/v and ke',0)
     !x	call grf3f1(zks,zf4,zke,ztemp,iy,ix,iy,
     !x     1		    'xksid','-descriminant','k.e.ion',' ',' ',' ',
     !x     2		    runid(1:lrunid)//'db'//iorblab//ifix1//ifix2,
     !x     3		    'max(0.0,-descr) vs. vp/v and ke',0)
     !
     goto 22000
     !
22240 continue
     call r8_grf3f1(zks,zf5,zke,ztemp,iy,ix,iy, &
          'xksid','|adrift|','k.e.ion',' ','cm/sec^2',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'adrift vs. vp/v and ke',0)
     !x	call grf3f1(zks,zf5,zke,ztemp,iy,ix,iy,
     !x     1		    'xksid','determinant','k.e.ion',' ',' ',' ',
     !x     2		    runid(1:lrunid)//'db'//iorblab//ifix1//ifix2,
     !x     3		    'determinant vs. vp/v and ke',0)
     !
     goto 22000
     !
22250 continue
     call r8_grf3f1(zks,zf6,zke,ztemp,iy,ix,iy, &
          'xksid','|vd|/|ad|','k.e.ion',' ','sec',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'vdrift/adrift vs. vp/v and ke',0)
     !x	call grf3f1(zks,zf6,zke,ztemp,iy,ix,iy,
     !x     1		    'xksid','vpphi-vrot','k.e.ion',' ',' ',' ',
     !x     2		    runid(1:lrunid)//'db'//iorblab//ifix1//ifix2,
     !x     3		    'vpphi-vrot vs. vp/v and ke',0)
     !
     goto 22000
     !
22260 continue
     call r8_grf3f1(zks,zf3,zke,ztemp,iy,ix,iy, &
          'xksid','dt','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dt vs. vp/v and ke',0)
     !
     goto 22000
     !
22270 continue
     call r8_grafx1(zke,zavdt,ix,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' <dt> vs. ke')
     !
     goto 22000
     !
22280 continue
     call r8_grafx1(zke,zavde,ix,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' <|de|/te> vs. ke')
     !
     goto 22000
     !
22290 continue
     call r8_grf3f1(zks,zf9,zke,ztemp,iy,ix,iy, &
          'xksid','zzdes','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'zzdes vs. vp/v and ke',0)
     !
     goto 22000
     !
22211 continue
     call r8_grf3f1(zks,zf11,zke,ztemp,iy,ix,iy, &
          'xksid','zdescr','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'zdescr vs. vp/v and ke',0)
     !
     goto 22000
     !
22212 continue
     call r8_grf3f1(zks,zf12,zke,ztemp,iy,ix,iy, &
          'xksid','zdpphi/mvr','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'zdpphi/mvr vs. vp/v and ke',0)
     !
     goto 22000
     !
22213 continue
     call r8_grafx1(zke,zavdp,ix,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' <|dpphi|/mvr> vs. ke')
     !
     goto 22000
     !
22214 continue
     call r8_grafx1(zke,zf7,ix,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' non-rotating critical xksid vs. ke')
     !
     goto 22000
     !
     !
22221 continue
     call r8_grf3f1(zks,zf8,zke,ztemp,iy,ix,iy, &
          'xksid','mrdvp/mvr','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'mrdvp/mvr vs. vp/v and ke',0)
     !
     goto 22000
     !
22222 continue
     call r8_grf3f1(zks,zf13,zke,ztemp,iy,ix,iy, &
          'xksid','d(xi)','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'd(xi) vs. vp/v and ke',0)
     !
     goto 22000
     !
22223 continue
     !
     !  nmodel=6 dt-code, identifying smallest timestep factor
     !   (cf zdtrob.fnb setting of odshrink:
     !      1st order               2nd order
     !         1        d(th)          5
     !         2        d(xi)          6
     !         3        d(vpll)        7
     !         4        d(xy-twiddle)  8  (near origin)  )
     !
     call r8_grf3f1(zks,zfdtcod,zke,ztemp,iy,ix,iy, &
          'xksid','dt-code','k.e.ion',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dt-code vs. vp/v and ke',0)
     !
     goto 22000
     !
     !	-----------------------------------------------
     !	3.0 compute & display things for fixed velocity
     !	-----------------------------------------------
     !
     !
30000 continue
     !
     call ureadr8('chkdt: enter fixed ion energy (ev):$$',zzefix)
     zzefix = max(zemin, min(zemax, zzefix))
     zvfix = sqrt(2.0_fp*zzefix/(1.04414e-12_fp*abeam))
     call ureadr8('chkdt: enter fixed xksid:$$',zksfix)
     zksfix = max(zminks, min(zmaxks, zksfix))
     !
     call ureadr8('chkdt: enter min xi (1.0e-5, 1.3):$$',zxi1)
     call ureadr8('chkdt: enter max xi (1.0e-5, 1.3):$$',zxi2)
     zximin = min(zxi1,zxi2)
     zximax = max(zxi1,zxi2)
     zximin = max(zminxi, min(zmaxxi, zximin))
     zximax = max(zminxi, min(zmaxxi, zximax))
     !
     call ureadr8('chkdt: enter min theta (-pi,+pi):$$',zth1)
     call ureadr8('chkdt: enter max theta (-pi,+pi):$$',zth2)
     zthmin = min(zth1,zth2)
     zthmax = max(zth1,zth2)
     zthmin = max(zminth,min(zmaxth,zthmin))
     zthmax = max(zminth,min(zmaxth,zthmax))
     !
     !
     !	    ---------------------------
     !	    3.1	loop over (xi,th) space
     !	    ---------------------------
     !
     !
31000 continue
     !
     do i = 1, ix
        zxi(i) = zximin + real(i-1,fp)*(zximax-zximin)/real(ix-1,fp)
        
     end do
     !
     do j = 1, iy
        zth(j) = zthmin + real(j-1,fp)*(zthmax-zthmin)/real(iy-1,fp)
        
     end do
     !
     do i = 1, ix
        do j = 1, iy
           !
           !   set up for orbg call
           !
           biv(ion) = zvfix
           biksid(ion) = zksfix
           bivpar(ion) = zvfix*zksfix
           bivprp(ion) = zvfix*sqrt(1.0_fp-zksfix**2)
           bixi(ion) = zxi(i)
           bithst(ion) = zth(j)
           call bintrp
           call jintrp
           !
           jx = int(bixi(ion)/dxi(lcentr,2)) + 1
           zh = bixi(ion) - real(jx-1,fp)*dxi(lcentr,2)
           !
           !  compute r,y, jacobean
           !
           call nbmomry(bixi(ion),bithst(ion),1,rmjion,yion,xbjaco)
           zr=rmjion
           !
           call orbg(ierr)
           !
           iloc = (i-1)*iy + j
           zm = zmp*abeam
           zpphimax = zm*biv(ion)*zr !max angular momentum for a neutral
           !
           zf1(iloc) = bidte(ion)/bite(ion)
           zf2(iloc) = bidpph(ion)/zpphimax  !relative error
           zf3(iloc) = bidt(ion)
           !
           zdrdt=xbjaco(1,1)*bixidt(ion)+xbjaco(1,2)*bithdt(ion)
           zdydt=xbjaco(2,1)*bixidt(ion)+xbjaco(2,2)*bithdt(ion)
           zdrdt2=xbjaco(1,1)*bixi2t(ion)+xbjaco(1,2)*bith2t(ion)
           zdydt2=xbjaco(2,1)*bixi2t(ion)+xbjaco(2,2)*bith2t(ion)
           !
           zvd=sqrt(zdrdt**2+zdydt**2)
           zad=sqrt(zdrdt2**2+zdydt2**2)
           !
           zf4(iloc)=zvd
           zf5(iloc)=zad
           !
           zf6(iloc)=zvd/max(1.0e-20_fp,zad)
           !
           zfdtcod(iloc) = odshrink
           
        end do
     end do
     !
     !
     !	    ----------------
     !	    3.2	plot results
     !	    ----------------
     !
     !
32000 continue
     !
     write(unit=iorblab,fmt=32100) ndtorb
32100 format(' ndtorb=',i2)
     write(unit=ifix1,fmt=32110) zzefix
32110 format(' eion=',1pe12.4,' ')
     write(unit=ifix2,fmt=32120) zksfix
32120 format(' vp/v=',1pe12.4,' ')
     !
     if(lxout(0)) write(lnout(0),32200)
32200 format(/,' chkdt: pick a function to plot',/, &
           ' 1 - de/t.e. vs. xi and theta',/, &
           ' 2 - dpphi/mvr vs. xi and theta',/, &
           ' 3 - dt vs. xi and theta',/, &
           ' 4 - |vdrift| vs. xi and theta'/ &
           ' 5 - |adrift| vs. xi and theta'/ &
           ' 6 - |vd|/|ad| vs. xi and theta'/ &
           ' t - min(dt) indicator code vs. xi and theta'/ &
           ' q - quit')
     !
     call uread('chkdt: enter 1-3 or q (quit):$')
     !
     if(lchar('1')) then
        goto 32210
     elseif(lchar('2')) then
        goto 32220
     elseif(lchar('3')) then
        goto 32230
     elseif(lchar('4')) then
        goto 32240
     elseif(lchar('5')) then
        goto 32250
     elseif(lchar('6')) then
        goto 32260
     else if(lchar('t')) then
        goto 32440
     elseif(lchar('q')) then
       cycle
     else
        call uermsg('chkdt: invalid option',' ',0)
        goto 32000
     end if
     !
     goto 32000
     !
32210 continue
     call r8_grf3f1(zth,zf1,zxi,ztemp,iy,ix,iy, &
          'theta','de/t.e.','xi',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'de/te vs. theta and xi',0)
     !
     goto 32000
     !
32220 continue
     call r8_grf3f1(zth,zf2,zxi,ztemp,iy,ix,iy, &
          'theta','dp-phi/mvr','xi',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dp-phi/mvr vs. theta and xi',0)
     !
     goto 32000
     !
32230 continue
     call r8_grf3f1(zth,zf3,zxi,ztemp,iy,ix,iy, &
          'theta','dt','xi',' ','sec',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dt vs. theta and xi',0)
     !
     goto 32000
     !
32240 continue
     call r8_grf3f1(zth,zf4,zxi,ztemp,iy,ix,iy, &
          'theta','|vdrift|','xi',' ','cm/sec',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'vdrift vs. theta and xi',0)
     !
     goto 32000
     !
32250 continue
     call r8_grf3f1(zth,zf5,zxi,ztemp,iy,ix,iy, &
          'theta','|adrift|','xi',' ','cm/sec^2',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'adrift vs. theta and xi',0)
     !
     goto 32000
     !
32260 continue
     call r8_grf3f1(zth,zf6,zxi,ztemp,iy,ix,iy, &
          'theta','|vd|/|ad|','xi',' ','sec',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'vdrift/adrift vs. theta and xi',0)
     !
     goto 32000
     !
32440 continue
     !  see comments ca line 22223 above (dmc 21 jul 94)
     call r8_grf3f1(zth,zfdtcod,zxi,ztemp,iy,ix,iy, &
          'theta','dt-code','xi',' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2, &
          'dt-code vs. theta and xi',0)
     !
     goto 32000
     !
     !
     !	--------------------------------------------------------------
     !	4.0 compute & display error scaling for a given ion & location
     !	--------------------------------------------------------------
     !
     !
40000 continue
     !
     call ureadr8('chkdt: enter fixed ion energy (ev):$$',zzefix)
     zzefix = max(zemin, min(zemax, zzefix))
     zvfix = sqrt(2.0_fp*zzefix/(1.04414e-12_fp*abeam))
     !
     call ureadr8('chkdt:enter fixed ion velocity (cm/sec):$$',zvfix)
     !
     call ureadr8('chkdt: enter fixed xksid:$$',zksfix)
     zksfix = max(zminks, min(zmaxks, zksfix))
     !
     call ureadr8('chkdt: enter fixed xi (1.0e-5, 1.3):$$',zxifix)
     zxifix = max(zminxi, min(zmaxxi, zxifix))
     !
     call ureadr8('chkdt: enter fixed th value:$$',zthfix)
     zthfix = max(zminth, min(zmaxth, zthfix))
     !
     jx = int(zxifix/dxi(lcentr,2)) + 1
     zh = zxifix - real(jx-1,fp)*dxi(lcentr,2)
     !
     !  compute r,y, jacobean
     !
     call nbmomry(zxifix,zthfix,1,rmjion,yion,xbjaco)
     zr=rmjion
     !
     !
     !	    -----------------------
     !	    4.1	loop over step size
     !	    -----------------------
     !
     !
41000 continue
     !
     zdtsav1 = bidtnv	  !save the current scalar multiplier
     zdtsav2 = odsmax	  !currently these are the two to remember
     !                                 !see zdtrob.fnb
     imid = (iscalmax+1)/2	!=33
     isteps = 2*imid-1   !multiply by 1.1**-n, 1 and 1.1**n, n=1,2,...32
     !			    !that is 1.1**-32 = 1/21.1 and 1.1**32 = 21.1
     do i = 1, isteps
        iexp = (i-imid)
        zfactor(i) = 1.1_fp**iexp
        
     end do
     !
     do i = 1, isteps
        bidtnv = zdtsav1*zfactor(i)	!modify the dt scale factor
        odsmax = zdtsav2*zfactor(i)	!for all methods (see zdtrob)
        !
        !	  ***set up for orbg call***
        !
        bixi(ion)   = zxifix
        bithst(ion) = zthfix
        biv(ion)    = zvfix
        biksid(ion) = zksfix
        bivpar(ion) = zksfix*zvfix
        bivprp(ion) = zvfix*sqrt(1.0_fp-zksfix**2)
        call bintrp
        call jintrp
        !
        call orbg(ierr)               !compute a single step
        !
        zdts(i) = bidt(ion)
        zdts3(i) = zdts(i)**3	!errors should scale with dt**3
        zdes(i) = bidte(ion)
        zdps(i) = bidpph(ion)
        zdxi(i) = bidxi(ion)
        zdth(i) = bidts(ion)
        zds(i) = sqrt(zdxi(i)**2 + (zxifix*zdth(i))**2) !arc length step size
        zds3(i) = zds(i)**3
        
     end do
     !
     !	put things back the way they were
     !
     bidtnv = zdtsav1
     odsmax = zdtsav2
     !
     !
     !	    ----------------
     !	    4.2	plot results
     !	    ----------------
     !
     !
42000 continue
     !
     write(unit=iorblab,fmt=42100) ndtorb
42100 format(' ndtorb=',i2)
     write(unit=ifix1,fmt=42110) zzefix
42110 format(' eion=',1pe12.4,' ')
     write(unit=ifix2,fmt=42120) zksfix
42120 format(' vp/v=',1pe12.4,' ')
     write(unit=ifix3,fmt=42130) zxifix
42130 format(' xi=',1pe12.4,' ')
     write(unit=ifix4,fmt=42140) zthfix
42140 format(' th=',1pe12.4,' ')
     !
     if(lxout(0)) write(lnout(0),42200)
42200 format(/,' chkdt: pick a function to plot',/, &
           ' 1 - de vs. ds**3',/, &
           ' 2 - dpphi vs. ds**3',/, &
           ' 3 - ds vs. dt',/, &
           ' 4 - de vs. dt**3',/, &
           ' 5 - dpphi vs. dt**3',/, &
           ' q - quit')
     !
     call uread('chkdt: enter 1-5 or q (quit):$')
     !
     if(lchar('1')) then
        goto 42210
     elseif(lchar('2')) then
        goto 42220
     elseif(lchar('3')) then
        goto 42230
     elseif(lchar('4')) then
        goto 42240
     elseif(lchar('5')) then
        goto 42250
     elseif(lchar('q')) then
       cycle
     else
        call uermsg('chkdt: invalid option',' ',0)
        goto 42000
     end if
     !
     goto 42000
     !
42210 continue
     call r8_grafx1(zds3,zdes,isteps,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' energy error vs. ds**3')
     !
     goto 42000
     !
42220 continue
     call r8_grafx1(zds3,zdps,isteps,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' momentum error vs. ds**3')
     !
     goto 42000
     !
42230 continue
     call r8_grafx1(zdts,zds,isteps,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' arc length (ds) vs. dt')
     !
     goto 42000
     !
42240 continue
     call r8_grafx1(zdts3,zdes,isteps,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' de vs. dt**3')
     !
     goto 42000
     !
42250 continue
     call r8_grafx1(zdts3,zdps,isteps,' ',' ',' ', &
          runid(1:lrunid)//'db'//iorblab//ifix1//ifix2//ifix3//ifix4, &
          ' dpphi vs. dt**3')
     !
     goto 42000
     !
     !
     !	----------------------------------
     !	5.0 change ndtorb, dtnv, or odsmax
     !	----------------------------------
     !
     !
50000 continue
     !
50100 continue
     call ureadi('chkdt: enter ndtorb (1-6):$',ndtorb)
     if(ndtorb.lt.1 .or. ndtorb.gt.6) then
        call uermsi( &
             '?chkdt: ndtorb must be btw 1 & 6',ndtorb,0)
        goto 50100
     end if
     !
     if    (ndtorb.eq.1 .or. ndtorb.eq.4 .or. ndtorb.eq.5) then
        !
        call ureadr8('chkdt: enter dtnv (models 1&4):$',dtnv)
        !
     elseif(ndtorb.eq.2 .or. ndtorb.eq.3) then
        !
        call ureadr8('chkdt: enter odsmax   (models 2&3):$',odsmax)
        call ureadr8('chkdt: enter odpowtip (models 2&3):$',odpowtip)
        call ureadr8('chkdt: enter odshrtip (models 2&3):$',odshrtip)
        call ureadr8('chkdt: enter odshrink (models 2&3):$',odshrink)
        !
     else if(ndtorb.eq.6) then
        call ureadr8('chkdt:  enter odthmx (ndtorb=6):$',odthmx)
        call ureadr8('chkdt:  enter odpsmx (ndtorb=6):$',odpsmx)
        call ureadr8('chkdt:  enter odvpmx (ndtorb=6):$',odvpmx)
        call ureadr8('chkdt:  enter odxymx (ndtorb=6):$',odxymx)
     end if
     !
     if(ndtorb.eq.3 .or. ndtorb.eq.5) then
        call ureadr8('chkdt: enter zrelpphi (dpphi/(mvr)):$',zrelpphi)
        call ureadi('chkdt: enter nlshrink (models 3&5):$',nlshrink)
     else
     end if
  enddo
  !
  !
  !
  !	--------------------
  !	7.0 cleanup & return
  !	--------------------
  !
  !
70000 continue	!come here for a normal return
  !
  !***
  !
79000 continue	!come here after handling a bad non-fatal error
  !
  return
end subroutine chkdt
! 11sep2004 fgtok
