diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index d4f94008ddc..081c65206ef 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1069,7 +1069,9 @@ class CConfig {
   bool Streamwise_Periodic_Temperature;              /*!< \brief Use real periodicity for Energy equation or otherwise outlet source term. */
   su2double Streamwise_Periodic_PressureDrop;        /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
   su2double Streamwise_Periodic_TargetMassFlow;      /*!< \brief Value of prescribed massflow [kg/s] which results in an delta p and therefore an artificial body force vector. */
+  su2double Streamwise_Periodic_ComputedMassFlow;      /*!< \brief Value of computed massflow [kg/s] corresponding to a prescribed delta p. */
   su2double Streamwise_Periodic_OutletHeat;          /*!< /brief Heatflux boundary [W/m^2] imposed at streamwise periodic outlet. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< /brief Exp coefficient for iso-thermal BCs Streamwise Periodic. */
 
   su2double *FreeStreamTurboNormal;     /*!< \brief Direction to initialize the flow in turbomachinery computation */
   su2double Restart_Bandwidth_Agg;      /*!< \brief The aggregate of the bandwidth for writing binary restarts (to be averaged later). */
@@ -3187,6 +3189,12 @@ class CConfig {
    */
   unsigned short GetnMarker_HeatFlux(void) const { return nMarker_HeatFlux; }
 
+  /*!
+   * \brief Get the total (local) number of isothermal markers.
+   * \return Total number of isothermal markers.
+   */
+  unsigned short GetnMarker_Isothermal(void) const { return nMarker_Isothermal; }
+
   /*!
    * \brief Get the total number of rough markers.
    * \return Total number of heat flux markers.
@@ -6219,6 +6227,17 @@ class CConfig {
    */
   void SetStreamwise_Periodic_PressureDrop(su2double Streamwise_Periodic_PressureDrop_) { Streamwise_Periodic_PressureDrop = Streamwise_Periodic_PressureDrop_; }
 
+  /*!
+   * \brief Get the value of the MassFlow from which body force vector is computed.
+   * \return MassFlow for body force computation.
+   */
+  su2double GetStreamwise_Periodic_ComputedMassFlow(void) const { return Streamwise_Periodic_ComputedMassFlow; }
+
+  /*!
+   * \brief Set the value of the MassFlow from which body force vector is computed. Necessary for Restart metadata??
+   */
+  void SetStreamwise_Periodic_ComputedMassFlow(su2double Streamwise_Periodic_MassFlow_) { Streamwise_Periodic_ComputedMassFlow = Streamwise_Periodic_MassFlow_; }
+
   /*!
    * \brief Get the value of the massflow from which body force vector is computed.
    * \return Massflow for body force computation.
@@ -9873,4 +9892,16 @@ class CConfig {
    */
   LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; }
 
+  /*!
+   * \brief Set Lambda L for Streamwise Periodic
+   * \return bool
+   */
+  void SetStreamwise_Periodic_LamdaL(su2double value) { Streamwise_Periodic_LambdaL = value; }
+
+  /*!
+   * \brief Get Lambda L for Streamwise Periodic
+   * \return su2double
+   */
+  su2double GetStreamwise_Periodic_LamdaL(void) const { return Streamwise_Periodic_LambdaL; }
+
 };
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 64c781a34f8..222f17a5c5b 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -1993,6 +1993,9 @@ enum ENUM_OBJECTIVE {
   TOPOL_DISCRETENESS = 63,      /*!< \brief Measure of the discreteness of the current topology. */
   TOPOL_COMPLIANCE = 64,        /*!< \brief Measure of the discreteness of the current topology. */
   STRESS_PENALTY = 65,          /*!< \brief Penalty function of VM stresses above a maximum value. */
+  STREAMWISE_LAMBDAL= 71, /*!< /brief temp. expo. coefficient for iso-thermal BCs Streamwise Periodic. */
+  STREAMWISE_PERIODIC_DP= 72, /*!< /brief pressure drop in Streamwise Periodic flow domain. */
+  STREAMWISE_PERIODIC_MASSFLOW= 73, /*!< /brief massflow rate in Streamwise Periodic flow domain. */
 };
 static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("DRAG", DRAG_COEFFICIENT)
@@ -2028,6 +2031,9 @@ static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("SURFACE_PRESSURE_DROP", SURFACE_PRESSURE_DROP)
   MakePair("SURFACE_SPECIES_0", SURFACE_SPECIES_0)
   MakePair("SURFACE_SPECIES_VARIANCE", SURFACE_SPECIES_VARIANCE)
