diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 43ce7ac1ed6..89773a41e2f 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -142,6 +142,17 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
+friction_coefficient = AnalysisSurfacePointwise(force_boundary_names,
+                                                SurfaceFrictionCoefficient(rho_inf(),
+                                                                           u_inf(equations),
+                                                                           l_inf()))
+
+pressure_coefficient = AnalysisSurfacePointwise(force_boundary_names,
+                                                SurfacePressureCoefficient(p_inf(),
+                                                                           rho_inf(),
+                                                                           u_inf(equations),
+                                                                           l_inf()))
+
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
                                      save_analysis = true,
@@ -149,7 +160,9 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      analysis_integrals = (drag_coefficient,
                                                            lift_coefficient,
                                                            drag_coefficient_shear_force,
-                                                           lift_coefficient_shear_force))
+                                                           lift_coefficient_shear_force),
+                                     analysis_pointwise = (friction_coefficient,
+                                                           pressure_coefficient))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index d9b9245918e..1b8cb5c5658 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -265,8 +265,10 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress
+       AnalysisSurfacePointwise, AnalysisSurfaceIntegral,
+       DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress,
+       SurfacePressureCoefficient, SurfaceFrictionCoefficient
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 06110d08d28..71d9ec82f12 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -9,11 +9,15 @@
 # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time
 """
     AnalysisCallback(semi; interval=0,
-                            save_analysis=false,
-                            output_directory="out",
-                            analysis_filename="analysis.dat",
-                            extra_analysis_errors=Symbol[],
-                            extra_analysis_integrals=())
+                           save_analysis=false,
+                           output_directory="out",
+                           analysis_filename="analysis.dat",
+                           analysis_errors = union(default_analysis_errors(equations),
+                                                   extra_analysis_errors),
+                           extra_analysis_integrals = (),
+                           analysis_integrals = union(default_analysis_integrals(equations),
+                                                      extra_analysis_integrals),
+                           analysis_pointwise = ())
 
 Analyze a numerical solution every `interval` time steps and print the
 results to the screen. If `save_analysis`, the results are also saved in
@@ -35,6 +39,16 @@ solution and integrated over the computational domain. Some examples for this ar
 [`entropy`](@ref), [`energy_kinetic`](@ref), [`energy_internal`](@ref), and [`energy_total`](@ref).
 You can also write your own function with the same signature as the examples listed above and
 pass it via `extra_analysis_integrals`.
+The default `analysis_integrals` is `(entropy_timederivative,)`.
+You can also request `extra_analysis_integrals` such as [`LiftCoefficientPressure`](@ref) or
+[`DragCoefficientPressure`](@ref) by constructing an [`AnalysisSurfaceIntegral`](@ref) with one of 
+the previously mentioned functions.
+
+Similarly, pointwise, i.e., per quadrature/interpolation point, quantities such at
+[`SurfacePressureCoefficient`](@ref) or [`SurfaceFrictionCoefficient`](@ref) can be computed.
+Instances of these need to be passed into [`AnalysisSurfacePointwise`](@ref) which is then in turn
+passed to `analysis_pointwise`.
+
 See the developer comments about `Trixi.analyze`, `Trixi.pretty_form_utf`, and
 `Trixi.pretty_form_ascii` for further information on how to create custom analysis quantities.
 
@@ -43,7 +57,8 @@ evaluating the computational performance, such as the total runtime, the perform
 (time/DOF/rhs!), the time spent in garbage collection (GC), or the current memory usage (alloc'd
 memory).
 """
-mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals,
+mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, AnalysisPointwise,
+                                InitialStateIntegrals,
                                 Cache}
     start_time::Float64
     start_time_last_analysis::Float64
@@ -56,6 +71,7 @@ mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegra
     analyzer::Analyzer
     analysis_errors::Vector{Symbol}
     analysis_integrals::AnalysisIntegrals
+    analysis_pointwise::AnalysisPointwise
     initial_state_integrals::InitialStateIntegrals
     cache::Cache
 end
@@ -80,6 +96,9 @@ function Base.show(io::IO, ::MIME"text/plain",
         for (idx, integral) in enumerate(analysis_callback.analysis_integrals)
             push!(setup, "│ integral " * string(idx) => integral)
         end
+        for (idx, quantity) in enumerate(analysis_callback.analysis_pointwise)
+            push!(setup, "│ pointwise " * string(idx) => quantity)
+        end
         push!(setup,
               "save analysis to file" => analysis_callback.save_analysis ? "yes" : "no")
         if analysis_callback.save_analysis
@@ -109,6 +128,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                           extra_analysis_integrals = (),
                           analysis_integrals = union(default_analysis_integrals(equations),
                                                      extra_analysis_integrals),
+                          analysis_pointwise = (),
                           RealT = real(solver),
                           uEltype = eltype(cache.elements),
                           kwargs...)
