!  debug plots of interpolated b field; comparisons of representations.
!
subroutine chkspb
  use nbi_com
  use nbi_dimensions
  !
  !============
  implicit none
  integer :: inxi,inth,io,imax
  integer :: ix,ians,ilinx,ilinth,ith,ixi,izon,ier
  !============
  real(fp) :: zxi1,zxi2,zth1,zth2,zpi,zximin,zximax,zthmin,zthmax
  real(fp) :: zth,zxi,zh,zzxi,zrpos,zypos,zbfs,zbpfs,zbtfs,zbinb
  real(fp) :: zdldth,ziot,zq,zg,zi,zval,zdthsdth,zdnudth,zbt,znudl
  real(fp) :: zthsdl,zbp,zb
  !============
  real(fp) :: zjac(2,2)
  real(fp), dimension(:), allocatable :: zxay,zthay
  real(fp), dimension(:,:), allocatable :: zb1ay,zb2ay,zwkay
  !
  integer :: isel3(6)
  !
  logical :: lchar,lxout
  !
  save zxi1,zxi2,zth1,zth2
  !
  character(len=20) :: zxilab,zthlab,zxiuns,zthuns,zbuns
  !
  data isel3/0,0,1,0,0,0/
  !
  data inxi/0/
  data inth/0/
  !
  data zxilab/'flux surf lbl (xi)'/
  data zthlab/'transp theta (th)'/
  data zxiuns/'xi'/
  data zthuns/'theta'/
  data zbuns/'tesla'/
  !
  data zpi/3.141593_fp/
  !
  !-------------------------------
  !
  !      call mombs_ngint  ! set ngint
  !
  io=nonlin
  !
