Repository: OceanParcels/parcels
Pull request #497 (2018-11-27 07:28:37)

Faster output writing by merging temporary files into one large netcdf file at end of run
This PR builds upon the code developed in #486 by @pjpetersik and Frederik Felius. It is a complete redesign of the ParticleFile writing routine in Parcels, and thus addresses #473. Instead of writing all particle data to a NetCDF file with two unlimited dimensions directly, this PR makes writing a two-step process:

1) During execution, every `ParticleFile.write()` creates a temporary `.npy` file, stored in a temporary directory
2) When the ParticleFile is closed (when the script terminates, or by issuing `ParticleFile.close()` or `ParticleFile.export()`), all temporary `.npy` files are combined in one large NetCDF file. Since dimensions are known by that time, the file creation is much faster and more efficient

Things still to do:

- [x] Confirm that `to_write=once` variables are aways written to correct Particle ID
- [x] Improve doc-strings of new functions/methods
- [x] Test in realistic runs, and confirm speed-up



Patch:
diff --git a/parcels/examples/parcels_tutorial.ipynb b/parcels/examples/parcels_tutorial.ipynb
index 443dcc8f..0407705e 100644
--- a/parcels/examples/parcels_tutorial.ipynb
+++ b/parcels/examples/parcels_tutorial.ipynb
@@ -195,10 +195,11 @@
     }
    ],
    "source": [
+    "output_file = pset.ParticleFile(name=\"EddyParticles.nc\", outputdt=timedelta(hours=1)) # the file name and the time step of the outputs\n",
     "pset.execute(AdvectionRK4,                 # the kernel (which defines how particles move)\n",
     "             runtime=timedelta(days=6),    # the total length of the run\n",
     "             dt=timedelta(minutes=5),      # the timestep of the kernel\n",
-    "             output_file=pset.ParticleFile(name=\"EddyParticles.nc\", outputdt=timedelta(hours=1)))  # the file name and the time step of the outputs"
+    "             output_file=output_file)"
    ]
   },
   {
@@ -250,7 +251,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The trajectory information of the particles is stored in the `EddyParticles.nc` file. It can be quickly plotted using the `plotTrajectoriesFile` function."
+    "The trajectory information of the particles can be written to the `EddyParticles.nc` file by using the `.export()` method on the output file. The trajectory can then be quickly plotted using the `plotTrajectoriesFile` function."
    ]
   },
   {
@@ -272,6 +273,7 @@
     }
    ],
    "source": [
+    "output_file.export()\n",
     "plotTrajectoriesFile('EddyParticles.nc');"
    ]
   },
@@ -1224,10 +1226,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "output_file = pset.ParticleFile(name=\"EddyParticles_Bwd.nc\", outputdt=timedelta(hours=1))  # the file name and the time step of the outputs\n",
     "pset.execute(AdvectionRK4,\n",
     "             dt=-timedelta(minutes=5),      # negative timestep for backward run\n",
     "             runtime=timedelta(days=6),     # the run time\n",
-    "             output_file=pset.ParticleFile(name=\"EddyParticles_Bwd.nc\", outputdt=timedelta(hours=1)))  # the file name and the time step of the outputs"
+    "             output_file=output_file)"
    ]
   },
   {
@@ -1326,10 +1329,11 @@
     "\n",
     "k_WestVel = pset.Kernel(WestVel)        # casting the WestVel function to a kernel object\n",
     "\n",
+    "output_file = pset.ParticleFile(name=\"EddyParticles_WestVel.nc\", outputdt=timedelta(hours=1))\n",
     "pset.execute(AdvectionRK4 + k_WestVel,  # simply add kernels using the + operator\n",
     "             runtime=timedelta(days=2),\n",
     "             dt=timedelta(minutes=5),\n",
-    "             output_file=pset.ParticleFile(name=\"EddyParticles_WestVel.nc\", outputdt=timedelta(hours=1)))"
+    "             output_file=output_file)"
    ]
   },
   {
@@ -1358,6 +1362,7 @@
     }
    ],
    "source": [
+    "output_file.export()\n",
     "plotTrajectoriesFile('EddyParticles_WestVel.nc');"
    ]
   },
@@ -1487,10 +1492,11 @@
     }
    ],
    "source": [
+    "output_file = pset.ParticleFile(name=\"GlobCurrentParticles.nc\", outputdt=timedelta(hours=6))\n",
     "pset.execute(AdvectionRK4,\n",
     "             runtime=timedelta(days=10),\n",
     "             dt=timedelta(minutes=5),\n",
-    "             output_file=pset.ParticleFile(name=\"GlobCurrentParticles.nc\", outputdt=timedelta(hours=6)))"
+    "             output_file=output_file)"
    ]
   },
   {
@@ -1526,6 +1532,7 @@
     }
    ],
    "source": [
+    "output_file.export()\n",
     "plotTrajectoriesFile('GlobCurrentParticles.nc',\n",
     "                     tracerfile='GlobCurrent_example_data/20020101000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc',\n",
     "                     tracerlon='lon',\n",
diff --git a/parcels/examples/tutorial_Argofloats.ipynb b/parcels/examples/tutorial_Argofloats.ipynb
index cb2dec69..520b53d1 100644
--- a/parcels/examples/tutorial_Argofloats.ipynb
+++ b/parcels/examples/tutorial_Argofloats.ipynb
@@ -123,9 +123,11 @@
     "# combine Argo vertical movement kernel with built-in Advection kernel\n",
     "kernels = ArgoVerticalMovement + pset.Kernel(AdvectionRK4)\n",
     "\n",
+    "# Create a ParticleFile object to store the output\n",
+    "output_file = pset.ParticleFile(name=\"argo_float\", outputdt=timedelta(minutes=30))\n",
+    "\n",
     "# Now execute the kernels for 30 days, saving data every 30 minutes\n",
-    "pset.execute(kernels, runtime=timedelta(days=30), dt=timedelta(minutes=5), \n",
-    "             output_file=pset.ParticleFile(name=\"argo_float\", outputdt=timedelta(minutes=30)))"
+    "pset.execute(kernels, runtime=timedelta(days=30), dt=timedelta(minutes=5), output_file=output_file)"
    ]
   },
   {
@@ -157,6 +159,8 @@
     "from mpl_toolkits.mplot3d import Axes3D\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
+    "\n",
     "nc = netCDF4.Dataset(\"argo_float.nc\")\n",
     "x = nc.variables[\"lon\"][:].squeeze()\n",
     "y = nc.variables[\"lat\"][:].squeeze()\n",
diff --git a/parcels/examples/tutorial_NestedFields.ipynb b/parcels/examples/tutorial_NestedFields.ipynb
index d2f9fd91..abf972c6 100644
--- a/parcels/examples/tutorial_NestedFields.ipynb
+++ b/parcels/examples/tutorial_NestedFields.ipynb
@@ -135,6 +135,8 @@
     "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[1000])\n",
     "output_file = pset.ParticleFile(name='NestedFieldParticle.nc', outputdt=50)\n",
     "pset.execute(AdvectionRK4, runtime=14000, dt=10, output_file=output_file)\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
+    "\n",
     "plt = plotTrajectoriesFile('NestedFieldParticle.nc', show_plt=False)\n",
     "plt.plot([0,2e3,2e3,0,0],[0,0,2e3,2e3,0], c='orange')\n",
     "plt.plot([-2e3,18e3,18e3,-2e3,-2e3],[-1e3,-1e3,3e3,3e3,-1e3], c='green');"
diff --git a/parcels/examples/tutorial_SummedFields.ipynb b/parcels/examples/tutorial_SummedFields.ipynb
index 7461f884..c02c5538 100644
--- a/parcels/examples/tutorial_SummedFields.ipynb
+++ b/parcels/examples/tutorial_SummedFields.ipynb
@@ -84,6 +84,7 @@
     "pset = ParticleSet(fieldset_flow, pclass=JITParticle, lon=[0], lat=[900])\n",
     "output_file = pset.ParticleFile(name='SummedFieldParticle_flow.nc', outputdt=1)\n",
     "pset.execute(AdvectionRK4, runtime=10, dt=1, output_file=output_file)\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('SummedFieldParticle_flow.nc');"
    ]
   },
