A time splitting numerical model is developed in this paper, which combined the explicit algorithm in the horizontal direction and the implicit algorithm in the vertical direction, and used to simulate the nonlinear propagation of the gravity wave (GW) in a 2-dimensional compressible atmo- sphere. Numerical simulation of the nonlinear propagation of the small amplitude GW is presented by using this model. The numerical results coincide well with linear theory and indicate that the numerical model is correct. The saturation and breakdown of the nonlinear propagation of the finite amplitude GW is simulated by using the model. The results show that: (1) The overturn appears before the saturation of GW, however, it requires a long time (nearly 1.5 periods of the GW) to induce the GW breakdown. The height (time) of instability predicted by linear saturation theory is higher (latter) than the corresponding results obtained from nonlinear numerical simulation because of the nonlinear wave-wave and wave-flow interaction. (2) The nonlinear wave-flow interaction produces energy transforring from the GW to mean flow prior to instability. The nonlinear wave-wave interaction is induced by the breakdown of GW directly. (3) The direction of the horizontal mean wind acceleration, the jet and the horizontal propagation of GW are consistent and indicate that the nonlinear wave-flow interaction accelerate the formation of horizontal mean wind shear and the development of the instability.