@@ -130,6 +150,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                                          analysis_filename,
                                          analyzer,
                                          analysis_errors, Tuple(analysis_integrals),
+                                         Tuple(analysis_pointwise),
                                          SVector(ntuple(_ -> zero(uEltype),
                                                         Val(nvariables(equations)))),
                                          cache_analysis)
@@ -156,7 +177,7 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
 
     analysis_callback = cb.affect!
     analysis_callback.initial_state_integrals = initial_state_integrals
-    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals = analysis_callback
+    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
 
     if save_analysis && mpi_isroot()
         mkpath(output_directory)
@@ -200,6 +221,8 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
             for quantity in analysis_integrals
                 @printf(io, "   %-14s", pretty_form_ascii(quantity))
             end
+            # Pointwise quantities are not saved in `analysis_filename`,
+            # i.e., `analysis.dat` but handle their own output.
 
             println(io)
         end
@@ -337,8 +360,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
         @notimeit timer() integrator.f(du_ode, u_ode, semi, t)
         u = wrap_array(u_ode, mesh, equations, solver, cache)
         du = wrap_array(du_ode, mesh, equations, solver, cache)
-        # Compute l2_error, linf_error
-        analysis_callback(io, du, u, u_ode, t, semi)
+        # Compute l2_error, linf_error among other quantities
+        analysis_callback(io, du, u, u_ode, t, semi, iter)
 
         mpi_println("─"^100)
         mpi_println()
@@ -365,8 +388,8 @@ end
 
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
-function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
-    @unpack analyzer, analysis_errors, analysis_integrals = analysis_callback
+function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi, iter)
+    @unpack analyzer, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
     cache_analysis = analysis_callback.cache
     _, equations, _, _ = mesh_equations_solver_cache(semi)
 
@@ -485,6 +508,8 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
 
     # additional integrals
     analyze_integrals(analysis_integrals, io, du, u, t, semi)
+    # additional pointwise quantities
+    analyze_pointwise(analysis_pointwise, du, u, t, semi, iter)
 
     return nothing
 end
@@ -570,6 +595,26 @@ function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
     nothing
 end
 
+# Iterate over tuples of pointwise analysis quantities in a type-stable way using "lispy tuple programming".
+function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
+                           semi, iter) where {N}
+
+    # Extract the first pointwise analysis quantity and process it; keep the remaining to be processed later
+    quantity = first(analysis_quantities)
+    remaining_quantities = Base.tail(analysis_quantities)
+
+    analyze(quantity, du, u, t, semi, iter)
+
+    # Recursively call this method with the unprocessed pointwise analysis quantities
+    analyze_pointwise(remaining_quantities, du, u, t, semi, iter)
+    return nothing
+end
+
+# terminate the type-stable iteration over tuples
+function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi, iter)
+    nothing
+end
+
 # used for error checks and EOC analysis
 function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
                                                                 Affect! <:
@@ -591,6 +636,13 @@ function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache)
 end
+# `analyze` function that passes also the iteration number `iter`along.
+# Required for callbacks that handle the write-off of results to disk themselves,
+# such as `AnalysisSurfacePointwise`.
+function analyze(quantity, du, u, t, semi::AbstractSemidiscretization, iter)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, iter)
+end
 function analyze(quantity, du, u, t, mesh, equations, solver, cache)
     integrate(quantity, u, mesh, equations, solver, cache, normalize = true)
 end
@@ -654,6 +706,7 @@ end # @muladd
 include("analysis_dg1d.jl")
 include("analysis_dg2d.jl")
 include("analysis_surface_integral_2d.jl")
