Estimates of the source parameters of gravitational-wave (GW) events produced by compact binary mergers rely on theoretical models for the GW signal. We present the first frequency-domain model for inspiral, merger and ringdown of the GW signal from precessing binary-black-hole systems that also includes multipoles beyond the leading-order quadrupole. Our model, {tt PhenomPv3HM}, is a combination of the higher-multipole non-precessing model {tt PhenomHM} and the spin-precessing model {tt PhenomPv3} that includes two-spin precession via a dynamical rotation of the GW multipoles. We validate the new model by comparing to a large set of precessing numerical-relativity simulations and find excellent agreement across the majority of the parameter space they cover. For mass ratios $<5$ the mismatch improves, on average, from $sim6%$ to $sim 2%$ compared to {tt PhenomPv3} when we include higher multipoles in the model. However, we find mismatches $sim8%$ for the mass-ratio $6$ and highly spinning simulation. As a first application of the new model we have analysed the binary black hole event GW170729. We find larger values for the primary black hole mass of $58.25^{+11.73}_{-12.53} , M_odot$ (90% credible interval). The lower limit ($sim 46 , M_odot$) is comparable to the proposed maximum black hole mass predicted by different stellar evolution models due to the pulsation pair-instability supernova (PPISN) mechanism. If we assume that the primary ac{BH} in GW170729 formed through a PPISN then out of the four PPISN models we considered only the model of Woosley (2017) is consistent with our mass measurements at the 90% level.