!-----------------------------------------------------------------------
!  neatom -- given input e/a find energy index in atomic physics table
!
!   .. dmc 14 oct 1991
!
!   the table energy grid is set up in subroutine ebatom (called in
!   nbatom)  -- this is for
!  the beam-target and beam-beam atomic physics tables set up in
!  nbatom and used in monte carlo deposition and recapture of fast
!  neutrals (cf dep and nutrav).
!
integer function neatom(zeova,isbeam,irange,istat)
  !
  use iso_c_binding, only: fp => c_double
  use nbatom_mod
  !
  !
  !  input -- zeova -- ev/nucleon -- energy per nucleon of the neutral
  !    particle.
  !           isbeam - beam species energy grid table to use.
  !  input/output -- irange -- initial guess as to the "energy range"
  !    on output the actual range into which the neutral falls
  !  output -- istat -- status/completion code:
  !    =0  normal
  !    =1  zeova is btw 0 and the 1st energy point in the table
  !    =2  **error** zeova is negative
  !    =3  **error** zeova exceeds the table max energy
  !
  !  output function value neatom --
  !    index pointing to the smallest table energy value that is
  !    *greater than* the energy zeova
  !    i.e. if neatom returns n the energy zeova lies btw table
  !    energy bnsves(n-1,isbeam) and bnsves(n,isbeam) -- but caution
  !    re istat=1:  n=1 is possible and must be handled specially
  !    by caller.
  !
  !  the species of the neutral is given by index "isbeam"
  !
  !	called from nbeloc and nbelc2.
  !

  implicit none
  real(fp) :: zeova
  integer :: isbeam,irange,istat

  real(fp) :: zeminr,zemaxr
  integer :: ir,ians,inc,inum,i0
  !-----------------------------------------------------------------------
  !
  !  start looking within the specified range -- use min/max in case the
  !  caller does not pass a valid range index
  !
  ir=max(1,min(nbnsver(isbeam),irange))
  istat=0
  !
  !  range search loop
  !
10 continue
  if(ir.eq.1) then
     zeminr=0.0_fp
  else
     zeminr=bnsves(lbnsve(ir-1,isbeam), isbeam)
  end if
  !
  zemaxr=bnsves(lbnsve(ir,isbeam), isbeam)
  !
  if(zeova.lt.zeminr) then
     !
     !  particle energy is in a lower range
     ir=ir-1
     if(ir.eq.0) then
        ians=0
        istat=2
        neatom=ians
        return
     end if
     goto 10                       ! look in lower range
     !
  else if(zeova.gt.zemaxr) then
     !
     !  particle energy is in a higher range
     ir=ir+1
     if(ir.gt.nbnsver(isbeam)) then
        ians=0
        istat=3
        neatom=ians
        return
     end if
     goto 10                       ! look in higher range
     !
  else
     !
     !  correct energy range has been found
     !
     irange=ir
     if(ir.eq.1) then
        i0=0
     else
        i0=lbnsve(ir-1,isbeam)
     end if
     !
     inum=lbnsve(ir,isbeam)-i0
     !
     !  find energy index -- energies evenly spaced within this range
     !
     inc=min((inum-1),int(inum*(zeova-zeminr)/(zemaxr-zeminr)))
     ians=i0+1+inc
     !
     if(ians.eq.1) then
        !  zeova is btw 0.0 and the first energy
        istat=1
     else
        istat=0
     end if
     !
  end if
  !
  !  exit
  !
  neatom=ians
  !
  return
end function neatom
