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:   allocate(dudtveryold(0:nx+1, 0:ny+1, layers))
    allocate(unew(0:nx+1, 0:ny+1, layers))
    allocate(dudt_bt(nx, 0:ny))
    allocate(ub(0:nx+1, 0:ny+1))
    allocate(uhalf(0:nx+1, 0:ny+1, layers))
    allocate(uav(0:nx+1, 0:ny+1, layers))
  
  
negative.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(nx, ny+1)
  
    vb = 0d0
  
  