diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index a3174a8b4e51f44b9b6597dced40553d094cdd95..8fdfe09127755231f51627ed5000d88c83ccf888 100644
--- a/src/GRHD/GRHD.jl
+++ b/src/GRHD/GRHD.jl
@@ -315,6 +315,33 @@ dg1d.@parameters GRHD begin
   doublecartoon_disable_derivatives = false
   @check doublecartoon_disable_derivatives isa Bool
 
+  """
+  If `true` use a time step that is tailored to the smallest grid spacing.
+  """
+  dt_min_grid_spacing = false
+  @check dt_min_grid_spacing isa Bool
+
+  """
+  If `true` enable the convex hull substep limiter.
+
+  Requires `Mesh.dg_kind = modal_bspline[1,2]`.
+  See also `limiter_limit_log_D, limiter_tci_log_D`.
+  """
+  limiter = false
+  @check limiter isa Bool
+
+  """
+  If `true` apply limiter to `log(D)` instead of `D`.
+  """
+  limiter_limit_log_D = false
+  @check limiter_limit_log_D isa Bool
+
+  """
+  If `true` apply TCI that is used to activate the limiter to `log(D)` instead of `D`.
+  """
+  limiter_tci_log_D = false
+  @check limiter_tci_log_D isa Bool
+
 end
 
 
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 66cf7115fb40a6568308a8ac4f8203c53e894e25..10db122ae2c76ca9005894ab6d0d8f056320b8ef 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -31,6 +31,21 @@ function make_callback(env, P)
 end
 
 
+function make_substep_callback(env, P)
+  if P.prms.limiter
+    if env.mesh.element.kind ∉ (:modal_bspline1, :modal_bspline2)
+      error("Substep limiter only works with Mesh.dg_kind = modal_bspline[1,2].")
+    end
+    tci = TCI.ModalDecayAverage(env.mesh, 1.0, 3.0)
+    cbfn_limiter(u, t) = callback_limiter(env, P, env.mesh, tci)
+    cb_limiter = FunctionCallback(cbfn_limiter, CallbackTiming(always_on=true))
+    return CallbackSet(cb_limiter)
+  else
+    return nothing
+  end
+end
+
+
 #######################################################################
 #                              Equation                               #
 #######################################################################
@@ -543,3 +558,47 @@ function callback_analyze0d(env, P::Project{:spherical1d}, mesh::Mesh1d)
     end
   end
 end
+
+
+#######################################################################
+#                               limiter                               #
+#######################################################################
+
+
+callback_limiter(env, P, mesh) = TODO()
+function callback_limiter(env, P::Project{:spherical1d}, mesh::Mesh1d{<:DGElement}, tci)
+  @unpack D, Sr = get_dynamic_variables(mesh)
+  @unpack D_modal, Sr_modal, log_D, c2p_reset_atm = get_static_variables(mesh)
+  @unpack flag = get_cell_variables(mesh)
+
+  broadcast_volume!(impose_atmosphere_spherical1d, P.equation, mesh)
+
+  @unpack limiter_tci_log_D, limiter_limit_log_D = P.prms
+
+  if limiter_tci_log_D || limiter_limit_log_D
+    @. log_D = log(D)
+  end
+
+  K = mesh.tree.dims[1]
+  dg1d.enter(P.prms.workspace) do
+    tmp_flag = dg1d.borrow(P.prms.workspace, length(flag))
+    TCI.compute_indicator!(tmp_flag, limiter_tci_log_D ? log_D : D, P.tci)
+    for (i,f) in enumerate(tmp_flag)
+      flag[i] = f > 0.9
+    end
+    TCI.compute_indicator!(tmp_flag, Sr, P.tci)
+    for (i,f) in enumerate(tmp_flag)
+      flag[i] = max(flag[i], f > 0.9)
+    end
+    tmp_flag .= (flag .+ flag[end:-1:1])./2
+    flag .= tmp_flag
+  end
+
+  dg1d.limit_slope_convex_hull!(limiter_limit_log_D ? log_D : D, flag, P.prms.workspace, mesh)
+  dg1d.limit_slope_convex_hull!(Sr, flag, P.prms.workspace, mesh)
+
+  if limiter_tci_log_D || limiter_limit_log_D
+    @. D = exp(log_D)
+  end
+
+end
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index 07a708d6b81c828d33ba8813f699cee9753cb97e..ad547ea134a2aa81af54a4350a624982dd9e38c2 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -395,3 +395,91 @@ end
   c2p_reset_atm = Float64(D<ρmin)
   @returns c2p_reset_atm
 end
