  We present an efficient and highly scalable geometric method for
two-dimensional ideal fluid dynamics on the sphere. The starting point is
Zeitlin's finite-dimensional model of hydrodynamics. The efficiency stems from
exploiting a tridiagonal splitting of the discrete spherical Laplacian combined
with highly optimized, scalable numerical algorithms. For time-stepping, we
adopt a recently developed isospectral integrator able to preserve the
geometric structure of Euler's equations, in particular conservation of the
Casimir functions. To overcome previous computational bottlenecks, we formulate
the matrix Lie algebra basis through a sequence of tridiagonal eigenvalue
problems, efficiently solved by well-established linear algebra libraries. The
same tridiagonal splitting allows for computation of the stream matrix,
involving the inverse Laplacian, for which we design an efficient parallel
implementation on distributed memory systems. The resulting overall
computational complexity is $\mathcal{O}(N^3)$ per time-step for $N^2$ spatial
degrees of freedom. The dominating computational cost is matrix-matrix
multiplication, carried out via the parallel library ScaLAPACK. Scaling tests
show approximately linear scaling up to around $2500$ cores for the matrix size
$N=4096$ with a computational time per time-step of about $0.55$ seconds. These
results allow for long-time simulations and the gathering of statistical
quantities while simultaneously conserving the Casimir functions. We illustrate
the developed algorithm for Euler's equations at the resolution $N=2048$.
