diff --git a/.gitignore b/.gitignore
index a2c1be2..dbd06b9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ BOUT.log.*
 *.pdf
 *~
 .BOUT.pid*
+.vscode/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 68759a2..8da6256 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,16 +13,16 @@ set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 
 # Update submodules
 # Adapted from https://cliutils.gitlab.io/modern-cmake/chapters/projects/submodule.html
-find_package(Git QUIET)
-if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
-  message(STATUS "Submodule update")
-  execute_process(COMMAND ${GIT_EXECUTABLE} -c submodule.recurse=false submodule update --init --recursive
-    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
-    RESULT_VARIABLE GIT_SUBMOD_RESULT)
-  if(NOT GIT_SUBMOD_RESULT EQUAL "0")
-    message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules")
-  endif()
-endif()
+# find_package(Git QUIET)
+# if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
+#   message(STATUS "Submodule update")
+#   execute_process(COMMAND ${GIT_EXECUTABLE} -c submodule.recurse=false submodule update --init --recursive
+#     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+#     RESULT_VARIABLE GIT_SUBMOD_RESULT)
+#   if(NOT GIT_SUBMOD_RESULT EQUAL "0")
+#     message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules")
+#   endif()
+# endif()
 
 # Get the Git revision
 include(GetGitRevisionDescription)
@@ -33,7 +33,8 @@ endif()
 message(STATUS "Git revision: ${HERMES_REVISION}")
 
 # BOUT++ is a dependency
-add_subdirectory(external/BOUT-dev)
+# add_subdirectory(external/BOUT-dev)
+find_package(bout++ REQUIRED)
 
 set(HERMES_SOURCES
     src/hermes-2.cxx
diff --git a/hermes-2 b/hermes-2
new file mode 100755
index 0000000..d77bd10
Binary files /dev/null and b/hermes-2 differ
diff --git a/include/hermes-2.hxx b/include/hermes-2.hxx
index 8d471cf..74141be 100644
--- a/include/hermes-2.hxx
+++ b/include/hermes-2.hxx
@@ -45,10 +45,18 @@ protected:
   
   int precon(BoutReal t, BoutReal gamma, BoutReal delta);
 private:
+  // grid region boundaries
+  int jyseps1_1;
+  int jyseps1_2;
+  int jyseps2_1;
+  int jyseps2_2;
+
   // Equilibrium current
   Field2D Jpar0;
 
+  BoutReal nelim_floor; // Global density floor
   BoutReal nesheath_floor; // Density floor used in sheath boundary conditions
+  BoutReal kappa_epar_floor; // Electron heat conductivity floor
 
   // Evolving variables
   Field3D Ne;         // Electron density
@@ -90,6 +98,10 @@ private:
   // Impurity radiation
   BoutReal fimp;             // Impurity fraction (of Ne)
   bool impurity_adas;        // True if using ImpuritySpecies, false if using
+  bool impurity_split_y;       // Split the fixed fraction impurities by region
+  BoutReal fimp_id;          // fimp in the inner divertor (0 -> jyseps1_1)
+  BoutReal fimp_od;          // fimp in the outer divertor (jyseps2_2 -> -1)
+  BoutReal fimp_us;          // fimp in the upstream (between jyseps1_1 and jyseps2_2)
   ImpuritySpecies *impurity; // Atomicpp impurity
   
   BoutReal carbon_fraction;
@@ -142,6 +154,7 @@ private:
   bool sheath_closure; // Sheath closure sink on vorticity (if sinks = true)
   bool drift_wave;     // Drift-wave closure (if sinks=true)
 
+  bool slab_radial_buffers; // an alternatice radial buffer region for use in slab geometry
   bool radial_buffers; // Radial buffer regions
   int radial_inner_width; // Number of points in the inner radial buffer
   int radial_outer_width; // Number of points in the outer radial buffer
@@ -181,14 +194,16 @@ private:
   BoutReal numdiff, hyper, hyperpar; ///< Numerical dissipation
   int low_pass_z; // Fourier filter in Z 
   BoutReal z_hyper_viscos, x_hyper_viscos, y_hyper_viscos; // 4th-order derivatives
-  bool low_n_diffuse; // Diffusion in parallel direction at low density
-  bool low_n_diffuse_perp; // Diffusion in perpendicular direction at low density
+  BoutReal low_n_diffuse; // Diffusion in parallel direction at low density
+  BoutReal low_n_diffuse_perp; // Diffusion in perpendicular direction at low density
   BoutReal ne_hyper_z, pe_hyper_z; // Hyper-diffusion
+  BoutReal nvi_hyper_z, vepsi_hyper_z; // hyper-diffusion on ion terms
   BoutReal scale_num_cs; // Scale numerical sound speed
   BoutReal floor_num_cs; // Apply a floor to the numerical sound speed
   bool vepsi_dissipation; // Dissipation term in VePsi equation
   bool vort_dissipation; // Dissipation term in Vorticity equation
-  bool phi_dissipation; // Dissipation term in Vorticity equation, depending on phi
+  BoutReal phi_dissipation; // Dissipation term in Vorticity equation, depending on phi
+  BoutReal delp2_dissipation; // Dissipation using delp2
   
   // Sources and profiles
   
@@ -222,6 +237,7 @@ private:
   
   // Curvature, Grad-B drift
   Vector3D Curlb_B; // Curl(b/B)
+  bool revCurlb_B; // Reverse direction of Curl(b/B) vector
   
   // Perturbed parallel gradient operators
   const Field3D Grad_parP(const Field3D &f);
diff --git a/include/mixed.hxx b/include/mixed.hxx
index 68926b0..f37c32e 100644
--- a/include/mixed.hxx
+++ b/include/mixed.hxx
@@ -33,6 +33,8 @@ private:
   Field3D Pnlim;  // Limited pressure, used to calculate pressure-driven diffusive flows
   Field3D Vn;
   Field3D Dnn;
+
+  Field3D Riz, Rrc, Rcx, Rex;
   
   bool sheath_ydown, sheath_yup;
   
diff --git a/include/neutral-model.hxx b/include/neutral-model.hxx
index ed4ad45..c5d300a 100644
--- a/include/neutral-model.hxx
+++ b/include/neutral-model.hxx
@@ -26,6 +26,7 @@ public:
     
     // Options for calculating rates
     OPTION(options, Eionize, 30); // Energy loss per ionisation [eV]
+    OPTION(options, h_excitation, true); // Add to plasma radiation sink
   }
   virtual ~NeutralModel() {}
   
