!******************** start file cxnsum.fnb ; group nubeam_lib  ********
!----------------------------------------------------------------
!  cxnsum
!
!  calculate ztscx =
!  sum n0*sigma*v for charge-exchange of beam ion with all
!  types of neutral densities present.  if beam ion is in "beam
!  region" include 1st generation beam neutrals and volume source
!  "halo" neutrals; include nth generation beam neutrals and wall
!  source neutrals in any case.
!
!
subroutine cxnsum(ztscx,isflag,zwght)
  !
  use nbatom_mod
  use physconst_mod, only : evptocm_sec
  use logmod
  !
  !----------------------------------------------------------------
  !
  !  -------> dmc 18 june 2001 ... convert to use ntcc `preact' module
  !           cf http://w3.pppl.gov/ntcc, catalog, `preact'...
  !           replacing presig, etc., for atomic physics data.
  !
  !  mod tbt jan 1992 - commented out ztemp debug statements.
  !  mod tbt dec 1991 - added ichrg=2 to call of presig to indicate that
  !                     the baem is charged and the background is neutral.
  !  mod tbt oct 1991 - change code to reflect nlvphi=t assumption.
  !
  !  mod tbt oct 1991 - do away with sigvbt with calls to presig.
  !                     this enables use of helium.
  !                     put in loop to do away with dn0out approximation.
  !
  !  mod dmc feb 1987 - zeovab computed locally
  !    beam-beam interaction velocities modified - both neutral and
  !  ion velocity vectors are evaluated in flux surface tang-normal
  !  frame instead of (r,y,z) frame
  !
  !  if flag isflag is set perform summing on monte-carlo arrays
  !  geometry generalization dmc apr 1985 -
  !   full 2d neutral density arrays for beam and halo neutrals assumed
  !  (old afac kluge deleted)
  !
  !  beam ion - 1st gen beam neutral interaction rewritten to use the
  !   neutral velocity vector output from dep
  !
  !----------------------------------------------------------------
  !
  !  zwght-- wght to be used in adding to cx thermal neutral sink
  !   profile or external cx beam ion loss rate-- these sums done only
  !   if isflag=1
  !
  !  for particles outside the plasma use only dn0out
  !  rjg 7/19/83

  use nbi_com
  !
  !%      use sglib_mod
  !
  !  ion velocity vector
  !============
  implicit none
  integer :: isflag,icx,ii,izb,irange,irange_bb,ig,iztarg,ind
  integer :: indp1,isb,ien,ib,iznb,ino,idir
  !============
  real(fp) :: zwght,ztscx,zeovab,zsumro,zintrp,znsigv,zt0a,zvz0
  real(fp) :: zdtscx,zwdtsx,zwdtsi,zebr,zerel,zhexfac,zsigv,zn02,zdtsii
  real(fp) :: zrati2,zratx2,zn0,zv2,zv0dir,zsigv_cx,zsvdtg,zdsncx
  real(fp) :: zsigv_ii,zdsnii
  real(fp) :: ze_neut,zksid_0
  integer :: lgyro_neut,ini,jba_0
  !============
  real(fp) :: zvi(3),a_rds
  !
  !%      real(fp) :: zzx(2),zzy(2),zr,zy,zdummy(2,2)
  !
  integer :: ichrg
  !
  integer :: i_fn_count,i_co_count, i_cnt_count,i_dir
  !
  integer :: jba, jbe !  pitch angle, energy indices
  real(fp) :: zep !ptcl. energy (ev)
  integer :: ndxfbe_r8 !integer function
  logical :: il_3dhalo
  !
  !-----------------------------------------------------------

  il_3dhalo=(levmod_halo.eq.2)
  ichrg = 2      ! = 2 indicates charged beam to hatom_btsigv
  icx = 1        ! indicates charge exchange to hatom_btsigv
  ii = 2         ! indicates impact ionization to hatom_btsigv

  izb = xzbeami + 0.1_fp

  !**************  sum dtgcx/tau-cx **************************
  !      taking into account all possible collisions
  !
  !  zero sum counter
  !
  !%      zzx=0; zzy=0
  ztscx=0.0_fp
  !
  zeovab=0.0_fp ! must be re-evaluated for each frame
  !
  !  deal with particles in the neutral halo outside the plasma
  !
  if(nlfprod(lsbeam)) then
     irange=3
  else
     irange=1
  end if
  irange_bb=2
  !
  if (nlpout) then
     zeovab=ep/evptocm_sec/evptocm_sec    !*5.2207e-13_fp

     !         ................................................ tbt 10/91
     zsumro = 0.0_fp
     do ig=1,ng
        zsumro = zsumro + rhob(ledge,ig,1) ! use for normalization.

     end do

     do ig=1,ng
        iztarg=backz(ig)+0.1_fp
        !  assume 1 ev cloud neutral temperature
        !dbg            call hatom_btsigv(icx,ichrg,zeovab*0.001,0.001,1,
        !dbg     >         izb,iztarg,zsigv,ier)
        !dbg            zsigv=1.0e6*zsigv
        !dbg            znsigv=dn0out*zsigv
        call cxnsum_intrp(zeovab,lsbeam,irange,ind,indp1,zintrp)
        znsigv=(1.0_fp-zintrp)*cxn_thcx_a(ig,ind,lsbeam) + &
             zintrp*cxn_thcx_a(ig,indp1,lsbeam)
        ztscx=ztscx + dtgcx * znsigv * (rhob(ledge,ig,1) / zsumro)

     end do
     !         ................................................

     !  sum for ptcl balance
     if(isflag.eq.1) then
        snbcxos(lsbeam)=snbcxos(lsbeam)+zwght*ztscx

        !%            if(nbplot.gt.0) then
        !%            zzx(1)=1.0
        !%            zzy(1)=1.0
        !%              if(nbplot.eq.1) then
        !%            call xmomry(xiflr,thflr,0,zr,zy,zdummy)
        !%            zzx(2)=zr
        !%            zzy(2)=zy
        !%              else
        !%            zzx(2)=vplflr
        !%            zzy(2)=vppflr
        !%              end if
        !%            call symbl(3)
        !%            call line(-1)
        !%            call r8_cplot(zzx,zzy)
        !%            end if

     end if
     goto 250      ! return
  end if
  !
  !---------------------------------------------------------------
  !  ion is in the plasma
  !
  !----------------------
  !
  !  ion velocity vector
  !   btiltc factors-- see bparms
  !
  !  toroidal direction:
  zvi(1)=btiltc(2,nflr2)*vplflr+btiltc(3,nflr2)*vppflr*cospho
  !  flux surface normal direction:
  zvi(2)=-vppflr*sinpho
  !  flux suface tangential direction in poloidal plane:
  zvi(3)=nsnccw*(vppflr*btiltc(2,nflr2)*cospho &
       -btiltc(3,nflr2)*vplflr)
  !
  !----------------------
  !
  !  wall (edge) source neutrals
  !
  do ig=1,ng
     iztarg=backz(ig)+0.1_fp
     if(walln0(ig,nflr).eq.0.0_fp) cycle
     zt0a=ewall0(ig,nflr)/aplasm(ig)
     zvz0=owall0(ig,nflr)*rmjflr
     zeovab=((zvi(1)-zvz0)**2+zvi(2)**2+zvi(3)**2)/evptocm_sec/evptocm_sec
     call cxnsum_intrp(zeovab,lsbeam,irange,ind,indp1,zintrp)
     !
     !      **d. mc cune 19 feb 1981** using 2d sigma*v lookup routines
     !  charge-exchange:
     !dbg        call hatom_btsigv(icx,ichrg,zeovab*0.001,zt0a*0.001,1,
     !dbg     >     izb,iztarg,zsigv,ier)
     !dbg        zsigv=1.0e6*zsigv
     !dbg        znsigv=walln0(ig,nflr)*zsigv
     znsigv=(1.0_fp-zintrp)*cxn_thcx_wa(ig,nflr,ind,lsbeam) + &
          zintrp*cxn_thcx_wa(ig,nflr,indp1,lsbeam)
     zdtscx=dtgcx*znsigv
     ztscx=ztscx+zdtscx
     if(isflag.ne.1) cycle
     !  charge-exchange
     zwdtsx=zwght*zdtscx
     snbxw0(lsbeam)=snbxw0(lsbeam)+zwdtsx
     sinbn0s(nflr,ig,lsbeam)=sinbn0s(nflr,ig,lsbeam)+zwdtsx
     snbn02(nflr2,ig,lsbeam)=snbn02(nflr2,ig,lsbeam)+zwdtsx
     !  impact-ionization:
     !dbg          call hatom_btsigv(ii,ichrg,zeovab*0.001,zt0a*0.001,1,
     !dbg     >       izb,iztarg,zsigv,ier)
     !dbg          zsigv=1.0e6*zsigv
     !dbg          znsigv=walln0(ig,nflr)*zsigv
     znsigv=(1.0_fp-zintrp)*cxn_thii_wa(ig,nflr,ind,lsbeam) + &
          zintrp*cxn_thii_wa(ig,nflr,indp1,lsbeam)

     zwdtsi=zwght*dtgcx*znsigv
     sinb0is(nflr,ig,lsbeam)=sinb0is(nflr,ig,lsbeam)+zwdtsi
     snb0i2(nflr2,ig,lsbeam)=snb0i2(nflr2,ig,lsbeam)+zwdtsi
     !
  end do
  !
  !  nth generation beam neutrals
  !
  do isb=1,nsbeam
     iznb=xzbeams(isb)
     iztarg = xzbeams(isb)+0.1_fp
     if(nlbbcx.and.(bn0x2p(nflr2,isb).gt.0.0_fp)) then
        zebr=((zvi(1)-bv0x2p(nflr2,isb))**2 &
             +(zvi(2)-bvrhox2p(nflr2,isb))**2 &
             +(zvi(3)-bvthx2p(nflr2,isb))**2)*v2toev
        !     mg            zerel=(zebr+be0x2p(nflr2,isb))/abeam !abeam here due to abeam in v2toev(see inispc.for)
        zerel=zebr/abeam    !abeam here due to abeam in v2toev(see inispc.for)
        if(nsigexc_use.gt.0) then
           zhexfac=max(one,hexbbxrp(nflr2,isb,lsbeam))
        else
           zhexfac=1.0_fp
        end if
        !dbg            call hatom_sigv(icx,0.001*zerel,1,iztarg,izb,zsigv,ier)
        !dbg            zsigv=1.0e6*zsigv
        call cxnsum_intrp(zerel,0,irange_bb,ind,indp1,zintrp)
        zsigv_cx=(1.0_fp-zintrp)*cxn_bbcx(iznb,ind,izb) + &
             zintrp*cxn_bbcx(iznb,indp1,izb)
        zsigv_ii=(1.0_fp-zintrp)*cxn_bbii(iznb,ind,izb) + &
             zintrp*cxn_bbii(iznb,indp1,izb)
        !     mg i think it is not consistent to use 2d look up table as a function of
        !     mg  z_rel and 'temperature' of neutrals
        !     mg 03/06/12            zsigv=(1.0_fp-zintrp)*cxn_bbcx_2d(isb,nflr2,ind,lsbeam) +
        !     mg 03/06/12      >         zintrp*cxn_bbcx_2d(isb,nflr2,indp1,lsbeam)
        !     mg
        zdtscx=dtgcx*zsigv_cx*zhexfac
        zdtsii=dtgcx*zsigv_ii*zhexfac
        ztscx=ztscx+zdtscx
        if(isflag.eq.1) then
           snbxbb1(lsbeam)=snbxbb1(lsbeam)+ &
                zdtscx*zwght*bn0x2p(nflr2,isb)/bmvol(nflr2)
           !     collect cx and ii nb * sig * v
           if((nflr2.gt.0).and.(nflr2.le.mimxbzf).and. &
                bn0x2p(nflr2,isb).gt.0.0_fp) then !inside plasma only and non zero neutral density
              ze_neut=bv0x2p(nflr2,isb)**2+bvrhox2p(nflr2,isb)**2+ &
                   bvthx2p(nflr2,isb)**2 ! energy of neutral in lab frame
              ze_neut=ze_neut*v2toev/abeam !abeam here due to abeam in v2toev(see inispc.for)
              ! find index on coarse energy grid for neutrals

              call cxnsum_intrp(ze_neut,-1,irange_bb,ind,indp1, &
                   zintrp)
              zksid_0=bvpllx2p(nflr2,isb) !pitch angle for neutral
              !  pitch angle,
              jba_0=int((1.0_fp+zksid_0)*0.5_fp*real(nznbma0))+1
              jba_0=min(nznbma0,max(1,jba_0))

              ! find neutral's lgyro in there
              if(.not.nlfbmflr) then
                 lgyro_neut=1
              else
                 call nbi_fbm_grid_indx(balpha_gyrop(nflr2,isb), &
                      lgyro_neut)
              end if
              if(il_3dhalo) then
                 if(iztarg.eq.1) then ! hydroden isotopes
                    !     cx nb * sig*v for  hydroden isotopes
                    svn3vcx_h(ind,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vcx_h(ind,jba_0,lgyro_neut,nflr2,lsbeam)+ &
                         (1.0_fp-zintrp)*zdtscx*zwght/bmvol(nflr2) ! [1/sec]
                    svn3vcx_h(indp1,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vcx_h(indp1,jba_0,lgyro_neut,nflr2,lsbeam) &
                         +zintrp*zdtscx*zwght/bmvol(nflr2) ! [1/sec]
                    ! ii nb * sig*v for  hydroden isotopes
                    svn3vii_h(ind,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vii_h(ind,jba_0,lgyro_neut,nflr2,lsbeam)+ &
                         (1.0_fp-zintrp)*zdtsii*zwght/bmvol(nflr2) ! [1/sec]
                    svn3vii_h(indp1,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vii_h(indp1,jba_0,lgyro_neut,nflr2,lsbeam) &
                         +zintrp*zdtsii*zwght/bmvol(nflr2) ! [1/sec]
                 elseif(iztarg.eq.2) then ! he isotopes
                    ! cx nb * sig*v for  he isotopes
                    svn3vcx_he(ind,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vcx_he(ind,jba_0,lgyro_neut,nflr2,lsbeam)+ &
                         (1.0_fp-zintrp)*zdtscx*zwght/bmvol(nflr2) ! [1/sec]
                    svn3vcx_he(indp1,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vcx_he(indp1,jba_0,lgyro_neut,nflr2,lsbeam) &
                         +zintrp*zdtscx*zwght/bmvol(nflr2) ! [1/sec]
                    ! ii nb * sig*v for  he isotopes
                    svn3vii_he(ind,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vii_he(ind,jba_0,lgyro_neut,nflr2,lsbeam)+ &
                         (1.0_fp-zintrp)*zdtsii*zwght/bmvol(nflr2) ! [1/sec]
                    svn3vii_he(indp1,jba_0,lgyro_neut,nflr2,lsbeam)= &
                         svn3vii_he(indp1,jba_0,lgyro_neut,nflr2,lsbeam) &
                         +zintrp*zdtsii*zwght/bmvol(nflr2) ! [1/sec]
                 end if
              end if !end il_3dhalo
           end if
        end if
     end if

  end do
  !
  !
  !  thermal volume source neutrals (most of these come from
  !   beam neutral - thermal ion cx)
  !
  !  rewritten dmc mar 1985 ***
  !   volume sce neutrals on 2d grid; cf subroutine parms
  !
  do ig=1,ng
     iztarg=backz(ig)+0.1_fp
     if(voln02(ig,nflr2).eq.0.0_fp) cycle
     zt0a=evol02(ig,nflr2)/aplasm(ig)
     zvz0=ovol02(ig,nflr2)*rmjflr
     zeovab=((zvi(1)-zvz0)**2+zvi(2)**2+zvi(3)**2)/evptocm_sec/evptocm_sec
     call cxnsum_intrp(zeovab,lsbeam,irange,ind,indp1,zintrp)
     !
     !        **d. mc cune 19 feb 1981  using 2d lookup routines
     !        charge-exchange:
     !dbg        call hatom_btsigv(icx,ichrg,zeovab*0.001,zt0a*0.001,1,
     !dbg     >     izb,iztarg,zsigv,ier)
     !dbg        zsigv=1.0e6*zsigv
     !dbg        znsigv=voln02(ig,nflr2)*zsigv
     znsigv=(1.0_fp-zintrp)*cxn_thcx_ha(ig,nflr2,ind,lsbeam) + &
          zintrp*cxn_thcx_ha(ig,nflr2,indp1,lsbeam)

     zdtscx=dtgcx*znsigv
     ztscx=ztscx+zdtscx
     if(isflag.ne.1) cycle
     !  charge-exchange
     zwdtsx=zwght*zdtscx
     snbxv0(lsbeam)=snbxv0(lsbeam)+zwdtsx
     sinbn0s(nflr,ig,lsbeam)=sinbn0s(nflr,ig,lsbeam)+zwdtsx
     snbn02(nflr2,ig,lsbeam)=snbn02(nflr2,ig,lsbeam)+zwdtsx
     !  impact-ionization:
     !dbg          call hatom_btsigv(ii,ichrg,zeovab*0.001,zt0a*0.001,1,
     !dbg     >       izb,iztarg,zsigv,ier)
     !dbg          zsigv=1.0e6*zsigv
     !dbg          znsigv=voln02(ig,nflr2) * zsigv
     znsigv=(1.0_fp-zintrp)*cxn_thii_ha(ig,nflr2,ind,lsbeam) + &
          zintrp*cxn_thii_ha(ig,nflr2,indp1,lsbeam)

     zwdtsi=zwght*dtgcx*znsigv

     sinb0is(nflr,ig,lsbeam)=sinb0is(nflr,ig,lsbeam)+zwdtsi
     snb0i2(nflr2,ig,lsbeam)=snb0i2(nflr2,ig,lsbeam)+zwdtsi

  end do
  !
  !  interact with first generation beam neutrals
  !  skip if there are none
  !
  if(.not.nlbbcx) goto 250
  if(xninj.le.0.0_fp) goto 250
  !
  !  do beam-beam charge-exchange and impact ionization
  !
  !  inside/outside index - suppressed dmc mar 1985 - passing to
  !  full 2d grid
  !
  ! collect density for fast ptcls. separately for fusion product
  ! and for beam ptcls. depend on their direction of injection.
  ! i_dir=1 co-injected ptcls
  ! i_dir=2 counter-injected ptcls
  ! i_dir=0 fusion product
  !
  !
  if(nlfprod(lsbeam)) then
     i_dir=0
  else
     i_dir=2
     if(nlco(nbsce)) then
        i_dir=1
     end if
  end if
  i_fn_count=0
  i_co_count=0
  i_cnt_count=0
  !
  !  loop over relative energies
  do ien=1,3
     !  loop over beamlines
     do ib=1,nbeam
        if(xninja(ib).eq.0.0_fp) cycle
        !  h/he beams only
        if(xzbeama(ib).gt.2.05_fp) cycle
        iznb=xzbeama(ib)+0.1_fp
        !  zero sums for 2d nb*sigma*v calculation (arrays have beamline index)
        zn02=0.0_fp
        zrati2=0.0_fp
        zratx2=0.0_fp
        a_rds=abeam*abeama(ib)/(abeam+abeama(ib)) !reduced mass
        !  loop over ingoing/outgoing
        do ino=1,2
           !  check for neutral density
           zn0=bn002(ib,ien,ino,nflr2)
           if(zn0.eq.0.0_fp) cycle
           zn0=zn0/bmvol(nflr2)
           !  evaluate relative energy from rel. velocity
           zv2=0.0_fp
           do idir=1,3
              !  directional compenent of avg 1st gen neutral velocity in zone
              zv0dir=viona(ien,ib)*vcxbn0(ib,idir,ino,nflr2)
              zv2=zv2+(zv0dir-zvi(idir))**2

           end do
           !  relative energy
           zerel=max(1.0_fp,(zv2*v2toev/abeam)) ![ev/amu]
           call cxnsum_intrp(zerel,0,irange_bb,ind,indp1,zintrp)
           !  (in ev/nucleon); get sigma*v's
           !dbg               call hatom_sigv(icx,0.001*zerel,1,iznb,izb,
           !dbg     >            zsigv_cx,ier)
           !dbg               zsigv_cx=1.0e6*zsigv_cx
           zsigv_cx=(1.0_fp-zintrp)*cxn_bbcx(iznb,ind,izb) + &
                zintrp*cxn_bbcx(iznb,indp1,izb)
           zsvdtg=dtgcx*zsigv_cx
           !  include excited states correction factor
           zdtscx=zn0*zsvdtg*hexbbcx0(ib,ien,ino,nflr2,lsbeam)
           ztscx=ztscx+zdtscx
           if(isflag.eq.1) then
              snbxbb0(lsbeam)=snbxbb0(lsbeam)+zdtscx*zwght
              !  nb*sigma*v sums - neutral density
              zn02=zn02+zn0
              !  reaction rates (no excited states correction)
              zdsncx=zsvdtg*zwght
              !dbg                  call hatom_sigv(ii,0.001*zerel,1,iznb,izb,
              !dbg     >               zsigv_ii,ier)
              !dbg                  zsigv_ii=1.0e6*zsigv_ii
              zsigv_ii=(1.0_fp-zintrp)*cxn_bbii(iznb,ind,izb) + &
                   zintrp*cxn_bbii(iznb,indp1,izb) !cm**3/sec
              zdsnii=dtgcx*zwght*zsigv_ii
              zrati2=zrati2+zn0*zdsnii
              zratx2=zratx2+zn0*zdsncx
              if((nflr2.gt.0).and.(nflr2.le.mimxbzf)) then !inside plasma only
                 if((ino.eq.1).and.(i_dir.eq.1)) then
                    ebb_rel_co(nflr2,ien,ib,lsbeam)= &
                         ebb_rel_co(nflr2,ien,ib,lsbeam)+zerel &
                         *dtgcx*a_rds*zwght !average e_rel  [ev/amu]
                    ebb_rel2_co(nflr2,ien,ib,lsbeam)= &
                         ebb_rel2_co(nflr2,ien,ib,lsbeam)+zerel &
                         *zerel*a_rds*a_rds*dtgcx*zwght !average  e_rel**2 [ev**2/amu**2]
                    if(i_co_count.eq.0) then
                       bnb_ien2(nflr2,i_dir,lsbeam)=bnb_ien2 &
                            (nflr2,i_dir,lsbeam)+dtgcx*zwght !weighting
                       i_co_count=i_co_count+1
                    end if
                 else if((ino.eq.2).and.(i_dir.eq.2)) then
                    ebb_rel_cnt(nflr2,ien,ib,lsbeam)= &
                         ebb_rel_cnt(nflr2,ien,ib,lsbeam)+zerel* &
                         dtgcx*a_rds*zwght !average e_rel  [ev/amu]
                    ebb_rel2_cnt(nflr2,ien,ib,lsbeam)= &
                         ebb_rel2_cnt(nflr2,ien,ib,lsbeam)+zerel &
                         *zerel*a_rds*a_rds*dtgcx*zwght !average  e_rel**2 [ev**2/amu**2]
                    if(i_cnt_count.eq.0) then
                       bnb_ien2(nflr2,i_dir,lsbeam)=bnb_ien2 &
                            (nflr2,i_dir,lsbeam)+dtgcx*zwght !weighting
                       i_cnt_count=i_cnt_count+1
                    end if
                 else if(i_dir.eq.0) then
                    ebb_rel_fn(nflr2,ien,ib,lsbeam)= &
                         ebb_rel_fn(nflr2,ien,ib,lsbeam)+zerel* &
                         dtgcx*a_rds*zwght/2 !average e_rel  [ev]
                    ebb_rel2_fn(nflr2,ien,ib,lsbeam)= &
                         ebb_rel2_fn(nflr2,ien,ib,lsbeam)+zerel &
                         *zerel*a_rds*a_rds*dtgcx*zwght/2. !average over ino  e_rel**2 [ev**2]
                    if(i_fn_count.le.0)then
                       bnb_fn(nflr2,lsbeam)=bnb_fn &
                            (nflr2,lsbeam)+dtgcx*zwght !weighting
                       i_fn_count=i_fn_count+1
                    end if
                 end if
              end if
           end if
        end do

        if((zn02.gt.0.0_fp).and.(isflag.eq.1)) then
           !  divide out neutral density for nb*sigma*v
           !  (excited states effect not included here; to be counted later)
           sn00i2(ib,ien,nflr2,lsbeam) = sn00i2(ib,ien,nflr2,lsbeam) &
                + zrati2/zn02/bmvol(nflr2)! 1/sec
           sn00x2(ib,ien,nflr2,lsbeam) = sn00x2(ib,ien,nflr2,lsbeam) &
                + zratx2/zn02/bmvol(nflr2)! 1/sec
        end if
     end do
  end do
  !
  !  exit
  !
250 continue
  if(isflag.eq.1) then
     snbxtot(lsbeam)=snbxtot(lsbeam)+ztscx*zwght
     snbx1d(ngc,lsbeam)=snbx1d(ngc,lsbeam)+ztscx*zwght
     !         snbx1d(nflr,lsbeam)=snbx1d(nflr,lsbeam)+ztscx*zwght
     !
     !
     !  pitch angle, energy indices (nflr2 = spatial index, cf nbiflr)
     jba=int((1.0_fp+xksid)*0.5_fp*real(nznbma))+1
     jba=min(nznbma,max(1,jba))
     zep = v2toev*ep
     !?mg or zep = v2toev*vion*vion
     jbe=0
     jbe=ndxfbe_r8(zep,lsbeam,jbe)
     if(jbe <= 0.or.jbe > mimxbe) then
        if(orbrzv_rf_option.ne.1.and. &
             (ntest_max.eq.0.or.abs(ntest_max).eq.10)) then
           call errorlog(' beamcx, ndxfbe_r8 : out of energy',zep)
        end if
        jbe = max(1,jbe)
        jbe = min(mimxbe,jbe)
     end if
     ! accumulate a guiding center oriented cx fast neutral source distribution vs. pitch, energy, and 2d position
     if(ngc2.gt.0.and.ngc2.le.mimxbz) then     !may include some zone rows beyond the plasma boundary
        ! collect cx distribution function if  ncx0 is more than zero
        if(ncx0.gt.0) then
           fbm_cx(jbe,jba,ngc2,lsbeam) = &
                fbm_cx(jbe,jba,ngc2,lsbeam) + ztscx*zwght
           !
        end if
     end if
  end if
  !
  return
end subroutine cxnsum
!
!---------------------------------------
!  internal subroutine for table interpolation
!  (see ebatom & nbatom subroutines for table setups)
!
subroutine cxnsum_intrp(zeova,isbe,irng,ie,iep,zifac)
  use nbatom_mod,only : nbnsve,nbnsver,bnsves,lbnsve
  use logmod
  use nbi_com,only : orbrzv_rf_option,ntest_max,fp
  real(fp), intent(in) :: zeova ! e/a for table
  integer, intent(in) :: isbe ! table id (0 or by beam specie) or -1 for coarse energy grid
  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 :: inume,inumr,imin,imax,inum,inc
  real(fp) :: zea,zemin,zemax

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

  inumr=nbnsver(isbe)       ! no. of energy range sectors
  inume=nbnsve(isbe)        ! tot. no. energy grid pts

  zea=zeova
  if(zea.gt.bnsves(inume,isbe)) then
     ! kjm-rfdbg : warning disabled when rf-option is on
     if(orbrzv_rf_option.ne.1.and. &
          (ntest_max.eq.0.or.abs(ntest_max).eq.10)) then
        call errorlog &
             (' cxnsum_intrp: energy interaction exceeds table limit')
        call errorlog('     energy/a = ',zea)
        call errorlog('     limit = ',bnsves(inume,isbe))
        call errorlog('     irange = ',irng)
        call errorlog('     isbe = ',isbe)
     end if
     zea=bnsves(inume,isbe)
     iep=inume
     ie=iep-1
     zifac=1.0_fp
     return
  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 cxnsum_intrp

!******************** end file cxnsum.for ; group beamcx ***************
