diff --git a/GPUSimulators/cuda/SWE2D_FORCE.cu b/GPUSimulators/cuda/SWE2D_FORCE.cu
index 1db7b17..42bb364 100644
--- a/GPUSimulators/cuda/SWE2D_FORCE.cu
+++ b/GPUSimulators/cuda/SWE2D_FORCE.cu
@@ -28,27 +28,18 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
 __device__ 
 void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float F[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
                   const float g_, const float dx_, const float dt_) {
-                      
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
     //Compute fluxes along the x axis
-    {
-        int j=ty;
-        const int l = j + 1; //Skip ghost cells
-        for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
-            const int k = i;
-            
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
             // Q at interface from the right and left
-            const float3 Qp = make_float3(Q[0][l][k+1],
-                                          Q[1][l][k+1],
-                                          Q[2][l][k+1]);
-            const float3 Qm = make_float3(Q[0][l][k],
-                                          Q[1][l][k],
-                                          Q[2][l][k]);
+            const float3 Qp = make_float3(Q[0][j][i+1],
+                                          Q[1][j][i+1],
+                                          Q[2][j][i+1]);
+            const float3 Qm = make_float3(Q[0][j][i],
+                                          Q[1][j][i],
+                                          Q[2][j][i]);
                                        
             // Computed flux
             const float3 flux = FORCE_1D_flux(Qm, Qp, g_, dx_, dt_);
@@ -65,27 +56,19 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
   */
 __device__ 
 void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float G[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
                   const float g_, const float dy_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
     //Compute fluxes along the y axis
-    for (int j=ty; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
-        const int l = j;
-        {
-            int i=tx;
-            const int k = i + 1; //Skip ghost cells
-            
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
             // Q at interface from the right and left
             // Note that we swap hu and hv
-            const float3 Qp = make_float3(Q[0][l+1][k],
-                                          Q[2][l+1][k],
-                                          Q[1][l+1][k]);
-            const float3 Qm = make_float3(Q[0][l][k],
-                                          Q[2][l][k],
-                                          Q[1][l][k]);
+            const float3 Qp = make_float3(Q[0][j+1][i],
+                                          Q[2][j+1][i],
+                                          Q[1][j+1][i]);
+            const float3 Qm = make_float3(Q[0][j][i],
+                                          Q[2][j][i],
+                                          Q[1][j][i]);
 
             // Computed flux
             // Note that we swap back
@@ -120,7 +103,7 @@ __global__ void FORCEKernel(
     const unsigned int vars = 3;
     
     __shared__ float Q[3][h+2][w+2];
-    __shared__ float F[3][h+1][w+1];
+    __shared__ float F[3][h+2][w+2];
     
     //Read into shared memory
     readBlock<w, h, gc>( h0_ptr_,  h0_pitch_, Q[0], nx_+2, ny_+2);
diff --git a/GPUSimulators/cuda/SWE2D_HLL.cu b/GPUSimulators/cuda/SWE2D_HLL.cu
index c05552c..6e2c4ff 100644
--- a/GPUSimulators/cuda/SWE2D_HLL.cu
+++ b/GPUSimulators/cuda/SWE2D_HLL.cu
@@ -31,24 +31,20 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
 __device__
 void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float F[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
                   const float g_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    {
-        const int j=ty;
-        const int l = j + 1; //Skip ghost cells     
-        for (int i=tx; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) { 
-            const int k = i;
-            
-            const float3 Q_l  = make_float3(Q[0][l][k  ], Q[1][l][k  ], Q[2][l][k  ]);
-            const float3 Q_r  = make_float3(Q[0][l][k+1], Q[1][l][k+1], Q[2][l][k+1]);
-            
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
+            // Q at interface from the right and left
+            const float3 Q_r = make_float3(Q[0][j][i+1],
+                                           Q[1][j][i+1],
+                                           Q[2][j][i+1]);
+            const float3 Q_l = make_float3(Q[0][j][i],
+                                           Q[1][j][i],
+                                           Q[2][j][i]);
+                       
+            // Computed flux
             const float3 flux = HLL_flux(Q_l, Q_r, g_);
-            
-            //Write to shared memory
             F[0][j][i] = flux.x;
             F[1][j][i] = flux.y;
             F[2][j][i] = flux.z;
@@ -65,27 +61,23 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
   */
 __device__
 void computeFluxG(float Q[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float G[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
                   const float g_) {
-    //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;
-        {
-            const int i=tx;
-            const int k = i + 1; //Skip ghost cells
-            
-            //NOte that hu and hv are swapped ("transposing" the domain)!
-            const float3 Q_l = make_float3(Q[0][l  ][k], Q[2][l  ][k], Q[1][l  ][k]);
-            const float3 Q_r = make_float3(Q[0][l+1][k], Q[2][l+1][k], Q[1][l+1][k]);
+    //Compute fluxes along the y axis
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
+            // Q at interface from the right and left
+            // Note that we swap hu and hv
+            const float3 Q_r = make_float3(Q[0][j+1][i],
+                                           Q[2][j+1][i],
+                                           Q[1][j+1][i]);
+            const float3 Q_l = make_float3(Q[0][j][i],
+                                           Q[2][j][i],
+                                           Q[1][j][i]);
                                        
             // Computed flux
-            const float3 flux = HLL_flux(Q_l, Q_r, g_);
-            
-            //Write to shared memory
             //Note that we here swap hu and hv back to the original
+            const float3 flux = HLL_flux(Q_l, Q_r, g_);
             G[0][j][i] = flux.x;
             G[1][j][i] = flux.z;
             G[2][j][i] = flux.y;
@@ -128,7 +120,7 @@ __global__ void HLLKernel(
     
     //Shared memory variables
     __shared__ float Q[3][h+2][w+2];
-    __shared__ float F[3][h+1][w+1];
+    __shared__ float F[3][h+2][w+2];
     
     //Read into shared memory
     readBlock<w, h, gc>( h0_ptr_,  h0_pitch_, Q[0], nx_+2, ny_+2);
diff --git a/GPUSimulators/cuda/SWE2D_HLL2.cu b/GPUSimulators/cuda/SWE2D_HLL2.cu
index 1b55c86..363710b 100644
--- a/GPUSimulators/cuda/SWE2D_HLL2.cu
+++ b/GPUSimulators/cuda/SWE2D_HLL2.cu
@@ -33,33 +33,26 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
 __device__
 void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float Qx[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, const float dx_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    {
-        const 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;
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
             // 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 float3 Q_rl = make_float3(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]);
-            const float3 Q_rr = make_float3(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]);
+            const float3 Q_rl = make_float3(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] - 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] - 0.5f*Qx[2][j][i+1]);
+            const float3 Q_rr = make_float3(Q[0][j][i+1] + 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] + 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] + 0.5f*Qx[2][j][i+1]);
                                          
-            const float3 Q_ll = make_float3(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]);
-            const float3 Q_lr = make_float3(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]);
+            const float3 Q_ll = make_float3(Q[0][j][i] - 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] - 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] - 0.5f*Qx[2][j][i]);
+            const float3 Q_lr = make_float3(Q[0][j][i] + 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] + 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] + 0.5f*Qx[2][j][i]);
                                         
             //Evolve half a timestep (predictor step)
             const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
@@ -85,34 +78,27 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
   */
 __device__
 void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float Qy[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, 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
+    for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) {
             // 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 float3 Q_rl = make_float3(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]);
-            const float3 Q_rr = make_float3(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]);
+            const float3 Q_rl = make_float3(Q[0][j+1][i] - 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] - 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] - 0.5f*Qy[1][j+1][i]);
+            const float3 Q_rr = make_float3(Q[0][j+1][i] + 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] + 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] + 0.5f*Qy[1][j+1][i]);
                                         
-            const float3 Q_ll = make_float3(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]);
-            const float3 Q_lr = make_float3(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]);
+            const float3 Q_ll = make_float3(Q[0][j][i] - 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] - 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] - 0.5f*Qy[1][j][i]);
+            const float3 Q_lr = make_float3(Q[0][j][i] + 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] + 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] + 0.5f*Qy[1][j][i]);
                                      
             //Evolve half a timestep (predictor step)
             const float3 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