@@ -74,11 +75,13 @@ protected:
   UpdatedRadiatedPower hydrogen; // Atomic rates (H.Willett)
   
   BoutReal Eionize;   // Energy loss per ionisation [eV]
+
+  bool h_excitation; // add excitation energy to the hydrogen radiation energy sink
   
   void neutral_rates(const Field3D &Ne, const Field3D &Te, const Field3D &Ti, const Field3D &Vi,    // Plasma quantities
                      const Field3D &Nn, const Field3D &Tn, const Field3D &Vnpar, // Neutral gas
                      Field3D &S, Field3D &F, Field3D &Qi, Field3D &R,  // Transfer rates
-                     Field3D &Riz, Field3D &Rrc, Field3D &Rcx);        // Rates
+                     Field3D &Riz, Field3D &Rrc, Field3D &Rcx, Field3D &Rex);        // Rates
 
 private:
   NeutralModel();
diff --git a/include/revision.hxx b/include/revision.hxx
new file mode 100644
index 0000000..4d4f17c
--- /dev/null
+++ b/include/revision.hxx
@@ -0,0 +1,26 @@
+/// Information about the version of Hermes
+///
+/// The build system will update this file on every commit, which may
+/// result in files that include it getting rebuilt. Therefore it
+/// should be included in as few places as possible
+
+#ifndef HERMES_REVISION_H
+#define HERMES_REVISION_H
+
+namespace hermes {
+namespace version {
+/// The git commit hash
+#ifndef HERMES_REVISION
+constexpr auto revision = "@HERMES_REVISION@";
+#else
+// Stringify value passed at compile time
+#define BUILDFLAG1_(x) #x
+#define BUILDFLAG(x) BUILDFLAG1_(x)
+constexpr auto revision = BUILDFLAG(HERMES_REVISION);
+#undef BUILDFLAG1
+#undef BUILDFLAG
+#endif
+} // namespace version
+} // namespace hermes
+
+#endif // HERMES_REVISION_H
diff --git a/src/full-velocity.cxx b/src/full-velocity.cxx
index 3b54947..a72925a 100644
--- a/src/full-velocity.cxx
+++ b/src/full-velocity.cxx
@@ -324,8 +324,8 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te,
   /////////////////////////////////////////////////////
   // Atomic processes
 
-  Field3D Riz, Rrc, Rcx;
-  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vnpar, S, F, Qi, Rp, Riz, Rrc, Rcx);
+  Field3D Riz, Rrc, Rcx, Rex;
+  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vnpar, S, F, Qi, Rp, Riz, Rrc, Rcx, Rex);
 
   Fperp = Rrc + Rcx; // Friction for vorticity
 
