!----------------- nbeloc.fnb --------- group nubeam -------------------
!
!  nbeloc - lookup interaction energy index in nb atomic physics table
!   formerly part of nutrav
!   broken out to facilitate generalization of code for rotating plasma
!

subroutine nbeloc(zeovab,irng,ind,zintri,indp1,zintrp)
  !
  use nbatom_mod
  !
  use nbi_com
  !
  !  called from: dep.fnb, nbain2.fnb, nbains.fnb, nutrav.fnb
  !
  !  arguments
  !   input:  zeovab - relative energy *per nucleon* of impinging neutral
  !           irng   - 1st guess of range # in energy grid.
  !
  !   output:
  !     ind     index to g.l.b. table entry rel. to zeovab
  !     zintri  interpolation factor
  !     indp1   index to l.u.b. table entry rel. to zeovab
  !     zintrp  interpolation factor = 1 - zintri
  !     irng    actual range # in energy grid.
  !
  !   use this output to form...
  !   answer = table(ind)*zintri + table(indp1)*zintrp
  !            unless ind=0, in which case
  !   answer = table(indp1)*zintrp
  !------------------------------------------------------------------
  !
  !  changes:
  !
  !  tbt nov 20, 1991 added lsbeam to call list to neatom.
  !  tbt oct 16, 1991 nbatom upgrade. added call to neatom. new grid.
  !
  !------------------------------------------------------------------
  !
  !============
  implicit none
  integer :: irng,ind,indp1,neatom,istat
  !============
  real(fp) :: zintri,zintrp,zeovab
  !============
  indp1 = neatom( zeovab, lsbeam, irng, istat)
  ind   = indp1-1
  !
  zintrp=0.0_fp
  !
  if(istat .eq. 0.or.istat .eq. 3) then
     if(istat .eq. 3) then
        !     .energy > maximum - error.

        write(nonlin,9901) bnsves(nbnsve(lsbeam),lsbeam), zeovab
9901    format(/'  ** error **'/ &
             ' nbeloc cross-section table max energy exceeded ** '/ &
             '   max energy in table (ev):  ',1pe11.4/ &
             '   fast neutral energy (ev):  ',1pe11.4/ &
             '   extrapolation will be taken'/)
        !
        ! mg 12/04/12 make linear extrapolation
        indp1 = lbnsve(nbnsver(lsbeam),lsbeam)
        ind   = indp1-1
     end if
     !     .normal status.
     zintrp=(zeovab              -bnsves(ind,lsbeam))/ &
          (bnsves(indp1,lsbeam)-bnsves(ind,lsbeam))

  else if(istat .eq. 1) then
     ! 	    .energy (zeovab) is between 0. and first energy grid point.
     zintrp = zeovab / bnsves(1,lsbeam)
     if (zintrp .gt. 1._fp.or. zintrp .lt. 0.0_fp)  then
        !     .error - zintrp out of range.

        write(nonlin,9903) zeovab, bnsves(1,lsbeam)
9903    format(/'  ** error **'/ &
             '     nbeloc:  zeovab > bnsves(1,lsbeam) - error ** '/ &
             '     values:  ',1pe11.4, 1pe11.4/)
     end if                 ! zintrp

  else if(istat .eq. 2) then
     !     .energy < 0.0 - error.

     write(nonlin,9902) zeovab
9902 format(/'  ** error **'/ &
          '   nbeloc input energy (zeovab) < 0.0 - error ** '/ &
          '   fast neutral energy (ev):  ',1pe11.4/)
     !
  else

     write(nonlin,9906) istat
9906 format(/'  ** error **'/ &
          '     nbeloc:  istat out of bounds - error ** '/ &
          '     value =:  ', i8 /)
  end if
  !
  if(istat .ne. 3) zintrp=max(0.0_fp,min(1.0_fp,zintrp))
  zintri=1.0_fp-zintrp
  !
  return
end subroutine nbeloc
