Magnetic helicity is robustly conserved in systems with large magnetic Reynolds numbers, including most systems of astrophysical interest. This plays a major role in suppressing the kinematic large scale dynamo and driving the large scale dynamo through the magnetic helicity flux. Numerical simulations of astrophysical systems typically lack sufficient resolution to enforce global magnetic helicity over several dynamical times. Errors in the internal distribution of magnetic helicity are equally serious and possibly larger. Here we propose an algorithm for enforcing strict local conservation of magnetic helicity in the Coulomb gauge in numerical simulations.