diff --git a/src/hermes-2.cxx b/src/hermes-2.cxx
index f9c9874..0d59a50 100644
--- a/src/hermes-2.cxx
+++ b/src/hermes-2.cxx
@@ -257,6 +257,11 @@ int Hermes::init(bool restarting) {
 
   OPTION(optsc, evolve_plasma, true);
 
+  mesh->get(jyseps1_1, "jyseps1_1");
+  mesh->get(jyseps1_2, "jyseps1_2");
+  mesh->get(jyseps2_1, "jyseps2_1");
+  mesh->get(jyseps2_2, "jyseps2_2");
+
   electromagnetic = optsc["electromagnetic"]
                         .doc("Include vector potential psi in Ohm's law?")
                         .withDefault<bool>(true);
@@ -306,7 +311,11 @@ int Hermes::init(bool restarting) {
 
   poloidal_flows = optsc["poloidal_flows"]
                        .doc("Include ExB flows in X-Y plane")
-                       .withDefault(true);
+                       .withDefault<bool>(true);
+
+  revCurlb_B = optsc["revCurlb_B"]
+                    .doc("Reverse the direction of Curl(b/B)")
+                    .withDefault<bool>(false);
   
   OPTION(optsc, ion_velocity, true);
 
@@ -368,13 +377,19 @@ int Hermes::init(bool restarting) {
   OPTION(optsc, vort_dissipation, false);
   phi_dissipation = optsc["phi_dissipation"]
     .doc("Add a dissipation term to vorticity, depending on reconstruction of potential?")
-    .withDefault<bool>(false);
+    .withDefault(-1.0);
+  delp2_dissipation = optsc["delp2_dissipation"]
+    .doc("Add a delp2 dissipation term to Ne, Pi, Pe, Vort")
+    .withDefault(-1.0);
+
 
+  OPTION(optsc, nvi_hyper_z, -1.0);
+  OPTION(optsc, vepsi_hyper_z, -1.0);
   OPTION(optsc, ne_hyper_z, -1.0);
   OPTION(optsc, pe_hyper_z, -1.0);
 
-  OPTION(optsc, low_n_diffuse, false);
-  OPTION(optsc, low_n_diffuse_perp, false);
+  OPTION(optsc, low_n_diffuse, -1.0);
+  OPTION(optsc, low_n_diffuse_perp, -1.0);
 
   OPTION(optsc, resistivity_multiply, 1.0);
 
@@ -387,7 +402,10 @@ int Hermes::init(bool restarting) {
   OPTION(optsc, sheath_ydown, true);     // Apply sheath at ydown?
   OPTION(optsc, test_boundaries, false); // Test boundary conditions
   
-  nesheath_floor = optsc["nesheath_floor"].doc("Ne sheath lower limit").withDefault(1e-5);
+  OPTION(optsc, nelim_floor, 1e-5);
+  nesheath_floor = optsc["nesheath_floor"].doc("Ne sheath lower limit").withDefault(nelim_floor);
+  kappa_epar_floor = optsc["kappa_epar_floor"].doc("kappa_epar lower limit").withDefault(1.0);
+
 
   // Fix profiles in SOL
   OPTION(optsc, sol_fix_profiles, false);
@@ -396,6 +414,8 @@ int Hermes::init(bool restarting) {
     sol_te = FieldFactory::get()->parse("sol_te", &optsc);
   }
 
+  OPTION(optsc, slab_radial_buffers, false);
+
   radial_buffers = optsc["radial_buffers"]
     .doc("Turn on radial buffer regions?").withDefault<bool>(false);
   OPTION(optsc, radial_inner_width, 4);
@@ -509,19 +529,86 @@ int Hermes::init(bool restarting) {
   // Get switches from each variable section
   auto& optne = opt["Ne"];
   NeSource = optne["source"].doc("Source term in ddt(Ne)").withDefault(Field3D{0.0});
-  NeSource /= Omega_ci;
-  Sn = DC(NeSource);
 
   // Inflowing density carries momentum
   OPTION(optne, density_inflow, false);
 
   auto& optpe = opt["Pe"];
   PeSource = optpe["source"].withDefault(Field3D{0.0});
-  PeSource /= Omega_ci;
-  Spe = DC(PeSource);
 
   auto& optpi = opt["Pi"];
   PiSource = optpi["source"].withDefault(Field3D{0.0});
+
+  // radial buffer setup
+  if (slab_radial_buffers || radial_buffers) {
+    // Need to set the sources in the radial buffer regions to zero
+
+    if ((mesh->getGlobalXIndex(mesh->xstart) - mesh->xstart) < radial_inner_width) {
+      // inner boundary
+      int imax = mesh->xstart + radial_inner_width - 1 - 
+                 (mesh->getGlobalXIndex(mesh->xstart) - mesh->xstart);
+      if (imax > mesh->xend) {
+        imax = mesh->xend;
+      }
+      
+      int imin = mesh->xstart;
+      if (!mesh->firstX()) {
+        --imin;
+      }
+
+      int ncz = mesh->LocalNz;
+
+      for (int i = imin; i <= imax; i++) {
+        for (int j = mesh->ystart; j <= mesh->yend; ++j) {
+          for (int k = 0; k < ncz; ++k) {
+            NeSource(i, j, k) = 0.0;
+            PiSource(i, j, k) = 0.0;
+            PeSource(i, j ,k) = 0.0;
+          }
+        }
+      }
+    }
+    
+    int nguard = mesh->LocalNx - mesh->xend - 1;
+    // this is the same as mesh->xstart
+
+    if (mesh->GlobalNx - nguard - mesh->getGlobalXIndex(mesh->xend) <=
+        radial_outer_width) {
+      // Outer boundary
+      int imin =
+          mesh->GlobalNx - nguard - radial_outer_width - mesh->getGlobalXIndex(0);
+      if (imin < mesh->xstart) {
+        imin = mesh->xstart;
+      }
+
+      int imax = mesh->xend;
+      if (!mesh->lastX()) {
+        ++imax;
+      }
+
+      int ncz = mesh->LocalNz;
+      for (int i = imin; i <= imax; ++i) {
+        for (int j = mesh->ystart; j <= mesh->yend; ++j) {
+          for (int k = 0; k < ncz; ++k) {
+            NeSource(i, j, k) = 0.0;
+            PiSource(i, j, k) = 0.0;
+            PeSource(i, j ,k) = 0.0;
+          }
+        }
+      }
+    }
+  }
+
+  NeSource /= Omega_ci;
+  Sn = DC(NeSource);
+
+  // auto& optpe = opt["Pe"];
+  // PeSource = optpe["source"].withDefault(Field3D{0.0});
+  PeSource /= Omega_ci;
+  Spe = DC(PeSource);
+
+  // auto& optpi = opt["Pi"];
+  // PiSource = optpi["source"].withDefault(Field3D{0.0});
   PiSource /= Omega_ci;
   Spi = DC(PiSource);
 
@@ -533,7 +620,7 @@ int Hermes::init(bool restarting) {
         for (int y = mesh->ystart; y <= mesh->yend; y++) {
           Sn(x, y) = 0.0;
           Spe(x, y) = 0.0;
-	  Spi(x, y) = 0.0;
+	        Spi(x, y) = 0.0;
         }
       }
     }
@@ -591,7 +678,7 @@ int Hermes::init(bool restarting) {
   }
 
   if (verbose) {
-    SAVE_REPEAT(Ti);
+    // SAVE_REPEAT(Ti);
     if (electron_ion_transfer) {
       SAVE_REPEAT(Wi);
     }
@@ -695,12 +782,22 @@ int Hermes::init(bool restarting) {
   impurity_adas = optsc["impurity_adas"]
                       .doc("Use Atomic++ interface to ADAS")
                       .withDefault<bool>(false);
+  
+  impurity_split_y = optsc["impurity_split_y"]
+                    .doc("Enable a region dependant fixed fraction")
+                    .withDefault<bool>(false);
 
   if (impurity_adas) {
 
-    fimp = optsc["impurity_fraction"]
+    if (impurity_split_y) {
+      OPTION(optsc, fimp_id, 0.0);
+      OPTION(optsc, fimp_od, 0.0);
+      OPTION(optsc, fimp_us, 0.0);
+    } else {
+      fimp = optsc["impurity_fraction"]
                .doc("Fixed fraction ADAS impurity, multiple of electron density")
                .withDefault(0.0);
+    }
 
     string impurity_species =
         optsc["impurity_species"]
@@ -837,6 +934,10 @@ int Hermes::init(bool restarting) {
       }
     }
   }
+
+  if (revCurlb_B) {
+    Curlb_B = -1.0 * Curlb_B;
+  }
   
   if (Options::root()["mesh"]["paralleltransform"].as<std::string>() == "shifted") {
     // Check if the gridfile was created for "shiftedmetric" or for "identity" parallel
@@ -925,7 +1026,7 @@ int Hermes::init(bool restarting) {
       if (!restarting) {
         // Start by setting to the sheath current = 0 boundary value
 
-        Field3D Nelim = floor(Ne, 1e-5);
+        Field3D Nelim = floor(Ne, nelim_floor);
         Telim = floor(Pe / Nelim, 0.1 / Tnorm);
         Tilim = floor(Pi / Nelim, 0.1 / Tnorm);
 
@@ -999,6 +1100,8 @@ int Hermes::init(bool restarting) {
 int Hermes::rhs(BoutReal t) {
   Coordinates *coord = mesh->getCoordinates();
   
+  Ne = floor(Ne, 0.0);
+
   if (!evolve_plasma) {
     Ne = 0.0;
     Pe = 0.0;
@@ -1021,7 +1124,7 @@ int Hermes::rhs(BoutReal t) {
     f->clearParallelSlices(); // Make sure no parallel slices
   }
 
-  Field3D Nelim = floor(Ne, 1e-5);
+  Field3D Nelim = floor(Ne, nelim_floor);
 
   Te = Pe / Nelim;
   Vi = NVi / Nelim;
@@ -1043,8 +1146,8 @@ int Hermes::rhs(BoutReal t) {
     for (int j = mesh->ystart; j <= mesh->yend; j++) {
       for (int k = 0; k < mesh->LocalNz; k++) {
         BoutReal ne_bndry = 0.5 * (Ne(1, j, k) + Ne(2, j, k));
-        if (ne_bndry < 1e-5)
-          ne_bndry = 1e-5;
+        if (ne_bndry < nelim_floor)
+          ne_bndry = nelim_floor;
         BoutReal pe_bndry = 0.5 * (Pe(1, j, k) + Pe(2, j, k));
         BoutReal pi_bndry = 0.5 * (Pi(1, j, k) + Pi(2, j, k));
 
@@ -1070,8 +1173,8 @@ int Hermes::rhs(BoutReal t) {
     for (int j = mesh->ystart; j <= mesh->yend; j++) {
       for (int k = 0; k < mesh->LocalNz; k++) {
         BoutReal ne_bndry = 0.5 * (Ne(n - 1, j, k) + Ne(n - 2, j, k));
-        if (ne_bndry < 1e-5)
-          ne_bndry = 1e-5;
+        if (ne_bndry < nelim_floor)
+          ne_bndry = nelim_floor;
         BoutReal pe_bndry = 0.5 * (Pe(n - 1, j, k) + Pe(n - 2, j, k));
         BoutReal pi_bndry = 0.5 * (Pi(n - 1, j, k) + Pi(n - 2, j, k));
 
@@ -1363,7 +1466,7 @@ int Hermes::rhs(BoutReal t) {
         // Solve Helmholtz equation for psi
 
         // aparSolver->setCoefA(-Ne*0.5*mi_me*beta_e);
-        Field2D NDC = DC(Ne);
+        Field2D NDC = DC(Nelim);
         aparSolver->setCoefs(1.0, -NDC * 0.5 * mi_me * beta_e);
         // aparSolver->setCoefs(1.0, -Ne*0.5*mi_me*beta_e);
         if (split_n0_psi) {
@@ -2284,8 +2387,8 @@ int Hermes::rhs(BoutReal t) {
 
         // Get density at the boundary
         BoutReal ne_bndry = 0.5 * (Ne(r.ind, jy, jz) + Ne(r.ind, jy + 1, jz));
-        if (ne_bndry < 1e-5)
-          ne_bndry = 1e-5;
+        if (ne_bndry < nelim_floor)
+          ne_bndry = nelim_floor;
 
         NVi(r.ind, jy, jz) = 2 * ne_bndry * vi_bndry - NVi(r.ind, jy + 1, jz);
         Pe(r.ind, jy, jz) = 2 * ne_bndry * te_bndry - Pe(r.ind, jy + 1, jz);
@@ -2308,8 +2411,8 @@ int Hermes::rhs(BoutReal t) {
 
         // Get density at the boundary
         BoutReal ne_bndry = 0.5 * (Ne(r.ind, jy, jz) + Ne(r.ind, jy - 1, jz));
-        if (ne_bndry < 1e-5)
-          ne_bndry = 1e-5;
+        if (ne_bndry < nelim_floor)
+          ne_bndry = nelim_floor;
 
         NVi(r.ind, jy, jz) = 2 * ne_bndry * vi_bndry - NVi(r.ind, jy - 1, jz);
         Pe(r.ind, jy, jz) = 2 * ne_bndry * te_bndry - Pe(r.ind, jy - 1, jz);
@@ -2349,7 +2452,7 @@ int Hermes::rhs(BoutReal t) {
   }
 
   // Ensure that Nelim, Telim and Pelim are calculated in guard cells
-  Nelim = floor(Ne, 1e-5);
+  Nelim = floor(Ne, nelim_floor);
   Telim = floor(Te, 0.1 / Tnorm);
   Pelim = Telim * Nelim;
   Tilim = floor(Ti, 0.1 / Tnorm);
@@ -2574,7 +2677,7 @@ int Hermes::rhs(BoutReal t) {
     Field3D qipar = -kappa_ipar * Grad_par(Tifree);
 
     // Limit the maximum value of tau_i
-    tau_i = ceil(tau_i, 1e4);
+    tau_i = ceil(tau_i, 5e4);
 
     // Square of total heat flux, parallel and perpendicular
     // The first Pi term cancels the parallel part of the second term
@@ -2713,7 +2816,7 @@ int Hermes::rhs(BoutReal t) {
 
     Dn = (Telim + Tilim) / (tau_e * mi_me * SQ(coord->Bxy));
     ddt(Ne) += FV::Div_a_Laplace_perp(Dn, Ne);
-    ddt(Ne) += FV::Div_a_Laplace_perp(Ne / (tau_e * mi_me * SQ(coord->Bxy)),
+    ddt(Ne) += FV::Div_a_Laplace_perp(Nelim / (tau_e * mi_me * SQ(coord->Bxy)),
                                       Ti - 0.5 * Te);
   }
   if (anomalous_D > 0.0) {
@@ -2762,13 +2865,13 @@ int Hermes::rhs(BoutReal t) {
   
   ddt(Ne) += NeSource;
   
-  if (low_n_diffuse) {
+  if (low_n_diffuse > 0.0) {
     // Diffusion which kicks in at very low density, in order to
     // help prevent negative density regions
-    ddt(Ne) += FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Ne);
+    ddt(Ne) += low_n_diffuse * FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Ne);
   }
-  if (low_n_diffuse_perp) {
-    ddt(Ne) += Div_Perp_Lap_FV_Index(1e-4 / Nelim, Ne, ne_bndry_flux);
+  if (low_n_diffuse_perp > 0.0) {
+    ddt(Ne) += low_n_diffuse_perp * Div_Perp_Lap_FV_Index(1e-4 / Nelim, Ne, ne_bndry_flux);
   }
 
   if (ne_hyper_z > 0.) {
@@ -2882,11 +2985,19 @@ int Hermes::rhs(BoutReal t) {
       ddt(Vort) -= FV::Div_par(Vort, 0.0, max_speed);
     }
 
-    if (phi_dissipation) {
+    if (phi_dissipation > 0.0) {
       // Adds dissipation term like in other equations, but depending on gradient of potential
       // Note: Doesn't seem to need faster than sound speed
-      ddt(Vort) -= FV::Div_par(-phi, 0.0, sound_speed);
+      ddt(Vort) -= phi_dissipation * FV::Div_par(-phi, 0.0, sound_speed);
     }
+
+    if (delp2_dissipation > 0.0) {
+      ddt(Ne) += delp2_dissipation * Delp2(Ne);
+      ddt(Pe) += delp2_dissipation * Delp2(Pe);
+      ddt(Pi) += delp2_dissipation * Delp2(Pi);
+      ddt(Vort) += delp2_dissipation * Delp2(Vort);
+    }
+
   }
 
   ///////////////////////////////////////////////////////////
@@ -2962,6 +3073,11 @@ int Hermes::rhs(BoutReal t) {
       // Adds dissipation term like in other equations
       ddt(VePsi) -= FV::Div_par(Ve - Vi, 0.0, max_speed);
     }
+
+    if (vepsi_hyper_z > 0.0) {
+      ddt(VePsi) -= vepsi_hyper_z * SQ(SQ(coord->dz)) * D4DZ4(VePsi);
+    }
+    
   }
 
   ///////////////////////////////////////////////////////////
@@ -2990,6 +3106,10 @@ int Hermes::rhs(BoutReal t) {
     if (pe_par) {
       ddt(NVi) -= Grad_parP(Pe + Pi);
     }
+
+    if (nvi_hyper_z > 0.) {
+      ddt(NVi) -= nvi_hyper_z * SQ(SQ(coord->dz)) * D4DZ4(NVi);
+    }
  
     if (ion_viscosity_par) {
       TRACE("NVi:ion viscosity parallel");
@@ -3049,13 +3169,13 @@ int Hermes::rhs(BoutReal t) {
       ddt(NVi) -= hyperpar * FV::D4DY4_Index(Vi) / mi_me;
     }
 
-    if (low_n_diffuse) {
+    if (low_n_diffuse > 0.0) {
       // Diffusion which kicks in at very low density, in order to
       // help prevent negative density regions
-      ddt(NVi) += FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, NVi);
+      ddt(NVi) += low_n_diffuse * FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, NVi);
     }
-    if (low_n_diffuse_perp) {
-      ddt(NVi) += Div_Perp_Lap_FV_Index(1e-4 / Nelim, NVi, ne_bndry_flux);
+    if (low_n_diffuse_perp > 0.0) {
+      ddt(NVi) += low_n_diffuse_perp * Div_Perp_Lap_FV_Index(1e-4 / Nelim, NVi, ne_bndry_flux);
     }
   }
 
@@ -3223,6 +3343,15 @@ int Hermes::rhs(BoutReal t) {
     ddt(Pe) -= (2. / 3) * 0.71 * Jpar * Grad_parP(Te);
   }
 
+  if (low_n_diffuse > 0.0) {
+    // Diffusion which kicks in at very low density, in order to
+    // help prevent negative density regions
+    ddt(Pe) += low_n_diffuse * FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Pe);
+  }
+  if (low_n_diffuse_perp > 0.0) {
+    ddt(Pe) += low_n_diffuse_perp * Div_Perp_Lap_FV_Index(1e-4 / Nelim, Pe, ne_bndry_flux);
+  }
+
   if (pe_par_p_term) {
     // This term balances energetically the pressure term
     // in Ohm's law
@@ -3338,6 +3467,15 @@ int Hermes::rhs(BoutReal t) {
     ddt(Pi) = 0.0;
   }
 
+  if (low_n_diffuse > 0.0) {
+    // Diffusion which kicks in at very low density, in order to
+    // help prevent negative density regions
+    ddt(Pi) += low_n_diffuse * FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Pi);
+  }
+  if (low_n_diffuse_perp > 0.0) {
+    ddt(Pi) += low_n_diffuse_perp * Div_Perp_Lap_FV_Index(1e-4 / Nelim, Pi, ne_bndry_flux);
+  }
+
   // Parallel flow
   if (parallel_flow_p_term) {
     ddt(Pi) -= FV::Div_par(Pi, Vi, sound_speed);
@@ -3761,6 +3899,173 @@ int Hermes::rhs(BoutReal t) {
     }
   }
 
+  if (slab_radial_buffers) {
+    /// radial buffer regions for slab geometry sims - same on both inner and outer boundaries
+    // output.write(" rad_buff zone 2 - slab \n");
+
+    // Calculate flux sZ averages
+    Field2D PeDC = DC(Pe);
+    Field2D PiDC = DC(Pi);
+    Field2D NeDC = DC(Ne);
+    Field2D VortDC = DC(Vort);
+
+    // Store Field3D Delp2 terms outside of loops for optimisation
+    Field3D Delp2_Pe = Delp2(Pe);
+    Field3D Delp2_Pi = Delp2(Pi);
+    Field3D Delp2_Ne = Delp2(Ne);
+    Field3D Delp2_NVi = Delp2(NVi);
+
+    if ((mesh->getGlobalXIndex(mesh->xstart) - mesh->xstart) < radial_inner_width) {
+      // This processor contains points inside the inner radial boundary
+
+      int imax = mesh->xstart + radial_inner_width - 1 -
+                 (mesh->getGlobalXIndex(mesh->xstart) - mesh->xstart);
+      if (imax > mesh->xend) {
+        imax = mesh->xend;
+      }
+
+      int imin = mesh->xstart;
+      if (!mesh->firstX()) {
+        --imin; // Calculate in guard cells, for radial fluxes
+      }
+      int ncz = mesh->LocalNz;
+
+      for (int i = imin; i <= imax; ++i) {
+        // position inside the boundary (0 = on boundary, 0.5 = first cell)
+        BoutReal pos =
+            static_cast<BoutReal>(mesh->getGlobalXIndex(i) - mesh->xstart) + 0.5;
+
+        // Diffusion coefficient which increases towards the boundary
+        BoutReal D = radial_buffer_D * (1. - pos / radial_inner_width);
+
+        for (int j = mesh->ystart; j <= mesh->yend; ++j) {
+          BoutReal dx = coord->dx(i, j);
+          BoutReal dx_xp = coord->dx(i + 1, j);
+          BoutReal J = coord->J(i, j);
+          BoutReal J_xp = coord->J(i + 1, j);
+
+          // Calculate metric factors for radial fluxes
+          BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp);
+          BoutReal x_factor = rad_flux_factor / (J * dx);
+          BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp);
+
+          for (int k = 0; k < ncz; ++k) {
+            // Relax towards constant value on flux surface
+            // output.write("00 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+            ddt(Pe)(i, j, k) -= D * (Pe(i, j, k) - PeDC(i, j));
+            // output.write("11 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+            ddt(Pi)(i, j, k) -= D * (Pi(i, j, k) - PiDC(i, j));
+            ddt(Ne)(i, j, k) -= D * (Ne(i, j, k) - NeDC(i, j));
+            ddt(Vort)(i, j, k) -= D * (Vort(i, j, k) - VortDC(i, j));
+            // ddt(NVi)(i, j, k) -= D * NVi(i, j, k);
+            
+            // Radial fluxes
+            BoutReal f = D * (Ne(i + 1, j, k) - Ne(i, j, k));
+            ddt(Ne)(i, j, k) += f * x_factor;
+            ddt(Ne)(i + 1, j, k) -= f * xp_factor;
+
+            f = D * (Pe(i + 1, j, k) - Pe(i, j, k));
+            ddt(Pe)(i, j, k) += f * x_factor;
+            ddt(Pe)(i + 1, j, k) -= f * xp_factor;
+            // output.write("22 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+
+            f = D * (Pi(i + 1, j, k) - Pi(i, j, k));
+            ddt(Pi)(i, j, k) += f * x_factor;
+            ddt(Pi)(i + 1, j, k) -= f * xp_factor;
+
+            f = D * (Vort(i + 1, j, k) - Vort(i, j, k));
+            ddt(Vort)(i, j, k) += f * x_factor;
+            ddt(Vort)(i + 1, j, k) -= f * xp_factor;
+
+            // finite difference dissipation applied for simplicity as 
+            // conservation not strictly necessary here (hopefully)
+            ddt(Pe)(i, j, k) += Delp2_Pe(i, j, k) * D;
+            ddt(Pi)(i, j, k) += Delp2_Pi(i, j, k) * D;
+            ddt(Ne)(i, j, k) += Delp2_Ne(i, j, k) * D;
+            ddt(NVi)(i, j, k) += Delp2_NVi(i, j, k) * D;
+
+          }
+        }
+      }
+    }
+    // Number of points in outer guard cells
+    int nguard = mesh->LocalNx - mesh->xend - 1;
+
+    if (mesh->GlobalNx - nguard - mesh->getGlobalXIndex(mesh->xend) <=
+        radial_outer_width) {
+
+      // Outer boundary
+      int imin =
+          mesh->GlobalNx - nguard - radial_outer_width - mesh->getGlobalXIndex(0);
+      if (imin < mesh->xstart) {
+        imin = mesh->xstart;
+      }
+      int imax = mesh->xend;
+      if (!mesh->lastX()) {
+        ++imax; // Calculate in guard cells, for radial fluxes
+      }
+      int ncz = mesh->LocalNz;
+      for (int i = imin; i <= imax; ++i) {
+
+        // position inside the boundary
+        BoutReal pos =
+            static_cast<BoutReal>(mesh->GlobalNx - nguard - mesh->getGlobalXIndex(i)) -
+            0.5;
+
+        // Diffusion coefficient which increases towards the boundary
+        // BoutReal D = radial_buffer_D * (1. - pos / radial_outer_width);
+        BoutReal D = radial_buffer_D * SQ((1. - pos / radial_outer_width));
+
+        for (int j = mesh->ystart; j <= mesh->yend; ++j) {
+          BoutReal dx = coord->dx(i, j);
+          BoutReal dx_xp = coord->dx(i + 1, j);
+          BoutReal J = coord->J(i, j);
+          BoutReal J_xp = coord->J(i + 1, j);
+
+          // Calculate metric factors for radial fluxes
+          BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp);
+          BoutReal x_factor = rad_flux_factor / (J * dx);
+          BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp);
+
+          for (int k = 0; k < ncz; ++k) {
+            // Relax towards constant value on flux surface
+            // output.write("00 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+            ddt(Pe)(i, j, k) -= D * (Pe(i, j, k) - PeDC(i, j));
+            // output.write("11 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+            ddt(Pi)(i, j, k) -= D * (Pi(i, j, k) - PiDC(i, j));
+            ddt(Ne)(i, j, k) -= D * (Ne(i, j, k) - NeDC(i, j));
+            ddt(Vort)(i, j, k) -= D * (Vort(i, j, k) - VortDC(i, j));
+            // ddt(NVi)(i, j, k) -= D * NVi(i, j, k);
+            
+            // Radial fluxes
+            BoutReal f = D * (Ne(i + 1, j, k) - Ne(i, j, k));
+            ddt(Ne)(i, j, k) += f * x_factor;
+            ddt(Ne)(i + 1, j, k) -= f * xp_factor;
+
+            f = D * (Pe(i + 1, j, k) - Pe(i, j, k));
+            ddt(Pe)(i, j, k) += f * x_factor;
+            ddt(Pe)(i + 1, j, k) -= f * xp_factor;
+            // output.write("22 ddt Pe = %e\n", ddt(Pe)(3,4,2));
+
+            f = D * (Pi(i + 1, j, k) - Pi(i, j, k));
+            ddt(Pi)(i, j, k) += f * x_factor;
+            ddt(Pi)(i + 1, j, k) -= f * xp_factor;
+
+            f = D * (Vort(i + 1, j, k) - Vort(i, j, k));
+            ddt(Vort)(i, j, k) += f * x_factor;
+            ddt(Vort)(i + 1, j, k) -= f * xp_factor;
+
+            // finite difference delp2 dissipation
+            ddt(Pe)(i, j, k) += Delp2(Pe)(i, j, k) * D;
+            ddt(Pi)(i, j, k) += Delp2(Pi)(i, j, k) * D;
+            ddt(Ne)(i, j, k) += Delp2(Ne)(i, j, k) * D;
+            ddt(NVi)(i, j, k) += Delp2(NVi)(i, j, k) * D;
+          }
+        }
+      }
+    }
+  }
+
   ///////////////////////////////////////////////////////////
   // Neutral gas
   if (neutrals) {
@@ -3892,6 +4197,16 @@ int Hermes::rhs(BoutReal t) {
         BoutReal Ne_C = Ne[i],
           Ne_L = 0.5 * (Ne[i.ym()] + Ne[i]),
           Ne_R = 0.5 * (Ne[i] + Ne[i.yp()]);
+
+        if (impurity_split_y) {
+          if (mesh->getGlobalYIndexNoBoundaries(mesh->yend) < jyseps1_1) {
+            fimp = fimp_id;
+          } else if (mesh->getGlobalYIndexNoBoundaries(mesh->yend) > jyseps2_2) {
+            fimp = fimp_od;
+          } else {
+            fimp = fimp_us;
+          }
+        }
         
         BoutReal Rz_L = computeRadiatedPower(*impurity,
                                              Te_L * Tnorm,        // electron temperature [eV]
@@ -4025,6 +4340,7 @@ int Hermes::precon(BoutReal t, BoutReal gamma, BoutReal delta) {
   }
   if (thermal_conduction) {
     // Set the coefficient in front of Grad2_par2
+    kappa_epar = floor(kappa_epar, kappa_epar_floor);
     inv->setCoefB(-(2. / 3) * gamma * kappa_epar);
     Field3D dT = ddt(Pe);
     dT.applyBoundary("neumann");
diff --git a/src/mixed.cxx b/src/mixed.cxx
index 530f171..188bedd 100644
--- a/src/mixed.cxx
+++ b/src/mixed.cxx
@@ -64,6 +64,8 @@ NeutralMixed::NeutralMixed(Solver *solver, Mesh *UNUSED(mesh), Options &options)
   F = 0;
   Qi = 0;
   Rp = 0;
+
+  SAVE_REPEAT4(Riz, Rrc, Rcx, Rex);
 }
 
 void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
@@ -168,8 +170,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
   // Atomic processes
   TRACE("Atomic processes");
 
-  Field3D Riz, Rrc, Rcx;
-  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vn, S, F, Qi, Rp, Riz, Rrc, Rcx);
+  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vn, S, F, Qi, Rp, Riz, Rrc, Rcx, Rex);
 
   // Neutral cross-field diffusion coefficient
   BoutReal neutral_lmax = 0.1 / Lnorm;
@@ -254,7 +255,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
       + FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity
       ;
 
-  Fperp = Rcx + Rrc;
+  Fperp = Rcx + Rrc + Rex;
 
   /////////////////////////////////////////////////////
   // Neutral pressure
diff --git a/src/neutral-model.cxx b/src/neutral-model.cxx
index 8722b4f..9764213 100644
--- a/src/neutral-model.cxx
+++ b/src/neutral-model.cxx
@@ -48,7 +48,7 @@ void NeutralModel::neutral_rates(
     const Field3D &Vi, // Plasma quantities
     const Field3D &Nn, const Field3D &Tn, const Field3D &Vnpar, // Neutral gas
     Field3D &S, Field3D &F, Field3D &Qi, Field3D &R, // Transfer rates
-    Field3D &Riz, Field3D &Rrc, Field3D &Rcx) {      // Rates
+    Field3D &Riz, Field3D &Rrc, Field3D &Rcx, Field3D &Rex) {      // Rates
 
   // Allocate output fields
   S = 0.0;
@@ -59,6 +59,7 @@ void NeutralModel::neutral_rates(
   Riz = 0.0;
   Rrc = 0.0;
   Rcx = 0.0;
+  Rex = 0.0;
 
   Coordinates *coord = mesh->getCoordinates();
 
@@ -184,6 +185,27 @@ void NeutralModel::neutral_rates(
         BoutReal R_iz_R =
             Ne_R * Nn_R * hydrogen.ionisation(Te_R * Tnorm) * Nnorm / Fnorm;
 
+        // Excitation
+        if (h_excitation) {
+          /////////////////////////////////////////////////////////
+          // Electron-neutral excitation
+          // Note: Rates need checking
+          // Currently assuming that quantity calculated is in [eV m^3/s]
+
+          BoutReal R_ex_L = Ne_L * Nn_L *
+                            hydrogen.excitation(Te_L * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+          BoutReal R_ex_C = Ne_C * Nn_C *
+                            hydrogen.excitation(Te_C * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+          BoutReal R_ex_R = Ne_R * Nn_R *
+                            hydrogen.excitation(Te_R * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+
+          Rex(i, j, k) = (J_L * R_ex_L + 4. * J_C * R_ex_C + J_R * R_ex_R) /
+                          (6. * J_C);
+        }
+
         // Neutral sink, plasma source
         S(i, j, k) -=
             (J_L * R_iz_L + 4. * J_C * R_iz_C + J_R * R_iz_R) / (6. * J_C);
@@ -200,10 +222,15 @@ void NeutralModel::neutral_rates(
                        (6. * J_C);
 
         // Ionisation and electron excitation energy
+
+
         R(i, j, k) += (Eionize / Tnorm) *
                       (J_L * R_iz_L + 4. * J_C * R_iz_C + J_R * R_iz_R) /
                       (6. * J_C);
 
+        // Add the hydrogen excitation rates to the overall total
+        R(i, j, k) += Rex(i, j, k); 
+
         // Cell-averaged rate
         Riz(i, j, k) =
             (J_L * R_iz_L + 4. * J_C * R_iz_C + J_R * R_iz_R) / (6. * J_C);