!******************** start file collid.for ; group collid ************
!
!	this subroutine performs fokker - planck collisions
!	transfering energy to the appropriate plasma species
!	and keeps track of various beam quantities such
!	as stored energy and current
!
!	upgrade to separate this section out as a subroutine
!	called by subroutine orbit; rob goldston  12-7-82
!
!	upgrade to include collisional classical transport
!	rob goldston 12-20-82
!
!	new version of beam code allows ions to travel outside
!	the plasma, but inside the limiters.  collid should not
!	be called for such ions, since we assume charge exchange is
!	much faster than any competing process in this region.
!	yion, new common variable, generated from zyion
!	for use in limiter loss program limlos.for
!	rjg july 19, 1983
!
!  dmc -- oct 1991 -- this routine also evaluates beam-target fusion
!    reactions -- support added for additional reactions to be
!    evaluated
!
!  dmc -- using a more exact formula for acceleration due to toroidal
!         e-field, to avoid trouble when the efield is large
!
!  dmc -- anomalous diffusion code removed to new subroutine abdifu
!    4 sep 1990
!
!   21 sep 1990 -- use energy dependent fokker planck coefficient
!      profiles; see subroutine fokker
!
!
!  dmc - full rotation upgrade
!   26 feb 1987
!
!       collisions in plasma frame
!       collisional torque profiles computed
!
!       algorithm modified to accumulate velocity changes -- these
!     are frame independent.  these are then applied to the g.c.
!     velocity and pll velocity
!--------------------------------------------------
!  algorithm... approximate transformation to plasma frame
!    using vpll-->vpll-rflr*omega*bt/b
!
!    note use of *flr* pll and perp velocities-- see nbiflr.fnb
!     rules:
!       rmjflr*vplflr = rmjion*vpll
!       v2toev*vflr*vflr + ephflr =
!                  v2toev*vion*vion + ephion
!
!    taking into account the electrostatic potential difference
!
!    collisions operate on vpll, vperp, as before, but now these
!    are local variables computed in the flr plasma frame
!
!    velocity changes dv=(dvpll,dvperp1,dvperp2) are accumulated
!
!    these are applicable in the lab frame and are used to
!    generate new g.c. vpll and v at the end of collid, assuring
!    change in momentum and energy of the ion consistent with
!    incremental heating and torques applied to plasma during the
!    collision step
!
!--------------------------------------------------
!       **classical transport disabled**
!     in order to compute the momentum effects of classical transport,
!     collid would have to be rewritten to work in the velocity coords
!     (vphi,vrho,vtheta)  rather than  (vpll,vperp1,vperp2), with
!
!       vphi= bt/b * vpll          + bp/b * cospho * vperp1
!                                      + nsnccw * bp/b * sinpho *vperp2
!
!       vrho=                      - sinpho * vperp1
!                                      + nsnccw * cospho *vperp2
!
!     vtheta= -nsnccw*bp/b *vpll   + nsnccw * bt/b * vperp1
!                                      + bt/b * sinpho *vperp2
! mg it is supposed to be written:
! mg     vtheta= -nsnccw*bp/b *vpll   + nsnccw * bt/b * vperp1 * cospho
! mg                                      + bt/b * sinpho *vperp2
!
!    where vperp2 (initially zero) is the perp velocity component
!    normal to the perp velocity on entering collid; sinpho and cospho
!    are the sin, cos of the gyrophase angle (cospho=1 indicates max
!    displacement normal to the flux surface out from the magnetic axis)
!    this is the +rho direction; the +theta direction is in the
!    poloidal plane normal to +rho (it traces the flux surface ccw
!    with theta=0 at the outer midplane crossing); the +phi direction
!    is normal to the poloidal plane in the direction of plasma current.
!
!    the above transformation is anti-orthoganal; the inverse transform
!    is the transpose of the 3 x 3 matrix of coefficients.
!
!    to include classical transport, one would have to transform to
!    (vphi,vrho,vtheta) noting the diffence dphi0 in gyro-avg mechanical
!    momentum (bt/b*vpll*rgc) and instantaneous momentum m*vphi*rflr
!
!    ...then transform to plasma frame by vphi-->vphi-omega*rflr
!    ...then have collisions; evaluate torques by sum delta(m*vphi*r)
!    ...then back transform vphi-->vphi+omega*rflr
!
!    ...then back transform to vpll,vperp1,vperp2 @flr
!    ...evaluate a new guiding ctr position and velocity vector
!       similar to what is done in depsub (when ionizing a neutral;
!       there is a velocity difference due to presence of the radial
!       electrostatic field)
!
!    getting a new difference dphi1 btw. gyro-avg and instantaneous
!    flr momentum.  the net jxb torque would arise from dphi1-dphi0
!    omega*this torque would match the electrostatic potential
!    difference btw old and new guiding ctr's
!
!    ***doing all this would result in a much more complicated and
!    costly algorithm, all to support a tiny effect which is
!    in any case of order (rho(flr)/(system size))**2 i.e. higher
!    order than what is supported elsewhere in the code
!
!    therefore the option to compute classical transport is disabled
!      (d mc 31 mar 1987)
!
!------------------------------------------------------------------
!
subroutine collid
  !
  use nbi_com
  use nbi_random_pll
  use nbi_aep_mod, only: aep_driver
  use nbi_aep_mod, only: rbq_driver
  use nbi_aep_com, only: ilaep_init
  !
  !============
  implicit none
  integer :: isteps,istep,iengy,iz,idum,ingc
  !============
  real(fp) :: zdtgfp,zdcwp,zcwp,zdwp,zbtilt,zangi,zangwi,zomega,zvp
  real(fp) :: zksid,zvprot,zeps,ztex,ztix,zep,ze1,zeinj,zcfa
  real(fp) :: zdfa,zefa,zdvpli,zdvpsi,zdvpdi,zdve,zdvi,zdpe,zdpi
  real(fp) :: zvppov,zdub,zdj,zdj2,zeforc,zdvef,zvpllp,zvpp
  real(fp) :: zdpoh,zepain,zdvs,zzfac,zdvdif,zdpdif,zdifrt,zdelv2
  real(fp) :: zdelv,zepold,zvold,zdves,ztqbe,zdvis,ztqbi,ztqbx
  real(fp) :: zdvdrg,zceff,zswsq,zdelth,zsndth,zcsdth,zrang,zsnphr
  real(fp) :: zcsphr,zxksol,zdvsam,zdvpll,zvpll,zvpp1,zvpp2,zvi
  real(fp) :: zfudg,zdummy
  !============
  logical :: ib_nbsce                 ! beam index nbsce.gt.0
  !
  real(fp) :: zvlab,zvdotb,z_dtgfp
  !-----------------------------------------------------------------------
  !  mod dmc 27 jan 1992 -- if collision processes are very rapid
  !   fpskip.lt.1.0 may be set in orbit.  then loop here over multiple
  !   collision events within one call to collid
  !
  isteps=1
  if(fpskip.lt.1.0_fp) then
     isteps=1+int(1.0_fp/fpskip)
  end if
  zdtgfp=dtgfp/isteps
  !
  !mg nclass dependence removed
  ib_nbsce = (nbsce.gt.0)
  !
  do istep=1,isteps
     !
     !**************  evaluate normalizing constants  ******************
     !
     !
     if(ilaep_init(lsbeam)) then
        if(maxval(aep_ascale).gt.0.0) then
           call aep_driver(bitime(1)-dtgfp* &
                (1.0_fp-real(istep,fp)/real(isteps,fp)))
        endif
        if(minval(aep_ascale).lt.0.0) then
           z_dtgfp=zdtgfp
           call rbq_driver(z_dtgfp)
        endif
        !
        !     if outside beam dist fcn zone region this is deemed an orbit loss!
        !
        !nbiflr called in aep_driver to get nflr, nfrl2 etc
     end if
     zdcwp=zdtgfp*v2toev*wght
     zcwp=v2toev*wght
     zdwp=zdtgfp*wght
     !
     if(nflr.gt.ledge) then
        bdenss(lep1,lsbeam)=bdenss(lep1,lsbeam)+zdwp
        viplasa=viplasa+vion*zdtgfp
         call gyr_larmor_update
         cycle
     end if
     !
     !  ! for torques et al
     !
     zbtilt=btiltc(2,ngc2)
     zangi=rmassb*rmjflr*zbtilt
     zangwi=zangi*wght
     !
     !
     !***********  terms due only to presence of beam particles  **********
     !
     !
     ! 1d outputs				lab frame
     bdenss(nflr,lsbeam)=bdenss(nflr,lsbeam)+zdwp
     ebapl2s(nflr,lsbeam)=ebapl2s(nflr,lsbeam)+vplflr*vplflr*zdcwp
     ebapp2s(nflr,lsbeam)=ebapp2s(nflr,lsbeam)+vppflr*vppflr*zdcwp
     !
     omegabs(nflr,lsbeam)=omegabs(nflr,lsbeam)+ &
          zdwp*zbtilt*vplflr/rmjflr
     !
     if(nlheavy(lsbeam)) then
        xzbmav(nflr,lsbeam)=xzbmav(nflr,lsbeam)+xzbeami*zdwp
     end if
     !
     ! 2d outputs				lab frame
     !   (dmc may 2003: moved to param_fbm using fbm(...))
     !xx      bdens2(nflr2,lsbeam)=bdens2(nflr2,lsbeam)+zdwp
     !xx      eba2pl(nflr2,lsbeam)=eba2pl(nflr2,lsbeam)+vplflr*vplflr*zdcwp
     !xx      eba2pp(nflr2,lsbeam)=eba2pp(nflr2,lsbeam)+vppflr*vppflr*zdcwp
     !
     !xx      vbtor2(nflr2,lsbeam)=vbtor2(nflr2,lsbeam)+zdwp*vplflr
     !
     if(nlorbo) then
        call gyr_larmor_update
        cycle
     endif
     !
     !  more lab frame stuff
     !
     zvdotb=vplflr*bimodb(1)
     zvlab=max(vmins(lsbeam),vflr)
     !
     !----------------------------------------------------------
     !
     !    euclidean transformation from lab frame to rotating plasma frame
     !
     !xx	zomega = fnbint(xiflr,omegag,1,mj,1,idummy,igc)
     zomega = omegag(nflr)+xinflr*(omegag(nflr+1)-omegag(nflr))
     !
     !    collid works on these *plasma frame* variables:
     !      zvp     --- plasma frame ion volocity
     !      viflr(1) -- plasma frame pll velocity
     !      viflr(2) -- plasma frame prp velocity component 1
     !      viflr(3) -- plasma frame prp velocity component 2
     !      zksid   --- viflr(1)/zvp
     !
     zvp=vflr
     viflr(1)=vplflr
     viflr(2)=vppflr
     viflr(3)=0.0_fp
     zksid=xksflr
     !
     !  plasma rotation speed
     zvprot = zomega*rmjflr*zbtilt
     !
     !  transform set (zksid, viflr(...)) to plasma frame
     !
     zeps = zvprot/zvp
     viflr(1)=viflr(1)-zvprot
     zvp=zvp*sqrt(max(0.0_fp,(1.0_fp+zeps*(zeps-2.0_fp*zksid))))
     !
     viflr(1)=max(-0.999999_fp*zvp,min(0.999999_fp*zvp,viflr(1)))
     !
     if(zvp.gt.0.0_fp) then
        zksid=viflr(1)/zvp
     else
        zksid=0.0_fp
     end if
     !
     !----------------------------------------------------------
     !
     ztex=te(nflr,2)
     ztix=ti(nflr,2)
     !
     !----------------------------------------------------------
     !
     zep=zvp*zvp
     ze1=zep*v2toev  ! ion energy, ev, plasma frame
     !
     !  compute energy index for fokker planck coefficients (cf
     !  subroutine fokker)
     !
     zeinj=einjs(lsbeam)
     if(ze1.ge.(fac_e*zeinj)) then
        iengy=1
     else
        if(ze1.ge.(fac_e*fac_e*zeinj)) then
           iengy=2
        else
           if(ze1.ge.(fac_e*fac_e*fac_e*zeinj)) then
              iengy=3
           else
              iengy=4
           end if
        end if
     end if
     !
     !  calculate beam-target neutron production.
     !
     call btfusn(ze1,zdwp)
     !
     !*******************  fokker-planck simulation  ********************
     !
     !  mod dmc jan 1992 -- read out fokker planck coefficients here;
     !  correct heavy ion f-p coeffs to take into account charge state
     !  of ion (cf subroutines fokker and hvsigma)
     !
     if(nlheavy(lsbeam)) then
        iz=xzbeami+0.5_fp
        zcfa=allame(nflr,iz,lsbeam)*cfa(nflr,lsbeam,iengy)
        zdfa=allamx(nflr,iz,lsbeam)*dfa(nflr,lsbeam,iengy)
        zefa=allamx(nflr,iz,lsbeam)*efa(nflr,lsbeam,iengy)
     else
        zcfa=cfa(nflr,lsbeam,iengy)
        zdfa=dfa(nflr,lsbeam,iengy)
        zefa=efa(nflr,lsbeam,iengy)
     end if
     if(minval(aep_ascale).lt.0.0) then
        zcfa=0.0_fp
        zdfa=0.0_fp
        zefa=0.0_fp
     endif
           !
     !  sum velocity changes to ion -- invariant btw frames  ***
     !
     zdvpli = 0.0_fp  ! in parallel/toroidal direction
     zdvpsi = 0.0_fp  ! perp/pldl "same" direction (as original vppflr)
     zdvpdi = 0.0_fp  ! perp/pldl "different" dir. (normal to orig. vppflr)
     !
     if(zvp.le.vmin(nflr,lsbeam)) then
        !  check for already - thermalized
        zdve=0.0_fp
        zdvi=0.0_fp
        zdpe=0.0_fp
        zdpi=0.0_fp
        zvppov=0.0_fp
        goto 90
     end if
     !
     !  j dot b sum (renorm - cf outptb for units conversion)
     zdub=viflr(1)*zdwp*bimodb(1)      ! v.b
     zdj=xzbeami*zdub
     pljbbms(nflr,lsbeam)=pljbbms(nflr,lsbeam)+zdj     ! unshielded
     !  j toroidal sum (beam grid; cf outptb)
     zdj2=xzbeami*viflr(1)*zdwp*zbtilt
     xcurb2(nflr2,lsbeam)=xcurb2(nflr2,lsbeam)+zdj2 ! unshielded
     !
     !  as above, for ions from *beam sources only*
     !  dmc sept 1990 -- nbsce=0 for fusion product fast ions
     if(nbsce.gt.0) then
        if(nlco(nbsce)) then
           pljbcct(nflr,lsbeam)=pljbcct(nflr,lsbeam)+zdj
           xcrb2c(nflr2,lsbeam)=xcrb2c(nflr2,lsbeam)+zdj2
        end if
        pljbbm_nb(nflr,nbsce)=pljbbm_nb(nflr,nbsce)+zdj
     end if
     !
     !  energy resolved sums of density and v.b; use lab frame info
     !
     bnbs_ien(nflr,nbien,lsbeam)=bnbs_ien(nflr,nbien,lsbeam)+zdwp
     upbs_ien(nflr,nbien,lsbeam)= &
          upbs_ien(nflr,nbien,lsbeam)+zvdotb*zdwp
     if(ib_nbsce) then
        bnb_ien(nflr,nbien,nbsce)=bnb_ien(nflr,nbien,nbsce)+zdwp
        upb_ien(nflr,nbien,nbsce)=upb_ien(nflr,nbien,nbsce)+zvdotb*zdwp
     end if
     !
     !  toroidal electric field effects
     !	xjbfac consistent with curb installed 12/27/82 rjg
     !
     !  local e field with shielding correction
     !
     !  mod dmc 24 oct 1997 -- use vpoh instead of vloop -- more
     !  stable for q analysis mode if user chooses nmodpoh=2 namelist option
     !
     if(nlbcoh) then
        zeforc=9.58e11_fp*xzbeami*vpoh(nflr)*xjbfac(nflr)*zbtilt/ &
             (twopi*rmjflr*abeam)
        !  acceleration
        zdvef=zeforc*zdtgfp
        zvpllp=viflr(1)
        viflr(1)=viflr(1)+zdvef
        zdvpli=zdvpli+zdvef
        !
        !  ohmic torque on beam ions
        tqohbs(nflr,lsbeam)=tqohbs(nflr,lsbeam)+zangwi*zdvef
        !
        !  dmc may 1991 -- exact formula for velocity change
        !   derived from vector sum; law of cosines; cos(alpha)=-zvpllp/zvpp
        !   as zdvef = change in velocity in pll direction only
        !
        zvpp=zvp
        zvp=sqrt(zvpp**2+zdvef**2+2.0_fp*zdvef*zvpllp)
        zdpoh=(zvp-zvpp)*(zvp+zvpp)
        pbohs(nflr,lsbeam)=pbohs(nflr,lsbeam)+zdpoh*zcwp
        !
        zksid=viflr(1)/zvp
        !
        zep=zvp*zvp
        !
     end if
     !
     zvppov=sqrt(1.0_fp-zksid**2)
     !
     !  electron drag (correct to first order in te/e)
     zepain=1.0_fp/(zep*abeam)
     zdve=zcfa*zdtgfp*zvp*(1._fp-ztex*1.915e12_fp*zepain)
     !
     !  ion drag (correct to first order in ti/e)
     zdvi=zdfa*zdtgfp*zepain*abeam* &
          (1._fp+ztix*9.577e11_fp*zepain)
     !
     !
     !  dmc feb 1994 -- guard against too much drag at thermalization time
     !  -change in velocity due to drag should not exceed the velocity itself
     !
     zdvs=zdve+zdvi
     if(zdvs.gt.zvp) then
        zzfac=zvp/zdvs    ! make it so drag results in zvp=0.0
        zdve=zzfac*zdve
        zdvi=zzfac*zdvi
        zdvs=zdve+zdvi
     end if
     !
     zdpe=(2._fp*zvp-zdve)*zdve
     zdpi=(2._fp*zvp-zdvi)*zdvi
     !
     zvp=zvp-zdvs
     zep=zvp*zvp
     !
     !  check for disabling of energy diffusion
     !   also disable if ptcl has thermalized
     !
