! 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
!---------------------------------------------------------------------
!  xjaset  look for collisionless orbit midplane intercepts
!    this means finding the roots of the function computed in xjasub
!  called via xjacon.  this routine finds the set of all such roots,
!  up to four:
!
!   co-outer --+
!  ctr-outer --+outer banana
!  ctr-inner
!   co-inner
!
subroutine xjaset(jdep,zvion)
  !
  !  input
  !    jdep	deposition zone
  !   zvion       velocity at ionization point
  !
  !  in xjasubs.blk common--
  !    zmass      ion mass
  !    zekin      ion kinetic energy
  !    zepot      ion (electrostatic) potential energy
  !    zmub       ion mu = vperp**2/mod(b)  ! constant of motion
  !    zpsito     ion canonical momentum
  !                      = m*vpll*r + e/c psi  ! constant of motion
  !
  !  output
  !    updated intercept information in xjasubs.blk common
  !
  !----------------------------
  !  dmc 7 nov 1988-- debugging
  !
  !    a numerical singularity can occur with f=0 and df/dx=0 at
  !  the same time, appearing numerically as too many f=0 intercepts.
  !  this is very rare.  fix it by compiling a list of regions of f
  !  not changing sign, and if too many intercepts occur, blow away
  !  the region of least measure as arbitrarily defined below ...
  !    jun 1990 -- patch also modified to deal with intercepts
  !  inferred by interpolation cf xjasng
  !
  !  dmc 14 feb 1989-- debugging
  !
  !   especially when depositing perp ions, all roots of f may fall
  !  between grid points and elude the search algorithm.  the roots of
  !  f must exist, so the fudge of subtracting a small constant from
  !  the numerical representation of f is introduced with a warning
  !  message and the search loop is re-entered.  the maximum that may
  !  be subtracted from f is 0.0002*zvion**2.
  !
  !  dmc 27 july 1992-- debugging -- expanding on 14 feb 1989 adjustment
  !
  !  if the user puts in a noisy rotation profile, with large changes
  !  spanning a single grid zone, the problem of search algorithm failure
  !  can be aggravated.  as we are seeing such data (on a few shots) i
  !  am expanding the correction fudge factor limit to 0.02*zvion**2
  !
  !-----------------------------------------------------------------
  !  note dmc oct 2003 -- this routine looks for midplane intercepts
  !  of orbits-- by seeking the roots of a function constructed on the
  !  basis of conservation of energy and momentum.  as has been known
  !  for some time, however, some orbits do not have midplane intercepts.
  !  this can happen e.g. if there is a magnetic well off the midplane
  !  as can happen in st mhd equilibria.  orbits lacking midplane
  !  intercepts will typically fail to find roots of the midplane
  !  function; the minimum value, instead, will be taken as an
  !  approximation.
  !
  !  warnings about this will only be issued a few times per run.
  !-----------------------------------------------------------------
  !
  use nbi_com
  !
  use xjasub_mod
  !
  !============
  implicit none
  integer :: jdep,ij,inumi,isigi,itmax,iter,istart,ilim,isrch
  integer :: ifound,iflagp,izap,ii,ini,inizap,izon00,iprev,il
  integer :: isign,isignl,iz001,iz002
  !============
  real(fp) :: zvion,zfmin,zzfa,xjasub,xjasbx,zxtarg,zrk,zpsik,zvk
  real(fp) :: zbk,zxif,zvpllf,zef,zbf,zrf,zvrotf,zfcorr,zxip,zvpllp
  real(fp) :: zep,zbp,zrp,zvrotp,zzap,zf1,zf2,zden,zfac,zmmx
  !============
  real(fp) :: zsumi(12)      ! measure of regions of same sign zfa
  integer :: indi(2,12)  ! indices to regions of same sign zfa
  !
  integer, save :: iwarn = 0
  integer, parameter :: iwarn_max = 5
  !
  !-----------------------------------------------------------------
  !  executable code --
  !
  !  start by looking for outer co intercept
  !
  !  look for outer co intercept first **
  !
  !  at outer bdy...
  !
  do ij=1,4
     iloc(ij)=-1
     zxidpi(ij)=0.0_fp
     zvplli(ij)=0.0_fp
     zei(ij)=0.0_fp
     zbi(ij)=0.0_fp
     zvroti(ij)=0.0_fp

  end do
  !
  !  evaluate intercept function
  !
  inumi=0
  isign=0
  !
  zfmin=epsinv
  !
  do ij=1,i2nzpx
     if(ij.le.i2nzp1) then
        zzfa=xjasub(ij,zxia(ij), &
             zvplla(ij),zea(ij),zba(ij),zvrota(ij))
     else
        !  beyond the plasma, large major radius side ...
        zvrota(ij)=zvrota(i2nzp1)
        zea(ij)=zea(i2nzp1)
        zxia(ij)=1.0_fp+real(ij-i2nzp1,fp)/nzones
        zzfa=xjasbx(zxia(ij),i2nzp1,zea(i2nzp1),zba(i2nzp1), &
             rmnmx(ij),psmnmx(ij),zvplla(ij),zba(ij))
        btbmnmx(ij)=btbmnmx(i2nzp1)
     end if
     !  ! convenient to bind away from zero (1 cm/sec is close enuf!)
     if(abs(zzfa).lt.1.0_fp) zzfa=1.0_fp
     zfa(ij)=zzfa
     zfmin=min(zzfa,zfmin)
     !  sign of region ...
     if(zzfa.gt.0.0_fp) then
        isigi=1
     else
        isigi=-1
     end if
     if(isigi.ne.isign) then
        !  sign change
        isign=isigi
        inumi=inumi+1
        zsumi(inumi)=abs(zzfa)
        indi(1,inumi)=ij
        !  close out previous sector?
        if(inumi.gt.1) indi(2,inumi-1)=ij-1
     else
        !  sum measure of this region of f of the same sign...
        zsumi(inumi)=zsumi(inumi)+abs(zzfa)
     end if

  end do
  !  close out last sector
  indi(2,inumi)=i2nzpx
  !
  !  f extrapolated one zone width at inner and outer bdys:
  !  inside
  zxtarg=1.0_fp+1.0_fp/nzones
  zfext1=xjasbx(zxtarg,1,zea(1),zba(1),zrk,zpsik,zvk,zbk)
  !  outside
  zxtarg=1.0_fp+real(i2nzpx-i2nzp1+1,fp)/nzones
  zfext2=xjasbx(zxtarg,i2nzp1,zea(i2nzp1),zba(i2nzp1), &
       zrk,zpsik,zvk,zbk)
  !bugfix-comment-this-line-dmc-feb-1997	zfmin=min(zfmin,zfext1,zfext2)
  !
  !  re-entry point:
  itmax=5
  iter=0
  !