+  MakePair("STREAMWISE_LAMBDAL", STREAMWISE_LAMBDAL)
+  MakePair("STREAMWISE_PERIODIC_DP", STREAMWISE_PERIODIC_DP)
+  MakePair("STREAMWISE_PERIODIC_MASSFLOW", STREAMWISE_PERIODIC_MASSFLOW)
   MakePair("CUSTOM_OBJFUNC", CUSTOM_OBJFUNC)
   MakePair("REFERENCE_GEOMETRY", REFERENCE_GEOMETRY)
   MakePair("REFERENCE_NODE", REFERENCE_NODE)
@@ -2655,6 +2661,8 @@ struct StreamwisePeriodicValues {
   su2double Streamwise_Periodic_InletTemperature;   /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
   su2double Streamwise_Periodic_BoundaryArea;       /*!< \brief Global Surface area of the streamwise periodic interface. */
   su2double Streamwise_Periodic_AvgDensity;         /*!< \brief Area avg density on the periodic interface. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< \brief Temperature Gradient in iso-thermal BCs. */
+  su2double Streamwise_Periodic_ThetaScaling;       /*!< \brief Scaled Temperature for iso-thermal BCs. */
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index bf2122f01a7..4c42b424841 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -5150,8 +5150,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
     if (nMarker_PerBound == 0)
       SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
-    if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
-      SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
+    // if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal > 0 && nMarker_HeatFlux > 0)
+    //   SU2_MPI::Error("MARKER_ISOTHERMAL and MARKER_HEATFLUX are not allowed simultaneously with STREAMWISE_PERIODIC_TEMPERATURE= YES, only one of these is allowed MARKER_ISOTHERMAL, MARKER_HEATFLUX in a fluid zone.", CURRENT_FUNCTION);
     if (Ref_Inc_NonDim != DIMENSIONAL)
       SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION);
     if (Axisymmetric)
diff --git a/SU2_CFD/include/interfaces/CInterface.hpp b/SU2_CFD/include/interfaces/CInterface.hpp
index 9ec0599827e..9be29da3a5c 100644
--- a/SU2_CFD/include/interfaces/CInterface.hpp
+++ b/SU2_CFD/include/interfaces/CInterface.hpp
@@ -219,4 +219,7 @@ class CInterface {
    */
   void GatherAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone);
 
+  inline virtual void SetLambdaL(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
+                          unsigned long Marker_Target) {};
+
 };
diff --git a/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp b/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp
index 919d4e8f7f1..49a957dfeb9 100644
--- a/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp
+++ b/SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp
@@ -70,4 +70,7 @@ class CConjugateHeatInterface : public CInterface {
    */
   void SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
                           unsigned long Marker_Target, unsigned long Vertex_Target, unsigned long Point_Target) override;
+
+  void SetLambdaL(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
+                          unsigned long Marker_Target) override;
 };
diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp
index ab9bedc8037..18ee562b414 100644
--- a/SU2_CFD/include/numerics/flow/flow_sources.hpp
+++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp
@@ -391,6 +391,8 @@ class CSourceIncStreamwise_Periodic final : public CSourceBase_Flow {
   bool turbulent; /*!< \brief Turbulence model used. */
   bool energy;    /*!< \brief Energy equation on. */
   bool streamwisePeriodic_temperature; /*!< \brief Periodicity in energy equation */
+  bool bool_heat_flux_bc; /*!< \brief HEAT Flux BC boolean. */
+  bool bool_isotherml_bc; /*!< \brief ISO THERMAL BC boolean. */
   su2double Streamwise_Coord_Vector[MAXNDIM] = {0.0}; /*!< \brief Translation vector between streamwise periodic surfaces. */
 
   su2double norm2_translation, /*!< \brief Square of distance between the 2 periodic surfaces. */
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 8fcb6d05bec..b22b9116925 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -2933,6 +2933,15 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
     case SURFACE_SPECIES_VARIANCE:
       objFun += weight * config.GetSurface_Species_Variance(0);
       break;
+    case STREAMWISE_LAMBDAL:
+      objFun += weight * config.GetStreamwise_Periodic_LamdaL();
+      break;
+    case STREAMWISE_PERIODIC_DP:
+      objFun += weight * config.GetStreamwise_Periodic_PressureDrop();
+      break;
+    case STREAMWISE_PERIODIC_MASSFLOW:
+      objFun += weight * config.GetStreamwise_Periodic_ComputedMassFlow();
+      break;
     case CUSTOM_OBJFUNC:
       objFun += weight * Total_Custom_ObjFunc;
       break;
diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp
index 627983ec096..5b0c853f03b 100644
--- a/SU2_CFD/include/solvers/CHeatSolver.hpp
+++ b/SU2_CFD/include/solvers/CHeatSolver.hpp
@@ -318,6 +318,20 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
                    CSolver **solver_container,
                    CConfig *config) override;
 
