The study of sound propagation in moving media is important in various fields, such as aeroacoustics and underwater acoustics. To address this problem, a three-dimensional Gaussian beam tracing model has been developed for subsonic moving media, based on the Helmholtz equation of velocity potential for high-frequency sound wave in a moving medium with arbitrary Mach numbers. By using the beam tracing method, the dynamic ray equations in the moving medium are derived, and the partial differential equation is replaced with ordinary differential equations, allowing for a more efficient and accurate calculation of the three-dimensional sound field in the moving medium. The Gaussian beam tracing method reveals that the expansion of the beam in a moving medium is more complex than in a static medium, and the energy in the ray tube is not necessarily conserved. The model has been applied to several problems, including point source sound propagation in semi-infinite homogeneous media, three-dimensional long-range sound propagation in horizontally layered atmospheres, and three-dimensional sound propagation in the Gulf Stream. The results of the point source sound propagation problem in the semi-infinite homogeneous medium verify the effectiveness and accuracy of the model. The results of the atmospheric sound propagation problem indicate that compared with the commonly used N@2D method, the three-dimensional Gaussian beam tracing in a moving medium fully considers the effect of medium motion, especially the effect of crosswind, and can calculate the sound pressure field more accurately. Although the Mach number of the ocean current is very small, its effect cannot be ignored and can quantitatively change the sound propagation mode and affect the convergence zone position. In some areas, the difference between considering and not considering the ocean current is more than 5dB. Moreover, the deviation of rays caused by lateral flow is much smaller and requires interface reflection to show more obvious deviation, and the impact of lateral flow on sound propagation is much smaller compared to the impact brought by flow velocity parallel to the propagation direction. In conclusion, the developed Gaussian beam tracing method provides an accurate and efficient approach to solve the sound propagation problem in subsonic moving media.