Multireservoir models are widely used for modeling and interpreting ultrafast magnetization dynamics. Here we introduce an alternative formulation to existing three-temperature models for the treatment of spin, electron, and lattice temperatures in magnetization dynamics simulations. In contrast to most existing models of calculations of energy transfer between reservoirs in these types of simulations, the heat distribution of the spin and lattice subsystems is evaluated during the simulation instead of being calculated a priori. The model is applied to investigate the demagnetization and remagnetization of fcc Ni, when subjected to a strong laser pulse. In particular, our model results in a fast interplay between the electron and spin subsystems which reproduces the main features of experimental data for fcc Ni significantly better than most reported three-temperature models. We also show that the way in which the electron, spin, and lattice heat capacities are described can have a significant impact on the simulated ultrafast dynamics. By introducing spin-lattice couplings in the simulation, it is shown that these explicit interactions only have a marginal impact on the magnetization dynamics of fcc Ni, albeit it is more pronounced for higher laser pulse powers.