diff --git a/src/lbm/lattice_model/CollisionModel.h b/src/lbm/lattice_model/CollisionModel.h
index bd521e042dedf772992dbf110d8751be30c95570..2aa4dd2c852809f1b0dc3aa40f8eeb76d2239bf3 100644
--- a/src/lbm/lattice_model/CollisionModel.h
+++ b/src/lbm/lattice_model/CollisionModel.h
@@ -109,6 +109,7 @@ public:
 
    real_t omega() const { return omega_; }
    inline real_t omega_bulk() const { return omega(); }
+   inline real_t omega_odd() const { return omega(); }
    real_t viscosity() const { return viscosity_; }
 
    real_t omega( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
@@ -274,6 +275,7 @@ public:
                  const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); }
 
    inline real_t omega_bulk() const { return omega(); }
+   inline real_t omega_odd() const { return lambda_d(); }
 
    static real_t lambda_e( const real_t _omega ) { return _omega; }
    static real_t lambda_d( const real_t _omega, const real_t _magicNumber = threeSixteenth )
@@ -484,6 +486,7 @@ public:
                  const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); }
 
    real_t omega_bulk() const { return s_[1]; }
+   real_t omega_odd() const { return s_[4]; }
 
    real_t viscosity() const { return viscosity_; }
    real_t viscosity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return viscosity_; }
diff --git a/src/lbm/lattice_model/ForceModel.h b/src/lbm/lattice_model/ForceModel.h
index 0c30441cae6c4bb32ebcb826c68909225e203863..365b64db1c56a16d1059e24749f839838267da3b 100644
--- a/src/lbm/lattice_model/ForceModel.h
+++ b/src/lbm/lattice_model/ForceModel.h
@@ -161,13 +161,13 @@ public:
    
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    { return DirectionIndependentTerms_T(); }
 
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & /*commonTerms*/, const real_t /*w*/,
-                     const real_t /*cx*/, const real_t /*cy*/, const real_t /*cz*/, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const { return real_t(0); }
+                     const real_t /*cx*/, const real_t /*cy*/, const real_t /*cz*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const { return real_t(0); }
 
    bool setConstantBodyForceIfPossible( const Vector3<real_t> &, const uint_t = uint_t(0) ) { return false; }
 };
@@ -214,13 +214,13 @@ public:
    
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    { return DirectionIndependentTerms_T(); }
 
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & /*commonTerms*/, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       return real_t(3.0) * w * ( cx * bodyForce_[0] + cy * bodyForce_[1] + cz * bodyForce_[2] );
    }
@@ -304,7 +304,7 @@ public:
 
    template< typename LatticeModel_T >
    real_t directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                     const Vector3<real_t> & /*velocity*/, const real_t rho, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                                     const Vector3<real_t> & /*velocity*/, const real_t rho, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       return DirectionIndependentTerm< LatticeModel_T >::get( rho );
    }
@@ -312,7 +312,7 @@ public:
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t rho,
                      const real_t commonTerm, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       const auto shiftedVelocity = ShiftedVelocity< LatticeModel_T >::get( velocity, force(x,y,z), commonTerm );
 
@@ -371,13 +371,13 @@ public:
    
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    { return DirectionIndependentTerms_T(); }
 
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & velocity, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & /*commonTerms*/, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       const Vector3<real_t> c( cx, cy, cz );
       return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ );
@@ -431,13 +431,13 @@ public:
    
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                                                          const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    { return DirectionIndependentTerms_T(); }
 
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & /*commonTerms*/, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       const Vector3<real_t> c( cx, cy, cz );
       return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) );
@@ -494,9 +494,9 @@ public:
    
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
-                                                          const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk ) const
+                                                          const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk, const real_t /*omega_odd*/ ) const
    {
-      if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value)
+      if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
       {
          const real_t one_over_d  = real_t(1) / real_t(LatticeModel_T::Stencil::D);
       
@@ -508,6 +508,7 @@ public:
       }
       else
       {
+         WALBERLA_ASSERT_FLOAT_EQUAL( omega, omega_bulk );
          return DirectionIndependentTerms_T();
       }
    }
@@ -515,16 +516,22 @@ public:
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & velocity, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & commonTerms, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t omega, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t omega, const real_t /*omega_bulk*/, const real_t omega_odd ) const
    {
       const Vector3<real_t> c( cx, cy, cz );
-      if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value)
+      if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
       {
          const real_t one_third  = real_t(1) / real_t(3);
       
          const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace();
          return real_t(3.0) * w * ( bodyForce_ * c + real_t(1.5) * common);
       }
+      else if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::TRT_tag >::value)
+      {
+         return real_t(3.0) * w * ( ( real_t(1) - real_t(0.5) * omega ) *
+                                    ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ ) +
+                                    ( omega - omega_odd ) * real_t(0.5) * (c * bodyForce_) );
+      }
       else
       {
          return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
@@ -580,9 +587,9 @@ public:
 
    template< typename LatticeModel_T >
    DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
-                                                          const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk ) const
+                                                          const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk, const real_t /*omega_odd*/ ) const
    {
-      if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value)
+      if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
       {
          const real_t one_over_d  = real_t(1) / real_t(LatticeModel_T::Stencil::D);
       
@@ -602,16 +609,22 @@ public:
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t /*rho*/,
                      const DirectionIndependentTerms_T & commonTerms, const real_t w,
-                     const real_t cx, const real_t cy, const real_t cz, const real_t omega, const real_t /*omega_bulk*/ ) const
+                     const real_t cx, const real_t cy, const real_t cz, const real_t omega, const real_t /*omega_bulk*/, const real_t omega_odd ) const
    {
       const Vector3<real_t> c( cx, cy, cz );
-      if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value)
+      if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
       {
          const real_t one_third  = real_t(1) / real_t(3);
       
          const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace();
          return real_t(3.0) * w * ( force(x,y,z) * c + real_t(1.5) * common);
       }
