diff --git a/src/fluid/calcRightHandSide.hpp b/src/fluid/calcRightHandSide.hpp
index c8738b4b..38ef76d5 100644
--- a/src/fluid/calcRightHandSide.hpp
+++ b/src/fluid/calcRightHandSide.hpp
@@ -149,47 +149,49 @@ struct Fluid_CorrectFluxFunctor {
         Flux(MX1+meanDir,k,j,i) += meanV * Flux(RHO,k,j,i);
       } // Fargo & Rotation corrections
 
-      real Ax = A(k,j,i);
+      //////////////////////////////////////////////
+      // Define correction factor for the fluxes
+      //////////////////////////////////////////////
 
-      for(int nv = 0 ; nv < Phys::nvar ; nv++) {
-        Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax;
-      }
+      real Ax[Phys::nvar]; // corrected Area
+      for(int nv = 0 ; nv < Phys::nvar ; nv++) Ax[nv] = A(k,j,i);
 
       // Curvature terms
-#if    (GEOMETRY == POLAR       && COMPONENTS >= 2) \
-    || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3)
-      if constexpr (dir==IDIR) {
-        // Conserve angular momentum, hence flux is R*Bphi
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i));
-        if constexpr(Phys::mhd) {
-          if(Ax<SMALL_NUMBER) Ax=SMALL_NUMBER;    //avoid singularity around poles
-          // No area for this one
-          Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) / Ax;
+      #if    (GEOMETRY == POLAR       && COMPONENTS >= 2) \
+          || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3)
+        if constexpr (dir==IDIR) {
+          // Conserve angular momentum, hence flux is R*Bphi
+          Ax[iMPHI] *= FABS(x1m(i));
+          if constexpr(Phys::mhd) {
+            Ax[iBPHI] = 1.0;    // corrected Flux is simply Bphi
+          }
         }
-      }
-#endif // GEOMETRY==POLAR OR CYLINDRICAL
+      #endif // GEOMETRY==POLAR OR CYLINDRICAL
 
-#if GEOMETRY == SPHERICAL
-      if constexpr(dir==IDIR) {
-  #if COMPONENTS == 3
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i));
-  #endif // COMPONENTS == 3
-        if constexpr(Phys::mhd) {
-          if(Ax<SMALL_NUMBER) Ax=SMALL_NUMBER;    // avoid singularity around poles
-          EXPAND(                                            ,
-              Flux(iBTH,k,j,i)  = Flux(iBTH,k,j,i) * x1m(i) / Ax;  ,
-              Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i) * x1m(i) / Ax; )
-        }
-      } else if constexpr (dir==JDIR) {
-  #if COMPONENTS == 3
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sinx2m(j));
-        if constexpr(Phys::mhd) {
-          if(Ax<SMALL_NUMBER) Ax=SMALL_NUMBER;    // avoid singularity around poles
-          Flux(iBPHI,k,j,i) = Flux(iBPHI,k,j,i)  / Ax;
+      #if GEOMETRY == SPHERICAL
+        if constexpr(dir==IDIR) {
+          #if COMPONENTS == 3
+            Ax[iMPHI] *= FABS(x1m(i));
+          #endif // COMPONENTS == 3
+          if constexpr(Phys::mhd) {
+            EXPAND(                                            ,
+                Ax[iBTH]  = x1m(i);  ,
+                Ax[iBPHI]  = x1m(i); )
+          }// Phys::mhd
+        } else if constexpr (dir==JDIR) {
+          #if COMPONENTS == 3
+              Ax[iMPHI] *= FABS(sinx2m(j));
+              if constexpr(Phys::mhd) {
+                Ax[iBPHI] = 1.0; // corrected Flux is simply Bphi
+              }
+          #endif // COMPONENTS = 3
         }
-  #endif // COMPONENTS = 3
+      #endif // GEOMETRY == SPHERICAL
+
+      // Finally correct the flux
+      for(int nv = 0 ; nv < Phys::nvar ; nv++) {
+        Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax[nv];
       }
-#endif // GEOMETRY == SPHERICAL
     }
 };
 
