Random fields on the sphere play a fundamental role in the natural sciences. This paper presents a simulation algorithm parenthetical to the spectral turning bands method used in Euclidean spaces, for simulating scalar- or vector-valued Gaussian random fields on the $d$-dimensional unit sphere. The simulated random field is obtained by a sum of Gegenbauer waves, each of which is variable along a randomly oriented arc and constant along the parallels orthogonal to the arc. Convergence criteria based on the Berry-Esséen inequality are proposed to choose suitable parameters for the implementation of the algorithm, which is illustrated through numerical experiments. A by-product of this work is a closed-form expression of the Schoenberg coefficients associated with the Chentsov and exponential covariance models on spheres of dimensions greater than or equal to 2.