diff --git a/doc/hermes-2-manual.tex b/doc/hermes-2-manual.tex
index 89c6653..273f1ae 100644
--- a/doc/hermes-2-manual.tex
+++ b/doc/hermes-2-manual.tex
@@ -103,6 +103,29 @@ \section{Model equations}
   \pi_{ci} &\simeq& \frac{3m_i}{4p_iT_i}\left[0.20q_{||i}^2 - 0.085q_i^2\right] + 0.96\frac{p_i}{\nu_i}\Bigg\{\kappa\cdot\left[\mathbf{V}_E + \mathbf{V}_{di} + 1.61\frac{\mathbf{b}\times\nabla T_i}{B}\right] \nonumber \\
   && - \frac{2}{\sqrt{B}}\partial_{||}\left(\sqrt{B}V_{||i}\right) - \frac{1.42}{p_i\sqrt{B}}\partial_{||}\left(\sqrt{B}q_{||i}\right) - \frac{0.49q_{||i}}{p_i}\left(2.27\partial_{||}\ln T_i - \partial_{||}\ln p_i\right)\Bigg\} \label{eq:pi_ci}
 \end{eqnarray}
+The non-Boussinesq version of the vorticity is
+\begin{equation}
+\omega = \nabla\cdot\left[\frac{1}{B^2}\left(n\nabla_\perp\phi + \nabla_\perp p_i\right)\right]
+\end{equation}
+and the corresponding evolution equation is:
+\begin{subequations}
+\begin{eqnarray}
+%
+% Vorticity
+%
+  \frac{\partial\omega}{\partial t} &=& \nabla\cdot\left[\left(p_e + p_i\right)\nabla\times\frac{\mathbf{b}}{B}\right] \label{eq:vort-nB_diamag} \\
+  && + \nabla_{||}j_{||} \label{eq:vort-nB_jpar} \\
+  && - \nabla\cdot\Big[ \frac{1}{2B^2}\nabla_\perp\left(\mathbf{V}_{E\times B}\cdot\nabla p_i\right) + \frac{\omega}{2}\mathbf{V}_{E\times B}  +\frac{n}{2B^2}\nabla_\perp^2\phi \left(\mathbf{V}_{E\times B} + \mathbf{V}_{di}\right) \Big] \label{eq:vort-nB_c} \\
+  && + \nabla\cdot\left[\left(\frac{1}{2B^2}\mathbf{V}_{E\times B}\cdot\nabla n\right)\nabla_\perp\phi\right] \\
+  && + \nabla\cdot\left(\frac{3T_i}{10\tau_iB^2}\nabla_\perp\omega\right) \label{eq:vort-nB_perpvis} \\
+  && + \nabla\cdot\left[\frac{\pi_{ci}}{2}\nabla\times\frac{\mathbf{b}}{B} - \frac{1}{3}\frac{\mathbf{b}\times\nabla\pi_{ci}}{B}\right] \\
+  && + \nabla\cdot\left(\nu_A\nabla_\perp\left<\omega\right>\right) \\
+  && - \nabla\cdot\left(\frac{F_\perp}{nB^2}\nabla_\perp\phi\right)
+\end{eqnarray}
+\end{subequations}
+where the ion diamagnetic velocity $\mathbf{V}_{di} = \frac{\mathbf{b}\times\nabla p_i}{n B}$ uses the full density $n$.
+In addition $n_0$ is replaced by the full $n$ in the ion pressure equation~\ref{eq:ion-pressure} and the conserved energy~\ref{eq:conserved-energy} for the non-Boussinesq version.
+
 The electron pressure equation is:
 \begin{subequations}
 \begin{eqnarray}
@@ -144,7 +167,7 @@ \section{Model equations}
 && + \frac{5}{2}\nabla\cdot\left(\frac{T_i}{T_e}\frac{\rho_e^2}{\tau_e}\left[ \nabla_\perp p_e + \nabla_\perp p_i - \frac{3}{2}n\nabla_\perp T_e \right]\right) \\
 && - \frac{3T_i}{10\tau_iB^2}\nabla_\perp\omega\cdot\nabla\left(\phi + \frac{p_i}{n_0}\right) - \left[\frac{\pi_{ci}}{2}\nabla\times\frac{\mathbf{b}}{B} - \frac{1}{3}\frac{\mathbf{b}\times\nabla\pi_{ci}}{B}\right]\cdot\nabla\left(\phi + \frac{p_i}{n_0}\right) \\
 && + S_{pi} - Q_i + W_i
