We develop a Spin-Lattice Dynamics (SLD) simulation model for ferromagnetic iron where atoms are treated as classical particles with spins. The atoms interact via many-body forces as well as via spin-orientation-dependent forces of the Heisenberg form. The coupling between the lattice and the spin degrees of freedom is described by a coordinate-dependent exchange function. An algorithm for integrating the spin-lattice dynamics equations of motion is based on the 2nd order Suzuki-Trotter decomposition for the non-commuting Liouville evolution operators for atomic coordinates and spins. The notions of the spin thermostat and the spin temperature are introduced through a combined application of the Langevin spin dynamics and the fluctuation-dissipation theorem. Several applications of the new method described in the paper illustrate the significant effect of the spin degrees of freedom on the dynamics of atomic motion in iron and iron-based alloys, and confirm that the Spin-Lattice Dynamics approach provides a viable framework for performing realistic large-scale simulations of magnetic materials.