module adas_sigfacz
!  use iso_c_binding, only: fp => c_double
  use physconst_mod, only : kev2ev=>zkev2ev, fp
  use adas_mod, only : adas_bms
  use adas_getsigs
  
  public :: sigfacz_adas
  interface sigfacz_adas
     module procedure  sigfacz_adas3d, sigfacz_adas1d, sigfacz_adas2d
  end interface
!
! subroutine sigfacz_adas computes excited states stopping enhancement 
! factor and total sig*v*n for excited state and ground state
!
!  
! enhancement factor can be computed as:
!                    enhanc_fact=bms_exc_tot/sigvnt
!
! formula for the free path time (t) in the ground state model is
! 1/t = <sig*v>_ei * n_e  +  sum_i (<sig*v>_cx+<sig*v>_ii)_i * n_i,
! 
! n_e - electron density,
! n_i - impurity density,
! < ... > - maxwellian averaged,
! <sig*v> - is a function of (eb,t_i)
! eb - beam energy,
! t_i - impurity temperature,
! ei- electron impact ionization,
! cx - charge exchange (non-neutralizing  and neutralizing when possible),
! ii - impact ionization.
! 
! in the excited states model (from  ref.1  p.794)
! 
! 1/t = s(eb,n_e,t_i) * sum_i(z_i*n_i),  where
! 
! z_i - impurity charge,
! s(eb,n_e,t_i) = sum_i(z_i * n_i *s_i(eb,n_equv,t_i)) / sum_i(z_i*n_i).
! 
! if all plasma species are counted then sum_i(z_i*n_i) = n_e, and
! 
! 1/t = s(eb,n_e,t_i) * n_e
! 
! 
! 
! enhancement factor is defined as
!  
! s(eb,n_e,t_i) * n_e / {<sig*v>_ei * n_e  +  sum_i (<sig*v>_cx+<sig*v>_ii)_i * n_i,}
! 
! 
! 1. h. anderson et al. (2000) "neutral beam stopping and emission in fusion plasmas i",
! plasma physics and controlled fusion vol. 42, pp 781-806
!
!

contains




