! dmc apr 2004 -- all xja routines use bmnmx instead of bmidp -- code
!  now in terms of min and max b on flux surfaces, not midplane values
!----------------------------------------------------------------------
!  xjasng  look for "singular" roots of orbit intercept function
!
!    called from xjacon; this routine deals with the situation where
!  f on the flux bdy grid approaches the axis but then moves away again
!  without crossing; there may be two roots btw. grid pts, rather than
!  none...
!
subroutine xjasng(ick,ifound,zxi,zvpll,ze,zb,zr,zvrot)
  !
  !  f(ick-1) thru f(ick+1) are all of same sign; axis crossing suspected
  !  close to f(ick)
  !
  use nbi_com
  !
  use xjasub_mod
  !
  !============
  implicit none
  integer :: isign,ifound,ick,ickp,ickm,idir,isid,i1,i2
  !============
  real(fp) :: zftmp,zsfmax,zfmm,zfpp,zfppa,zfmma,zfsma,zfspa,zs1
  real(fp) :: ztes1,zs4,ztes4,zxf2,zxf1,zwtol,zpsi,zbtb
  !============
  real(fp) :: zxi(2),zvpll(2),ze(2),zb(2),zr(2),zvrot(2)
  !
  real(fp) :: zfwrk(4)
  !
  !------------------
  !
  ifound=0
  !
  isign=sign(1.0_fp,zfa(ick))
  ickp=ick+1
  ickm=ick-1
  ! test code
  zftmp=1.e-10_fp*zfa(ick)
  if((zftmp*zfa(ickp).le.0.0_fp).or. &
       (zftmp*zfa(ickm).le.0.0_fp)) then
     write(nonlin,*) ' ?nubeam xjasng -- call/algorithm error.'
     return
  end if
  ! end test code
  !
  !  test if singularity is local maximum modulus instead of local minimum
  zsfmax=max(isign*zfa(ickp),isign*zfa(ickm))
  if(isign*zfa(ick).gt.zsfmax) return
  !
  !  ok have local minimum.  fetch values 2 pts away in each direction
  !  careful:  can go off end of array...
  !
  if(ickm.gt.1) then
     zfmm=zfa(ickm-1)
  else
     zfmm=zfext1
  end if
  !
  if(ickp.lt.i2nzpx) then
     zfpp=zfa(ickp+1)
  else
     zfpp=zfext2
  end if
  !
  !  look for root towards side with points with strongest tendancy
  !  to zero first...
  !
  zfppa=max(0.1_fp*isign*zfa(ickp),isign*(zfpp-zfa(ickp)))
  zfmma=max(0.1_fp*isign*zfa(ickm),isign*(zfmm-zfa(ickm)))
  !
  zfsma=isign*zfa(ickm)/zfmma
  zfspa=isign*zfa(ickp)/zfppa
  !
  if(zfsma.lt.zfspa) then
     idir=1
  else
     idir=-1
  end if
  !
  do isid=1,2
     idir=-idir
     !
     !  mod dmc oct 1989 -- in this loop, look on both sides of singular
     !  pt for a numeric intercept-- instead of only on one preferred side
     !  as before
     !
     if(idir.eq.-1) then
        zfwrk(4)=zfa(ickp)
        zfwrk(3)=zfa(ick)
        zfwrk(2)=zfa(ickm)
        !
        i1=ickm
        i2=ick
        !
        zfwrk(1)=zfmm
        !
     else
        zfwrk(1)=zfa(ickm)
        zfwrk(2)=zfa(ick)
        zfwrk(3)=zfa(ickp)
        !
        i1=ick
        i2=ickp
        !
        zfwrk(4)=zfpp
        !
     end if
     !
     !  array zfwrk now contains 4 consecutive pts from f
     !   which look like:
     !
     !   *            *
     !       *    *                     i1  i2
     ! -------------------  or  ----------------------
     !       i1   i2                    *   *
     !                             *           *
     !
     !
     !  if this is a valid numerical singularity,
     !  changes just beyond interval must substantially exceed values on
     !  interval
     !
     zs1=isign*(zfwrk(1)-zfwrk(2))
     ztes1=zs1/(isign*zfwrk(2))
     if(ztes1.lt.0.1_fp) cycle
     !
     zs4=isign*(zfwrk(4)-zfwrk(3))
     ztes4=zs4/(isign*zfwrk(3))
     if(ztes4.lt.0.1_fp) cycle
     !
     if(max(ztes1,ztes4).lt.10.0_fp) then
        if(min(ztes1,ztes4).lt.1.0_fp) cycle
     end if
     !
     zxf2=1.0_fp-isign*zfwrk(3)/zs4
     zxf1=isign*zfwrk(2)/zs1
     !
     zwtol=0.2_fp
     if(zxf1.lt.zwtol) then
        zxf2=max(zxf1,zxf2)
     else if(zxf2.gt.(1.0_fp-zwtol)) then
        zxf1=min(zxf1,zxf2)
     else if(zxf2.lt.(zxf1+zwtol)) then
        cycle
     end if
  !  intercepts taken as found
  !
     ifound=2
  !
     zxi(1)=zxia(i1)+zxf1*(zxia(i2)-zxia(i1))
     ze(1)=zea(i1)+zxf1*(zea(i2)-zea(i1))
     zb(1)=zba(i1)+zxf1*(zba(i2)-zba(i1))
     zvrot(1)=zvrota(i1)+zxf1*(zvrota(i2)-zvrota(i1))
     zr(1)=rmnmx(i1)+zxf1*(rmnmx(i2)-rmnmx(i1))
     zpsi=xzbeami*(psmnmx(i1)+zxf1*(psmnmx(i2)-psmnmx(i1)))
     zbtb=btbmnmx(i1)+zxf1*(btbmnmx(i2)-btbmnmx(i1))
     zvpll(1)=1.e7_fp*xe*(zpsi-zpsito)/(zmass*zr(1))/zbtb
     !
     zxi(2)=zxia(i1)+zxf2*(zxia(i2)-zxia(i1))
     ze(2)=zea(i1)+zxf2*(zea(i2)-zea(i1))
     zb(2)=zba(i1)+zxf2*(zba(i2)-zba(i1))
     zvrot(2)=zvrota(i1)+zxf2*(zvrota(i2)-zvrota(i1))
     zr(2)=rmnmx(i1)+zxf2*(rmnmx(i2)-rmnmx(i1))
     zpsi=xzbeami*(psmnmx(i1)+zxf2*(psmnmx(i2)-psmnmx(i1)))
     zbtb=btbmnmx(i1)+zxf2*(btbmnmx(i2)-btbmnmx(i1))
     zvpll(2)=1.e7_fp*xe*(zpsi-zpsito)/(zmass*zr(2))/zbtb
     return  ! intercepts found, leave now
  end do
  !
  return      ! intercepts not found (fell thru loop)
  !
end subroutine xjasng
