Spin dynamics (SD) methods have been developed to compute NMR paramagnetic relaxation enhancements (NMR-PRE) produced by solutes with electron spin S greater than or equal to 1 in solution. The spin dynamics algorithms, which are implemented as the computer program SpinDyn.f, are similar in spirit to molecular dynamics calculations in statistical mechanics, except that the spin motion is propagated numerically in time using quantum mechanical equations of motion of the spin operators, rather than Newtonian equations of motion of the molecular degrees of freedom as in MD simulations. SD simulations as implemented in SpinDyn.f provide accurate, flexible, and rapid calculations of NMR-PRE phenomena with few of the assumptions or limitations of previous analytical theories. The program calculates inter- and intramolecular NMR-PRE phenomena for both integer and half-integer spin systems processing under arbitrary Zeeman and zfs Hamiltonians in the presence of Brownian reorientation. Isotropic Brownian reorientation is simulated by means of a finite-step algorithm with adjustable step size. Simulations computed by SpinDyn.f have been used in a systematic study aimed at better understanding the influence of Brownian reorientation on the NMR-PRE of an S=1 ion in a non-Zeeman-limit physical situation. Conditions required for the validity of zfs-limit analytical theory an given. SpinDyn.f has also been used to assess quantitatively the effects of molecular reorientation on a prior analysis of NMR-PRE data for the model S=2 complex ion [tris(acetylacetonato)manganese(III)] in acetone solution; this system was found to be well described by zfs-limit analytical theory.