We present an adaptive algorithm for bivariate interpolation of irregularly distributed points. The adaptive scheme is based on the radial basis function (RBF) partition of unity method, which is particularly suitable for the interpolation of a large number of points. By this technique we can decompose the entire domain into several subdomains of variable size, so that the algorithm can be applied successfully on highly varying point distributions. The use of efficient searching procedures together with a maximum likelihood estimation criterion enables us a fast and accurate computation of the RBF partition of unity interpolant. Numerical results show good performance of this adaptive algorithm.