diff --git a/CHANGELOG.md b/CHANGELOG.md
index a59eebab5..1535d376c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -33,6 +33,7 @@ Attention: The newest changes should be on top -->
 ### Added
 
 - ENH: Support for the RSE file format has been added to the library [#798](https://github.com/RocketPy-Team/RocketPy/pull/798)
+- ENH: Introduce Net Thrust with pressure corrections [#789](https://github.com/RocketPy-Team/RocketPy/pull/789)
 
 ### Changed
 
diff --git a/rocketpy/motors/empty_motor.py b/rocketpy/motors/empty_motor.py
index 4e844472b..ff7128c11 100644
--- a/rocketpy/motors/empty_motor.py
+++ b/rocketpy/motors/empty_motor.py
@@ -19,6 +19,7 @@ def __init__(self):
             reshape_thrust_curve=False,
             interpolation_method="linear",
             coordinate_system_orientation="nozzle_to_combustion_chamber",
+            reference_pressure=0,
         )
 
         # Mass properties
diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py
index ea01e686f..1de9f435c 100644
--- a/rocketpy/motors/hybrid_motor.py
+++ b/rocketpy/motors/hybrid_motor.py
@@ -29,6 +29,8 @@ class HybridMotor(Motor):
         "combustion_chamber_to_nozzle".
     HybridMotor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    HybridMotor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     HybridMotor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See
@@ -147,7 +149,11 @@ class HybridMotor(Motor):
         HybridMotor.propellant_I_22 and HybridMotor.propellant_I_33 for
         more information.
     HybridMotor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from thrust source, in Newtons, as a
+        function of time.
+    HybridMotor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     HybridMotor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     HybridMotor.max_thrust : float
@@ -167,7 +173,10 @@ class HybridMotor(Motor):
         Total motor burn duration, in seconds. It is the difference between the
         ``burn_out_time`` and the ``burn_start_time``.
     HybridMotor.exhaust_velocity : Function
-        Propulsion gases exhaust velocity, assumed constant, in m/s.
+        Effective exhaust velocity of the propulsion gases in m/s. Computed
+        as the thrust divided by the mass flow rate. This corresponds to the
+        actual exhaust velocity only when the nozzle exit pressure equals the
+        atmospheric pressure.
     HybridMotor.burn_area : Function
         Total burn area considering all grains, made out of inner
         cylindrical burn area and grain top and bottom faces. Expressed
@@ -181,6 +190,9 @@ class HybridMotor(Motor):
         Method of interpolation used in case thrust curve is given
         by data set in .csv or .eng, or as an array. Options are 'spline'
         'akima' and 'linear'. Default is "linear".
+    HybridMotor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     def __init__(  # pylint: disable=too-many-arguments
@@ -203,6 +215,7 @@ def __init__(  # pylint: disable=too-many-arguments
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Initialize Motor class, process thrust curve and geometrical
         parameters and store results.
@@ -298,6 +311,8 @@ class Function. Thrust units are Newtons.
             positions specified. Options are "nozzle_to_combustion_chamber" and
             "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
 
         Returns
         -------
@@ -314,6 +329,7 @@ class Function. Thrust units are Newtons.
             reshape_thrust_curve=reshape_thrust_curve,
             interpolation_method=interpolation_method,
             coordinate_system_orientation=coordinate_system_orientation,
+            reference_pressure=reference_pressure,
         )
         self.liquid = LiquidMotor(
             thrust_source,
@@ -326,6 +342,7 @@ class Function. Thrust units are Newtons.
             reshape_thrust_curve,
             interpolation_method,
             coordinate_system_orientation,
+            reference_pressure,
         )
         self.solid = SolidMotor(
             thrust_source,
@@ -346,6 +363,7 @@ class Function. Thrust units are Newtons.
             reshape_thrust_curve,
             interpolation_method,
             coordinate_system_orientation,
+            reference_pressure,
         )
 
         self.positioned_tanks = self.liquid.positioned_tanks
@@ -371,6 +389,11 @@ def exhaust_velocity(self):
         -------
         self.exhaust_velocity : Function
             Gas exhaust velocity of the motor.
+
+        Notes
+        -----
+        This corresponds to the actual exhaust velocity only when the nozzle
+        exit pressure equals the atmospheric pressure.
         """
         return Function(
             self.total_impulse / self.propellant_initial_mass
@@ -676,6 +699,7 @@ def from_dict(cls, data):
             grains_center_of_mass_position=data["grains_center_of_mass_position"],
             nozzle_position=data["nozzle_position"],
             throat_radius=data["throat_radius"],
+            reference_pressure=data["reference_pressure"],
         )
 
         for tank in data["positioned_tanks"]:
diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py
index 7e763dd70..5edaae3b8 100644
--- a/rocketpy/motors/liquid_motor.py
+++ b/rocketpy/motors/liquid_motor.py
@@ -29,6 +29,8 @@ class LiquidMotor(Motor):
         "combustion_chamber_to_nozzle".
     LiquidMotor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    LiquidMotor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     LiquidMotor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See
@@ -122,7 +124,11 @@ class LiquidMotor(Motor):
         LiquidMotor.propellant_I_22 and LiquidMotor.propellant_I_33 for
         more information.
     LiquidMotor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from thrust source, in Newtons, as a
+        function of time.
+    LiquidMotor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     LiquidMotor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     LiquidMotor.max_thrust : float
@@ -142,7 +148,13 @@ class LiquidMotor(Motor):
         Total motor burn duration, in seconds. It is the difference between the
         burn_out_time and the burn_start_time.
     LiquidMotor.exhaust_velocity : Function
-        Propulsion gases exhaust velocity in m/s.
+        Effective exhaust velocity of the propulsion gases in m/s. Computed
+        as the thrust divided by the mass flow rate. This corresponds to the
+        actual exhaust velocity only when the nozzle exit pressure equals the
+        atmospheric pressure.
+    LiquidMotor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     def __init__(
@@ -157,6 +169,7 @@ def __init__(
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Initialize LiquidMotor class, process thrust curve and geometrical
         parameters and store results.
@@ -230,6 +243,8 @@ class Function. Thrust units are Newtons.
             positions specified. Options are "nozzle_to_combustion_chamber"
             and "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         """
         super().__init__(
             thrust_source=thrust_source,
@@ -242,6 +257,7 @@ class Function. Thrust units are Newtons.
             reshape_thrust_curve=reshape_thrust_curve,
             interpolation_method=interpolation_method,
             coordinate_system_orientation=coordinate_system_orientation,
+            reference_pressure=reference_pressure,
         )
 
         self.positioned_tanks = []
@@ -264,7 +280,8 @@ def exhaust_velocity(self):
         -----
         The exhaust velocity is computed as the ratio of the thrust and the
         mass flow rate. Therefore, this will vary with time if the mass flow
-        rate varies with time.
+        rate varies with time. This corresponds to the actual exhaust velocity
+        only when the nozzle exit pressure equals the atmospheric pressure.
         """
         times, thrusts = self.thrust.source[:, 0], self.thrust.source[:, 1]
         mass_flow_rates = self.mass_flow_rate(times)
@@ -511,6 +528,7 @@ def from_dict(cls, data):
             nozzle_position=data["nozzle_position"],
             interpolation_method=data["interpolate"],
             coordinate_system_orientation=data["coordinate_system_orientation"],
+            reference_pressure=data["reference_pressure"],
         )
 
         for tank in data["positioned_tanks"]:
diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py
index 6a12665dd..ebc39b05a 100644
--- a/rocketpy/motors/motor.py
+++ b/rocketpy/motors/motor.py
@@ -29,6 +29,8 @@ class Motor(ABC):
         "combustion_chamber_to_nozzle".
     Motor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    Motor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     Motor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See :ref:`positions` for more information.
@@ -122,7 +124,11 @@ class Motor(ABC):
         e_3 axes in kg*m^2, as a function of time. See Motor.propellant_I_22
         and Motor.propellant_I_33 for more information.
     Motor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from the thrust source, in Newtons, as a
+        function of time.
+    Motor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     Motor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     Motor.max_thrust : float
@@ -142,11 +148,17 @@ class Motor(ABC):
         Total motor burn duration, in seconds. It is the difference between
         the burn_out_time and the burn_start_time.
     Motor.exhaust_velocity : Function
-        Propulsion gases exhaust velocity in m/s.
+        Effective exhaust velocity of the propulsion gases in m/s. Computed
+        as the thrust divided by the mass flow rate. This corresponds to the
+        actual exhaust velocity only when the nozzle exit pressure equals the
+        atmospheric pressure.
     Motor.interpolate : string
         Method of interpolation used in case thrust curve is given
         by data set in .csv or .eng, or as an array. Options are 'spline'
         'akima' and 'linear'. Default is "linear".
+    Motor.reference_pressure : int, float, None
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     # pylint: disable=too-many-statements
@@ -162,6 +174,7 @@ def __init__(
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Initialize Motor class, process thrust curve and geometrical
         parameters and store results.
@@ -237,6 +250,8 @@ class Function. Thrust units are Newtons.
             positions specified. Options are "nozzle_to_combustion_chamber" and
             "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
 
         Returns
         -------
@@ -258,7 +273,9 @@ class Function. Thrust units are Newtons.
         self.interpolate = interpolation_method
         self.nozzle_position = nozzle_position
         self.nozzle_radius = nozzle_radius
+        self.nozzle_area = np.pi * nozzle_radius**2
         self.center_of_dry_mass_position = center_of_dry_mass_position
+        self.reference_pressure = reference_pressure
 
         # Inertia tensor setup
         inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia
@@ -314,10 +331,6 @@ class Function. Thrust units are Newtons.
         self.burn_out_time = self.burn_time[1]
         self.burn_duration = self.burn_time[1] - self.burn_time[0]
 
-        # Define motor attributes
-        self.nozzle_radius = nozzle_radius
-        self.nozzle_position = nozzle_position
-
         # Compute thrust metrics
         self.max_thrust = np.amax(self.thrust.y_array)
         max_thrust_index = np.argmax(self.thrust.y_array)
@@ -410,7 +423,7 @@ def total_impulse(self):
     @property
     @abstractmethod
     def exhaust_velocity(self):
-        """Exhaust velocity of the motor gases.
+        """Effective exhaust velocity of the motor gases.
 
         Returns
         -------
@@ -429,6 +442,9 @@ def exhaust_velocity(self):
         - The ``LiquidMotor`` class favors the more accurate data from the
           Tanks's mass flow rates. Therefore the exhaust velocity is generally
           variable, being the ratio of the motor thrust by the mass flow rate.
+
+        This corresponds to the actual exhaust velocity only when the nozzle
+        exit pressure equals the atmospheric pressure.
         """
 
     @funcify_method("Time (s)", "Total mass (kg)")
@@ -1127,6 +1143,45 @@ def import_eng(file_name):
         # Return all extract content
         return comments, description, data_points
 
+    @cached_property
+    def vacuum_thrust(self):
+        """Calculate the vacuum thrust from the raw thrust and the reference
+        pressure at which the thrust data was recorded.
+
+        Returns
+        -------
+        vacuum_thrust : Function
+            The rocket's thrust in a vaccum.
+        """
+        if self.reference_pressure is None:
+            warnings.warn(
+                "Reference pressure not set. Returning thrust instead.",
+                UserWarning,
+            )
+            return self.thrust
+
+        return self.thrust + self.reference_pressure * self.nozzle_area
+
+    def pressure_thrust(self, pressure):
+        """Computes the contribution to thrust due to the difference between
+        the atmospheric pressure and the reference pressure at which the
+        thrust data was recorded.
+
+        Parameters
+        ----------
+        pressure : float
+            Atmospheric pressure in Pa.
+
+        Returns
+        -------
+        pressure_thrust : float
+            Thrust component resulting from the pressure difference.
+        """
+        if self.reference_pressure is None:
+            return 0
+
+        return (self.reference_pressure - pressure) * self.nozzle_area
+
     def export_eng(self, file_name, motor_name):
         """Exports thrust curve data points and motor description to
         .eng file format. A description of the format can be found
@@ -1193,17 +1248,20 @@ def to_dict(self, include_outputs=False):
             "dry_I_13": self.dry_I_13,
             "dry_I_23": self.dry_I_23,
             "nozzle_radius": self.nozzle_radius,
+            "nozzle_area": self.nozzle_area,
             "center_of_dry_mass_position": self.center_of_dry_mass_position,
             "dry_mass": self.dry_mass,
             "nozzle_position": self.nozzle_position,
             "burn_time": self.burn_time,
             "interpolate": self.interpolate,
             "coordinate_system_orientation": self.coordinate_system_orientation,
+            "reference_pressure": self.reference_pressure,
         }
 
         if include_outputs:
             data.update(
                 {
+                    "vacuum_thrust": self.vacuum_thrust,
                     "total_mass": self.total_mass,
                     "propellant_mass": self.propellant_mass,
                     "mass_flow_rate": self.mass_flow_rate,
@@ -1264,6 +1322,7 @@ class GenericMotor(Motor):
     therefore for more accurate results, use the ``SolidMotor``, ``HybridMotor``
     or ``LiquidMotor`` classes."""
 
+    # pylint: disable=too-many-arguments
     def __init__(
         self,
         thrust_source,
@@ -1280,6 +1339,7 @@ def __init__(
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Initialize GenericMotor class, process thrust curve and geometrical
         parameters and store results.
@@ -1369,6 +1429,8 @@ def __init__(
             positions specified. Options are "nozzle_to_combustion_chamber" and
             "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
         """
         super().__init__(
             thrust_source=thrust_source,
@@ -1381,6 +1443,7 @@ def __init__(
             reshape_thrust_curve=reshape_thrust_curve,
             interpolation_method=interpolation_method,
             coordinate_system_orientation=coordinate_system_orientation,
+            reference_pressure=reference_pressure,
         )
 
         self.chamber_radius = chamber_radius
@@ -1418,6 +1481,11 @@ def exhaust_velocity(self):
         -------
         self.exhaust_velocity : Function
             Gas exhaust velocity of the motor.
+
+        Notes
+        -----
+        This corresponds to the actual exhaust velocity only when the nozzle
+        exit pressure equals the atmospheric pressure.
         """
         return Function(
             self.total_impulse / self.propellant_initial_mass
@@ -1542,6 +1610,7 @@ def load_from_eng_file(
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Loads motor data from a .eng file and processes it.
 
@@ -1603,6 +1672,8 @@ def load_from_eng_file(
             positions specified. Options are "nozzle_to_combustion_chamber" and
             "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
 
         Returns
         -------
@@ -1652,6 +1723,7 @@ def load_from_eng_file(
             reshape_thrust_curve=reshape_thrust_curve,
             interpolation_method=interpolation_method,
             coordinate_system_orientation=coordinate_system_orientation,
+            reference_pressure=reference_pressure,
         )
 
     @staticmethod
@@ -1824,4 +1896,5 @@ def from_dict(cls, data):
             ),
             nozzle_position=data["nozzle_position"],
             interpolation_method=data["interpolate"],
+            reference_pressure=data["reference_pressure"],
         )
diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py
index 37d89130e..e736831bb 100644
--- a/rocketpy/motors/solid_motor.py
+++ b/rocketpy/motors/solid_motor.py
@@ -29,6 +29,8 @@ class SolidMotor(Motor):
         "combustion_chamber_to_nozzle".
     SolidMotor.nozzle_radius : float
         Radius of motor nozzle outlet in meters.
+    SolidMotor.nozzle_area : float
+        Area of motor nozzle outlet in square meters.
     SolidMotor.nozzle_position : float
         Motor's nozzle outlet position in meters, specified in the motor's
         coordinate system. See
@@ -147,7 +149,11 @@ class SolidMotor(Motor):
         See SolidMotor.propellant_I_22 and SolidMotor.propellant_I_33 for more
         information.
     SolidMotor.thrust : Function
-        Motor thrust force, in Newtons, as a function of time.
+        Motor thrust force obtained from thrust source, in Newtons, as a
+        function of time.
+    SolidMotor.vacuum_thrust : Function
+        Motor thrust force when the rocket is in a vacuum. In Newtons, as a
+        function of time.
     SolidMotor.total_impulse : float
         Total impulse of the thrust curve in N*s.
     SolidMotor.max_thrust : float
@@ -167,7 +173,10 @@ class SolidMotor(Motor):
         Total motor burn duration, in seconds. It is the difference between the
         ``burn_out_time`` and the ``burn_start_time``.
     SolidMotor.exhaust_velocity : Function
-        Propulsion gases exhaust velocity, assumed constant, in m/s.
+        Effective exhaust velocity of the propulsion gases in m/s. Computed
+        as the thrust divided by the mass flow rate. This corresponds to the
+        actual exhaust velocity only when the nozzle exit pressure equals the
+        atmospheric pressure.
     SolidMotor.burn_area : Function
         Total burn area considering all grains, made out of inner
         cylindrical burn area and grain top and bottom faces. Expressed
@@ -181,6 +190,9 @@ class SolidMotor(Motor):
         Method of interpolation used in case thrust curve is given
         by data set in .csv or .eng, or as an array. Options are 'spline'
         'akima' and 'linear'. Default is "linear".
+    SolidMotor.reference_pressure : int, float
+        Atmospheric pressure in Pa at which the thrust data was recorded.
+        It will allow to obtain the net thrust in the Flight class.
     """
 
     # pylint: disable=too-many-arguments
@@ -204,6 +216,7 @@ def __init__(
         reshape_thrust_curve=False,
         interpolation_method="linear",
         coordinate_system_orientation="nozzle_to_combustion_chamber",
+        reference_pressure=None,
     ):
         """Initialize Motor class, process thrust curve and geometrical
         parameters and store results.
@@ -299,6 +312,8 @@ class Function. Thrust units are Newtons.
             positions specified. Options are "nozzle_to_combustion_chamber" and
             "combustion_chamber_to_nozzle". Default is
             "nozzle_to_combustion_chamber".
+        reference_pressure : int, float, optional
+            Atmospheric pressure in Pa at which the thrust data was recorded.
 
         Returns
         -------
@@ -315,6 +330,7 @@ class Function. Thrust units are Newtons.
             reshape_thrust_curve=reshape_thrust_curve,
             interpolation_method=interpolation_method,
             coordinate_system_orientation=coordinate_system_orientation,
+            reference_pressure=reference_pressure,
         )
         # Nozzle parameters
         self.throat_radius = throat_radius
@@ -379,6 +395,11 @@ def exhaust_velocity(self):
         -------
         self.exhaust_velocity : Function
             Gas exhaust velocity of the motor.
+
+        Notes
+        -----
+        This corresponds to the actual exhaust velocity only when the nozzle
+        exit pressure equals the atmospheric pressure.
         """
         return Function(
             self.total_impulse / self.propellant_initial_mass
@@ -800,4 +821,5 @@ def from_dict(cls, data):
             throat_radius=data["throat_radius"],
             interpolation_method=data["interpolate"],
             coordinate_system_orientation=data["coordinate_system_orientation"],
+            reference_pressure=data["reference_pressure"],
         )
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 6c0e14312..3b857600a 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -358,6 +358,11 @@ class Flight:
         as array.
         Direction 3 is in the rocket's body axis and points in the
         direction of cylindrical symmetry.
+    Flight.net_thrust : Function
+        Rocket's engine net thrust as a function of time in Newton.
+        This is the actual thrust force experienced by the rocket.
+        It may be corrected with the atmospheric pressure if a reference
+        pressure is defined. Can be called or accessed as array.
     Flight.aerodynamic_lift : Function
         Resultant force perpendicular to rockets axis due to
         aerodynamic effects as a function of time. Units in N.
@@ -1376,12 +1381,17 @@ def udot_rail1(self, t, u, post_processing=False):
         drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
 
         # Calculate Forces
-        thrust = self.rocket.motor.thrust.get_value_opt(t)
+        pressure = self.env.pressure.get_value_opt(z)
+        net_thrust = max(
+            self.rocket.motor.thrust.get_value_opt(t)
+            + self.rocket.motor.pressure_thrust(pressure),
+            0,
+        )
         rho = self.env.density.get_value_opt(z)
         R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff)
 
         # Calculate Linear acceleration
-        a3 = (R3 + thrust) / total_mass_at_t - (
+        a3 = (R3 + net_thrust) / total_mass_at_t - (
             e0**2 - e1**2 - e2**2 + e3**2
         ) * self.env.gravity.get_value_opt(z)
         if a3 > 0:
@@ -1450,6 +1460,8 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u
         # Determine lift force and moment
         R1, R2, M1, M2, M3 = 0, 0, 0, 0, 0
+        # Thrust correction parameters
+        pressure = self.env.pressure.get_value_opt(z)
         # Determine current behavior
         if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time:
             # Motor burning
@@ -1467,10 +1479,14 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
             mass_flow_rate_at_t = self.rocket.motor.mass_flow_rate.get_value_opt(t)
             propellant_mass_at_t = self.rocket.motor.propellant_mass.get_value_opt(t)
             # Thrust
-            thrust = self.rocket.motor.thrust.get_value_opt(t)
+            net_thrust = max(
+                self.rocket.motor.thrust.get_value_opt(t)
+                + self.rocket.motor.pressure_thrust(pressure),
+                0,
+            )
             # Off center moment
-            M1 += self.rocket.thrust_eccentricity_y * thrust
-            M2 -= self.rocket.thrust_eccentricity_x * thrust
+            M1 += self.rocket.thrust_eccentricity_y * net_thrust
+            M2 -= self.rocket.thrust_eccentricity_x * net_thrust
         else:
             # Motor stopped
             # Inertias
@@ -1483,7 +1499,7 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
             # Mass
             mass_flow_rate_at_t, propellant_mass_at_t = 0, 0
             # thrust
-            thrust = 0
+            net_thrust = 0
 
         # Retrieve important quantities
         # Inertias
@@ -1686,7 +1702,7 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
                 + 2 * c * mass_flow_rate_at_t * omega1
             )
             / total_mass_at_t,
-            (R3 - b * propellant_mass_at_t * (alpha2 - omega1 * omega3) + thrust)
+            (R3 - b * propellant_mass_at_t * (alpha2 - omega1 * omega3) + net_thrust)
             / total_mass_at_t,
         ]
         ax, ay, az = K @ Vector(L)
@@ -1711,7 +1727,22 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
 
         if post_processing:
             self.__post_processed_variables.append(
-                [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3]
+                [
+                    t,
+                    ax,
+                    ay,
+                    az,
+                    alpha1,
+                    alpha2,
+                    alpha3,
+                    R1,
+                    R2,
+                    R3,
+                    M1,
+                    M2,
+                    M3,
+                    net_thrust,
+                ]
             )
 
         return u_dot
@@ -1789,8 +1820,15 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         free_stream_mach = free_stream_speed / speed_of_sound
 
         if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time:
+            pressure = self.env.pressure.get_value_opt(z)
+            net_thrust = max(
+                self.rocket.motor.thrust.get_value_opt(t)
+                + self.rocket.motor.pressure_thrust(pressure),
+                0,
+            )
             drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach)
         else:
+            net_thrust = 0
             drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach)
         R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff
         for air_brakes in self.rocket.air_brakes:
@@ -1853,14 +1891,13 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
             M3 += L
 
         # Off center moment
-        thrust = self.rocket.motor.thrust.get_value_opt(t)
         M1 += (
             self.rocket.cp_eccentricity_y * R3
-            + self.rocket.thrust_eccentricity_y * thrust
+            + self.rocket.thrust_eccentricity_y * net_thrust
         )
         M2 -= (
             self.rocket.cp_eccentricity_x * R3
-            + self.rocket.thrust_eccentricity_x * thrust
+            + self.rocket.thrust_eccentricity_x * net_thrust
         )
         M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1
 
@@ -1871,7 +1908,7 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         T00 = total_mass * r_CM
         T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot
         T04 = (
-            Vector([0, 0, thrust])
+            Vector([0, 0, net_thrust])
             - total_mass * r_CM_ddot
             - 2 * total_mass_dot * r_CM_dot
             + total_mass_ddot * (r_NOZ - r_CM)
@@ -1915,7 +1952,7 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
 
         if post_processing:
             self.__post_processed_variables.append(
-                [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3]
+                [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust]
             )
 
         return u_dot
@@ -1988,7 +2025,7 @@ def u_dot_parachute(self, t, u, post_processing=False):
 
         if post_processing:
             self.__post_processed_variables.append(
-                [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0]
+                [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0, 0]
             )
 
         return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
@@ -2208,6 +2245,13 @@ def M3(self):
         Sometimes referred to as roll moment."""
         return self.__evaluate_post_process[:, [0, 12]]
 
+    @funcify_method("Time (s)", "Net Thrust (N)", "linear", "zero")
+    def net_thrust(self):
+        """Net thrust of the rocket as a Function of time. This is the
+        actual thrust force experienced by the rocket. It may be corrected
+        with the atmospheric pressure if a reference pressure is defined."""
+        return self.__evaluate_post_process[:, [0, 13]]
+
     @funcify_method("Time (s)", "Pressure (Pa)", "spline", "constant")
     def pressure(self):
         """Air pressure felt by the rocket as a Function of time."""
@@ -2647,13 +2691,10 @@ def total_energy(self):
         return self.kinetic_energy + self.potential_energy
 
     # thrust Power
-    @funcify_method("Time (s)", "thrust Power (W)", "spline", "zero")
+    @funcify_method("Time (s)", "Thrust Power (W)", "spline", "zero")
     def thrust_power(self):
         """Thrust power as a Function of time."""
-        thrust = deepcopy(self.rocket.motor.thrust)
-        thrust = thrust.set_discrete_based_on_model(self.speed)
-        thrust_power = thrust * self.speed
-        return thrust_power
+        return self.net_thrust * self.speed
 
     # Drag Power
     @funcify_method("Time (s)", "Drag Power (W)", "spline", "zero")
@@ -3059,7 +3100,7 @@ def __evaluate_post_process(self):
         np.array
             An array containing all post-processed variables evaluated at each
             time step. Each element of the array is a list containing:
-            [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3]
+            [t, ax, ay, az, alpha1, alpha2, alpha3, R1, R2, R3, M1, M2, M3, net_thrust]
         """
         self.__post_processed_variables = []
         for phase_index, phase in self.time_iterator(self.flight_phases):