  We discuss recent advances made in modelling the complex magnetohydrodynamics
of the Sun using our anelastic spherical harmonics (ASH) code. We have
conducted extensive 3--D simulations of compressible convection in rotating
spherical shells with and without magnetic fields, to study the coupling
between global-scale convection and rotation in seeking to understand how the
solar differential rotation is established and maintained. Such simulations
capable of studying fairly turbulent convection have been enabled by massively
parallel supercomputers. The resulting convection within domains that capture a
good fraction of the bulk of the convection zone is highly time dependent and
intricate, and is dominated by intermittent upflows and networks of strong
downflows. A high degree of coherent structures involving downflowing plumes
can be embedded in otherwise chaotic flow fields. These vortical structures
play a significant role in yielding Reynolds stresses that serve to
redistribute angular momentum, leading to differential rotation profiles with
pole to equator contrasts of about 30% in angular velocity $\Omega$ and some
constancy along radial line at mid latitudes, thereby making good contact with
deductions from helioseismology. When a magnetic field is introduced, a dynamo
regime can be found that does not destroy the strong differential rotation
achieved in pure hydrodynamics cases. The magnetic fields are found to
concentrate around the downflowing networks and to have significant north-south
asymmetry and helicity.
