A Bayesian Linear Dynamical Mode (LDM) decomposition method is applied to isolate robust modes of climate variability in the observed surface air temperature (SAT) field. This decomposition finds the optimal number of internal modes characterized by their own time scales, which enter the cost function through a specific choice of prior probabilities. The forced climate response, with time dependence estimated from state-of-the-art climate-model simulations, is also incorporated in the present LDM decomposition and shown to increase its optimality from a Bayesian standpoint. On top of the forced signal, the decomposition identifies five distinct LDMs of internal climate variability. The first three modes exhibit multidecadal scales, while the remaining two modes are attributable to interannual-to-decadal variability associated with El Niño-Southern oscillation; all of these modes contribute to the secular climate signal-the so-called global stadium wave-missing in the climate-model simulations. One of the multidecadal LDMs is associated with Atlantic multidecadal oscillation. The two remaining slow modes have secular time scales and patterns exhibiting regional-to-global similarities to the forced-signal pattern. These patterns have a global scale and contribute significantly to SAT variability over the Southern and Pacific Oceans. In combination with low-frequency modulation of the fast LDMs, they explain the vast majority of the variability associated with interdecadal Pacific oscillation. The global teleconnectivity of the secular climate modes and their possible crucial role in shaping the forced climate response are the two key dynamical questions brought about by the present analysis.