**Abstract** : This paper investigates the three-dimensional stability of a Lamb-Chaplygin columnar vertical vortex pair as a function of the vertical wavenumber k(z), horizontal Froude number F(h), Reynolds number Re and Schmidt number Sc. The horizontal Froude number F(h) (F(h) = U/NR, where U is the dipole travelling velocity, R the dipole radius and N the Brunt-Vaisala frequency) is varied in the range [0.033, 8[ and three set of Reynolds-Schmidt numbers are investigated: {Re = 10 000, Sc = 1}, {Re = 1000, Sc = 1}, {Re = 200, Sc = 637}. In the whole range of F(h) and Re, the dominant mode is always antisymmetric with respect to the middle plane between the vortices but its physical nature and properties change when F(h) is varied. An elliptic instability prevails for F(h) > 0.25, independently of the Reynolds number. It manifests itself by the bending of each vortex core in the opposite direction to the vortex periphery. The growth rate of the elliptic instability is reduced by stratification effects but its spatial structure is almost unaffected. In the range 0.2 < F(h) < 0.25, a continuous transition occurs from the elliptic instability to a different instability called zigzag instability. The transitional range F(hc) = 0.2-0.25 is in good agreement with the value F(h) = 0.22 at which the elliptic instability of an infinite uniform vortex is suppressed by the stratification. The zigzag instability dominates for F(h) = 0.2 and corresponds to a vertically modulated bending and twisting of the whole vortex pair. The experimental evidence for this zigzag instability in a strongly stratified fluid reported in the first part of this study (Billant and Chomaz 2000a) are therefore confirmed and extended. The numerically calculated wavelength and growth rate for low Reynolds number compare well with experimental measurements. The present numerical stability analysis fully agrees with the inviscid asymptotic analysis carried out in the second part of this investigation (Billant and Chomaz 2000b) for small Froude number F(h) and long wavelength. This confirms that the zigzag instability is related to the breaking of translational and rotational invariances. As predicted, the growth rate of the zigzag instability is observed to be self-similar with respect to the variable F(h)k(z), implying that the maximum growth rate is independent of F(h) while the most amplified dimensional wavenumber varies with N/U. The numerically computed eigenmode and dispersion relation are in striking agreement with the analytical results.