!-----------------------------------------------------------------------
subroutine nbi_tausl
! nbi_tausl -- compute a volume average tau slowing down, using power 
! weighted average e_inj over beam lines (full energy fraction only), and 
! ignoring any plasma frame correction due to rotation velocity
  use nbi_com
  use physconst_mod, only : evptocm_sec
  implicit none
  integer j,ib,ireact,iz
  real(fp) ze,zvrel,zcf,ztausl,volzc_tot,zpw
!
! initial energy (ev) of fusion product so produced:
  real(fp) zereact(mireact)
  real(fp) znmax,zn
  logical ilreact(mireact)
  real(fp) ::  evptocm_sec2
!
!-----------------------------------------------------------------------
!
!  only do this calculation on the first call to depall
!
!
  evptocm_sec2= evptocm_sec*evptocm_sec
  if(.not.nlfprod(lsbeam)) then
!finding   power weighted average e_inj over beam lines (full energy fraction only)
     ze=0.0_fp
     zpw=0.0_fp
     do ib=1,nbeam
        if((abs(abeama(ib)-abeams(lsbeam)).lt.0.05_fp).and.&
             &         (abs(xzbeama(ib)-xzbeams(lsbeam)).lt.0.05_fp)) then
           ze=ze+einja(ib)*pinja(ib)
           zpw=zpw+pinja(ib) !total power for species
        endif
     enddo
     if(zpw.eq.0.0_fp) then
        ze=0.0_fp
     else
        ze=ze/zpw ![ev]
     endif
     zvrel=sqrt(ze*evptocm_sec2/abeams(lsbeam))  ! taken in the plasma frame
                                                    ! ignoring any plasma frame correction 
                                                    ! due to rotation velocity
  else !for fusion product
     call inispc !for fusion product
     call fuslis(lsbeam,ilreact,zereact)
!  choose the energy from the "dominant" fusion reaction
!
     znmax=-epsinv
!
     do ireact=1,nreact
        if(.not.ilreact(ireact)) cycle
!
!  total source from this reaction
!
        zn=0.0_fp
        do iz=1,nfbzns
           zn=zn+sfusion(iz,ireact)
        enddo
        if(zn.gt.znmax) then
           znmax=zn
           ze=zereact(ireact)
        endif
     enddo
     zvrel=sqrt(ze/v2toev)  ! taken in the plasma frame
  endif
!  volume averaging
!  the cf profiles are evaluated in subroutine fokker
!  
  volzc_tot=0.0_fp !total volume
  tausl_avg(lsbeam)=0.0_fp
  do  j=lcentr,ledge
     zcf=0.5_fp*(cfa(j,lsbeam,1)+cfa(j,lsbeam,2))
     ztausl=log((zvrel**3+vcrit3(j,lsbeam))/&
          &(vmin3(j,lsbeam)+vcrit3(j,lsbeam)))/(3.0_fp*zcf)
     tausl_avg(lsbeam)=tausl_avg(lsbeam)+ztausl*dvolzc(j)
     volzc_tot=volzc_tot+dvolzc(j)
  enddo
  if(volzc_tot.eq.0.0_fp) then
     tausl_avg(lsbeam)=0.0_fp
  else
     tausl_avg(lsbeam)=tausl_avg(lsbeam)/volzc_tot
  endif
!
end subroutine nbi_tausl
!-----------------------------------------------------------------------
subroutine nbi_tausl_vol_avg
! nbi_tausl_vol_avg -- compute a volume average tau slowing down to be used for next timestep,  
! restricted to full energy ions, and average over co and ctr for beam ions
  use nbi_com
  implicit none
  integer j
  real(fp) ztausl,volzc_tot,tausl_avg_prev
!
  tausl_avg_prev=tausl_avg(lsbeam) !save previous value
  volzc_tot=0.0_fp 
  tausl_avg(lsbeam)=0.0_fp
  do  j=lcentr,ledge
     if(.not.nlfprod(lsbeam)) then
        ztausl=0.5_fp*(taubsl(j,1,lsbeam)+taubsl(j,4,lsbeam))    ! average over co and ctr for beam ions,
     else
        ztausl= taubsl(j,1,lsbeam)    ! fusion product
                                      ! full energy only
     endif
     tausl_avg(lsbeam)=tausl_avg(lsbeam)+ztausl*dvolzc(j)
     volzc_tot=volzc_tot+dvolzc(j) !total volume
  enddo
  if(volzc_tot.eq.0.0_fp) then
     tausl_avg(lsbeam)=0.0_fp
  else
     tausl_avg(lsbeam)=tausl_avg(lsbeam)/volzc_tot
  endif
  if(tausl_avg(lsbeam).le.0.0_fp) then 
     tausl_avg(lsbeam)=tausl_avg_prev !return previous value, if 
                                      !taubsl... profiles havn't
                                      !been computed
  endif
end subroutine nbi_tausl_vol_avg
