Numerical simulations of isotropic, homogeneous, forced and dissipative two-dimensional (2D) turbulence in the energy transfer subrange are complicated by the inverse cascade that continuously propagates energy to the large scale modes. To avoid energy condensation in the lowest modes, an energy sink, or a large scale drag is usually introduced. With a few exceptions, simulations with different formulations of the large scale drag reveal the development of strong coherent vortices and steepening of energy and enstrophy spectra that lead to erosion and eventual destruction of Kolmogorov–Batchelor–Kraichnan (KBK) statistical laws. Being attributed to the intrinsic anomalous fluctuations independent of the large scale drag formulation, these coherent vortices have prompted conjectures that KBK 2D turbulence in the energy subrange is irreproducible in long term simulations. Here, we advance a different point of view, according to which the emergence of coherent vortices is triggered by the inverse energy cascade distortion directly attributable to the choice of a large scale drag formulation. We subdivide the computational modes into explicit and implicit, or supergrid scale (SPGS), which are the few lowest wave numbers modes that adhere to KBK statistics. Then, we introduce a new concept of the large scale drag—rather than being an energy sink, it accounts for the energy and enstrophy exchange between the explicit and SPGS modes. The new SPGS parameterization was used in both direct numerical simulations (DNS) and large eddy simulations (LES) in a doubly periodic box setting. It was found that the new technique enables both DNS and LES to reach a steady state preserved for many large scale eddy turnover times. For the entire time of integration, the flow field remained structureless and in good agreement with the KBK statistical laws. We conclude that homogeneous, isotropic, forced, dissipative 2D turbulence in the energy subrange is statistically stable, does not produce coherent structures, and obeys the KBK statistical laws for as long as its inverse energy cascade remains undisturbed. The proposed new technique of computing the intermediate modes while the statistics of the largest scales is known may find a wide range of applications.