diff --git a/src/mesh/mesh_vars.f90 b/src/mesh/mesh_vars.f90
index cf565b9..377a2e6 100644
--- a/src/mesh/mesh_vars.f90
+++ b/src/mesh/mesh_vars.f90
@@ -83,8 +83,10 @@ MODULE MOD_Mesh_Vars
   REAL      :: pvar(NVAR)           ! Primitive Variables
   REAL      :: cvar(NVAR)           ! Conservative Variables
   REAL      :: cvar_stage(NVAR)     ! Conservative Variables at initial RK stage (needed for 2nd order RK time stepping)
-  REAL      :: u_x(NVAR)            ! x-gradient of pvar
-  REAL      :: u_y(NVAR)            ! y-gradient of pvar
+  REAL      :: u_x(NVAR)            ! x-gradient of pvar (limited for convection)
+  REAL      :: u_y(NVAR)            ! y-gradient of pvar (limited for convection)
+  REAL      :: u_x_unlim(NVAR)      ! unlimited x-gradient (for viscous terms)
+  REAL      :: u_y_unlim(NVAR)      ! unlimited y-gradient (for viscous terms)
   REAL      :: u_t(NVAR)            ! t-gradient of pvar
   REAL      :: source(NVAR)         ! source term
   REAL      :: dt                   ! Cell timestep
diff --git a/src/navierstokes/flux/fluxcalculation.f90 b/src/navierstokes/flux/fluxcalculation.f90
index d562ed4..da1ff5c 100644
--- a/src/navierstokes/flux/fluxcalculation.f90
+++ b/src/navierstokes/flux/fluxcalculation.f90
@@ -83,8 +83,13 @@ SUBROUTINE FluxCalculation()
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! Extract left and right gradient
   state_mean(:)  = 0.5*(aSide%connection%pvar(:) + aSide%pvar(:) )
-  grad_Uxmean(:) = 0.5*(aSide%Elem%u_x(:)+aSide%connection%Elem%u_x(:))
-  grad_Uymean(:) = 0.5*(aSide%Elem%u_y(:)+aSide%connection%Elem%u_y(:))
+  IF (ASSOCIATED(aSide%connection%BC)) THEN ! the side is a boundary side
+      grad_Uxmean(:) = aSide%Elem%u_x_unlim(:)
+      grad_Uymean(:) = aSide%Elem%u_y_unlim(:)
+  ELSE
+    grad_Uxmean(:) = 0.5*(aSide%Elem%u_x_unlim(:)+aSide%connection%Elem%u_x_unlim(:))
+    grad_Uymean(:) = 0.5*(aSide%Elem%u_y_unlim(:)+aSide%connection%Elem%u_y_unlim(:))
+  END IF
   BaryBary(:)    = aSide%BaryBaryVec(:)/aSide%BaryBaryDist
   correction(:)  = (grad_Uxmean(:)*BaryBary(1) + grad_Uymean(:)*BaryBary(2) -       &
                    (aSide%connection%Elem%pvar(:)-aSide%Elem%pvar(:))/aSide%BaryBaryDist)
diff --git a/src/reconstruction/reconstruction.f90 b/src/reconstruction/reconstruction.f90
index bc0c9c5..200c16f 100644
--- a/src/reconstruction/reconstruction.f90
+++ b/src/reconstruction/reconstruction.f90
@@ -56,6 +56,8 @@ SUBROUTINE SpatialReconstruction(time)
     aSide => aElem%firstSide
     aElem%u_x(:) = 0
     aElem%u_y(:) = 0
+    aElem%u_x_unlim(:) = 0
+    aElem%u_y_unlim(:) = 0
     aElem%u_t(:) = 0
     DO WHILE(ASSOCIATED(aSide))
       aSide%pVar(:) = aElem%pVar(:)
@@ -74,6 +76,8 @@ SUBROUTINE SpatialReconstruction(time)
     aElem => Elems(iElem)%Elem
     aElem%u_x(:) = 0
     aElem%u_y(:) = 0
+    aElem%u_x_unlim(:) = 0
+    aElem%u_y_unlim(:) = 0
     aElem%u_t(:) = 0
   END DO
   !$omp end parallel do
@@ -125,6 +129,15 @@ SUBROUTINE SpatialReconstruction(time)
 !  END DO
 !  !$omp end parallel do
 
+! Save unlimited gradients before limiting
+  !$omp parallel do private(aElem)
+  DO iElem = 1, nElems
+    aElem => Elems(iElem)%Elem
+    aElem%u_x_unlim = aElem%u_x
+    aElem%u_y_unlim = aElem%u_y
+  END DO
+  !$omp end parallel do
+
 !-----------------------------------------------------------------------------------------------------------------------------------
   !$omp parallel do private(aElem,aSide,dx,dy)
 
@@ -151,4 +164,3 @@ SUBROUTINE SpatialReconstruction(time)
 END SUBROUTINE SpatialReconstruction
 
 END MODULE MOD_Reconstruction
-