anchor.diff:    double precision :: vb(nx, ny+1)
     double precision :: etastar(0:nx+1, 0:ny+1)
   
  +  ! External solver variables
  +  integer   :: offsets(2,5)
  +  integer*8 :: stencil
  +  integer :: ierr
  +  integer :: i, j ! loop variables
  +  integer*8 :: hypre_grid
  +  integer*8 :: hypre_A
  +  integer*8 :: hypre_b
  +  integer*8 :: hypre_x
  +  integer   ::  nentries, nvalues, stencil_indices(5)
  +  double precision, dimension(:), allocatable :: values
  +
  +  integer :: MPI_COMM_WORLD
  +  integer :: num_procs
  +  integer :: ilower(0:num_procs-1,2), iupper(0:num_procs-1,2)
  +
     ! Calculate the barotropic velocities
     call calc_baro_u(ub, unew, hnew, eta, freesurfFac, nx, ny, layers)
     call calc_baro_v(vb, vnew, hnew, eta, freesurfFac, nx, ny, layers)
  
positive.code: subroutine barotropic_correction(hnew, unew, vnew, eta, etanew, depth, a, &
      dx, dy, wetmask, hfacW, hfacS, dt, &
      maxits, eps, rjac, freesurfFac, thickness_error, &
      g_vec, nx, ny, layers, n)
    implicit none
  
    double precision, intent(inout) :: hnew(0:nx+1, 0:ny+1, layers)
  
negative.code:   ! Prevent barotropic signals from bouncing around outside the
    ! wet region of the model.
    ! etastar = etastar*wetmask
  #ifndef useExtSolver
    call SOR_solver(a, etanew, etastar, freesurfFac, nx, ny, &
        dt, rjac, eps, maxits, n)
    ! print *, maxval(abs(etanew))
  #elseifdef useExtSolver
    call HYPRE_StructVectorCreate(MPI_COMM_WORLD, hypre_grid, hypre_b, ierr)
  #endif
  
    etanew = etanew*wetmask
    call wrap_fields_2D(etanew, nx, ny)
  