diff --git a/src/olbProcessors/AdeConcBoundaries2D.h b/src/olbProcessors/AdeConcBoundaries2D.h
new file mode 100644
index 0000000..8e5e754
--- /dev/null
+++ b/src/olbProcessors/AdeConcBoundaries2D.h
@@ -0,0 +1,57 @@
+#pragma once
+
+#include <olb2D.h>
+#include <olb2D.hh>
+
+namespace arch {
+
+// Forward declared dependencies
+template<typename T>
+class Opening;
+
+template <typename T>
+class Node;
+
+}
+
+namespace olb {
+
+template<typename T>
+class AdeConcBC : public AnalyticalF2D<T,T> {
+
+private:
+  T width;                ///< Physical width of the opening.
+  T xc;                   ///< The physical x coordinate, xc, of the opening node (center).
+  T yc;                   ///< The physical y coordinate, yc, of the opening node (center).
+  T x0;                   ///< The physical x coordinate, x0, at which the opening "starts" (bottom for right facing opening).
+  T y0;                   ///< The physical y coordinate, y0, at which the opening "starts" (bottom for right facing opening).
+  std::vector<T> array;   ///< Vector that contains the array of concentration values along the opening.
+
+public:
+
+  AdeConcBC(const arch::Opening<T>& opening, std::vector<T> concentrationField);
+
+  bool operator()(T output[], const T input[]) override;
+
+};
+
+template<typename T>
+class AdeConc1D : public SuperPlaneIntegralF2D<T,T> {
+
+private:
+  T width;                ///< Physical width of the opening.
+  T xc;                   ///< The physical x coordinate, xc, of the opening node (center).
+  T yc;                   ///< The physical y coordinate, yc, of the opening node (center).
+  T x0;                   ///< The physical x coordinate, x0, at which the opening "starts" (bottom for right facing opening).
+  T y0;                   ///< The physical y coordinate, y0, at which the opening "starts" (bottom for right facing opening).
+  std::vector<T> array;   ///< Vector that contains the array of concentration values along the opening.
+
+public:
+
+  AdeConc1D(const arch::Opening<T>& opening, std::vector<T> concentrationField);
+
+  bool operator()(T output[], const T input[]) override;
+
+};
+
+}
diff --git a/src/olbProcessors/AdeConcBoundaries2D.hh b/src/olbProcessors/AdeConcBoundaries2D.hh
new file mode 100644
index 0000000..71eef5a
--- /dev/null
+++ b/src/olbProcessors/AdeConcBoundaries2D.hh
@@ -0,0 +1,37 @@
+#include <AdeConcBoundary2D.h>
+
+namespace olb {
+
+template<typename T>
+AdeConcBoundary2D<T>::AdeConcBoundary2D(const arch::Opening<T>& opening, std::vector<T> concentrationField) :
+  width(opening.width),
+  xc(opening.node->getPosition()[0]),
+  yc(opening.node->getPosition()[0]),
+  x0(xc - opening.tangent[0]*0.5*width),
+  y0(yc - opening.tangent[1]*0.5*width),
+  array(concentrationField) { }
+
+template<typename T>
+bool AdeConcBoundary2D<T>::operator()(T output[], const T input[]) {
+  T s   = std::sqrt((input[0]-x0)*(input[0]-x0) + (input[1]-y0)*(input[1]-y0));
+  T d   = std::sqrt((input[0]-xc)*(input[0]-xc) + (input[1]-yc)*(input[1]-yc));
+  T res = width / array.size();
+  if (s < width) {
+    output[0] = array[std::floor(s/res)];
+  }
+  if (0.5*width - d < 0.0) {
+    output[0] = T();
+  }
+  return true;
+}
+
+template<typename T>
+AdeConc1D<T>::AdeConc1D(const arch::Opening<T>& opening, std::vector<T> concentrationField) 
+  : SuperPlaneIntegralF2D<T>()
+
+template<typename T>
+bool AdeConc1D<T>::operator()(T output[], const T input[]) {
+
+}
+
+} // namespace olb
\ No newline at end of file
diff --git a/src/simulation/CFDSim.hh b/src/simulation/CFDSim.hh
index 79b546c..77ec6e3 100644
--- a/src/simulation/CFDSim.hh
+++ b/src/simulation/CFDSim.hh
@@ -47,7 +47,7 @@ namespace sim {
     }
 
     template<typename T>
