subroutine r8_nbi_datchkf(nonlin,nlfsp,nlfsp14,nlfst,nlfhe3,nlfhe4, &
       mibs,nsbeam,mimxbe,nlfprod, &
       ebdmaxa,nznbmea,einjs,abeams,xzbeams, &
       nfprot,nfprot14,nftrit,nfhel3,nfhel4,ierct)
 
  use iso_c_binding, only: fp => c_double
  implicit none
 
  !  using no use or include statements, get fusion product species indices
  !  from fusion product request switches
  !  the names are from transp common, originally...
 
  !-------------------------------------------
  !  caution -- dmc jan 2007 -- order of fast 
  !  species is important -- p, t, he3, he4 -- must match
  !  order setting in source/namelist_adapt/namelist_adapt.f90 ****
  !-------------------------------------------

  integer,intent(in) :: nonlin  ! lun for messages
 
  ! switches for fusion product protons, tritons, he3, and he4 respectively:
  logical,intent(in) :: nlfsp,nlfsp14,nlfst,nlfhe3,nlfhe4
 
  integer, intent(in) :: mibs      ! max no. of distinct species
  integer, intent(inout) :: nsbeam   ! no. of distinct species found
  !  ... on input: no. of beam species;
  !  ... on output: no. of beam + fusion product species
 
  integer, intent(in) :: mimxbe    ! no. of energy zones for fusion products
 
  ! the following are "inout" because elements (1:nsbeam) already set
 
  logical, intent(inout) :: nlfprod(mibs)  ! .true. for fusion products
  real(fp), intent(inout) :: ebdmaxa(mibs) ! max energy of distribution fcn
  integer, intent(inout) :: nznbmea(mibs) ! no. of energy zones (copied)
  real(fp), intent(inout) :: einjs(mibs)   ! "injection energy"
  real(fp), intent(inout) :: abeams(mibs),xzbeams(mibs) ! a & z of species
 
  integer, intent(out) :: nfprot   ! index to h fusion product, or zero
  integer, intent(out) :: nfprot14   ! index to h fusion product, or zero
  integer, intent(out) :: nftrit   ! index to t fusion product, or zero
  integer, intent(out) :: nfhel3   ! index to he3 fusion product, or zero
  integer, intent(out) :: nfhel4   ! index to he4 fusion product, or zero
 
  integer,intent(inout) :: ierct  ! incremented if specie classification fails
                                  ! of if max no. species exceeded
  real(fp), parameter ::  elex_ah = 1.0079_fp
 
  !---------------------------------------
  ! local...
 
  integer ifct
  !---------------------------------------
 
  nfprot=0
  nfprot14=0
  nftrit=0
  nfhel3=0
  nfhel4=0
 
  ifct=0
  if(nlfsp) ifct=ifct+1
  if(nlfsp14) ifct=ifct+1
  if(nlfst) ifct=ifct+1
  if(nlfhe3) ifct=ifct+1
  if(nlfhe4) ifct=ifct+1
 
  if(nsbeam+ifct.gt.mibs) then
     write(nonlin,*) ' ?nbi_datchkf:  too many distinct fast ion species.'
     write(nonlin,*) '  max no. of distinct species: ',mibs
     write(nonlin,*) '  no. of beam species requested: ',nsbeam
     write(nonlin,*) '  no. of fusion products requested: ',ifct
     ierct=ierct+1
     return
  endif
 
  if(nlfsp) then
 
     ! dd -> p,t reaction (p product)
 
     nsbeam=nsbeam+1
     nfprot=nsbeam
 
     nlfprod(nsbeam)=.true.
     call r8_nbfusn_emax(1,1,ebdmaxa(nsbeam))
     nznbmea(nsbeam)=mimxbe
     einjs(nsbeam)=3.02e+6_fp
     abeams(nsbeam)=1.0_fp
     xzbeams(nsbeam)=1.0_fp
 
  endif

  if(nlfsp14) then
 
     ! dhe3 -> p,he4 reaction (p product)
 
     nsbeam=nsbeam+1
     nfprot14=nsbeam
 
     nlfprod(nsbeam)=.true.
     nznbmea(nsbeam)=mimxbe
     einjs(nsbeam)=14.7e+6_fp
     ebdmaxa(nsbeam)=14.7e+6_fp*1.25
     abeams(nsbeam)=1.0_fp
     xzbeams(nsbeam)=1.0_fp
 
  endif
 

  if(nlfst) then
 
     ! dd -> p,t reaction (t product)
 
     nsbeam=nsbeam+1
     nftrit=nsbeam
 
     nlfprod(nsbeam)=.true.
     call r8_nbfusn_emax(1,3,ebdmaxa(nsbeam))
     nznbmea(nsbeam)=mimxbe
     einjs(nsbeam)=1.01e+6_fp
     abeams(nsbeam)=3.016049_fp/elex_ah !convert from chemical definition to proton mass
     xzbeams(nsbeam)=1.0_fp
 
  endif
 
  if(nlfhe3) then
 
     ! dd -> n,he3 reaction (he3 product)
 
     nsbeam=nsbeam+1
     nfhel3=nsbeam
 
     nlfprod(nsbeam)=.true.
     call r8_nbfusn_emax(2,3,ebdmaxa(nsbeam))
     nznbmea(nsbeam)=mimxbe
     einjs(nsbeam)=0.82e+6_fp
     abeams(nsbeam)=3.016029_fp/elex_ah !convert from chemical definition to proton mass
     xzbeams(nsbeam)=2.0_fp
 
  endif
 
  if(nlfhe4) then
 
     ! dt -> n,he4 reaction (alpha product)
     ! tt -> 2n,he4 also.
 
     nsbeam=nsbeam+1
     nfhel4=nsbeam
 
     nlfprod(nsbeam)=.true.
     call r8_nbfusn_emax(2,4,ebdmaxa(nsbeam))
     nznbmea(nsbeam)=mimxbe
     einjs(nsbeam)=3.50e+6_fp
     abeams(nsbeam)=4.0026_fp/elex_ah !convert from chemical definition to proton mass
     xzbeams(nsbeam)=2.0_fp
 
  endif
 
