!-----------------------------------------------------------------------
!  beamiz0 -- promote ionization level of orbitting heavy beam ion
!    subroutine of beamiz
!
subroutine beamiz0(isw,ievent,iflri,zdtiz0,zangph)
  !
  use nbi_com
  use nbi_random_pll
  !
  !%      use sglib_mod
  !-----------------------------------------------------------------------
  !
  !  see comments in the caller routine -- beamiz
  !
  !  the purpose of this routine is to compute the effect of stripping
  !  an electron from an orbitting heavy ion-- this involves changing the
  !  ion charge, and its gyro radius and hence its guiding center.
  !
  !  passed arguments
  !    isw -- output set =1 if the orbit changes
  !    ievent -- output set =1 if a stripping event occurs
  !    iflri -- input -- if =1, integrate the gyrophase angle instead
  !      of computing a new gyrophase angle randomly
  !    zdtiz0 -- input -- the timestep
  !    zangph -- input/output -- the gyrophase angle, if iflri=1.
  !
  !-----------------------------------------------------------------------
  !  local:
  !
  !============
  implicit none
  integer :: ievent,iflri,isw,iz,ictiz,idum,inoflr,ilim
  !============
  real(fp) :: zdtiz0,zangph,zdtiz,ztaui,zdtsiz,ztscxp,zfact,zdtnow
  real(fp) :: zxi0,zth0,zbmod,zbtflr,zr,zy,zxi,zepogc,zphigc,zepoflr
  real(fp) :: zphiflr,zbtgc,zthea,zphia,zlndep,zeflr,zegc,zdummy
  real(fp) :: zrmji,zyi,zrl,zsinph,zcosph,zsinth,zcosth
  !============
  integer :: iwarn=50
  !
  real(fp) :: zbdir(4)  ! flr gyro point field data
  !
  !%      real(fp) :: zrmjpl(100),zypl(100)
  !%      integer :: inpl
  !-----------------------------------------------------------------------
  !
  !%      inpl=1
  !%      zrmjpl=0; zypl=0
  !%      zrmjpl(inpl+1)=rmjion
  !%      zypl(inpl+1)=yion
  !
  !  current ionization level
  !
  iz=xzbeami+0.5_fp
  !
  !  timestep
  !
  zdtiz=zdtiz0
  !
  !  event loop -- .gt. 1 event / timestep is possible
  !
  ictiz=0
