diff --git a/CREDITS b/CREDITS
index 448267ee..1489656a 100644
--- a/CREDITS
+++ b/CREDITS
@@ -10,6 +10,7 @@ The development team consists of (in alphabetic order):
   Sylwester Arabas 
     (core code, library design, OOP-formulae concepts [4], 
     concurrency, gnuplot output)
+  Piotr Dziekan
   Dorota Jarecka 
     (shallow-water systems, OOP-formulae concepts [4])
   Anna Jaruga 
diff --git a/libmpdata++/output/detail/output_common.hpp b/libmpdata++/output/detail/output_common.hpp
index 47920259..b3258c18 100644
--- a/libmpdata++/output/detail/output_common.hpp
+++ b/libmpdata++/output/detail/output_common.hpp
@@ -29,7 +29,7 @@ namespace libmpdataxx
 
         int do_record_cnt = 0;
         typename parent_t::real_t record_time;
-        const typename parent_t::advance_arg_t outfreq;
+        const typename parent_t::advance_arg_t outfreq, outstart;
         const int outwindow;
         const std::string outdir;
 
@@ -158,7 +158,7 @@ namespace libmpdataxx
               record_time = this->time;
               for (int t = 0; t < outwindow; ++t)
               {
-                if ((this->timestep - t) % static_cast<int>(outfreq) == 0) record_all();
+                if ((this->timestep - t) % static_cast<int>(outfreq) == 0 && this->timestep >= static_cast<int>(outstart)) record_all();
               }
             }
           }
@@ -170,7 +170,8 @@ namespace libmpdataxx
 
         struct rt_params_t : parent_t::rt_params_t
         {
-          typename parent_t::advance_arg_t outfreq = 1;
+          typename parent_t::advance_arg_t outfreq = 1,
+                                           outstart = 0; // TODO: make it work for var_dt
           int outwindow = 1;
           std::map<int, info_t> outvars;
           std::string outdir;
@@ -184,11 +185,15 @@ namespace libmpdataxx
         ) :
           parent_t(args, p),
           outfreq(p.outfreq),
+          outstart(p.outstart),
           outwindow(p.outwindow),
           outvars(p.outvars),
           outdir(p.outdir),
           intrp_vars(args.mem->tmp[__FILE__][0])
         {
+          assert(int(outstart) % int(outfreq) == 0);
+          assert(int(outstart)==0 || this->var_dt==0);
+
           // default value for outvars
           if (this->outvars.size() == 0 && parent_t::n_eqns == 1)
             outvars = {{0, {"", ""}}};
diff --git a/libmpdata++/output/gnuplot.hpp b/libmpdata++/output/gnuplot.hpp
index 9ffd8f4d..432e1163 100644
--- a/libmpdata++/output/gnuplot.hpp
+++ b/libmpdata++/output/gnuplot.hpp
@@ -100,6 +100,7 @@ namespace libmpdataxx
 
           for (int t = 0; t <= nt; t+=p.outfreq)
           {
+            if(t > 0 && t < p.outstart) continue;
             for (const auto &v : this->outvars)
             {
               *gp << ", '-'";
diff --git a/libmpdata++/output/hdf5.hpp b/libmpdata++/output/hdf5.hpp
index 7e890916..a9e07ae2 100644
--- a/libmpdata++/output/hdf5.hpp
+++ b/libmpdata++/output/hdf5.hpp
@@ -198,12 +198,21 @@ namespace libmpdataxx
 
             // T
             {
-              const hsize_t
-                nt_out = nt / this->outfreq + 1; // incl. t=0
+              // incl. t=0 and t=outstart
+              const hsize_t nt_out = this->outstart == 0 ? 
+                (nt - this->outstart) / this->outfreq + 1 :
+                (nt - this->outstart) / this->outfreq + 2;  // for outstart>0 we still want to store t=0
+
               float dt = this->dt;
 
               blitz::Array<typename solver_t::real_t, 1> coord(nt_out);
-              coord = (this->var_dt ? this->outfreq : this->outfreq * this->dt) * blitz::firstIndex();
+              if(this->outstart == 0)
+                coord(blitz::Range(0,nt_out-1)) = (this->var_dt ? this->outfreq : this->outfreq * this->dt) * (blitz::firstIndex() + this->outstart/this->outfreq);
+              else
+              {
+                coord(blitz::Range(1,nt_out-1)) = (this->var_dt ? this->outfreq : this->outfreq * this->dt) * (blitz::firstIndex() + this->outstart/this->outfreq);
+                coord(0)=0;
+              }
 
               auto curr_dim = (*hdfp).createDataSet("T", flttype_output, H5::DataSpace(1, &nt_out));
 
@@ -365,6 +374,22 @@ namespace libmpdataxx
         aux.write(data, flttype_solver, H5::DataSpace(parent_t::n_dims, shape.data()), space, dxpl_id);
       }
 
+      // for 1-D arrays
+      void record_aux_hlpr(const std::string &name, typename solver_t::real_t *data, hsize_t size, H5::H5File hdf)
+      {
+        assert(this->rank == 0);
+
+        auto aux = hdf.createDataSet(
+          name,
+          flttype_output,
+          H5::DataSpace(1, &size)
+        );
+
+        auto space = aux.getSpace();
+        space.selectHyperslab(H5S_SELECT_SET, &size, &zero);
+        aux.write(data, flttype_solver, H5::DataSpace(1, &size), space, dxpl_id);
+      }
+
       // for discontiguous array with halos
       void record_aux_dsc_hlpr(const std::string &name, const typename solver_t::arr_t &arr, H5::H5File hdf, bool srfc = false)
       {
@@ -516,10 +541,24 @@ namespace libmpdataxx
       // has to be called after const file was created (i.e. after start())
       void record_aux_const(const std::string &name, typename solver_t::real_t *data)
       {
-        H5::H5File hdfcp(const_file, H5F_ACC_RDWR); // reopen the const file
+        H5::H5File hdfcp(const_file, H5F_ACC_RDWR
+#if defined(USE_MPI)
+          , H5P_DEFAULT, fapl_id
+#endif
+        ); // reopen the const file
         record_aux_hlpr(name, data, hdfcp);
       }
 
+      void record_aux_const(const std::string &name, typename solver_t::real_t *data, const int &size)
+      {
+        H5::H5File hdfcp(const_file, H5F_ACC_RDWR
+#if defined(USE_MPI)
+          , H5P_DEFAULT, fapl_id
+#endif
+        ); // reopen the const file
+        record_aux_hlpr(name, data, size, hdfcp);
+      }
+
       void record_aux_const(const std::string &name, const std::string &group_name, typename solver_t::real_t data)
       {
         H5::H5File hdfcp(const_file, H5F_ACC_RDWR