A high-order numerical method for three-dimensional hydrodynamics is presented. The present method applies high-order compact schemes in space and a Runge-Kutta scheme in time to solve the Reynolds-averaged Navier-Stokes equations with the $k-ϵ$ turbulence model in an orthogonal curvilinear coordinate system. In addition, a two-dimensional equation is derived from the depth-averaged momentum equations to predict the water level. The proposed method is first validated by its application to simulate flow in a $180^◦$ curved laboratory flume. It is found that the simulated results agree with measurements and are better than those from SIMPLEC algorithm. Then the method is applied to study three-dimensional hydrodynamics in a natural river, and the simulated results are in accordance with measurements.