!------------------------------ nbain2.fnb ----- group nubeam ----------
!
!  nbain2 - interpolate beam-beam cross-section data from nbatom
!           lookup table
!

subroutine nbain2(zbbnsva,izone,zvtor,zvr,zvy,irng)
  !
  use nbatom_mod
  use nbi_com
  !
  implicit none
  real(fp) :: zbbnsva(2,mibs)
  integer :: irng(mibs)

  integer :: izone
  real(fp) :: zvtor,zvr,zvy
  !
  !  input interpolation factors ind,... computed in subroutine nbeloc
  !
  !  input izone gives plasma zone index
  !  input: irng - 1st guess to range # of energy grid.
  !
  !  output: zbbnsva
  !          array of numbers
  !          interpolated from nbatom bbnsvcx & bbnsvii tables.
  !           zbbnsva(1,isb) -- ionization by fast species isb
  !           zbbnsva(2,isb) -- cx by fast species isb
  !
  !  output: irng - actual range # in energy grid.
  !
  !  updates:
  !
  !  dmc may 1993 -- output argument changed, prep work for excited
  !                  states model.
  !
  !  tbt nov 1991 -- put in call to nbelc2 (instead of nbeloc). test to
  !		   see is one of the species is fusion product and if
  !                  it is use the appropriate energy grid.
  !
  !  tbt oct 1991 -- nbatom upgrade. eliminate bbnsig. define bbnsvcx and
  !                  bbnsvii. new output arguments zbbnsvi and zbbnsvt.
  !		   defined irng.
  !
  !  dmc aug 1988 -- generalized for multispecies fast ion populations
  !   this means multiple average velocity frames (cf vbtr2p array...)
  !
  integer :: isb,isbeam,ind,indp1
  real(fp) :: zv2,zeovab,zintrp,zintri
  !
  !	-------------------------------------------------------------
  !
  do isb=1,nsbeam
     if((nlfprod(isb).and..not.nlfatom).or.nlheavy(isb)) then
        zbbnsva(1:2,isb) = 0.0_fp
        cycle
     endif
     !	  .choose the grid that has the maximum energy per nucleon. tbt
     if (bnsves( nbnsve(isb),    isb) .gt. &
          bnsves( nbnsve(lsbeam), lsbeam)) then
        isbeam = isb
     else
        isbeam = lsbeam
     end if                          ! bnsves
     
     zv2=(zvtor-vbtr2p(izone,isb))**2+zvr**2+zvy**2
     zeovab=v2toev*zv2/abeam
     !
     call nbelc2(zeovab,irng(isbeam),isbeam,ind,zintri,indp1,zintrp)
     !
     if (ind .eq. 0)  then
        !	    .assume energy at grid pt to the left of first grid pt is 0.
        
        zbbnsva(1,isb) = bbnsvii(izone,indp1,lsbeam,isb)*zintrp
        zbbnsva(2,isb) = bbnsvcx(izone,indp1,lsbeam,isb)*zintrp
        !
     else                           ! normal interpolation.
        
        zbbnsva(1,isb) = &
             bbnsvii(izone,ind,  lsbeam,isb)*zintri + &
             bbnsvii(izone,indp1,lsbeam,isb)*zintrp
        !
        zbbnsva(2,isb) = &
             bbnsvcx(izone,ind,  lsbeam,isb)*zintri + &
             bbnsvcx(izone,indp1,lsbeam,isb)*zintrp
        
        zbbnsva(1,isb)=max(zbbnsva(1,isb),0.0_fp)
        zbbnsva(2,isb)=max(zbbnsva(2,isb),0.0_fp)
        !
     end if                         ! ind
     cycle
  end do
  !
  return
end subroutine nbain2
