diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp
index e3658b3fd..2c043ef2c 100644
--- a/src/M2ulPhyS.cpp
+++ b/src/M2ulPhyS.cpp
@@ -600,7 +600,7 @@ void M2ulPhyS::initVariables() {
 
   gradUp_A = new GradNonLinearForm(gradUpfes, intRules, dim, num_equation);
   gradUp_A->AddInteriorFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation));
-  gradUp_A->AddBdrFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation, bcIntegrator));
+  gradUp_A->AddBdrFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation, bcIntegrator, config.useBCinGrad));
 
   rhsOperator =
       new RHSoperator(iter, dim, num_equation, order, eqSystem, max_char_speed, intRules, intRuleType, d_fluxClass,
@@ -2874,6 +2874,8 @@ void M2ulPhyS::parseReactionInputs() {
 }
 
 void M2ulPhyS::parseBCInputs() {
+  tpsP->getInput("boundaryConditions/useBCinGrad", config.useBCinGrad, false);
+
   // number of BC regions defined
   int numWalls, numInlets, numOutlets;
   tpsP->getInput("boundaryConditions/numWalls", numWalls, 0);
diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp
index 790a43660..0f0061edd 100644
--- a/src/faceGradientIntegration.cpp
+++ b/src/faceGradientIntegration.cpp
@@ -34,8 +34,8 @@
 
 // Implementation of class FaceIntegrator
 GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation,
-                                       BCintegrator *bc)
-    : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc) {}
+                                       BCintegrator *bc, bool useBCinGrad)
+    : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc), useBCinGrad_(useBCinGrad) {}
 
 void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2,
                                             FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) {
@@ -96,18 +96,22 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
     if (Tr.Elem2No < 0) {
       assert(bc_ != NULL);
 
-      const int attr = Tr.Attribute;
-      std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = bc_->inletBCmap.find(attr);
-      std::unordered_map<int, BoundaryCondition *>::const_iterator obc = bc_->outletBCmap.find(attr);
-      std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = bc_->wallBCmap.find(attr);
-      if (ibc != bc_->inletBCmap.end()) {
-        ibc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
-      }
-      if (obc != bc_->outletBCmap.end()) {
-        obc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
-      }
-      if (wbc != bc_->wallBCmap.end()) {
-        wbc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+      if (useBCinGrad_) {
+        const int attr = Tr.Attribute;
+        std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = bc_->inletBCmap.find(attr);
+        std::unordered_map<int, BoundaryCondition *>::const_iterator obc = bc_->outletBCmap.find(attr);
+        std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = bc_->wallBCmap.find(attr);
+        if (ibc != bc_->inletBCmap.end()) {
+          ibc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+        if (obc != bc_->outletBCmap.end()) {
+          obc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+        if (wbc != bc_->wallBCmap.end()) {
+          wbc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+      } else {
+        iUp2 = iUp1;
       }
     } else {
       elfun2_mat.MultTranspose(shape2, iUp2);
diff --git a/src/faceGradientIntegration.hpp b/src/faceGradientIntegration.hpp
index 61ebfdef4..f2a16f0a3 100644
--- a/src/faceGradientIntegration.hpp
+++ b/src/faceGradientIntegration.hpp
@@ -45,10 +45,13 @@ class GradFaceIntegrator : public NonlinearFormIntegrator {
   const int dim;
   const int num_equation;
   IntegrationRules *intRules;
+
   BCintegrator *bc_;  // NB: GradFaceIntegrator is a friend of BCintegrator
+  const bool useBCinGrad_;
 
  public:
-  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation, BCintegrator *bc = NULL);
+  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation, BCintegrator *bc = NULL,
+                     bool useBCinGrad = false);
 
   virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr,
                                   const Vector &elfun, Vector &elvect);
diff --git a/src/run_configuration.hpp b/src/run_configuration.hpp
index 19e9c3930..46d6cc401 100644
--- a/src/run_configuration.hpp
+++ b/src/run_configuration.hpp
@@ -183,6 +183,9 @@ class RunConfiguration {
   std::vector<WallData> wallBC;
   std::vector<pair<int, WallType>> wallPatchType;
 
+  // Apply BCs in gradient calculation
+  bool useBCinGrad;
+
   // Passive scalar data
   Array<passiveScalarData*> arrayPassiveScalar;