subroutine r8_nbi_efbm(ilfprod,ilfi_mcrf,edbfac,inum,zemax,ilin,zelin,efbmb,efbm)
  !
  !  compute the energy grid for a mc fast ion specie dist. fcn.
  !  (code taken from old trcore/datckb.for -- dmc 10/04/01)
 
  use iso_c_binding, only: fp => c_double
  implicit none
 
  logical,intent(in) :: ilfprod    ! .true. if this is a fusion product specie
  logical,intent(in) :: ilfi_mcrf  ! .true. if this is a rf kick off operator product specie
                                   ! or maxwellian test
  integer,intent(in) ::  inum      ! total no. of zones
  real(fp),intent(in) ::  edbfac     ! max energy factor, for rf kick off operator
  real(fp),intent(in) ::  zemax      ! max energy
 
  integer,intent(out) ::  ilin   ! no. of linearly spaced zones; ilin.ge.1
  real(fp),intent(out) ::  zelin     ! max. linearly spaced energy (min is zero)
 
  real(fp),intent(out) ::  efbmb(inum+1)  ! energy grid bdys
  real(fp),intent(out) ::  efbm(inum)     ! energy grid zone ctrs

  !  beam species -- linear grid
  !  fusion products of rf kick -- linear at low energies, otherwise logarithmic
 
  !$r8real_input: zelin,zemax
  !$r8real_output: efbmb,efbm
 
  ! note if ilin.eq.inum, zelin.eq.zemax is expected; no checks here!
 
 
  ! local...
 
  integer il,ilp1
  real(fp) zfac
 
  !---------------------------------
  if(ilfprod) then
     ilin=inum/10
     zelin=zemax/100
  else
     if(ilfi_mcrf) then
        ilin=inum*2/10 
        zelin=zemax/edbfac
     else
        ilin=inum
        zelin=zemax
     endif
  endif
 
  efbmb(1)=0.0_fp
 
!  linear part
  do il=1,ilin
     ilp1=il+1
     efbmb(ilp1)=il*zelin/ilin
  enddo
 
!  log part
  if(inum.gt.ilin) then
     zfac=exp(log(zemax/zelin)/(inum-ilin))
     do il=ilin+1,inum
        ilp1=il+1
        efbmb(ilp1)=zfac*efbmb(il)
     enddo
  endif
 
!  fill in zone ctr values
  do il=1,inum
     efbm(il)=(efbmb(il)+efbmb(il+1))/2
  enddo
 
end subroutine r8_nbi_efbm
