subroutine orbgdb
#ifdef __defunct

  use nbi_com
  use nbi_dimensions
  !
  !  routine for debugging subroutine orbg -- orbit integrator
  !
  !xxx      include 'oepcom'
  implicit none
  LOGICAL, EXTERNAL :: LCHAR,LXOUT,LYES,LYESC                !
  !
  !============
  integer :: ifmax,idsiz,luntrm,lnout,ilorbgdb,istart,inmax,izgdt
  integer :: izgerre,izgerrp,izgdxi,izgdth,izgde,izgdk,izgxd1
  integer :: izgtd1,izgxd2,izgtd2,izgzgb,izgzdgbdt,izgbibdts
  integer :: izgz2btst,izggcej,izgdxdps,izgmuxp,i,ii,idcod,iert
  integer :: idim1,icfx,j,iadr_offset,ifun,iadr
  !============
  real(fp) :: zxsave,ztsave,zvsave,zesave,zksave,zxuse,ztuse,zvuse
  real(fp) :: zkuse,zmin,zmax,zzmin,zzmax,zf,zfi,zeuse,zparm1,zparm2
  !============
  integer :: indx(2)
  character(len=10) :: ixlbl(4),ixuns(4)
  !
  real(fp) :: zdim1(idim1),zdim2(idim1),zwk(idim1)
  !
  real(fp) :: zextrem(2,4)
  !
  parameter (ifmax=40)
  !
  character(len=20) :: zflbl(0:ifmax)
  character(len=10) :: zfuns(0:ifmax)
  !
  character(len=35) :: zfxlbl
  integer :: iccfx(2),iccfx2(2)
  !
  real(fp), dimension(:), allocatable :: zdatbuf
  !-------------
  !
  data zextrem/0.0_fp,1.3_fp, &
       -3.1415926535897931_fp,3.1415926535897931_fp, &
       1.0_fp,1.0e7_fp,  -1.0_fp,1.0_fp/
  !
  data iccfx/1,19/
  data iccfx2/17,35/
  !
  data ixlbl/'xi        ', &
       'th        ', &
       'e         ', &
       'vpll/v    '/
  !
  data ixuns/'xi        ', &
       'radians   ', &
       'ev        ', &
       'vpll/v    '/
  !
  !--------------------------------------------------
  !
  write(6,*) ' ?? orbgdb: disabled:  cannot use datbuf, mc code'
  write(6,*) '    is being modularized...'
  !
  allocate(zdatbuf((1+ifmax)*idsiz + 1)); zdatbuf=0
  !
  luntrm=lnout(0)
  !
  ilorbgdb=.true.
  !
  zxsave=xion
  ztsave=th
  zvsave=vion
  zesave=vion*vion*v2toev
  zksave=xksid
  !
  write(6,1000) zxsave,ztsave,zesave,zksave
