!  debug plots of mod(grad(xi)), mod(grad(theta)), mod(grad(b)),
!                 mod(grad(philb)),
!                   xx = xihat dot ( phihat ^ lhat)
!                     where xihat is unit vector in grad(xi) direction
!                          phihat is unit vector in toroidal direction
!                            lhat is unit vector in grad(b) direction
!
!                   philb' = d/dxi ( philb ) | fixed b
!                      = grad(philb) dot (phihat ^ lhat)
!                        -------------------------------
!                        grad(xi)    dot (phihat ^ lhat)
!
subroutine chkspd
  use nbi_com
  use nbi_com, only: zpi=>xpi
  use nbi_dimensions
  !
  !============
  implicit none
  integer :: inxi,inth,io,imax
  integer :: ix,it,ians,ith,ixi,izon,ioff,ia,ia2
  !============
  real(fp) :: zxi1,zxi2,zth1,zth2,zximin,zximax,zthmin,zthmax
  real(fp) :: zth,zxi,zh,zr,zy,zdxir,zdxiy,zans,zdthr,zdthy,zbgauss
  real(fp) :: zdbdxi,zdbdth,zdbr,zdby,zgbext,zph,zdphdxi,zdphdth
  real(fp) :: zdphdr,zdphdy,znorm,zxhatr,zxhaty,zgradb,zmr,zmy
  real(fp) :: zdenom,znum
  !============
  logical :: lchar,lxout
  !
  save zxi1,zxi2,zth1,zth2
  !
  character(len=20) :: zxilab,zthlab,zxiuns,zthuns,zbuns
  !
  real(fp), dimension(:), allocatable :: zxay,zthay
  real(fp), dimension(:,:), allocatable :: zwkay,zf1ay,zf2ay
  !
  data inxi/0/
  data inth/0/
  !
  data zxilab/'flux surf lbl (xi)'/
  data zthlab/'transp theta (th)'/
  data zxiuns/'xi'/
  data zthuns/'theta'/
   !
  !-------------------------------
  !
  io=nonlin
  !
10 continue
  !
  write(io,'(/'' plotting:  miscellaneous quantities:  '')')
  !
  call ureadi('enter inxi= no. of xi points (0 to quit)',inxi)
  if(inxi.eq.0) return
  inxi=max(3,min(1001,inxi))
  call ureadi('enter inth= no. of th points (0 to quit)',inth)
  if(inth.eq.0) return
  inth=max(3,min(1001,inth))

20 continue
  call ureadr8('enter xi(min) (btw 0 and 1):$$',zxi1)
  call ureadr8('enter xi(max) (btw 0 and 1):$$',zxi2)
  zximin=max(0.0_fp,min(1.0_fp,min(zxi1,zxi2)))
  zximax=max(0.0_fp,min(1.0_fp,max(zxi1,zxi2)))

  if(zximax.le.zximin) then
     call uermsg(' % xi range specification error',' ',0)
     goto 20
  end if

30 continue
  call ureadr8('enter th(min) (btw -pi and pi):$$',zth1)
  call ureadr8('enter th(max) (btw -pi and pi):$$',zth2)
  zthmin=max(-zpi,min(zpi,min(zth1,zth2)))
  zthmax=max(-zpi,min(zpi,max(zth1,zth2)))

  if(zthmax.le.zthmin) then
     call uermsg(' % th range specification error',' ',0)
     goto 30
  end if
  !
  !----------------------
  !
  imax=max(inxi,inth)
  allocate(zxay(inxi),zthay(inth))
  allocate(zwkay(imax,2),zf1ay(inxi,inth),zf2ay(inxi,inth))
  zwkay=0.0
  !
  do ix=1,inxi
     zxay(ix)=zximin+ &
          real(ix-1,fp)*(zximax-zximin)/real(inxi-1,fp)
  end do
  !
  do it=1,inth
     zthay(it)=zthmin+ &
          real(it-1,fp)*(zthmax-zthmin)/real(inth-1,fp)
  end do
  !
  !  interactive menu
  !
100 continue
  if(lxout(0)) then
     write(io,1100)