diff --git a/src/rkl/rkl.hpp b/src/rkl/rkl.hpp
index a2f6734f..c8dd8cb9 100644
--- a/src/rkl/rkl.hpp
+++ b/src/rkl/rkl.hpp
@@ -762,32 +762,42 @@ void RKLegendre<Phys>::CalcParabolicRHS(real t) {
              data->beg[IDIR],data->end[IDIR]+ioffset,
     KOKKOS_LAMBDA (int n, int k, int j, int i) {
       real Ax = A(k,j,i);
-
-#if GEOMETRY != CARTESIAN
-      if(Ax<SMALL_NUMBER)
-        Ax=SMALL_NUMBER;    // Essentially to avoid singularity around poles
-#endif
-
       const int nv = varList(n);
 
-      Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax;
-
       // Curvature terms
-#if    (GEOMETRY == POLAR       && COMPONENTS >= 2) \
-    || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3)
-      if(dir==IDIR && nv==iMPHI) {
-        // Conserve angular momentum, hence flux is R*Vphi
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i));
-      }
-#endif // GEOMETRY==POLAR OR CYLINDRICAL
+      #if    (GEOMETRY == POLAR       && COMPONENTS >= 2) \
+      || (GEOMETRY == CYLINDRICAL && COMPONENTS == 3)
+        if(dir==IDIR) {
+          if(nv==iMPHI) {
+            Ax *= FABS(x1m(i));
+          }
+          if constexpr(Phys::mhd) {
+            if(nv==iBPHI) {
+              // No area for this one
+              Ax = 1;
+            }
+          }
+        }
+      #endif // GEOMETRY==POLAR OR CYLINDRICAL
+
+      #if GEOMETRY == SPHERICAL
+        if(dir==IDIR && nv==VX3) {
+          Ax *= FABS(x1m(i));
+        } else if(dir==JDIR && nv==VX3) {
+          Ax *= FABS(sm(j));
+        }
+
+        if constexpr(Phys::mhd) {
+          if(dir == IDIR  && (nv==BX3 || nv == BX2)) {
+            Ax = x1m(i);
+          }
+          if(dir==JDIR && nv==BX3) {
+            Ax = 1.0;
+          }
+        }
+      #endif // GEOMETRY == SPHERICAL
 
-#if GEOMETRY == SPHERICAL && COMPONENTS == 3
-      if(dir==IDIR && nv==iMPHI) {
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(x1m(i));
-      } else if(dir==JDIR && nv==iMPHI) {
-        Flux(iMPHI,k,j,i) = Flux(iMPHI,k,j,i) * FABS(sm(j));
-      }
-#endif // GEOMETRY == SPHERICAL && COMPONENTS == 3
+      Flux(nv,k,j,i) = Flux(nv,k,j,i) * Ax;
     }
   );
 
@@ -818,17 +828,42 @@ void RKLegendre<Phys>::CalcParabolicRHS(real t) {
       }
 
 #if GEOMETRY != CARTESIAN
-  #ifdef iMPHI
-      if((dir==IDIR) && (nv == iMPHI)) {
+  if(dir==IDIR) {
+    #ifdef iMPHI
+      if((nv == iMPHI)) {
         rhs /= x1(i);
       }
+    #endif // iMPHI
+    real dx_ = dx(i);
+    real x1_ = x1(i);
+
+    if constexpr(Phys::mhd) {
+      #if (GEOMETRY == POLAR || GEOMETRY == CYLINDRICAL) &&  (defined iBPHI)
+        if(nv==iBPHI) rhs = - 1 / dx_ * (Flux(iBPHI, k, j, i+1) - Flux(iBPHI, k, j, i) );
+
+      #elif (GEOMETRY == SPHERICAL)
+        real q = 1 / (x1_*dx_);
+        if(nv == BX2 || nv == BX3) {
+          rhs = -q * ((Flux(nv, k, j, i+1)  - Flux(nv, k, j, i) ));
+        }
+      #endif // GEOMETRY
+    } // MHD
+  } // dir==IDIR
+  if(dir==JDIR) {
     #if (GEOMETRY == SPHERICAL) && (COMPONENTS == 3)
-      if((dir==JDIR) && (nv == iMPHI)) {
+      if(nv == iMPHI) {
         rhs /= FABS(s(j));
       }
+      real dx_ = dx(j);
+      real rt_ = rt(i);
+      if constexpr(Phys::mhd) {
+        if(nv == iBPHI) {
+          rhs = - 1 / (rt_*dx_) * (Flux(nv, k, j+1, i) - Flux(nv, k, j, i));
+        }
+      }
     #endif // GEOMETRY
-      // Nothing for KDIR
-  #endif  // iMPHI
+  } // dir==JDIR
+  // Nothing for KDIR
 #endif // GEOMETRY != CARTESIAN
 
       // store the field components