In overloaded Massive MIMO (mMIMO) systems, wherein the number $K$ of user equipments (UEs) exceeds the number of base station antennas $M$, it has recently been shown that non-orthogonal multiple access (NOMA) can increase the sum spectral efficiency. This paper aims at identifying cases where code-domain NOMA can improve the spectral efficiency of mMIMO in the classical regime where $K < M$. Novel spectral efficiency expressions are provided for the uplink and downlink with arbitrary spreading signatures and spatial correlation matrices. Particular attention is devoted to the planar arrays that are currently being deployed in pre-5G and 5G networks (in sub$-6$ GHz bands), which are characterized by limited spatial resolution. Numerical results show that mMIMO with such planar arrays can benefit from NOMA in scenarios where the UEs are spatially close to each other. A two-step UE grouping scheme is proposed for NOMA-aided mMIMO systems that is applicable to the spatial correlation matrices of the UEs that are currently active in each cell. Numerical results are used to investigate the performance of the algorithm under different operating conditions and types of spreading signatures (orthogonal, sparse and random sets). The analysis reveals that orthogonal signatures provide the highest average spectral efficiency.