@@ -163,8 +149,8 @@ __global__ void HLL2Kernel(
             
     //Shared memory variables
     __shared__ float  Q[3][h+4][w+4];
-    __shared__ float Qx[3][h+2][w+2];
-    __shared__ float  F[3][h+1][w+1];
+    __shared__ float Qx[3][h+4][w+4];
+    __shared__ float  F[3][h+4][w+4];
     
     //Read into shared memory
     readBlock<w, h, gc>( h0_ptr_,  h0_pitch_, Q[0], nx_+2, ny_+2);
diff --git a/GPUSimulators/cuda/SWE2D_KP07.cu b/GPUSimulators/cuda/SWE2D_KP07.cu
index 3740776..4565cdf 100644
--- a/GPUSimulators/cuda/SWE2D_KP07.cu
+++ b/GPUSimulators/cuda/SWE2D_KP07.cu
@@ -96,6 +96,39 @@ void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
 
 
 
+__device__ void minmodSlopeX(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
+                  const float theta_) {
+    //Reconstruct slopes along x axis
+    for (int p=0; p<3; ++p) {
+        {
+            const int j = threadIdx.y+2;
+            for (int i=threadIdx.x+1; i<BLOCK_WIDTH+3; i+=BLOCK_WIDTH) {
+                Qx[p][j-2][i-1] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_);
+            }
+        }
+    }
+}
+
+
+/**
+  * Reconstructs a minmod slope for a whole block along the ordinate
+  */
+__device__ void minmodSlopeY(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
+                  const float theta_) {
+    //Reconstruct slopes along y axis
+    for (int p=0; p<3; ++p) {
+        const int i = threadIdx.x + 2;
+        for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+3; j+=BLOCK_HEIGHT) {
+            {
+                Qy[p][j-1][i-2] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_);
+            }
+        }
+    }
+}
+
+
 /**
   * This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
   */
