diff --git a/WENO.F90 b/WENO.F90
index 36608d072e00ce30bb9df3c403faf27939a6aef5..0b4bbb03b52bf977e316be6aced36717a9174d58 100644
--- a/WENO.F90
+++ b/WENO.F90
@@ -30,6 +30,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
   real, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
 
   real flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
+
   real generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
   real generalized_states_frozen(nvars, -2:3)
   real generalized_fluxes_frozen(nvars, ndim, -2:3)
@@ -38,7 +40,6 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
   real lambda(nvars)
   integer i, j, k, l, m
   real delta_xi, delta_eta, delta_zeta
-  real flux(nvars)
   real characteristic_fluxes_pos(nvars, -2:3)
   real characteristic_fluxes_neg(nvars, -2:3)
   real metrics_frozen(ndim, ndim)
@@ -50,6 +51,7 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
   delta_eta = 1.0
   delta_zeta = 1.0
 
+  !$loopy begin tagged: init_flux
   do k=1,nz
     do j=1,ny
       do i=1,nx
@@ -61,11 +63,16 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
       end do
     end do
   end do
+  !$loopy end tagged: init_flux
 
+  !$loopy begin tagged: to_generalized
   call convert_to_generalized(nvars, ndim, nx, ny, nz, &
     fluxes, metrics, metric_jacobians, generalized_fluxes)
+  !$loopy end tagged: to_generalized
 
   ! for the next three loops, note carefully that the inner loop indices have a different range
+
+  !$loopy begin tagged: flux_x_compute
   do k=1,nz
     do j=1,ny
       do i=0,nx
@@ -96,18 +103,27 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
         call weno_flux(nvars, &
                         generalized_fluxes(:, 1, i-2:i+3, j, k), &
                         characteristic_fluxes_pos, characteristic_fluxes_neg, &
-                        combined_frozen_metrics(1), R, flux)
+                        combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k))
+      end do
+    end do
+  end do
+  !$loopy end tagged: flux_x_compute
 
+  !$loopy begin tagged: flux_x_diff
+  do k=1,nz
+    do j=1,ny
+      do i=1,nx
         do v=1,nvars
           flux_derivatives_generalized(v, 1, i, j, k) = &
-            flux_derivatives_generalized(v, 1, i, j, k) + flux(v) / delta_xi
-          flux_derivatives_generalized(v, 1, i+1, j, k) = &
-            flux_derivatives_generalized(v, 1, i+1, j, k) - flux(v) / delta_xi
+            flux_derivatives_generalized(v, 1, i, j, k) &
+            + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) / delta_xi
         end do
       end do
     end do
   end do
+  !$loopy end tagged: flux_x_diff
 
+  !$loopy begin tagged: flux_y_compute
   do k=1,nz
     do i=1,nx
       do j=0,ny
@@ -138,18 +154,28 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
         call weno_flux(nvars, &
                         generalized_fluxes(:, 2, i, j-2:j+3, k), &
                         characteristic_fluxes_pos, characteristic_fluxes_neg, &
-                        combined_frozen_metrics(2), R, flux)
+                        combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k))
 
+      end do
+    end do
+  end do
+  !$loopy end tagged: flux_y_compute
+
+  !$loopy begin tagged: flux_y_diff
+  do k=1,nz
+    do j=1,ny
+      do i=1,nx
         do v=1,nvars
           flux_derivatives_generalized(v, 2, i, j, k) = &
-            flux_derivatives_generalized(v, 2, i, j, k) + flux(v) / delta_eta
-          flux_derivatives_generalized(v, 2, i, j+1, k) = &
-            flux_derivatives_generalized(v, 2, i, j+1, k) - flux(v) / delta_eta
+            flux_derivatives_generalized(v, 2, i, j, k) &
+            + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) / delta_eta
         end do
       end do
     end do
   end do
+  !$loopy end tagged: flux_y_diff
 
+  !$loopy begin tagged: flux_z_compute
   do j=1,ny
     do i=1,nx
       do k=0,nz
@@ -180,23 +206,33 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
         call weno_flux(nvars, &
                         generalized_fluxes(:, 3, i, j, k-2:k+3), &
                         characteristic_fluxes_pos, characteristic_fluxes_neg, &
-                        combined_frozen_metrics(3), R, flux)
+                        combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k))
+      end do
+    end do
+  end do
+  !$loopy end tagged: flux_z_compute
 
+  !$loopy begin tagged: flux_z_diff
+  do k=1,nz
+    do j=1,ny
+      do i=1,nx
         do v=1,nvars
           flux_derivatives_generalized(v, 3, i, j, k) = &
-            flux_derivatives_generalized(v, 3, i, j, k) + flux(v) / delta_zeta
-          flux_derivatives_generalized(v, 3, i, j, k+1) = &
-            flux_derivatives_generalized(v, 3, i, j, k+1) - flux(v) / delta_zeta
+            flux_derivatives_generalized(v, 3, i, j, k) + &
+            (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) / delta_zeta
         end do
       end do
     end do
   end do
+  !$loopy end tagged: flux_z_diff
 
+  !$loopy begin tagged: from_generalized
   call convert_from_generalized(nvars, ndim, nx, ny, nz, &
                                 flux_derivatives_generalized, &
                                 metrics, &
                                 metric_jacobians, &
                                 flux_derivatives)
+  !$loopy end tagged: from_generalized
 end subroutine
 
 subroutine pointwise_eigenvalues(nvars, d, states, lambda_pointwise)