-\end{eqnarray}
+\end{eqnarray}\label{eq:ion-pressure}
 \end{subequations}
 where the terms are a) cross-field E$\times$B and magnetic drifts including compression; b) parallel flow including compression;
 c) energy exchange with diamagnetic flows; d) parallel viscous heating; e) parallel and perpendicular collisional heat conduction; f) collisional resistive drift; g) heating due to perpendicular viscosity; h)
@@ -207,6 +230,7 @@ \section{Model equations}
 The plasma model is cast in conservative form, conserving particle number, and an energy
 \begin{equation}
 E = \int dv \frac{m_in_0}{2}\left|\frac{\nabla_\perp\phi}{B} + \frac{\nabla_\perp p_i}{en_0B}\right|^2 + \frac{1}{2}m_inV_{||i}^2 + \frac{3}{2}p_e + \frac{3}{2}p_i + \frac{1}{4}\beta_e\left|\nabla_\perp\psi\right|^2 + \frac{m_e}{m_i}\frac{1}{2}\frac{j_{||}^2}{n}
+\label{eq:conserved-energy}
 \end{equation}
 where the terms correspond to the ion perpendicular kinetic energy ($E\times B$ and diamagnetic flow), ion parallel kinetic energy, electron and ion thermal energy, and electromagnetic field energy. Details are given in section~\ref{sec:conservation}. Differential operators are discretised using flux-conservative Finite Volume methods, details of which are given in section~\ref{sec:operators}. 
 
@@ -224,7 +248,7 @@ \subsection{Input options}
 These switches control the terms included in the plasma
 evolution equations. They are grouped so that switching on
 or off a term maintains energy conservation, though the
-form of the conserved energy may change. See section~\ref{sec:rhsfunc} for details.
+form of the conserved energy may change. See section~\ref{sec:conservation} for details.
 \begin{center}
 \begin{tabular}{l c c}
   Setting & Default & Description \\
diff --git a/external/BOUT-dev b/external/BOUT-dev
index 639d581..2b2223b 160000
--- a/external/BOUT-dev
+++ b/external/BOUT-dev
@@ -1 +1 @@
-Subproject commit 639d581279fab03ad7906f191c6c12c2bc437fbc
+Subproject commit 2b2223b51cb4c1f2013950a4586619f8a917d32c
diff --git a/src/hermes-2.cxx b/src/hermes-2.cxx
index f9c9874..5c77219 100644
--- a/src/hermes-2.cxx
+++ b/src/hermes-2.cxx
@@ -905,7 +905,14 @@ int Hermes::init(bool restarting) {
         // Use older Laplacian solver
         phiSolver = Laplacian::create(&opt["phiSolver"]);
         // Set coefficients for Boussinesq solve
-        phiSolver->setCoefC(1. / SQ(coord->Bxy));
+        if (boussinesq) {
+          phiSolver->setCoefC(1. / SQ(coord->Bxy));
+        } else {
+          if (not phiSolver.uses3DCoefs()) {
+            throw BoutException("Non-Boussinesq solve requires phiSolver to support 3D "
+                                "coefficients.")
+          }
+        }
       }
     }
     phi = 0.0;
@@ -1335,7 +1342,51 @@ int Hermes::rhs(BoutReal t) {
         ////////////////////////////////////////////
         // Non-Boussinesq
         //
-        throw BoutException("Non-Boussinesq not implemented yet");
+        if (split_n0) {
+          throw BoutException("split_n0 not compaible with non-Boussinesq solve");
+        }
+        if (newXZsolver) {
+          throw BoutException("Non-Boussinesq not implemented yet for newXZsolver");
+        }
+        // Update boundary conditions.
+        // The INVERT_SET flag takes the value in the guard (boundary) cell and sets the
+        // boundary between cells to this value.  This shift by 1/2 grid cell is
+        // important.
+
+        if (mesh->firstX()) {
+          for (int j = mesh->ystart; j <= mesh->yend; j++) {
+            for (int k = 0; k < mesh->LocalNz; k++) {
+              // Average phi + Pi at the boundary, and set the boundary cell
+              // to this value. The phi solver will then put the value back
+              // onto the cell mid-point
+              phi_boundary3d(mesh->xstart - 1, j, k) =
+                  0.5
+                  * (phi_boundary3d(mesh->xstart - 1, j, k) +
+                     phi_boundary3d(mesh->xstart, j, k));
+            }
+          }
+        }
+
+        if (mesh->lastX()) {
+          for (int j = mesh->ystart; j <= mesh->yend; j++) {
+            for (int k = 0; k < mesh->LocalNz; k++) {
+              phi_boundary3d(mesh->xend + 1, j, k) =
+                  0.5
+                  * (phi_boundary3d(mesh->xend + 1, j, k) +
+                     phi_boundary3d(mesh->xend, j, k));
+            }
+          }
+        }
+        auto n_B2 = Ne/SQ(coord->Bxy);
+        // Use older Laplacian solver
+        // Need to subtract Div(B^-2 Grad_perp(pi)) from rhs for non-Boussinesq solve, as
+        // coefficients of Grad_perp(phi) and Grad_perp(Pi) terms are not the same in
+        // non-Boussinesq vorticity
+        phiSolver->setCoefC(n_B2); // Set when initialised
+        // phi_boundary3d here is equivalent to withBoundary(phi, phi_boundary3d) because
+        // of the way phi_boundary3d is initialised
+        phi = phiSolver->solve((Vort - FV::Div_a_Laplace_perp(1/SQ(coord->Bxy),Pi))/n_B2,
+                               phi_boundary3d);
       }
     }
     phi.applyBoundary(t);
