diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
index c485a848..86cf2b21 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
@@ -103,12 +103,16 @@ species::species
     partMolarConc_(NULL),
     loopCounter_(-1),
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
+    couplingTimestep_(0.0),
     massSourceCurr_(0.0),
     massSourceTot_(0.0),
     initialized_(false)
 {
     particleCloud_.checkCG(false);
     allocateMyArrays();
+    scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
+    scalar dtCFD = mesh_.time().deltaTValue();
+    couplingTimestep_ = max(dtDEM,dtCFD);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -300,7 +304,6 @@ void species::execute()
 
     // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
     {
-        scalar timestep = mesh_.time().deltaTValue();
         changeOfGasMassField_.primitiveFieldRef() = 0.0;
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
@@ -322,17 +325,17 @@ void species::execute()
 
             // take care for implementation in LIGGGHTS: species produced from particles defined positive
             // changeOf...Fields need to be mass per volume per timestep
-            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep);
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * couplingTimestep_);
             changeOfSpeciesMassFields_[i].correctBoundaryConditions();
             changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
 
             if (verbose_)
             {
                 Info << "total conversion of species" << speciesNames_[i] << " = "
-                     << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep) << endl;
+                     << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * couplingTimestep_) << endl;
             }
         }
-        massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * timestep);
+        massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * couplingTimestep_);
         massSourceTot_ += massSourceCurr_;
 
         if (verbose_)
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
index 2b9529a0..fa411506 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H
@@ -111,6 +111,8 @@ private:
 
     label Nevery_;
 
+    scalar couplingTimestep_;
+
     scalar massSourceCurr_;
 
     scalar massSourceTot_;
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
index a5f3e2e1..3b0cfe0c 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C
@@ -63,9 +63,15 @@ reactionHeat::reactionHeat
         ),
         mesh_,
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0),0.0)
-    )
+    ),
+    loopCounter_(-1),
+    Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
+    couplingTimestep_(0.0)
 {
     allocateMyArrays();
+    scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
+    scalar dtCFD = mesh_.time().deltaTValue();
+    couplingTimestep_ = max(dtDEM,dtCFD);
 
     if(propsDict_.found("maxsource"))
     {
@@ -94,6 +100,12 @@ void reactionHeat::allocateMyArrays() const
 
 void reactionHeat::calcEnergyContribution()
 {
+    loopCounter_++;
+    if (loopCounter_ % Nevery_ != 0)
+    {
+        return;
+    }
+
    // realloc the arrays
     allocateMyArrays();
 
@@ -118,7 +130,7 @@ void reactionHeat::calcEnergyContribution()
         NULL
     );
 
-    reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V());
+    reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V() * Nevery_ * couplingTimestep_);
 
     forAll(reactionHeatField_,cellI)
     {
@@ -130,6 +142,12 @@ void reactionHeat::calcEnergyContribution()
         }
     }
 
+    if (verbose_)
+    {
+        Info << "reaction heat per unit time = "
+                     << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
+    }
+
     reactionHeatField_.correctBoundaryConditions();
 }
 
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
index 4ebd7a02..a827fb20 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H
@@ -60,6 +60,12 @@ protected:
 
     volScalarField reactionHeatField_;
 
+    label loopCounter_;
+
+    label Nevery_;
+
+    scalar couplingTimestep_;
+
     void allocateMyArrays() const;
 
 public: