anchor.diff:      call wrap_fields_2D(etanew, nx, ny)
   
       if (nit.gt.1.and.norm.lt.eps*norm0) then
  -      ! print 10014,  eps,  nit
  -      ! 10014   format('Res less than eps of original. eps, iteration: ', F5.4, I7)
  +
         return
  -      ! elseif (nit.gt.1.and.abs(norm - norm_old).lt.norm_old*eps) then
  -      !     print *, 'change in residual less than eps of previous residual. iteration:', eps, nit
  -      !     return
  -      ! This was not a good way of exiting the solver - it would occasionally leave after 2 or 3 iterations.
  +
       end if
     end do
   
  
positive.code:   double precision rjac, eps
    double precision rhs(nx, ny)
    double precision res(nx, ny)
    double precision norm, norm0, norm_old
    double precision relax_param
    character(30) Format
  
    rhs = -etastar(1:nx,1:ny)/dt**2
    ! first guess for etanew
  
negative.code:   double precision :: H0 ! default depth in no file specified
    ! Pressure solver variables
    double precision :: a(5, nx, ny)
  
    ! Bernoulli potential and relative vorticity
    double precision :: b(0:nx+1, 0:ny+1, layers)
  