In this study the authors apply a recently developed clustering method for the systematic identification of metastable atmospheric regimes in high-dimensional datasets generated by atmospheric models. The novelty of this approach is that it decomposes the phase space in, possibly, overlapping clusters and simultaneously estimates the most likely switching sequence among the clusters. The parameters of the clustering and switching are estimated by a finite element approach. The switching among the clusters can be described by a Markov transition matrix. Possible metastable regime behavior is assessed by inspecting the eigenspectrum of the associated transition probability matrix. The recently introduced metastable data-analysis method is applied to high-dimensional datasets produced by a barotropic model and a comprehensive atmospheric general circulation model (GCM). Significant and dynamically relevant metastable regimes are successfully identified in both models. The metastable regimes in the barotropic model correspond to blocked and zonal states. Similar regime states were already previously identified in highly reduced phase spaces of just one and two dimensions in the same model. Next, the clustering method is applied to a comprehensive atmospheric GCM in which seven significant flow regimes are identified. The spatial structures of the regimes correspond to, among others, both phases of the Northern Annular Mode and Pacific blocking. The regimes are maintained predominantly by transient eddy fluxes of low-pass-filtered anomalies. It is demonstrated how the dynamical description of the slow process switching between the regimes can be acquired from the analysis results, and an investigation of the resulting simplified dynamical model with respect to predictability is performed. A predictability study shows that a simple Markov model is able to predict the regimes up to six days ahead, comparable to the ability of high-resolution state-of-the-art numerical weather prediction models to accurately predict the onset and decay of blockings. The implications of the results for derivation of reduced models for extended-range predictability are discussed.