diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp
index a51a7f46b..2801f7af0 100644
--- a/src/oneD/Boundary1D.cpp
+++ b/src/oneD/Boundary1D.cpp
@@ -191,6 +191,8 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
             // The third flow residual is for T, where it is set to T(0).  Subtract
             // the local temperature to hold the flow T to the inlet T.
             rb[c_offset_T] -= m_temp;
+        } else {
+            rb[c_offset_T] -= m_flow->T_fixed(0);
         }
 
         if (m_flow->isFree()) {
@@ -221,6 +223,8 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
         rb[c_offset_V] -= m_V0;
         if (m_flow->doEnergy(m_flow->nPoints() - 1)) {
             rb[c_offset_T] -= m_temp; // T
+        } else {
+            rb[c_offset_T] -= m_flow->T_fixed(m_flow->nPoints() - 1);
         }
         rb[c_offset_U] += m_mdot; // u
         for (size_t k = 0; k < m_nsp; k++) {
@@ -380,6 +384,11 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg,
         for (size_t k = c_offset_Y; k < nc; k++) {
             rb[k] = xb[k] - xb[k + nc];
         }
+        if (m_flow_left->doEnergy(0)) {
+            rb[c_offset_T] = xb[c_offset_T + nc] - xb[c_offset_T]; // zero T gradient
+        } else {
+            rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(0);
+        }
     }
 
     if (m_flow_left) {
@@ -389,8 +398,11 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg,
         double* rb = r - nc;
         int* db = diag - nc;
 
-        if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
+        size_t last = m_flow_left->nPoints() - 1;
+        if (m_flow_left->doEnergy(last)) {
             rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
+        } else {
+            rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(last);
         }
         size_t kSkip = c_offset_Y + m_flow_left->rightExcessSpecies();
         for (size_t k = c_offset_Y; k < nc; k++) {
@@ -467,13 +479,11 @@ void OutletRes1D::eval(size_t jg, double* xg, double* rg,
         double* xb = x;
         double* rb = r;
 
-        // this seems wrong...
-        // zero Lambda
-        rb[c_offset_U] = xb[c_offset_L];
-
         if (m_flow_right->doEnergy(0)) {
             // zero gradient for T
             rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc];
+        } else {
+            rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(0);
         }
 
         // specified mass fractions
@@ -488,11 +498,11 @@ void OutletRes1D::eval(size_t jg, double* xg, double* rg,
         double* rb = r - nc;
         int* db = diag - nc;
 
-        if (!m_flow_left->usesLambda()) {
-            rb[c_offset_L] = xb[c_offset_L]; // zero Lambda
-        }
-        if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
-            rb[c_offset_T] = xb[c_offset_T] - m_temp; // zero dT/dz
+        size_t last = m_flow_left->nPoints() - 1;
+        if (m_flow_left->doEnergy(last)) {
+            rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
+        } else {
+            rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(last);
         }
         size_t kSkip = m_flow_left->rightExcessSpecies();
         for (size_t k = c_offset_Y; k < nc; k++) {
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 3c495b1fe..f0a7af653 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -480,16 +480,12 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
             // a result, these residual equations will force the solution
             // variables to the values for the boundary object
             rsd[index(c_offset_V,0)] = V(x,0);
-            if (doEnergy(0)) {
-                rsd[index(c_offset_T,0)] = T(x,0);
-            } else {
-                rsd[index(c_offset_T,0)] = T(x,0) - T_fixed(0);
-            }
-            if (m_isFree) {
-                rsd[index(c_offset_L, 0)] = lambda(x, 0);
-            } else {
-                // used to set mass flow rate (both axisymmetric and unstrained)
+            rsd[index(c_offset_T,0)] = T(x,0);
+            if (m_usesLambda) {
                 rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
+            } else {
+                rsd[index(c_offset_L, 0)] = lambda(x, 0);
+                diag[index(c_offset_L, 0)] = 0;
             }
 
             // The default boundary condition for species is zero flux. However,
@@ -574,10 +570,10 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
                 diag[index(c_offset_T, j)] = 0;
             }
 
-            if (m_isFree) {
-                rsd[index(c_offset_L, j)] = lambda(x, j);
-            } else {
+            if (m_usesLambda) {
                 rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
+            } else {
+                rsd[index(c_offset_L, j)] = lambda(x, j);
             }
             diag[index(c_offset_L, j)] = 0;
         }
@@ -1015,8 +1011,6 @@ void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
 
     rsd[index(c_offset_V,j)] = V(x,j);
     double sum = 0.0;
-    rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
-    diag[index(c_offset_L, j)] = 0;
     // set residual of poisson's equ to zero
     rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
     for (size_t k = 0; k < m_nsp; k++) {
@@ -1025,24 +1019,16 @@ void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
     }
     rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
     diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
-    if (m_isFree) {
-        rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
-        rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
-    } else if (m_usesLambda) {
-        rsd[index(c_offset_U,j)] = rho_u(x,j);
-        if (m_do_energy[j]) {
-            rsd[index(c_offset_T,j)] = T(x,j);
-        } else {
-            rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
-        }
+    if (m_usesLambda) {
+        rsd[index(c_offset_U, j)] = rho_u(x, j);
+        rsd[index(c_offset_L, j)] = lambda(x, j);
     } else {
-        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
-        if (m_do_energy[j]) {
-            rsd[index(c_offset_T, j)] = T(x, j);
-        } else {
-            rsd[index(c_offset_T, j)] = T(x, j) - T_fixed(j);
-        }
+        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1);
+        rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1);
     }
+
+    diag[index(c_offset_L, j)] = 0;
+    rsd[index(c_offset_T, j)] = T(x, j);
 }
 
 void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
@@ -1079,6 +1065,7 @@ void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double
                 - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
         }
     } else {
+        // unstrained
         rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
     }
 }