diff --git a/src/fluid/RiemannSolver/riemannSolver.hpp b/src/fluid/RiemannSolver/riemannSolver.hpp
index abb70eed..b6c1e6f1 100644
--- a/src/fluid/RiemannSolver/riemannSolver.hpp
+++ b/src/fluid/RiemannSolver/riemannSolver.hpp
@@ -70,7 +70,6 @@ class RiemannSolver {
 
   IdefixArray4D<real> Vc;
   IdefixArray4D<real> Vs;
-  IdefixArray4D<real> Flux;
   IdefixArray3D<real> cMax;
   Fluid<Phys>* hydro;
   DataBlock *data;
@@ -90,7 +89,6 @@ class RiemannSolver {
 template <typename Phys>
 RiemannSolver<Phys>::RiemannSolver(Input &input, Fluid<Phys>* hydro) : Vc{hydro->Vc},
                                       Vs{hydro->Vs},
-                                      Flux{hydro->FluxRiemann},
                                       cMax{hydro->cMax},
                                       hydro{hydro},
                                       data{hydro->data}
@@ -148,8 +146,6 @@ RiemannSolver<Phys>::RiemannSolver(Input &input, Fluid<Phys>* hydro) : Vc{hydro-
     mySolver = HLL_DUST;
   }
 
-
-
   // Shock flattening
   this->haveShockFlattening = input.CheckEntry(std::string(Phys::prefix),"shockFlattening")>=0;
   // Init shock flattening
diff --git a/src/fluid/addNonIdealMHDFlux.hpp b/src/fluid/addNonIdealMHDFlux.hpp
index 52c281d2..9ae6af30 100644
--- a/src/fluid/addNonIdealMHDFlux.hpp
+++ b/src/fluid/addNonIdealMHDFlux.hpp
@@ -19,8 +19,7 @@ void Fluid<Phys>::AddNonIdealMHDFlux(const real t) {
 
   if constexpr(Phys::mhd) {
     int ioffset,joffset,koffset;
-
-    IdefixArray4D<real> Flux = this->FluxRiemann;
+    IdefixArray4D<real> Flux = this->FluxRiemann[dir];
     IdefixArray4D<real> Vc   = this->Vc;
     IdefixArray4D<real> Vs   = this->Vs;
     IdefixArray3D<real> dMax = this->dMax;
diff --git a/src/fluid/calcParabolicFlux.hpp b/src/fluid/calcParabolicFlux.hpp
index a8c6761c..8713f910 100644
--- a/src/fluid/calcParabolicFlux.hpp
+++ b/src/fluid/calcParabolicFlux.hpp
@@ -42,7 +42,7 @@ void Fluid<Phys>::CalcParabolicFlux(const real t) {
     if(data->haveFargo && viscosityStatus.isExplicit) {
       data->fargo->AddVelocityFluid(t,this);
     }
-    this->viscosity->AddViscousFlux(dir,t, this->FluxRiemann);
+    this->viscosity->AddViscousFlux(dir,t, this->FluxRiemann[dir]);
 
     // Remove back Fargo velocity
     if(data->haveFargo && viscosityStatus.isExplicit) {
@@ -53,18 +53,18 @@ void Fluid<Phys>::CalcParabolicFlux(const real t) {
   // Add thermal diffusion
   if( (thermalDiffusionStatus.isExplicit && (!data->rklCycle))
     || (thermalDiffusionStatus.isRKL && data->rklCycle))  {
-    this->thermalDiffusion->AddDiffusiveFlux(dir,t, this->FluxRiemann);
+    this->thermalDiffusion->AddDiffusiveFlux(dir,t, this->FluxRiemann[dir]);
   }
 
   if( (bragViscosityStatus.isExplicit && (!data->rklCycle))
     || (bragViscosityStatus.isRKL && data->rklCycle))  {
-    this->bragViscosity->AddBragViscousFlux(dir,t, this->FluxRiemann);
+    this->bragViscosity->AddBragViscousFlux(dir,t, this->FluxRiemann[dir]);
   }
 
   // Add braginskii thermal diffusion
   if( (bragThermalDiffusionStatus.isExplicit && (!data->rklCycle))
     || (bragThermalDiffusionStatus.isRKL && data->rklCycle))  {
-    this->bragThermalDiffusion->AddBragDiffusiveFlux(dir,t, this->FluxRiemann);
+    this->bragThermalDiffusion->AddBragDiffusiveFlux(dir,t, this->FluxRiemann[dir]);
   }
 
   idfx::popRegion();
diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp
index 38ef76d5..357fcaac 100644
--- a/src/fluid/calcRightHandSide.hpp
+++ b/src/fluid/calcRightHandSide.hpp
@@ -21,7 +21,7 @@ struct Fluid_CorrectFluxFunctor {
   explicit Fluid_CorrectFluxFunctor (Fluid<Phys> *hydro, real dt) {
     Uc   = hydro->Uc;
     Vc   = hydro->Vc;
-    Flux = hydro->FluxRiemann;
+    Flux = hydro->FluxRiemann[dir];
     A    = hydro->data->A[dir];
     dV   = hydro->data->dV;
     x1m  = hydro->data->xl[IDIR];
@@ -205,7 +205,7 @@ struct Fluid_CalcRHSFunctor {
   explicit Fluid_CalcRHSFunctor (Fluid<Phys> *hydro, real dt) {
     Uc   = hydro->Uc;
     Vc   = hydro->Vc;
-    Flux = hydro->FluxRiemann;
+    Flux = hydro->FluxRiemann[dir];
     A    = hydro->data->A[dir];
     dV   = hydro->data->dV;
     x1m  = hydro->data->xl[IDIR];
diff --git a/src/fluid/evolveStage.hpp b/src/fluid/evolveStage.hpp
index 31dae5b7..b0a467ba 100644
--- a/src/fluid/evolveStage.hpp
+++ b/src/fluid/evolveStage.hpp
@@ -14,20 +14,20 @@ template<typename Phys>
 template<int dir>
 void Fluid<Phys>::LoopDir(const real t, const real dt) {
     // Step 2: compute the intercell flux with our Riemann solver, store the resulting InvDt
-    this->rSolver->template CalcFlux<dir>(this->FluxRiemann);
+    this->rSolver->template CalcFlux<dir>(this->FluxRiemann[dir]);
 
     // Step 2.5: compute intercell parabolic flux when needed
     if(haveExplicitParabolicTerms) CalcParabolicFlux<dir>(t);
 
     // If we have tracers, compute the tracer intercell flux
     if(haveTracer) {
-      this->tracer->template CalcFlux<dir, Phys>(this->FluxRiemann);
+      this->tracer->template CalcFlux<dir, Phys>(this->FluxRiemann[dir]);
     }
 
     // Step 3: compute the resulting evolution of the conserved variables, stored in Uc
     CalcRightHandSide<dir>(t,dt);
     if(haveTracer) {
-      this->tracer->template CalcRightHandSide<dir, Phys>(this->FluxRiemann,t ,dt);
+      this->tracer->template CalcRightHandSide<dir, Phys>(this->FluxRiemann[dir],t ,dt);
     }
 
     // Recursive: do next dimension
diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp
index 44cace6b..9e0615de 100644
--- a/src/fluid/fluid.hpp
+++ b/src/fluid/fluid.hpp
@@ -67,7 +67,7 @@ class Fluid {
   void EvolveStage(const real, const real);
   void ResetStage();
   void ShowConfig();
-  IdefixArray4D<real> GetFlux() {return this->FluxRiemann;}
+  IdefixArray4D<real> GetFlux(int dir) {return this->FluxRiemann.at(dir);}
   int CheckNan();
 
   // Our boundary conditions
@@ -182,7 +182,7 @@ class Fluid {
   // Required by time integrator
   IdefixArray3D<real> InvDt;
 
-  IdefixArray4D<real> FluxRiemann;
+  std::array<IdefixArray4D<real>,DIMENSIONS> FluxRiemann;
   IdefixArray3D<real> dMax;    // Maximum diffusion speed
 
   std::unique_ptr<RiemannSolver<Phys>> rSolver;
@@ -539,8 +539,12 @@ Fluid<Phys>::Fluid(Grid &grid, Input &input, DataBlock *datain, int n) {
                               data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]);
   dMax = IdefixArray3D<real>(prefix+"_dMax",
                               data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]);
-  FluxRiemann =  IdefixArray4D<real>(prefix+"_FluxRiemann", Phys::nvar+nTracer,
-                                     data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]);
+
+  for(int i = 0 ; i < DIMENSIONS ; i++) {
+    FluxRiemann[i] = IdefixArray4D<real>(prefix+"_FluxRiemann_X"+std::to_string(i),
+                                        Phys::nvar+nTracer,
+                                        data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]);
+  }
 
   if constexpr(Phys::mhd) {
     Vs = IdefixArray4D<real>(prefix+"_Vs", DIMENSIONS,
diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp
index c8dd8cb9..a2636eda 100644
--- a/src/rkl/rkl.hpp
+++ b/src/rkl/rkl.hpp
@@ -30,7 +30,7 @@ class RKLegendre {
   RKLegendre(Input &, Fluid<Phys>*);
   void Cycle();
   void ResetStage();
-  void ResetFlux();
+  void ResetFlux(int dir);
   void EvolveStage(real);
   template <int> void CalcParabolicRHS(real);
   void ComputeDt();
@@ -555,9 +555,9 @@ void RKLegendre<Phys>::Cycle() {
 
 
 template<typename Phys>
-void RKLegendre<Phys>::ResetFlux() {
+void RKLegendre<Phys>::ResetFlux(int dir) {
   idfx::pushRegion("RKLegendre::ResetFlux");
-  IdefixArray4D<real> Flux = hydro->FluxRiemann;
+  IdefixArray4D<real> Flux = hydro->FluxRiemann[dir];
   IdefixArray1D<int> vars = this->varList;
   idefix_for("RKL_ResetFlux",
              0,nvarRKL,
@@ -576,7 +576,6 @@ template<typename Phys>
 struct RKLegendre_ResetStageFunctor {
   explicit RKLegendre_ResetStageFunctor(RKLegendre<Phys> *rkl) {
     dU = rkl->dU;
-    Flux = rkl->hydro->FluxRiemann;
     vars = rkl->varList;
     stage = rkl->stage;
     nvar = rkl->nvarRKL;
@@ -596,7 +595,6 @@ struct RKLegendre_ResetStageFunctor {
   }
 
   IdefixArray4D<real> dU;
-  IdefixArray4D<real> Flux;
   IdefixArray1D<int> vars;
   IdefixArray4D<real> dA, dB;
   IdefixArray3D<real> ex,ey,ez;
@@ -680,7 +678,7 @@ void RKLegendre<Phys>::ComputeDt() {
 template<typename Phys>
 template<int dir>
 void RKLegendre<Phys>::LoopDir(real t) {
-    ResetFlux();
+    ResetFlux(dir);
 
     // CalcParabolicFlux
     hydro->template CalcParabolicFlux<dir>(t);
@@ -724,7 +722,7 @@ template <int dir>
 void RKLegendre<Phys>::CalcParabolicRHS(real t) {
   idfx::pushRegion("RKLegendre::CalcParabolicRHS");
 
-  IdefixArray4D<real> Flux = hydro->FluxRiemann;
+  IdefixArray4D<real> Flux = hydro->FluxRiemann[dir];
   IdefixArray3D<real> A    = data->A[dir];
   IdefixArray3D<real> dV   = data->dV;
   IdefixArray1D<real> x1m  = data->xl[IDIR];
diff --git a/test/MHD/AmbipolarWind/setup.cpp b/test/MHD/AmbipolarWind/setup.cpp
index 78620f6e..bde4ada9 100644
--- a/test/MHD/AmbipolarWind/setup.cpp
+++ b/test/MHD/AmbipolarWind/setup.cpp
@@ -335,7 +335,7 @@ void EmfBoundary(DataBlock& data, const real t) {
 }
 
 void FluxBoundary(DataBlock & data, int dir, BoundarySide side, const real t) {
-    IdefixArray4D<real> Flux = data.hydro->FluxRiemann;
+    IdefixArray4D<real> Flux = data.hydro->FluxRiemann[dir];
     if( dir==IDIR && side == left) {
         int iref = data.beg[IDIR];
 
diff --git a/test/SelfGravity/UniformCollapse/setup.cpp b/test/SelfGravity/UniformCollapse/setup.cpp
index 6f773a15..22f4a34d 100644
--- a/test/SelfGravity/UniformCollapse/setup.cpp
+++ b/test/SelfGravity/UniformCollapse/setup.cpp
@@ -65,7 +65,7 @@ void FluxBoundary(Fluid<DefaultPhysics> *hydro, int dir, BoundarySide side, cons
     if((dir==IDIR) && (side == left)) {
       // Loading needed data
       DataBlock &data = *hydro->data;
-      IdefixArray4D<real> Flux = hydro->FluxRiemann;
+      IdefixArray4D<real> Flux = hydro->FluxRiemann[dir];
       real halfDt = data.dt/2.; // RK2, dt is actually half at each flux calculation
       int iref = data.nghost[IDIR];
       real rin = data.xbeg[IDIR];