diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst
index a4db823369..c1526c14e2 100644
--- a/manual/sphinx/user_docs/time_integration.rst
+++ b/manual/sphinx/user_docs/time_integration.rst
@@ -165,6 +165,9 @@ nonlinear solvers:
   `CVodeSetEpsLin
   <https://sundials.readthedocs.io/en/latest/cvodes/Usage/SIM.html#c.CVodeSetEpsLin>`_.
 
+The linear solver type can be set using the ``linear_solver`` option.
+Valid choices include ``gmres`` (the default), ``fgmres``, ``tfqmr``, ``bcgs``.
+
 IMEX-BDF2
 ---------
 
diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx
index f0d42d39bc..8b81a66f6d 100644
--- a/src/solver/impls/cvode/cvode.cxx
+++ b/src/solver/impls/cvode/cvode.cxx
@@ -51,6 +51,10 @@
 #include <cvode/cvode_bbdpre.h>
 #include <cvode/cvode_ls.h>
 
+#include <sunlinsol/sunlinsol_spbcgs.h>
+#include <sunlinsol/sunlinsol_spfgmr.h>
+#include <sunlinsol/sunlinsol_sptfqmr.h>
+
 #include <algorithm>
 #include <iterator>
 #include <numeric>
@@ -59,6 +63,8 @@
 BOUT_ENUM_CLASS(positivity_constraint, none, positive, non_negative, negative,
                 non_positive);
 
+BOUT_ENUM_CLASS(linear_solver, gmres, fgmres, tfqmr, bcgs);
+
 // NOLINTBEGIN(readability-identifier-length)
 namespace {
 int cvode_linear_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data);
@@ -344,7 +350,24 @@ int CvodeSolver::init() {
 
     const auto prectype =
         use_precon ? (rightprec ? SUN_PREC_RIGHT : SUN_PREC_LEFT) : SUN_PREC_NONE;
-    sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, maxl);
+
+    switch ((*options)["linear_solver"]
+                .doc("Set linear solver type. Default is gmres.")
+                .withDefault(linear_solver::gmres)) {
+    case linear_solver::gmres:
+      sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, maxl);
+      break;
+    case linear_solver::fgmres:
+      sun_solver = callWithSUNContext(SUNLinSol_SPFGMR, suncontext, uvec, prectype, maxl);
+      break;
+    case linear_solver::tfqmr:
+      sun_solver =
+          callWithSUNContext(SUNLinSol_SPTFQMR, suncontext, uvec, prectype, maxl);
+      break;
+    case linear_solver::bcgs:
+      sun_solver = callWithSUNContext(SUNLinSol_SPBCGS, suncontext, uvec, prectype, maxl);
+      break;
+    };
     if (sun_solver == nullptr) {
       throw BoutException("Creating SUNDIALS linear solver failed\n");
     }