This paper presents an algorithm for fitting a bivariate spline function to a set of scattered data. The number of knots and their position are determined automatically. However, the user has to provide a non-negative constant to control the tradeoff between closeness of fit and smoothness of fit.