Models of planetary, atmospheric and oceanic circulation involve eddy viscosity and eddy diffusivity, Km and K-H, that account for unresolved turbulent mixing and diffusion. The most sophisticated turbulent closure models used today for geophysical applications belong in the family of the Reynolds stress models. These models are formulated for the physical space variables; they consider a hierarchy of turbulent correlations and employ a rational way of its truncation. In the process, unknown correlations are related to the known ones via "closure assumptions" that are based upon physical plausibility, preservation of tensorial properties, and the principle of the invariant modeling according to which the constants in the closure relationships are universal. Although a great deal of progress has been achieved with Reynolds stress closure models over the years, there are still situations in which these models fail. The most difficult flows for the Reynolds stress modeling are those with anisotropy and waves because these processes are scale-dependent and cannot be included in the closure assumptions that pertain to ensemble-averaged quantities. Here, we develop an alternative approach of deriving expressions for KM and KH using the spectral space representation and employing a self-consistent, quasi-normal scale elimination QNSE) algorithm. More specifically, the QNSE procedure is based upon the quasi-Gaussian mapping of the velocity and temperature fields using the Langevin equations. Turbulence and waves are treated as one entity and the effect of the internal waves is easily identifiable. This model implies partial averaging and, thus, is scale-dependent; it allows one to easily introduce into consideration such parameters as the grid resolution, the degree of the anisotropy, and spectral characteristics, among others. Applied to turbulent flows affected by anisotropy and waves, the method traces turbulence anisotropization and shows how the dispersion relationships for linear waves are modified by turbulence. In addition, one can derive the internal wave frequency shift and the threshold criterion of internal wave generation in the presence of turbulence. The spectral method enables one to derive analytically various one-dimensional and three-dimensional spectra that reflect the effects of waves and anisotropy. When averaging is extended to all scales, the method yields a Reynolds-averaged, Navier-Stokes equations based model (RANS). This RANS model shows that there exists a range of Ri, approximately between 0.1 and 1, in which turbulence undergoes remarkable anisotropization; the vertical mixing becomes suppressed while the horizontal mixing is enhanced. Although KH decreases at large Ri and tends to its molecular value, KM remains finite and larger than its molecular value. This behavior is attributable to the effect of internal waves that mix the momentum but do not mix a scalar. In the Reynolds stress models, this feature is not replicated; instead, all Reynolds stress models predict K-M -> 0 at some value of Ri <= 1 which varies from one model to another. The presented spectral model indicates that there is no a single-valued critical Richardson number Ri at which turbulence is fully suppressed by stable stratification. This result is in agreement with large volume of atmospheric, oceanic and laboratory data. The new spectral model has been implemented in the K-epsilon format and tested in simulations of the stably stratified atmospheric boundary layers. The results of these simulations are in good greement with the data collected in BASE, SHEBA and CASES99 campaigns. Implementation of the QNSE-derived K-M and K-H in the high-resolution weather prediction system HIRLAM results in significant improvement of its predictive skills.