anchor.diff:      type(StructuredMesh), pointer :: m
       type(Particle),       pointer :: p
   
  -    integer :: i, j
  -    integer :: ijk0(3)
  -    integer :: ijk1(3)
  -    integer :: ijkd_cross(4) ! (i,j,k) and direction of crossing
  -    integer :: n_cross
  -    real(8) :: uvw(3)
  -    real(8) :: xyz0(3)
  -    real(8) :: xyz1(3)
  -    real(8) :: xyz_cross(3)
  -    real(8) :: d(3)
  -    real(8) :: distance
  -    logical :: start_in_mesh
  -    logical :: end_in_mesh
  -    logical :: x_same
  -    logical :: y_same
  -    logical :: z_same
  +    integer :: i             ! loop indices
  +    integer :: j             ! loop indices
  +    integer :: ijk0(3)       ! indices of starting coordinates
  +    integer :: ijk1(3)       ! indices of ending coordinates
  +    integer :: n_cross       ! number of surface crossings
  +    integer :: surface       ! surface/direction of crossing, e.g. IN_RIGHT
  +    real(8) :: uvw(3)        ! cosine of angle of particle
  +    real(8) :: xyz0(3)       ! starting/intermediate coordinates
  +    real(8) :: xyz1(3)       ! ending coordinates of particle
  +    real(8) :: xyz_cross(3)  ! coordinates of bounding surfaces
  +    real(8) :: d(3)          ! distance to each bounding surface
  +    real(8) :: distance      ! actual distance traveled
  +    logical :: start_in_mesh ! particle's starting xyz in mesh?
  +    logical :: end_in_mesh   ! particle's ending xyz in mesh?
  +    logical :: x_same        ! same starting/ending x index (i)
  +    logical :: y_same        ! same starting/ending y index (j)
  +    logical :: z_same        ! same starting/ending z index (k)
   
       ! Copy starting and ending location of particle
       xyz0 = p % last_xyz
  
positive.code:        if (uvw(1) > 0) then
            do i = ijk0(1), ijk1(1) - 1
               ijk0(1) = i
               print *, ijk0, "OUT_TOP"
            end do
         else
            do i = ijk0(1) - 1, ijk1(1), -1
               ijk0(1) = i
               print *, ijk0, "IN_TOP"
            end do
         end if
         return
  
negative.code:              ! current in (i-1,j,k)
               ijk0(1) = ijk0(1) - 1
               xyz_cross(1) = xyz_cross(1) - m % width(1)
               if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) print *, ijk0, "IN_RIGHT"
  
            end if
         elseif (distance == d(2)) then
            if (uvw(2) > 0) then
               ! Crossing into front mesh cell -- this is treated as outgoing
               ! current in (i,j,k)
               if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) print *, ijk0, "OUT_FRONT"
               ijk0(2) = ijk0(2) + 1
               xyz_cross(2) = xyz_cross(2) + m % width(2)
            else
  