@@ -149,6 +150,7 @@
     "pset = ParticleSet(fieldset_stokes, pclass=JITParticle, lon=[0], lat=[900])\n",
     "output_file = pset.ParticleFile(name='SummedFieldParticle_stokes.nc', outputdt=1)\n",
     "pset.execute(AdvectionRK4, runtime=10, dt=1, output_file=output_file)\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('SummedFieldParticle_stokes.nc');"
    ]
   },
@@ -202,6 +204,7 @@
     "pset = ParticleSet(fieldset_sum, pclass=JITParticle, lon=[0], lat=[900])\n",
     "output_file = pset.ParticleFile(name='SummedFieldParticle_sum.nc', outputdt=1)\n",
     "pset.execute(AdvectionRK4, runtime=10, dt=1, output_file=output_file)\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('SummedFieldParticle_sum.nc');"
    ]
   },
diff --git a/parcels/examples/tutorial_delaystart.ipynb b/parcels/examples/tutorial_delaystart.ipynb
index 6a720548..b3e31f09 100644
--- a/parcels/examples/tutorial_delaystart.ipynb
+++ b/parcels/examples/tutorial_delaystart.ipynb
@@ -102,9 +102,10 @@
     }
    ],
    "source": [
-    "output_file = pset.ParticleFile(name=\"DelayParticle_time\", outputdt=delta(hours=1))\n",
+    "output_file = pset.ParticleFile(name=\"DelayParticle_time.nc\", outputdt=delta(hours=1))\n",
     "pset.execute(AdvectionRK4, runtime=delta(hours=24), dt=delta(minutes=5),\n",
-    "             output_file=output_file)"
+    "             output_file=output_file)\n",
+    "output_file.export()  # export the trajectory data to a netcdf file"
    ]
   },
   {
@@ -1162,6 +1163,7 @@
     }
    ],
    "source": [
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('DelayParticle_releasedt.nc', mode='movie2d_notebook')"
    ]
   },
@@ -1718,6 +1720,7 @@
     "pset.execute(AdvectionRK4, runtime=delta(hours=15), dt=delta(minutes=5),\n",
     "             output_file=output_file)\n",
     "\n",
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('DelayParticle_releasedt_9hrs.nc', mode='movie2d_notebook')"
    ]
   },
diff --git a/parcels/examples/tutorial_nemo_curvilinear.ipynb b/parcels/examples/tutorial_nemo_curvilinear.ipynb
index 7a07100c..942a2ffa 100644
--- a/parcels/examples/tutorial_nemo_curvilinear.ipynb
+++ b/parcels/examples/tutorial_nemo_curvilinear.ipynb
@@ -167,6 +167,7 @@
     }
    ],
    "source": [
+    "pfile.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile(\"nemo_particles.nc\");"
    ]
   }
diff --git a/parcels/examples/tutorial_periodic_boundaries.ipynb b/parcels/examples/tutorial_periodic_boundaries.ipynb
index cbc0b280..98c96ed6 100644
--- a/parcels/examples/tutorial_periodic_boundaries.ipynb
+++ b/parcels/examples/tutorial_periodic_boundaries.ipynb
@@ -156,6 +156,7 @@
     }
    ],
    "source": [
+    "output_file.export()  # export the trajectory data to a netcdf file\n",
     "plotTrajectoriesFile('PeriodicParticle.nc');"
    ]
   },
diff --git a/parcels/particle.py b/parcels/particle.py
index 30edeba3..9e4b8cd7 100644
--- a/parcels/particle.py
+++ b/parcels/particle.py
@@ -165,7 +165,6 @@ class ScipyParticle(_Particle):
     depth = Variable('depth', dtype=np.float32)
     time = Variable('time', dtype=np.float64)
     id = Variable('id', dtype=np.int32)
-    fileid = Variable('fileid', dtype=np.int32, to_write=False)
     dt = Variable('dt', dtype=np.float32, to_write=False)
     state = Variable('state', dtype=np.int32, initial=ErrorCode.Success, to_write=False)
 
@@ -179,7 +178,6 @@ class ScipyParticle(_Particle):
         type(self).time.initial = time
         type(self).id.initial = lastID
         lastID += 1
-        type(self).fileid.initial = -1  # -1 means particle is not written yet
         type(self).dt.initial = None
         super(ScipyParticle, self).__init__()
 
diff --git a/parcels/particlefile.py b/parcels/particlefile.py
index fe48eb87..4e6e5f90 100644
--- a/parcels/particlefile.py
+++ b/parcels/particlefile.py
@@ -3,12 +3,21 @@ import numpy as np
 import netCDF4
 from datetime import timedelta as delta
 from parcels.tools.loggers import logger
+import os
+import shutil
+import string
+import random
 from parcels.tools.error import ErrorCode
-from os import path
 try:
     from parcels._version import version as parcels_version
 except:
     raise EnvironmentError('Parcels version can not be retrieved. Have you run ''python setup.py install''?')
+try:
+    from os import getuid
+except:
+    # Windows does not have getuid(), so define to simply return 'tmp'
+    def getuid():
+        return 'tmp'
 
 
 __all__ = ['ParticleFile']