100 continue
  do
  !
     if(nflr.gt.ledge) exit
  !
     ictiz=ictiz+1
     !
     iz=iz+1
     !
     ztaui=hvsigi(nflr,iz,lsbeam)
     !
     zdtsiz=zdtiz*ztaui
     !
     !  test for event
     !
     ztscxp=tscx
     tscx=tscx+zdtsiz
     if(tscx.le.ttcx) goto 1000
  !
  !-----------------------------------------------------------------------
  !  event detected!
  !
  !  time location w/in step
  !
     zfact=(ttcx-ztscxp)/(tscx-ztscxp)
     zdtnow=zfact*zdtiz
     !
     !  time left
     !
     zdtiz=(1.0_fp-zfact)*zdtiz
     !
     !  rel. time to subsequent event
     !
     ttcx=-log(nbrandm_pll(idum))
     tscx=0.0_fp
     !
  !-----------------------------------------------------------------------
  !  process event
  !
     ievent=1
     isw=1
  !
  !  update radial current sum
  !
     call radcur(xistate,xiflr,wght,tirorbs)
     xistate=xiflr
  !
  !  get field data out on the gyro
  !
     zxi0=xiflr
     zth0=thflr
  !
     call nblocb(rmjflr,yflr,xiflr,thflr,zbmod,zbdir,2,0,idum)
     zbtflr=1.e-4_fp*zbmod*zbdir(1)
  !
  !  prepare to call cxline -- as in nutrav -- get info out on gyro
  !
     zcosph=cospho
     zsinph=sinpho
     !
     zr=rmjflr
     zy=yflr
     zxi=xiflr
     !
     zepogc=ephion
     zphigc=bphion
     !
     !  get direction angles of ion out on the gyro and it momentum and
     !  potential energy
     !
     call cxline(zcosph,zsinph,zr,zy,zxi,zepogc,zphigc,zbdir, &
          zepoflr,zphiflr)
     !
     !  jxb gc -> flr; zepo argumenter order reversed because these are
     !  potential rather than kinetic energies
     !
     call nb0jxb(wght,zepoflr,zepogc,zphigc,zphiflr,zbtgc,zbtflr, &
          ngc,nflr,xion,xiflr,2,.false.)
     !
     !  trajectory angle info
     !
     zthea=degrad*thea
     zphia=degrad*phia
     !
     zsinth=sin(zthea)
     zcosth=cos(zthea)
     !
     zsinph=sin(zphia)
     zcosph=cos(zphia)
     !
     !-----------------------------------------------------------------------
     !
     !  promote the ionization level and then "recapture" the ion to a
     !  gyro, similar to nutrav, except that the neutral track is of
     !  zero length
     !
     xzbeami=iz                    ! promote ionization level
     sorbes(nflr,lsbeam)=sorbes(nflr,lsbeam)+wght  ! supply electron source
     !
     zlndep=0.0_fp
     !
     !  get new initial orbit parameters
     !
     if(nlbflr) then
        inoflr=0
     else
        inoflr=1
     end if
     !
     call depsub(zlndep,zsinth,zcosth,zsinph,zcosph,inoflr, &
          zphiflr,zphigc,zeflr,zegc,zbtflr,zbtgc)
     !
     !  jxb:  flr -> new gc
     !
     call nb0jxb(wght,zeflr,zegc,zphiflr,zphigc,zbtflr,zbtgc, &
          nflr,ngc,xiflr,xion,3,.false.)
     !
     !  update gc common parameters
     !
     call nbigc(1)
     !
     !%      inpl=inpl+1
     !%      zrmjpl(inpl+1)=rmjion
     !%      zypl(inpl+1)=yion
     !
     !  new random gyrophase
     !
     if(nlbflr) then
        if(iflri.eq.0) then
           call gphase
        else
           zangph=zangph+zdtnow*vperp/rlar
           cospho=cos(zangph)
           sinpho=sin(zangph)
        end if
     end if
     !
     !  update flr common parameters
     !
     call nbiflr(zdummy,0)
     !
     !  check for additional events within this timestep
     !
  enddo
    !
     !-----------------------------------------------------------------------
  !
  !  if iflri is set, do limiter loss checking here instead of in
  !  subroutine orbit.
  !
  if(iflri.eq.1) then
     !
     zrmji=rmjion
     zyi=yion
     zrl=rlar
     zxi=xion
     !
     rmjion=rmjflr
     yion=yflr
     rlar=0.0_fp
     xion=xiflr
     !
     if(nlbgflr.and.(iwarn.gt.0)) then
        iwarn=iwarn-1
        write(nonlin,*) ' %beamiz0 -- nlbgiflr=.true. ignored.'
        write(nonlin,*) '  no generalized flr for heavy beam ions.'
     end if

     call xlimlos(ilim)
     if(ilim.eq.1) ievent=2
     !
     rmjion=zrmji
     yion=zyi
     rlar=zrl
     xion=zxi
     !
  end if
  !
1000 continue
  !
  if(iflri.eq.1) then
     zangph=zangph+zdtiz*vperp/rlar
     cospho=cos(zangph)
     sinpho=sin(zangph)
     call nbiflr(zdummy,0)
  end if
  !
  !%      if((isw.eq.1).and.(nbplot.gt.0)) then
  !%        zrmjpl(1)=inpl
  !%        zypl(1)=inpl
  !%        call symbl(2)
  !%        call line(-1)
  !%        call r8_cplot(zrmjpl,zypl)
  !%        call tsend
  !%      end if
  !
  return
end subroutine beamiz0
