!-----------------------------------------------------------------------
!  fusion -- compute fusion product spatially 2d ion source functions
!
subroutine fusion
  !
  use nbi_com
  !
  !-------------------------------------------------------------
  !  modified -- jfs 11 feb 94
  !
  !  added protons
  !  moved number and power summations into separate blocks
  !
  !-------------------------------------------------------------
  !  initial version -- dmc sept 1990
  !
  !  mod dmc 7 jul 1992 -- implement source renormalization option
  !  -- magnitude of fusion ion sources may be renormalized by
  !  -- input data
  !
  !  the source is added up from thermonuclear, beam-target, and
  !  beam-beam reactions
  !
  !-------------------------------------------------------------
  !
  !============
  implicit none
  integer :: ireact,iz,ibs
  !============
  real(fp) :: zfac
  !============
  if(.not.nlfusn) return
  !
  do ireact=1,nreact
     !
     !  the reaction index comes originally from hht's nuclear fusion
     !  cross section routine sigma.for:
     !	  =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
     !            (reaction 6 not supported in transp dmc sept 1990)
     !
     !  with modification:
     !         =7:  reversed d-t reaction
     !         =8:  reversed d-he3 reaction
     !
     !  and for "burnup" simulations:
     !         =9:  reversed d-t reaction -- t is a fusion product
     !         =10: reversed d-he3 reaction -- he3 is a fusion product
     !
     !  and for rf tail ion fusion simulations:
     !         =11:  reversed d-t reaction -- t is rf tail ion, d is **
     !         =12:  reversed d-he3 reaction -- he3 is rf tail ion, d is **
     !         =13:  forward d-t reaction -- d is rf tail ion, t is **
     !         =14:  forward d-he3 reaction -- d is rf tail ion he3 is **
     !         =15:  forward d+(d,p) reaction
     !         =16:  forward d+(d,n) reaction
     !
     !  ** beam or thermal
     !
     !  note (dmc oct 22 1997):  the following fast ion -- fast ion fusion
     !          rates are not computed:  (fusion product - fusion product)
     !                                   (rf tail - fusion product)
     !                                   (rf tail - rf tail)
     !
     !
     !  thermonuclear contribution
     !
     call map12(thntns(1,ireact),sfusion(1,ireact),0)
     !
     sftot(ireact)=0.0_fp
     !
     do iz=1,nfbzns
        !
        if(ireact.eq.5) then
           zfac=0.5_fp  ! only one ion per two neutrons for t-t
           sfusion(iz,ireact)=0.5_fp*sfusion(iz,ireact)
        else
           zfac=1.0_fp  ! one ion per neutron
        end if
        !
        !  beam-target contribution
        !
        sfusion(iz,ireact)=sfusion(iz,ireact)+zfac*btntn2(iz,ireact)

        !
        !  beam-beam contribution
        !
        sfusion(iz,ireact)=sfusion(iz,ireact)+zfac*bbnvol(iz,ireact)
        !
        !  add to total
        !
        sftot(ireact)=sftot(ireact)+sfusion(iz,ireact)
        !

     end do
     !

  end do
  !
  !  renormalize source rates if enabled
  !
  if(nlfst) then
     if(nlusft) call fusnrm(nftrit,xdatsft)
  end if
  !
  if(nlfsp) then
     if(nlusfp) call fusnrm(nfprot,xdatsfp)
  end if
  !
  if(nlfsp14) then
     if(nlusfp14) call fusnrm(nfprot14,xdatsp3)
  end if
  !
  if(nlfhe3) then
     if(nlusf3) call fusnrm(nfhel3,xdatsf3)
  end if
  !
  if(nlfhe4) then
     if(nlusfa) call fusnrm(nfhel4,xdatsfa)
  end if
  !
  !  add up summary outputs
  !   for fusion product species to be followed in the monte carlo code
  !
  do ibs=1,mibs
     do ireact=1,mireact
        nlreact(ireact,ibs)=.false.
        ereact(ireact,ibs)=0.0_fp
     end do
  end do
  !
  sftota=0.0_fp
  pftota=0.0_fp
  !
  if(nlfst) then
     call fuspwr(nftrit)
     sftota=sftota+xninjs(nftrit)
     pftota=pftota+pinjs(nftrit)
     call fuslis(nftrit,nlreact(1,nftrit),ereact(1,nftrit))
  end if
  !
  if(nlfsp) then
     call fuspwr(nfprot)
     sftota=sftota+xninjs(nfprot)
     pftota=pftota+pinjs(nfprot)
     call fuslis(nfprot,nlreact(1,nfprot),ereact(1,nfprot))
  end if
  !
  if(nlfsp14) then
     call fuspwr(nfprot14)
     sftota=sftota+xninjs(nfprot14)
     pftota=pftota+pinjs(nfprot14)
     call fuslis(nfprot14,nlreact(1,nfprot14),ereact(1,nfprot14))
  end if
  !
  if(nlfhe3) then
     call fuspwr(nfhel3)
     sftota=sftota+xninjs(nfhel3)
     pftota=pftota+pinjs(nfhel3)
     call fuslis(nfhel3,nlreact(1,nfhel3),ereact(1,nfhel3))
  end if
  !
  if(nlfhe4) then
     call fuspwr(nfhel4)
     sftota=sftota+xninjs(nfhel4)
     pftota=pftota+pinjs(nfhel4)
     call fuslis(nfhel4,nlreact(1,nfhel4),ereact(1,nfhel4))
  end if
  !
  return
end subroutine fusion
