!******************** start file nbatom.fnb ; group nbatom *************
!-----------------------------------------------------------------------
!
subroutine nbatom(ierr)
  !
  use nbi_com
  use physconst_mod, only : me_kg,mp_kg  
  use adas_mod, only : adas_bms,adas_btsigv,adas_zstop_sigvz, &
       adas_sigvte_ioniz,adas_sigv
  use adas_getsigs
  use adas_sigfacz
  use nbatom_mod
  !
  !============
  implicit none
  integer :: iabs
  !============
  real(fp) :: stat
  !============
  integer, intent(out) :: ierr   ! status code (0=ok)
  !
  !------------------------------------------------------------------
  !  nbatom - atomic physics for stopping neutral beam
  !
  !------> dmc 18 june 2001 -- convert to `preact' modular atomic physics
  !        package... hatom_btsigv instead of presig calls
  !        cf http://w3.pppl.gov/ntcc, catalog, `preact' module
  !
  !        dmc 28 june 2001 -- also, use nbatom_mod to hold nbi
  !        internal tables, and, calculate tables for cxnsum and
  !        btfusn, not just for dep/nutrav...
  !
  !........................
  !
  !            created from old nbatom.fnb 10/15/91  tbt
  !
  !    (former parts of dep and nutrav combined here - dmc mar 1985)
  !
  !  mod dmc may 1994 -- replace hexnb with sigfacz (sigmas) call
  !       new fitting approximation for excited states deposition effect
  !
  !  mod dmc jan 1992 -- output suppressed for heavy neutral beam species
  !
  !  mod tbt dec 10, 1991 - added ichrg=1, to argument list of presig to
  !                         indicate charged population and neutral
  !			  particle.
  !  mod tbt nov 20, 1991 - added ibmgrd array to store index of the beam
  !                         grid with the higher max energy for beam-beam.
  !  mod tbt oct  6, 1991 - added test for nlfatom.
  !  mod tbt sep 18, 1991 - switch backz and xzbeams in argument to presig
  !                         call.
  !  mod tbt may 14, 1991 - redo - change sigvbt calls to presig calls.
  !
  !
  !  mod dmc jan 1987 - take into acct possibility of rotating plasma
  !   neutral beam energy (in plasma frame for atomic physics collisions)
  !   becomes a function of plasma position, so that the beam indexed
  !   part of the output tables re. interaction with thermal plasma are
  !   not useful and are not computed.  however the beam-beam information
  !   is useful and should be used.  cf routines dep and nutrav
  !
  !  bugfix aug 1988 -- needed for proper modeling of rotation
  !   bnsigv and bbnsig now give sums of n*sigma*v(rel), not n*sigma
  !   when dep and nutrav use these quantities they should divide by
  !   the appropriate lab frame velocity, not known here, to get
  !   inverse mean free paths.
  !
  !  inputs:
  !   beam mass, injection energies
  !   plasma parameters (temperature and density profiles)
  !   beam-beam interaction data from previous beam timestep
  !
  !  outputs:
  !   atomic physics summary tables
  !     bnsviif, bnsvcxf, bnsvtot, bbnsvcx and bbnsvii.
  !   used by dep and nutrav.
  !
  !
  !----------------------
  !
  !  modif dmc dec 1984-- support excited state deposition (post-boley
  !    code hexnb).
  !    n*sigma*v rate calculation results now stored in common array
  !    bnsigv.  this includes a section of tabulated precalculated rates
  !    for use in nutrav subroutine for cx recapture
  !
  !------------------------------------
  !  organization of bnsviif, bnsvcxf, and bnsvtot.
  !
  !   bnsvief(i,ie,isb)   *fraction* of stopping rate due to electron
  !			impact ionizations
  !   bnsvizf(i,ie,isb)   *fraction* of stopping rate due to impurity
  !                       impact ionizations + cxs
  !   bnsviif(i,ie,isb)   *fraction* of stopping rate due to impact
  !                       ionizations by the thermal plasma ion species
  !   bnsvcxf(i,ig,ie,isb)*fraction* of neutral stopping rate which results
  !                       in species 1 or 2 hydrogenic source (i.e. cx
  !                       with background ions)
  !		        dimensionless units.
  !   bnsvtot(i,ie,isb)   total stopping rate (units:  n*sigma*v, sec**-1)
  !
  !     i   - plasma zone index
  !     ig  - (1 to ngmax):  thermal ion species
  !     ie  - energy grid point
  !     isb - species of beam.
  !
  !
  !  the energy range for nutrav is
  !  defined in subroutine ebatom, which is called in nbinit0.
  !  the table of energies for nutrav are in common array bnsves
  !
  !-----------------
  !  bbnsig(iz,j,isb,k)  is for beam-beam reactions for nutrav only
  !   units n*sigma*v, sec**-1
  !
  !   isb is the beam species index (new aug 88)
  !   iz is the 2d beam dist fcn zone spatial index
  !  rearranged dmc 8/88 for multispecies beam support
  !   j = 1:  beam beam impact ionization (with electron source)
  !   j = 2:  beam beam cx
  !            -- interaction with nth generation species isb beam
  !   neutrals *only*; i.e. used in nutrav only!.
  !   1st generation (dep) neutrals are handled via sn00i2/sn00x2
  !
  !   k  is the energy index (just as with bnsigv)
  !            ... d. mc cune
  !
  !-----------------
  ! lev_nbidep = 1	! beam "deposision" model level:
  ! 1 -- legacy of preact atomic physics
  ! 2 -- adas atomic physics
  ! 3 -- preact ground state atomic physics
  !      adas interaction with impurities
  ! 4 -- adas ground state atomic physics
  !      preact interaction with impurities
  ! nsigexc=0 means ground state model.
  ! nsigexc=1 means "appropriate applicable" excited states model:
  ! adas for lev_nbidep=2 and
  ! janev-boley-post for other values of lev_nbidep=1
  ! nsigexc=2 means janev-boley-post applied to any model including adas.
  ! janev/boley/post article, nucl. fusion vol.29, no. 12 (1989).
  ! nsigexc=3 means adas excited states model (adas310 library):
  ! adas for lev_nbidep=2
  !-----------------
  integer :: ibmgrd(mibs, mibs)          ! ibmgrd(i,j) = i if
  !                          ! bnsves(nbnsve(i),i) > bnsves(nbnsve(j),j)
  !                          !                   else = j.
  !
  integer :: ichrg      ! 1 for neutral particle and charged backgrnd
  !
  !------------------
  integer :: isb0,izb,iel,i,ig,iztarg,ier,ien,ij,inumi,inumi0
  integer :: ijx,isb,izbi,isbe,iz,ifrac,ib0,ix,ib,iz0,izi
  real(fp) :: ztia,ztix,zsigv,zfac,ze,zeova,zvi,zv2eva,zerel,zte
  real(fp) :: zsgvnt,zetran,zmemp,zdum1,zdum2
  real(fp) :: sigfacz,psigls_olson
  !
  real(fp), dimension(:), allocatable :: zearg,ztarg,zsigv_cx,zsigv_ii
  real(fp), dimension(:,:), allocatable :: zearg_adas
  real(fp), dimension(:,:,:), allocatable :: bmsexc2d
  real(fp), dimension(:,:,:), allocatable :: sigvnt2d
  real(fp), dimension(:,:), allocatable :: bmsexc1d
  real(fp), dimension(:,:), allocatable :: sigvnt1d
  real(fp), dimension(:), allocatable :: te_wrk,ne_wrk
  !
  integer, parameter :: icx=1, ii=2
  integer, parameter :: i_neutral=1, i_ion=2
  !
  real(fp), parameter :: zli = 3.0_fp
  !
  !------------------
  !  for sigfacz and sigmas excited states fit
  !
  logical :: ilhydr,nl_fprod,nl_heavy
  !
  integer :: imaxn,nneut
  real(fp),dimension(:), allocatable :: zzi,zti,zni,zza
  real(fp),dimension(:,:), allocatable :: zt_i,zn_i,zni_2d
  integer :: istat
  real(fp) :: z_beam,a_beam !current particle z and a
  integer :: i_nzone_min,i_nzone_max,i_zones,ik,ik2d,i_nzone_cl
  integer :: n_fast,i_fn
  integer,dimension(:), allocatable :: indx1d_2d
  integer :: i1,i2,n_cl
  real(fp) :: a_rds
  integer :: ix1,nrhix_s
  real(fp) :: rhix_s,xzimpx_s