subroutine sigfacz_adas1d(erel_vec,ie_vec,z_beam,a_beam,zne,zte,zni,zzi,&
     &                     zza,ztia_vec,inumi,i_nsbeam,ing, &
     &                     bms_exc_tot,sigvnt,il_excited,ier)
  implicit none
  integer, intent(in) :: ie_vec !vector size
  real(fp),intent(in) :: erel_vec(ie_vec,inumi) !interaction energy [kev/amu] beam -> target
  real(fp),intent(in) :: z_beam !beam charge
  real(fp),intent(in) :: a_beam !beam amu
  real(fp),intent(in) :: zne(ie_vec) !electron density [1/cm3]
  real(fp),intent(in) :: zte(ie_vec) !electron temperature [kev]
  real(fp),intent(in) :: zza(inumi) ! ion species atomic mass
  real(fp),intent(in) :: ztia_vec(inumi,ie_vec) !ion species ti/ai temperature [kev/amu]
  real(fp),intent(in) :: zni(inumi,ie_vec) !ion species density [1/cm3]
  integer,intent(in) :: inumi ! number of distinct ion species
  integer,intent(in) :: i_nsbeam ! number of fast ion species, if positive or zero do not include them
                                 ! to calculation of the <sig*v>*n for ground state and
                                 ! bms for excited state
                                 ! if negative, then fast ion species will be taken into
                                 ! account
  integer,intent(in) :: ing ! number of h isotopes ng + 3
  real(fp),intent(in) :: zzi(inumi) ! ion species charge
  real(fp),intent(out) :: sigvnt(ie_vec,3) ! ground state sig*v*n : 1-- sigvne, 2 -- sigviin3, -- sigvcxn 
  real(fp),intent(out) :: bms_exc_tot(ie_vec,ing) 
                                              ! il_excited=.false. adas310 bms data: 1 - total excited state  sig*v*n 
                                              ! il_excited=.true. - adas310_fortran_driver:
                                              ! excited state model sig*v*n  1 -- ie ,2 -- ii, 3 -- ii for impuritie
                                              ! cx for ng h isotopes
  logical,intent(in) :: il_excited !.true. -- adas excited state model will be used to calculate sig*v*n 
                                   !.false.  -- adas beam stopping coeffitient will be used to
                                   ! estimate enhancement factor
  real(fp),dimension(:,:),allocatable :: conc_imp !target species fraction   z_imp * n_imp/sum(z_imp*n_imp)
  real(fp),dimension(:),allocatable :: sum_conc !sum(z_imp*n_imp)
  real(fp),dimension(:), allocatable :: sigv
  real(fp),dimension(:), allocatable :: zne_equiv  !n_e*sum(z_imp*z_imp*n_imp)/sum(z_imp * n_imp)/z_imp

  integer ii,ij,ier,ierr,inumz,iel,inumz1,&
       & inumz_min,inumz_max
  logical :: incl_ei !include ei reaction into ground state model if it is not beam-beam interactions
  !
  !ground state stopping
  inumz1=inumi-abs(i_nsbeam )
  if(i_nsbeam.ge.0) then !only thermal + impurities
     inumz_min=1
     inumz_max=inumz1  
     inumz=inumz1
  else !only only beam ion species
     inumz_min=inumz1+1
     inumz_max=inumi
     inumz=abs(i_nsbeam )
  endif

  !excited state stopping
  bms_exc_tot=0.0_fp
  !ground state stopping
  sigvnt=0.0_fp
  if(il_excited) then
     call getsigs_adas310_fast(erel_vec(1:ie_vec,inumz_min:inumz_max)*kev2ev,ie_vec,&
          & z_beam,a_beam,zte*kev2ev,&
          & ztia_vec(inumz_min:inumz_max,1:ie_vec)*kev2ev,zne,&
          ! t_i/amu [ev]      
          & zni(inumz_min:inumz_max,1:ie_vec),&
          & zzi(inumz_min:inumz_max),zza(inumz_min:inumz_max),&
          & inumz,bms_exc_tot(1:ie_vec,1:ing),ing,ier)
     return
  endif

  allocate(conc_imp(inumi,ie_vec),sum_conc(ie_vec))
  conc_imp=0.0_fp
  allocate(zne_equiv(ie_vec))
  zne_equiv=0._fp
  !
  ! target species fraction is z_imp * n_imp/sum(z_imp*n_imp)
  !
  !calculate  z_i*n_i 
  do ii=1,inumi !thermal + impurities and beam ion densities for quasi-neutrality
     conc_imp(ii,:)=zzi(ii)*zni(ii,:)
  enddo
  
  sum_conc(:)=sum(conc_imp(1:inumi,:),dim=1) !sum_i(z_i*n_i)
  
  do ij=1,ie_vec
     if( sum_conc(ij).eq.0.0_fp) then
        conc_imp(:,ij)=0.0_fp
     else
        do ii=1,inumi
           conc_imp(ii,ij)=conc_imp(ii,ij)/sum_conc(ij)  !z_i*n_i/sum_i(z_i*n_i)
        enddo
     endif
  enddo
  !
  !calculate electron equivalent density sum_i(z_i*z_i*n_i)
  do ii=1,inumi !thermal + impurities and beam ion densities
     zne_equiv(:)=zne_equiv(:)+zzi(ii)*zzi(ii)*zni(ii,:) !sum_i(z_i**2*n_i)
  enddo
  
  do ij=1,ie_vec
     if( sum_conc(ij).eq.0.0_fp) then
        zne_equiv(ij)=0.0_fp
     else
        if(i_nsbeam.ge.0) then
           zne_equiv(ij)=zne_equiv(ij)*zne(ij)/sum_conc(ij) !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)
           ! else assume that n_e = sum_i(z_i*n_i)
        endif
     endif
  enddo
  !
  ! ground state <sig*v>*n
  
  incl_ei=(i_nsbeam.ge.0)
  
  call getsigs_adas(erel_vec(1:ie_vec,inumz_min:inumz_max),ie_vec,z_beam,zte,zne,inumz,&
       & zzi(inumz_min:inumz_max),&
       & zni(inumz_min:inumz_max,1:ie_vec),&
       & ztia_vec(inumz_min:inumz_max,1:ie_vec),&
       & sigvnt,incl_ei,ier)
  
  !excited state beam stopping coefficient
  
  bms_exc_tot=0.0_fp
  
  allocate(sigv(ie_vec))
  sigv=0._fp
  
  
  do ii=inumz_min,inumz_max
     sigv=0._fp
     zne_equiv(:)=zne_equiv(:)/zzi(ii) !equivalent electron density for each impurity
     !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)/z_i
     
     call adas_bms(erel_vec(:,ii),ie_vec,z_beam,conc_imp(ii,:),&
             & ztia_vec(ii,:)*zza(ii),&
!mg 12/05/11 i found that bms tables were created with timp [kev] not [kev/amu] as described in fiels
          & zzi(ii:ii),zne_equiv,sigv,ierr) !only diagonal elements will be computed
     
     bms_exc_tot(:,1)=bms_exc_tot(:,1) + sigv(:) * sum_conc(:)*1.e+6 !<sig*v>*n_e - includes excited states
     !if all plasma species are counted
     ! then sum_conc = zne
     
     zne_equiv(:)=zne_equiv(:)*zzi(ii) !resume for next impurity
  enddo
  
  
  deallocate(conc_imp,sum_conc,sigv,zne_equiv)