90   continue
     if((.not.nlbcde).or.(zvp.le.vmin(nflr,lsbeam))) then
        zdvdif=0.0_fp
        zdpdif=0.0_fp
        zdifrt=0.0_fp
        goto 100
     end if
     !
     !  energy diffusion on both electrons and ions
     !
     zdelv2=zdtgfp*(zcfa*ztex*1.915e12_fp/abeam+ &
          ztix*1.915e12_fp/(zep*zvp*abeam)*zdfa)
     !
     zdelv=sqrt(zdelv2)
     !
     zepold=zep
     zvold=zvp
     zvp=max(0.5_fp*vmin(nflr,lsbeam),nbgauss_pll(zdelv,zvp))
     !
     zep=zvp*zvp
     !
     zdpdif=zepold-zep
     zdvdif=zvold-zvp
     !
     zdifrt=zcfa*ztex*zepold*zvold/ &
          (zdfa*ztix)
     zdifrt=zdifrt/(zdifrt+1._fp)
     !
     !  accumulate drag and diffusion energy transfer to electrons and ions,
     !    add up torques to thermal plasma species
     !
100  continue
     zdpe=zdifrt*zdpdif+zdpe
     zdves=zdifrt*zdvdif+zdve
     pbes(nflr,lsbeam)=pbes(nflr,lsbeam)+zdpe*zcwp
     ztqbe=zdves*zksid*zangwi
     !
     zdpi=(1._fp-zdifrt)*zdpdif+zdpi
     zdvis=(1._fp-zdifrt)*zdvdif+zdvi
     pbis(nflr,lsbeam)=pbis(nflr,lsbeam)+zdpi*zcwp
     ztqbi=zdvis*zksid*zangwi
     ztqbx=pbxfrc(nflr,lsbeam)*ztqbi
     !
     !  nclass...
     pbes_ien(nflr,nbien,lsbeam)=pbes_ien(nflr,nbien,lsbeam)+zdpe*zcwp
     pbis_ien(nflr,nbien,lsbeam)=pbis_ien(nflr,nbien,lsbeam)+zdpi*zcwp
     if(ib_nbsce) then
        pbe_ien(nflr,nbien,nbsce)=pbe_ien(nflr,nbien,nbsce)+zdpe*zcwp
        pbi_ien(nflr,nbien,nbsce)=pbi_ien(nflr,nbien,nbsce)+zdpi*zcwp
     end if
     !
     !  velocity change thus far...
     zdvdrg=zdves+zdvis
     zdvpli=zdvpli-zksid*zdvdrg
     zdvpsi=zdvpsi-zvppov*zdvdrg
     !
     !  check for thermalization - set flag
     if(zvp.le.vmin(nflr,lsbeam)) then
        nbflg=2
        goto 150
     end if
     !
     !  pitch angle scattering
     !	using technique of shanny et al. phys. flu. 10 (1967) 1281
     !	in order to be able to include classical transport
     !
     if(nlbcpa) then
        !
        zceff=2._fp*zefa/(zep*zvp)*zdtgfp
        zswsq=(1._fp-zksid*zksid)*zceff
        !
        !	old technique:
        !
        !x	zsw=sqrt(zswsq)
        !x	zksid=nbgauss_pll(zsw,zksid*(1.-zceff))
        !x	if(zksid.gt..999) zksid=.999
        !x	if(zksid.lt.-.999) zksid=-.999
        !
        !	new technique:
        !
        zdelth=sqrt(-2.0_fp*zceff*log(nbrandm_pll(idum)))
        zsndth=sin(zdelth)
        zcsdth=cos(zdelth)
        !
        !	compute the sine and cosine of the (random) angle between
        !	the scatter and going towards vpar = v
        !
        !	this section of code evaluates cos and sin of a random
        !	angle, using different techniques on the vax and on the
        !	cray
        !	rob goldston 12/15/82
        !
        zrang=nbrandm_pll(idum)*twopi		!@
        zsnphr=sin(zrang)			!@
        zcsphr=cos(zrang)			!@
        !
        !@ 120	continue
        !@	zxposi=2.0*nbrandm_pll(idum)-1.0
        !@	zyposi=2.0*nbrandm_pll(idum)-1.0
        !@	zsqrad=zxposi*zxposi+zyposi*zyposi
        !@	if((zsqrad.gt.1.0).or.(zsqrad.lt.1.0e-15)) goto 120
        !@	zsqrad=1.0/zsqrad
        !@	zsnphr=2.0*zxposi*zyposi*zsqrad
        !@	zcsphr=(zxposi*zxposi-zyposi*zyposi)*zsqrad
        !
        !	do actual pitch angle scattering
        !
        zxksol=zksid
        zksid=zsndth*zvppov*zcsphr+zcsdth*zxksol
        !
        zdvsam=zvp*((zcsdth-1.0_fp)*zvppov-zsndth*zcsphr*zxksol)
        zdvpsi=zdvpsi+zdvsam
        !
        zdvdif=zvp*zsndth*zsnphr
        zdvpdi=zdvpdi+zdvdif
        !
        zksid=max(-1.0_fp,min(1.0_fp,zksid))
        !
        !  accumulate torques
        !
        zdvpll=zvp*(zxksol-zksid)
        zdvpli=zdvpli-zdvpll
        !
        ztqbi=ztqbi+zangwi*zdvpll
        ztqbx=ztqbx+paxfrc(nflr,lsbeam)*zangwi*zdvpll
        !
     else
        zceff=0.0_fp
        zswsq=0.0_fp
     end if
     !
     !  accumulate expected drag and pitch angle scattering for goose
     !	calculation
     !
     viplasa=viplasa+zvp*zdtgfp
     dvg=zdve+zdvi+dvg
     dksgc=dksgc+zceff*zksid
     dksgd=dksgd+zswsq
     !
     !-------------------- end of collisions with plasma -----------------
     !  deposit torques
     !  switch back to lab frame
     !
