subroutine bphip(ierr)
  !
  !  at many flux surfaces:  evaluate b(theta), find all b singularities,
  !  and evaluate phi(p)
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: ifree,i,iaoff,isurf,inumth,iath,iab,iaph,ith,ito,iuse
  integer :: iaps,iathb,iaflag
  !============
  real(fp) :: zxi,zth,zb,zdbdth,zd2bdth2,zth1,zb1,zthb,zp,zphib
  !============
  integer,intent(out) :: ierr       ! status code on exit (0=ok)
  !
  character(len=60) :: zlbl                 ! for debugging plot
  !
  !  to get debugging going...
  !
  if(nphipzc.eq.0) nphipzc=2
  !
  !-----------------------------------------
  !
  ierr=0
  !
  nphipc=nzones*nphipzc
  !
  !  ndatspl(...) gives information about each surface, and where
  !    surface data is stored in datspl.
  !
  !  ndatspl(1 thru nphipc): ptrs to each surface's ndatspl data block.
  !    i.e. block for surface k starts at ndatspl(k).
  !
  !  within each ndatspl block, i.e. for the kth surface:
  !
  !  item#1:  number of theta pts.
  !  item#2:  address (in datspl) of surface's theta grid
  !  item#3:  address of b field strength saved at grid pts
  !  item#4:  address of phi(b) calc. and saved at grid pts
  !  item#5:  number of maxima in b(theta) fcn, this surface
  !  item#6:  address of singularity description data
  !  item#7:  address of auxilliary singularity data
  !  item#8:  address of phi(p) calc. and saved at grid pts
  !  item#9:  address of theta(b) -- poloidal arc length, radians
  !  item#10 and up:  theta-sorted list of singularities.
  !
  !  surfaces to be considered:  xi=i*1.0_fp/nphipc, i=1 to nphipc.
  !
  ifree=nphipc+1
  !
  do i=1,nphipc
     !
     ndatspl(i)=ifree
     iaoff=ifree-1
     !
     isurf=i
     zxi=real(i,fp)/real(nphipc,fp)
     !
     inumth=nthpmin+zxi*(nthpmax-nthpmin)
     ndatspl(iaoff+1)=inumth
     !
     ndatspl(iaoff+2)=lspuse
     iath=lspuse
     !
     ndatspl(iaoff+3)=lspuse+inumth
     iab=lspuse+inumth
     !
     ndatspl(iaoff+4)=lspuse+2*inumth
     iaph=lspuse+2*inumth
     !
     lspuse=lspuse+3*inumth
     !
     do ith=1,inumth
        ito=ith-1
        zth=-xpi + (real(ith-1,fp)/real(inumth-1,fp))*(xpi+xpi)
        datspl(iath+ito)=zth
        !
        call bphip_gb(zxi,zth,zb,zdbdth,zd2bdth2)
        datspl(iab+ito)=zb
        !
        !  temporarily store db/dth
        !
        datspl(iaph+ito)=zdbdth
        if(nlsym2b) then
           if((ith.eq.1).or.(ith.eq.inumth)) datspl(iaph+ito)=0.0_fp
        end if
        !
     end do
     !
     !  analyze singularities
     !
     call bphip_sng(zxi,isurf,iaoff,iuse,ierr)
     if(ierr.ne.0) return
     ifree=iuse+1
     !
     !  evaluate p(psi) for phi(p), and phi(b)
     !
     iaps=ndatspl(iaoff+8)
     iathb=ndatspl(iaoff+9)
     !
     !  check dimensions
     !
     if(iaps.gt.midspl-1000) then
        ierr=ierr+1
        write(nonlin,*) ' ?bphip:  datspl array too small, rebuild!'
     end if
     if(ifree+inumth.gt.minspl) then
        ierr=ierr+1
        write(nonlin,*) ' ?bphip: ndatspl array too small, rebuild!'
     end if
     if(ierr.gt.0) return
     !
     do ith=1,inumth
        ito=ith-1
        zth1=datspl(iath+ito)
        zb1 =datspl(iab+ito)
        call bphip_phif(isurf,zth1,zb1,zthb,zp,zphib,iaflag,ierr)
        if(ierr.ne.0) return
        ndatspl(ifree+ito)=iaflag
        datspl(iaph+ito)=zphib
        datspl(iaps+ito)=zp
        datspl(iathb+ito)=zthb
     end do
     !
     !  assymptote patch on p(psi)
     !
     call bphip_minp(datspl(iath), &
          datspl(iaps),ndatspl(ifree),inumth,iaoff)
     !
     !dbg         write(zlbl,'('' surface zxi='',1pe11.4)') zxi
     !dbg         call grafx1(datspl(iath),datspl(iaph),inumth,
     !dbg     >     'radians','radians',runid,zlbl,
     !dbg     >     'phi(b) vs. bounce angle')
     !dbg         call grafx1(datspl(iath),datspl(iaps),inumth,
     !dbg     >     'radians','radians/rho',runid,zlbl,
     !dbg     >     'p(psi) vs. bounce angle')
     !
  end do
  !
  !  debug option
  !
  !dbg      call chkspl
  !
  ! (not working: modularization) c%      call chkphip
  !
  return
end subroutine bphip
!
!---------------------
!
subroutine bphip_gb(zxi,zth,zb,zdbdth,zd2bdth2)
  !
  !  evaluate b field & 1st & 2nd derivatives w.r.t. theta(transp).
  !
  use nbi_com
  !
  !============
  implicit none
  real(fp) :: zth,zb,zdbdth,zd2bdth2,zxi,zbgaus,zdbgaus,zdbg2
  integer :: ierr
  !============
  bixi(1)=zxi
  bithst(1)=zth
  call bintrp
  zbgaus=1.0_fp/bi1ovb(1)
  zdbgaus=-zbgaus*zbgaus*bibdts(1)
  zdbg2=zbgaus*zbgaus*(2.0_fp*zbgaus*bibdts(1)*bibdts(1)-bib2ts(1))
  !
  !  b, db/dth in tesla...
  !
  zb=1.e-4_fp*zbgaus
  zdbdth=1.e-4_fp*zdbgaus
  zd2bdth2=1.e-4_fp*zdbg2
  !
  return
