subroutine r8_nbi_datchkb(nonlin,nbeam,abeama,xzbeama, &
       mibs,abeams,xzbeams,nsbeam, &
       nbhydr,nbdeut,nbtrit,nbhel3,nbhel4,nlheavy,ierct)
 
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zel, vc, ps_mp=>mp_kg, zmtocm
  use periodic_table_mod
  implicit none
 
  !  using no use or include statements, get beam species indices
  !  from neutral beam species a & z lists
  !  the names are from transp common, originally...
 
integer,intent(in) :: nonlin  ! lun for messages
 
  integer,intent(in) :: nbeam      ! no. of beams
  real(fp),intent(inout) :: abeama(nbeam)
  real(fp),intent(in) ::  xzbeama(nbeam)  ! a & z of beams
 
  integer, intent(in) :: mibs      ! max no. of distinct species
 
  real(fp), intent(out) :: abeams(mibs),xzbeams(mibs) ! list of distinct species
  integer, intent(out) :: nsbeam   ! no. of distinct species found
 
  integer, intent(out) :: nbhydr   ! index to h beam specie, or zero
  integer, intent(out) :: nbdeut   ! index to d beam specie, or zero
  integer, intent(out) :: nbtrit   ! index to t beam specie, or zero
  integer, intent(out) :: nbhel3   ! index to he3 beam specie, or zero
  integer, intent(out) :: nbhel4   ! index to he4 beam specie, or zero
 
  logical, intent(out) :: nlheavy(mibs)  ! flag species other than h/he
 
  integer,intent(inout) :: ierct  ! incremented if specie classification fails
                                  ! of if max no. species exceeded
 
  !$r8real_input: abeama,xzbeama
  !$r8real_output: abeams,xzbeams
 
  !---------------------------------------
  ! local...
 
  integer ib,isb,imatch
  integer ier
  integer :: iz, ia
  real(fp) :: r_mass, q_atom, q_ch
 
  !---------------------------------------
 
  nsbeam=0
  abeams=0.0
  xzbeams=0.0
  nlheavy=.false.
 
  nbhydr=0
  nbdeut=0
  nbtrit=0
  nbhel3=0
  nbhel4=0
 
  do ib=1,nbeam
 
     ! have we seen this specie before...
 
     imatch=0
     do isb=1,nsbeam
        if((abs(abeama(ib)-abeams(isb)).le.0.05_fp).and. &
             (abs(xzbeama(ib)-xzbeams(isb)).le.0.05_fp)) then
           imatch=isb
           exit
        endif
     enddo
 
     if(imatch.gt.0) cycle
 
     ! new specie
 
     if(nsbeam.ge.mibs) then
        write(nonlin,*) ' ?nbi_datchkb:  too many distinct beam species.'
        write(nonlin,*) '  max no. of distinct species: ',mibs
        ierct=ierct+1
        return
     endif
 
     nsbeam=nsbeam+1
     abeams(nsbeam)=abeama(ib)
     xzbeams(nsbeam)=xzbeama(ib)
 
     ! classification...
 
     if((abs(abeams(nsbeam)-1.0_fp).le.0.05_fp).and. &
          (abs(xzbeams(nsbeam)-1.0_fp).le.0.05_fp)) then
        nbhydr=nsbeam
     else if((abs(abeams(nsbeam)-2.0_fp).le.0.05_fp).and. &
          (abs(xzbeams(nsbeam)-1.0_fp).le.0.05_fp)) then
        nbdeut=nsbeam
     else if((abs(abeams(nsbeam)-3.0_fp).le.0.05_fp).and. &
          (abs(xzbeams(nsbeam)-1.0_fp).le.0.05_fp)) then
        nbtrit=nsbeam
     else if((abs(abeams(nsbeam)-3.0_fp).le.0.05_fp).and. &
          (abs(xzbeams(nsbeam)-2.0_fp).le.0.05_fp)) then
        nbhel3=nsbeam
     else if((abs(abeams(nsbeam)-4.0_fp).le.0.05_fp).and. &
          (abs(xzbeams(nsbeam)-2.0_fp).le.0.05_fp)) then
        nbhel4=nsbeam
     else if((xzbeams(nsbeam).ge.2.95_fp).and. &
          (abeams(nsbeam).ge.1.95_fp*xzbeams(nsbeam))) then
        nlheavy(nsbeam)=.true.
     else
        ierct=ierct+1
        write(nonlin,*) '?nbi_datchkb:  could not classify beam specie:'
        write(nonlin,*) '   abeams(',nsbeam,')=',abeams(nsbeam)
        write(nonlin,*) '    xzbeams(',nsbeam,')=',xzbeams(nsbeam)
     endif
 
  enddo
! convert from chemical definition to proton mass
  do ib=1,nbeam
     iz=xzbeama(ib) + 0.1_fp
     ia=abeama(ib) + 0.5_fp
     call tr_species_convert(iz,iz,abeama(ib),q_atom,q_ch,r_mass,ier)
     abeama(ib)=r_mass/ps_mp !convert from kg to proton mass
  enddo
  do ib=1,nsbeam
     iz=xzbeams(ib) + 0.1_fp
     ia=abeams(ib) + 0.5_fp
     call tr_species_convert(iz,iz,abeams(ib),q_atom,q_ch,r_mass,ier)
     abeams(ib)=r_mass/ps_mp !convert from kg to proton mass
  enddo
  

 
end subroutine r8_nbi_datchkb
