!----------------- nbelc2.fnb --------- group nubeam -------------------
!
!  nbelc2 - lookup interaction energy index in nb atomic physics table
!   formerly part of nutrav
!   called from nbain2 only for beam-beam cross section data.
!   broken out to facilitate generalization of code for rotating plasma
!

subroutine nbelc2(zeovab,irange,isbeam,ind,zintri,indp1,zintrp)
  use iso_c_binding, only: fp => c_double
  !
  use nbatom_mod
  !
  !
  !  called from: dep.fnb, nbain2.fnb, nbains.fnb, nutrav.fnb
  !
  !  arguments
  !   input:  zeovab - relative energy *per nucleon* of impinging neutral
  !           irange   - 1st guess of range # in energy grid.
  !	    isbeam - # of beam species to use in look up table.
  !
  !   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
  !     irange    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
  !
  implicit none
  real(fp) :: zeovab,zintri,zintrp
  integer :: irange,isbeam,ind,indp1

  integer :: istat

  integer :: neatom

  !------------------------------------------------------------------
  !
  !  changes:
  !
  !  tbt nov 18, 1991 created from nbeloc. added argument isbeam.
  !                   added isbeam to call of neatom.
  !  tbt oct 16, 1991 nbatom upgrade. added call to neatom. new grid.
  !
  !------------------------------------------------------------------

  indp1 = neatom( zeovab, isbeam, irange, 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(lun_nbatom_mod,9901) bnsves(nbnsve(isbeam),isbeam), &
             zeovab
9901    format(/'  ** error **'/ &
             ' nbelc2 cross-section table max energy exceeded ** '/ &
             '   max energy in table (ev):  ',1pe11.4/ &
             '   fast neutral energy (ev):  ',1pe11.4/ &
             '   extrapolation will be taken'/)
        indp1 = lbnsve(nbnsver(isbeam),isbeam)
        ind   = indp1-1
     end if
     !
     !     .normal status.
     zintrp=(zeovab-bnsves(ind,isbeam))/ &
          (bnsves(indp1,isbeam)-bnsves(ind,isbeam))

  else if(istat .eq. 1) then
     ! 	    .energy (zeovab) is between 0. and first energy grid point.
     zintrp = zeovab / bnsves(1,isbeam)
     if (zintrp .gt. 1._fp.or. zintrp .lt. 0.0_fp)  then
        !	      .error - zintrp out of range.
        write(lun_nbatom_mod,9903) zeovab, bnsves(1,isbeam)
9903    format(/'  ** error **'/ &
             '     nbelc2:  zeovab > bnsves(1,isbeam) - error ** '/ &
             '     values:  ',1pe11.4, 1pe11.4/)
        zintrp=max(0.0_fp,min(1.0_fp,zintrp))
     end if                         ! zintrp

  else if(istat .eq. 2) then
     !	    .energy < 0.0 - error.
     write(lun_nbatom_mod,9902) zeovab
9902 format(/'  ** error **'/ &
          '   nbelc2 input energy (zeovab) < 0.0 - error ** '/ &
          '   fast neutral energy (ev):  ',1pe11.4/)
     !
  else
     write(lun_nbatom_mod,9906) istat
9906 format(/'  ** error **'/ &
          '     nbelc2:  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 nbelc2