end subroutine bphip_gb
!
!----------------------
!
subroutine bphip_sng(zxi,isurf,iaoff,iuse,ierr)
  !
  !  find b field singularities (i.e. db/dtheta=0 points, corresponding
  !  to local maxima or minima)
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: isurf,iaoff,iuse,imax_sng,inumth,iath,iab,iadb,iammx
  integer :: iammx0,imax,il1,il2,iskip,iskipn,ith,im,iaa,iammn
  integer :: iammn0,imin,insing,iaux,inper,iphip,isinga,icm,iic
  integer :: iauxmn,ii,icthwm,ioffhwm,iauxhwm,imx2,idir,i2x,iaminp
  integer :: imingot,imaxgot,icmin,icmax,ithfndp,iauxmx,iaminm
  integer :: ithfndm,ilast,ia,ifirst,ins,itype,is,itypep,is2,ifnd
  integer :: iasave,inums,iam,iap,imaxp,inhwm,imaxfrom,iside
  integer :: idegen,inumdg,ia0,idegen2,iad,ioff,ierr
  !============
  real(fp) :: zxi,zdb1,zdbn,zdbcut,zdb2,zth1,zth2,zthmax,zbmax
  real(fp) :: zb2sng,bphip_sngk,zthtmp,zbtmp,zb2tmp,zthmin,zbmin
  real(fp) :: zbphlm,bphip_plim,zbphtmp,zphib,zphipmn,zthmax2,zbmax2
  real(fp) :: zthbminp,zbminp,zthbmaxp,zbmaxp,zthmax2a,zbmax2a
  real(fp) :: zthfndp,zdelthp,zphibp,zthbminm,zbminm,zthbmaxm,zbmaxm
  real(fp) :: zthfndm,zdelthm,zphibm,zthlim,zthlow,zthlimp,zthloc
  real(fp) :: nbi_d1mach,zthmx,zk,zfac,zthmk,zk2,ztol
  !============
  parameter (imax_sng=50)
  integer :: iccmin(2,imax_sng),iccmax(2,imax_sng)
  !
  !--------------------------------------------
  !
  ierr=0
  !
  inumth=ndatspl(iaoff+1)
  iath=ndatspl(iaoff+2)
  iab=ndatspl(iaoff+3)
  iadb=ndatspl(iaoff+4)
  !
  ndatspl(iaoff+6)=lspuse   ! address of singularity data
  iammx=lspuse
  iammx0=iammx
  !
  !  1.  look for maxima
  !
  imax=0
  !
  !     special treatment at the branch cut
  !
  zdb1=datspl(iadb)
  zdbn=datspl(iadb+inumth-1)
  zdbcut=0.5_fp*(zdb1+zdbn)
  datspl(iadb)=zdbcut
  datspl(iadb+inumth-1)=zdbcut
  !
  !  main search loop
  !
  il1=1
  il2=inumth-il1
  !
  iskip=0
  iskipn=0
  !
  zdb2=datspl(iadb+il1-1)
  do ith=il1,il2
     !
     if(ith.eq.il2) then
        if(iskipn.eq.1) then
           iskip=1
        end if
     end if
     !
     zdb1=zdb2
     zdb2=datspl(iadb+ith)
     if((zdb1.ge.0.0_fp).and.(zdb2.le.0.0_fp).and.(iskip.eq.0)) then
        !  found one
        imax=imax+1
        zth1=datspl(iath+ith-1)
        zth2=datspl(iath+ith)
        call bphip_sfnd(zxi,zth1,zth2,zdb1,zdb2,zthmax,zbmax,ierr)
        !mg  bug          if(ierr.ne.1) return
        if(ierr.ne.0) return
        zb2sng=bphip_sngk(zxi,zthmax,ierr)
        if(ierr.ne.1) return
        if(imax.gt.1) then
           !  keep descending order
           do im=1,imax-1
              iaa=iammx0+3*(im-1)
              if(zbmax.gt.datspl(iaa+1)) then
                 zthtmp=datspl(iaa)
                 zbtmp=datspl(iaa+1)
                 zb2tmp=datspl(iaa+2)
                 datspl(iaa)=zthmax
                 datspl(iaa+1)=zbmax
                 datspl(iaa+2)=zb2sng
                 zthmax=zthtmp
                 zbmax=zbtmp
                 zb2sng=zb2tmp
              end if
           end do
        end if  ! imax.gt.1
        datspl(iammx)=zthmax
        datspl(iammx+1)=zbmax
        datspl(iammx+2)=zb2sng
        iammx=iammx+3               ! leave space for possible later use
        iskip=1                     ! prevent dbl detection of end pt max
        if(ith.eq.1) iskipn=1
     else
        iskip=0
     end if  ! local maxima detected
  end do  ! search loop
  !
  !  save no. of maxima found
  !
  ndatspl(iaoff+5)=imax
  !
  !  update free pointer
  !
  lspuse=iammx
  !
  !  2.  look for minima
  !
  iammn=lspuse
  iammn0=iammn
  !
  imin=0
  !
  !  main search loop
  !
  il1=1
  il2=inumth-il1
  !
  iskip=0
  iskipn=0
  !
  zdb2=datspl(iadb+il1-1)
  do ith=il1,il2
     !
     if(ith.eq.il2) then
        if(iskipn.eq.1) then
           iskip=1
        end if
     end if
     !
     zdb1=zdb2
     zdb2=datspl(iadb+ith)
     if((zdb1.le.0.0_fp).and.(zdb2.ge.0.0_fp).and.(iskip.eq.0)) then
        !  found one
        imin=imin+1
        zth1=datspl(iath+ith-1)
        zth2=datspl(iath+ith)
        call bphip_sfnd(zxi,zth1,zth2,zdb1,zdb2,zthmin,zbmin,ierr)
        !mg bug            if(ierr.ne.1) return
        if(ierr.ne.0) return
        zbphlm=bphip_plim(zxi,zthmin,ierr)
        if(ierr.ne.0) return
        if(imin.gt.1) then
           !  keep ascending order in theta
           do im=1,imin-1
              iaa=iammn0+3*(im-1)
              if(zthmin.lt.datspl(iaa)) then
                 zthtmp=datspl(iaa)
                 zbtmp=datspl(iaa+1)
                 zbphtmp=datspl(iaa+2)
                 datspl(iaa)=zthmin
                 datspl(iaa+1)=zbmin
                 datspl(iaa+2)=zbphlm
                 zthmin=zthtmp
                 zbmin=zbtmp
                 zbphlm=zbphtmp
              end if
           end do
        end if  ! imin.gt.1
        datspl(iammn)=zthmin
        datspl(iammn+1)=zbmin
        !  evaluate limiting p(psi) value
        datspl(iammn+2)=zbphlm
        iammn=iammn+3
        iskip=1                     ! prevent dbl detection of end pt min
        if(ith.eq.1) iskipn=1
     else
        iskip=0
     end if  ! local minima detected
  end do  ! search loop
  !
  !  update free pointer -- allow space for high water mark points
  !  allows space for interval characterization
  !
  insing=4*imax   ! add space for 2 addl for branch cut handling
  iaux=iammn + 2*(imax - 1)
  ndatspl(iaoff+7)=iaux    ! addr auxilliary singularity data
  inper=6                           ! #data items / singularity
  iphip=iaux+inper*(insing-2)
  ndatspl(iaoff+8)=iphip   ! addr phip(p) at grid pts
  ndatspl(iaoff+9)=iphip+inumth  ! addr theta bounce data
  lspuse=iphip + 2*inumth
  !
  isinga=iaoff+10
  iuse=iaoff+9+insing
  !
  if(imin.ne.imax) then
     write(nonlin,'('' ?bphip_sng:  imax.ne.imin'')')
     ierr=1; return
  end if
  if(imax.gt.imax_sng) then
     write(nonlin,'('' ?bphip_sng:  too many singularities'')')
     ierr=1; return
  end if
  !
  do icm=1,imax
     iccmin(1,icm)=0
     iccmin(2,icm)=0
     iccmax(1,icm)=0
     iccmax(2,icm)=0
  end do
  !
  !  poloidal and toroidal orbit arclength limits at surface global maximum
  datspl(iaux)=-twopi
  datspl(iaux+1)=+twopi
  iic = 1
  call bphip_phib(isurf,-xpi,xpi,zphib,iic,0.0_fp,ierr)
  if(ierr.ne.0) return
  datspl(iaux+2)=zphib
  datspl(iaux+3)=zphib
  !
  !  ...and at all minima...
  !
  iauxmn=iaux+inper*imax
  do ii=1,imin
     datspl(iauxmn+inper*(ii-1))=0.0_fp
     datspl(iauxmn+inper*(ii-1)+1)=0.0_fp
     datspl(iauxmn+inper*(ii-1)+2)=0.0_fp
     datspl(iauxmn+inper*(ii-1)+3)=0.0_fp
     zphipmn=datspl(iammn0+3*(ii-1)+2)
     datspl(iauxmn+inper*(ii-1)+4)=zphipmn
     datspl(iauxmn+inper*(ii-1)+5)=zphipmn
  end do
  !
  !  3.  "high water" marks --  for local maxima that are not the highest
  !      maximum, there exist points which define the other side of the
  !      magnetic well delimitted by the local maxima:  find these points.
  !
  icthwm=0
  ioffhwm=iaux+inper*(imin+imax)
  iauxhwm=ioffhwm
  do imx2=2,imax
     iaa=iammx0+3*(imx2-1)
     !
     zthmax2=datspl(iaa)
     zbmax2=datspl(iaa+1)
     !
     !  find monotonic search range in +theta direction
     idir=1
     i2x=int(2.0_fp+0.5_fp*idir)
     call bphip_mmx(isurf,zthmax2,zbmax2,idir, &
          zthbminp,zbminp,iaminp,imingot,zthbmaxp,zbmaxp,imaxgot, &
          zthmax2a,zbmax2a)
     !  find exact match in +theta direction
     icmin=iccmin(i2x,imingot)+1
     iccmin(i2x,imingot)=icmin
     icmax=iccmax(i2x,imaxgot)+1
     iccmax(i2x,imaxgot)=icmax
     call bphip_bfnd(isurf,zbmax2,zthbminp,zbminp,zthbmaxp,zbmaxp, &
          zthfndp,ithfndp,icmin,icmax,ierr)
     if(ierr.ne.0) return
     !
     datspl(iammn+icthwm)=zthfndp
     icthwm=icthwm+1
     !
     !  singularity data:  delta(theta), delta(phi)
     zdelthp=zthfndp-zthmax2
     iic=1
     call bphip_phib(isurf,zthmax2,zthfndp,zphibp,iic,0.0_fp,ierr)
     if(ierr.ne.0) return
     datspl(ioffhwm)=-zdelthp
     datspl(ioffhwm+2)=zphibp
     iauxmx=iaux+inper*(imx2-1)
     datspl(iauxmx+1)=zdelthp
     datspl(iauxmx+3)=zphibp
     !
     !  find monotonic search range in -theta direction
     idir=-1
     i2x=int(2.0_fp+0.5_fp*idir)
     call bphip_mmx(isurf,zthmax2,zbmax2,idir, &
          zthbminm,zbminm,iaminm,imingot,zthbmaxm,zbmaxm,imaxgot, &
          zthmax2a,zbmax2a)
     !  find exact match in -theta direction
     icmin=iccmin(i2x,imingot)+1
     iccmin(i2x,imingot)=icmin
     icmax=iccmax(i2x,imaxgot)+1
     iccmax(i2x,imaxgot)=icmax
     call bphip_bfnd(isurf,zbmax2,zthbminm,zbminm,zthbmaxm,zbmaxm, &
          zthfndm,ithfndm,icmin,icmax,ierr)
     if(ierr.ne.0) return
     !
     datspl(iammn+icthwm)=zthfndm
     icthwm=icthwm+1
     !  singularity data:
     zdelthm=zthmax2-zthfndm
     iic=-1
     call bphip_phib(isurf,zthmax2,zthfndm,zphibm,iic,0.0_fp,ierr)
     if(ierr.ne.0) return
     datspl(ioffhwm+inper+1)=zdelthm
     datspl(ioffhwm+inper+3)=zphibm
     datspl(iauxmx)=-zdelthm
     datspl(iauxmx+2)=zphibm
     !
     !  span both gaps, just above the high water marks.
     datspl(ioffhwm+1)=-zdelthm-zdelthp
     datspl(ioffhwm+3)=zphibm+zphibp
     datspl(ioffhwm+inper)=zdelthm+zdelthp
     datspl(ioffhwm+inper+2)=zphibm+zphibp
     !
     !dbg         write(nonlin,8888)
     !dbg     >      zthmax2,zbmax2,zthfndm,zthfndp,ithfndm,ithfndp,
     !dbg     >      inumth,zxi
     !
     ioffhwm=ioffhwm+2*inper
  end do
8888 format(/' th(local_max)=',1pe11.4,' b(local_max)=',1pe11.4, &
       ' th- = ',1pe11.4,' th+ = ',1pe11.4/ &
       '    ithfndm,p=',i4,1x,i4,' inumth=',i4,' .... zxi = ',1pe11.4)
  !
  !--------------------------
  !  4.  singularity pointers sorted in order of increasing theta
  !      allow space for one extra ptr at each end, for "wrapping" the
  !      periodic data set.
  !
  ilast=iuse
  ia=isinga
  ifirst=ia
  !
  ins=insing-2
  zthlim=-epsinv
  zthlow=-epsinv
  itype=0
  !
  do is=1,ins
     ia=ia+1
     zthlimp=zthlim
     zthlim=epsinv
     itypep=itype
     !  scan maxima
     do is2=1,imax
        zthloc=datspl(iammx0+3*(is2-1))
        if((zthloc.lt.zthlim).and.(zthloc.gt.zthlow)) then
           ifnd=is2
           zthlim=zthloc
           itype=1
        end if
     end do
     !  scan minima
     do is2=1,imin
        zthloc=datspl(iammn0+3*(is2-1))
        if(zthloc.lt.zthlim) then
           if((zthloc.gt.zthlow).or. &
                ((itypep.eq.1).and.(zthloc.eq.zthlow))) then
              ifnd=is2+imax
              zthlim=zthloc
              itype=2
           end if
        end if
     end do
     !  scan high water marks
     do is2=1,2*(imax-1)
        zthloc=datspl(iammn+is2-1)
        if(zthloc.lt.zthlim) then
           if((zthloc.gt.zthlow).or. &
                ((itypep.eq.1).and.(zthloc.eq.zthlow))) then
              ifnd=is2+2*imax
              zthlim=zthloc
              itype=3
           end if
        end if
     end do
     if(zthlim.lt.zthlimp) then
        ierr=1; return  ! must find *new* theta glb each time thru
     else if((zthlim.lt.-xpi).or.(zthlim.gt.xpi)) then
        ierr=1; return  ! not in range (-pi,pi)
     else
        zthlow=zthlim
        ndatspl(ia)=ifnd
        if(ifnd.eq.1) iasave=ia   ! global maximum
     end if
  end do
  !  wrap ends
  ndatspl(ilast)=ndatspl(ifirst+1)
  ndatspl(ifirst)=ndatspl(ilast-1)
  !
  !  now place the singularity constants for phi(p) for the maxima and
  !  high water mark points.  phi(p) singularities look like
  !     k*ln(delta(theta)) locally where k is related to d2b/dtheta2.
  !  if two neighbouring maxima are equal in magnitude, the k btw them
  !  doubles as the phi(p) integral is singular at both ends.
  !
  !  high-water-mark points normally see k/2 due to db/dtheta .ne. 0
  !  there.  see notes, dmc nov 1996
  !
  !  set standard values for maxima -- correct if neighbouring maxima
  !  is not of the same height.
  !
  ztol=3.0_fp*nbi_d1mach(4)*xpi
  !
  ia=iasave-1
  do ii=1,ins
     ia=ia+1
     if(ia.ge.ilast) ia=ifirst+1
     inums=ndatspl(ia)
     if(inums.le.imax) then
        !  maxima
        zthmx=datspl(iammx0+3*(inums-1))
        zk=datspl(iammx0+3*(inums-1)+2)
        if(inums.eq.1) then
           zfac=2.0_fp  ! global
        else
           zfac=1.0_fp  ! local
        end if
        iam=iaux+(inums-1)*inper+4
        !  set if not set already due to hwm processing
        datspl(iam)=zfac*zk
        iap=iaux+(inums-1)*inper+5
        datspl(iap)=zfac*zk
        imaxp=inums
     end if
  end do
  !
  !  set high water marks and check for "degeneracy":  two neighbouring
  !  maxima with the same value.
  !
  ia=iasave-1
  do ii=1,ins
     ia=ia+1
     if(ia.ge.ilast) ia=ifirst+1
     inums=ndatspl(ia)
     if(inums.gt.(2*imax)) then
        !  high water mark
        inhwm=inums-2*imax
        zthmk=datspl(iammn+inhwm-1)
        imaxfrom=1 + (1+inhwm)/2   ! of which maxima
        zk=datspl(iammx0+3*(imaxfrom-1)+2)
        iside=2*(imaxfrom-1)-inhwm ! 1: + side, 0: - side of maxima
        !  check for degeneracy -- theta match with neighbouring maxima
        call bphip_dck(ia,ifirst,ilast,iammx0, &
             imax,zthmk,ztol,idegen,inumdg)
        !
        iam=iauxhwm+(inhwm-1)*inper+4
        iap=iam+1
        if(idegen.eq.0) then
           !  normal case
           if(iside.eq.1) then
              datspl(iam)=0.5_fp*zk
              datspl(iap)=zk
           else
              datspl(iam)=zk
              datspl(iap)=0.5_fp*zk
           end if
        else
           !  degenerate case:  this means two successive equal maxima, that
           !  need to be connected.
           zk2=datspl(iammx0+3*(inumdg-1)+2)
           if(iside.eq.1) then
              datspl(iam)=zk+zk2
              datspl(iap)=0.0_fp  ! don't know yet
              datspl(iaux+(imaxfrom-1)*inper+5)=zk+zk2
              datspl(iaux+(inumdg-1)*inper+4)=zk+zk2
           else
              datspl(iap)=zk+zk2
              datspl(iam)=0.0_fp  ! don't know yet
              datspl(iaux+(imaxfrom-1)*inper+4)=zk+zk2
              datspl(iaux+(inumdg-1)*inper+5)=zk+zk2
           end if
        end if  ! degeneracy
     end if  ! is a high water mark
  end do
  !
  !  ok now can pick up unknown pts
  ia=iasave-1
  do ii=1,ins
     ia=ia+1
     if(ia.ge.ilast) ia=ifirst+1
     inums=ndatspl(ia)
     if(inums.gt.(2*imax)) then
        !  high water mark
        inhwm=inums-2*imax
        zthmk=datspl(iammn+inhwm-1)
        ia0=iauxhwm+(inhwm-1)*inper
        iam=ia0+4
        iap=iam+1
        if((datspl(iam).eq.0.0_fp).or.(datspl(iap).eq.0.0_fp)) then
           idegen2=1
        else
           idegen2=0
        end if
        call bphip_dck(ia,ifirst,ilast,iammx0, &
             imax,zthmk,ztol,idegen,inumdg)
        if(idegen.ne.idegen2) then
           ierr=1; return  ! this loop not consistent w/prev. loop
        end if
        !  copy aux data from singularity to which hwm is degenerate
        if(idegen.eq.1) then
           iad=iaux+(inumdg-1)*inper
           do ioff=0,5
              datspl(ia0+ioff)=datspl(iad+ioff)
              datspl(ia0+ioff)=datspl(iad+ioff)
           end do
        end if
     end if  ! is a hwm
  end do
  !
  return
