diff --git a/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.cpp b/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.cpp
index da82b2f5ea8cbf62a77413ab50a124603e1f9ce6..33d2dded3aba94e98b97daa0482fa60cb9317d8d 100755
--- a/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.cpp
+++ b/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.cpp
@@ -326,6 +326,7 @@ namespace Nektar
                       const int n,
                       const NekDouble time,
                       int &cnt,
+                      Array<OneD, Array<OneD, NekDouble> > &Fwd,
                       Array<OneD, Array<OneD, NekDouble> > &inarray)
     {
         std::string varName;
@@ -335,48 +336,48 @@ namespace Nektar
         {
             if(boost::iequals(userDefStr,"Wall"))
             {
-                WallBC(n, cnt, inarray);
+                WallBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"WallViscous") ||
                     boost::iequals(userDefStr,"WallAdiabatic"))
             {
                 // Wall Boundary Condition
-                WallViscousBC(n, cnt, inarray);
+                WallViscousBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"Symmetry"))
             {
                 // Symmetric Boundary Condition
-                SymmetryBC(n, cnt, inarray);
+                SymmetryBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"RiemannInvariant"))
             {
                 // Riemann invariant characteristic Boundary Condition
-                RiemannInvariantBC(n, cnt, inarray);
+                RiemannInvariantBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"PressureOutflowNonReflective"))
             {
                 // Pressure outflow non-reflective Boundary Condition
-                PressureOutflowNonReflectiveBC(n, cnt, inarray);
+                PressureOutflowNonReflectiveBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"PressureOutflow"))
             {
                 // Pressure outflow Boundary Condition
-                PressureOutflowBC(n, cnt, inarray);
+                PressureOutflowBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"PressureOutflowFile"))
             {
                 // Pressure outflow Boundary Condition from file 
-                PressureOutflowFileBC(n, cnt, inarray);
+                PressureOutflowFileBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"PressureInflowFile"))
             {
                 // Pressure inflow Boundary Condition from file
-                PressureInflowFileBC(n, cnt, inarray);
+                PressureInflowFileBC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"ExtrapOrder0"))
             {
                 // Extrapolation of the data at the boundaries
-                ExtrapOrder0BC(n, cnt, inarray);
+                ExtrapOrder0BC(n, cnt, Fwd, inarray);
             }
             else if(boost::iequals(userDefStr,"TimeDependent"))
             {
@@ -401,6 +402,7 @@ namespace Nektar
     void CompressibleFlowSystem::WallBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i;
@@ -410,14 +412,6 @@ namespace Nektar
         const Array<OneD, const int> &traceBndMap
             = m_fields[0]->GetTraceBndMap();
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Adjust the physical values of the trace to take
         // user defined boundaries into account
         int e, id1, id2, nBCEdgePts, eMax;
@@ -498,6 +492,7 @@ namespace Nektar
     void CompressibleFlowSystem::WallViscousBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i;
@@ -507,14 +502,6 @@ namespace Nektar
         const Array<OneD, const int> &traceBndMap
             = m_fields[0]->GetTraceBndMap();
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Take into account that for PDE based shock capturing, eps = 0 at the
         // wall. Adjust the physical values of the trace to take user defined
         // boundaries into account
@@ -573,6 +560,7 @@ namespace Nektar
     void CompressibleFlowSystem::SymmetryBC(
         int                                      bcRegion,
         int                                      cnt,
+        Array<OneD, Array<OneD, NekDouble> >    &Fwd,
         Array<OneD, Array<OneD, NekDouble> >    &physarray)
     {
         int i;
@@ -582,14 +570,6 @@ namespace Nektar
         const Array<OneD, const int> &traceBndMap
             = m_fields[0]->GetTraceBndMap();
 
-        // Get physical values of the forward trace (from exp to phys)
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Take into account that for PDE based shock capturing, eps = 0 at the
         // wall.
         int e, id1, id2, nBCEdgePts, eMax;
@@ -669,6 +649,7 @@ namespace Nektar
     void CompressibleFlowSystem::RiemannInvariantBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -706,14 +687,6 @@ namespace Nektar
             Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
         }
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Computing the normal velocity for characteristics coming
         // from inside the computational domain
         Array<OneD, NekDouble > Vn (nTracePts, 0.0);
@@ -926,6 +899,7 @@ namespace Nektar
     void CompressibleFlowSystem::PressureOutflowNonReflectiveBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -962,14 +936,6 @@ namespace Nektar
             Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
         }
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Computing the normal velocity for characteristics coming
         // from inside the computational domain
         Array<OneD, NekDouble > Vn (nTracePts, 0.0);
@@ -1086,6 +1052,7 @@ namespace Nektar
     void CompressibleFlowSystem::PressureOutflowBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -1122,14 +1089,6 @@ namespace Nektar
             Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
         }
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Computing the normal velocity for characteristics coming
         // from inside the computational domain
         Array<OneD, NekDouble > Vn (nTracePts, 0.0);
@@ -1247,6 +1206,7 @@ namespace Nektar
     void CompressibleFlowSystem::PressureOutflowFileBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -1283,15 +1243,6 @@ namespace Nektar
             Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
         }
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Computing the normal velocity for characteristics coming
         // from inside the computational domain
         Array<OneD, NekDouble > Vn (nTracePts, 0.0);
@@ -1411,6 +1362,7 @@ namespace Nektar
     void CompressibleFlowSystem::PressureInflowFileBC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -1447,15 +1399,6 @@ namespace Nektar
             Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
         }
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         // Computing the normal velocity for characteristics coming
         // from inside the computational domain
         Array<OneD, NekDouble > Vn (nTracePts, 0.0);
