Many geophysical and astrophysical phenomena such as magnetic fields generation, or the differential rotation observed in the atmospheres of the major planets are studied by means of numerical simulations of the Navier-Stokes equations in rotating spherical shells. Two different computational codes, spatially discretized using spherical harmonics in the angular variables, are presented. The first code, PARODY, solves the magneto-hydrodynamic anelastic convective equations with finite a difference discretization in the radial direction. This allows the parallelization on distributed memory computers to run massive numerical simulations of second order in time. It is mainly designed to perform direct numerical simulations. The second code, SPHO, solves the fully spectral Boussinesq convective equations, and its variationals, parallelized on shared memory architectures and it uses optimized linear algebra libraries. High-order time integration methods are implemented to allow the use of dynamical systems tools for the study of complex dynamics.