  In this study we have developed a flexible and efficient numerical scheme for
the simulation of three-dimensional incompressible flows in spherical
coordinates. The main idea, inspired by a similar strategy as (Verzicco, R.,
Orlandi, P., 1996, A Finite-Difference Scheme for Three-Dimensional
Incompressible Flows in Cylindrical Coordinates) for cylindrical coordinates,
consists of a change of variables combined with a discretization on a staggered
mesh and the special treatment of few discrete terms that remove the
singularities of the Navier-Stokes equations at the sphere centre and along the
polar axis. This new method alleviates also the time step restrictions
introduced by the discretization around the polar axis while the sphere centre
still yields strong limitations, although only in very unfavourable flow
configurations.
  The scheme is second-order accurate in space and is verified and validated by
computing numerical examples that are compared with similar results produced by
other codes or available from the literature.
  The method can cope with flows evolving in the whole sphere, in a spherical
shell and in a sector without any change and, thanks to the flexibility of
finite-differences, it can employ generic mesh stretching (in two of the three
directions) and complex boundary conditions.
