This study proposes an algorithm for modeling compressible flows in spherical shells in nearly incompressible and weakly compressible regimes based on an implicit direction splitting approach. The method retains theoretically expected convergence rates and remains stable for extremely small values of the characteristic Mach number. The staggered spatial discretization on the MAC stencil, commonly used in numerical methods for incompressible Navier-Stokes equations, was found to be convenient for the discretization of the compressible Navier-Stokes equations written in the non-conservative form in terms of the primitive variables. This approach helped to avoid the high-frequency oscillations without any artificial stabilization terms. Nonlinear Picard iterations with the splitting error reduction were also implemented to allow one to obtain a solution of the fully nonlinear system of equations. These results, alongside excellent parallel performance, prove the viability of the direction splitting approach in large-scale high-resolution high-performance simulations of atmospheric and oceanic flows.