+
+
+#######################################################################
+#                     primitives to conservatives                     #
+#######################################################################
+
+
+@with_signature function prim2cons_spherical1d(eq::Equation)
+  @accepts ρ, p, ϵ, vr, γrr
+
+  @unpack atm_density, atm_threshold, cold_eos = eq
+  ρatm = atm_density
+  patm = cold_eos(Pressure, Density, ρatm)
+  ϵatm = cold_eos(InternalEnergy, Density, ρatm)
+  ρmin = atm_density * atm_threshold
+
+  γuurr = 1/γrr
+  v2 = γuurr*vr*vr
+  h  = 1.0+ϵ+p/ρ
+  W2 = 1.0/(1.0-v2)
+  W  = sqrt(W2)
+  ρhW2 = ρ*h*W2
+  D  = W*ρ
+  Sr = ρhW2*vr
+  τ  = ρhW2 - p - D
+
+  q = Ï„ / D
+  rr = Sr / D
+  r2 = γuurr*rr*rr
+  r  = sqrt(r2)
+  k = r / (1.0+q)
+  if D < ρmin || !(0.0 ≤ k ≤ 1.0)
+    ρ  = ρatm
+    vr = 0.0
+    W  = 1.0
+    p  = patm
+    ϵ  = ϵatm
+
+    v2 = γuurr*vr*vr
+    h  = 1.0+ϵ+p/ρ
+    W2 = 1.0/(1.0-v2)
+    W  = sqrt(W2)
+    ρhW2 = ρ*h*W2
+    D  = W*ρ
+    Sr = ρhW2*vr
+    τ  = ρhW2 - p - D
+  end
+
+  @returns D, Sr, τ, ρhW2, ρ, p, ϵ, vr
+end
+
+
+@with_signature function impose_atmosphere_spherical1d(eq::Equation)
+  @accepts D, Sr, τ, γrr, c2p_reset_atm_v2
+
+  @unpack atm_density, atm_threshold, cold_eos = eq
+  ρatm = atm_density
+  patm = cold_eos(Pressure, Density, ρatm)
+  ϵatm = cold_eos(InternalEnergy, Density, ρatm)
+  ρmin = atm_density * atm_threshold
+  γuurr = 1/γrr
+
+  q = Ï„ / D
+  rr = Sr / D
+  r2 = γuurr*rr*rr
+  r  = sqrt(r2)
+  k = r / (1.0+q)
+  c2p_reset_atm_v2 = 0.0
+  if D < ρmin || !(0.0 ≤ k ≤ 1.0)
+    c2p_reset_atm_v2 = 1.0
+    ρ  = ρatm
+    vr = 0.0
+    W  = 1.0
+    p  = patm
+    ϵ  = ϵatm
+
+    v2 = γuurr*vr*vr
+    h  = 1.0+ϵ+p/ρ
+    W2 = 1.0/(1.0-v2)
+    W  = sqrt(W2)
+    ρhW2 = ρ*h*W2
+    D  = W*ρ
+    Sr = ρhW2*vr
+    τ  = ρhW2 - p - D
+  end
+
+  @returns D, Sr, Ï„, c2p_reset_atm_v2
+end
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 6a1b86368d6ddb77521157338e0fc236ea76ca01..ad1cafa6fc8b387aa6879048f64bdb58ed1fc5de 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -34,8 +34,12 @@ function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
   vmax = get_maxspeed(mesh, P.equation)
   L, = dg1d.widths(mesh)
   K, = mesh.tree.dims
