diff --git a/src/ScalarEq/ScalarEq.jl b/src/ScalarEq/ScalarEq.jl
index 528b764df7f149fd0a0f6d4f1838be66f8c12fc6..5557e6eb558508520486567d287584dcb6bfd193 100644
--- a/src/ScalarEq/ScalarEq.jl
+++ b/src/ScalarEq/ScalarEq.jl
@@ -3,6 +3,7 @@ module ScalarEq
 
 
 using dg1d
+using LinearAlgebra
 using LoopVectorization
 
 
diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 75100c42b3a8369bd1470fe804dc7319af3a2d23..0e07391b1c0d6ae0c0cbcf749180cbdf308a16d6 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -267,7 +267,18 @@ end
 #######################################################################
 
 
-function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+function rhs!(env, mesh::Mesh1d{DGElement}, P::Project, hrsc::Nothing)
+  if mesh.element.kind === :nodal_lagrange
+    nodal_lagrange_rhs!(env, mesh, P, hrsc)
+  elseif mesh.element.kind === :modal_bernstein
+    modal_bernstein_rhs!(env, mesh, P, hrsc)
+  else
+    TODO(mesh.element.kind)
+  end
+end
+
+
+function nodal_lagrange_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
 
   @unpack cache, mesh = env
   @unpack equation = P
@@ -288,6 +299,34 @@ function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
 end
 
 
+function modal_bernstein_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+
+  @unpack cache, mesh = env
+  @unpack equation = P
+
+  @unpack u_modal, flx_u, src_u = get_static_variables(cache)
+  @unpack u                     = get_dynamic_variables(cache)
+  @unpack rhs_u                 = get_rhs_variables(cache)
+  @unpack nflx_u, bdry_u        = get_bdry_variables(cache)
+
+  u_modal .= u
+  for (v_u, v_u_modal) in zip(eachcell(mesh,u),eachcell(mesh,u_modal))
+    mul!(v_u, P.prms.V_bernstein, v_u_modal)
+  end
+
+  dg1d.interpolate_face_data!(mesh, u, bdry_u)
+  broadcast_volume!(flux_source, equation, mesh)
+  broadcast_faces!(llf_1d, equation, mesh)
+  # broadcast_bdry!(bdryllf_1d, equation, P.bdrycond, mesh)
+
+  compute_rhs_weak_form!(rhs_u, flx_u, src_u, nflx_u, mesh)
+
+  u .= u_modal
+
+  return
+end
+
+
 function rhs!(env, mesh::Mesh2d, P::Project, hrsc::Nothing)
 
   @unpack cache = env
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index f92c621ec692f6e37359fd2bc050993e38fb87b4..f8bc0adf34fd46b19f0921e6345198704acb2ffc 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -13,9 +13,11 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   bernstein = BernsteinReconstruction(mesh)
   analyze_error = prms["ScalarEq"]["analyze_error"]
   analyze_error_norm = prms["ScalarEq"]["analyze_error_norm"]
+  V_bernstein = dg1d.Bernstein.vandermonde_matrix(mesh.element.z)
   fixedprms = (; bernstein, av_derivative_scheme, av_drag, hrsc, muscl_omega,
                  slope_limiter_method, slope_limiter_tvb_M,
                  av_recompute_substeps,
+                 V_bernstein,
                  analyze_error, analyze_error_norm)
 
   # construct Project
@@ -123,6 +125,10 @@ function register_equation!(mesh::Mesh1d, P::Project)
     else
       TODO(P.prms.hrsc)
     end
+  elseif mesh.element isa DGElement
+    if mesh.element.kind === :modal_bernstein
+      register_variables!(mesh.cache, static_variablenames=(:u_modal,))
+    end
   end
 end
 function register_equation!(mesh::Mesh2d, P::Project)
diff --git a/src/dgelement.jl b/src/dgelement.jl
index f9de0ac97e7a7b03d29d73398ca7b5c9f9cfcf9a..dcb1914e575def18f029dd0ee7a264112533697f 100644
--- a/src/dgelement.jl
+++ b/src/dgelement.jl
@@ -6,7 +6,7 @@ struct DGElement
   # quadrature rule
   # quadr_w are such that data sampled on `z` can be directly passed to the integrate
   # functions below. In particular, any interpolation from `z` onto `quadr_z`
-  # (need kind === :nodal_bernstein) is already accounted for.
+  # (need kind === :modal_bernstein) is already accounted for.
   quadr::Symbol             # :lgl, :lgl_mass_lumped, :glgl
   quadr_z::Vector{Float64}  # quadrature points
   quadr_w::Vector{Float64}  # quadrature weights
@@ -61,7 +61,7 @@ struct DGElement
       l_rhs = zeros(Float64, N+1)
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
-    elseif kind === :nodal_bernstein
+    elseif kind === :modal_bernstein
       if quadr === :lgl || quadr === :lgl_lumped
         quadr_z, quadr_w = LGL.rule(Npts)
       elseif quadr === :glgl
@@ -82,7 +82,7 @@ struct DGElement
       l_lhs[1]   = 1.0
       l_rhs[end] = 1.0
     else
