We consider linear waves and baroclinic instability in a multilayer primitive equation model for ocean dynamics and thermodynamics. The model is a generalization to an arbitrary number of layers of Ripa's single-layer model with variable velocity shear and stratification. In addition to vary arbitrarily in horizontal position and time, the model's horizontal velocity and buoyancy fields are allowed to vary linearly with depth within each layer. The model enjoys several properties which make it very attractive. For instance: unlike slab layer models, in which all fields are set to be depth independent, the model represents explicitly within each layer the thermal-wind balance which dominates at low frequency; unlike homogeneous-density layered models, the model can incorporate thermodynamic processes, which are particularly important in the upper part of the ocean; and, in the absence of forcing and dissipation, the model has several integrals of motion including volume, mass, buoyancy variance, energy, and momentum. Our preliminary results on waves and instabilities suggest, in particular, that a configuration involving a few layers can be used as a basis for a quite accurate and, at the same time, numerically efficient ocean model.