[eng] Gravitational-waveobservations of binary black holes currently rely on theoretical modelsthat predict the dominantmultipoles (l = 2,|m| = 2) oftheradiationduringinspiral,merger,andringdown.Weintroducea simple method to include the subdominant multipoles to binary black hole gravitational waveforms, given a frequency-domain model for the dominant multipoles. The amplitude and phase of the original model are appropriately stretched and rescaled using post-Newtonian results (for the inspiral), perturbation theory (for the ringdown), and a smooth transition between the two. No additional tuning to numerical-relativity simulations is required. We apply avariant of this method to the nonprecessing PhenomD model. The result, PhenomHM ,constitutes the first higher-multipole model ofspinning andcoalescing black-holebinaries,and currently includes the (l.|m|)= 2 (2,2) (3,3) (4,4) (2,1) (4,2)(4,3)radiative moments. Comparisons with numerical-relativity waveforms demonstrate that PhenomHM is more accurate than dominant-multipole- only models for all binary configurations, and typically improves the measurement of binary properties