!---------------------------------------------------------------------
!  sigfacz -- interface to excited states stopping enhancement fit
!
!    sigfacz parcels the z>1 densities into incremental zeff's and
!    calls sigfacx which evaluates a stopping enhancement based on
!    the fit coded in sigmas.for
!
!  mod jec 14oct 2003 -- suppress erel error if erel .lt. previous
!                        high value
!
!  mod dmc 6 may 1994 -- the fit is weak for low energies especially
!    at high zeff's.  resetting emin to 20kev ...
!
real(fp) function sigfacz(nonlin, &
     erel,amu,zne,zte,zni,zzi,ztia,inumi, &
     lambdgs,lambdexc)
  !
  !  input:
  !  nonlin -- i/o unit for messages
  !    erel -- interaction energy (ev) beam -> target
  !     amu -- beam mass (amu)
  !     zne -- electron density /cm3
  !     zte -- electron temperature ev
  !     zni -- ion densities /cm3
  !     zzi -- ion species z numbers
  !    ztia -- ion species ti/ai
  !   inumi -- number of distinct ion species
  !
  !  note the fit does not depend on ion temperature!
  !  i am using ti/ai as a correction to the e/amu "average"
  !  interaction energy...
  !
  !============
  use iso_c_binding, only: fp => c_double
  implicit none
  integer :: inumi,iwarneh,iwarnel,iwarnt,iwarnn,i,iz,izp1,nonlin
  !============
  real(fp) :: ztiav,znito,zzil,zefi,zw2,zw1,zcorr,zrel,zminr,zmaxr
  real(fp) :: zzne,zminn,zmaxn,zzte,zmint,zmaxt,zfac,sigfacx
  !============
  real(fp) :: erel,amu,zne,zte,zni(inumi),zzi(inumi),ztia(inumi)
  !
  !  output:  sigfacz = enhancement factor (over ground state)
  !     of total beam stopping cross section, due to excitation
  !     effects (estimated from fit data)
  !
  !  addl output:
  !   lambdgs -- ground state model beam mean free path, cm
  !   lambdexc-- excited states model beam mean free path, cm
  !
  real(fp) :: lambdgs,lambdexc
  !
  !   (sigfacz = lambdgs/lambdexc)
  !
  !----------------
  !   local -- to use the fit routine the ion densities must be expressed
  !     as incremental zeff's at z=2,6,8,26.  input z's .gt. 26 are in
  !     error.  code here parses the densities accordingly
  !
  real(fp) :: zf(4),zz(4)
  integer :: iwarn,jwarn
  !
  data iwarn/0/
  data jwarn/4/
  data iwarneh/0/
  data iwarnel/0/
  data iwarnt/0/
  data iwarnn/0/
  !
  data zz/2.0_fp,6.0_fp,8.0_fp,26.0_fp/
  !----------------
  !
  zf(1:4)=0.0_fp
  !
  ztiav=0.0_fp
  znito=0.0_fp
  do i=1,inumi
     ztiav=ztiav+zni(i)*ztia(i)
     znito=znito+zni(i)
     if(zzi(i).gt.1.0_fp) then
        zzil=zzi(i)
        if((zzil.lt.2.0_fp).or.(zzil.gt.26.0_fp)) then
           if((zzil.lt.1.99_fp).or.(zzil.gt.26.01_fp)) then
              iwarn=iwarn+1
              if(iwarn.le.jwarn) write(nonlin,1001) zzil
           end if
           zzil=min(26.0_fp,max(2.0_fp,zzil))
        end if
1001    format(' %sigfacz:  ion z out of range, z=',1pe11.4, &
             ' (patched).')
        do iz=1,3
           izp1=iz+1
           if((zz(iz).le.zzil).and.(zzil.le.zz(izp1))) exit
        end do
        !
        !  distribute the zeff increment
        !
        zefi=zni(i)*zzi(i)*zzi(i)/zne
        zw2=(zzil-zz(iz))/(zz(izp1)-zz(iz))
        zw1=1.0_fp-zw2
        zf(iz)=zf(iz)+zw1*zefi
        zf(izp1)=zf(izp1)+zw2*zefi
        !
     end if  ! zi > 1
     !
  end do
  !
  ztiav=ztiav/znito
  !
  !  compute enhancement factor (sigfacx function, calls sigmas function)
  !
  zcorr=1.0_fp
  !
  !  dmc-- the atomic physics energy tables from preact go up to 25.0e6 ev!
  !        these tables only go up to 10e6 ev!  warn a few times about this
  !        but not incessantly.  the enhancement factor for 10 mev is used.
  !        this is not going to affect any known neutral beam deposition
  !        calculations.
  !
  zrel=erel/amu+1.5_fp*ztiav
  zminr=20.0e3_fp  ! reset to this from 10.0e3 due to weakness of fit
  zmaxr=10000.0e3_fp
  if((zrel.lt.zminr).or.(zrel.gt.zmaxr)) then
     if(zrel.lt.zminr) then
        iwarnel=iwarnel+1
        if(iwarnel.le.jwarn) then
           write(nonlin,2001) zrel,zminr,zmaxr
        end if
        zcorr=zcorr*(zrel/zminr)
     else
        iwarneh=iwarneh+1
        if(iwarneh.le.jwarn) then
           write(nonlin,2001) zrel,zminr,zmaxr
        end if
     end if
2001 format( &
          ' %sigfacz (exc.states) erel=',1pe11.4,' out of range'/ &
          '   min=',1pe11.4,' max=',1pe11.4,' ev/amu; fixup taken.')
     if(zrel.lt.zminr) then
     end if
     zrel=min(zmaxr,max(zminr,zrel))
  end if
  !
  zzne=zne
  zminn=1.0e12_fp
  zmaxn=1.0e15_fp
  if((zzne.lt.zminn).or.(zzne.gt.zmaxn)) then
     iwarnn=iwarnn+1
     if((iwarnn.le.jwarn).or.(zzne.gt.1.05_fp*zmaxn)) then
        write(nonlin,2002) zzne,zminn,zmaxn
2002    format( &
             ' %sigfacz (exc.states) ne=',1pe11.4,' out of range'/ &
             '   min=',1pe11.4,' max=',1pe11.4,' /cm3; fixup taken.')
     end if
     if(zzne.lt.zminn) then
        zcorr=zcorr*(zzne/zminn)
     end if
     zzne=min(zmaxn,max(zminn,zzne))
  end if
  !
  zzte=zte
  zmint=1000.0_fp
  zmaxt=50000.0_fp
  if((zzte.lt.zmint).or.(zzte.gt.zmaxt)) then
     iwarnt=iwarnt+1
     if((iwarnt.le.jwarn).or.(zzte.gt.1.05_fp*zmaxt)) then
        write(nonlin,2003) zzte,zmint,zmaxt
2003    format( &
             ' %sigfacz (exc.states) te=',1pe11.4,' out of range'/ &
             '   min=',1pe11.4,' max=',1pe11.4,' ev; fixup taken.')
     end if
     zzte=min(zmaxt,max(zmint,zzte))
  end if
  !
  zfac=sigfacx(zrel,1.0_fp,zzne,zzte,zf,lambdgs,lambdexc)
  sigfacz=1.0_fp+max(0.0_fp,zcorr*(zfac-1.0_fp))
  !
  return
end function sigfacz
