!----------------------------------------------------------------------
module adas_getsigs
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: me_kg,mp_kg
  real(fp),parameter :: zmemp=me_kg/mp_kg
  private :: fp,zmemp 
   
  interface getsigs_adas
     module procedure getsigs_adas1d, getsigs_adas2d, getsigs_adas3d
  end interface getsigs_adas

contains

  subroutine getsigs_adas3d(zeova,ie_vec,z_beam,zte,it_vec,zne, &
       ine_vec,inumz,zz,zn_imp,ztxova, &
       sigvnt,incl_ei,ierr)
    !     program getsigs_adas calculates  ground state stopping coeffisient
    !     electron impact ionization (ei)+ charge exchane (cx)+ ion ionization (ii)
    !     beam particles on plasma species.
    use adas_mod, only : adas_sigvte_ioniz, adas_btsigv
    use adas_mod, only : adas_zstop_sigvz
    !============
    implicit none
    integer, intent(in) :: ie_vec,ine_vec,it_vec !vector size
    real(fp),intent(in) :: zeova(ie_vec,inumz)!interaction energy [kev/amu] beam -> target
    real(fp),intent(in) :: z_beam !beam charge
    real(fp),intent(in) :: zte(it_vec)!electron temperature [kev]
    real(fp),intent(in) :: zne(ine_vec)!electron density [1/cm*3]
    integer, intent(in) :: inumz !# of species (impurities+bakground plasma+fast particles)
    real(fp), intent(in) :: zz(inumz)! z  each species
    real(fp), intent(in) :: zn_imp(inumz,ine_vec)!  density for each impurity [1/cm*3]
    real(fp), intent(in) :: ztxova(inumz,it_vec) ! t_imp/ax for impurities [kev/amu]
    logical, intent(in) :: incl_ei
    !
    real(fp),intent(out),dimension(ie_vec,ine_vec,it_vec,3) :: sigvnt ! <sig*v>*n: 1 -- ei, 2 -- ii, 3 -- cx
    integer,intent(out) ::ierr !control for errors, if =0 means sigvnt.ne.0
    integer, parameter :: neutral=1, ion=2
    !------------------
    ! local
    integer :: ier,i,j,ij
    real(fp),dimension(:), allocatable ::sigvt,sigve,sigvi,zt_vec

    sigvnt=0.0_fp
    !------------------
    allocate(sigvt(it_vec),sigve(ie_vec),sigvi(ie_vec),zt_vec(ie_vec))
    sigvt=0.0_fp
    sigve=0.0_fp
    sigvi=0.0_fp
    zt_vec=0.0_fp
    ierr=1
    if(incl_ei) then
       do ij=1,it_vec
          do i=1,ie_vec
             zt_vec(i)=zte(ij)+0.66666667_fp*zeova(i,1)*zmemp
          end do
          sigve=0.0_fp
          call adas_sigvte_ioniz(z_beam,zt_vec,ie_vec,sigve,ier) ! <sigma*v> in mks
          ierr=ierr*ier
          do j=1,ine_vec
             sigvnt(1:ie_vec,j,ij,1)=sigvnt(1:ie_vec,j,ij,1)+ &
                  sigve(1:ie_vec)*zne(j)*1.e+6 ! cgs conversion * add density
          end do
       end do
    end if
    do i=1,inumz
       do j=1,it_vec
          if(ztxova(i,j).eq.0.0_fp) cycle
          zt_vec(1:ie_vec)=ztxova(i,j)
          sigve=0.0_fp
          sigvi=0.0_fp
          call adas_btsigv(1,neutral,zeova(:,i),zt_vec,ie_vec, &
               z_beam,zz(i),sigve, ier) !cx
          ierr=ierr*ier
          do ij=1,ine_vec
             sigvnt(1:ie_vec,ij,j,3)=sigvnt(1:ie_vec,ij,j,3)+ &
                  sigve(1:ie_vec)*zn_imp(i,ij)*1.e+6
          end do
          call adas_btsigv(2,neutral,zeova(:,i),zt_vec,ie_vec, &
               z_beam,zz(i),sigvi,ier) ! ii
          ierr=ierr*ier
          do ij=1,ine_vec
             sigvnt(1:ie_vec,ij,j,2)=sigvnt(1:ie_vec,ij,j,2)+ &
                  sigvi(1:ie_vec)*zn_imp(i,ij)*1.e+6
          end do
       end do
    end do


    deallocate(sigvt,sigve,sigvi,zt_vec)

  end subroutine getsigs_adas3d

  !
  !
  !

  subroutine getsigs_adas1d(zeova,ie_vec,z_beam,zte,zne, &
       inumz,zz,zn_imp,ztxova, &
       sigvnt,incl_ei,ierr)
    !     program getsigs_adas calculates  ground state stopping coeffisient
    !     electron impact ionization (ei)+ charge exchane (cx)+ ion ionization (ii)
    !     beam particles on plasma species.
    use adas_mod, only : adas_sigvte_ioniz, adas_btsigv
    use adas_mod, only : adas_zstop_sigvz
    !============
    implicit none
    integer, intent(in) :: ie_vec !vector size
    real(fp),intent(in) :: zeova(ie_vec,inumz)!interaction energy [kev/amu] beam -> target
    real(fp),intent(in) :: z_beam !beam charge
    real(fp),intent(in) :: zte(ie_vec)!electron temperature [kev]
    real(fp),intent(in) :: zne(ie_vec)!electron density [1/cm*3]
    integer, intent(in) :: inumz !# of species (impurities+bakground plasma+fast particles)
    real(fp), intent(in) :: zz(inumz)! z  each species
    real(fp), intent(in) :: zn_imp(inumz,ie_vec)!  density for each impurity [1/cm*3]
    real(fp), intent(in) :: ztxova(inumz,ie_vec) ! t_imp/ax for impurities [kev/amu]
    logical, intent(in) :: incl_ei
    !
    real(fp),intent(out),dimension(ie_vec,3) :: sigvnt !<sig*v>*n: 1 -- ei, 2 -- ii, 3 -- cx
    integer,intent(out) ::ierr !control for errors, if =0 means sigvnt.ne.0
    integer, parameter :: neutral=1, ion=2
    !------------------
    ! local
    integer :: ier,i,ij
    real(fp),dimension(:), allocatable ::sigvt,sigve,sigvi,zt_vec

    sigvnt=0.0_fp
    !------------------
    allocate(sigvt(ie_vec),sigve(ie_vec),sigvi(ie_vec),zt_vec(ie_vec))
    sigvt=0.0_fp
    sigve=0.0_fp
    sigvi=0.0_fp
    zt_vec=0.0_fp
    ierr=1
    if(incl_ei) then
       do i=1,ie_vec
          zt_vec(i)=zte(i)+0.66666667_fp*zeova(i,1)*zmemp
       end do
       sigve=0.0_fp
       call adas_sigvte_ioniz(z_beam,zt_vec,ie_vec,sigve,ier) ! <sigma*v> in mks
       ierr=ierr*ier
       sigvnt(:,1)=sigvnt(:,1)+ &
            sigve(:)*zne(:)*1.e+6 ! cgs conversion * add density
    end if
    do i=1,inumz
       if(sum(ztxova(i,:)).eq.0.0_fp) cycle
       zt_vec(1:ie_vec)=ztxova(i,1:ie_vec)
       sigve=0.0_fp
       sigvi=0.0_fp
       call adas_btsigv(1,neutral,zeova(1:ie_vec,i),zt_vec,ie_vec, &
            z_beam,zz(i),sigve, ier) !cx
       ierr=ierr*ier
       sigvnt(1:ie_vec,3)=sigvnt(1:ie_vec,3)+ &
            sigve(1:ie_vec)*zn_imp(i,1:ie_vec)*1.e+6
       call adas_btsigv(2,neutral,zeova(1:ie_vec,i),zt_vec,ie_vec, &
            z_beam,zz(i),sigvi,ier) ! ii
       ierr=ierr*ier
       sigvnt(1:ie_vec,2)=sigvnt(1:ie_vec,2)+ &
            sigvi(1:ie_vec)*zn_imp(i,1:ie_vec)*1.e+6
    end do


    deallocate(sigvt,sigve,sigvi,zt_vec)

  end subroutine getsigs_adas1d
  !
  !
  !


  subroutine getsigs_adas2d(zeova,ie_vec,z_beam,zte,it_vec,zne, &
       inumz,zz,zn_imp,ztxova, &
       sigvnt,incl_ei,ierr)
    !     program getsigs_adas calculates  ground state stopping coeffisient
    !     electron impact ionization (ei)+ charge exchane (cx)+ ion ionization (ii)
    !     beam particles on plasma species.
    use adas_mod, only : adas_sigvte_ioniz, adas_btsigv
    use adas_mod, only : adas_zstop_sigvz
    !============
    implicit none
    integer, intent(in) :: ie_vec,it_vec !vector size
    real(fp),intent(in) :: zeova(ie_vec,inumz)!interaction energy [kev/amu] beam -> target
    real(fp),intent(in) :: z_beam !beam charge
    real(fp),intent(in) :: zte(it_vec)!electron temperature [kev]
    real(fp),intent(in) :: zne(it_vec)!electron density [1/cm*3]
    integer, intent(in) :: inumz !# of species (impurities+bakground plasma+fast particles)
    real(fp), intent(in) :: zz(inumz)! z  each species
    real(fp), intent(in) :: zn_imp(inumz,it_vec)!  density for each impurity [1/cm*3]
    real(fp), intent(in) :: ztxova(inumz,it_vec) ! t_imp/ax for impurities [kev/amu]
    logical, intent(in) :: incl_ei
    !
    real(fp),intent(out),dimension(ie_vec,it_vec,3) :: sigvnt !<sig*v>*n: 1-- ei, 2-- ii, 3 --cx
    integer,intent(out) ::ierr !control for errors, if =0 means sigvnt.ne.0
    integer, parameter :: neutral=1, ion=2
    !------------------
    ! local
    integer :: ier,i,j,ij
    real(fp),dimension(:), allocatable ::sigvt,sigve,sigvi,zt_vec

    sigvnt=0.0_fp
    !------------------
    allocate(sigvt(it_vec),sigve(ie_vec),sigvi(ie_vec),zt_vec(ie_vec))
    sigvt=0.0_fp
    sigve=0.0_fp
    sigvi=0.0_fp
    zt_vec=0.0_fp
    ierr=1
    if(incl_ei) then
       do ij=1,it_vec
          do i=1,ie_vec
             zt_vec(i)=zte(ij)+0.66666667_fp*zeova(i,1)*zmemp
          end do
          sigve=0.0_fp
          call adas_sigvte_ioniz(z_beam,zt_vec,ie_vec,sigve,ier) ! <sigma*v> in mks
          ierr=ierr*ier
          sigvnt(1:ie_vec,ij,1)=sigvnt(1:ie_vec,ij,1)+ &
               sigve(1:ie_vec)*zne(ij)*1.e+6 ! cgs conversion * add density
       end do
    end if
    do i=1,inumz
       do j=1,it_vec
          if(ztxova(i,j).eq.0.0_fp) cycle
          zt_vec(1:ie_vec)=ztxova(i,j)
          sigve=0.0_fp
          sigvi=0.0_fp
          if(zne(j).ne.0.0_fp) then
             call adas_btsigv(1,neutral,zeova(:,i),zt_vec,ie_vec, &
                  z_beam,zz(i),sigve, ier) !cx
             ierr=ierr*ier
             sigvnt(1:ie_vec,j,3)=sigvnt(1:ie_vec,j,3)+ &
                  sigve(1:ie_vec)*1.e+6*zn_imp(i,j)
             !     >              sigve(1:ie_vec)*1.e+6*zn_imp(i,j)*zz(i)/zne(j) !?mg what is *zz(i)/zne(j)
             call adas_btsigv(2,neutral,zeova(:,i),zt_vec,ie_vec, &
                  z_beam,zz(i),sigvi,ier) ! ii
             ierr=ierr*ier
             sigvnt(1:ie_vec,j,2)=sigvnt(1:ie_vec,j,2)+ &
                  sigvi(1:ie_vec)*1.e+6*zn_imp(i,j)
             !     >              sigvi(1:ie_vec)*1.e+6*zn_imp(i,j)*zz(i)/zne(j)!?mg what is *zz(i)/zne(j)
          end if
       end do
    end do


    deallocate(sigvt,sigve,sigvi,zt_vec)

  end subroutine getsigs_adas2d

end module adas_getsigs
