  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.
