Theories of gravity other than general relativity (GR) can explain the observed cosmic acceleration without a cosmological constant. One such class of theories of gravity is f(R). Metric f(R) theories have been proven to be equivalent to Brans-Dicke (BD) scalar-tensor gravity without a kinetic term. Using this equivalence and a 3+1 decomposition of the theory it has been shown that metric f(R) gravity admits a well-posed initial value problem. However, it has not been proven that the 3+1 evolution equations of metric f(R) gravity preserve the (hamiltonian and momentum) constraints. In this paper we show that this is indeed the case. In addition, we show that the mathematical form of the constraint propagation equations in BD-equilavent f(R) gravity and in f(R) gravity in both the Jordan and Einstein frames, is exactly the same as in the standard ADM 3+1 decomposition of GR. Finally, we point out that current numerical relativity codes can incorporate the 3+1 evolution equations of metric f(R) gravity by modifying the stress-energy tensor and adding an additional scalar field evolution equation. We hope that this work will serve as a starting point for relativists to develop fully dynamical codes for valid f(R) models.