  We study the thermal evolution of axisymmetric rotating neutron stars in full
general relativity. To this aim we develop "NSCool 2D Rot", a major upgrade of
the 1D neutron stars thermal evolution code "NSCool" by D. Page. As a first
application of our new code we address the standard cooling of isolated neutron
stars with rotation frequencies up to the mass shedding limit. We investigate
the effects of the equation of state (EOS) by considering different
combinations of core and crust EOSs. The results indicate complex
time-dependent evolution of temperature distribution throughout the whole
volume of the star and, in particular, in the crust. We show that most of that
complexity can be attributed to the formation of a "heat blob" in the crust and
to the latitude dependence of the heat diffusion timescale through the crust.
