diff --git a/CHANGES.rst b/CHANGES.rst
index 0f4cb12cd..f88226eb0 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -41,6 +41,7 @@ Breaking changes
 * `black` formatting style has been updated to the 2024 stable conventions. `isort` has been added to the `dev` installation recipe. (:pull:`1626`).
 * The indice and indicator for ``winter_storm`` has been removed (deprecated since `xclim` v0.46.0 in favour of ``snd_storm_days``). (:pull:`1565`).
 * `xclim` has dropped support for `scipy` version below v1.9.0 and `numpy` versions below v1.20.0. (:pull:`1565`).
+* For generic function ``select_resample_op`` and ``core.units.to_agg_units``, operation "sum" will now return the same units as the input, and not implicitly be translated to an "integral". (:issue:`1645`, :pull:`1649`).
 
 Bug fixes
 ^^^^^^^^^
diff --git a/tests/test_atmos.py b/tests/test_atmos.py
index 214e0dbfc..1d5a6e6fc 100644
--- a/tests/test_atmos.py
+++ b/tests/test_atmos.py
@@ -280,7 +280,7 @@ def test_wind_power_potential_from_3h_series():
     power = out * 100
     power.attrs["units"] = "MW"
     annual_power = convert_units_to(
-        select_resample_op(power, op="sum", freq="D"), "MWh"
+        select_resample_op(power, op="integral", freq="D"), "MWh"
     )
     assert (annual_power == 100 * 24).all()
 
diff --git a/tests/test_units.py b/tests/test_units.py
index 5dac52545..febcd27be 100644
--- a/tests/test_units.py
+++ b/tests/test_units.py
@@ -22,6 +22,7 @@
     pint_multiply,
     rate2amount,
     str2pint,
+    to_agg_units,
     units,
     units2pint,
 )
@@ -317,3 +318,28 @@ def index(data: xr.DataArray, thresh: Quantified, dthreshdt: Quantified):
 
     with pytest.raises(ValidationError):
         index_full_mm("1 mm", "2 Pa", "3 mm/s")
+
+
+@pytest.mark.parametrize(
+    "in_u,opfunc,op,exp,exp_u",
+    [
+        ("m/h", "sum", "integral", 8760, "m"),
+        ("m/h", "sum", "sum", 365, "m/h"),
+        ("K", "mean", "mean", 1, "K"),
+        ("", "sum", "count", 365, "d"),
+        ("", "sum", "count", 365, "d"),
+        ("kg m-2", "var", "var", 0, "kg2 m-4"),
+        ("°C", "argmax", "doymax", 0, ""),
+    ],
+)
+def test_to_agg_units(in_u, opfunc, op, exp, exp_u):
+    da = xr.DataArray(
+        np.ones((365,)),
+        dims=("time",),
+        coords={"time": xr.cftime_range("1993-01-01", periods=365, freq="D")},
+        attrs={"units": in_u},
+    )
+
+    out = to_agg_units(getattr(da, opfunc)(), da, op)
+    np.testing.assert_allclose(out, exp)
+    assert out.attrs["units"] == exp_u
diff --git a/xclim/core/units.py b/xclim/core/units.py
index c6a73c91d..2d58a33f8 100644
--- a/xclim/core/units.py
+++ b/xclim/core/units.py
@@ -471,8 +471,8 @@ def to_agg_units(
         The original array before the aggregation operation,
         used to infer the sampling units and get the variable units.
     op : {'min', 'max', 'mean', 'std', 'var', 'doymin', 'doymax',  'count', 'integral', 'sum'}
-        The type of aggregation operation performed. The special "delta_*" ops are used
-        with temperature units needing conversion to their "delta" counterparts (e.g. degree days)
+        The type of aggregation operation performed. "integral" is mathematically equivalent to "sum",
+        but the units are multiplied by the timestep of the data (requires an inferrable frequency).
     dim : str
         The time dimension along which the aggregation was performed.
 
@@ -521,7 +521,7 @@ def to_agg_units(
     >>> degdays.units
     'K d'
     """
-    if op in ["amin", "min", "amax", "max", "mean", "std"]:
+    if op in ["amin", "min", "amax", "max", "mean", "std", "sum"]:
         out.attrs["units"] = orig.attrs["units"]
 
     elif op in ["var"]:
@@ -532,7 +532,7 @@ def to_agg_units(
             units="", is_dayofyear=np.int32(1), calendar=get_calendar(orig)
         )
 
-    elif op in ["count", "integral", "sum"]:
+    elif op in ["count", "integral"]:
         m, freq_u_raw = infer_sampling_units(orig[dim])
         orig_u = str2pint(orig.units)
         freq_u = str2pint(freq_u_raw)
@@ -540,8 +540,14 @@ def to_agg_units(
 
         if op == "count":
             out.attrs["units"] = freq_u_raw
-        elif op in ["integral", "sum"]:
-            out.attrs["units"] = pint2cfunits(orig_u * freq_u)
+        elif op == "integral":
+            if "[time]" in orig_u.dimensionality:
+                # We need to simplify units after multiplication
+                out_units = (orig_u * freq_u).to_reduced_units()
+                out = out * out_units.magnitude
+                out.attrs["units"] = pint2cfunits(out_units)
+            else:
+                out.attrs["units"] = pint2cfunits(orig_u * freq_u)
     else:
         raise ValueError(
             f"Unknown aggregation op {op}. "