#ifdef __DEBUG
  real(fp), dimension(:), allocatable :: zwk
#endif
  !
  !------------------------------------------------------------------
  !
  ierr=0
  !
  zmemp=me_kg/mp_kg ! me/mp
  !
  ichrg = 1          ! for neutral beam and charged background.

  call ebatom        ! set up energy grid.

  do isb=1,nsbeam
     if(nlheavy(isb)) cycle

     do isb0=1,nsbeam
        if(nlheavy(isb0)) cycle

        if ( bnsves(nbnsve(isb ),isb ) .gt. &
             bnsves(nbnsve(isb0),isb0))     then
           ibmgrd(isb,isb0) = isb
        else
           ibmgrd(isb,isb0) = isb0
        end if                       ! bnsves

     end do
  end do

  !  dmc oct 84-- loop over all beamlines plus all nutrav energies
  do ib=1,nbeam
     !         3 energy fractions for deposition
     ifrac=3
     !
     do ien=1,ifrac
        !
        !  depo stuff
        zfac=1.0_fp/real(ien)
        !  modify depo energy for trace element beam injection
        if((ien.gt.1).and.(ntrace(ib).ne.0)) then
           ib0=iabs(ntrace(ib))
           zfac=abeama(ib)/(abeama(ib)+(ien-1)*abeama(ib0))
        end if
        !
        ze=einja(ib)*zfac
        if(ze.eq.0.0_fp) ze=einj*zfac
        !
        eiona(ien,ib)=ze
        zeova=ze/abeama(ib) !mass per nucleon
        zvi=1.384e6_fp*sqrt(zeova) !v=sqrt(2*e/m)
        viona(ien,ib)=zvi
        !
     end do

  end do

  !        for nutrav-- loop over fast ion species
  !
  if(allocated(zearg)) deallocate(zearg)
  allocate(zearg(nbnsvmx))
  if(allocated(ztarg)) deallocate(ztarg)
  allocate(ztarg(nbnsvmx))
  if(allocated(zsigv_cx)) deallocate(zsigv_cx)
  if(allocated(zsigv_ii)) deallocate(zsigv_ii)
  allocate(zsigv_cx(nbnsvmx),zsigv_ii(nbnsvmx))
  !
  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
  !     find # impurity elements
  nrhix_s=0
  do ix = 1, nrhix
     ix1=min(ix+1,nrhix)
     if(ix.eq.ix1) then
        nrhix_s=nrhix_s+1
     end if
     if(aimpx(ix).eq.aimpx(ix1))then
        cycle
     else
        nrhix_s=nrhix_s+1
     end if
  end do
  !
  do isb0=1,nneut
     !mg         nsbeam

     !
     nl_fprod=.false.
     nl_heavy=.false.
     if(isb0.le.nsbeam) then
        nl_fprod=nlfprod(isb0)
        nl_heavy=nlheavy(isb0)
        izb = xzbeams(isb0)+0.1_fp
        a_beam = abeams(isb0)
     else
        izb = backz(isb0-nsbeam)+0.1_fp
        a_beam = aplasm(isb0-nsbeam)
     end if
     if(nl_fprod .and. .not. nlfatom) goto 900 ! tbt 11/91
     if(nl_heavy) goto 900
     !
     ilhydr = (izb.eq.1)            ! hydrogenic isotope?
     !-----------------------------------------------------------------
     !  plasma zone index - 1d profile data
     !
     !
     write(nonlin,*)' %nbatom deposition model: lev_nbidep=', &
          lev_nbidep
     write(nonlin,*)' %nbatom deposition model: nsigexc=', &
          nsigexc
     if(nsigexc.eq.3.and.ilhydr) goto 500 !n*sig*v for all reactions will be computed by adas310_fortran_driver
     !in excited state model, for hydrogenic isotopes only

     do i=lcentr,ledge

        !           interact with plasma ions
        do ig=1,ng
           iztarg=backz(ig)+0.1_fp
           ztia=ti(i,2)/aplasm(ig)
           !
           !     loop over all energy grid pts.
           do iel=1,nbnsve(isb0)
              zeova=bnsves(iel,isb0)
              zearg(iel)=0.001_fp*zeova
              ztarg(iel)=0.001_fp*ztia
           end do
           !
           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              !     .charge-exchange:
              call hatom_btsigv(icx,i_neutral,zearg,ztarg, &
                   nbnsve(isb0),izb,iztarg,zsigv_cx,ier)
              !     .impact-ionization
              call hatom_btsigv(ii,i_neutral,zearg,ztarg, &
                   nbnsve(isb0),izb,iztarg,zsigv_ii,ier)
           elseif(lev_nbidep.eq.2.or.lev_nbidep.eq.4) then !adas ground state atomic physics
              !     .charge-exchange:
              call adas_btsigv(icx,i_neutral,zearg,ztarg, &
                   nbnsve(isb0),izb,iztarg,zsigv_cx,ier)
              if(ier.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              !     .impact-ionization
              call adas_btsigv(ii,i_neutral,zearg,ztarg, &
                   nbnsve(isb0),izb,iztarg,zsigv_ii,ier)
              if(ier.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if
           !
           do iel=1,nbnsve(isb0)
              bnsvcxf(i,ig,iel,isb0) = &
                   rhob(i,ig,2) * 1.0e6_fp*zsigv_cx(iel)
              bnsviif(i,iel,isb0) = bnsviif(i,iel,isb0) + &
                   rhob(i,ig,2) * 1.0e6_fp*zsigv_ii(iel)
           end do


        end do
     end do                 ! i=lcentr,ledge

     !  interact with impurities
     do i=lcentr,ledge
        if (nlminsv.and.(nrhix.gt.1)) then
           ! -- treat each impurity separately --
           bnsvizf(i,1:nbnsve(isb0),isb0) = 0.0_fp
           rhix_s=0.0_fp
           xzimpx_s=0.0_fp
           do ix = 1, nrhix
              ix1=min(ix+1,nrhix)
              rhix_s=rhix_s+rhix(i,2,ix) !summed impurity density for all charge state
              xzimpx_s=max(xzimpx_s,xzimpx(ix)) !max charge state (corresponds fully stripped ion)
              if(aimpx(ix).eq.aimpx(ix1))then
                 if(ix.ne.ix1) then
                    cycle
                 end if
              end if
              do iel=1,nbnsve(isb0)
                 zeova=bnsves(iel,isb0)
                 zearg(iel)=0.001_fp* &
                      (zeova+1.5_fp*ti(i,2)/aimpx(ix)) ! kev/amu
              end do
              if((izb.eq.1).and.(nmsigx.eq.1)) then
                 do iel=1,nbnsve(isb0)
                    !     option for ancient olson cross section... watch units...
                    zsigv_ii(iel) = 1.384e6_fp* &
                         sqrt(1000.0_fp*zearg(iel))* &
                         psigls_olson(1000.0_fp*zearg(iel), &
                         xzimpx_s)*1.0e-6_fp
                 end do
              else
                 !     standard impurity stopping cross section
                 if((lev_nbidep.eq.1).or.(lev_nbidep.eq.4)) then !preact impurity stopping rates
                    call zstop_sigvz(izb,xzimpx_s,zearg, &
                         nbnsve(isb0),zsigv_ii)
                 elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.3)) &
                      then !adas impurity stopping rates
                    call adas_zstop_sigvz(izb,xzimpx_s,zearg, &
                         nbnsve(isb0),zsigv_ii,ier)
                    if(ier.eq.10) then
                       write(nonlin,*) &
                            '?nbatom->adas data is requested but not found'
                       goto 999
                    end if
                 end if

              end if
              do iel=1,nbnsve(isb0)
                 bnsvizf(i,iel,isb0) = bnsvizf(i,iel,isb0) + &
                      rhix_s*1.0e6_fp*zsigv_ii(iel)
              end do
              rhix_s=0.0_fp
              xzimpx_s=0.0_fp
           end do

        else
           !
           !     treat as one impurity...
           !
           do iel=1,nbnsve(isb0)
              zeova=bnsves(iel,isb0)
              zearg(iel)=0.001_fp*(zeova+1.5_fp*ti(i,2)/aimpj(i))
           end do
           if((izb.eq.1).and.(nmsigx.eq.1)) then
              do iel=1,nbnsve(isb0)
                 zsigv_ii(iel) = 1.384e6_fp* &
                      sqrt(1000.0_fp*zearg(iel))* &
                      psigls_olson(1000.0_fp*zearg(iel), &
                      xzimpj(i))*1.0e-6_fp
              end do
           else
              if((lev_nbidep.eq.1).or.(lev_nbidep.eq.4)) then !preact impurity stopping rates
                 call zstop_sigvz(izb,xzimpj(i),zearg, &
                      nbnsve(isb0),zsigv_ii)
              elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.3)) then !adas impurity stopping rates
                 call adas_zstop_sigvz(izb,xzimpj(i),zearg, &
                      nbnsve(isb0),zsigv_ii,ier)
                 if(ier.eq.10) then
                    write(nonlin,*) &
                         '?nbatom->adas data is requested but not found'
                    goto 999
                 end if
              end if
           end if
           do iel=1,nbnsve(isb0)
              bnsvizf(i,iel,isb0) = rhi(i,2)*1.0e6_fp* &
                   zsigv_ii(iel)
           end do

        end if
     end do                 ! i=lcentr,ledge


     !  interact with electrons.
     !  zmemp=me/mp -- approx correction to <sigma*v> due to beam energy
     !  now available -- dmc mar 1995
     do i=lcentr,ledge
        do iel=1,nbnsve(isb0)
           zeova=bnsves(iel,isb0)
           !
           if(nlebei) then
              zte=te(i,2)+0.6666667_fp*zmemp*zeova
           else
              zte=te(i,2)
           end if
           ztarg(iel)=0.001_fp*zte !kev
        end do
        !
        if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact e ii  atomic physics
           call sigvte_ioniz(izb,ztarg,nbnsve(isb0),zsigv_ii)
        elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas e ii atomic physics
           call adas_sigvte_ioniz(izb,ztarg,nbnsve(isb0),zsigv_ii, &
                ier)
           if(ier.eq.10) then
              write(nonlin,*) &
                   '?nbatom->adas data is requested but not found'
              goto 999
           end if
        end if
        !
        do iel=1,nbnsve(isb0)
           bnsvief(i,iel,isb0) = rhoel(i,2)*1.0e6_fp*zsigv_ii(iel)
        end do
     end do                 ! i loop
     !
     !  have all thermal reactions now...
     !  calculate (sigma*v*n) total and normalize ratios
     do i=lcentr,ledge
        do iel=1,nbnsve(isb0)
           zsgvnt=bnsviif(i,iel,isb0)+bnsvief(i,iel,isb0)+ &
                bnsvizf(i,iel,isb0)

           do ig=1,ng
              zsgvnt=zsgvnt+bnsvcxf(i,ig,iel,isb0)

           end do
           !
           !           .ratio making electron source
           bnsviif(i,iel,isb0) = bnsviif(i,iel,isb0)/zsgvnt
           bnsvief(i,iel,isb0) = bnsvief(i,iel,isb0)/zsgvnt
           bnsvizf(i,iel,isb0) = bnsvizf(i,iel,isb0)/zsgvnt

           do ig=1,ng
              !             .ratio making thermal neutral source (cx except beam-beam)
              bnsvcxf(i,ig,iel,isb0) =bnsvcxf(i,ig,iel,isb0)/zsgvnt

           end do
           !
           bnsvtot(i,iel,isb0) = zsgvnt
        end do               ! iel loop
        !
     end do                 ! i loop
     !
     !  mod dmc nov 1984 -- electron state excitation correction in
     !    neutral beams-- see also nutrav.for
     !
     !  mod dmc may 1994 -- hexnb replaced by sigmas fit --
     !    sigfacz -> sigfacx -> sigmas.  sigmas implements janev/post/boley
     !    fit described in 1989 nucl. fusion article, vol. 29., no. 12,
     !    pp 2125-2139 -- coefficients modified to extend the effective fit
     !    range down to the lower beam energies occurring in tftr and other
     !    real tokamaks.  the effect of the fast ion distributions in the
     !    target are handled in a very rough manner.
     !
     !    note -- only good for hydrogenic beams --
     !
