diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl
index dbf067e3f..c391e3743 100644
--- a/src/Blocks/Blocks.jl
+++ b/src/Blocks/Blocks.jl
@@ -6,8 +6,9 @@ using ModelingToolkit, Symbolics
 import IfElse: ifelse
 import ..@symcheck
 using ModelingToolkit: getdefault
+using ..DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export RealInput, RealOutput, SISO
diff --git a/src/Blocks/math.jl b/src/Blocks/math.jl
index 86045ec52..91280a91b 100644
--- a/src/Blocks/math.jl
+++ b/src/Blocks/math.jl
@@ -38,19 +38,19 @@ Output the product of a gain matrix with the input signal vector.
   - `output`
 """
 @mtkmodel MatrixGain begin
-    @parameters begin
-        K, [description = "Matrix gain"]
+    @structural_parameters begin
+        K
     end
     begin
-        nout = size(getdefault(K), 1)
-        nin = size(getdefault(K), 2)
+        nout = size(K, 1)
+        nin = size(K, 2)
     end
     @components begin
-        input = RealInput(; nin = size(K, 2))
-        output = RealOutput(; nout = size(K, 1))
+        input = RealInput(; nin = nin)
+        output = RealOutput(; nout = nout)
     end
     @equations begin
-        [(@info i, j; output.u[i] ~ sum(getdefault(K)[i, j] * input.u[j])) for j in 1:nin
+        [(output.u[i] ~ sum(K[i, j] * input.u[j])) for j in 1:nin
          for i in 1:nout]...
     end
 end
diff --git a/src/Blocks/nonlinear.jl b/src/Blocks/nonlinear.jl
index bab526a83..416166cc4 100644
--- a/src/Blocks/nonlinear.jl
+++ b/src/Blocks/nonlinear.jl
@@ -97,7 +97,6 @@ Initial value of state `Y` can be set with `int.y`
         rising = 1.0, [description = "Maximum rising slew rate of SlewRateLimiter"]
         falling = -rising, [description = "Derivative time constant of SlewRateLimiter"]
         Td = 0.001, [description = "Derivative time constant"]
-        y_start
     end
     begin
         getdefault(rising) ≥ getdefault(falling) ||
@@ -105,7 +104,7 @@ Initial value of state `Y` can be set with `int.y`
         getdefault(Td) > 0 ||
             throw(ArgumentError("Time constant `Td` must be strictly positive"))
     end
-    @extend u, y = siso = SISO(y_start = y_start)
+    @extend u, y = siso = SISO(; y_start)
     @equations begin
         D(y) ~ max(min((u - y) / Td, rising), falling)
     end
diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl
index 0b7b278a1..4a476d4a7 100644
--- a/src/Blocks/sources.jl
+++ b/src/Blocks/sources.jl
@@ -76,7 +76,7 @@ Generate constant signal.
 """
 @mtkmodel Constant begin
     @components begin
-        output = RealOutput()
+        output = RealOutput(; unit)
     end
     @parameters begin
         k = 0.0, [description = "Constant output value of block"]
@@ -101,7 +101,7 @@ The input variable `t` can be changed by passing a different variable as the key
         f
     end
     @components begin
-        output = RealOutput()
+        output = RealOutput(; unit)
     end
     @equations begin
         output.u ~ first(getdefault(f))(t)
@@ -135,8 +135,9 @@ Generate sine signal.
     phase = 0,
     offset = 0,
     start_time = 0,
-    smooth = false)
-    @named output = RealOutput()
+    smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
     equation = if smooth == false
         offset + ifelse(t < start_time, 0,
@@ -178,8 +179,9 @@ Cosine signal.
     phase = 0,
     offset = 0,
     start_time = 0,
-    smooth = false)
-    @named output = RealOutput()
+    smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
     equation = if smooth == false
         offset + ifelse(t < start_time, zero(t),
@@ -209,8 +211,8 @@ Generate current time signal.
 
   - `output`
 """
-@component function ContinuousClock(; name, offset = 0, start_time = 0)
-    @named output = RealOutput()
+@component function ContinuousClock(; name, offset = 0, start_time = 0, output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters offset=offset start_time=start_time
     eqs = [
         output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time),
@@ -242,8 +244,9 @@ Generate ramp signal.
     duration = 1,
     offset = 0,
     start_time = 0,
-    smooth = false)
-    @named output = RealOutput()
+    smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters offset=offset start_time=start_time height=height duration=duration
     equation = if smooth == false
         offset + ifelse(t < start_time, 0,
@@ -280,8 +283,9 @@ Generate smooth square signal.
   - `output`
 """
 @component function Square(; name, frequency = 1.0, amplitude = 1.0,
-    offset = 0.0, start_time = 0.0, smooth = false)
-    @named output = RealOutput()
+    offset = 0.0, start_time = 0.0, smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters begin
         frequency = frequency
         amplitude = amplitude
@@ -322,8 +326,8 @@ Generate step signal.
   - `output`
 """
 @component function Step(; name, height = 1, offset = 0, start_time = 0, duration = Inf,
-    smooth = 1e-5)
-    @named output = RealOutput()
+    smooth = 1e-5, output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     duration_numeric = duration
     pars = @parameters offset=offset start_time=start_time height=height duration=duration
     equation = if smooth == false # use comparison in case smooth is a float
@@ -372,8 +376,9 @@ Exponentially damped sine signal.
     phase = 0,
     offset = 0,
     start_time = 0,
-    smooth = false)
-    @named output = RealOutput()
+    smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping
 
     equation = if smooth == false
@@ -413,8 +418,9 @@ Generate smooth triangular signal for frequencies less than or equal to 25 Hz
   - `output`
 """
 @component function Triangular(; name, amplitude = 1.0, frequency = 1.0,
-    offset = 0.0, start_time = 0.0, smooth = false)
-    @named output = RealOutput()
+    offset = 0.0, start_time = 0.0, smooth = false,
+    output__unit = nothing)
+    @named output = RealOutput(; unit = output__unit)
     pars = @parameters begin
         amplitude = amplitude
         frequency = frequency
@@ -604,13 +610,13 @@ data input component.
 # Connectors:
   - `output`
 """
-@component function SampledData(; name, buffer)
+@component function SampledData(; name, buffer, unit = nothing)
     pars = @parameters begin
         buffer = buffer
     end
     vars = []
     systems = @named begin
-        output = RealOutput()
+        output = RealOutput(; unit)
     end
     eqs = [
         output.u ~ get_sampled_data(t, buffer),
diff --git a/src/Blocks/utils.jl b/src/Blocks/utils.jl
index 59db3655b..9254db35a 100644
--- a/src/Blocks/utils.jl
+++ b/src/Blocks/utils.jl
@@ -1,14 +1,30 @@
-@connector function RealInput(; name, nin = 1, u_start = nin > 1 ? zeros(nin) : 0.0)
+@connector function RealInput(; name, nin = 1, u_start = nin > 1 ? zeros(nin) : 0.0, unit = nothing)
     if nin == 1
-        @variables u(t)=u_start [
-            input = true,
-            description = "Inner variable in RealInput $name",
-        ]
+        if unit !== nothing
+            @variables u(t)=u_start [
+                input = true,
+                description = "Inner variable in RealInput $name",
+                unit = unit
+            ]
+        else
+            @variables u(t)=u_start [
+                input = true,
+                description = "Inner variable in RealInput $name",
+            ]
+        end
     else
-        @variables u(t)[1:nin]=u_start [
-            input = true,
-            description = "Inner variable in RealInput $name",
-        ]
+        if unit !== nothing
+            @variables u(t)[1:nin]=u_start [
+                input = true,
+                description = "Inner variable in RealInput $name",
+                unit = unit
+            ]
+        else
+            @variables u(t)[1:nin]=u_start [
+                input = true,
+                description = "Inner variable in RealInput $name",
+            ]
+        end
         u = collect(u)
     end
     ODESystem(Equation[], t, [u...], []; name = name)
@@ -26,18 +42,34 @@ Connector with one input signal of type Real.
 - `u`: Value of the connector; if nin=1 this is a scalar
 """ RealInput
 
-@connector function RealOutput(; name, nout = 1, u_start = nout > 1 ? zeros(nout) : 0.0)
+@connector function RealOutput(; name, nout = 1, u_start = nout > 1 ? zeros(nout) : 0.0, unit = nothing)
     if nout == 1
-        @variables u(t)=u_start [
-            output = true,
-            description = "Inner variable in RealOutput $name",
-        ]
+        if unit !== nothing
+            @variables u(t)=u_start [
+                output = true,
+                description = "Inner variable in RealOutput $name",
+                unit = unit
+            ]
+        else
+            @variables u(t)=u_start [
+                output = true,
+                description = "Inner variable in RealOutput $name",
+            ]
+        end
     else
-        @variables u(t)[1:nout]=u_start [
-            output = true,
-            description = "Inner variable in RealOutput $name",
-        ]
-        u = collect(u)
+        if unit !== nothing
+            @variables u(t)[1:nout]=u_start [
+                output = true,
+                description = "Inner variable in RealOutput $name",
+                unit = unit
+            ]
+        else
+            @variables u(t)[1:nout]=u_start [
+                output = true,
+                description = "Inner variable in RealOutput $name",
+                ]
+            end
+            u = collect(u)
     end
     ODESystem(Equation[], t, [u...], []; name = name)
 end
diff --git a/src/Electrical/Analog/ideal_components.jl b/src/Electrical/Analog/ideal_components.jl
index 82a787e90..68d446339 100644
--- a/src/Electrical/Analog/ideal_components.jl
+++ b/src/Electrical/Analog/ideal_components.jl
@@ -38,7 +38,7 @@ See [OnePort](@ref)
 @mtkmodel Resistor begin
     @extend v, i = oneport = OnePort()
     @parameters begin
-        R, [description = "Resistance"]
+        R, [description = "Resistance", unit = u"Ω"]
     end
     @equations begin
         v ~ i * R
@@ -66,7 +66,7 @@ See [OnePort](@ref)
 @mtkmodel Conductor begin
     @extend v, i = oneport = OnePort()
     @parameters begin
-        G, [description = "Conductance"]
+        G, [description = "Conductance", unit = u"S"]
     end
     @equations begin
         i ~ v * G
@@ -93,13 +93,10 @@ See [OnePort](@ref)
   - `C`: [`F`] Capacitance
 """
 @mtkmodel Capacitor begin
+    @extend v, i = oneport = OnePort(; v = 0.0)
     @parameters begin
-        C, [description = "Capacitance"]
+        C, [description = "Capacitance",unit = u"F"]
     end
-    @variables begin
-        v
-    end
-    @extend v, i = oneport = OnePort(; v = v)
     @equations begin
         D(v) ~ i / C
     end
@@ -125,13 +122,10 @@ See [OnePort](@ref)
   - `L`: [`H`] Inductance
 """
 @mtkmodel Inductor begin
+    @extend v, i = oneport = OnePort(; i = 0.0)
     @parameters begin
-        L, [description = "Inductance"]
-    end
-    @variables begin
-        i
+        L, [description = "Inductance", unit = u"H"]
     end
-    @extend v, i = oneport = OnePort(; i = i)
     @equations begin
         D(i) ~ 1 / L * v
     end
@@ -211,9 +205,9 @@ Temperature dependent electrical resistor
         heat_port = HeatPort()
     end
     @parameters begin
-        R_ref = 1.0, [description = "Reference resistance"]
-        T_ref = 300.15, [description = "Reference temperature"]
-        alpha = 0, [description = "Temperature coefficient of resistance"]
+        R_ref = 1.0, [description = "Reference resistance", unit = u"Ω"]
+        T_ref = 300.15, [description = "Reference temperature", unit = u"K"]
+        alpha = 0, [description = "Temperature coefficient of resistance", unit = u"1/K"]
     end
     @variables begin
         R(t) = R_ref
