A special place in climatology is taken by the so-called conceptual climate models. These, relatively simple, sets of differential equations can successfully describe single mechanisms of the climate. We focus on one family of such models based on the global energy balance. This gives rise to a degenerate nonlocal parabolic nonlinear partial differential equation for the zonally averaged temperature. We construct a fully discrete numerical method that has an optimal spectral accuracy in space and second order in time. Our scheme is based on Galerkin formulation of the Legendre basis expansion which is particularly convenient for this setting. By using extrapolation the numerical scheme is linear even though the original equation is strongly nonlinear. We also test our theoretical result during various numerical simulations that confirm the aforementioned accuracy of the scheme. All implementations are coded in Julia programming language with the use of parallelization (multi-threading).