!-----------------------------------------------------------------------
!  chkspl -- check out spline functions & fourier splines
!
subroutine chkspl
  use nbi_com
  use nbi_com, only : zpi=>xpi
  use nbi_dimensions
  !
  !	5/15/90  tbt  added code for the derivatives of j.
  !
  implicit none
  LOGICAL, EXTERNAL :: LCHAR,LXOUT,LYES,LYESC                !
  !
  !============
  integer :: lnout,ixi,ith,itot,i,ifun2,ix,it,idim,iloc
  !============
  real(fp) :: zf1d,zzh,zf0,zf1,zf2,zf3,zxi1,zxi2,zximin,zximax
  real(fp) :: zth1,zth2,zthmin,zthmax,zxil,zthl,zans,zr,zy
  !============
  integer :: imax1,imax2
  parameter (imax1=1024)
  parameter (imax2=32768)
  !
  real(fp) :: zxi(imax1),zth(imax1),zfunc(imax2),ztemp(imax1)
  real(fp) :: zjac(2,2)
  !
  character(len=16) :: zzlbl
  !
  !  executable code
  !  look at 1d spline or 2d spline ?
  !
10 continue
  call neupag
  !
  if(lxout(0)) write(lnout(0),1001) imax1,imax2
1001 format(/' chkspl:  plot 1d splines vs xi or '/ &
       '          2d fourier splines vs xi, theta'/ &
       ' choose dimensionality, range and no. of points first'/ &
       ' choice of function offerred last'/ &
       '  max 1d spline evaluation:  imax1 = ',i5,' points'/ &
       '  max 2d spline evaluation:  imax2 = ',i5,' points'/ &
       '  or enter "b" for b field representation comparisons'/ &
       '  or enter "c" for bicubic spline evaluations'/ &
       '  or enter "d" for derived quantities plots')
  !
  call uread( &
       'chkspl: enter 1, 2, b or c, or q for 1d, 2d spline or quit:$')
  !
  if(lchar('q')) then
     return
  else if(lchar('1')) then
     idim=1
  else if(lchar('2')) then
     idim=2
  else if(lchar('c')) then
15   continue
     if((lspphilb.eq.0).or.(lspthets.eq.0)) then
        call uermsg( &
             'chkspl:  no nu/phi/theta boozer splines.',' ',0)
        goto 10
     end if
     if(lxout(0)) write(lnout(0),1002)
1002 format(/' bicubic spline vs. (xi,theta) selection:'/ &
          ' n -- nu (boozer coord. conversion, zeta = phi + nu'/ &
          ' p -- phi arc length along field line'/ &
          ' t -- theta (boozer coord.)'/ &
          ' q -- quit')
     call uread('chkspl:  bicubic spline option p/t/q:$')
     if(lchar('q')) then
        goto 10
     else if(lchar('p')) then
        call bpsplot(-1,'field line phi arc')
     else if(lchar('n')) then
        call bpsplot(lspphilb,'boozer nu')
     else if(lchar('t')) then
        call bpsplot(lspthets,'boozer theta')
     else
        call uermsg('chkspl:  invalid bicubic selection',' ',0)
     end if
     goto 15
  else if(lchar('b')) then
     call chkspb
     goto 10
  else if(lchar('d')) then
     call chkspd
     goto 10
  else
     call uermsg('chkspl:  invalid option',' ',0)
     goto 10
  end if
  !
20 continue
  call ureadi('chkspl:  no. of "xi" points (.le.imax1):$',ixi)
  if((ixi.lt.2).or.(ixi.gt.imax1)) then
     call uermsi( &
          '?chkspl:  no. of "xi" pts must be btw 2 and imax1',ixi,0)
     goto 20
  end if
  !
  if(idim.eq.1) goto 40
  !
30 continue
  call ureadi('chkspl:  no. of "th" points (.le.imax1):$',ith)
  if((ith.lt.2).or.(ith.gt.imax1)) then
     call uermsi( &
          '?chkspl:  no. of "th" pts must be btw 2 and imax1',ith,0)
     goto 30
  end if
  !
  itot=ixi*ith
  if(itot.gt.imax2) then
     call uermsi( &
          '?chkspl:  tot. no. "xi" pts x "th" pts exceeds imax2', &
          itot,0)
     goto 20
  end if
  !
