  In this paper we present a next-to-next-to-leading order (NNLO) calculation
of the process $pp\rightarrow \gamma\gamma$ that we have implemented into the
parton level Monte Carlo code MCFM. We do not find agreement with the previous
calculation of this process in the literature. In addition to the
$\mathcal{O}(\alpha_s^2)$ corrections present at NNLO, we include some effects
arising at $\mathcal{O}(\alpha_s^3)$, namely those associated with
gluon-initiated closed fermion loops. We investigate the role of this process
in the context of studies of QCD at colliders and as a background for searches
for new physics, paying particular attention to the diphoton invariant mass
spectrum. We demonstrate that the NNLO QCD prediction for the shape of this
spectrum agrees well with functional forms used in recent data-driven fits.
