anchor.diff:  
     ! norm_old = norm0
   
  -  do nit = 1,maxits
  +  do nit = 1, maxits
       ! norm_old = norm
       norm=0.d0
  -    do i=1,nx
  -      do j=1,ny
  +    do i=1, nx
  +      do j=1, ny
           res(i,j) = &
               a(1,i,j)*etanew(i+1,j) &
               + a(2,i,j)*etanew(i,j+1) &
  
positive.code: ! -----------------------------------------------------------------------------
  !> Evaluate relative vorticity at lower left grid boundary (du/dy
  !! and dv/dx are at lower left corner as well)
  subroutine evaluate_zeta(zeta,u,v,nx,ny,layers,dx,dy)
    implicit none
  
    integer nx,ny,layers
    integer i,j,k
    double precision h(0:nx,0:ny,layers)
    double precision u(0:nx+1,0:ny+1,layers)
    double precision v(0:nx+1,0:ny+1,layers)
    double precision zeta(0:nx+1,0:ny+1,layers)
    double precision dx, dy
  
    zeta = 0d0
  
    do k = 1,layers
      do j=1,ny+1
        do i=1,nx+1
          zeta(i,j,k)=(v(i,j,k)-v(i-1,j,k))/dx-(u(i,j,k)-u(i,j-1,k))/dy
        end do
      end do
  
negative.code:   ! - hfacW and hfacS are zero where the transition between
    !   wet and dry cells occurs.
    ! - wetmask is 1 in wet cells, and zero in dry cells.
    do k = 1, layers
      u(:, :, k) = u(:, :, k) * hfacW * wetmask(:, :)
      v(:, :, k) = v(:, :, k) * hfacS * wetmask(:, :)
    end do
  
    ! Wrap fields around for periodic simulations
    call wrap_fields_3D(u, nx, ny, layers)
    call wrap_fields_3D(v, nx, ny, layers)
    call wrap_fields_3D(h, nx, ny, layers)
  
    ! Now the model is ready to start.
    ! - We have h, u, v at the zeroth time step, and the tendencies at
  