@@ -248,25 +242,19 @@ Electromotoric force (electric/mechanic transformer)
   - `k`: [`N⋅m/A`] Transformation coefficient
 """
 @mtkmodel EMF begin
+    @extend OnePort()
     @components begin
-        p = Pin()
-        n = Pin()
         flange = Flange()
         support = Support()
     end
     @parameters begin
-        k, [description = "Transformation coefficient"]
+        k, [description = "Transformation coefficient", unit = u"N*m/A"]
     end
     @variables begin
-        v(t) = 0.0
-        i(t) = 0.0
-        phi(t) = 0.0
-        w(t) = 0.0
+        phi(t) = 0.0, [description = "Rotation angle", unit = u"rad"]
+        w(t) = 0.0, [description = "Angular velocity", unit = u"rad/s"]
     end
     @equations begin
-        v ~ p.v - n.v
-        0 ~ p.i + n.i
-        i ~ p.i
         phi ~ flange.phi - support.phi
         D(phi) ~ w
         k * w ~ v
diff --git a/src/Electrical/Analog/sensors.jl b/src/Electrical/Analog/sensors.jl
index 5c6e7f16d..94af63487 100644
--- a/src/Electrical/Analog/sensors.jl
+++ b/src/Electrical/Analog/sensors.jl
@@ -19,7 +19,7 @@ an ideal ammeter.
         n = Pin()
     end
     @variables begin
-        i(t)
+        i(t), [description = "Measured current", unit = u"A"]
     end
     @equations begin
         p.v ~ n.v
@@ -46,7 +46,7 @@ Creates a circuit component which measures the potential at a pin.
         p = Pin()
     end
     @variables begin
-        phi(t)
+        phi(t), [description = "Measured potential", unit = u"V"]
     end
     @equations begin
         p.i ~ 0
@@ -74,7 +74,7 @@ Creates a circuit component that measures the voltage across it. Analogous to an
         n = Pin()
     end
     @variables begin
-        v(t)
+        v(t), [description = "Voltage difference from positive to negative pin", unit = u"V"]
     end
     @equations begin
         p.i ~ 0
@@ -112,7 +112,7 @@ consumed by a circuit.
         current_sensor = CurrentSensor()
     end
     @variables begin
-        power(t)
+        power(t), [description = "Power being consumed", unit = u"W"]
     end
     @equations begin
         connect(voltage_sensor.p, pv)
@@ -150,8 +150,8 @@ Combines a [`VoltageSensor`](@ref) and a [`CurrentSensor`](@ref).
         current_sensor = CurrentSensor()
     end
     @variables begin
-        i(t) = 1.0
-        v(t) = 1.0
+        i(t) = 1.0, [description = "Current", unit = u"A"]
+        v(t) = 1.0, [description = "Voltage", unit = u"V"]
     end
     @equations begin
         connect(voltage_sensor.p, pv)
diff --git a/src/Electrical/Analog/sources.jl b/src/Electrical/Analog/sources.jl
index 807fa91f1..4d3f20b3e 100644
--- a/src/Electrical/Analog/sources.jl
+++ b/src/Electrical/Analog/sources.jl
@@ -16,7 +16,7 @@ See [OnePort](@ref)
 @mtkmodel Voltage begin
     @extend v, i = oneport = OnePort()
     @components begin
-        V = RealInput()
+        V = RealInput(unit = u"V")
     end
     @equations begin
         v ~ V.u
@@ -41,7 +41,7 @@ See [OnePort](@ref)
 @mtkmodel Current begin
     @extend v, i = oneport = OnePort()
     @components begin
-        I = RealInput()
+        I = RealInput(unit = u"A")
     end
     @equations begin
         i ~ I.u
diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl
index 687f92b20..720febb03 100644
--- a/src/Electrical/Electrical.jl
+++ b/src/Electrical/Electrical.jl
@@ -8,8 +8,13 @@ using ModelingToolkit, Symbolics, IfElse
 using ..Thermal: HeatPort
 using ..Mechanical.Rotational: Flange, Support
 using ..Blocks: RealInput, RealOutput
+using ..DynamicQuantities
+using ..DynamicQuantities: @u_str
 
-@parameters t
+import ..rad
+import ..S
+
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export Pin, OnePort
diff --git a/src/Electrical/utils.jl b/src/Electrical/utils.jl
index ba7d0f541..32151812a 100644
--- a/src/Electrical/utils.jl
+++ b/src/Electrical/utils.jl
@@ -1,6 +1,6 @@
 @connector Pin begin
-    v(t)                    # Potential at the pin [V]
-    i(t), [connect = Flow]    # Current flowing into the pin [A]
+    v(t), [description = "Voltage", unit = u"V"]                  # Potential at the pin [V]
+    i(t), [description = "Current", unit = u"A", connect = Flow]    # Current flowing into the pin [A]
 end
 @doc """
     Pin(; name)
@@ -33,8 +33,8 @@ Component with two electrical pins `p` and `n` and current `i` flows from `p` to
         n = Pin()
     end
     @variables begin
-        v(t) = 0.0
-        i(t) = 0.0
+        v(t) = 0.0, [description = "Voltage", unit = u"V"]
+        i(t) = 0.0, [description = "Current", unit = u"A", connect = Flow]
     end
     @equations begin
         v ~ p.v - n.v
@@ -70,10 +70,10 @@ Current `i1` flows from `p1` to `n1` and `i2` from `p2` to `n2`.
         n2 = Pin()
     end
     @variables begin
-        v1(t) = 0.0
-        i1(t) = 0.0
-        v2(t) = 0.0
-        i2(t) = 0.0
+        v1(t) = 0.0, [description = "Voltage", unit = u"V"]
+        i1(t) = 0.0, [description = "Current", unit = u"A"]
+        v2(t) = 0.0, [description = "Voltage", unit = u"V"]
+        i2(t) = 0.0, [description = "Current", unit = u"A"]
     end
     @equations begin
         v1 ~ p1.v - n1.v
diff --git a/src/Magnetic/FluxTubes/FluxTubes.jl b/src/Magnetic/FluxTubes/FluxTubes.jl
index c082b5f0a..4d39a3a9f 100644
--- a/src/Magnetic/FluxTubes/FluxTubes.jl
+++ b/src/Magnetic/FluxTubes/FluxTubes.jl
@@ -1,8 +1,10 @@
 module FluxTubes
 using ModelingToolkit
 using ...Electrical: Pin
+import ...Wb
+using ...DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export PositiveMagneticPort, NegativeMagneticPort, TwoPort
diff --git a/src/Magnetic/FluxTubes/basic.jl b/src/Magnetic/FluxTubes/basic.jl
index e19ab198c..e8d1a97f9 100644
--- a/src/Magnetic/FluxTubes/basic.jl
+++ b/src/Magnetic/FluxTubes/basic.jl
@@ -20,7 +20,7 @@ Idle running branch.
 @mtkmodel Idle begin
     @extend (Phi,) = two_port = TwoPort()
     @equations begin
-        Phi ~ 0
+        Phi ~ 0.0
     end
 end
 
@@ -68,7 +68,7 @@ Constant permeance.
 @mtkmodel ConstantPermeance begin
     @extend V_m, Phi = two_port = TwoPort()
     @parameters begin
-        G_m = 1.0, [description = "Magnetic permeance"]
+        G_m = 1.0, [description = "Magnetic permeance", unit = u"H"]
     end
     @equations begin
         Phi ~ G_m * V_m
@@ -87,7 +87,7 @@ Constant reluctance.
 @mtkmodel ConstantReluctance begin
     @extend V_m, Phi = two_port = TwoPort(; Phi = 0.0)
     @parameters begin
-        R_m = 1.0, [description = "Magnetic reluctance"]
+        R_m = 1.0, [description = "Magnetic reluctance", unit = u"H^-1"]
     end
     @equations begin
         V_m ~ Phi * R_m
@@ -114,9 +114,10 @@ Initial magnetic flux flowing into the port_p can be set with `Phi` ([Wb])
         N, [description = "Number of turns"]
     end
     @variables begin
-        v(t)
-        i(t)
-        Phi
+        v(t),
+        [description = "Voltage difference from positive to negative pin", unit = u"V"]
+        i(t), [description = "Current", unit = u"A"]
+        Phi, [description = "Magnetic flux", unit = Wb]
     end
     @extend V_m, Phi = two_port = TwoPort(; Phi = Phi)
     @components begin
@@ -147,13 +148,14 @@ Initial magnetic flux flowing into the port_p can be set with `Phi` ([`Wb`])
 """
 @mtkmodel EddyCurrent begin
     @variables begin
-        Phi
+        Phi, [description = "Magnetic flux", unit = Wb]
     end
     @parameters begin
-        rho = 0.098e-6, [description = "Resistivity of flux tube material"]
-        l = 1, [description = "Average length of eddy current path"]
-        A = 1, [description = "Cross sectional area of eddy current path"]
-        R = rho * l / A # Electrical resistance of eddy current path
+        rho = 0.098e-6, [description = "Resistivity of flux tube material", unit = u"Ω*m"]
+        l = 1, [description = "Average length of eddy current path", unit = u"m"]
+        A = 1, [description = "Cross sectional area of eddy current path", unit = u"m^2"]
+        R = rho * l / A,
+        [description = "Electrical resistance of eddy current path", unit = u"Ω"]
     end
     @extend (V_m, Phi) = two_port = TwoPort(; Phi = Phi)
     @equations begin
diff --git a/src/Magnetic/FluxTubes/sources.jl b/src/Magnetic/FluxTubes/sources.jl
index d2c5f8dfe..dff4d228d 100644
--- a/src/Magnetic/FluxTubes/sources.jl
+++ b/src/Magnetic/FluxTubes/sources.jl
@@ -7,16 +7,17 @@ Parameters:
 
   - `V_m`: [A] Magnetic potential difference
 """
+
 @mtkmodel ConstantMagneticPotentialDifference begin
     @components begin
         port_p = PositiveMagneticPort()
         port_n = NegativeMagneticPort()
     end
     @parameters begin
-        V_m = 0.0, [description = "Magnetic potential difference"]
+        V_m = 0.0, [description = "Magnetic potential difference", unit = u"A"]
     end
     @variables begin
-        Phi(t)
+        Phi(t), [description = "Magnetic flux", unit = u"Wb"]
     end
     @equations begin
         V_m ~ port_p.V_m - port_n.V_m
@@ -40,10 +41,10 @@ Parameters:
         port_n = NegativeMagneticPort()
     end
     @parameters begin
-        Phi = 0.0, [description = "Magnetic flux"]
+        Phi = 0.0, [description = "Magnetic flux", unit = u"Wb"]
     end
     @variables begin
-        V_m(t)
+        V_m(t), [description = "Magnetic potential difference", unit = u"A"]
     end
     @equations begin
         V_m ~ port_p.V_m - port_n.V_m
diff --git a/src/Magnetic/FluxTubes/utils.jl b/src/Magnetic/FluxTubes/utils.jl
index f6154d81f..15ee132d6 100644
--- a/src/Magnetic/FluxTubes/utils.jl
+++ b/src/Magnetic/FluxTubes/utils.jl
@@ -1,6 +1,7 @@
 @connector MagneticPort begin
-    V_m(t), [description = "Magnetic potential at the port"]
-    Phi(t), [connect = Flow, description = "Magnetic flux flowing into the port"]
+    V_m(t), [description = "Magnetic potential at the port", unit = u"A"]
+    Phi(t),
+    [connect = Flow, description = "Magnetic flux flowing into the port", unit = u"Wb"]
 end
 Base.@doc "Port for a Magnetic system." MagneticPort
 
@@ -30,8 +31,12 @@ Partial component with magnetic potential difference between two magnetic ports
         port_n = NegativeMagneticPort()
     end
     @variables begin
-        V_m(t) = 0.0
-        Phi(t) = 0.0
+        V_m(t) = 0.0, [description = "Magnetic potential at the port", unit = u"A"]
+        Phi(t) = 0.0, [
+            connect = Flow,
+            description = "Magnetic flux flowing into the port",
+            unit = u"Wb",
+        ]
     end
     @equations begin
         V_m ~ port_p.V_m - port_n.V_m
diff --git a/src/Mechanical/MultiBody2D/MultiBody2D.jl b/src/Mechanical/MultiBody2D/MultiBody2D.jl
index ed9ccb80d..53dcf7ae0 100644
--- a/src/Mechanical/MultiBody2D/MultiBody2D.jl
+++ b/src/Mechanical/MultiBody2D/MultiBody2D.jl
@@ -2,8 +2,9 @@ module MultiBody2D
 
 using ModelingToolkit, Symbolics, IfElse
 using ..Translational
