  We present the first numerical simulations of the solar interior to exhibit a
tachocline consistent with the Gough and McIntyre (1998) model. We find
nonlinear, axisymmetric, steady-state numerical solutions in which: (1) a
large-scale primordial field is confined within the radiation zone by
downwelling meridional flows that are gyroscopically pumped in the convection
zone (2) the radiation zone is in almost-uniform rotation, with a rotation rate
consistent with observations (3) the bulk of the tachocline is magnetic free,
in thermal-wind balance and in thermal equilibrium and (4) the interaction
between the field and the flows takes place within a very thin magnetic
boundary layer, the tachopause, located at the bottom of the tachocline. We
show that the thickness of the tachocline scales with the amplitude of the
meridional flows exactly as predicted by Gough and McIntyre. We also determine
the parameter conditions under which such solutions can be obtained, and
provide a simple explanation for the failure of previous numerical attempts at
reproducing the Gough and McIntyre model. Finally, we discuss the implications
of our findings for future numerical models of the solar interior, and for
future observations of the Sun and other stars.