end subroutine sigfacz_adas1d


subroutine sigfacz_adas2d(erel_vec,ie_vec,z_beam,a_beam,zne,it_vec,zte,zni,zzi,&
  &                     zza,ztia_vec,inumi,i_nsbeam,ing, &
  &                     bms_exc_tot,sigvnt,il_excited,ier)
  implicit none
  integer, intent(in) :: ie_vec,it_vec !vector size
  real(fp),intent(in) :: erel_vec(ie_vec,inumi) !interaction energy [kev/amu] beam -> target
  real(fp),intent(in) :: z_beam !beam charge
  real(fp),intent(in) :: a_beam !beam amu
  real(fp),intent(in) :: zne(it_vec) !electron density [1/cm3]
  real(fp),intent(in) :: zte(it_vec) !electron temperature [kev]
  real(fp),intent(in) :: ztia_vec(inumi,it_vec) !ion species ti/ai temperature [kev/amu]
  real(fp),intent(in) :: zni(inumi,it_vec) !ion species density [1/cm3]
  integer,intent(in) :: inumi ! number of distinct ion species
  real(fp),intent(in) :: zza(inumi) ! ion species atomic mass
  integer,intent(in) :: i_nsbeam ! number of fast ion species, if positive or zero do not include them
                                 ! to calculation of the <sig*v>*n for ground state and
                                 ! bms for excited state
                                 ! if negative, then fast ion species will be taken into
                                 ! account
  real(fp),intent(in) :: zzi(inumi) ! ion species charge
  logical,intent(in) :: il_excited !.true. -- adas excited state model will be used to calculate  sig*v*n 
                                   !.false.  -- adas beam stopping coeffitient will be used to
                                   ! estimate enhancement factor
  integer,intent(in) :: ing ! number of h isotopes ng + 3
  real(fp),intent(out) :: sigvnt(ie_vec,it_vec,3) !  ground state sig*v*n : 1 -- sigvne,2 -- sigviin,
                                                !                         3 -- sigvcxn 
  real(fp),intent(out) :: bms_exc_tot(ie_vec,it_vec,ing) 
                                              ! il_excited=.false. adas310 bms data: 1 - total excited state  sig*v*n 
                                              ! il_excited=.true. - adas310_fortran_driver:
                                              ! excited state model sig*v*n  1 -- ie ,2 -- ii, 3 -- ii for impuritie
                                              ! cx for ng h isotopes!
  real(fp),dimension(:,:),allocatable :: conc_imp !target species fraction   z_imp * n_imp/sum(z_imp*n_imp)
  real(fp),dimension(:),allocatable :: sum_conc !sum(z_imp*n_imp)
  real(fp),dimension(:,:), allocatable :: sigv
  real(fp),dimension(:), allocatable :: zne_equiv  !n_e*sum(z_imp*z_imp*n_imp)/sum(z_imp * n_imp)/z_imp
  real(fp),dimension(:), allocatable ::zt_i !ion temperature
  integer ii,ij,ier,ierr,inumz,iel,inumz1,&
       & inumz_min,inumz_max
  logical :: incl_ei !include ei reaction into ground state model if it is not beam-beam interactions
  
  inumz1=inumi-abs(i_nsbeam )
  
  if(i_nsbeam.ge.0) then !only thermal + impurities
     inumz_min=1
     inumz_max=inumz1  
     inumz=inumz1
  else !only beam ion species,fast ion population
     inumz_min=inumz1+1
     inumz_max=inumi
     inumz=abs(i_nsbeam )
  endif
  !excited state stopping
  bms_exc_tot=0.0_fp
  !ground state stopping
  sigvnt=0.0_fp

  if(il_excited) then
     !adas310 library will be used 
     call getsigs_adas310(erel_vec(1:ie_vec,inumz_min:inumz_max)*kev2ev,ie_vec,&
          & z_beam,a_beam,zte*kev2ev,&
          & ztia_vec(inumz_min:inumz_max,1:it_vec)*kev2ev,it_vec,zne,&
          ! t_i/amu [ev]      
          & zni(inumz_min:inumz_max,1:it_vec),&
          & zzi(inumz_min:inumz_max),zza(inumz_min:inumz_max),&
          & inumz,bms_exc_tot(1:ie_vec,1:it_vec,1:ing),ing,ier)
     return
  endif
     

  allocate(conc_imp(inumi,it_vec),sum_conc(it_vec))
  conc_imp=0.0_fp
  allocate(zne_equiv(it_vec))
  zne_equiv=0._fp
  !
  ! target species fraction is z_imp * n_imp/sum(z_imp*n_imp)
  !
  !calculate  z_i*n_i 
  do ii=1,inumi !thermal + impurities and beam ion densities for quasi-neutrality
     conc_imp(ii,:)=zzi(ii)*zni(ii,:)
  enddo
  
  sum_conc(:)=sum(conc_imp(1:inumi,:),dim=1) !sum_i(z_i*n_i)
  
  do ij=1,it_vec
     if( sum_conc(ij).eq.0.0_fp) then
        conc_imp(:,ij)=0.0_fp
     else
        do ii=1,inumi
           conc_imp(ii,ij)=conc_imp(ii,ij)/sum_conc(ij)  !z_i*n_i/sum_i(z_i*n_i)
        enddo
     endif
  enddo
  !
  !calculate electron equivalent density sum_i(z_i*z_i*n_i)
  do ii=1,inumi !thermal + impurities and beam ion densities
     zne_equiv(:)=zne_equiv(:)+zzi(ii)*zzi(ii)*zni(ii,:) !sum_i(z_i**2*n_i)
  enddo
  
  do ij=1,it_vec
     if( sum_conc(ij).eq.0.0_fp) then
        zne_equiv(ij)=0.0_fp
     else
        if(i_nsbeam.ge.0) then
           zne_equiv(ij)=zne_equiv(ij)*zne(ij)/sum_conc(ij) !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)
           ! else assume that n_e = sum_i(z_i*n_i)
        endif
     endif
  enddo
  !
  
  
  ! ground state <sig*v>*n
  
  incl_ei=(i_nsbeam.ge.0)

  call getsigs_adas(erel_vec(1:ie_vec,inumz_min:inumz_max),ie_vec,z_beam,zte,it_vec,zne,inumz,&
       & zzi(inumz_min:inumz_max),&
       & zni(inumz_min:inumz_max,1:it_vec),&
       & ztia_vec(inumz_min:inumz_max,1:it_vec),&
       & sigvnt,incl_ei,ier)
  
  !excited state beam stopping coefficient
  
  if(ier.ne.0) then
     if(allocated(conc_imp))deallocate(conc_imp)
     if(allocated(sum_conc))deallocate(sum_conc)
     if(allocated(zne_equiv))deallocate(zne_equiv)
     return
  endif

  allocate(sigv(ie_vec,it_vec))
  sigv=0._fp
  
  
  do ii=inumz_min,inumz_max
     sigv=0._fp
     zne_equiv(:)=zne_equiv(:)/zzi(ii) !equivalent electron density for each impurity
     !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)/z_i
     
     call adas_bms(erel_vec(:,ii),ie_vec,z_beam,conc_imp(ii,:),&
          !mg             & ztia_vec(ii,:),it_vec,&
          & ztia_vec(ii,:)*zza(ii),it_vec,&
          !mg 12/05/11 i found that bms tables were created with timp [kev] not [kev/amu] as described in fiels
          & zzi(ii:ii),zne_equiv,sigv,ierr) !only diagonal elements will be computed
     
     do ij=1,it_vec
        bms_exc_tot(:,ij,1)=bms_exc_tot(:,ij,1)+sigv(:,ij)*sum_conc(ij)*1.e+6 !<sig*v>*n_e - includes excited states
        !if all plasma species are counted
        ! then sum_conc = zne
     enddo
     
     zne_equiv(:)=zne_equiv(:)*zzi(ii) !resume for next impurity
  enddo
  if(allocated(sigv))deallocate(sigv)
  
  
  if(allocated(conc_imp))deallocate(conc_imp)
  if(allocated(sum_conc))deallocate(sum_conc)
  if(allocated(zne_equiv))deallocate(zne_equiv)
  