+using ...DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export Link
diff --git a/src/Mechanical/Rotational/Rotational.jl b/src/Mechanical/Rotational/Rotational.jl
index 9a6fda31f..4e0546263 100644
--- a/src/Mechanical/Rotational/Rotational.jl
+++ b/src/Mechanical/Rotational/Rotational.jl
@@ -6,8 +6,11 @@ module Rotational
 using ModelingToolkit, Symbolics, IfElse
 using ...Blocks: RealInput, RealOutput
 import ...@symcheck
+using ...DynamicQuantities: @u_str
+import ...Wb
+import ...rad
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export Flange, Support
@@ -22,4 +25,5 @@ include("sources.jl")
 export AngleSensor, SpeedSensor, TorqueSensor, RelSpeedSensor
 include("sensors.jl")
 
+
 end
diff --git a/src/Mechanical/Rotational/components.jl b/src/Mechanical/Rotational/components.jl
index 6df5d3b1c..b667e257f 100644
--- a/src/Mechanical/Rotational/components.jl
+++ b/src/Mechanical/Rotational/components.jl
@@ -16,7 +16,7 @@ Flange fixed in housing at a given angle.
         flange = Flange()
     end
     @parameters begin
-        phi0 = 0.0, [description = "Fixed offset angle of flange"]
+        phi0 = 0.0, [description = "Fixed offset angle of flange", unit = u"rad"]
     end
     @equations begin
         flange.phi ~ phi0
@@ -45,7 +45,7 @@ end
 """
 @mtkmodel Inertia begin
     @parameters begin
-        J, [description = "Moment of inertia"]
+        J, [description = "Moment of inertia", unit = u"kg*m^2"]
     end
     @components begin
         flange_a = Flange()
@@ -55,9 +55,9 @@ end
         @symcheck J > 0 || throw(ArgumentError("Expected `J` to be positive"))
     end
     @variables begin
-        phi(t) = 0.0, [description = "Absolute rotation angle"]
-        w(t) = 0.0, [description = "Absolute angular velocity"]
-        a(t) = 0.0, [description = "Absolute angular acceleration"]
+        phi(t) = 0.0, [description = "Absolute rotation angle", unit = u"rad"]
+        w(t) = 0.0, [description = "Absolute angular velocity", unit = u"rad*s^-1"]
+        a(t) = 0.0, [description = "Absolute angular acceleration", unit = u"rad*s^-2"]
     end
     @equations begin
         phi ~ flange_a.phi
@@ -94,8 +94,8 @@ Linear 1D rotational spring
         @symcheck c > 0 || throw(ArgumentError("Expected `c` to be positive"))
     end
     @parameters begin
-        c, [description = "Spring constant"]
-        phi_rel0 = 0.0, [description = "Unstretched spring angle"]
+        c, [description = "Spring constant", unit = u"N*m"]
+        phi_rel0 = 0.0, [description = "Unstretched spring angle", unit = u"rad"]
     end
     @equations begin
         tau ~ c * (phi_rel - phi_rel0)
@@ -129,7 +129,7 @@ Linear 1D rotational damper
         @symcheck d > 0 || throw(ArgumentError("Expected `d` to be positive"))
     end
     @parameters begin
-        d, [description = "Damping constant"]
+        d, [description = "Damping constant", unit = u"N*m*s*rad^-1"]
     end
     @equations begin
         tau ~ d * w_rel
@@ -159,10 +159,10 @@ Linear 1D rotational spring and damper
   - `phi_rel0`: [`rad`] Unstretched spring angle
 """
 @mtkmodel SpringDamper begin
-    @extend phi_rel, w_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
+    @extend phi_rel, w_rel, a_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
     @variables begin
-        tau_c(t), [description = "Spring torque"]
-        tau_d(t), [description = "Damper torque"]
+        tau_c(t), [description = "Spring torque", unit = u"N*m"]
+        tau_d(t), [description = "Damper torque", unit = u"N*m"]
     end
     @parameters begin
         d, [description = "Damping constant"]
@@ -208,8 +208,8 @@ This element characterizes any type of gear box which is fixed in the ground and
         ratio, [description = "Transmission ratio"]
     end
     @variables begin
-        phi_a(t) = 0.0, [description = "Relative angle between shaft a and the support"]
-        phi_b(t) = 0.0, [description = "Relative angle between shaft b and the support"]
+        phi_a(t) = 0.0, [description = "Relative angle between shaft a and the support", unit = u"rad"]
+        phi_b(t) = 0.0, [description = "Relative angle between shaft b and the support", unit = u"rad"]
     end
     @equations begin
         phi_a ~ flange_a.phi - phi_support
@@ -247,14 +247,13 @@ Friction model: "Armstrong, B. and C.C. de Wit, Friction Modeling and Compensati
   - `tau_brk`: [`N⋅m`] Breakaway friction torque
 """
 @mtkmodel RotationalFriction begin
-    @extend w_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
+    @extend phi_rel, w_rel, a_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
     @parameters begin
-        f, [description = "Viscous friction coefficient"]
-        tau_c, [description = "Coulomb friction torque"]
-        w_brk, [description = "Breakaway friction velocity"]
-        tau_brk, [description = "Breakaway friction torque"]
+        f, [description = "Viscous friction coefficient", unit = u"N*m*s/rad"]
+        tau_c, [description = "Coulomb friction torque", unit = u"N*m"]
+        w_brk, [description = "Breakaway friction velocity", unit = u"rad*s^-1"]
+        tau_brk, [description = "Breakaway friction torque", unit = "N*m"]
     end
-
     begin
         str_scale = sqrt(2 * exp(1)) * (tau_brk - tau_c)
         w_st = w_brk * sqrt(2)
diff --git a/src/Mechanical/Rotational/sensors.jl b/src/Mechanical/Rotational/sensors.jl
index d3be89feb..b90c5271c 100644
--- a/src/Mechanical/Rotational/sensors.jl
+++ b/src/Mechanical/Rotational/sensors.jl
@@ -11,7 +11,7 @@ Ideal sensor to measure the absolute flange angle
 @mtkmodel AngleSensor begin
     @components begin
         flange = Flange()
-        phi = RealOutput()
+        phi = RealOutput(unit = u"rad")
     end
     @equations begin
         phi.u ~ flange.phi
@@ -32,7 +32,7 @@ Ideal sensor to measure the absolute flange angular velocity
 @mtkmodel SpeedSensor begin
     @components begin
         flange = Flange()
-        w = RealOutput()
+        w = RealOutput(unit = u"rad/s")
     end
     @equations begin
         D(flange.phi) ~ w.u
@@ -55,7 +55,7 @@ Ideal sensor to measure the torque between two flanges (`= flange_a.tau`)
     @components begin
         flange_a = Flange()
         flange_b = Flange()
-        tau = RealOutput()
+        tau = RealOutput(unit = u"N*m")
     end
     @equations begin
         flange_a.phi ~ flange_b.phi
@@ -72,16 +72,16 @@ Ideal sensor to measure the relative angular velocity
 
   - `flange_a`: [Flange](@ref) Flange of shaft from which sensor information shall be measured
   - `flange_b`: [Flange](@ref) Flange of shaft from which sensor information shall be measured
-  - `w`: [RealOutput](@ref) Absolute angular velocity of flange
+  - `w_rel`: [RealOutput](@ref) Relative angular velocity
 """
 @mtkmodel RelSpeedSensor begin
     @components begin
         flange_a = Flange()
         flange_b = Flange()
-        w_rel = RealOutput()
+        w_rel = RealOutput(unit = u"rad/s")
     end
     @variables begin
-        phi_rel(t) = 0.0
+        phi_rel(t) = 0.0, [description = "Relative rotation angle", unit = u"rad"]
     end
     @equations begin
         0 ~ flange_a.tau + flange_b.tau
diff --git a/src/Mechanical/Rotational/sources.jl b/src/Mechanical/Rotational/sources.jl
index 35abf2c52..8ac70636a 100644
--- a/src/Mechanical/Rotational/sources.jl
+++ b/src/Mechanical/Rotational/sources.jl
@@ -5,9 +5,8 @@
     @extend flange, phi_support = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
     @variables begin
         phi(t),
-        [
-            description = "Angle of flange with respect to support (= flange.phi - support.phi)",
-        ]
+        [description = "Angle of flange with respect to support",
+            unit = u"rad"]
     end
     @equations begin
         phi ~ flange.phi - phi_support
@@ -33,12 +32,9 @@ Input signal acting as external torque on a flange
   - `use_support`
 """
 @mtkmodel Torque begin
-    @parameters begin
-        use_support
-    end
-    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(; use_support)
     @components begin
-        tau = RealInput()
+        tau = RealInput(unit = u"N*m")
     end
     @equations begin
         flange.tau ~ -tau.u
@@ -63,19 +59,20 @@ Constant torque source
 - `tau_constant`: The constant torque applied by the source
 - `use_support`: Whether or not an internal support flange is added.
 """
-@mtkmodel ConstantTorque begin #(; name, tau_constant, use_support = false)
+@mtkmodel ConstantTorque begin
     @parameters begin
         tau_constant,
         [
             description = "Constant torque (if negative, torque is acting as load in positive direction of rotation)",
-        ]
+            unit = u"N*m"]
         use_support
     end
     @extend flange, phi = partial_element = PartialTorque(; use_support = use_support)
     @variables begin
-        tau(t), [description = "Accelerating torque acting at flange (= -flange.tau)"]
-        w(t),
-        [description = "Angular velocity of flange with respect to support (= der(phi))"]
+        tau(t), [description = "Accelerating torque acting at flange (= -flange.tau)",
+            unit = u"N*m"]
+        w(t), [description = "Angular velocity of flange with respect to support",
+            unit = u"rad*s^-1"]
     end
     @equations begin
         w ~ D(phi)
@@ -108,7 +105,9 @@ Forced movement of a flange according to a reference angular velocity signal
     @named partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
     @unpack flange, phi_support = partial_element
     @named w_ref = RealInput()
-    @variables phi(t)=0.0 w(t)=0.0 a(t)=0.0
+    @variables phi(t) = 0.0 [description = "Angle of flange with respect to support", unit = u"rad"]
+    @variables w(t) = 0.0 [description = "Angular velocity", unit = u"rad*s^-1"]
+    @variables a(t) = 0.0 [description = "Angular acceleration", unit = u"rad*s^-2"]
     eqs = [phi ~ flange.phi - phi_support
         D(phi) ~ w]
     if exact
@@ -116,10 +115,13 @@ Forced movement of a flange according to a reference angular velocity signal
         push!(eqs, w ~ w_ref.u)
         push!(eqs, a ~ 0)
     else
-        pars = @parameters tau_filt = tau_filt
+        pars = @parameters tau_filt=tau_filt [
+            description = "Time constant of low-pass filter",
+            unit = u"rad*s^-1",
+        ]
         push!(eqs, D(w) ~ a)
         push!(eqs, a ~ (w_ref.u - w) * tau_filt)
     end
