Linear stability and the secondary flow pattern of the rectangular cell flow, $\ensuremath{\Psi}=\mathrm{sin}\mathrm{kx}\mathrm{sin}y$ $(0lkl\ensuremath{\infty}),$ are investigated in an infinitely long array of the $x$ direction $[(\ensuremath{-}\ensuremath{\infty},\ensuremath{\infty})\ifmmode\times\else\texttimes\fi{}[0,\ensuremath{\pi}]]$ or various finite $M$ arrays $([0,M\ensuremath{\pi}/k]\ifmmode\times\else\texttimes\fi{}[0,\ensuremath{\pi}])$ on the assumption of a stress-free boundary condition on the lateral walls. The numerical results of the eigenvalue problems on the infinite array show that a mode representing a global circulating vortex in the whole region $(\ensuremath{\psi}\ensuremath{\approx}\mathrm{sin}y)$ appears in the $y$-elongated cases $(kg1),$ which confirm the secondary flow observed in Tabeling et al. [J. Fluid Mech. 213, 511 (1990)], while a mode representing quasiperiodic arrays of counter-rotating vortices appears in the $x$-elongated cases $(kl1)$ at large critical Reynolds number. In large finite arrays the mode connected with those of the case $M=\ensuremath{\infty}$ appears for most cases while another (oscillatory) mode appears for vortices elongated in the $y$ direction. The parameter region of the oscillatory modes becomes wider when the system size $(M)$ becomes smaller. For a pair of counter-rotating vortices $(M=2)$ at the point ${k}_{0}$ between the regions of the two modes the critical Reynolds number takes an extreme large value. Analysis of a finite nonlinear system obtained by the Galerkin method shows the nonlinear saturation of the critical modes, though its results are in quantitative agreement with those of the linear stability in a limited region of $k.$