of $\eta_i$, $\epsilon_n$ and $k_{\parallel}$. The results are also contrasted to previous results obtained in models excluding parallel ion motion.

Results and discussion
======================

In the present section the solutions to the zonal flow dispersion relation are discussed and compared with previous findings. The zonal flow dispersion relation Eq. 18 is solved numerically. In the numerical analysis it is assumed that the mode coupling saturation level is reached for the drift wave turbulence [@a45] $$\begin{aligned}
\tilde{\phi}_{DW} = \frac{\gamma}{\omega_{\star}}\frac{1}{k_y L_n}.\end{aligned}$$ Of particular interest in the present work is the effect of toroidicity on the zonal flow generation. In this analysis this effect is included in the parameter $\epsilon_n$ where the slab like ITG mode is found for $\epsilon_n = 0$. The results are presented both in direct scalings and as the quality of zonal flows. The quality of zonal flows is defined as the ratio of the growth rate and the real frequency ($Q = \Omega^{IM}/\Omega^{RE}$).

First the effects of $\epsilon_n$ on zonal flow generation will be discussed. The quality of the zonal flow (Figure 1) with the parallel wave number ($k_{\parallel}^2$) as a parameter is investigated. The parameters are $\tau = 1$, $k_x = k_y = q_x = 0.3$, $\mu 