-    bool conductADSimulation(const std::unordered_map<int, std::unique_ptr<CFDSimulator<T>>>& cfdSimulators) {
+    bool conductADSimulation(const std::unordered_map<int, std::unique_ptr<CFDSimulator<T>>>& cfdSimulators, bool fieldValues) {
 
         bool allConverge = true;
 
@@ -61,7 +61,7 @@ namespace sim {
             assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::ESS_LBM);
             throw std::runtime_error("Simulation of Advection Diffusion not defined for ESS LBM.");
             #endif
-            cfdSimulator.second->adSolve();
+            cfdSimulator.second->adSolve(fieldValues);
 
             if (!cfdSimulator.second->hasAdConverged()) {
                 allConverge = false;
diff --git a/src/simulation/MixingModels.hh b/src/simulation/MixingModels.hh
index 01fbb44..cea9fcf 100644
--- a/src/simulation/MixingModels.hh
+++ b/src/simulation/MixingModels.hh
@@ -469,14 +469,18 @@ void DiffusionMixingModel<T>::generateInflows(T timeStep, arch::Network<T>* netw
             for (auto& [channelId, channel] : network->getChannels()) { // There is only one channel connected to a node at a CFD hybrid node connection
                 if (cfdSimulator->getFlowRates().at(nodeId) > 0.0 && (channel->getNodeA() == nodeId || channel->getNodeB() == nodeId)) {
                     std::unordered_map<int, std::tuple<std::function<T(T)>,std::vector<T>,T>> newDistributions;
-                    std::cout << "Getting here, nodeId is " + std::to_string(nodeId) << std::endl;
+                    std::cout << "Getting here 1, nodeId is " + std::to_string(nodeId) << std::endl;
                     for (const auto& [specieId, concentrationField] : concentrationFieldsIn.at(nodeId)) {
-                        std::cout << "Getting here 2" << std::endl;
+                        std::cout << "Getting here 1.1" << std::endl;
                         T dx = channel->getWidth() / concentrationField.size();
+                        std::cout << "Getting here 1.2" << std::endl;
                         T pecletNr = (std::abs(channel->getFlowRate()) / channel->getHeight()) / (sim->getSpecie(specieId))->getDiffusivity();
+                        std::cout << "Getting here 1.3" << std::endl;
                         std::tuple<std::function<T(T)>, std::vector<T>, T> analyticalResult = getAnalyticalSolutionHybrid(
                             channel->getLength(), channel->getFlowRate(), channel->getWidth(), noFourierTerms, pecletNr, concentrationField, dx);
+                        std::cout << "Getting here 1.4" << std::endl;
                         newDistributions.try_emplace(specieId, analyticalResult);
+                        std::cout << "Getting here 1.5" << std::endl;
                     }
                     //Create new DiffusiveMixture
                     DiffusiveMixture<T>* newMixture = dynamic_cast<DiffusiveMixture<T>*>(sim->addDiffusiveMixture(newDistributions));
@@ -484,9 +488,11 @@ void DiffusionMixingModel<T>::generateInflows(T timeStep, arch::Network<T>* netw
                     this->injectMixtureInEdge(newMixture->getId(), channelId);
                 }
                 else if (cfdSimulator->getFlowRates().at(nodeId) < 0.0 && (channel->getNodeA() == nodeId || channel->getNodeB() == nodeId)) {
+                    std::cout << "Getting here 2, nodeId is " + std::to_string(nodeId) << std::endl;
                     int resolutionIntoCfd = cfdSimulator->getResolution(nodeId);
                     concentrationFieldsOut.try_emplace(nodeId, std::unordered_map<int, std::vector<T>>());
                     for (const auto& [specieId, speciePtr] : sim->getSpecies()) {
+                        std::cout << "Getting here 2.1" << std::endl;
                         std::vector<T> concentrationFieldForCfdInput = defineConcentrationNodeFieldForCfdInput(resolutionIntoCfd, specieId, channelId, channel->getWidth(), mixtures, noFourierTerms);
                         concentrationFieldsOut.at(nodeId).try_emplace(specieId, concentrationFieldForCfdInput);
                     }                    
diff --git a/src/simulation/Simulation.hh b/src/simulation/Simulation.hh
index e02aaaa..a971911 100644
--- a/src/simulation/Simulation.hh
+++ b/src/simulation/Simulation.hh
@@ -639,7 +639,7 @@ namespace sim {
             // Obtain overal steady-state concentration results
             bool concentrationConverged = false;
             while (!concentrationConverged) {
-                concentrationConverged = conductADSimulation(cfdSimulators);
+                concentrationConverged = conductADSimulation(cfdSimulators, mixingModel->isDiffusive());
                 this->mixingModel->propagateSpecies(network, this);
             }
         }
@@ -682,7 +682,7 @@ namespace sim {
             // Obtain overal steady-state concentration results
             bool concentrationConverged = false;
             while (!concentrationConverged) {
-                concentrationConverged = conductADSimulation(cfdSimulators);
+                concentrationConverged = conductADSimulation(cfdSimulators, mixingModel->isDiffusive());
                 this->mixingModel->propagateSpecies(network, this);
             }
         }
diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h
index 47ca506..bccf7b2 100755
--- a/src/simulation/simulators/cfdSimulator.h
+++ b/src/simulation/simulators/cfdSimulator.h
@@ -256,7 +256,7 @@ class CFDSimulator {
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    virtual void setBoundaryValues(int iT) = 0;
+    virtual void setBoundaryValues(int iT, bool fieldValues=false) = 0;
 
     /**
      * @brief Prepare the LBM geometry of this simulator.
@@ -285,7 +285,7 @@ class CFDSimulator {
     /**
      * @brief Conducts the collide and stream operations of the AD lattice(s).
     */
-    virtual void adSolve() 
+    virtual void adSolve(bool fieldValues) 
     {
         throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
     }
diff --git a/src/simulation/simulators/essContinuous.h b/src/simulation/simulators/essContinuous.h
index eb500b1..5a7bf7d 100644
--- a/src/simulation/simulators/essContinuous.h
+++ b/src/simulation/simulators/essContinuous.h
@@ -83,7 +83,7 @@ namespace sim {
              * @brief Set the boundary values on the lattice at the module nodes.
              * @param[in] iT Iteration step.
             */
-            void setBoundaryValues(int iT);
+            void setBoundaryValues(int iT, bool fieldValues=false);
 
             /**
              * @brief Conducts the collide and stream operations of the lattice.
diff --git a/src/simulation/simulators/essContinuous.hh b/src/simulation/simulators/essContinuous.hh
index 873b8bc..9d363f3 100644
--- a/src/simulation/simulators/essContinuous.hh
+++ b/src/simulation/simulators/essContinuous.hh
@@ -17,7 +17,7 @@ namespace sim{
     }
 
     template<typename T>
-    void essLbmSimulator<T>::setBoundaryValues(int iT)
+    void essLbmSimulator<T>::setBoundaryValues(int iT, bool fieldValues)
     {
         setFlowRates(flowRates);
         setPressures(pressures);
diff --git a/src/simulation/simulators/essMixing.h b/src/simulation/simulators/essMixing.h
index eb500b1..5a7bf7d 100644
--- a/src/simulation/simulators/essMixing.h
+++ b/src/simulation/simulators/essMixing.h
@@ -83,7 +83,7 @@ namespace sim {
              * @brief Set the boundary values on the lattice at the module nodes.
              * @param[in] iT Iteration step.
             */
-            void setBoundaryValues(int iT);
+            void setBoundaryValues(int iT, bool fieldValues=false);
 
             /**
              * @brief Conducts the collide and stream operations of the lattice.
diff --git a/src/simulation/simulators/essMixing.hh b/src/simulation/simulators/essMixing.hh
index 28c0319..e012bae 100644
--- a/src/simulation/simulators/essMixing.hh
+++ b/src/simulation/simulators/essMixing.hh
@@ -18,7 +18,7 @@ namespace sim{
     }
 
     template<typename T>
-    void essLbmSimulator<T>::setBoundaryValues(int iT)
+    void essLbmSimulator<T>::setBoundaryValues(int iT, bool fieldValues)
     {
         setFlowRates(flowRates);
         setPressures(pressures);
diff --git a/src/simulation/simulators/olbContinuous.h b/src/simulation/simulators/olbContinuous.h
index 9ff6572..29bd03e 100644
--- a/src/simulation/simulators/olbContinuous.h
+++ b/src/simulation/simulators/olbContinuous.h
@@ -174,7 +174,7 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    void setBoundaryValues(int iT) override;
+    void setBoundaryValues(int iT, bool fieldValues=false) override;
 
     /**
      * @brief Conducts the collide and stream operations of the lattice.
diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh
index 5d49fcb..a1e05ac 100644
--- a/src/simulation/simulators/olbContinuous.hh
+++ b/src/simulation/simulators/olbContinuous.hh
@@ -74,7 +74,7 @@ void lbmSimulator<T>::prepareLattice () {
 }
 
 template<typename T>
-void lbmSimulator<T>::setBoundaryValues (int iT) {
+void lbmSimulator<T>::setBoundaryValues (int iT, bool fieldValues) {
 
     for (auto& [key, Opening] : this->moduleOpenings) {
         if (this->groundNodes.at(key)) {
@@ -389,15 +389,15 @@ template<typename T>
 void lbmSimulator<T>::setFlowProfile2D (int openingKey, T openingWidth)  {
     T maxVelocity = (3./2.)*(flowRates[openingKey]/(openingWidth));
     T distance2Wall = getConverter().getConversionFactorLength()/2.;
-    this->flowProfiles.at(openingKey) = std::make_shared<olb::Poiseuille2D<T>>(getGeometry(), openingKey+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall);
-    getLattice().defineU(getGeometry(), openingKey+3, *this->flowProfiles.at(openingKey));
+    flowProfiles.at(openingKey) = std::make_shared<olb::Poiseuille2D<T>>(getGeometry(), openingKey+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall);
+    getLattice().defineU(getGeometry(), openingKey+3, *flowProfiles.at(openingKey));
 }
 
 template<typename T>
 void lbmSimulator<T>::setPressure2D (int openingKey)  {
     T rhoV = getConverter().getLatticeDensityFromPhysPressure((pressures[openingKey]));
-    this->densities.at(openingKey) = std::make_shared<olb::AnalyticalConst2D<T,T>>(rhoV);
-    getLattice().defineRho(getGeometry(), openingKey+3, *this->densities.at(openingKey));
+    densities.at(openingKey) = std::make_shared<olb::AnalyticalConst2D<T,T>>(rhoV);
+    getLattice().defineRho(getGeometry(), openingKey+3, *densities.at(openingKey));
 }
 
 template<typename T>
diff --git a/src/simulation/simulators/olbMixing.h b/src/simulation/simulators/olbMixing.h
index 82d9400..f021515 100644
--- a/src/simulation/simulators/olbMixing.h
+++ b/src/simulation/simulators/olbMixing.h
@@ -44,8 +44,8 @@ using ADDynamics = olb::AdvectionDiffusionBGKdynamics<T,ADDESCRIPTOR>;
 using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
 
 protected:
-    std::unordered_map<int, std::unordered_map<int, T>> concentrations;   ///< Vector of concentration values at module nodes. <nodeId, <speciId, conc>>
-    std::unordered_map<int, std::unordered_map<int, std::vector<T>>> nodeConcentrationFields;   ///< Vector of concentration values at module nodes. <nodeId, <speciId, <conc>>>
+    std::unordered_map<int, std::unordered_map<int, T>> concentrations;   ///< Vector of concentration values at module nodes. <nodeId, <specieId, conc>>
+    std::unordered_map<int, std::unordered_map<int, std::vector<T>>> nodeConcentrationFields;   ///< Vector of concentration values at module nodes. <nodeId, <specieId, <conc>>>
 
     std::unordered_map<int, Specie<T>*> species;
 
@@ -61,7 +61,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     std::unordered_map<int, T*> fluxWall;
     T zeroFlux = 0.0;
 
-    std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::AnalyticalConst2D<T,T>>>> concentrationProfiles;
+    std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::AdeConcBoundary2D<T,T>>>> concentrationProfiles;   ///< Map of concentration field in olb format at module nodes. <nodeId, <specieId, <conc>>>
     std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::SuperPlaneIntegralFluxPressure2D<T>>>> meanConcentrations;       ///< Map of mean pressure values at module nodes.
 
     auto& getAdConverter(int key) {
@@ -94,13 +94,15 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     */
     void executeCoupling() override;
 
-    void setConcentration2D(int key);
+    void setConcentrationConst2D(int key);
+
+    void setConcentrationField2D(int key);
 
     /**
      * @brief Update the values at the module nodes based on the simulation result after stepIter iterations.
      * @param[in] iT Iteration step.
     */
-    void storeCfdResults(int iT);
+    void storeCfdResults(int iT, bool fieldValues);
 
 public:
     /**
@@ -159,7 +161,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    void setBoundaryValues(int iT) override;
+    void setBoundaryValues(int iT, bool fieldValues) override;
 
     /**
      * @brief Conducts the collide and stream operations of the lattice.
@@ -174,7 +176,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     /**
      * @brief Conducts the collide and stream operations of the AD lattice(s).
     */
-    void adSolve() override;
+    void adSolve(bool fieldValues) override;
 
     /**
      * @brief Write the vtk file with results of the CFD simulation to file system.
diff --git a/src/simulation/simulators/olbMixing.hh b/src/simulation/simulators/olbMixing.hh
index 7e6cb0c..d48b638 100644
--- a/src/simulation/simulators/olbMixing.hh
+++ b/src/simulation/simulators/olbMixing.hh
@@ -84,7 +84,7 @@ void lbmMixingSimulator<T>::prepareLattice () {
 }
 
 template<typename T>
-void lbmMixingSimulator<T>::setBoundaryValues (int iT) {
+void lbmMixingSimulator<T>::setBoundaryValues (int iT, bool fieldValues) {
 
     for (auto& [key, Opening] : this->moduleOpenings) {
         if (this->groundNodes.at(key)) {
@@ -92,12 +92,19 @@ void lbmMixingSimulator<T>::setBoundaryValues (int iT) {
         } else {
             this->setPressure2D(key);
         }
-        setConcentration2D(key);
+        // Constant concentration values
+        if (!fieldValues) {
+            setConcentrationConst2D(key);
+        }
+        // Constant field values
+        if (fieldValues) {
+            setConcentrationField2D(key);
+        }
     }    
 }
 
 template<typename T>
-void lbmMixingSimulator<T>::storeCfdResults (int iT) {
+void lbmMixingSimulator<T>::storeCfdResults (int iT, bool fieldValues) {
     int input[1] = { };
     T output[10];
     
@@ -115,13 +122,22 @@ void lbmMixingSimulator<T>::storeCfdResults (int iT) {
     for (auto& [key, Opening] : this->moduleOpenings) {
         // If the node is an outflow, write the concentration value
         if (this->flowRates.at(key) < 0.0) {
-            for (auto& [speciesId, adLattice] : adLattices) {
-                this->meanConcentrations.at(key).at(speciesId)->operator()(output,input);
-                this->concentrations.at(key).at(speciesId) = output[0];
-                if (iT % 1000 == 0) {
-                    this->meanConcentrations.at(key).at(speciesId)->print();
+            // Constant concentration values
+            if (!fieldValues) {
+                for (auto& [speciesId, adLattice] : adLattices) {
+                    this->meanConcentrations.at(key).at(speciesId)->operator()(output,input);
+                    this->concentrations.at(key).at(speciesId) = output[0];
+                    if (iT % 1000 == 0) {
+                        this->meanConcentrations.at(key).at(speciesId)->print();
+                    }
                 }
             }
+            // Field concentration values
+            if (fieldValues) {
+                /** TODO:
+                 * 
+                 */
+            }
         }
     }
     
@@ -226,9 +242,9 @@ void lbmMixingSimulator<T>::nsSolve() {
 }
 
 template<typename T>
-void lbmMixingSimulator<T>::adSolve() {
+void lbmMixingSimulator<T>::adSolve(bool fieldValues) {
     // theta = 10
-    this->setBoundaryValues(this->step);
+    this->setBoundaryValues(this->step, fieldValues);
     for (int iT = 0; iT < 100; ++iT){
         for (auto& [speciesId, adLattice] : adLattices) {
             adLattice->collideAndStream();
@@ -236,7 +252,7 @@ void lbmMixingSimulator<T>::adSolve() {
         writeVTK(this->step);
         this->step += 1;
     }
-    storeCfdResults(this->step);
+    storeCfdResults(this->step, fieldValues);
 }
 
 template<typename T>
@@ -370,7 +386,7 @@ void lbmMixingSimulator<T>::prepareCoupling() {
 }
 
 template<typename T>
-void lbmMixingSimulator<T>::setConcentration2D (int key) {
+void lbmMixingSimulator<T>::setConcentrationConst2D (int key) {
     // Set the boundary concentrations for inflows and outflows
     // If the boundary is an inflow
     if (this->flowRates.at(key) >= 0.0) {
@@ -393,6 +409,30 @@ void lbmMixingSimulator<T>::setConcentration2D (int key) {
     }
 }
 
+template<typename T>
+void lbmMixingSimulator<T>::setConcentrationField2D (int key) {
+
+    // Set the boundary concentrations for inflows and outflows
+    // If the boundary is an inflow
+    if (this->flowRates.at(key) >= 0.0) {
+        for (auto& [speciesId, adLattice] : adLattices) {
+            olb::setAdvectionDiffusionTemperatureBoundary<T,ADDESCRIPTOR>(*adLattice, getAdConverter(speciesId).getLatticeRelaxationFrequency(), this->getGeometry(), key+3);
+            adLattice->defineRho(this->getGeometry(), key+3, *concentrationProfiles.at(key).at(speciesId));
+        }
+    }
+    // If the boundary is an outflow
+    else if (this->flowRates.at(key) < 0.0) {
+        for (auto& [speciesId, adLattice] : adLattices) {
+            olb::setRegularizedHeatFluxBoundary<T,ADDESCRIPTOR>(*adLattice, getAdConverter(speciesId).getLatticeRelaxationFrequency(), this->getGeometry(), key+3, &zeroFlux);
+        }
+    }
+
+    else {
+        std::cerr << "Error: Invalid Flow Type Boundary Node." << std::endl;
+        exit(1);
+    }
+}
+
 template<typename T>
 void lbmMixingSimulator<T>::storeConcentrations(std::unordered_map<int, std::unordered_map<int, T>> concentrations_) {
     this->concentrations = concentrations_;
@@ -400,6 +440,11 @@ void lbmMixingSimulator<T>::storeConcentrations(std::unordered_map<int, std::uno
 
 template<typename T>
 void lbmMixingSimulator<T>::storeNodeConcentrationFields(std::unordered_map<int, std::unordered_map<int, std::vector<T>>> concentrationFieldsOut_) {
+    for (auto& [nodeId, fields] : concentrationFieldsOut_) {
+        for (auto& [specieId, concField] : fields) {
+            concentrationProfiles.at(nodeId).at(specieId) = std::make_shared<olb::AdeConcBoundary2D>(this->moduleOpenings.at(nodeId), concentrationFieldsOut_.at(nodeId).at(specieId));
+        }
+    }
     this->nodeConcentrationFields = concentrationFieldsOut_;
 }