10 continue
  !
  write(io,'(/'' plotting:  b field comparisons:  '')')
  !
  call ureadi('enter inxi= no. of xi points (0 to quit)',inxi)
  if(inxi.eq.0) goto 1000
  inxi=max(3,min(1001,inxi))
  call ureadi('enter inth= no. of th points (0 to quit)',inth)
  if(inth.eq.0) goto 1000
  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),zb1ay(inxi,inth),zb2ay(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 ith=1,inth
     zthay(ith)=zthmin+ &
          real(ith-1,fp)*(zthmax-zthmin)/real(inth-1,fp)
  end do
  !
  !  interactive menu
  !
100 continue
  if(lxout(0)) then
     write(io,1100)
1100 format(/ &
          ' chkspb options:'/ &
          ' r -- reset no. of pts and/or (xi,theta) plotting range'/ &
          ' q -- quit'/ &
          ' 1 -- bpol:  boozer representation vs. fourier-spline formula'/ &
          ' 2 -- mod(b):  boozer representation vs. 1/(1/b) fourier-spline'/ &
          ' 3 -- mod(b):  boozer rep. vs. mombs fourier-spline formula'/ &
          ' 4 -- mod(b):  1/(1/b) f-s vs. mombs f-s formula')
     !
     if(lspthets.eq.0) then
        write(io,*) ' '
        write(io,*) &
             ' %nb -- no boozer data; options 1-3 unavailable!'
        write(io,*) ' '
     end if
     !
  end if
  !
  call uread('chkspb:  choose option:$$')
  !
  ians=0
  if(lchar('q')) then
     deallocate(zxay,zthay,zwkay,zb1ay,zb2ay)
     goto 1000
  else if(lchar('r')) then
     deallocate(zxay,zthay,zwkay,zb1ay,zb2ay)
     goto 10
  end if

  if(lspthets.ne.0) then
     if(lchar('1')) then
        ians=1
     else if(lchar('2')) then
        ians=2
     else if(lchar('3')) then
        ians=3
     end if
  end if
  if(lchar('4')) then
     ians=4
  end if
  !
  if(ians.eq.0) then
     call uermsg('?chkspb:  illegal option code.',' ',0)
     goto 100
  end if
  !
  !  compute b field representations
  !
  ilinx=1
  ilinth=1
  !
  do ith=1,inth
     zth=zthay(ith)
     do ixi=1,inxi
        zxi=zxay(ixi)
        !
        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 f-s b and components; this also fills 2x2 xbjaco
        !  with (dr/dxi, dr/dth / dy/dxi, dy/dth )
        !
        zzxi = max(.001_fp,zxi)  ! avoid singularity
        call mombs_nb(zzxi,zth,zrpos,zypos,zjac,zbfs,zbpfs,zbtfs)
        !
        if((ians.eq.2).or.(ians.eq.4)) then
           bixi(1)=zxi
           bithst(1)=zth
           call bintrp
           zbinb=1.0e-4_fp/bi1ovb(1)
        end if
        !
        if(ians.le.3) then
           zdldth=sqrt(zjac(1,2)**2+zjac(2,2)**2) ! dlp/dtheta
           ziot=spliot(1,ix)+zh*(spliot(2,ix) + zh*( &
                spliot(3,ix)+zh*spliot(4,ix)))
           zq=1.0_fp/ziot
           zg=1.0e-4_fp*(splg(1,ix)+zh*(splg(2,ix) + zh*( &
                splg(3,ix)+zh*splg(4,ix))))
           zi=1.0e-4_fp*(splimhd(1,ix)+zh*(splimhd(2,ix) + zh*( &
                splimhd(3,ix)+zh*splimhd(4,ix))))
           !
           call bcspeval(zxi,zth,isel3,zval, &
                xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
                ilinx,ilinth, &
                datspl(lspthets),nzp1,ier)
           zdthsdth=zval            ! dtheta*/dtheta
           !
           call bcspeval(zxi,zth,isel3,zval, &
                xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
                ilinx,ilinth, &
                datspl(lspphilb),nzp1,ier)
           zdnudth=zval             ! dnu/dtheta, zeta = phi + nu
           !
           !  b(vector) = g*grad(zeta) + i*grad(theta*) + correction to keep
           !                                              field in flux surface
           !
           zbt = zg/zrpos     ! toroidal component
           !
           znudl = zdnudth/zdldth
           zthsdl = zdthsdth/zdldth
           zbp = zi*zthsdl + zg*znudl  ! poloidal component
           !
           zb = sqrt(zbt*zbt + zbp*zbp)
           !
        end if
        !
        if(ians.eq.1) then
           zb1ay(ixi,ith)=zbpfs
           zb2ay(ixi,ith)=zbp
        else if(ians.eq.2) then
           zb1ay(ixi,ith)=zb
           zb2ay(ixi,ith)=zbinb
        else if(ians.eq.3) then
           zb1ay(ixi,ith)=zb
           zb2ay(ixi,ith)=zbfs
        else if(ians.eq.4) then
           zb1ay(ixi,ith)=zbinb
           zb2ay(ixi,ith)=zbfs
        end if
        !
     end do
  end do
  !
  !  menus for comparitive & standalone plots
  !
  if(ians.eq.1) then
     call r8_grf3f2(zxay,zb1ay,zb2ay, &
          zthay,zwkay,inxi,inth,inxi, &
          zxilab,'bpol(boozer)','bpol(mombs)',zthlab, &
          zxiuns,zbuns,zthuns,'poloidal field comparison', &
          runid,0)
  else if(ians.eq.2) then
     call r8_grf3f2(zxay,zb1ay,zb2ay, &
          zthay,zwkay,inxi,inth,inxi, &
          zxilab,'mod(b)(boozer)','mod(b)(1/(1/b)f-s)',zthlab, &
          zxiuns,zbuns,zthuns,'mod(b) field comparison', &
          runid,0)
  else if(ians.eq.3) then
     call r8_grf3f2(zxay,zb1ay,zb2ay, &
          zthay,zwkay,inxi,inth,inxi, &
          zxilab,'mod(b)(boozer)','mod(b)(mombs)',zthlab, &
          zxiuns,zbuns,zthuns,'mod(b) field comparison', &
          runid,0)
  else if(ians.eq.4) then
     call r8_grf3f2(zxay,zb1ay,zb2ay, &
          zthay,zwkay,inxi,inth,inxi, &
          zxilab,'mod(b)(1/(1/b)f-s)','mod(b)(mombs)',zthlab, &
          zxiuns,zbuns,zthuns,'mod(b) field comparison', &
          runid,0)
  end if
  !
  goto 100
  !
1000 continue
  !      call mombs_ngint                  ! restore ngint
  return
  !
end subroutine chkspb
! 11sep2004 fgtok
