  We present a detailed numerical study of the Gough & McIntyre model for the
solar tachocline. This model explains the uniformity of the rotation profile
observed in the bulk of the radiative zone by the presence of a large-scale
primordial magnetic field confined below the tachocline by flows originating
from within the convection zone. We attribute the failure of previous numerical
attempts at reproducing even qualitatively Gough & McIntyre's idea to the use
of boundary conditions which inappropriately model the radiative--convective
interface. We emphasize the key role of flows downwelling from the convection
zone in confining the assumed internal field. We carefully select the range of
parameters used in the simulations to guarantee a faithful representation of
the hierarchy of expected lengthscales. We then present, for the first time, a
fully nonlinear and self-consistent numerical solution of the Gough & McIntyre
model which qualitatively satisfies the following set of observational
constraints: (i) the quenching of the large-scale differential rotation below
the tachocline - including in the polar regions - as seen by helioseismology
(ii) the confinement of the large-scale meridional flows to the uppermost
layers of the radiative zone as required by observed light element abundances
and suggested by helioseismic sound-speed data.
