subroutine nbmomry(zxii,zthi,ijac,zr,zz,zjac)
 
  use iso_c_binding, only: fp => c_double
  use nbspline_mod
  use physconst_mod, only: twopi
  !  use bicubic splines to evaluate what was done by the legacy
  !  momry routine
 
  implicit none
  real(fp), intent(in) :: zxii,zthi        ! mag. coordinates in
  integer, intent(in) :: ijac        ! =0: no jacobian; .gt.0: jacobian
  real(fp), intent(out) :: zr,zz         ! r,z coordinates out
  real(fp), intent(out) :: zjac(2,2)     ! 2x2 jacobian out
 
  !
  !  zjac(1,1) = dr/dxi     zjac(1,2) = dr/dth
  !  zjac(2,1) = dz/dxi     zjac(2,2) = dz/dth
  !
 
  real(fp) zxi,zth,zval(3)
  integer iwant(6),iwantj(6),ierr
  real(fp) zthi8,zth8
 
  data iwant /1,0,0,0,0,0/
  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))
 
  if(ixi.ge.nbrk) then
 
     ! external spline
     ixi = ixi-nbrk+1
 
     if(ijac.le.0) then
 
        ! values only
 
        nbsii=ith
        nbsjj=ixi

        call fvbicub(iwant,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,rsplx,nth,next)
        zr=zval(1)
 
        call fvbicub(iwant,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,zsplx,nth,next)
        zz=zval(1)
 
     else
 
        ! values and jacobian

        nbsii=ith
        nbsjj=ixi
 
        call fvbicub(iwantj,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,rsplx,nth,next)
        zr=zval(1)
        zjac(1,1)=zval(3)
        zjac(1,2)=zval(2)
 
        call fvbicub(iwantj,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,zsplx,nth,next)
        zz=zval(1)
        zjac(2,1)=zval(3)
        zjac(2,2)=zval(2)
 
     endif
 
  else
 
     ! internal spline
 
 
     if(ijac.le.0) then
 
        ! values only

        nbsii=ith
        nbsjj=ixi
 
        call fvbicub(iwant,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,rspl,nth,nint)
        zr=zval(1)
 
        call fvbicub(iwant,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,zspl,nth,nint)
        zz=zval(1)
 
     else
 
        ! values and jacobian

        nbsii=ith
        nbsjj=ixi
 
        call fvbicub(iwantj,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,rspl,nth,nint)
        zr=zval(1)
        zjac(1,1)=zval(3)
        zjac(1,2)=zval(2)
 
        call fvbicub(iwantj,1,1,zval,nbsii,nbsjj,thparam,xiparam, &
             hth,hthi,hxi,hxii,zspl,nth,nint)
        zz=zval(1)
        zjac(2,1)=zval(3)
        zjac(2,2)=zval(2)
 
     endif
 
  endif
 
end subroutine nbmomry
! 11sep2004 fgtok -s r8_precision.sub nubeam_test_r8.sub "r8con.csh conversion"
