In this work we present an extension of the time domain phenomenological model IMRPhenomT for gravitational wave signals from binary black hole coalescences to include subdominant harmonics, specifically the $(l=2, m=pm 1)$, $(l=3, m=pm 3)$, $(l=4, m=pm 4)$ and $(l=5, m=pm 5)$ spherical harmonics. We also improve our model for the dominant $(l=2, m=pm 2)$ mode and discuss mode mixing for the $(l=3, m=pm 2)$ mode. The model is calibrated to numerical relativity solutions of the full Einstein equations up to mass ratio 18, and to numerical solutions of the Teukolsky equations for higher mass ratios. This work complements the latest generation of traditional frequency domain phenomenological models (IMRPhenomX), and provides new avenues to develop computationally efficient models for gravitational wave signals from generic compact binaries.