end subroutine sigfacz_adas2d

  subroutine sigfacz_adas3d(erel_vec,ie_vec,z_beam,a_beam,zne,ine_vec,zte,zni,zzi,&
       &                     zza,ztia_vec,it_vec,inumi,i_nsbeam,ing, &
       &                     bms_exc_tot,sigvnt,il_excited,ier)
  implicit none
  integer, intent(in) :: ie_vec,ine_vec,it_vec !vector size
  real(fp),intent(in) :: erel_vec(ie_vec,inumi) !interaction energy [kev/amu] beam -> target
  real(fp),intent(in) :: z_beam !beam charge
  real(fp),intent(in) :: a_beam !beam amu
  real(fp),intent(in) :: zne(ine_vec) !electron density [1/cm3]
  real(fp),intent(in) :: zte(it_vec) !electron temperature [kev]
  real(fp),intent(in) :: ztia_vec(inumi,it_vec) !ion species ti/ai temperature [kev/amu]
  real(fp),intent(in) :: zni(inumi,ine_vec) !ion species density [1/cm3]
  integer,intent(in) :: inumi ! number of distinct ion species
  real(fp),intent(in) :: zza(inumi) ! ion species atomic mass
  integer,intent(in) :: i_nsbeam ! number of fast ion species, if positive or zero do not include them
                                 ! to calculation of the <sig*v>*n for ground state and
                                 ! bms for excited state
                                 ! if negative, then fast ion species will be taken into
                                 ! account
  real(fp),intent(in) :: zzi(inumi) ! ion species charge
  logical,intent(in) :: il_excited !.true. -- adas excited state model will be used to calculate
                                   ! enhancement factor
                                   !.false.  -- adas beam stopping coeffitient will be used to
                                 ! estimate enhancement factor
  integer,intent(in) :: ing ! number of h isotopes ng + 3
  real(fp),intent(out) :: sigvnt(ie_vec,ine_vec,it_vec,3) ! ground state sig*v*n : 1 -- sigvne, 2 -- sigviin,3 -- sigvcxn 
  real(fp),intent(out) :: bms_exc_tot(ie_vec,ine_vec,it_vec,ing)
                                              ! il_excited=.false. adas310 bms data: 1 - total excited state  sig*v*n 
                                              ! il_excited=.true. - adas310_fortran_driver:
                                              ! excited state model sig*v*n  1 -- ie ,2 -- ii, 3 -- ii for impuritie
                                              ! cx for ng h isotopes
  real(fp),dimension(:,:),allocatable :: conc_imp !target species fraction   z_imp * n_imp/sum(z_imp*n_imp)
  real(fp),dimension(:),allocatable :: sum_conc !sum(z_imp*n_imp)
  real(fp),dimension(:,:,:), allocatable :: sigv
  real(fp),dimension(:), allocatable :: zne_equiv  !n_e*sum(z_imp*z_imp*n_imp)/sum(z_imp * n_imp)/z_imp

  integer ii,ij,ier,ierr,inumz,iel,inumz1,&
       & inumz_min,inumz_max
  logical :: incl_ei !include ei reaction into ground state model if it is not beam-beam interactions

  inumz1=inumi-abs(i_nsbeam )
  
  if(i_nsbeam.ge.0) then !only thermal + impurities
     inumz_min=1
     inumz_max=inumz1  
     inumz=inumz1
  else !only only beam ion species
     inumz_min=inumz1+1
     inumz_max=inumi
     inumz=abs(i_nsbeam )
  endif
  
  !excited_state stopping
  bms_exc_tot=0.0_fp
  !ground state stopping
  sigvnt=0.0_fp
  if(il_excited) then
     !adas310 library will be used 
     call getsigs_adas310_all(erel_vec(1:ie_vec,inumz_min:inumz_max)*kev2ev,ie_vec,&
          & z_beam,a_beam,zte*kev2ev,&
          & ztia_vec(inumz_min:inumz_max,1:it_vec)*kev2ev,it_vec,zne,&
          ! t_i/amu [ev]      
          & zni(inumz_min:inumz_max,1:ine_vec),ine_vec,&
          & zzi(inumz_min:inumz_max),zza(inumz_min:inumz_max),&
          & inumz,bms_exc_tot(1:ie_vec,1:ine_vec,1:it_vec,1:ing),ing,ier)
     return
  endif
  allocate(conc_imp(inumi,ine_vec),sum_conc(ine_vec))
  conc_imp=0.0_fp
  allocate(zne_equiv(ine_vec))
  zne_equiv=0._fp