-    return extend(ODESystem(eqs, t, [phi, w, a], pars; name = name, systems = [w_ref]),
+    return extend(ODESystem(eqs, t, vars, pars; name = name, systems = [w_ref]),
         partial_element)
 end
diff --git a/src/Mechanical/Rotational/utils.jl b/src/Mechanical/Rotational/utils.jl
index b54bfd06d..1d6b23fd0 100644
--- a/src/Mechanical/Rotational/utils.jl
+++ b/src/Mechanical/Rotational/utils.jl
@@ -1,6 +1,6 @@
 @connector Flange begin
-    phi(t), [description = "Rotation angle of flange"]
-    tau(t), [connect = Flow, description = "Cut torque in flange"]
+    phi(t), [description = "Rotation angle of flange", unit = u"rad"]
+    tau(t), [connect = Flow, description = "Cut torque in flange", unit = u"N*m"]
 end
 
 Base.@doc """
@@ -14,8 +14,8 @@ Base.@doc """
 """ Flange
 
 @connector Support begin
-    phi(t), [description = "Rotation angle of flange"]
-    tau(t), [connect = Flow, description = "Cut torque in flange"]
+    phi(t), [description = "Rotation angle of flange", unit = u"rad"]
+    tau(t), [connect = Flow, description = "Cut torque in flange", unit = u"N*m"]
 end
 
 # Base.@doc """
@@ -71,8 +71,8 @@ Partial model for the compliant connection of two rotational 1-dim. shaft flange
         flange_b = Flange()
     end
     @variables begin
-        phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges"]
-        tau(t) = 0.0, [description = "Torque between flanges"]
+        phi_rel(t)= 0.0, [description = "Relative rotation angle between flanges", unit = u"rad"]
+        tau(t) = 0.0, [description = "Torque between flanges", unit = u"N*m"]
     end
     @equations begin
         phi_rel ~ flange_b.phi - flange_a.phi
@@ -104,10 +104,13 @@ Partial model for the compliant connection of two rotational 1-dim. shaft flange
         flange_b = Flange()
     end
     @variables begin
-        phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges"]
-        w_rel(t) = 0.0, [description = "Relative angular velocity between flanges"]
-        a_rel(t) = 0.0, [description = "Relative angular acceleration between flanges"]
-        tau(t) = 0.0, [description = "Torque between flanges"]
+        phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges", unit = u"rad"]
+        w_rel(t) = 0.0, [description = "Relative angular velocity between flanges", unit = u"rad*s^-1"]
+        a_rel(t) = 0.0, [
+                description = "Relative angular acceleration between flanges",
+                unit = u"rad*s^-2",
+            ]
+        tau(t) = 0.0, [description = "Torque between flanges", unit = u"N*m"]
     end
     @equations begin
         phi_rel ~ flange_b.phi - flange_a.phi
@@ -138,7 +141,10 @@ Partial model for a component with one rotational 1-dim. shaft flange and a supp
 @component function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
     @named flange = Flange()
     sys = [flange]
-    @variables phi_support(t)=0.0 [description = "Absolute angle of support flange"]
+    @variables phi_support(t)=0.0 [
+        description = "Absolute angle of support flange",
+        unit = u"rad",
+    ]
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support
@@ -173,7 +179,10 @@ Partial model for a component with two rotational 1-dim. shaft flanges and a sup
     @named flange_a = Flange()
     @named flange_b = Flange()
     sys = [flange_a, flange_b]
-    @variables phi_support(t)=0.0 [description = "Absolute angle of support flange"]
+    @variables phi_support(t)=0.0 [
+        description = "Absolute angle of support flange",
+        unit = u"rad",
+    ]
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support
diff --git a/src/Mechanical/Translational/Translational.jl b/src/Mechanical/Translational/Translational.jl
index d20961a43..2560b002e 100644
--- a/src/Mechanical/Translational/Translational.jl
+++ b/src/Mechanical/Translational/Translational.jl
@@ -8,8 +8,9 @@ using ModelingToolkit: getdefault
 
 using ModelingToolkitStandardLibrary.Blocks: RealInput, RealOutput
 using IfElse: ifelse
+using ...DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export MechanicalPort
diff --git a/src/Mechanical/Translational/components.jl b/src/Mechanical/Translational/components.jl
index 959c2aa5e..299b44786 100644
--- a/src/Mechanical/Translational/components.jl
+++ b/src/Mechanical/Translational/components.jl
@@ -12,7 +12,7 @@ Use to close a system that has un-connected `MechanicalPort`'s where the force s
         flange = MechanicalPort()
     end
     @variables begin
-        f(t) = 0.0
+        f(t) = 0.0, [description = "Force", unit = u"N"]
     end
     @equations begin
         flange.f ~ f
@@ -60,21 +60,20 @@ Sliding mass with inertia
 """
 @component function Mass(; name, v = 0.0, m, s = nothing, g = nothing)
     pars = @parameters begin
-        m = m
+        m = m,  [description = "Mass of sliding body", unit = u"kg"]
     end
     @named flange = MechanicalPort(; v = v)
 
-    vars = @variables begin
-        v(t) = v
-        f(t) = 0
-    end
+    @variables v(t) = v [description = "Absolute linear velocity of sliding mass", unit = u"m/s"]
+    @variables f(t) = 0 [description = "Force", unit = u"N"]
+    vars = [v, f]
 
     eqs = [flange.v ~ v
         flange.f ~ f]
 
     # gravity option
     if g !== nothing
-        @parameters g = g
+        @parameters g = g [description = "Gravitational force", unit = u"N"]
         push!(pars, g)
         push!(eqs, D(v) ~ f / m + g)
     else
@@ -83,7 +82,7 @@ Sliding mass with inertia
 
     # position option
     if s !== nothing
-        @variables s(t) = s
+        @variables s(t) = s [description = "Absolute position of sliding mass", unit = u"m"]
         push!(vars, s)
 
         push!(eqs, D(s) ~ v)
@@ -103,7 +102,7 @@ Linear 1D translational spring
 # Parameters:
 
   - `k`: [N/m] Spring constant
-  - `delta_s`: initial spring stretch
+  - `delta_s`: [m] initial spring stretch
   - `va`: [m/s] Initial value of absolute linear velocity at flange_a (default 0 m/s)
   - `v_b_0`: [m/s] Initial value of absolute linear velocity at flange_b (default 0 m/s)
 
@@ -119,12 +118,10 @@ end # default
 @component function Spring(::Val{:relative}; name, k, delta_s = 0.0, flange_a__v = 0.0,
     flange_b__v = 0.0)
     pars = @parameters begin
-        k = k
-    end
-    vars = @variables begin
-        delta_s(t) = delta_s
-        f(t) = 0
+        k = k, [description = "Spring constant", unit = u"N/m"]
     end
+    @variables delta_s(t) = delta_s [description = "Initial spring stretch", unit = u"m"]
+    @variables f(t) = 0 [description = "Force", unit = u"N"]
 
     @named flange_a = MechanicalPort(; v = flange_a__v)
     @named flange_b = MechanicalPort(; v = flange_b__v)
@@ -133,7 +130,7 @@ end # default
         f ~ k * delta_s
         flange_a.f ~ +f
         flange_b.f ~ -f]
-    return compose(ODESystem(eqs, t, vars, pars; name = name),
+    return compose(ODESystem(eqs, t, [delta_s, f], pars; name = name),
         flange_a,
         flange_b) #flange_a.f => +k*delta_s, flange_b.f => -k*delta_s
 end
@@ -142,14 +139,12 @@ const ABS = Val(:absolute)
 @component function Spring(::Val{:absolute}; name, k, sa = 0, sb = 0, flange_a__v = 0.0,
     flange_b__v = 0.0, l = 0)
     pars = @parameters begin
-        k = k
+        k = k, [description = "Spring constant", unit = u"N/m"]
         l = l
     end
-    vars = @variables begin
-        sa(t) = sa
-        sb(t) = sb
-        f(t) = 0
-    end
+    @variables sa(t) = sa [unit = u"m"]
+    @variables sb(t) = sb [unit = u"m"]
+    @variables f(t) = 0 [description = "Force", unit = u"N"]
 
     @named flange_a = MechanicalPort(; v = flange_a__v)
     @named flange_b = MechanicalPort(; v = flange_b__v)
@@ -159,7 +154,7 @@ const ABS = Val(:absolute)
         f ~ k * (sa - sb - l) #delta_s
         flange_a.f ~ +f
         flange_b.f ~ -f]
-    return compose(ODESystem(eqs, t, vars, pars; name = name),
+    return compose(ODESystem(eqs, t, [sa, sb, f], pars; name = name),
         flange_a,
         flange_b) #, flange_a.f => k * (flange_a__s - flange_b__s - l)
 end
@@ -180,11 +175,11 @@ Linear 1D translational damper
 """
 @mtkmodel Damper begin
     @parameters begin
-        d
+        d, [description = "Damping constant", unit = u"N*s/m"]
     end
     @variables begin
-        v(t)
-        f(t) = 0.0
+        v(t), [description = "Velocity", unit = u"m/s"]
+        f(t) = 0.0, [description = "Force", unit = u"N"]
     end
 
     @components begin
diff --git a/src/Mechanical/Translational/sensors.jl b/src/Mechanical/Translational/sensors.jl
index ec778ee58..1dab8adf5 100644
--- a/src/Mechanical/Translational/sensors.jl
+++ b/src/Mechanical/Translational/sensors.jl
@@ -11,7 +11,7 @@ Linear 1D force input sensor.
 @mtkmodel ForceSensor begin
     @components begin
         flange = MechanicalPort()
-        output = RealOutput()
+        output = RealOutput(unit = u"N")
     end
 
     @equations begin
@@ -36,11 +36,11 @@ Linear 1D position input sensor.
 @mtkmodel PositionSensor begin
     @components begin
         flange = MechanicalPort()
-        output = RealOutput()
+        output = RealOutput(unit = u"m")
     end
 
     @variables begin
-        s(t) = 0.0
+        s(t) = 0.0, [description = "Absolute position", unit = u"m"]
     end
 
     @equations begin
diff --git a/src/Mechanical/Translational/sources.jl b/src/Mechanical/Translational/sources.jl
index 71d4f62f4..07a4fd21c 100644
--- a/src/Mechanical/Translational/sources.jl
+++ b/src/Mechanical/Translational/sources.jl
@@ -11,7 +11,7 @@ Linear 1D force input source
 @mtkmodel Force begin
     @components begin
         flange = MechanicalPort(; v = 0.0)
-        f = RealInput()
+        f = RealInput(unit = u"N")
     end
 
     @equations begin
@@ -22,35 +22,33 @@ end
 """
     Position(solves_force = true; name)
 
-Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the position is given, the respective force needed is already provided elsewhere in the model).  
+Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the position is given, the respective force needed is already provided elsewhere in the model).
 
 # Connectors:
 
   - `flange`: 1-dim. translational flange
   - `s`: real input
 """
-@component function Position(solves_force = true; name)
-    vars = []
+@mtkmodel Position begin
+    @structural_parameters begin
+        solves_force = true
+    end
 
-    systems = @named begin
+    @components begin
         flange = MechanicalPort(; v = 0)
-        s = RealInput()
+        s = RealInput(unit = u"m")
     end
 
-    eqs = [
-        D(s.u) ~ flange.v,
-    ]
-
-    !solves_force && push!(eqs, 0 ~ flange.f)
-
-    ODESystem(eqs, t, vars, [];
-        name, systems)
+    @equations begin
+        D(s.u) ~ flange.v
+        0 ~ flange.f
+    end
 end
 
 """
     Velocity(solves_force = true; name)
 
-Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the velocity is given, the respective force needed is already provided elsewhere in the model).  
+Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the velocity is given, the respective force needed is already provided elsewhere in the model).
 
 # Connectors:
 
@@ -60,7 +58,7 @@ Linear 1D position input source.  Set `solves_force=false` to force input force
 @component function Velocity(solves_force = true; name)
     systems = @named begin
         flange = MechanicalPort(; v = 0)
-        v = RealInput()
+        v = RealInput(unit = u"m/s")
     end
 
     eqs = [
@@ -75,25 +73,30 @@ end
 """
 Acceleration(solves_force = true; name)
 
-Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the acceleration is given, the respective force needed is already provided elsewhere in the model).  
+Linear 1D position input source.  Set `solves_force=false` to force input force to 0 (i.e. only the acceleration is given, the respective force needed is already provided elsewhere in the model).
 
 # Connectors:
 
   - `flange`: 1-dim. translational flange
   - `a`: real input
 """
-@component function Acceleration(solves_force = true; name)
-    systems = @named begin
+
+@mtkmodel Acceleration begin#(solves_force = true; name)
+    @structural_parameters begin
+        solves_force = true
+    end
+    @variables begin
+        v(t) = 0, [unit = u"m/s"]
+    end
+    @components begin
         flange = MechanicalPort(; v = 0)
-        a = RealInput()
+        a = RealInput(unit = u"m/s^2")
+    end
+    @equations begin
+        v ~ flange.v
+        D(v) ~ a.u
+        if solves_force
+            0 ~ flange.f
+        end
     end
-
-    vars = @variables v(t) = 0
-
-    eqs = [v ~ flange.v
-        D(v) ~ a.u]
-
-    !solves_force && push!(eqs, 0 ~ flange.f)
-
-    ODESystem(eqs, t, vars, []; name, systems)
 end
diff --git a/src/Mechanical/Translational/utils.jl b/src/Mechanical/Translational/utils.jl
index 7184d6002..085aa3687 100644
--- a/src/Mechanical/Translational/utils.jl
+++ b/src/Mechanical/Translational/utils.jl
@@ -1,6 +1,6 @@
 @connector MechanicalPort begin
-    v(t) = 0.0
-    f(t) = 0.0, [connect = Flow]
+    v(t) = 0.0, [description = "Velocity of the node", unit = u"m/s"]
+    f(t) = 0.0, [connect = Flow, description = "Force entering the node", unit = u"N"]
 end
 Base.@doc """
     MechanicalPort(;name)