500  continue
     nsigexc_use=nsigexc
     bnsvexc(:,:,isb0)=1.0_fp

     if((nsigexc.eq.1.or.nsigexc.eq.3).and. &
          ilhydr.and.(lev_nbidep.eq.2)) then !adas excited state data only for hydrogen isotopes
        !

        inumi0=ng+nrhix_s

        if((nbsteps.eq.0).or.(nl_fprod)) then
           n_fast=nsbeam
           imaxn=ng+nrhix_s+n_fast

           call contribute_step0(ierr)

        else if((nbsteps.gt.0).and.(.not.nl_fprod)) then
           i_fn=0
           do ij=1,nsbeam
              if(nlfprod(ij)) i_fn=i_fn+1
           end do

           n_fast=(nsbeam-i_fn)*2 !beam ptcls populations by direction of injection
           n_fast=n_fast+i_fn !fast ptcls including fusion products

           imaxn=ng+nrhix_s+n_fast

           call contribute_thermal(ierr)

        end if
        !
        if(ierr.ne.0) then
           write(nonlin,*) &
                ' ?? nbatom: no adas data for enhancement factor', &
                '            forced nsigexc_use = 0'
           nsigexc_use=0
        else
           if(nsigexc.eq.1.and.isb0.le.nsbeam) then !there is no needs to compute enhancement factor
              !for nsigexc.eq.3 - if excited state model are used
              if(nlbbcx.and.(nbbcx_bb.eq.0).and.(nbsteps.ne.0).and. &
                   (.not.nl_fprod)) then ! include contribution from beam-beam reactions
                 ! to the enhancement factor

                 if(nbbcx_avg.eq.0) then !use non average value for e_rel and t
                    i_nzone_min=1
                    i_nzone_cl=nfbznsi
                 else if(nbbcx_avg.gt.0) then
                    i_nzone_min=lcentr
                    i_nzone_cl=ledge
                 end if
                 i_zones=(i_nzone_cl-i_nzone_min+1)*nbeam*3
                 i_nzone_max=i_nzone_min+i_zones-1
                 !     allocate arrays
                 call alloc_wrk

                 if(nbbcx_avg.eq.0) then

                    allocate(indx1d_2d(nfbznsi))
                    indx1d_2d=0

                    call corresp1d_2d


                    !
                    !     expand density array from 1d radial grid to 2d mc grid. it is possible
                    !     because there is one-to-one correspondence between 2d mc grid index and 1d radial index,
                    do ik2d=1,nfbznsi
                       ik=indx1d_2d(ik2d)
                       if(ik.eq.0) cycle
                       ne_wrk(ik2d)=rhoel(ik,2)
                       te_wrk(ik2d)=te(ik,2)*0.001_fp
                       do ij=1,ijx
                          zni_2d(ij,ik2d)=zn_i(ij,ik)
                       end do
                    end do
                    deallocate(indx1d_2d)

                 else if(nbbcx_avg.gt.0)  then
                    ne_wrk(i_nzone_min:i_nzone_cl)= &
                         rhoel(i_nzone_min:i_nzone_cl,2)
                    te_wrk(i_nzone_min:i_nzone_cl)= &
                         te(i_nzone_min:i_nzone_cl,2)*0.001_fp
                    zni_2d(:,i_nzone_min:i_nzone_cl)= &
                         zn_i(:,i_nzone_min:i_nzone_cl)
                 end if



                 call contribute_fast

                 i1=i_nzone_min
                 i2=i_nzone_cl
                 n_cl=i_nzone_cl-i_nzone_min+1

                 do ib=1,nbeam

                    if(lsbeama(ib).ne.isb0) cycle

                    do ien=1,3
                       !mg                           if(nsigexc.eq.1)then
                       bms_exc_2d(i_nzone_min:i_nzone_cl,ien,ib)= &
                            bmsexc1d(i1:i2,1)
                       !mg                           elseif(nsigexc.eq.3)then
                       !mg                              bms_exc_2d(i_nzone_min:i_nzone_cl,ien,ib)=
                       !mg     >                             sum(bmsexc1d(i1:i2,:),dim=2)-
                       !mg     >                             bmsexc1d(i1:i2,1) !beam-beam total - ei
                       !mg                          end if
                       sigvnt_2d(i_nzone_min:i_nzone_cl,ien,ib)= &
                            sigvnt1d(i1:i2,1)+sigvnt1d(i1:i2,2)+ &
                            sigvnt1d(i1:i2,3)
                       i1=i2+1
                       i2=i2+n_cl
                    end do
                 end do
                 deallocate(zni_2d)
              else          !do not contribute to the enhancement factor from  beam-beam reactions
                 !if nsigexc_use.eq.3 then n*sig*v are calculated in excited state model using adas310_fortran_driver
                 ! and there is  no needs to calculated enhancement factor
                 do i=lcentr,ledge
                    do iel=1,nbnsve(isb0)
                       if(sigvnt_1d(iel,i,isb0).ne.0.0_fp) then
                          bnsvexc(i,iel,isb0)=bms_exc_1d(iel,i,isb0) &
                               /sigvnt_1d(iel,i,isb0)
                       end if
                    end do
                 end do
              end if
           end if            !if(nlbbcx)

        end if               !if(ierr.ne.0) then

     end if

     call nbatom_free_memory

     !if((nsigexc.eq.1).and.ilhydr.and.(lev_nbidep.eq.2))
     !
     if(((nsigexc.eq.1.or.nsigexc.eq.3).and.ilhydr.and. &
          (lev_nbidep.ne.2)).or.((nsigexc.eq.2).and.ilhydr)) then
        imaxn=ng+1+nsbeam
        allocate(zzi(imaxn),zni(imaxn),zti(imaxn))
        do i=lcentr,ledge
           do iel=1,nbnsve(isb0)
              inumi=ng+1
              do ij=1,ng
                 zni(ij)=rhob(i,ij,2)
                 zzi(ij)=backz(ij)
                 zti(ij)=ti(i,2)/aplasm(ij)
              end do
              zni(inumi)=rhi(i,2)
              zzi(inumi)=xzimpj(i)
              zti(inumi)=ti(i,2)/aimpj(i)
              !     fast ions
              inumi0=inumi
              inumi=inumi+nsbeam
              do ij=1,nsbeam
                 ijx=inumi0+ij
                 zni(ijx)=rhbs(i,2,ij)
                 if(nlheavy(ij)) then
                    zzi(ijx)=max(zli,xzbmav(i,ij))
                 else
                    zzi(ijx)=xzbeams(ij)
                 end if
                 zti(ijx)=0.6666667_fp*ebavps(i,ij)/abeams(ij)
              end do
              !
              zeova=bnsves(iel,isb0)
              bnsvexc(i,iel,isb0)=sigfacz(nonlin, &
                   zeova,1.0_fp,rhoel(i,2), &
                   te(i,2),zni,zzi,zti,inumi,zdum1,zdum2)
              !
              !
           end do            ! iel loop
           !
        end do              ! i loop
        if (allocated(zzi)) deallocate(zzi,zni,zti)
     end if
