diff --git a/div_ops.cxx b/div_ops.cxx
index f1ce897..e93970c 100644
--- a/div_ops.cxx
+++ b/div_ops.cxx
@@ -23,10 +23,10 @@
 #include "div_ops.hxx"
 
 #include <bout/mesh.hxx>
-#include <globals.hxx>
-#include <derivs.hxx>
-#include <output.hxx>
-#include <utils.hxx>
+#include <bout/globals.hxx>
+#include <bout/derivs.hxx>
+#include <bout/output.hxx>
+#include <bout/utils.hxx>
 #include <bout/assert.hxx>
 
 #include <cmath>
diff --git a/div_ops.hxx b/div_ops.hxx
index 0c08a73..3bc7a23 100644
--- a/div_ops.hxx
+++ b/div_ops.hxx
@@ -24,7 +24,7 @@
 #ifndef __DIV_OPS_H__
 #define __DIV_OPS_H__
 
-#include <field3d.hxx>
+#include <bout/field3d.hxx>
 
 /*!
  * Parallel diffusion (in y)
diff --git a/external/BOUT-dev b/external/BOUT-dev
index dfef9f6..fa83d8e 160000
--- a/external/BOUT-dev
+++ b/external/BOUT-dev
@@ -1 +1 @@
-Subproject commit dfef9f688e9500cf4ea65110f4d1809ee4ee3bff
+Subproject commit fa83d8eda4e8bd890ca8b12e7989de23855842bb
diff --git a/loadmetric.cxx b/loadmetric.cxx
index adcd56f..4b1fcf2 100644
--- a/loadmetric.cxx
+++ b/loadmetric.cxx
@@ -18,10 +18,10 @@
     along with SD1D.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <globals.hxx>
-#include <output.hxx>
-#include <utils.hxx>
-#include <field2d.hxx>
+#include <bout/globals.hxx>
+#include <bout/output.hxx>
+#include <bout/utils.hxx>
+#include <bout/field2d.hxx>
 #include <bout/mesh.hxx>
 
 #include "loadmetric.hxx"
diff --git a/loadmetric.hxx b/loadmetric.hxx
index 5fa3acc..260bf47 100644
--- a/loadmetric.hxx
+++ b/loadmetric.hxx
@@ -21,7 +21,7 @@
 #ifndef __LOADMETRIC_H__
 #define __LOADMETRIC_H__
 
-#include <bout_types.hxx>
+#include <bout/bout_types.hxx>
 
 void LoadMetric(BoutReal Lnorm, BoutReal Bnorm);
 
diff --git a/radiation.cxx b/radiation.cxx
index 9709c80..3d69c55 100644
--- a/radiation.cxx
+++ b/radiation.cxx
@@ -1,8 +1,8 @@
 
-#include <globals.hxx> // for mesh object
-#include <output.hxx>
-#include <boutexception.hxx>
-#include <utils.hxx>
+#include <bout/globals.hxx> // for mesh object
+#include <bout/output.hxx>
+#include <bout/boutexception.hxx>
+#include <bout/utils.hxx>
 
 #include "radiation.hxx"
 
diff --git a/radiation.hxx b/radiation.hxx
index cbbd1c2..36aeb4a 100644
--- a/radiation.hxx
+++ b/radiation.hxx
@@ -3,8 +3,8 @@
 #ifndef __RADIATION_H__
 #define __RADIATION_H__
 
-#include <field3d.hxx>
-#include <bout_types.hxx>
+#include <bout/field3d.hxx>
+#include <bout/bout_types.hxx>
 
 #include <vector>
 #include <cmath>
diff --git a/sd1d.cxx b/sd1d.cxx
index eb0fd70..a662553 100644
--- a/sd1d.cxx
+++ b/sd1d.cxx
@@ -40,8 +40,8 @@
 
 #include <bout/constants.hxx>
 #include <bout/physicsmodel.hxx>
-#include <derivs.hxx>
-#include <field_factory.hxx>
+#include <bout/derivs.hxx>
+#include <bout/field_factory.hxx>
 #include <bout/invert_pardiv.hxx>
 #include <bout/snb.hxx>
 #include <bout/fv_ops.hxx>
@@ -1850,7 +1850,7 @@ class SD1D : public PhysicsModel {
                         {"long_name", "neutral atom pressure source"}});
       }
 
-      if (mesh->lastY()) {
+      if (mesh->lastY(0)) {
         set_with_attrs(state["flux_ion"], flux_ion,
                        {{"units", "m^-2 s^-1"},
                         {"conversion", Nnorm * Cs0},