anchor.diff:    print *, 'iupper (x) = ', iupper(:,1)
     print *, 'iupper (y) = ', iupper(:,2)
   end if
  -    call HYPRE_IJMatrixCreate(mpi_comm, ilower(myid,1), & 
  -           iupper(myid,1), ilower(myid,2), iupper(myid,2), A, ierr)
  -
  -  ! ! Choose a parallel csr format storage (apparently his is the only one supported)
  -  call HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR, ierr)
  -
  -  ! ! Initialize before setting coefficients
  -  call HYPRE_IJMatrixInitialize(A, ierr)
  +  !   call HYPRE_IJMatrixCreate(mpi_comm, ilower(myid,1), & 
  +  !          iupper(myid,1), ilower(myid,2), iupper(myid,2), A, ierr)
  +
  +  ! ! ! Choose a parallel csr format storage (apparently his is the only one supported)
  +  ! call HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR, ierr)
  +
  +  ! ! ! Initialize before setting coefficients
  +  ! call HYPRE_IJMatrixInitialize(A, ierr)
  +
  +  call Hypre_StructGridCreate(MPI_COMM_WORLD, 2, hypre_grid, ierr)
  +
  +  do i = 0, num_procs-1
  +    call HYPRE_StructGridSetExtents(hypre_grid, ilower(i,:),iupper(i,:), ierr)
  +  end do
  +
  +  call HYPRE_StructGridAssemble(hypre_grid, ierr)
   #endif
   
   
  
positive.code:   double precision :: pi
    integer :: nwrite, avwrite
    double precision :: rjac
  
    ! Time step loop variable
    integer :: n
  
negative.code:   allocate(iupper(0:num_procs-1, 2))
  
  
    do i = 0, nProcX - 1
      ilower(i * nProcY:(i+1)*nProcY - 1,1) = i * nx / nProcX
      iupper(i * nProcY:(i+1)*nProcY - 1,1) = ((i+1) * nx / nProcX) - 1
    end do
    ! correct first ilower value to exclude the global halo
    ilower(0,1) = 1
  
    do j = 0, nProcY - 1
      ilower(j * nProcX:(j+1)*nProcX - 1,2) = j * ny / nProcY
      iupper(j * nProcX:(j+1)*nProcX - 1,2) = ((j+1) * ny / nProcY) - 1
    end do
    ! correct first ilower value to exclude the global halo
    ilower(0,2) = 1
  
  if (myid .eq. 0) then
    print *, 'ilower (x) = ', ilower(:,1)
  