!-------------------------------------------------------
subroutine r8_nbsums_errck( &
     zpinja,zeinja,zffulla,zfhalfa, &
     ibeama,zninja,zminja, &
     zpinj,zeinj,zminj,zninj,zcofrac, &
     zpinjs,zeinjs,zminjs,zninjs,zcofracs, &
     iwarn,ierr)
  !
  use nbi_com
  !
  !
  !  using common information, compute injected power, max voltage, etc.,
  !  per specie and overall, given power and voltage of individual beams;
  !  also compute momentum and ptcls/sec for each beam.
  !
  !  data for each of the beams...
  !
  implicit none
  real(fp), intent(in) :: &
       zpinja(nbeam),zeinja(nbeam),zffulla(nbeam),zfhalfa(nbeam)
  !
  !  no. of beams with non-zero power
  !
  integer, intent(out) :: ibeama
  !
  !  momentum and #/sec injected, per beamline
  !
  real(fp), intent(out) :: zninja(nbeam),zminja(nbeam)
  !
  !  voltages... remember the past
  !
  real(fp), intent(inout) :: zeinj,zeinjs(nsbeam)  ! in volts
  !
  !  data (summed over all beams)
  !
  real(fp), intent(out) :: &
       zpinj,zminj,zninj,zcofrac
  !
  !  power,voltage,momentum,#ptcls,power-fraction-co
  !  watts, volts, g*cm*cm/sec, #/sec, ratio
  !
  !  data, by species (summed over beams for each given species)
  !
  real(fp), intent(out) :: &
       zpinjs(nsbeam),zminjs(nsbeam),zninjs(nsbeam), &
       zcofracs(nsbeam)
  !
  integer, intent(out) :: iwarn   ! >0 if einj < 10 kev
  integer, intent(out) :: ierr    ! >0 if ffull <= 0.0
  !
  !$r8real_input: zpinja,zeinja,zffulla,zfhalfa
  !$r8real_output: zninja,zminja,zpinj,zminj,zninj,zcofrac
  !$r8real_output: zpinjs,zminjs,zninjs,zcofracs
  !
  real(fp), parameter :: zzeinj_min = 10000.0_fp
  !
  integer :: ib
  real(fp) :: zzeinja(nbeam)
  !-------------------------------------
  !
  iwarn = 0
  ierr = 0
  !
  do ib=1,nbeam
     zzeinja(ib) = 0.0_fp
     if(zpinja(ib).le.0.0_fp) cycle

     zzeinja(ib) = max(zzeinj_min,zeinja(ib))
     if(zeinja(ib).lt.zzeinj_min) then
        write(nonlin,1001) ib,zzeinj_min,zeinja(ib)
        iwarn = iwarn + 1
     end if
     if(zffulla(ib).le.0.0_fp) then
        write(nonlin,1002) ib,zffulla(ib)
        ierr = ierr + 1
     end if
     if(zfhalfa(ib).lt.0.0_fp) then
        write(nonlin,1003) ib,zfhalfa(ib)
        ierr = ierr + 1
     end if
     if((zffulla(ib)+zfhalfa(ib)).gt.1.000001_fp) then
        write(nonlin,1004) ib,zffulla(ib)+zfhalfa(ib)
        ierr = ierr + 1
     end if
  end do
  !
1001 format(' %r8_nbsums_errck: beam #',i3,'  einj_min=', &
       1pe11.4,'  einj(actual) = ',1pe11.4,' (ev)')
  !
1002 format(' ?r8_nbsums_errck: beam #',i3,'  ffull = ',1pe11.4, &
       '; positive full energy fraction expected.')
  !
1003 format(' ?r8_nbsums_errck: beam #',i3,'  fhalf = ',1pe11.4, &
       '; non-negative half energy fraction expected.')
  !
1004 format(' ?r8_nbsums_errck: beam #',i3,'  ffull+fhalf = ',1pe11.4, &
       '; energy fraction sum <= 1.0 expected.')
  !
  if(ierr.ne.0) then
     ibeama=nbeam

     zninja=0.0_fp
     zminja=0.0_fp

     zpinj=0.0_fp
     zminj=0.0_fp
     zninj=0.0_fp
     zcofrac=0.0_fp

     zpinjs=0.0_fp
     zminjs=0.0_fp
     zninjs=0.0_fp
     zcofracs=0.0_fp

     return
  end if

  call r8_nbsums(zpinja,zzeinja,zffulla,zfhalfa, &
       ibeama,zninja,zminja, &
       zpinj,zeinj,zminj,zninj,zcofrac, &
       zpinjs,zeinjs,zminjs,zninjs,zcofracs)

  return
