<p>With the continually accumulated magnetic measurements and the gradually reliable global models of the lithospheric magnetic field by several advanced satellites (such as CHAMP, <i>Swarm</i>, CSES-1 and MSS-1), now present a requirement and also a challenge to develop realistic forward modelling methods for magnetic fields (i.e., the magnetic potential and its derivatives) that take the curvature of the Earth into account. The spatial discretization by a set of elementary tesseroids is generally utilized to approximate the complex magnetized source in spherical domain by the principle of superposition and saturate the source volume without “holes”. Since there is no analytic solution for magnetic fields of the tesseroid (except for the special points located on the polar axis), the numerical solution is the efficient way, where the Gauss–Legendre quadrature (GLQ) is usually employed. However, the required computation becomes notably time-consuming when the geometric sizes of the tesseroids are very large or the distances between the tesseroids and the observation points are very close, that is, the distance-to-size ratio (DSR) is quite small. Moreover, in an actual application, the DSRs vary with relative distances between source locations and observation points and hence are often non-uniform. Therefore, in order to reduce the computational time while maintaining a desired accuracy (i.e., relative percentage error) of each observation point, an efficient forward modelling scheme is employed. The key point of this scheme is the adoption of a new simple and efficient adaptive subdivision method. It is an equidistant subdivision method based on the longest side length, rather than recursion or stacking. By comparing the number of subdivided tesseroids, this method demonstrates its ability to avoid over-subdivision and perform more efficient calculations compared to the recursive method, because it adopts a new priori termination condition for subdivision rather than the traditional posteriori way. We obtain the required DSRs with errors of 0.1% and 0.01% through numerical simulation. At the same time, we package this scheme and release the open-source forward calculation software written by the C++ programming. Then, the analytical solution of the global homogeneous spherical shell using Runcorn’s theorem is utilized to test our newly proposed method. As a practical application, the impacts of Earth’s curvature on forward modelling of the magnetic fields are investigated.</p>