The computation of higher-order corrections to cross-sections relevant at LHC involves the evaluation of phase-space integrals that exhibit soft and collinear divergences. The subtraction of these divergences is a key ingredient to obtain fully-differential predictions for physical observables. We discuss a subtraction method to handle these divergences based on the construction of universal local counterterms. The integration of the counterterms is carried out analytically, giving a strong control on the numerical stability of our predictions. We implement our method in a numerical program, that we dub NNLOCAL, and validate it by computing the fully-differential NNLO cross-section for Higgs boson production in gluon-gluon fusion.