diff --git a/docs/notebook/shallow_water_adjoint.ipynb b/docs/notebook/shallow_water_adjoint.ipynb
new file mode 100644
index 000000000..d0a8c0062
--- /dev/null
+++ b/docs/notebook/shallow_water_adjoint.ipynb
@@ -0,0 +1,768 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Sensitivities computation via automated adjoint for a 2D shallow water model\n",
+    "\n",
+    "This notebook demonstrates how to compute the sensitivities (gradient) of a cost function with respect to control parameters using automated adjoints, a tool provided by the Firedrake library. The model used in this notebook is a shallow water model on a sphere. \n",
+    "\n",
+    "## Cost function and model\n",
+    "The cost function measures the summation of the kinetic and the potential energy of the fluid, which is modelled by the following expression:\n",
+    "\\begin{equation}\n",
+    "J = \\int_{\\Omega} \\left( \\frac{1}{2} \\textbf{u}(\\mathbf{x}, t_f) \\cdot \\textbf{u}(\\mathbf{x}, t_f) + \\frac{1}{2} g D^2(\\mathbf{x}, t_f) \\right) \\, dx,\n",
+    "\\end{equation}\n",
+    "where $\\textbf{u} (\\mathbf{x}, t_f)$ is the velocity of the fluid at the final step $t_f$, $D(\\mathbf{x}, t_f)$ is the fluid depth at the final step, and $g$ is the gravitational acceleration. The model is governed by the following set of equations:\n",
+    "\n",
+    "\\begin{align}\n",
+    "    \\textbf{u}_t + (\\textbf{u} \\cdot \\nabla) \\textbf{u} + f \\textbf{u}^{\\perp} + g \\nabla (D + b)  &= 0, \\tag{2}\\\\\n",
+    "    D_t + \\nabla \\cdot (D \\textbf{u}) &= 0, \\tag{3}\n",
+    "\\end{align}\n",
+    "where $f$ is the Coriolis parameter, the fluid depth is given by $D = H+h-b$, $H$ is the mean fluid depth, $h$ is the free surface height and $b$ is the topography. No boundary conditions are required as we are solving on a spherical domain, $\\Omega$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Solve the model with Gusto\n",
+    "We first import the libraries required for the computations: ``gusto`` and ``firedrake``."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 80,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from firedrake import *\n",
+    "from gusto import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We then define the shallow water parameter $H$ (The fluid depth)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 81,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H = 5960.0\n",
+    "parameters = ShallowWaterParameters(H=H)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We will use one of the spherical meshes provided by Firedrake: the ``IcosahedralSphereMesh``. As the spherical domain we are solving over is the Earth we specify the radius as 6371220m. The refinement level, ``ref_level``, specifies the number of times the base icosahedron is refined. The argument ``degree`` specifies the polynomial degree of the function space used to represent the coordinates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 82,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "R = 6371220. # Earth radius in meters\n",
+    "mesh = IcosahedralSphereMesh(radius=R, refinement_level=3, degree=2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A domain is created with the ``Domain`` Gusto object, which holds the mesh and the function spaces defined on it. It also holds the model’s time interval."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 83,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dt = 900.\n",
+    "domain = Domain(mesh, dt, 'BDM', 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now set up the finite element form of the shallow water equations by passing ``ShallowWaterEquations`` Gusto class the following arguments: ``domain``, ``parameters``,  the expressions ``fexpr`` for the Coriolis parameter and ``bexpr`` for the bottom surface of the fluid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 84,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = SpatialCoordinate(mesh)\n",
+    "Omega = parameters.Omega  # Spatial domain\n",
+    "fexpr = 2*Omega*x[2]/R  # Expression for the Coriolis parameter\n",
+    "# ``lonlatr_from_xyz`` is returning the required spherical longitude ``lamda``, latitude ``theta``.\n",
+    "lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2])\n",
+    "R0 = pi/9. # radius of mountain (rad)\n",
+    "R0sq = R0**2  \n",
+    "lamda_c = -pi/2. # longitudinal centre of mountain (rad)\n",
+    "lsq = (lamda - lamda_c)**2 \n",
+    "theta_c = pi/6.  # latitudinal centre of mountain (rad)\n",
+    "thsq = (theta - theta_c)**2\n",
+    "rsq = min_value(R0sq, lsq+thsq)\n",
+    "r = sqrt(rsq)\n",
+    "bexpr = 2000 * (1 - r/R0)  # An expression for the bottom surface of the fluid.\n",
+    "eqn = ShallowWaterEquations(domain, parameters, fexpr=fexpr, bexpr=bexpr)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We specify instructions regarding the output using the ``OutputParameters`` class. The directory name, ``dirname``, must be specified. To prevent losing hard-earned simulation data, Gusto will not allow existing files to be overwritten. Hence, if one wishes to re-run a simulation with the same output filename, the existing results file needs to be moved or deleted first. This is also the place to specify the output frequency (in timesteps) of VTU files. The default is ``dumpfreq=1``, which outputs VTU files at every timestep (very useful when first setting up a problem!). Below we set ``dumpfreq=5``.\n",
+    "\n",
+    "We can specify which diagnostics to record over a simulation. The list of available diagnostics in the gusto source code: [diagnostics](https://github.com/firedrakeproject/gusto/blob/main/gusto/diagnostics/diagnostics.py). Since this flow should be in a steady state, it is also instructive to output the steady state error fields for both $\\textbf{u}$ and $D$ as this will tell us how close the simulation is to be correct. The errors should not grow in time. They should reduce as the mesh and timestep are refined. We pass these diagnostics into the ``IO`` class, which controls the input and output and stores the fields which will be updated at each timestep."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 85,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "output = OutputParameters(dirname=\"adjoint_sw\", log_courant=False)\n",
+    "io = IO(domain, output)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now start to build a setup for time discretisation. We use the ``SemiImplicitQuasiNewton`` approach, which splits the equation into `transport` terms and 'forcing' terms (i.e. everything that does not transport) and solves each separately. That allows for different time-steppers used to transport the velocity and depth fields. We are employing an Implicit Midpoint method for the velocity and an explicit strong stability preserving RK3 (SSPRK3) method for the depth. Since the Courant number for a stable SSPRK3 scheme is lower than that for the Implicit Midpoint method, we do two subcycles of the SSPRK3 scheme per timestep, allowing us to use a longer timestep overall. The full list of time-stepping methods is available at the [time discretisation](https://github.com/firedrakeproject/gusto/blob/main/gusto/time_discretisation.py) python file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 86,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO     Physical parameters that take non-default values:\n",
+      "INFO:gusto:Physical parameters that take non-default values:\n",
+      "INFO     H: 5960.0\n",
+      "INFO:gusto:H: 5960.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Transport schemes\n",
+    "transported_fields = [TrapeziumRule(domain, \"u\"), SSPRK3(domain, \"D\")]\n",
+    "transport_methods = [DGUpwind(eqn, \"u\"), DGUpwind(eqn, \"D\")]\n",
+    "\n",
+    "# Time stepper\n",
+    "stepper = SemiImplicitQuasiNewton(eqn, io, transported_fields, transport_methods)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We are now ready to specify the initial conditions:\n",
+    "\\begin{align}\n",
+    "    \\textbf{u}_0 &= \\frac{u_{max}}{R} [-y,x,0], \\tag{4}\\\\\n",
+    "    D_0 &= H - \\frac{\\Omega u_{max} z^2}{g R} \\tag{5}\n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 87,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Coefficient(WithGeometry(IndexedProxyFunctionSpace(<firedrake.mesh.MeshTopology object at 0x1682cee90>, FiniteElement('Discontinuous Lagrange', triangle, 1), name='L2', index=1, component=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 2), dim=3), 175601)), 248029)"
+      ]
+     },
+     "execution_count": 87,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "u0 = stepper.fields('u')\n",
+    "D0 = stepper.fields('D')\n",
+    "u_max = 20.   # Maximum amplitude of the zonal wind (m/s)\n",
+    "uexpr = as_vector([-u_max*x[1]/R, u_max*x[0]/R, 0.0])\n",
+    "g = parameters.g  # acceleration due to gravity (m/s^2)\n",
+    "Rsq = R**2\n",
+    "Dexpr = H - ((R * Omega * u_max + 0.5*u_max**2)*x[2]**2/Rsq)/g - bexpr\n",
+    "\n",
+    "u0.project(uexpr)\n",
+    "D0.interpolate(Dexpr)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When using the `SemiImplicitQuasiNewton` time ``stepper``, we also have to set up any non-zero reference profiles."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 88,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Dbar = Function(D0.function_space()).assign(H)\n",
+    "stepper.set_reference_profiles([('D', Dbar)])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We are almost ready to run the model. Before we do so, we have to tape the model to compute the sensitivities via automated adjoint. Hence, at this point, we import ``firedrake.adjoint`` and enable the tape of the forward solver (Shallow Water model) with ``continue_annotation()``."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 89,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from firedrake.adjoint import *\n",
+    "continue_annotation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We are finally ready to run the model!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 90,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "DEBUG:gusto:Dumping output to disk\n",
+      "INFO     \n",
+      "INFO:gusto:\n",
+      "INFO     ========================================\n",
+      "INFO:gusto:========================================\n",
+      "INFO     at start of timestep 1, t=0.0, dt=900.0\n",
+      "INFO:gusto:at start of timestep 1, t=0.0, dt=900.0\n",
+      "INFO     Semi-implicit Quasi Newton: Explicit forcing\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Explicit forcing\n",
+      "DEBUG:gusto:    Residual norms for ExplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 3.717272816438e+05 true resid norm 3.717272816438e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.620107299083e-10 true resid norm 1.620107299083e-10 ||r(i)||/||b|| 4.358322294556e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.585397212708e+05 true resid norm 6.127594999762e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.064757082310e+04 true resid norm 6.297866206201e+03 ||r(i)||/||b|| 1.027787607772e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 1.984167553564e+02 true resid norm 2.154172949430e+02 ||r(i)||/||b|| 3.515527624645e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 4.896138420464e+00 true resid norm 5.372370990865e+00 ||r(i)||/||b|| 8.767503386031e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 4.896138420921e+00 true resid norm 5.372370992372e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.213180858366e-01 true resid norm 1.505771328497e-01 ||r(i)||/||b|| 2.802805931748e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.386110502201e-03 true resid norm 7.233829275614e-03 ||r(i)||/||b|| 1.346487293205e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.888081698314e-04 true resid norm 2.155980867138e-04 ||r(i)||/||b|| 4.013090068052e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.508958609457e-06 true resid norm 8.412262624625e-06 ||r(i)||/||b|| 1.565837995285e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.901594446440e-07 true resid norm 5.901594446440e-07 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.718314622419e-22 true resid norm 3.718314622419e-22 ||r(i)||/||b|| 6.300525487076e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 3.478820869067e+05 true resid norm 3.478820869067e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.994472980877e-10 true resid norm 1.994472980877e-10 ||r(i)||/||b|| 5.733186777771e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.587569440754e+05 true resid norm 6.129105611037e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.065740127148e+04 true resid norm 6.297890541694e+03 ||r(i)||/||b|| 1.027538264368e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 1.985356265347e+02 true resid norm 2.156115172251e+02 ||r(i)||/||b|| 3.517830021347e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 4.881602606427e+00 true resid norm 5.364127741803e+00 ||r(i)||/||b|| 8.751893150843e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 4.881602606361e+00 true resid norm 5.364127745118e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.189409731723e-01 true resid norm 1.489319454693e-01 ||r(i)||/||b|| 2.776442928767e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.287887155869e-03 true resid norm 7.132861396961e-03 ||r(i)||/||b|| 1.329733693135e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.876689010863e-04 true resid norm 2.147829601784e-04 ||r(i)||/||b|| 4.004061245072e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.484304727953e-06 true resid norm 8.385391455195e-06 ||r(i)||/||b|| 1.563234854507e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.375704065564e+05 true resid norm 3.375200844565e+05 ||r(i)||/||b|| 2.453435247486e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 5.515968938792e-11 true resid norm 5.515968938792e-11 ||r(i)||/||b|| 4.009560687406e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.066019019371e+05 true resid norm 2.423054418141e+05 ||r(i)||/||b|| 2.272993609037e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.961061610235e-11 true resid norm 3.961061610235e-11 ||r(i)||/||b|| 3.715751349890e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "DEBUG:gusto:Leaving Semi-implicit Quasi-Newton timestep method\n",
+      "INFO     \n",
+      "INFO:gusto:\n",
+      "INFO     ========================================\n",
+      "INFO:gusto:========================================\n",
+      "INFO     at start of timestep 2, t=900.0, dt=900.0\n",
+      "INFO:gusto:at start of timestep 2, t=900.0, dt=900.0\n",
+      "INFO     Semi-implicit Quasi Newton: Explicit forcing\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Explicit forcing\n",
+      "DEBUG:gusto:    Residual norms for ExplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.362120374912e+04 true resid norm 3.859973858259e+05 ||r(i)||/||b|| 7.198596055991e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.893877795969e-11 true resid norm 1.893877795969e-11 ||r(i)||/||b|| 3.531956881890e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.604766785161e+05 true resid norm 6.152090559168e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.068399186592e+04 true resid norm 6.329441734408e+03 ||r(i)||/||b|| 1.028827790088e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.000708822109e+02 true resid norm 2.162420601441e+02 ||r(i)||/||b|| 3.514936232886e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 4.931342828052e+00 true resid norm 5.396506875620e+00 ||r(i)||/||b|| 8.771826135717e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 4.931342827824e+00 true resid norm 5.396506873935e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.227791166505e-01 true resid norm 1.521523930563e-01 ||r(i)||/||b|| 2.819460747677e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.471748032847e-03 true resid norm 7.301465459939e-03 ||r(i)||/||b|| 1.352998454464e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.898529502076e-04 true resid norm 2.167792649846e-04 ||r(i)||/||b|| 4.017029349701e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.555933266040e-06 true resid norm 8.456006918931e-06 ||r(i)||/||b|| 1.566940822363e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.922290117006e-07 true resid norm 1.066019019371e+05 ||r(i)||/||b|| 1.800011479190e+11\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.792661951956e-22 true resid norm 3.792661951956e-22 ||r(i)||/||b|| 6.404046200075e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 4.082629514187e+05 true resid norm 4.082629514187e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.613042924564e-10 true resid norm 1.613042924564e-10 ||r(i)||/||b|| 3.950990210987e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.613736253633e+05 true resid norm 6.164185032864e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.069411765753e+04 true resid norm 6.344717297526e+03 ||r(i)||/||b|| 1.029287288376e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.006967024553e+02 true resid norm 2.165291243542e+02 ||r(i)||/||b|| 3.512696701994e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 4.955181626652e+00 true resid norm 5.408670055491e+00 ||r(i)||/||b|| 8.774347341385e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 4.955181626602e+00 true resid norm 5.408670055861e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.250540665393e-01 true resid norm 1.538410757549e-01 ||r(i)||/||b|| 2.844342031702e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.571771905522e-03 true resid norm 7.390528689363e-03 ||r(i)||/||b|| 1.366422542517e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.910313160878e-04 true resid norm 2.176203600240e-04 ||r(i)||/||b|| 4.023546597897e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.593703303884e-06 true resid norm 8.487476707684e-06 ||r(i)||/||b|| 1.569235434964e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.547025377642e+05 true resid norm 4.944980624256e+05 ||r(i)||/||b|| 3.196444412432e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 6.244438529493e-11 true resid norm 6.244438529493e-11 ||r(i)||/||b|| 4.036416350848e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.129613821142e+05 true resid norm 2.644949396295e+05 ||r(i)||/||b|| 2.341463380486e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 5.188974790380e-11 true resid norm 5.188974790380e-11 ||r(i)||/||b|| 4.593582951327e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "DEBUG:gusto:Leaving Semi-implicit Quasi-Newton timestep method\n",
+      "INFO     \n",
+      "INFO:gusto:\n",
+      "INFO     ========================================\n",
+      "INFO:gusto:========================================\n",
+      "INFO     at start of timestep 3, t=1800.0, dt=900.0\n",
+      "INFO:gusto:at start of timestep 3, t=1800.0, dt=900.0\n",
+      "INFO     Semi-implicit Quasi Newton: Explicit forcing\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Explicit forcing\n",
+      "DEBUG:gusto:    Residual norms for ExplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 4.918470774046e+04 true resid norm 4.603785515452e+04 ||r(i)||/||b|| 9.360196953382e-01\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.709363306585e-11 true resid norm 1.709363306585e-11 ||r(i)||/||b|| 3.475395880372e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.628475674062e+05 true resid norm 6.207601980794e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.066580831543e+04 true resid norm 6.385116455001e+03 ||r(i)||/||b|| 1.028596304137e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.011618065183e+02 true resid norm 2.163257559986e+02 ||r(i)||/||b|| 3.484852229056e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.009107162475e+00 true resid norm 5.422766781313e+00 ||r(i)||/||b|| 8.735686982011e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.009107162336e+00 true resid norm 5.422766781943e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.310238285189e-01 true resid norm 1.585661985134e-01 ||r(i)||/||b|| 2.924082943073e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.804915566986e-03 true resid norm 7.593813246366e-03 ||r(i)||/||b|| 1.400357705157e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.933081298615e-04 true resid norm 2.185050909137e-04 ||r(i)||/||b|| 4.029402327263e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.653093695177e-06 true resid norm 8.546658817820e-06 ||r(i)||/||b|| 1.576069774249e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.845374237257e-07 true resid norm 1.129613821142e+05 ||r(i)||/||b|| 1.932491873560e+11\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.760932244576e-22 true resid norm 3.760932244576e-22 ||r(i)||/||b|| 6.434031581082e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 4.068447998955e+05 true resid norm 4.068447998955e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.302395903596e-10 true resid norm 1.302395903596e-10 ||r(i)||/||b|| 3.201210643298e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.637275978130e+05 true resid norm 6.223345449248e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.068196732817e+04 true resid norm 6.400945806866e+03 ||r(i)||/||b|| 1.028537763019e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.017322034986e+02 true resid norm 2.166513722980e+02 ||r(i)||/||b|| 3.481268620949e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.031640760550e+00 true resid norm 5.431779556009e+00 ||r(i)||/||b|| 8.728070135757e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.031640760269e+00 true resid norm 5.431779557485e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.327606111123e-01 true resid norm 1.598490164410e-01 ||r(i)||/||b|| 2.942848006795e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.868316502526e-03 true resid norm 7.641905040180e-03 ||r(i)||/||b|| 1.406887919384e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.945661057793e-04 true resid norm 2.190233858132e-04 ||r(i)||/||b|| 4.032258369384e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.688507454708e-06 true resid norm 8.579874687314e-06 ||r(i)||/||b|| 1.579569751775e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.199022954527e+05 true resid norm 4.055284066275e+05 ||r(i)||/||b|| 3.382157156344e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 5.290535120365e-11 true resid norm 5.290535120365e-11 ||r(i)||/||b|| 4.412371840247e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 9.656092509913e+04 true resid norm 2.151199852505e+05 ||r(i)||/||b|| 2.227816117437e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 4.368178127942e-11 true resid norm 4.368178127942e-11 ||r(i)||/||b|| 4.523753395545e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "DEBUG:gusto:Leaving Semi-implicit Quasi-Newton timestep method\n",
+      "INFO     \n",
+      "INFO:gusto:\n",
+      "INFO     ========================================\n",
+      "INFO:gusto:========================================\n",
+      "INFO     at start of timestep 4, t=2700.0, dt=900.0\n",
+      "INFO:gusto:at start of timestep 4, t=2700.0, dt=900.0\n",
+      "INFO     Semi-implicit Quasi Newton: Explicit forcing\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Explicit forcing\n",
+      "DEBUG:gusto:    Residual norms for ExplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 3.623731520079e+04 true resid norm 6.791995697676e+04 ||r(i)||/||b|| 1.874309854370e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.377006062130e-11 true resid norm 1.377006062130e-11 ||r(i)||/||b|| 3.799967118149e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.655756934086e+05 true resid norm 6.263209139889e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.071506957247e+04 true resid norm 6.441284811339e+03 ||r(i)||/||b|| 1.028432017433e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.023667345222e+02 true resid norm 2.177420973958e+02 ||r(i)||/||b|| 3.476526051303e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.037205816478e+00 true resid norm 5.414858593667e+00 ||r(i)||/||b|| 8.645501807023e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.037205816945e+00 true resid norm 5.414858596158e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.308305953126e-01 true resid norm 1.583788999825e-01 ||r(i)||/||b|| 2.924894476374e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.800830603278e-03 true resid norm 7.535497309851e-03 ||r(i)||/||b|| 1.391633258013e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.951814065919e-04 true resid norm 2.176322724719e-04 ||r(i)||/||b|| 4.019168157527e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.684913337251e-06 true resid norm 8.591881993703e-06 ||r(i)||/||b|| 1.586723243299e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.930426159714e-07 true resid norm 9.656092509912e+04 ||r(i)||/||b|| 1.628229110331e+11\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.788205038540e-22 true resid norm 3.788205038540e-22 ||r(i)||/||b|| 6.387745056626e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 2.870374044092e+05 true resid norm 2.870374044092e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.282656149014e-10 true resid norm 1.282656149014e-10 ||r(i)||/||b|| 4.468602799884e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.664936934029e+05 true resid norm 6.275513396511e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.075136327118e+04 true resid norm 6.442870975745e+03 ||r(i)||/||b|| 1.026668348653e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.030257661540e+02 true resid norm 2.186361917364e+02 ||r(i)||/||b|| 3.483957055337e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.034072905720e+00 true resid norm 5.418192805198e+00 ||r(i)||/||b|| 8.633863817756e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.034072905562e+00 true resid norm 5.418192806697e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.298184557140e-01 true resid norm 1.575351913962e-01 ||r(i)||/||b|| 2.907522803572e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.767180921667e-03 true resid norm 7.495866514393e-03 ||r(i)||/||b|| 1.383462490506e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.951345461690e-04 true resid norm 2.173722948979e-04 ||r(i)||/||b|| 4.011896635152e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.682088429556e-06 true resid norm 8.582316804832e-06 ||r(i)||/||b|| 1.583981432743e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.346881528775e+05 true resid norm 2.754043729863e+05 ||r(i)||/||b|| 2.044755734655e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 4.386840848005e-11 true resid norm 4.386840848005e-11 ||r(i)||/||b|| 3.257035421665e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.032755124426e+05 true resid norm 2.367036499603e+05 ||r(i)||/||b|| 2.291962967425e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.628393263694e-11 true resid norm 3.628393263694e-11 ||r(i)||/||b|| 3.513314219293e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "DEBUG:gusto:Leaving Semi-implicit Quasi-Newton timestep method\n",
+      "INFO     \n",
+      "INFO:gusto:\n",
+      "INFO     ========================================\n",
+      "INFO:gusto:========================================\n",
+      "INFO     at start of timestep 5, t=3600.0, dt=900.0\n",
+      "INFO:gusto:at start of timestep 5, t=3600.0, dt=900.0\n",
+      "INFO     Semi-implicit Quasi Newton: Explicit forcing\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Explicit forcing\n",
+      "DEBUG:gusto:    Residual norms for ExplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 2.993856955300e+04 true resid norm 3.914845825554e+04 ||r(i)||/||b|| 1.307626210606e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 1.117982053676e-11 true resid norm 1.117982053676e-11 ||r(i)||/||b|| 3.734253407454e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.712284756826e+05 true resid norm 6.318020945328e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.092236059309e+04 true resid norm 6.475907500247e+03 ||r(i)||/||b|| 1.024989875198e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.068099444904e+02 true resid norm 2.232405514586e+02 ||r(i)||/||b|| 3.533393659033e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.037550826001e+00 true resid norm 5.449947978305e+00 ||r(i)||/||b|| 8.626036579279e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.037550825972e+00 true resid norm 5.449947976511e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.262184219078e-01 true resid norm 1.548699263179e-01 ||r(i)||/||b|| 2.841677149679e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.660094330537e-03 true resid norm 7.385726291473e-03 ||r(i)||/||b|| 1.355192072164e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.943513776045e-04 true resid norm 2.167003730878e-04 ||r(i)||/||b|| 3.976191589750e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.676621445170e-06 true resid norm 8.553356571439e-06 ||r(i)||/||b|| 1.569438205337e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 0: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 5.802708268023e-07 true resid norm 1.032755124426e+05 ||r(i)||/||b|| 1.779781227530e+11\n",
+      "DEBUG:gusto:    1 KSP no resid norm 3.767179856699e-22 true resid norm 3.767179856699e-22 ||r(i)||/||b|| 6.492106241941e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (0, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 2.027107463164e+05 true resid norm 2.027107463164e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 7.869384266651e-11 true resid norm 7.869384266651e-11 ||r(i)||/||b|| 3.882075523697e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (0, 1)\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: u\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: u\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.733314345061e+05 true resid norm 6.333470343240e+05 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 1.098965232923e+04 true resid norm 6.488475328056e+03 ||r(i)||/||b|| 1.024473941838e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 2.084150171367e+02 true resid norm 2.252619676010e+02 ||r(i)||/||b|| 3.556690967084e-04\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 5.035195935951e+00 true resid norm 5.465094804057e+00 ||r(i)||/||b|| 8.628910388585e-06\n",
+      "DEBUG:gusto:    Residual norms for uTrapeziumRule_ solve\n",
+      "DEBUG:gusto:    0 KSP preconditioned resid norm 5.035195935774e+00 true resid norm 5.465094802197e+00 ||r(i)||/||b|| 1.000000000000e+00\n",
+      "DEBUG:gusto:    1 KSP preconditioned resid norm 2.252681932078e-01 true resid norm 1.541017486928e-01 ||r(i)||/||b|| 2.819745206082e-02\n",
+      "DEBUG:gusto:    2 KSP preconditioned resid norm 8.634772911632e-03 true resid norm 7.362233497288e-03 ||r(i)||/||b|| 1.347137380733e-03\n",
+      "DEBUG:gusto:    3 KSP preconditioned resid norm 2.938027804860e-04 true resid norm 2.166120739369e-04 ||r(i)||/||b|| 3.963555652316e-05\n",
+      "DEBUG:gusto:    4 KSP preconditioned resid norm 9.671700592498e-06 true resid norm 8.527053629568e-06 ||r(i)||/||b|| 1.560275519125e-06\n",
+      "INFO     Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Transport 1: D\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 0)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.401792544219e+05 true resid norm 2.581791465204e+05 ||r(i)||/||b|| 1.841778568342e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 5.654885791223e-11 true resid norm 5.654885791223e-11 ||r(i)||/||b|| 4.034038998527e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 0)\n",
+      "INFO     Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Implicit forcing (1, 1)\n",
+      "DEBUG:gusto:    Residual norms for ImplicitForcingSolver_ solve\n",
+      "DEBUG:gusto:    0 KSP no resid norm 1.099362982095e+05 true resid norm 2.494150460533e+05 ||r(i)||/||b|| 2.268723343568e+00\n",
+      "DEBUG:gusto:    1 KSP no resid norm 4.166161929637e-11 true resid norm 4.166161929637e-11 ||r(i)||/||b|| 3.789614529042e-16\n",
+      "INFO     Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "INFO:gusto:Semi-implicit Quasi Newton: Mixed solve (1, 1)\n",
+      "DEBUG:gusto:Leaving Semi-implicit Quasi-Newton timestep method\n",
+      "INFO     TIMELOOP complete. t=4500.00000, tmax=4500.00000\n",
+      "INFO:gusto:TIMELOOP complete. t=4500.00000, tmax=4500.00000\n"
+     ]
+    }
+   ],
+   "source": [
+    "timesteps = 5 \n",
+    "stepper.run(0., timesteps*dt)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, the final solution of the shallow water solver is used to compute the cost functional (1) as shown in the following cell code."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 91,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u_tf = stepper.fields('u')  # Final velocity field\n",
+    "D_tf = stepper.fields('D')  # Final depth field\n",
+    "\n",
+    "J = assemble(0.5*inner(u_tf, u_tf)*dx + 0.5*g*D_tf**2*dx)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We define the control variable: initial fluid depth $D_0$ and the velocity $\\textbf{u}_0$. Finally, we run the sensitivity computation: gradient of the cost function with respect $D_0$ and $\\textbf{u}_0$. The gradient computation is achieved with the ``compute_gradient`` function from the `firedrake.adjoint` module."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 92,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "control = [Control(D0), Control(u0)] # Control variables\n",
+    "grad = compute_gradient(J, control)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can execute a gradient verification with second order [Taylor's test](https://www.dolfinadjoint.org/en/latest/documentation/verification.html). To do so, we define the reduced functional as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 93,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "J_hat = ReducedFunctional(J, control)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We then verify the gradient computation using the ``taylor_test``. We already taped the forward model. Any additional operation as below, for  Taylor's test, does not require the tape to be enabled. Thus, we can stop annotating."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 94,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running Taylor test\n",
+      "Computed residuals: [151981167509504.0, 37995189125120.0, 9498745905152.0, 2374616748032.0]\n",
+      "Computed convergence rates: [np.float64(2.0000039015457785), np.float64(2.0000078031232102), np.float64(2.0000423626811754)]\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "with stop_annotating():\n",
+    "    # Stop annotation to perform the Taylor test\n",
+    "    h0 = Function(D0.function_space())\n",
+    "    h1 = Function(u0.function_space())\n",
+    "    h0.assign(D0 * np.random.rand())\n",
+    "    h1.assign(u0 * np.random.rand())\n",
+    "    taylor_test(J_hat, [D0, u0], [h0, h1])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "That is great! We have a sensitivity computation via automated adjoint with Gusto! In addition, we have verified the gradient computation with a second-order Taylor's test, which is a good practice to ensure the correctness of the gradient computation. See the results above are close to the expected values, 2.0."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For a notebook demonstration of adjoints, it is good practice to clear the tape. This is because if we do not clear the tape, the tape will keep growing every time we run the model."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 95,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tape = get_working_tape()\n",
+    "tape.clear_tape()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "firedrake",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/boussinesq/skamarock_klemp_compressible.py b/examples/boussinesq/skamarock_klemp_compressible.py
index c15daff84..5fb9ce075 100644
--- a/examples/boussinesq/skamarock_klemp_compressible.py
+++ b/examples/boussinesq/skamarock_klemp_compressible.py
@@ -66,7 +66,7 @@ def skamarock_klemp_compressible_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters(cs=cs)
+    parameters = BoussinesqParameters(mesh, cs=cs)
     eqns = BoussinesqEquations(domain, parameters)
 
     # I/O
diff --git a/examples/boussinesq/skamarock_klemp_incompressible.py b/examples/boussinesq/skamarock_klemp_incompressible.py
index b7bb5fd5e..8aa3d7265 100644
--- a/examples/boussinesq/skamarock_klemp_incompressible.py
+++ b/examples/boussinesq/skamarock_klemp_incompressible.py
@@ -65,7 +65,7 @@ def skamarock_klemp_incompressible_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters()
+    parameters = BoussinesqParameters(mesh)
     eqns = BoussinesqEquations(domain, parameters, compressible=False)
 
     # I/O
diff --git a/examples/boussinesq/skamarock_klemp_linear.py b/examples/boussinesq/skamarock_klemp_linear.py
index e8145548d..2d430b58b 100644
--- a/examples/boussinesq/skamarock_klemp_linear.py
+++ b/examples/boussinesq/skamarock_klemp_linear.py
@@ -63,7 +63,7 @@ def skamarock_klemp_linear_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters(cs=cs)
+    parameters = BoussinesqParameters(mesh, cs=cs)
     eqns = LinearBoussinesqEquations(domain, parameters)
 
     # I/O
diff --git a/examples/compressible_euler/dcmip_3_1_gravity_wave.py b/examples/compressible_euler/dcmip_3_1_gravity_wave.py
index 6862dd94f..15cf22c1d 100644
--- a/examples/compressible_euler/dcmip_3_1_gravity_wave.py
+++ b/examples/compressible_euler/dcmip_3_1_gravity_wave.py
@@ -43,16 +43,10 @@ def dcmip_3_1_gravity_wave(
     # Test case parameters
     # ------------------------------------------------------------------------ #
 
-    parameters = CompressibleParameters()
     a_ref = 6.37122e6               # Radius of the Earth (m)
     X = 125.0                       # Reduced-size Earth reduction factor
     a = a_ref/X                     # Scaled radius of planet (m)
-    g = parameters.g                # Acceleration due to gravity (m/s^2)
     N = 0.01                        # Brunt-Vaisala frequency (1/s)
-    p_0 = parameters.p_0            # Reference pressure (Pa, not hPa)
-    c_p = parameters.cp             # SHC of dry air at constant pressure (J/kg/K)
-    R_d = parameters.R_d            # Gas constant for dry air (J/kg/K)
-    kappa = parameters.kappa        # R_d/c_p
     T_eq = 300.0                    # Isothermal atmospheric temperature (K)
     p_eq = 1000.0 * 100.0           # Reference surface pressure at the equator
     u_max = 20.0                    # Maximum amplitude of the zonal wind (m/s)
@@ -81,6 +75,7 @@ def dcmip_3_1_gravity_wave(
         extrusion_type="radial"
     )
     domain = Domain(mesh, dt, "RTCF", element_order)
+    parameters = CompressibleParameters(mesh)
 
     # Equation
     eqns = CompressibleEulerEquations(
@@ -137,6 +132,13 @@ def dcmip_3_1_gravity_wave(
     # Initial conditions with u0
     uexpr = as_vector([-u_max*y/a, u_max*x/a, 0.0])
 
+    # Parameters required for initialising
+    g = parameters.g                # Acceleration due to gravity (m/s^2)
+    p_0 = parameters.p_0            # Reference pressure (Pa, not hPa)
+    c_p = parameters.cp             # SHC of dry air at constant pressure (J/kg/K)
+    R_d = parameters.R_d            # Gas constant for dry air (J/kg/K)
+    kappa = parameters.kappa        # R_d/c_p
+
     # Surface temperature
     G = g**2/(N**2*c_p)
     Ts_expr = (
diff --git a/examples/compressible_euler/dry_bryan_fritsch.py b/examples/compressible_euler/dry_bryan_fritsch.py
index dc4c92e9e..380f182fd 100644
--- a/examples/compressible_euler/dry_bryan_fritsch.py
+++ b/examples/compressible_euler/dry_bryan_fritsch.py
@@ -66,7 +66,7 @@ def dry_bryan_fritsch(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, params, u_transport_option=u_eqn_type,
         no_normal_flow_bc_ids=[1, 2]
diff --git a/examples/compressible_euler/straka_bubble.py b/examples/compressible_euler/straka_bubble.py
index fed55a1ec..d33eef5d5 100644
--- a/examples/compressible_euler/straka_bubble.py
+++ b/examples/compressible_euler/straka_bubble.py
@@ -70,8 +70,8 @@ def straka_bubble(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    parameters = CompressibleParameters()
-    diffusion_params = DiffusionParameters(kappa=kappa, mu=mu0/delta)
+    parameters = CompressibleParameters(mesh)
+    diffusion_params = DiffusionParameters(mesh, kappa=kappa, mu=mu0/delta)
     diffusion_options = [("u", diffusion_params), ("theta", diffusion_params)]
     eqns = CompressibleEulerEquations(
         domain, parameters, u_transport_option=u_eqn_type,
diff --git a/examples/compressible_euler/unsaturated_bubble.py b/examples/compressible_euler/unsaturated_bubble.py
index 59a911cb9..f7af29759 100644
--- a/examples/compressible_euler/unsaturated_bubble.py
+++ b/examples/compressible_euler/unsaturated_bubble.py
@@ -84,7 +84,7 @@ def unsaturated_bubble(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     tracers = [WaterVapour(), CloudWater(), Rain()]
     eqns = CompressibleEulerEquations(
         domain, params, active_tracers=tracers, u_transport_option=u_eqn_type
diff --git a/examples/shallow_water/linear_thermal_galewsky_jet.py b/examples/shallow_water/linear_thermal_galewsky_jet.py
index a4fd68631..a5622c055 100644
--- a/examples/shallow_water/linear_thermal_galewsky_jet.py
+++ b/examples/shallow_water/linear_thermal_galewsky_jet.py
@@ -64,7 +64,7 @@ def linear_thermal_galewsky_jet(
     domain = Domain(mesh, dt, 'BDM', element_order)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H)
+    parameters = ShallowWaterParameters(mesh, H=H)
     Omega = parameters.Omega
     fexpr = 2*Omega*xyz[2]/R
     eqns = LinearThermalShallowWaterEquations(domain, parameters, fexpr=fexpr)
diff --git a/examples/shallow_water/linear_williamson_2.py b/examples/shallow_water/linear_williamson_2.py
index 9d022bf11..a50b68697 100644
--- a/examples/shallow_water/linear_williamson_2.py
+++ b/examples/shallow_water/linear_williamson_2.py
@@ -57,7 +57,7 @@ def linear_williamson_2(
     domain = Domain(mesh, dt, 'BDM', element_order)
 
     # Equation
-    parameters = ShallowWaterParameters(H=mean_depth)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
     eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr)
diff --git a/examples/shallow_water/moist_convective_williamson_2.py b/examples/shallow_water/moist_convective_williamson_2.py
index 07fc7e96c..12c585128 100644
--- a/examples/shallow_water/moist_convective_williamson_2.py
+++ b/examples/shallow_water/moist_convective_williamson_2.py
@@ -76,7 +76,7 @@ def moist_convect_williamson_2(
     _, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equations
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 
diff --git a/examples/shallow_water/moist_thermal_equivb_gw.py b/examples/shallow_water/moist_thermal_equivb_gw.py
index 912e9deba..f2928461c 100644
--- a/examples/shallow_water/moist_thermal_equivb_gw.py
+++ b/examples/shallow_water/moist_thermal_equivb_gw.py
@@ -59,7 +59,8 @@ def moist_thermal_gw(
     xyz = SpatialCoordinate(mesh)
 
     # Equation parameters
-    parameters = ShallowWaterParameters(H=mean_depth, q0=q0, nu=nu, beta2=beta2)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, q0=q0,
+                                        nu=nu, beta2=beta2)
     Omega = parameters.Omega
     fexpr = 2*Omega*xyz[2]/radius
 
diff --git a/examples/shallow_water/moist_thermal_williamson_5.py b/examples/shallow_water/moist_thermal_williamson_5.py
index f17c278c0..044719a7e 100644
--- a/examples/shallow_water/moist_thermal_williamson_5.py
+++ b/examples/shallow_water/moist_thermal_williamson_5.py
@@ -86,7 +86,7 @@ def moist_thermal_williamson_5(
     lamda, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equation: coriolis
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 
diff --git a/examples/shallow_water/shallow_water_1d_wave.py b/examples/shallow_water/shallow_water_1d_wave.py
index 3e0b1582b..057a1f157 100644
--- a/examples/shallow_water/shallow_water_1d_wave.py
+++ b/examples/shallow_water/shallow_water_1d_wave.py
@@ -62,9 +62,9 @@ def shallow_water_1d_wave(
 
     # Diffusion
     delta = domain_length / ncells
-    u_diffusion_opts = DiffusionParameters(kappa=kappa)
-    v_diffusion_opts = DiffusionParameters(kappa=kappa, mu=10/delta)
-    D_diffusion_opts = DiffusionParameters(kappa=kappa, mu=10/delta)
+    u_diffusion_opts = DiffusionParameters(mesh, kappa=kappa)
+    v_diffusion_opts = DiffusionParameters(mesh, kappa=kappa, mu=10/delta)
+    D_diffusion_opts = DiffusionParameters(mesh, kappa=kappa, mu=10/delta)
     diffusion_options = [
         ("u", u_diffusion_opts),
         ("v", v_diffusion_opts),
@@ -72,7 +72,7 @@ def shallow_water_1d_wave(
     ]
 
     # Equation
-    parameters = ShallowWaterParameters(H=1/epsilon, g=1/epsilon)
+    parameters = ShallowWaterParameters(mesh, H=1/epsilon, g=1/epsilon)
     eqns = ShallowWaterEquations_1d(
         domain, parameters, fexpr=Constant(1/epsilon),
         diffusion_options=diffusion_options
diff --git a/examples/shallow_water/thermal_williamson_2.py b/examples/shallow_water/thermal_williamson_2.py
index 7adc305dc..318f44962 100644
--- a/examples/shallow_water/thermal_williamson_2.py
+++ b/examples/shallow_water/thermal_williamson_2.py
@@ -68,7 +68,7 @@ def thermal_williamson_2(
     x, y, z = SpatialCoordinate(mesh)
 
     # Equations
-    params = ShallowWaterParameters(H=mean_depth, g=g)
+    params = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = params.Omega
     fexpr = 2*Omega*z/radius
     eqns = ThermalShallowWaterEquations(
diff --git a/examples/shallow_water/williamson_2.py b/examples/shallow_water/williamson_2.py
index 8c5f60917..f362605bd 100644
--- a/examples/shallow_water/williamson_2.py
+++ b/examples/shallow_water/williamson_2.py
@@ -62,7 +62,7 @@ def williamson_2(
     xyz = SpatialCoordinate(mesh)
 
     # Equation
-    parameters = ShallowWaterParameters(H=mean_depth)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth)
     Omega = parameters.Omega
     _, lat, _ = rotated_lonlatr_coords(xyz, rotate_pole_to)
     e_lon, _, _ = rotated_lonlatr_vectors(xyz, rotate_pole_to)
diff --git a/examples/shallow_water/williamson_5.py b/examples/shallow_water/williamson_5.py
index efe41b183..13f48a691 100644
--- a/examples/shallow_water/williamson_5.py
+++ b/examples/shallow_water/williamson_5.py
@@ -65,7 +65,7 @@ def williamson_5(
     lamda, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equation: coriolis
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 
diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py
index 40852033a..46ee2c557 100644
--- a/gusto/core/configuration.py
+++ b/gusto/core/configuration.py
@@ -1,7 +1,7 @@
 """Some simple tools for configuring the model."""
-from abc import ABCMeta, abstractproperty
+from abc import ABCMeta, abstractmethod
 from enum import Enum
-from firedrake import sqrt, Constant
+from firedrake import sqrt, Function, FunctionSpace
 
 
 __all__ = [
@@ -50,11 +50,14 @@ class TransportEquationType(Enum):
 class Configuration(object):
     """A base configuration object, for storing aspects of the model."""
 
-    def __init__(self, **kwargs):
+    mesh = None
+
+    def __init__(self, mesh=None, **kwargs):
         """
         Args:
             **kwargs: attributes and their values to be stored in the object.
         """
+        self.mesh = mesh
         for name, value in kwargs.items():
             self.__setattr__(name, value)
 
@@ -64,7 +67,7 @@ def __setattr__(self, name, value):
 
         When attributes are provided as floats or integers, these are converted
         to Firedrake :class:`Constant` objects, other than a handful of special
-        integers (dumpfreq, pddumpfreq, chkptfreq and log_level).
+        integers.
 
         Args:
             name: the attribute's name.
@@ -77,14 +80,18 @@ def __setattr__(self, name, value):
         if not hasattr(self, name):
             raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name))
 
-        # Almost all parameters should be Constants -- but there are some
-        # specific exceptions which should be kept as integers
+        # Almost all parameters should be functions on the real space
+        # -- but there are some specific exceptions which should be
+        # kept as integers
+        if self.mesh is not None:
+            R = FunctionSpace(self.mesh, 'R', 0)
         non_constants = [
             'dumpfreq', 'pddumpfreq', 'chkptfreq',
             'fixed_subcycles', 'max_subcycles', 'subcycle_by_courant'
         ]
         if type(value) in [float, int] and name not in non_constants:
-            object.__setattr__(self, name, Constant(value))
+            print(name, type(value))
+            object.__setattr__(self, name, Function(R, val=float(value)))
         else:
             object.__setattr__(self, name, value)
 
@@ -167,7 +174,7 @@ class ShallowWaterParameters(Configuration):
 class WrapperOptions(Configuration, metaclass=ABCMeta):
     """Base class for specifying options for a transport scheme."""
 
-    @abstractproperty
+    @abstractmethod
     def name(self):
         pass
 
diff --git a/gusto/core/domain.py b/gusto/core/domain.py
index 27addd9fa..2caedfab3 100644
--- a/gusto/core/domain.py
+++ b/gusto/core/domain.py
@@ -65,15 +65,16 @@ def __init__(self, mesh, dt, family, degree=None,
         # -------------------------------------------------------------------- #
 
         # Store central dt for use in the rest of the model
+        R = FunctionSpace(mesh, "R", 0)
         if type(dt) is Constant:
-            self.dt = dt
+            self.dt = Function(R, val=float(dt))
         elif type(dt) in (float, int):
-            self.dt = Constant(dt)
+            self.dt = Function(R, val=dt)
         else:
             raise TypeError(f'dt must be a Constant, float or int, not {type(dt)}')
 
         # Make a placeholder for the time
-        self.t = Constant(0.0)
+        self.t = Function(R, val=0.0)
 
         # -------------------------------------------------------------------- #
         # Build compatible function spaces
diff --git a/gusto/core/io.py b/gusto/core/io.py
index bd80f86d5..1e69a8651 100644
--- a/gusto/core/io.py
+++ b/gusto/core/io.py
@@ -8,7 +8,8 @@
 from gusto.diagnostics import Diagnostics, CourantNumber
 from gusto.core.meshes import get_flat_latlon_mesh
 from firedrake import (Function, functionspaceimpl, Constant, COMM_WORLD,
-                       DumbCheckpoint, FILE_CREATE, FILE_READ, CheckpointFile)
+                       DumbCheckpoint, FILE_CREATE, FILE_READ, CheckpointFile,
+                       MeshGeometry)
 from firedrake.output import VTKFile
 from pyop2.mpi import MPI
 import numpy as np
@@ -257,7 +258,7 @@ def log_parameters(self, equation):
         """
         if hasattr(equation, 'parameters') and equation.parameters is not None:
             logger.info("Physical parameters that take non-default values:")
-            logger.info(", ".join("%s: %s" % (k, float(v)) for (k, v) in vars(equation.parameters).items()))
+            logger.info(", ".join("%s: %s" % (k, float(v)) for (k, v) in vars(equation.parameters).items() if type(v) is not MeshGeometry))
 
     def setup_log_courant(self, state_fields, name='u', component="whole",
                           expression=None):
diff --git a/gusto/equations/boussinesq_equations.py b/gusto/equations/boussinesq_equations.py
index b9823f110..d94e2a038 100644
--- a/gusto/equations/boussinesq_equations.py
+++ b/gusto/equations/boussinesq_equations.py
@@ -168,10 +168,12 @@ def __init__(self, domain, parameters,
         # -------------------------------------------------------------------- #
         if compressible:
             cs = parameters.cs
+            # On assuming ``cs`` as a constant, it is right keep it out of the
+            # integration.
             linear_div_form = divergence(subject(
-                prognostic(cs**2 * phi * div(u_trial) * dx, 'p'), self.X))
+                prognostic(cs**2 * (phi * div(u_trial) * dx), 'p'), self.X))
             divergence_form = divergence(linearisation(
-                subject(prognostic(cs**2 * phi * div(u) * dx, 'p'), self.X),
+                subject(prognostic(cs**2 * (phi * div(u) * dx), 'p'), self.X),
                 linear_div_form))
         else:
             # This enforces that div(u) = 0
diff --git a/gusto/equations/compressible_euler_equations.py b/gusto/equations/compressible_euler_equations.py
index 73bc76d45..fa805ae93 100644
--- a/gusto/equations/compressible_euler_equations.py
+++ b/gusto/equations/compressible_euler_equations.py
@@ -17,7 +17,6 @@
 )
 from gusto.equations.active_tracers import Phases, TracerVariableType
 from gusto.equations.prognostic_equations import PrognosticEquationSet
-
 __all__ = ["CompressibleEulerEquations", "HydrostaticCompressibleEulerEquations"]
 
 
@@ -45,7 +44,7 @@ def __init__(self, domain, parameters, sponge_options=None,
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
                 mesh and the compatible function spaces.
-            parameters (:class:`Configuration`, optional): an object containing
+            x (:class:`Configuration`, optional): an object containing
                 the model's physical parameters.
             sponge_options (:class:`SpongeLayerParameters`, optional): any
                 parameters for applying a sponge layer to the upper boundary.
@@ -109,6 +108,7 @@ def __init__(self, domain, parameters, sponge_options=None,
         u_trial = split(self.trials)[0]
         _, rho_bar, theta_bar = split(self.X_ref)[0:3]
         zero_expr = Constant(0.0)*theta
+        # Check this for adjoints
         exner = exner_pressure(parameters, rho, theta)
         n = FacetNormal(domain.mesh)
 
diff --git a/gusto/equations/shallow_water_equations.py b/gusto/equations/shallow_water_equations.py
index fd61cedf7..0e9aebe2a 100644
--- a/gusto/equations/shallow_water_equations.py
+++ b/gusto/equations/shallow_water_equations.py
@@ -163,8 +163,9 @@ def _setup_residual(self, fexpr, topog_expr, u_transport_option):
         # -------------------------------------------------------------------- #
         # Pressure Gradient Term
         # -------------------------------------------------------------------- #
+        # On assuming ``g``, it is right to keep it out of the integral.
         pressure_gradient_form = pressure_gradient(
-            subject(prognostic(-g*div(w)*D*dx, 'u'), self.X))
+            subject(prognostic(-g*(div(w)*D*dx), 'u'), self.X))
 
         residual = (mass_form + adv_form + pressure_gradient_form)
 
@@ -418,6 +419,7 @@ def _setup_residual(self, fexpr, topog_expr, u_transport_option):
         # provide linearisation
         if self.equivalent_buoyancy:
             beta2 = self.parameters.beta2
+
             qsat_expr = self.compute_saturation(self.X)
             qv = conditional(qt < qsat_expr, qt, qsat_expr)
             qvbar = conditional(qtbar < qsat_expr, qtbar, qsat_expr)
diff --git a/gusto/physics/shallow_water_microphysics.py b/gusto/physics/shallow_water_microphysics.py
index 29b3637bd..f30685fa0 100644
--- a/gusto/physics/shallow_water_microphysics.py
+++ b/gusto/physics/shallow_water_microphysics.py
@@ -3,7 +3,8 @@
 """
 
 from firedrake import (
-    conditional, Function, dx, min_value, max_value, Constant, assemble, split
+    conditional, Function, dx, min_value, max_value, FunctionSpace,
+    assemble, split
 )
 from firedrake.__future__ import interpolate
 from firedrake.fml import subject
@@ -97,13 +98,14 @@ def __init__(self, equation, saturation_curve,
         self.source_expr = split(self.source)[self.Vv_idx]
         self.source_int = self.source.subfunctions[self.Vv_idx]
 
+        R = FunctionSpace(equation.domain.mesh, "R", 0)
         # tau is the timescale for conversion (may or may not be the timestep)
         if tau is not None:
             self.set_tau_to_dt = False
-            self.tau = tau
+            self.tau = Function(R).assign(tau)
         else:
             self.set_tau_to_dt = True
-            self.tau = Constant(0)
+            self.tau = Function(R)
             logger.info("Timescale for rain conversion has been set to dt. If this is not the intention then provide a tau parameter as an argument to InstantRain.")
 
         if self.time_varying_saturation:
@@ -280,12 +282,13 @@ def __init__(self, equation, saturation_curve,
             V_idxs.append(self.Vb_idx)
 
         # tau is the timescale for condensation/evaporation (may or may not be the timestep)
+        R = FunctionSpace(equation.domain.mesh, "R", 0)
         if tau is not None:
             self.set_tau_to_dt = False
-            self.tau = tau
+            self.tau = Function(R).assign(tau)
         else:
             self.set_tau_to_dt = True
-            self.tau = Constant(0)
+            self.tau = Function(R)
             logger.info("Timescale for moisture conversion between vapour and cloud has been set to dt. If this is not the intention then provide a tau parameter as an argument to SWSaturationAdjustment.")
 
         if self.time_varying_saturation:
diff --git a/gusto/spatial_methods/diffusion_methods.py b/gusto/spatial_methods/diffusion_methods.py
index 768b7532e..bc33da4df 100644
--- a/gusto/spatial_methods/diffusion_methods.py
+++ b/gusto/spatial_methods/diffusion_methods.py
@@ -161,4 +161,5 @@ def __init__(self, equation, variable, diffusion_parameters):
             kappa = diffusion_parameters.kappa
             self.form = diffusion(kappa * self.test.dx(0) * self.field.dx(0) * dx)
         else:
-            raise NotImplementedError("CG diffusion only implemented in 1D")
+            kappa = diffusion_parameters.kappa
+            self.form = diffusion(kappa * inner(grad(self.test), grad(self.field)) * dx)
diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py
index db71412ca..7b5eee9a6 100644
--- a/gusto/time_discretisation/time_discretisation.py
+++ b/gusto/time_discretisation/time_discretisation.py
@@ -9,8 +9,8 @@
 import math
 
 from firedrake import (Function, TestFunction, TestFunctions, DirichletBC,
-                       Constant, NonlinearVariationalProblem,
-                       NonlinearVariationalSolver)
+                       NonlinearVariationalProblem, NonlinearVariationalSolver,
+                       FunctionSpace)
 from firedrake.fml import (replace_subject, replace_test_function, Term,
                            all_terms, drop)
 from firedrake.formmanipulation import split_form
@@ -88,10 +88,10 @@ def __init__(self, domain, field_name=None, subcycling_options=None,
         self.domain = domain
         self.field_name = field_name
         self.equation = None
-
-        self.dt = Constant(0.0)
+        R = FunctionSpace(domain.mesh, "R", 0)
+        self.dt = Function(R, val=0.0)
         self.dt.assign(domain.dt)
-        self.original_dt = Constant(0.0)
+        self.original_dt = Function(R, val=0.0)
         self.original_dt.assign(self.dt)
         self.options = options
         self.limiter = limiter
diff --git a/gusto/timestepping/timestepper.py b/gusto/timestepping/timestepper.py
index 24a2058f7..f4f5bba0a 100644
--- a/gusto/timestepping/timestepper.py
+++ b/gusto/timestepping/timestepper.py
@@ -235,7 +235,7 @@ def run(self, t, tmax, pick_up=False):
 
             self.timestep()
 
-            self.t.assign(self.t + self.dt)
+            self.t.assign(float(self.t) + float(self.dt))
             self.step += 1
 
             with timed_stage("Dump output"):
diff --git a/integration-tests/adjoints/test_diffusion_sensitivity.py b/integration-tests/adjoints/test_diffusion_sensitivity.py
new file mode 100644
index 000000000..7c353eef9
--- /dev/null
+++ b/integration-tests/adjoints/test_diffusion_sensitivity.py
@@ -0,0 +1,78 @@
+import pytest
+import numpy as np
+
+from firedrake import *
+from firedrake.adjoint import *
+from pyadjoint import get_working_tape
+from gusto import *
+
+
+@pytest.fixture(autouse=True)
+def handle_taping():
+    yield
+    tape = get_working_tape()
+    tape.clear_tape()
+
+
+@pytest.fixture(autouse=True, scope="module")
+def handle_annotation():
+    from firedrake.adjoint import annotate_tape, continue_annotation
+    if not annotate_tape():
+        continue_annotation()
+    yield
+    # Ensure annotation is paused when we finish.
+    annotate = annotate_tape()
+    if annotate:
+        pause_annotation()
+
+
+@pytest.mark.parametrize("nu_is_control", [True, False])
+def test_diffusion_sensitivity(nu_is_control, tmpdir):
+    assert get_working_tape()._blocks == []
+    n = 30
+    mesh = PeriodicUnitSquareMesh(n, n)
+    output = OutputParameters(dirname=str(tmpdir))
+    dt = 0.01
+    domain = Domain(mesh, 10*dt, family="BDM", degree=1)
+    io = IO(domain, output)
+
+    V = VectorFunctionSpace(mesh, "CG", 2)
+    domain.spaces.add_space("vecCG", V)
+
+    R = FunctionSpace(mesh, "R", 0)
+    # We need to define nu as a function in order to have a control variable.
+    nu = Function(R, val=0.0001)
+    diffusion_params = DiffusionParameters(mesh, kappa=nu)
+    eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params)
+
+    diffusion_scheme = BackwardEuler(domain)
+    diffusion_methods = [CGDiffusion(eqn, "f", diffusion_params)]
+    timestepper = Timestepper(eqn, diffusion_scheme, io, spatial_methods=diffusion_methods)
+
+    x = SpatialCoordinate(mesh)
+    fexpr = as_vector((sin(2*pi*x[0]), cos(2*pi*x[1])))
+    timestepper.fields("f").interpolate(fexpr)
+
+    end = 0.1
+    timestepper.run(0., end)
+
+    u = timestepper.fields("f")
+    J = assemble(inner(u, u)*dx)
+
+    if nu_is_control:
+        control = Control(nu)
+        h = Function(R, val=0.0001)  # the direction of the perturbation
+    else:
+        control = Control(u)
+        # the direction of the perturbation
+        h = Function(V).interpolate(fexpr * np.random.rand())
+
+    # the functional as a pure function of nu
+    Jhat = ReducedFunctional(J, control)
+
+    if nu_is_control:
+        assert np.allclose(J, Jhat(nu))
+        assert taylor_test(Jhat, nu, h) > 1.95
+    else:
+        assert np.allclose(J, Jhat(u))
+        assert taylor_test(Jhat, u, h) > 1.95
diff --git a/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py b/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py
new file mode 100644
index 000000000..17e527d6f
--- /dev/null
+++ b/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py
@@ -0,0 +1,218 @@
+"""
+The moist thermal form of Test Case 5 (flow over a mountain) of Williamson et
+al, 1992:
+``A standard test set for numerical approximations to the shallow water
+equations in spherical geometry'', JCP.
+
+The initial conditions are taken from Zerroukat & Allen, 2015:
+``A moist Boussinesq shallow water equations set for testing atmospheric
+models'', JCP.
+
+Three moist variables (vapour, cloud liquid and rain) are used. This set of
+equations involves an active buoyancy field.
+
+The example here uses the icosahedral sphere mesh and degree 1 spaces. An
+explicit RK4 timestepper is used.
+"""
+import pytest
+import numpy as np
+
+from firedrake import (
+    SpatialCoordinate, as_vector, pi, sqrt, min_value, exp, cos, sin, assemble, dx, inner, Function
+)
+from firedrake.adjoint import *
+from pyadjoint import get_working_tape
+from gusto import (
+    Domain, IO, OutputParameters, Timestepper, RK4, DGUpwind,
+    ShallowWaterParameters, ThermalShallowWaterEquations, lonlatr_from_xyz,
+    InstantRain, SWSaturationAdjustment, WaterVapour, CloudWater,
+    Rain, GeneralIcosahedralSphereMesh
+)
+
+
+@pytest.fixture(autouse=True)
+def handle_taping():
+    yield
+    tape = get_working_tape()
+    tape.clear_tape()
+
+
+@pytest.fixture(autouse=True, scope="module")
+def handle_annotation():
+    from firedrake.adjoint import annotate_tape, continue_annotation
+    if not annotate_tape():
+        continue_annotation()
+    yield
+    # Ensure annotation is paused when we finish.
+    annotate = annotate_tape()
+    if annotate:
+        pause_annotation()
+
+
+def test_moist_thermal_williamson_5_sensitivity(
+        tmpdir, ncells_per_edge=8, dt=600, tmax=50.*24.*60.*60.,
+        dumpfreq=2880
+):
+    assert get_working_tape()._blocks == []
+    # ------------------------------------------------------------------------ #
+    # Parameters for test case
+    # ------------------------------------------------------------------------ #
+    radius = 6371220.           # planetary radius (m)
+    mean_depth = 5960           # reference depth (m)
+    g = 9.80616                 # acceleration due to gravity (m/s^2)
+    u_max = 20.                 # max amplitude of the zonal wind (m/s)
+    epsilon = 1/300             # linear air expansion coeff (1/K)
+    theta_SP = -40*epsilon      # value of theta at south pole (no units)
+    theta_EQ = 30*epsilon       # value of theta at equator (no units)
+    theta_NP = -20*epsilon      # value of theta at north pole (no units)
+    mu1 = 0.05                  # scaling for theta with longitude (no units)
+    mu2 = 0.98                  # proportion of qsat to make init qv (no units)
+    q0 = 135                    # qsat scaling, gives init q_v of ~0.02, (kg/kg)
+    beta2 = 10*g                # buoyancy-vaporisation factor (m/s^2)
+    nu = 20.                    # qsat factor in exponent (no units)
+    qprecip = 1e-4              # cloud to rain conversion threshold (kg/kg)
+    gamma_r = 1e-3              # rain-coalescence implicit factor
+    mountain_height = 2000.     # height of mountain (m)
+    R0 = pi/9.                  # radius of mountain (rad)
+    lamda_c = -pi/2.            # longitudinal centre of mountain (rad)
+    phi_c = pi/6.               # latitudinal centre of mountain (rad)
+
+    # ------------------------------------------------------------------------ #
+    # Our settings for this set up
+    # ------------------------------------------------------------------------ #
+
+    element_order = 1
+    u_eqn_type = 'vector_invariant_form'
+
+    # ------------------------------------------------------------------------ #
+    # Set up model objects
+    # ------------------------------------------------------------------------ #
+
+    # Domain
+    mesh = GeneralIcosahedralSphereMesh(radius, ncells_per_edge, degree=2)
+    domain = Domain(mesh, dt, "BDM", element_order)
+    x, y, z = SpatialCoordinate(mesh)
+    lamda, phi, _ = lonlatr_from_xyz(x, y, z)
+
+    # Equation: coriolis
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
+    Omega = parameters.Omega
+    fexpr = 2*Omega*z/radius
+
+    # Equation: topography
+    rsq = min_value(R0**2, (lamda - lamda_c)**2 + (phi - phi_c)**2)
+    r = sqrt(rsq)
+    tpexpr = mountain_height * (1 - r/R0)
+
+    # Equation: moisture
+    tracers = [
+        WaterVapour(space='DG'), CloudWater(space='DG'), Rain(space='DG')
+    ]
+    eqns = ThermalShallowWaterEquations(
+        domain, parameters, fexpr=fexpr, topog_expr=tpexpr,
+        active_tracers=tracers, u_transport_option=u_eqn_type
+    )
+
+    # I/O
+    output = OutputParameters(
+        dirname=str(tmpdir), dumplist_latlon=['D'], dumpfreq=dumpfreq,
+        dump_vtus=True, dump_nc=False, log_courant=False,
+        dumplist=['D', 'b', 'water_vapour', 'cloud_water']
+    )
+    io = IO(domain, output)
+
+    # Physics ------------------------------------------------------------------
+    # Saturation function -- first define simple expression
+    def q_sat(b, D):
+        return (q0/(g*D + g*tpexpr)) * exp(nu*(1 - b/g))
+
+    # Function to pass to physics (takes mixed function as argument)
+    def phys_sat_func(x_in):
+        D = x_in.sub(1)
+        b = x_in.sub(2)
+        return q_sat(b, D)
+
+    # Feedback proportionality is dependent on D and b
+    def gamma_v(x_in):
+        D = x_in.sub(1)
+        b = x_in.sub(2)
+        return 1.0 / (1.0 + nu*beta2/g*q_sat(b, D))
+
+    SWSaturationAdjustment(
+        eqns, phys_sat_func, time_varying_saturation=True,
+        parameters=parameters, thermal_feedback=True,
+        beta2=beta2, gamma_v=gamma_v, time_varying_gamma_v=True
+    )
+
+    InstantRain(
+        eqns, qprecip, vapour_name="cloud_water", rain_name="rain",
+        gamma_r=gamma_r
+    )
+
+    transport_methods = [
+        DGUpwind(eqns, field_name) for field_name in eqns.field_names
+    ]
+
+    # Timestepper
+    stepper = Timestepper(
+        eqns, RK4(domain), io, spatial_methods=transport_methods
+    )
+
+    # ------------------------------------------------------------------------ #
+    # Initial conditions
+    # ------------------------------------------------------------------------ #
+
+    u0 = stepper.fields("u")
+    D0 = stepper.fields("D")
+    b0 = stepper.fields("b")
+    v0 = stepper.fields("water_vapour")
+    c0 = stepper.fields("cloud_water")
+    r0 = stepper.fields("rain")
+
+    uexpr = as_vector([-u_max*y/radius, u_max*x/radius, 0.0])
+
+    Dexpr = (
+        mean_depth - tpexpr
+        - (radius * Omega * u_max + 0.5*u_max**2)*(z/radius)**2/g
+    )
+
+    # Expression for initial buoyancy - note the bracket around 1-mu
+    theta_expr = (
+        2/(pi**2) * (
+            phi*(phi - pi/2)*theta_SP
+            - 2*(phi + pi/2) * (phi - pi/2)*(1 - mu1)*theta_EQ
+            + phi*(phi + pi/2)*theta_NP
+        )
+        + mu1*theta_EQ*cos(phi)*sin(lamda)
+    )
+    bexpr = g * (1 - theta_expr)
+
+    # Expression for initial vapour depends on initial saturation
+    vexpr = mu2 * q_sat(bexpr, Dexpr)
+
+    # Initialise (cloud and rain initially zero)
+    u0.project(uexpr)
+    D0.interpolate(Dexpr)
+    b0.interpolate(bexpr)
+    v0.interpolate(vexpr)
+    c0.assign(0.0)
+    r0.assign(0.0)
+
+    # ----------------------------------------------------------------- #
+    # Run
+    # ----------------------------------------------------------------- #
+    stepper.run(t=0, tmax=dt)
+
+    u_tf = stepper.fields('u')  # Final velocity field
+    D_tf = stepper.fields('D')  # Final depth field
+
+    J = assemble(0.5*inner(u_tf, u_tf)*dx + 0.5*g*D_tf**2*dx)
+
+    Jhat = ReducedFunctional(J, Control(D0))
+
+    assert np.allclose(Jhat(D0), J)
+    with stop_annotating():
+        # Stop annotation to perform the Taylor test
+        h0 = Function(D0.function_space())
+        h0.assign(D0 * np.random.rand())
+        assert taylor_test(Jhat, D0, h0) > 1.95
diff --git a/integration-tests/adjoints/test_shallow_water_sensitivity.py b/integration-tests/adjoints/test_shallow_water_sensitivity.py
new file mode 100644
index 000000000..22ffb9681
--- /dev/null
+++ b/integration-tests/adjoints/test_shallow_water_sensitivity.py
@@ -0,0 +1,100 @@
+import pytest
+import numpy as np
+
+from firedrake import *
+from firedrake.adjoint import *
+from pyadjoint import get_working_tape
+from gusto import *
+
+
+@pytest.fixture(autouse=True)
+def handle_taping():
+    yield
+    tape = get_working_tape()
+    tape.clear_tape()
+
+
+@pytest.fixture(autouse=True, scope="module")
+def handle_annotation():
+    from firedrake.adjoint import annotate_tape, continue_annotation
+    if not annotate_tape():
+        continue_annotation()
+    yield
+    # Ensure annotation is paused when we finish.
+    annotate = annotate_tape()
+    if annotate:
+        pause_annotation()
+
+
+def test_shallow_water(tmpdir):
+    assert get_working_tape()._blocks == []
+    # setup shallow water parameters
+    R = 6371220.
+    H = 5960.
+    dt = 900.
+
+    # Domain
+    mesh = IcosahedralSphereMesh(radius=R, refinement_level=3, degree=2)
+    x = SpatialCoordinate(mesh)
+    domain = Domain(mesh, dt, 'BDM', 1)
+    parameters = ShallowWaterParameters(mesh, H=H)
+
+    # Equation
+    Omega = parameters.Omega
+    fexpr = 2*Omega*x[2]/R
+    lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2])
+    R0 = pi/9.
+    R0sq = R0**2
+    lamda_c = -pi/2.
+    lsq = (lamda - lamda_c)**2
+    theta_c = pi/6.
+    thsq = (theta - theta_c)**2
+    rsq = min_value(R0sq, lsq+thsq)
+    r = sqrt(rsq)
+    bexpr = 2000 * (1 - r/R0)
+    eqn = ShallowWaterEquations(domain, parameters, fexpr=fexpr, topog_expr=bexpr)
+
+    # I/O
+    output = OutputParameters(dirname=str(tmpdir), log_courant=False)
+    io = IO(domain, output)
+
+    # Transport schemes
+    transported_fields = [TrapeziumRule(domain, "u"), SSPRK3(domain, "D")]
+    transport_methods = [DGUpwind(eqn, "u"), DGUpwind(eqn, "D")]
+
+    # Time stepper
+    stepper = SemiImplicitQuasiNewton(
+        eqn, io, transported_fields, transport_methods
+    )
+
+    u0 = stepper.fields('u')
+    D0 = stepper.fields('D')
+    u_max = 20.   # Maximum amplitude of the zonal wind (m/s)
+    uexpr = as_vector([-u_max*x[1]/R, u_max*x[0]/R, 0.0])
+    g = parameters.g
+    Rsq = R**2
+    Dexpr = H - ((R * Omega * u_max + 0.5*u_max**2)*x[2]**2/Rsq)/g - bexpr
+
+    u0.project(uexpr)
+    D0.interpolate(Dexpr)
+
+    Dbar = Function(D0.function_space()).assign(H)
+    stepper.set_reference_profiles([('D', Dbar)])
+
+    stepper.run(0., 5*dt)
+
+    u_tf = stepper.fields('u')  # Final velocity field
+    D_tf = stepper.fields('D')  # Final depth field
+
+    J = assemble(0.5*inner(u_tf, u_tf)*dx + 0.5*g*D_tf**2*dx)
+
+    control = [Control(D0), Control(u0)]  # Control variables
+    J_hat = ReducedFunctional(J, control)
+    assert np.isclose(J_hat([D0, u0]), J, rtol=1e-10)
+    with stop_annotating():
+        # Stop annotation to perform the Taylor test
+        h0 = Function(D0.function_space())
+        h1 = Function(u0.function_space())
+        h0.assign(D0 * np.random.rand())
+        h1.assign(u0 * np.random.rand())
+        assert taylor_test(J_hat, [D0, u0], [h0, h1]) > 1.95
diff --git a/integration-tests/balance/test_boussinesq_balance.py b/integration-tests/balance/test_boussinesq_balance.py
index 3fae46620..21e55d9e4 100644
--- a/integration-tests/balance/test_boussinesq_balance.py
+++ b/integration-tests/balance/test_boussinesq_balance.py
@@ -31,7 +31,7 @@ def setup_balance(dirname):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    parameters = BoussinesqParameters()
+    parameters = BoussinesqParameters(mesh)
     eqns = BoussinesqEquations(domain, parameters)
 
     # I/O
diff --git a/integration-tests/balance/test_compressible_balance.py b/integration-tests/balance/test_compressible_balance.py
index a9726dd5a..29edcfa8d 100644
--- a/integration-tests/balance/test_compressible_balance.py
+++ b/integration-tests/balance/test_compressible_balance.py
@@ -31,7 +31,7 @@ def setup_balance(dirname):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(domain, parameters)
 
     # I/O
diff --git a/integration-tests/balance/test_saturated_balance.py b/integration-tests/balance/test_saturated_balance.py
index 2c89da089..42124bcb6 100644
--- a/integration-tests/balance/test_saturated_balance.py
+++ b/integration-tests/balance/test_saturated_balance.py
@@ -39,7 +39,7 @@ def setup_saturated(dirname, recovered):
         u_transport_option = "vector_advection_form"
     else:
         u_transport_option = "vector_invariant_form"
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, parameters, u_transport_option=u_transport_option, active_tracers=tracers)
 
diff --git a/integration-tests/balance/test_unsaturated_balance.py b/integration-tests/balance/test_unsaturated_balance.py
index c6e3f33a9..d3a6463f8 100644
--- a/integration-tests/balance/test_unsaturated_balance.py
+++ b/integration-tests/balance/test_unsaturated_balance.py
@@ -42,7 +42,7 @@ def setup_unsaturated(dirname, recovered):
         u_transport_option = "vector_advection_form"
     else:
         u_transport_option = "vector_invariant_form"
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, parameters, u_transport_option=u_transport_option, active_tracers=tracers)
 
diff --git a/integration-tests/diffusion/test_diffusion.py b/integration-tests/diffusion/test_diffusion.py
index fd570107d..15ea87be4 100644
--- a/integration-tests/diffusion/test_diffusion.py
+++ b/integration-tests/diffusion/test_diffusion.py
@@ -32,7 +32,7 @@ def test_scalar_diffusion(tmpdir, DG, tracer_setup):
 
     mu = 5.
 
-    diffusion_params = DiffusionParameters(kappa=kappa, mu=mu)
+    diffusion_params = DiffusionParameters(mesh, kappa=kappa, mu=mu)
     eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params)
     diffusion_scheme = BackwardEuler(domain)
     diffusion_methods = [InteriorPenaltyDiffusion(eqn, "f", diffusion_params)]
diff --git a/integration-tests/equations/test_advection_diffusion.py b/integration-tests/equations/test_advection_diffusion.py
index 778eb9f73..aac3ab2ee 100644
--- a/integration-tests/equations/test_advection_diffusion.py
+++ b/integration-tests/equations/test_advection_diffusion.py
@@ -22,7 +22,7 @@ def run_advection_diffusion(tmpdir):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    diffusion_params = DiffusionParameters(kappa=0.75, mu=5)
+    diffusion_params = DiffusionParameters(mesh, kappa=0.75, mu=5)
     V = domain.spaces("DG")
     Vu = VectorFunctionSpace(mesh, "CG", 1)
 
diff --git a/integration-tests/equations/test_boussinesq.py b/integration-tests/equations/test_boussinesq.py
index 2655d87d3..8de65aef6 100644
--- a/integration-tests/equations/test_boussinesq.py
+++ b/integration-tests/equations/test_boussinesq.py
@@ -31,7 +31,7 @@ def run_boussinesq(tmpdir, compressible):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    parameters = BoussinesqParameters()
+    parameters = BoussinesqParameters(mesh)
     eqn = BoussinesqEquations(domain, parameters, compressible=compressible)
 
     # I/O
diff --git a/integration-tests/equations/test_dry_compressible.py b/integration-tests/equations/test_dry_compressible.py
index 846d38bbf..282667f97 100644
--- a/integration-tests/equations/test_dry_compressible.py
+++ b/integration-tests/equations/test_dry_compressible.py
@@ -30,7 +30,7 @@ def run_dry_compressible(tmpdir):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     eqn = CompressibleEulerEquations(domain, parameters)
 
     # I/O
diff --git a/integration-tests/equations/test_linear_sw_wave.py b/integration-tests/equations/test_linear_sw_wave.py
index 95eecdca4..e7cd59d74 100644
--- a/integration-tests/equations/test_linear_sw_wave.py
+++ b/integration-tests/equations/test_linear_sw_wave.py
@@ -28,7 +28,7 @@ def run_linear_sw_wave(tmpdir):
     domain = Domain(mesh, dt, 'BDM', 1)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H, g=g)
+    parameters = ShallowWaterParameters(mesh, H=H, g=g)
     fexpr = Constant(1)
     eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr)
 
diff --git a/integration-tests/equations/test_moist_compressible.py b/integration-tests/equations/test_moist_compressible.py
index f2909f75e..d02c0ced0 100644
--- a/integration-tests/equations/test_moist_compressible.py
+++ b/integration-tests/equations/test_moist_compressible.py
@@ -30,7 +30,7 @@ def run_moist_compressible(tmpdir):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     tracers = [WaterVapour(name='vapour_mixing_ratio'), CloudWater(name='cloud_liquid_mixing_ratio')]
     eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers)
 
diff --git a/integration-tests/equations/test_sw_fplane.py b/integration-tests/equations/test_sw_fplane.py
index 3d0acff5f..711adb7bb 100644
--- a/integration-tests/equations/test_sw_fplane.py
+++ b/integration-tests/equations/test_sw_fplane.py
@@ -23,7 +23,7 @@ def run_sw_fplane(tmpdir):
     # Equation
     H = 2
     g = 50
-    parameters = ShallowWaterParameters(H=H, g=g)
+    parameters = ShallowWaterParameters(mesh, H=H, g=g)
     f0 = 10
     fexpr = Constant(f0)
     eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr)
diff --git a/integration-tests/equations/test_sw_linear_triangle.py b/integration-tests/equations/test_sw_linear_triangle.py
index d533b3286..e34bb122f 100644
--- a/integration-tests/equations/test_sw_linear_triangle.py
+++ b/integration-tests/equations/test_sw_linear_triangle.py
@@ -30,7 +30,7 @@ def setup_sw(dirname):
     domain = Domain(mesh, dt, "BDM", degree=1)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H)
+    parameters = ShallowWaterParameters(mesh, H=H)
     Omega = parameters.Omega
     fexpr = 2*Omega*x[2]/R
     eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr)
diff --git a/integration-tests/equations/test_sw_triangle.py b/integration-tests/equations/test_sw_triangle.py
index 455595c31..bc25b2e63 100644
--- a/integration-tests/equations/test_sw_triangle.py
+++ b/integration-tests/equations/test_sw_triangle.py
@@ -36,7 +36,7 @@ def setup_sw(dirname, dt, u_transport_option):
     x = SpatialCoordinate(mesh)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H)
+    parameters = ShallowWaterParameters(mesh, H=H)
     Omega = parameters.Omega
     fexpr = 2*Omega*x[2]/R
     eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr,
diff --git a/integration-tests/equations/test_thermal_sw.py b/integration-tests/equations/test_thermal_sw.py
index f201e78c6..09b706fe9 100644
--- a/integration-tests/equations/test_thermal_sw.py
+++ b/integration-tests/equations/test_thermal_sw.py
@@ -37,7 +37,7 @@ def setup_sw(dirname, dt, u_transport_option):
     x = SpatialCoordinate(mesh)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H, g=g)
+    parameters = ShallowWaterParameters(mesh, H=H, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*x[2]/R
     eqns = ThermalShallowWaterEquations(domain, parameters, fexpr=fexpr,
diff --git a/integration-tests/model/test_checkpointing.py b/integration-tests/model/test_checkpointing.py
index dc72fb1de..09bb5ff16 100644
--- a/integration-tests/model/test_checkpointing.py
+++ b/integration-tests/model/test_checkpointing.py
@@ -15,7 +15,7 @@ def set_up_model_objects(mesh, dt, output, stepper_type, ref_update_freq):
 
     domain = Domain(mesh, dt, "CG", 1)
 
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(domain, parameters)
 
     # Have two diagnostic fields that depend on initial values -- check if
diff --git a/integration-tests/model/test_passive_tracer.py b/integration-tests/model/test_passive_tracer.py
index 41bd5608e..2482722b0 100644
--- a/integration-tests/model/test_passive_tracer.py
+++ b/integration-tests/model/test_passive_tracer.py
@@ -26,7 +26,7 @@ def run_tracer(setup):
 
     x = SpatialCoordinate(mesh)
     H = 0.1
-    parameters = ShallowWaterParameters(H=H)
+    parameters = ShallowWaterParameters(mesh, H=H)
     Omega = parameters.Omega
     g = parameters.g
     umax = setup.umax
diff --git a/integration-tests/model/test_simultaneous_SIQN.py b/integration-tests/model/test_simultaneous_SIQN.py
index ea3a12327..384062823 100644
--- a/integration-tests/model/test_simultaneous_SIQN.py
+++ b/integration-tests/model/test_simultaneous_SIQN.py
@@ -62,7 +62,7 @@ def run_simult_SIQN(tmpdir, order):
                           density_name='rho')]
 
     # Equation
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, params, active_tracers=tracers, u_transport_option=u_eqn_type
     )
diff --git a/integration-tests/model/test_split_timestepper.py b/integration-tests/model/test_split_timestepper.py
index 5f675cfd6..3f3eff94c 100644
--- a/integration-tests/model/test_split_timestepper.py
+++ b/integration-tests/model/test_split_timestepper.py
@@ -25,7 +25,7 @@ def run_split_timestepper_adv_diff_physics(tmpdir, timestepper):
     domain = Domain(mesh, dt, "CG", 1)
 
     # Equation
-    diffusion_params = DiffusionParameters(kappa=0.75, mu=5)
+    diffusion_params = DiffusionParameters(mesh, kappa=0.75, mu=5)
     V = domain.spaces("DG")
     Vu = VectorFunctionSpace(mesh, "CG", 1)
 
diff --git a/unit-tests/diagnostic_tests/test_cumulative_sum.py b/unit-tests/diagnostic_tests/test_cumulative_sum.py
index 295e522ca..918447fb9 100644
--- a/unit-tests/diagnostic_tests/test_cumulative_sum.py
+++ b/unit-tests/diagnostic_tests/test_cumulative_sum.py
@@ -18,7 +18,7 @@ def test_dewpoint():
     mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers)
 
     domain = Domain(mesh, 0.1, 'CG', 1)
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     active_tracers = [WaterVapour()]
     eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers)
     prog_fields = TimeLevelFields(eqn)
diff --git a/unit-tests/diagnostic_tests/test_dewpoint_temperature.py b/unit-tests/diagnostic_tests/test_dewpoint_temperature.py
index eaf1fbe55..1b4fb2980 100644
--- a/unit-tests/diagnostic_tests/test_dewpoint_temperature.py
+++ b/unit-tests/diagnostic_tests/test_dewpoint_temperature.py
@@ -18,7 +18,7 @@ def test_dewpoint():
     mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers)
 
     domain = Domain(mesh, 0.1, 'CG', 1)
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     active_tracers = [WaterVapour()]
     eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers)
     prog_fields = TimeLevelFields(eqn)
diff --git a/unit-tests/diagnostic_tests/test_wetbulb_temperature.py b/unit-tests/diagnostic_tests/test_wetbulb_temperature.py
index 5909fe710..1410a7b30 100644
--- a/unit-tests/diagnostic_tests/test_wetbulb_temperature.py
+++ b/unit-tests/diagnostic_tests/test_wetbulb_temperature.py
@@ -18,7 +18,7 @@ def test_wetbulb_temperature():
     mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers)
 
     domain = Domain(mesh, 0.1, 'CG', 1)
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     active_tracers = [WaterVapour()]
     eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers)
     prog_fields = TimeLevelFields(eqn)