@@ -31,122 +40,148 @@ def _set_calendar(origin_calendar):
 
 
 class ParticleFile(object):
-    """Initialise netCDF4.Dataset for trajectory output.
-
-    The output follows the format outlined in the Discrete Sampling Geometries
-    section of the CF-conventions:
-    http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#discrete-sampling-geometries
-
-    The current implementation is based on the NCEI template:
-    http://www.nodc.noaa.gov/data/formats/netcdf/v2.0/trajectoryIncomplete.cdl
-
-    Developer note: We cannot use xray.Dataset here, since it does not yet allow
-    incremental writes to disk: https://github.com/pydata/xarray/issues/199
-
+    """Initialise trajectory output.
     :param name: Basename of the output file
     :param particleset: ParticleSet to output
     :param outputdt: Interval which dictates the update frequency of file output
                      while ParticleFile is given as an argument of ParticleSet.execute()
                      It is either a timedelta object or a positive double.
     :param write_ondelete: Boolean to write particle data only when they are deleted. Default is False
-    :param chunksizes: 2d vector of sizes for NetCDF chunking. Lower values means smaller files, but also slower IO.
-                       See e.g. https://www.unidata.ucar.edu/blogs/developer/entry/chunking_data_choosing_shapes
+    :param tempwritedir: directory to write temporary files to during executing.
+                        Default is out-XXXXXX where Xs are random capitals
+    :param pset_info: dictionary of info on the ParticleSet, stored in tempwritedir/pset_info.npy,
+                      used to create NetCDF file from npy-files.
     """
 
-    def __init__(self, name, particleset, outputdt=np.infty, write_ondelete=False, chunksizes=None):
+    def __init__(self, name, particleset, outputdt=np.infty, write_ondelete=False, tempwritedir=None, pset_info=None):
 
-        self.name = name
-        self.chunksizes = [max([len(particleset), 1]), 1] if chunksizes is None else chunksizes
         self.write_ondelete = write_ondelete
         self.outputdt = outputdt
-        self.lasttraj = 0  # id of last particle written
         self.lasttime_written = None  # variable to check if time has been written already
 
         self.dataset = None
-        self.particleset = particleset
-
-    def open_dataset(self):
-        extension = path.splitext(str(self.name))[1]
+        self.metadata = {}
+        if pset_info is not None:
+            for v in pset_info.keys():
+                setattr(self, v, pset_info[v])
+        else:
+            self.name = name
+            self.particleset = particleset
+            self.parcels_mesh = self.particleset.fieldset.gridset.grids[0].mesh
+            self.time_origin = self.particleset.time_origin
+            self.lonlatdepth_dtype = self.particleset.lonlatdepth_dtype
+            self.var_names = []
+            self.var_names_once = []
+            for v in self.particleset.ptype.variables:
+                if v.to_write == 'once':
+                    self.var_names_once += [v.name]
+                elif v.to_write is True:
+                    self.var_names += [v.name]
+            if len(self.var_names_once) > 0:
+                self.written_once = []
+                self.file_list_once = []
+
+            self.file_list = []
+            self.time_written = []
+            self.maxid_written = -1
+        self.to_export = False
+
+        if tempwritedir is None:
+            self.tempwritedir = os.path.join(os.path.dirname(str(self.name)), "out-%s" % ''.join(random.choice(string.ascii_uppercase) for _ in range(8)))
+        else:
+            self.tempwritedir = tempwritedir
+
+        self.delete_tempwritedir()
+
+    def open_netcdf_file(self, data_shape):
+        """Initialise NetCDF4.Dataset for trajectory output.
+        The output follows the format outlined in the Discrete Sampling Geometries
+        section of the CF-conventions:
+        http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#discrete-sampling-geometries
+        The current implementation is based on the NCEI template:
+        http://www.nodc.noaa.gov/data/formats/netcdf/v2.0/trajectoryIncomplete.cdl
+        :param data_shape: shape of the variables in the NetCDF4 file
+        """
+        extension = os.path.splitext(str(self.name))[1]
         fname = self.name if extension in ['.nc', '.nc4'] else "%s.nc" % self.name
+        if os.path.exists(str(fname)):
+            os.remove(str(fname))
         self.dataset = netCDF4.Dataset(fname, "w", format="NETCDF4")
-        self.dataset.createDimension("obs", None)
-        self.dataset.createDimension("traj", None)
+        self.dataset.createDimension("obs", data_shape[1])
+        self.dataset.createDimension("traj", data_shape[0])
         coords = ("traj", "obs")
         self.dataset.feature_type = "trajectory"
         self.dataset.Conventions = "CF-1.6/CF-1.7"
         self.dataset.ncei_template_version = "NCEI_NetCDF_Trajectory_Template_v2.0"
         self.dataset.parcels_version = parcels_version
-        self.dataset.parcels_mesh = self.particleset.fieldset.gridset.grids[0].mesh
+        self.dataset.parcels_mesh = self.parcels_mesh
 
         # Create ID variable according to CF conventions
-        self.id = self.dataset.createVariable("trajectory", "i4", coords, chunksizes=self.chunksizes)
+        self.id = self.dataset.createVariable("trajectory", "i4", coords)
         self.id.long_name = "Unique identifier for each particle"
         self.id.cf_role = "trajectory_id"
 
         # Create time, lat, lon and z variables according to CF conventions:
-        self.time = self.dataset.createVariable("time", "f8", coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.time = self.dataset.createVariable("time", "f8", coords, fill_value=np.nan)
         self.time.long_name = ""
         self.time.standard_name = "time"
-        if self.particleset.time_origin.calendar is None:
+        if self.time_origin.calendar is None:
             self.time.units = "seconds"
         else:
-            self.time.units = "seconds since " + str(self.particleset.time_origin)
-            self.time.calendar = _set_calendar(self.particleset.time_origin.calendar)
+            self.time.units = "seconds since " + str(self.time_origin)
+            self.time.calendar = 'standard' if self.time_origin.calendar == 'np_datetime64' else self.time_origin.calendar
         self.time.axis = "T"
 
-        if self.particleset.lonlatdepth_dtype is np.float64:
+        if self.lonlatdepth_dtype is np.float64:
             lonlatdepth_precision = "f8"
         else:
             lonlatdepth_precision = "f4"
 
-        self.lat = self.dataset.createVariable("lat", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.lat = self.dataset.createVariable("lat", lonlatdepth_precision, coords, fill_value=np.nan)
         self.lat.long_name = ""
         self.lat.standard_name = "latitude"
         self.lat.units = "degrees_north"
         self.lat.axis = "Y"
 
-        self.lon = self.dataset.createVariable("lon", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.lon = self.dataset.createVariable("lon", lonlatdepth_precision, coords, fill_value=np.nan)
         self.lon.long_name = ""
         self.lon.standard_name = "longitude"
         self.lon.units = "degrees_east"
         self.lon.axis = "X"
 
-        self.z = self.dataset.createVariable("z", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.z = self.dataset.createVariable("z", lonlatdepth_precision, coords, fill_value=np.nan)
         self.z.long_name = ""
         self.z.standard_name = "depth"
         self.z.units = "m"
         self.z.positive = "down"
 
-        self.user_vars = []
-        self.user_vars_once = []
-        """
-        :user_vars: list of additional user defined particle variables to write for all particles and all times
-        :user_vars_once: list of additional user defined particle variables to write for all particles only once at initial time.
-        """
+        for vname in self.var_names:
+            if vname not in ['time', 'lat', 'lon', 'depth', 'id']:
+                setattr(self, vname, self.dataset.createVariable(vname, "f4", coords, fill_value=np.nan))
+                getattr(self, vname).long_name = ""
+                getattr(self, vname).standard_name = vname
+                getattr(self, vname).units = "unknown"
+
+        for vname in self.var_names_once:
+            setattr(self, vname, self.dataset.createVariable(vname, "f4", "traj", fill_value=np.nan))
+            getattr(self, vname).long_name = ""
+            getattr(self, vname).standard_name = vname
+            getattr(self, vname).units = "unknown"
 
-        for v in self.particleset.ptype.variables:
-            if v.name in ['time', 'lat', 'lon', 'depth', 'z', 'id']:
-                continue
-            if v.to_write:
-                if v.to_write is True:
-                    setattr(self, v.name, self.dataset.createVariable(v.name, "f4", coords, fill_value=np.nan, chunksizes=self.chunksizes))
-                    self.user_vars += [v.name]
-                elif v.to_write == 'once':
-                    setattr(self, v.name, self.dataset.createVariable(v.name, "f4", "traj", fill_value=np.nan, chunksizes=[self.chunksizes[0]]))
-                    self.user_vars_once += [v.name]
-                getattr(self, v.name).long_name = ""
-                getattr(self, v.name).standard_name = v.name
-                getattr(self, v.name).units = "unknown"
-
-        self.idx = np.empty(shape=0)
+        for name, message in self.metadata.items():
+            setattr(self.dataset, name, message)
 
     def __del__(self):
-        if self.dataset:
-            self.dataset.close()
+        if self.to_export:
+            self.close()
 
-    def sync(self):
-        """Write all buffered data to disk"""
-        self.dataset.sync()
+    def close(self):
+        """Close the ParticleFile object by exporting and then deleting
+        the temporary npy files"""
+        self.export()
+        self.delete_tempwritedir()
+        self.dataset.close()
+        self.to_export = False
 
     def add_metadata(self, name, message):
         """Add metadata to :class:`parcels.particleset.ParticleSet`