diff --git a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
index 2aec491b1..f2cf06b5a 100644
--- a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
+++ b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
@@ -5,8 +5,9 @@ module TranslationalModelica
 
 using ModelingToolkit, Symbolics, IfElse
 using ...Blocks: RealInput, RealOutput
+using ...DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export Flange
diff --git a/src/Mechanical/TranslationalModelica/components.jl b/src/Mechanical/TranslationalModelica/components.jl
index ad8c92b82..f1ea11817 100644
--- a/src/Mechanical/TranslationalModelica/components.jl
+++ b/src/Mechanical/TranslationalModelica/components.jl
@@ -13,7 +13,7 @@ Flange fixed in housing at a given position.
 """
 @mtkmodel Fixed begin
     @parameters begin
-        s0
+        s0 = 0.0, [description = "Absolute position of sliding mass", unit = u"m"]
     end
 
     @components begin
@@ -44,15 +44,14 @@ Sliding mass with inertia
   - `flange: 1-dim. translational flange of mass`
 """
 @mtkmodel Mass begin
+    @extend flange_a, flange_b, s = pr = PartialRigid(; L = 0.0, s = 0.0)
     @parameters begin
-        m = 0.0, [description = "Mass of sliding mass [kg]"]
+        m = 0.0, [description = "Mass of sliding mass", unit = u"kg"]
     end
     @variables begin
-        s
-        v(t) = 0.0, [description = "Absolute linear velocity of sliding mass [m/s]"]
-        a(t) = 0.0, [description = "Absolute linear acceleration of sliding mass [m/s^2]"]
+        v(t) = 0.0, [description = "Absolute linear velocity of sliding mass", unit = u"m/s"]
+        a(t) = 0.0, [description = "Absolute linear acceleration of sliding mass", unit = u"m/s^2"]
     end
-    @extend flange_a, flange_b, s = pr = PartialRigid(; L = 0.0, s = s)
     @equations begin
         v ~ D(s)
         a ~ D(v)
@@ -78,8 +77,8 @@ Linear 1D translational spring
 @mtkmodel Spring begin
     @extend flange_a, flange_b, s_rel, f = pc = PartialCompliant()
     @parameters begin
-        c = 0.0, [description = "Spring constant [N/m]"]
-        s_rel0 = 0.0, [description = "Unstretched spring length [m]"]
+        c = 0.0, [description = "Spring constant", unit = u"N/m"]
+        s_rel0 = 0.0, [description = "Unstretched spring length", unit = u"m"]
     end
 
     @equations begin
@@ -104,10 +103,10 @@ Linear 1D translational damper
 @mtkmodel Damper begin
     @extend flange_a, flange_b, v_rel, f = pc = PartialCompliantWithRelativeStates()
     @parameters begin
-        d = 0.0, [description = "Damping constant [Ns/m]"]
+        d = 0.0, [description = "Damping constant", unit = u"N*s/m"]
     end
     @variables begin
-        lossPower(t) = 0.0, [description = "Power dissipated by the damper [W]"]
+        lossPower(t) = 0.0, [description = "Power dissipated by the damper", unit = u"W"]
     end
     @equations begin
         f ~ d * v_rel
diff --git a/src/Mechanical/TranslationalModelica/sources.jl b/src/Mechanical/TranslationalModelica/sources.jl
index 0110ed2d2..04c5a49dd 100644
--- a/src/Mechanical/TranslationalModelica/sources.jl
+++ b/src/Mechanical/TranslationalModelica/sources.jl
@@ -4,12 +4,9 @@
 Input signal acting as external force on a flange
 """
 @mtkmodel Force begin
-    @parameters begin
-        use_support
-    end
-    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(; use_support)
     @components begin
-        f = RealInput() # Accelerating force acting at flange (= -flange.tau)
+        f = RealInput(unit = u"N") # Accelerating force acting at flange (= -flange.tau)
     end
     @equations begin
         flange.f ~ -f.u
diff --git a/src/Mechanical/TranslationalModelica/utils.jl b/src/Mechanical/TranslationalModelica/utils.jl
index a12016087..c445b8b81 100644
--- a/src/Mechanical/TranslationalModelica/utils.jl
+++ b/src/Mechanical/TranslationalModelica/utils.jl
@@ -1,6 +1,6 @@
 @connector Flange begin
-    s(t)
-    f(t), [connect = Flow]
+    s(t), [description = "Absolute position of sliding mass", unit = u"m"]
+    f(t), [description = "Cut force into the flange", unit = u"N", connect = Flow]
 end
 Base.@doc """
     Flange(;name)
@@ -13,8 +13,8 @@ Base.@doc """
 """ Flange
 
 @connector Support begin
-    s(t)
-    f(t), [connect = Flow]
+    s(t), [description = "Absolute position of sliding mass", unit = u"m"]
+    f(t), [description = "Cut force into the flange", unit = u"N", connect = Flow]
 end
 Base.@doc """
     Support(;name)
@@ -46,8 +46,8 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
 @mtkmodel PartialCompliant begin
     @extend (flange_a, flange_b) = pt = PartialTwoFlanges()
     @variables begin
-        s_rel(t) = 0.0, [description = "Relative distance between flanges"]
-        f(t) = 0.0, [description = "Force between flanges"]
+        s_rel(t) = 0.0, [description = "Relative distance between flanges", unit = u"m"]
+        f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
     end
 
     @equations begin
@@ -71,9 +71,9 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
 @mtkmodel PartialCompliantWithRelativeStates begin
     @extend flange_a, flange_b = pt = PartialTwoFlanges()
     @variables begin
-        s_rel(t) = 0.0, [description = "Relative distance between flanges"]
-        v_rel(t) = 0.0, [description = "Relative linear velocity))"]
-        f(t) = 0.0, [description = "Forces between flanges"]
+        s_rel(t) = 0.0, [description = "Relative distance between flanges", unit = u"m"]
+        v_rel(t) = 0.0, [description = "Relative linear velocity", unit = u"m/s"]
+        f(t) = 0.0, [description = "Forces between flanges", unit = u"N"]
     end
 
     @equations begin
@@ -99,9 +99,11 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
 """
 function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
     @named flange = Flange()
-    @variables s_support(t) [description = "Absolute position of support flange"]
+    @variables s_support(t) [description = "Absolute position of support flange", unit = u"m"]
     @variables s(t) [
         description = "Distance between flange and support (= flange.s - support.s)",
+        unit = u"m"
+
     ]
     eqs = [s ~ flange.s - s_support]
     if use_support
@@ -130,9 +132,9 @@ Partial model for a component with two translational 1-dim. flanges and a suppor
 function PartialElementaryTwoFlangesAndSupport2(; name, use_support = false)
     @named flange = Flange()
 
-    @variables s_a(t) [description = "Distance between left flange and support"]
-    @variables s_b(t) [description = "Distance between right flange and support"]
-    @variables s_support(t) [description = "Absolute position of support flange"]
+    @variables s_a(t) [description = "Distance between left flange and support", unit = u"m"]
+    @variables s_b(t) [description = "Distance between right flange and support", unit = u"m"]
+    @variables s_support(t) [description = "Absolute position of support flange", unit = u"m"]
 
     eqs = [s_a ~ flange_a.s - s_support
         s_b ~ flange_b.s - s_support]
@@ -149,10 +151,10 @@ end
 @mtkmodel PartialRigid begin
     @extend flange_a, flange_b = ptf = PartialTwoFlanges()
     @variables begin
-        s(t) = 0.0, [description = "Absolute position of center of component"]
+        s(t) = 0.0, [description = "Absolute position of center of component", unit = u"m"]
     end
     @parameters begin
-        L = 0.0, [description = "Length of component, from left flange to right flange"]
+        L = 0.0, [description = "Length of component, from left flange to right flange", unit = u"m"]
     end
     @equations begin
         flange_a.s ~ s - L / 2
diff --git a/src/Mechanical/TranslationalPosition/TranslationalPosition.jl b/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
index 72c91835c..eb6a58153 100644
--- a/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
+++ b/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
@@ -5,8 +5,9 @@ module TranslationalPosition
 
 using ModelingToolkit, Symbolics, IfElse
 using ...Blocks: RealInput, RealOutput
+using ...DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export Flange
diff --git a/src/Mechanical/TranslationalPosition/components.jl b/src/Mechanical/TranslationalPosition/components.jl
index e1901cb24..2b5c7e495 100644
--- a/src/Mechanical/TranslationalPosition/components.jl
+++ b/src/Mechanical/TranslationalPosition/components.jl
@@ -13,7 +13,7 @@ Flange fixed in housing at a given position.
 """
 @mtkmodel Fixed begin
     @parameters begin
-        s_0
+        s_0, [description = "Fixed offset position of housing", unit = u"m"]
     end
     @components begin
         flange = Flange(; s = s_0)
@@ -31,13 +31,12 @@ Sliding mass with inertia
 # Parameters:
 
   - `m`: [kg] Mass of sliding mass
-  - `s_0`: [m] Initial value of absolute position of sliding mass
-  - `v_0`: [m/s] Initial value of absolute linear velocity of sliding mass
 
 # States:
 
-  - `s`: [m] Absolute position of sliding mass
-  - `v`: [m/s] Absolute linear velocity of sliding mass (= der(s))
+  - `s`: [m] Absolute position of sliding mass. It accepts an initial value, which defaults to 0.0.
+  - `v`: [m/s] Absolute linear velocity of sliding mass (= der(s)). It accepts an initial value, which defaults to 0.0.
+  - `f`: [N] Force. It accepts an initial value, which defaults to 0.0.
 
 # Connectors:
 
@@ -45,12 +44,12 @@ Sliding mass with inertia
 """
 @mtkmodel Mass begin
     @parameters begin
-        m
+        m, [description = "Mass of sliding mass", unit = u"kg"]
     end
     @variables begin
-        s(t) = 0.0
-        v(t) = 0.0
-        f(t) = 0.0
+        s(t) = 0.0, [description = "Absolute position of sliding mass", unit = u"m"]
+        v(t) = 0.0, [description = "Absolute linear velocity of sliding mass", unit = u"m*s^-1"]
+        f(t) = 0.0, [description = "Force", unit = u"N"]
     end
     @components begin
         flange = Flange(; s = s)
@@ -105,7 +104,7 @@ Linear 1D translational spring
 # Parameters:
 
   - `k`: [N/m] Spring constant
-  - `l`: Unstretched spring length
+  - `l`: [m] Unstretched spring length
   - `flange_a__s`: [m] Initial value of absolute position of flange_a
   - `flange_b__s`: [m] Initial value of absolute position of flange_b
 
@@ -122,12 +121,10 @@ end #default function
     name, k, flange_a__s = 0,
     flange_b__s = 0, l = 0)
     pars = @parameters begin
-        k = k
-        l = l
-    end
-    vars = @variables begin
-        f(t) = k * (flange_a__s - flange_b__s - l)
+        k = k, [description = "Spring constant", unit = u"N/m"]
+        l = l, [description = "Unstretched spring length", unit = u"m"]
     end
+    @variables f(t) = k * (flange_a__s - flange_b__s - l) [description = "Force", unit = u"N"]
 
     @named flange_a = Flange(; s = flange_a__s, f = k * (flange_a__s - flange_b__s - l))
     @named flange_b = Flange(; s = flange_a__s, f = -k * (flange_a__s - flange_b__s - l))
@@ -137,7 +134,7 @@ end #default function
         f ~ k * (flange_a.s - flange_b.s - l) #delta_s
         flange_a.f ~ +f
         flange_b.f ~ -f]
-    return compose(ODESystem(eqs, t, vars, pars; name = name), flange_a, flange_b)
+    return compose(ODESystem(eqs, t, [f], pars; name = name), flange_a, flange_b)
 end
 
 """
@@ -158,12 +155,12 @@ Linear 1D translational damper
 """
 @mtkmodel Damper begin
     @parameters begin
-        d
+        d, [description = "Damping constant", unit = u"N*s/m"]
     end
     @variables begin
-        va(t) = 0.0
-        vb(t) = 0.0
-        f(t) = +(va - vb) * d
+        va(t) = 0.0, [description = "Velocity of flage a", unit = u"m/s"]
+        vb(t) = 0.0, [description = "Velocity of flage b", unit = u"m/s"]
+        f(t) = +(va - vb) * d, [description = "Force", unit = u"N"]
     end
 
     @components begin
