subroutine nbsplj(zxii,zthi,zjac,zjacth,zjacxi)
 
  use iso_c_binding, only: fp => c_double
  use nbspline_mod
 
  !  use bicubic splines to evaluate what was done by the legacy
  !  jintrp routine
 
  implicit none
  real(fp), intent(in) :: zxii,zthi        ! mag. coordinates in
  real(fp), intent(out) :: zjac,zjacth,zjacxi  ! j,dj/dth,dj/dxi
 
  !
 
  real(fp) zxi,zth,zval(3),twopi
  integer iwantj(6),ierr
  real(fp) zthi8,zth8
 
  data iwantj /1,1,1,0,0,0/
 
  integer ith,ixi
  real(fp) xiparam(1),thparam(1)
 
  !---------------------------------------------
  ! normalize range
 
  zxi=max(xi(1),min(xi(nxi),zxii))
  if((zthi.ge.th(1)).and.(zthi.le.th(nth))) then
     zth=zthi
  else
     if(zthi.lt.th(1)) zth=zthi+twopi
     if(zthi.gt.th(nth)) zth=zthi-twopi
     if((zth.lt.th(1)).or.(zth.gt.th(nth))) then
        zthi8=zthi
        zth8=mod(zthi8,thmax)
        if(zthi8.lt.0) zth8=zth8+thmax
        zth=zth8
        zth=max(th(1),min(th(nth),zth))
     endif
  endif
 
  ixi=min( (nxi-1), int(1+xspani*(zxi-xi(1))*(nxi-1)) )
  if(xi(ixi).gt.zxi) then
     ixi=ixi-1
  endif
  if(xi(ixi+1).lt.zxi) then
     ixi=min((nxi-1),(ixi+1))
  endif
  xiparam=(zxi-xi(ixi))/(xi(ixi+1)-xi(ixi))
 
  ith=min( (nth-1), int(1+thspani*(zth-th(1))*(nth-1)) )
  if(th(ith).gt.zth) then
     ith=max(1,(ith-1))
  endif
  if(th(ith+1).lt.zth) then
     ith=min((nth-1),(ith+1))
  endif
  thparam=(zth-th(ith))/(th(ith+1)-th(ith))
 
  nbsii=ith
  nbsjj=ixi

  call fvbicub(iwantj,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
       hth,hthi,hxi,hxii,jspl,nth,nxi)
 
  zjac=zval(1)
  zjacth=zval(2)
  zjacxi=zval(3)
 
end subroutine nbsplj
! 11sep2004 fgtok -s r8_precision.sub nubeam_test_r8.sub "r8con.csh conversion"
