diff --git a/WRF.md b/WRF.md
index 0a196e680bb1443bfe7d194219329b9c944c3dc4..de2eb2d9011c5f7caddfea943ddac0f7418f4eda 100644
--- a/WRF.md
+++ b/WRF.md
@@ -12,6 +12,7 @@
     * [Running WRF in a software container](#)
     * [Running an idealized simulation](#)
     * [Running a real-case simulation](#)
+    * [Output and restart files](#)
     * [Suggested workflow](#)
     * [Analysing model output](#)
     * [Important namelist settings](#)
@@ -344,12 +345,60 @@ then you have a problem, and there is no unique solution. Take a closer look at
 
 ## Running a real-case simulation
 
+## Output and restart files
+
+incl. how to modify output paths
+
 ## Suggested workflow
 
 ## Analysing model output
 
+Things to remember:
+
+* staggered grid (Arakawa-C)
+* mass-based vertical coordinate (level height AGL is time-dependent)
+* terrain-following coordinate system (curvilinear)
+* in the model output, some variables are split into base state + perturbation
+
 [Python interface to WRF](https://wrf-python.readthedocs.io/en/latest/)  
 
+Example of a very basic Python class to create an object from a WRF run, initialized with only some basic information:
+```
+class wrfrun:
+    def __init__(self, filename):
+        self.filename = filename
+        self.nc = netCDF4.Dataset(filename)
+        self.dx = self.nc.DX
+        self.dy = self.nc.DY
+        self.nx = self.nc.dimensions['west_east'].size
+        self.ny = self.nc.dimensions['south_north'].size
+        self.x = np.arange(0,self.nx*self.dx,self.dx)
+        self.y = np.arange(0,self.ny*self.dy,self.dy)
+        self.valid_times = self.nc['XTIME'][:]*60
+        self.current_time = 0
+
+    def set_time(self,step):
+        self.current_time = step
+
+    def add_wind(self):
+        udum = self.nc['U'][self.current_time,:,:,:]
+        vdum = self.nc['V'][self.current_time,:,:,:]
+        wdum = self.nc['W'][self.current_time,:,:,:]
+        self.u = 0.5*(udum[:,:,:-1]+udum[:,:,1:])
+        self.v = 0.5*(vdum[:,:-1,:]+vdum[:,1:,:])
+        self.w = 0.5*(wdum[:-1,:,:]+wdum[1:,:,:])
+        del udum,vdum,wdum
+```
+The last function adds 3D wind variables at a specific time, after destaggering.
+
+The `wrfrun` class is then used as follows:
+```
+wrf = wrfrun("./wrfout_d01_0001-01-01_00:00:00")
+wrf.set_time(36)
+wrf.add_wind()
+```
+Variables are then accessible as `wrf.u`, `wrf.v` etc.
+
 ## Important namelist settings
 
 # Advanced usage