diff --git a/src/Mechanical/TranslationalPosition/sources.jl b/src/Mechanical/TranslationalPosition/sources.jl
index 67a94a152..63c9eab89 100644
--- a/src/Mechanical/TranslationalPosition/sources.jl
+++ b/src/Mechanical/TranslationalPosition/sources.jl
@@ -4,12 +4,9 @@
 Input signal acting as external force on a flange
 """
 @mtkmodel Force begin
-    @parameters begin
-        use_support
-    end
-    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+    @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(; use_support)
     @components begin
-        f = RealInput() # Accelerating force acting at flange (= -flange.tau)
+        f = RealInput(unit = u"N") # Accelerating force acting at flange (= -flange.tau)
     end
     @equations begin
         flange.f ~ -f.u
diff --git a/src/Mechanical/TranslationalPosition/utils.jl b/src/Mechanical/TranslationalPosition/utils.jl
index 4e3dd86d9..107d7be82 100644
--- a/src/Mechanical/TranslationalPosition/utils.jl
+++ b/src/Mechanical/TranslationalPosition/utils.jl
@@ -1,6 +1,6 @@
 @connector Flange begin
-    s(t)
-    f(t), [connect = Flow]
+    s(t), [description = "Absolute position of flange", unit = u"m"]
+    f(t), [connect = Flow, description = " Cut force into the flange", unit = u"N"]
 end
 Base.@doc """
     Flange(;name)
@@ -13,8 +13,8 @@ Base.@doc """
 """ Flange
 
 @connector Support begin
-    s(t)
-    f(t), [connect = Flow]
+    s(t), [description = "Absolute position of the support/housing", unit = u"m"]
+    f(t), [connect = Flow, description = "Cut force into the flange", unit = u"N"]
 end
 Base.@doc """
     Support(;name)
@@ -36,16 +36,16 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
   - `s_rel`: [m] Relative distance (= flange_b.s - flange_a.s). It accepts an initial value, which defaults to 0.0.
   - `f`: [N] Force between flanges (= flange_b.f). It accepts an initial value, which defaults to 0.0.
 """
-@mtkmodel PartialCompliant begin#(; name, s_rel_start = 0.0, f_start = 0.0)
+@mtkmodel PartialCompliant begin
     @components begin
         flange_a = Flange()
         flange_b = Flange()
     end
     @variables begin
-        v_a(t) = 0.0
-        v_b(t) = 0.0
-        s_rel(t) = 0.0
-        f(t) = 0.0
+        v_a(t) = 0.0, [description = "Velocity", unit = u"m/s"]
+        v_b(t) = 0.0, [description = "Velocity", unit = u"m/s"]
+        s_rel(t) = 0.0, [description = "Relative distance ", unit = u"m"]
+        f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
     end
     @equations begin
         D(flange_a.s) ~ v_a
@@ -61,18 +61,9 @@ end
 
 Partial model for the compliant connection of two translational 1-dim. flanges.
 
-# Parameters:
-
-  - `s_rel_start`: [m] Initial relative distance
-  - `v_rel_start`: [m/s] Initial relative linear velocity (= der(s_rel))
-  - `a_rel_start`: [m/s²] Initial relative linear acceleration (= der(v_rel))
-  - `f_start`: [N] Initial force between flanges
-
-# States:
+  # States:
 
-  - `s_rel`: [m] Relative distance (= flange_b.phi - flange_a.phi)
-  - `v_rel`: [m/s] Relative linear velocity (= der(s_rel))
-  - `a_rel`: [m/s²] Relative linear acceleration (= der(v_rel))
+  - `delta_s`: [m]
   - `f`: [N] Force between flanges (= flange_b.f)
 """
 @mtkmodel PartialCompliantWithRelativeStates begin
@@ -81,8 +72,8 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
         flange_b = Flange()
     end
     @variables begin
-        delta_s(t) = 0.0
-        f(t) = 0.0
+        delta_s(t) = 0.0, [unit = u"m"]
+        f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
     end
     @equations begin
         delta_s ~ flange_a.s - flange_b.s
@@ -107,7 +98,8 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
 @component function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
     @named flange = Flange()
     sys = [flange]
-    @variables s_support(t)
+    @variables s_support(t),
+    [description = "Absolute position of support flange", unit = u"m"]
     if use_support
         @named support = Support()
         eqs = [support.s ~ s_support
@@ -120,7 +112,7 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
 end
 
 """
-    PartialElementaryTwoFlangesAndSupport2(;name, use_support=false)
+    PartialElementaryTwoFlangesAndSupport2(; name, use_support = false)
 
 Partial model for a component with two translational 1-dim. flanges and a support used for textual modeling, i.e., for elementary models
 
@@ -136,7 +128,7 @@ Partial model for a component with two translational 1-dim. flanges and a suppor
     @named flange_a = Flange()
     @named flange_b = Flange()
     sys = [flange_a, flange_b]
-    @variables s_support(t) = 0.0
+    @variables s_support(t) = 0.0 [description = "Absolute position of support flange", unit = u"m"]
     if use_support
         @named support = Support()
         eqs = [support.s ~ s_support
diff --git a/src/ModelingToolkitStandardLibrary.jl b/src/ModelingToolkitStandardLibrary.jl
index d764f8d81..bff0ab19b 100644
--- a/src/ModelingToolkitStandardLibrary.jl
+++ b/src/ModelingToolkitStandardLibrary.jl
@@ -1,9 +1,10 @@
 module ModelingToolkitStandardLibrary
 import Symbolics: unwrap
+using DynamicQuantities
 
 """
   @symcheck J > 0 || throw(ArgumentError("Expected `J` to be positive"))
-  
+
 Omits the check expression if the argument `J` is symbolic.
 """
 macro symcheck(ex)
@@ -16,6 +17,7 @@ macro symcheck(ex)
     end
 end
 
+include("utils.jl")
 include("Blocks/Blocks.jl")
 include("Mechanical/Mechanical.jl")
 include("Thermal/Thermal.jl")
diff --git a/src/Thermal/HeatTransfer/ideal_components.jl b/src/Thermal/HeatTransfer/ideal_components.jl
index f4e0e19c0..027e1e045 100644
--- a/src/Thermal/HeatTransfer/ideal_components.jl
+++ b/src/Thermal/HeatTransfer/ideal_components.jl
@@ -21,11 +21,11 @@ Lumped thermal element storing heat
         port = HeatPort()
     end
     @parameters begin
-        C, [description = "Heat capacity of element"]
+        C, [description = "Heat capacity of element", unit = u"J/K"]
     end
     @variables begin
-        T(t) = 273.15 + 20
-        der_T(t) = 0.0
+        T(t) = 273.15 + 20, [description = "Temperature of element", unit = u"K"]
+        der_T(t) = 0.0, [description = "Time derivative of temperature", unit = u"K/s"]
     end
 
     @equations begin
@@ -56,7 +56,7 @@ see [`Element1D`](@ref)
 @mtkmodel ThermalConductor begin
     @extend Q_flow, dT = element1d = Element1D()
     @parameters begin
-        G
+        G, [description = "Constant thermal conductance of material", unit = u"W/K"]
     end
     @equations begin
         Q_flow ~ G * dT
@@ -85,7 +85,7 @@ Lumped thermal element transporting heat without storing it.
 @mtkmodel ThermalResistor begin
     @extend Q_flow, dT = element1d = Element1D()
     @parameters begin
-        R
+        R, [description = "Constant thermal resistance of material", unit = u"K/W"]
     end
     @equations begin
         dT ~ R * Q_flow
@@ -114,7 +114,7 @@ Lumped thermal element for heat convection.
 @mtkmodel ConvectiveConductor begin
     @extend Q_flow, dT = convective_element1d = ConvectiveElement1D()
     @parameters begin
-        G
+        G, [description = "Convective thermal conductance", unit = u"W/K"]
     end
     @equations begin
         Q_flow ~ G * dT
@@ -143,7 +143,7 @@ Lumped thermal element for heat convection.
 @mtkmodel ConvectiveResistor begin
     @extend Q_flow, dT = convective_element1d = ConvectiveElement1D()
     @parameters begin
-        R
+        R, [description = "Constant thermal resistance of material", unit = u"K/W"]
     end
     @equations begin
         dT ~ R * Q_flow
@@ -176,7 +176,7 @@ Lumped thermal element for radiation heat transfer.
 
     @extend Q_flow, dT, port_a, port_b = element1d = Element1D()
     @parameters begin
-        G
+        G, [description = "Net radiation conductance between two surfaces", unit = "m^2"]
     end
     @equations begin
         Q_flow ~ G * sigma * (port_a.T^4 - port_b.T^4)
diff --git a/src/Thermal/HeatTransfer/sensors.jl b/src/Thermal/HeatTransfer/sensors.jl
index 2282af693..3df90b753 100644
--- a/src/Thermal/HeatTransfer/sensors.jl
+++ b/src/Thermal/HeatTransfer/sensors.jl
@@ -20,7 +20,7 @@ lags are associated with this sensor model.
         port = HeatPort()
     end
     @variables begin
-        T(t)
+        T(t), [description = "Absolute temperature", unit = u"K"]
     end
     @equations begin
         T ~ port.T
@@ -51,7 +51,7 @@ output signal in kelvin.
         port_b = HeatPort()
     end
     @variables begin
-        T(t)
+        T(t), [description = "Relative temperature", unit = u"K"]
     end
     @equations begin
         T ~ port_a.T - port_b.T
@@ -85,7 +85,7 @@ The output signal is positive, if the heat flows from `port_a` to `port_b`.
         port_b = HeatPort()
     end
     @variables begin
-        Q_flow(t)
+        Q_flow(t), [connect = Flow, description = "Heat flow rate", unit = u"W"]
     end
     @equations begin
         port_a.T ~ port_b.T
diff --git a/src/Thermal/HeatTransfer/sources.jl b/src/Thermal/HeatTransfer/sources.jl
index 1ee8d36f7..74863af4e 100644
--- a/src/Thermal/HeatTransfer/sources.jl
+++ b/src/Thermal/HeatTransfer/sources.jl
@@ -19,9 +19,9 @@ the component FixedHeatFlow is connected, if parameter `Q_flow` is positive.
 """
 @mtkmodel FixedHeatFlow begin
     @parameters begin
-        Q_flow = 1.0, [description = "Fixed heat flow rate at port"]
-        T_ref = 293.15, [description = "Reference temperature"]
-        alpha = 0.0, [description = "Temperature coefficient of heat flow rate"]
+        Q_flow = 1.0, [description = "Fixed heat flow rate at port", unit = u"W"]
+        T_ref = 293.15, [description = "Reference temperature", unit = u"K"]
+        alpha = 0.0, [description = "Temperature coefficient of heat flow rate", unit = u"1/K"]
     end
     @components begin
         port = HeatPort()
@@ -52,7 +52,7 @@ This model defines a fixed temperature `T` at its port in kelvin, i.e., it defin
         port = HeatPort()
     end
     @parameters begin
-        T, [description = "Fixed temperature boundary condition"]
+        T, [description = "Fixed temperature boundary condition", unit = u"K"]
     end
     @equations begin
         port.T ~ T
@@ -82,12 +82,12 @@ dependent losses (which are given a reference temperature T_ref).
 """
 @mtkmodel PrescribedHeatFlow begin
     @parameters begin
-        T_ref = 293.15, [description = "Reference temperature"]
-        alpha = 0.0, [description = "Temperature coefficient of heat flow rate"]
+        T_ref = 293.15, [description = "Reference temperature", unit = u"K"]
+        alpha = 0.0, [description = "Temperature coefficient of heat flow rate", unit = u"1/K"]
     end
     @components begin
         port = HeatPort()
-        Q_flow = RealInput()
+        Q_flow = RealInput(; unit = u"W")
     end
     @equations begin
         port.Q_flow ~ -Q_flow.u * (1 + alpha * (port.T - T_ref))
@@ -111,7 +111,7 @@ the temperature at the specified value.
 @mtkmodel PrescribedTemperature begin
     @components begin
         port = HeatPort()
-        T = RealInput()
+        T = RealInput(; unit = u"K")
     end
     @equations begin
         port.T ~ T.u
diff --git a/src/Thermal/Thermal.jl b/src/Thermal/Thermal.jl
index 1a2958fd9..ce0900ac5 100644
--- a/src/Thermal/Thermal.jl
+++ b/src/Thermal/Thermal.jl
@@ -3,9 +3,10 @@ Library of thermal system components to model heat transfer.
 """
 module Thermal
 using ModelingToolkit, Symbolics, IfElse
-using ...Blocks: RealInput, RealOutput
+using ..Blocks: RealInput, RealOutput
+using ..DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 export HeatPort, Element1D
diff --git a/src/Thermal/utils.jl b/src/Thermal/utils.jl
index fe50d4c7d..ba38a9731 100644
--- a/src/Thermal/utils.jl
+++ b/src/Thermal/utils.jl
@@ -1,6 +1,6 @@
 @connector HeatPort begin
-    T(t) = 273.15 + 20.0
-    Q_flow(t) = 0.0, [connect = Flow]
+    T(t) = 273.15 + 20.0, [description = "Temperature", unit = u"K"]
+    Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
 end
 Base.@doc """
     HeatPort(; name, T = 273.15 + 20.0, Q_flow = 0.0)
@@ -35,8 +35,8 @@ flow rate through the element from `port_a` to `port_b`, `Q_flow`.
         port_b = HeatPort()
     end
     @variables begin
-        dT(t) = 0.0
-        Q_flow(t) = 0.0
+        dT(t) = 0.0, [description = "Temperature difference across the component", unit = u"K"]
+        Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
     end
     @equations begin
         dT ~ port_a.T - port_b.T
@@ -69,8 +69,8 @@ flow rate through the element from `solid` to `fluid`, `Q_flow`.
         fluid = HeatPort()
     end
     @variables begin
-        dT(t) = 0.0
-        Q_flow(t) = 0.0
+        dT(t) = 0.0, [description = "Temperature difference across the component", unit = u"K"]
+        Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
     end
     @equations begin
         dT ~ solid.T - fluid.T
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 000000000..c98db13f3
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,6 @@
+# Define units
+
+@register_unit Wb u"m^2*kg*s^-2*A^-1"
+@register_unit H u"kg*m^2*s^−2*A^−2"
+@register_unit rad u"1"
+@register_unit S u"1/Ω"
diff --git a/test/Blocks/continuous.jl b/test/Blocks/continuous.jl
index bdd4739c8..8b778c518 100644
--- a/test/Blocks/continuous.jl
+++ b/test/Blocks/continuous.jl
@@ -1,8 +1,9 @@
 using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 
 #=
 Testing strategy:
diff --git a/test/Blocks/math.jl b/test/Blocks/math.jl
index 1ad2317fc..e6945696e 100644
--- a/test/Blocks/math.jl
+++ b/test/Blocks/math.jl
@@ -3,8 +3,9 @@ using ModelingToolkit, OrdinaryDiffEq, Test
 using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone
 using ModelingToolkit: inputs, unbound_inputs, bound_inputs
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 
 @testset "Gain" begin
     @named c = Constant(; k = 1)
diff --git a/test/Blocks/nonlinear.jl b/test/Blocks/nonlinear.jl
index e13a45047..ca14763ab 100644
--- a/test/Blocks/nonlinear.jl
+++ b/test/Blocks/nonlinear.jl
@@ -2,8 +2,9 @@ using ModelingToolkit, OrdinaryDiffEq
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 
 @testset "Limiter" begin
     @testset "Constant" begin
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index b2df09461..4bd864ba0 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -4,8 +4,9 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp
     smooth_square, smooth_step, smooth_ramp,
     smooth_triangular, triangular, square
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 @testset "Constant" begin
@@ -408,14 +409,15 @@ end
         amplitude, damping, phase,
         offset, start_time) atol=1e-3
 end
