anchor.diff:  
         ! Now update the velocities using the barotropic tendency due to
         ! the pressure gradient.
  -      call update_velocities_for_barotropic_tendency(unew, etanew, g_vec, 1, 0, dx, dt, nx, ny, layers)
  -      call update_velocities_for_barotropic_tendency(vnew, etanew, g_vec, 0, 1, dy, dt, nx, ny, layers)
  +      call update_velocities_for_barotropic_tendency(unew, etanew, g_vec, &
  +          1, 0, dx, dt, nx, ny, layers)
  +      call update_velocities_for_barotropic_tendency(vnew, etanew, g_vec, &
  +          0, 1, dy, dt, nx, ny, layers)
   
         ! We now have correct velocities at the next time step, but the
         ! layer thicknesses were updated using the velocities without
         ! the barotropic pressure contribution. Force consistency
         ! between layer thicknesses and ocean depth by scaling
         ! thicknesses to agree with free surface.
  -      call enforce_moderate_free_surface(hnew, etanew, depth, freesurfFac, thickness_error, nx, ny, layers)
  +      call enforce_moderate_free_surface(hnew, etanew, depth, &
  +          freesurfFac, thickness_error, nx, ny, layers)
   
         ! Apply the boundary conditions
         call apply_boundary_conditions(unew, hfacW, wetmask, nx, ny, layers)
  
positive.code:   return
  end subroutine evaluate_zeta
  
  ! -----------------------------------------------------------------------------
  !> Calculate the tendency of layer thickness for each of the active layers
  !! dh/dt is in the centre of each grid point.
  subroutine evaluate_dhdt(dhdt, h, u, v, ah, dx, dy, nx, ny, layers, &
  
negative.code:   return
  end subroutine calc_eta_star
  
  ! ---------------------------------------------------------------------------
  !> Use the successive over-relaxation algorithm to solve the backwards
  !! Euler timestepping for the free surface anomaly, or for the surface
  !! pressure required to keep the barotropic flow nondivergent.
  