@@ -2805,35 +2856,46 @@ int Hermes::rhs(BoutReal t) {
     }
 
     // Advection of vorticity by ExB
-    if (boussinesq) {
-      TRACE("Vort:boussinesq");
-      // Using the Boussinesq approximation
-
-      ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_bndry_flux,
-                                         poloidal_flows);
+    ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_bndry_flux,
+                                       poloidal_flows);
+    // V_ExB dot Grad(Pi)
+    Field3D vEdotGradPi = bracket(phi, Pi, BRACKET_ARAKAWA);
+    vEdotGradPi.applyBoundary("free_o2");
 
-      // V_ExB dot Grad(Pi)
-      Field3D vEdotGradPi = bracket(phi, Pi, BRACKET_ARAKAWA);
-      vEdotGradPi.applyBoundary("free_o2");
-      // delp2(phi) term
-      Field3D DelpPhi_2B2 = 0.5 * Delp2(phi) / SQ(coord->Bxy);
-      DelpPhi_2B2.applyBoundary("free_o2");
+    // delp2(phi) term
+    Field3D DelpPhi_2B2 = 0.5 * Delp2(phi) / SQ(coord->Bxy);
+    DelpPhi_2B2.applyBoundary("free_o2");
 
-      mesh->communicate(vEdotGradPi, DelpPhi_2B2);
+    mesh->communicate(vEdotGradPi, DelpPhi_2B2);
 
-      ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / SQ(coord->Bxy), vEdotGradPi);
+    ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / SQ(coord->Bxy), vEdotGradPi);
 
+    if (boussinesq) {
+      TRACE("Vort:boussinesq_terms");
       if (j_pol_terms){
-	ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi, vort_bndry_flux,
+        // Using the Boussinesq approximation
+        ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi, vort_bndry_flux,
+                                           poloidal_flows);
+      }
+    } else {
+      TRACE("Vort:non-boussinesq_terms");
+      if (j_pol_terms){
+        // When the Boussinesq approximation is not made,
+        // then the changing ion density introduces a number
+        // of other terms.
+
+	ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Ne*DelpPhi_2B2, phi, vort_bndry_flux,
+					   poloidal_flows);
+	ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, Pi, vort_bndry_flux,
 					   poloidal_flows);
       }
 
-    } else {
-      // When the Boussinesq approximation is not made,
-      // then the changing ion density introduces a number
-      // of other terms.
+      // V_ExB dot Grad(n)
+      Field3D vEdotGradNe = bracket(phi, Ne, BRACKET_ARAKAWA);
+      vEdotGradNe.applyBoundary("free_o2");
+      mesh->communicate(vEdotGradNe);
 
-      throw BoutException("Hot ion non-Boussinesq not implemented yet\n");
+      ddt(Vort) += FV::Div_a_Laplace_perp(0.5 * vEdotGradNe / SQ(coord->Bxy), phi);
     }
 
     if (classical_diffusion) {
@@ -3352,7 +3414,11 @@ int Hermes::rhs(BoutReal t) {
     // in the vorticity equation
     ddt(Pi) -= j_diamag_scale * (2. / 3) * Pi * (Curlb_B * Grad(phi));
 
-    ddt(Pi) += j_diamag_scale * Pi * Div((Pe + Pi) * Curlb_B);
+    if (boussinesq) {
+      ddt(Pi) += j_diamag_scale * Pi * Div((Pe + Pi) * Curlb_B);
+    } else {
+      ddt(Pi) += j_diamag_scale * Pi / Nelim * Div((Pe + Pi) * Curlb_B);
+    }
   }
 
   if (j_par) {