!
! target species fraction is z_imp * n_imp/sum(z_imp*n_imp)
!
!calculate  z_i*n_i 
  do ii=1,inumi !thermal + impurities and beam ion densities for quasi-neutrality
     conc_imp(ii,:)=zzi(ii)*zni(ii,:)
  enddo

  sum_conc(:)=sum(conc_imp(1:inumi,:),dim=1) !sum_i(z_i*n_i)
  
  do ij=1,ine_vec
       if( sum_conc(ij).eq.0.0_fp) then
          conc_imp(:,ij)=0.0_fp
       else
          do ii=1,inumi
             conc_imp(ii,ij)=conc_imp(ii,ij)/sum_conc(ij)  !z_i*n_i/sum_i(z_i*n_i)
          enddo
       endif
  enddo
!
!calculate electron equivalent density sum_i(z_i*z_i*n_i)
  do ii=1,inumi !thermal + impurities and beam ion densities
     zne_equiv(:)=zne_equiv(:)+zzi(ii)*zzi(ii)*zni(ii,:) !sum_i(z_i**2*n_i)
  enddo

  do ij=1,ine_vec
     if( sum_conc(ij).eq.0.0_fp) then
        zne_equiv(ij)=0.0_fp
     else
        if(i_nsbeam.ge.0) then
           zne_equiv(ij)=zne_equiv(ij)*zne(ij)/sum_conc(ij) !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)