@@ -1574,6 +1517,7 @@ namespace Nektar
     void CompressibleFlowSystem::ExtrapOrder0BC(
         int                                   bcRegion,
         int                                   cnt,
+        Array<OneD, Array<OneD, NekDouble> > &Fwd,
         Array<OneD, Array<OneD, NekDouble> > &physarray)
     {
         int i, j;
@@ -1586,14 +1530,6 @@ namespace Nektar
         const Array<OneD, const int> &traceBndMap
             = m_fields[0]->GetTraceBndMap();
 
-        // Get physical values of the forward trace
-        Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
-        for (i = 0; i < nVariables; ++i)
-        {
-            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
-            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
-        }
-
         int eMax;
 
         eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
diff --git a/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h b/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h
index a08268e99a406fb45a7f212f711a56ae4f63216b..6db0e95ffb135d0599ebab5e38f17af2adff835d 100755
--- a/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h
+++ b/solvers/CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h
@@ -164,43 +164,53 @@ namespace Nektar
                          const int n,
                          const NekDouble time,
                          int &cnt,
+                         Array<OneD, Array<OneD, NekDouble> > &Fwd,
                          Array<OneD, Array<OneD, NekDouble> > &inarray);
         
         void WallBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void WallViscousBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void SymmetryBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void RiemannInvariantBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void PressureOutflowNonReflectiveBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void PressureOutflowBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void PressureOutflowFileBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void PressureInflowFileBC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void ExtrapOrder0BC(
             int                                                 bcRegion,
             int                                                 cnt,
+            Array<OneD, Array<OneD, NekDouble> > &Fwd,
             Array<OneD, Array<OneD, NekDouble> >               &physarray);
         void GetVelocityVector(
             const Array<OneD,       Array<OneD, NekDouble> > &physfield,
diff --git a/solvers/CompressibleFlowSolver/EquationSystems/EulerADCFE.cpp b/solvers/CompressibleFlowSolver/EquationSystems/EulerADCFE.cpp
index 09c730beaa0880afdde15825ac0f4e40fdf27f63..8deec8469b031525552d17e215a901ec7f4a9cde 100755
--- a/solvers/CompressibleFlowSolver/EquationSystems/EulerADCFE.cpp
+++ b/solvers/CompressibleFlowSolver/EquationSystems/EulerADCFE.cpp
@@ -275,6 +275,15 @@ namespace Nektar
         std::string varName;
         int cnt        = 0;
 
+        int nTracePts = GetTraceTotPoints();
+
+        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
+        for (int i = 0; i < nvariables; ++i)
+        {
+            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
+            m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
+        }
+
         // loop over Boundary Regions
         for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
         {
@@ -289,12 +298,12 @@ namespace Nektar
             }
             else
             {
-                SetCommonBC(type,n,time, cnt,inarray);
+                SetCommonBC(type, n, time, cnt, inarray);
             }
 
-            // no User Defined conditions provided so skip cnt 
-            // this line is left in case solver specific condition is added. 
-            cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize(); 
+            // no User Defined conditions provided so skip cnt
+            // this line is left in case solver specific condition is added.
+            cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
         }
     }
 }
diff --git a/solvers/CompressibleFlowSolver/EquationSystems/EulerCFE.cpp b/solvers/CompressibleFlowSolver/EquationSystems/EulerCFE.cpp
index 593c9f050cd158e6c16619d7b6aa410e56f165ff..728551536f52cb6afdf86ea4ca0eda853a2472b3 100755
--- a/solvers/CompressibleFlowSolver/EquationSystems/EulerCFE.cpp
+++ b/solvers/CompressibleFlowSolver/EquationSystems/EulerCFE.cpp
@@ -223,6 +223,14 @@ namespace Nektar
     {
         std::string varName;
         int cnt        = 0;
+        int nTracePts = GetTraceTotPoints();
+        
+        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
+        for (int i = 0; i < nvariables; ++i)
+        {
+            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
+            m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
+        }
         
         std::string userDefStr;
         int nreg = m_fields[0]->GetBndConditions().num_elements();
@@ -249,7 +257,7 @@ namespace Nektar
                 else
                 {
                     // set up userdefined BC common to all solvers
-                    SetCommonBC(userDefStr,n,time, cnt,inarray);
+                    SetCommonBC(userDefStr, n, time, cnt, Fwd, inarray);
                 }
             }
 
diff --git a/solvers/CompressibleFlowSolver/EquationSystems/NavierStokesCFE.cpp b/solvers/CompressibleFlowSolver/EquationSystems/NavierStokesCFE.cpp
index 77d347d70634d2ab40808a6df7d6ba764bb95471..305fbb706ec547eb414274523fabd3084a244fc5 100755
--- a/solvers/CompressibleFlowSolver/EquationSystems/NavierStokesCFE.cpp
+++ b/solvers/CompressibleFlowSolver/EquationSystems/NavierStokesCFE.cpp
@@ -251,12 +251,20 @@ namespace Nektar
     {
         std::string varName;
         int cnt        = 0;
+        int nTracePts = GetTraceTotPoints();
+
+        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
+        for (int i = 0; i < nvariables; ++i)
+        {
+            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
+            m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
+        }
 
         // loop over Boundary Regions
         for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
         {
             std::string type = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
-            SetCommonBC(type,n,time,cnt,inarray);
+            SetCommonBC(type, n, time, cnt, Fwd, inarray);
         }
     }
 }