end subroutine bphip_sng
!
!--------------------------------------------------
subroutine bphip_dck(ia,ifirst,ilast,iammx0, &
     imax,zthmk,ztol,idegen,inumdg)
  !
  !  "degeneracy detector" subroutine for bphip_sng
  !
  use nbi_com
  !
  !------------------
  !
  !  look at neighbouring singularities
  !============
  implicit none
  integer :: ifirst,ilast,iammx0,imax,idegen,inumdg,ia,iapeek1
  integer :: iapeek2,ins1,ins2
  !============
  real(fp) :: zthmk,zthmx,ztol
  !============
  iapeek1=ia-1
  if(iapeek1.lt.ifirst) iapeek1=ilast
  iapeek2=ia+1
  if(iapeek2.gt.ilast) iapeek2=ifirst
  idegen=0
  ins1=ndatspl(iapeek1)
  ins2=ndatspl(iapeek2)
  !
  !  are they maxima?  do their locations coincide with this hwm
  !  location?
  if(ins1.le.imax) then
     zthmx=datspl(iammx0+3*(ins1-1))
     if(iapeek1.eq.ifirst) zthmx=zthmx-twopi
     if(iapeek1.eq.ilast) zthmx=zthmx+twopi
     if(abs(zthmx-zthmk).le.ztol) then
        idegen=1
        inumdg=ins1
     end if
  end if
  if(ins2.le.imax) then
     zthmx=datspl(iammx0+3*(ins2-1))
     if(iapeek2.eq.ifirst) zthmx=zthmx-twopi
     if(iapeek2.eq.ilast) zthmx=zthmx+twopi
     if(abs(zthmx-zthmk).le.ztol) then
        idegen=1
        inumdg=ins2
     end if
  end if
  !
  return
