The production yield of Z bosons is measured in the electron and muon decay channels in Pb+Pb collisions at √sNN=5.02 TeV with the ATLAS detector. Data from the 2015 LHC run corresponding to an integrated luminosity of 0.49 nb^-1 are used for the analysis. The Z boson yield, normalised by the total number of minimum-bias events and the mean nuclear thickness function, is measured as a function of dilepton rapidity and event centrality. The measurements in Pb+Pb collisions are compared with similar measurements made in proton–proton collisions at the same centre-of-mass energy. The nuclear modification factor is found to be consistent with unity for all centrality intervals. The results are compared with theoretical predictions obtained at next-to-leading order using nucleon and nuclear parton distribution functions. The normalised Z boson yields in Pb+Pb collisions lie 1–3σ above the predictions. The nuclear modification factor measured as a function of rapidity agrees with unity and is consistent with a next-to-leading-order QCD calculation including the isospin effect.