diff --git a/bfps/PP.py b/bfps/PP.py
new file mode 100644
index 0000000000000000000000000000000000000000..ffa0b2da49b77a45ecd7d5a3e07a787a059c3fc8
--- /dev/null
+++ b/bfps/PP.py
@@ -0,0 +1,676 @@
+#######################################################################
+#                                                                     #
+#  Copyright 2015 Max Planck Institute                                #
+#                 for Dynamics and Self-Organization                  #
+#                                                                     #
+#  This file is part of bfps.                                         #
+#                                                                     #
+#  bfps 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.                             #
+#                                                                     #
+#  bfps 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 bfps.  If not, see <http://www.gnu.org/licenses/>       #
+#                                                                     #
+# Contact: Cristian.Lalescu@ds.mpg.de                                 #
+#                                                                     #
+#######################################################################
+
+
+
+import os
+import sys
+import shutil
+import subprocess
+import argparse
+import h5py
+import math
+import numpy as np
+import warnings
+
+import bfps
+from ._code import _code
+from bfps import tools
+
+class PP(_code):
+    """This class is meant to stitch together the C++ code into a final source file,
+    compile it, and handle all job launching.
+    """
+    def __init__(
+            self,
+            work_dir = './',
+            simname = 'test'):
+        _code.__init__(
+                self,
+                work_dir = work_dir,
+                simname = simname)
+        self.host_info = {'type'        : 'cluster',
+                          'environment' : None,
+                          'deltanprocs' : 1,
+                          'queue'       : '',
+                          'mail_address': '',
+                          'mail_events' : None}
+        self.generate_default_parameters()
+        return None
+    def set_precision(
+            self,
+            fluid_dtype):
+        if fluid_dtype in [np.float32, np.float64]:
+            self.fluid_dtype = fluid_dtype
+        elif fluid_dtype in ['single', 'double']:
+            if fluid_dtype == 'single':
+                self.fluid_dtype = np.dtype(np.float32)
+            elif fluid_dtype == 'double':
+                self.fluid_dtype = np.dtype(np.float64)
+        self.rtype = self.fluid_dtype
+        if self.rtype == np.float32:
+            self.ctype = np.dtype(np.complex64)
+            self.C_field_dtype = 'float'
+            self.fluid_precision = 'single'
+        elif self.rtype == np.float64:
+            self.ctype = np.dtype(np.complex128)
+            self.C_field_dtype = 'double'
+            self.fluid_precision = 'double'
+        return None
+    def write_src(self):
+        self.version_message = (
+                '/***********************************************************************\n' +
+                '* this code automatically generated by bfps\n' +
+                '* version {0}\n'.format(bfps.__version__) +
+                '***********************************************************************/\n\n\n')
+        self.include_list = [
+                '"base.hpp"',
+                '"scope_timer.hpp"',
+                '"fftw_interface.hpp"',
+                '"full_code/main_code.hpp"',
+                '<cmath>',
+                '<iostream>',
+                '<hdf5.h>',
+                '<string>',
+                '<cstring>',
+                '<fftw3-mpi.h>',
+                '<omp.h>',
+                '<cfenv>',
+                '<cstdlib>',
+                '"full_code/{0}.hpp"\n'.format(self.dns_type)]
+        self.main = """
+            int main(int argc, char *argv[])
+            {{
+                bool fpe = (
+                    (getenv("BFPS_FPE_OFF") == nullptr) ||
+                    (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+                return main_code< {0} >(argc, argv, fpe);
+            }}
+            """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+        self.includes = '\n'.join(
+                ['#include ' + hh
+                 for hh in self.include_list])
+        with open(self.name + '.cpp', 'w') as outfile:
+            outfile.write(self.version_message + '\n\n')
+            outfile.write(self.includes + '\n\n')
+            outfile.write(self.main + '\n')
+        return None
+    def generate_default_parameters(self):
+        # these parameters are relevant for all DNS classes
+        self.parameters['dealias_type'] = int(1)
+        self.parameters['dkx'] = float(1.0)
+        self.parameters['dky'] = float(1.0)
+        self.parameters['dkz'] = float(1.0)
+        self.parameters['nu'] = float(0.1)
+        self.parameters['fmode'] = int(1)
+        self.parameters['famplitude'] = float(0.5)
+        self.parameters['fk0'] = float(2.0)
+        self.parameters['fk1'] = float(4.0)
+        self.parameters['forcing_type'] = 'linear'
+        self.pp_parameters = {}
+        self.pp_parameters['iteration_list'] = np.zeros(1).astype(np.int)
+        return None
+    def get_data_file_name(self):
+        return os.path.join(self.work_dir, self.simname + '.h5')
+    def get_data_file(self):
+        return h5py.File(self.get_data_file_name(), 'r')
+    def get_particle_file_name(self):
+        return os.path.join(self.work_dir, self.simname + '_particles.h5')
+    def get_particle_file(self):
+        return h5py.File(self.get_particle_file_name(), 'r')
+    def get_postprocess_file_name(self):
+        return os.path.join(self.work_dir, self.simname + '_postprocess.h5')
+    def get_postprocess_file(self):
+        return h5py.File(self.get_postprocess_file_name(), 'r')
+    def compute_statistics(self, iter0 = 0, iter1 = None):
+        """Run basic postprocessing on raw data.
+        The energy spectrum :math:`E(t, k)` and the enstrophy spectrum
+        :math:`\\frac{1}{2}\omega^2(t, k)` are computed from the
+
+        .. math::
+
+            \sum_{k \\leq \\|\\mathbf{k}\\| \\leq k+dk}\\hat{u_i} \\hat{u_j}^*, \\hskip .5cm
+            \sum_{k \\leq \\|\\mathbf{k}\\| \\leq k+dk}\\hat{\omega_i} \\hat{\\omega_j}^*
+
+        tensors, and the enstrophy spectrum is also used to
+        compute the dissipation :math:`\\varepsilon(t)`.
+        These basic quantities are stored in a newly created HDF5 file,
+        ``simname_postprocess.h5``.
+        """
+        if len(list(self.statistics.keys())) > 0:
+            return None
+        self.read_parameters()
+        with self.get_data_file() as data_file:
+            if 'moments' not in data_file['statistics'].keys():
+                return None
+            iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
+                         self.parameters['niter_stat']-1),
+                        iter0)
+            if type(iter1) == type(None):
+                iter1 = data_file['iteration'].value
+            else:
+                iter1 = min(data_file['iteration'].value, iter1)
+            ii0 = iter0 // self.parameters['niter_stat']
+            ii1 = iter1 // self.parameters['niter_stat']
+            self.statistics['kshell'] = data_file['kspace/kshell'].value
+            self.statistics['kM'] = data_file['kspace/kM'].value
+            self.statistics['dk'] = data_file['kspace/dk'].value
+            computation_needed = True
+            pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
+            if 'ii0' in pp_file.keys():
+                computation_needed =  not (ii0 == pp_file['ii0'].value and
+                                           ii1 == pp_file['ii1'].value)
+                if computation_needed:
+                    for k in pp_file.keys():
+                        del pp_file[k]
+            if computation_needed:
+                pp_file['iter0'] = iter0
+                pp_file['iter1'] = iter1
+                pp_file['ii0'] = ii0
+                pp_file['ii1'] = ii1
+                pp_file['t'] = (self.parameters['dt']*
+                                self.parameters['niter_stat']*
+                                (np.arange(ii0, ii1+1).astype(np.float)))
+                pp_file['energy(t, k)'] = (
+                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 0, 0] +
+                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 1, 1] +
+                    data_file['statistics/spectra/velocity_velocity'][ii0:ii1+1, :, 2, 2])/2
+                pp_file['enstrophy(t, k)'] = (
+                    data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 0, 0] +
+                    data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 1, 1] +
+                    data_file['statistics/spectra/vorticity_vorticity'][ii0:ii1+1, :, 2, 2])/2
+                pp_file['vel_max(t)'] = data_file['statistics/moments/velocity']  [ii0:ii1+1, 9, 3]
+                pp_file['renergy(t)'] = data_file['statistics/moments/velocity'][ii0:ii1+1, 2, 3]/2
+            for k in ['t',
+                      'energy(t, k)',
+                      'enstrophy(t, k)',
+                      'vel_max(t)',
+                      'renergy(t)']:
+                if k in pp_file.keys():
+                    self.statistics[k] = pp_file[k].value
+            self.compute_time_averages()
+        return None
+    def compute_time_averages(self):
+        """Compute easy stats.
+
+        Further computation of statistics based on the contents of
+        ``simname_postprocess.h5``.
+        Standard quantities are as follows
+        (consistent with [Ishihara]_):
+
+        .. math::
+
+            U_{\\textrm{int}}(t) = \\sqrt{\\frac{2E(t)}{3}}, \\hskip .5cm
+            L_{\\textrm{int}}(t) = \\frac{\pi}{2U_{int}^2(t)} \\int \\frac{dk}{k} E(t, k), \\hskip .5cm
+            T_{\\textrm{int}}(t) =
+            \\frac{L_{\\textrm{int}}(t)}{U_{\\textrm{int}}(t)}
+
+            \\eta_K = \\left(\\frac{\\nu^3}{\\varepsilon}\\right)^{1/4}, \\hskip .5cm
+            \\tau_K = \\left(\\frac{\\nu}{\\varepsilon}\\right)^{1/2}, \\hskip .5cm
+            \\lambda = \\sqrt{\\frac{15 \\nu U_{\\textrm{int}}^2}{\\varepsilon}}
+
+            Re = \\frac{U_{\\textrm{int}} L_{\\textrm{int}}}{\\nu}, \\hskip
+            .5cm
+            R_{\\lambda} = \\frac{U_{\\textrm{int}} \\lambda}{\\nu}
+
+        .. [Ishihara] T. Ishihara et al,
+                      *Small-scale statistics in high-resolution direct numerical
+                      simulation of turbulence: Reynolds number dependence of
+                      one-point velocity gradient statistics*.
+                      J. Fluid Mech.,
+                      **592**, 335-366, 2007
+        """
+        for key in ['energy', 'enstrophy']:
+            self.statistics[key + '(t)'] = (self.statistics['dk'] *
+                                            np.sum(self.statistics[key + '(t, k)'], axis = 1))
+        self.statistics['Uint(t)'] = np.sqrt(2*self.statistics['energy(t)'] / 3)
+        self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi /
+                                       (2*self.statistics['Uint(t)']**2)) *
+                                      np.nansum(self.statistics['energy(t, k)'] /
+                                                self.statistics['kshell'][None, :], axis = 1))
+        for key in ['energy',
+                    'enstrophy',
+                    'vel_max',
+                    'Uint',
+                    'Lint']:
+            if key + '(t)' in self.statistics.keys():
+                self.statistics[key] = np.average(self.statistics[key + '(t)'], axis = 0)
+        for suffix in ['', '(t)']:
+            self.statistics['diss'    + suffix] = (self.parameters['nu'] *
+                                                   self.statistics['enstrophy' + suffix]*2)
+            self.statistics['etaK'    + suffix] = (self.parameters['nu']**3 /
+                                                   self.statistics['diss' + suffix])**.25
+            self.statistics['tauK'    + suffix] =  (self.parameters['nu'] /
+                                                    self.statistics['diss' + suffix])**.5
+            self.statistics['Re' + suffix] = (self.statistics['Uint' + suffix] *
+                                              self.statistics['Lint' + suffix] /
+                                              self.parameters['nu'])
+            self.statistics['lambda' + suffix] = (15 * self.parameters['nu'] *
+                                                  self.statistics['Uint' + suffix]**2 /
+                                                  self.statistics['diss' + suffix])**.5
+            self.statistics['Rlambda' + suffix] = (self.statistics['Uint' + suffix] *
+                                                   self.statistics['lambda' + suffix] /
+                                                   self.parameters['nu'])
+            self.statistics['kMeta' + suffix] = (self.statistics['kM'] *
+                                                 self.statistics['etaK' + suffix])
+            if self.parameters['dealias_type'] == 1:
+                self.statistics['kMeta' + suffix] *= 0.8
+        self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
+        self.statistics['Taylor_microscale'] = self.statistics['lambda']
+        return None
+    def set_plt_style(
+            self,
+            style = {'dashes' : (None, None)}):
+        self.style.update(style)
+        return None
+    def convert_complex_from_binary(
+            self,
+            field_name = 'vorticity',
+            iteration = 0,
+            file_name = None):
+        """read the Fourier representation of a vector field.
+
+        Read the binary file containing iteration ``iteration`` of the
+        field ``field_name``, and write it in a ``.h5`` file.
+        """
+        data = np.memmap(
+                os.path.join(self.work_dir,
+                             self.simname + '_{0}_i{1:0>5x}'.format('c' + field_name, iteration)),
+                dtype = self.ctype,
+                mode = 'r',
+                shape = (self.parameters['ny'],
+                         self.parameters['nz'],
+                         self.parameters['nx']//2+1,
+                         3))
+        if type(file_name) == type(None):
+            file_name = self.simname + '_{0}_i{1:0>5x}.h5'.format('c' + field_name, iteration)
+            file_name = os.path.join(self.work_dir, file_name)
+        f = h5py.File(file_name, 'a')
+        f[field_name + '/complex/{0}'.format(iteration)] = data
+        f.close()
+        return None
+    def job_parser_arguments(
+            self,
+            parser):
+        parser.add_argument(
+                '--ncpu',
+                type = int,
+                dest = 'ncpu',
+                default = -1)
+        parser.add_argument(
+                '--np', '--nprocesses',
+                metavar = 'NPROCESSES',
+                help = 'number of mpi processes to use',
+                type = int,
+                dest = 'nb_processes',
+                default = 4)
+        parser.add_argument(
+                '--ntpp', '--nthreads-per-process',
+                type = int,
+                dest = 'nb_threads_per_process',
+                metavar = 'NTHREADS_PER_PROCESS',
+                help = 'number of threads to use per MPI process',
+                default = 1)
+        parser.add_argument(
+                '--no-submit',
+                action = 'store_true',
+                dest = 'no_submit')
+        parser.add_argument(
+                '--environment',
+                type = str,
+                dest = 'environment',
+                default = None)
+        parser.add_argument(
+                '--minutes',
+                type = int,
+                dest = 'minutes',
+                default = 5,
+                help = 'If environment supports it, this is the requested wall-clock-limit.')
+        parser.add_argument(
+               '--njobs',
+               type = int, dest = 'njobs',
+               default = 1)
+        return None
+    def simulation_parser_arguments(
+            self,
+            parser):
+        parser.add_argument(
+                '--simname',
+                type = str, dest = 'simname',
+                default = 'test')
+        parser.add_argument(
+               '-n', '--grid-size',
+               type = int,
+               dest = 'n',
+               default = 32,
+               metavar = 'N',
+               help = 'code is run by default in a grid of NxNxN')
+        for coord in ['x', 'y', 'z']:
+            parser.add_argument(
+                   '--L{0}'.format(coord), '--box-length-{0}'.format(coord),
+                   type = float,
+                   dest = 'L{0}'.format(coord),
+                   default = 2.0,
+                   metavar = 'length{0}'.format(coord),
+                   help = 'length of the box in the {0} direction will be `length{0} x pi`'.format(coord))
+        parser.add_argument(
+                '--wd',
+                type = str, dest = 'work_dir',
+                default = './')
+        parser.add_argument(
+                '--precision',
+                choices = ['single', 'double'],
+                type = str,
+                default = 'single')
+        parser.add_argument(
+                '--src-wd',
+                type = str,
+                dest = 'src_work_dir',
+                default = '')
+        parser.add_argument(
+                '--src-simname',
+                type = str,
+                dest = 'src_simname',
+                default = '')
+        parser.add_argument(
+                '--src-iteration',
+                type = int,
+                dest = 'src_iteration',
+                default = 0)
+        parser.add_argument(
+               '--kMeta',
+               type = float,
+               dest = 'kMeta',
+               default = 2.0)
+        parser.add_argument(
+               '--dtfactor',
+               type = float,
+               dest = 'dtfactor',
+               default = 0.5,
+               help = 'dt is computed as DTFACTOR / N')
+        return None
+    def particle_parser_arguments(
+            self,
+            parser):
+        parser.add_argument(
+               '--particle-rand-seed',
+               type = int,
+               dest = 'particle_rand_seed',
+               default = None)
+        parser.add_argument(
+               '--pclouds',
+               type = int,
+               dest = 'pclouds',
+               default = 1,
+               help = ('number of particle clouds. Particle "clouds" '
+                       'consist of particles distributed according to '
+                       'pcloud-type.'))
+        parser.add_argument(
+                '--pcloud-type',
+                choices = ['random-cube',
+                           'regular-cube'],
+                dest = 'pcloud_type',
+                default = 'random-cube')
+        parser.add_argument(
+               '--particle-cloud-size',
+               type = float,
+               dest = 'particle_cloud_size',
+               default = 2*np.pi)
+        return None
+    def add_parser_arguments(
+            self,
+            parser):
+        subparsers = parser.add_subparsers(
+                dest = 'DNS_class',
+                help = 'type of simulation to run')
+        subparsers.required = True
+        parser_native_binary_to_hdf5 = subparsers.add_parser(
+                'native_binary_to_hdf5',
+                help = 'convert native binary to hdf5')
+        self.simulation_parser_arguments(parser_native_binary_to_hdf5)
+        self.job_parser_arguments(parser_native_binary_to_hdf5)
+        self.parameters_to_parser_arguments(parser_native_binary_to_hdf5)
+        return None
+    def prepare_launch(
+            self,
+            args = []):
+        """Set up reasonable parameters.
+
+        With the default Lundgren forcing applied in the band [2, 4],
+        we can estimate the dissipation, therefore we can estimate
+        :math:`k_M \\eta_K` and constrain the viscosity.
+
+        In brief, the command line parameter :math:`k_M \\eta_K` is
+        used in the following formula for :math:`\\nu` (:math:`N` is the
+        number of real space grid points per coordinate):
+
+        .. math::
+
+            \\nu = \\left(\\frac{2 k_M \\eta_K}{N} \\right)^{4/3}
+
+        With this choice, the average dissipation :math:`\\varepsilon`
+        will be close to 0.4, and the integral scale velocity will be
+        close to 0.77, yielding the approximate value for the Taylor
+        microscale and corresponding Reynolds number:
+
+        .. math::
+
+            \\lambda \\approx 4.75\\left(\\frac{2 k_M \\eta_K}{N} \\right)^{4/6}, \\hskip .5in
+            R_\\lambda \\approx 3.7 \\left(\\frac{N}{2 k_M \\eta_K} \\right)^{4/6}
+
+        """
+        opt = _code.prepare_launch(self, args = args)
+        self.set_precision(opt.precision)
+        self.dns_type = opt.DNS_class
+        self.name = self.dns_type + '-' + self.fluid_precision + '-v' + bfps.__version__
+        # merge parameters if needed
+        for k in self.pp_parameters.keys():
+             self.parameters[k] = self.pp_parameters[k]
+        if len(opt.src_work_dir) == 0:
+            opt.src_work_dir = os.path.realpath(opt.work_dir)
+        if type(opt.dkx) == type(None):
+            opt.dkx = 2. / opt.Lx
+        if type(opt.dky) == type(None):
+            opt.dky = 2. / opt.Ly
+        if type(opt.dkx) == type(None):
+            opt.dkz = 2. / opt.Lz
+        if type(opt.nx) == type(None):
+            opt.nx = opt.n
+        if type(opt.ny) == type(None):
+            opt.ny = opt.n
+        if type(opt.nz) == type(None):
+            opt.nz = opt.n
+        self.pars_from_namespace(opt)
+        return opt
+    def launch(
+            self,
+            args = [],
+            **kwargs):
+        opt = self.prepare_launch(args = args)
+        self.launch_jobs(opt = opt, **kwargs)
+        return None
+    def get_checkpoint_0_fname(self):
+        return os.path.join(
+                    self.work_dir,
+                    self.simname + '_checkpoint_0.h5')
+    def generate_tracer_state(
+            self,
+            rseed = None,
+            species = 0):
+        with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
+            dset = data_file[
+                'tracers{0}/state/0'.format(species)]
+            if not type(rseed) == type(None):
+                np.random.seed(rseed)
+            nn = self.parameters['nparticles']
+            cc = int(0)
+            batch_size = int(1e6)
+            while nn > 0:
+                if nn > batch_size:
+                    dset[cc*batch_size:(cc+1)*batch_size] = np.random.random(
+                            (batch_size, 3))*2*np.pi
+                    nn -= batch_size
+                else:
+                    dset[cc*batch_size:cc*batch_size+nn] = np.random.random(
+                            (nn, 3))*2*np.pi
+                    nn = 0
+                cc += 1
+        return None
+    def generate_vector_field(
+            self,
+            rseed = 7547,
+            spectra_slope = 1.,
+            amplitude = 1.,
+            iteration = 0,
+            field_name = 'vorticity',
+            write_to_file = False,
+            # to switch to constant field, use generate_data_3D_uniform
+            # for scalar_generator
+            scalar_generator = tools.generate_data_3D):
+        """generate vector field.
+
+        The generated field is not divergence free, but it has the proper
+        shape.
+
+        :param rseed: seed for random number generator
+        :param spectra_slope: spectrum of field will look like k^(-p)
+        :param amplitude: all amplitudes are multiplied with this value
+        :param iteration: the field is written at this iteration
+        :param field_name: the name of the field being generated
+        :param write_to_file: should we write the field to file?
+        :param scalar_generator: which function to use for generating the
+            individual components.
+            Possible values: bfps.tools.generate_data_3D,
+            bfps.tools.generate_data_3D_uniform
+        :type rseed: int
+        :type spectra_slope: float
+        :type amplitude: float
+        :type iteration: int
+        :type field_name: str
+        :type write_to_file: bool
+        :type scalar_generator: function
+
+        :returns: ``Kdata``, a complex valued 4D ``numpy.array`` that uses the
+            transposed FFTW layout.
+            Kdata[ky, kz, kx, i] is the amplitude of mode (kx, ky, kz) for
+            the i-th component of the field.
+            (i.e. x is the fastest index and z the slowest index in the
+            real-space representation).
+        """
+        np.random.seed(rseed)
+        Kdata00 = scalar_generator(
+                self.parameters['nz']//2,
+                self.parameters['ny']//2,
+                self.parameters['nx']//2,
+                p = spectra_slope,
+                amplitude = amplitude).astype(self.ctype)
+        Kdata01 = scalar_generator(
+                self.parameters['nz']//2,
+                self.parameters['ny']//2,
+                self.parameters['nx']//2,
+                p = spectra_slope,
+                amplitude = amplitude).astype(self.ctype)
+        Kdata02 = scalar_generator(
+                self.parameters['nz']//2,
+                self.parameters['ny']//2,
+                self.parameters['nx']//2,
+                p = spectra_slope,
+                amplitude = amplitude).astype(self.ctype)
+        Kdata0 = np.zeros(
+                Kdata00.shape + (3,),
+                Kdata00.dtype)
+        Kdata0[..., 0] = Kdata00
+        Kdata0[..., 1] = Kdata01
+        Kdata0[..., 2] = Kdata02
+        Kdata1 = tools.padd_with_zeros(
+                Kdata0,
+                self.parameters['nz'],
+                self.parameters['ny'],
+                self.parameters['nx'])
+        if write_to_file:
+            Kdata1.tofile(
+                    os.path.join(self.work_dir,
+                                 self.simname + "_c{0}_i{1:0>5x}".format(field_name, iteration)))
+        return Kdata1
+    def copy_complex_field(
+            self,
+            src_file_name,
+            src_dset_name,
+            dst_file,
+            dst_dset_name,
+            make_link = True):
+        # I define a min_shape thingie, but for now I only trust this method for
+        # the case of increasing/decreasing by the same factor in all directions.
+        # in principle we could write something more generic, but i'm not sure
+        # how complicated that would be
+        dst_shape = (self.parameters['nz'],
+                     self.parameters['ny'],
+                     (self.parameters['nx']+2) // 2,
+                     3)
+        src_file = h5py.File(src_file_name, 'r')
+        if (src_file[src_dset_name].shape == dst_shape):
+            if make_link and (src_file[src_dset_name].dtype == self.ctype):
+                dst_file[dst_dset_name] = h5py.ExternalLink(
+                        src_file_name,
+                        src_dset_name)
+            else:
+                dst_file.create_dataset(
+                        dst_dset_name,
+                        shape = dst_shape,
+                        dtype = self.ctype,
+                        fillvalue = 0.0)
+                for kz in range(src_file[src_dset_name].shape[0]):
+                    dst_file[dst_dset_name][kz] = src_file[src_dset_name][kz]
+        else:
+            print('aloha')
+            min_shape = (min(dst_shape[0], src_file[src_dset_name].shape[0]),
+                         min(dst_shape[1], src_file[src_dset_name].shape[1]),
+                         min(dst_shape[2], src_file[src_dset_name].shape[2]),
+                         3)
+            print(self.ctype)
+            dst_file.create_dataset(
+                    dst_dset_name,
+                    shape = dst_shape,
+                    dtype = np.dtype(self.ctype),
+                    fillvalue = complex(0))
+            for kz in range(min_shape[0]):
+                dst_file[dst_dset_name][kz,:min_shape[1], :min_shape[2]] = \
+                        src_file[src_dset_name][kz, :min_shape[1], :min_shape[2]]
+        return None
+    def launch_jobs(
+            self,
+            opt = None,
+            particle_initial_condition = None):
+        self.rewrite_par(
+                group = self.dns_type,
+                parameters = self.pp_parameters)
+        self.run(
+                nb_processes = opt.nb_processes,
+                nb_threads_per_process = opt.nb_threads_per_process,
+                njobs = opt.njobs,
+                hours = opt.minutes // 60,
+                minutes = opt.minutes % 60,
+                no_submit = opt.no_submit)
+        return None
+
diff --git a/bfps/__main__.py b/bfps/__main__.py
index f1ee7278da0ea555bab5307d4e5d7785352030fa..fa30289b430c7db9508c2e9fae796f6fe47b8eb6 100644
--- a/bfps/__main__.py
+++ b/bfps/__main__.py
@@ -29,6 +29,7 @@ import argparse
 
 import bfps
 from .DNS import DNS
