diff --git a/GPUSimulators/HLL2Euler.py b/GPUSimulators/EE2D_KP07_dimsplit.py
similarity index 76%
rename from GPUSimulators/HLL2Euler.py
rename to GPUSimulators/EE2D_KP07_dimsplit.py
index 9346d9c..c98136f 100644
--- a/GPUSimulators/HLL2Euler.py
+++ b/GPUSimulators/EE2D_KP07_dimsplit.py
@@ -20,8 +20,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 #Import packages we need
-import numpy as np
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
         
@@ -34,7 +34,7 @@ from GPUSimulators import Simulator, Common
 """
 Class that solves the SW equations using the Forward-Backward linear scheme
 """
-class HLL2Euler (Simulator.BaseSimulator):
+class EE2D_KP07_dimsplit (Simulator.BaseSimulator):
 
     """
     Initialization routine
@@ -47,31 +47,40 @@ class HLL2Euler (Simulator.BaseSimulator):
     dx: Grid cell spacing along x-axis
     dy: Grid cell spacing along y-axis
     dt: Size of each timestep 
-    g: Gravitational accelleration
+    gamma: Gas constant
+    p: pressure
     """
     def __init__(self, \
                  context, \
                  rho, rho_u, rho_v, E, \
                  nx, ny, \
                  dx, dy, dt, \
-                 g, \
-                 theta=1.8, \
+                 gamma, \
+                 theta=1.3, \
                  block_width=16, block_height=16):
                  
         # Call super constructor
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
-            block_width, block_height);
+            block_width, block_height)
+        self.gamma = np.float32(gamma)
             
         self.theta = np.float32(theta)
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("HLL2_kernel.cu", "HLL2Kernel", \
-                                        "iifffffiPiPiPiPiPiPi", \
-                                        BLOCK_WIDTH=self.local_size[0], \
-                                        BLOCK_HEIGHT=self.local_size[1])
+        #Get kernels
+        self.kernel = context.get_prepared_kernel("cuda/EE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \
+                                        "iifffffiPiPiPiPiPiPiPiPi", \
+                                        defines={
+                                            'BLOCK_WIDTH': self.block_size[0], 
+                                            'BLOCK_HEIGHT': self.block_size[1]
+                                        }, \
+                                        compile_args={
+                                            'no_extern_c': True,
+                                            'options': ["--use_fast_math"], 
+                                        }, \
+                                        jit_compile_args={})
         
         #Create data by uploading to device
         self.u0 = Common.ArakawaA2D(self.stream, \
@@ -90,10 +99,10 @@ class HLL2Euler (Simulator.BaseSimulator):
         return self.stepDimsplitXY(dt)
                 
     def stepDimsplitXY(self, dt):
-        self.kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
+        self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \
                 self.nx, self.ny, \
                 self.dx, self.dy, dt, \
-                self.g, \
+                self.gamma, \
                 self.theta, \
                 np.int32(0), \
                 self.u0[0].data.gpudata, self.u0[0].data.strides[0], \
@@ -102,16 +111,16 @@ class HLL2Euler (Simulator.BaseSimulator):
                 self.u0[3].data.gpudata, self.u0[3].data.strides[0], \
                 self.u1[0].data.gpudata, self.u1[0].data.strides[0], \
                 self.u1[1].data.gpudata, self.u1[1].data.strides[0], \
-                self.u1[2].data.gpudata, self.u1[2].data.strides[0]
+                self.u1[2].data.gpudata, self.u1[2].data.strides[0], \
                 self.u1[3].data.gpudata, self.u1[3].data.strides[0])
         self.u0, self.u1 = self.u1, self.u0
         self.t += dt
             
     def stepDimsplitYX(self, dt):