40 continue
  !
  call ureadr8('chkspl:  enter "xi" limit no. 1:$$',zxi1)
  call ureadr8('chkspl:  enter "xi" limit no. 2:$$',zxi2)
  zximin=min(zxi1,zxi2)
  zximax=max(zxi1,zxi2)
  !
  zximin=max(0.0_fp,min(2.0_fp,zximin))
  zximax=max(0.0_fp,min(2.0_fp,zximax))
  !
  if(zximin.ge.zximax) then
     write(lnout(0),9001) zximin,zximax
9001 format(' ?chkspl:  zximin=',1pe11.4,' zximax=',1pe11.4/ &
          '  input "xi" limits must be distinct and btw 0.0 and 2.0')
     goto 40
  end if
  !
  do i=1,ixi
     zxi(i)=((ixi-i)*zximin+(i-1)*zximax)/real(ixi-1,fp)

  end do
  !
50 continue
  if(idim.eq.1) goto 60
  !
  call ureadr8('chkspl:  enter "th" limit no. 1:$$',zth1)
  call ureadr8('chkspl:  enter "th" limit no. 2:$$',zth2)
  zthmin=min(zth1,zth2)
  zthmax=max(zth1,zth2)
  !
  zthmin=max(-zpi,min(zpi,zthmin))
  zthmax=max(-zpi,min(zpi,zthmax))
  !
  if(zthmin.ge.zthmax) then
     write(lnout(0),9002) zthmin,zthmax
9002 format(' ?chkspl:  zthmin=',1pe11.4,' zthmax=',1pe11.4/ &
          '  input "th" limits must be distinct and btw -pi and pi')
     goto 50
  end if
  !
  do i=1,ith
     zth(i)=((ith-i)*zthmin+(i-1)*zthmax)/real(ith-1)

  end do
  !
60 continue
  !
  if(idim.eq.1) then
     if(lxout(0)) write(lnout(0),2001)
2001 format(/' options:'/ &
          ' i - plot splimhd - i-mhd'/ &
          ' j - plot iota(bar) spline'/ &
          ' g - plot splg - g-mhd'/ &
          ' v - plot splphi - phi-r - radial electric potential'/ &
          ' 0 - quit'/)
     call uread('chkspl:  enter 1d function choice:$')
     !
     if(lchar('0')) then
        goto 10
     else if(lchar('i')) then
        call chksp1(zxi,zfunc,ixi,splimhd,mtbl,'i-mhd')
     else if(lchar('j')) then
        call chksp1(zxi,zfunc,ixi,spliot,mtbl,'iota(bar)')
     else if(lchar('g')) then
        call chksp1(zxi,zfunc,ixi,splg,mtbl,'g-mhd')
     else if(lchar('v')) then
        call chksp1(zxi,zfunc,ixi,splphi,mtbl,'v-pot')
     else
        call uermsg('?chkspl -- unrecognized function code',' ',0)
        goto 60   !  tbt  5/90
     end if
     !
  else
     !
     !   add ploting of dj/d(xi), dj/d(th), d(1/b)/d(xi), d2(1/b)/d2(xi), etc.
     !
     if(lxout(0)) write(lnout(0),2002)
