diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f3a8f16..936b50c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -88,7 +88,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().relaxSources(td.cloud().cloudCopy());
         }
@@ -99,7 +99,7 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        if (solution_.coupled())
+        if (solution_.couplingRegime() > 1.5)
         {
             td.cloud().scaleSources();
         }
@@ -165,7 +165,7 @@ template<class CloudType>
 template<class TrackData>
 void Foam::KinematicCloud<CloudType>::evolveCloud(TrackData& td)
 {
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         td.cloud().resetSourceTerms();
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 72cad26..2dbe962 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -555,7 +555,7 @@ Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
             << max(UCoeff()).value() << endl;
     }
 
-    if (solution_.coupled())
+    if (solution_.couplingRegime() > 1.5)
     {
         if (solution_.semiImplicit("U"))
         {
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index fd7bf31..7fab1df 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
@@ -43,6 +43,7 @@ Foam::cloudSolution::cloudSolution
     iter_(1),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(1.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(true),
@@ -69,6 +70,7 @@ Foam::cloudSolution::cloudSolution
     iter_(cs.iter_),
     trackTime_(cs.trackTime_),
     coupled_(cs.coupled_),
+    couplingRegime_(cs.couplingRegime_),
     cellValueSourceCorrection_(cs.cellValueSourceCorrection_),
     maxTrackTime_(cs.maxTrackTime_),
     resetSourcesOnStartup_(cs.resetSourcesOnStartup_),
@@ -90,6 +92,7 @@ Foam::cloudSolution::cloudSolution
     iter_(0),
     trackTime_(0.0),
     coupled_(false),
+    couplingRegime_(0.0),
     cellValueSourceCorrection_(false),
     maxTrackTime_(0.0),
     resetSourcesOnStartup_(false),
@@ -126,6 +129,8 @@ void Foam::cloudSolution::read()
 
     if (coupled_)
     {
+        couplingRegime_ = 2.0;
+
         const dictionary&
             schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
index e155abc..5423258 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H
@@ -88,6 +88,14 @@ class cloudSolution
             //  carrier phase
             Switch coupled_;
 
+            //- Coupling regime of parcels and carrier phase:
+            //  0-way coupling: Only non-coupled forces are calculated
+            //  1-way coupling: Coupled and non-coupled forces are calculated,
+            //                  but no contribution to the source term
+            //  2-way coupling: Coupled and non-coupled forces are calculated
+            //                  and changes are contributed to the source term
+            scalar couplingRegime_;
+
             //- Flag to correct cell values with latest transfer information
             //  during the lagrangian timestep
             Switch cellValueSourceCorrection_;
@@ -174,8 +182,11 @@ public:
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
-            //- Return non-const access to the coupled flag
-            inline Switch& coupled();
+            //- Return const access to the coupling regime
+            inline scalar couplingRegime() const;
+
+            //- Return non-const access to the coupling regime
+            inline scalar& couplingRegime();
 
             //- Return const access to the cell value correction flag
             inline const Switch cellValueSourceCorrection() const;
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
index 6280b97..f2f5d1b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H
@@ -101,15 +101,21 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
-inline Foam::Switch& Foam::cloudSolution::coupled()
+inline const Foam::Switch Foam::cloudSolution::coupled() const
 {
     return coupled_;
 }
 
 
-inline const Foam::Switch Foam::cloudSolution::coupled() const
+inline Foam::scalar& Foam::cloudSolution::couplingRegime()
 {
-    return coupled_;
+    return couplingRegime_;
+}
+
+
+inline Foam::scalar Foam::cloudSolution::couplingRegime() const
+{
+    return couplingRegime_;
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index 7abbc9f..6e399ed 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -105,7 +105,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
     volScalarField& Yi
 ) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("Yi"))
         {
@@ -180,7 +180,7 @@ Foam::ReactingCloud<CloudType>::Srho(const label i) const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& rhoi = tRhoi();
         rhoi = rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
@@ -217,7 +217,7 @@ Foam::ReactingCloud<CloudType>::Srho() const
         )
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         scalarField& sourceField = trhoTrans();
         forAll(rhoTrans_, i)
@@ -236,7 +236,7 @@ template<class CloudType>
 inline Foam::tmp<Foam::fvScalarMatrix>
 Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
 {
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         tmp<volScalarField> trhoTrans
         (
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index e97b6cf..9f1600b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -51,7 +51,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         this->subModelProperties().lookup("radiation") >> radiation_;
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index ea1b79f..0b9a45a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -249,7 +249,7 @@ Foam::ThermoCloud<CloudType>::Sh(volScalarField& hs) const
             << max(hsCoeff()).value() << endl;
     }
 
-    if (this->solution().coupled())
+    if (this->solution().couplingRegime() > 1.5)
     {
         if (this->solution().semiImplicit("h"))
         {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 5d9e4ef..d64b433 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -132,14 +132,11 @@ void Foam::KinematicParcel<ParcelType>::calc
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
     // Update momentum transfer
     td.cloud().UTrans()[cellI] += np0*dUTrans;
 
     // Update momentum transfer coefficient
     td.cloud().UCoeff()[cellI] += np0*Spu;
-    }
 }
 
 
@@ -166,7 +163,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     // Momentum source due to particle forces
     const parcelType& p = static_cast<const parcelType&>(*this);
-    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    forceSuSp Fcp(vector::zero, 0.0);
+    if (td.cloud().solution().couplingRegime() > 0.5)
+    {
+        Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
+    }
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
     const forceSuSp Feff = Fcp + Fncp;
     const scalar massEff = forces.massEff(p, mass);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 72ad743..78e9548 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -345,8 +345,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
-        {
         scalar dm = np0*mass0;
 
         // Absorb parcel into carrier phase
@@ -373,7 +371,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
 
         td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
-        }
 
         return;
     }
@@ -420,8 +417,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
-    {
     // Transfer mass lost to carrier mass, momentum and enthalpy sources
     forAll(YGas_, i)
     {
@@ -479,7 +474,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 47a1310..96ad1ba 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -375,8 +375,6 @@ void Foam::ReactingParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        if (td.cloud().solution().coupled())
-        {
         scalar dm = np0*mass0;
 
         // Absorb parcel into carrier phase
@@ -392,7 +390,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().UTrans()[cellI] += dm*U0;
 
         td.cloud().phaseChange().addToPhaseChangeMass(np0*mass1);
-        }
 
         return;
     }
@@ -438,8 +435,6 @@ void Foam::ReactingParcel<ParcelType>::calc
     // 4. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (td.cloud().solution().coupled())
-    {
     // Transfer mass lost to carrier mass, momentum and enthalpy sources
     forAll(dMass, i)
     {
@@ -469,7 +464,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 74ece70..e98f928 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -237,8 +237,6 @@ void Foam::ThermoParcel<ParcelType>::calc
 
     //  Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
-    {
     // Update momentum transfer
     td.cloud().UTrans()[cellI] += np0*dUTrans;
 
@@ -260,7 +258,6 @@ void Foam::ThermoParcel<ParcelType>::calc
         td.cloud().radT4()[cellI] += dt*np0*T4;
         td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
     }
-    }
 }
 
 
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index 22d9811..ed96db0 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -68,13 +68,13 @@ void Foam::SprayParcel<ParcelType>::calc
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
-    const bool coupled = td.cloud().solution().coupled();
+    const scalar coupled = td.cloud().solution().couplingRegime();
 
     // check if parcel belongs to liquid core
     if (liquidCore() > 0.5)
     {
-        // liquid core parcels should not interact with the gas
-        td.cloud().solution().coupled() = false;
+        // no drag calculation for liquid core parcels
+        td.cloud().solution().couplingRegime() = 0.0;
     }
 
     // get old mixture composition
@@ -131,7 +131,7 @@ void Foam::SprayParcel<ParcelType>::calc
     }
 
     // restore coupled
-    td.cloud().solution().coupled() = coupled;
+    td.cloud().solution().couplingRegime() = coupled;
 }