diff --git a/examples/example_manual_aggregate.py b/examples/example_manual_aggregate.py
new file mode 100755
index 000000000..48d84ba4c
--- /dev/null
+++ b/examples/example_manual_aggregate.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+"""
+Manual aggregate
+==================
+"""
+
+from datetime import datetime, timedelta
+import pandas as pd
+import matplotlib.pyplot as plt
+from opendrift.readers import reader_netCDF_CF_generic
+from opendrift.readers import open_mfdataset_overlap
+from opendrift.models.oceandrift import OceanDrift
+
+#%
+# Create manual aggregate from individual URLs, for NorKyst ocean model initialized at 00 hours every day
+ds = open_mfdataset_overlap(
+    'https://thredds.met.no/thredds/dodsC/fou-hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.%Y%m%d%H.nc',
+    time_series=pd.date_range(datetime.now().date()-timedelta(days=1), datetime.now().date(), freq='1D'))
+
+#%
+# Create reader from Xarray dataset
+rm = reader_netCDF_CF_generic.Reader(ds, name='NorKyst manual aggregate')
+print(rm)
+
+om = OceanDrift()
+om.add_reader(rm)
+om.seed_elements(lon=3.4, lat=60.23, number=1000, radius=100, time=rm.start_time)
+om.run(end_time=rm.end_time)
+
+#%
+# Second simulation using ready made aggregate from thredds
+ot = OceanDrift()
+ot.add_readers_from_list(['https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be'])
+ot.seed_elements(lon=3.4, lat=60.23, number=1000, radius=100, time=rm.start_time)
+ot.run(end_time=rm.end_time)
+
+#%
+# Simulation should be identical, but we see that manual aggregate is significantly slower than using thredds aggregate
+om.animation(compare=ot,
+             legend=[f'NorKyst manual aggregate {om.timing["total time"]}',
+                     f'NorKyst thredds aggregate {ot.timing["total time"]}'])
+
+#%%
+# .. image:: /gallery/animations/example_manual_aggregate_0.gif
diff --git a/opendrift/readers/__init__.py b/opendrift/readers/__init__.py
index 7a8f52d60..9f6f268c1 100644
--- a/opendrift/readers/__init__.py
+++ b/opendrift/readers/__init__.py
@@ -20,11 +20,26 @@
     See :class:`.basereader.BaseReader` for how readers work internally.
 """
 
+from datetime import datetime, timedelta
 import importlib
 import logging; logger = logging.getLogger(__name__)
 import glob
 import json
 import opendrift
+import xarray as xr
+
+def open_mfdataset_overlap(url_base, time_series=None, start_time=None, end_time=None, freq=None, timedim='time'):
+    if time_series is None:
+        construct_from_times
+    urls = [t.strftime(url_base) for t in time_series]
+    time_step = time_series[1] - time_series[0]
+    print('Opening individual URLs...')
+    datasets = [xr.open_dataset(u, chunks='auto').sel({timedim: slice(t, t+time_step-timedelta(seconds=1))})
+                for u,t in zip(urls, time_series)]
+    print('Concatenating...')
+    ds = xr.concat(datasets, dim=timedim,
+                   compat='override', combine_attrs='override', join='override', coords='all')
+    return ds
 
 def reader_from_url(url, timeout=10):
     '''Make readers from URLs or paths to datasets'''