We present a local but fully nonlinear model of the solar tachocline, using three-dimensional direct numerical simulations. The tachocline forms naturally as a statistically steady balance between Coriolis, pressure, buoyancy and Lorentz forces beneath a turbulent convection zone. Uniform rotation is maintained in the radiation zone by a primordial magnetic field, which is confined by meridional flows in the tachocline and convection zone. Such balanced dynamics has previously been found in idealised laminar models, but never in fully self-consistent numerical simulations.