  Solving inverse problems and achieving statistical rigour in landscape
evolution models requires running many model realizations. Parallel computation
is necessary to achieve this in a reasonable time. However, no previous
algorithm is well-suited to leveraging modern parallelism. Here, I describe an
algorithm that can utilize the parallel potential of GPUs, many-core
processors, and SIMD instructions, in addition to working well in serial. The
new algorithm runs 43x faster (70s vs. 3,000s on a 10,000x10,000 input) than
the previous state of the art and exhibits sublinear scaling with input size. I
also identify methods for using multidirectional flow routing and quickly
eliminating landscape depressions and local minima. Tips for parallelization
and a step-by-step guide to achieving it are given to help others achieve good
performance with their own code. Complete, well-commented, easily adaptable
source code for all versions of the algorithm is available as a supplement and
on Github.
