Subcooled boiling flow is widely used for heat and mass transfer in industries such as nuclear reactor, electronic equipment, and spacecraft, where the bubble dynamics is essential for understanding the gas-liquid flow behavior. Current studies focused on the gas void fraction or bubble size distribution at a global scale, but the dynamics of single bubbles are still obscure, which inhibits the fundamental understanding of gas-liquid flow and the precise prediction of heat transfer. Here, we develop a numerical method that captures the motion and condensation of bubbles in subcooled boiling flow using Euler-Lagrange method, with a proposed model for spatiotemporal lift-off bubble distribution on the heated wall as boundary condition. Our numerical method is validated by the agreement between the calculated gas void fraction distributions with the existing experimental data in vertical subcooled boiling flow. We further investigate the distribution of bubble position and size in the tube and find that the bubble size distribution shows different modes at different radial positions. Analysis of the single bubble dynamics shows that small and large bubbles suffer from centrifugal and centripetal lift forces, and thus tend to accumulate at a radial position near the wall and near the tube center, respectively, which reveals the mechanism of bubble size distribution. Besides, we find that the centripetal wall lubrication force is essential for the path instability of small bubbles near the wall. These findings promote the fundamental understanding of the gas-liquid phase distribution in the subcooled boiling flow and provide a numerical strategy for the precise prediction of mass and heat transfer in subcooled boiling flow.