+      else if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::TRT_tag >::value)
+      {
+         return real_t(3.0) * w * ( ( real_t(1) - real_t(0.5) * omega ) *
+                                    ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) ) +
+                                    ( omega - omega_odd ) * real_t(0.5) * (c * force(x,y,z)) );
+      }
       else
       {
          return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
@@ -659,7 +672,7 @@ public:
    
    template< typename LatticeModel_T >
    Vector3<real_t> directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
-                                              const Vector3<real_t> & velocity, const real_t rho, const real_t omega, const real_t /*omega_bulk*/ ) const
+                                              const Vector3<real_t> & velocity, const real_t rho, const real_t omega, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       const auto rhoVelocity = rho * velocity;
       Vector3<real_t> & previousRhoVelocity = previousRhoVelocity_->get(x,y,z);
@@ -670,7 +683,7 @@ public:
 
    template< typename LatticeModel_T >
    real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/,
-                     const Vector3<real_t> & commonTerm, const real_t w, const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/ ) const
+                     const Vector3<real_t> & commonTerm, const real_t w, const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
    {
       return w * ( cx * commonTerm[0] + cy * commonTerm[1] + cz * commonTerm[2] );
    }
diff --git a/src/lbm/mrt/CellwiseSweep.impl.h b/src/lbm/mrt/CellwiseSweep.impl.h
index b9a3628ca16c3d82ee5bc1d1f28a487e6a48f240..fc1b50b2788deac48ad2a3204aaddf1c06056ce3 100644
--- a/src/lbm/mrt/CellwiseSweep.impl.h
+++ b/src/lbm/mrt/CellwiseSweep.impl.h
@@ -185,9 +185,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
          
          if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false)
          {
-            const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk() );
+            const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
             for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
-               dst->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk() );
+               dst->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
          }
       }
 
@@ -324,9 +324,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
          
          if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false)
          {
-            const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk() );
+            const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
             for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
-               src->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk() );
+               src->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
          }
       }
 
diff --git a/src/lbm/srt/CellwiseSweep.impl.h b/src/lbm/srt/CellwiseSweep.impl.h
index 37e156f57dcf1a150946db8c1ae96ad707bbee63..ec6d7fa4ba77200c5f8a99e590cd3ef4f4e9de02 100644
--- a/src/lbm/srt/CellwiseSweep.impl.h
+++ b/src/lbm/srt/CellwiseSweep.impl.h
@@ -1965,13 +1965,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
 
          const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
 
-         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega );
+         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega, omega );
 
          // collide
          for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
          {
             const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega );
+                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega, omega );
 
             dst->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * dst->get( x, y, z, d.toIdx() ) +
                                                              omega   * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
@@ -1996,12 +1996,12 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
 
          const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
 
-         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega );
+         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega, omega );
 
          for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
          {
             const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega );
+                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega, omega );
 
             src->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * src->get( x, y, z, d.toIdx() ) +
                                                              omega   * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
diff --git a/src/lbm/trt/CellwiseSweep.impl.h b/src/lbm/trt/CellwiseSweep.impl.h
index 241aeedfce15834f6d6cc8bac42fde42f27c5de4..4780803d80873db706e808b91452586f303ea41c 100644
--- a/src/lbm/trt/CellwiseSweep.impl.h
+++ b/src/lbm/trt/CellwiseSweep.impl.h
@@ -1630,13 +1630,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
 
          this->densityVelocityOut( x, y, z, lm, velocity, rho );
 
-         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e );
+         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e, lambda_d );
 
          // collide
          for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
          {
             const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e );
+                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e, lambda_d );
 
             const real_t fsym  = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho );
             const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
@@ -1681,7 +1681,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
 
          this->densityVelocityOut( x, y, z, lm, velocity, rho );
 
-         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e );
+         const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e, lambda_d );
 
          for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
             pdfs[ d.toIdx() ] = src->get( x, y, z, d.toIdx() );
@@ -1689,7 +1689,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
          for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
          {
             const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e );
+                                                                                           real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e, lambda_d );
 
             const real_t fsym  = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho );
             const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
diff --git a/src/pe_coupling/partially_saturated_cells_method/PSMSweep.h b/src/pe_coupling/partially_saturated_cells_method/PSMSweep.h
index 2b326f116aac44498e61010ae8a2987e8c212505..2c512c9b58bce7b4c2fa6b8c01b7b851d8d98e83 100644
--- a/src/pe_coupling/partially_saturated_cells_method/PSMSweep.h
+++ b/src/pe_coupling/partially_saturated_cells_method/PSMSweep.h
@@ -187,7 +187,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
 
                // possible external forcing on fluid
-               const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk() );
+               const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                // check if body is present
                if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
@@ -252,7 +252,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                   {
                      // external forcing
                      const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() );
+                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                      dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT
                                                       + omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
@@ -265,7 +265,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                   {
                      // external forcing
                      const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() );
+                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                      dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) + forceTerm;
                   }
@@ -340,7 +340,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
 
                // possible external forcing on fluid
-               const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk() );
+               const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                // check if a body is present
                if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
@@ -406,7 +406,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                   {
                      // external forcing
                      const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() );
+                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                      src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT
                                                       + omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
@@ -419,7 +419,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
                   {
                      // external forcing
                      const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
-                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() );
+                                                                                                    real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
 
                      src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] )  + forceTerm;
                   }