-  dL = L/K
-  dt = dL / (vmax * dtfactor(mesh))
+  dl = if P.prms.dt_min_grid_spacing
+    min_grid_spacing(mesh)
+  else
+    L/K
+  end
+  dt = dl / (vmax * dtfactor(mesh))
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
@@ -1197,6 +1201,15 @@ end
 
 
 function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
+  if mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
+    modal_rhs!(mesh, P, hrsc)
+  else
+    nodal_rhs!(mesh, P, hrsc)
+  end
+end
+
+
+function nodal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
 
   @unpack cache = mesh
   @unpack equation = P
@@ -1285,6 +1298,146 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
 end
 
 
+function modal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
+
+  @unpack cache = mesh
+  @unpack equation = P
+
+  @unpack D, Sr, Ï„                 = get_dynamic_variables(cache)
+  @unpack rhs_D, rhs_Sr, rhs_Ï„     = get_rhs_variables(cache)
+  @unpack flr_D, flr_Sr, flr_Ï„,
+          max_v, vr, p,
+          src_D, src_Sr, src_Ï„,
+          D_modal, Sr_modal, Ï„_modal,
+          flr_D_modal, flr_Sr_modal, flr_Ï„_modal,
+          src_D_modal, src_Sr_modal, src_Ï„_modal = get_static_variables(cache)
+  @unpack nflr_D, nflr_Sr, nflr_Ï„,
+          bdry_D, bdry_Sr, bdry_Ï„,
+          bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
+  @unpack flag = get_cell_variables(mesh)
+  @unpack c2p_reset_atm, c2p_reset_atm_v2, log_D = get_static_variables(mesh)
+
+  # evaluate nodal values
+  D_modal .= D
+  Sr_modal .= Sr
+  Ï„_modal .= Ï„
+  for (var,var_modal) in ( (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal) )
+    for (v_var,v_var_modal) in zip(eachcell(mesh,var),eachcell(mesh,var_modal))
+      mul!(v_var, P.prms.V_bspline, v_var_modal)
+    end
+  end
+
+  if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
+    broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
+  else
+    broadcast_volume!(cons2prim_spherical1d, equation, mesh)
+  end
+  broadcast_volume!(maxspeed_spherical1d, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, Sr,    bdry_Sr)
+  dg1d.interpolate_face_data!(mesh, Ï„,     bdry_Ï„)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, vr,    bdry_vr)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+  impose_symmetry_sources!(P, mesh)
+  broadcast_faces!(llf_spherical1d, equation, mesh)
+  id = P.prms.id
+  bc = P.prms.bc
+  if id == "bondi_accretion" || id == "bondi_accretion_infall"
+    if bc == "from_id"
+      # already accounts for inflow boundary
+      broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
+  elseif id == "tov"
+    if bc == "tov_symmetric_domain"
+      # nothing todo, atmosphering should take care of outer bdrys
+    elseif bc == "tov_reflective_origin"
+      broadcast_bdry!(bdryllf_spherical1d_tov, equation, mesh)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
+  else
+    TODO("unknown boundary conditions for initial data $id: $bc")
+  end
+
+
+  # evaluate modal coefficients
+  for (var,var_modal) in ( (flr_D,flr_D_modal), (flr_Sr,flr_Sr_modal), (flr_Ï„,flr_Ï„_modal),
+                           (src_D,src_D_modal), (src_Sr,src_Sr_modal), (src_Ï„,src_Ï„_modal),
+                           # (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal)
+                          )
+    for (v_var,v_var_modal) in zip(eachcell(mesh,var),eachcell(mesh,var_modal))
+      mul!(v_var_modal, P.prms.invV_bspline, v_var)
+    end
+  end
+
+  if :D ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_D,  flr_D_modal,  src_D_modal,  nflr_D,  mesh)
+  end
+  if :Sr ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Sr, flr_Sr_modal, src_Sr_modal, nflr_Sr, mesh)
+  end
+  if :τ ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„_modal,  src_Ï„_modal,  nflr_Ï„,  mesh)
+  end
+
+  if !P.prms.atm_evolve
+    broadcast_volume!(determine_atmosphere, P.equation, mesh)
+    broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
+  end
+
+  if bc == "tov_symmetric_domain"
+    # Need to symmetrize data around origin to avoid asymmetries to blow up there.
+    # This is needed, because near the surface the solution becomes asymmetric,
+    # and those errors are amplified when propagating inwards.
+    # Not sure if this is caused by the LDG/AV method or by the DG method itself.
+    # korigin = find_cell_index(0.0, mesh)
+    # dg1d.enter(P.prms.workspace) do
+    #   buf = dg1d.borrow(P.prms.workspace, mesh.element.Npts)
+    #   for (k,(v_rhs_D, v_rhs_Sr, v_rhs_Ï„)) in enumerate(zip(eachcell(mesh,rhs_D),
+    #                                                         eachcell(mesh,rhs_Sr),
+    #                                                         eachcell(mesh,rhs_Ï„)))
+    #     if korigin == k
+    #       # @views buf .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+    #       # v_rhs_D .= buf
+    #       # @views buf .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+    #       # v_rhs_Sr .= buf
+    #       # @views buf .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+    #       # v_rhs_Ï„ .= buf
+    #       # @views v_rhs_D  .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+    #       # @views v_rhs_Sr .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+    #       # @views v_rhs_Ï„  .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+    #       break
+    #     end
+    #   end
+    # end
+    K = mesh.tree.dims[1]
+    dg1d.enter(P.prms.workspace) do
+      buf = dg1d.borrow(P.prms.workspace, mesh.element.Npts*K)
+      @views begin
+        @. buf  = (rhs_D + rhs_D[end:-1:1])/2
+        rhs_D .= buf
+        @. buf = (rhs_Sr - rhs_Sr[end:-1:1])/2
+        rhs_Sr .= buf
+        @. buf  = (rhs_Ï„ + rhs_Ï„[end:-1:1])/2
+        rhs_Ï„ .= buf
+      end
+    end
+  end
+
+  D .= D_modal
+  Sr .= Sr_modal
+  Ï„ .= Ï„_modal
+
+  return
+end
+
+
 function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
 
   @unpack cache = mesh
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 3c76c8bb3320a4cf0c4dbaa76d53108b3e526a44..7638364e5c0ca3bbb515dcd9c932261f637adcdf 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -40,15 +40,35 @@ function Project(env::Environment, prms)
   c2p_enforce_causal_atm = Bool(prms["GRHD"]["c2p_enforce_causal_atm"])
   atm_evolve = Bool(prms["GRHD"]["atm_evolve"])
   hrsc = Symbol(prms["GRHD"]["hrsc"])
