  A method is presented for constructing energy-conserving Galerkin
approximations in the vertical coordinate of the full quasigeostrophic model
with active surface buoyancy. The derivation generalizes the approach of Rocha
\emph{et al.} (2016) to allow for general bases. Details are then presented for
a specific set of bases: Legendre polynomials for potential vorticity and a
recombined Legendre basis from Shen (1994) for the streamfunction. The method
is tested in the context of linear baroclinic instability calculations, where
it is compared to the standard second-order finite-difference method and to a
Chebyshev collocation method. The Galerkin scheme is quite accurate even for a
small number of degrees of freedom $N$, and growth rates converge much more
quickly with increasing $N$ for the Galerkin scheme than for the
finite-difference scheme. The Galerkin scheme is at least as accurate as finite
differences and can in some cases achieve the same accuracy as the finite
difference scheme with ten times fewer degrees of freedom. The
energy-conserving Galerkin scheme is of comparable accuracy to the Chebyshev
collocation scheme in most linear stability calculations, but not in the Eady
problem where the Chebyshev scheme is significantly more accurate. Finally the
three methods are compared in the context of a simplified version of the
nonlinear equations: the two-surface model with zero potential vorticity. The
Chebyshev scheme is the most accurate, followed by the Galerkin scheme and then
the finite difference scheme. All three methods conserve energy with similar
accuracy, despite not having any a priori guarantee of energy conservation for
the Chebyshev scheme. Further nonlinear tests with non-zero potential vorticity
to assess the merits of the methods will be performed in a future work.
