In this paper, we investigate a numerical solution of the problem of viscous fluid convective flow in a twodimensional region with a free surface. The presence of a free surface poses a big problem for the numerical simulation since the surface position is unknown. It is necessary to determine the interface position both for the non-stationary problem and for the stationary case. The various numerical methods are used nowadays for solving such problems. We use the finite difference method and formulate the problem in terms of eddycurrent variables. This approach has several advantages when we consider a two-dimensional flow because, in this case, the continuity equation is fulfilled exactly. This is especially important when there is no outflow or inflow of liquid in the fluid flow domain. However, in the case of stream function and vorticity the form of the boundary conditions on the free surface is much more complicated. It is difficult to construct the efficient numerical algorithms. We propose the method, in which the free boundary is found by solving the boundary value problem for the ordinary differential equations of third order with nonlocal condition. The test calculations are carried out. The results demonstrate the effectiveness of the proposed method.