+/*!
+   * \brief Compute the viscous residuals for the turbulent equation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics_container - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
+   */
+  void Source_Residual(CGeometry *geometry,
+                        CSolver **solver_container,
+                        CNumerics **numerics_container,
+                        CConfig *config,
+                        unsigned short iMesh) override;
   /*!
    * \brief Get value of the heat load (integrated heat flux).
    * \return Value of the heat load (integrated heat flux).
diff --git a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
index 54ad481b488..89463b47dc1 100644
--- a/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
+++ b/SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp
@@ -171,3 +171,47 @@ void CConjugateHeatInterface::SetTarget_Variable(CSolver *target_solution, CGeom
                                               target_config->GetRelaxation_Factor_CHT(), Target_Variable[3]);
   }
 }
+
+void CConjugateHeatInterface::SetLambdaL(CSolver *target_solution, CGeometry *target_geometry,
+                                                 const CConfig *target_config, unsigned long Marker_Target) {
+  unsigned short iPoint, iMarker, iVertex;
+  su2double term1=0.0, term2=0.0, term3=0.0;
+  auto nDim = target_geometry->GetnDim();
+
+    /*--- Loop over all heatflux Markers ---*/
+  for (iMarker = 0; iMarker < target_config->GetnMarker_All(); iMarker++) {
+
+    if (target_config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) {
+
+    /*--- Identify the boundary by string name and retrive heatflux from config ---*/
+    const auto Marker_StringTag = target_config->GetMarker_All_TagBound(iMarker);
+    const su2double Wall_HeatFlux = target_config->GetWall_HeatFlux(Marker_StringTag);
+
+      for (iVertex = 0ul; iVertex < target_geometry->nVertex[iMarker]; iVertex++) {
+
+      auto iPoint = target_geometry->vertex[iMarker][iVertex]->GetNode();
+
+      if (!target_geometry->nodes->GetDomain(iPoint)) continue;
+
+      const auto AreaNormal = target_geometry->vertex[iMarker][iVertex]->GetNormal();
+
+      const su2double FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
+
+      // compute the constant
+      const auto GradT = target_solution->GetNodes()->GetGradient_Primitive(iPoint)[0]; //TODO: check index
+       
+      term1 += GeometryToolbox::DotProduct(nDim, GradT, AreaNormal);
+
+      // compute coefficient of lambda
+      term2 += target_solution->GetNodes()->GetTemperature(iPoint) * FaceArea;
+      }
+    } // loop Vertices
+  }
+
+  // compute coefficient of lambda^2
+for (iPoint = 0; iPoint<=target_geometry->GetnPointDomain(); iPoint++) {
+  term3 += target_solution->GetNodes()->GetTemperature(iPoint) * target_geometry->nodes->GetVolume(iPoint);
+}
+
+  // Compute LambdaL using the quadratic equation
+}
\ No newline at end of file
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index 9ed1305f3ce..3f8d9f3b830 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -723,6 +723,8 @@ CSourceIncStreamwise_Periodic::CSourceIncStreamwise_Periodic(unsigned short val_
   turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
   energy    = config->GetEnergy_Equation();
   streamwisePeriodic_temperature = config->GetStreamwise_Periodic_Temperature();
+  bool_heat_flux_bc = (config->GetnMarker_HeatFlux() > 0);
+  bool_isotherml_bc = (config->GetnMarker_Isothermal() > 0);
 
   for (unsigned short iDim = 0; iDim < nDim; iDim++)
     Streamwise_Coord_Vector[iDim] = config->GetPeriodic_Translation(0)[iDim];
@@ -737,6 +739,7 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
 
   /* Value of prescribed pressure drop which results in an artificial body force vector. */
   const su2double delta_p = SPvals.Streamwise_Periodic_PressureDrop;
+  auto implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
 
   for (unsigned short iVar = 0; iVar < nVar; iVar++) residual[iVar] = 0.0;
 
@@ -749,23 +752,65 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
   /*--- Compute the periodic temperature contribution to the energy equation, if energy equation is considered ---*/
   if (energy && streamwisePeriodic_temperature) {
 
-    scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
+    /*--- Compute scalar factors, if there is a HEAT FLUX BC ---*/
+    if (bool_heat_flux_bc) {
+      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
 
-    /*--- Compute scalar-product dot_prod(v*t) ---*/
-    dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+      /*--- Compute scalar-product dot_prod(v*t) ---*/
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+    }
+
+    /*--- Compute scalar factors, if there is an ISOTHERMAL BC ---*/
+    if (bool_isotherml_bc) {
+
+      /*--- Compute three terms associated with the source terms for iso-thermal BCs ---*/
+      // Reference Eq(20) Stalio et. al, doi:10.1115/1.2717235
+
+      // Displacement with temperature
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, PrimVar_Grad_i[nDim+1]);
+      
+      su2double u_i = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+      su2double temp_i = V_i[nDim+1];
+      su2double term1 = 2 * Thermal_Conductivity_i * dot_product;
+      su2double term2 = - SPvals.Streamwise_Periodic_LambdaL * Thermal_Conductivity_i * temp_i;
+      su2double term3 = - temp_i * u_i * DensityInc_i *  config->GetSpecific_Heat_Cp();
+      scalar_factor = SPvals.Streamwise_Periodic_LambdaL * (term1 + term2 + term3);
+      dot_product = 1.0;
+    }
 
     residual[nDim+1] = Volume * scalar_factor * dot_product;
 
-    /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
-    if(turbulent) {
+  if (implicit) {
+
+    // /*--- Jacobian is set to zero on initialization. ---*/
 
-      /*--- Compute a scalar factor ---*/
-      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+    // jacobian[nDim+1][nDim+1] = Volume * scalar_factor * dot_product;
+
+  }
 
+    /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
+    if(turbulent) {
+      if (bool_heat_flux_bc) {
+        /*--- Compute a scalar factor ---*/
+        scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+      }
+      
+      if (bool_isotherml_bc) {
+        scalar_factor = (-V_i[nDim + 1] * SPvals.Streamwise_Periodic_LambdaL) * config->GetSpecific_Heat_Cp() / Prandtl_Turb;
+      }
+      
       /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
       dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
 
       residual[nDim+1] -= Volume * scalar_factor * dot_product;
+
+        if (implicit) {
+
+        // /*--- Jacobian is set to zero on initialization. ---*/
+
+        // jacobian[nDim+1][nDim+1] -= Volume * scalar_factor * dot_product;
+
+        }
     } // if turbulent
   } // if energy
 
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
index cd275b8f51e..1770a9963f8 100644
--- a/SU2_CFD/src/output/CFlowIncOutput.cpp
+++ b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -178,6 +178,8 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){
     AddHistoryOutput("STREAMWISE_MASSFLOW", "SWMassflow", ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Massflow in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_DP",       "SWDeltaP",   ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Pressure drop in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_HEAT",     "SWHeat",     ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Integrated heat for streamwise periodic flow");
+    if (config->GetnMarker_Isothermal() > 0)
+      AddHistoryOutput("STREAMWISE_LAMBDAL", "SWLambdaL", ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Exponential factor for iso-thermal temperature gradient");
   }
 
   AddAnalyzeSurfaceOutput(config);
@@ -258,10 +260,12 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
 
   LoadHistoryDataScalar(config, solver);
 
-  if(streamwisePeriodic) {
+  if (streamwisePeriodic) {
     SetHistoryOutputValue("STREAMWISE_MASSFLOW", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_MassFlow);
     SetHistoryOutputValue("STREAMWISE_DP", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_PressureDrop);
     SetHistoryOutputValue("STREAMWISE_HEAT", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_IntegratedHeatFlow);
+    if (config->GetnMarker_Isothermal() > 0)
+      SetHistoryOutputValue("STREAMWISE_LAMBDAL", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_LambdaL);
   }
 
   /*--- Set the analyse surface history values --- */
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index 5ffee0d23f5..5ba4b0dd1af 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -2399,6 +2399,7 @@ void CFlowOutput::WriteMetaData(const CConfig *config){
 
     if(config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
       meta_file << "STREAMWISE_PERIODIC_PRESSURE_DROP=" << GetHistoryFieldValue("STREAMWISE_DP") << endl;
+      meta_file << "STREAMWISE_PERIODIC_LAMBDAL=" << GetHistoryFieldValue("STREAMWISE_LAMBDAL") << endl;
     }
   }
 
diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
index 1b1989fed7d..c6304149235 100644
--- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
+++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
@@ -87,7 +87,7 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di
 
   /*--- Allocate extra solution variables, if any are in use. ---*/
 
-  const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config);
+  const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config); // TODO: Do we need this? I think no. Nitish
   nodes->AllocateAdjointSolutionExtra(nVarExtra);
 
   switch(KindDirect_Solver){
@@ -162,7 +162,7 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) {
   /*--- Boolean true indicates that an input is registered ---*/
   direct_solver->GetNodes()->RegisterSolution(true);
 
-  direct_solver->RegisterSolutionExtra(true, config);
+//   direct_solver->RegisterSolutionExtra(true, config);
 
   if (time_n_needed)
     direct_solver->GetNodes()->RegisterSolution_time_n();
@@ -261,6 +261,8 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo
     config->SetIncPressureOut_BC(BPressure);
     config->SetIncTemperature_BC(Temperature);
 
+    direct_solver->RegisterSolutionExtra(reset, config);
+
   }
 
   /*--- Register incompressible radiation values as input ---*/
@@ -299,7 +301,7 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) {
 
   direct_solver->GetNodes()->RegisterSolution(false);
 
-  direct_solver->RegisterSolutionExtra(false, config);
+//   direct_solver->RegisterSolutionExtra(false, config);
 }
 
 void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) {
diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp
index 64ddb898e74..c4a3a1e0712 100644
--- a/SU2_CFD/src/solvers/CHeatSolver.cpp
+++ b/SU2_CFD/src/solvers/CHeatSolver.cpp
@@ -706,6 +706,42 @@ void CHeatSolver::Heat_Fluxes(CGeometry *geometry, CSolver **solver_container, C
   Total_HeatFlux = AllBound_HeatFlux;
 }
 
+void CHeatSolver::Source_Residual(CGeometry *geometry,
+                      CSolver **solver_container,
+                      CNumerics **numerics_container,
+                      CConfig *config,
+                      unsigned short iMesh) {
+auto HeatGen = 0.0; //config->GetSolid_Heat_Generation();
+auto streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE);
+unsigned long nPointDomain = geometry->GetnPointDomain(), iPoint;
+unsigned short iVar;
+
+if (HeatGen != 0.0) {
+    /*--- loop over points ---*/
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+
+      /*--- Get control volume ---*/
+      su2double Volume = geometry->nodes->GetVolume(iPoint);
+
+      /*--- Get volumetric heat generation term and add to residual ---*/
+      for (iVar = 0; iVar < nVar; iVar++) {
+        LinSysRes(iPoint,iVar) -= Volume * HeatGen;
+      }
+    }
+    END_SU2_OMP_FOR
+  }
+if (streamwise_periodic) {
+    su2double lambda_l = 0.0, term1= 0.0, term2=0.0, grad_T=0.0, Volume=0.0, T = 0.0;
+    term1 = lambda_l * 2 * grad_T;
+    term2 = lambda_l * lambda_l * T;
+
+    for (iVar = 0; iVar < nVar; iVar++) {
+      LinSysRes(iPoint,iVar) -= Volume * (term1 - term2);
+    }
+  }
+}
+
 void CHeatSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config,
                                unsigned short iMesh, unsigned long Iteration) {
 
@@ -731,7 +767,7 @@ void CHeatSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container,
   } END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
   /*--- Loop domain points. ---*/
