The accuracy of advection has a great influence on the visual effect of fluid simulation. Constrained interpolation profile (CIP) method has been an important advection scheme because of its third-order accuracy and the fact that it only needs to be performed over a compact stencil, but extending it to high-dimensional advection equations is not easy, because it involves complex calculations and large memory overheads, and is usually unstable. In this article, we propose a stable and efficient three-dimensional (3D) CIP scheme which can maintain high accuracy but requires low computation and memory cost. We first construct an efficient two-dimensional (2D) CIP scheme based on dimensional splitting and local Taylor expansions, and then propose an effective way to extend it for 3D applications without decreasing the computational accuracy or affecting the stability. The experimental results show the advantages of our method over the state-of-the-art advection schemes. K E Y W O R D S constrained interpolation profile, fluid simulation, high-order advection, Taylor expansion 1 INTRODUCTION Various kinds of methods 1 have been proposed to solve the advection step in Eulerian fluid simulation, but most of them can only reach the second-order accuracy, such as back and forth error compensation and correction (BFECC) 2,3 and MacCormack, 4 and so forth. The advection methods with higher order accuracy usually need to be performed over a wide stencil, such as the essentially nonoscillatory (ENO) 5 and weighted ENO (WENO) 5 methods, which results in a large amount of calculation and cannot be extended to nonuniform grids.