Turbulent convection is a phenomenon relevant to both stellar structure and accretion disks. In the latter, a basic parameter such as the turbulent viscosity nu(t) is still treated phenomenologically; in the case of stellar structure, most of the work still relies on the mixing-length theory (MLT) which assumes homogeneity and thus lacks diffusion terms (divergence of third-order moments like w2-thetaBAR, w-theta-2BAR, q 2wBAR). To include them, one needs a new formalism. We review and discuss the Reynolds stress approach (proven successful in other fields) which provides a set of coupled differential equations that yield all the turbulent quantities of interest. Although the system can only be solved numerically, some features can be listed: 1. The convective flux F(c) = c(p)rho-w-thetaBAR is not given simply by (kappa(t) is the turbulent conductivity) F(c) = F(c)MLT almost-equal-to kappa(t)(del - del(ad)). 2. Inclusion of the diffusion terms related to w2-thetaBAR and w-theta-2BAR contributes a countergradient term GAMMA(c), which may carry heat from cold to hot regions, F(c) almost-equal-to kappa(t)(del - del(ad) + GAMMA(c)). The term GAMMA(c) was first discussed by Deardorff in the context of atmospheric turbulence. 3. Inclusion of the diffusion term related to 1/2q2wBAR (turbulent kinetic energy flux) contributes an additional term (first discussed in atmospheric turbulence by Tennekes) F(c) almost-equal-to kappa(t)(del - del(ad) + GAMMA(c)) + F(c)diff, which is responsible for overshooting. In addition to the convective flux, we also derive a model expression for nu(t) as a function of both shear and buoyency: it is needed in the numerical simulation of stellar convection and in accretion disks to replace the phenomenological expressions used thus far.