anchor.diff:      ! Find number of cells
       geo%n_cells = size(x)
       geo%n_masked = geo%n_cells
  -
  -    print *,geo%n_cells
   
       ! Allocate cells
       allocate(geo%cells(geo%n_cells))
  
positive.code:     integer :: i, n_neighbors
      type(vector3d_dp) :: r0, ri, rmid, n
      real(dp) :: t, tx, ty, tz
  
      icell => geo%cells(p%icell%ic)
  
  
negative.code:   type(grid_geometry_desc),public,target :: geo
  
    integer :: n_filled = 0
  
    ! Declare here otherwise pointer becomes corrupt
    real(dp),allocatable, target :: points(:,:)
  
  contains
  
  