Thurbulent thermal convection between two plates, bottom plate is at higher temperature $T_{h}$ and the upper plate is at lower temperature $T_{i}$ is numerically investigated. Model equations are abridged Reynolds stress equations; full Reynolds stress equations are simplified to yield algebraic relations in case of mean square velocity fluctuations in vertical and horizontal directions. Boundary conditions for turbulent kinetic energy k and mean square temperature variance .thera.$^{2}$oner bar at the plate surfaces are set to be zero and those of dissipation rate of turbulent kinetic energy .epsilon. and dissipation rate of mean square temperature variance .epsilon.$_{ heta}$ are assumed at first grid point nearest to the boundary surfaces, whose values are approximated by inviscid estimates. Results show that temperature profiles are in good agreement with experimental data except transition region, in which temperature is over-predicted. Such discrepancy becomes larger as the Rayleigh number becomes smaller. Nusselt numbers, which are calculated from the temperature gradients at the boundary surfaces, are also in good agreement with experimental data.a.a.