-
+#### TODO
+#=
 @testset "SampledData" begin
     using DataInterpolations
 
     dt = 4e-4
     t_end = 10.0
-    time = 0:dt:t_end
-    x = @. time^2 + 1.0
+    time_span = 0:dt:t_end
+    x = @. time_span^2 + 1.0
 
     vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
     @named src = SampledData(Float64)
@@ -441,3 +443,4 @@ end
     @test sol[dy][end]≈2 * time[end] atol=1e-3
     @test sol[ddy][end]≈2 atol=1e-3
 end
+=#
\ No newline at end of file
diff --git a/test/Blocks/test_analysis_points.jl b/test/Blocks/test_analysis_points.jl
index 39cc222ef..cdbd33e32 100644
--- a/test/Blocks/test_analysis_points.jl
+++ b/test/Blocks/test_analysis_points.jl
@@ -4,6 +4,7 @@ using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 using ModelingToolkit: get_eqs, vars, @set!, get_iv
 using ControlSystemsBase
+using DynamicQuantities: @u_str
 
 @named P = FirstOrder(k = 1, T = 1)
 @named C = Gain(; k = -1)
@@ -231,10 +232,11 @@ Si = ss(matrices...)
 @test tf(So) ≈ tf(Si)
 
 ## A simple multi-level system with loop openings
-@parameters t
+@parameters t [unit = u"s"]
 @named P_inner = FirstOrder(k = 1, T = 1)
 @named feedback = Feedback()
 @named ref = Step()
+### TODO
 @named sys_inner = ODESystem([connect(P_inner.output, :y, feedback.input2)
         connect(feedback.output, :u, P_inner.input)
         connect(ref.output, :r, feedback.input1)], t,
diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl
index 35fce60d5..46e427806 100644
--- a/test/Electrical/analog.jl
+++ b/test/Electrical/analog.jl
@@ -4,13 +4,15 @@ using ModelingToolkitStandardLibrary.Blocks: Step,
     Square, Triangular
 using ModelingToolkitStandardLibrary.Blocks: square, triangular
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-# using Plots
+@parameters t [unit = u"s"]
+D = Differential(t)
 
-@parameters t
+# using Plots
 
 @testset "sensors" begin
-    @named source = Sine(offset = 1, amplitude = 10, frequency = 5)
+    @named source = Sine(offset = 1, amplitude = 10, frequency = 5, output__unit = u"V")
     @named voltage = Voltage()
     @named resistor = Resistor(R = 1)
     @named capacitor = Capacitor(C = 1, v = 0.0)
@@ -43,7 +45,7 @@ using OrdinaryDiffEq: ReturnCode.Success
             power_sensor,
         ])
     sys = structural_simplify(model)
-    prob = ODAEProblem(sys, [], (0.0, 10.0))
+    prob = ODEProblem(sys, [], (0.0, 10.0))
     sol = solve(prob, Tsit5())
 
     # Plots.plot(sol; vars=[capacitor.v, voltage_sensor.v])
@@ -57,7 +59,7 @@ end
 
 # simple voltage divider
 @testset "voltage divider with a short branch" begin
-    @named source = Constant(k = 10)
+    @named source = Constant(k = 10, output.unit = u"V")
     @named voltage = Voltage()
     @named R0 = Resistor(R = 1e3)
     @named R1 = Resistor(R = 1e3)
@@ -86,7 +88,7 @@ end
 
 # simple RC
 @testset "RC" begin
-    @named source = Constant(k = 10)
+    @named source = Constant(k = 10, output.unit = u"V")
     @named voltage = Voltage()
     @named resistor = Resistor(R = 1)
     @named capacitor = Capacitor(C = 1, v = 0.0)
@@ -100,7 +102,7 @@ end
     @named model = ODESystem(connections, t;
         systems = [resistor, capacitor, source, voltage, ground])
     sys = structural_simplify(model)
-    prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
     sol = solve(prob, Tsit5())
 
     # Plots.plot(sol; vars=[source.v, capacitor.v])
@@ -110,7 +112,7 @@ end
 
 # simple RL
 @testset "RL" begin
-    @named source = Constant(k = 10)
+    @named source = Constant(k = 10, output.unit = u"V")
     @named voltage = Voltage()
     @named resistor = Resistor(R = 1)
     @named inductor = Inductor(L = 1.0, i = 0.0)
@@ -124,7 +126,7 @@ end
     @named model = ODESystem(connections, t;
         systems = [resistor, inductor, source, voltage, ground])
     sys = structural_simplify(model)
-    prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
     sol = solve(prob, Tsit5())
 
     # Plots.plot(sol; vars=[inductor.i, inductor.i])
@@ -135,15 +137,15 @@ end
 @testset "RC with voltage sources" begin
     R, C = 1, 1
     @named voltage = Voltage()
-    @named source_const = Constant(k = 10)
+    @named source_const = Constant(k = 10, output.unit = u"V")
     @named source_sin = Sine(offset = 1, amplitude = 10, frequency = 2, start_time = 0.5,
-        phase = 0)
-    @named source_step = Step(offset = 1, height = 10, start_time = 0.5)
+        phase = 0, output__unit = u"V")
+    @named source_step = Step(offset = 1, height = 10, start_time = 0.5, output__unit = u"V")
     @named source_tri = Triangular(offset = 1, start_time = 0.5, amplitude = 10,
-        frequency = 2)
+        frequency = 2, output__unit = u"V")
     @named source_dsin = ExpSine(offset = 1, amplitude = 10, frequency = 2,
-        start_time = 0.5, phase = 0, damping = 0.5)
-    @named source_ramp = Ramp(offset = 1, height = 10, start_time = 0.5, duration = 1)
+        start_time = 0.5, phase = 0, damping = 0.5, output__unit = u"V")
+    @named source_ramp = Ramp(offset = 1, height = 10, start_time = 0.5, duration = 1, output__unit = u"V")
     sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp]
 
     @named resistor = Resistor(; R)
@@ -159,7 +161,7 @@ end
         @named model = ODESystem(connections, t;
             systems = [resistor, capacitor, source, ground, voltage])
         sys = structural_simplify(model)
-        prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
         sol = solve(prob, Tsit5())
         @test sol.retcode == Success
         sol = solve(prob, Rodas4())
@@ -173,7 +175,7 @@ end
 @testset "RC with current sources" begin
     start_time = 2
     @named current = Current()
-    @named source = Step(start_time = 2)
+    @named source = Step(start_time = 2, output__unit = u"A")
     @named resistor = Resistor(R = 1)
     @named capacitor = Capacitor(C = 1, v = 0.0)
     @named ground = Ground()
@@ -186,7 +188,7 @@ end
     @named model = ODESystem(connections, t;
         systems = [ground, resistor, current, capacitor, source])
     sys = structural_simplify(model)
-    prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
     sol = solve(prob, Tsit5())
     y(x, st) = (x .> st) .* abs.(collect(x) .- st)
     @test sol.retcode == Success
@@ -202,7 +204,7 @@ end
     @named R2 = Resistor(R = 100 * R)
     @named C1 = Capacitor(C = 1 / (2 * pi * f * R), v = 0.0)
     @named opamp = IdealOpAmp()
-    @named square_source = Square(amplitude = Vin)
+    @named square_source = Square(amplitude = Vin, output__unit = u"V")
     @named voltage = Voltage()
     @named sensor = VoltageSensor()
 
@@ -251,15 +253,19 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
     @named ground = Ground()
     @named voltage = Voltage()
     @named voltage_sensor = VoltageSensor()
-    @named step = Step(start_time = st, offset = o, height = h)
+    @named step = Step(start_time = st, offset = o, height = h, output__unit = u"V")
     @named cosine = Cosine(offset = o, amplitude = A, frequency = f, start_time = st,
-        phase = ϕ)
-    @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st, phase = ϕ)
+        phase = ϕ, output__unit = u"V")
+    @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st,
+        phase = ϕ, output__unit = u"V")
     @named damped_sine = ExpSine(offset = o, amplitude = A, frequency = f, start_time = st,
-        phase = ϕ, damping = d)
-    @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h)
-    @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f)
-    @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f)
+        phase = ϕ, damping = d, output__unit = u"V")
+    @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h,
+        output__unit = u"V")
+    @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f,
+        output__unit = u"V")
+    @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f,
+        output__unit = u"V")
     # @named vsawtooth = SawTooth(amplitude=A, start_time=st, frequency=f, offset=o)
 
     sources = [step, cosine, sine, damped_sine, ramp, tri, vsquare] #, vsawtooth]
@@ -294,7 +300,7 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
 
         u0 = [cap.v => 0.0]
 