-
+ 
   SU2_OMP_FOR_DYN(omp_chunk_size)
   for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
 
@@ -1023,3 +1059,28 @@ void CHeatSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
 
   AD::EndNoSharedReading();
 }
+// su2double CHeatSolver::GetLambdaL(CGeometry *geometry, CSolver **solver_container, CConfig *config) {
+
+
+//       /*--- Loop over all heatflux Markers ---*/
+//     for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) {
+
+//       const su2double volume = geometry->nodes->GetVolume(iPoint);
+
+//       const su2double Temp = nodes->GetTemperature(iPoint);
+
+//       Volume_Local += volume;
+
+//       Volume_TempS_Local += volume * Temp;
+
+//       Volume_Temp_Local += volume * Temp * nodes->GetThermalConductivity(iPoint);
+
+//       // coeff_b1 for turbulence
+//       if (turbulent && (config->GetnMarker_Isothermal() != 0))
+//         turb_b1_coeff_Local += Temp * nodes->GetAuxVarGradient(iPoint, 0, 0) * config->GetSpecific_Heat_Cp() * volume  / config->GetPrandtl_Turb();
+
+//       Volume_VTemp_Local += volume * Temp * nodes->GetVelocity(iPoint, 2) * nodes->GetDensity(iPoint);
+
+//     } // points
+
+// }
\ No newline at end of file
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 5b81394bd9a..b4404f005ed 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -1693,6 +1693,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
     }
     else {
       numerics->SetStreamwisePeriodicValues(SPvals);
+      config->SetStreamwise_Periodic_ComputedMassFlow(SPvals.Streamwise_Periodic_MassFlow);
     }
 
     AD::StartNoSharedReading();
