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:   ! Wind
    double precision, dimension(:,:),   allocatable :: wind_x
    double precision, dimension(:,:),   allocatable :: wind_y
  
    allocate(dhdt(0:nx+1, 0:ny+1, layers))
    allocate(dhdtold(0:nx+1, 0:ny+1, layers))
  
negative.code: 
    integer, dimension(:,:), allocatable :: ilower, iupper
    integer, dimension(:,:), allocatable :: jlower, jupper
    integer :: i, j
  
    integer*8 :: parcsr_A
    integer*8 :: par_b
    integer*8 :: par_x
    integer*8 :: solver
    integer*8 :: precond
    integer*8 :: hypre_grid
    integer   :: offsets(2,5)
    ! Set default values here
  
    ! TODO Possibly wait until the model is split into multiple files,
  