subroutine philb(zxi,zth,zph,zdphdxi,zdphdth)
  !
  !  philb -- return the toroidal arc spanned by the field line from
  !    theta=-pi to the passed theta value (zth) on the specified flux
  !    surface (zxi); also return the derivatives w.r.t. xi and theta.
  !
  !  input:  (zxi,zth) -- mag. coordinates zxi and theta
  !
  ! output:  zph = arc length
  !          zdphdxi = d(arc length)/dxi
  !          zdphdth = d(arc length)/dth
  !
  !  method:  the precomputed bicubic spline functions theta*
  !           ("theta-star") and nu are related to the desired quantities
  !           as follows:
  !
  !   philb(zxi,zth) = q*(theta* + pi) - nu
  !
  !           where q = 1/iota(bar) is evaluated at zxi
  !           and theta* and nu at (zxi,zth)
  !
  !  and the derivatives are formed from this (see trdebug/bpsplot.fdb also).
  !
  !------------------------
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: intvl,jx,ilinx,ilinth,ier,ii,jj
  !============
  real(fp) :: zth,zph,zdphdxi,zdphdth,zxi,zf1d,zzh,zf0,zf1,zf2,zf3
  real(fp) :: zdf1d,zx0,zh,ziotb,zq,zdiot,zdqdxi,zdxi,zdth,zths
  real(fp) :: zdthsdxi,zdthsdth,znu,zdnudxi,zdnudth
  !============
  real(fp) :: zvalnu(6),zvalth(6)  ! place to return spline evals & derivs
  !
  integer :: isel(6)
  data isel/1,1,1,0,0,0/
  !
  intvl=11       ! spline data assumed to be on evenly spaced grid
  !
  jx = int(zxi/dxi(lcentr,2)) + 1
  zx0 = real(jx-1,fp)*dxi(lcentr,2)
  zh = zxi - zx0
  !
  ziotb = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx), &
       spliot(4,jx))
  !
  !  q = 1/iota(bar)
  !
  zq=1.0_fp/ziotb
  !
  !  dq/dxi = -q*q*d(iota(bar))/dxi
  !
  zdiot=zdf1d(zh,spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zdqdxi=-zq*zq*zdiot
  !
  !  when calling splines get value & derivatives
  !
  !  return elements (1) - value, (2) - xi derivative, (3) - th derivative
  !
  ilinx=1
  ilinth=1
  ier=0
  !
  call bcspevxy(zxi,zth,xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth,ii,jj,zdxi,zdth,ier)
  !
  if(ier.gt.0) then
     write(nonlin,*) ' ?philb:  spline evaluation error.'
  end if
  !
  call bcspevfn(isel,1,1,zvalth,ii,jj,zdxi,zdth, &
       datspl(lspthets),nzp1,nthspl)
  !
  call bcspevfn(isel,1,1,zvalnu,ii,jj,zdxi,zdth, &
       datspl(lspphilb),nzp1,nthspl)
  !
  zths=zvalth(1)
  zdthsdxi=zvalth(2)
  zdthsdth=zvalth(3)
  !
  znu=zvalnu(1)
  zdnudxi=zvalnu(2)
  zdnudth=zvalnu(3)
  !
  !  phi(lb):
  !
  zph=zq*(zths+xpi) - znu
  !
  !  derivatives
  !
  zdphdxi=zdqdxi*(zths+xpi) + zq*zdthsdxi - zdnudxi
  !
  zdphdth=zq*zdthsdth - zdnudth
  !
  return
end subroutine philb
! 11sep2004 fgtok