@@ -154,54 +189,149 @@ class ParticleFile(object):
         :param message: message to be written
         """
         if self.dataset is None:
-            self.open_dataset()
-        setattr(self.dataset, name, message)
-
-    def write(self, pset, time, sync=True, deleted_only=False):
-        """Write :class:`parcels.particleset.ParticleSet` data to file
+            self.metadata[name] = message
+        else:
+            setattr(self.dataset, name, message)
 
+    def convert_pset_to_dict(self, pset, time, deleted_only=False):
+        """Convert all Particle data from one time step to a python dictionary.
         :param pset: ParticleSet object to write
         :param time: Time at which to write ParticleSet
-        :param sync: Optional argument whether to write data to disk immediately. Default is True
-
+        :param deleted_only: Flag to write only the deleted Particles
+        returns two dictionaries: one for all variables to be written each outputdt,
+         and one for all variables to be written once
         """
-        if self.dataset is None:
-            self.open_dataset()
-        if isinstance(time, delta):
-            time = time.total_seconds()
-        if self.lasttime_written != time and \
-           (self.write_ondelete is False or deleted_only is True):
-            if pset.size > 0:
+        data_dict = {}
+        data_dict_once = {}
 
-                first_write = [p for p in pset if (p.fileid < 0 or len(self.idx) == 0) and _is_particle_started_yet(p, time)]  # len(self.idx)==0 in case pset is written to new ParticleFile
-                for p in first_write:
-                    p.fileid = self.lasttraj  # particle id in current file
-                    self.lasttraj += 1
+        time = time.total_seconds() if isinstance(time, delta) else time
 
-                self.idx = np.append(self.idx, np.zeros(len(first_write)))
+        if self.lasttime_written != time and \
+           (self.write_ondelete is False or deleted_only is True):
+            if pset.size == 0:
+                logger.warning("ParticleSet is empty on writing as array at time %g" % time)
+            else:
+                for var in self.var_names:
+                    data_dict[var] = np.nan * np.zeros(len(pset))
 
+                i = 0
                 for p in pset:
-                    if _is_particle_started_yet(p, time):
-                        i = p.fileid
-                        self.id[i, self.idx[i]] = p.id
-                        self.time[i, self.idx[i]] = p.time
-                        self.lat[i, self.idx[i]] = p.lat
-                        self.lon[i, self.idx[i]] = p.lon
-                        self.z[i, self.idx[i]] = p.depth
-                        for var in self.user_vars:
-                            getattr(self, var)[i, self.idx[i]] = getattr(p, var)
+                    if p.dt*p.time <= p.dt*time:
+                        for var in self.var_names:
+                            data_dict[var][i] = getattr(p, var)
                         if p.state != ErrorCode.Delete and not np.allclose(p.time, time):
                             logger.warning_once('time argument in pfile.write() is %g, but a particle has time %g.' % (time, p.time))
-
-                for p in first_write:
-                    for var in self.user_vars_once:
-                        getattr(self, var)[p.fileid] = getattr(p, var)
-            else:
-                logger.warning("ParticleSet is empty on writing as array at time %g" % time)
+                        self.maxid_written = np.max([self.maxid_written, p.id])
+                        i += 1
+
+                save_ind = np.isfinite(data_dict["id"])
+                for key in self.var_names:
+                    data_dict[key] = data_dict[key][save_ind]
+
+                if time not in self.time_written:
+                    self.time_written.append(time)
+
+                if len(self.var_names_once) > 0:
+                    first_write = [p for p in pset if (p.id not in self.written_once) and _is_particle_started_yet(p, time)]
+                    data_dict_once['id'] = np.nan * np.zeros(len(first_write))
+                    for var in self.var_names_once:
+                        data_dict_once[var] = np.nan * np.zeros(len(first_write))
+
+                    i = 0
+                    for p in first_write:
+                        self.written_once.append(p.id)
+                        data_dict_once['id'][i] = p.id
+                        for var in self.var_names_once:
+                            data_dict_once[var][i] = getattr(p, var)
+                        i += 1
 
             if not deleted_only:
