!-----------------------------------------------------------------------
!  exc_sig -- deal with simplified sigmas excited states depo model
!    dmc 4 may 1994
!
subroutine exc_sig(zexcfac,zbbxfac,zwkfac, &
     zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto, &
     zbbnsva,ierr)
  !
  !  ------------------------------------------------------------
  !     as of 4 may 1994 ...
  !  this routine simply applies the same factor, zexcfac, to all
  !  components of deposition; the thermal target fractions are
  !  not affected.
  !
  !     mod sep 2009 -- mg & dmc -- zexcfac gives the net enhancement
  !  but we assume that ionization reaction rates only, not charge exchange,
  !  are enhanced.  (note: for now, impurity rates are enhanced as if all
  !  ionization even though some portion of stopping is non-neutralizing
  !  charge exchange).  the ionization enhancement factor has to be solved
  !  for; it is output in elements of "zwkfac(...)".
  !
  !  ------------------------------------------------------------
  !
  use nbi_com
  use logmod
  !
  !============
  implicit none
  integer :: i,isb,ii,ierr
  !============
  real(fp) :: zbnsvii,zbnsvie,zbnsviz,zbnsvto !will be enhanced at the output by factor zexcfac_ii
  real(fp) :: zexcfac
  !============
  real(fp) :: zbnsvcx(mig)
  real(fp) :: zbbnsva(2,mibs)!zbbnsva(1,mibs) will be enhanced at the output zexcfac_ii
  !
  !  output:
  !
  !     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
  !     note: beam-fast ion reaction rate will not be changed.
  !
  !     if nsigexc_use.ne.3
  !    zbnsvto:  sum ( n*sigma*v ) for all stopping reactions on thermal
  !              plasma species  -- 1/sec
  !
  !  the following are dimensionless fractions:
  !    zbnsvcx(ig) = portion of zbnsvto due to cx with species ig
  !              i.e. n*sigma*v (cx w/ig) is zbnsvto*zbnsvcx(ig)
  !    zbnsvii = portion of zbnsvto due to ii on thermal ion species
  !    zbnsvie = portion of zbnsvto due to ii on electrons
  !    zbnsviz = portion of zbnsvto due to ii on impurity species
  !
  !  output:  beam-beam corrected reaction rate coefficients:
  !    zbbnsva(1,isb) = ionization n*sigma*v on fast species isb
  !    zbbnsva(2,isb) = cx n*sigma*v on fast species isb
  !
  !  output:  zwkfac(...) enhancement factors:
  !
  real(fp) :: zbbxfac(mibs)
  !
  !    zbbxfac(isb) = cx enhancement factor on fast species isb
  !
  real(fp) :: zwkfac(7)       !  output excited states enhancements
  !                   1 -> electron impact ionization
  !                   2 -> impurity impact ionization
  !                   3 -> th. plasma species impact ionization
  !                   4 -> th. plasma species charge exchange
  !                   5 -> beam beam impact ionization
  !                   6 -> beam beam charge exchange
  !                   7 -> net effect on stopping cross section
  !
  !-------------------------------------
  real(fp) :: znsv_iitherm
  real(fp) :: znsv_cxtherm
  real(fp) :: znsv_iifast
  real(fp) :: znsv_cxfast
  real(fp) :: zbnsvto_orig
  real(fp) :: zexcfac_ii
  !-------------------------------------
  !
  !  see also:  related modules:  exc_pre.fnb  exc_do.fnb
  !
  ierr=0
  if(nsigexc_use.eq.3) then
     !     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

     zwkfac(1:7)=1.0_fp
     zbbxfac(1:nsbeam)=1.0_fp
     return
  end if
  !---------------------------------------------
  if((zexcfac.le.0.99_fp).and.(nbfallgr.ge.1)) then
     call errorlog('exc_sig: error enhancement factor < 0.99')
     call errorlog('   zexcfac = ',zexcfac)
     ierr=1
     return
  end if
  if(zbnsvto.le.0.0_fp) return
  !---------------------------------------------
  !
  !  ground state sums:
  zbnsvto_orig = zbnsvto
  znsv_iitherm = zbnsvto*(zbnsvii+zbnsvie+zbnsviz)
  znsv_iifast  = sum(zbbnsva(1,1:nsbeam))

  if((znsv_iitherm+znsv_iifast).le.0.0_fp) return

  znsv_cxtherm = zbnsvto*(sum(zbnsvcx(1:ng)))
  znsv_cxfast  = sum(zbbnsva(2,1:nsbeam))

  zexcfac = max(1.00_fp,zexcfac)
  !
  !  zexcfac= [zexcfac_ii*(znsv_iitherm+znsv_iifast) + znsv_cxtherm+znsv_cxfast]/
  !             [znsv_iitherm+znsv_iifast + znsv_cxtherm+znsv_cxfast]
  !
  !  zexcfac is input, so:
  !
  !  zexcfac_ii =
  !     (zexcfac*[znsv_iitherm+znsv_iifast + znsv_cxtherm+znsv_cxfast] -
  !      (znsv_cxtherm+znsv_cxfast))/(znsv_iitherm+znsv_iifast)
  !
  !  total thermal target stopping enhancement (new zbnsvto):

  !mg (this is not consistent with a model)
  !mg      zbnsvto = zexcfac*zbnsvto_orig

  zexcfac_ii = ( &
       zexcfac*(znsv_iitherm+znsv_iifast+znsv_cxtherm+znsv_cxfast) &
       - (znsv_cxtherm+znsv_cxfast))/(znsv_iitherm+znsv_iifast)

  zwkfac(1)=zexcfac_ii
  zwkfac(2)=zexcfac_ii
  zwkfac(3)=zexcfac_ii
  zwkfac(4)=1.0_fp
  zwkfac(5)=zexcfac_ii
  zwkfac(6)=1.0_fp
  zwkfac(7)=zexcfac
  !
  do i=1,nsbeam
     zbbxfac(i)=1.0_fp
  end do
  !     recalculate total thermal bms
  zbnsvto = zexcfac_ii*znsv_iitherm+znsv_cxtherm
  !     recalculate ionization fraction
  zbnsvii = zbnsvii*zexcfac_ii*zbnsvto_orig/zbnsvto
  zbnsvie = zbnsvie*zexcfac_ii*zbnsvto_orig/zbnsvto
  zbnsviz = zbnsviz*zexcfac_ii*zbnsvto_orig/zbnsvto
  !mg 06/04/2018 bug fixed, re-calculate fraction of cx reactions that contributes into fast ptcl stopping coefficient
  zbnsvcx(1:ng)=zbnsvcx(1:ng)*zbnsvto_orig/zbnsvto
  !
  !     beam beam enhancement -- ionization only
  !
  do isb=1,nsbeam
     zbbnsva(1,isb)=zbbnsva(1,isb)*zexcfac_ii
  end do
  !
  return
end subroutine exc_sig