! else assume that n_e = sum_i(z_i*n_i)
        endif
     endif
  enddo
!
  ! ground state <sig*v>*n

  incl_ei=(i_nsbeam.ge.0)

  call getsigs_adas(erel_vec(1:ie_vec,inumz_min:inumz_max),ie_vec,z_beam,zte,it_vec,zne,ine_vec,inumz,&
       & zzi(inumz_min:inumz_max),&
       & zni(inumz_min:inumz_max,1:ine_vec),&
       & ztia_vec(inumz_min:inumz_max,1:it_vec),&
       & sigvnt,incl_ei,ier)
  
  !excited state beam stopping coefficient
  
  if(ier.ne.0) then
     if(allocated(conc_imp))deallocate(conc_imp)
     if(allocated(sum_conc))deallocate(sum_conc)
     if(allocated(zne_equiv))deallocate(zne_equiv)
     return
  endif
  
  if(allocated(sigv))deallocate(sigv)
  allocate(sigv(ie_vec,ine_vec,it_vec))
  sigv=0._fp
  do ii=inumz_min,inumz_max
     sigv=0._fp
     zne_equiv(:)=zne_equiv(:)/zzi(ii) !equivalent electron density for each impurity
     !n_e*sum_i(z_i**2*n_i)/sum_i(z_i*n_i)/z_i
     
     call adas_bms(erel_vec(:,ii),ie_vec,z_beam,conc_imp(ii,:),&
          !mg             & ztia_vec(ii,:),it_vec,&
          & ztia_vec(ii,:)*zza(ii),it_vec,&
          !mg 12/05/11 i found that bms tables were created with timp [kev] not [kev/amu] as described in fiels
          & zzi(ii:ii),zne_equiv,ine_vec,sigv,ierr) !only diagonal elements will be computed
     
     do ij=1,ine_vec
        bms_exc_tot(:,ij,:,1)=bms_exc_tot(:,ij,:,1)+sigv(:,ij,:)*sum_conc(ij)*1.e+6 !<sig*v>*n_e - includes excited states
        !if all plasma species are counted
        !then sum_conc = zne
     enddo
     
     zne_equiv(:)=zne_equiv(:)*zzi(ii) !resume for next impurity
  enddo
  
  if(allocated(sigv))deallocate(sigv)
  if(allocated(conc_imp))deallocate(conc_imp)
  if(allocated(sum_conc))deallocate(sum_conc)
  if(allocated(zne_equiv))deallocate(zne_equiv)
  
end subroutine sigfacz_adas3d

end module adas_sigfacz

