!-----------------------------------------------------------------------
!  ndxfbe -- find index to fast ion dist. fcn energy zone
!
integer function ndxfbe_r8(zenergy,isbeam,istart)
  !
  !  input
  !    zenergy = energy, ev
  !    isbeam  = nubeam fast ion species index
  !    istart  = where to start search
  !
  !  output
  !    ndxfbe_r8  = index to fast ion dist fcn zone containing the energy
  !              zenergy, or 0 if zenergy is negative or .gt. maximum
  !
  use nbi_com
  use logmod
  !
  implicit none
  real(fp), intent(in) :: zenergy     ! energy, ev
  integer, intent(in) :: isbeam     ! nubeam specie index
  integer, intent(in) :: istart     ! search index (initial guess)
  !
  !-----------------------------------------------------------------------
  integer :: is,in,inp1,idir,ism1,isp1,isp,ideln
  integer :: ipmin, ipmax
  real(fp) :: ze1,ze2
  integer, save :: iwarn = 50
  integer, save :: jwarn = 99
  !-----------------------------------------------------------------------
  !
  is=istart
  in=nznbmea(isbeam)
  inp1=in+1
  !
  if((is.gt.0).and.(is.le.in)) then
     !  see if still in same zone as before
     idir=0
     if(efbmb(is,isbeam).gt.zenergy) then
        idir=-1    ! energy is below range of starting zone
     else
        isp1=is+1
        if(efbmb(isp1,isbeam).lt.zenergy) then
           idir=1   ! energy is above range of starting zone
           is=isp1
        end if
     end if
     if(idir.eq.0) then
        !  in same zone already
        ndxfbe_r8=is
        return
     end if
  else
     !  start index out of range; will start from bottom up
     is=1
     idir=1
  end if
  !
  !  not in same zone, must search...
  !
  !  first check neighbouring zones
  !
  !  down 1
  if(idir.eq.-1) then
     if(is.gt.1) then
        ism1=is-1
        if((efbmb(ism1,isbeam).le.zenergy).and. &
             (zenergy.le.efbmb(is,isbeam))) then
           ndxfbe_r8=ism1
           return
        end if
     end if
  end if
  !
  !  up 1
  if(idir.eq.1) then
     if(is.lt.inp1) then
        isp1=is+1
        if((efbmb(is,isbeam).le.zenergy).and. &
             (zenergy.le.efbmb(isp1,isbeam))) then
           ndxfbe_r8=is
           return
        end if
     end if
  end if
  !
  !  check if energy is in range at all...
  !
  if((zenergy.lt.efbmb(1,isbeam)).or. &
       (zenergy.gt.efbmb(inp1,isbeam))) then
     !  energy completely out of range
     if(zenergy.lt.efbmb(1,isbeam)) ndxfbe_r8=0
     if(zenergy.gt.efbmb(inp1,isbeam)) ndxfbe_r8=inp1
     if(iwarn.gt.0) then
        jwarn=jwarn+1
        if(jwarn.eq.100) then
           jwarn=0
           iwarn=iwarn-1
           !              write(nonlin,*) ' !ndxfbe warning: energy out of range.'
           !              write(nonlin,*) '   species, e, emin, emax: ',isbeam,
           !     >             zenergy, efbmb(1,isbeam), efbmb(inp1,isbeam)
           call errorlog( ' !ndxfbe warning: energy out of range.')
           call errorlog('   species: ',isbeam)
           call errorlog('   e: ',zenergy)
           call errorlog('   emin: ',efbmb(1,isbeam))
           call errorlog('   emax: ',efbmb(inp1,isbeam))
        end if
     end if
     return
  end if
  !
  !  ok, use generalized interpolation search
  !
  ! rga, decreasing bracket   ipmin<=is<=ipmax  maintained
  !      use the above condition that the nearby intervals do not include zenergy
  !
  if (idir==1) then
     ipmin=min(in,is+1)
     ipmax=in
  else
     ipmin=1
     ipmax=max(1,is-2)
  end if
  !
10 continue
  if (ipmin>=ipmax) then
     if (zenergy>=efbmb(ipmin,isbeam) &
          .and. zenergy<=efbmb(ipmin+1,isbeam)) then
        ndxfbe_r8=ipmin
        return
     else
        call errorlog( ' ?ndxfbe: logic error by rga')
        !            write(nonlin, *) '?ndxfbe: logic error by rga'
        call bad_exit
     end if
  end if

  ze1 = efbmb(ipmin,isbeam)
  ze2 = efbmb(ipmax+1,isbeam)
  ideln = ipmax+1-ipmin
  is  = ipmin +(zenergy-ze1)*ideln/(ze2-ze1)
  is  = max(ipmin,min(is,ipmax))

  if (zenergy<efbmb(is,isbeam)) then
     ipmax=is-1
  else if (zenergy>efbmb(is+1,isbeam)) then
     ipmin=is+1
  else
     ndxfbe_r8=is
     return
  end if
  !
  goto 10
  !
end function ndxfbe_r8
!--------------------------------------
