diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp
index 023bc0422..913df9bf9 100644
--- a/src/dg/strong_dg.cpp
+++ b/src/dg/strong_dg.cpp
@@ -211,8 +211,8 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_derivatives(
 
     std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< ADArray > source_at_q(n_quad_pts);
-    std::vector< ADArray > physical_source_at_q(n_quad_pts);
+    std::vector< ADArray > source_at_q;
+    std::vector< ADArray > physical_source_at_q;
 
     // AD variable
     std::vector< FadType > soln_coeff(n_dofs_cell);
@@ -242,12 +242,14 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_derivatives(
         diss_phys_flux_at_q[iquad] = this->pde_physics_fad->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
         if(this->pde_physics_fad->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> real_quad_points = fe_values_vol.quadrature_point(iquad);
             dealii::Point<dim,FadType> ad_points;
             for (int d=0;d<dim;++d) { ad_points[d] = real_quad_points[d]; }
             physical_source_at_q[iquad] = this->pde_physics_fad->physical_source_term (ad_points, soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> real_quad_point = fe_values_vol.quadrature_point(iquad);
             dealii::Point<dim,FadType> ad_point;
             for (int d=0;d<dim;++d) { ad_point[d] = real_quad_point[d]; }
@@ -550,8 +552,8 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
 
     std::vector< realArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< realArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< realArray > source_at_q(n_quad_pts);
-    std::vector< realArray > physical_source_at_q(n_quad_pts);
+    std::vector< realArray > source_at_q;
+    std::vector< realArray > physical_source_at_q;
 
     // AD variable
     std::vector< realtype > soln_coeff(n_dofs_cell);
@@ -578,10 +580,12 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
         // Evaluate physical convective flux and source term
         conv_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
-        if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+        if(this->pde_physics_double->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             physical_source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->physical_source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad], current_cell_index);
         }
     }
@@ -637,7 +641,7 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
 
             // Physical source
-            if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+            if(this->pde_physics_double->has_nonzero_physical_source){
                 rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * physical_source_at_q[iquad][istate] * JxW[iquad];
             }
             // Source
diff --git a/src/dg/weak_dg.cpp b/src/dg/weak_dg.cpp
index 054651c34..c9c13c531 100644
--- a/src/dg/weak_dg.cpp
+++ b/src/dg/weak_dg.cpp
@@ -577,8 +577,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
 
     std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< doubleArray > source_at_q(n_quad_pts);
-    std::vector< doubleArray > physical_source_at_q(n_quad_pts);
+    std::vector< doubleArray > source_at_q;
+    std::vector< doubleArray > physical_source_at_q;
 
     std::vector< real > soln_coeff(n_soln_dofs_int);
     for (unsigned int idof = 0; idof < n_soln_dofs_int; ++idof) {
@@ -628,7 +628,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
         conv_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
-        if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+        if(this->pde_physics_double->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> points = fe_values_vol.quadrature_point(iquad);
             physical_source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->physical_source_term (points,soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
@@ -639,6 +640,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             }
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> point = fe_values_vol.quadrature_point(iquad);
             source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->source_term (point, soln_at_q[iquad], current_cell_index);
             //std::array<real,nstate> artificial_source_at_q = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->artificial_source_term (artificial_diss_coeff, point, soln_at_q[iquad]);
@@ -678,7 +680,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             //// Note that for diffusion, the negative is defined in the physics_double
             rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
             // Physical source
-            if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+            if(this->pde_physics_double->has_nonzero_physical_source){
                 rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * physical_source_at_q[iquad][istate] * JxW[iquad];
             }
             // Source
@@ -3623,8 +3625,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
 
     std::vector< ArrayTensor > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ArrayTensor > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< Array > source_at_q(n_quad_pts);
-    std::vector< Array > physical_source_at_q(n_quad_pts);
+    std::vector< Array > source_at_q;
+    std::vector< Array > physical_source_at_q;
 
     for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
         for (int istate=0; istate<nstate; istate++) {
@@ -3642,6 +3644,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
         diss_phys_flux_at_q[iquad] = physics.dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
         if(physics.has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             dealii::Point<dim,real2> ad_points;
             for (int d=0;d<dim;++d) { ad_points[d] = 0.0;}
             for (unsigned int idof = 0; idof < n_metric_dofs; ++idof) {
@@ -3661,6 +3664,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
         }
 
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             dealii::Point<dim,real2> ad_point;
             for (int d=0;d<dim;++d) { ad_point[d] = 0.0;}
             for (unsigned int idof = 0; idof < n_metric_dofs; ++idof) {
diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt
index d05d6e8da..9bae894d7 100644
--- a/src/physics/CMakeLists.txt
+++ b/src/physics/CMakeLists.txt
@@ -15,7 +15,7 @@ SET(SOURCE
     model.cpp
     large_eddy_simulation.cpp
     reynolds_averaged_navier_stokes.cpp
-    negative_spalart_allmaras.cpp
+    negative_spalart_allmaras_rans_model.cpp
     model_factory.cpp)
 
 foreach(dim RANGE 1 3)
diff --git a/src/physics/convection_diffusion.h b/src/physics/convection_diffusion.h
index 4ddcafc22..d11525ce5 100644
--- a/src/physics/convection_diffusion.h
+++ b/src/physics/convection_diffusion.h
@@ -62,9 +62,7 @@ class ConvectionDiffusion : public PhysicsBase <dim, nstate, real>
             diffusion_scaling_coeff(input_diffusion_coefficient),
             hasConvection(convection), 
             hasDiffusion(diffusion)
-    {
-        //static_assert(nstate<=5, "Physics::ConvectionDiffusion() should be created with nstate<=5");
-    };
+    {};
 
     /// Destructor
     ~ConvectionDiffusion () {};
diff --git a/src/physics/model.cpp b/src/physics/model.cpp
index 21e05b0c9..089eb30f7 100644
--- a/src/physics/model.cpp
+++ b/src/physics/model.cpp
@@ -44,37 +44,37 @@ template <int dim, int nstate, typename real>
 void ModelBase<dim, nstate, real>
 ::boundary_face_values (
    const int /*boundary_type*/,
-   const dealii::Point<dim, real> &pos,
-   const dealii::Tensor<1,dim,real> &normal_int,
-   const std::array<real,nstate> &soln_int,
-   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   const dealii::Point<dim, real> &/*pos*/,
+   const dealii::Tensor<1,dim,real> &/*normal_int*/,
+   const std::array<real,nstate> &/*soln_int*/,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
    std::array<real,nstate> &soln_bc,
    std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
 {
-    std::array<real,nstate> boundary_values;
-    std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
-    for (int s=0; s<nstate; s++) {
-        boundary_values[s] = this->manufactured_solution_function->value (pos, s);
-        boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
-    }
+    //std::array<real,nstate> boundary_values;
+    //std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
+    //for (int s=0; s<nstate; s++) {
+    //    boundary_values[s] = this->manufactured_solution_function->value (pos, s);
+    //    boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
+    //}
 
     for (int istate=0; istate<dim+2; ++istate) {
         soln_bc[istate] = 0.0;
         soln_grad_bc[istate] = 0.0;
     }
-    for (int istate=dim+2; istate<nstate; ++istate) {
-
-        std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
-        const bool inflow = (characteristic_dot_n[istate] <= 0.);
-
-        if (inflow) { // Dirichlet boundary condition
-            soln_bc[istate] = boundary_values[istate];
-            soln_grad_bc[istate] = soln_grad_int[istate];
-        } else { // Neumann boundary condition
-            soln_bc[istate] = soln_int[istate];
-            soln_grad_bc[istate] = soln_grad_int[istate];
-        }
-    }
+    //for (int istate=dim+2; istate<nstate; ++istate) {
+//
+    //    std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
+    //    const bool inflow = (characteristic_dot_n[istate] <= 0.);
+//
+    //    if (inflow) { // Dirichlet boundary condition
+    //        soln_bc[istate] = boundary_values[istate];
+    //        soln_grad_bc[istate] = soln_grad_int[istate];
+    //    } else { // Neumann boundary condition
+    //        soln_bc[istate] = soln_int[istate];
+    //        soln_grad_bc[istate] = soln_grad_int[istate];
+    //    }
+    //}
 }
 //----------------------------------------------------------------
 template <int dim, int nstate, typename real>
diff --git a/src/physics/model.h b/src/physics/model.h
index bd1f00f24..6744a9c82 100644
--- a/src/physics/model.h
+++ b/src/physics/model.h
@@ -78,8 +78,7 @@ class ModelBase
         std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;
 
     /// Returns current vector solution to be used by PhysicsPostprocessor to output current solution.
-    /** The implementation in this Model base class simply returns the stored solution.
-     */
+    /** The implementation in this Model base class simply returns the stored solution. */
     virtual dealii::Vector<double> post_compute_derived_quantities_vector (
         const dealii::Vector<double>              &uh,
         const std::vector<dealii::Tensor<1,dim> > &/*duh*/,
@@ -88,13 +87,11 @@ class ModelBase
         const dealii::Point<dim>                  &/*evaluation_points*/) const;
 
     /// Returns DataComponentInterpretation of the solution to be used by PhysicsPostprocessor to output current solution.
-    /** Treats every solution state as an independent scalar.
-     */
+    /** Treats every solution state as an independent scalar. */
     virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const;
 
     /// Returns names of the solution to be used by PhysicsPostprocessor to output current solution.
-    /** The implementation in this Model base class simply returns "state(dim+2+0), state(dim+2+1), etc.".
-     */
+    /** The implementation in this Model base class simply returns "state(dim+2+0), state(dim+2+1), etc.". */
     virtual std::vector<std::string> post_get_names () const;
 
     // Quantities needed to be updated by DG for the model -- accomplished by DGBase update_model_variables()
diff --git a/src/physics/model_factory.cpp b/src/physics/model_factory.cpp
index ee098233c..999f19d9e 100644
--- a/src/physics/model_factory.cpp
+++ b/src/physics/model_factory.cpp
@@ -9,7 +9,7 @@
 #include "manufactured_solution.h"
 #include "large_eddy_simulation.h"
 #include "reynolds_averaged_navier_stokes.h"
-#include "negative_spalart_allmaras.h"
+#include "negative_spalart_allmaras_rans_model.h"
 
 namespace PHiLiP {
 namespace Physics {
diff --git a/src/physics/negative_spalart_allmaras.cpp b/src/physics/negative_spalart_allmaras_rans_model.cpp
similarity index 91%
rename from src/physics/negative_spalart_allmaras.cpp
rename to src/physics/negative_spalart_allmaras_rans_model.cpp
index 32402204a..403e0ddb1 100644
--- a/src/physics/negative_spalart_allmaras.cpp
+++ b/src/physics/negative_spalart_allmaras_rans_model.cpp
@@ -6,7 +6,7 @@
 
 #include "model.h"
 #include "reynolds_averaged_navier_stokes.h"
-#include "negative_spalart_allmaras.h"
+#include "negative_spalart_allmaras_rans_model.h"
 
 namespace PHiLiP {
 namespace Physics {
@@ -586,6 +586,67 @@ ::compute_production_dissipation_cross_term (
 }
 //----------------------------------------------------------------
 template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_wall (
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    // Wall boundary condition for nu_tilde (working variable of negative SA model)
+    // nu_tilde = 0
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = 0.0;
+        soln_grad_bc[istate] = 0.0;
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_outflow (
+   const std::array<real,nstate> &soln_int,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = soln_int[istate];
+        soln_grad_bc[istate] = soln_grad_int[istate];
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_inflow (
+   const std::array<real,nstate> &/*soln_int*/,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    const real density_bc = this->navier_stokes_physics->density_inf;
+    const real dynamic_viscosity_coefficient_bc = this->navier_stokes_physics->viscosity_coefficient_inf;
+    const real kinematic_viscosity_coefficient_bc = dynamic_viscosity_coefficient_bc/density_bc;
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = density_bc*3.0*kinematic_viscosity_coefficient_bc;
+        soln_grad_bc[istate] = soln_grad_int[istate];
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_farfield (
+   std::array<real,nstate> &soln_bc) const
+{
+   const real density_bc = this->navier_stokes_physics->density_inf;
+   const real dynamic_viscosity_coefficient_bc = this->navier_stokes_physics->viscosity_coefficient_inf;
+   const real kinematic_viscosity_coefficient_bc = dynamic_viscosity_coefficient_bc/density_bc;
+   
+   // Farfield boundary condition for nu_tilde (working variable of negative SA model)
+   // nu_tilde = 3.0*kinematic_viscosity_inf to 5.0*kinematic_viscosity_inf
+   for (int istate=dim+2; istate<nstate; ++istate) {
+      soln_bc[istate] = density_bc*3.0*kinematic_viscosity_coefficient_bc;
+   }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
 dealii::Vector<double> ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
 ::post_compute_derived_quantities_vector (
     const dealii::Vector<double>              &uh,
diff --git a/src/physics/negative_spalart_allmaras.h b/src/physics/negative_spalart_allmaras_rans_model.h
similarity index 91%
rename from src/physics/negative_spalart_allmaras.h
rename to src/physics/negative_spalart_allmaras_rans_model.h
index cec51ecde..ee66fe30f 100644
--- a/src/physics/negative_spalart_allmaras.h
+++ b/src/physics/negative_spalart_allmaras_rans_model.h
@@ -1,5 +1,5 @@
-#ifndef __NEGATIVE_SPALART_ALLMARAS__
-#define __NEGATIVE_SPALART_ALLMARAS__
+#ifndef __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__
+#define __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__
 
 #include "euler.h"
 #include "navier_stokes.h"
@@ -115,6 +115,35 @@ class ReynoldsAveragedNavierStokes_SAneg : public ReynoldsAveragedNavierStokesBa
     template<typename real2> real2 scale_coefficient(
         const real2 coefficient) const;
 
+    /// Wall boundary condition
+    /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."
+     *  eq.(7)
+     */ 
+    void boundary_wall (
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Inflow boundary condition
+    void boundary_outflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Inflow boundary condition
+    void boundary_inflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Farfield boundary conditions based on freestream values
+    /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."
+     *  eq.(8)
+     */ 
+    void boundary_farfield (
+        std::array<real,nstate> &soln_bc) const;
+
 private:
     /// Templated nondimensionalized eddy viscosity for the negative SA model.
     /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."
diff --git a/src/physics/physics_model.cpp b/src/physics/physics_model.cpp
index 1f56b109e..b6d8936b8 100644
--- a/src/physics/physics_model.cpp
+++ b/src/physics/physics_model.cpp
@@ -248,6 +248,11 @@ ::boundary_face_values (
         std::array<real,nstate_baseline_physics> baseline_soln_bc;
         std::array<dealii::Tensor<1,dim,real>,nstate_baseline_physics> baseline_soln_grad_bc;
 
+        for (int istate=0; istate<nstate_baseline_physics; istate++) {
+            baseline_soln_bc[istate]      = 0;
+            baseline_soln_grad_bc[istate] = 0;
+        }
+
         physics_baseline->boundary_face_values(
                 boundary_type, pos, normal_int, baseline_soln_int, baseline_soln_grad_int, 
                 baseline_soln_bc, baseline_soln_grad_bc);
diff --git a/src/physics/reynolds_averaged_navier_stokes.cpp b/src/physics/reynolds_averaged_navier_stokes.cpp
index 8df0760bb..1c41f4daf 100644
--- a/src/physics/reynolds_averaged_navier_stokes.cpp
+++ b/src/physics/reynolds_averaged_navier_stokes.cpp
@@ -701,6 +701,86 @@ ::physical_source_source_term (
     return physical_source_source_term;
 }
 //----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokesBase<dim,nstate,real>
+::boundary_manufactured_solution (
+    const dealii::Point<dim, real> &pos,
+    const dealii::Tensor<1,dim,real> &normal_int,
+    const std::array<real,nstate> &soln_int,
+    const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+    std::array<real,nstate> &soln_bc,
+    std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    // Manufactured solution
+    std::array<real,nstate> conservative_boundary_values;
+    std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
+    for (int s=0; s<nstate; s++) {
+        conservative_boundary_values[s] = this->manufactured_solution_function->value (pos, s);
+        boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
+    }
+
+    for (int istate=dim+2; istate<nstate; ++istate) {
+
+        std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int);
+        const bool inflow = (characteristic_dot_n[istate] <= 0.);
+
+        if (inflow) { // Dirichlet boundary condition
+            soln_bc[istate] = conservative_boundary_values[istate];
+            soln_grad_bc[istate] = soln_grad_int[istate];
+        } else { // Neumann boundary condition
+            soln_bc[istate] = soln_int[istate];
+            soln_grad_bc[istate] = soln_grad_int[istate];
+        }
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokesBase<dim,nstate,real>
+::boundary_face_values (
+   const int boundary_type,
+   const dealii::Point<dim, real> &pos,
+   const dealii::Tensor<1,dim,real> &normal_int,
+   const std::array<real,nstate> &soln_int,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    if (boundary_type == 1000) {
+        // Manufactured solution boundary condition
+        boundary_manufactured_solution (pos, normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1001) {
+        // Wall boundary condition for working variables of RANS turbulence model
+        boundary_wall (soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1002) {
+        // Outflow boundary condition 
+        boundary_outflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1003) {
+        // Inflow boundary condition
+        boundary_inflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1004) {
+        // Riemann-based farfield boundary condition
+        std::cout << "Riemann boundary condition is not implemented!" << std::endl;
+        std::abort();
+    } 
+    else if (boundary_type == 1005) {
+        // Simple farfield boundary condition
+        boundary_farfield(soln_bc);
+    } 
+    else if (boundary_type == 1006) {
+        // Slip wall boundary condition
+        std::cout << "Slip wall boundary condition is not implemented!" << std::endl;
+        std::abort();
+    } 
+    else {
+        std::cout << "Invalid boundary_type: " << boundary_type << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
 //----------------------------------------------------------------
 //----------------------------------------------------------------
 // Instantiate explicitly
diff --git a/src/physics/reynolds_averaged_navier_stokes.h b/src/physics/reynolds_averaged_navier_stokes.h
index 4d3f0b006..3d82a212b 100644
--- a/src/physics/reynolds_averaged_navier_stokes.h
+++ b/src/physics/reynolds_averaged_navier_stokes.h
@@ -120,6 +120,16 @@ class ReynoldsAveragedNavierStokesBase : public ModelBase <dim, nstate, real>
         const std::array<real,nstate> &conservative_solution,
         const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const = 0;
 
+    /// Boundary condition handler
+    void boundary_face_values (
+        const int /*boundary_type*/,
+        const dealii::Point<dim, real> &/*pos*/,
+        const dealii::Tensor<1,dim,real> &/*normal*/,
+        const std::array<real,nstate> &/*soln_int*/,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
+        std::array<real,nstate> &/*soln_bc*/,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const override;
+
 protected:
     /// Returns the square of the magnitude of the vector 
     template<typename real2> 
@@ -236,6 +246,38 @@ class ReynoldsAveragedNavierStokesBase : public ModelBase <dim, nstate, real>
     std::array<real,nstate> physical_source_source_term (
         const dealii::Point<dim,real> &pos,
         const dealii::types::global_dof_index cell_index) const;
+
+    /// Evaluate the manufactured solution boundary conditions.
+    void boundary_manufactured_solution (
+        const dealii::Point<dim, real> &pos,
+        const dealii::Tensor<1,dim,real> &normal_int,
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Wall boundary condition
+    virtual void boundary_wall (
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const = 0;
+
+    /// Outflow Boundary Condition 
+    virtual void boundary_outflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const = 0;
+
+    /// Inflow boundary conditions 
+    virtual void boundary_inflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const = 0;
+
+    /// Farfield boundary conditions based on freestream values
+    virtual void boundary_farfield (
+        std::array<real,nstate> &soln_bc) const = 0;
 };
 
 } // Physics namespace
diff --git a/tests/unit_tests/navier_stokes_unit_test/reynolds_averaged_navier_stokes_sa_neg_manufactured_solution_source.cpp b/tests/unit_tests/navier_stokes_unit_test/reynolds_averaged_navier_stokes_sa_neg_manufactured_solution_source.cpp
index 340d7cab7..0df554eb4 100644
--- a/tests/unit_tests/navier_stokes_unit_test/reynolds_averaged_navier_stokes_sa_neg_manufactured_solution_source.cpp
+++ b/tests/unit_tests/navier_stokes_unit_test/reynolds_averaged_navier_stokes_sa_neg_manufactured_solution_source.cpp
@@ -13,7 +13,7 @@
 #include "assert_compare_array.h"
 #include "parameters/parameters.h"
 #include "physics/reynolds_averaged_navier_stokes.h"
-#include "physics/negative_spalart_allmaras.h"
+#include "physics/negative_spalart_allmaras_rans_model.h"
 
 const double TOLERANCE = 1E-5;