We study the statistical properties of spherical harmonic modes of temperature maps of the cosmic microwave background. Unlike other studies, which focus mainly on properties of the amplitudes of these modes, we look instead at their phases. In particular, we present a simple measure of phase correlation that can be diagnostic of departures from the standard assumption that primordial density fluctuations constitute a statistically homogeneous and isotropic Gaussian random field, which should possess phases that are uniformly random on the unit circle. The method we discuss checks for the uniformity of the distribution of phase angles using a non-parametric descriptor based on the use order statistics, which is known as Kuiper's statistic. The particular advantage of the method we present is that, when coupled to the judicious use of Monte Carlo simulations, it can deliver very interesting results from small data samples. In particular, it is useful for studying the properties of spherical harmonics at low l for which there are only small number of independent values of m and which therefore furnish only a small number of phases for analysis. We apply the method to the COBE-DMR and WMAP sky maps, and find departures from uniformity in both. In the case of WMAP, our results probably reflect Galactic contamination or the known variation of signal-to-noise across the sky rather than primordial non-Gaussianity.