+include("analysis_surface_pointwise_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
@@ -661,16 +714,19 @@ include("analysis_dg3d_parallel.jl")
 # This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires
 # `semi` to be passed along to retrieve the current boundary indices, which are non-static 
 # in the case of AMR.
-function analyze(quantity::AnalysisSurfaceIntegral, du, u, t,
-                 semi::AbstractSemidiscretization)
+function analyze(quantity::AnalysisSurfaceIntegral{Variable},
+                 du, u, t,
+                 semi::AbstractSemidiscretization) where {Variable}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache, semi)
 end
 
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
-# Note that this needs to be included after `analysis_surface_integral_2d.jl` to
-# have `VariableViscous` available.
+# Note that this needs to be defined after `analysis_surface_integral_2d.jl` is included 
+# to have `VariableViscous` available.
+# We require `semi` to be passed along to retrieve the current boundary indices, which are non-static 
+# in the case of AMR.
 function analyze(quantity::AnalysisSurfaceIntegral{Variable},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
@@ -682,3 +738,31 @@ function analyze(quantity::AnalysisSurfaceIntegral{Variable},
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
             cache_parabolic)
 end
+
+# This version of `analyze` is used for `AnalysisSurfacePointwise` such as `SurfacePressureCoefficient`.
+# We need the iteration number `iter` to be passed in here 
+# as for `AnalysisSurfacePointwise` the writing to disk is handled by the callback itself.
+function analyze(quantity::AnalysisSurfacePointwise{Variable},
+                 du, u, t,
+                 semi::AbstractSemidiscretization,
+                 iter) where {Variable}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, semi, iter)
+end
+# Special analyze for `SemidiscretizationHyperbolicParabolic` such that
+# precomputed gradients are available. Required for `AnalysisSurfacePointwise` equipped 
+# with `VariableViscous` such as `SurfaceFrictionCoefficient`.
+# As for the inviscid version, we need to pass in the iteration number `iter` as 
+# for `AnalysisSurfacePointwise` the writing to disk is handled by the callback itself.
+function analyze(quantity::AnalysisSurfacePointwise{Variable},
+                 du, u, t,
+                 semi::SemidiscretizationHyperbolicParabolic,
+                 iter) where {
+                              Variable <:
+                              VariableViscous}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    equations_parabolic = semi.equations_parabolic
+    cache_parabolic = semi.cache_parabolic
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
+            cache_parabolic, iter)
+end
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 2b5790e9fb0..71aa95df598 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -18,7 +18,7 @@ the boundary/boundaries associated with particular name(s) given in `boundary_sy
 or `boundary_symbols`.
 For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or
 drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries
   where the quantity of interest is computed
@@ -34,7 +34,7 @@ struct AnalysisSurfaceIntegral{Variable, NBoundaries}
     end
 end
 
-struct ForceState{RealT <: Real}
+struct FlowStateDirectional{RealT <: Real}
     psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream
     rhoinf::RealT
     uinf::RealT
@@ -42,11 +42,11 @@ struct ForceState{RealT <: Real}
 end
 
 struct LiftCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 # Abstract base type used for dispatch of `analyze` for quantities
@@ -54,11 +54,11 @@ end
 abstract type VariableViscous end
 
 struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 """
@@ -85,7 +85,7 @@ function LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientPressure(FlowStateDirectional(psi_lift, rhoinf, uinf, linf))
 end
 
 """
@@ -108,7 +108,7 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientPressure(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientPressure(FlowStateDirectional(psi_drag, rhoinf, uinf, linf))
 end
 
 """
@@ -135,7 +135,8 @@ function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientShearStress(FlowStateDirectional(psi_lift, rhoinf, uinf,
+                                                           linf))
 end
 
 """
@@ -158,13 +159,14 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientShearStress(FlowStateDirectional(psi_drag, rhoinf, uinf,
+                                                           linf))
 end
 
 function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -173,7 +175,7 @@ end
 function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -215,29 +217,31 @@ function viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                    gradients_1, gradients_2)
 
     # Viscous stress vector: Stress tensor * normal vector
-    visc_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2]
-    visc_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2]
+    viscous_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2]
+    viscous_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2]
 
-    return (visc_stress_vector_1, visc_stress_vector_2)
+    return (viscous_stress_vector_1, viscous_stress_vector_2)
 end
 
 function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
-                                               gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
-    return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
+    return (viscous_stress_vector_[1] * psi[1] + viscous_stress_vector_[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
 function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
-                                               gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
-    return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
+    return (viscous_stress_vector_[1] * psi[1] + viscous_stress_vector_[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
new file mode 100644
index 00000000000..145bba56374
--- /dev/null
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -0,0 +1,303 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+# This file contains callbacks that are performed on the surface like computation of
+# pointwise surface forces.
+
+"""
+    AnalysisSurfacePointwise{Variable, NBoundaries}(boundary_symbol_or_boundary_symbols,
+                                                    variable, output_directory = "out")
+
+This struct is used to compute pointwise surface values of a quantity of interest `variable`
+alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`
+or `boundary_symbols`.
+For instance, this can be used to compute the surface pressure coefficient [`SurfacePressureCoefficient`](@ref) or
+surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
+symbol `:Airfoil` in 2D.
+
+- `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries
+  where the quantity of interest is computed
+- `variable::Variable`: Quantity of interest, like lift or drag
+- `output_directory = "out"`: Directory where the pointwise value files are stored.
+"""
+struct AnalysisSurfacePointwise{Variable, NBoundaries}
+    variable::Variable # Quantity of interest, like lift or drag
+    boundary_symbols::NTuple{NBoundaries, Symbol} # Name(s) of the boundary/boundaries
+    output_directory::String
+
+    function AnalysisSurfacePointwise(boundary_symbols::NTuple{NBoundaries, Symbol},
+                                      variable,
+                                      output_directory = "out") where {NBoundaries}
+        return new{typeof(variable), NBoundaries}(variable, boundary_symbols,
+                                                  output_directory)
+    end
+end
+
+struct FlowState{RealT <: Real}
+    rhoinf::RealT
+    uinf::RealT
+    linf::RealT
+end
+
+struct SurfacePressureCoefficient{RealT <: Real}
+    pinf::RealT # Free stream pressure
+    flow_state::FlowState{RealT}
+end
+
+struct SurfaceFrictionCoefficient{RealT <: Real} <: VariableViscous
+    flow_state::FlowState{RealT}
+end
+
+"""
+    SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+
+Compute the surface pressure coefficient
+```math
+C_p \\coloneqq \\frac{p - p_{p_\\infty}}
+                     {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the pressure distribution along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurfacePointwise`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `pinf::Real`: Free-stream pressure
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+    return SurfacePressureCoefficient(pinf, FlowState(rhoinf, uinf, linf))
+end
+
+"""
+SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+
+Compute the surface skin friction coefficient
+```math
+C_f \\coloneqq \\frac{\\boldsymbol \\tau_w  \\boldsymbol n^\\perp}
+                     {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the wall shear stress vector ``\\tau_w`` along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurfacePointwise`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+    return SurfaceFrictionCoefficient(FlowState(rhoinf, uinf, linf))
+end
+
+# Compute local pressure coefficient.
+# Works for both purely hyperbolic and hyperbolic-parabolic systems.
+# C_p(x) = (p(x) - p_inf) / (0.5 * rho_inf * u_inf^2 * l_inf)
+function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
+    p = pressure(u, equations)
+    @unpack pinf = pressure_coefficient
+    @unpack rhoinf, uinf, linf = pressure_coefficient.flow_state
+    return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
+end
+
+# Compute local friction coefficient.
+# Works only in conjunction with a hyperbolic-parabolic system.
+# C_f(x) = (tau_w(x) * n_perp(x)) / (0.5 * rho_inf * u_inf^2 * l_inf)
+function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
+                                                        equations_parabolic,
+                                                        gradients_1, gradients_2)
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
+    @unpack rhoinf, uinf, linf = surface_friction.flow_state
+
+    # Normalize as `normal_direction` is not necessarily a unit vector
+    n = normal_direction / norm(normal_direction)
+    # Tangent vector = perpendicular vector to normal vector
+    t = (-n[2], n[1])
+    return (viscous_stress_vector_[1] * t[1] +
+            viscous_stress_vector_[2] * t[2]) /
+           (0.5 * rhoinf * uinf^2 * linf)
+end
+
+# Compute and save to disk a space-dependent `surface_variable`.
+# For the purely hyperbolic, i.e., non-parabolic case, this is for instance 
+# the pressure coefficient `SurfacePressureCoefficient`.
+# The boundary/boundaries along which this quantity is to be integrated is determined by
+# `boundary_symbols`, which is retrieved from `surface_variable`.
+function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
+                 mesh::P4estMesh{2},
+                 equations, dg::DGSEM, cache, semi, iter)
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coordinates = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+
+    index_range = eachnode(dg)
+
+    global_node_counter = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
+            value = variable(u_node, equations)
+
+            coordinates[global_node_counter, 1] = x[1]
+            coordinates[global_node_counter, 2] = x[2]
+            values[global_node_counter] = value
+
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_counter += 1
+        end
+    end
+
+    # Save to disk
+    save_pointwise_file(surface_variable.output_directory, varname(variable),
+                        coordinates, values, t, iter)
+end
+
+# Compute and save to disk a space-dependent `surface_variable`.
+# For the purely hyperbolic-parabolic case, this is may be for instance 
+# the surface skin fricition coefficient `SurfaceFrictionCoefficient`.
+# The boundary/boundaries along which this quantity is to be integrated is determined by
+# `boundary_symbols`, which is retrieved from `surface_variable`.
+function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
+                 du, u, t, mesh::P4estMesh{2},
+                 equations, equations_parabolic,
+                 dg::DGSEM, cache, semi,
+                 cache_parabolic, iter) where {Variable <: VariableViscous}
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coordinates = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+
+    # Additions for parabolic
+    @unpack viscous_container = cache_parabolic
+    @unpack gradients = viscous_container
+
+    gradients_x, gradients_y = gradients
+
+    index_range = eachnode(dg)
+    global_node_counter = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
+
+            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg,
+                                        i_node, j_node, element)
+            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg,
+                                        i_node, j_node, element)
+
+            # Extract normal direction at nodes which points from the 
+            # fluid cells *outwards*, i.e., *into* the structure.
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node, element)
+
+            # Integral over whole boundary surface
+            value = variable(u_node, normal_direction, x, t, equations_parabolic,
+                             gradients_1, gradients_2)
+
+            coordinates[global_node_counter, 1] = x[1]
+            coordinates[global_node_counter, 2] = x[2]
+            values[global_node_counter] = value
+
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_counter += 1
+        end
+    end
+
+    # Save to disk
+    save_pointwise_file(surface_variable.output_directory, varname(variable),
+                        coordinates, values, t, iter)
+end
+
+varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
+varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
+varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
+
+# Helper function that saves a space-dependent quantity `values`
+# at every solution/quadrature point `coordinates` at 
+# time `t` and iteration `iter` to disk.
+# The file is written to the `output_directory` with name `varname` in HDF5 (.h5) format.
+# The latter two are retrieved from the `surface_variable`,
+# an instantiation of `AnalysisSurfacePointwise`.
+function save_pointwise_file(output_directory, varname, coordinates, values, t, iter)
+    n_points = length(values)
+
+    filename = joinpath(output_directory, varname) * @sprintf("_%06d.h5", iter)
+
+    h5open(filename, "w") do file
+        # Add context information as attributes
+        attributes(file)["n_points"] = n_points
+        attributes(file)["variable_name"] = varname
+
+        file["time"] = t
+        file["timestep"] = iter
+        file["point_coordinates"] = coordinates
+        file["point_data"] = values
+    end
+end
+
+function pretty_form_ascii(::AnalysisSurfacePointwise{<:SurfacePressureCoefficient{<:Any}})
+    "CP(x)"
+end
+function pretty_form_utf(::AnalysisSurfacePointwise{<:SurfacePressureCoefficient{<:Any}})
+    "CP(x)"
+end
+
+function pretty_form_ascii(::AnalysisSurfacePointwise{<:SurfaceFrictionCoefficient{<:Any}})
+    "CF(x)"
+end
+function pretty_form_utf(::AnalysisSurfacePointwise{<:SurfaceFrictionCoefficient{<:Any}})
+    "CF(x)"
+end
+end # muladd
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index ceefb65e99b..d1547f23915 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -729,7 +729,7 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 3000
     end
 
     u_ode = copy(sol.u[end])
@@ -755,6 +755,35 @@ end
 
     @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13)
     @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13)
+
+    # Check cp, cf written to outfiles
+
+    cp_vals = read(Trixi.h5open("out/CP_x_000007.h5"))["point_data"]
+    cf_vals = read(Trixi.h5open("out/CF_x_000007.h5"))["point_data"]
+    @test sort(cp_vals[1:10]) ≈ [
+        1.7642563337066723,
+        1.8309926967660561,
+        1.8676797245949768,
+        1.8715139312836404,
+        1.8742476977595561,
+        1.9596647512001575,
+        1.9749658702816537,
+        2.0400855926847923,
+        2.095692020352621,
+        2.1591959225932738,
+    ]
+    @test sort(cf_vals[1:10]) ≈ [
+        -1.2838932978862392,
+        -1.1230998392551437,
+        -1.1167735453217076,
+        -0.9947974567347279,
+        -0.6878311452704875,
+        -0.45370785698055793,
+        -0.37053727238603945,
+        -0.23515678100990403,
+        0.045526882648679205,
+        0.20383549655613004,
+    ]
 end
 end