Abstract. The computational cost of a spectral model using spherical harmonics (SH) increases significantly at high resolution because the transform method with SH requires O(N3) operations, where N is the truncation wavenumber. One way to solve this problem is to use double Fourier series (DFS) instead of SH, which requires O(N2 log N) operations. This paper proposes a new DFS method that improves the numerical stability of the model compared with the conventional DFS methods by adopting the following two improvements: a new expansion method that employs the least-squares method (or the Galerkin method) to calculate the expansion coefficients in order to minimize the error caused by wavenumber truncation, and new basis functions that satisfy the continuity of both scalar and vector variables at the poles. In the semi-implicit semi-Lagrangian shallow water model using the new DFS method, the Williamson test cases 2 and 5 and the Galewsky test case give stable results without the appearance of high-wavenumber noise near the poles, even without using horizontal diffusion. The new DFS model is faster than the SH model, especially at high resolutions, and gives almost the same results.