The implementation of the free surface boundary condition into a composite spectral Chebyshev-Fourier method for solving the three dimensional equations of motion is presented. In this method spatial derivatives with respect to the vertical direction are calculated with a Chebyshev spectral method, while differencing with respect to the horizontal directions is performed with the Fourier method. The technique can handle free surface boundary conditions in a more rigorous manner than the ordinary Fourier method. It therefore appears well-suited for realistic full wave modeling, in particular of near surface layer problems. Isotropic media and transversely isotropic media with a vertical axis of symmetry are considered. A comparison of the isotropic modeling results with the analytic solution for Lamb's problem shows the high accuracy of the alogorithm. Modeling examples are presented for a 3D thin layer overlaying an elastic halfspace and a 3D transversely isotropic halfspace.