-      error("Invalid quadrature rule specified")
+      error("Invalid DG formulation specified: $kind. Use on of nodal_lagrange, modal_bernstein.")
     end
 
     MDM    = invM * transpose(S)
diff --git a/src/parameters.jl b/src/parameters.jl
index c13ec135e015face3dc637511955a7165168dd59..31b2078fdd837de5a16429f64db4afbe7de53dca 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -510,12 +510,14 @@ end
   - `nodal_lagrange` ... Nodal approximation using Lagrange polynomials collocated with
     the quadrature nodes provided by `basis`. This is the standard method discussed
     in Hesthaven, Warburton, 2007.
-  - `nodal_bernstein` ... (Quasi-)Nodal approximation using Bernstein polynomials.
-    The coefficients correspond to the nodal values sampled on an equidistant grid.
+  - `modal_bernstein` ... A modal approximation using Bernstein polynomials.
+    The coefficients correspond to the nodal values sampled on an equidistant grid,
+    but the Bernstein polynomials are not interpolating.
+    Hence, this Ansatz is also refered to as quasi-nodal.
     See Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
   """
   dg_kind = "nodal_lagrange"
-  @check scheme in [ "nodal_lagrange", "nodal_bernstein" ]
+  @check scheme in [ "nodal_lagrange", "modal_bernstein" ]
 
   """
   Periodicity in cartesian directions `x,y`
diff --git a/test/IntegrationTests/refs/modal_bernstein_advection_sine/interpolate.h5 b/test/IntegrationTests/refs/modal_bernstein_advection_sine/interpolate.h5
new file mode 100644
index 0000000000000000000000000000000000000000..bf79b76dcc1db77644baf904aa72eeb8aa845e5f
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bernstein_advection_sine/interpolate.h5 differ
diff --git a/test/IntegrationTests/refs/modal_bernstein_advection_sine/modal_bernstein_advection_sine.toml b/test/IntegrationTests/refs/modal_bernstein_advection_sine/modal_bernstein_advection_sine.toml
new file mode 100644
index 0000000000000000000000000000000000000000..6b41c05f43261490bdd4a4db8702a2671847c0d2
--- /dev/null
+++ b/test/IntegrationTests/refs/modal_bernstein_advection_sine/modal_bernstein_advection_sine.toml
@@ -0,0 +1,24 @@
+[Evolution]
+cfl = 0.5
+tend = 2.2
+
+[ScalarEq]
+analyze_error = false
+equation = "advection"
+id = "sine"
+
+[Output]
+substep_every_iteration = 200
+interpolate_variables = ["u"]
+substep_variables = ["u"]
+interpolate_aligned_ts = [0.3, 0.6, 0.8]
+aligned_ts = [0.05, 0.1, 0.15, 0.2]
+interpolate_nodes = "$(collect(range(0.0,1.0,length=10)))"
+variables1d = ["u"]
+
+[Mesh]
+dg_kind = "modal_bernstein"
+range = [0.0, 1.0]
+k = 40
+basis = "lgl"
+n = 4
diff --git a/test/IntegrationTests/refs/modal_bernstein_advection_sine/output1d.h5 b/test/IntegrationTests/refs/modal_bernstein_advection_sine/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..6515693b408270ee7095a7d995b1b9d80881ac4e
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bernstein_advection_sine/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/modal_bernstein_advection_sine/substep_output.h5 b/test/IntegrationTests/refs/modal_bernstein_advection_sine/substep_output.h5
new file mode 100644
index 0000000000000000000000000000000000000000..f67e822f1758cd4bd521039a5362ae5668e170a4
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bernstein_advection_sine/substep_output.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index 5b254ec5ecda84e0d509d04a29d1ce52af4c0721..b5cdbccf82d11f4257b10818dcd97077dd2bfe0b 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -7,6 +7,9 @@ default_abstol1d = 1e-8
 [advection_sine]
 variables1d = [ "u", "u_err" ]
 
+[modal_bernstein_advection_sine]
+variables1d = [ "u" ]
+
 [advection_sine_analysis]
 variables0d = [ "u_err_norm", "u_norm2" ]
 variables1d = [ "u" ]
diff --git a/test/UnitTests/src/test_dgelement.jl b/test/UnitTests/src/test_dgelement.jl
index cf43f55b56bf7f5969a20a588c5402392d77fe5d..028a9d33df9eb60e94f169addea43be4b80b0660 100644
--- a/test/UnitTests/src/test_dgelement.jl
+++ b/test/UnitTests/src/test_dgelement.jl
@@ -209,7 +209,7 @@ end
     @test all(MD.≈S)
   end
 
-  el = dg1d.DGElement(N, :lgl, :nodal_bernstein)
+  el = dg1d.DGElement(N, :lgl, :modal_bernstein)
   T, invT = dg1d.Bernstein.bernstein_legendre_transformation(el.Npts)
   M_legendre = dg1d.Legendre.mass_matrix(el.Npts)
   TMT = T*M_legendre*transpose(T) |> dg1d.purge_zeros