!--------------------------------------------------------------------
!  addtrap -- compute trapping fraction profile sum -- for the
!     slowing down population
!
subroutine addtrap(zwght,itrap,zeiav)
  !
  use nbi_com
  use xjasub_mod, only : zxidpi,impin,impout
  !
  !  auxilliary values returned from xjadep calls...
  !============
  implicit none
  integer :: itrap,ixmin,ixmax,ixiav,ie
  !============
  real(fp) :: zeiav,zwght,zxmin,zxmax
  real(fp) :: zxspani,zxwght,ze1,ze2
  real(fp) :: zxmp1,zxmp2
  !
  !  zxmp1 -- xi(midplane) intercept #1 (xjadep common)
  !  zxmp2 -- xi(midplane) intercept #2 (xjadep common)
  !  zwght -- ptcl weight
  !  itrap -- =1 if this is a trapped orbit; 0 otherwise.
  !  zeiav -- avg energy of orbit
  !
  !  output:  total orbit weights profile sum and trapped orbits weights
  !    profile sums are incremented.
  !
  !  if nerngfi is set, energy binned sums are also incremented.
  !
  !  this routine is called from orball after orbit timesteps --
  !  sums applied to slowing down mc ion populations.
  !
  !--------------------------------------------------------------------
  !
  zxmp1=zxidpi(impin)
  zxmp2=zxidpi(impout)
  zxmin=max(0.0_fp,min(zxmp1,zxmp2))
  zxmax=max(0.0_fp,max(zxmp1,zxmp2))
  !
  ixmin=min(lep1,int(lcentr+zxmin*nzones))
  ixmax=min(lep1,int(lcentr+zxmax*nzones))
  !
  if(ixmin.ne.ixmax) zxspani=1.0_fp/(zxmax-zxmin)
  !
  do ixiav=ixmin,ixmax
     if(ixmin.eq.ixmax) then
        !  orbit confined within a single zone
      	zxwght=1.0_fp
     else if(ixiav.eq.ixmin) then
      	zxwght=(xi(ixmin+1,1)-zxmin)*zxspani
     else if(ixiav.eq.ixmax) then
      	zxwght=(zxmax-xi(ixmax,1))*zxspani
     else
      	zxwght=(xi(ixiav+1,1)-xi(ixiav,1))*zxspani
     end if
     !  trapping fraction profile, entire slowing down population
     wtrap(ixiav,lsbeam)=wtrap(ixiav,lsbeam)+zwght*zxwght
     if(itrap.eq.1) then
        btrap(ixiav,lsbeam)=btrap(ixiav,lsbeam)+zwght*zxwght
     end if
     !
     !  trapping fraction profiles, within energy ranges
     do ie=1,nerngfi
        ze1=erngfi(ie,lsbeam)
        if(ie.lt.nerngfi) then
           ze2=erngfi(ie+1,lsbeam)
        else
           ze2=epsinv
        end if
        if((ze1.le.zeiav).and.(zeiav.le.ze2)) then
           wtrape(ixiav,ie,lsbeam)= &
                wtrape(ixiav,ie,lsbeam)+zwght*zxwght
           if(itrap.eq.1) then
              btrape(ixiav,ie,lsbeam)= &
                   btrape(ixiav,ie,lsbeam)+zwght*zxwght
           end if
        end if
     end do
  end do
  !
  return
end subroutine addtrap
