A subharmonic method was proposed based on the multi-order frequency grids dividing for the generation of the turbulent phase screen to compensate the low frequency components. Frequency grids were divided by multiple orders in the frequency domain of phase screen, and the frequency grids of each order were further finely subdivided, and then the frequency components in low frequency range which contribute the most to the turbulent spectrum could be compensated sufficiently. The phase structure function and the power spectrum of atmospheric turbulence were used to analyze the accuracy of the generated phase screen. Simulation results show that the discrepancy between the phase structure function derived from the generated phase screen and the theoretical value directly calculated from the turbulent spectrum is decreased by the increasing of the total number of frequency grid orders M, and frequency grid point number of each order N′xand N′y. For every group of M, N′xand N′y, there is an optimum value for frequency space before subdividing Δf1. By selecting M = 3, N′x=N′y=3 and Δf1=2 m-1, the phase structure function and the power spectrum of the screen can be quite closely to theory. The overall slope in the generated random phase screen reflects an accurate sampling of the low frequency components.