A theoretical model is presented describing harmonic heat flow in a two layer system heated by a modulated Gaussian laser beam. Amplitude and phase of the modulated temperature rise in the layers, as well as in the backing substrate and adjacent atmosphere, are calculated by solving the three-dimensional heat conduction equation with a source term including exponential absorption of the laser light in one or two layers. Heat conduction is assumed to be isotropic throughout the system, however, a thermal contact resistance between the two layers can be taken into account. Results are presented for single and double layer systems of gold and various dielectric thin film materials on glass substrates. Amplitude and phase of the harmonic temperature variation are calculated either as a function of position in the sample system or at the surface as a function of the laser beam modulation frequency. It is found that both amplitude and phase of the calculated temperature rise exhibit typical thin film features in their frequency dependence, however, the phase is more sensitive to thin film phenomena than the amplitude. The phase shows typical extrema in that frequency region, where the thermal diffusion length in the film is equal to the film thickness. Based on these findings, we demonstrate how these calculations can be utilized for the interpretation of thin film thermal parameter measurements. The influence of thermal wave interference is demonstrated, and it is discussed how the main thermal parameters like conductivity, effusivity, and thermal contact resistance of the thin film system can be extracted from measurements by a fit of theoretical curves to experimental data. Applying a simple one-dimensional thermal expansion model, surface displacements for thin film systems are calculated and the applicability of photothermal surface displacement for thin film conductivity measurements is discussed.