1100 format(/ &
          ' chkspd options:'/ &
          ' r -- reset no. of pts and/or (xi,theta) plotting range'/ &
          ' q -- quit'/ &
          ' 1 -- mod(grad(xi))'/ &
          ' 2 -- xi*mod(grad(theta))'/ &
          ' 3 -- mod(grad(b))'/ &
          ' 4 -- xi*mod(grad(philb))'/ &
          '      philb = toroidal arc length of field line'/ &
          ' 5 -- component of xihat along field contour'/ &
          '      xihat = unit vector perp to flux surface'/ &
          ' 6 -- philb#:  d(philb)/dxi along field contour')
     !
     if(lspthets.le.0) then
        write(io,*) ' '
        write(io,*) ' %warning: no boozer data, philb=0 assumed.'
        write(io,*) ' '
     end if
     !
  end if
  !
  call uread('chkspd:  choose option:$$')
  !
  ians=0
  if(lchar('q')) then
     deallocate(zxay,zthay,zwkay,zf1ay,zf2ay)
     return
  else if(lchar('r')) then
     deallocate(zxay,zthay,zwkay,zf1ay,zf2ay)
     goto 10
  else if(lchar('1')) then
     ians=1
  else if(lchar('2')) then
     ians=2
  else if(lchar('3')) then
     ians=3
  else if(lchar('4')) then
     ians=4
  else if(lchar('5')) then
     ians=5
  else if(lchar('6')) then
     ians=6
  else
     call uermsg('?chkspd:  illegal option code.',' ',0)
  end if
  !
  if(ians.eq.0) goto 100
  !
  !  compute stuff
  !
  do ith=1,inth
     zth=zthay(ith)
     do ixi=1,inxi
        zxi=zxay(ixi)
        zxi=min(0.999999_fp,max(1.e-5_fp,zxi))
        !
        izon=min(ledge, (lcentr + int(nzones*zxi)))
        if(zxi.lt.xi(izon,1)) izon=izon-1
        if(zxi.gt.xi(izon+1,1)) izon=izon+1
        ix=izon-lcentr+1
        zh=zxi-xi(izon,1)
        !
        !  compute the r-y location and jacobian (xbjaco)
        !
        call nbmomry(zxi,zth,1,zr,zy,xbjaco)
        !
        !  grad(xi)
        call rygrad(1.0_fp,0.0_fp,zdxir,zdxiy)
        if(ians.eq.1) then
           zans=sqrt(zdxir**2+zdxiy**2)
           goto 150
        end if
        !  xi*grad(th)
        call rygrad(0.0_fp,1.0_fp,zdthr,zdthy)
        if(ians.eq.2) then
           zans=zxi*sqrt(zdthr**2+zdthy**2)
           goto 150
        end if
        !  grad(b)
        bixi(1)=zxi
        bithst(1)=zth
        call bintrp
        zbgauss=1.0_fp/bi1ovb(1)
        zdbdxi=-1.0e-4_fp*zbgauss*zbgauss*bibdxi(1)
        zdbdth=-1.0e-4_fp*zbgauss*zbgauss*bibdts(1)
        call rygrad(zdbdxi,zdbdth,zdbr,zdby)
        if(ians.eq.3) then
           zans=sqrt(zdbr**2+zdby**2)
           zgbext=bzxr/(zr*zr)
           goto 150
        end if
        !  grad(philb)
        if(lspthets.gt.0) then
           call philb(zxi,zth,zph,zdphdxi,zdphdth)
        else
           zph=0.0_fp
           zdphdxi=0.0_fp
           zdphdth=0.0_fp
        end if
        call rygrad(zdphdxi,zdphdth,zdphdr,zdphdy)
        if(ians.eq.4) then
           zans=zxi*sqrt(zdphdr**2+zdphdy**2)
           goto 150
        end if
        !  xihat dot (unit vector in direction of field contour)
        !
        !  first the xihat vector, form from grad(xi) vector
        znorm=sqrt(zdxir**2+zdxiy**2)
        zxhatr=zdxir/znorm
        zxhaty=zdxiy/znorm
        !  next use the grad(b) vector:
        zgradb=sqrt(zdbr**2+zdby**2)
        if(zgradb.eq.0.0_fp) then
           zdbr=-1.0_fp
           zdby=0.0_fp
           zgradb=1.0_fp
        end if
        !  cross with phihat and normalize:
        zmr=zdby/zgradb
        zmy=-zdbr/zgradb
        !  eval the dot product
        if(ians.eq.5) then
           zans=zxhatr*zmr+zxhaty*zmy
           goto 150
        end if
        !  d(philb)/dxi along const(b) contour (perp to grad(b))
        !  m vector is unit vector in const(b) contour direction
        !
        zdenom=zdxir*zmr+zdxiy*zmy
        znum=zdphdr*zmr+zdphdy*zmy
        !
        if(zdenom.eq.0.0_fp) then
           if(znum.eq.0.0_fp) then
              zans=1.0e-8_fp
           else if(znum.gt.0.0_fp) then
              zans=1.0e8_fp
           else if(znum.lt.0.0_fp) then
              zans=-1.0e8_fp
           end if
        else
           zans=znum/zdenom
        end if
        !
        !--------------------------
        !
150     continue
        zf1ay(ixi,ith)=zans
        if(ians.eq.3) zf2ay(ixi,ith)=zgbext
        !
     end do
  end do
  !
  !  plots
  !
  if(ians.eq.1) then
     call r8_grf3f1( &
          zxay,zf1ay,zthay,zwkay, &
          inxi,inth,inxi,zxilab,'mod(grad(xi))',zthlab, &
          zxiuns,'1/cm',zthuns,'chkspd debug plot',runid,0)
  else if(ians.eq.2) then
     call r8_grf3f1( &
          zxay,zf1ay,zthay,zwkay, &
          inxi,inth,inxi,zxilab,'xi*mod(grad(th))',zthlab, &
          zxiuns,'radians/cm',zthuns,'chkspd debug plot',runid,0)
  else if(ians.eq.3) then
     call r8_grf3f2(zxay,zf1ay,zf2ay, &
          zthay,zwkay,inxi,inth,inxi, &
          zxilab,'mod(grad(b))','mod(grad(bext))',zthlab, &
          zxiuns,'t/cm',zthuns,'chkspd debug plot',runid,0)
  else if(ians.eq.4) then
     call r8_grf3f1( &
          zxay,zf1ay,zthay,zwkay, &
          inxi,inth,inxi,zxilab,'xi*mod(grad(philb))',zthlab, &
          zxiuns,'radians/cm',zthuns,'chkspd debug plot',runid,0)
  else if(ians.eq.5) then
     call r8_grf3f1( &
          zxay,zf1ay,zthay,zwkay, &
          inxi,inth,inxi,zxilab,'xihat.(phihat^gbhat)',zthlab, &
          zxiuns,' ',zthuns, &
          'chkspd: component of perp unit vector along const(b)', &
          runid,0)
  else if(ians.eq.6) then
     call r8_grf3f1( &
          zxay,zf1ay,zthay,zwkay, &
          inxi,inth,inxi,zxilab,'d(philb)/d(xi)|b',zthlab, &
          zxiuns,'radians',zthuns,'chkspd debug plot',runid,0)
  end if
  !
  goto 100
  !
end subroutine chkspd
! 11sep2004 fgtok