150  continue
     !
     tqbes(nflr,lsbeam)=tqbes(nflr,lsbeam)+ztqbe
     tqbis(nflr,lsbeam)=tqbis(nflr,lsbeam)+ztqbi
     !tq
     if(ib_nbsce) then
        tqbco_nb(nflr,nbsce)=tqbco_nb(nflr,nbsce)+ztqbe+ztqbi
        ! collisional torque is a sum tqbe tqbi separated by beam
     end if

     wtqbxfs(nflr,lsbeam)=wtqbxfs(nflr,lsbeam)+abs(ztqbi)
     tqbxfs(nflr,lsbeam)=tqbxfs(nflr,lsbeam)+abs(ztqbx)
     !
     !
     !    lab frame...
     !       apply velocity changes calculated in plasma frame
     !--------------------------------------------------------------------
     !  reconstitute the guiding ctr velocities - lab frame
     !
     zvpll=vpll+zdvpli
     zvpp1=vperp+zdvpsi
     zvpp2=zdvpdi
     zvi=sqrt(zvpll*zvpll+zvpp1*zvpp1+zvpp2*zvpp2)
     zvi=max(1.0_fp,zvi)  ! at least 1 cm/sec, 1/v2toev << 1 ev
     !
     vion=zvi
     ep=vion*vion
     !
     zfudg=0.999999_fp  ! bind v away from bdy ...
     vpll=max(-zfudg*vion,min(zfudg*vion,zvpll))
     xksid=vpll/vion
     !
     vperp=sqrt(1.0_fp-xksid*xksid)*vion
     !
     bphion=rmassb*vpll*rmjion*zbtilt
     !
     if(nbflg.ne.2) then
        call gyr_larmor_update
        cycle
     endif
     !
     !-------------------------------------------------------------------
     !
     !  ** thermalization **
     !      deposit at old g.c. to avoid jxb correction
     !
     ingc=min(ledge,ngc)
     if(nbs.gt.0) sorbhs(ingc,nbs,lsbeam)= &
          sorbhs(ingc,nbs,lsbeam)+wght
     sorbths(ingc,lsbeam)=sorbths(ingc,lsbeam)+wght
     tqbths(ingc,lsbeam)=tqbths(ingc,lsbeam)+bphion*wght
     ! thermalization torque
     if(ib_nbsce) then
        tqbth_nb(ingc,nbsce)=tqbth_nb(ingc,nbsce)+bphion*wght
     end if
     !
     bphion=0.0_fp
     pbths(ingc,lsbeam)=pbths(ingc,lsbeam)+ &
          (zep+(ephflr-ephion)/v2toev)*zcwp
     !beam thermalization power
     if(ib_nbsce) then
        pbth_nb(ingc,nbsce)=pbth_nb(ingc,nbsce)+ &
             (zep+(ephflr-ephion)/v2toev)*zcwp
     end if
     ep=0.0_fp
     zomega=0.5_fp*(omegag(ingc)+omegag(ingc+1))
     bpthsfs(lsbeam)=bpthsfs(lsbeam)-zomega*rmjion*zomega*rmjion*zcwp
     pbthr2s(ingc,lsbeam)=pbthr2s(ingc,lsbeam)+wght*rmjion*rmjion
     return
  end do  ! end of isteps loop
  return

contains

  subroutine gyr_larmor_update
    if(istep.lt.isteps) then
       !  multiple steps -- new larmour radius and gyrophase angle
       call nbigc(0)
       call gphase
       call nbiflr(zdummy,0)
    end if
  end subroutine gyr_larmor_update
  !
end subroutine collid
!******************** end file collid.for ; group collid ***************
