anchor.diff:    ! add free surface elevation to the upper layer
     h_temp(:, :, 1) = h(:, :, 1) + eta*freesurfFac
   
  -  do i = 1, nx
  +  do i = 1, nx+1
       do j = 1, ny
         do k = 1, layers
           ub(i,j) = ub(i,j) + u(i,j,k)*(h_temp(i,j,k)+h_temp(i-1,j,k))/2d0
         end do
       end do
     end do
  -
  -  call wrap_fields_2D(ub, nx, ny)
   
     return
   end subroutine calc_baro_u
  
positive.code:   double precision eta(0:nx+1, 0:ny+1)
    double precision freesurfFac
    double precision v(0:nx+1, 0:ny+1, layers)
    double precision vb(0:nx+1, 0:ny+1)
  
    vb = 0d0
  
  
negative.code:   integer i, j
    double precision eta(0:nx+1, 0:ny+1)
    double precision etastar(0:nx+1, 0:ny+1)
    double precision ub(nx+1, ny)
    double precision vb(nx, ny+1)
    double precision freesurfFac, dx, dy, dt
  
    etastar = 0d0
  