-                self.idx += 1
                 self.lasttime_written = time
 
-        if sync:
-            self.sync()
+        return data_dict, data_dict_once
+
+    def dump_dict_to_npy(self, data_dict, data_dict_once):
+        """Buffer data to set of temporary numpy files, using np.save"""
+
+        if not os.path.exists(self.tempwritedir):
+            os.mkdir(self.tempwritedir)
+
+        if len(data_dict) > 0:
+            tmpfilename = os.path.join(self.tempwritedir, str(len(self.file_list)) + ".npy")
+            with open(tmpfilename, 'wb') as f:
+                np.save(f, data_dict)
+            self.file_list.append(tmpfilename)
+
+        if len(data_dict_once) > 0:
+            tmpfilename = os.path.join(self.tempwritedir, str(len(self.file_list)) + '_once.npy')
+            with open(tmpfilename, 'wb') as f:
+                np.save(f, data_dict_once)
+            self.file_list_once.append(tmpfilename)
+
+        self.to_export = True
+
+    def dump_psetinfo_to_npy(self):
+        pset_info = {}
+        attrs_to_dump = ['name', 'var_names', 'var_names_once', 'time_origin', 'lonlatdepth_dtype',
+                         'file_list', 'file_list_once', 'maxid_written', 'time_written', 'parcels_mesh',
+                         'metadata']
+        for a in attrs_to_dump:
+            if hasattr(self, a):
+                pset_info[a] = getattr(self, a)
+        with open(os.path.join(self.tempwritedir, 'pset_info.npy'), 'wb') as f:
+            np.save(f, pset_info)
+
+    def write(self, pset, time, deleted_only=False):
+        """Write all data from one time step to a temporary npy-file
+        using a python dictionary. The data is saved in the folder 'out'.
+        :param pset: ParticleSet object to write
+        :param time: Time at which to write ParticleSet
+        :param deleted_only: Flag to write only the deleted Particles
+        """
+
+        data_dict, data_dict_once = self.convert_pset_to_dict(pset, time, deleted_only=deleted_only)
+        self.dump_dict_to_npy(data_dict, data_dict_once)
+        self.dump_psetinfo_to_npy()
+
+    def read_from_npy(self, file_list, time_steps, var):
+        """Read NPY-files for one variable using a loop over all files.
+        :param file_list: List that  contains all file names in the output directory
+        :param time_steps: Number of time steps that were written in out directory
+        :param var: name of the variable to read
+        """
+
+        data = np.nan * np.zeros((self.maxid_written+1, time_steps))
+        time_index = np.zeros(self.maxid_written+1, dtype=int)
+        t_ind_used = np.zeros(time_steps, dtype=int)
+
+        # loop over all files
+        for npyfile in file_list:
+            data_dict = np.load(npyfile, allow_pickle=True).item()
+            id_ind = np.array(data_dict["id"], dtype=int)
+            t_ind = time_index[id_ind] if 'once' not in file_list[0] else 0
+            t_ind_used[t_ind] = 1
+            data[id_ind, t_ind] = data_dict[var]
+            time_index[id_ind] = time_index[id_ind] + 1
+
+        # remove rows and columns that are completely filled with nan values
+        tmp = data[time_index > 0, :]
+        return tmp[:, t_ind_used == 1]
+
+    def export(self):
+        """Exports outputs in temporary NPY-files to NetCDF file"""
+
+        for var in self.var_names:
+            data = self.read_from_npy(self.file_list, len(self.time_written), var)
+            if var == self.var_names[0]:
+                self.open_netcdf_file(data.shape)
+            varout = 'z' if var == 'depth' else var
+            getattr(self, varout)[:, :] = data
+
+        if len(self.var_names_once) > 0:
+            for var in self.var_names_once:
+                getattr(self, var)[:] = self.read_from_npy(self.file_list_once, 1, var)
+
+    def delete_tempwritedir(self):
+        """Deleted all temporary npy files"""
+        if os.path.exists(self.tempwritedir):
+            shutil.rmtree(self.tempwritedir)
diff --git a/parcels/particleset.py b/parcels/particleset.py
index 58a3e40b..9002abfc 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -411,6 +411,10 @@ class ParticleSet(object):
         while (time < endtime and dt > 0) or (time > endtime and dt < 0) or dt == 0:
             if verbose_progress is None and time_module.time() - walltime_start > 10:
                 # Showing progressbar if runtime > 10 seconds
+                if output_file:
+                    logger.info('Temporary output files are stored in %s.' % output_file.tempwritedir)
+                    logger.info('You can use "parcels_convert_npydir_to_netcdf %s" to convert these '
+                                'to a NetCDF file during the run.' % output_file.tempwritedir)
                 pbar = progressbar.ProgressBar(max_value=abs(endtime - _starttime)).start()
                 verbose_progress = True
             if dt > 0:
diff --git a/parcels/scripts/convert_npydir_to_netcdf.py b/parcels/scripts/convert_npydir_to_netcdf.py
new file mode 100644
index 00000000..0413b1c3
--- /dev/null
+++ b/parcels/scripts/convert_npydir_to_netcdf.py
@@ -0,0 +1,35 @@
+from parcels import ParticleFile
+import numpy as np
+from os import path
+from argparse import ArgumentParser
+
+
+def convert_npydir_to_netcdf(tempwritedir):
+    """Convert npy files in tempwritedir to a NetCDF file
+    :param tempwritedir: directory where the temporary npy files are stored (can be obtained from ParticleFile.tempwritedir attribute)
+    """
+
+    pyset_file = path.join(tempwritedir, 'pset_info.npy')
+    if not path.isdir(tempwritedir):
+        raise ValueError('Output directory "%s" does not exist' % tempwritedir)
+    if not path.isfile(pyset_file):
+        raise ValueError('Output directory "%s" does not contain a pset_info.npy file' % tempwritedir)
+
+    pset_info = np.load(pyset_file, allow_pickle=True).item()
+    pfile = ParticleFile(None, None, pset_info=pset_info)
+    pfile.export()
+    pfile.dataset.close()
+
+
+def main(tempwritedir=None):
+    if tempwritedir is None:
+        p = ArgumentParser(description="""Script to convert temporary npy output files to NetCDF""")
+        p.add_argument('tempwritedir', help='Name of directory where temporary npy files are stored')
+        args = p.parse_args()
+        tempwritedir = args.tempwritedir
+
+    convert_npydir_to_netcdf(tempwritedir)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/setup.py b/setup.py
index 6c1047f8..d51bf82f 100644
--- a/setup.py
+++ b/setup.py
@@ -15,5 +15,6 @@ setup(name='parcels',
       package_data={'parcels': ['include/*',
                                 'examples/*']},
       entry_points={'console_scripts': [
-          'parcels_get_examples = parcels.scripts.get_examples:main']}
+          'parcels_get_examples = parcels.scripts.get_examples:main',
+          'parcels_convert_npydir_to_netcdf = parcels.scripts.convert_npydir_to_netcdf:main']}
       )


