In this paper, we propose a model reduction method for linear systems evolving on large-scale complex networks, called dynamical networks. In this method, we construct a set of clusters (i.e., disjoint subsets of state variables) based on a notion of cluster controllability that characterizes local controllability of the state-space of the dynamical networks. We aggregate the constructed clusters to obtain a reduced model that preserves connection topology of the original system as well as the stability and some particular properties, such as steady-state characteristic and system positivity. In addition, we derive an H_∞-error bound of the state discrepancy caused by the aggregation. The efficiency of the proposed method is shown through a numerical example including a large-scale complex network.