Accurate modeling of turbulence is required for investigating interactions between turbulence, blooming, and environmental and laser beam characteristics. A layered model propagates the optic field across each layer by Fresnel diffraction and a phase screen. Turbulence is modeled statistically with the Kolmogorov spectrum which goes to infinity as spatial frequency decreases. Below a cut-off frequency, the spectrum transitions to zero at a point corresponding to the largest eddies in the turbulence. Two different approaches for computing phase screens with this spectrum are considered. In the first approach, the phase screens are computed directly from the spectrum by inverse transformation. Gaussian random numbers are molded to the spectrum. A discrete 2D Fourier transform provides sample random phase screens. We show that the discrete Fourier transform cannot provide an accurate transform because of the character of the Kolmogorov spectrum. In the second approach, a covariance matrix is developed using structure functions. This avoids the need to compute a Fourier transform. The eigenvalues and eigenvectors are computed for the matrix. Eigenvalue weighted Gaussian random variables are premultiplied by the eigenvector array to generate phase screens. The advantages, disadvantages, and computational effort are discussed.