+  limiter = prms["GRHD"]["limiter"]
+  limiter_limit_log_D = prms["GRHD"]["limiter_limit_log_D"]
+  limiter_tci_log_D = prms["GRHD"]["limiter_tci_log_D"]
   freeze_vars = Symbol.(prms["GRHD"]["freeze_vars"])
   fv_numerical_flux = prms["GRHD"]["fv_numerical_flux"]
   id = prms["GRHD"]["id"]
   bc = prms["GRHD"]["bc"]
+  dt_min_grid_spacing = prms["GRHD"]["dt_min_grid_spacing"]
+  workspace = dg1d.Workspace()
+  if mesh.element isa DGElement
+    if mesh.element.kind === :modal_bspline2
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
+      invV_bspline = pinv(V_bspline)
+    elseif mesh.element.kind === :modal_bspline1
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs1)
+      invV_bspline = pinv(V_bspline)
+    else
+      V_bspline, invV_bspline = nothing, nothing
+    end
+  else
+    V_bspline, invV_bspline = nothing, nothing
+  end
   fixedprms = (; av_regularization=:covariant, id_smooth=true,
                  bernstein, slope_limiter_method, slope_limiter_tvb_M,
                  c2p_dynamic_atm, atm_evolve, atm_equalize_on_interface,
                  c2p_set_atmosphere_on_failure, c2p_enforce_causal_atm, c2p_enforce_atm,
                  av_drag, av_sensor_abslog_D, problem, freeze_vars,
+                 V_bspline, invV_bspline, workspace, dt_min_grid_spacing,
+                 limiter, limiter_limit_log_D, limiter_tci_log_D,
                  id, bc, hrsc, fv_numerical_flux)
 
   # construct Project
