diff --git a/post_process/.gitignore b/post_process/.gitignore
index e2e7327..40906ba 100644
--- a/post_process/.gitignore
+++ b/post_process/.gitignore
@@ -1 +1,2 @@
+local.hdf
 /out
diff --git a/post_process/post_process/__main__.py b/post_process/post_process/__main__.py
new file mode 100644
index 0000000..674394d
--- /dev/null
+++ b/post_process/post_process/__main__.py
@@ -0,0 +1,115 @@
+from pprint import pp
+import logging
+import configparser
+from pathlib import Path
+import re
+from multiprocessing.pool import ThreadPool
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+import fluidfoam
+
+logging.basicConfig(level='INFO')
+log = logging.getLogger('post_process')
+
+config = configparser.ConfigParser()
+config.read('config.ini')
+
+root = Path(config.get('main', 'root')).expanduser()
+cases = config.get('main', 'cases').split(',')
+case_root = root.joinpath(cases[0])
+
+if config.getboolean('main', 'loaddata', fallback=False):
+    mesh = fluidfoam.readmesh(str(case_root))
+
+    t = []
+    for item in case_root.iterdir():
+        if item.is_dir() and re.fullmatch(r'[0-9]+(\.[0-9]+)?', item.name):
+            t.append(item.name)
+    timesteps = pd.Series(t, dtype='float', name='t', index=t)\
+        .sort_values()
+    timesteps.drop('0', inplace=True)
+    
+    data = {}
+    with pd.HDFStore(
+            Path(config.get('main', 'savefile')),
+            mode='w',
+            complib='blosc',
+        ) as hdf:
+        hdf.put('timesteps', timesteps)
+        hdf.put('mesh', pd.DataFrame(mesh, index=('x', 'y', 'z')).T)
+        for case in cases:
+            with ThreadPool(
+                        config.getint('main', 'nthreads', fallback=1)
+                    ) as pool:
+                sensor = dict(zip(timesteps, pool.map(
+                    lambda t:np.apply_along_axis(
+                        lambda x:np.sqrt((x**2).sum()),
+                        0,
+                        fluidfoam.readvector(
+                            root.joinpath(case),
+                            t,
+                            config.get('sensor', 'value')
+                        ),
+                    ),
+                    timesteps.index
+                )))
+        
+                hdf.put(case, pd.DataFrame(sensor))
+                del sensor
+
+with pd.HDFStore(
+            Path(config.get('bathy', 'source')),
+            mode='r',
+        ) as bathy:
+    rubble_line = bathy.get('rubble')
+
+
+with pd.HDFStore(
+            Path(config.get('main', 'savefile')),
+            mode='r',
+        ) as data:
+    timesteps = data.get('timesteps')
+    mesh = data.get('mesh')
+    rawxs = config.get('sensor', 'x').split(',')
+    fig, ax = plt.subplots(len(rawxs))
+    figmax = config.getfloat('figure', 'max')
+    for axi, rawx in zip(ax, rawxs):
+        x0 = float(rawx)
+        y0 = rubble_line.iloc[np.abs(x0-rubble_line.index).argmin()]
+        pp((x0,y0))
+
+        sensor_id = ((mesh.x-x0)**2+(mesh.y-y0)**2).argmin()
+        pp(mesh.iloc[sensor_id])
+
+        sensor_data = pd.DataFrame(index=timesteps)
+        for case, title, color, ls \
+                in zip(
+                    cases,
+                    config.get('main', 'case_titles').split(';'),
+                    config.get('main', 'case_colors').split(','),
+                    config.get('main', 'case_ls').split(','),
+                ):
+            sensor_data = pd.concat((
+                sensor_data,
+                data.get(case).iloc[sensor_id].rename(case),
+            ), axis=1)
+            axi.plot(
+                sensor_data.index,
+                sensor_data[case],
+                label=title,
+                color=color,
+                ls=ls,
+            )
+        axi.legend()
+        axi.grid()
+        axi.set(
+            xlabel='t (s)',
+            ylabel='U (m/s)',
+            xlim=(0,timesteps.max()),
+            ylim=(0,figmax),
+        )
+    fig.savefig(config.get('figure', 'save'))
+    plt.show()