diff --git a/CHANGELOG.md b/CHANGELOG.md
index fe957ca..107ba3b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,10 @@
 * [Issue 282](https://github.com/MassimoCimmino/pygfunction/issues/282) - Enabled the use of negative mass flow rates in `Pipe` and `Network` classes to model reversed flow direction.
 * [Pull Request 308](https://github.com/MassimoCimmino/pygfunction/pull/308) - Introduced a new `borefield` module. The new `Borefield` class replaces lists of `Borehole` objects as the preferred way to configure bore fields. The `Borefield.evaluate_g_function` method evaluates g-functions using the 'UHTR' and 'UBWT' boundary conditions. Deprecated bore field creation functions in the `boreholes` module (e.g. `boreholes.rectangle_field()`). These functions are replaced by methods of the new `Borefield` class. They will be removed in `v3.0.0`.
 
+### Bug fixes
+
+* [Issue 305](https://github.com/MassimoCimmino/pygfunction/issues/305) - Fixed `ClaessonJaved` to return a float when the *g*-function is a vector (i.e. when there is only one heat source). This is required for compatibility with `numpy` version `2.x`.
+
 ### Other changes
 
 * [Issue 312](https://github.com/MassimoCimmino/pygfunction/issues/312) - The installation of `matplotlib` is now optional. Using `pip install pygfunction` will not install `matplotlib` and `pip install pygfunction[plot]` should be used instead.
diff --git a/pygfunction/load_aggregation.py b/pygfunction/load_aggregation.py
index 9416ece..b95a8fb 100644
--- a/pygfunction/load_aggregation.py
+++ b/pygfunction/load_aggregation.py
@@ -93,6 +93,10 @@ def initialize(self, g_d):
             If nSources=1, g_d can be 1 dimensional.
 
         """
+        if g_d.ndim == 1:
+            self._scalar_output = True
+        else:
+            self._scalar_output = False
         if self.nSources==1:
             g_d = g_d.reshape(1, 1, -1)
         # Build matrix of thermal response factor increments
@@ -187,6 +191,8 @@ def temporal_superposition(self):
         #        deltaT += (self.dg[:,:,i]).dot(self.q_b[:,i])
 
         deltaT = np.einsum('ijk,jk', self.dg, self.q_b)
+        if self._scalar_output:
+            deltaT = deltaT.item()
         return deltaT
 
     def _build_cells(self, dt, tmax, nSources, cells_per_level):