@@ -77,7 +97,9 @@ function Project(env::Environment, prms)
 
   initialdata!(env, P)
   projectcb = make_callback(env, P)
-  append!(env.callbacks, CallbackSet(projectcb.callbacks))
+  append!(env.callbacks, projectcb)
+  substep_projectcb = make_substep_callback(env, P)
+  append!(env.callbacks_substep, substep_projectcb)
 
   return P, nothing
 end
@@ -169,11 +191,13 @@ function register_evolution!(mesh, P::Project{:spherical1d})
     rhs_variablenames     = (:rhs_D, :rhs_Sr, :rhs_Ï„), # rhs of conservatives
     static_variablenames  = (:flr_D, :flr_Sr, :flr_Ï„,  # conservatives' fluxes
                              :src_D, :src_Sr, :src_Ï„,  # sources
-                             :ρ, :vr, :p, :ϵ,          # primitives
+                             :ρ, :vr, :p, :ϵ,                 # primitives
+                             :raw_ρ, :raw_vr, :raw_p, :raw_ϵ, # unlimited primitives
                              :max_v, :r,               # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
                              :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
+                             :log_D,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_Ï„,
                              :bdry_ρ, :bdry_ϵ, :bdry_vr, :bdry_p,
@@ -184,6 +208,12 @@ function register_evolution!(mesh, P::Project{:spherical1d})
   )
   @unpack r, x = get_static_variables(mesh)
   r .= x
+  register_variables!(mesh,
+    static_variablenames = (:D_modal, :Sr_modal, :Ï„_modal,
+                            :flr_D_modal, :flr_Sr_modal, :flr_Ï„_modal,
+                            :src_D_modal, :src_Sr_modal, :src_Ï„_modal,
+                            :c2p_failed)
+  )
 end
 function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
   hrsc = P.prmsdb["GRHD"]["hrsc"]
@@ -480,8 +510,10 @@ end
 function register_analysis!(mesh::Mesh1d, P)
   register_variables!(mesh,
       static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
-                              :c2p_freeze_atm, :c2p_init_admissible, :v),
-      bdry_variablenames = (:bdry_c2p_reset_atm,)
+                              :c2p_freeze_atm, :c2p_init_admissible, :v,
+                              :c2p_reset_atm_v2),
+      bdry_variablenames = (:bdry_c2p_reset_atm,),
+      cell_variablenames = (:flag,),
   )
   if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
     register_variables!(mesh, global_variablenames = (:Mtot,))
@@ -576,8 +608,8 @@ register_tci!(mesh, P, tci::TCI.AbstractTCI) = TODO(tci)
 
 
 function register_tci!(mesh, P::Project{:spherical1d}, tci::TCI.AbstractTCI)
-  TODO(tci)
-  register_variables!(mesh, cell_variablenames = (:flag,:sD_flag,:sSr_flag,:stau_flag))
+  return
+  register_variables!(mesh, cell_variablenames = (:flag,))
 end
 
 function register_tci!(mesh, P::Project{:spherical1d}, tci::TCI.EntropyProduction)