Test patch:
diff --git a/tests/test_particle_file.py b/tests/test_particle_file.py
index 2980fdd3..61f52bc3 100644
--- a/tests/test_particle_file.py
+++ b/tests/test_particle_file.py
@@ -1,7 +1,215 @@
+from parcels import (FieldSet, ParticleSet, ScipyParticle, JITParticle,
+                     Variable, ErrorCode)
 from parcels.particlefile import _set_calendar
 from parcels.tools.converters import _get_cftime_calendars, _get_cftime_datetimes
+import numpy as np
+import pytest
+import os
+from netCDF4 import Dataset
 import cftime
 
+ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
+
+
+@pytest.fixture
+def fieldset(xdim=40, ydim=100):
+    U = np.zeros((ydim, xdim), dtype=np.float32)
+    V = np.zeros((ydim, xdim), dtype=np.float32)
+    lon = np.linspace(0, 1, xdim, dtype=np.float32)
+    lat = np.linspace(-60, 60, ydim, dtype=np.float32)
+    depth = np.zeros(1, dtype=np.float32)
+    data = {'U': np.array(U, dtype=np.float32), 'V': np.array(V, dtype=np.float32)}
+    dimensions = {'lat': lat, 'lon': lon, 'depth': depth}
+    return FieldSet.from_data(data, dimensions)
+
+
+def close_and_compare_netcdffiles(filepath, ofile, assystemcall=False):
+    if assystemcall:
+        os.system('parcels_convert_npydir_to_netcdf %s' % ofile.tempwritedir)
+    else:
+        import parcels.scripts.convert_npydir_to_netcdf as convert
+        convert.convert_npydir_to_netcdf(ofile.tempwritedir)
+
+    ncfile1 = Dataset(filepath, 'r', 'NETCDF4')
+
+    ofile.name = filepath + 'b.nc'
+    ofile.export()
+    ofile.dataset.close()
+    ncfile2 = Dataset(filepath + 'b.nc', 'r', 'NETCDF4')
+
+    for v in ncfile2.variables.keys():
+        assert np.allclose(ncfile1.variables[v][:], ncfile2.variables[v][:])
+
+    for a in ncfile2.ncattrs():
+        assert getattr(ncfile1, a) == getattr(ncfile2, a)
+
+    ncfile2.close()
+    return ncfile1
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_pfile_array_remove_particles(fieldset, mode, tmpdir, npart=10):
+    filepath = tmpdir.join("pfile_array_remove_particles.nc")
+    pset = ParticleSet(fieldset, pclass=ptype[mode],
+                       lon=np.linspace(0, 1, npart),
+                       lat=0.5*np.ones(npart))
+    pfile = pset.ParticleFile(filepath)
+    pfile.write(pset, 0)
+    pset.remove(3)
+    pfile.write(pset, 1)
+    ncfile = close_and_compare_netcdffiles(filepath, pfile)
+    ncfile.close()
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_pfile_array_remove_all_particles(fieldset, mode, tmpdir, npart=10):
+
+    filepath = tmpdir.join("pfile_array_remove_particles.nc")
+    pset = ParticleSet(fieldset, pclass=ptype[mode],
+                       lon=np.linspace(0, 1, npart),
+                       lat=0.5*np.ones(npart))
+    pfile = pset.ParticleFile(filepath)
+    pfile.write(pset, 0)
+    for _ in range(npart):
+        pset.remove(-1)
+    pfile.write(pset, 1)
+    pfile.write(pset, 2)
+    ncfile = close_and_compare_netcdffiles(filepath, pfile)
+    ncfile.close()
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+@pytest.mark.parametrize('assystemcall', [True, False])
+def test_variable_written_ondelete(fieldset, mode, tmpdir, assystemcall, npart=3):
+    filepath = tmpdir.join("pfile_on_delete_written_variables.nc")
+
+    def move_west(particle, fieldset, time):
+        tmp = fieldset.U[time, particle.depth, particle.lat, particle.lon]  # to trigger out-of-bounds error
+        particle.lon -= 0.1 + tmp
+
+    def DeleteP(particle, fieldset, time):
+        particle.delete()
+
+    lon = np.linspace(0.05, 0.95, npart)
+    lat = np.linspace(0.95, 0.05, npart)
+
+    (dt, runtime) = (0.1, 0.8)
+    lon_end = lon - runtime/dt*0.1
+    noutside = len(lon_end[lon_end < 0])
+
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
+
+    outfile = pset.ParticleFile(name=filepath, write_ondelete=True)
+    outfile.add_metadata('runtime', runtime)
+    pset.execute(move_west, runtime=runtime, dt=dt, output_file=outfile,
+                 recovery={ErrorCode.ErrorOutOfBounds: DeleteP})
+
+    ncfile = close_and_compare_netcdffiles(filepath, outfile, assystemcall=assystemcall)
+    assert ncfile.runtime == runtime
+    lon = ncfile.variables['lon'][:]
+    assert (lon.size == noutside)
+    ncfile.close()
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_variable_write_double(fieldset, mode, tmpdir):
+    filepath = tmpdir.join("pfile_variable_write_double.nc")
+
+    def Update_lon(particle, fieldset, time):
+        particle.lon += 0.1
+
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
+    ofile = pset.ParticleFile(name=filepath, outputdt=0.1)
+    pset.execute(pset.Kernel(Update_lon), endtime=1, dt=0.1, output_file=ofile)
+
+    ncfile = close_and_compare_netcdffiles(filepath, ofile)
+    lons = ncfile.variables['lon'][:]
+    assert (isinstance(lons[0, 0], np.float64))
+    ncfile.close()
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+@pytest.mark.parametrize('npart', [1, 2, 5])
+def test_variable_written_once(fieldset, mode, tmpdir, npart):
+    filepath = tmpdir.join("pfile_once_written_variables.nc")
+
+    def Update_v(particle, fieldset, time):
+        particle.v_once += 1.
+        particle.age += particle.dt
+
+    class MyParticle(ptype[mode]):
+        v_once = Variable('v_once', dtype=np.float64, initial=0., to_write='once')
+        age = Variable('age', dtype=np.float32, initial=0.)
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
+    time = np.arange(0, npart/10., 0.1, dtype=np.float64)
+    pset = ParticleSet(fieldset, pclass=MyParticle, lon=lon, lat=lat, time=time, v_once=time)
+    ofile = pset.ParticleFile(name=filepath, outputdt=0.1)
+    pset.execute(pset.Kernel(Update_v), endtime=1, dt=0.1, output_file=ofile)
+
+    assert np.allclose([p.v_once - vo - p.age*10 for p, vo in zip(pset, time)], 0, atol=1e-5)
+    ncfile = close_and_compare_netcdffiles(filepath, ofile)
+    vfile = np.ma.filled(ncfile.variables['v_once'][:], np.nan)
+    assert (vfile.shape == (npart, ))
+    assert np.allclose(vfile, time)
+    ncfile.close()
+
+
+@pytest.mark.parametrize('type', ['repeatdt', 'timearr'])
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+@pytest.mark.parametrize('repeatdt', range(1, 3))
+@pytest.mark.parametrize('dt', [-1, 1])
+@pytest.mark.parametrize('maxvar', [2, 4, 10])
+def test_pset_repeated_release_delayed_adding_deleting(type, fieldset, mode, repeatdt, tmpdir, dt, maxvar, runtime=10):
+    fieldset.maxvar = maxvar
+
+    class MyParticle(ptype[mode]):
+        sample_var = Variable('sample_var', initial=0.)
+    if type == 'repeatdt':
+        pset = ParticleSet(fieldset, lon=[0], lat=[0], pclass=MyParticle, repeatdt=repeatdt)
+    elif type == 'timearr':
+        pset = ParticleSet(fieldset, lon=np.zeros(runtime), lat=np.zeros(runtime), pclass=MyParticle, time=list(range(runtime)))
+    outfilepath = tmpdir.join("pfile_repeated_release.nc")
+    pfile = pset.ParticleFile(outfilepath, outputdt=abs(dt))
+
+    def IncrLon(particle, fieldset, time):
+        particle.sample_var += 1.
+        if particle.sample_var > fieldset.maxvar:
+            particle.delete()
+    for i in range(runtime):
+        pset.execute(IncrLon, dt=dt, runtime=1., output_file=pfile)
+
+    ncfile = close_and_compare_netcdffiles(outfilepath, pfile)
+    samplevar = ncfile.variables['sample_var'][:]
+    if type == 'repeatdt':
+        assert samplevar.shape == (runtime // repeatdt+1, min(maxvar+1, runtime)+1)
+        assert np.allclose([p.sample_var for p in pset], np.arange(maxvar, -1, -repeatdt))
+    elif type == 'timearr':
+        assert samplevar.shape == (runtime, min(maxvar + 1, runtime) + 1)
+    # test whether samplevar[:, k] = k
+    for k in range(samplevar.shape[1]):
+        assert np.allclose([p for p in samplevar[:, k] if np.isfinite(p)], k)
+    filesize = os.path.getsize(str(outfilepath))
+    assert filesize < 1024 * 65  # test that chunking leads to filesize less than 65KB
+    ncfile.close()
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_write_timebackward(fieldset, mode, tmpdir):
+    outfilepath = tmpdir.join("pfile_write_timebackward.nc")
+
+    def Update_lon(particle, fieldset, time):
+        particle.lon -= 0.1 * particle.dt
+
+    pset = ParticleSet(fieldset, pclass=JITParticle, lat=np.linspace(0, 1, 3), lon=[0, 0, 0],
+                       time=[1, 2, 3])
+    pfile = pset.ParticleFile(name=outfilepath, outputdt=1.)
+    pset.execute(pset.Kernel(Update_lon), runtime=4, dt=-1.,
+                 output_file=pfile)
+    ncfile = close_and_compare_netcdffiles(outfilepath, pfile)
+    trajs = ncfile.variables['trajectory'][:, 0]
+    assert np.all(np.diff(trajs) > 0)  # all particles written in order of traj ID
+
 
 def test_set_calendar():
     for calendar_name, cf_datetime in zip(_get_cftime_calendars(), _get_cftime_datetimes()):
diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py
index cb8e43a2..8e43f52b 100644
--- a/tests/test_particle_sets.py
+++ b/tests/test_particle_sets.py
@@ -2,8 +2,6 @@ from parcels import (FieldSet, ParticleSet, Field, ScipyParticle, JITParticle,
                      Variable, ErrorCode, CurvilinearZGrid)
 import numpy as np
 import pytest
-import os
-from netCDF4 import Dataset
 
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
@@ -132,44 +130,6 @@ def test_pset_repeated_release(fieldset, mode, npart=10):
     assert np.allclose([p.lon for p in pset], np.arange(npart, 0, -1))
 
 
-@pytest.mark.parametrize('type', ['repeatdt', 'timearr'])
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-@pytest.mark.parametrize('repeatdt', range(1, 3))
-@pytest.mark.parametrize('dt', [-1, 1])
-@pytest.mark.parametrize('maxvar', [2, 4, 10])
-def test_pset_repeated_release_delayed_adding_deleting(type, fieldset, mode, repeatdt, tmpdir, dt, maxvar, runtime=10):
-    fieldset.maxvar = maxvar
-
-    class MyParticle(ptype[mode]):
-        sample_var = Variable('sample_var', initial=0.)
-    if type == 'repeatdt':
-        pset = ParticleSet(fieldset, lon=[0], lat=[0], pclass=MyParticle, repeatdt=repeatdt)
-    elif type == 'timearr':
-        pset = ParticleSet(fieldset, lon=np.zeros(runtime), lat=np.zeros(runtime), pclass=MyParticle, time=list(range(runtime)))
-    outfilepath = tmpdir.join("pfile_repeated_release")
-    pfile = pset.ParticleFile(outfilepath, outputdt=abs(dt))
-
-    def IncrLon(particle, fieldset, time):
-        particle.sample_var += 1.
-        if particle.sample_var > fieldset.maxvar:
-            particle.delete()
-    for i in range(runtime):
-        pset.execute(IncrLon, dt=dt, runtime=1., output_file=pfile)
-    ncfile = Dataset(outfilepath+".nc", 'r', 'NETCDF4')
-    samplevar = ncfile.variables['sample_var'][:]
-    ncfile.close()
-    if type == 'repeatdt':
-        assert samplevar.shape == (runtime // repeatdt+1, min(maxvar+1, runtime)+1)
-        assert np.allclose([p.sample_var for p in pset], np.arange(maxvar, -1, -repeatdt))
-    elif type == 'timearr':
-        assert samplevar.shape == (runtime, min(maxvar + 1, runtime) + 1)
-    # test whether samplevar[:, k] = k
-    for k in range(samplevar.shape[1]):
-        assert np.allclose([p for p in samplevar[:, k] if np.isfinite(p)], k)
-    filesize = os.path.getsize(str(outfilepath+".nc"))
-    assert filesize < 1024 * 65  # test that chunking leads to filesize less than 65KB
-
-
 def test_pset_repeatdt_check_dt(fieldset):
     pset = ParticleSet(fieldset, lon=[0], lat=[0], pclass=ScipyParticle, repeatdt=5)
 
@@ -217,8 +177,7 @@ def test_pset_custom_ptype(fieldset, mode, npart=100):
                        lon=np.linspace(0, 1, npart),
                        lat=np.linspace(1, 0, npart))
     assert(pset.size == 100)
-    # FIXME: The float test fails with a conversion error of 1.e-8
-    # assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), rtol=1e-12)
+    assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), rtol=1e-12)
     assert np.allclose([p.n - 2 for p in pset], np.zeros(npart), rtol=1e-12)
 
 