+from .PP import PP
 from .NavierStokes import NavierStokes
 from .NSVorticityEquation import NSVorticityEquation
 from .FluidResize import FluidResize
@@ -65,7 +66,7 @@ def main():
                'NSManyParticles-double']
     parser.add_argument(
             'base_class',
-            choices = ['DNS'] +
+            choices = ['DNS', 'PP'] +
                       NSoptions +
                       NSVEoptions +
                       FRoptions +
@@ -81,6 +82,10 @@ def main():
         c = DNS()
         c.launch(args = sys.argv[2:])
         return None
+    if opt.base_class == 'PP':
+        c = PP()
+        c.launch(args = sys.argv[2:])
+        return None
     if 'double' in opt.base_class:
         precision = 'double'
     else:
diff --git a/bfps/_code.py b/bfps/_code.py
index 37a996f855a82db8ce4154a5ce10a18ae9b55dfb..22bcd9101ff6591e00f0455c1de1af2698c5f842 100644
--- a/bfps/_code.py
+++ b/bfps/_code.py
@@ -249,6 +249,8 @@ class _code(_base):
             assert(self.compile_code() == 0)
             if self.work_dir != os.path.realpath(os.getcwd()):
                 shutil.copy(self.name, self.work_dir)
+        if 'niter_todo' not in self.parameters.keys():
+            self.parameters['niter_todo'] = 1
         current_dir = os.getcwd()
         os.chdir(self.work_dir)
         os.chdir(current_dir)
diff --git a/bfps/cpp/full_code/native_binary_to_hdf5.cpp b/bfps/cpp/full_code/native_binary_to_hdf5.cpp
index 40d64e99e2bccf9dff77025e9c107f35fea7e6f6..78a29addbb18dcb0f162ce4e95d7cc7733a5684d 100644
--- a/bfps/cpp/full_code/native_binary_to_hdf5.cpp
+++ b/bfps/cpp/full_code/native_binary_to_hdf5.cpp
@@ -50,6 +50,20 @@ int native_binary_to_hdf5<rnumber>::finalize(void)
     return EXIT_SUCCESS;
 }
 
