subroutine bpsplot(iaddrf,zlabel)
  use nbi_com
  use nbi_com, only: zpi=>xpi
  use nbi_dimensions
  !
  !  test plot bicubic spline function and/or derivatives
  !
  ! input:
  !  iaddrf -- address of bicubic spline data in datspl
  !      or negative number for calculated functions
  !         -1;  phi arc length
  !  zlabel -- label for data
  !
  !============
  implicit none
  integer :: inx,inth,imx,ix
  integer :: ith,lnout,ibad,ideriv,ilinx,ilinth,ier
  !============
  real(fp) :: zx1,zx2,zth1,zth2,zxmin,zxmax,zthmin,zthmax,zth,zx
  real(fp) :: zval
  !============
  integer :: iaddrf
  character(len=*) zlabel
  !
  character(len=20) :: zlocal
  !
  integer :: isela(6,0:5)
  !
  logical :: lchar,lxout
  !
  real(fp), dimension(:), allocatable :: zxay,zthay,zwkay
  real(fp), dimension(:,:), allocatable :: zdatay
  !
  data isela/ &
       1,0,0,0,0,0, &
       0,1,0,0,0,0, &
       0,0,1,0,0,0, &
       0,0,0,1,0,0, &
       0,0,0,0,1,0, &
       0,0,0,0,0,1/
  !------------------------------------------------------------------
  !
  zx1=0.0_fp
  zx2=1.0_fp
  zth1=-zpi
  zth2=zpi
  !
