diff --git a/src/GksGpu/FlowStateData/FlowStateData.cuh b/src/GksGpu/FlowStateData/FlowStateData.cuh
index b304edfd62df7b643f268ca16130ac073090d9e2..1ae9f6225b3c4ca6b535b64c67ec6816bb2a1c1c 100644
--- a/src/GksGpu/FlowStateData/FlowStateData.cuh
+++ b/src/GksGpu/FlowStateData/FlowStateData.cuh
@@ -127,7 +127,7 @@ __host__ __device__ inline ConservedVariables toConservedVariables( const Primit
                              ,prim.U * prim.rho
                              ,prim.V * prim.rho
                              ,prim.W * prim.rho
-                             ,( K + two ) / ( four * prim.lambda ) * prim.rho + c1o2 * prim.rho * ( prim.U * prim.U + prim.V * prim.V + prim.W * prim.W )
+                             ,( K + three ) / ( four * prim.lambda ) * prim.rho + c1o2 * prim.rho * ( prim.U * prim.U + prim.V * prim.V + prim.W * prim.W )
     #ifdef USE_PASSIVE_SCALAR
                              ,prim.S * prim.rho
     #endif
@@ -144,7 +144,7 @@ __host__ __device__ inline PrimitiveVariables toPrimitiveVariables( const Conser
 						     ,cons.rhoU / cons.rho
 						     ,cons.rhoV / cons.rho
 						     ,cons.rhoW / cons.rho
-						     ,( K + two ) * cons.rho / ( four * ( cons.rhoE - c1o2 * ( cons.rhoU * cons.rhoU + cons.rhoV * cons.rhoV + cons.rhoW * cons.rhoW ) / cons.rho ) )
+						     ,( K + three ) * cons.rho / ( four * ( cons.rhoE - c1o2 * ( cons.rhoU * cons.rhoU + cons.rhoV * cons.rhoV + cons.rhoW * cons.rhoW ) / cons.rho ) )
     #ifdef USE_PASSIVE_SCALAR
                              ,cons.rhoS / cons.rho
     #endif
diff --git a/src/GksGpu/FluxComputation/Reconstruction.cuh b/src/GksGpu/FluxComputation/Reconstruction.cuh
index 999ee2e45cb5cbfffef71301894cbf852f050cfb..ef852525510eaed3a22664a6c5e711565d3e9541 100644
--- a/src/GksGpu/FluxComputation/Reconstruction.cuh
+++ b/src/GksGpu/FluxComputation/Reconstruction.cuh
@@ -100,6 +100,81 @@ __host__ __device__ inline void computeGradN( const Parameters& parameters,
 #endif // USE_PASSIVE_SCALAR
 }
 
+__host__ __device__ inline void computeGradT( const DataBaseStruct& dataBase,
+                                              const Parameters& parameters,
+                                              const uint posCellIndexT[2],
+                                              const uint negCellIndexT[2],
+                                              const PrimitiveVariables& facePrim,
+                                              ConservedVariables& gradN )
+{
+    ConservedVariables cons;
+
+    //////////////////////////////////////////////////////////////////////////
+    {
+        readCellData(posCellIndexT[0], dataBase, cons);
+
+        gradN.rho  += c1o2 * cons.rho;
+        gradN.rhoU += c1o2 * cons.rhoU;
+        gradN.rhoV += c1o2 * cons.rhoV;
+        gradN.rhoW += c1o2 * cons.rhoW;
+        gradN.rhoE += c1o2 * cons.rhoE;
+    #ifdef USE_PASSIVE_SCALAR
+        gradN.rhoS += c1o2 * cons.rhoS;
+    #endif // USE_PASSIVE_SCALAR
+    }
+    {
+        readCellData(posCellIndexT[1], dataBase, cons);
+
+        gradN.rho  += c1o2 * cons.rho;
+        gradN.rhoU += c1o2 * cons.rhoU;
+        gradN.rhoV += c1o2 * cons.rhoV;
+        gradN.rhoW += c1o2 * cons.rhoW;
+        gradN.rhoE += c1o2 * cons.rhoE;
+    #ifdef USE_PASSIVE_SCALAR
+        gradN.rhoS += c1o2 * cons.rhoS;
+    #endif // USE_PASSIVE_SCALAR
+    }
+    //////////////////////////////////////////////////////////////////////////
+    {
+        readCellData(negCellIndexT[0], dataBase, cons);
+
+        gradN.rho  -= c1o2 * cons.rho;
+        gradN.rhoU -= c1o2 * cons.rhoU;
+        gradN.rhoV -= c1o2 * cons.rhoV;
+        gradN.rhoW -= c1o2 * cons.rhoW;
+        gradN.rhoE -= c1o2 * cons.rhoE;
+    #ifdef USE_PASSIVE_SCALAR
+        gradN.rhoS -= c1o2 * cons.rhoS;
+    #endif // USE_PASSIVE_SCALAR
+    }
+    {
+        readCellData(negCellIndexT[1], dataBase, cons);
+
+        gradN.rho  -= c1o2 * cons.rho;
+        gradN.rhoU -= c1o2 * cons.rhoU;
+        gradN.rhoV -= c1o2 * cons.rhoV;
+        gradN.rhoW -= c1o2 * cons.rhoW;
+        gradN.rhoE -= c1o2 * cons.rhoE;
+    #ifdef USE_PASSIVE_SCALAR
+        gradN.rhoS -= c1o2 * cons.rhoS;
+    #endif // USE_PASSIVE_SCALAR
+    }
+    //////////////////////////////////////////////////////////////////////////
+    {
+        gradN.rho  /= two * parameters.dx * facePrim.rho;
+        gradN.rhoU /= two * parameters.dx * facePrim.rho;
+        gradN.rhoV /= two * parameters.dx * facePrim.rho;
+        gradN.rhoW /= two * parameters.dx * facePrim.rho;
+        gradN.rhoE /= two * parameters.dx * facePrim.rho;
+    #ifdef USE_PASSIVE_SCALAR
+        gradN.rhoS /= two * parameters.dx * facePrim.rho;
+    #endif // USE_PASSIVE_SCALAR
+    }
+    //////////////////////////////////////////////////////////////////////////
+}
+
+
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -141,69 +216,18 @@ __host__ __device__ inline void reconstructFiniteDifferences( const uint faceInd
         if( direction == 'y' ) getCellIndicesTZ(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT1, negCellIndexT1);
         if( direction == 'z' ) getCellIndicesTX(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT1, negCellIndexT1);
 
-        ConservedVariables cons;
-
-        //////////////////////////////////////////////////////////////////////////
-        {
-            readCellData(posCellIndexT1[0], dataBase, cons);
-
-            gradN.rho  += c1o2 * cons.rho;
-            gradN.rhoU += c1o2 * cons.rhoU;
-            gradN.rhoV += c1o2 * cons.rhoV;
-            gradN.rhoW += c1o2 * cons.rhoW;
-            gradN.rhoE += c1o2 * cons.rhoE;
-        #ifdef USE_PASSIVE_SCALAR
-            gradN.rhoS += c1o2 * cons.rhoS;
-        #endif // USE_PASSIVE_SCALAR
-        }
-        {
-            readCellData(posCellIndexT1[1], dataBase, cons);
-
-            gradN.rho  += c1o2 * cons.rho;
-            gradN.rhoU += c1o2 * cons.rhoU;
-            gradN.rhoV += c1o2 * cons.rhoV;
-            gradN.rhoW += c1o2 * cons.rhoW;
-            gradN.rhoE += c1o2 * cons.rhoE;
-        #ifdef USE_PASSIVE_SCALAR
-            gradN.rhoS += c1o2 * cons.rhoS;
-        #endif // USE_PASSIVE_SCALAR
-        }
-        //////////////////////////////////////////////////////////////////////////
-        {
-            readCellData(negCellIndexT1[0], dataBase, cons);
-
-            gradN.rho  -= c1o2 * cons.rho;
-            gradN.rhoU -= c1o2 * cons.rhoU;
-            gradN.rhoV -= c1o2 * cons.rhoV;
-            gradN.rhoW -= c1o2 * cons.rhoW;
-            gradN.rhoE -= c1o2 * cons.rhoE;
-        #ifdef USE_PASSIVE_SCALAR
-            gradN.rhoS -= c1o2 * cons.rhoS;
-        #endif // USE_PASSIVE_SCALAR
-        }
-        {
-            readCellData(negCellIndexT1[1], dataBase, cons);
-
-            gradN.rho  -= c1o2 * cons.rho;
-            gradN.rhoU -= c1o2 * cons.rhoU;
-            gradN.rhoV -= c1o2 * cons.rhoV;
-            gradN.rhoW -= c1o2 * cons.rhoW;
-            gradN.rhoE -= c1o2 * cons.rhoE;
-        #ifdef USE_PASSIVE_SCALAR
-            gradN.rhoS -= c1o2 * cons.rhoS;
-        #endif // USE_PASSIVE_SCALAR
-        }
-        //////////////////////////////////////////////////////////////////////////
-        {
-            gradN.rho  /= two * parameters.dx * facePrim.rho;
-            gradN.rhoU /= two * parameters.dx * facePrim.rho;
-            gradN.rhoV /= two * parameters.dx * facePrim.rho;
-            gradN.rhoW /= two * parameters.dx * facePrim.rho;
-            gradN.rhoE /= two * parameters.dx * facePrim.rho;
-        #ifdef USE_PASSIVE_SCALAR
-            gradN.rhoS /= two * parameters.dx * facePrim.rho;
-        #endif // USE_PASSIVE_SCALAR
-        }
+        computeGradT( dataBase, parameters, posCellIndexT1, negCellIndexT1, facePrim, gradT1 );
+    }
+
+    {
+        uint posCellIndexT2[2];
+        uint negCellIndexT2[2];
+    
+        if( direction == 'x' ) getCellIndicesTZ(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT2, negCellIndexT2);
+        if( direction == 'y' ) getCellIndicesTX(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT2, negCellIndexT2);
+        if( direction == 'z' ) getCellIndicesTY(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT2, negCellIndexT2);
+
+        computeGradT( dataBase, parameters, posCellIndexT2, negCellIndexT2, facePrim, gradT2 );
     }
 }