2002 format(/' options:'/ &
          ' b - plot 1/mod(b)'/ &
          ' s - plot d(1/b)/d(theta)'/ &
          ' t - plot d(1/b)/d(xi)'/ &
          ' u - plot d2(1/b)/d2(theta)'/ &
          ' v - plot d2(1/b)/d2(xi)'/ &
          ' w - plot d2(1/b)/d(theta)d(xi)'/ &
          '                               '/ &
          ' j - plot det(jacobean) from orbit code'/ &
          ' d - plot d(j)/d(theta) from orbit code'/ &
          ' e - plot d(j)/d(xi)    from orbit code'/ &
          ' k - plot det(jacobean) from nbmomry code -- spline interp.'/ &
          ' r - plot major radius (nbmomry)'/ &
          ' y - plot elevation (nbmomry)'/)
     !
     call uread('chkspl:  enter 2d function choice:$')
     !
     if(lchar('0')) then
        goto 10
     else if(lchar('q')) then
        goto 10
     else if(lchar('b')) then
        ifun2=1
        zzlbl='1/mod(b) orbg'
     else if(lchar('s')) then
        ifun2=8
        zzlbl='d(1/b)/d(th) '
     else if(lchar('t')) then
        ifun2=9
        zzlbl='d(1/b)/d(xi) '
     else if(lchar('u')) then
        ifun2=10
        zzlbl='d2(1/b)/d2(th)'
     else if(lchar('v')) then
        ifun2=11
        zzlbl='d2(1/b)/d2(xi)'
     else if(lchar('w')) then
        ifun2=12
        zzlbl='d2(1/b)/dthdxi'

     else if(lchar('j')) then
        ifun2=2
        zzlbl='det[j] orbg'
     else if(lchar('d')) then
        ifun2=66
        zzlbl='d(j)/d(theta)'
     else if(lchar('e')) then
        ifun2=7
        zzlbl='d(j)/d(xi)'
     else if(lchar('k')) then
        ifun2=3
        zzlbl='det[j] nbmomry(s)'
     else if(lchar('l')) then
        ifun2=6
        zzlbl='det[j] nbmomry(l)'
     else if(lchar('r')) then
        ifun2=4
        zzlbl='r nbmomry'
     else if(lchar('y')) then
        ifun2=5
        zzlbl='y nbmomry'
     else
        call uermsg('?chkspl -- unrecognized function code',' ',0)
        goto 60   ! tbt 5/90
     end if
     !
     do ix=1,ixi
        !
        zxil=zxi(ix)
        do it=1,ith
           !
           zthl=zth(it)
           !
           if(ifun2.eq.1) then
              !  1/b orbg
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bi1ovb(1)
              end if
              !
           else if(ifun2.eq.8) then
              !  d(1/b)/d(theta)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bibdts(1)
              end if

              !
           else if(ifun2.eq.9) then
              !  d(1/b)/d(xi)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bibdxi(1)
              end if

              !
           else if(ifun2.eq.10) then
              !  d2(1/b)/d2(theta)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bib2ts(1)
              end if

              !
           else if(ifun2.eq.11) then
              !  d2(1/b)/d2(xi)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bib2xi(1)
              end if

              !
           else if(ifun2.eq.12) then
              !  d2(1/b)/d(theta)d(xi)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call bintrp
                 zans=bibtsx(1)
              end if

           else if(ifun2.eq.2) then
              !  jacobean orbg
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call jintrp
                 zans=bijac(1)
              end if

           else if(ifun2.eq.66) then
              !  d(j)/d(theta)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call jintrp
                 zans=bijdts(1)
              end if

           else if(ifun2.eq.7) then
              !  d(j)/d(xi)
              if(zxil.gt.(xbmbnd+0.5_fp*dxizc(ledge))) then
                 zans=0.0_fp
              else
                 bixi(1)=zxil
                 bithst(1)=zthl
                 call jintrp
                 zans=bijdxi(1)
              end if

           else if(ifun2.eq.3) then
              !  jacobean nbmomry -- spline
              if(zxil.gt.xiblo(nznblo+lcm1)) then
                 zans=0.0_fp
              else
                 call nbmomry(zxil,zthl,1,zr,zy,zjac)
                 zans=(zjac(1,1)*zjac(2,2)-zjac(1,2)*zjac(2,1))
              end if
           else if(ifun2.eq.6) then
              !  jacobean nbmomry -- linear
              if(zxil.gt.xiblo(nznblo+lcm1)) then
                 zans=0.0_fp
              else
                 call nbmomry(zxil,zthl,3,zr,zy,zjac)
                 zans=(zjac(1,1)*zjac(2,2)-zjac(1,2)*zjac(2,1))
              end if
           else if(ifun2.eq.4) then
              !  r nbmomry
              if(zxil.gt.xiblo(nznblo+lcm1)) then
                 zans=0.0_fp
              else
                 call nbmomry(zxil,zthl,0,zr,zy,zjac)
                 zans=zr
              end if
           else if(ifun2.eq.5) then
              !  y nbmomry
              if(zxil.gt.xiblo(nznblo+lcm1)) then
                 zans=0.0_fp
              else
                 call nbmomry(zxil,zthl,0,zr,zy,zjac)
                 zans=zy
              end if
           end if
           !
           iloc=(ix-1)*ith+it
           zfunc(iloc)=zans
           !

        end do

     end do
     !
     call r8_grf3f1(zth,zfunc,zxi,ztemp,ith,ixi,ith, &
          'theta',zzlbl,'xi',' ',' ',' ', &
          'transp run '//runid(1:lrunid)//' debug plot', &
          'run '//runid(1:lrunid),0)
     !
  end if
  !
  goto 60
  !
end subroutine chkspl
! 11sep2004 fgtok