-        self.kernel.prepared_async_call(self.global_size, self.local_size, self.stream, \
+        self.kernel.prepared_async_call(self.grid_size, self.block_size, self.stream, \
                 self.nx, self.ny, \
                 self.dx, self.dy, dt, \
-                self.g, \
+                self.gamma, \
                 self.theta, \
                 np.int32(1), \
                 self.u0[0].data.gpudata, self.u0[0].data.strides[0], \
@@ -120,7 +129,7 @@ class HLL2Euler (Simulator.BaseSimulator):
                 self.u0[3].data.gpudata, self.u0[3].data.strides[0], \
                 self.u1[0].data.gpudata, self.u1[0].data.strides[0], \
                 self.u1[1].data.gpudata, self.u1[1].data.strides[0], \
-                self.u1[2].data.gpudata, self.u1[2].data.strides[0]
+                self.u1[2].data.gpudata, self.u1[2].data.strides[0], \
                 self.u1[3].data.gpudata, self.u1[3].data.strides[0])
         self.u0, self.u1 = self.u1, self.u0
         self.t += dt
diff --git a/GPUSimulators/FORCE.py b/GPUSimulators/FORCE.py
index e88224c..bd8074c 100644
--- a/GPUSimulators/FORCE.py
+++ b/GPUSimulators/FORCE.py
@@ -22,6 +22,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #Import packages we need
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
 
@@ -62,11 +63,11 @@ class FORCE (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
+        self.g = np.float32(g) 
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_FORCE.cu", "FORCEKernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_FORCE.cu", "FORCEKernel", \
                                         "iiffffPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/HLL.py b/GPUSimulators/HLL.py
index bc77ff8..11e5219 100644
--- a/GPUSimulators/HLL.py
+++ b/GPUSimulators/HLL.py
@@ -21,6 +21,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #Import packages we need
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
 
@@ -57,11 +58,11 @@ class HLL (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
+        self.g = np.float32(g) 
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_HLL.cu", "HLLKernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL.cu", "HLLKernel", \
                                         "iiffffPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/HLL2.py b/GPUSimulators/HLL2.py
index f4e141b..e45291d 100644
--- a/GPUSimulators/HLL2.py
+++ b/GPUSimulators/HLL2.py
@@ -20,8 +20,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 #Import packages we need
-import numpy as np
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
         
@@ -61,13 +61,13 @@ class HLL2 (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
+        self.g = np.float32(g) 
             
         self.theta = np.float32(theta)
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_HLL2.cu", "HLL2Kernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_HLL2.cu", "HLL2Kernel", \
                                         "iifffffiPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/KP07.py b/GPUSimulators/KP07.py
index 46ad8f4..75e4a11 100644
--- a/GPUSimulators/KP07.py
+++ b/GPUSimulators/KP07.py
@@ -25,8 +25,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 #Import packages we need
-import numpy as np
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
 
@@ -62,13 +62,12 @@ class KP07 (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
-            
+        self.g = np.float32(g)             
         self.theta = np.float32(theta)
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_KP07.cu", "KP07Kernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07.cu", "KP07Kernel", \
                                         "iifffffiPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/KP07_dimsplit.py b/GPUSimulators/KP07_dimsplit.py
index 92e7fa5..4fd080e 100644
--- a/GPUSimulators/KP07_dimsplit.py
+++ b/GPUSimulators/KP07_dimsplit.py
@@ -25,8 +25,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 #Import packages we need
-import numpy as np
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
 
@@ -62,13 +62,12 @@ class KP07_dimsplit (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
-            block_width, block_height);
-            
+            block_width, block_height)
+        self.g = np.float32(g)
         self.theta = np.float32(theta)
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_KP07_dimsplit.cu", "KP07DimsplitKernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \
                                         "iifffffiPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/LxF.py b/GPUSimulators/LxF.py
index 33ab080..c8aa554 100644
--- a/GPUSimulators/LxF.py
+++ b/GPUSimulators/LxF.py
@@ -22,6 +22,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #Import packages we need
 from GPUSimulators import Simulator, Common, CudaContext
+import numpy as np
 
 
 
@@ -58,11 +59,11 @@ class LxF (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
+        self.g = np.float32(g) 
 
         # Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_LxF.cu", "LxFKernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_LxF.cu", "LxFKernel", \
                                         "iiffffPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/Simulator.py b/GPUSimulators/Simulator.py
index 28fcd7a..3bf7abf 100644
--- a/GPUSimulators/Simulator.py
+++ b/GPUSimulators/Simulator.py
@@ -44,13 +44,11 @@ class BaseSimulator:
     dx: Grid cell spacing along x-axis (20 000 m)
     dy: Grid cell spacing along y-axis (20 000 m)
     dt: Size of each timestep (90 s)
-    g: Gravitational accelleration (9.81 m/s^2)
     """
     def __init__(self, \
                  context, \
                  nx, ny, \
                  dx, dy, dt, \
-                 g, \
                  block_width, block_height):
         #Get logger
         self.logger = logging.getLogger(__name__ + "." + self.__class__.__name__)
@@ -64,7 +62,6 @@ class BaseSimulator:
         self.dx = np.float32(dx)
         self.dy = np.float32(dy)
         self.dt = np.float32(dt)
-        self.g = np.float32(g) 
         
         #Handle autotuning block size
         if (self.context.autotuner):
diff --git a/GPUSimulators/WAF.py b/GPUSimulators/WAF.py
index c22c7f6..d9ae132 100644
--- a/GPUSimulators/WAF.py
+++ b/GPUSimulators/WAF.py
@@ -21,8 +21,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 #Import packages we need
-import numpy as np
 from GPUSimulators import Simulator, Common
+import numpy as np
 
 
 
@@ -57,11 +57,11 @@ class WAF (Simulator.BaseSimulator):
         super().__init__(context, \
             nx, ny, \
             dx, dy, dt, \
-            g, \
             block_width, block_height);
+        self.g = np.float32(g) 
 
         #Get kernels
-        self.kernel = context.get_prepared_kernel("cuda/SWE_WAF.cu", "WAFKernel", \
+        self.kernel = context.get_prepared_kernel("cuda/SWE2D_WAF.cu", "WAFKernel", \
                                         "iiffffiPiPiPiPiPiPi", \
                                         defines={
                                             'BLOCK_WIDTH': self.block_size[0], 
diff --git a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
new file mode 100644
index 0000000..c404ed2
--- /dev/null
+++ b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
@@ -0,0 +1,244 @@
+ /*
+This kernel implements the Central Upwind flux function to
+solve the Euler equations 
+
+Copyright (C) 2018  SINTEF Digital
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+
+#include "common.h"
+#include "EulerCommon.h"
+#include "limiters.h"
+
+
+__device__
+void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float Qx[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
+                  float F[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  const float gamma_, const float dx_, const float dt_) {
+    //Index of thread within block
+    const int tx = threadIdx.x;
+    const int ty = threadIdx.y;
+    
+    {
+        int j=ty;
+        const int l = j + 2; //Skip ghost cells
+        for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
+            const int k = i + 1;
+            // Reconstruct point values of Q at the left and right hand side 
+            // of the cell for both the left (i) and right (i+1) cell 
+            const float4 Q_rl = make_float4(Q[0][l][k+1] - 0.5f*Qx[0][j][i+1],
+                                            Q[1][l][k+1] - 0.5f*Qx[1][j][i+1],
+                                            Q[2][l][k+1] - 0.5f*Qx[2][j][i+1],
+                                            Q[4][l][k+1] - 0.5f*Qx[4][j][i+1]);
+            const float4 Q_rr = make_float4(Q[0][l][k+1] + 0.5f*Qx[0][j][i+1],
+                                            Q[1][l][k+1] + 0.5f*Qx[1][j][i+1],
+                                            Q[2][l][k+1] + 0.5f*Qx[2][j][i+1],
+                                            Q[4][l][k+1] + 0.5f*Qx[4][j][i+1]);
+                                         
+            const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qx[0][j][i],
+                                            Q[1][l][k] - 0.5f*Qx[1][j][i],
+                                            Q[2][l][k] - 0.5f*Qx[2][j][i],
+                                            Q[4][l][k] - 0.5f*Qx[4][j][i]);
+            const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qx[0][j][i],
+                                            Q[1][l][k] + 0.5f*Qx[1][j][i],
+                                            Q[2][l][k] + 0.5f*Qx[2][j][i],
+                                            Q[4][l][k] + 0.5f*Qx[4][j][i]);
+                                    
+            //Evolve half a timestep (predictor step)
+            const float4 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
+            const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
+
+            // Compute flux based on prediction
+            const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
+            
+            //Write to shared memory
+            F[0][j][i] = flux.x;
+            F[1][j][i] = flux.y;
+            F[2][j][i] = flux.z;
+            F[3][j][i] = flux.w;
+        }
+    }    
+}
+
+__device__
+void computeFluxG(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float Qy[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
+                  float G[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  const float gamma_, const float dy_, const float dt_) {
+    //Index of thread within block
+    const int tx = threadIdx.x;
+    const int ty = threadIdx.y;
+    
+    for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
+        const int l = j + 1;
+        {
+            int i=tx;
+            const int k = i + 2; //Skip ghost cells
+            // Reconstruct point values of Q at the left and right hand side 
+            // of the cell for both the left (i) and right (i+1) cell 
+            //NOte that hu and hv are swapped ("transposing" the domain)!
+            const float4 Q_rl = make_float4(Q[0][l+1][k] - 0.5f*Qy[0][j+1][i],
+                                            Q[2][l+1][k] - 0.5f*Qy[2][j+1][i],
+                                            Q[1][l+1][k] - 0.5f*Qy[1][j+1][i],
+                                            Q[3][l+1][k] - 0.5f*Qy[3][j+1][i]);
+            const float4 Q_rr = make_float4(Q[0][l+1][k] + 0.5f*Qy[0][j+1][i],
+                                            Q[2][l+1][k] + 0.5f*Qy[2][j+1][i],
+                                            Q[1][l+1][k] + 0.5f*Qy[1][j+1][i],
+                                            Q[3][l+1][k] + 0.5f*Qy[3][j+1][i]);
+                                       
+            const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qy[0][j][i],
+                                            Q[2][l][k] - 0.5f*Qy[2][j][i],
+                                            Q[1][l][k] - 0.5f*Qy[1][j][i],
+                                            Q[3][l][k] - 0.5f*Qy[3][j][i]);
+            const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qy[0][j][i],
+                                            Q[2][l][k] + 0.5f*Qy[2][j][i],
+                                            Q[1][l][k] + 0.5f*Qy[1][j][i],
+                                            Q[3][l][k] + 0.5f*Qy[3][j][i]);
+                                     
+            //Evolve half a timestep (predictor step)
+            const float4 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
+            const float4 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
+            
+            // Compute flux based on prediction
+            const float4 flux = make_float4(0.01, 0.01, 0.01, 0.01);//CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
+            
+            //Write to shared memory
+            //Note that we here swap hu and hv back to the original
+            G[0][j][i] = flux.x;
+            G[1][j][i] = flux.z;
+            G[2][j][i] = flux.y;
+            G[3][j][i] = flux.w;
+        }
+    }
+}
+
+
+
+
+/**
+  * This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
+  */
+extern "C" {
+__global__ void KP07DimsplitKernel(
+        int nx_, int ny_,
+        float dx_, float dy_, float dt_,
+        float gamma_,
+        
+        float theta_,
+        
+        int step_,
+        
+        //Input h^n
+        float* rho0_ptr_, int rho0_pitch_,
+        float* rho_u0_ptr_, int rho_u0_pitch_,
+        float* rho_v0_ptr_, int rho_v0_pitch_,
+        float* E0_ptr_, int E0_pitch_,
+        
+        //Output h^{n+1}
+        float* rho1_ptr_, int rho1_pitch_,
+        float* rho_u1_ptr_, int rho_u1_pitch_,
+        float* rho_v1_ptr_, int rho_v1_pitch_,
+        float* E1_ptr_, int E1_pitch_) {
+        
+    const unsigned int w = BLOCK_WIDTH;
+    const unsigned int h = BLOCK_HEIGHT;
+    const unsigned int gc = 2;
+        
+    //Shared memory variables
+    __shared__ float  Q[4][h+4][w+4];
+    __shared__ float Qx[4][h+2][w+2];
+    __shared__ float  F[4][h+1][w+1];
+    
+    
+    
+    //Read into shared memory
+    readBlock<w, h, gc>(  rho0_ptr_,   rho0_pitch_, Q[0], nx_+2, ny_+2);
+    readBlock<w, h, gc>(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_+2, ny_+2);
+    readBlock<w, h, gc>(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_+2, ny_+2);
+    readBlock<w, h, gc>(    E0_ptr_,     E0_pitch_, Q[3], nx_+2, ny_+2);
+    __syncthreads();
+    
+    
+    //Fix boundary conditions
+    noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
+    noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+    noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
+    noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
+    __syncthreads();
+    
+    
+    //Step 0 => evolve x first, then y
+    if (step_ == 0) {
+        //Compute fluxes along the x axis and evolve
+        minmodSlopeX(Q, Qx, theta_);
+        __syncthreads();
+        computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
+        __syncthreads();
+        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        __syncthreads();
+        
+        //Set boundary conditions
+        noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
+        noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
+        __syncthreads();
+        
+        //Compute fluxes along the y axis and evolve
+        minmodSlopeY(Q, Qx, theta_);
+        __syncthreads();
+        computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
+        __syncthreads();
+        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        __syncthreads();
+    }
+    //Step 1 => evolve y first, then x
+    else {
+        //Compute fluxes along the y axis and evolve
+        minmodSlopeY(Q, Qx, theta_);
+        __syncthreads();
+        computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
+        __syncthreads();
+        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        __syncthreads();
+        
+        //Set boundary conditions
+        noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
+        noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
+        __syncthreads();
+        
+        //Compute fluxes along the x axis and evolve
+        minmodSlopeX(Q, Qx, theta_);
+        __syncthreads();
+        computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
+        __syncthreads();
+        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        __syncthreads();
+    }
+    
+    
+    // Write to main memory for all internal cells
+    writeBlock<w, h, gc>(  rho1_ptr_,   rho1_pitch_, Q[0], nx_, ny_);
+    writeBlock<w, h, gc>(rho_u1_ptr_, rho_u1_pitch_, Q[1], nx_, ny_);
+    writeBlock<w, h, gc>(rho_v1_ptr_, rho_v1_pitch_, Q[2], nx_, ny_);
+    writeBlock<w, h, gc>(    E1_ptr_,     E1_pitch_, Q[3], nx_, ny_);
+}
+
+} // extern "C"
\ No newline at end of file
diff --git a/GPUSimulators/cuda/EulerCommon.h b/GPUSimulators/cuda/EulerCommon.h
index 47b1eb2..a748e0f 100644
--- a/GPUSimulators/cuda/EulerCommon.h
+++ b/GPUSimulators/cuda/EulerCommon.h
@@ -66,15 +66,15 @@ __device__ float4 CentralUpwindFlux(const float4 Qm, float4 Qp, const float gamm
     const float Pp = pressure(Qp, gamma);
     const float4 Fp = F_func(Qp, Pp);
     const float up = Qp.y / Qp.x;   // rho*u / rho
-    const float cp = sqrt(gamma*P*Qp.x); // sqrt(gamma*P/rho)
+    const float cp = sqrt(gamma*Pp*Qp.x); // sqrt(gamma*P/rho)
 
     const float Pm = pressure(Qm, gamma);
-    const float3 Fm = F_func(Qm, Pm);
+    const float4 Fm = F_func(Qm, Pm);
     const float um = Qm.y / Qm.x;   // rho*u / rho
-    const float cm = sqrt(gamma*P/Qm.x); // sqrt(g*h)
+    const float cm = sqrt(gamma*Pm/Qm.x); // sqrt(g*h)
     
     const float am = min(min(um-cm, up-cp), 0.0f); // largest negative wave speed
     const float ap = max(max(um+cm, up+cp), 0.0f); // largest positive wave speed
     
-    return ((ap*Fm - am*Fp) + ap*am*(Qp-Qm))/(ap-am);
+    return  ((ap*Fm - am*Fp) + ap*am*(Qp-Qm))/(ap-am);
 }
\ No newline at end of file
diff --git a/GPUSimulators/cuda/SWE_FORCE.cu b/GPUSimulators/cuda/SWE2D_FORCE.cu
similarity index 94%
rename from GPUSimulators/cuda/SWE_FORCE.cu
rename to GPUSimulators/cuda/SWE2D_FORCE.cu
index adab937..0762029 100644
--- a/GPUSimulators/cuda/SWE_FORCE.cu
+++ b/GPUSimulators/cuda/SWE2D_FORCE.cu
@@ -136,7 +136,10 @@ __global__ void FORCEKernel(
     //Compute flux along x, and evolve
     computeFluxF(Q, F, g_, dx_, dt_);
     __syncthreads();
-    evolveF1(Q, F, nx_, ny_, dx_, dt_);
+    
+    evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+    evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+    evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
     __syncthreads();
     
     //Set boundary conditions
@@ -148,7 +151,10 @@ __global__ void FORCEKernel(
     //Compute flux along y, and evolve
     computeFluxG(Q, F, g_, dy_, dt_);
     __syncthreads();
-    evolveG1(Q, F, nx_, ny_, dy_, dt_);
+    
+    evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+    evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+    evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
     __syncthreads();
     
     //Write to main memory
diff --git a/GPUSimulators/cuda/SWE_HLL.cu b/GPUSimulators/cuda/SWE2D_HLL.cu
similarity index 94%
rename from GPUSimulators/cuda/SWE_HLL.cu
rename to GPUSimulators/cuda/SWE2D_HLL.cu
index e34658c..8b0c2cb 100644
--- a/GPUSimulators/cuda/SWE_HLL.cu
+++ b/GPUSimulators/cuda/SWE2D_HLL.cu
@@ -144,7 +144,10 @@ __global__ void HLLKernel(
     //Compute F flux
     computeFluxF(Q, F, g_);
     __syncthreads();
-    evolveF1(Q, F, nx_, ny_, dx_, dt_);
+    
+    evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+    evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+    evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
     __syncthreads();
     
     //Set boundary conditions
@@ -156,7 +159,10 @@ __global__ void HLLKernel(
     //Compute G flux
     computeFluxG(Q, F, g_);
     __syncthreads();
-    evolveG1(Q, F, nx_, ny_, dy_, dt_);
+    
+    evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+    evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+    evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
     __syncthreads();
     
     // Write to main memory for all internal cells
diff --git a/GPUSimulators/cuda/SWE_HLL2.cu b/GPUSimulators/cuda/SWE2D_HLL2.cu
similarity index 93%
rename from GPUSimulators/cuda/SWE_HLL2.cu
rename to GPUSimulators/cuda/SWE2D_HLL2.cu
index dc63d66..40910cc 100644
--- a/GPUSimulators/cuda/SWE_HLL2.cu
+++ b/GPUSimulators/cuda/SWE2D_HLL2.cu
@@ -184,7 +184,10 @@ __global__ void HLL2Kernel(
         __syncthreads();
         computeFluxF(Q, Qx, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
         
         //Set boundary conditions
@@ -198,7 +201,10 @@ __global__ void HLL2Kernel(
         __syncthreads();
         computeFluxG(Q, Qx, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
     }
     //Step 1 => evolve y first, then x
@@ -208,7 +214,10 @@ __global__ void HLL2Kernel(
         __syncthreads();
         computeFluxG(Q, Qx, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
         
         //Set boundary conditions
@@ -222,7 +231,10 @@ __global__ void HLL2Kernel(
         __syncthreads();
         computeFluxF(Q, Qx, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
     }
     
diff --git a/GPUSimulators/cuda/SWE_KP07.cu b/GPUSimulators/cuda/SWE2D_KP07.cu
similarity index 100%
rename from GPUSimulators/cuda/SWE_KP07.cu
rename to GPUSimulators/cuda/SWE2D_KP07.cu
diff --git a/GPUSimulators/cuda/SWE_KP07_dimsplit.cu b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
similarity index 91%
rename from GPUSimulators/cuda/SWE_KP07_dimsplit.cu
rename to GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
index 74381d0..2957dcc 100644
--- a/GPUSimulators/cuda/SWE_KP07_dimsplit.cu
+++ b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
@@ -181,7 +181,9 @@ __global__ void KP07DimsplitKernel(
         __syncthreads();
         computeFluxF(Q, Qx, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
         
         //Set boundary conditions
@@ -190,12 +192,18 @@ __global__ void KP07DimsplitKernel(
         noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
         __syncthreads();
         
+        
+        
         //Compute fluxes along the y axis and evolve
         minmodSlopeY(Q, Qx, theta_);
         __syncthreads();
+        
         computeFluxG(Q, Qx, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
     }
     //Step 1 => evolve y first, then x
@@ -205,7 +213,10 @@ __global__ void KP07DimsplitKernel(
         __syncthreads();
         computeFluxG(Q, Qx, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
         
         //Set boundary conditions
@@ -219,15 +230,17 @@ __global__ void KP07DimsplitKernel(
         __syncthreads();
         computeFluxF(Q, Qx, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
     }
     
     
     // Write to main memory for all internal cells
-    writeBlock<w, h, 2>( h1_ptr_,  h1_pitch_, Q[0], nx_, ny_);
-    writeBlock<w, h, 2>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_);
-    writeBlock<w, h, 2>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_);
+    writeBlock<w, h, gc>( h1_ptr_,  h1_pitch_, Q[0], nx_, ny_);
+    writeBlock<w, h, gc>(hu1_ptr_, hu1_pitch_, Q[1], nx_, ny_);
+    writeBlock<w, h, gc>(hv1_ptr_, hv1_pitch_, Q[2], nx_, ny_);
 }
 
 } // extern "C"
\ No newline at end of file
diff --git a/GPUSimulators/cuda/SWE_LxF.cu b/GPUSimulators/cuda/SWE2D_LxF.cu
similarity index 100%
rename from GPUSimulators/cuda/SWE_LxF.cu
rename to GPUSimulators/cuda/SWE2D_LxF.cu
diff --git a/GPUSimulators/cuda/SWE_WAF.cu b/GPUSimulators/cuda/SWE2D_WAF.cu
similarity index 90%
rename from GPUSimulators/cuda/SWE_WAF.cu
rename to GPUSimulators/cuda/SWE2D_WAF.cu
index 62a82e2..1b80e4d 100644
--- a/GPUSimulators/cuda/SWE_WAF.cu
+++ b/GPUSimulators/cuda/SWE2D_WAF.cu
@@ -161,7 +161,10 @@ __global__ void WAFKernel(
         //Compute fluxes along the x axis and evolve
         computeFluxF(Q, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
         
         //Fix boundary conditions
@@ -173,7 +176,10 @@ __global__ void WAFKernel(
         //Compute fluxes along the y axis and evolve
         computeFluxG(Q, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
     }
     //Step 1 => evolve y first, then x
@@ -181,7 +187,10 @@ __global__ void WAFKernel(
         //Compute fluxes along the y axis and evolve
         computeFluxG(Q, F, g_, dy_, dt_);
         __syncthreads();
-        evolveG2(Q, F, nx_, ny_, dy_, dt_);
+        
+        evolveG<w, h, gc>(Q[0], F[0], dy_, dt_);
+        evolveG<w, h, gc>(Q[1], F[1], dy_, dt_);
+        evolveG<w, h, gc>(Q[2], F[2], dy_, dt_);
         __syncthreads();
         
         //Fix boundary conditions
@@ -193,7 +202,10 @@ __global__ void WAFKernel(
         //Compute fluxes along the x axis and evolve
         computeFluxF(Q, F, g_, dx_, dt_);
         __syncthreads();
-        evolveF2(Q, F, nx_, ny_, dx_, dt_);
+        
+        evolveF<w, h, gc>(Q[0], F[0], dx_, dt_);
+        evolveF<w, h, gc>(Q[1], F[1], dx_, dt_);
+        evolveF<w, h, gc>(Q[2], F[2], dx_, dt_);
         __syncthreads();
     }
 
diff --git a/GPUSimulators/cuda/common.h b/GPUSimulators/cuda/common.h
index faf961d..9fe08ff 100644
--- a/GPUSimulators/cuda/common.h
+++ b/GPUSimulators/cuda/common.h
@@ -44,6 +44,28 @@ inline __device__ float3 operator+(const float3 a, const float3 b) {
     return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
 }
 
+/**
+  * Float4 operators 
+  */
+inline __device__ float4 operator*(const float a, const float4 b) {
+    return make_float4(a*b.x, a*b.y, a*b.z, a*b.w);
+}
+
+inline __device__ float4 operator/(const float4 a, const float b) {
+    return make_float4(a.x/b, a.y/b, a.z/b, a.w/b);
+}
+
+inline __device__ float4 operator-(const float4 a, const float4 b) {
+    return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
+}
+
+inline __device__ float4 operator+(const float4 a, const float4 b) {
+    return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
+}
+
+
+
+
 inline __device__ __host__ float clamp(const float f, const float a, const float b) {
     return fmaxf(a, fminf(f, b));
 }
@@ -223,62 +245,23 @@ __device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2
 
 
 
-
-
-
-/**
-  * Evolves the solution in time along the x axis (dimensional splitting)
-  */
-__device__ void evolveF1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-              float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-              const int nx_, const int ny_,
+template<int block_width, int block_height, int ghost_cells>
+__device__ void evolveF(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells],
+              float F[block_height+1][block_width+1],
               const float dx_, const float dt_) {
     //Index of thread within block
     const int tx = threadIdx.x;
     const int ty = threadIdx.y;
+        
+    const int i = tx + ghost_cells; //Skip local ghost cells
+    const int j = ty + ghost_cells;
     
     //Index of cell within domain
-    const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1
-    const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1;
+    //const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells; //Skip global ghost cells, i.e., +1
+    //const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells;
+    //if (ti > ghost_cells-1 && ti < nx_+ghost_cells && tj > ghost_cells-1 && tj < ny_+ghost_cells) {
+    Q[j][i] = Q[j][i] + (F[ty][tx] - F[ty][tx+1]) * dt_ / dx_;
     
-    if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) {
-        const int i = tx + 1; //Skip local ghost cells, i.e., +1
-        const int j = ty + 1;
-        
-        Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_;
-        Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_;
-        Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_;
-    }
-}
-
-
-
-
-
-
-/**
-  * Evolves the solution in time along the x axis (dimensional splitting)
-  */
-__device__ void evolveF2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-              float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-              const int nx_, const int ny_,
-              const float dx_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    //Index of cell within domain
-    const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
-    const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
-    
-    if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) {
-        const int i = tx + 2; //Skip local ghost cells, i.e., +1
-        const int j = ty + 2;
-        
-        Q[0][j][i] = Q[0][j][i] + (F[0][ty][tx] - F[0][ty][tx+1]) * dt_ / dx_;
-        Q[1][j][i] = Q[1][j][i] + (F[1][ty][tx] - F[1][ty][tx+1]) * dt_ / dx_;
-        Q[2][j][i] = Q[2][j][i] + (F[2][ty][tx] - F[2][ty][tx+1]) * dt_ / dx_;
-    }
 }
 
 
@@ -289,57 +272,18 @@ __device__ void evolveF2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
 /**
   * Evolves the solution in time along the y axis (dimensional splitting)
   */
-__device__ void evolveG1(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-              float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-              const int nx_, const int ny_,
+template<int block_width, int block_height, int ghost_cells>
+__device__ void evolveG(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells],
+              float G[block_height+1][block_width+1],
               const float dy_, const float dt_) {
     //Index of thread within block
     const int tx = threadIdx.x;
     const int ty = threadIdx.y;
     
-    //Index of cell within domain
-    const int ti = blockDim.x*blockIdx.x + threadIdx.x + 1; //Skip global ghost cells, i.e., +1
-    const int tj = blockDim.y*blockIdx.y + threadIdx.y + 1;
+    const int i = tx + ghost_cells; //Skip local ghost cells, i.e., +1
+    const int j = ty + ghost_cells;
     
-    if (ti > 0 && ti < nx_+1 && tj > 0 && tj < ny_+1) {
-        const int i = tx + 1; //Skip local ghost cells, i.e., +1
-        const int j = ty + 1;
-        
-        Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_;
-        Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_;
-        Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_;
-    }
-}
-
-
-
-
-
-
-
-/**
-  * Evolves the solution in time along the y axis (dimensional splitting)
-  */
-__device__ void evolveG2(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-              float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-              const int nx_, const int ny_,
-              const float dy_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    //Index of cell within domain
-    const int ti = blockDim.x*blockIdx.x + threadIdx.x + 2; //Skip global ghost cells, i.e., +2
-    const int tj = blockDim.y*blockIdx.y + threadIdx.y + 2;
-    
-    if (ti > 1 && ti < nx_+2 && tj > 1 && tj < ny_+2) {
-        const int i = tx + 2; //Skip local ghost cells, i.e., +2
-        const int j = ty + 2;
-        
-        Q[0][j][i] = Q[0][j][i] + (G[0][ty][tx] - G[0][ty+1][tx]) * dt_ / dy_;
-        Q[1][j][i] = Q[1][j][i] + (G[1][ty][tx] - G[1][ty+1][tx]) * dt_ / dy_;
-        Q[2][j][i] = Q[2][j][i] + (G[2][ty][tx] - G[2][ty+1][tx]) * dt_ / dy_;
-    }
+    Q[j][i] = Q[j][i] + (G[ty][tx] - G[ty+1][tx]) * dt_ / dy_;
 }