In this paper, new vertically discrete versions of the surface quasigeostrophic (SQG) model with two boundaries are formulated. For any number of partition levels, the equations of the discrete model are written in the form of conservation laws for two Lagrangian invariants, which have the meaning of buoyancy distributions at the horizontal boundaries of the fluid layer. The values of the invariants are expressed in terms of the values of the stream function at two internal levels and contain higher order elliptic operators. The use of discrete models greatly simplifies the solution of problems of the linear theory of hydrodynamic stability and provides high accuracy even with a small number of vertical discrete levels. Using the two-level version of the SQG model, which is similar to the classical two–layer Phillips model, we investigated the linear stability of jet flows induced by piecewise constant boundary distributions of buoyancy. For these flows, analytical expressions for the growth rate of perturbations have been obtained and it is shown that the most unstable perturbation has a wavelength of the order of the Rossby baroclinic radius of deformation. Flows with vertical shear induced by smooth and slowly varying boundary buoyancy distributions are also considered. The instability of these flows is found to be absolute, that is, independent of the velocity profile horizontal structure.