diff --git a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py b/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py
index 5ce8e5e49..7fad18ca5 100644
--- a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py
+++ b/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py
@@ -26,7 +26,8 @@
     logger, SUPGOptions, Perturbation, CompressibleParameters,
     CompressibleEulerEquations, HydrostaticCompressibleEulerEquations,
     compressible_hydrostatic_balance, RungeKuttaFormulation, CompressibleSolver,
-    SubcyclingOptions, hydrostatic_parameters
+    SubcyclingOptions, LinearVariationalProblem, LinearVariationalSolver,
+    TestFunction, TrialFunction, dx
 )
 
 skamarock_klemp_nonhydrostatic_defaults = {
@@ -67,6 +68,8 @@ def skamarock_klemp_nonhydrostatic(
     # ------------------------------------------------------------------------ #
 
     element_order = 1
+    alpha = 0.5
+    u_eqn_type = 'vector_advection_form'
 
     # ------------------------------------------------------------------------ #
     # Set up model objects
@@ -78,9 +81,11 @@ def skamarock_klemp_nonhydrostatic(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    parameters = CompressibleParameters(mesh)
+    parameters = CompressibleParameters()
     if hydrostatic:
-        eqns = HydrostaticCompressibleEulerEquations(domain, parameters)
+        eqns = HydrostaticCompressibleEulerEquations(
+            domain, parameters, u_transport_option=u_eqn_type
+        )
     else:
         eqns = CompressibleEulerEquations(domain, parameters)
 
@@ -135,17 +140,14 @@ def skamarock_klemp_nonhydrostatic(
 
     # Linear solver
     if hydrostatic:
-        linear_solver = CompressibleSolver(
-            eqns, solver_parameters=hydrostatic_parameters,
-            overwrite_solver_parameters=True
-        )
+        linear_solver = CompressibleSolver(eqns, alpha=alpha)
     else:
         linear_solver = CompressibleSolver(eqns)
 
     # Time stepper
     stepper = SemiImplicitQuasiNewton(
         eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver
+        linear_solver=linear_solver, alpha=alpha, num_outer=2, num_inner=2
     )
 
     # ------------------------------------------------------------------------ #
@@ -175,12 +177,23 @@ def skamarock_klemp_nonhydrostatic(
     # Calculate hydrostatic exner
     compressible_hydrostatic_balance(eqns, theta_b, rho_b)
 
+    # Define initial theta
     theta_pert = (
         deltaTheta * sin(pi*z/domain_height)
         / (1 + (x - domain_width/2)**2 / pert_width**2)
     )
     theta0.interpolate(theta_b + theta_pert)
-    rho0.assign(rho_b)
+
+    # find perturbed rho
+    gamma = TestFunction(Vr)
+    rho_trial = TrialFunction(Vr)
+    dx_qp = dx(degree=domain.max_quad_degree)
+    lhs = gamma * rho_trial * dx_qp
+    rhs = gamma * (rho_b * theta_b / theta0) * dx_qp
+    rho_problem = LinearVariationalProblem(lhs, rhs, rho0)
+    rho_solver = LinearVariationalSolver(rho_problem)
+    rho_solver.solve()
+
     u0.project(as_vector([wind_initial, 0.0]))
 
     stepper.set_reference_profiles([('rho', rho_b), ('theta', theta_b)])
diff --git a/examples/compressible_euler/test_compressible_euler_examples.py b/examples/compressible_euler/test_compressible_euler_examples.py
index d74d5dade..b060f2403 100644
--- a/examples/compressible_euler/test_compressible_euler_examples.py
+++ b/examples/compressible_euler/test_compressible_euler_examples.py
@@ -65,8 +65,6 @@ def test_skamarock_klemp_nonhydrostatic_parallel():
     test_skamarock_klemp_nonhydrostatic()
 
 
-# Hydrostatic equations not currently working
-@pytest.mark.xfail
 def test_hyd_switch_skamarock_klemp_nonhydrostatic():
     from skamarock_klemp_nonhydrostatic import skamarock_klemp_nonhydrostatic
     test_name = 'hyd_switch_skamarock_klemp_nonhydrostatic'
@@ -81,8 +79,6 @@ def test_hyd_switch_skamarock_klemp_nonhydrostatic():
     )
 
 
-# Hydrostatic equations not currently working
-@pytest.mark.xfail
 @pytest.mark.parallel(nprocs=2)
 def test_hyd_switch_skamarock_klemp_nonhydrostatic_parallel():
     test_hyd_switch_skamarock_klemp_nonhydrostatic()
diff --git a/gusto/core/labels.py b/gusto/core/labels.py
index b49bbefad..3e4aec99a 100644
--- a/gusto/core/labels.py
+++ b/gusto/core/labels.py
@@ -114,3 +114,5 @@ def __call__(self, target, value=None):
 hydrostatic = DynamicsLabel("hydrostatic")
 incompressible = DynamicsLabel("incompressible")
 sponge = DynamicsLabel("sponge")
+horizontal_prognostic = Label("horizontal_prognostic")
+vertical_prognostic = Label("vertical_prognostic")
diff --git a/gusto/equations/compressible_euler_equations.py b/gusto/equations/compressible_euler_equations.py
index d512432cf..bbb14cfa3 100644
--- a/gusto/equations/compressible_euler_equations.py
+++ b/gusto/equations/compressible_euler_equations.py
@@ -4,10 +4,11 @@
     sin, pi, inner, dx, div, cross, FunctionSpace, FacetNormal, jump, avg, dS_v,
     conditional, SpatialCoordinate, split, Constant, as_vector
 )
-from firedrake.fml import subject, replace_subject
+from firedrake.fml import subject, drop, keep
 from gusto.core.labels import (
     time_derivative, transport, prognostic, hydrostatic, linearisation,
-    pressure_gradient, coriolis, gravity, sponge
+    pressure_gradient, coriolis, gravity, sponge, implicit,
+    horizontal_prognostic
 )
 from gusto.equations.thermodynamics import exner_pressure
 from gusto.equations.common_forms import (
@@ -335,49 +336,62 @@ def __init__(self, domain, parameters, sponge_options=None,
                          no_normal_flow_bc_ids=no_normal_flow_bc_ids,
                          active_tracers=active_tracers)
 
-        # Replace
-        self.residual = self.residual.label_map(
-            lambda t: t.has_label(time_derivative),
-            map_if_true=lambda t: self.hydrostatic_projection(t, 'u')
-        )
-
-        # Add an extra hydrostatic term
         u_idx = self.field_names.index('u')
         u = split(self.X)[u_idx]
+        w = self.tests[u_idx]
         k = self.domain.k
-        self.residual += hydrostatic(
-            subject(
-                prognostic(
-                    -inner(k, self.tests[u_idx]) * inner(k, u) * dx, "u"),
-                self.X
-            )
+        u_vert = k*inner(k, u)
+        u_hori = u - u_vert
+
+        # -------------------------------------------------------------------- #
+        # Add hydrostatic term
+        # -------------------------------------------------------------------- #
+        # Term to appear in both explicit and implicit forcings
+        # For the explicit forcing, this will cancel out the u^n part of the
+        # time derivative
+        self.residual += hydrostatic(subject(prognostic(
+            inner(k, w)*inner(k, u)/domain.dt*dx, 'u'), self.X)
         )
 
-    def hydrostatic_projection(self, term, field_name):
-        """
-        Performs the 'hydrostatic' projection.
+        # Term that appears only in implicit forcing
+        # For the implicit forcing, in combination with the term above, the
+        # u^{n+1} term will be cancelled out
+        self.residual += implicit(hydrostatic(subject(prognostic(
+            Constant(-1.0)*inner(k, w)*inner(k, u)/domain.dt*dx, 'u'), self.X))
+        )
 
-        Takes a term involving a vector prognostic variable and replaces the
-        prognostic with only its horizontal components. It also adds the
-        'hydrostatic' label to that term.
+        # # Add Euler-Poincare term
+        # self.residual += hydrostatic(subject(prognostic(
+        #     Constant(0.5)*div(w)*inner(u_hori, u)*dx, 'u'), self.X)
+        # )
 
-        Args:
-            term (:class:`Term`): the term to perform the projection upon.
-            field_name (str): the prognostic field to make hydrostatic.
+        # -------------------------------------------------------------------- #
+        # Only transport horizontal wind
+        # -------------------------------------------------------------------- #
+        # Drop wind transport term, it needs replacing with a version that
+        # includes only the horizontal components for the transported field
+        self.residual = self.residual.label_map(
+            lambda t: t.has_label(transport) and t.get(prognostic) == 'u',
+            map_if_true=drop,
+            map_if_false=keep
+        )
 
-        Returns:
-            :class:`LabelledForm`: the labelled form containing the new term.
-        """
+        # u_term = prognostic(
+        #     advection_equation_circulation_form(domain, w, u_hori, u), 'u'
+        # )
 
-        f_idx = self.field_names.index(field_name)
-        k = self.domain.k
-        X = term.get(subject)
-        field = split(X)[f_idx]
-
-        new_subj = field - inner(field, k) * k
-        # In one step:
-        # - set up the replace_subject routine (which returns a function)
-        # - call that function on the supplied `term` argument,
-        #   to replace the subject with the new hydrostatic subject
-        # - add the hydrostatic label
-        return replace_subject(new_subj, old_idx=f_idx)(term)
+        # Velocity transport term -- depends on formulation
+        if u_transport_option == "vector_invariant_form":
+            u_term = prognostic(vector_invariant_form(domain, w, u_hori, u), 'u')
+        elif u_transport_option == "vector_advection_form":
+            u_term = prognostic(advection_form(w, u_hori, u), 'u')
+        elif u_transport_option == "circulation_form":
+            circ_form = prognostic(
+                advection_equation_circulation_form(domain, w, u_hori, u), 'u'
+            )
+            ke_form = prognostic(kinetic_energy_form(w, u_hori, u), 'u')
+            u_term = ke_form + circ_form
+        else:
+            raise ValueError("Invalid u_transport_option: %s" % u_transport_option)
+
+        self.residual += horizontal_prognostic(subject(u_term, self.X))
diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py
index 213d09df2..57775ab8e 100644
--- a/gusto/solvers/linear_solvers.py
+++ b/gusto/solvers/linear_solvers.py
@@ -23,8 +23,8 @@
     logger, DEBUG, logging_ksp_monitor_true_residual,
     attach_custom_monitor
 )
-from gusto.core.labels import linearisation, time_derivative, hydrostatic
-from gusto.equations import thermodynamics
+from gusto.core.labels import linearisation, time_derivative
+from gusto.equations import thermodynamics, HydrostaticCompressibleEulerEquations
 from gusto.recovery.recovery_kernels import AverageWeightings, AverageKernel
 from abc import ABCMeta, abstractmethod, abstractproperty
 
@@ -135,7 +135,7 @@ class CompressibleSolver(TimesteppingSolver):
     (3) Reconstruct theta
     """
 
-    solver_parameters = {'mat_type': 'matfree',
+    hybrid_parameters = {'mat_type': 'matfree',
                          'ksp_type': 'preonly',
                          'pc_type': 'python',
                          'pc_python_type': 'firedrake.SCPC',
@@ -152,8 +152,33 @@ class CompressibleSolver(TimesteppingSolver):
                                                            'pc_type': 'bjacobi',
                                                            'sub_pc_type': 'ilu'}}}
 
-    def __init__(self, equations, alpha=0.5, tau_values=None,
-                 solver_parameters=None, overwrite_solver_parameters=False):
+    full_parameters = {
+        'pc_type': 'fieldsplit',
+        'pc_fieldsplit_type': 'schur',
+        'ksp_type': 'gmres',
+        'ksp_max_it': 100,
+        'ksp_gmres_restart': 50,
+        'pc_fieldsplit_schur_fact_type': 'FULL',
+        'pc_fieldsplit_schur_precondition': 'selfp',
+        'fieldsplit_0': {'ksp_type': 'preonly',
+                         'pc_type': 'bjacobi',
+                         'sub_pc_type': 'ilu'},
+        'fieldsplit_1': {'ksp_type': 'preonly',
+                         'pc_type': 'gamg',
+                         'mg_levels': {'ksp_type': 'chebyshev',
+                                       'ksp_chebyshev_esteig': True,
+                                       'ksp_max_it': 1,
+                                       'pc_type': 'bjacobi',
+                                       'sub_pc_type': 'ilu'}}
+    }
+
+    solver_parameters = None
+
+    def __init__(
+            self, equations, alpha=0.5, tau_values=None,
+            solver_parameters=None, overwrite_solver_parameters=False,
+            formulation='hybridized'
+    ):
         """
         Args:
             equations (:class:`PrognosticEquation`): the model's equation.
@@ -168,9 +193,20 @@ def __init__(self, equations, alpha=0.5, tau_values=None,
                 `solver_parameters` that have been passed in. If False then
                 update the default parameters with the `solver_parameters`
                 passed in. Defaults to False.
+            formulation (str, optional): the formulation to use. Valid options
+                are 'hybridized' and 'full'. Defaults to 'hybridized'.
         """
+        assert formulation in ['hybridized', 'full'], \
+            f'Invalid solver formulation: {formulation}'
+
         self.equations = equations
         self.quadrature_degree = equations.domain.max_quad_degree
+        self.formulation = formulation
+
+        if formulation == 'hybridized':
+            self.solver_parameters = self.hybrid_parameters
+        elif formulation == 'full':
+            self.solver_parameters = self.full_parameters
 
         super().__init__(equations, alpha, tau_values, solver_parameters,
                          overwrite_solver_parameters)
@@ -178,43 +214,74 @@ def __init__(self, equations, alpha=0.5, tau_values=None,
     @timed_function("Gusto:SolverSetup")
     def _setup_solver(self):
 
+        # Declare constants and relaxation parameters --------------------------
         equations = self.equations
         cp = equations.parameters.cp
         dt = self.dt
+
         # Set relaxation parameters. If an alternative has not been given, set
         # to semi-implicit off-centering factor
         beta_u = dt*self.tau_values.get("u", self.alpha)
         beta_t = dt*self.tau_values.get("theta", self.alpha)
         beta_r = dt*self.tau_values.get("rho", self.alpha)
 
+        # Specify degree for some terms as estimated degree is too large -------
+        dx_qp = dx(degree=(equations.domain.max_quad_degree))
+        dS_v_qp = dS_v(degree=(equations.domain.max_quad_degree))
+        dS_h_qp = dS_h(degree=(equations.domain.max_quad_degree))
+        ds_v_qp = ds_v(degree=(equations.domain.max_quad_degree))
+        ds_tb_qp = (ds_t(degree=(equations.domain.max_quad_degree))
+                    + ds_b(degree=(equations.domain.max_quad_degree)))
+
+        # Split up the rhs vector (symbolically) -------------------------------
+        self.xrhs = Function(self.equations.function_space)
+        u_in, rho_in, theta_in = split(self.xrhs)[0:3]
+
+        # Get the function spaces ----------------------------------------------
         Vu = equations.domain.spaces("HDiv")
-        Vu_broken = FunctionSpace(equations.domain.mesh, BrokenElement(Vu.ufl_element()))
         Vtheta = equations.domain.spaces("theta")
         Vrho = equations.domain.spaces("DG")
 
-        h_deg = Vrho.ufl_element().degree()[0]
-        v_deg = Vrho.ufl_element().degree()[1]
-        Vtrace = FunctionSpace(equations.domain.mesh, "HDiv Trace", degree=(h_deg, v_deg))
+        if self.formulation == 'hybridized':
+            h_deg = Vrho.ufl_element().degree()[0]
+            v_deg = Vrho.ufl_element().degree()[1]
+            Vtrace = FunctionSpace(equations.domain.mesh, "HDiv Trace", degree=(h_deg, v_deg))
+            Vu_broken = FunctionSpace(equations.domain.mesh, BrokenElement(Vu.ufl_element()))
 
-        # Split up the rhs vector (symbolically)
-        self.xrhs = Function(self.equations.function_space)
-        u_in, rho_in, theta_in = split(self.xrhs)[0:3]
+            # Build the function space for "broken" u, rho, and pressure trace
+            M = MixedFunctionSpace((Vu_broken, Vrho, Vtrace))
+            w, phi, dl = TestFunctions(M)
+            u, rho, l0 = TrialFunctions(M)
 
-        # Build the function space for "broken" u, rho, and pressure trace
-        M = MixedFunctionSpace((Vu_broken, Vrho, Vtrace))
-        w, phi, dl = TestFunctions(M)
-        u, rho, l0 = TrialFunctions(M)
-
-        n = FacetNormal(equations.domain.mesh)
+        elif self.formulation == 'full':
+            # Mixed function space is just for velocity and density
+            M = MixedFunctionSpace((Vu, Vrho))
+            w, phi = TestFunctions(M)
+            u, rho = TrialFunctions(M)
 
-        # Get background fields
+        # Get background fields ------------------------------------------------
         _, rhobar, thetabar = split(equations.X_ref)[0:3]
         exnerbar = thermodynamics.exner_pressure(equations.parameters, rhobar, thetabar)
         exnerbar_rho = thermodynamics.dexner_drho(equations.parameters, rhobar, thetabar)
         exnerbar_theta = thermodynamics.dexner_dtheta(equations.parameters, rhobar, thetabar)
 
-        # Analytical (approximate) elimination of theta
+        # Set up elimination of theta  =========================================
         k = equations.domain.k             # Upward pointing unit vector
+        n = FacetNormal(equations.domain.mesh)
+
+        # Vertical projection
+        def V(u):
+            return k*inner(u, k)
+
+        # Hydrostatic projection
+        h_project = lambda u: u - k*inner(u, k)
+
+        if isinstance(self.equations, HydrostaticCompressibleEulerEquations):
+            u_mass = inner(w, (h_project(u) - u_in))*dx
+        else:
+            u_mass = inner(w, (u - u_in))*dx
+
+        # Analytical (approximate) elimination of theta
         theta = -dot(k, u)*dot(k, grad(thetabar))*beta_t + theta_in
 
         # Only include theta' (rather than exner') in the vertical
@@ -224,21 +291,6 @@ def _setup_solver(self):
         # for linear perturbations)
         exner = exnerbar_theta*theta + exnerbar_rho*rho
 
-        # Vertical projection
-        def V(u):
-            return k*inner(u, k)
-
-        # hydrostatic projection
-        h_project = lambda u: u - k*inner(u, k)
-
-        # Specify degree for some terms as estimated degree is too large
-        dx_qp = dx(degree=(equations.domain.max_quad_degree))
-        dS_v_qp = dS_v(degree=(equations.domain.max_quad_degree))
-        dS_h_qp = dS_h(degree=(equations.domain.max_quad_degree))
-        ds_v_qp = ds_v(degree=(equations.domain.max_quad_degree))
-        ds_tb_qp = (ds_t(degree=(equations.domain.max_quad_degree))
-                    + ds_b(degree=(equations.domain.max_quad_degree)))
-
         # Add effect of density of water upon theta, using moisture reference profiles
         # TODO: Explore if this is the right thing to do for the linear problem
         if equations.active_tracers is not None:
@@ -258,65 +310,92 @@ def V(u):
             theta_w = theta
             thetabar_w = thetabar
 
-        _l0 = TrialFunction(Vtrace)
-        _dl = TestFunction(Vtrace)
-        a_tr = _dl('+')*_l0('+')*(dS_v_qp + dS_h_qp) + _dl*_l0*ds_v_qp + _dl*_l0*ds_tb_qp
-
-        def L_tr(f):
-            return _dl('+')*avg(f)*(dS_v_qp + dS_h_qp) + _dl*f*ds_v_qp + _dl*f*ds_tb_qp
-
-        cg_ilu_parameters = {'ksp_type': 'cg',
-                             'pc_type': 'bjacobi',
-                             'sub_pc_type': 'ilu'}
-
-        # Project field averages into functions on the trace space
-        rhobar_avg = Function(Vtrace)
-        exnerbar_avg = Function(Vtrace)
-
-        rho_avg_prb = LinearVariationalProblem(a_tr, L_tr(rhobar), rhobar_avg,
-                                               constant_jacobian=True)
-        exner_avg_prb = LinearVariationalProblem(a_tr, L_tr(exnerbar), exnerbar_avg,
-                                                 constant_jacobian=True)
-
-        self.rho_avg_solver = LinearVariationalSolver(rho_avg_prb,
-                                                      solver_parameters=cg_ilu_parameters,
-                                                      options_prefix='rhobar_avg_solver')
-        self.exner_avg_solver = LinearVariationalSolver(exner_avg_prb,
-                                                        solver_parameters=cg_ilu_parameters,
-                                                        options_prefix='exnerbar_avg_solver')
-
-        # "broken" u, rho, and trace system
-        # NOTE: no ds_v integrals since equations are defined on
-        # a periodic (or sphere) base mesh.
-        if any([t.has_label(hydrostatic) for t in self.equations.residual]):
-            u_mass = inner(w, (h_project(u) - u_in))*dx
-        else:
-            u_mass = inner(w, (u - u_in))*dx
-
-        eqn = (
-            # momentum equation
-            u_mass
-            - beta_u*cp*div(theta_w*V(w))*exnerbar*dx_qp
-            # following does nothing but is preserved in the comments
-            # to remind us why (because V(w) is purely vertical).
-            # + beta*cp*jump(theta_w*V(w), n=n)*exnerbar_avg('+')*dS_v_qp
-            + beta_u*cp*jump(theta_w*V(w), n=n)*exnerbar_avg('+')*dS_h_qp
-            + beta_u*cp*dot(theta_w*V(w), n)*exnerbar_avg*ds_tb_qp
-            - beta_u*cp*div(thetabar_w*w)*exner*dx_qp
-            # trace terms appearing after integrating momentum equation
-            + beta_u*cp*jump(thetabar_w*w, n=n)*l0('+')*(dS_v_qp + dS_h_qp)
-            + beta_u*cp*dot(thetabar_w*w, n)*l0*(ds_tb_qp + ds_v_qp)
-            # mass continuity equation
-            + (phi*(rho - rho_in) - beta_r*inner(grad(phi), u)*rhobar)*dx
-            + beta_r*jump(phi*u, n=n)*rhobar_avg('+')*(dS_v + dS_h)
-            # term added because u.n=0 is enforced weakly via the traces
-            + beta_r*phi*dot(u, n)*rhobar_avg*(ds_tb + ds_v)
-            # constraint equation to enforce continuity of the velocity
-            # through the interior facets and weakly impose the no-slip
-            # condition
-            + dl('+')*jump(u, n=n)*(dS_v + dS_h)
-            + dl*dot(u, n)*(ds_t + ds_b + ds_v)
-        )
+        cg_ilu_parameters = {
+            'ksp_type': 'cg',
+            'pc_type': 'bjacobi',
+            'sub_pc_type': 'ilu'
+        }
+
+        # Hybridization formulation ============================================
+        if self.formulation == 'hybridized':
+            _l0 = TrialFunction(Vtrace)
+            _dl = TestFunction(Vtrace)
+            a_tr = _dl('+')*_l0('+')*(dS_v_qp + dS_h_qp) + _dl*_l0*ds_v_qp + _dl*_l0*ds_tb_qp
+
+            def L_tr(f):
+                return _dl('+')*avg(f)*(dS_v_qp + dS_h_qp) + _dl*f*ds_v_qp + _dl*f*ds_tb_qp
+
+            # Project field averages into functions on the trace space
+            rhobar_avg = Function(Vtrace)
+            exnerbar_avg = Function(Vtrace)
+
+            rho_avg_prb = LinearVariationalProblem(
+                a_tr, L_tr(rhobar), rhobar_avg, constant_jacobian=True
+            )
+            exner_avg_prb = LinearVariationalProblem(
+                a_tr, L_tr(exnerbar), exnerbar_avg, constant_jacobian=True
+            )
+
+            self.rho_avg_solver = LinearVariationalSolver(
+                rho_avg_prb, solver_parameters=cg_ilu_parameters,
+                options_prefix='rhobar_avg_solver'
+            )
+            self.exner_avg_solver = LinearVariationalSolver(
+                exner_avg_prb, solver_parameters=cg_ilu_parameters,
+                options_prefix='exnerbar_avg_solver'
+            )
+
+            # Function for the hybridized solutions
+            self.urhol0 = Function(M)
+
+            # "broken" u, rho, and trace system
+            # NOTE: no ds_v integrals since equations are defined on
+            # a periodic (or sphere) base mesh.
+
+            eqn = (
+                # momentum equation
+                u_mass
+                - beta_u*cp*div(theta_w*V(w))*exnerbar*dx_qp
+                # following does nothing but is preserved in the comments
+                # to remind us why (because V(w) is purely vertical).
+                # + beta*cp*jump(theta_w*V(w), n=n)*exnerbar_avg('+')*dS_v_qp
+                + beta_u*cp*jump(theta_w*V(w), n=n)*exnerbar_avg('+')*dS_h_qp
+                + beta_u*cp*dot(theta_w*V(w), n)*exnerbar_avg*ds_tb_qp
+                - beta_u*cp*div(thetabar_w*w)*exner*dx_qp
+                # trace terms appearing after integrating momentum equation
+                + beta_u*cp*jump(thetabar_w*w, n=n)*l0('+')*(dS_v_qp + dS_h_qp)
+                + beta_u*cp*dot(thetabar_w*w, n)*l0*(ds_tb_qp + ds_v_qp)
+                # mass continuity equation
+                + (phi*(rho - rho_in) - beta_r*inner(grad(phi), u)*rhobar)*dx
+                + beta_r*jump(phi*u, n=n)*rhobar_avg('+')*(dS_v + dS_h)
+                # term added because u.n=0 is enforced weakly via the traces
+                + beta_r*phi*dot(u, n)*rhobar_avg*(ds_tb + ds_v)
+                # constraint equation to enforce continuity of the velocity
+                # through the interior facets and weakly impose the no-slip
+                # condition
+                + dl('+')*jump(u, n=n)*(dS_v + dS_h)
+                + dl*dot(u, n)*(ds_t + ds_b + ds_v)
+            )
+
+        # Full formulation =====================================================
+        elif self.formulation == 'full':
+            # Function to store result of u-rho solve
+            self.urho = Function(M)
+            one = Constant(1.0)
+            eqn = (
+                u_mass
+                - beta_u*cp*div(theta_w*V(w))*exnerbar*dx_qp
+                # following does nothing but is preserved in the comments
+                # to remind us why (because V(w) is purely vertical.
+                # + beta*cp*jump(theta*V(w),n)*avg(pibar)*dS_v_qp
+                - beta_u*cp*div(thetabar_w*w)*exner*dx_qp
+                - (one - beta_u)*cp*div(thetabar_w*V(w))*exner*dx_qp
+                + beta_u*cp*jump(thetabar_w*w, n)*avg(exner)*dS_v_qp
+                + (phi*(rho - rho_in) - beta_r*inner(grad(phi), u)*rhobar)*dx
+                + beta_r*jump(phi*u, n)*avg(rhobar)*(dS_v + dS_h)
+            )
+
+        # Add additional terms =================================================
         # TODO: can we get this term using FML?
         # contribution of the sponge term
         if hasattr(self.equations, "mu"):
@@ -329,68 +408,99 @@ def L_tr(f):
         aeqn = lhs(eqn)
         Leqn = rhs(eqn)
 
-        # Function for the hybridized solutions
-        self.urhol0 = Function(M)
-
-        hybridized_prb = LinearVariationalProblem(aeqn, Leqn, self.urhol0,
-                                                  constant_jacobian=True)
-        hybridized_solver = LinearVariationalSolver(hybridized_prb,
-                                                    solver_parameters=self.solver_parameters,
-                                                    options_prefix='ImplicitSolver')
-        self.hybridized_solver = hybridized_solver
-
-        # Project broken u into the HDiv space using facet averaging.
-        # Weight function counting the dofs of the HDiv element:
-        self._weight = Function(Vu)
-        weight_kernel = AverageWeightings(Vu)
-        weight_kernel.apply(self._weight)
-
-        # Averaging kernel
-        self._average_kernel = AverageKernel(Vu)
-
-        # HDiv-conforming velocity
-        self.u_hdiv = Function(Vu)
+        # Set up the problem ===================================================
+        if self.formulation == 'hybridized':
+            hybridized_prb = LinearVariationalProblem(
+                aeqn, Leqn, self.urhol0, constant_jacobian=True
+            )
+            hybridized_solver = LinearVariationalSolver(
+                hybridized_prb, solver_parameters=self.solver_parameters,
+                options_prefix='HybridImplicitSolver'
+            )
+            self.hybridized_solver = hybridized_solver
+
+            # Project broken u into the HDiv space using facet averaging.
+            # Weight function counting the dofs of the HDiv element:
+            self._weight = Function(Vu)
+            weight_kernel = AverageWeightings(Vu)
+            weight_kernel.apply(self._weight)
+
+            # Averaging kernel
+            self._average_kernel = AverageKernel(Vu)
+
+            # HDiv-conforming velocity
+            self.u_hdiv = Function(Vu)
+            u_hdiv = self.u_hdiv
+
+            # Store boundary conditions for the div-conforming velocity to apply
+            # post-solve
+            self.bcs = self.equations.bcs['u']
 
-        # Reconstruction of theta
+        else:
+            # Boundary conditions (assumes extruded mesh)
+            self.bcs = [
+                DirichletBC(M.sub(0), 0.0, "bottom"),
+                DirichletBC(M.sub(0), 0.0, "top")
+            ]
+
+            # Solver for u, rho
+            urho_problem = LinearVariationalProblem(
+                aeqn, Leqn, self.urho, bcs=self.bcs
+            )
+            self.urho_solver = LinearVariationalSolver(
+                urho_problem, solver_parameters=self.solver_parameters,
+                options_prefix='ImplicitSolver'
+            )
+            # Velocity to appear in theta reconstruction
+            u_hdiv = self.urho.subfunctions[0]
+
+        # Reconstruction of theta ==============================================
         theta = TrialFunction(Vtheta)
         gamma = TestFunction(Vtheta)
 
         self.theta = Function(Vtheta)
         theta_eqn = gamma*(theta - theta_in
-                           + dot(k, self.u_hdiv)*dot(k, grad(thetabar))*beta_t)*dx
-
-        theta_problem = LinearVariationalProblem(lhs(theta_eqn), rhs(theta_eqn), self.theta,
-                                                 constant_jacobian=True)
-        self.theta_solver = LinearVariationalSolver(theta_problem,
-                                                    solver_parameters=cg_ilu_parameters,
-                                                    options_prefix='thetabacksubstitution')
+                           + dot(k, u_hdiv)*dot(k, grad(thetabar))*beta_t)*dx
 
-        # Store boundary conditions for the div-conforming velocity to apply
-        # post-solve
-        self.bcs = self.equations.bcs['u']
+        theta_problem = LinearVariationalProblem(
+            lhs(theta_eqn), rhs(theta_eqn), self.theta, constant_jacobian=True
+        )
+        self.theta_solver = LinearVariationalSolver(
+            theta_problem, solver_parameters=cg_ilu_parameters,
+            options_prefix='thetabacksubstitution'
+        )
 
-        # Log residuals on hybridized solver
-        self.log_ksp_residuals(self.hybridized_solver.snes.ksp)
-        # Log residuals on the trace system too
-        python_context = self.hybridized_solver.snes.ksp.pc.getPythonContext()
-        attach_custom_monitor(python_context, logging_ksp_monitor_true_residual)
+        if self.formulation == 'hybridized':
+            # Log residuals on hybridized solver
+            self.log_ksp_residuals(self.hybridized_solver.snes.ksp)
+            # Log residuals on the trace system too
+            python_context = self.hybridized_solver.snes.ksp.pc.getPythonContext()
+            attach_custom_monitor(python_context, logging_ksp_monitor_true_residual)
+        elif self.formulation == 'full':
+            # Log residuals on mixed solver
+            self.log_ksp_residuals(self.urho_solver.snes.ksp)
 
     @timed_function("Gusto:UpdateReferenceProfiles")
     def update_reference_profiles(self):
-        with timed_region("Gusto:HybridProjectRhobar"):
-            logger.info('Compressible linear solver: rho average solve')
-            self.rho_avg_solver.solve()
 
-        with timed_region("Gusto:HybridProjectExnerbar"):
-            logger.info('Compressible linear solver: Exner average solve')
-            self.exner_avg_solver.solve()
+        if self.formulation == 'hydridized':
+            with timed_region("Gusto:HybridProjectRhobar"):
+                logger.info('Compressible linear solver: rho average solve')
+                self.rho_avg_solver.solve()
+
+            with timed_region("Gusto:HybridProjectExnerbar"):
+                logger.info('Compressible linear solver: Exner average solve')
+                self.exner_avg_solver.solve()
 
-        # Because the left hand side of the hybridised problem depends
-        # on the reference profile, the Jacobian matrix should change
-        # when the reference profiles are updated. This call will tell
-        # the hybridized_solver to reassemble the Jacobian next time
-        # `solve` is called.
-        self.hybridized_solver.invalidate_jacobian()
+            # Because the left hand side of the hybridised problem depends
+            # on the reference profile, the Jacobian matrix should change
+            # when the reference profiles are updated. This call will tell
+            # the hybridized_solver to reassemble the Jacobian next time
+            # `solve` is called.
+            self.hybridized_solver.invalidate_jacobian()
+
+        elif self.formulation == 'full':
+            self.urho_solver.invalidate_jacobian()
 
     @timed_function("Gusto:LinearSolve")
     def solve(self, xrhs, dy):
@@ -405,23 +515,30 @@ def solve(self, xrhs, dy):
         """
         self.xrhs.assign(xrhs)
 
-        # Solve the hybridized system
-        logger.info('Compressible linear solver: hybridized solve')
-        self.hybridized_solver.solve()
+        # MIXED SOLVE ==========================================================
+        if self.formulation == 'hybridized':
+            # Solve the hybridized system
+            logger.info('Compressible linear solver: hybridized solve')
+            self.hybridized_solver.solve()
+
+            broken_u, rho1, _ = self.urhol0.subfunctions
+            u1 = self.u_hdiv
 
-        broken_u, rho1, _ = self.urhol0.subfunctions
-        u1 = self.u_hdiv
+            # Project broken_u into the HDiv space
+            u1.assign(0.0)
 
-        # Project broken_u into the HDiv space
-        u1.assign(0.0)
+            with timed_region("Gusto:HybridProjectHDiv"):
+                logger.info('Compressible linear solver: restore continuity')
+                self._average_kernel.apply(u1, self._weight, broken_u)
 
-        with timed_region("Gusto:HybridProjectHDiv"):
-            logger.info('Compressible linear solver: restore continuity')
-            self._average_kernel.apply(u1, self._weight, broken_u)
+            # Reapply bcs to ensure they are satisfied
+            for bc in self.bcs:
+                bc.apply(u1)
 
-        # Reapply bcs to ensure they are satisfied
-        for bc in self.bcs:
-            bc.apply(u1)
+        elif self.formulation == 'full':
+            logger.info('Compressible linear solver: mixed solve')
+            self.urho_solver.solve()
+            u1, rho1 = self.urho.subfunctions
 
         # Copy back into u and rho cpts of dy
         u, rho, theta = dy.subfunctions[0:3]
diff --git a/gusto/spatial_methods/spatial_methods.py b/gusto/spatial_methods/spatial_methods.py
index 6ad214dcd..bc0c9a1e1 100644
--- a/gusto/spatial_methods/spatial_methods.py
+++ b/gusto/spatial_methods/spatial_methods.py
@@ -3,9 +3,11 @@
 spatial discretisation of some term.
 """
 
-from firedrake import split
+from firedrake import split, dot
 from firedrake.fml import Term, keep, drop
-from gusto.core.labels import prognostic
+from gusto.core.labels import (
+    prognostic, horizontal_prognostic, vertical_prognostic
+)
 
 __all__ = ['SpatialMethod']
 
@@ -47,6 +49,15 @@ def __init__(self, equation, variable, term_label):
         assert num_terms == 1, f'Unable to find {term_label.label} term ' \
             + f'for {variable}. {num_terms} found'
 
+        # If specified, replace field with only horizontal or vertical part
+        if self.original_form.terms[0].has_label(horizontal_prognostic):
+            k = self.equation.domain.k
+            self.field = self.field - k*dot(k, self.field)
+
+        if self.original_form.terms[0].has_label(vertical_prognostic):
+            k = self.equation.domain.k
+            self.field = k*dot(k, self.field)
+
     def replace_form(self, equation):
         """
         Replaces the form for the term in the equation with the form for the
diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py
index 7b5eee9a6..95e2e03ab 100644
--- a/gusto/time_discretisation/time_discretisation.py
+++ b/gusto/time_discretisation/time_discretisation.py
@@ -423,7 +423,7 @@ def update_subcycling(self):
                         'Adaptive subcycling: capping number of subcycles at '
                         f'{max_subcycles}'
                     )
-                    self.ncycles = max_subcycles
+                self.ncycles = min(self.ncycles, max_subcycles)
 
             logger.debug(f'Performing {self.ncycles} subcycles')
             self.dt.assign(self.original_dt/self.ncycles)
diff --git a/gusto/timestepping/semi_implicit_quasi_newton.py b/gusto/timestepping/semi_implicit_quasi_newton.py
index c0ed20fb1..d172fc6fa 100644
--- a/gusto/timestepping/semi_implicit_quasi_newton.py
+++ b/gusto/timestepping/semi_implicit_quasi_newton.py
@@ -5,15 +5,17 @@
 
 from firedrake import (
     Function, Constant, TrialFunctions, DirichletBC, div, assemble,
-    LinearVariationalProblem, LinearVariationalSolver, FunctionSpace
+    LinearVariationalProblem, LinearVariationalSolver, FunctionSpace,
+    interpolate, dot
 )
 from firedrake.fml import drop, replace_subject
-from firedrake.__future__ import interpolate
+import numpy as np
 from pyop2.profiling import timed_stage
 from gusto.core import TimeLevelFields, StateFields
-from gusto.core.labels import (transport, diffusion, time_derivative,
-                               hydrostatic, physics_label, sponge,
-                               incompressible)
+from gusto.core.labels import (
+    transport, diffusion, time_derivative, hydrostatic, physics_label, sponge,
+    incompressible, implicit
+)
 from gusto.solvers import LinearTimesteppingSolver, mass_parameters
 from gusto.core.logging import logger, DEBUG, logging_ksp_monitor_true_residual
 from gusto.time_discretisation.time_discretisation import ExplicitTimeDiscretisation
@@ -398,6 +400,12 @@ def finish_spinup(self):
         # This only needs doing once, so update the flag
         self.spinup_done = True
 
+    def log_w(self, field, time_stage):# TODO TB: remove this function
+        k = self.equation.domain.k
+        w = Function(self.equation.domain.spaces('theta'))
+        w.interpolate(dot(k, field))
+        logger.info(f'{time_stage} -- w: min {np.min(w.dat.data)}, max {np.max(w.dat.data)}')
+
     def timestep(self):
         """Defines the timestep"""
         xn = self.x.n
@@ -433,6 +441,7 @@ def timestep(self):
             logger.info('Semi-implicit Quasi Newton: Explicit forcing')
             # Put explicit forcing into xstar
             self.forcing.apply(x_after_slow, xn, xstar(self.field_name), "explicit")
+            self.log_w(xstar('u'), 'Explicit forcing')
 
         # set xp here so that variables that are not transported have
         # the correct values
@@ -446,6 +455,7 @@ def timestep(self):
                 self.io.log_courant(self.fields, 'transporting_velocity',
                                     message=f'transporting velocity, outer iteration {outer}')
                 self.transport_fields(outer, xstar, xp)
+                self.log_w(xp('u'), 'Transport')
 
             # Fast physics -----------------------------------------------------
             x_after_fast(self.field_name).assign(xp(self.field_name))
@@ -467,6 +477,7 @@ def timestep(self):
                     if (inner > 0 and self.accelerator):
                         # Zero implicit forcing to accelerate solver convergence
                         self.forcing.zero_forcing_terms(self.equation, xp, xrhs, self.transported_fields)
+                    self.log_w(xrhs.subfunctions[0], 'Implicit forcing')
 
                 xrhs -= xnp1(self.field_name)
                 xrhs += xrhs_phys
@@ -478,6 +489,7 @@ def timestep(self):
 
                 xnp1X = xnp1(self.field_name)
                 xnp1X += dy
+                self.log_w(xnp1('u'), 'After solver')
 
             # Update xnp1 values for active tracers not included in the linear solve
             self.copy_active_tracers(x_after_fast, xnp1)
@@ -547,7 +559,7 @@ def __init__(self, equation, alpha):
         """
 
         self.field_name = equation.field_name
-        implicit_terms = [incompressible, sponge]
+        implicit_terms = [incompressible, sponge, implicit]
         dt = equation.domain.dt
 
         W = equation.function_space