diff --git a/src/projects/GRHD/cons2prim_implementations.jl b/src/projects/GRHD/cons2prim_implementations.jl
index 6c83868733f63f18f14b51520c6e8b5a8d0a781c..12825c4170993705bca9bcec24bf9a20934bbbf4 100644
--- a/src/projects/GRHD/cons2prim_implementations.jl
+++ b/src/projects/GRHD/cons2prim_implementations.jl
@@ -18,8 +18,11 @@
   #   vr   = Formulae.vr(B, D, Sr, tau, p)
   #   eps  = Formulae.ϵ(B, D, Sr, tau, p)
   # end
+
+  D, Sr, tau = conservative_fixing(D, Sr, tau, B, equation)
   rho, vr, eps, p = cons2prim_kastaun((D, Sr, tau, p, r, B), equation)
-  @returns rho, vr, eps, p
+  # rho, vr, eps, p, D, Sr, tau = atmosphering(rho, vr, eps, p, D, Sr, tau)
+  @returns rho, vr, eps, p, D, Sr, tau
 end
 
 
@@ -153,8 +156,8 @@ end
   function fn_What(mu)
     _vhat = min(mu * s, v0)
     if _vhat < 0
-      println((; q, _vhat, s, mu, v0))
       error()
+      println((; q, _vhat, s, mu, v0))
     end
     return 1 / sqrt(1 - _vhat^2)
   end
@@ -218,7 +221,8 @@ end
       # enforce EOS bounds
       _epsmin, _epsmax = fn_epsmin(rhohat), fn_epsmax(rhohat)
       if rhohat < rhomin || epshat < _epsmin
-        rhohat, epshat, What = rhomin, _epsmin
+        # TODO Should we reset What here?
+        rhohat, epshat = rhomin, _epsmin
       elseif rhohat > rhomax || epshat > _epsmax
         rhohat, epshat = rhomax, _epsmax
       end
@@ -276,3 +280,40 @@ end
 
   @returns rho, vr, eps, p
 end
+
+
+#######################################################################
+#                    Conservative variable fixing                     #
+#######################################################################
+
+
+function conservative_fixing(D, Sr, tau, B, equation)
+
+  @unpack eos, atmosphere = equation
+  rhomin = atmosphere.atm_density * atmosphere.threshold
+  rhoatm = atmosphere.atm_density
+
+  if D < rhomin
+    D = rhoatm
+  end
+
+  if tau < 0
+    tau = 1e-12
+  end
+
+  # (48)-(50) in Deppe paper where we put B=0
+  # that = tau / D
+  # muhat = 0
+  # Solution to (52) when B=0
+  What = 1 + tau / D
+
+  # in spherical symmetry: S^2 = S_r^2 / B^2
+  # Hence, Smax^2/S^2 = Srmax^2 / Sr^2
+  S2 = Sr^2 / B^2
+  factor = sqrt( (1-1e-12) * (What^2-1) * D^2 / (S2 + D^2 * 1e-16) )
+  if factor > 1
+    Sr /= factor * (1 + 1e-8)
+  end
+
+  return D, Sr, tau
+end
diff --git a/src/projects/GRHD/equation.jl b/src/projects/GRHD/equation.jl
index 6cc35897c43f4d6c5b041984606aed741068200d..d56760f1af37f005fec691658c032113f9341dab 100644
--- a/src/projects/GRHD/equation.jl
+++ b/src/projects/GRHD/equation.jl
@@ -118,10 +118,16 @@ end
 
 @with_signature function flux(equation::Equation)
   @accepts State(sD, sSr, stau), p, rho, vr, eps, A, B, r
-  rAB3 = r * A * B^3
-  flx_sD   = rAB3 * Formulae.fD(B, rho, vr, eps, p)
-  flx_sSr  = rAB3 * Formulae.fSr(B, rho, vr, eps, p)
-  flx_stau = rAB3 * Formulae.fτ(B, rho, vr, eps, p)
+  rAB3     = r * A * B^3
+  flx_sD, flx_sSr, flx_stau = 0, 0, 0
+  try
+    flx_sD   = rAB3 * Formulae.fD(B, rho, vr, eps, p)
+    flx_sSr  = rAB3 * Formulae.fSr(B, rho, vr, eps, p)
+    flx_stau = rAB3 * Formulae.fτ(B, rho, vr, eps, p)
+  catch e
+    println("failed @ r = $r")
+    rethrow(e)
+  end
   @returns flx_sD, flx_sSr, flx_stau
 end
 
diff --git a/src/projects/GRHD/rhs.jl b/src/projects/GRHD/rhs.jl
index 457208bb50bb96fc9ce7f3df5a867b9fbfb01620..c40fd0b802133f131f25866412b470519d0eb096 100644
--- a/src/projects/GRHD/rhs.jl
+++ b/src/projects/GRHD/rhs.jl
@@ -226,7 +226,7 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
     bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh, equation, rsolver,
-          ldg_rsolver, av_rsolver                      = rhs
+          ldg_rsolver, av_rsolver                = rhs
   @unpack sD, sSr, stau                          = get_dynamic_variables(cache)
   @unpack flx_sD, flx_sSr, flx_stau,
           ldg_sD, ldg_sSr, ldg_stau,
@@ -239,7 +239,12 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
           lhs_numflx_ldg_sSr,  rhs_numflx_ldg_sSr,
           lhs_numflx_ldg_stau, rhs_numflx_ldg_stau = get_bdry_variables(cache)
 
-  # broadcast_volume!(cons2prim, equation, cache)
+  broadcast_volume!(scaledcons2cons, equation, cache)
+  postprocess_scaledcons2cons!(cache, equation, mesh)
+  # broadcast_volume!(conservative_fixing, equation, cache)
+  broadcast_volume!(cons2prim, equation, cache)
+  # broadcast_volume!(atmosphering, equation, cache)
+  broadcast_volume!(cons2scaledcons, equation, cache)
 
   ## solve auxiliary equation: q + ∂x u = 0
   broadcast_volume!(ldg_flux, equation, cache)
@@ -254,8 +259,6 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
   broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
   broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
 
-  broadcast_volume!(scaledcons2cons, equation, cache)
-  postprocess_scaledcons2cons!(cache, equation, mesh)
   broadcast_volume!(flux, equation, cache)
   broadcast_volume!(sources, equation, cache)
   broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)