The production yield of $Z$ bosons is measured in the electron and muon decay channels in Pb+Pb collisions at $\sqrt{s_{\textrm{NN}}}$ = 5.02 TeV with the ATLAS detector. Data from the 2015 LHC run corresponding to an integrated luminosity of 0.49 $\textrm{nb}^{-1}$ are used for the analysis. The $Z$ boson yield, normalised by the total 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$$\sigma$ 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.