anchor.diff:            + a(3,i,j)*etanew(i-1,j) &
             + a(4,i,j)*etanew(i,j-1) &
             + a(5,i,j)*etanew(i,j)   &
  -          - freesurfFac*etanew(i,j)/dt**2 &
             - rhs(i,j)
         norm0 = norm0 + abs(res(i,j))
         etanew(i,j) = etanew(i,j)-relax_param*res(i,j)/a(5,i,j)
  
positive.code:     a(2, i, ny) = 0.0
      a(4, i, 1) = 0.0
    end do
    do j = 1, ny
      do i = 1, nx
        a(5,i,j) = -a(1,i,j)-a(2,i,j)-a(3,i,j)-a(4,i,j)
      end do
    end do
  
  
negative.code:             + a(3,i,j)*etanew(i-1,j) &
              + a(4,i,j)*etanew(i,j-1) &
              + a(5,i,j)*etanew(i,j)   &
              - rhs(i,j)
          norm = norm + abs(res(i,j))
          etanew(i,j) = etanew(i,j)-relax_param*res(i,j)/(a(5,i,j))
  