diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index 918f35748..cf01af3f6 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -1409,6 +1409,34 @@ public:
         if (enableAquifers_)
             aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
 
+        // Add source term from deck
+        const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
+        if (source.size() > 0) {
+            std::array<int,3> ijk;
+                this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
+            RateVector massRate(0.0);
+            if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = source.rate({ijk, SourceComponent::OIL}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = source.rate({ijk, SourceComponent::GAS}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = source.rate({ijk, SourceComponent::WATER}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if constexpr (enableSolvent) {
+                massRate[Indices::solventSaturationIdx] = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
+            rate.setMassRate(massRate, pvtRegionIdx);
+            if constexpr (enablePolymer) {
+                rate[Indices::polymerConcentrationIdx] = source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if constexpr (enableEnergy) {
+                rate[Indices::contiEnergyEqIdx] = source.hrate(ijk) / this->model().dofTotalVolume(globalDofIdx);
+            }
+        }
+
         // if requested, compensate systematic mass loss for cells which were "well
         // behaved" in the last time step
         // Note that we don't allow for drift compensation if there are no active wells.
diff --git a/tests/test_ParallelSerialization.cpp b/tests/test_ParallelSerialization.cpp
index 1b90eb07a..f5e8d9c8f 100644
--- a/tests/test_ParallelSerialization.cpp
+++ b/tests/test_ParallelSerialization.cpp
@@ -52,6 +52,7 @@
 #include <opm/input/eclipse/EclipseState/IOConfig/IOConfig.hpp>
 #include <opm/input/eclipse/Schedule/RPTConfig.hpp>
 #include <opm/input/eclipse/Schedule/RSTConfig.hpp>
+#include <opm/input/eclipse/Schedule/Source.hpp>
 #include <opm/input/eclipse/Schedule/SummaryState.hpp>
 #include <opm/input/eclipse/Schedule/Action/ActionResult.hpp>
 #include <opm/input/eclipse/Schedule/Action/State.hpp>
@@ -281,6 +282,7 @@ TEST_FOR_TYPE(SkprpolyTable)
 TEST_FOR_TYPE(SkprwatTable)
 TEST_FOR_TYPE(SICD)
 TEST_FOR_TYPE(SolventDensityTable)
+TEST_FOR_TYPE(Source)
 TEST_FOR_TYPE(SummaryConfig)
 TEST_FOR_TYPE(SummaryConfigNode)
 TEST_FOR_TYPE(SummaryState)