In this paper we employ renormalized viscosity and thermal diffusivity to construct a subgrid-scale model for large eddy simulation (LES) of turbulent thermal convection. For LES, we add $\nu_\mathrm{ren} \propto \Pi_u^{1/3} (\pi/\Delta)^{-1/3}$ to the kinematic viscosity; here $\Pi_u$ is the turbulent kinetic energy flux, and $\Delta$ is the grid spacing. In our model, the turbulent Prandtl number is unity. We performed LES of turbulent thermal convection on a $128^3$ grid and compare the results with direct numerical simulation (DNS) on a $512^3$ grid. There is a good agreement between the LES and DNS results on the evolution of kinetic energy and entropy, spectra and fluxes of velocity and temperature fields, and the isosurfaces of temperature. We also show the capability of our LES to simulate thermal convection at very high Rayleigh numbers and exhibit some results for $\mathrm{Ra=10^{18}}$.