@@ -1714,6 +1715,9 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
       if(streamwise_periodic_temperature && turbulent)
         numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
 
+      /*--- Load the Prim variable gradient that we already computed. ---*/
+      numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr);
+
       /*--- Compute the streamwise periodic source residual and add to the total ---*/
       auto residual = numerics->ComputeResidual(config);
       LinSysRes.AddBlock(iPoint, residual);
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index d8327f1297a..0162620d67f 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -55,9 +55,11 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short
 
   /*--- Set the initial Streamwise periodic pressure drop value. ---*/
 
-  if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE)
+  if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE){
     // Note during restarts, the flow.meta is read first. But that sets the cfg-value so we are good here.
     SPvals.Streamwise_Periodic_PressureDrop = config->GetStreamwise_Periodic_PressureDrop();
+    SPvals.Streamwise_Periodic_LambdaL = config->GetStreamwise_Periodic_LamdaL();
+  }
 }
 
 void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh,
@@ -120,7 +122,7 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container
 void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
                                                      CConfig *config,
                                                      const unsigned short iMesh) {
-
+    const bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
   /*---------------------------------------------------------------------------------------------*/
   // 1. Evaluate massflow, area avg density & Temperature and Area at streamwise periodic outlet.
   // 2. Update delta_p is target massflow is chosen.
@@ -191,6 +193,9 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
   SPvals.Streamwise_Periodic_BoundaryArea = Area_Global;
   SPvals.Streamwise_Periodic_AvgDensity = Average_Density_Global;
 
+  su2double HeatFlow_Local = 0.0, HeatFlow_Global = 0.0;
+  su2double dTdn_Local = 0.0, dTdn_Global = 0.0;
+
   if (config->GetEnergy_Equation()) {
     /*---------------------------------------------------------------------------------------------*/
     /*--- 3. Compute the integrated Heatflow [W] for the energy equation source term, heatflux  ---*/
@@ -198,8 +203,6 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
     /*---    Here the Heatflux from all Boundary markers in the config-file is used.            ---*/
     /*---------------------------------------------------------------------------------------------*/
 
-    su2double HeatFlow_Local = 0.0, HeatFlow_Global = 0.0;
-
     /*--- Loop over all heatflux Markers ---*/
     for (auto iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
 
@@ -222,13 +225,99 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
           HeatFlow_Local += FaceArea * (-1.0) * Wall_HeatFlux/config->GetHeat_Flux_Ref();
         } // loop Vertices
       } // loop Heatflux marker
+
+      if (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) {
+        /*--- Identify the boundary by string name and retrive ISOTHERMAL from config ---*/
+
+        for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+
+          const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+
+          if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+          /*--- Compute wall heat flux (normal to the wall) based on computed temperature gradient ---*/
+          const auto AreaNormal = geometry->vertex[iMarker][iVertex]->GetNormal();
+
+          const auto GradT = nodes->GetGradient_Primitive(iPoint)[prim_idx.Temperature()];
+
+          dTdn_Local += nodes->GetThermalConductivity(iPoint) * GeometryToolbox::DotProduct(nDim, GradT, AreaNormal);
+        } // loop Vertices
+      } // loop Isothermal Marker
     } // loop AllMarker
 
     /*--- MPI Communication sum up integrated Heatflux from all processes ---*/
     SU2_MPI::Allreduce(&HeatFlow_Local, &HeatFlow_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&dTdn_Local, &dTdn_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    
+    su2double Volume_Temp_Local = 0.0, Volume_Temp_Global = 0.0;
+    su2double Volume_TempS_Local = 0.0, Volume_TempS_Global = 0.0;
+    su2double Volume_Local = 0.0, Volume_Global = 0.0;
+    su2double Volume_VTemp_Local = 0.0, Volume_VTemp_Global = 0.0;
+    su2double turb_b1_coeff_Local = 0.0, turb_b1_coeff_Global = 0.0;
+
+    /*--- Loop over all heatflux Markers ---*/
+    for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) {
+
+      const su2double volume = geometry->nodes->GetVolume(iPoint);
+
+      const su2double Temp = nodes->GetTemperature(iPoint);
+
+      Volume_Local += volume;
+
+      Volume_TempS_Local += volume * Temp;
+
+      Volume_Temp_Local += volume * Temp * nodes->GetThermalConductivity(iPoint);
+
+      // coeff_b1 for turbulence
+      if (turbulent && (config->GetnMarker_Isothermal() != 0)) {
+        su2double dot_product= 0.0;
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+          dot_product += config->GetPeriodic_Translation(0)[iDim]* nodes->GetAuxVarGradient(iPoint, 0, iDim);
+        }
+        // su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), nodes->GetAuxVarGradient(iPoint, 0));
+        turb_b1_coeff_Local += Temp * dot_product * config->GetSpecific_Heat_Cp() * volume  / config->GetPrandtl_Turb();
+      }
+        su2double dot_product_vel= 0.0;
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+          dot_product_vel += config->GetPeriodic_Translation(0)[iDim]* nodes->GetVelocity(iPoint, iDim);
+        }
+
+      Volume_VTemp_Local += volume * Temp * dot_product_vel * nodes->GetDensity(iPoint) * config->GetSpecific_Heat_Cp();
+
+    } // points
+
+    /*--- MPI Communication sum up integrated Heatflux from all processes ---*/
+    SU2_MPI::Allreduce(&Volume_Temp_Local, &Volume_Temp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_VTemp_Local, &Volume_VTemp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_Local, &Volume_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_TempS_Local, &Volume_TempS_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    if (turbulent && (config->GetnMarker_Isothermal() != 0))
+      SU2_MPI::Allreduce(&turb_b1_coeff_Local, &turb_b1_coeff_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
 
     /*--- Set the solver variable Integrated Heatflux ---*/
-    SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global;
+    if (config->GetnMarker_HeatFlux() > 0)
+      SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global;
+
+    if (config->GetnMarker_Isothermal() > 0) {
+      SPvals.Streamwise_Periodic_ThetaScaling = Volume_TempS_Global/Volume_Global;
+      
+      for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++)
+        nodes->SetStreamwise_Periodic_RecoveredTemperature(iPoint, nodes->GetTemperature(iPoint)/SPvals.Streamwise_Periodic_ThetaScaling);
+
+      /*--- Set the solver variable Lambda_L for iso-thermal BCs ---*/
+      const su2double b0_coeff =  Volume_Temp_Global; 
+      const su2double b1_coeff = Volume_VTemp_Global  + turb_b1_coeff_Global;
+      const su2double b2_coeff = -dTdn_Global;
+
+      /*--- Find the value of Lambda L by solving the quadratic equation ---*/
+      const su2double pred_lambda = (- b1_coeff + sqrt(b1_coeff * b1_coeff - 4 * b0_coeff * b2_coeff))/(2 * b0_coeff);
+      if (!config->GetDiscrete_Adjoint())
+        SPvals.Streamwise_Periodic_LambdaL -= 0.01 * (SPvals.Streamwise_Periodic_LambdaL - pred_lambda);
+      else
+        SPvals.Streamwise_Periodic_LambdaL = pred_lambda;
+        
+      config->SetStreamwise_Periodic_LamdaL(SPvals.Streamwise_Periodic_LambdaL);
+    } // if isothermal
   } // if energy
 }
 
