anchor.diff:  !! Euler timestepping for the free surface anomaly, or for the surface
   !! pressure required to keep the barotropic flow nondivergent.
   
  -subroutine SOR_solver(a, etanew, etastar, freesurfFac, nx, ny, dt, rjac, eps, maxits, n)
  +subroutine SOR_solver(a, etanew, etastar, freesurfFac, nx, ny, dt, &
  +    rjac, eps, maxits, n)
     implicit none
   
     double precision a(5, nx, ny)
  
positive.code: 
    ! Now calculate d/dt of u, v, h and store as dhdtveryold, dudtveryold
    ! and dvdtveryold
    call evaluate_dhdt(dhdtveryold, hhalf, uhalf, vhalf, ah, dx, dy, nx, ny, layers, &
        spongeHTimeScale, spongeH, wetmask, RedGrav)
  
    call evaluate_dudt(dudtveryold, hhalf, uhalf, vhalf, b, zeta, &
        wind_x, fu, au, ar, slip, dx, dy, hfacN, hfacS, nx, ny, layers, rho0, &
  
negative.code:   call evaluate_zeta(zeta, uhalf, vhalf, nx, ny, layers, dx, dy)
  
    ! Now calculate d/dt of u, v, h and store as dhdtold, dudtold and dvdtold
    call evaluate_dhdt(dhdtold, hhalf, uhalf, vhalf, ah, dx, dy, &
        nx, ny, layers, spongeHTimeScale, spongeH, wetmask, RedGrav)
    call evaluate_dudt(dudtold, hhalf, uhalf, vhalf, b, zeta, wind_x, &
        fu, au, ar, slip, dx, dy, hfacN, hfacS, nx, ny, layers, rho0, &
        spongeUTimeScale, spongeU, RedGrav, botDrag)
  