+template <typename rnumber>
+int native_binary_to_hdf5<rnumber>::read_parameters(void)
+{
+    this->postprocess::read_parameters();
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string(".h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    this->iteration_list = hdf5_tools::read_vector<int>(
+            parameter_file,
+            "/native_binary_to_hdf5/iteration_list");
+    return EXIT_SUCCESS;
+}
+
 template class native_binary_to_hdf5<float>;
 template class native_binary_to_hdf5<double>;
 
diff --git a/bfps/cpp/full_code/postprocess.cpp b/bfps/cpp/full_code/postprocess.cpp
index 3c6eed59921afa17c3b7c4fa048df41f42eef350..193dafeb2ab4b6a892be11bcf913f766ecaaa531 100644
--- a/bfps/cpp/full_code/postprocess.cpp
+++ b/bfps/cpp/full_code/postprocess.cpp
@@ -26,3 +26,66 @@ int postprocess::main_loop(void)
     return EXIT_SUCCESS;
 }
 
+
+int postprocess::read_parameters()
+{
+    hid_t parameter_file;
+    hid_t dset, memtype, space;
+    char fname[256];
+    hsize_t dims[1];
+    char *string_data;
+    sprintf(fname, "%s.h5", this->simname.c_str());
+    parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+    dset = H5Dopen(parameter_file, "/parameters/dealias_type", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dealias_type);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/dkx", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkx);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/dky", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dky);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/dkz", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkz);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/dt", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dt);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/famplitude", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->famplitude);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/fk0", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->fk0);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/fk1", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->fk1);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/fmode", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->fmode);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/forcing_type", H5P_DEFAULT);
+    space = H5Dget_space(dset);
+    memtype = H5Dget_type(dset);
+    string_data = (char*)malloc(256);
+    H5Dread(dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &string_data);
+    sprintf(this->forcing_type, "%s", string_data);
+    free(string_data);
+    H5Sclose(space);
+    H5Tclose(memtype);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/nu", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nu);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/nx", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nx);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/ny", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->ny);
+    H5Dclose(dset);
+    dset = H5Dopen(parameter_file, "/parameters/nz", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nz);
+    H5Dclose(dset);
+    H5Fclose(parameter_file);
+    return 0;
+}
+
diff --git a/bfps/cpp/full_code/postprocess.hpp b/bfps/cpp/full_code/postprocess.hpp
index 5d3e578d7dad0975fe47848e1f6224442a219b80..c80fc3f2dfdc35691d9e69442fa3ad7b6e592891 100644
--- a/bfps/cpp/full_code/postprocess.hpp
+++ b/bfps/cpp/full_code/postprocess.hpp
@@ -40,6 +40,15 @@ class postprocess: public code_base
         std::vector<int> iteration_list;
         hid_t stat_file;
 
