diff --git a/analysis/python/jupyter/notebook_example.ipynb b/analysis/python/jupyter/notebook_example.ipynb
index b968762..feb0a05 100644
--- a/analysis/python/jupyter/notebook_example.ipynb
+++ b/analysis/python/jupyter/notebook_example.ipynb
@@ -7,7 +7,8 @@
    "outputs": [],
    "source": [
     "import astar.data as ad\n",
-    "import astar.visual as vis"
+    "import astar.visual as vis\n",
+    "import numpy as np"
    ]
   },
   {
@@ -23,7 +24,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "meshdir  = \"/scratch/data/mvaisala/iotest/\""
+    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\""
    ]
   },
   {
@@ -32,10 +33,26 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "imesh = 0\n",
+    "imesh = 10000\n",
     "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Auxuliary variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -51,7 +68,11 @@
    },
    "outputs": [],
    "source": [
-    "vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = False, fname = 'uux')"
+    "vis.slices.plot_3(mesh, uu_tot, title = r'$\\|u\\|$')\n",
+    "vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$')\n",
+    "vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$')\n",
+    "vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$')\n",
+    "vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\\ln_\\rho$')"
    ]
   },
   {
@@ -86,6 +107,13 @@
    "outputs": [],
    "source": []
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/analysis/python/samples/readtest.py b/analysis/python/samples/readtest.py
index 81f7512..a2263c7 100644
--- a/analysis/python/samples/readtest.py
+++ b/analysis/python/samples/readtest.py
@@ -45,7 +45,7 @@ print("sys.argv", sys.argv)
 #meshdir = "/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/"
 #meshdir  = "/tiara/ara/data/mvaisala/tmp/astaroth-code/astaroth_2.0/build/"
 #meshdir = "/tiara/ara/data/mvaisala/asth_testbed_double/"
-meshdir  = "/scratch/data/mvaisala/iotest/"
+meshdir  = "/scratch/data/mvaisala/forcingtest/"
 
 if "xtopbound" in sys.argv: 
     for i in range(0, 171):
@@ -169,7 +169,7 @@ if '1d' in sys.argv:
 if 'sl' in sys.argv:
     #maxfiles = 200002
     #stride = 10000
-    maxfiles = 1001
+    maxfiles = 5001
     stride = 100
     for i in range(0, maxfiles, stride):
         mesh = ad.read.Mesh(i, fdir=meshdir) 
diff --git a/src/core/kernels/kernels.cuh b/src/core/kernels/kernels.cuh
index b4d748c..97d0803 100644
--- a/src/core/kernels/kernels.cuh
+++ b/src/core/kernels/kernels.cuh
@@ -545,10 +545,12 @@ normalized(const AcReal3& vec)
 
 // Sinusoidal forcing
 // https://arxiv.org/pdf/1704.04676.pdf
+// NOTE: This method of forcing is depracated. However, it will remain in here
+// until a corresponding scheme exists in the new code. 
 __constant__ AcReal3 forcing_vec;
 __constant__ AcReal forcing_phi;
 static __device__ __forceinline__ AcReal3
-forcing(const int i, const int j, const int k)
+DEPRECATED_forcing(const int i, const int j, const int k)
 {
 #define DOMAIN_SIZE_X (DCONST_INT(AC_nx) * DCONST_REAL(AC_dsx))
 #define DOMAIN_SIZE_Y (DCONST_INT(AC_ny) * DCONST_REAL(AC_dsy))