!-----------------------------------------------------------------------
! btfusn -- evaluate beam target fusion reactions
!
!   dmc oct 1991 -- this routine is called from collid (collid.fnb)
!     this code extracted from collid and made a separate subroutine
!     2 oct 1991
!
!  this has been generalized from the old beam-target neutron
!  calculation -- now more fusion reactions are supported, even
!  those that do not produce neutrons.  the convention for output
!  is:  if the reaction produces neutrons, output the contribution
!  to the neutron production rate (i.e. include a factor of 2 for
!  the t+t-->2n,he4 reaction); but if the reaction does not
!  produce neutrons just output the contribution to the total
!  no. of reactions
!
!  note -- some reactions are not supported:
!     t+he3
!     he3+he3
!     t+t where one of the t's is a fusion product; for fusion
!      product burnup, burnup on background d *only* is computed
!      in this version of the routine (dmc 3 oct 1991)
!
!  this routine also does fusion product - target burnup reactions
!  e.g. burnup of fusion tritons from the d+d-->p+t reaction, if
!  we are following these fusion ions in the monte carlo code
!    the fast ion index to fusion product tritons is "nftrit";
!  nftrit is 0 if we are not following these fusion product ions.
!
!-----------------------------------------------------------------------
!
!
subroutine btfusn(ze1,zdwp)
  !
  !  ze1 -- beam ion energy (ev) **in the plasma frame**
  !  ztix -- target temperature (ev)
  !
  !  zdwp -- weight factor, from s.r. collid
  !
  !-----------------------------------------------------------------------
  !
  use nbatom_mod
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: irange,ie1,ie2
  !============
  real(fp) :: ztix,zdwp,ze1,zintrp,zsigv,zsigvd,zdneut
  !-----------------------------------------------------------------------
  !
  !
  if(nlfprod(lsbeam)) then
     irange=3
  else
     irange=2
  end if
  !
  ztix=tmjsm(nflr)
  !
  call btfusn_intrp(ze1,irange,ie1,ie2,zintrp)
  !
  !  calculate beam-target neutron production.
  !
  if(lsbeam.eq.nbdeut) then
     !
     !  beam ion is deuterium
     !
     !	d. mc cune 20 oct 1982 -- new improved <sigma*v> d-d neutrons
     !	ze1 is energy of beam ion in ev, plasma frame
     !
     !
     !  dt beam target fusion
     !
     zsigvd=0.0_fp
     if(ntrit.gt.0) then
        zsigv=(1.0_fp-zintrp)*btfus_dt(nflr,ie1)+ &
             zintrp*btfus_dt(nflr,ie2)
        zsigvd=rhob(nflr,ntrit,2)*zsigv
        zdneut=zsigvd*zdwp
        btntns(nflr,1)=btntns(nflr,1)+zdneut
        btntn2(nflr2,1)=btntn2(nflr2,1)+zdneut
        btntrs(nflr,1)=btntrs(nflr,1)+zsigv*zdwp
        btntws(nflr,1)=btntws(nflr,1)+zdwp
     end if
     !
     !  d beam rf min t target fusion
     !
     if(nrftrit.gt.0) then
        if(nthrf.ne.1)  then
           call fusn_btsigv_dt(ze1*0.001_fp, &
                timin(nflr,2,mrftrit)*0.001_fp,zsigv,1)
           zsigvd=rhmin(nflr,2,mrftrit)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           btntns(nflr,11)=btntns(nflr,11)+zdneut
           btntn2(nflr2,11)=btntn2(nflr2,11)+zdneut
           btntrs(nflr,11)=btntrs(nflr,11)+zsigv*zdwp
           btntws(nflr,11)=btntws(nflr,11)+zdwp
        else
           zsigvd=rhmin(nflr,2,mrftrit)*zsigv
           zdneut=zsigvd*zdwp
           btntns(nflr,1)=btntns(nflr,1)+zdneut
           btntn2(nflr2,1)=btntn2(nflr2,1)+zdneut
           btntrs(nflr,1)=btntrs(nflr,1)+zsigv*zdwp
           btntws(nflr,1)=btntws(nflr,1)+zdwp
        end if
     end if
     !
     !  d-he3 beam target fusion
     !
     if(nhel3.gt.0) then
        zsigv=(1.0_fp-zintrp)*btfus_d3(nflr,ie1)+ &
             zintrp*btfus_d3(nflr,ie2)
        zsigvd=rhob(nflr,nhel3,2)*zsigv
        zdneut=zsigvd*zdwp
        btntns(nflr,2)=btntns(nflr,2)+zdneut
        btntn2(nflr2,2)=btntn2(nflr2,2)+zdneut
        btntrs(nflr,2)=btntrs(nflr,2)+zsigv*zdwp
        btntws(nflr,2)=btntws(nflr,2)+zdwp
     end if
     !
     !  d beam rf min he3 target fusion
     !
     if(nrfhel3.gt.0) then
        if(nthrf.ne.1) then
           call fusn_btsigv_d3(ze1*0.001_fp, &
                timin(nflr,2,mrfhel3)*0.001_fp,zsigv,1)
           zsigvd=rhmin(nflr,2,mrfhel3)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           btntns(nflr,12)=btntns(nflr,12)+zdneut
           btntn2(nflr2,12)=btntn2(nflr2,12)+zdneut
           btntrs(nflr,12)=btntrs(nflr,12)+zsigv*zdwp
           btntws(nflr,12)=btntws(nflr,12)+zdwp
        else
           zsigvd=rhmin(nflr,2,mrfhel3)*zsigv
           zdneut=zsigvd*zdwp
           btntns(nflr,2)=btntns(nflr,2)+zdneut
           btntn2(nflr2,2)=btntn2(nflr2,2)+zdneut
           btntrs(nflr,2)=btntrs(nflr,2)+zsigv*zdwp
           btntws(nflr,2)=btntws(nflr,2)+zdwp
        end if
     end if

     !
     !  dd beam target fusion
     !
     if(ndeut.gt.0) then
        call fusn_btsigv_ddn(ze1*0.001_fp,ztix*0.001_fp,zsigv,1)
        zsigvd=rhob(nflr,ndeut,2)*zsigv*1.0e6_fp
        zdneut=zsigvd*zdwp
        btntns(nflr,4)=btntns(nflr,4)+zdneut
        btntn2(nflr2,4)=btntn2(nflr2,4)+zdneut
        btntrs(nflr,4)=btntrs(nflr,4)+zsigv*zdwp
        btntws(nflr,4)=btntws(nflr,4)+zdwp
        !     dd(p) fusion
        if(nlddp) then
           call fusn_btsigv_ddp(ze1*0.001_fp,ztix*0.001_fp,zsigv,1)
           zsigvd=rhob(nflr,ndeut,2)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           btntns(nflr,3)=btntns(nflr,3)+zdneut
           btntn2(nflr2,3)=btntn2(nflr2,3)+zdneut
           btntrs(nflr,3)=btntrs(nflr,3)+zsigv*zdwp
           btntws(nflr,3)=btntws(nflr,3)+zdwp
        end if
     end if

     !  d beam rf min d target fusion
     if(nrfdeut.gt.0) then
        if(nthrf.ne.1) then
           call fusn_btsigv_ddn(ze1*0.001_fp, &
                timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
           zsigvd=rhmin(nflr,2,mrfdeut)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           btntns(nflr,16)=btntns(nflr,16)+zdneut
           btntn2(nflr2,16)=btntn2(nflr2,16)+zdneut
           btntrs(nflr,16)=btntrs(nflr,16)+zsigv*zdwp
           btntws(nflr,16)=btntws(nflr,16)+zdwp
        else
           zsigvd=rhmin(nflr,2,mrfdeut)*zsigv
           zdneut=zsigvd*zdwp
           btntns(nflr,4)=btntns(nflr,4)+zdneut
           btntn2(nflr2,4)=btntn2(nflr2,4)+zdneut
           btntrs(nflr,4)=btntrs(nflr,4)+zsigv*zdwp
           btntws(nflr,4)=btntws(nflr,4)+zdwp
        end if

        if(nlddp) then
           !     dd(p) fusion
           call fusn_btsigv_ddp(ze1*0.001_fp, &
                timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
           zsigvd=rhmin(nflr,2,mrfdeut)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           !
           btntns(nflr,15)=btntns(nflr,15)+zdneut
           btntn2(nflr2,15)=btntn2(nflr2,15)+zdneut
           btntrs(nflr,15)=btntrs(nflr,15)+zsigv*zdwp
           btntws(nflr,15)=btntws(nflr,15)+zdwp
        end if
     end if
     !
     !-----------------------------------------------------------------------
     !
  else if(lsbeam.eq.nbtrit) then
     !
     !  beam ion is tritium
     !
     !
     !	d. mc cune 20 oct 1982 -- new improved <sigma*v> d-d neutrons
     !	ze1 is energy of beam ion in ev, plasma frame
     !       td beam target fusion
     !
     zsigvd=0.0_fp
     if(ndeut.gt.0) then
        zsigv=(1.0_fp-zintrp)*btfus_td(nflr,ie1)+ &
             zintrp*btfus_td(nflr,ie2)
        zsigvd=rhob(nflr,ndeut,2)*zsigv
        zdneut=zsigvd*zdwp
        btntns(nflr,7)=btntns(nflr,7)+zdneut
        btntn2(nflr2,7)=btntn2(nflr2,7)+zdneut
        btntrs(nflr,7)=btntrs(nflr,7)+zsigv*zdwp
        btntws(nflr,7)=btntws(nflr,7)+zdwp
     end if
     if(nrfdeut.gt.0) then
        if(nthrf.ne.1) then
           call fusn_btsigv_td(ze1*0.001_fp, &
                timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
           zsigvd=rhmin(nflr,2,mrfdeut)*zsigv*1.0e6_fp
           zdneut=zsigvd*zdwp
           btntns(nflr,13)=btntns(nflr,13)+zdneut
           btntn2(nflr2,13)=btntn2(nflr2,13)+zdneut
           btntrs(nflr,13)=btntrs(nflr,13)+zsigv*zdwp
           btntws(nflr,13)=btntws(nflr,13)+zdwp
        else
           zsigvd=rhmin(nflr,2,mrfdeut)*zsigv
           btntns(nflr,7)=btntns(nflr,7)+zdneut
           btntn2(nflr2,7)=btntn2(nflr2,7)+zdneut
           btntrs(nflr,7)=btntrs(nflr,7)+zsigv*zdwp
           btntws(nflr,7)=btntws(nflr,7)+zdwp
        end if
     end if
     !
     !  tt beam target fusion
     zsigvd=0.0_fp
     if(ntrit.gt.0.or.nrftrit.gt.0) then
        zsigv=(1.0_fp-zintrp)*btfus_tt(nflr,ie1)+ &
             zintrp*btfus_tt(nflr,ie2)
        if(ntrit.gt.0) then
           zsigvd=rhob(nflr,ntrit,2)*zsigv
        end if
        if(nrftrit.gt.0) then
           if(nthrf.ne.1) then
              call fusn_btsigv_tt(ze1*0.001_fp, &
                   timin(nflr,2,mrftrit)*0.001_fp,zsigv,1)
              zsigv=zsigv*1.0e6_fp
           end if
           zsigvd=zsigvd+rhmin(nflr,2,mrftrit)*zsigv
        end if
        !     get 2 neutrons per reaction
        zdneut=zsigvd*zdwp*2.0_fp
        btntns(nflr,5)=btntns(nflr,5)+zdneut
        btntn2(nflr2,5)=btntn2(nflr2,5)+zdneut
        btntrs(nflr,5)=btntrs(nflr,5)+zsigv*zdwp
        btntws(nflr,5)=btntws(nflr,5)+zdwp
     end if
     !
  else if(lsbeam.eq.nbhel3) then
     !
     !  beam ion is helium-3
     !
     !
     !  he3-d beam target fusion
     zsigvd=0.0_fp
     if(ndeut.gt.0.or.nrfdeut.gt.0) then
        zsigv=(1.0_fp-zintrp)*btfus_3d(nflr,ie1)+ &
             zintrp*btfus_3d(nflr,ie2)
        if(ndeut.gt.0) then
           zsigvd=rhob(nflr,ndeut,2)*zsigv
           zdneut=zsigvd*zdwp
           btntns(nflr,8)=btntns(nflr,8)+zdneut
           btntn2(nflr2,8)=btntn2(nflr2,8)+zdneut
           btntrs(nflr,8)=btntrs(nflr,8)+zsigv*zdwp
           btntws(nflr,8)=btntws(nflr,8)+zdwp
        end if
        if(nrfdeut.gt.0) then
           if(nthrf.eq.1) then
              zsigvd=rhmin(nflr,2,mrfdeut)*zsigv
              zdneut=zsigvd*zdwp
              btntns(nflr,8)=btntns(nflr,8)+zdneut
              btntn2(nflr2,8)=btntn2(nflr2,8)+zdneut
              btntrs(nflr,8)=btntrs(nflr,8)+zsigv*zdwp
              btntws(nflr,8)=btntws(nflr,8)+zdwp
           else
              call fusn_btsigv_3d(ze1*0.001_fp, &
                   timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
              zsigvd=rhmin(nflr,2,mrfdeut)*zsigv*1.0e6_fp
              zdneut=zsigvd*zdwp
              btntns(nflr,14)=btntns(nflr,14)+zdneut
              btntn2(nflr2,14)=btntn2(nflr2,14)+zdneut
              btntrs(nflr,14)=btntrs(nflr,14)+zsigv*zdwp
              btntws(nflr,14)=btntws(nflr,14)+zdwp
           end if
        end if
     end if
     !
  else if(lsbeam.eq.nftrit) then
     !
     !  **fusion product ion**  burnup calculation  t + d
     !
     zsigvd=0.0_fp
     if(ndeut.gt.0.or.(nrfdeut.gt.0.and.nthrf.eq.1)) then
        zsigv=(1.0_fp-zintrp)*btfus_td(nflr,ie1)+ &
             zintrp*btfus_td(nflr,ie2)
        if(ndeut.gt.0) then
           zsigvd=rhob(nflr,ndeut,2)*zsigv
        end if
     end if
     if(nrfdeut.gt.0) then
        if(nthrf.ne.1) then
           call fusn_btsigv_td(ze1*0.001_fp, &
                timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
           zsigv=zsigv*1.0e6_fp
        end if
        zsigvd=zsigvd+rhmin(nflr,2,mrfdeut)*zsigv
     end if
     zdneut=zsigvd*zdwp
     btntns(nflr,9)=btntns(nflr,9)+zdneut
     btntn2(nflr2,9)=btntn2(nflr2,9)+zdneut
     btntrs(nflr,9)=btntrs(nflr,9)+zsigv*zdwp
     btntws(nflr,9)=btntws(nflr,9)+zdwp
     !
  else if(lsbeam.eq.nfhel3) then
     !
     !  **fusion product ion**  burnup calculation  he3 + d
     !
     zsigvd=0.0_fp
     if(ndeut.gt.0.or.(nrfdeut.gt.0.and.nthrf.eq.1)) then
        zsigv=(1.0_fp-zintrp)*btfus_3d(nflr,ie1)+ &
             zintrp*btfus_3d(nflr,ie2)
        if(ndeut.gt.0) then
           zsigvd=rhob(nflr,ndeut,2)*zsigv
        end if
     end if
     if(nrfdeut.gt.0) then
        if(nthrf.ne.1) then
           call fusn_btsigv_3d(ze1*0.001_fp, &
                timin(nflr,2,mrfdeut)*0.001_fp,zsigv,1)
           zsigv=1.0e6_fp*zsigv
        end if
        zsigvd= zsigvd+rhmin(nflr,2,mrfdeut)*zsigv
     end if
     zdneut=zsigvd*zdwp
     btntns(nflr,10)=btntns(nflr,10)+zdneut
     btntn2(nflr2,10)=btntn2(nflr2,10)+zdneut
     btntrs(nflr,10)=btntrs(nflr,10)+zsigv*zdwp
     btntws(nflr,10)=btntws(nflr,10)+zdwp
  end if
  !
  !
10000 return
  !---------------------------------------
  !  internal subroutine for table interpolation
  !  (see ebatom & nbatom subroutines for table setups)
  !
contains
  subroutine btfusn_intrp(ze,irng,ie,iep,zifac)
    use iso_c_binding, only: fp => c_double


    real(fp), intent(in) :: ze         ! e for table
    integer, intent(inout) :: irng ! table energy range sector
    integer, intent(out) :: ie,iep ! bracketing energy indices
    real(fp), intent(out) :: zifac     ! linear interpolation factor

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

    integer :: isbe,inume,inumr,imin,imax,inum,inc
    real(fp) :: zea,zemin,zemax

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

    isbe=0                         ! always "0" for btfusn
    inumr=nbnsver(isbe)            ! no. of energy range sectors
    inume=nbnsve(isbe)             ! tot. no. energy grid pts

    zea=ze
    if(zea.gt.bnsves(inume,isbe)) then
       if(orbrzv_rf_option.ne.1.and. &
            (ntest_max.eq.0.or.abs(ntest_max).eq.10)) then ! mg: warning disabled when rf-option is on
          write(lun_nbatom_mod,9001) zea,bnsves(inume,isbe)
9001      format( &
               ' ?btfusn_intrp:  interaction energy/a exceeds'/ &
               '  table limit:  energy/a = ',1pe11.4,' limit = ', &
               1pe10.4,' ev/nucleon')
       end if
       zea = bnsves(inume,isbe)
    end if
    !
    !  guard low energy limit
    !
    zea=max(bnsves(1,isbe),zea)
    !
    !  guard sector range
    !
    irng=max(1,min(inumr,irng))
    !
    do
       if(irng.eq.1) then
          imin=1
       else
          imin=lbnsve(irng-1,isbe)
       end if
       imax=lbnsve(irng,isbe)
       !
       zemin=bnsves(imin,isbe)
       zemax=bnsves(imax,isbe)
       !
       if(zea.lt.zemin) then
          irng=irng-1
          cycle
          !
       else if(zea.gt.zemax) then
          irng=irng+1
          cycle
          !
       else
          !  in range...
          inum=imax-imin
          inc=min((inum-1),int(inum*(zea-zemin)/(zemax-zemin)))
          !
          ie=imin+inc
          iep=ie+1
          !
          zifac=(zea-bnsves(ie,isbe))/ &
               (bnsves(iep,isbe)-bnsves(ie,isbe))
          !
          exit
          !
       end if
    end do
    !
  end subroutine btfusn_intrp
end subroutine btfusn

!------------------------------------------------------------------------
!  dmc: removing lsbeamr calculation from subroutine btfusn

subroutine btfusn_init
  ! dmc -- lsbeamr calculation now done from nbstart, not inside particle
  ! loop.  this avoids problem with mpi reduce (inconsistent logic btw
  ! processors, as nuse -> 0 so that on some processors btfusn was never
  ! called.

  use nbi_com

  !---------------------------------------------
  integer :: i_lsbeam
  !---------------------------------------------
  !
  !
  do i_lsbeam=1,nsbeam
     if(nlfbmfpp(i_lsbeam)) cycle
     !
     !  calculate beam-target neutron production.
     !
     if(i_lsbeam.eq.nbdeut) then
        !
        !  beam ion is deuterium
        !
        !  dt beam target fusion
        if(ntrit.gt.0) then
           lsbeamr(1)=i_lsbeam
        end if

        !  d beam rf min t target fusion
        if(nrftrit.gt.0) then
           lsbeamr(11)=i_lsbeam
        end if
        !
        !  d-he3 beam target fusion
        if(nhel3.gt.0) then
           lsbeamr(2)=i_lsbeam
        end if

        !  d beam rf min he3 target fusion
        if(nrfhel3.gt.0) then
           lsbeamr(12)=i_lsbeam
        end if
        !
        !  dd beam target fusion
        !  dd(n) fusion; d thermal target always there if d beam ions exist
        lsbeamr(4)=i_lsbeam

        !  d beam rf min d target fusion
        if(nrfdeut.gt.0) then
           lsbeamr(16)=i_lsbeam
        end if

        if(nlddp) then
           !  dd(p) fusion
           lsbeamr(3)=i_lsbeam
           if(nrfdeut.gt.0) then
              lsbeamr(15)=i_lsbeam
           end if
        end if
        !
        !-----------------------------------------------------------------------
        !
     else if(i_lsbeam.eq.nbtrit) then
        !
        !  beam ion is tritium
        !
        if(ndeut.gt.0) then
           lsbeamr(7)=i_lsbeam
        end if

        if(nrfdeut.gt.0) then
           lsbeamr(13)=i_lsbeam
        end if
        !
        !  tt beam target fusion  --  t thermal ions always there
        !
        lsbeamr(5)=i_lsbeam
        !
     else if(i_lsbeam.eq.nbhel3) then
        !
        !  beam ion is helium-3
        !
        !  he3->d beam target
        if(ndeut.gt.0) then
           lsbeamr(8)=i_lsbeam
        end if

        if(nrfdeut.gt.0) then
           lsbeamr(14)=i_lsbeam
        end if

        !
     else if(i_lsbeam.eq.nftrit) then
        !
        !  **fusion product ion**  burnup calculation  t + d
        !
        if((ndeut.gt.0).or.(nrfdeut.gt.0)) then
           lsbeamr(9)=i_lsbeam
        end if
        !
     else if(i_lsbeam.eq.nfhel3) then
        !
        !  **fusion product ion**  burnup calculation  he3 + d
        !
        if((ndeut.gt.0).or.(nrfdeut.gt.0)) then
           lsbeamr(10)=i_lsbeam
        end if
        !
     end if
     !
  end do

end subroutine btfusn_init
