High-order numerical methods are attracting growing attention for numerical weather prediction, yet several challenges must be addressed before they can be considered viable alternatives. Atmospheric models often deal with slow-moving weather systems that involve very low wind speeds. In this regime, numerical schemes for compressible flows commonly encounter challenges related to convergence and accuracy. This issue mainly arises because convective propagation speeds are much slower than acoustic wave speeds, and because the dissipation introduced by upwind fluxes prevents the discrete compressible equations from properly approaching the incompressible limit. Another challenge commonly faced during the study of many atmospheric phenomena is maintaining a non-trivial stationary solution, known as the hydrostatic steady state. Classical numerical schemes often fail to preserve the steady state over long periods and also lead to erroneous solutions when trying to compute small perturbations around the stationary solution. This is particularly important because most atmospheric flow problems of practical interest are characterized by small perturbations from the hydro-static steady state. Both of these challenges typically require very fine meshes to achieve satisfactory accuracy. However, using such high-resolution grids is often impractical for large-scale three-dimensional simulations. In this work, we incorporate a low-Mach number flux function and a deviation-based well-balanced technique into the global atmospheric dynamical core WxFactory, which utilizes direct flux reconstruction and exponential time integration on a rotated cubed-sphere grid. We demonstrate the impact of this choice over a series of numerical examples, including convergence and validation studies.