! 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
!----------------------------------------------------------------------
!  insert intercept point in xjasubs.blk common
!
!  perform extrapolation if necessary
!   this is done by extrapolating the f fcn (cf subroutine xjacon,
!   xjasub) assuming e independent of position, b varying as 1/r,
!   and flux varying as (r-r0)/a
!
subroutine xjastl(ij,iloci,zxi,zvpll,ze,zb,zr,zvrot)
  !
  use nbi_com
  !
  use xjasub_mod
  !
  !--------------------------------
  !
  !============
  implicit none
  integer :: iloci,ij,isav,iter,imax
  !============
  real(fp) :: zxi,zvpll,ze,zb,zr,zvrot,zxp,zfp,zvpllp,zrp,zbp,zfacx
  real(fp) :: zx,zf,xjasbx,zrk,zpsik,zvpllk,zbk,ztmp,zxl,zfaci
  !============
  iloc(ij)=iloci
  !
  zei(ij)=ze
  !
  if(iloci.gt.0) then
     !  mark the interior point
     zxidpi(ij)=zxi
     zvplli(ij)=zvpll
     zbi(ij)=zb
     zri(ij)=zr
     zvroti(ij)=zvrot
  else
     !
     if(ij.eq.1) then
        isav=i2nzp1
     else if(ij.eq.4) then
        isav=1
     else
        write(nonlin,*) &
             ' ?nubeam xjastl: unexpected call, interior intercept.'
        return
     end if
     !
     !  extrapolation
     iter=0
     !
     zxp=zxi
     zfp=zfa(isav)
     zvpllp=zvpll
     zrp=zr
     zbp=zb
     !
     zfacx=1.2_fp
     !
     imax=10
     !
10   continue
     iter=iter+1
     if(iter.gt.imax) then
        write(nonlin,*) ' ?nubeam xjastl: too many iterations.'
        return
     end if
     !
     zx=zxp*zfacx
     !
     zf=xjasbx(zx,isav,ze,zb,zrk,zpsik,zvpllk,zbk)
     !
     ! careful of overflows here (dmc 4/91):
     ztmp=1.e-10_fp*zfp
     if(((ztmp*zf).le.0.0_fp).or.(iter.ge.imax)) then
        if((ztmp*zf).gt.0.0_fp) then
           zxl=zx
        else
           zxl=max(zxp,min(zx,(zf*zxp-zfp*zx)/(zf-zfp)))
        end if
        zxidpi(ij)=zxl
        zfaci=(zxl-zxp)/(zx-zxp)
        zvplli(ij)=zvpllp+zfaci*(zvpllk-zvpllp)
        zbi(ij)=zbp+zfaci*(zbk-zbp)
        zri(ij)=zrp+zfaci*(zrk-zrp)
        zvroti(ij)=zvrot
     else
        zfp=zf
        zxp=zx
        zvpllp=zvpllk
        zbp=zbk
        zrp=zrk
        goto 10
     end if
     !
  end if
  !
  return
end subroutine xjastl