end subroutine bphip_dck
!---------------------------------
real(fp) function bphip_sngk(zxi,ztheta,ierr)
  !
  use nbi_com
  !
  !  given the (zxi,ztheta) coordinates of a b point that is a local
  !  maximum in the direction of variation of ztheta, compute the quantity
  !
  !     (r*bz)ext/(dpsi/dxi) * (i + g*q) / 2*b**2 * (db/dxi)
  !        * dtheta(gwb)/dtheta(transp) * sqrt( -b/d2b/dtheta2 )
  !
  !  which is useful in characterizing the precesion phi(p) singularity
  !  at a field maximum in the zero banana width approximation.
  !
  !---------------------------------
  !
  !============
  implicit none
  integer :: jx,ilinx,ilinth,ierr
  !============
  real(fp) :: ztheta,zxi,zxip,ztmp,zfac,zi,zg,zq,zf1d
  real(fp) :: zx0,zh,ziot,zcv,zdpsidx,zb,zdbdth,zd2bdth2,zbg
  real(fp) :: zdbdxi,zarg,zsqrt,zdthcon,zans
  !============
  integer :: isel(6)
  !
  data isel/0,0,1,0,0,0/
  !
  data zxip/1.0e10_fp/
  data ztmp/-1.0e10_fp/
  !
  save zxip,ztmp,zfac,zi,zg,zq
  !
  ierr=0
  !
  if((zxi.ne.zxip).or.(ztmp.ne.tbm1)) then
     zxip=zxi
     ztmp=tbm1
     !
     !  flux surface constants:  i, g, q, (r*bt)vacuum/(dpsi/dxi)
     !  eval once on this surface at this time.
     !
     jx=int(zxi/dxi(lcentr,2)) + 1
     zx0=real(jx-1,fp)*dxi(lcentr,2)
     zh=zxi-zx0
     !
     ziot=zf1d(zh, &
          spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
     zq=1.0_fp/ziot
     !
     zcv=1.0e-4_fp  !t/gauss
     !
     zg=zcv*zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
     !
     zi=zcv*zf1d(zh,splimhd(1,jx), &
          splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
     !
     zdpsidx=(zxi*1.e4_fp*tflux)/(zq*xpi)
     zfac=bzxr/zdpsidx
     !
  end if
  !
  !  b and derivatives at singularity
  !
  call bphip_gb(zxi,ztheta,zb,zdbdth,zd2bdth2)
  zbg=1.e4_fp*zb
  zdbdxi=-1.e-4_fp*zbg*zbg*bibdxi(1)
  !
  !  sqrt(b/(d2b/dth2))
  !
  zarg=-zd2bdth2/zb
  zsqrt=sqrt(1.0_fp/max(1.0e-8_fp,zarg))
  !
  !  dtheta(gwb)/dtheta(transp)
  !
  ilinx=1
  ilinth=1
  !
  call bcspeval(zxi,ztheta,isel,zdthcon, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspthets),nzp1,ierr)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?bphip_sngk: bcspeval (spline) error.'
     return
  end if
  !
  zans=zfac*(zi+zg*zq)*zdbdxi*zdthcon*zsqrt/(2.0_fp*zb*zb)
  if(zans.eq.0.0_fp) then
     zans=epslon                    ! never return 0.0
  end if
  !
  bphip_sngk=zans
  !
  return
end function bphip_sngk
!---------------------------------
!
subroutine bphip_sfnd(zxi,zth1,zth2,zdb1,zdb2, &
     zthmax,zbmax,ierr)
  !
  !  find singularity on zxi surface btw zth1 and zth2; singularity is
  !  defined as the point where db/dth = 0.0.
  !
  !  zdb1, zdb2 are the db/dth values at the endpoints and are of
  !  oposite sign; db/dth is continuous.
  !
  use nbi_com
  !
  !----------------------------
  !
  !============
  implicit none
  integer :: iter,itrap,ierr
  !============
  real(fp) :: zth1,zth2,zdb1,zdb2,zthmax,zbmax,zxi,zdthtol,zdb1a
  real(fp) :: zdb2a,zth1a,zth2a,zthlim1,zdblim1,zthlim2,zdblim2,zth
  real(fp) :: zd2b,zthnex,zdb,zbg
  !============
  ierr=0
  !
  zdthtol=0.01_fp*(zth2-zth1)
  !
  bixi(1)=zxi
  !
  if(zth2.le.zth1) then
     write(nonlin,*) 'bphip_sfnd: zth1,zth2 args out of order: ', &
          zth1,zth2
     ierr=1
     return
  end if
  !
  if(zdb1.eq.0.0_fp) then
     zthmax=zth1
     goto 900
  else if(zdb2.eq.0.0_fp) then
     zthmax=zth2
     goto 900
  else
     zdb1a=zdb1
     zdb2a=zdb2
     zth1a=zth1
     zth2a=zth2
     !
     iter=0
     !
100  continue
     iter=iter+1
     if(iter.gt.10) then
        write(nonlin,*) ' %lots of iterations in bphip_sfnd'
     end if
     if(iter.gt.100) then
        write(nonlin,*) ' ?too many iterations in bphip_sfnd'
        ierr=1; return                !  bug check
     end if
     if((zdb1a*zdb2a).lt.0.0_fp) then
        itrap=1
        !  next guess -- trapzoidal rule
        zthlim1=zth1a
        zdblim1=zdb1a
        zthlim2=zth2a
        zdblim2=zdb2a
        !
        zth=(zdb2a*zth1a-zdb1a*zth2a)/(zdb2a-zdb1a)
        !
     else
        itrap=0
        !  both previous estimates to one side (and iter>1); try a newton step
        if(zd2b.eq.0.0_fp) then
           zthnex=epsinv
        else
           zthnex=zth-zdb/zd2b
        end if
        if((zthnex.lt.zthlim1).or.(zthnex.gt.zthlim2)) then
           !  out of range; try to reset
           if((zdb*zdblim1).lt.0.0_fp) then
              zth=0.5_fp*(zth+zthlim1)
           else
              zth=0.5_fp*(zth+zthlim2)
           end if
        else
           !  accept the newton step
           zth=zthnex
        end if
     end if
     !  evaluate at this guess
     bithst(1)=zth
     call bintrp
     zbg=1.0_fp/bi1ovb(1)
     zdb=-1.0e-4_fp*zbg*zbg*bibdts(1)
     if(zdb.eq.0.0_fp) goto 200
     zd2b=1.0e-4_fp*zbg*zbg*(2.0_fp*zbg*(bibdts(1))**2-bib2ts(1))
     !  do another guess
     if(itrap.eq.1) then
        if(abs(zdb).lt.0.1_fp*max(abs(zdb1a),abs(zdb2a))) then
           !  go by modulus if reduction was achieved
           itrap=0
        end if
     end if
     if(itrap.eq.1) then
        !  keep db/dtheta bracketed
        if((zdb1a*zdb).lt.0.0_fp) then
           zdb2a=zdb
           zth2a=zth
        else
           zdb1a=zdb
           zth1a=zth
        end if
     else
        if(abs(zdb1a).lt.abs(zdb2a)) then
           zdb2a=zdb
           zth2a=zth
        else
           zdb1a=zdb
           zth1a=zth
        end if
     end if
     !  otherwise keep
     if(abs(zth2a-zth1a).lt.zdthtol) goto 200
     goto 100
  end if
  !
200 continue
  !
  zthmax=zth
  if(zthmax.gt.zth2) then
     if((zthmax-zth2).gt.zdthtol) then
        write(nonlin,*) ' ?bphip_sfnd: zth2 range error.'
        ierr=1; return
     end if
     zthmax=zth2
     goto 900
  else if(zthmax.lt.zth1) then
     if((zth1-zthmax).gt.zdthtol) then
        write(nonlin,*) ' ?bphip_sfnd: zth1 range error.'
        ierr=1; return
     end if
     zthmax=zth1
     goto 900
  else
     zbmax=1.e-4_fp*zbg !tesla
     goto 1000
  end if
  !
900 continue
  !
  !  endpoint match:  evaluate b there.
  !
  bithst(1)=zthmax
  call bintrp
  zbmax=1.0e-4_fp/bi1ovb(1) !tesla
  !
  !  exit
  !
1000 continue
  return
end subroutine bphip_sfnd
!
!---------------------------------
subroutine bphip_mmx(isurf,zthstart,zbwant,idir, &
     zthbmin,zbmin,iamin,imingot,zthbmax,zbmax,imaxgot, &
     zthbmax2,zbmax2)
  !
  use nbi_com
  !
  !  find search range for remote reflection point of orbit that has
  !  bounced at a particular b field value; i.e. find a b field maximum
  !  on this surface that is greater than the b field value given, and
  !  then find the next closer minimum; the search range is the monotonic
  !  region between the minimum and the maximum.
  !
  !  also find the index of the minimum found (iamin)
  !  and the bounding maximum in the -idir direction from this minimum
  !  (zthbmax2,zbmax2)
  !
  !  input:
  !============
  implicit none
  integer :: iaoff,imax,iammx0,iammn0,im,ia0,imxsmal,imxfnd,imn
  integer :: imnfnd,imx
  !============
  real(fp) :: zbtol,nbi_d1mach,zb0,zthmx,zthbmaxp,zthmn,zthbminp
  !============
  integer :: isurf   ! surface index
  real(fp) :: zthstart   ! starting theta location
  real(fp) :: zbwant     ! b field value sought
  integer :: idir    ! search direction +1 / -1
  !     ** but if idir=0 on input, compute idir here and idir is output **
  !
  !  output:
  real(fp) :: zthbmin     ! theta of bmin
  real(fp) :: zbmin      ! bmin value
  integer :: iamin    ! address of minima data
  integer :: imingot ! minima index
  !
  real(fp) :: zthbmax     ! theta of bmax
  real(fp) :: zbmax      ! bmax value
  integer :: imaxgot  ! maxima index
  !
  real(fp) :: zthbmax2    ! theta of bmax on other side of minimum
  real(fp) :: zbmax2     ! value of bmax on other side (-idir side).
  !
  !------------------------------------------
  !
  if(idir.eq.0) then
     call bphip_dir(isurf,zthstart,iaoff,imax,iammx0,iammn0,idir)
  else
     iaoff=ndatspl(isurf)-1
     !
     imax=ndatspl(iaoff+5)
     iammx0=ndatspl(iaoff+6)
     iammn0=iammx0+3*imax
  end if
  !
  !  find smallest maximum .gt. zbwant
  !
  zbtol=zbwant*3.0_fp*nbi_d1mach(4)
  do im=1,imax
     ia0=iammx0+3*(im-1)
     zb0=datspl(ia0+1)
     if(zb0.le.zbwant-zbtol) goto 10
  end do
  im=imax+1
10 continue
  imxsmal=max(1,(im-1))
  !
  !  of the maxima .gt. zbwant find the one that is closest in the `idir'
  !  direction
  !
  zthbmax=idir*epsinv
  do im=1,imxsmal
     ia0=iammx0+3*(im-1)
     zthmx=datspl(ia0)
     zthbmaxp=zthbmax
     if(idir.eq.1) then
        if(zthmx.gt.zthstart) then
           zthbmax=min(zthbmax,zthmx)
        else
           zthbmax=min(zthbmax,(zthmx+twopi))
        end if
     else
        if(zthmx.lt.zthstart) then
           zthbmax=max(zthbmax,zthmx)
        else
           zthbmax=max(zthbmax,(zthmx-twopi))
        end if
     end if
     if(zthbmax.ne.zthbmaxp) imxfnd=im
  end do
  zbmax=datspl(iammx0+3*(imxfnd-1)+1)
  imaxgot=imxfnd
  !
  !  find the first minimum that is in the -`idir' direction from the
  !  first bounding maximum just found.
  !
  zthbmin=-idir*epsinv
  do imn=1,imax
     ia0=iammn0+3*(imn-1)
     zthmn=datspl(ia0)
     if(idir.eq.1) then
        if(zthmn.lt.zthstart) then
           zthmn=zthmn+twopi
        end if
     else
        if(zthmn.gt.zthstart) then
           zthmn=zthmn-twopi
        end if
     end if
     zthbminp=zthbmin
     if(idir.eq.-1) then
        if(zthmn.gt.zthbmax) then
           zthbmin=min(zthbmin,zthmn)
        end if
     else
        if(zthmn.lt.zthbmax) then
           zthbmin=max(zthbmin,zthmn)
        end if
     end if
     if(zthbmin.ne.zthbminp) imnfnd=imn
  end do
  iamin=iammn0+3*(imnfnd-1)
  zbmin=datspl(iamin+1)
  imingot=imnfnd
  !
  !  find the first maximum that is in the -`idir' direction from the
  !  bounding minimum just found.
  !
  zthbmax2=-idir*epsinv
  do imx=1,imax
     ia0=iammx0+3*(imx-1)
     zthmx=datspl(ia0)
     if(idir.eq.1) then
        if(zthmx.gt.zthbmin) then
           zthmx=zthmx-twopi
        end if
     else
        if(zthmx.lt.zthbmin) then
           zthmx=zthmx+twopi
        end if
     end if
     zthbmaxp=zthbmax2
     if(idir.eq.-1) then
        if(zthmx.gt.zthbmin) then
           zthbmax2=min(zthbmax2,zthmx)
        end if
     else
        if(zthmx.lt.zthbmin) then
           zthbmax2=max(zthbmax2,zthmx)
        end if
     end if
     if(zthbmax2.ne.zthbmaxp) imxfnd=imx
  end do
  zbmax2=datspl(iammx0+3*(imxfnd-1)+1)
  !
  return
end subroutine bphip_mmx
!-------------------------------------------------------------------
subroutine bphip_bfnd(isurf,zbwant,zthbmin,zbmin,zthbmax,zbmax, &
     zthfnd,ithfnd,icmin,icmax,ierr)
  !
  use nbi_com
  !
  !  find matching b (=zbwant) on specified surface, searching the
  !  monotonic region from zthbmin to zthbmax; note that zthbmin might
  !  be larger than zthbmax; also one or the other might be past the
  !  (-pi,+pi) branch cut (but by no more than 2pi).
  !
  !  zbwant is bracketed by zbmin and zbmax.
  !   the theta value returned is btw zthbmin and zthbmax but will be
  !   normalized into the range (-pi,pi).  for reasons of supporting
  !   the "sorting" logic of calling routines, the output theta value will
  !   be bound away from the end points of the theta range (zthbmin,zthbmax).
  !
  !  input:
  !
  !============
  implicit none
  integer :: iaoff,inumth,iath,iab,idir,inzn,iznbmin,ith_bphip
  integer :: iznbmax,izn1,izn2,izntest,izz,iter,ierr
  !============
  real(fp) :: zbmax,zztol0,nbi_d1mach,zbtol,zdthtol,zthinc,zxi,zthbmin2
  real(fp) :: zthbmax2,zth1,zth2,zb1,zb2,zb1a,zb2a,zth1a,zth2a,ztmp
  real(fp) :: zthlim1,zblim1,zthlim2,zblim2,zth,zdb,zthnex,zb,zbg
  !============
  integer :: isurf                     ! the surface to search on
  real(fp) :: zbwant                       ! target b value
  real(fp) :: zthbmin                      ! theta of local bmin
  real(fp) :: zbmin                        ! bmin
  real(fp) :: zthbmax                      ! theta of local bmax
  !
  integer :: icmin             ! # of del(th) increments from min to bind
  integer :: icmax             ! # of del(th) increments from max to bind
  !
  !  output:
  !
  real(fp) :: zthfnd                       ! theta where b value is matched
  integer :: ithfnd                    ! index to theta zone
  !
  !----------------------------
  ierr = 0
  !
  zztol0=nbi_d1mach(4)
  zbtol=10.0_fp*zztol0*zbwant
  !
  zdthtol=twopi*1.0e-5_fp
  !
  iaoff=ndatspl(isurf)-1
  !
  inumth=ndatspl(iaoff+1)
  iath=ndatspl(iaoff+2)
  iab=ndatspl(iaoff+3)
  !
  !  check endpts; done if already w/in tolerance to these...
  !
  if(zthbmin.lt.zthbmax) then
     idir=1
  else
     idir=-1
  end if
  !
  if(zbwant.gt.zbmax-zbtol) then
     if(zbwant.gt.zbmax+zbtol) write(nonlin,8902) zbwant,zbmax
     zthinc=idir*icmax*zztol0*abs(zthbmax)
     zthfnd=zthbmax-zthinc
     goto 900
  else if(zbwant.lt.zbmin+zbtol) then
     if(zbwant.lt.zbmin-zbtol) write(nonlin,8901) zbwant,zbmin
     zthinc=idir*icmin*zztol0*abs(zthbmin)
     zthfnd=zthbmin+zthinc
     goto 900
  end if
8901 format(' %bphip_bfnd warning:  zbwant .lt. zbmin: ',2(1x,1pe11.4))
8902 format(' %bphip_bfnd warning:  zbwant .gt. zbmax: ',2(1x,1pe11.4))
  !
  zxi=real(isurf,fp)/nphipc
  bixi(1)=zxi
  !
  !----------------------------------------------------
  !  find the theta zone within which to search more precisely...
  !
  !  get indices of zones containing end pts...
  !   number indices as: 0,1,2,...,(inumth-2)
  !    there are (inumth) bdys and (inumth-1) zones defined.
  !
  inzn=inumth-1
  !
  !  (zthbmin2,zthbmax2) are forced into range (-pi,pi)
  !
  zthbmin2=zthbmin
  iznbmin=ith_bphip(zthbmin2,iath,inumth)-1
  if(zthbmin.lt.-xpi) iznbmin=iznbmin-inzn
  if(zthbmin.gt.xpi) iznbmin=iznbmin+inzn
  !
  zthbmax2=zthbmax
  iznbmax=ith_bphip(zthbmax2,iath,inumth)-1
  if(zthbmax.lt.-xpi) iznbmax=iznbmax-inzn
  if(zthbmax.gt.xpi) iznbmax=iznbmax+inzn
  !
  izn1=iznbmin
  izn2=iznbmax
  !
10 continue
  izntest=(izn1+izn2)/2
  if(izntest.lt.0) then
     izz=izntest+inzn
  else if(izntest.ge.inzn) then
     izz=izntest-inzn
  else
     izz=izntest
  end if
  !
  if(izntest.eq.iznbmin) then
     if(idir.eq.1) then
        zth1=zthbmin2
        zth2=datspl(iath+izz+1)
        zb1=zbmin
        zb2=datspl(iab+izz+1)
     else
        zth1=datspl(iath+izz)
        zth2=zthbmin2
        zb1=datspl(iab+izz)
        zb2=zbmin
     end if
  else if(izntest.eq.iznbmax) then
     if(idir.eq.1) then
        zth1=datspl(iath+izz)
        zth2=zthbmax2
        zb1=datspl(iab+izz)
        zb2=zbmax
     else
        zth1=zthbmax2
        zth2=datspl(iath+izz+1)
        zb1=zbmax
        zb2=datspl(iab+izz+1)
     end if
  else
     zth1=datspl(iath+izz)
     zth2=datspl(iath+izz+1)
     zb1=datspl(iab+izz)
     zb2=datspl(iab+izz+1)
  end if
  if((zb1-zbwant)*(zb2-zbwant).le.0.0_fp) then
     !  zone found...
     ithfnd=izz+1
     goto 50
  end if
  !
  !  not in range -- employ binary search
  !
  if(izn1.eq.izn2) then
     write(nonlin,*) ' ?bphip_bfnd algorithm failure.'
     ierr=1
     return
  end if
  !
  if(zbwant.gt.max(zb1,zb2)) then
     izn1=izntest+idir
  else
     izn2=izntest-idir
  end if
  !
  goto 10
  !
  !----------------------------------------------------
  !  now do a combined newton / trapezoidal-rule search to
  !  find the matching b precisely.
  !
50 continue
  !
  if(zb1.eq.zbwant) then
     zthfnd=zth1
     goto 900
  else if(zb2.eq.zbwant) then
     zthfnd=zth2
     goto 900
  else
     zb1a=zb1
     zb2a=zb2
     zth1a=zth1
     zth2a=zth2
     !
     iter=0
     !
100  continue
     !
     !  check ordering
     if(zth1a.gt.zth2a) then
        ztmp=zth1a
        zth1a=zth2a
        zth2a=ztmp
        ztmp=zb1a
        zb1a=zb2a
        zb2a=ztmp
     end if
     !
     iter=iter+1
     if(iter.gt.10) then
        write(nonlin,*) ' %lots of iterations in bphip_bfnd'
     end if
     if(iter.gt.100) then
        write(nonlin,*) ' ?too many iterations in bphip_bfnd'
        ierr=1; return                 !  bug check
     end if
     if((zb1a-zbwant)*(zb2a-zbwant).lt.0.0_fp) then
        !  next guess -- trapzoidal rule -- 1st guess always goes here...

        zthlim1=zth1a
        zblim1=zb1a
        zthlim2=zth2a
        zblim2=zb2a
        !
        zth=((zb2a-zbwant)*zth1a-(zb1a-zbwant)*zth2a)/(zb2a-zb1a)
        !
     else
        !  both previous estimates to one side (and iter>1); try a newton step
        if(zdb.eq.0.0_fp) then
           zthnex=epsinv
        else
           zthnex=zth-(zb-zbwant)/zdb
        end if
        if((zthnex.lt.zthlim1).or.(zthnex.gt.zthlim2)) then
           !
           !  newton step went
           !  out of range; try to reset to a bracketing pair of guesses
           !
           if((zb-zbwant)*(zblim1-zbwant).lt.0.0_fp) then
              if(abs(zth1a-zthlim1).lt.abs(zth2a-zthlim1)) then
                 zth2a=zth1a
                 zb2a=zb1a
              end if
              zth1a=zthlim1
              zb1a=zblim1
           else
              if(abs(zth2a-zthlim2).lt.abs(zth1a-zthlim2)) then
                 zth1a=zth2a
                 zb1a=zb2a
              end if
              zth2a=zthlim2
              zb2a=zblim2
           end if
           goto 100
           !
        else
           !  accept the newton step
           zth=zthnex
        end if
     end if
     !  evaluate at this guess
     bithst(1)=zth
     call bintrp
     zbg=1.0_fp/bi1ovb(1)
     zb=1.0e-4_fp*zbg
     if(abs(zb-zbwant).lt.zbtol) goto 200
     zdb=-1.0e-4_fp*zbg*zbg*bibdts(1)
     !  do another guess
     if(abs(zb1a-zbwant).lt.abs(zb2a-zbwant)) then
        zb2a=zb
        zth2a=zth
     else
        zb1a=zb
        zth1a=zth
     end if
     if(abs(zth2a-zth1a).lt.zdthtol) goto 200
     goto 100
  end if
  !
200 continue
  !
  zthfnd=zth
  goto 1000
  !
900 continue
  !
  !  endpoint match:
  !
  ithfnd=ith_bphip(zthfnd,iath,inumth)
  !
  !  exit
  !
1000 continue
  return
end subroutine bphip_bfnd
!-------------------------
!
integer function ith_bphip(zth,iath,inumth)
  !
  !  find zone index, for passed theta value.
  !   ** side effect ** -- zth value is normalized to range
  !    (-pi,pi)
  !
  !  value returned is in range 1 to (inumth-1)
  !   inumth=# of zone bdys.
  !
  use nbi_com
  !
  !----------------
  !
  !============
  implicit none
  integer :: iath,inumth,ithzon0,ithfnd
  !============
  real(fp) :: zth,zthz1,zthz2
  !============
  if(zth.lt.-xpi) zth=zth+twopi
  if(zth.gt.xpi) zth=zth-twopi
  !
  ithzon0=(inumth-2)*(zth+xpi)/twopi
  ithzon0=max(0,min((inumth-2),ithzon0))
  !
  zthz1=datspl(iath+ithzon0)
  if(zthz1.gt.zth) then
     ithzon0=ithzon0-1
     if(ithzon0.lt.0) ithzon0=inumth-2
  end if
  zthz2=datspl(iath+ithzon0+1)
  if(zthz2.lt.zth) then
     ithzon0=ithzon0+1
     if(ithzon0.gt.(inumth-2)) ithzon0=0
  end if
  !
  ithfnd=ithzon0+1
  ith_bphip=ithfnd
  !
  return
end function ith_bphip
!
!------------------------------------------
subroutine bphip_phif(isurf,zth1,zb1,zthb,zp,zphib,iaflag,ierr)
  !
  use nbi_com
  !
  !  on the selected surface,
  !  find p(psi) and the poloidal arc length theta(b) and the
  !  toroidal arc length phi(b)
  !
  !  bphip_phip subroutine finds the other side of the
  !  orbit and p(psi); bphip_phib subroutine used to compute
  !  the phi(b) of this orbit.
  !
  !
  ! input:
  !============
  implicit none
  integer :: idir,ierr
  !============
  real(fp) :: zthb,zxi,zdb,zd2b,nbi_d1mach,zth1a,zth2a,ztol
  !============
  integer :: isurf                     ! surface index input
  real(fp) :: zth1                         ! theta angle input
  real(fp) :: zb1                          ! b at zth1 (output if 0.0)
  ! output:
  real(fp) :: zth2                         ! theta at other side of robit
  real(fp) :: zp                           ! p(psi) output
  real(fp) :: zphib                        ! phi(b) output
  integer :: iaflag                    ! bphip_phip asymptotic flag
  !
  !-------------------------------
  !  find direction
  !
  idir=0  ! this tells bphip_phip to compute idir as well
  !
  !  find other side of orbit
  !
  zxi=real(isurf,fp)/real(nphipc,fp)
  if(zb1.eq.0.0_fp) then
     call bphip_gb(zxi,zth1,zb1,zdb,zd2b)
  end if
  !
  !  compute p(psi) and find other bounce point of orbit
  !
  call bphip_phip(isurf,zth1,zb1,zth2,zp,idir,iaflag,ierr)
  if(ierr.ne.0) return
  !
  !  ok now compute phib
  !
  if(iaflag.eq.0) then
     ztol=10.0_fp*nbi_d1mach(4)*xpi
  else
     ztol=0.0_fp
  end if
  call bphip_phib(isurf,zth1,zth2,zphib,idir,ztol,ierr)
  if(ierr.ne.0) return
  !
  !  standardize scan direction locally
  !
  if(idir.eq.1) then
     zth1a=zth1
     zth2a=zth2
  else
     zth1a=zth2
     zth2a=zth1
  end if
  !
  !  check for crossing branch cut
  !
  if(zth2a.lt.zth1a+ztol) then
     zth2a=zth2a+twopi
  end if
  !
  !  poloidal arclength
  !
  zthb=idir*(zth2a-zth1a)
  !
  return
end subroutine bphip_phif
!------------------------------------------
subroutine bphip_phib(isurf,zth1,zth2,zphib,idir,ztol,ierr)
  !
  use nbi_com
  !
  !  on the specified surface, find the toroidal arc length phi(b)
  !  along a field line spanning poloidal arc from zth1 to zth2.
  !
  !  scan in the direction specified by idir, unless idir=0.
  !
  !  if idir=0:
  !  note there are actually two arcs spanning this range; choose the one
  !  on the small field side, determined by looking at db/dth at zth1;
  !  if db/dth > 0, go in -theta direction from zth1 otherwise go in
  !  +theta direction.
  !
  !  the output arclength is always positive, regardless of the order
  !  of (zth1,zth2).  the output idir indicates the direction followed
  !  from zth1 to define the arc.
  !
  !  input:
  !
  !============
  implicit none
  integer :: iaoff,imax,iammx0,iammn0,jx,ibranch,ilinx,ilinth,ier
  integer :: ii,jj,ier1,ierr
  !============
  real(fp) :: zf1d,zxi,zx0,zh,ziot,zq,zth1a
  real(fp) :: zth2a,zphi2add,zdxi,zdthi,zths1,znu1,zphi1,zths2,znu2
  real(fp) :: zphi2,ztol
  !============
  integer :: isurf                     ! surface index
  real(fp) :: zth1                         ! zth1 in range (-pi,pi)
  real(fp) :: zth2                         ! zth2 in range (-pi,pi)
  !
  !  output:
  !
  real(fp) :: zphib                        ! toroidal arc length along b
  integer :: idir                      ! +1 if arc starts in +theta
  ! direction from zth1, -1 o.w.
  !
  integer :: iselth(6)
  !
  data iselth/1,0,0,0,0,0/
  !
  ierr = 0
  !
  !  addresses; determine direction of orbit reflection, from
  !  zth1
  !
  if(idir.eq.0) then
     call bphip_dir(isurf,zth1,iaoff,imax,iammx0,iammn0,idir)
  end if
  !
  !  local q value --
  !
  zxi=real(isurf,fp)/real(nphipc,fp)
  !
  jx=int(zxi/dxi(lcentr,2)) + 1
  zx0=(jx-1)*dxi(lcentr,2)
  zh=zxi-zx0
  ziot=zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zq=1.0_fp/ziot
  !
  ibranch=0
  !
  !  standardize scan direction locally
  !
  if(idir.eq.1) then
     zth1a=zth1
     zth2a=zth2
  else
     zth1a=zth2
     zth2a=zth1
  end if
  !
  !  check for crossing branch cut
  !
  if(zth2a.lt.zth1a+ztol) then
     ibranch=1
     zphi2add=twopi*zq    ! fieldline lengh for one poloidal circuit
  else
     ibranch=0
     zphi2add=0.0_fp
  end if
  !
  !  evaluate bicubic splines to get phi location
  !
  ilinx=1
  ilinth=0
  !
  ier=0
  !
  call bcspevxy(zxi,zth1a,xizb(lcentr),nzp1,datspl(lsptheta), &
       nthspl,ilinx,ilinth,ii,jj,zdxi,zdthi,ier1)
  ier=max(ier,ier1)
  !
  call bcspevfn(iselth,1,1,zths1,ii,jj,zdxi,zdthi, &
       datspl(lspthets),nzp1,nthspl)
  !
  call bcspevfn(iselth,1,1,znu1,ii,jj,zdxi,zdthi, &
       datspl(lspphilb),nzp1,nthspl)
  !
  zphi1=zq*(xpi+zths1)-znu1
  !
  call bcspevxy(zxi,zth2a,xizb(lcentr),nzp1,datspl(lsptheta), &
       nthspl,ilinx,ilinth,ii,jj,zdxi,zdthi,ier1)
  ier=max(ier,ier1)
  !
  call bcspevfn(iselth,1,1,zths2,ii,jj,zdxi,zdthi, &
       datspl(lspthets),nzp1,nthspl)
  !
  call bcspevfn(iselth,1,1,znu2,ii,jj,zdxi,zdthi, &
       datspl(lspphilb),nzp1,nthspl)
  !
  zphi2=zphi2add+zq*(xpi+zths2)-znu2
  !
  if(ier.gt.0) then
     write(nonlin,*) 'bphip_phib spline eval error.'
     ierr=1
  else
     zphib=zphi2-zphi1
  end if
  !
  return
end subroutine bphip_phib
!--------------------------------------------
subroutine bphip_phip(isurf0,zth1,zb1,zth2,zp,idir,iaflag,ierr)
  use bphip_mod
  !
  !  on the specified surface, find information related to the toroidal
  !  precession of a fast ion trapped orbit.
  !
  !  input isurf0 -- surface index (datspl scheme)
  !        zth1  -- start angle -- reflection point of trapped orbit
  !         zb1  -- b field at the reflection point
  !
  ! output zth2  -- stop angle -- reflection point at other end of
  !                 trapped orbit (where b=zb1 also)
  !
  !          zp  -- p(psi) integral, eqn 11,
  !                 r. b. white, r. j. goldston, "ripple-induced energetic
  !                 particle loss in tokamaks"
  !
  !                 toroidal precession angle is phi(p) = rho*p(psi), where
  !                 rho = v/omega0/r0, omega0 = cyc.freq on axis, r0=
  !                 major radius of axis
  !
  !                 p(psi) is indep. of energy; all energy dependence of
  !                 phi(p) comes via rho.
  !
  !                 caution -- expression ignores exb (rotation related)
  !                     effects!
  !
  !        idir  -- direction from zth1 of orbit
  !
  !        iaflag - if non-zero this is the index of the local b minimum
  !                 from which an asymptotic value of p(psi) was taken
  !                 this is done if the orbit is very deeply trapped
  !                 and there is insufficient variation in b to do the
  !                 numerical p(psi) integral.
  !
  !  the output integral value p(psi) is always evaluated by integrating
  !  from a smaller theta value to a larger theta value, regardless of
  !  the ordering of zth1 and zth2.  (zth1 input, zth2 found here).
  !  i.e. a call bphip_phip(isurf,zth2,...) (were zth2 known)
  !  would return the same positive value of p(psi) as a call
  !  bphip_phip(isurf,zth1,...); only idir would be reversed for the two
  !  situations.
  !
  !-------
  !
  use nbi_com
  !
  !-------
  !  for p(psi) integrand
  !
  !============
  implicit none
  integer :: idir,iaflag,isurf0,isurf,iamin,idum1,idum2,ith2
  integer :: jx,inumi,it,ierr
  !============
  real(fp) :: zth1,zb1,zth2,zp
  real(fp) :: zf1d,zthbmin
  real(fp) :: zthbmax,zthbmax2,zbmax2,zbmaxb,zbdeep,zth1a,zth2a
  real(fp) :: nbi_d1mach,zepsth,zdthmin,zthav,zx0,zh,ziot,zdiot,zcv
  real(fp) :: zbdum,zsum,zf,zdth,zterm0,ztermm,ztermp,zd2
  real(fp) :: ziterm,zdf1d
  !
  ierr = 0
  !
  isurf=isurf0
  !
  !  determine direction of orbit reflection, and find zth2, given
  !  zth1
  !
  !   ..first find the search range; btw which two singularities to search
  !
  idir=0  ! this tells bphip_mmx to compute idir as well
  !
  call bphip_mmx(isurf,zth1,zb1,idir, &
       zthbmin,zbmin,iamin,idum1,zthbmax,zbmax,idum2,zthbmax2,zbmax2)
  !
  !   ..then find zth2 where b=zb1 occurs on the other side of the orbit
  !
  call bphip_bfnd(isurf,zb1,zthbmin,zbmin,zthbmax,zbmax, &
       zth2,ith2,1,1,ierr)
  if(ierr.ne.0) return
  !
  !  determine if orbit is very deeply trapped (bmax is < 1.004*bmin)
  !  adjustment factor determined empirically -- due to singularity there
  !  start to be problems with the evaluation of the integral
  !
  iaflag=0
  zbmaxb=min(zbmax,zbmax2)          ! bounds the well at iamin
  zbdeep=min(1.004_fp*zbmin, &
       (0.6666667_fp*zbmin+0.3333333_fp*zbmaxb))
  if(zb1.le.zbdeep) then
     iaflag=iamin
     zp=datspl(iamin+2)             ! asymptotic p(psi) value
     goto 1000
  end if
  !
  !  ok numeric integration will be done
  !
  !  at this point, both zth1 and zth2 are in the range (-pi,pi); zth2
  !  less than zth1 is possible, depending on idir
  !
  !  before doing the integration we want to standardize the ordering of
  !  the limits; this may include crossing the branch cut at +pi.
  !
  if(idir.eq.1) then
     zth1a=zth1
     zth2a=zth2
  else
     zth1a=zth2
     zth2a=zth1
  end if
  !
  !  check for crossing branch cut
  !
  zeps=nbi_d1mach(4)    ! machine precision
  !
  !  allow for some slop here.  if zth1a and zth2a were really close we
  !  would have taken the deep-trap exit (iaflag=iamin) above.
  !
  zepsth=zeps*10.0_fp*xpi
  if(zth2a.lt.zth1a+zepsth) then
     zth2a=zth2a+twopi
  end if
  !
  !  if limits are too close we are near a b minimum; spread the limits
  !  slightly to avoid -> 0/0 situation in the integration
  !
  !  this code may never be entered, due to use of asymptotic value
  !  in deep trapping limit
  !
  zdthmin=1.0e-3_fp
  if((zth2a-zth1a).lt.zdthmin) then
     zthav=0.5_fp*(zth1a+zth2a)
     zth1a=zthav-0.5_fp*zdthmin
     zth2a=zthav+0.5_fp*zdthmin
  end if
  !
  !  ok limits are set.
  !  evaluate quantities need for numeric evaluation of integral
  !
  zbmax=zb1
  !
  zxi=real(isurf)/real(nphipc)
  !
  jx=int(zxi/dxi(lcentr,2)) + 1
  zx0=real(jx-1)*dxi(lcentr,2)
  zh=zxi-zx0
  !
  ziot=zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zq=1.0_fp/ziot
  zdiot=zdf1d(zh,spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zdqdxi=-zq*zq*zdiot
  !
  zcv=1.0e-4_fp !t/gauss
  !
  zg=zcv*zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  zdgdxi=zcv*zdf1d(zh,splg(2,jx),splg(3,jx),splg(4,jx))
  !
  zi=zcv*zf1d(zh,splimhd(1,jx), &
       splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  zdidxi=zcv*zdf1d(zh,splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  !
  zthimin=zth1a
  zthimax=zth2a
  !
  !  b derivatives at end pts
  call bphip_gb(zxi,zthimin,zbdum,zmnbth1,zmnbth2)
  call bphip_gb(zxi,zthimax,zbdum,zmxbth1,zmxbth2)
  !
  zsum=0.0_fp
  !
  !  theta zone ctrs
  !
  inumi=intot
  do it=1,inumi
     zf=(real(it)-0.5_fp)/real(inumi)
     ztha(it)=zthimin+zf*(zthimax-zthimin)
  end do
  !
  call bphip_pig2(ierr)
  if(ierr.ne.0) return
  !
  zdth=(zthimax-zthimin)/inumi
  !
  !  for non edge points, a numerically evaluated 2nd derivative based
  !  correction is applied... this helps the accuracy of the integral in
  !  the singular edge region
  !
  ziterm=zterm(1)+zterm(inumi)
  !
  do it=2,inumi-1
     !  integrand
     zterm0=zterm(it)
     ztermm=zterm(it-1)
     ztermp=zterm(it+1)
     !  2nd order correction
     zd2=(ztermp+ztermm-2*zterm0)/24.0_fp
     !  sums:     zdth not yet factored in
     ziterm=ziterm+zterm0+zd2
  end do
  !
  zp=ziterm
  !
  zp=zp*zdth
  !
1000 continue
  return
end subroutine bphip_phip
!--------------------------------------------
real(fp) function bphip_plim(zxi,zthmin,ierr)
  !
  use nbi_com
  !
  !  evaluate deeply trapped limit value of p(psi) integral,
  !  at the indicated local minimum bmin @thetamin on the zxi flux
  !  surface.
  !
  !  this is a limiting expression for the same integrand evaluated
  !  in subroutine bphip_pig2.
  !
  !  using calculus (cf notes dmc 1 aug 1996) it can be shown that the
  !  limiting expression is:
  !
  !  p(psi) = pi**2*q*bzxr/(zxi*tflux) * dtheta(gwb)/dtheta(transp)
  !           *(i + g*q)*(db/dxi)/(2*bmin**2*sqrt((d2b/dtheta2)/2*bmin))
  !
  !  the extra factor of pi comes from the analytic evaluation of
  !    2 * integral[th = 0 to delth](dth/sqrt(delth**2 - th**2))
  !  which appears when the original p(psi) expression is taken to the
  !  limit.
  !
  !============
  implicit none
  integer :: ier1,jx,ierr
  !============
  real(fp) :: zf1d,zdthcon,zx0,zh,ziot,zq,zcv,zg
  real(fp) :: zi,zbmin,zbg,zdbdxi,zd2bdt2,zfac,zarg,zans
  !============
  real(fp) :: zxi                          ! flux surface (sqrt norm tor.flx)
  real(fp) :: zthmin                       ! theta of bmin
  !
  !  the caller is presumed to have found the correct theta such that
  !  b(zthmin) is a local minimum on this zxi flux surface.
  !
  integer :: ideriv(6),ilinx,ilinth
  !
  data ilinx/1/
  data ilinth/1/
  !
  data ideriv/0,0,1,0,0,0/
  !
  !--------------------------------------------
  !
  ierr=0
  !
  !  dtheta(gwb)/dtheta(transp) conversion factor
  !
  call bcspeval(zxi,zthmin,ideriv,zdthcon, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       ilinx,ilinth, &
       datspl(lspthets),nzp1,ier1)
  if(ier1.ne.0) then
     write(nonlin,*) 'bphip_plim: spline eval error.'
     ierr=1
     return
  end if
  !
  !  evaluate q, g, i
  !
  jx=int(zxi/dxi(lcentr,2)) + 1
  zx0=real(jx-1)*dxi(lcentr,2)
  zh=zxi-zx0
  !
  ziot=zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zq=1.0_fp/ziot
  !
  zcv=1.0e-4_fp !t/gauss
  !
  zg=zcv*zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  !
  zi=zcv*zf1d(zh,splimhd(1,jx), &
       splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  !
  !  evaluate b and derivatives
  !
  bixi(1)=zxi
  bithst(1)=zthmin
  !
  call bintrp
  !
  zbmin=bimodb(1)
  zbg=1.0e4_fp*zbmin                   ! b in gauss
  zdbdxi=-1.0e-4_fp*zbg*zbg*bibdxi(1)  ! db/dxi @bmin
  !  d2b/dth2:
  zd2bdt2=1.0e-4_fp*zbg*zbg* &
       (2.0_fp*zbg*bibdts(1)*bibdts(1) - bib2ts(1))
  !
  !  ok:
  !
  !dbg      zdthcon=1.0
  zfac=zdthcon*zq*xpi*bzxr/(zxi*1.0e4_fp*tflux) ! cm**-1
  !
  zarg=max(1.0e-8_fp,(zd2bdt2/(2.0_fp*zbmin)))
  zans=-xpi*zfac*(zi+zg*zq)*zdbdxi/(2.0_fp*zbmin*zbmin*sqrt(zarg))
  !
  bphip_plim=zans
  !
  return
end function bphip_plim
!
!
!--------------------------------------------
!
subroutine bphip_pig2(ier)
  use bphip_mod
  !
  !  phi(p) p(psi) integrand
  !  eqn 11, "ripple-induced energetic particle loss in tokamaks",
  !  white, goldston, redi, budny
  !
  !  dmc -- jul-sep 1996 -- trying various integration methods
  !    and general testing...
  !
  !  evaluate at an array of theta points
  !
  !  dnc sep 1996 -- part of integrand is singular at endpoint with
  !   1/sqrt(delta(theta)) type singularity; trying integration by
  !   parts method to avoid problems due to this.
  !
  !  can transform by parts int[ dtheta*h(theta)/sqrt(delta(theta)) ] to
  !     the form 2k -int[ dtheta*phi(theta)*h'(theta)*sqrt(delta(theta)) ]
  !     where phi(theta)= -1 if theta < theta(av), +1 if theta > theta(av),
  !       h(theta)=a(theta)*b(theta)*c(theta)
  !
  !       a(theta)=sqrt( delta(theta) / (1 - b/bmax) )
  !       b(theta)=[ 1/b**2 - 1/(2*b*bmax) ]
  !       c(theta)=[ db/dxi*dthgwb/dtheta - db/dtheta*dthgwb/dxi ]
  !
  !    midpoint:
  !       theta(av)=0.5*(theta(max)+theta(min))
  !    distance from nearest endpoint:
  !       delta(theta)=
  !         0.5*(theta(max)-theta(min)) - (theta - theta(av))
  !             if theta > theta(av),
  !         0.5*(theta(max)-theta(min)) - (theta(av) - theta)
  !             if theta < theta(av)
  !
  !       bmax = b(theta(max)) = b(theta(min))
  !          1 - b/bmax = 0 (singularity) at endpoints only.
  !       b(xi,theta) is known      ! b field
  !       thgwb(xi,theta) is known  ! transformation to gwb coordinates
  !
  !       k = h(theta(av))*sqrt(0.5*(theta(max)-theta(min)))
  !
  !   a(theta) -> 0/0 at endpoints but can be shown to have well behaved
  !   limits... so, endpoint singularities are handled.
  !
  !     derivation consists of breaking the original integral into
  !     two pieces each with one singularity:  from theta(av) to theta(max)
  !     and from theta(min) to theta(av).  an integration-by-parts is then
  !     done to remove the singularity from under the integral sign; this
  !     yields a constant part k for each piece that is of equal magnitude
  !     giving the factor of 2.
  !
  !     h'(theta) =a'*b*c + a*b'*c + a*b*c'
  !
  !-------------------------
  !
  use nbi_com
  !
  !  common block set up in caller (bphip_phip)
  !
  !
  !============
  implicit none
  integer :: ier,ier1,i,ii,jj
  !============
  real(fp) :: zthavg,zthrng,ztharg,zdelth,zsqrtdel,zdthdxi,zdthcon,zb,zbg
  real(fp) :: zdbdxi,zdbdth,zmbb,zafa,zbfa,zcfa,zha,z2k,z2kth
  real(fp) :: zdpsidx,zfac,ztheta,zdelpr,zphi,zdxi,zdthi,zd2thdt2
  real(fp) :: zd2thdxt,zdnu2,zd2bdxth,zd2bdth2,zrhr,zthdif,zhprime
  real(fp) :: ziterm
  !-----------------
  !  debugging
  !
  real(fp) :: zaf(intot),zdaf(intot)      ! a(theta) and deriv.
  real(fp) :: zbf(intot),zdbf(intot)      ! b(theta) and deriv.
  real(fp) :: zcf(intot),zdcf(intot)      ! c(theta) and deriv.
  !
  !-----------------
  !  local...
  !
  integer :: ilinx,ilinth,iselec0(6),iselect(6),ideriv2(6)
  !
  data ilinx/1/
  data ilinth/1/
  !
  data iselec0/0,1,1,0,0,0/
  data iselect/0,1,1,0,1,1/
  data ideriv2/0,0,0,0,0,1/
  !
  real(fp) :: zderiv(6)
  !
  !-----------------
  !
  ier=0
  zaf=0.0_fp
  zbf=0.0_fp
  zcf=0.0_fp
  zdaf=0.0_fp
  zdbf=0.0_fp
  zdcf=0.0_fp
  !
  zthavg=0.5_fp*(zthimin+zthimax)
  zthrng=0.5_fp*(zthimax-zthimin)
  !
  ztharg=zthavg
  if(ztharg.gt.xpi) ztharg=ztharg-twopi
  if(ztharg.lt.-xpi) ztharg=ztharg+twopi
  !
  !  evaluate h(zthavg)
  !
  zdelth=zthrng
  zsqrtdel=sqrt(zdelth)
  !
  call bcspeval(zxi,ztharg,iselec0,zderiv, &
       xizb(lcentr),nzp1,datspl(lsptheta),nthspl,ilinx,ilinth, &
       datspl(lspthets),nzp1,ier1)
  ier=max(ier,ier1)
  !
  zdthdxi=zderiv(1)                 ! dth(gwb)/dxi
  zdthcon=zderiv(2)                 ! dth(gwb)/dth(transp)
  !
  !  b and db/dxi and db/dth information
  !
  bixi(1)=zxi
  bithst(1)=ztharg
  !
  call bintrp
  !
  zb=bimodb(1)
  zbg=1.0e4_fp*bimodb(1)
  zdbdxi=-1.0e-4_fp*zbg*zbg*bibdxi(1)
  zdbdth=-1.0e-4_fp*zbg*zbg*bibdts(1)
  !
  !  sqrt(1-b/bmax)/b = roscoe's rho/rho||
  !
  zmbb=max(zeps,(1.0_fp-zb/zbmax))
  !
  !  a,b,c functions
  zafa=sqrt(zdelth/zmbb)
  zbfa=(1.0_fp/(zb*zb) - 1.0_fp/(2.0_fp*zb*zbmax))
  zcfa=(zdbdxi*zdthcon-zdbdth*zdthdxi)
  
  !  h(theta(av))
  zha=zafa*zbfa*zcfa
  !
  z2k=2.0_fp*(zi+zg*zq)*zha*zsqrtdel
  z2kth=z2k/(zthimax-zthimin)       ! add into integrand
  !
  !------------------------
  !
  !  common factor (includes dxi/dpsi = q*pi/(zxi*tflux) to transform
  !  d/dxi to d/dpsi in terms)
  !
  zdpsidx=(zxi*1.e4_fp*tflux)/(zq*xpi)
  zfac=bzxr/zdpsidx
  !
  do i=1,intot
     ztheta=ztha(i)
     ztharg=ztheta
     if(ztharg.gt.xpi) ztharg=ztharg-twopi
     if(ztharg.lt.-xpi) ztharg=ztharg+twopi
     !
     zdelth=zthrng-abs(ztheta-zthavg)
     zsqrtdel=sqrt(zdelth)
     !
     if(ztheta.gt.zthavg) then
        zdelpr=-1.0_fp  ! delta' if ztheta > zthavg
        zphi=1.0_fp
     else
        zdelpr=1.0_fp  ! delta' if ztheta < zthavg
        zphi=-1.0_fp
     end if
     !
     !  dtheta(gwb)/dtheta(transp) conversion factor
     !
     call bcspevxy(zxi,ztharg, &
          xizb(lcentr),nzp1,datspl(lsptheta),nthspl,ilinx,ilinth, &
          ii,jj,zdxi,zdthi,ier1)
     ier=max(ier,ier1)
     !
     call bcspevfn(iselect,1,1,zderiv,ii,jj,zdxi,zdthi, &
          datspl(lsptheta),nzp1,nthspl)
     !
     zdthdxi=zderiv(1)  ! dth(gwb)/dxi
     zdthcon=zderiv(2)  ! dth(gwb)/dth(transp)
     zd2thdt2=zderiv(3) ! d2th(gwb)/dth2(transp)
     zd2thdxt=zderiv(4) ! d2th(gwb)/dxi/dth(transp)
     !
     !  d2nu/dxi/dtheta(transp)
     !
     call bcspevfn(ideriv2,1,1,zdnu2,ii,jj,zdxi,zdthi, &
          datspl(lspphilb),nzp1,nthspl)
     ier=max(ier,ier1)
     !
     !  b and db/dxi and db/dth information
     !
     bixi(1)=zxi
     bithst(1)=ztharg
     !
     call bintrp
     !
     zb=bimodb(1)
     zbg=1.0e4_fp*bimodb(1)
     zdbdxi=-1.0e-4_fp*zbg*zbg*bibdxi(1)
     zdbdth=-1.0e-4_fp*zbg*zbg*bibdts(1)
     zd2bdxth = 1.0e-4_fp*zbg*zbg* &
          (2.0_fp*zbg*bibdxi(1)*bibdts(1)-bibtsx(1))
     zd2bdth2 = 1.0e-4_fp*zbg*zbg* &
          (2.0_fp*zbg*bibdts(1)*bibdts(1)-bib2ts(1))
     !
     !  sqrt(1-b/bmax)/b = roscoe's rho/rho||
     !
     zmbb=max(zeps,(1.0_fp-zb/zbmax))
     zrhr=sqrt(zmbb)/zb
     !
     !  if too close to singularity, use expansion for 1-b/bmax
     if(zmbb.lt.10.0_fp*zeps) then
        if(ztheta.gt.zthavg) then
           zthdif=ztheta-zthimax
           zmbb=-zthdif*(zmxbth1+0.5_fp*zmxbth2*zthdif)/zbmax
        else
           zthdif=ztheta-zthimin
           zmbb=-zthdif*(zmnbth1+0.5_fp*zmnbth2*zthdif)/zbmax
        end if
     end if
     zmbb=max(1.0e-8_fp,zmbb)
     !
     !  a,b,c functions
     zaf(i)=sqrt(zdelth/zmbb)
     zbf(i)=(1.0_fp/(zb*zb) - 1.0_fp/(2.0_fp*zb*zbmax))
     zcf(i)=(zdbdxi*zdthcon-zdbdth*zdthdxi)
     !
     !  a,b,c derivatives
     zdaf(i)=(zdelpr/zaf(i) + zaf(i)*zdbdth/zbmax)/(2.0_fp*zmbb)
     zdbf(i)=(-2.0_fp/zb + 1.0_fp/(2.0_fp*zbmax))*zdbdth/(zb*zb)
     zdcf(i)=(zd2bdxth*zdthcon + zdbdxi*zd2thdt2) - &
          (zd2bdth2*zdthdxi + zdbdth*zd2thdxt)
     !
     !  h', from h=a*b*c
     !
     zhprime=zdaf(i)*zbf(i)*zcf(i) &
          +zaf(i)*zdbf(i)*zcf(i) &
          +zaf(i)*zbf(i)*zdcf(i)
     !
     !  terms  -- all have units cm (zi,zg & derivs are t*cm)
     !
     !  grad(b) + curvature(b) term
     zterm(i) = -2.0_fp*(zphi*(zi+zg*zq)*zhprime*zsqrtdel + z2kth)
     !
     !  curl(b) term  (as per roscoe)
     zterm(i) = zterm(i) + zrhr*(zdidxi+zq*zdgdxi)*zdthcon
     !
     !  sheer term  (as per roscoe with d2nu/dnu/dxi term added)
     zterm(i) = zterm(i) + zg*zrhr*(zdqdxi*zdthcon-zdnu2)
     !
     zterm(i)=zfac*zterm(i)
     !
  end do
  !
  ! debug
  !dbg      call grafx2(ztha,zaf,zdaf,intot,'radians','a',
  !dbg     >     runid,'bphip_pig2 debug plot','function and derivative')
  !dbg      call grafx2(ztha,zbf,zdbf,intot,'radians','b',
  !dbg     >     runid,'bphip_pig2 debug plot','function and derivative')
  !dbg      call grafx2(ztha,zcf,zdcf,intot,'radians','c',
  !dbg     >     runid,'bphip_pig2 debug plot','function and derivative')
  !
  if(ier.gt.0) write(nonlin,*) ' ?bphip_pig2: spline eval error.'
  !
  return
end subroutine bphip_pig2
!
!--------------------------------------------
!
subroutine bphip_dir(isurf,zth1,iaoff,imax,iammx0,iammn0,idir)
  !
  !  on surface isurf from starting angle zth1 (btw -pi and pi), find
  !  the direction to the small field side.  also return some addressing
  !  information re. singularity data in datspl(...).
  !
  !  idir=1 if db/dth < 0 (+theta => smaller field)
  !  idir=-1 if db/dth > 0
  !
  !  determination is made by looking at nearby singularities in b(theta)
  !  on the indicated surface.
  !
  !  input:  isurf, the surface
  !           zth1, the start angle
  !
  ! output:  iaoff, imax, iammx0, iammn0
  !          addressing info; no. of maxima = no. of minima (since b is
  !          periodic in theta).
  !
  !          idir -- the direction as defined above.
  !
  !-------
  !
  use nbi_com
  !
  !-------
  !
  !============
  implicit none
  integer :: iaoff,imax,iammx0,iammn0,idir,isurf,im,ia0
  !============
  real(fp) :: zth1,zthsing,zth
  !============
  iaoff=ndatspl(isurf)-1
  !
  imax=ndatspl(iaoff+5)
  iammx0=ndatspl(iaoff+6)
  iammn0=iammx0+3*imax
  !
  !
  !  find nearest singularity at theta > zth1
  !   if nearest singularity is maxima, idir=-1
  !   if nearest singularity is minima, idir=+1
  !
  zthsing=epsinv
  !
  !  check maxima
  !
  do im=1,imax
     ia0=iammx0+3*(im-1)
     zth=datspl(ia0)
     if(zth.lt.zth1) zth=zth+twopi
     if(zth.lt.zthsing) then
        zthsing=zth
        idir=-1
     end if
  end do
  !
  !  check minima
  !
  do im=1,imax
     ia0=iammn0+3*(im-1)
     zth=datspl(ia0)
     if(zth.lt.zth1) zth=zth+twopi
     if(zth.lt.zthsing) then
        zthsing=zth
        idir=1
     end if
  end do
  !
  return
end subroutine bphip_dir
!------------------------------------------
subroutine bphip_minp(ztha,zpha,ia,inum,iaoff)
  !
  use nbi_com
  !
  !  patch phi(p) in deeply trapped regions
  !
  !============
  implicit none
  integer :: inum,imax,inslim,isinga,ising0,ising1,iskip,it,is
  integer :: inums,iaa,it2,iflag,iflag2,it3,itx,ict,itx1
  !============
  real(fp) :: zthtol,nbi_d1mach,zdthl,zthbmin,zthglb,zthlub,zthsing
  real(fp) :: zthtes,zth2,za2,zth3,zth,zthp
  !============
  real(fp) :: ztha(inum)                   ! theta array
  real(fp) :: zpha(inum)                   ! phi(p) vs. theta
  integer :: ia(inum)                  ! deep trap marks
  !
  integer :: iaoff                     ! ptr into ndatspl data
  !
  !  asymptotic values are known for most-deeply-trapped orbits;
  !  numerical evaluation of p(psi) is unreliable in that vicinity.
  !  so, use the asymptotic value and a cubic fit in these regions
  !  which are in the vicinity of each local minimum in b.
  !
  !  the deeply trapped regions {j} are marked by ia({j}).ne.0
  !  zpha({j}) currently contains the asymptotic value.
  !
  !  at the minimum, p(psi) is known and d/dtheta (p(psi)) is zero.
  !
  !  find two good points beyond the deeply trapped region and fit a
  !  cubic.
  !
  !  if two good points are not available, do nothing (leave the asymptotic
  !  value alone; the resulting function may not be smooth).
  !
  !  the ztha array spans the region (-pi,pi) with the endpoints being
  !  -pi,pi.  p(psi) is periodic but careful:  the endpoints at +/-pi are
  !  duplicates.
  !
  !  the ztha grid spacing is assumed to be monotonic increasing and
  !  fairly evenly spaced
  !
  !-------------------
  !
  real(fp) :: za,zb,zd,zx,zx2,zx3,zx22,zx23,zx32,zx33
  real(fp) :: zp,zppp,zpp2,zpp3
  !-------------------
  !
  imax=ndatspl(iaoff+5)
  inslim=4*imax
  isinga=iaoff+10
  ising0=ndatspl(iaoff+6)
  ising1=ising0+3*2*imax
  !
  zthtol=10.0_fp*nbi_d1mach(4)*xpi
  !
  !  outer loop:  look for deeply trapped regions
  !
  zdthl=0.5_fp*min((ztha(2)-ztha(1)),(ztha(inum)-ztha(inum-1)))
  iskip=0
  do it=1,inum
     if((it.gt.iskip).and.(ia(it).ne.0)) then
        !
        !  the theta value and asymptotic p(psi) value;
        !  find bounding singularities around the minimum
        !
        zthbmin=datspl(ia(it))
        zppp=zpha(it)
        !
        zthglb=-epsinv              ! greatest lower bound
        zthlub=epsinv               ! least upper bound
        do is=1,inslim
           inums=ndatspl(isinga+is-1)
           if(inums.le.2*imax) then
              iaa=ising0+3*(inums-1)
           else
              iaa=ising1+(inums-2*imax)-1
           end if
           zthsing=datspl(iaa)
           if(is.eq.1) zthsing=zthsing-twopi
           if(is.eq.inslim) zthsing=zthsing+twopi
           !
           if(abs(zthsing-zthbmin).gt.zthtol) then
              if(zthsing.lt.zthbmin) then
                 zthglb=max(zthglb,zthsing)
              else
                 zthlub=min(zthlub,zthsing)
              end if
           end if
        end do
        !
        it2=it
10      continue
        !
        !  find small theta edge of trapped region
        !
        it2=it2-1
        if(it2.eq.0) it2=inum-1
        if(ia(it2).ne.0) then
           goto 10
        else
           !
           !  found end of deep trapped region -- ok?
           !
           zthtes=ztha(it2)
           if(zthtes.gt.zthbmin) zthtes=zthtes-twopi
           if(zthtes.lt.zthglb+zthtol) then
              iflag=0               ! past or too close to singularity
              iflag2=0
           else
              iflag2=1              ! there is one good point
              !
              it3=it2-1
              if(it3.eq.0) it3=inum-1
              if(ia(it3).ne.0) then
                 iflag=0            ! not enough good points
              else
                 zthtes=ztha(it3)
                 if(zthtes.gt.zthbmin) zthtes=zthtes-twopi
                 if(zthtes.lt.zthglb+zthtol) then
                    iflag=0         ! past or too close to singularity
                 else
                    iflag=1         ! ok
                 end if
              end if
           end if
           !
           zth2=ztha(it2)
           if(zth2.gt.zthbmin) zth2=zth2-twopi
           if(iflag.eq.1) then
              zx2=zthbmin-zth2
              zpp2=zpha(it2)
              za2=(zpp2-zppp)/(zx2*zx2)
              !
              zth3=ztha(it3)
              if(zth3.gt.zthbmin) zth3=zth3-twopi
              zx3=zthbmin-zth3
              zpp3=zpha(it3)
              !
              !  cubic for small theta side
              !
              zx22=zx2*zx2
              zx23=zx22*zx2
              zx32=zx3*zx3
              zx33=zx32*zx3
              !
              zd=zx22*zx32*(zx3-zx2)
              !
              zb=(zx22*(zpp3-zppp)-zx32*(zpp2-zppp))/zd
              za=(zx33*(zpp2-zppp)-zx23*(zpp3-zppp))/zd
              !
           end if
           !
           !  zp = zppp + zx*zx*(za+zx*zb)
           !
           itx=it2
           zth=zth2
           !
           ict=0
20         continue
           !
           zthp=zth
           itx=itx+1
           if(itx.gt.inum) itx=1
           zth=ztha(itx)
           if(zth.gt.zthbmin) zth=zth-twopi
           if(zth.gt.(zthp-0.1_fp*zdthl)) then
              !  zth will drop by ~2pi when zthbmin is crossed
              if(iflag.eq.1) then
                 zx=zthbmin-zth
                 zp = zppp + zx*zx*(za+zx*zb)
                 zpha(itx)=zp
                 ia(itx)=0
              end if
              ict=ict+1
              if(ict.eq.1) itx1=itx
              goto 20  ! haven't reached zthbmin yet
           end if
           !
           if((iflag.eq.0).and.(iflag2.eq.1).and.(ict.gt.1)) then
              zpha(itx1)=0.9_fp*zpha(itx1)+0.1_fp*zpha(it2)
           end if
           !
        end if  ! end of small side region was found
        !--------
        !  now the large theta side
        !
        it2=itx
30      continue
        !
        !  find large theta edge of trapped region
        !
        it2=it2+1
        if(it2.gt.inum) it2=2
        if(ia(it2).ne.0) then
           goto 30
        else
           !
           !  found end of deep trapped region
           !
           iskip=it2                ! for the main loop
           !
           zthtes=ztha(it2)
           if(zthtes.lt.zthbmin) zthtes=zthtes+twopi
           if(zthtes.gt.zthlub-zthtol) then
              iflag=0               ! past or too close to singularity
              iflag2=0
           else
              iflag2=1              ! there is one good point
              !
              it3=it2+1
              if(it3.gt.inum) it3=2
              if(ia(it3).ne.0) then
                 iflag=0            ! not enough good points
              else
                 zthtes=ztha(it3)
                 if(zthtes.lt.zthbmin) zthtes=zthtes+twopi
                 if(zthtes.gt.zthlub-zthtol) then
                    iflag=0         ! past or too close to singularity
                 else
                    iflag=1         ! ok
                 end if
              end if
           end if
           !
           zth2=ztha(it2)
           if(zth2.lt.zthbmin) zth2=zth2+twopi
           if(iflag.eq.1) then
              zx2=zth2-zthbmin
              zpp2=zpha(it2)
              za2=(zpp2-zppp)/(zx2*zx2)
              !
              zth3=ztha(it3)
              if(zth3.lt.zthbmin) zth3=zth3+twopi
              zx3=zth3-zthbmin
              zpp3=zpha(it3)
              !
              !  cubic for large theta side
              !
              zx22=zx2*zx2
              zx23=zx22*zx2
              zx32=zx3*zx3
              zx33=zx32*zx3
              !
              zd=zx22*zx32*(zx3-zx2)
              !
              zb=(zx22*(zpp3-zppp)-zx32*(zpp2-zppp))/zd
              za=(zx33*(zpp2-zppp)-zx23*(zpp3-zppp))/zd
              !
           end if
           !
           !  zp = zppp + zx*zx*(za+zx*zb)
           !
           itx=it2
           zth=zth2
           !
           ict=0
40         continue
           !
           zthp=zth
           itx=itx-1
           if(itx.eq.0) itx=inum
           zth=ztha(itx)
           if(zth.lt.zthbmin) zth=zth+twopi
           if(zth.lt.(zthp+0.1_fp*zdthl)) then
              !  zth will jump by ~2pi when zthbmin is crossed
              if(iflag.eq.1) then
                 zx=zth-zthbmin
                 zp = zppp + zx*zx*(za+zx*zb)
                 zpha(itx)=zp
                 ia(itx)=0
              end if
              ict=ict+1
              if(ict.eq.1) itx1=itx
              goto 40  ! haven't reached zthbmin yet
           end if
           !
           if((iflag.eq.0).and.(iflag2.eq.1).and.(ict.gt.1)) then
              zpha(itx1)=0.9_fp*zpha(itx1)+0.1_fp*zpha(it2)
           end if
           !
        end if  ! end of large side region was found
        !
     end if  ! 1st clause in main loop
  end do  ! main loop
  !
  return
end subroutine bphip_minp

! 11sep2004 fgtok
! 11sep2004 fgtok
