!--------------------  nbains.fnb  --  group nubeam  ----------------
!----------------------------------------------------------------------
!  nbains - interpolate thermal plasma cross-section data from nbatom
!           lookup table
!

subroutine nbains( &
     zbnsvcx, zbnsvii, zbnsvie, zbnsviz, zbnsvto, zexcfac, &
     zbmsexc1d, zsigvnt1d,izone,ind,zintri,indp1,zintrp)
  !
  use nbatom_mod
  use nbi_com
  !
  !
  implicit none
  real(fp) :: zbnsvcx(mig)
  real(fp) :: zbnsvii
  real(fp) :: zbnsvie
  real(fp) :: zbnsviz
  real(fp) :: zbnsvto
  !
  real(fp) :: zexcfac
  real(fp) :: zbmsexc1d, zsigvnt1d
  !
  real(fp) :: zintri,zintrp
  integer :: izone,ind,indp1
  logical :: ilbms1d !ture - beam-beam rections are included and for adas beam
  !stopping coefficient and enhancement factor calculation
  !turn on
  !
  !
  !  input interpolation factors ind,... computed in subroutine nbeloc
  !
  !  input izone gives plasma zone index
  !
  !  output: zbnsvcx(ig), zbnsvii, zbnsvie, zbnsviz,
  !                                  zbnsvto - mig+4 numbers interpolated
  !                                          from nbatom tables.
  !       these numbers are fractions (cx, th.ion impact ionization,
  !          electron and impurity impact ionizations (dimensionless)
  !          and the total stopping rate 1/sec (sum of n*sigma*v's))
  !
  !
  !	--------------------------------------------------------------
  !
  ! updates:
  !
  !  dmc   5 may 1994 -- added zexcfac -- beam excitation effect factor
  !         set if nsigexc.gt.0
  !
  !  dmc  13 may 1993 -- added zbnsvie and zbnsviz to output set
  !
  !  tbt  oct. 17, 1991 - nbatom upgrade: got rid of bnsigv.
  !		        use instead bnsviif, bnsvcxf and bnsvtot.
  !		        took out znsigv & idim arguments.
  !     		        use instead zbnsvcx(mig), zbnsvii & zbnsvto
  !
  integer :: ig
  !
  integer :: ict=1                  ! set to 0 for initial debug output
  save ict
  !	---------------------------------------------------------------
  zbmsexc1d=0.0_fp
  zsigvnt1d=0.0_fp
  ilbms1d=(nlbbcx.and.(nsigexc_use.eq.1).and. &
       (lev_nbidep.eq.2))
  ! if nsigexc_use.eq.3 then n*sig*v are calculated in excited state model using adas310_fortran_driver (nbatom.for)
  if (ind .eq. 0) then
     !
     !	  .energy is between 0 and first grid point.
     !	  .fractions (zbnsvcx&zbnsvii) set to constant
     !                                           (value at 1st grid poiunt.)
     !         .total n*sigma*v interpolation from 0 to value at 1st grid pt

     do ig=1,ng
        zbnsvcx(ig) = bnsvcxf(izone,ig,indp1,lsbeam)
     end do

     zbnsvii = bnsviif(izone,indp1,lsbeam)
     zbnsviz = bnsvizf(izone,indp1,lsbeam)
     zbnsvie = bnsvief(izone,indp1,lsbeam)
     zbnsvto = bnsvtot(izone,indp1,lsbeam)*zintrp

     zexcfac = bnsvexc(izone,indp1,lsbeam)
     if(ilbms1d) then    !use excited state model and contribution from  beam-beam rection will be included
        zbmsexc1d = bms_exc_1d(indp1,izone,lsbeam)
        zsigvnt1d = sigvnt_1d(indp1,izone,lsbeam)
     end if
  else                              ! normal interpolation.

     do ig=1,ng
        zbnsvcx(ig) = bnsvcxf(izone,ig,ind,  lsbeam)*zintri + &
             bnsvcxf(izone,ig,indp1,lsbeam)*zintrp
        zbnsvcx(ig) =max(zbnsvcx(ig),0.0_fp)
     end do

     zbnsvii = bnsviif(izone,ind,  lsbeam)*zintri + &
          bnsviif(izone,indp1,lsbeam)*zintrp
     zbnsviz = bnsvizf(izone,ind,  lsbeam)*zintri + &
          bnsvizf(izone,indp1,lsbeam)*zintrp
     zbnsvie = bnsvief(izone,ind,  lsbeam)*zintri + &
          bnsvief(izone,indp1,lsbeam)*zintrp
     zbnsvto = bnsvtot(izone,ind,  lsbeam)*zintri + &
          bnsvtot(izone,indp1,lsbeam)*zintrp

     zbnsvii=max(zbnsvii,0.0_fp)
     zbnsviz =max(zbnsviz,0.0_fp)
     zbnsvie =max(zbnsvie,0.0_fp)
     zbnsvto =max(zbnsvto,0.0_fp)
     zexcfac = bnsvexc(izone,ind  ,lsbeam)*zintri + &
          bnsvexc(izone,indp1,lsbeam)*zintrp
     zexcfac = max(zexcfac,1.0_fp)

     if(ilbms1d) then       !use excited state model and contribution from  beam-beam rection will be included
        zbmsexc1d = bms_exc_1d(ind,izone,lsbeam)*zintri + &
             bms_exc_1d(indp1,izone,lsbeam)*zintrp

        zsigvnt1d = sigvnt_1d(ind,izone,lsbeam)*zintri + &
             sigvnt_1d(indp1,izone,lsbeam)*zintrp

        zbmsexc1d =max(zbmsexc1d,0.0_fp)
        zsigvnt1d =max(zsigvnt1d,0.0_fp)
     end if

  end if  ! ind

  ict=ict+1
  if(ict.eq.1) then
     write(nonlin,1001) izone,ind,lsbeam,zintri,zintrp, &
          zbnsvii,zbnsviz,zbnsvie,zbnsvto,zexcfac
1001 format(/'nbains first call:'/ &
          ' izone=',i4,' ind=',i4,' lsbeam=',i2, &
          ' zintri=',1pe11.4,' zintrp=',1pe11.4/ &
          ' bnsvii=',1pe11.4,' bnsviz=',1pe11.4, &
          ' bnsvie=',1pe11.4,' bnsvtot=',1pe11.4/ &
          ' excfac=',1pe11.4/)
  end if

  return
end subroutine nbains