@@ -236,9 +325,6 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
 void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, const CGeometry *geometry,
                                                                 const unsigned short iMesh) {
 
-  const bool energy = (config->GetEnergy_Equation() && config->GetStreamwise_Periodic_Temperature());
-  const auto InnerIter = config->GetInnerIter();
-
   /*--- Reference node on inlet periodic marker to compute relative distance along periodic translation vector. ---*/
   const auto ReferenceNode = geometry->GetStreamwise_Periodic_RefNode();
 
@@ -260,12 +346,12 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config,
     nodes->SetStreamwise_Periodic_RecoveredPressure(iPoint, Pressure_Recovered);
 
     /*--- InnerIter > 0 as otherwise MassFlow in the denominator would be zero ---*/
-    if (energy && InnerIter > 0) {
-      su2double Temperature_Recovered = nodes->GetTemperature(iPoint);
-      Temperature_Recovered += SPvals.Streamwise_Periodic_IntegratedHeatFlow /
-                              (SPvals.Streamwise_Periodic_MassFlow * nodes->GetSpecificHeatCp(iPoint) * norm2_translation) * dot_product;
-      nodes->SetStreamwise_Periodic_RecoveredTemperature(iPoint, Temperature_Recovered);
-    }
+    // if (energy && InnerIter > 0) {
+    //   su2double Temperature_Recovered = nodes->GetTemperature(iPoint);
+    //   Temperature_Recovered += SPvals.Streamwise_Periodic_IntegratedHeatFlow /
+    //                           (SPvals.Streamwise_Periodic_MassFlow * nodes->GetSpecificHeatCp(iPoint) * norm2_translation) * dot_product;
+    //   nodes->SetStreamwise_Periodic_RecoveredTemperature(iPoint, Temperature_Recovered);
+    // }
   } // for iPoint
   END_SU2_OMP_FOR
 
