In the past, linear quasi-geostrophic theory has proven successful in modeling the vertical and meridional propagation of stationary planetary waves in the stratosphere. Since in such models the wave solution does not sensitively depend on the wave damping, the latter was usually implemented as relaxation with a simple damping coefficient. This is likely to underestimate the amount of damping in regions of wave breaking. In the present study, a recently developed parameterization for Rossby wave breaking (Garcia, 1991) is applied to obtain an improved representation of wave damping in the context of a linear model. As a result the divergence of the Eliassen-Palm flux calculated from the wave solution represents a meaningful diagnostic for how strongly the mean flow is affected by the waves, which is generally not the case in the standard Matsuno-model. A residual circulation is computed with quantifies the non-local effect of the wave forcing on long time scales. For illustration the new method is applied to different basic states representing the different phases of an idealized QBO. The divergence of the EP flux and the strength of the diagnosed residual circulation sensitively depend on small differences in the basic state. The dependences are qualitatively in agreement with the results from previous investigations. The current model allows the attribution of the differences in wave mean flow interaction to differences in the basic flow. Furthermore, it facilitates the distinction between wave propagation, which affects wave amplitude, and wave dissipation, which feeds back on the mean flow.