  The demand for substantial increases in the spatial resolution of global
weather- and climate- prediction models makes it necessary to use numerically
efficient and highly scalable algorithms to solve the equations of large scale
atmospheric fluid dynamics. For stability and efficiency reasons several of the
operational forecasting centres, in particular the Met Office and the ECMWF in
the UK, use semi-implicit semi-Lagrangian time stepping in the dynamical core
of the model. The additional burden with this approach is that a three
dimensional elliptic partial differential equation (PDE) for the pressure
correction has to be solved at every model time step and this often constitutes
a significant proportion of the time spent in the dynamical core. To run within
tight operational time scales the solver has to be parallelised and there seems
to be a (perceived) misconception that elliptic solvers do not scale to large
processor counts and hence implicit time stepping can not be used in very high
resolution global models. After reviewing several methods for solving the
elliptic PDE for the pressure correction and their application in atmospheric
models we demonstrate the performance and very good scalability of Krylov
subspace solvers and multigrid algorithms for a representative model equation
with more than $10^{10}$ unknowns on 65536 cores on HECToR, the UK's national
supercomputer. For this we tested and optimised solvers from two existing
numerical libraries (DUNE and hypre) and implemented both a Conjugate Gradient
solver and a geometric multigrid algorithm based on a tensor-product approach
which exploits the strong vertical anisotropy of the discretised equation. We
study both weak and strong scalability and compare the absolute solution times
for all methods; in contrast to one-level methods the multigrid solver is
robust with respect to parameter variations.