end subroutine r8_nbi_datchkf

subroutine nbfusn_emax(iz,ia,emax)

  use iso_c_binding, only: fp => c_double
  implicit none

  integer, intent(in) :: iz  ! z of fusion product
  integer, intent(in) :: ia  ! a of fusion product

  real, intent(out) :: emax  ! a number safely above the birth energy

  !-------------------------------------

  emax = 0.0

  if(iz.eq.1) then
     if(ia.eq.1) then
        emax = 5.00e+6_fp
     else if(ia.eq.3) then
        emax = 1.25e+6_fp
     endif
  else if(iz.eq.2) then
     if(ia.eq.3) then
        emax = 1.00e+6_fp
     else if(ia.eq.4) then
        emax = 5.00e+6_fp
     endif
  endif

end subroutine nbfusn_emax

subroutine r8_nbfusn_emax(iz,ia,emax)

  use iso_c_binding, only: fp => c_double
  implicit none

  integer, intent(in) :: iz  ! z of fusion product
  integer, intent(in) :: ia  ! a of fusion product

  real(fp), intent(out) :: emax  ! a number safely above the birth energy

  !-------------------------------------

  emax = 0.0_fp

  if(iz.eq.1) then
     if(ia.eq.1) then
        emax = 5.00e+6_fp
     else if(ia.eq.3) then
        emax = 1.25e+6_fp
     endif
  else if(iz.eq.2) then
     if(ia.eq.3) then
        emax = 1.00e+6_fp
     else if(ia.eq.4) then
        emax = 5.00e+6_fp
     endif
  endif
end subroutine r8_nbfusn_emax