@@ -487,8 +573,8 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
 
       /*--- Apply a weak boundary condition for the energy equation.
       Compute the residual due to the prescribed heat flux. ---*/
-
-      LinSysRes(iPoint, nDim+1) -= thermal_conductivity*dTdn*Area;
+      // Test with relaxation factor with a factor of 0.5
+      LinSysRes(iPoint, nDim+1) -= 0.5*thermal_conductivity*dTdn*Area;
 
       /*--- Jacobian contribution for temperature equation. ---*/
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index de0ec914b73..41c3900a5ee 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -3329,6 +3329,7 @@ void CSolver::Read_SU2_Restart_Metadata(CGeometry *geometry, CConfig *config, bo
   su2double dCMy_dCL_ = config->GetdCMy_dCL();
   su2double dCMz_dCL_ = config->GetdCMz_dCL();
   su2double SPPressureDrop_ = config->GetStreamwise_Periodic_PressureDrop();
+  su2double SPLambdaL_ = config->GetStreamwise_Periodic_LamdaL();
   string::size_type position;
   unsigned long InnerIter_ = 0;
   ifstream restart_file;
@@ -3415,6 +3416,12 @@ void CSolver::Read_SU2_Restart_Metadata(CGeometry *geometry, CConfig *config, bo
         text_line.erase (0,34); SPPressureDrop_ = atof(text_line.c_str());
       }
 
+      position = text_line.find ("STREAMWISE_PERIODIC_LAMBDAL=",0);
+      if (position != string::npos) {
+              // Erase the name from the line, 'STREAMWISE_PERIODIC_LAMBDAL=' has 28 chars.
+        text_line.erase (0,28); SPLambdaL_ = atof(text_line.c_str());
+      }
+
     }
 
     /*--- Close the restart meta file. ---*/