50 continue
  !
  istart=0
  !
  ilim=4
  if(zfa(i2nzpx).lt.0.0_fp) then
     !  its outside the plasma/ mapped region
     call xjastl(1,0, &
          zxia(i2nzp1),zvplla(i2nzp1),zea(i2nzp1),zba(i2nzp1), &
          rmnmx(i2nzp1),zvrota(i2nzp1))
  else
     isrch=1
     call xjafnd(istart,ifound,1,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     if(ifound.eq.0) then
        !  dmc feb 1989 -- permit retry with fudge on root not found
        iter=iter+1
        if(iter.le.itmax) then
           zfcorr=min(zfmin, (2.0_fp*zvion**2/10.0_fp**(itmax-iter)) )
           zfext1=zfext1-zfcorr
           if(abs(zfext1).lt.1.0_fp) zfext1=1.0_fp
           zfext2=zfext2-zfcorr
           if(abs(zfext2).lt.1.0_fp) zfext2=1.0_fp
           if(itmax.ge.4) then
              if(iwarn.lt.iwarn_max) then
                 iwarn=iwarn+1
                 write(nonlin,*) ' **xjaset: .ge.3 iterations...'
                 write(nonlin,*) ' zfmin,zfcorr=',zfmin,zfcorr
                 write(nonlin,*) ' zfext1,zfext2=',zfext1,zfext2
              end if
           end if
           do ij=1,i2nzpx
              zfa(ij)=zfa(ij)-zfcorr
              if(abs(zfa(ij)).lt.1.0_fp) zfa(ij)=1.0_fp
              if(ij.eq.1) then
                 zfmin=zfa(ij)
              else
                 zfmin=min(zfmin,zfa(ij))
              end if

           end do
           zfmin=max(100.0_fp,zfmin)

           if(iwarn.lt.iwarn_max) write(nonlin,9901) iter
9901       format( &
                ' %xjaset -- root not found, trying standard fixup iteration=',i2)
           if(iter.eq.itmax) then

              if(iwarn.lt.iwarn_max) write(nonlin,9911)
9911          format(' xjaset ** max iteration reached **')
           end if
           goto 50
        else
           if(iwarn.lt.iwarn_max) write(nonlin,9902)
9902       format(' ?nubeam xjaset-- no orbit intercepts found')
        end if
     else
        call xjastl(1,istart,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     end if
  end if
  !
  !---------------------------
  !  outer co intercept located; find next intercept in; this is either
  !  the outer ctr intercept or the inner co intercept
  !
  isrch=2
  call xjafnd(istart,ifound,-1,zxif,zvpllf,zef,zbf,zrf,zvrotf)
  if(ifound.eq.0) then
     call xjastl(4,0,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     return
  else
     !  nature of intercept not yet known
     !
     zxip=zxif
     zvpllp=zvpllf
     zep=zef
     zbp=zbf
     zrp=zrf
     zvrotp=zvrotf
     iflagp=istart
     !
  end if
  !
300 continue
  call xjafnd(istart,ifound,1,zxif,zvpllf,zef,zbf,zrf,zvrotf)
  if(ifound.gt.0) then
     if(iflagp.gt.0) then
        ! previously found intercept was outer ctr
        !
        call xjastl(2,iflagp,zxip,zvpllp,zep,zbp,zrp,zvrotp)
        !
     end if
     !
     call xjastl(3,istart,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     !
  else
     !  no more intercepts
     !
     if(iflagp.eq.0) then
        !  none found at all...
        write(nonlin,*) &
             ' ?nubeam xjaset -- no orbit midplane intercepts.'
     else
        !
        !  assume trapped if pll velocity .lt. rotation velocity, passing
        !  w/o trap intercepts otherwise
        !
        if(zvpllp.lt.zvrotp) then
           call xjastl(2,iflagp,zxip,zvpllp,zep,zbp,zrp,zvrotp)
        else
           call xjastl(4,iflagp,zxip,zvpllp,zep,zbp,zrp,zvrotp)
        end if
        return
     end if
     !
  end if
  !
  !  look for 4th intercept...
  call xjafnd(istart,ifound,-1,zxif,zvpllf,zef,zbf,zrf,zvrotf)
  if(ifound.gt.0) then
     call xjastl(4,istart,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     !  check for error-- too many intercepts
     call xjafnd(istart,ifound,1,zxif,zvpllf,zef,zbf,zrf,zvrotf)
     if(ifound.gt.0) then
        !===================================================
        !  f=0 and df/dx=0 at same time singularity?  try blowing away
        !  a small glitch in f...
        izap=0
        zzap=epsinv
        do ii=1,inumi
           !  to qualify, interval must be of small measure
           if(zsumi(ii).lt.0.003_fp*zekin) then
              !  and must be of small extent
              ini=indi(2,ii)-indi(1,ii)+1
              if(ini.le.2) then
                 if(zsumi(ii).lt.zzap) then
                    zzap=zsumi(ii)  ! ok, mark this zone
                    izap=ii
                    inizap=ini
                 end if
              end if
           end if

        end do
        if(izap.gt.0) then
           write(nout,9904)
9904       format( &
                ' %xjaset-- too many orbit intercepts, recovery attempt')
           !  blow away the marked interval-- find surrounding f values
           if(izap.eq.1) then
              zf1=zfa(indi(2,izap)+1)
              zfext1=zf1
           else
              zf1=zfa(indi(1,izap)-1)
              zsumi(izap-1)=epsinv
           end if
           if(izap.eq.inumi) then
              zf2=zfa(indi(1,izap)-1)
              zfext2=zf2
           else
              zf2=zfa(indi(2,izap)+1)
              zsumi(izap+1)=epsinv
           end if
           zsumi(izap)=epsinv
           !
           zden=1.0_fp/real(inizap+1,fp)
           do ij=indi(1,izap),indi(2,izap)
              zfac=zden*real(ij+1-indi(1,izap),fp)
              zfa(ij)=zf1*(1.0_fp-zfac)+zf2*zfac

           end do
           !  try again:
           goto 50
        else
           !  dmc june 1990 -- patching
           izon00=0
           iprev=ifound
           do il=4,1,-1
              if(iloc(il).eq.iprev) izon00=iloc(il)
              iprev=iloc(il)
           end do
           if(izon00.gt.0) then
              !  pair of intercepts found by xjasng -- remove by
              !  patching local zfa values
              write(nout,9904)
              if(zfa(izon00).gt.0.0_fp) then
                 isignl=1
              else
                 isignl=-1
              end if
              zmmx=0.0_fp
              iz001=max(1,izon00-1)
              iz002=izon00+2
              do il=iz001,iz002
                 zmmx=max(zmmx,isignl*zfa(il))
              end do
              zmmx=isignl*zmmx
              if(zfa(izon00)*zfa(max(1,izon00-1)).ge.0.0_fp) &
                   zfa(izon00)=zmmx
              if(zfa(izon00+1)*zfa(izon00+2).ge.0.0_fp) &
                   zfa(izon00+1)=zmmx
              !     try again:
              goto 50
           else
              write(nonlin,9903) jdep
9903          format( &
                   ' %xjaset-- too many orbit intercepts found; jdep=',i3/ &
                   ' ...intercept(s) at inside major radius ignored!')
           end if
        end if
     end if
     !===================================================
     !  4th intercept is beyond plasma on small major radius side...
     !
  else
     call xjastl(4,0,zxif,zvpllf,zef,zbf,zrf,zvrotf)
  end if
  return
end subroutine xjaset
!======================================================================
! 11sep2004 fgtok
