We introduce \TEOBiResumSM{}, an improved version of the effective-one-body (EOB) waveform model \TEOBResumS{} for spin-aligned, coalescing black hole binaries, that includes subdominant gravitational waveform modes completed through merger and ringdown. Beyond the dominant $(\ell,|m|)=(2,2)$ one, the more robust multipoles all over the parameter space are: $(2,1)$, $(3,3)$, $(3,2)$, $(4,4)$ and $(5,5)$. The multipolar ringdown EOB waveform stems from suitably fitting many numerical relativity (NR) waveform data from the Simulating eXtreme Spacetimes (SXS) collaboration together with test-mass waveform data. Mode-mixing effects are not incorporated. The orbital (nonspinning) part of the multipolar waveform amplitudes includes test-mass results up to (relative) 6PN order and, for most modes,is Pad\'e resummed. The $m$=odd waveform multipoles (up to $\ell=5$) incorporate most of the currently available spin-dependent analytical information. Improving on previous work, we confirm that certain $m=\text{odd}$ modes, e.g. the $(2,1)$, and even the $(3,1)$, may develop a zero (or a minimum) in the amplitude for nearly equal-mass binaries and for several combinations of the individual spins. A remarkable EOB/NR agreement around such zero is found for these modes. The new waveform, and radiation reaction, prompts a new NR-calibration of the spinning sector of the model, done with only $32$ datasets.The maximum $(2,2)$ EOB/NR unfaithfulness $\bar{F}$ with Advanced LIGO noise against the SXS catalog ($\sim 595$ datasets) is always below $0.5\%$ for binaries with total mass $M$ as $10M_\odot\leq M \leq 200M_\odot$, except for a single outlier with $\max{(\bar{F})}\sim 0.85\%$. When $(2,1)$, $(3,3)$ and $(4,4)$ modes are included, one finds an excellent EOB/NR agreement up to $M\sim 120M_\odot$, above which the performance degrades slightly and moves above $3\%$