-        prob = ODAEProblem(vsys, u0, (0, 10.0))
+        prob = ODEProblem(vsys, u0, (0, 10.0))
         sol = solve(prob, dt = 0.1, Tsit5())
 
         @test sol.retcode == Success
@@ -314,15 +320,19 @@ end
     @named cap = Capacitor(C = 1, v = 0.0)
     @named current_sensor = CurrentSensor()
     @named current = Current()
-    @named step = Step(start_time = st, offset = o, height = h)
+    @named step = Step(start_time = st, offset = o, height = h, output__unit = u"A")
     @named cosine = Cosine(offset = o, amplitude = A, frequency = f, start_time = st,
-        phase = ϕ)
-    @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st, phase = ϕ)
+        phase = ϕ, output__unit = u"A")
+    @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st,
+        phase = ϕ, output__unit = u"A")
     @named damped_sine = ExpSine(offset = o, amplitude = A, frequency = f, start_time = st,
-        phase = ϕ, damping = d)
-    @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h)
-    @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f)
-    @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f)
+        phase = ϕ, damping = d, output__unit = u"A")
+    @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h,
+        output__unit = u"A")
+    @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f,
+        output__unit = u"A")
+    @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f,
+        output__unit = u"A")
     # @named isawtooth = SawTooth(amplitude=A, start_time=st, frequency=f, offset=o)
 
     sources = [step, cosine, sine, damped_sine, ramp, tri, vsquare] #, idamped_sine]
@@ -358,7 +368,7 @@ end
 
         u0 = [cap.v => 0.0]
 
-        prob = ODAEProblem(isys, u0, (0, 10.0))
+        prob = ODEProblem(isys, u0, (0, 10.0))
         sol = solve(prob, dt = 0.1, Tsit5())
 
         @test sol.retcode == Success
diff --git a/test/Mechanical/multibody.jl b/test/Mechanical/multibody.jl
index c00203c84..46d2a0f53 100644
--- a/test/Mechanical/multibody.jl
+++ b/test/Mechanical/multibody.jl
@@ -5,8 +5,6 @@ using DifferentialEquations
 # using Setfield
 using Test
 
-@parameters t
-
 @named link1 = Link(; m = 1, l = 10, I = 84, g = -9.807)
 @named link2 = Link(; m = 1, l = 10, I = 84, g = -9.807, x1_0 = 10)
 @named cart = Mass(; m = 1, s_0 = 0)
diff --git a/test/Mechanical/rotational.jl b/test/Mechanical/rotational.jl
index 79bbbf3ae..6198d68d8 100644
--- a/test/Mechanical/rotational.jl
+++ b/test/Mechanical/rotational.jl
@@ -6,7 +6,7 @@ using OrdinaryDiffEq: ReturnCode.Success
 
 # using Plots
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 @testset "two inertias" begin
diff --git a/test/Mechanical/translational.jl b/test/Mechanical/translational.jl
index 2af419ea3..34f09e6b1 100644
--- a/test/Mechanical/translational.jl
+++ b/test/Mechanical/translational.jl
@@ -1,17 +1,19 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
 using ModelingToolkitStandardLibrary.Blocks
+import ModelingToolkitStandardLibrary: Mechanical
 import ModelingToolkitStandardLibrary.Mechanical.Translational as TV
 import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition as TP
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 @testset "Free" begin
     function System(; name)
         systems = @named begin
             acc = TV.Acceleration()
-            a = Constant(; k = -10)
+            a = Constant(; k = -10, output.unit = u"m/s^2")
             mass = TV.Mass(; m = 100)
             free = TV.Free()
         end
@@ -85,7 +87,7 @@ end
     @named fv = TV.Force()
     @named fp = TP.Force(use_support = false)
 
-    @named source = Sine(frequency = 3, amplitude = 2)
+    @named source = Sine(frequency = 3, amplitude = 2, unit = u"N")
 
     function System(damping, spring, body, ground, f, source)
         eqs = [connect(f.f, source.output)
@@ -126,11 +128,11 @@ end
 
             spring = TV.Spring(; k = 1000)
 
-            src1 = Sine(frequency = 100, amplitude = 2)
-            src2 = Sine(frequency = 100, amplitude = -1)
+            src1 = Sine(frequency = 100, amplitude = 2, output__unit = u"m")
+            src2 = Sine(frequency = 100, amplitude = -1, output__unit = u"N")
 
-            pos_value = RealInput()
-            force_output = RealOutput()
+            pos_value = RealInput(unit = u"m")
+            force_output = RealOutput(unit = u"N")
         end
 
         eqs = [connect(pos.s, src1.output)
diff --git a/test/Mechanical/translational_modelica.jl b/test/Mechanical/translational_modelica.jl
index 015b36572..8cd466890 100644
--- a/test/Mechanical/translational_modelica.jl
+++ b/test/Mechanical/translational_modelica.jl
@@ -1,16 +1,17 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
 using ModelingToolkitStandardLibrary.Blocks
-import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as TP
+import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as TM
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 @testset "spring damper mass fixed" begin
-    @named damper = TP.Damper(; d = 1)
-    @named spring = TP.Spring(; c = 1, s_rel0 = 1)
-    @named mass = TP.Mass(; m = 1, v = 1)
-    @named fixed = TP.Fixed(s0 = 1)
+    @named damper = TM.Damper(; d = 1)
+    @named spring = TM.Spring(; c = 1, s_rel0 = 1)
+    @named mass = TM.Mass(; m = 1, v = 1)
+    @named fixed = TM.Fixed(s0 = 1)
 
     eqs = [connect(spring.flange_a, mass.flange_a, damper.flange_a)
         connect(spring.flange_b, damper.flange_b, fixed.flange)]
@@ -29,12 +30,12 @@ D = Differential(t)
 end
 
 @testset "driven spring damper mass" begin
-    @named damper = TP.Damper(; d = 1)
-    @named spring = TP.Spring(; c = 1, s_rel0 = 1)
-    @named mass = TP.Mass(; m = 1, v = 1)
-    @named fixed = TP.Fixed(; s0 = 1)
-    @named force = TP.Force(use_support = false)
-    @named source = Sine(frequency = 3, amplitude = 2)
+    @named damper = TM.Damper(; d = 1)
+    @named spring = TM.Spring(; c = 1, s_rel0 = 1)
+    @named mass = TM.Mass(; m = 1, v = 1)
+    @named fixed = TM.Fixed(; s0 = 1)
+    @named force = TM.Force(use_support = false)
+    @named source = Sine(frequency = 3, amplitude = 2, output__unit = u"N")
 
     eqs = [connect(force.f, source.output)
         connect(force.flange, mass.flange_a)
diff --git a/test/Thermal/demo.jl b/test/Thermal/demo.jl
index 0cc950c56..e27d7925b 100644
--- a/test/Thermal/demo.jl
+++ b/test/Thermal/demo.jl
@@ -1,7 +1,8 @@
 using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 # Modelica example
diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl
index d8f69fd41..1c705272d 100644
--- a/test/Thermal/thermal.jl
+++ b/test/Thermal/thermal.jl
@@ -1,8 +1,9 @@
 using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
 using ModelingToolkitStandardLibrary.Blocks: Constant, Step
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 #=
 # Test HeatCapacitor, TemperatureSensor, RelativeTemperatureSensor, FixedTemperature
@@ -210,8 +211,8 @@ end
     @named T_core = TemperatureSensor()
     @named convection = ConvectiveConductor(G = 25)
     @named environment = PrescribedTemperature()
-    @named amb = Constant(k = T_amb)
-    @named core_losses_const = Constant(k = 500)
+    @named amb = Constant(k = T_amb, output.unit = u"K")
+    @named core_losses_const = Constant(k = 500, output.unit = u"")
     @named winding_losses = Step(height = 900, offset = 100, start_time = 360,
         duration = Inf, smooth = false)
     connections = [connect(windingLosses.port, winding.port)
diff --git a/test/chua_circuit.jl b/test/chua_circuit.jl
index c12896d3c..15aa05c42 100644
--- a/test/chua_circuit.jl
+++ b/test/chua_circuit.jl
@@ -4,22 +4,25 @@ using ModelingToolkitStandardLibrary.Electrical: OnePort
 using OrdinaryDiffEq
 using OrdinaryDiffEq: ReturnCode.Success
 using IfElse: ifelse
+using DynamicQuantities: @u_str
 
 @testset "Chua Circuit" begin
-    @parameters t
+    @parameters t [unit = u"s"]
 
-    @component function NonlinearResistor(; name, Ga, Gb, Ve)
-        @named oneport = OnePort()
-        @unpack v, i = oneport
-        pars = @parameters Ga=Ga Gb=Gb Ve=Ve
-        eqs = [
+    @mtkmodel NonlinearResistor begin
+        @extend OnePort()
+        @parameters begin
+            Ga, [description = "Conductance in inner voltage range", unit = u"1/Ω"]
+            Gb, [description = "Conductance in outer voltage range", unit = u"1/Ω"]
+            Ve, [description = "Inner voltage range limit", unit = u"V"]
+        end
+        @equations begin
             i ~ ifelse(v < -Ve,
                 Gb * (v + Ve) - Ga * Ve,
                 ifelse(v > Ve,
                     Gb * (v - Ve) + Ga * Ve,
-                    Ga * v)),
-        ]
-        extend(ODESystem(eqs, t, [], pars; name = name), oneport)
+                    Ga * v))
+        end
     end
 
     @named L = Inductor(L = 18, i = 0.0)
diff --git a/test/multi_domain.jl b/test/multi_domain.jl
index 6b113c698..6027b45c4 100644
--- a/test/multi_domain.jl
+++ b/test/multi_domain.jl
@@ -6,9 +6,10 @@ using ModelingToolkitStandardLibrary.Thermal
 import ModelingToolkitStandardLibrary
 using ModelingToolkit, OrdinaryDiffEq, Test
 using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
 # using Plots
 
-@parameters t
+@parameters t [unit = u"s"]
 D = Differential(t)
 
 @testset "DC motor" begin
@@ -21,13 +22,13 @@ D = Differential(t)
     tau_L_step = -3
     @named ground = Ground()
     @named source = Voltage()
-    @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+    @named voltage_step = Blocks.Step(height = V_step, start_time = 0, output__unit = u"V")
     @named R1 = Resistor(R = R)
     @named L1 = Inductor(L = L, i = 0.0)
     @named emf = EMF(k = k)
     @named fixed = Fixed()
     @named load = Torque(use_support = false)
-    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3, output__unit = u"N*m")
     @named inertia = Inertia(J = J)
     @named friction = Damper(d = f)
 
@@ -72,7 +73,7 @@ D = Differential(t)
     @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
     @test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
 
-    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+    prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 6.0))
     sol = solve(prob, DFBDF())
     @test sol.retcode == Success
     # EMF equations
@@ -103,13 +104,13 @@ end
     tau_L_step = -3
     @named ground = Ground()
     @named source = Voltage()
-    @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+    @named voltage_step = Blocks.Step(height = V_step, start_time = 0, output__unit = u"V")
     @named R1 = Resistor(R = R)
     @named L1 = Inductor(L = L, i = 0.0)
     @named emf = EMF(k = k)
     @named fixed = Fixed()
     @named load = Torque(use_support = false)
-    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3, output__unit = u"N*m")
     @named inertia = Inertia(J = J)
     @named friction = Damper(d = f)
     @named speed_sensor = SpeedSensor()
@@ -161,7 +162,7 @@ end
 
     @test all(sol[inertia.w] .== sol[speed_sensor.w.u])
 
-    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+    prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 6.0))
     sol = solve(prob, DFBDF())
 
     @test sol.retcode == Success
@@ -183,7 +184,7 @@ end
 @testset "El. Heating Circuit" begin
     @named ground = Ground()
     @named source = Voltage()
-    @named voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1)
+    @named voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1, output__unit = u"V")
     @named heating_resistor = HeatingResistor(R_ref = 100, alpha = 1e-3, T_ref = 293.15)
     @named thermal_conductor = ThermalConductor(G = 50)
     @named env = FixedTemperature(T = 273.15 + 20)
@@ -207,5 +208,6 @@ end
     prob = ODEProblem(sys, Pair[], (0, 6.0))
     sol = solve(prob, Rodas4())
     @test sol.retcode == Success
+    ### TODO
     @test sol[source.v * source.i] == -sol[env.port.Q_flow]
 end