+        /* parameters that are read in read_parameters */
+        double dt;
+        double famplitude;
+        double fk0;
+        double fk1;
+        int fmode;
+        char forcing_type[512];
+        double nu;
+
         postprocess(
                 const MPI_Comm COMMUNICATOR,
                 const std::string &simulation_name):
@@ -53,6 +62,7 @@ class postprocess: public code_base
         virtual int finalize(void) = 0;
 
         int main_loop(void);
+        virtual int read_parameters(void);
 };
 
 #endif//POSTPROCESS_HPP
diff --git a/bfps/cpp/hdf5_tools.cpp b/bfps/cpp/hdf5_tools.cpp
index 6158607dc829db4af434d203b1cf7639b9b5f3c7..abd894aa15c5cc7a5f7ad3f983bf7a35182d87e3 100644
--- a/bfps/cpp/hdf5_tools.cpp
+++ b/bfps/cpp/hdf5_tools.cpp
@@ -111,32 +111,29 @@ std::vector<dtype> hdf5_tools::read_vector_with_single_rank(
     return data;
 }
 
-namespace hdf5_tools
-{
-template <>
-std::vector<int> read_vector(
+template
+std::vector<int> hdf5_tools::read_vector<int>(
         const hid_t,
         const std::string);
 
-template <>
-std::vector<double> read_vector(
+template
+std::vector<double> hdf5_tools::read_vector<double>(
         const hid_t,
         const std::string);
 
-template <>
-std::vector<int> read_vector_with_single_rank(
+template
+std::vector<int> hdf5_tools::read_vector_with_single_rank<int>(
         const int myrank,
         const int rank_to_use,
         const MPI_Comm COMM,
         const hid_t file_id,
         const std::string dset_name);
 
-template <>
-std::vector<double> read_vector_with_single_rank(
+template
+std::vector<double> hdf5_tools::read_vector_with_single_rank<double>(
         const int myrank,
         const int rank_to_use,
         const MPI_Comm COMM,
         const hid_t file_id,
         const std::string dset_name);
-}