@@ -372,101 +331,6 @@ def test_density(fieldset, mode, area_scale):
             assert np.allclose(fieldset.U.lat[inds[0][i]], pset[i].lat, atol=fieldset.U.lat[1]-fieldset.U.lat[0])
 
 
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_pfile_array_remove_particles(fieldset, mode, tmpdir, npart=10):
-    filepath = tmpdir.join("pfile_array_remove_particles")
-    pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart),
-                       lat=0.5*np.ones(npart))
-    pfile = pset.ParticleFile(filepath)
-    pfile.write(pset, 0)
-    pset.remove(3)
-    pfile.write(pset, 1)
-
-
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_pfile_array_remove_all_particles(fieldset, mode, tmpdir, npart=10):
-
-    filepath = tmpdir.join("pfile_array_remove_particles")
-    pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart),
-                       lat=0.5*np.ones(npart))
-    pfile = pset.ParticleFile(filepath)
-    pfile.write(pset, 0)
-    for _ in range(npart):
-        pset.remove(-1)
-    pfile.write(pset, 1)
-    pfile.write(pset, 2)
-
-
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-@pytest.mark.parametrize('npart', [1, 2, 5])
-def test_variable_written_once(fieldset, mode, tmpdir, npart):
-    filepath = tmpdir.join("pfile_once_written_variables")
-
-    def Update_v(particle, fieldset, time):
-        particle.v_once += 1.
-        particle.age += particle.dt
-
-    class MyParticle(ptype[mode]):
-        v_once = Variable('v_once', dtype=np.float32, initial=0., to_write='once')
-        age = Variable('age', dtype=np.float32, initial=0.)
-    lon = np.linspace(0, 1, npart)
-    lat = np.linspace(1, 0, npart)
-    pset = ParticleSet(fieldset, pclass=MyParticle, lon=lon, lat=lat, repeatdt=0.1)
-    pset.execute(pset.Kernel(Update_v), endtime=1, dt=0.1,
-                 output_file=pset.ParticleFile(name=filepath, outputdt=0.1))
-    assert np.allclose([p.v_once - p.age * 10 for p in pset], 0, atol=1e-5)
-    ncfile = Dataset(filepath+".nc", 'r', 'NETCDF4')
-    vfile = ncfile.variables['v_once'][:]
-    assert (vfile.shape == (npart*11, ))
-    assert [v == 0 for v in vfile]
-
-
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_variable_write_double(fieldset, mode, tmpdir):
-    filepath = tmpdir.join("pfile_variable_write_double")
-
-    def Update_lon(particle, fieldset, time):
-        particle.lon += 0.1
-
-    pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
-    pset.execute(pset.Kernel(Update_lon), endtime=1, dt=0.1,
-                 output_file=pset.ParticleFile(name=filepath, outputdt=0.1))
-    ncfile = Dataset(filepath+".nc", 'r', 'NETCDF4')
-    lons = ncfile.variables['lon'][:]
-    assert (isinstance(lons[0, 0], np.float64))
-
-
-@pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_variable_written_ondelete(fieldset, mode, tmpdir, npart=3):
-    filepath = tmpdir.join("pfile_on_delete_written_variables")
-
-    def move_west(particle, fieldset, time):
-        tmp = fieldset.U[time, particle.depth, particle.lat, particle.lon]  # to trigger out-of-bounds error
-        particle.lon -= 0.1 + tmp
-
-    def DeleteP(particle, fieldset, time):
-        particle.delete()
-
-    lon = np.linspace(0.05, 0.95, npart)
-    lat = np.linspace(0.95, 0.05, npart)
-
-    (dt, runtime) = (0.1, 0.8)
-    lon_end = lon - runtime/dt*0.1
-    noutside = len(lon_end[lon_end < 0])
-
-    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
-
-    outfile = pset.ParticleFile(name=filepath, write_ondelete=True)
-    outfile.add_metadata('runtime', runtime)
-    pset.execute(move_west, runtime=runtime, dt=dt, output_file=outfile,
-                 recovery={ErrorCode.ErrorOutOfBounds: DeleteP})
-    ncfile = Dataset(filepath+".nc", 'r', 'NETCDF4')
-    lon = ncfile.variables['lon'][:]
-    assert (lon.size == noutside)
-
-
 @pytest.mark.parametrize('staggered_grid', ['Agrid', 'Cgrid'])
 def test_from_field_exact_val(staggered_grid):
     xdim = 4
diff --git a/tests/test_scripts.py b/tests/test_scripts.py
index e10ac314..ff5d9acf 100644
--- a/tests/test_scripts.py
+++ b/tests/test_scripts.py
@@ -24,12 +24,14 @@ def create_outputfiles(dir):
     output_file = pset.ParticleFile(name=fp, outputdt=delaytime)
 
     for t in range(npart):
-        pset.add(JITParticle(lon=x, lat=lat[t], fieldset=fieldset))
+        time = 0 if len(pset) == 0 else pset[0].time
+        pset.add(JITParticle(lon=x, lat=lat[t], fieldset=fieldset, time=time))
         pset.execute(AdvectionRK4, runtime=delaytime, dt=delta(minutes=5),
                      output_file=output_file)
 
     pset.execute(AdvectionRK4, runtime=endtime-npart*delaytime,
                  dt=delta(minutes=5), output_file=output_file)
+    output_file.close()
 
     return fp