Two simple and yet sufficiently accurate approaches to predict surface temperature of a vaporizing droplet (higher order polynomial approximation and the heat balance integral method) are proposed. Being computationally inexpensive these approaches are tested as candidates for use in high-resolution LES spray modelling. Robust and efficient numerical algorithm for solving inherently stiff equations of droplet heating and evaporation has been developed. Robustness and computational efficiency of the proposed algorithm is achieved by use of unconditionally stable strongly implicit integration scheme and appropriate adaptation of the time step. The above methodology has been implemented in CFD spray model included in the in-house Fire3D code. The spray model has been applied to replicate three essentially different experimental scenarios in which turbulent sprays of water, acetone, and diesel fuel were investigated. Reasonable agreement has been demonstrated for predicted and measured droplet sizes and velocities as well as for the spray tip penetration dynamics. New numerical algorithms used to calculate surface temperatures of evaporating droplets with non-uniform internal temperature did not incur observable increase of CPU time in turbulent spray simulations.