@@ -159,8 +192,8 @@ __global__ void KP07Kernel(
     
     
     //Reconstruct slopes along x and axis
-    minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_);
-    minmodSlopeY<w, h, gc, vars>(Q, Qy, theta_);
+    minmodSlopeX(Q, Qx, theta_);
+    minmodSlopeY(Q, Qy, theta_);
     __syncthreads();
     
     
diff --git a/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
index 385e8c8..be50059 100644
--- a/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
+++ b/GPUSimulators/cuda/SWE2D_KP07_dimsplit.cu
@@ -31,77 +31,65 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 __device__
 void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qx[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float Qx[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, 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 float3 Q_rl = make_float3(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]);
-        const float3 Q_rr = make_float3(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]);
-                                     
-        const float3 Q_ll = make_float3(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]);
-        const float3 Q_lr = make_float3(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]);
-                                
-        //Evolve half a timestep (predictor step)
-        const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
-        const float3 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
+            // 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 float3 Q_rl = make_float3(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] - 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] - 0.5f*Qx[2][j][i+1]);
+            const float3 Q_rr = make_float3(Q[0][j][i+1] + 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] + 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] + 0.5f*Qx[2][j][i+1]);
 
-        // Compute flux based on prediction
-        const float3 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, g_);
-        
-        //Write to shared memory
-        F[0][j][i] = flux.x;
-        F[1][j][i] = flux.y;
-        F[2][j][i] = flux.z;
+            const float3 Q_ll = make_float3(Q[0][j][i] - 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] - 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] - 0.5f*Qx[2][j][i]);
+            const float3 Q_lr = make_float3(Q[0][j][i] + 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] + 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] + 0.5f*Qx[2][j][i]);
+                                    
+            //Evolve half a timestep (predictor step)
+            const float3 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
+            const float3 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, g_) - F_func(Q_lr, g_));
+
+            // Compute flux based on prediction
+            const float3 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, g_);
+            
+            //Write to shared memory
+            F[0][j][i] = flux.x;
+            F[1][j][i] = flux.y;
+            F[2][j][i] = flux.z;
+        }
     }    
 }
 
 __device__
 void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qy[3][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float Qy[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, 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
+    for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) {
             // 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 float3 Q_rl = make_float3(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]);
-            const float3 Q_rr = make_float3(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]);
+            const float3 Q_rl = make_float3(Q[0][j+1][i] - 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] - 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] - 0.5f*Qy[1][j+1][i]);
+            const float3 Q_rr = make_float3(Q[0][j+1][i] + 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] + 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] + 0.5f*Qy[1][j+1][i]);
                                         
