diff --git a/python/amun.py b/python/amun.py
new file mode 100644
index 0000000..0e139cb
--- /dev/null
+++ b/python/amun.py
@@ -0,0 +1,204 @@
+"""
+================================================================================
+
+ This file is part of the AMUN source code, a program to perform
+ Newtonian or relativistic magnetohydrodynamical simulations on uniform or
+ adaptive mesh.
+
+ Copyright (C) 2018 Grzegorz Kowal <grzegorz@amuncode.org>
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+================================================================================
+
+ module: AMUN
+
+ Python module with subroutines to read AMUN code HDF5 files.
+
+ The only requirements for this module are:
+ - h5py
+ - numpy
+
+--------------------------------------------------------------------------------
+"""
+import h5py as h5
+import numpy as np
+import os.path as op
+
+def amun_attribute(fname, aname):
+ '''
+ Subroutine to reads global attributes from AMUN HDF5 snapshots.
+
+ Arguments:
+
+ fname - the AMUN HDF5 file name;
+ aname - the attribute name;
+
+ Return values:
+
+ ret - the value read from the attribute;
+
+ Examples:
+
+ time = amun_attribute('p000010_00000.h5', 'time')
+
+ '''
+ # open the file
+ #
+ f = h5.File(fname, 'r')
+
+ # check if the file is written in the AMUN format
+ #
+ if f.attrs.get('code')[0].astype(str) != "AMUN" or \
+ not 'attributes' in f or \
+ not 'coordinates' in f or \
+ not 'variables' in f:
+ print('It seems this HDF5 file is corrupted or not compatible with the AMUN format!')
+ f.close()
+ return False
+
+ # open the group of attributes
+ #
+ g = f['attributes'].attrs
+
+ # get attribute's value
+ #
+ ret = g.get(aname)
+ if len(ret) == 1:
+ ret = ret[0]
+
+ # close the file
+ #
+ f.close()
+
+ # return the value of attribute
+ #
+ return ret
+
+def amun_dataset(fname, vname):
+ '''
+ Subroutine to reads dataset from AMUN HDF5 snapshots.
+
+ Arguments:
+
+ fname - the AMUN HDF5 file name;
+ vname - the variable name;
+
+ Return values:
+
+ ret - the array of values for the variable;
+
+ Examples:
+
+ dn = amun_dataset('p000010_00000.h5', 'dens')
+
+ '''
+ # open the file
+ #
+
+ f = h5.File(fname, 'r')
+
+ # check if the file is written in the AMUN format
+ #
+ if f.attrs.get('code')[0].astype(str) != "AMUN" or \
+ not 'attributes' in f or \
+ not 'coordinates' in f or \
+ not 'variables' in f:
+ print('It seems this HDF5 file is corrupted or not compatible with the AMUN format!')
+ return False
+
+ # build the list of supported variables
+ #
+ variables = []
+ for var in f['variables'].keys():
+ variables.append(var)
+
+ # check if the requested variable is in the variable list
+ #
+ if not vname in variables:
+ print('The requested variable cannot be extracted from the file datasets!')
+ return False
+
+ # get attributes necessary to reconstruct the domain
+ #
+ g = f['attributes'].attrs
+
+ # get the snapshot number, the number of domain files, and the number of blocks
+ #
+ nr = g.get('isnap')[0]
+ nc = g.get('nprocs')[0]
+ nl = g.get('nleafs')[0]
+
+ # prepare array to hold data
+ #
+ bm = g.get('dims')
+ if bm[2] > 1:
+ ndims = 3
+ else:
+ ndims = 2
+ rm = g.get('rdims')
+ ng = g.get('nghosts')
+ ml = g.get('maxlev')[0]
+ dm = rm[0:ndims] * bm[0:ndims] * 2**(ml - 1)
+ ret = np.zeros(dm[::-1])
+
+ f.close()
+
+ # iterate over all subdomain files
+ #
+ dname = op.dirname(fname)
+ for n in range(nc):
+ lname = op.join(dname, 'p%06d_%05d.h5' % (nr, n))
+ f = h5.File(lname, 'r')
+ g = f['attributes'].attrs
+ dblocks = g.get('dblocks')[0]
+ if dblocks > 0:
+ g = f['coordinates']
+ levels = g['levels'][()]
+ coords = g['coords'][()]
+ g = f['variables']
+ if ndims == 3:
+ dataset = g[vname][ng:-ng,ng:-ng,ng:-ng,:]
+ for l in range(dblocks):
+ nn = 2**(ml - levels[l])
+ cm = bm * nn
+ ibeg = coords[:,l] * cm[:]
+ iend = ibeg + cm[:]
+ ib, jb, kb = ibeg[0], ibeg[1], ibeg[2]
+ ie, je, ke = iend[0], iend[1], iend[2]
+ ret[kb:ke,jb:je,ib:ie] = np.repeat(np.repeat(np.repeat(dataset[:,:,:,l], nn, axis = 0), nn, axis = 1), nn, axis = 2)
+ else:
+ dataset = g[vname][0,ng:-ng,ng:-ng,:]
+ for l in range(dblocks):
+ nn = 2**(ml - levels[l])
+ cm = bm[0:ndims] * nn
+ ibeg = coords[:,l] * cm[:]
+ iend = ibeg + cm[:]
+ ib, jb = ibeg[0], ibeg[1]
+ ie, je = iend[0], iend[1]
+ ret[jb:je,ib:ie] = np.repeat(np.repeat(dataset[:,:,l], nn, axis = 0), nn, axis = 1)
+ f.close()
+
+ return ret
+
+if __name__ == "__main__":
+ fname = './p000030_00000.h5'
+
+ ret = amun_attribute(fname, 'time')
+ print(ret)
+ ret = amun_attribute(fname, 'dims')
+ print(ret)
+
+ ret = amun_dataset(fname, 'dens')
+ print(ret.shape, ret.min(), ret.max())