As detailed terrain data becomes available. GIS terrain applications target larger geographic areas at finer resolutions. Processing the massive datasets involved in such applications presents significant challenges to GIS systems and demands algorithms that are optimized for both data movement and computation. In this paper we present efficient algorithms for flow routing on massive grid terrain datasets, extending our previous work on flow accumulation. Our algorithms are developed in the framework of external memory algorithms and use I/O-techniques to achieve efficiency. We have implemented the algorithms in the Terraflow system, which is the first comprehensive terrain flow software system designed and optimized for massive data. We compare the performance of Terraflow with that of state-of-the-art commercial and open-source GIS systems. On large terrains, Terraflow outperforms existing systems by a factor of 2 to 1,000, and is capable of solving problems no system was previously able to solve.