Acoustical shock waves can be generated by numerous atmospheric sources, either natural – like thunder and volcanoes – or anthropic – like explosions, sonic boom or buzz saw noise. The prediction of their long-range propagation remains a numerical challenge at 3D because of the large propagation distance to wavelength ratio, and of the high frequency / small wavelength content associated to shocks. In this paper, an original numerical method for propagating acoustical shock waves in three-dimensional heterogeneous media is proposed. Heterogeneities can result from temperature or density gradients and also from atmospheric shear and turbulent flows. The method called FLHOWARD (for FLow and Heterogeneities in a One-Way Approximation of the nonlineaR wave equation in 3D) is based on a one-way solution of a generalized nonlinear wave equation. Even though backscattering is neglected, it does not suffer from the limitations of classical ray theory nor from the angular limitations of the popular parabolic methods. The numerical approach is based on a split-step method, which has the advantage of splitting the original equation into simpler ones associated with specific physical mechanisms: diffraction, flows, heterogeneities, nonlinearities, absorption and relaxation. The method has been developed on parallel architecture for very high demanding 3D configurations using the Single Method Multiple Data paradigm. The method is validated through several test cases. A study of the lateral cut-off of the sonic boom finally illustrates the potentialities of the method for realistic cases.