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: #endif
        write(coll_scatterUnit,'(1x,i2,2x,i4,2x,i5,2x,i1,3(2x,e18.10))') icoll,iturn,ipart,nabs_tmp,(p-pBef)/pBef,xp-xpBef,zp-zpBef
  #ifdef HDF5
        endif
  #endif
      end if
      return
    end if
  
  !++  Otherwise do multi-coulomb scattering.
  !++  REGULAR STEP IN ITERATION LOOP
    call k2coll_mcs(s)
  
  !++  Check if particle is outside of collimator (X.LT.0) after
  !++  MCS. If yes, calculate output longitudinal position (s),
  !++  reduce momentum (output as dpop) and return.
  !++  PARTICLE LEFT COLLIMATOR BEFORE ITS END.
    if(x.le.zero) then
      s=(zlm-rlen)+s
  
  #ifndef MERLINSCATTER
      call k2coll_calcIonLoss(mat,p,rlen,m_dpodx)
      p=p-m_dpodx*s
  #endif
  #ifdef MERLINSCATTER
      call merlinscatter_calc_ion_loss(p,edens(mat),pleng(mat),exenergy(mat),s,m_dpodx)
      p=p-m_dpodx
  #endif
      dpop=(p-p0)/p0
  
      if(dowrite_impact) then
        ! write coll_scatter.dat for complete scattering histories
  
negative.code: #endif
      write(coll_scatterUnit,'(1x,i2,2x,i4,2x,i5,2x,i1,3(2x,e14.6))') icoll,iturn,ipart,nabs_tmp,-one,zero,zero
  #ifdef HDF5
      end if
  #endif
      return
    end if
  
    ! Now treat the other types of interaction, as determined by ICHOIX:
  
    ! Nuclear-Elastic:          inter = 2
    ! pp Elastic:               inter = 3
    ! Single-Diffractive:       inter = 4    (changes momentum p)
    ! Coulomb:                  inter = 5
  
    ! As the single-diffractive interaction changes the momentum, save input momentum in p1.
    p1 = p
  
    ! Gettran returns some monte carlo number, that, as I believe, gives the rms transverse momentum transfer.
    t = k2coll_gettran(inter,mat,p)
  
    ! Tetat calculates from the rms transverse momentum transfer in
    ! monte-carlo fashion the angle changes for x and z planes. The
    ! angle change is proportional to SQRT(t) and 1/p, as expected.
    call k2coll_tetat(t,p,dxp,dzp)
  
    ! Apply angle changes
    xp = xp+dxp
    zp = zp+dzp
  
    ! Treat single-diffractive scattering.
    if(inter == 4) then
  
      ! added update for s
      s    = (zlm-rlen)+zlm1
      xpsd = dxp
      zpsd = dzp
      psd  = p1
  
      ! Add this code to get the momentum transfer also in the calling routine
      dpop = (p-p0)/p0
    end if
  
    if(dowrite_impact) then
  