!******************** start file bbneut.for ; group bbneut *************
!
subroutine bbneut(ierr)
  !
  !  dmc jan 1995 -- adding support for calculation of burnup of rf
  !   tail ions from the dist. fcn. computed in the spruce/fpp package.
  !
  !  dmc oct 1991 -- adding support for d-he3 reaction; adding support
  !   for fusion product burnup:  d-t and d-he3.
  !   for fusion reaction indexing info:  see comments, fusion.fnb
  !
  !  dmc sept 1991 -- the restriction that all fast ion species
  !   distribution functions have the same energy grid is being lifted.
  !   also, the distribution function energy grid no longer has to be
  !   constructed out of steps of constant d(energy).
  !     the inputs delebd, delebda, ebdmax, and ebdmaxa are replaced
  !   with the energy zone bdy array efbmb:
  !     for fast ion species isb, energy zone j covers the range
  !       efbmb(j,isb) to efbmb(j+1,isb)
  !
  !  dmc sept 1991 -- beam dist fcn modified to cover only the plasma;
  !   zones beyond the plasma are no longer computed
  !
  !  dmc june 1990 -- modified to support fpmain (nlbfpp=.true.)
  !  monte-carlo-ize the integral for a continuous fbm (not sparse)
  !  computed from the fpp code
  !  use new routines:
  !    bbnpre -- decide no. of "events" / zone and no. of bins into
  !      which to discretize the distribution function(s).
  !    bbnpik -- do the distribution function discretization; renormalize
  !      to the total no. of ions in the zone.
  !
  !	bob mccann, oct. 1982
  !	rob goldston, oct. 1982 (if then else, elseif, rng)
  !
  !	mod dmc dec. 1986 - output 1d profile of flux surface
  !	averaged beambeam neutrons, for comparison to other neutron
  !       sources
  !
  !       mod dmc aug 1988 - restored old functionality (d-d fus. only)
  !   with new common structures supporting future upgrade for multiple
  !   fusion reactions.
  !
  !   *** correction june 2001 *** using preact now; transp common
  !       (datbuf) no longer used for rate data!
  !xx
  !xx       mod dmc oct 1988 - use trdat common (datbuf) as source of
  !xx   <sigma*v> data instead of defunct bbsigv array.  use flipr output
  !xx   array nfbmpt giving no. of non-zero fbm elements / zone, instead of
  !xx   counting them up again here.
  !
  !	please note: rtm sept 1988 - this calculates the reaction <sigma*v>
  !   integrated over the distribution functions for fast ions.  for reactions
  !   which produce only one neutron, this gives the fusion neutron production
  !   rate.  it also gives the production rate for the other reaction products.
  !   for example, for d(d,n)he3, this gives the production rate for both n and
  !   he3, but, for t(t,2n)he4, it gives the production rate for he4 as well as
  !   half the production rate for n!
  !
  !	this calculates the neutron production from beam-beam
  !   collision reactions using the beam distribution function
  !   which is a function of beam particle energy and pitch
  !   angle.  it also uses the collision cross-section already
  !   integrated over phase angle.  the arguments of this
  !   function are a = vp1**2 + vp2**2 + (vpar2-vpar1)**2 and
  !   b = 2*vp1*vp2, were vp1 and vp2 are the perpendicular
  !   components of the first and second particle velocity and
  !   vpar1 and vpar2 are the parallel components.  the distribu-
  !   tion function may equivalently be thought of as having
  !   v**2 and chi as arguments instead of e and chi.  in these
  !   terms we have the following
  !
  !	vp1**2 = v1**2(1 - chi1**2)
  !	vpar1 = v1*chi1
  !
  !	vp2**2 = v2**2(1 - chi2**2)
  !	vpar2 = v2*chi2
  !
  !	a = v1**2*(1-chi1**2) + v2**2*(1-chi2**2) +
  !		(v2*chi2 - v1*chi1)**2
  !   or
  !	a = v1**2 + v2**2 - 2.0*v1*chi1*v2*chi2
  !
  !	b = 2*sqrt(v1**2*(1-chi1**2)*v2**2*(1-chi2**2) )
  !
  !   in computing the gyrophase integral, e = a-b*cos(phi), but
  !   for convenience i have used  e = a*(1-bf*cos(phi)) and it
  !   is a and bf which are written into the data file where
  !   bf varies from 0.0 to 1.0.
  !   (a is called bba in trans.com and bf is called bbb)
  !
  !	for cases where the target and projectile species are different,
  !   there is a minor change in a, which is basically the kinetic energy
  !   of the projectile in the target's rest frame.  let v1 be the velocity
  !   of the target in the lab frame, and let v2 be the velocity of the
  !   projectile in the lab frame.  also let r = sqrt(m2/m1), where m1 and m2
  !   are the masses (atomic numbers) of the target and projectile species
  !   respectively.  then, in the target's rest frame,
  !
  !	a = v2*v2 + r*v1*r*v1 - 2.0*r*v1*cos1*v2*cos2
  !
  !	b = 2.0*v1*r*sin1*v2*sin2
  !
  !   that is, we replace v1 by r*v1.  b is still an evenly spaced fraction.
  !
  !	the idea of this subroutine is to identify all the a's
  !   and bf's for a given spatial zone and then to sum them over
  !   the paritally integrated collision cross-section and to
  !   then multiply by the appropriate units factor to give the
  !   rate of neutron production.
  !
  !	we assume that fbm has been "flipped" to indexing
  !   by (energy, angle, zone) in flipr.for.  this subroutine
  !   also provides nfbmpt data
  !
  !-----------------------------------------------
  !
  use nbi_com
  use nbi_random
  !
  !  only set in case of "severe" error:
  !============
  implicit none
  integer :: itreact,isb,i,j,ireact,izone,isb1,isb2,ireacl
  integer :: iz,idslas,id1,id2,inum,igot,iproj,ik,idum,ja
  integer :: ie,ilower,itarg,isblcl,izr,ij1,ij2,ij,izz2,izz1,iz1
  integer :: iz2,jzntot,iprod,ilike,jtot2,ii,jmax
  !============
  real(fp) :: zs,zdel,zr2,zr,zebdmax,zdoit,zdoit2,zdelebd,zvol,zjmax
  real(fp) :: zit,zievents
  !============
  integer,intent(out) :: ierr  ! status code, 0=ok
  !
  !-----------------------------------------------
  !	input variables from trans.com
  !-----------------------------------------------
  !
  !	fbm(e,chi,s,species)	beam dist fcn.  by fast ion species
  !                       vs energy(ev), pitch angle, space
  !			units = ptcls*(4*pi/d(omega))/(ev*cm**3)
  !			n.b. fbm has no gyrophase dependence
  !       nfbmpt(s,species)  no. of non-zero elements of fbm in zone s
  !                       for specified species -- from "flipr"
  !
  !	nbbcal		number of beam-beam collisions to calculate
  !			(namelist control)
  !
  !	nznbmr		# beam radial zones
  !       lcm1            lcentr-1 (dummy zone offset)
  !       nthzsm(...)     total number of beam zones to zone radial row...
  !
  !       nfbzns          grand total number of beam zones
  !       ..nfbznsi         no. of beam zones *inside the plasma**
  !          dmc sept 1991 -- use nfbznsi instead of nfbzns
  !
  !	bmvol       	volume of each deposition zone
  !
  !	nznbmea(isb)		# beam particle energies (v**2's)
  !                         now a function of fast ion species
  !
  !	nznbma		# beam pitch angles
  !			(-1 to +1) in nznbma zones
  !       efbmb(j,isb)    lower bound of energy grid zone j for
  !                         fast ion species isb (ev)
  !
  !	abeams(lsbeam)		!atomic weight of beam species
  !	nout		!tty unit number
  !
  !-----------------------------------------------
  !	output information to trans.com
  !-----------------------------------------------
  !
  !	bbnvol(nfbznsi,ireact)	neutron production/2d zone by reaction
  !       bbntns(j,ireact)       neutron production/1d zone by reaction
  !       bbnta(ireact)           beam-beam neutron production by reaction
  !       bbntnt(j)              total neutron production/ 1d zone
  !	bbntot		total beam-beam neutron production
  !
  !-----------------------------------------------
  !	local variable declaration
  !-----------------------------------------------
  !
  integer, save :: izntot
  !			!this is ((#pitch angles)*(#energies))**2
  !
  integer :: izntot_tmp
  real(fp) ::	zdvol		!integration volume element
  !				![delebd/nznbma]**2	!needs c-m fix
  real(fp) ::	zunits		!zdvol*units (*0.5 for self-self collisions)
  !
  real(fp), dimension(:), allocatable, save :: zfpk1,zfpk2, &
       zcos,zsin
  real(fp) :: zcspk1,zsnpk1,zvvpk1,zvpk1,zcspk2,zsnpk2,zvvpk2,zvpk2, &
       za,zb
  !
  !  (converted to allocatable arrays...)
  !x      real	zfpk1(izntot)	!fbm(,,,i1) without zeroes
  !x      real	zfpk2(izntot)	!fbm(,,,i2) without zeroes
  !x      real	za	! a's used
  !x      real	zb	! b's used
  !x      real	zcos(mimxba)	!array of cos(pitch angle)
  !x      real	zsin(mimxba)	!array of sin(pitch angle)
  !
  !x      real	zcspk1	!cos(pitch angle) used in fbm(,,,i1)
  !x      real	zsnpk1)	!ditto - sin(pitch angle) used
  !x      real	zvvpk1	!ditto - v**2's
  !x      real	zvpk1	!ditto - v's
  !
  !x      real	zcspk2	!cos(pitch angle) used in fbm(,,,i2)
  !x      real	zsnpk2	!ditto - sin(pitch angle) used
  !x      real	zvvpk2	!ditto - v**2's
  !x      real	zvpk2	!ditto - v's
  !
  !   the following are used for the interpolation
  !
  !  note dmc june 1990 -- ial and zbbal are also used to hold random
  !  no. sequences when bbneut is called from fpmain -- to monte-carlo-ize
  !  the beam beam neutron integral on a continuous distribution function
  !  see the bbnpik subroutine
  !
  !  dmc -- interpolation code replace with `preact' call
  !   old local arrays removed...
  !
  integer, dimension(:), allocatable, save :: ial, &
       jchpk1,ivvpk1,idist1,jchpk2,ivvpk2,idist2
  real(fp), dimension(:), allocatable, save :: zbbal, &
       zevals,zevars,zsigva,zcont
  !
  !x      integer ::	ial(izntot)	!index of a.le.za -- no:  rng seq info.
  !
  !x      real	zbbal(izntot)	!max a .le.za  -- no:  rng seq info.
  !
  !
  !x      real zevals(izntot),zevars(izntot),zsigva(izntot)
  !
  !   end of interpolation stuf
  !
  !x      real	zcont(izntot)	!neutron production from zsigv
  !
  !x      integer ::	jchpk1(izntot)	!indices for v**2's used
  !x      integer ::	ivvpk1(izntot)	!ditto - vvi**2's
  !x      integer ::	idist1(mimxbz)	!# nonzero elements in zfpk1
  !				!jchpk1, and ivvpk1
  !x      integer ::	jchpk2(izntot)	!indices for v**2's used
  !x      integer ::	ivvpk2(izntot)	!ditto - vvi**2's
  !x      integer ::	idist2(mimxbz)	!# nonzero elements in zfpk2
  !				!jchpk2, and ivvpk2
  logical	ilsame		!t=> target=projectile d(d,n)he3
  !
  !
  !--------------------------------------------
  real(fp) :: fbm_flr !function to choose distribution function
  ! depending on nlfbmflr, fbm if nlfbmflr=.false.
  ! fbm_ptcl if nlfbmflr=.true.
  !
  !  table of neutron producting fusion reactions
  !  see the comments in routine fusion.fnb for description of
  !  reaction indexing
  !
  parameter (itreact=14)
  integer :: ineuts(itreact)
  !
  data ineuts/1,0,0,1,1,1,1,0,0,0,0,0,0,0/
  !
  !  0's for computed fusion reactions that do not produce neutrons
  !  also a 0 for the fusion-product burnup neutrons category (#9)
  !  and the rf tail burnup neutrons categories (#11 and #13)
  !  which are counted separately cf sumntx.for (dmc 17 oct 1991, & 1997)
  !
  !--------------------------------------------
  !
  !===============================================
  !
  !
  !	----------------------
  !	0.1	initialization
  !	----------------------
  !
  ierr = 0
  !

  izntot=mimxba*mimxbe
  izntot_tmp= min(int(sqrt(real(nbbcal,fp)))+1,100000)
  izntot=max(izntot,izntot_tmp)
  nbbcal=izntot*izntot

  !
  !  allocate and initialize bbneut arrays
  !

  call dealloc_arr
  allocate(zfpk1(izntot),zfpk2(izntot))
  allocate(zcos(mimxba),zsin(mimxba))
  zfpk1=0
  zfpk2=0
  za=0; zb=0
  zcos=0; zsin=0
  zcspk1=0; zsnpk1=0; zvvpk1=0; zvpk1=0
  zcspk2=0; zsnpk2=0; zvvpk2=0; zvpk2=0

  allocate(ial(izntot))
  allocate(jchpk1(izntot),ivvpk1(izntot),idist1(mimxbz))
  allocate(jchpk2(izntot),ivvpk2(izntot),idist2(mimxbz))
  allocate(zevals(izntot))
  allocate(zevars(izntot),zsigva(izntot),zcont(izntot))
  ial=0;
  jchpk1=0; ivvpk1=0; idist1=0
  jchpk2=0; ivvpk2=0; idist2=0
  zevals=0
  zevars=0; zsigva=0; zcont=0
  !
  !
  !---------------------------------
  !
  !   now define constant part of volume element (dv)**6
  !   note that i have already integrated once over relative
  !   gyrophase when i created <sigma-v>.  the asolute phase
  !   then integrates to another 2*pi.
  !
  !   so, initially for the usual distribution function, what i
  !   want is
  !
  !	[ 2*pi*v**2*dv*d(cos(th)) ]**2
  !   or
  !	[ pi*v*d(v*v)*d(cos(th)) ]**2
  !
  !   however, rob has used fancy units for fbm, so zdvol is
  !   much simpler, and i don't have to have the v in (dv)**3
  !   also, i don't need the pi in (dv)**3.  that just leaves
  !
  !	[d(v*v)*d(cos(th))]**2   ...   i hope   ...
  !
  zdvol = 1.0_fp/(real(nznbma,fp)**2)	!needs c-m fix
  !
  zdvol = zdvol
  !
  !   need array of sqrt(energies) for all species
  !
  !
  zs = -1.0_fp- 1.0_fp/real(nznbma,fp)
  zdel = 2.0_fp/real(nznbma,fp)
  do i = 1, nznbma
     zcos(i) = zs + zdel*i
     zsin(i) = sqrt(1.0_fp- zcos(i)**2)
  end do
  !
  !  dmc aug 1988 -- also zero output arrays indexed by fusion reaction
  !    nreact = number of known reactions.
  !  dmc oct 1991 -- nreact=12 -- see comments in fusion.fnb
  !
  bbntnt(:)=0.0_fp
  bbntns=0.0_fp
  !
  bbnvol(1:nfbznsi,1:nreact) = 0.0_fp !neutron production for zone i
  !
  bbntot = 0.0_fp		!total neutron production
  !
  bbnta(1:nreact)=0.0_fp
  !
  !	  	at this point i need to have two pointer vectors
  !		for the nonzero entries in fbm for each species
  !
  !
  !
  !	-----------------------------------------
  !	a.0	loop over all beam-beam reactions
  !	-----------------------------------------
  !
  !
  do  ireact = 1, nreact
     if(.not.nlbba(ireact)) cycle !only do the ones we want
     if(ireact.ge.11) cycle ! neutrons from rf species are calculated in btfusn
     isb1 = nbbsp(ireact,1)	!transp species index for target
     isb2 = nbbsp(ireact,2)	!transp species index for projectile
     if(norbtot.eq.0) then
        !     skip if no mc orbits for mc species
        if(.not.nlfbmfpp(isb1)) cycle
        if(.not.nlfbmfpp(isb2)) cycle
     end if
     !
     !
     ireacl=0               ! <sigma*v> data selector
     !
     !				!check for equality of target & projectile
     !
     if (isb1.eq.isb2) then
        !
        !  symmetric reaction
        !
        ilsame = .true.
        zr2 = 1.0_fp
        zr = 1.0_fp
        zunits = zdvol*0.5_fp !beam species = target species only
        !
        if(isb1.eq.nbtrit) then !t(t,2n)he3
           !
           !  ** tt reaction ** 0.5 [self collision] * 2.0 [ntns/reaction] = 1.0
           !
           zunits = zdvol   !(beam=target, but twice as many neutrons)
        end if
        !
        !  dmc setup to access tables
        if(isb1.eq.nbdeut) then
           if(ireact.eq.3) then
              !  dd proton branch
              ireacl=1
           else
              !  dd neutron branch
              ireacl=2
           end if
        else if(isb1.eq.nbtrit) then
           ireacl=3         ! tt reaction
        end if
        !
     else
        !
        !  asymmetric reaction
        !
        ilsame = .false.
        zunits = zdvol
        !
        zr2 = abeams(isb2)/abeams(isb1)
        zr = sqrt(zr2)      !ratio for v1*r used in computing "a"
        !
        !  dmc setup to access tables
        if((isb2.eq.nbdeut).and.(isb1.eq.nbtrit)) then
           !  dt beam-beam
           ireacl=4
        else if((isb2.eq.nbdeut).and.(isb1.eq.nbhel3)) then
           !  d-he3 beam-beam
           ireacl=5
        else if((isb2.eq.nbdeut).and.(isb1.eq.nftrit)) then
           !  dt fusion-product-on-beam
           ireacl=4
        else if((isb2.eq.nbdeut).and.(isb1.eq.nfhel3)) then
           !  d-he3 fusion-product-on-beam
           ireacl=5
        else if((isb2.eq.nbdeut).and.(isb1.eq.nrftrit)) then
           !  dt t-rf-tail-on-d-beam
           ireacl=4
        else if((isb2.eq.nbdeut).and.(isb1.eq.nrfhel3)) then
           !  d-he3 he3-rf-tail-on-d-beam
           ireacl=5
        else if((isb2.eq.nbtrit).and.(isb1.eq.nrfdeut)) then
           !  dt d-rf-tail-on-t-beam
           ireacl=4
        else if((isb2.eq.nbhel3).and.(isb1.eq.nrfdeut)) then
           !  d-he3 d-rf-tail-on-he3-beam
           ireacl=5
        end if
     end if
     !
     !  dmc -- set up bba and bbb arrays for *this* reaction
     !
     if(ireacl.eq.0) then

        write(nonlin,9901) ireact
9901    format(' ? bbneut -- reaction rate data not known, ireact=',i2)
        ierr=1
        call dealloc_arr
        return
     else
        !
        if(isb1.eq.isb2) then
           zebdmax=4.0_fp*efbmb(nznbmea(isb1),isb1)
        else
           zebdmax= &
                efbmb(nznbmea(isb1),isb1)+efbmb(nznbmea(isb2),isb2)+ &
                2.0_fp*sqrt(efbmb(nznbmea(isb1),isb1)* &
                efbmb(nznbmea(isb2),isb2))
        end if
        !
     end if
     !
     !   issue an error if nbbcal is large.  when i wrote this, nbbcal=1.0e5
     !   and izntot**2 = 50**4 = 6.25e6, so this warning shouldn't happen unless
     !   someone gets greedy with nbbcal!
     !
     if(nbbcal.gt.izntot**2) then
        jtot2 = izntot**2
        write(nout,80501) nbbcal, jtot2
        zjmax = 0
        do ii = 1, nfbznsi
           zit = nfbmpt(ii,isb1)*nfbmpt(ii,isb2)
           if(zit.gt.zjmax) then
              zjmax = zit
           end if
        end do
        !
        jmax = zdoit*sqrt(zjmax) + 1
        write(nout,80511) jmax
        ierr=1
        call dealloc_arr
        return
     endif
     !
     !
     ! calculate bdens2 and nfbmpt for isb1, isb2
     call bbnpre_bdens2(isb1,isb2)

     !   count the number of events over the whole plasma
     !
     if(nlfbmfpp(isb1).or.nlfbmfpp(isb2)) then
        !
        !  from fpp -- decide on fbm discretization in each zone, in advance
        !    bbnpre returns answer in common array nfbmpt(*,isb1)
        !
        call bbnpre_nevents(isb1,isb2)
        !
        zdoit=0.0_fp
        zdoit2=0.0_fp
        !
     else
        !
        !  from monte carlo beam code -- fbm is already discrete and sparse
        !  count the no. of non-zero points
        !
        zievents = 0.0_fp
        do iz = 1, nfbznsi
           zievents = zievents + nfbmpt(iz,isb1)*nfbmpt(iz,isb2)

        end do
        !
        !	set up the control for the monte-carlo integration
        !	and adjust zunits to take m-c into account
        !
        !       bail out if there are no events
        !
        if(zievents.le.0.0_fp) cycle
        !
        zdoit2 = real(nbbcal,fp)/zievents !events_used/total_events
        zdoit2 = min(zdoit2, 1.0_fp) !1.0 means "use them all"
        zdoit  = sqrt(zdoit2) !sqrt(events) = # fbm elements
        !			  !zdoit is the fraction of nonzero fbm elements to use
        zunits=zunits/(zdoit2)
     end if
     !
     !	-----------------------------------
     !	a.1	loop over all spatial zones
     !	-----------------------------------
     !
     !
     do izone = 1, nfbznsi
        !	  	!# nonzero elements of fbm/zone actually used
        !               !hope that this is .le. izntot
        idist1(izone) = 0   !target loops
        idist2(izone) = 0   !projectile loops
        !
        !
        !	--------------------------------------
        !	1.0	create the list of a's and b's
        !		for this spatial zone
        !	--------------------------------------
        !
        !
        !
        !		initialize the pointers, etc.  (this is to help debugging)
        !
        if(izone.eq.1) then
           idslas=izntot	!initialize all the first time through
        else                !re-initialize only those last used
           id1 = idist1(izone-1)
           id2 = idist2(izone-1)
           idslas=max(id1,id2)
        end if
        !
        do i = 1, idslas
           ivvpk1(i) = 0    !points to v**2's used in fbm(,,,i1)
           jchpk1(i) = 0    !points to chi's used
           zfpk1(i) = 0.0_fp
           !
           ivvpk2(i) = 0    !points to v**2's used in fbm(,,,i2)
           jchpk2(i) = 0    !points to chi's used
           zfpk2(i) = 0.0_fp
           !
        end do
        !
        !		squeeze the zeroes out of the dist. function
        !
        if(ilsame) then     !identical species
           !  symmetric reaction
           !
           if(nlfbmfpp(isb1)) then
              !  fokker planck beam code:
              !  make a discretized approximation of the local dist. fcn.
              inum=nfbmpt(izone,isb1)
              if(inum.eq.0) then
                 idist1(izone)=0
                 idist2(izone)=0
                 cycle ! dist fcn empty
              end if
              call bbnpik(izone,isb1,inum, &
                   ial,zbbal,zfpk1,ivvpk1,jchpk1,igot)
              if(igot.eq.-1) then
                 write(nonlin,*) ' bbneut -- error exit, bbnpik error.'
                 write(nonlin,*) '   ireact = ',ireact,'  izone = ',izone
                 write(nonlin,*) '   isb = ',isb1 ! isb1 = isb2, symm. reaction.
                 ierr=1
                 call dealloc_arr
                 return
              end if
              idist1(izone)=igot
              !  and copy (symmetric reaction)
              idist2(izone)=idist1(izone)
              do iproj=1,idist2(izone)
                 zfpk2(iproj)=zfpk1(iproj)
                 ivvpk2(iproj)=ivvpk1(iproj)
                 jchpk2(iproj)=jchpk1(iproj)
              end do
              !
           else
              !
              !  symmetric reaction
              !  monte carlo beam code:
              !  pack zeroes out of the m.c. distribution fcn; use russian roulette
              !  factor as well
              do j = 1, nznbma
                 do ik = 1, nznbmea(isb1)
                    if (fbm_flr(ik,j,izone,isb1,0).eq.0.0_fp) then
                       cycle !don't do anything else
                    else if(nbrandm(idum).gt.zdoit) then
                       cycle !don't do anything else
                    else
                       idist1(izone) = idist1(izone) + 1 !target loop
                       idist2(izone) = idist1(izone) !projectile loop
                       if(idist1(izone).gt.izntot) then
                          jzntot = izntot
                          iprod = nfbmpt(izone,isb1)*nfbmpt(izone,isb2)
                          ilike = iprod*zdoit2
                          write(nout,80401) jzntot, izone, isb1, isb2, &
                               nfbmpt(izone,isb1), nfbmpt(izone,isb2), &
                               iprod, ilike, nbbcal, zievents, zdoit2
                          ierr=1
                          call dealloc_arr
                          return       
                       endif
                       !
                       zdelebd=efbmb(ik+1,isb1)-efbmb(ik,isb1)
                       zfpk1(idist1(izone)) = fbm_flr(ik,j,izone,isb1,0)*zdelebd
                       zfpk2(idist1(izone)) = fbm_flr(ik,j,izone,isb1,0)*zdelebd
                       !
                       ivvpk1(idist1(izone)) = ik !pointers to e's
                       ivvpk2(idist1(izone)) = ik
                       !
                       jchpk1(idist1(izone)) = j !pointers to v_par/v's
                       jchpk2(idist1(izone)) = j
                    end if
                    !
                 end do
              end do
           end if
           !
           if(idist1(izone).eq.0) cycle!nothing to integrate in this zone
           !
        else                !different species
           !
           !  assymetric reaction
           !
           if(nlfbmfpp(isb1).or.nlfbmfpp(isb2)) then
              !  fokker planck beam code:
              !  make a discretized approximation of the local dist. fcns.
              inum=nfbmpt(izone,isb1)
              if(inum.eq.0) then
                 idist1(izone)=0
                 idist2(izone)=0
                 cycle ! dist fcn empty
              end if
              !
              call bbnpik(izone,isb1,inum, &
                   ial,zbbal,zfpk1,ivvpk1,jchpk1,igot)
              if(igot.eq.-1) then
                 write(nonlin,*) ' bbneut -- error exit, bbnpik error.'
                 write(nonlin,*) '   ireact = ',ireact,'  izone = ',izone
                 write(nonlin,*) '   isb1 = ',isb1
                 ierr=1
                 call dealloc_arr
                 return
              end if
              idist1(izone)=igot
              !
              call bbnpik(izone,isb2,inum, &
                   ial,zbbal,zfpk2,ivvpk2,jchpk2,igot)
              if(igot.eq.-1) then
                 write(nonlin,*) ' bbneut -- error exit, bbnpik error.'
                 write(nonlin,*) '   ireact = ',ireact,'  izone = ',izone
                 write(nonlin,*) '   isb2 = ',isb2
                 ierr=1
                 call dealloc_arr
                 return
              end if
              idist2(izone)=igot
           else
              !
              !  monte carlo beam code
              !  pack non-zero elements of dist. fcn together; use russian roulette
              !  to control statistics - runtime
              do ja = 1, nznbma
                 do ie = 1, nznbmea(isb1)
                    if(fbm_flr(ie,ja,izone,isb1,0).eq.0.0_fp) then !do nothing
                       !
                    else if(nbrandm(idum).gt.zdoit) then !throw away enough, i hope
                       !						       !don't do anything else
                    else    !record fbm element
                       idist1(izone) = idist1(izone) + 1
                       if(idist1(izone).gt.izntot) then
                          jzntot = izntot
                          iprod = nfbmpt(izone,isb1)*nfbmpt(izone,isb2)
                          ilike = iprod*zdoit2
                          write(nout,80401) jzntot, izone, isb1, isb2, &
                               nfbmpt(izone,isb1), nfbmpt(izone,isb2), &
                               iprod, ilike, nbbcal, zievents, zdoit2
                          ierr=1
                          call dealloc_arr
                          return       
                       endif
                       zdelebd=efbmb(ie+1,isb1)-efbmb(ie,isb1)
                       zfpk1(idist1(izone)) = fbm_flr(ie,ja,izone,isb1,0)*zdelebd
                       ivvpk1(idist1(izone)) = ie
                       jchpk1(idist1(izone)) = ja
                    end if
                    !
                 end do
              end do
              !
              if(idist1(izone).eq.0) cycle !nothing to integrate in this zone
              !
              do ja = 1, nznbma
                 do ie = 1, nznbmea(isb2)
                    if(fbm_flr(ie,ja,izone,isb2,0).eq.0.0_fp) then !do nothing
                       !
                    else if(nbrandm(idum).gt.zdoit) then !throw away enough, i hope
                       !						      !don't do anything else
                    else    !record fbm element
                       idist2(izone) = idist2(izone) + 1
                       if(idist2(izone).gt.izntot) then
                          jzntot = izntot
                          iprod = nfbmpt(izone,isb1)*nfbmpt(izone,isb2)
                          ilike = iprod*zdoit2
                          write(nout,80401) jzntot, izone, isb1, isb2, &
                               nfbmpt(izone,isb1), nfbmpt(izone,isb2), &
                               iprod, ilike, nbbcal, zievents, zdoit2
                          ierr=1
                          call dealloc_arr
                          return       
                       endif
                       zdelebd=efbmb(ie+1,isb2)-efbmb(ie,isb2)
                       zfpk2(idist2(izone)) = fbm_flr(ie,ja,izone,isb2,0)*zdelebd
                       ivvpk2(idist2(izone)) = ie
                       jchpk2(idist2(izone)) = ja
                    end if
                    !
                 end do
              end do
              !
              if(idist2(izone).eq.0) cycle !nothing to integrate in this zone
              !
           end if            ! fpp / monte carlo code
           !
        end if               ! symmetric / assymetric reaction
        !
        !		store the v**2's and chi's
        !		this is a "gather" operation
        !		so use the omnilib routine on the cray
        !
        !
        !	-----------------------------------------
        !	b.0	now double loop over the dist fnc
        !		to convolve it with itself
        !	-----------------------------------------
        !
        !   i'm assuming a very dumb vector optimizer
        !   so i'll parse the expressions myself so that
        !   it all vectorizes
        !
        do iproj = 1, idist2(izone) !projectile loop
           !
           !		now form the args for the gyrophase
           !		averaged production cross-section
           !
           if(ilsame) then  !same species
              ilower = iproj
           else             !different species
              ilower = 1
           end if

           zvvpk2 = efbm(ivvpk2(iproj),isb2)
           zvpk2  = sqrt(zvvpk2)
           zcspk2 = zcos(jchpk2(iproj))
           zsnpk2 = zsin(jchpk2(iproj))
           !
           !		***note:(itarg=iproj symmetry,
           !			 so do half & double off diag part)***
           !
           do itarg = ilower, idist1(izone) !target loop

              zvvpk1 = efbm(ivvpk1(itarg),isb1)
              zvpk1  = sqrt(zvvpk1)
              zcspk1 = zcos(jchpk1(itarg))
              zsnpk1 = zsin(jchpk1(itarg))
              !
              !   a = zr2*v1**2 + v2**2 - 2*zr*v1*cos(th1)*v2*cos(th2)
              !
              za = (-2.0_fp)*zvpk2 !pk2 is iproj(projectils)
              za = za*zvpk1*zr !pk1 is itarg (target)
              za = za*zcspk2
              za = za*zcspk1
              za = za + zvvpk2
              za = za + zvvpk1*zr2
              !
              !   b = 2*zr*v1*sin(th1)*v2*sin(th2)
              !
              zb = 2.0_fp*zvpk2
              zb = zb*zsnpk2
              zb = zb*zvpk1*zr
              zb = zb*zsnpk1
              !
              !   now normalize zb to za since bbb is a fraction (0 to 1)
              !   also force zb to be .lt. 1.0
              !   ***n.b. here i assume za.ne.0!  that should be true given
              !           the error testing in "0.1".
              !
              zevars(itarg)=min(1.0_fp,(zb/za))
              !
              !  convert to kev
              !
              zevals(itarg)=za*0.001_fp ! kev
           end do
           !
           !
           !	---------------------------------------------
           !	2.0	linearly interpolate the phase avg cs
           !	---------------------------------------------
           !
           !
           zsigva=0.0_fp
           if(ireacl.eq.1) then
              call fusn_bbsigv_ddp(zevals(ilower:idist1(izone)), &
                   zevars(ilower:idist1(izone)),zsigva(ilower:idist1(izone)), &
                   idist1(izone)-ilower+1)
           else if(ireacl.eq.2) then
              call fusn_bbsigv_ddn(zevals(ilower:idist1(izone)), &
                   zevars(ilower:idist1(izone)),zsigva(ilower:idist1(izone)), &
                   idist1(izone)-ilower+1)
           else if(ireacl.eq.3) then
              if(lev_ttnuk.eq.0) then
                 call fusn_bbsigv_tt(zevals(ilower:idist1(izone)), &
                      zevars(ilower:idist1(izone)), &
                      zsigva(ilower:idist1(izone)), &
                      idist1(izone)-ilower+1)
              elseif(lev_ttnuk.eq.1) then
                 call fusn_bbsigv_dpctt &
                      (zevals(ilower:idist1(izone)), &
                      zevars(ilower:idist1(izone)), &
                      zsigva(ilower:idist1(izone)), &
                      idist1(izone)-ilower+1)
              end if
           else if(ireacl.eq.4) then
              call fusn_bbsigv_dt(zevals(ilower:idist1(izone)), &
                   zevars(ilower:idist1(izone)),zsigva(ilower:idist1(izone)), &
                   idist1(izone)-ilower+1)
           else if(ireacl.eq.5) then
              call fusn_bbsigv_d3(zevals(ilower:idist1(izone)), &
                   zevars(ilower:idist1(izone)),zsigva(ilower:idist1(izone)), &
                   idist1(izone)-ilower+1)
           else if(ireacl.eq.6) then
              call fusn_bbsigv_t3(zevals(ilower:idist1(izone)), &
                   zevars(ilower:idist1(izone)),zsigva(ilower:idist1(izone)), &
                   idist1(izone)-ilower+1)
           end if
           zsigva(ilower:idist1(izone))=1.e6_fp* &
                zsigva(ilower:idist1(izone))
           !
           !	-------------------------------------------
           !	3.0	sum the contributions for this zone
           !	-------------------------------------------
           !
           !
           !
           if(ilsame) then  !same species
              !
              !   now add the diagonal contribution
              !
              zcont(iproj) = zsigva(iproj)*zfpk1(iproj)**2
              bbnvol(izone,ireact) = bbnvol(izone,ireact) + &
                   zcont(iproj)
              !
              !   next add twice the off diagonal contribution to take
              !   care of the symmetry
              !
              if(iproj.ne.idist1(izone)) then	!some compilers aren't ansi
              !
                 do itarg = iproj+1, idist1(izone) !target loop
                    zcont(itarg) = zsigva(itarg)*zfpk1(iproj)
                    zcont(itarg) = zcont(itarg)*zfpk1(itarg)
                    zcont(itarg) = 2.0_fp*zcont(itarg)
                    !
                    !   here i'm using a lowest order integration scheme
                    !   by just adding rectangular parallelpipeds
                    !
                    bbnvol(izone,ireact) = bbnvol(izone,ireact) + &
                         zcont(itarg)
                 end do
              endif
              !
           else             !different species
              !
              do itarg = 1, idist1(izone) !target loop
                 zcont(itarg) = zsigva(itarg)*zfpk1(itarg)
                 zcont(itarg) = zcont(itarg)*zfpk2(iproj)
                 bbnvol(izone, ireact) = bbnvol(izone, ireact) &
                      + zcont(itarg)
              end do
              !
           end if
           !
        end do            !end outer loop of convolution
        !
     end do               !end loop over zones
     !
     !
     !	------------------------------
     !	4.0	sum over all the zones and fix units
     !	------------------------------
     !
     !
     !   first take care of the deposition volumes
     !
     do izone = 1, nfbznsi
        bbnvol(izone,ireact) = &
             (bbnvol(izone,ireact)*zunits)*bmvol(izone)
        !
        !   now add all the zones together
        !
        bbnta(ireact) = bbnta(ireact) + bbnvol(izone,ireact)
     end do
     !
     !
     !	---------------------------------------------
     !	6.0	end the loop over beam-beam reactions
     !	---------------------------------------------
  end do                  !ireact=1,nreact loop
  !
  !
  !	------------------------
  !	7.0	cleanup & return
  !	------------------------
  !
  !
  !
  !  dmc dec 1986:
  !  bb neuts from each radial zone (orig 1d code grid)
  !
  do ireact=1,nreact
     if(.not.nlbba(ireact)) cycle !only do the ones we need
     !
     do izr=1,nznbmr
        !
        ij1=lcentr+(izr-1)*nxskpb
        ij2=ij1+nxskpb-1
        !
        if(ij1.gt.ledge) cycle
        !  volume
        zvol=0.0_fp
        do ij=ij1,ij2
           zvol=zvol+dvolzc(ij)
        end do
     !
     !  sum bbntns (zero initialized above)
     !  beam zone indices
        izz2=izr+lcm1
        izz1=izz2-1
        iz1=nthzsm(izz1)+1
        iz2=nthzsm(izz2)
     !
        do iz=iz1,iz2
           do ij=ij1,ij2
              bbntns(ij,ireact)=bbntns(ij,ireact) + &
                   bbnvol(iz,ireact)*dvolzc(ij)/zvol
           end do
        end do
  !
     end do               ! izr loop
!
  end do                  ! ireact loop
!
!  dmc aug 1988 -- sum over reactions
!    sums zeroed above at top of routine
!
  do ireact=1,itreact
     if(.not.nlbba(ireact)) cycle !only do the ones we need
     bbntot=bbntot+ineuts(ireact)*bbnta(ireact)
     do j=lcentr,ledge
        bbntnt(j)=bbntnt(j)+ineuts(ireact)*bbntns(j,ireact)
     end do
  end do
!
!
  call dealloc_arr
  return
!
!
!	-----------------------------------
!	8.0	error handling and messages
!	-----------------------------------
!
!
!--------------------------
!
80401 format(' bbneut error: too many collisions. izntot=',i5, &
           ' dimension exceeded izone=',i5,/, &
           ' target species=',i2,' projectile species=',i2,/, &
           ' #nonzero target elements=',i6,/, &
           ' #nonzero projectile elements=',i6,/, &
           ' product=max # collisions to compute=',i9,/, &
           ' most likely # collisions to compute=',i9, &
           ' nbbcal (# to use over whole plasma)=',i9,/, &
           ' zievents (# potential events overall)=',e12.4,/, &
           ' zdoit2 = nbbcal/zievents = ',1p,e12.4)
!

!
80501 format(' bbneut warning: nbbcal=',i9,'.gt.izntot**2=',i9,/, &
           ' subscript out of bounds possible')
!
80511 format(' bbneut: maximum expected # fbm elements used=',i5)
  !
  !
contains
  
  subroutine dealloc_arr
    
    if(allocated(zfpk1)) then
       deallocate(zfpk1,zfpk2,zcos,zsin)
       deallocate(ial,jchpk1,ivvpk1,idist1)
       deallocate(jchpk2,ivvpk2,idist2)
       deallocate(zevals,zevars,zsigva,zcont)
    end if

end subroutine dealloc_arr

end subroutine bbneut
!******************** end file bbneut.for ; group bbneut ***************
! 11sep2004 fgtok
