anchor.diff:  !!     RB: adding as input arguments to jaw variables icoll,iturn,ipart
   !!         these are only used for the writeout of particle histories
   !!
  -!!++  Input:   ZLM is interaction length
  -!!++           MAT is choice of material
  +!! Input:   ZLM is interaction length
  +!!          MAT is choice of material
   !!
  -!!++  Output:  nabs = 1   Particle is absorped
  -!!++           nabs = 4   Single-diffractive scattering
  -!!++           dpop       Adjusted for momentum loss (dE/dx)
  -!!++           s          Exit longitudinal position
  +!! Output:  nabs = 1   Particle is absorped
  +!!          nabs = 4   Single-diffractive scattering
  +!!          dpop       Adjusted for momentum loss (dE/dx)
  +!!          s          Exit longitudinal position
   !!
  -!!++  Physics:  If monte carlo interaction length greater than input
  -!!++            interaction length, then use input interaction length
  -!!++            Is that justified???
  -!!
  -!!     nabs=1....absorption
  -!!
  +!! Physics:  If monte carlo interaction length greater than input
  +!!           interaction length, then use input interaction length
  +!!           Is that justified???
   !<
  -subroutine k2coll_jaw(s,nabs,icoll,iturn,ipart)
  +subroutine k2coll_jaw(s, nabs, icoll, iturn, ipart)
   
     use mod_ranlux
     use coll_common
  
positive.code:   use hdf5_output
  #endif
  
    implicit none
  
    integer nabs,inter,iturn,icoll,ipart,nabs_tmp ! RB: added variables icoll,iturn,ipart for writeout
    real(kind=fPrec) m_dpodx     !CT, RB, DM
    real(kind=fPrec) p,rlen,s,t,dxp,dzp,p1,zpBef,xpBef,pBef
  !...cne=1/(sqrt(b))
  !...dpodx=dE/(dx*c)
  
  !++  Note that the input parameter is dpop. Here the momentum p is
  !++  constructed out of this input.
  
    p=p0*(one+dpop)
    nabs=0
    nabs_tmp=nabs
  
    if(mat.eq.nmat) then
  !++  Collimator treated as black absorber
      nabs=1
      nabs_tmp=nabs
      s=zero
  
      if(dowrite_impact) then
        ! write coll_scatter.dat for complete scattering histories
  
negative.code:     return
    end if
  
    ! Initialize the interaction length to input interaction length
    rlen = zlm
  
    ! Do a step for a point-like interaction.
    ! Get monte-carlo interaction length.
  10 continue
    zlm1     = (-one*xintl(mat))*log_mb(real(rndm4(),fPrec))
    nabs_tmp = 0  ! type of interaction reset before following scattering process
    xpBef    = xp ! save angles and momentum before scattering
    zpBef    = zp
    pBef     = p
  
    ! If the monte-carlo interaction length is longer than the
    ! remaining collimator length, then put it to the remaining
    ! length, do multiple coulomb scattering and return.
    ! LAST STEP IN ITERATION LOOP
    if(zlm1 > rlen) then
      zlm1 = rlen
      call k2coll_mcs(s)
      s = (zlm-rlen)+s
  #ifdef MERLINSCATTER
      call merlinscatter_calc_ion_loss(p,edens(mat), pleng(mat),exenergy(mat),s,m_dpodx)
      p = p-m_dpodx
  #else
      call k2coll_calcIonLoss(mat,p,rlen,m_dpodx)  ! DM routine to include tail
      p = p-m_dpodx*s
  #endif
  
      dpop = (p-p0)/p0
      if(dowrite_impact) then
        ! write coll_scatter.dat for complete scattering histories
  #ifdef HDF5
  