A method is presented for the accurate and efficient calculation of the dynamic Green's functions for a layered viscoelastic solid under plane strain conditions. The method is based on the classical wavenumber integral representation and a delta matrix formulation of the elastodynamic field. The high efficiency is achieved through the introduction of a new quadrature scheme in which the kernels of the wavenumber integrals are represented by means of polynomials in finite and semi-infinite panels, and the resulting oscillatory integrals are evaluated analytically. The accuracy of the results are controlled through the use of an adaptive procedure whereby the number of panels are increased successively until the desired accuracy is reached, with no previous function evaluations wasted. The computer program currently runs on IBM PCs for all multi-layered structures at all finite frequencies.