anchor.diff:              obj % ntot = obj % ntot * numCells(j)
           enddo
           
  -        allocate(obj % dimProd, size = obj % ndims)
  -        
  -        obj % dimProd(ndims) = 1
  +        allocate(obj % dimProd(obj % ndims))
  +        
  +        obj % dimProd(obj % ndims) = 1
           do j = obj % ndims - 1, 1, -1
               obj % dimProd(j) =  obj % dimProd(j + 1) * obj % numCells(j + 1)
           enddo
           
  -        allocate(obj % f, size = obj % ntot)
  +        allocate(obj % f(obj % ntot))
           
           ! initialize
           obj % f = 0
  
positive.code:         integer, intent(out) :: res(:)
      
          integer :: i
          do i = 1:obj % ndims
              res(i) = mod(div((flatIndex - 1), obj % dimProd[i]), obj % numCells[i]) + 1
          enddo
      end subroutine
      
  
negative.code:         integer iunit, i
          
          open(newunit = iunit, file = filename, status = 'new')
          write(iunit, *) '# vtk DataFile Version 2.0'
          write(iunit, *) 'upwind.cxx'
          write(iunit, *) 'ASCII'
          write(iunit, *) 'DATASET RECTILINEAR_GRID'
          write(iunit, *) 'DIMENSIONS'
          ! in VTK the first dimension varies fastest so need 
          ! to invert the order of the dimensions
          if (obj % ndims > 2) then
              write(iunit, *) ' ', obj % numCells(3) + 1
          else
              write(iunit, *) ' 1'
          endif 
          
          if (obj % ndims > 1) then
              write(iunit, *) ' ', obj % numCells(2) + 1
          else
              write(iunit, *) ' 1'
          endif
          
          write(iunit, *) ' ', obj % numCells(1) + 1
          write(iunit, *) 'X_COORDINATES '
          if (obj % ndims > 2) then
              write(iunit, *) obj % numCells(3) + 1, ' double'
              do i = 1, obj % numCells(3)
                  write(iunit, *) ' ', 0.0 + obj % deltas(3) * i
              enddo     
    
          else 
              write(iunit, *) '1 double'
              write(iunit, *) '0.0'
          endif
      
          write(iunit, *) 'Y_COORDINATES '
          if (obj % ndims > 1) then
              write(iunit, *) obj % numCells(2) + 1, ' double'
              do i = 1, obj % numCells(2)
                  write(iunit, *) ' ', 0.0 + obj % deltas(2) * i
              enddo      
          else
              write(iunit, *) '1 double'
              write(iunit, *) '0.0'
          endif
        
          write(iunit, *) 'Z_COORDINATES '
          write(iunit, *) obj % numCells(1) + 1, ' double'
          do i = 1, obj % numCells(1)
              write(iunit, *) ' ', 0.0 + obj % deltas(1) * i
          enddo
          write(iunit, *) 'CELL_DATA ', obj % ntot
          write(iunit, *) 'SCALARS f double 1'
          write(iunit, *) 'LOOKUP_TABLE default'
          do i = 1, obj % ntot
              write(iunit, *) obj % f(i), ' '
          enddo
      
          close(iunit)
      
      end subroutine
  
      subroutine upwind_print(obj)
  