The study of sound propagation in moving medium is important in various fields, such as atmospheric sound and underwater acoustics. To address this problem, a three-dimensional Gaussian beam tracing model is developed for subsonic moving medium, based on the Helmholtz equation of velocity potential for high-frequency sound wave in a moving medium with arbitrary Mach numbers. The dynamic ray equations in the moving medium are derived by using the beam tracing method, and further the partial differential equation is transformed into ordinary differential equations, so as to be able to more efficiently and accurately calculate 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 is applied to several problems, including point source sound propagation in a semi-infinite homogeneous medium, 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. The ocean current can quantitatively change the sound propagation mode and affect the convergence zone position. In some areas, the difference between calculation results with and without considering the ocean current is more than 5 dB. Moreover, the deviation of rays caused by lateral flow is much smaller, and even in the areas with complex terrain, the deviation becomes more obvious only after being reflected by the interface. Moreover, the influence of lateral flow on sound propagation is much smaller than that caused by flow velocity parallel to the propagation direction. In conclusion, the developed Gaussian beam tracing method provides an accurate and efficient approach for solving the sound propagation problem in subsonic moving media.