  A three-dimensional finite-difference solver has been developed and
implemented for Boussinesq convection in a spherical shell. The solver
transforms any complex curvilinear domain into an equivalent Cartesian domain
using Jacobi transformation and solves the governing equations in the latter.
This feature enables the solver to account for the effects of the non-spherical
shape of the convective regions of planets and stars. Apart from
parallelization using MPI, implicit treatment of the viscous terms using a
pipeline alternating direction implicit scheme and HYPRE multigrid accelerator
for pressure correction makes the solver efficient for high-fidelity direct
numerical simulations. We have performed simulations of Rayleigh-B\'enard
convection at three Rayleigh numbers $Ra=10^{5}, 10^{7}$ and $10^{8}$ while
keeping the Prandtl number fixed at unity ($Pr=1$). The average radial
temperature profile and the Nusselt number match very well, both qualitatively
and quantitatively, with the existing literature. Closure of the turbulent
kinetic energy budget, apart from the relative magnitude of the grid spacing
compared to the local Kolmogorov scales, assures sufficient spatial resolution.
