subroutine rdbbsg
  use nbi_com

  !
  !  read the beam-beam fusion neutron gyro-phase avg sigma*v file
  !   (of use only if beam dist fcn is evaluated)
  !   note dmc -- july 1990 -- reaction rate data is now read in
  !   through the trdat interface -- set up indexing in rdbbsg
  !
  !   mod dmc oct 1991 -- support addl fusion reactions in bbneut
  !   support fusion triton and he3 burnup calculations as options
  !
  !============
  implicit none
  integer :: ir
  !============
  call catfus_num(nreact)
  !
  if(nlfbm) then
     nlbba(1:mireact)=.false.
     nlbba(1:6)=.true.
     nlbba(9:nreact)=.true.
     !  ... don't want d->d (p,t) unless following fusion tritons
     if(.not.nlddp) nlbba(3)=.false.
     call rdbbsg0
  end if

  return
end subroutine rdbbsg
!******************** start file rdbbsg.for ; group bbneut *************
!
!----------------------------------------------------------
!  modified dmc 13 oct 1988
!    file functionality is superseded -- trdat loads transp common with
!    beam-beam sigma*v gyro integral table data, so the file reading
!    previously done in this routine is no longer needed.  only the
!    function of setting up the bbneut indexing is preserved
!
subroutine rdbbsg0
  !
  !
  !    set up reaction indexing info for bbneut -- beam-beam fusion
  !  calculation.  mod dmc oct 1991 -- also support fusion products
  !  burnup.  special nreact values 9 and 10 are used for t->d and
  !  he3->d fusion product burnup reactions, respectively -- see
  !  comments in fusion.fnb
  !
  !
  !	the general case will be to input a whole set of <sigma*v>'s for
  !   the beam-beam ion reactions listed in harry towner's sigma routine:
  !
  !...n - tells what reaction you want
  !	  =1 you get the cross section for the reaction  t(d,n)he4
  !	  =2 you get the cross section for the reaction  he3(d,p)he4
  !	  =3 you get the cross section for the reaction  d(d,p)t
  !	  =4 you get the cross section for the reaction  d(d,n)he3
  !	  =5 you get the cross section for the reaction  t(t,2n)he4
  !	  =6 you get the cross section for the reaction  t(he3,x)y
  !
  !   nbbsp(reaction,1) = target species in fbm(energy,pitch,zone,species)
  !   nbbsp(reaction,2) = impacting species in fbm(energy,pitch,zone,species)
  !
  !   for example: if species=1 is h and species=2 is d, and species=3 is t, then
  !		 for the first reaction, t(d,n)he3,
  !
  !			nbbsp(1,1) = 3	!species=1 is t which is the target
  !			nbbsp(1,2) = 2	!species=2 is d which is the projectile
  !
  !		and for the fourth reaction, d(d,n)he3,
  !
  !			nbbsp(4,1) = 2
  !			nbbsp(4,2) = 2
  !
  !   to figure out which species is which, we need both z and a for each species
  !
  !	species		z	a	standard order subscript
  !
  !	   h		1	1		1
  !	   d		1	2		2
  !	   t		1	3		3
  !	  he3		2	3		4
  !	  he4		2	4		5
  !
  use nbi_com
  !
  !  cray only -- filename
  !@      character*8 zfile
  !@      data zfile/'allsigv'/
  !
  !--------------------------------------------
  !	input variables from trans.com
  !--------------------------------------------
  !
  !	lunphd	unit number for reading gyroavg beam**2 cs.
  !
  !--------------------------------------------
  !	output variables to trans.com
  !--------------------------------------------
  !
  !
  !	npbba	#of bba values read
  !	npbbb	#of bbb values read
  !	bba	vector of relative energies for collisions
  !	bbb	fraction from 0 to 1 (sin(pitch angle))
  !	bbsigv(bba,bbb)	beam-beam neutron production cross
  !			section*v, gyrophase integrated
  !			units=barns*cm/sec
  !
  !
  !	---------------
  !	local variables
  !	---------------
  !
  !============
  implicit none
  integer :: inspec,ireact,instd,irstd,i,j
  !============
  parameter	(inspec=5)	!number of species in std reactions
  parameter	(ireact=6)	!number of standard reactions
  !
  real(fp) ::		zstdz(inspec)	!z for standard order list
  real(fp) ::		zstda(inspec)	!a for standard order list
  integer ::		iperm(inspec)	!permutation from standard to transp
  !				!if transp species 1 is t, then iperm(3)=1
  integer ::		itarg(ireact)	!std target species
  integer ::		iproj(ireact)	!std projectile species
  !
  character(len=20) ::	zsreact(ireact)	!standard reaction names
  !
  data	zstdz/1.0_fp,1.0_fp,1.0_fp,2.0_fp,2.0_fp/!sorted into std order (h,d,t,he3,he4)
  data	zstda/1.0_fp,2.0_fp,3.0_fp,3.0_fp,4.0_fp/
  data	itarg/3,4,2,2,3,3/	!std target species
  data	iproj/2,2,2,2,3,4/	!std projectile species
  !
  data	zsreact/'t(d,n)he4','he3(d,p)he4','d(d,p)t', & !standard reaction &
       'd(d,n)he3','t(t,2n)he4','t(he3,x)y'/           !names
  !
  !==========================================================
  !
  !
  !	----------------------
  !	0.1	initialization
  !	----------------------
  !
  !
  instd = inspec		!this is to make it possible to look at these
  irstd = ireact		!during debugging
  !
  !
  !	--------------------------------------
  !	1.0	open & read the <sigma*v> file
  !	--------------------------------------
  !       deleted dmc oct 1988
  !
  !
  !
  !
  !	--------------------------------------------------------------------
  !	2.0	create the permutation from transp species to standard order
  !	--------------------------------------------------------------------
  !
  iperm(1:instd) = 0
  !
  do i = 1, nsbeam	!loop over the transp beam species
     if(nlfprod(i)) cycle  !ignore fusion product ions!
     do j = 1, instd	!loop over the std list
        !
        if(abs(xzbeams(i)-zstdz(j)) .gt. 0.05_fp) cycle
        if(abs(abeams(i)-zstda(j)) .gt. 0.05_fp) cycle
        !
        iperm(j) = i
     end do
  end do
  !
  !
  !	-----------------------------------------
  !	3.0	use the permutation to fill nbbsp
  !	-----------------------------------------
  !
  nbbsp(1:irstd,1:2) = 0
  do i = 1, irstd                  !loop over standard reactions
     if(nlbba(i)) then	    ! only if reaction is under consideration
        nbbsp(i,1) = iperm(itarg(i))	!the transp species target
        nbbsp(i,2) = iperm(iproj(i))	!the transp species projectile
     end if
  end do
  !
  do i = 1, irstd
     if(nbbsp(i,1)*nbbsp(i,2).eq.0) then
        nlbba(i) = .false.
     else
        nlbba(i) = .true.
     end if
  end do
  !
  !
  !  indexing for fusion product burnup reactions...
  !    burnup of fusion product on beam ions -- fast ion + fast ion
  !    fusion
  !
  if((nftrit.gt.0).and.(nbdeut.gt.0)) then
     nlbba(9)=.true.
     nbbsp(9,1)=nftrit
     nbbsp(9,2)=nbdeut
  else
     nlbba(9)=.false.
  end if
  !
  if((nfhel3.gt.0).and.(nbdeut.gt.0)) then
     nlbba(10)=.true.
     nbbsp(10,1)=nfhel3
     nbbsp(10,2)=nbdeut
  else
     nlbba(10)=.false.
  end if
  !
  !  indexing for rf minority burnup reactions...
  !    burnup of rf tail ions on beam ions -- fast ion + fast ion
  !    fusion
  !
  if((nbdeut.gt.0).and.(nrftrit.gt.0)) then
     nlbba(11)=.true.
     nbbsp(11,1)=nrftrit
     nbbsp(11,2)=nbdeut
  else
     nlbba(11)=.false.
  end if
  !
  if((nbdeut.gt.0).and.(nrfhel3.gt.0)) then
     nlbba(12)=.true.
     nbbsp(12,1)=nrfhel3
     nbbsp(12,2)=nbdeut
  else
     nlbba(12)=.false.
  end if
  !
  if((nbtrit.gt.0).and.(nrfdeut.gt.0)) then
     nlbba(13)=.true.
     nbbsp(13,1)=nrfdeut
     nbbsp(13,2)=nbtrit
  else
     nlbba(13)=.false.
  end if
  !
  if((nbhel3.gt.0).and.(nrfdeut.gt.0)) then
     nlbba(14)=.true.
     nbbsp(14,1)=nrfdeut
     nbbsp(14,2)=nbhel3
  else
     nlbba(14)=.false.
  end if
  return
end subroutine rdbbsg0
!******************** end file rdbbsg.for ; group bbneut ***************
