From 0431b19504527ca5e662f4e5e87d8bc6f4b9864a Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Fri, 16 Dec 2011 10:28:56 -0600
Subject: [PATCH] Added out-of-bounds check on lattice elements in find_cell.

---
 src/geometry.F90 | 88 ++++++++++++++++++++++++++++--------------------
 1 file changed, 52 insertions(+), 36 deletions(-)

diff --git a/src/geometry.F90 b/src/geometry.F90
index 3ac71e768c..29ad60626c 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -171,26 +171,42 @@ contains
              x = ceiling((xyz(1) - lat % x0)/lat % width_x)
              y = ceiling((xyz(2) - lat % y0)/lat % width_y)
 
-             ! Create new level of coordinates
-             p % in_lower_universe = .true.
-             allocate(p % coord % next)
-             
-             ! adjust local position of particle
-             p % coord % next % xyz(1) = p % coord % xyz(1) - &
-                  (lat%x0 + (x-0.5_8)*lat%width_x)
-             p % coord % next % xyz(2) = p % coord % xyz(2) - &
-                  (lat%y0 + (y-0.5_8)*lat%width_y)
-             p % coord % next % xyz(3) = p % coord % xyz(3)
-             p % coord % next % uvw    = p % coord % uvw
+             ! Check if lattice coordinates are within bounds
+             if (x < 1 .or. x > lat % n_x .or. &
+                  y < 1 .or. y > lat % n_y) then
 
-             ! Move particle to next level
-             p % coord => p % coord % next
+                ! This condition should only get hit in rare circumstances where
+                ! a neutron hits the corner of a lattice. In this case, the
+                ! neutron may need to be moved diagonally across the lattice. To
+                ! do so, we remove all lower coordinate levels and then search
+                ! from universe 0.
+
+                p % coord => p % coord0
+                call deallocate_coord(p % coord % next)
+
+             else
+
+                ! Create new level of coordinates
+                p % in_lower_universe = .true.
+                allocate(p % coord % next)
              
-             ! set particle lattice indices
-             p % coord % lattice   = c % fill
-             p % coord % lattice_x = x
-             p % coord % lattice_y = y
-             p % coord % universe  = lat % element(x,y)
+                ! adjust local position of particle
+                p % coord % next % xyz(1) = p % coord % xyz(1) - &
+                     (lat%x0 + (x-0.5_8)*lat%width_x)
+                p % coord % next % xyz(2) = p % coord % xyz(2) - &
+                     (lat%y0 + (y-0.5_8)*lat%width_y)
+                p % coord % next % xyz(3) = p % coord % xyz(3)
+                p % coord % next % uvw    = p % coord % uvw
+
+                ! Move particle to next level
+                p % coord => p % coord % next
+             
+                ! set particle lattice indices
+                p % coord % lattice   = c % fill
+                p % coord % lattice_x = x
+                p % coord % lattice_y = y
+                p % coord % universe  = lat % element(x,y)
+             end if
 
              call find_cell(p, found)
              if (.not. found) exit
@@ -216,19 +232,19 @@ contains
     type(Particle), pointer :: p
     integer, intent(in)     :: last_cell  ! last cell particle was in
 
-    real(8)                 :: x          ! x-x0 for sphere
-    real(8)                 :: y          ! y-y0 for sphere
-    real(8)                 :: z          ! z-z0 for sphere
-    real(8)                 :: R          ! radius of sphere
-    real(8)                 :: u          ! x-component of direction
-    real(8)                 :: v          ! y-component of direction
-    real(8)                 :: w          ! z-component of direction
-    real(8)                 :: n1         ! x-component of surface normal
-    real(8)                 :: n2         ! y-component of surface normal
-    real(8)                 :: n3         ! z-component of surface normal
-    real(8)                 :: dot_prod   ! dot product of direction and normal
-    real(8)                 :: norm       ! "norm" of surface normal
-    logical                 :: found      ! particle found in universe?
+    real(8) :: x        ! x-x0 for sphere
+    real(8) :: y        ! y-y0 for sphere
+    real(8) :: z        ! z-z0 for sphere
+    real(8) :: R        ! radius of sphere
+    real(8) :: u        ! x-component of direction
+    real(8) :: v        ! y-component of direction
+    real(8) :: w        ! z-component of direction
+    real(8) :: n1       ! x-component of surface normal
+    real(8) :: n2       ! y-component of surface normal
+    real(8) :: n3       ! z-component of surface normal
+    real(8) :: dot_prod ! dot product of direction and normal
+    real(8) :: norm     ! "norm" of surface normal
+    logical :: found    ! particle found in universe?
     type(Surface),  pointer :: surf => null()
 
     surf => surfaces(abs(p % surface))
@@ -429,11 +445,11 @@ contains
     type(Particle), pointer :: p
     integer, intent(in) :: lattice_crossed
 
-    integer        :: i_x      ! x index in lattice
-    integer        :: i_y      ! y index in lattice
-    real(8)        :: x0       ! half the width of lattice element
-    real(8)        :: y0       ! half the height of lattice element
-    logical        :: found    ! particle found in cell?
+    integer :: i_x   ! x index in lattice
+    integer :: i_y   ! y index in lattice
+    real(8) :: x0    ! half the width of lattice element
+    real(8) :: y0    ! half the height of lattice element
+    logical :: found ! particle found in cell?
     type(Lattice),  pointer :: lat => null()
 
     lat => lattices(p % coord % lattice)