#ifdef __DEBUG
     !         write(nonlin,*) 'no. of energy pts: ',nbnsve(isb0),' isb=',isb0
     !         if((nsigexc_use.gt.0).and.ilhydr) then
     !            allocate(zwk(max(nzones,nbnsve(isb0))))
     !            zwk=0.              ! for debugger
     !            call r8_grf3f1(xi(lcentr:ledge,2),
     !     >         bnsvexc(lcentr:ledge,1:nbnsve(isb0),isb0),
     !     >         bnsves(1:nbnsve(isb0),isb0),zwk,
     !     >         nzones,nbnsve(isb0),nzones,
     !     >         'xi','excitation factor','e/a',
     !     >         ' ',' ','ev/amu','excited states correction',runid,0)
     !            deallocate(zwk)
     !         end if
#endif
     !
     !---------------------------------------------------------------------
     !  beam-beam interaction stuff-- 2d beam dist fcn zone geometry
     !  ** no excited state correction as of mar 1985 **
     !

     zv2eva=1.0_fp/(1.384e6_fp**2)        ! v**2 --> ev/a
     !
     if(.not.nlbbcx.or.(isb0.gt.nsbeam)) goto 900      ! skip if no beam-beam

     do isb=1,nsbeam
        if(nlfprod(isb) .and. .not. nlfatom) goto 800 ! tbt 11/06/91
        if(nlheavy(isb)) goto 800
        !
        izbi=xzbeams(isb)+0.1_fp
        !
        !	 .select energy grid with larger range.
        isbe = ibmgrd(isb,isb0)
        !
        !             beam-beam interactions
        !
        !             nth gen (for nutrav)
        !
        !     .loop over beam zones
        do iz=1,nfbznsi

           !	 .loop over energy grid points
           if(bdns2p(iz,isb).gt.0.0_fp) then
              do iel=1,nbnsve(isbe)
                 zeova = bnsves(iel,isbe)
                 zetran=vbtr2p(iz,isb)**2*zv2eva
                 zerel=zeova+(eb2ppp(iz,isb)+eb2plp(iz,isb))/ &
                      abeams(isb)-zetran
                 zearg(iel)=0.001_fp*zerel
              end do
              !
              if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
                 call hatom_sigv(ii,zearg,nbnsve(isbe),izb,izbi, &
                      zsigv_ii,ierr)
                 call hatom_sigv(icx,zearg,nbnsve(isbe),izb,izbi, &
                      zsigv_cx,ierr)
              elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
                 call adas_sigv(ii,zearg,nbnsve(isbe),izb,izbi, &
                      zsigv_ii,ierr)
                 if(ierr.eq.10) then
                    write(nonlin,*) &
                         '?nbatom->adas data is requested but not found'
                    goto 999
                 end if
                 call adas_sigv(icx,zearg,nbnsve(isbe),izb,izbi, &
                      zsigv_cx,ierr)
                 if(ierr.eq.10) then
                    write(nonlin,*) &
                         '?nbatom->adas data is requested but not found'
                    goto 999
                 end if
              end if

              !
              do iel=1,nbnsve(isbe)
                 bbnsvii(iz,iel,isb0,isb) = &
                      bdns2p(iz,isb)/bmvol(iz)*1.0e6_fp*zsigv_ii(iel)
                 bbnsvcx(iz,iel,isb0,isb) = &
                      bdns2p(iz,isb)/bmvol(iz)*1.0e6_fp*zsigv_cx(iel)
              end do
           else
              do iel=1,nbnsve(isbe)
                 bbnsvii(iz,iel,isb0,isb) = 0.0_fp
                 bbnsvcx(iz,iel,isb0,isb) = 0.0_fp
              end do
           end if                    ! bdns2p


        end do
        !
        !---------------------------------------------------------------------
        !
        !        end of loop over beam ion species (isb)

800     continue
     end do
     !
     !	end of loop over beam neutral species (isb0)

