diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 3e4679f5af..e000c8f556 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -15,7 +15,7 @@ using Trixi
 #   Structured and Unstructured Grid Methods (2016)
 #   [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623)
 # - Deep Ray, Praveen Chandrashekar (2017)
-#   An entropy stable finite volume scheme for the 
+#   An entropy stable finite volume scheme for the
 #   two dimensional Navier–Stokes equations on triangular grids
 #   [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020)
 
@@ -142,6 +142,17 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   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 b8364eef44..9a477bcc3a 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -263,8 +263,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 860e3fa21d..dfb7efde89 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -43,7 +43,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 +57,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 +82,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 +114,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 +136,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 +163,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)
@@ -201,6 +208,10 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
                 @printf(io, "   %-14s", pretty_form_ascii(quantity))
             end
 
+            for quantity in analysis_pointwise
+                @printf(io, "   %-14s", pretty_form_ascii(quantity))
+            end
+
             println(io)
         end
     end
@@ -337,8 +348,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 +376,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)
 
@@ -484,6 +495,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
@@ -569,6 +582,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! <:
@@ -590,6 +623,10 @@ 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
+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
@@ -653,6 +690,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")
@@ -660,8 +698,8 @@ 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
@@ -681,3 +719,23 @@ function analyze(quantity::AnalysisSurfaceIntegral{Variable},
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
             cache_parabolic)
 end
+
+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
+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 d92ac5a207..1a95db66d8 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -9,16 +9,16 @@
 # surface forces
 
 """
-    AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi,
-                                                          boundary_symbol_or_boundary_symbols,
-                                                          variable)
+    AnalysisSurfaceIntegral{Variable, NBoundaries}(semi,
+                                                   boundary_symbols::NTuple{NBoundaries, Symbol},
+                                                   variable)
 
 This struct is used to compute the surface integral 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 lift [`LiftCoefficientPressure`](@ref) or
 drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries
@@ -36,7 +36,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
@@ -44,11 +44,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
@@ -56,11 +56,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
 
 """
@@ -87,7 +87,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
 
 """
@@ -110,7 +110,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
 
 """
@@ -137,7 +137,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
 
 """
@@ -160,13 +161,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)
@@ -175,7 +177,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)
@@ -217,29 +219,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 0000000000..7e7de2ef0e
--- /dev/null
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -0,0 +1,280 @@
+# 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
+
+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
+
+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)
+    n_perp = (-n[2], n[1])
+    return (viscous_stress_vector_[1] * n_perp[1] +
+            viscous_stress_vector_[2] * n_perp[2]) /
+           (0.5 * rhoinf * uinf^2 * linf)
+end
+
+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_index = 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_index, 1] = x[1]
+            coordinates[global_node_index, 2] = x[2]
+            values[global_node_index] = value
+
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+
+    save_pointwise_file(surface_variable.output_directory, varname(variable),
+                        coordinates,
+                        values, t, iter)
+end
+
+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_index = 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)
+            # Extract normal direction at nodes which points from the elements outwards,
+            # i.e., *into* the structure.
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    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)
+
+            # Integral over whole boundary surface
+            value = variable(u_node, normal_direction, x, t, equations_parabolic,
+                             gradients_1, gradients_2)
+
+            coordinates[global_node_index, 1] = x[1]
+            coordinates[global_node_index, 2] = x[2]
+            values[global_node_index] = value
+
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+
+    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"
+
+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 d038354f88..ed84322b7c 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,31 @@ 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.5152278762705806
+           1.5182860925755697
+           1.7417814818119175
+           1.871513931283643
+           1.874247697759559
+           2.0354491762572806
+           2.0400855926847936
+           2.0956920203526193
+           2.098639791879171
+           2.1591959225932777]
+    @test sort(cf_vals[1:10]) ≈ [-0.4537078569805583
+           -0.37053727238603884
+           -0.23515678100990495
+           0.045526882648679115
+           0.2038354965561312
+           0.2634810993266428
+           0.42919006623618894
+           0.6611394672752093
+           0.7458959144321556
+           0.8113349266282898]
 end
 end