@@ -3510,6 +3517,9 @@ void CSolver::Read_SU2_Restart_Metadata(CGeometry *geometry, CConfig *config, bo
     if ((config->GetStreamwise_Periodic_PressureDrop() != SPPressureDrop_) && (rank == MASTER_NODE))
       cout <<"WARNING: SU2 will use the STREAMWISE_PERIODIC_PRESSURE_DROP provided in the direct solution file: " << std::setprecision(16) << SPPressureDrop_ << endl;
     config->SetStreamwise_Periodic_PressureDrop(SPPressureDrop_);
+    if ((config->GetStreamwise_Periodic_LamdaL() != SPLambdaL_) && (rank == MASTER_NODE))
+      cout <<"WARNING: SU2 will use the STREAMWISE_PERIODIC_LAMBDAL provided in the direct solution file: " << std::setprecision(16) << SPLambdaL_ << endl;
+    config->SetStreamwise_Periodic_LamdaL(SPLambdaL_);
   }
   else {
     if ((config->GetStreamwise_Periodic_PressureDrop() != SPPressureDrop_) && (rank == MASTER_NODE))
diff --git a/externals/mel b/externals/mel
index 46205ab019e..2484cd3258e 160000
--- a/externals/mel
+++ b/externals/mel
@@ -1 +1 @@
-Subproject commit 46205ab019e5224559091375a6d71aabae6bc5b9
+Subproject commit 2484cd3258ef800a10e361016cb341834ee7930b
diff --git a/su2preconfig.timestamp b/su2preconfig.timestamp
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/subprojects/CoolProp b/subprojects/CoolProp
index bafdea1f39e..0ce42fcf3bb 160000
--- a/subprojects/CoolProp
+++ b/subprojects/CoolProp
@@ -1 +1 @@
-Subproject commit bafdea1f39ee873a6bb9833e3a21fe41f90b85e8
+Subproject commit 0ce42fcf3bb2c373512bc825a4f0c1973a78f307
diff --git a/subprojects/MLPCpp b/subprojects/MLPCpp
index c19c53ea2b8..665c45b7d35 160000
--- a/subprojects/MLPCpp
+++ b/subprojects/MLPCpp
@@ -1 +1 @@
-Subproject commit c19c53ea2b85ccfb185f1c6c87044dc0b5bc7ae0
+Subproject commit 665c45b7d3533c977eb1f637918d5b8b75c07d3b