900  continue
  end do
  !
  !-----------------------------------------------------
  !  now prepare tables for subroutine "cxnsum"
  !
  do isb=1,nsbeam  ! loop over beam ion species (isb)
     !
     if(nlfprod(isb) .and. .not. nlfatom) cycle
     if(nlheavy(isb)) cycle
     !
     izb = xzbeams(isb)+0.1_fp
     !
     !  beyond bdy:
     !
     do ig=1,ng
        iz0=backz(ig)+0.1_fp
        !
        do iel=1,nbnsve(isb)
           zearg(iel)=0.001_fp*bnsves(iel,isb)
        end do
        if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
           call hatom_sigv(icx,zearg,nbnsve(isb),iz0,izb, &
                zsigv_cx,ierr)
           if(ierr.ne.0) then
              write(nonlin,*) '?nbatom->hatom code error.'
              goto 999
           end if
        elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
           call adas_sigv(icx,zearg,nbnsve(isb),iz0,izb, &
                zsigv_cx,ierr)
           if(ierr.eq.10) then
              write(nonlin,*) &
                   '?nbatom->adas data is requested but not found'
              goto 999
           end if
        end if
        do iel=1,nbnsve(isb)
           cxn_thcx_a(ig,iel,isb)=dn0out*1.0e6_fp*zsigv_cx(iel)
        end do
     end do
     !
     !  wall neutrals:
     !
     do i=lcentr,ledge
        do ig=1,ng
           iz0=backz(ig)+0.1_fp
           !
           do iel=1,nbnsve(isb)
              zearg(iel)=0.001_fp*bnsves(iel,isb)
              ztarg(iel)=0.001_fp*ewall0(ig,i)/aplasm(ig)
           end do
           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              call hatom_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
              call hatom_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
           elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
              call adas_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              call adas_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if
           do iel=1,nbnsve(isb)
              cxn_thcx_wa(ig,i,iel,isb)= &
                   walln0(ig,i)*1.0e6_fp*zsigv_cx(iel)
              cxn_thii_wa(ig,i,iel,isb)= &
                   walln0(ig,i)*1.0e6_fp*zsigv_ii(iel)
           end do
        end do
     end do
     !
     !  halo neutrals  sigma*v*n for fast ion + thermal halo neutrals, it could be updated
     !  after 3d halo model was called  in depbox
     !  by  calling fast_halo_sigvn(isb)
     do iz=1,nfbznsi
        do ig=1,ng
           iz0=backz(ig)+0.1_fp
           !
           do iel=1,nbnsve(isb)
              zearg(iel)=0.001_fp*bnsves(iel,isb)
              ztarg(iel)=0.001_fp*evol02(ig,iz)/aplasm(ig)
           end do
           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              call hatom_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
              call hatom_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
           elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
              call adas_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              call adas_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if

           do iel=1,nbnsve(isb)
              cxn_thcx_ha(ig,iz,iel,isb)= &
                   voln02(ig,iz)*1.0e6_fp*zsigv_cx(iel)
              cxn_thii_ha(ig,iz,iel,isb)= &
                   voln02(ig,iz)*1.0e6_fp*zsigv_ii(iel)
           end do
        end do
     end do

     !
     ilhydr = (izb.eq.1)            ! hydrogenic isotope?
     !
     !
     ! fast ion intereaction with nth generation neutrals
     !
     !
     do isb0=1,nsbeam !loop over neutrals species
        iz0=xzbeams(isb0)+0.1_fp

        do iz=1,nfbznsi
           do iel=1,nbnsve(0)
              zearg(iel)=0.001_fp*bnsves(iel,0)
              ztarg(iel)=0.001_fp*0.6666667_fp*be0x2p(iz,isb0) &
                   /abeams(isb0) !'temperature' of neutrals species population
           end do
           if( bn0x2p(iz,isb0).le.0.0_fp) cycle

           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              call hatom_btsigv(icx,i_ion,zearg,ztarg,nbnsve(0), &
                   izb,iz0,zsigv_cx,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
              call hatom_btsigv(ii,i_ion,zearg,ztarg,nbnsve(0), &
                   izb,iz0,zsigv_ii,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
           elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
              call adas_btsigv(icx,i_ion,zearg,ztarg,nbnsve(0), &
                   izb,iz0,zsigv_cx,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              call adas_btsigv(ii,i_ion,zearg,ztarg,nbnsve(0), &
                   izb,iz0,zsigv_ii,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if

           do iel=1,nbnsve(0)
              cxn_bbcx_2d(isb0,iz,iel,isb)=1.0e6_fp*zsigv_cx(iel) &
                   *bn0x2p(iz,isb0)/bmvol(iz) ![1/sec]
              cxn_bbii_2d(isb0,iz,iel,isb)=1.0e6_fp*zsigv_ii(iel) &
                   *bn0x2p(iz,isb0)/bmvol(iz) ![1/sec]
           end do
        end do
     end do
     !
  end do !end loop over fast ion species
  !
  !  sigma*v tables for beam-beam calculations in cxnsum
  !
  if(cxn_zmax.gt.0) then
     do izi=cxn_zmin,cxn_zmax
        do iz0=cxn_zmin,cxn_zmax
           do iel=1,nbnsve(0)
              zearg(iel)=0.001_fp*bnsves(iel,0)
           end do
           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              call hatom_sigv(icx,zearg,nbnsve(0),iz0,izi, &
                   zsigv_cx,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
              call hatom_sigv(ii,zearg,nbnsve(0),iz0,izi, &
                   zsigv_ii,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
           elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
              call adas_sigv(icx,zearg,nbnsve(0),iz0,izi, &
                   zsigv_cx,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              call adas_sigv(ii,zearg,nbnsve(0),iz0,izi, &
                   zsigv_ii,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if

           do iel=1,nbnsve(0)
              cxn_bbcx(iz0,iel,izi)=zsigv_cx(iel)*1.0e6_fp !(cm**3)/sec
              cxn_bbii(iz0,iel,izi)=zsigv_ii(iel)*1.0e6_fp !(cm**3)/sec
           end do
        end do
     end do
  end if
  !-------------------------------------------
  !  beam target fusion tables
  !   ...for use in subroutine btfusn
  !
  do i=lcentr,ledge
     !
     ztix=tmjsm(i)*.001_fp
     !
     do iel=1,nbnsve(0)
        zearg(iel)=0.001_fp*bnsves(iel,0)
        ztarg(iel)=ztix
     end do
     !
     !  do the fusion reactions... cf ebatom.for allocated only arrays for
     !  reactions that are needed...
     !
     if(allocated(btfus_ddp)) then
        call fusn_btsigv_ddp(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_ddp(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_ddn)) then
        call fusn_btsigv_ddn(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_ddn(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_dt)) then
        call fusn_btsigv_dt(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_dt(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_d3)) then
        call fusn_btsigv_d3(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_d3(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_tt)) then
        if(lev_ttnuk.eq.0) then
           call fusn_btsigv_tt(zearg,ztarg,zsigv_ii,nbnsve(0))
        elseif(lev_ttnuk.eq.1) then
           call fusn_btsigv_dpctt(zearg,ztarg,zsigv_ii,nbnsve(0))
        end if
        do iel=1,nbnsve(0)
           btfus_tt(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_td)) then
        call fusn_btsigv_td(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_td(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
     if(allocated(btfus_3d)) then
        call fusn_btsigv_3d(zearg,ztarg,zsigv_ii,nbnsve(0))
        do iel=1,nbnsve(0)
           btfus_3d(i,iel)=1.0e6_fp*zsigv_ii(iel)
        end do
     end if
     !
  end do
  !
999 continue
  deallocate(zearg,ztarg,zsigv_cx,zsigv_ii, &
       stat=istat)
  !
  return

contains
  !********************
  subroutine alloc_wrk
    ! allocates an working arrays

    allocate(te_wrk(i_nzone_min:i_nzone_max), &
         ne_wrk(i_nzone_min:i_nzone_max), &
         sigvnt1d(i_nzone_min:i_nzone_max,3))
    ne_wrk=0.0_fp
    te_wrk=0.0_fp
    sigvnt1d=0.0_fp

    allocate(zni_2d(1:imaxn,i_nzone_min:i_nzone_cl))
    zni_2d=0.0_fp

    if(allocated(zt_i)) deallocate(zt_i)
    allocate(zt_i(1:ijx,i_nzone_min:i_nzone_max))
    zt_i=0.0_fp

  end subroutine  alloc_wrk
  !
  !-----------------------------
  !
  subroutine corresp1d_2d
    integer :: i,l_edge
    ! write the correspondence between 2d and 1d grid indecies
    l_edge=min(ledge,size(nthzsm))
    l_edge=max(lcentr,l_edge)
    do i=lcentr,l_edge
       if((nthzsm(i-1)+1.le.0).or.(nthzsm(i).gt.nfbznsi)) cycle
       indx1d_2d(nthzsm(i-1)+1:nthzsm(i))=i
    end do

  end subroutine corresp1d_2d
  !
  !-----------------------------
  !
  subroutine contribute_step0(irr)
    integer :: irr,kk,iij,n_fast_in,i3ng
    integer :: ix, ij, ij1
    real(fp) :: bms_tot, xzimpx_s
    logical :: il_excited
    irr=0

    il_excited=(nsigexc.eq.3)

    if(allocated(zzi)) deallocate(zzi)
    allocate(zzi(imaxn))
    if(allocated(zn_i)) deallocate(zn_i)
    allocate(zn_i(imaxn,mj))
    if(allocated(zt_i)) deallocate(zt_i)
    allocate(zt_i(imaxn,mj))
    if(allocated(zza)) deallocate(zza) !atomic mass
    allocate(zza(imaxn))
    !      allocate(enhanc_fact(nbnsve(isb0),lcentr:ledge,
    !     >     lcentr:ledge))
    if(allocated(zearg_adas)) deallocate(zearg_adas)
    allocate(zearg_adas(nbnsvmx,imaxn))
    zn_i=0.0_fp
    zt_i=0.0_fp
    zzi=0.0_fp
    zza=0.0_fp
    zearg_adas=0.0_fp
    !      enhanc_fact=1.0_fp
    do kk=1,nbnsve(isb0)
       zearg_adas(kk,1:imaxn)=0.001_fp*bnsves(kk,isb0) !energy per nucleon [kev/amu]
    end do
    !     background plasma
    do ij=1,ng
       zn_i(ij,:)=rhob(:,ij,2)
       zzi(ij)=backz(ij)
       zza(ij)=aplasm(ij)
       zt_i(ij,:)=ti(:,2)*0.001_fp/aplasm(ij)
    end do
    !     impurities
    ix=ng+1
    zn_i(ng+1:ng+nrhix_s,:)=0.0_fp
    xzimpx_s=0.0_fp
    do ij=1,nrhix
       ij1=min(ij+1,nrhix)
       zn_i(ix,:)=zn_i(ix,:)+rhix(:,2,ij) !summed impurity density for all charge state
       xzimpx_s=max(xzimpx_s,xzimpx(ij)) !max charge state (corresponds fully stripped ion)
       if(aimpx(ij).eq.aimpx(ij1))then
          if(ij.ne.ij1) then
             cycle
          end if
       end if
       zzi(ix)=xzimpx_s
       zza(ix)=aimpx(ij)
       zt_i(ix,:)=ti(:,2)*0.001_fp /aimpx(ij) !converted to [kev/amu]
       xzimpx_s=0.0_fp
       ix=ix+1
    end do
    !     fast ions
    iij=0
    do ij=1,nsbeam
       iij=iij+1
       ijx=inumi0+iij         !total # of plasma species as a target (thermal+impurities+beam ion)
       zzi(ijx)=xzbeams(ij)
       zn_i(ijx,lcentr:ledge)=rhbs(lcentr:ledge,2,ij) ! [1/cm**3]
       zza(ijx)=abeams(ij)
       zt_i(ijx,lcentr:ledge)=0.6666667_fp*ebavps(lcentr:ledge,ij) &
            /abeams(ij)*0.001_fp !converted to [kev/amu]
    end do
    n_fast_in=n_fast
    !if(nlbbcx) n_fast_in=0 !include contributions from fast prtcls
    !to compute enhancement factor
    !
    z_beam=real(izb,fp)
    !     compute enhancement factor, sig*v*n for excited state and ground state
    !
    !     get enhancement factor <sig*v>*n for ground state and for excited state without
    !     beam-beam reactions, only density for beam ion species is included for quasineutrality
    !
    if(il_excited)then
       allocate(bmsexc2d(1:nbnsve(isb0),lcentr:ledge,3+ng))
       i3ng=3+ng
    else
       allocate(bmsexc2d(1:nbnsve(isb0),lcentr:ledge,3))
       i3ng=3
    end if
    bmsexc2d=0.0_fp

    allocate(sigvnt2d(1:nbnsve(isb0),lcentr:ledge,3))
    sigvnt2d=0.0_fp

    call sigfacz_adas(zearg_adas(1:nbnsve(isb0),1:ijx),nbnsve(isb0), &
         z_beam,a_beam,rhoel(lcentr:ledge,2), &
         nzones,te(lcentr:ledge,2)*0.001_fp,zn_i(1:ijx, &
         lcentr:ledge), &
         zzi,zza,zt_i(1:ijx,lcentr:ledge),ijx, &
         n_fast_in,i3ng, &
         bmsexc2d(1:nbnsve(isb0),lcentr:ledge,1:i3ng), &
         sigvnt2d(1:nbnsve(isb0),lcentr:ledge,1:3),il_excited,irr)

    do i=lcentr,ledge
       do iel=1,nbnsve(isb0)
          if(il_excited) then
             bms_tot=sum(bmsexc2d(iel,i,:))
             if(bms_tot.le.0.0_fp) cycle
             bnsvtot(i,iel,isb0)=bms_tot
             bnsvief(i,iel,isb0)=bmsexc2d(iel,i,1)/bms_tot
             bnsviif(i,iel,isb0)=bmsexc2d(iel,i,2)/bms_tot
             bnsvizf(i,iel,isb0)=bmsexc2d(iel,i,3)/bms_tot
             bnsvcxf(i,1:ng,iel,isb0)=bmsexc2d(iel,i,4:i3ng)/bms_tot
             bnsvexc(i,iel,isb0)=1.0_fp
          else
             bms_exc_1d(iel,i,isb0)=bmsexc2d(iel,i,1) !total bms
             sigvnt_1d(iel,i,isb0)=sigvnt2d(iel,i,1)+sigvnt2d(iel,i,2) &
                  +sigvnt2d(iel,i,3) !total sigvnt
          end if
       end do
    end do
    deallocate(bmsexc2d,sigvnt2d,zearg_adas)

  end subroutine contribute_step0
  !
  !-----------------------------
  !
  subroutine contribute_thermal(irr)
    integer :: irr,kk,iij,ino,n_fast_in,i3ng
    integer :: ix, ij, ij1
    real(fp) :: bms_tot,xzimpx_s
    logical :: il_excited
    irr=0

    il_excited=(nsigexc.eq.3)

    if(allocated(zzi)) deallocate(zzi)
    allocate(zzi(imaxn))
    if(allocated(zn_i)) deallocate(zn_i)
    allocate(zn_i(imaxn,mj))
    if(allocated(zt_i)) deallocate(zt_i)
    allocate(zt_i(imaxn,mj))
    if(allocated(zza)) deallocate(zza) !atomic mass
    allocate(zza(imaxn))
    !      allocate(enhanc_fact(nbnsve(isb0),lcentr:ledge,
    !     >     lcentr:ledge))
    if(allocated(zearg_adas)) deallocate(zearg_adas)
    allocate(zearg_adas(nbnsvmx,imaxn))
    zn_i=0.0_fp
    zt_i=0.0_fp
    zzi=0.0_fp
    zza=0.0_fp
    zearg_adas=0.0_fp
    !      enhanc_fact=1.0_fp
    do kk=1,nbnsve(isb0)
       zearg_adas(kk,1:imaxn)=0.001_fp*bnsves(kk,isb0) !energy per nucleon [kev/amu]
    end do
    !     background plasma
    do ij=1,ng
       zn_i(ij,:)=rhob(:,ij,2)
       zzi(ij)=backz(ij)
       zza(ij)=aplasm(ij)
       zt_i(ij,:)=ti(:,2)*0.001_fp/aplasm(ij)
    end do
    !     impurities
    ix=ng+1
    zn_i(ng+1:ng+nrhix_s,:)=0.0_fp
    xzimpx_s=0.0_fp

    do ij=1,nrhix
       ij1=min(ij+1,nrhix)
       zn_i(ix,:)=zn_i(ix,:)+rhix(:,2,ij) !summed impurity density for all charge state
       xzimpx_s=max(xzimpx_s,xzimpx(ij)) !max charge state (corresponds fully stripped ion)
       if(aimpx(ij).eq.aimpx(ij1))then
          if(ij.ne.ij1) then
             cycle
          end if
       end if
       zzi(ix)=xzimpx_s
       zza(ix)=aimpx(ij)
       zt_i(ix,:)=ti(:,2)*0.001_fp /aimpx(ij) !converted to [kev/amu]
       xzimpx_s=0.0_fp
       ix=ix+1
    end do
    !     fast ions, beam ptcls
    iij=0
    do ij=1,nsbeam
       if(nlfprod(ij)) cycle
       do ino=1,2
          iij=iij+1
          ijx=inumi0+iij      !total # of plasma species as a target (thermal+impurities+beam ion)
          zzi(ijx)=xzbeams(ij)
          zza(ijx)=abeams(ij)
          zn_i(ijx,lcentr:ledge)=bnb_ienx(lcentr:ledge,ino,ij)* &
               1.e-6_fp       !to[1/cm**3]
       end do
    end do
    !     fast ions fusion product
    do ij=1,nsbeam
       if(nlfprod(ij)) then
          ijx=ijx+1
          zzi(ijx)=xzbeams(ij)
          zza(ijx)=abeams(ij)
          zn_i(ijx,lcentr:ledge)=bnb_fnx(lcentr:ledge,ij)* &
               1.e-6_fp       !to[1/cm**3]
       end if
    end do

    !
    n_fast_in=n_fast !contribution from thermal + impurities only

    z_beam=real(izb,fp)
    !     compute enhancement factor, sig*v*n for excited state and ground state
    !
    !     get enhancement factor <sig*v>*n for ground state and for excited state without
    !     beam-beam reactions, only density for beam ion species is included for quasineutrality
    !
    if(il_excited)then
       allocate(bmsexc2d(1:nbnsve(isb0),lcentr:ledge,3+ng))
       i3ng=3+ng
    else
       allocate(bmsexc2d(1:nbnsve(isb0),lcentr:ledge,3))
       i3ng=3
    end if
    bmsexc2d=0.0_fp

    allocate(sigvnt2d(1:nbnsve(isb0),lcentr:ledge,3))
    sigvnt2d=0.0_fp

    call sigfacz_adas(zearg_adas(1:nbnsve(isb0),1:ijx),nbnsve(isb0), &
         z_beam,a_beam,rhoel(lcentr:ledge,2), &
         nzones,te(lcentr:ledge,2)*0.001_fp,zn_i(1:ijx, &
         lcentr:ledge), &
         zzi,zza,zt_i(1:ijx,lcentr:ledge),ijx, &
         n_fast_in,i3ng,&
    !     >     enhanc_fact(1:nbnsve(isb0),lcentr:ledge,
    !     >     lcentr:ledge), &
         bmsexc2d(1:nbnsve(isb0),lcentr:ledge,1:i3ng), &
         sigvnt2d(1:nbnsve(isb0),lcentr:ledge,1:3),il_excited,irr)

    do i=lcentr,ledge
       do iel=1,nbnsve(isb0)
          if(il_excited) then
             bms_tot=sum(bmsexc2d(iel,i,:))
             if(bms_tot.le.0.0_fp) cycle
             bnsvtot(i,iel,isb0)=bms_tot
             bnsvief(i,iel,isb0)=bmsexc2d(iel,i,1)/bms_tot
             bnsviif(i,iel,isb0)=bmsexc2d(iel,i,2)/bms_tot
             bnsvizf(i,iel,isb0)=bmsexc2d(iel,i,3)/bms_tot
             bnsvcxf(i,1:ng,iel,isb0)=bmsexc2d(iel,i,4:i3ng)/bms_tot
             bnsvexc(i,iel,isb0)=1.0_fp
          else
             bms_exc_1d(iel,i,isb0)=bmsexc2d(iel,i,1) !total bms
             sigvnt_1d(iel,i,isb0)=sigvnt2d(iel,i,1)+sigvnt2d(iel,i,2) &
                  +sigvnt2d(iel,i,3) !total sigvn
          end if
       end do
    end do
    deallocate(bmsexc2d,sigvnt2d,zearg_adas)

  end subroutine contribute_thermal
  !
  !-----------------------------
  !
  subroutine contribute_fast
    integer :: kk,iij,ino,ijx,i1,i2,ib0,ien0,n_cl,i3ng
    logical :: il_excited

    il_excited=(nsigexc.eq.3)
    !for nsigexc.eq.3 - excited state model will be used
    if(allocated(zearg_adas)) deallocate(zearg_adas)
    allocate(zearg_adas(i_nzone_min:i_nzone_max,imaxn))
    zearg_adas=0.0_fp
    if(allocated(zn_i)) deallocate(zn_i)
    allocate(zn_i(imaxn,i_nzone_min:i_nzone_max))
    iij=0
    do ij=1,nsbeam
       if(nlfprod(ij)) cycle
       do ino=1,2
          iij=iij+1
          ijx=inumi0+iij
          i1=i_nzone_min
          i2=i_nzone_cl
          n_cl=i_nzone_cl-i_nzone_min+1
          do ib0=1,nbeam
             if(lsbeama(ib0).ne.isb0) cycle
             a_rds=abeama(ib0)*abeams(ij)/(abeama(ib0)+abeams(ij))
             do ien0=1,3
                zn_i(ijx,i1:i2)=zni_2d(ijx,i_nzone_min:i_nzone_cl)
                ne_wrk(i1:i2)=ne_wrk(i_nzone_min:i_nzone_cl)
                te_wrk(i1:i2)=te_wrk(i_nzone_min:i_nzone_cl)
                if(ino.eq.1)then
                   if(nbbcx_avg.eq.0) then
                      zt_i(ijx,i1:i2)= &
                           rms_ebb_co(i_nzone_min:i_nzone_cl, &
                           ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                      zearg_adas(i1:i2,ijx)=ebb_rel_co &
                           (i_nzone_min:i_nzone_cl,ien0,ib0,ij)/ &
                           a_rds*0.001_fp !converted to [kev/amu]
                   else if(nbbcx_avg.gt.0) then
                      zt_i(ijx,i1:i2)= &
                           rms_ebbx_co(i_nzone_min:i_nzone_cl, &
                           ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                      zearg_adas(i1:i2,ijx)=ebb_relx_co &
                           (i_nzone_min:i_nzone_cl,ien0,ib0,ij)/ &
                           a_rds*0.001_fp !converted to [kev/amu]
                   end if
                else
                   if(nbbcx_avg.eq.0) then
                      zt_i(ijx,i1:i2)= &
                           rms_ebb_cnt(i_nzone_min:i_nzone_cl, &
                           ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                      zearg_adas(i1:i2,ijx)=ebb_rel_cnt &
                           (i_nzone_min:i_nzone_cl,ien0,ib0,ij)/ &
                           a_rds*0.001_fp !converted to [kev/amu]
                   else if(nbbcx_avg.gt.0) then
                      zt_i(ijx,i1:i2)= &
                           rms_ebbx_cnt(i_nzone_min:i_nzone_cl, &
                           ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                      zearg_adas(i1:i2,ijx)=ebb_relx_cnt &
                           (i_nzone_min:i_nzone_cl,ien0,ib0,ij)/ &
                           a_rds*0.001_fp !converted to [kev/amu]
                   end if
                end if
                i1=i2+1
                i2=i2+n_cl
             end do            !ien0
          end do               !ib0
       end do
    end do
    ! fast ions fusion product
    do ij=1,nsbeam
       if(nlfprod(ij)) then
          ijx=ijx+1
          i1=i_nzone_min
          i2=i_nzone_cl
          n_cl=i_nzone_cl-i_nzone_min+1
          do ib0=1,nbeam
             if(lsbeama(ib0).ne.isb0) cycle
             a_rds=abeama(ib0)*abeams(ij)/(abeama(ib0)+abeams(ij))
             do ien0=1,3
                ne_wrk(i1:i2)=ne_wrk(i_nzone_min:i_nzone_cl)
                te_wrk(i1:i2)=te_wrk(i_nzone_min:i_nzone_cl)
                zn_i(ijx,i1:i2)=zni_2d(ijx,i_nzone_min:i_nzone_cl)
                if(nbbcx_avg.eq.0) then
                   zt_i(ijx,i1:i2)= &
                        rms_ebb_fn(i_nzone_min:i_nzone_cl, &
                        ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                   zearg_adas(i1:i2,ijx)=ebb_rel_fn &
                        (i_nzone_min:i_nzone_cl,ien0,ib0,ij) &
                        /a_rds*0.001_fp !converted to [kev/amu]
                else if(nbbcx_avg.gt.0) then
                   zt_i(ijx,i1:i2)= &
                        rms_ebbx_fn(i_nzone_min:i_nzone_cl, &
                        ien0,ib0,ij)/a_rds*0.001_fp !converted to [kev/amu]
                   zearg_adas(i1:i2,ijx)=ebb_relx_fn &
                        (i_nzone_min:i_nzone_cl,ien0,ib0,ij) &
                        /a_rds*0.001_fp !converted to [kev/amu]
                end if
                i1=i2+1
                i2=i2+n_cl
             end do            !ien0
          end do               !ib0
       end if

    end do

    do ij=1,inumi0
       i1=i_nzone_min
       i2=i_nzone_cl
       n_cl=i_nzone_cl-i_nzone_min+1
       do ib0=1,nbeam
          if(lsbeama(ib0).ne.isb0) cycle
          do ien0=1,3
             ne_wrk(i1:i2)=ne_wrk(i_nzone_min:i_nzone_cl)
             te_wrk(i1:i2)=te_wrk(i_nzone_min:i_nzone_cl)
             zn_i(ij,i1:i2)=zni_2d(ij,i_nzone_min:i_nzone_cl)
             i1=i2+1
             i2=i2+n_cl
          end do               !ien0
       end do                  !ib0
    end do
    i3ng=3+n_fast
    allocate(bmsexc1d(i_nzone_min:i_nzone_max,i3ng))
    bmsexc1d=0.0_fp
    z_beam=real(izb,fp)

    call sigfacz_adas( &
         zearg_adas(i_nzone_min:i_nzone_max,1:ijx),i_zones, &
         z_beam,a_beam,ne_wrk(i_nzone_min:i_nzone_max), &
         te_wrk(i_nzone_min:i_nzone_max), &
         zn_i(1:ijx,i_nzone_min:i_nzone_max), &
         zzi(1:ijx),zza(1:ijx),zt_i(1:ijx, &
         i_nzone_min:i_nzone_max), &
         ijx,-n_fast,i3ng, &
         bmsexc1d(i_nzone_min:i_nzone_max,1:i3ng), &
         sigvnt1d(i_nzone_min:i_nzone_max,1:3),il_excited, &
         ierr)
    deallocate(zearg_adas)

  end subroutine contribute_fast

  subroutine nbatom_free_memory

    if (allocated(zzi)) deallocate(zzi)
    if (allocated(zt_i)) deallocate(zt_i)
    if (allocated(zn_i)) deallocate(zn_i)
    if (allocated(te_wrk)) deallocate(te_wrk)
    if (allocated(ne_wrk)) deallocate(ne_wrk)
    if (allocated(bmsexc1d)) deallocate(bmsexc1d)
    if (allocated(sigvnt1d)) deallocate(sigvnt1d)

  end subroutine nbatom_free_memory

end subroutine nbatom

subroutine fast_halo_sigvn
  !     halo neutrals  sigma*v*n for fast ion + thermal halo neutrals, 
  !     it could be updated in depbox
  use iso_c_binding, only: fp => c_double
  use adas_mod, only : adas_btsigv
  use nbatom_mod, only : nbnsve,bnsves,nbnsvmx
  use nbatom_mod, only : cxn_thii_ha,cxn_thcx_ha
  use nbi_com, only : nfbznsi,ng, backz,aplasm,voln02,evol02,xzbeams
  use nbi_com, only : lev_nbidep,nsbeam,nonlin
  !
  !============
  implicit none
  integer, parameter :: icx=1, ii=2
  integer, parameter :: i_neutral=1, i_ion=2
  integer :: iz, ig, izb, iz0,iel,isb,ierr
  real(fp), dimension(:), allocatable :: zearg,ztarg,zsigv_cx,zsigv_ii

  if(allocated(zearg)) deallocate(zearg)
  allocate(zearg(nbnsvmx))
  if(allocated(ztarg)) deallocate(ztarg)
  allocate(ztarg(nbnsvmx))
  if(allocated(zsigv_cx)) deallocate(zsigv_cx)
  if(allocated(zsigv_ii)) deallocate(zsigv_ii)
  allocate(zsigv_cx(nbnsvmx),zsigv_ii(nbnsvmx))
  zsigv_cx=0.0_fp
  zsigv_ii=0.0_fp
  zearg=0.0_fp
  ztarg=0.0_fp

  do isb=1,nsbeam
     izb = xzbeams(isb)+0.1_fp
     do iz=1,nfbznsi
        do ig=1,ng
           iz0=backz(ig)+0.1_fp
           !
           do iel=1,nbnsve(isb)
              zearg(iel)=0.001_fp*bnsves(iel,isb)
              ztarg(iel)=0.001_fp*evol02(ig,iz)/aplasm(ig)
           end do
           if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics
              call hatom_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
              call hatom_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.ne.0) then
                 write(nonlin,*) '?nbatom->hatom error.'
                 goto 999
              end if
           elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
              call adas_btsigv(icx,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_cx,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
              call adas_btsigv(ii,i_ion,zearg,ztarg, &
                   nbnsve(isb),izb,iz0,zsigv_ii,ierr)
              if(ierr.eq.10) then
                 write(nonlin,*) &
                      '?nbatom->adas data is requested but not found'
                 goto 999
              end if
           end if

           do iel=1,nbnsve(isb)
              cxn_thcx_ha(ig,iz,iel,isb)= &
                   voln02(ig,iz)*1.0e6_fp*zsigv_cx(iel)
              cxn_thii_ha(ig,iz,iel,isb)= &
                   voln02(ig,iz)*1.0e6_fp*zsigv_ii(iel)
           end do
        end do
     end do
  end do
999 if(allocated(zearg)) deallocate(zearg)
  if(allocated(ztarg)) deallocate(ztarg)
  if(allocated(zsigv_cx)) deallocate(zsigv_cx)
  if(allocated(zsigv_ii)) deallocate(zsigv_ii)

  return

end subroutine fast_halo_sigvn

!******************** end file nbatom.fnb ; group nbatom ***************
! 11sep2004 fgtok
