!-----------------------------------------------------------------------
!  fuspwr -- add up source power and no. of fusion product fast ions
!
subroutine fuspwr(ident)
  !
  !  changed for uda zone volumes - pam stubberfield 14 july 1994
  !
  !  ident -- identifies the fusion product fast ion species
  !    (see subroutine fuslis)
  !
  use xstraln_calls
  use nbi_com
  !
  !  local:
  !
  !============
  implicit none
  integer :: ident,ireact,jz2,iz2,izr,iz,jz1,jz2p1,idummy,igc,iz1
  integer :: inth,iwarn
  !============
  real(fp) :: zomega,zxi,fnbint,zdth,zerot,zth,zr,zy,zvrot
  !============
  real(fp) :: zjac(2,2)
  !
  !  set .true. for reactions that produce a desired fusion product:
  logical ilreact(mireact)
  !
  !  initial energy (ev) of fusion product so produced:
  real(fp) :: zereact(mireact)
  !
  !-------------------------------------------------------------
  !
  !  get the list of reactions with this fusion product, and the product
  !  energy from these reactions
  !
  call fuslis(ident,ilreact,zereact)
  !
  !  the source rate per reaction has already been computed; do sums here.
  !
  pinjs(ident)=0.0_fp
  xninjs(ident)=0.0_fp
  !
  if(nlorbo) return  ! if debug switch for orbitting only is set.
  !
  do ireact=1,nreact
     if(ilreact(ireact)) then
        !
        jz2=lcm1
        iz2=nthzsm(lcm1)
        izr=lcm1
        !
        do iz=1,nfbzns
           !
           !  ptcl sum
           !
           xninjs(ident)=xninjs(ident)+sfusion(iz,ireact)
           !
           !  power sum to take a rotation correction into account
           !
           if(iz.gt.iz2) then
              !  new zone row, new rotation velocity:
              jz1=min(lep1,(jz2+1))
              jz2=min(lep1,(jz1+nxskpb-1))
              if(jz1.eq.lep1) then
                 zomega=0.0_fp
              else
                 jz2p1=min(lep1,(jz2+1))
                 zxi=0.5_fp*(xi(jz1,1)+xi(jz2p1,1))
                 zomega = fnbint(zxi,omegag,1,mj,1,idummy,igc)
              end if
              !  new theta indexing
              izr=izr+1
              iz1=iz2+1
              iz2=nthzsm(izr)
              inth=iz2-iz1+1
              zdth=(thbdy1-thbdy0)/inth
           end if
           !
           if(zomega.eq.0.0_fp) then
              zerot=0.0_fp
           else
              zth=zdth*(0.5_fp+(iz-iz1)) + thbdy0
              call xmomry(zxi,zth,0,zr,zy,zjac)
!              call xs_rzget(0.0_fp,zth,zxi,100.0_fp,zr,zy,iwarn)
              zvrot=zr*zomega
              zerot=v2toev*zvrot*zvrot
           end if
           !
           pinjs(ident)=pinjs(ident)+ &
                sfusion(iz,ireact)*xe*(zereact(ireact)+zerot)

        end do
     end if

  end do
  !
  if(xninjs(ident).gt.0.0_fp) then
     einjs(ident)=pinjs(ident)/(xe*xninjs(ident))
  end if
  !
  return
end subroutine fuspwr
