Spherical harmonics power spectrum of the geopotential field of Gaussian-bell type on the sphere was investigated using integral formula that is associated with Legendre polynomials. The geopotential field of Gaussian-bell type is defined as a function of sine of angular distance from the bell's center in order to guarantee the continuity on the global domain. Since the integral-formula associated with the Legendre polynomials was represented with infinite series of polynomial, an estimation method was developed to make the procedure computationally efficient while preserving the accuracy. The spherical harmonics power spectrum was shown to vary significantly depending on the scale parameter of the Gaussian bell. Due to the accurate procedure of the new method, the power (degree variance) spanning over orders that were far higher than machine roundoff was well explored. When the scale parameter (or width) of the Gaussian bell is large, the spectrum drops sharply with the total wavenumber. On the other hand, in case of small scale parameter the spectrum tends to be flat, showing very slow decaying with the total wavenumber. The accuracy of the new method was compared with theoretical values for various scale parameters. The new method was found advantageous over discrete numerical methods, such as Gaussian quadrature and Fourier method, in that it can produce the power spectrum with accuracy and computational efficiency for all range of total wavenumber. The results of present study help to determine the allowable maximum scale parameter of the geopotential field when a Gaussian-bell type is adopted as a localized function.