The numerical implementation of the Green’s function parabolic equation (GFPE) method for atmospheric sound propagation is discussed. Four types of numerical errors are distinguished: (i) errors in the forward Fourier transform; (ii) errors in the inverse Fourier transform; (iii) errors in the refraction factor; and (iv) errors caused by the split-step approximation. The sizes of the errors depend on the choice of the numerical parameters, in particular the range step and the vertical grid spacing. It is shown that this dependence is related to the stationary phase point of the inverse Fourier integral. The errors of type (i) can be reduced by increasing the range step and/or decreasing the vertical grid spacing, but can be reduced much more efficiently by using an improved approximation for the forward Fourier integral. The errors of type (ii) can be reduced by using a numerical filter in the inverse Fourier integral. The errors of type (iii) can be reduced slightly by using an improved refraction factor. The errors of type (iv) can be reduced only by reducing the range step. The reduction of the four types of errors is illustrated for realistic test cases, by comparison with analytic solutions and results of the Crank–Nicholson PE (CNPE) method. Further, optimized values are presented for the parameters that determine the computational speed of the GFPE method. The computational speed difference between GFPE and CNPE is discussed in terms of numbers of floating point operations required by both methods.