diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp
index 216d40088e9b3ae0ee14b4accd20eec3eaf4d3a3..d59676359c3b208ca797a2cd73c707a638e096c8 100644
--- a/src/scheme/Order2AcousticSolver.cpp
+++ b/src/scheme/Order2AcousticSolver.cpp
@@ -511,8 +511,31 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
         const std::shared_ptr<const DiscreteFunctionVariant>& p,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
+    std::shared_ptr m_app = getCommonMesh({rho, c, u, p});
+    std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
+    std::shared_ptr<const DiscreteFunctionVariant> u_app = u;
+    std::shared_ptr<const DiscreteFunctionVariant> E_app = E;
+
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
-    return apply_fluxes(dt, rho, u, E, ur, Fjr);
+    std::tie(m_app,rho_app,u_app,E_app) = apply_fluxes(dt/2, rho, u, E, ur, Fjr);
+
+    const double& gamma = 3;
+
+    const DiscreteScalarFunction& rho_d = rho_app->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u_d   = u_app->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& E_d   = E_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& eps_d = E_d - 0.5 * dot(u_d, u_d);
+    const DiscreteScalarFunction& p_d   = (gamma - 1) * rho_d * eps_d;
+    const DiscreteScalarFunction& c_d   = sqrt(gamma * p_d / rho_d);
+
+    const std::shared_ptr<const DiscreteFunctionVariant>& p_app =
+      std::make_shared<const DiscreteFunctionVariant>(p_d);
+    const std::shared_ptr<const DiscreteFunctionVariant>& c_app = 
+      std::make_shared<const DiscreteFunctionVariant>(c_d);
+
+    auto [ur_app, Fjr_app] = compute_fluxes(solver_type,rho_app,c_app,u_app,p_app,bc_descriptor_list);
+    return apply_fluxes(dt, rho, u, E, ur_app, Fjr_app);
+
   }
 
   Order2AcousticSolver()                 = default;