-            const float3 Q_ll = make_float3(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]);
-            const float3 Q_lr = make_float3(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]);
+            const float3 Q_ll = make_float3(Q[0][j][i] - 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] - 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] - 0.5f*Qy[1][j][i]);
+            const float3 Q_lr = make_float3(Q[0][j][i] + 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] + 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] + 0.5f*Qy[1][j][i]);
                                      
             //Evolve half a timestep (predictor step)
             const float3 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, g_) - F_func(Q_rr, g_));
@@ -153,8 +141,8 @@ __global__ void KP07DimsplitKernel(
         
     //Shared memory variables
     __shared__ float  Q[3][h+4][w+4];
-    __shared__ float Qx[3][h+2][w+2];
-    __shared__ float  F[3][h+1][w+1];
+    __shared__ float Qx[3][h+4][w+4];
+    __shared__ float  F[3][h+4][w+4];
     
     
     
diff --git a/GPUSimulators/cuda/SWE2D_WAF.cu b/GPUSimulators/cuda/SWE2D_WAF.cu
index 7a56d8e..7c096a9 100644
--- a/GPUSimulators/cuda/SWE2D_WAF.cu
+++ b/GPUSimulators/cuda/SWE2D_WAF.cu
@@ -34,23 +34,15 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
 __device__
 void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float F[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float F[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, 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;
-            
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
             // Q at interface from the right and left
-            const float3 Ql2 = make_float3(Q[0][l][k-1], Q[1][l][k-1], Q[2][l][k-1]);
-            const float3 Ql1 = make_float3(Q[0][l][k  ], Q[1][l][k  ], Q[2][l][k  ]);
-            const float3 Qr1 = make_float3(Q[0][l][k+1], Q[1][l][k+1], Q[2][l][k+1]);
-            const float3 Qr2 = make_float3(Q[0][l][k+2], Q[1][l][k+2], Q[2][l][k+2]);
+            const float3 Ql2 = make_float3(Q[0][j][i-1], Q[1][j][i-1], Q[2][j][i-1]);
+            const float3 Ql1 = make_float3(Q[0][j][i  ], Q[1][j][i  ], Q[2][j][i  ]);
+            const float3 Qr1 = make_float3(Q[0][j][i+1], Q[1][j][i+1], Q[2][j][i+1]);
+            const float3 Qr2 = make_float3(Q[0][j][i+2], Q[1][j][i+2], Q[2][j][i+2]);
 
             // Computed flux
             const float3 flux = WAF_1D_flux(Ql2, Ql1, Qr1, Qr2, g_, dx_, dt_);
@@ -73,24 +65,16 @@ void computeFluxF(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
   */
 __device__
 void computeFluxG(float Q[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float G[3][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
+                  float G[3][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
                   const float g_, const float dy_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    //Compute fluxes along the y axis
-    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
+    for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) {
             // Q at interface from the right and left
             // Note that we swap hu and hv
-            const float3 Ql2 = make_float3(Q[0][l-1][k], Q[2][l-1][k], Q[1][l-1][k]);
-            const float3 Ql1 = make_float3(Q[0][l  ][k], Q[2][l  ][k], Q[1][l  ][k]);
-            const float3 Qr1 = make_float3(Q[0][l+1][k], Q[2][l+1][k], Q[1][l+1][k]);
-            const float3 Qr2 = make_float3(Q[0][l+2][k], Q[2][l+2][k], Q[1][l+2][k]);
+            const float3 Ql2 = make_float3(Q[0][j-1][i], Q[2][j-1][i], Q[1][j-1][i]);
+            const float3 Ql1 = make_float3(Q[0][j  ][i], Q[2][j  ][i], Q[1][j  ][i]);
+            const float3 Qr1 = make_float3(Q[0][j+1][i], Q[2][j+1][i], Q[1][j+1][i]);
+            const float3 Qr2 = make_float3(Q[0][j+2][i], Q[2][j+2][i], Q[1][j+2][i]);
             
             // Computed flux
             // Note that we swap back
@@ -138,7 +122,7 @@ __global__ void WAFKernel(
          
     //Shared memory variables
     __shared__ float Q[3][h+4][w+4];
-    __shared__ float F[3][h+1][w+1];
+    __shared__ float F[3][h+4][w+4];