We provide a method for fast and exact simulation of Gaussian random fields on the sphere having isotropic covariance functions. The method proposed is then extended to Gaussian random fields defined over the sphere cross time and having covariance functions that depend on geodesic distance in space and on temporal separation. The crux of the method is in the use of block circulant matrices obtained working on regular grids defined over longitude and latitude.