A two layer version of HPM method was developed to model ocean flows above and below the thermocline, using a rigid lid approximation. In the figure is a simulation of a baroclinic instability, showing: the baroclinic and barotropic vorticities (top), the velocity fields in both layers (middle), the exchange of potential for kinetic energies (bottom left), and the error in total energy conservation (bottom right).