end subroutine r8_nbsums_errck
!-------------------------------------------------------
subroutine r8_nbsums(zpinja,zeinja,zffulla,zfhalfa, &
     ibeama,zninja,zminja, &
     zpinj,zeinj,zminj,zninj,zcofrac, &
     zpinjs,zeinjs,zminjs,zninjs,zcofracs)
  !
  use nbi_com
  use physconst_mod, only : evptocm_sec
  !
  !
  !  using common information, compute injected power, max voltage, etc.,
  !  per specie and overall, given power and voltage of individual beams;
  !  also compute momentum and ptcls/sec for each beam.
  !
  !  data for each of the beams...
  !
  implicit none
  real(fp), intent(in) :: &
       zpinja(nbeam),zeinja(nbeam),zffulla(nbeam),zfhalfa(nbeam)
  !
  !  no. of beams with non-zero power
  !
  integer, intent(out) :: ibeama
  !
  !  momentum and #/sec injected, per beamline
  !
  real(fp), intent(out) :: zninja(nbeam),zminja(nbeam)
  !
  !  voltages... remember the past
  !
  real(fp), intent(inout) :: zeinj,zeinjs(nsbeam)  ! in volts
  !
  !  data (summed over all beams)
  !
  real(fp), intent(out) :: &
       zpinj,zminj,zninj,zcofrac
  !
  !  power,voltage,momentum,#ptcls,power-fraction-co
  !  watts, volts, g*cm*cm/sec, #/sec, ratio
  !
  !  data, by species (summed over beams for each given species)
  !
  real(fp), intent(out) :: &
       zpinjs(nsbeam),zminjs(nsbeam),zninjs(nsbeam), &
       zcofracs(nsbeam)
  !
  !$r8real_input: zpinja,zeinja,zffulla,zfhalfa
  !$r8real_output: zninja,zminja,zpinj,zminj,zninj,zcofrac
  !$r8real_output: zpinjs,zminjs,zninjs,zcofracs
  !-------------------------------------
  !  this routine used to be part of datcor/nbpdat.for
  !
  !  separated as part of modularization effort -- dmc 10/01/01
  !
  !  rev-- may 1980 d. mc cune; pinj,einj, and cofrac are now
  !  defined each timestep from the input variables (arrays)
  !  nbeam,pinja,einja,nlco,pdelta,tbona,tboffa
  !
  !  einj is approximate-- used as normalizer in fokker
  !   precise energies and fractions of each beamline are used in dep
  !
  !  mod dmc jan 1990 -- added ibeams array and check on einjs(is) in
  !  410 loop to fix bug reported by jet in multispecies beam runs where
  !  one of the species goes away before the beams are fully off
  !
  !-----------------------------------------------
  !
  real(fp) :: zcofs(nsbeam),zcof
  real(fp) :: zeinjsp(nsbeam),zeinjp
  integer :: ibeams(nsbeam)
  real(fp) :: zbmcurs(nsbeam)
  integer :: is,ib,ib0,isb
  real(fp) :: zv0,zva,zfthrd,zcurnr,zincb
  real(fp) :: zvsinh,zvcosh,zrmass,zeinjl
  real(fp) :: zfac2,zfac3
  !
  !-----------------------------------------------
  !
  zninja=0.0_fp
  zminja=0.0_fp
  zpinj=0.0_fp
  zminj=0.0_fp
  zninj=0.0_fp
  !
  zeinjp=einj
  zeinj=0.0_fp
  !
  zcofrac=0.0_fp
  zcof=0.0_fp
  !
  ibeama=0
  !
  !  specie dependent powers, etc...
  !
  do is=1,nsbeam
     ibeams(is)=0
     zcofs(is)=0.0_fp
     zcofracs(is)=0.0_fp
     zeinjsp(is)=zeinjs(is)
     if(.not.nlfprod(is)) zeinjs(is)=0.0_fp
     zpinjs(is)=0.0_fp
     zminjs(is)=0.0_fp
     zninjs(is)=0.0_fp
     zbmcurs(is)=0.0_fp

  end do

  !
  do ib=1,nbeam
     !
     !  get species indices
     !
     is=0
     do isb=1,nsbeam
        if(nlfprod(isb)) cycle
        if((abs(abeama(ib)-abeams(isb)).lt.0.05_fp).and. &
             (abs(xzbeama(ib)-xzbeams(isb)).lt.0.05_fp)) is=isb
     end do
     if(is.eq.0) cycle
     !
     zninja(ib)=0.0_fp
     zminja(ib)=0.0_fp
     zeinj=max(zeinj,zeinja(ib))
     zeinjs(is)=max(zeinjs(is),zeinja(ib))
     !
     if(zpinja(ib).le.0.0_fp) cycle
     !
     ibeama=ibeama+1
     ibeams(is)=ibeams(is)+1
     zpinj=zpinj+zpinja(ib)
     zpinjs(is)=zpinjs(is)+zpinja(ib)
     !
     zv0=evptocm_sec*sqrt(zeinja(ib)/abeama(ib)) ! cm/sec
     !
     if(nlco(ib)) then
        zcof=zcof+zpinja(ib)
        zcofs(is)=zcofs(is)+zpinja(ib)
     else
        zv0=-zv0                    !
     end if
     !
     !  support trace beams...
     !
     ib0=abs(ntrace(ib))
     if(ib0.eq.0) ib0=ib
     !
     !  usually: 0.5:
     zfac2=abeama(ib0)/(abeama(ib0)+abeama(ib))
     !
     !  usually: 0.333333:
     zfac3=abeama(ib0)/(abeama(ib0)+2*abeama(ib))
     !
     zfthrd=1.0_fp-zffulla(ib)-zfhalfa(ib)
     zva=zv0* &
          (zffulla(ib)+zfhalfa(ib)*sqrt(zfac2)+ &
          zfthrd*sqrt(zfac3))
     !
     zcurnr=1.0_fp+ zfac2*zfhalfa(ib)/zffulla(ib) + &
          zfac3*zfthrd/zffulla(ib)
     !
     zincb=zpinja(ib)/(zeinja(ib)*zcurnr*zffulla(ib))
     !
     zbmcurs(is)=zbmcurs(is)+zincb
     !
     zninja(ib)=zincb/xe
     !
     !  for trace beams use parent beam index (ib0) for geometry information
     !  trace beam index (ib) for z,a,power,voltage,current,momentum source
     !
     zvsinh=(xybsca(ib0)-xybapa(ib0))/xlbapa(ib0)
     zvcosh=sqrt(max(0.0_fp,(1.0_fp-zvsinh**2)))
     zrmass=abeama(ib)*1.6726e-24_fp  ! grams
     zminja(ib)=zrmass*zva*zvcosh*rtcena(ib0)*zninja(ib)
     !
     zminj=zminj+zminja(ib)
     zminjs(is)=zminjs(is)+zminja(ib)
     !
  end do
  !
  !
  do is=1,nsbeam
     if(nlfprod(is)) cycle
     zninjs(is)=zbmcurs(is)/xe
     zninj=zninj+zninjs(is)
  end do
  !
  !  if nbeama=0, no more injection, but may have a few orbits left
  !  if all beams are off set einj=value from prev step for fokker
  !   mod dmc jan 1990 -- do this individually for each fast ion species!
  !
  !   mod dmc sept 1990 -- add code with zeinjl to make sure einj is
  !  positive even if nbpdat is called before any beams have been turned
  !  on.
  !
  !   the point of all this is to give sr fokker in the beam code a
  !  reference energy for evaluating log(lamda) effects in the collision
  !  operator.  a non-zero energy must be provided, but only order-of-
  !  magnitude accuracy is required
  !
  zeinjl=zeinj
  if(ibeama.eq.0) zeinj=zeinjp
  if(zeinj.le.0.0_fp) zeinj=zeinjp
  if(zeinj.le.0.0_fp) zeinj=max(1000.0_fp,zeinjl)
  do is=1,nsbeam
     if(nlfprod(is)) cycle
     if(ibeams(is).eq.0) zeinjs(is)=zeinjsp(is)
     if(zeinjs(is).le.0.0_fp) zeinjs(is)=max(1000.0_fp,zeinjsp(is))
  end do
  !
  if(ibeama.gt.0) then
     if(zpinj.gt.0.0_fp) zcofrac=zcof/zpinj
     do is=1,nsbeam
        if(nlfprod(is)) cycle
        if(zpinjs(is).gt.0.0_fp) zcofracs(is)=zcofs(is)/zpinjs(is)
     end do
  end if
  !
  return
end subroutine r8_nbsums