1000 format(' orbgdb call:  xion=',1pe11.4,' th=',1pe11.4, &
       ' e=',1pe11.4,' ksid=',1pe11.4)
  !
  zxuse=xion
  ztuse=th
  zvuse=vion
  zkuse=xksid
  !
  istart=0
  !
  inmax=-1
  !
  zflbl(izgdt)='time step'
  zfuns(izgdt)='seconds'
  inmax=inmax+1
  !
  zflbl(izgerre)='rel. energy error'
  zfuns(izgerre)=' '
  inmax=inmax+1
  !
  zflbl(izgerrp)='rel. momentum error'
  zfuns(izgerrp)=' '
  inmax=inmax+1
  !
  zflbl(izgdxi)='delta(xi)'
  zfuns(izgdxi)=' '
  inmax=inmax+1
  !
  zflbl(izgdth)='delta(th)'
  zfuns(izgdth)='radians'
  inmax=inmax+1
  !
  zflbl(izgde)='delta(e)'
  zfuns(izgde)='ev'
  inmax=inmax+1
  !
  zflbl(izgdk)='delta(vpll/v)'
  zfuns(izgdk)=' '
  inmax=inmax+1
  !
  zflbl(izgxd1)='d(xi)/dt'
  zfuns(izgxd1)='sec**-1'
  inmax=inmax+1
  !
  zflbl(izgtd1)='xi*d(th)/dt'
  zfuns(izgtd1)='sec**-1'
  inmax=inmax+1
  !
  zflbl(izgxd2)='d2(xi)/dt2'
  zfuns(izgxd2)='sec**-2'
  inmax=inmax+1
  !
  zflbl(izgtd2)='xi*d2(th)/dt2'
  zfuns(izgtd2)='sec**-2'
  inmax=inmax+1
  !
  zflbl(izgzgb)='grad(b) drift cf'
  zfuns(izgzgb)='g/sec'
  inmax=inmax+1
  !
  zflbl(izgzdgbdt)='d/dt(grad(b) drift cf)'
  zfuns(izgzdgbdt)='g/sec**2'
  inmax=inmax+1
  !
  zflbl(izgbibdts)='d(1/b)/dth'
  zfuns(izgbibdts)='g**-1'
  inmax=inmax+1
  !
  zflbl(izgz2btst)='d/dt(d(1/b)/dth)'
  zfuns(izgz2btst)='(g*sec)**-1'
  inmax=inmax+1
  !
  zflbl(izggcej)='zg*zc/(ze*zj)'
  zfuns(izggcej)='(???)'
  inmax=inmax+1
  !
  zflbl(izgdxdps)='d(xi)/d(psi)'
  zfuns(izgdxdps)='(g*cm)**-2'
  inmax=inmax+1
  !
  zflbl(izgmuxp)='mu+m*vpll**2/b'
  zfuns(izgmuxp)='g*(cm/s)**2/g'
  inmax=inmax+1
  !
  !--------------------------------------------------
10 continue
  !
  indx(1)=0
  indx(2)=0
  !
  if(lxout(0)) then
     write(luntrm,'('' orbgdb:  can vary two of the following:'')')
     do i=1,4
        write(luntrm,1001) i,ixlbl(i)
1001    format(3x,i1,'.  ',a)
     end do
  end if
  !
  call uread('orbgdb:  1st dimension of variation (1/2/3/4/q):$')
  if(lchar('q')) goto 10
  !
  ii=idcod(iert)
  if(iert.ne.0) goto 10
  if((ii.lt.1).or.(ii.gt.4)) then
     call uermsg('?orbgdb:  dim code not in range 1 to 4',' ',0)
     goto 10
  end if
  !
  indx(1)=ii
  !
  call uread('orbgdb:  2nd dimension of variation (1/2/3/4/q):$')
  if(lchar('q')) goto 10
  !
  ii=idcod(iert)
  if(iert.ne.0) goto 10
  if((ii.lt.1).or.(ii.gt.4)) then
     call uermsg('?orbgdb:  dim code not in range 1 to 4',' ',0)
     goto 10
  end if
  !
  if(ii.eq.indx(1)) then
     call uermsg('?orbgdb:  that dim already chosen.',' ',0)
     goto 10
  end if
  !
  indx(2)=ii
  !
20 continue
  if(lxout(0)) write(luntrm,1020) zextrem(1,indx(1)), &
       zextrem(2,indx(1))
