!-----------------------------------------------------------------------
!  ebatom -- set up energy grid for nbatom atomic physics tables
!
!      .. dmc 14 oct 1991
!
!  mod dmc 27 jun 2001 -- using nbatom_mod module instead of trcom
!                         for output tables
!
!
!  note -- as with most atomic physics code in transp -- we are
!  dealing with energy *per nucleon*, ev/a, not actual particle
!  energies -- i.e. a 100 kev deuteron is 50 kev/nucleon -- 50 "kev"
!  in the bnsves table units -- this is the relevant "energy" for
!  doing the reaction rate table lookup.
!
!  this routine sets up a species dependent energy grid to be used by
!  nbatom to set up atomic physics (ionization + cx) reaction rate data
!  for use in the beam code deposition and charge exchange neutral
!  recapture calculations.
!
!  the energy grid is set up in multiple ranges of linearly spaced
!  points.  the maximum energy and number of energy ranges is dependent
!  on the fast ion species -- generally, fusion product species will
!  have a higher maximum energy and a greater number of energy ranges
!  than a beam species.
!
!  if we ever implement a monte carlo rf resonant beam ion species, we
!  may want to give it the energy range treatment of a fusion product...
!
!  output:
!   for each fast ion species isb,
!
!    n=nbnsve(isb) = no. of table energies
!    bnsves(1...n,isb) = the table energies **per nucleon**
!    nr=nbnsver(isb) = no. of table energy ranges
!    lbnsve(1...nr,isb) = pointer to maximum energy in the indicated
!                    energy range
!
!    for a given energy range jr, the minimum energy is 0.0 if jr=1, or
!    bnsves(lbnsve(jr-1,isb),isb) if jr.gt.1;
!
!    this range is subtended by lbnsve(jr,isb) pts if jr=1
!    or lbnsve(jr,isb)-lbnsve(jr-1,isb)+1 pts if jr.gt.1
!
!
!
subroutine ebatom
  !
  use nbatom_mod
  !
  use nbi_com
  !
  !
  implicit none
  real(fp) :: zerange(10)
  integer :: irange(10)
  !
  integer :: ii,ie,ir,iz,izh,ihandle,idum1,idum2
  integer :: indmin,inrange,infrac,isb,indmax
  real(fp) :: zelo,zehi,zelimit,zeova,zemaxa,zemaxl,zeminl,zfac,zfaci
  real(fp) :: zehi_fusn
  integer :: ntbb_pnts,nneut
  logical :: nl_fprod
  !-----------------------------------------------------------------------
  !
  !  the minimum of the energy limit in the ionization reaction rate data
  !  and that in the fusion tables is used as a constraining energy limit
  !  in this routine

  ii=2
  izh=1
  call hatom_handle(ii,izh,izh,ihandle)   ! h0,h ii "handle"
  call ii_coldtarget_erange(ihandle,zelo,zehi,idum1,idum2)
  zelimit=1.0e3_fp*zehi                ! ev
  !
  call fusn_handle_ddp(ihandle)
  call fs_coldtarget_erange(ihandle,zelo,zehi,idum1,idum2)
  zelimit=min(zelimit,1.0e3_fp*zehi)   ! ev
  !
  zehi_fusn=1.0e3_fp*zehi
  !
  !---------------------------------
  !  nbatom/nutrav/dep stuff
  !  allocate the arrays...
  if(allocated(nbnsve)) deallocate(nbnsve)
  if(allocated(lbnsve)) deallocate(lbnsve)
  if(allocated(nbnsver)) deallocate(nbnsver)
  !
  !  sizes, subrange information
  !
  nbnsvmx=400
  nbnsvmx_f=nznbme0
  inrange=4
  if(levmod_halo.eq.2) then
     nneut=nsbeam+ng !we need data for interaction of thermal neutrals with all plasma components
  else
     nneut=nsbeam    !only data for  beam neutrals  are needed
  end if
  !
  allocate(nbnsve(-1:nneut))
  nbnsve=0
  !
  allocate(lbnsve(inrange,-1:nneut))
  lbnsve=0
  !
  allocate(nbnsver(-1:nneut))
  nbnsver=0
  !
  !  energies and velocities of neutral beams (all fractions)
  !
  if(allocated(viona)) deallocate(viona)
  if(allocated(eiona)) deallocate(eiona)
  !
  infrac=3                          ! no. of energy fractions
  allocate(viona(infrac,nbeam))
  viona=0.0_fp
  allocate(eiona(infrac,nbeam))
  eiona=0.0_fp
  !
  !  the actual energy grid
  !
  if(allocated(bnsves)) deallocate(bnsves)
  allocate(bnsves(nbnsvmx,-1:nneut))
  bnsves=0.0_fp
  !
  !  total stopping n*sigma*v summed over thermal species
  !  excitation correction
  !
  if(allocated(bnsvtot)) deallocate(bnsvtot)
  if(allocated(bnsvexc)) deallocate(bnsvexc)
  allocate(bnsvtot(lep1,nbnsvmx,nneut))
  allocate(bnsvexc(lep1,nbnsvmx,nneut))
  bnsvtot=0.0_fp
  bnsvexc=0.0_fp
  if((nsigexc.eq.1.or.nsigexc.eq.3).and.(lev_nbidep.eq.2)) then
     if(allocated(bms_exc_1d)) deallocate(bms_exc_1d)
     allocate(bms_exc_1d(nbnsvmx,lep1,nneut))
     bms_exc_1d=0.0_fp
     !
     if(allocated(sigvnt_1d)) deallocate(sigvnt_1d)
     allocate(sigvnt_1d(nbnsvmx,lep1,nneut))
     sigvnt_1d=0.0_fp
     !
     if(nlbbcx) then
        if(nbbcx_avg.eq.0) then !use non average over flux surface ebb_rel_avg as energy grid

           ntbb_pnts=mimxbzf
        else                !use average over flux surface ebb_rel_avgx as energy grid

           ntbb_pnts=lep1
        end if
        if(allocated(bms_exc_2d)) deallocate(bms_exc_2d)
        allocate(bms_exc_2d(ntbb_pnts,3,nbeam))
        bms_exc_2d=0.0_fp
        !
        if(allocated(sigvnt_2d)) deallocate(sigvnt_2d)
        allocate(sigvnt_2d(ntbb_pnts,3,nbeam))
        sigvnt_2d=0.0_fp

     end if
  end if

  !
  !  fraction of stopping, by thermal reaction type:
  !
  if(allocated(bnsviif)) deallocate(bnsviif)
  if(allocated(bnsvief)) deallocate(bnsvief)
  if(allocated(bnsvizf)) deallocate(bnsvizf)
  if(allocated(bnsvcxf)) deallocate(bnsvcxf)
  allocate(bnsviif(lep1,nbnsvmx,nneut))
  allocate(bnsvief(lep1,nbnsvmx,nneut))
  allocate(bnsvizf(lep1,nbnsvmx,nneut))
  allocate(bnsvcxf(lep1,ng,nbnsvmx,nneut))
  bnsviif=0.0_fp
  bnsvief=0.0_fp
  bnsvizf=0.0_fp
  bnsvcxf=0.0_fp
  !
  !  beam-beam n*sigma*v's ... these are on the mc codes 2d irregular
  !  (rho,theta) grid
  !
  if(allocated(bbnsvcx)) deallocate(bbnsvcx)
  if(allocated(bbnsvii)) deallocate(bbnsvii)
  allocate(bbnsvcx(nfbznsi,nbnsvmx,nsbeam,nsbeam))
  allocate(bbnsvii(nfbznsi,nbnsvmx,nsbeam,nsbeam))
  bbnsvcx=0.0_fp
  bbnsvii=0.0_fp
  !
  !---------------------------------
  !  cxnsum stuff...
  !
  !  find min and max(z) of any normal fast species
  cxn_zmin=0  ! integer ::
  cxn_zmax=0  ! integer ::
  do isb=1,nsbeam
     if(.not.nlheavy(isb)) then
        if(cxn_zmin.eq.0) then
           cxn_zmin=xzbeams(isb)+0.1_fp
           cxn_zmax=cxn_zmin
        else
           iz=xzbeams(isb)+0.1_fp
           cxn_zmin=min(iz,cxn_zmin)
           cxn_zmax=max(iz,cxn_zmax)
        end if
     end if
  end do
  !
  !  for beyond-the-edge neutrals (dn0out)
  if(allocated(cxn_thcx_a)) deallocate(cxn_thcx_a)
  allocate(cxn_thcx_a(ng,nbnsvmx,nsbeam))
  cxn_thcx_a=0
  !
  !  for "wall" neutrals
  if(allocated(cxn_thcx_wa)) deallocate(cxn_thcx_wa)
  allocate(cxn_thcx_wa(ng,lep1,nbnsvmx,nsbeam))
  cxn_thcx_wa=0
  !
  if(allocated(cxn_thii_wa)) deallocate(cxn_thii_wa)
  allocate(cxn_thii_wa(ng,lep1,nbnsvmx,nsbeam))
  cxn_thii_wa=0
  !
  !  for "halo" neutrals
  if(allocated(cxn_thcx_ha)) deallocate(cxn_thcx_ha)
  allocate(cxn_thcx_ha(ng,nfbznsi,nbnsvmx,nsbeam))
  cxn_thcx_ha=0
  !
  if(allocated(cxn_thii_ha)) deallocate(cxn_thii_ha)
  allocate(cxn_thii_ha(ng,nfbznsi,nbnsvmx,nsbeam))
  cxn_thii_ha=0
  !
  ! for nth generation neutrals
  !
  if(allocated(cxn_bbcx_2d)) deallocate(cxn_bbcx_2d) !cx n*sig*v table for plasma fast neutral of nth generation
  allocate(cxn_bbcx_2d(nsbeam,nfbznsi,nbnsvmx,nsbeam))
  cxn_bbcx_2d=0
  !
  if(allocated(cxn_bbii_2d)) deallocate(cxn_bbii_2d)!cx n*sig*v table for plasma fast neutral of nth generation
  allocate(cxn_bbii_2d(nsbeam,nfbznsi,nbnsvmx,nsbeam))
  cxn_bbii_2d=0
  !
  !  for fast neutrals
  !
  if(allocated(cxn_bbcx)) deallocate(cxn_bbcx)
  if(allocated(cxn_bbii)) deallocate(cxn_bbii)
  if(cxn_zmin.gt.0) then
     allocate(cxn_bbcx(cxn_zmin:cxn_zmax,nbnsvmx,cxn_zmin:cxn_zmax))
     cxn_bbcx=0
     allocate(cxn_bbii(cxn_zmin:cxn_zmax,nbnsvmx,cxn_zmin:cxn_zmax))
     cxn_bbii=0
  end if
  !
  !---------------------------------
  !  btfusn stuff...
  !
  if(allocated(btfus_dt)) deallocate(btfus_dt)
  if(allocated(btfus_d3)) deallocate(btfus_d3)
  if(allocated(btfus_ddp)) deallocate(btfus_ddp)
  if(allocated(btfus_ddn)) deallocate(btfus_ddn)
  if(allocated(btfus_td)) deallocate(btfus_td)
  if(allocated(btfus_tt)) deallocate(btfus_tt)
  if(allocated(btfus_3d)) deallocate(btfus_3d)
  !
  if(nbdeut.gt.0) then
     !  d beam present
     allocate(btfus_ddp(lep1,nbnsvmx))
     btfus_ddp=0.0_fp
     !
     allocate(btfus_ddn(lep1,nbnsvmx))
     btfus_ddn=0.0_fp
     !
     if(max(ntrit,nrftrit).gt.0) then
        !  thermal t may be present
        allocate(btfus_dt(lep1,nbnsvmx))
        btfus_dt=0.0_fp
        !
     end if
     !
     if(max(nhel3,nrfhel3).gt.0) then
        !  thermal he3 may be present
        allocate(btfus_d3(lep1,nbnsvmx))
        btfus_d3=0.0_fp
        !
     end if
     !
  end if
  !
  if(max(nbtrit,nftrit).gt.0) then
     !  t beam or fusion product present
     allocate(btfus_tt(lep1,nbnsvmx))
     btfus_tt=0.0_fp
     !
     if(max(ndeut,nrfdeut).gt.0) then
        !  thermal d may be present
        allocate(btfus_td(lep1,nbnsvmx))
        btfus_td=0.0_fp
        !
     end if
     !
  end if
  !
  if(max(nbhel3,nfhel3).gt.0) then
     !  he3 beam or fusion product present
     !
     if(max(ndeut,nrfdeut).gt.0) then
        !  thermal d may be present
        allocate(btfus_3d(lep1,nbnsvmx))
        btfus_3d=0.0_fp
        !
     end if
     !
  end if
  !
  !---------------------------------
  !
  do isb=0,nneut

     nl_fprod=.false.
     if(isb.le.nsbeam.and.isb.ne.0) nl_fprod=nlfprod(isb)
     !
     if(isb.eq.0) then
        !     energy grid for beam target fusion tables...
        !
        inrange=4
        !
        irange(1)=nbnsvmx/inrange
        zerange(1)=zehi_fusn*0.001_fp
        irange(2)=(2*nbnsvmx)/inrange
        zerange(2)=zehi_fusn*0.01_fp
        irange(3)=(3*nbnsvmx)/inrange
        zerange(3)=zehi_fusn*0.1_fp
        irange(inrange)=nbnsvmx
        zerange(inrange)=zehi_fusn*0.9999999_fp
        !

     else if(nl_fprod) then
        !
        !  fusion product species:  will have four ranges covering the
        !  entire raw reaction table energy range as follows: 1/4 from
        !  0 to 0.1%, 1/4 from 0.1% to 1%, 1/4 from 1% to 10% and 1/4 of
        !  the points will cover the range 10% to 100% of the table limit.
        !
        !
        inrange=4
        !
        irange(1)=nbnsvmx/inrange
        zerange(1)=zelimit*0.001_fp
        irange(2)=(2*nbnsvmx)/inrange
        zerange(2)=zelimit*0.01_fp
        irange(3)=(3*nbnsvmx)/inrange
        zerange(3)=zelimit*0.1_fp
        irange(inrange)=nbnsvmx
        zerange(inrange)=zelimit*0.9999999_fp
        !
     else if(isb.gt.nsbeam) then
        !
        !  thermal  species:  will have two linearly spaced grid ranges, one from
        !  0 energy to an energy / nucleon of order the injection energy,
        !  and a second range to cover substantially higher energies.
        zemaxa= 3.5_fp*maxval(ti)/aplasm(isb-nsbeam)
        inrange=2
        !
        irange(1)=3*nbnsvmx/5
        zerange(1)=zemaxa
        irange(2)=nbnsvmx
        zerange(2)= &
             min(max(10.0_fp,edbfac)*zemaxa,zelimit*0.9999999_fp)
     else
        !
        !  beam species:  will have two linearly spaced grid ranges, one from
        !  0 energy to an energy / nucleon of order the injection energy,
        !  and a second range to cover substantially higher energies.
        !
        !  kjm-rfdbg
        if(orbrzv_rf_option.eq.1.or.abs(ntest_max).eq.1) then
           zeova=18.0_fp*einj/abeams(isb)
        else
           zeova=1.5_fp*einj/abeams(isb)
        end if
        !zeova=1.5_fp*einj/abeams(isb)
        zemaxa=max(75.0e3_fp,zeova)
        !
        inrange=2
        !
        irange(1)=3*nbnsvmx/5
        zerange(1)=zemaxa
        irange(2)=nbnsvmx
        zerange(2)= &
             min(max(10.0_fp,edbfac)*zemaxa,zelimit*0.9999999_fp)
        !
     end if
     !
     !  insert in common
     !
     nbnsve(isb)=nbnsvmx
     nbnsver(isb)=inrange
     !
     indmax=0
     zemaxl=0.0_fp
     do ir=1,inrange
        !
        zeminl=zemaxl
        zemaxl=zerange(ir)
        !
        indmin=indmax+1
        indmax=irange(ir)
        !
        lbnsve(ir,isb) = indmax
        !
        do ie=indmin,indmax
           zfac=real(ie+1-indmin)/real(indmax+1-indmin)
           zfaci=1.0_fp-zfac
           bnsves(ie,isb)=zeminl*zfaci+zemaxl*zfac

        end do
        !

     end do                        ! energy range loop
     !

  end do                          ! species loop



  !     energy grid for  fast ion fast neutral of generation>'0'
  inrange=4
  !
  irange(1)=nbnsvmx_f/inrange
  zerange(1)=zehi_fusn*0.001_fp
  irange(2)=(2*nbnsvmx_f)/inrange
  zerange(2)=zehi_fusn*0.01_fp
  irange(3)=(3*nbnsvmx_f)/inrange
  zerange(3)=zehi_fusn*0.1_fp
  irange(inrange)=nbnsvmx_f
  zerange(inrange)=zehi_fusn*0.9999999_fp
  nbnsve(-1)=nbnsvmx_f
  nbnsver(-1)=inrange
  zemaxl=0.0_fp
  indmax=0
  do ir=1,inrange
     !
     zeminl=zemaxl
     zemaxl=zerange(ir)
     !
     indmin=indmax+1
     indmax=irange(ir)
     !
     lbnsve(ir,-1) = indmax
     !
     do  ie=indmin,indmax
        zfac=real(ie+1-indmin)/real(indmax+1-indmin)
        zfaci=1.0_fp-zfac
        bnsves(ie,-1)=zeminl*zfaci+zemaxl*zfac
     end do
     !
  end do                     ! energy range loop
  !
  !
  lun_nbatom_mod=nonlin
  !
  return

end subroutine ebatom
! 11sep2004 fgtok
