  In this paper we present a next-to-next-to-leading order (NNLO) QCD
calculation of the processes $pp\rightarrow l^+l^-\gamma$ and $pp\rightarrow
\nu\bar\nu\gamma$ that we have implemented in MCFM. Our calculation includes
QCD corrections at NNLO both for the Standard Model (SM) and additionally in
the presence of $Z\gamma\gamma$ and $ZZ\gamma$ anomalous couplings. We compare
our implementation, obtained using the jettiness slicing approach, with a
previous SM calculation and find broad agreement. Focusing on the sensitivity
of our results to the slicing parameter, we show that using our setup we are
able to compute NNLO cross sections with numerical uncertainties of about
$0.1\%$, which is small compared to residual scale uncertainties of a few
percent. We study potential improvements using two different jettiness
definitions and the inclusion of power corrections. At $\sqrt{s}=13$ TeV we
present phenomenological results and consider $Z\gamma$ as a background to
$H\to Z\gamma$ production. We find that, with typical cuts, the inclusion of
NNLO corrections represents a small effect and loosens the extraction of limits
on anomalous couplings by about $10\%$.
