anchor.diff:    call HYPRE_ParCSRPCGSolve(hypre_solver, hypre_A, hypre_b, &
                               hypre_x, ierr)
   
  -  ! code for printing out results form the external solver
  +  ! code for printing out results from the external solver
  +  ! Not being used, but left here since the manual isn't very helpful
  +  ! and this may be useful in the future.
     ! call HYPRE_ParCSRPCGGetNumIterations(hypre_solver, & 
  -  !   temp(1), ierr)
  -  ! print *, 'num iterations = ', temp(1)
  +  !   hypre_out(1), ierr)
  +  ! print *, 'num iterations = ', hypre_out(1)
   
     ! call HYPRE_ParCSRPCGGetFinalRelative(hypre_solver, & 
  -  !   temp(2), ierr)
  -  ! print *, 'final residual norm = ', temp(2)
  +  !   hypre_out(2), ierr)
  +  ! print *, 'final residual norm = ', hypre_out(2)
   
     do i = 0, num_procs-1
       call HYPRE_StructVectorGetBoxValues(hypre_x, & 
         ilower(i,:), iupper(i,:), values, ierr)
     end do
   
  -  ! print *, 'values = ', values
     do i = 1, nx ! loop over every grid point
       do j = 1, ny
     ! the 2D array is being laid out like
     ! [x1y1, x1y2, x1y3, x2y1, x2y2, x2y3, x3y1, x3y2, x3y3]
       etanew(i,j) = values( ((i-1)*ny + j) )
  -    ! print *, etanew(i,j)
  -    end do
  -  end do
  -
  -  ! print *, 'etanew = ', etanew
  -
  -  ! debuggin commands from hypre library - dump out a single copy of these two variables. Can be used to check that the values have been properly allocated.
  -  call HYPRE_StructVectorPrint(hypre_x, ierr)
  +    end do
  +  end do
  +
  +  ! debugging commands from hypre library - dump out a single 
  +  ! copy of these two variables. Can be used to check that the 
  +  ! values have been properly allocated.
  +  ! call HYPRE_StructVectorPrint(hypre_x, ierr)
     ! call HYPRE_StructMatrixPrint(hypre_A, ierr)
   
  -  ! do i = 1, nx !0, num_procs-1
  -  !   do j = 1, ny
  -  !   temp(1) = i
  -  !   temp(2) = j
  -  !   ! call HYPRE_StructVectorSetBoxValues(hypre_b, & 
  -  !   !   ilower(i,:), iupper(i,:), values, ierr)
  -
  -  !   call HYPRE_StructVectorGetValues(hypre_x, & 
  -  !     temp, etanew(i,j), ierr)
  -  ! end do
  -  ! end do
  -  !   call HYPRE_StructVectorSetBoxValues(hypre_x, & 
  -  !     ilower(myid,:), iupper(myid,:), values, ierr)
  -
  -  ! print *, 'values from hypre_x', values  
  -  ! print *, 'etastar = ', etastar
  -  ! print *, 'etanew = ', etanew
  -  
     call HYPRE_StructPCGDestroy(hypre_solver, ierr)
   
  -  !   do i = 1, nx !0, num_procs-1
  -  !   do j = 1, ny
  -  !   temp(1) = i
  -  !   temp(2) = j
  -  !   ! call HYPRE_StructVectorSetBoxValues(hypre_b, & 
  -  !   !   ilower(i,:), iupper(i,:), values, ierr)
  -
  -  !   call HYPRE_StructVectorGetValues(hypre_x, & 
  -  !     temp, etastar(i,j), ierr)
  -  ! end do
  -  ! end do
  -
  -  ! do i = 0, num_procs-1
  -  !   call hypre_StructVectorGetBoxValues(hypre_x, & 
  -  !     ilower(i,:), iupper(i,:), values, ierr)
  -  ! end do
  -  ! ! possibly this is an alternative routine for getting the values? Manual is not very helpful
  -
  -
  -  ! ! Move the values back into the 2D array etanew
  -  ! if (myid .eq. 0) then
  -  ! do i = 1, nx ! loop over every grid point
  -  !   do j = 1, ny
  -  ! ! the 2D array is laid out like
  -  ! ! [x1y1, x1y2, x1y3, x2y1, x2y2, x2y3, x3y1, x3y2, x3y3]
  -  !     etanew(i,j) =  values( ((i-1)*ny + j) )
  -  !   end do
  -  ! end do
  -  ! end if
  -  
  -
  -
  -
  -
   #endif
   
  -
  -    ! print *, 'here'
     etanew = etanew*wetmask
  -  ! print *, 'etanew after wetmask', etanew
  +
     call wrap_fields_2D(etanew, nx, ny)
   
     ! Now update the velocities using the barotropic tendency due to
  
positive.code: 
      call HYPRE_StructMatrixInitialize(hypre_A, ierr)
  
      ! nentries = 5 ! a five point stncil means five entries for each grid point
      ! nvalues = nx * ny * 5 ! total number of grid points to simulate * five stencil entries per grid point
      ! allocate(values(nvalues))
  
      ! do i = 1, nx !values, nentries ! loop over every grid point, skipping by the number of points in the stencil. This 
      ! ! selects only the central points of the stencil
      !   do j = 1, ny
  
  
      ! ! the 2D array is being laid out like
      ! ! [x1y1, x1y2, x1y3, x2y1, x2y2, x2y3, x3y1, x3y2, x3y3]
      !   values( ((i-1)*ny + j)*5 )    = a(5,i,j)
      !   values( ((i-1)*ny + j)*5 + 1) = a(3,i,j)
      !   values( ((i-1)*ny + j)*5 + 2) = a(1,i,j)
      !   values( ((i-1)*ny + j)*5 + 3) = a(4,i,j)
      !   values( ((i-1)*ny + j)*5 + 4) = a(2,i,j)
  
      !   end do
      ! end do
  
      ! do i = 0, num_procs-1
      !   call HYPRE_StructMatrixSetBoxValues(hypre_A, & 
      !     ilower(i,:), iupper(i,:), 5, stencil_indices, & 
      !     values, ierr)
      ! end do
  
    do i = 1, nx !0, num_procs-1
      do j = 1, ny
      temp(1) = i
      temp(2) = j
  
      
      call HYPRE_StructMatrixSetValues(hypre_A, & 
          temp, 1, 0, & 
          a(5,i,j) - freesurfFac/dt**2, ierr)
      call HYPRE_StructMatrixSetValues(hypre_A, & 
          temp, 1, 1, & 
          a(3,i,j), ierr)
      call HYPRE_StructMatrixSetValues(hypre_A, & 
          temp, 1, 2, & 
          a(1,i,j), ierr)
      call HYPRE_StructMatrixSetValues(hypre_A, & 
          temp, 1, 3, & 
          a(4,i,j), ierr)
      call HYPRE_StructMatrixSetValues(hypre_A, & 
          temp, 1, 4, & 
          a(2,i,j), ierr)
  
    end do
    end do
  
  
  
      call HYPRE_StructMatrixAssemble(hypre_A, ierr)
  
        call MPI_Barrier(  MPI_COMM_WORLD, ierr)
  
  
  #endif
  
  
negative.code: 
  #ifdef useExtSolver
  
    ! Create the rhs vector, b
    call HYPRE_StructVectorCreate(MPI_COMM_WORLD, hypre_grid, hypre_b, ierr)
    call HYPRE_StructVectorInitialize(hypre_b, ierr)
  
    ! set rhs values (vector b)
    do i = 1, nx ! loop over every grid point
  