10 continue
  !
  write(6,'(/'' plotting:  spline of:  '',a/)') zlabel
  !
  call ureadi('no. of x points:$',inx)
  inx=max(3,min(1001,inx))
  call ureadi('no. of theta points:$',inth)
  inth=max(3,min(1001,inth))
  !
  call ureadr8('minimum x (btw 0 and 1):$',zx1)
  call ureadr8('maximum x (btw 0 and 1):$',zx2)
  !
  zxmin=max(0.0_fp,min(1.0_fp,min(zx1,zx2)))
  zxmax=max(0.0_fp,min(1.0_fp,max(zx1,zx2)))
  !
  if(zxmax.le.zxmin) then
     call uermsg('% x range or separation error',' ',0)
     goto 10
  end if
  !
  call ureadr8('minimum th (btw -pi and pi):$',zth1)
  call ureadr8('maximum th (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 or separation error',' ',0)
     goto 10
  end if
  !
  !-----------------------
  !
  imx=max(inx,inth)
  allocate(zxay(inx),zthay(inth),zwkay(imx)); zwkay = 0.0
  allocate(zdatay(inx,inth))
  !
  do ix=1,inx
     zxay(ix)=zxmin+real(ix-1,fp)*(zxmax-zxmin)/real(inx-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(6,'(/'' plotting:  spline of f = '',a/)') zlabel
     write(lnout(0),1101)
1101 format(/ &
          ' options:'/ &
          '  r  -- reset (xi,theta) plotting range'/ &
          '  q  -- quit'/ &
          '  a  -- plot f'/ &
          '  b  -- plot df/dxi'/ &
          '  c  -- plot df/dtheta')
     if(iaddrf.gt.0) write(lnout(0),1102)
1102 format( &
          '  d  -- plot d2f/dxi2'/ &
          '  e  -- plot d2f/dtheta2'/ &
          '  f  -- plot d2f/dxidth')
  end if
  call uread(' bpsplot:  enter opcode:$')
  !
  ibad=0
  if(lchar('q')) then
     deallocate(zxay,zthay,zwkay,zdatay)
     return
  else if(lchar('r')) then
     goto 10
  else if(lchar('a')) then
     ideriv=0
     zlocal='f'
  else if(lchar('b')) then
     ideriv=1
     zlocal='df/dxi'
  else if(lchar('c')) then
     ideriv=2
     zlocal='df/dth'
  else if(lchar('d')) then
     ideriv=3
     zlocal='d2f/dxi2'
     if(iaddrf.lt.0) ibad=1
  else if(lchar('e')) then
     ideriv=4
     zlocal='d2f/dth2'
     if(iaddrf.lt.0) ibad=1
  else if(lchar('f')) then
     ideriv=5
     zlocal='d2f/dxidth'
     if(iaddrf.lt.0) ibad=1
  else
     ibad=1
  end if
  if(ibad.eq.1) then
     call uermsg(' ?bpsplot:  invalid option.',' ',0)
     goto 100
  end if
  !
  ilinx=1
  ilinth=1
  do ith=1,inth
     zth=zthay(ith)
     do ix=1,inx
        zx=zxay(ix)
        if(iaddrf.gt.0) then
           call bcspeval(zx,zth,isela(1,ideriv),zval, &
                xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
                ilinx,ilinth, &
                datspl(iaddrf),nzp1,ier)
        else
           call bpsplot_cal(zx,zth,iaddrf,ideriv,isela,ilinx,ilinth, &
                zval,ier)
        end if
        if(ier.ne.0) then
           call uermsi(' %bcspeval error occurred, ier=',ier,0)
           goto 100
        end if
        zdatay(ix,ith)=zval
     end do
  end do
  !
  call r8_grf3f1(zxay,zdatay,zthay,zwkay, &
       inx,inth,inx,'xi',zlocal//' : '//zlabel,'theta', &
       'xi',zlocal,'th',runid//' bicubic spline plot', &
       ' ',0)
  !
  goto 100
end subroutine bpsplot
!-----------------
!  bpsplot_cal -- calculations
!
subroutine bpsplot_cal(zxi,zth,icode,ideriv,isela,ilinx,ilinth, &
     zval,ier)
  !
  use nbi_com
  use nbi_com, only: zpi=>xpi
  use nbi_dimensions
  !
  !============
  implicit none
  integer :: icode,ideriv,ilinx,ilinth,ier,jx,i,ierl
  !============
  real(fp) :: zth,zval,zxi,zf1d,zzh,zf0,zf1,zf2,zf3,zdf1d,zx0,zh
  real(fp) :: ziotb,zq,zdiot,zdqdxi
  !============
  real(fp) :: zfths(6),zfnus(6)
  integer :: isela(6,0:5)
  !
  integer :: iselnu(6),iselth(6)
   !
  ier=0
  if(icode.ne.-1) then
     ier=-1
     goto 900
  end if
  if(ideriv.gt.2) then
     ier=-2
     goto 900
  end if
  !
  ier=0
  !
  !  calculate function or derivative vs. transp xi,theta:
  !
  !    splines:  theta*(xi,theta) and nu(xi,theta)
  !
  !    philb(xi,theta) = phi-arclength of fieldline on xi surface starting
  !                      at th=-pi and going to th=theta
  !
  !               = q*(theta*(xi,theta) + pi) - nu(xi,theta)
  !
  !    return this for ideriv=0
  !
  !  derivatives:
  !
  !    ideriv=2:
  !   d(philb)/dtheta = q*dtheta*/dtheta - dnu/dtheta
  !    ideriv=1:
  !   d(philb)/dxi = dq/dxi*(theta*(xi,theta) + pi) +q*dtheta*/dxi - dnu/dxi
  !
  jx = int(zxi/dxi(lcentr,2)) + 1
  zx0 = real(jx-1,fp)*dxi(lcentr,2)
  zh = zxi - zx0
  !
  do i=1,6
     iselnu(i)=isela(i,ideriv)
     iselth(i)=0
  end do
  !
  ziotb = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx), &
       spliot(4,jx))
  zq=1.0_fp/ziotb
  !
  if(ideriv.eq.0) then
     iselth(1)=1
  else if(ideriv.eq.1) then
     iselth(1)=1
     iselth(2)=1
     zdiot=zdf1d(zh,spliot(2,jx),spliot(3,jx),spliot(4,jx))
     zdqdxi=-zq*zq*zdiot
  else if(ideriv.eq.2) then
     iselth(3)=1
  end if
  !
  call bcspeval(zxi,zth,iselth,zfths, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspthets),nzp1,ierl)
  ier=max(ier,ierl)
  !
  call bcspeval(zxi,zth,iselnu,zfnus, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspphilb),nzp1,ierl)
  ier=max(ier,ierl)
  !
  !  zfnus(1) always contains the proper nu derivative
  !  zfths(1) contains the th derivative, or else
  !   zfths(1) contains th and zfths(2) contains the derivative.
  !
  if(ideriv.eq.0) then
     zval = zq*(zfths(1)+zpi) - zfnus(1)
  else if(ideriv.eq.1) then
     zval = zdqdxi*(zfths(1)+zpi) + zq*zfths(2) - zfnus(1)
  else if(ideriv.eq.2) then
     zval = zq*zfths(1) - zfnus(1)
  end if
  !
900 continue
  return
end subroutine bpsplot_cal
! 11sep2004 fgtok