1020 format(' -> range limit:  ',1pe11.4,' to ',1pe11.4,'.')
  call ureadr8('orbgdb:  enter min '//ixlbl(1)//':$',zmin)
  call ureadr8('orbgdb:  enter max '//ixlbl(1)//':$',zmax)
  if(zmin.eq.zmax) then
     call uermsg(' min and max values must be distinct',' ',0)
     goto 20
  end if
  zzmin=min(zmin,zmax)
  zzmax=max(zmin,zmax)
  if((zzmin.lt.zextrem(1,indx(1))).or. &
       (zzmax.gt.zextrem(2,indx(1)))) then
     call uermsg(' %orbgdb:  limit range error',' ',0)
     goto 20
  end if
  do ii=1,idim1
     zf=real(ii-1,fp)/real(idim1-1,fp)
     zfi=1.0_fp-zfi
     zdim1(ii)=zfi*zzmin+zf*zzmax
  end do
  !
30 continue
  if(lxout(0)) write(luntrm,1020) zextrem(1,indx(2)), &
       zextrem(2,indx(2))
  call ureadr8('orbgdb:  enter min '//ixlbl(2)//':$',zmin)
  call ureadr8('orbgdb:  enter max '//ixlbl(2)//':$',zmax)
  if(zmin.eq.zmax) then
     call uermsg(' min and max values must be distinct',' ',0)
     goto 30
  end if
  zzmin=min(zmin,zmax)
  zzmax=max(zmin,zmax)
  if((zzmin.lt.zextrem(1,indx(2))).or. &
       (zzmax.gt.zextrem(2,indx(2)))) then
     call uermsg(' %orbgdb:  limit range error',' ',0)
     goto 20
  end if
  do ii=1,idim1
     zf=real(ii-1,fp)/real(idim1-1,fp)
     zfi=1.0_fp-zfi
     zdim2(ii)=zfi*zzmin+zf*zzmax
  end do
  !
  zfxlbl=' '
  icfx=0
  if((indx(1).ne.1).and.(indx(2).ne.1)) then
42   continue
     if(lxout(0)) write(luntrm,7701) zxuse
7701 format(' fixed:  xi=',1pe11.4)
     call ureadr8('enter new value of xi or "," not to change:$$', &
          zxuse)
     if((zxuse.lt.zextrem(1,1).or.(zxuse.gt.zextrem(2,1)))) then
        call uermsr8('range error on xi input',zxuse,0)
        zxuse=zxsave
        goto 42
     end if
     icfx=icfx+1
     write(zfxlbl(iccfx(icfx):iccfx2(icfx)),'('' xi= '',1pe10.3)') &
          zxuse
  end if
  !
  if((indx(1).ne.2).and.(indx(2).ne.2)) then
44   continue
     if(lxout(0)) write(luntrm,7702) ztuse
7702 format(' fixed:  th=',1pe11.4)
     call ureadr8('enter new value of th or "," not to change:$$', &
          ztuse)
     if((ztuse.lt.zextrem(1,2).or.(ztuse.gt.zextrem(2,2)))) then
        call uermsr8('range error on th input',ztuse,0)
        ztuse=ztsave
        goto 44
     end if
     icfx=icfx+1
     write(zfxlbl(iccfx(icfx):iccfx2(icfx)),'('' th= '',1pe10.3)') &
          ztuse
  end if
  !
  if((indx(1).ne.3).and.(indx(2).ne.3)) then
46   continue
     zeuse=zvuse*zvuse*v2toev
     if(lxout(0)) write(luntrm,7703) zeuse
7703 format(' fixed:  e=',1pe11.4)
     call ureadr8('enter new value of e or "," not to change:$$', &
          zeuse)
     if((zeuse.lt.zextrem(1,3).or.(zeuse.gt.zextrem(2,3)))) then
        call uermsr8('range error on e input',zeuse,0)
        zeuse=zvuse*zvuse*v2toev
        goto 46
     end if
     zvuse=sqrt(zeuse/v2toev)
     icfx=icfx+1
     write(zfxlbl(iccfx(icfx):iccfx2(icfx)),'('' e= '',1pe10.3)') &
          zeuse
  end if
  !
  if((indx(1).ne.4).and.(indx(2).ne.4)) then
48   continue
     if(lxout(0)) write(luntrm,7704) zkuse
7704 format(' fixed:  vpll/v=',1pe11.4)
     call ureadr8( &
          'enter new value of vpll/v or "," not to change:$$', &
          zkuse)
     if((zkuse.lt.zextrem(1,4).or.(zkuse.gt.zextrem(2,4)))) then
        call uermsr8('range error on vpll/v input',zvuse,0)
        zvuse=zvsave
        goto 48
     end if
     icfx=icfx+1
     write(zfxlbl(iccfx(icfx):iccfx2(icfx)),'(''vpll/v='',1pe10.3)') &
          zkuse
  end if
  !
  !--------------------------------------------------
  !
  !  calculate...
  !
  do j=1,idim1
     do i=1,idim1
        !
        !  address offset -- in oepcom common:
        !
        iadr_offset = (j-1)*idim1 + (i-1)
        !
        xion=zxuse
        th=ztuse
        vion=zvuse
        xksid=zkuse
        !
        bidt(1)=0.0_fp
        bimodb(1)=0.0_fp
        !
        zparm1=zdim1(i)
        zparm2=zdim2(j)
        !
        !  first parm substitution
        !
        if(indx(1).eq.1) then
           xion=zparm1
        else if(indx(1).eq.2) then
           th=zparm1
        else if(indx(1).eq.3) then
           vion=sqrt(zparm1/v2toev)
        else
           xksid=zparm1
        end if
        !
        !  second parm substitution
        !
        if(indx(2).eq.1) then
           xion=zparm2
        else if(indx(2).eq.2) then
           th=zparm2
        else if(indx(2).eq.3) then
           vion=sqrt(zparm2/v2toev)
        else
           xksid=zparm2
        end if
        !
        vpll=vion*xksid
        vperp=vion*sqrt(1.0_fp-xksid*xksid)
        !
        !  do an orbit step
        !
        call nbigc(1)
        !
        call orbpre
        call orbg(ierr)
        call orbout
        !
     end do
  end do
  !
  !  restore state
  !
  !
  xion=zxsave
  th=ztsave
  vion=zvsave
  xksid=zksave
  vpll=vion*xksid
  vperp=vion*sqrt(1.0_fp-xksid*xksid)
  !
  bidt(1)=0.0_fp
  bimodb(1)=0.0_fp
  !
  call nbigc(1)
  !
  !------------------------
  !
50 continue
  if(lxout(0)) then
     write(luntrm,1050) ixlbl(indx(1)),zdim1(1),zdim1(idim1), &
          ixuns(indx(1)), &
          ixlbl(indx(2)),zdim2(1),zdim2(idim1),ixuns(indx(2))
1050 format(/ &
          ' orbgdb:  plot axes:'/ &
          '   ',a,' from ',1pe10.4,' to ',1pe10.4,1x,a/ &
          '   ',a,' from ',1pe10.4,' to ',1pe10.4,1x,a/ &
          ' choose index number for function plot;'/ &
          ' "c" for chkspl call'/ &
          ' "l" for a function list'/ &
          ' "x" to change plot axes'/ &
          ' "q" to exit this routine')
  end if
  call uread('orbgdb:  enter plot option:$$')
  !
  if(lchar('q')) goto 2000
  !
  if(lchar('x')) then
     goto 10
  else if(lchar('c')) then
     call chkspl
     goto 50
  else if(lchar('l')) then
     write(luntrm,1052)
1052 format(/'  code#  label                 units'/ &
          '  -----  --------------------  ----------')
     do ii=0,inmax
        write(luntrm,'(4x,i2,3x,a20,2x,a10)') ii,zflbl(ii),zfuns(ii)
     end do
     goto 50
  end if
  !
  ifun=idcod(iert)
  if(iert.ne.0) goto 50
  !
  if((ifun.lt.0).or.(ifun.gt.inmax)) then
     call uermsi('?orbgdb:  function code out of range:  ',ifun,0)
     goto 50
  end if
  !
  iadr=1+ifun*idsiz
  call r8_grf3f1(zdim1,zdatbuf(iadr),zdim2,zwk,idim1,idim1,idim1, &
       ixlbl(indx(1)),zflbl(ifun),ixlbl(indx(2)), &
       ixuns(indx(1)),zfuns(ifun),ixuns(indx(2)), &
       runid//' orbgdb plot',zfxlbl,0)
  !
  goto 50
  !
2000 continue
  !
  ilorbgdb=.false.
  return
  !
#else
  write(6,*) ' %nbdebug: orbgdb subroutine has been disabled.'
#endif
end subroutine orbgdb
! 11sep2004 fgtok
