We used the Green's function method to describe the multiple scattering of surface waves generated by interactions with the complex 3-D earth structures. The basic process involves an integral equation of convolutional type. An efficient multilevel fast multipole method is utilized to accelerate the matrix-vector products that correspond to the integrals on the horizontal plane. This new algorithm is shown to be quite successful when compared to its more traditional complement (e.g. direct integration method, DIM) with particularly the large number of data points. The fast execution is achieved through well-organized truncated multipole expansions, functions grouping and translation operators. In general, the new algorithm has a fruitful logarithmic time complexity as opposed to an uncomfortable exponential time complexity attainable with the traditional algorithms. This algorithm requires the user provide some important parameters to operate. One of them is the truncation number of the infinite series associated with the multipole expansions for which two linear relationships are derived for an automatic determination. The other central parameter is the clustering number that is used to group data points in the data structure. We showed that the clustering number around 5 is mostly an optimum value providing a minimum run time. However, greater clustering numbers (i.e. similar to 10) become necessary for an optimum operation when the number of data points gets real large.