Repository: OceanParcels/Parcels (language: python)

Issue:
Documentation revamp
I'll use this issue to outline the ideas I have for the documentation, so that discussion can be centralised here with PRs linking back to this. Any of these points are up for discussion. I'm happy to work on all of these (although may need some additional help with config). If anyone else would like to contribute to a bulletpoint, drop a comment 🙂 (although there's no expectation on my end).

As these changes get implemented, I will edit this description adding the PR numbers to the bullet points.

# Changes
## Major changes
- [x] Documentation restructure
  - Add sectioning to the documentation, breaking apart the documentation into its conceptual components so that it isn't one large page.
- [x] Migrate to PyData theme (#1321)
  - The PyData theme is a much cleaner theme, with more features.
- [x] Create contribution guidelines for `parcels`
  - Having guidelines that contributors can follow would be good to communicate the standards adopted by the codebase.
  - The contributing [guidelines from geopandas](https://geopandas.org/en/stable/community/contributing.html) might be good to copy and slightly adapt to our purposes (ie. we only support Python 3.8+ (@erikvansebille it would be good to confirm this, I'm going off the `environment_*.yml` files), we don't use `black`, and [commit message conventions](https://geopandas.org/en/stable/community/contributing.html#commit-message-conventions) might be a bit extra for the size of this project).
  - The contribution guidelines can be its own page on the documentation website.
  - Perhaps [oceanparcels.org](https://oceanparcels.org) can be updated to link to the contributing page for dev installs.
- [x] Use [`flake8-docstrings`](https://pypi.org/project/flake8-docstrings/) to monitor docstring quality (#1321)
- [x] Migrate docstring format to the [Numpy Docstring Standard](https://numpydoc.readthedocs.io/en/latest/format.html#docstring-standard)
  - The `rst` format is quite difficult to read when displayed in IDEs.
  - Due to this ^, the Numpy Docstring Standard has been adopted by many other scientific packages (numpy, pandas, geopandas, xarray, to name a few).
  - Existing docstrings can be easily converted from `rst` format into Numpy format using [`pyment`](https://pypi.org/project/pyment/).
  - Compiling and [validation](https://numpydoc.readthedocs.io/en/latest/validation.html) of Numpy docstrings is possible in the build process through Sphinx with the `numpydoc` package. This validation is in extension to `flake8-docstrings`, making sure that the docstring is in the Numpy standard format.
  - The validation step is particularly important to make sure breaking documentation changes fail CI/CD on pull requests, allowing them to be changed before being merged into main (the ReadTheDocs integration feature just below will greatly help with viewing documentation changes, and failing builds that break validation rules).
- [x] ReadTheDocs build on pull requests
  - ReadTheDocs can [build on pull requests, and report build status back to GitHub](https://docs.readthedocs.io/en/stable/pull-requests.html#features). This can be [set up in the ReadTheDocs panel](https://docs.readthedocs.io/en/stable/guides/pull-requests.html) (@erikvansebille your help here would be great).
- [x] Type annotations and static typechecking
  - Having type annotations will greatly help with documentation, and reduce ambiguity amongst users on the types of different parameters.
  - Using `mypy` to do static typechecking will validate type annotations, making sure internally `parcels` is correctly typed.
- [x] Improve support for documentation notebooks
  - Automatically link functions/classes to example notebooks in which they are called.
  - Generally increase the visibility of the example notebooks (eg. perhaps creating a gallery; although [oceanparcels.org](https://oceanparcels.org) does a good job of that already).
  - Still need to explore the viability of this feature, and how it will be implemented 🤔

## Minor changes
- [x] Have "view source" links in documentation link directly back to GitHub (#1321)
- [x] Disable search engine indexing of old documentation

## Bugfixes
- [x] ReadTheDocs + Pydata theme compatability
  - Currently on [my ReadTheDocs (`docs` branch)](https://veckoparcels.readthedocs.io/en/docs/) the PyData theme makes the dropdown (to switch between the different versions) completely disappear.

# Notes
PRs will be made into the `new_docs` branch to act as a staging area before they are rolled out to the community. Once this issue is completed, we'll merge `new_docs` to `master`.

---
The type annotation and docstring specific changes will mainly be applied to the following modules which are currently part of the documentation website as defined in the `.. automodule::` calls in `docs/index.rst`:

- `parcels.particleset.particlesetsoa`
- `parcels.fieldset`
- `parcels.field`
- `parcels.gridset`
- `parcels.grid`
- `parcels.particle`
- `parcels.application_kernels.advection`
- `parcels.application_kernels.advectiondiffusion`
- `parcels.application_kernels.EOSseawaterproperties`
- `parcels.application_kernels.TEOSseawaterdensity`
- `parcels.compilation`
- `parcels.interaction`
- `parcels.interaction.neighborsearch`
- `parcels.kernel.basekernel`
- `parcels.kernel.kernelaos`
- `parcels.kernel.kernelsoa`
- `parcels.particlefile.baseparticlefile`
- `parcels.particlefile.particlefileaos`
- `parcels.particlefile.particlefilesoa`
- `parcels.rng`
- `parcels.particleset.baseparticleset`
- `parcels.particleset.particlesetaos`
- `parcels.tools.statuscodes`
- `parcels.tools.converters`
- `parcels.tools.loggers`
- `parcels.plotting`
- `parcels.scripts.plottrajectoriesfile`
- `parcels.scripts.get_examples`



Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5405bdd93..79c46106f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -10,6 +10,7 @@ repos:
     rev: '6.0.0'
     hooks:
     -   id: flake8
+        additional_dependencies: ["flake8-docstrings"]
 -   repo: https://github.com/pycqa/isort
     rev: 5.12.0
     hooks:
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index e6d49aa27..05ea1a4fa 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -3,5 +3,9 @@ build:
   os: ubuntu-22.04
   tools:
     python: mambaforge-4.10
+
+sphinx:
+  configuration: docs/conf.py
+
 conda:
   environment: environment_py3_linux.yml
diff --git a/docs/_static/parcelslogo-inverted.png b/docs/_static/parcelslogo-inverted.png
new file mode 100644
index 000000000..89aa520a3
Binary files /dev/null and b/docs/_static/parcelslogo-inverted.png differ
diff --git a/docs/parcelslogo.png b/docs/_static/parcelslogo.png
similarity index 100%
rename from docs/parcelslogo.png
rename to docs/_static/parcelslogo.png
diff --git a/docs/conf.py b/docs/conf.py
index d304cfd89..abbdd8f86 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -11,8 +11,10 @@
 # All configuration values have a default; values that are commented out
 # serve to show the default.
 
+import inspect
 import os
 import sys
+import warnings
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
@@ -30,7 +32,7 @@
 extensions = [
     'sphinx.ext.autodoc',
     'sphinx.ext.todo',
-    'sphinx.ext.viewcode',
+    "sphinx.ext.linkcode",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -57,9 +59,11 @@
 # built documents.
 #
 # The short X.Y version.
-version = '2.4.0'
+import parcels  # noqa: E402
+
+version = parcels.__version__
 # The full version, including alpha/beta/rc tags.
-release = '2.4.0'
+release = version
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -110,7 +114,7 @@
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'bizstyle'
+html_theme = 'pydata_sphinx_theme'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
@@ -127,9 +131,6 @@
 # A shorter title for the navigation bar.  Default is the same as html_title.
 # html_short_title = None
 
-# The name of an image file (relative to this directory) to place at the top
-# of the sidebar.
-html_logo = 'parcelslogo.png'
 
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
@@ -140,6 +141,89 @@
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
 html_static_path = ['_static']
+html_theme_options = {
+    "logo": {
+        "image_light": "parcelslogo.png",
+        "image_dark": "parcelslogo-inverted.png",
+    },
+    "use_edit_page_button": True,
+    "github_url": "https://github.com/OceanParcels/parcels",
+    "icon_links": [
+        {
+            "name": "Conda Forge",
+            "url": "https://anaconda.org/conda-forge/parcels",  # required
+            "icon": "fa-solid fa-box",
+            "type": "fontawesome",
+        }
+    ]
+}
+
+html_context = {
+    "github_user": "OceanParcels",
+    "github_repo": "parcels",
+    "github_version": "master",
+    "doc_path": "docs",
+}
+
+
+# based on pandas doc/source/conf.py
+def linkcode_resolve(domain, info):
+    """Determine the URL corresponding to Python object."""
+    if domain != "py":
+        return None
+
+    modname = info["module"]
+    fullname = info["fullname"]
+
+    submod = sys.modules.get(modname)
+    if submod is None:
+        return None
+
+    obj = submod
+    for part in fullname.split("."):
+        try:
+            with warnings.catch_warnings():
+                # Accessing deprecated objects will generate noisy warnings
+                warnings.simplefilter("ignore", FutureWarning)
+                obj = getattr(obj, part)
+        except AttributeError:
+            return None
+
+    try:
+        fn = inspect.getsourcefile(inspect.unwrap(obj))
+    except TypeError:
+        try:  # property
+            fn = inspect.getsourcefile(inspect.unwrap(obj.fget))
+        except (AttributeError, TypeError):
+            fn = None
+    if not fn:
+        return None
+
+    try:
+        source, lineno = inspect.getsourcelines(obj)
+    except TypeError:
+        try:  # property
+            source, lineno = inspect.getsourcelines(obj.fget)
+        except (AttributeError, TypeError):
+            lineno = None
+    except OSError:
+        lineno = None
+
+    if lineno:
+        linespec = f"#L{lineno}-L{lineno + len(source) - 1}"
+    else:
+        linespec = ""
+
+    fn = os.path.relpath(fn, start=os.path.dirname(parcels.__file__))
+
+    if "-" in parcels.__version__:
+        return f"https://github.com/OceanParcels/parcels/blob/master/parcels/{fn}{linespec}"
+    else:
+        return (
+            f"https://github.com/OceanParcels/parcels/blob/"
+            f"{parcels.__version__}/parcels/{fn}{linespec}"
+        )
+
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
diff --git a/environment_py3_linux.yml b/environment_py3_linux.yml
index ee8d8f812..6b178a6a4 100644
--- a/environment_py3_linux.yml
+++ b/environment_py3_linux.yml
@@ -38,7 +38,10 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
+  - flake8-docstrings
   - pre_commit
 
   # Docs
   - sphinx<6
+  - pydata-sphinx-theme
+  - sphinx-autobuild
diff --git a/environment_py3_osx.yml b/environment_py3_osx.yml
index 9bf9c1028..09287a2f6 100644
--- a/environment_py3_osx.yml
+++ b/environment_py3_osx.yml
@@ -38,7 +38,10 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
+  - flake8-docstrings
   - pre_commit
 
   # Docs
   - sphinx<6
+  - pydata-sphinx-theme
+  - sphinx-autobuild
diff --git a/environment_py3_win.yml b/environment_py3_win.yml
index 9030b1b65..795f9c25b 100644
--- a/environment_py3_win.yml
+++ b/environment_py3_win.yml
@@ -35,7 +35,10 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
+  - flake8-docstrings
   - pre_commit
 
   # Docs
   - sphinx<6
+  - pydata-sphinx-theme
+  - sphinx-autobuild
diff --git a/parcels/application_kernels/EOSseawaterproperties.py b/parcels/application_kernels/EOSseawaterproperties.py
index baf2ff254..86e8142aa 100644
--- a/parcels/application_kernels/EOSseawaterproperties.py
+++ b/parcels/application_kernels/EOSseawaterproperties.py
@@ -1,4 +1,4 @@
-"""Collection of pre-built eos sea water property kernels"""
+"""Collection of pre-built eos sea water property kernels."""
 import math
 
 __all__ = ['PressureFromLatDepth', 'AdiabticTemperatureGradient', 'PtempFromTemp', 'TempFromPtemp', 'UNESCODensity']
@@ -106,7 +106,6 @@ def PtempFromTemp(particle, fieldset, time):
        Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6
 
     """
-
     s = fieldset.psu_salinity[time, particle.depth, particle.lat, particle.lon]
     t = fieldset.temperature[time, particle.depth, particle.lat, particle.lon]
     pres, pr = particle.pressure, fieldset.refpressure
diff --git a/parcels/application_kernels/TEOSseawaterdensity.py b/parcels/application_kernels/TEOSseawaterdensity.py
index e02d162bf..fa34abc31 100644
--- a/parcels/application_kernels/TEOSseawaterdensity.py
+++ b/parcels/application_kernels/TEOSseawaterdensity.py
@@ -1,12 +1,12 @@
-"""Collection of pre-built sea water density kernels"""
+"""Collection of pre-built sea water density kernels."""
 import math
 
 __all__ = ['PolyTEOS10_bsq']
 
 
 def PolyTEOS10_bsq(particle, fieldset, time):
-    '''
-    calculates density based on the polyTEOS10-bsq algorithm from Appendix A.2 of
+    """
+    Calculates density based on the polyTEOS10-bsq algorithm from Appendix A.2 of
     https://www.sciencedirect.com/science/article/pii/S1463500315000566
     requires fieldset.abs_salinity and fieldset.cons_temperature Fields in the fieldset
     and a particle.density Variable in the ParticleSet
@@ -19,8 +19,7 @@ def PolyTEOS10_bsq(particle, fieldset, time):
     Accurate and computationally efficient algorithms for potential
     temperature and density of seawater.  Journal of Atmospheric and
     Oceanic Technology, 20, 730-741.
-    '''
-
+    """
     Z = - math.fabs(particle.depth)  # Z needs to be negative
     SA = fieldset.abs_salinity[time, particle.depth, particle.lat, particle.lon]
     CT = fieldset.cons_temperature[time, particle.depth, particle.lat, particle.lon]
diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py
index 800e0f02e..c2a441c6c 100644
--- a/parcels/application_kernels/advection.py
+++ b/parcels/application_kernels/advection.py
@@ -1,4 +1,4 @@
-"""Collection of pre-built advection kernels"""
+"""Collection of pre-built advection kernels."""
 import math
 
 from parcels.tools.statuscodes import OperationCode
@@ -10,7 +10,8 @@
 def AdvectionRK4(particle, fieldset, time):
     """Advection of particles using fourth-order Runge-Kutta integration.
 
-    Function needs to be converted to Kernel object before execution"""
+    Function needs to be converted to Kernel object before execution.
+    """
     (u1, v1) = fieldset.UV[particle]
     lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt)
     (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
@@ -25,7 +26,8 @@ def AdvectionRK4(particle, fieldset, time):
 def AdvectionRK4_3D(particle, fieldset, time):
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity.
 
-    Function needs to be converted to Kernel object before execution"""
+    Function needs to be converted to Kernel object before execution.
+    """
     (u1, v1, w1) = fieldset.UVW[particle]
     lon1 = particle.lon + u1*.5*particle.dt
     lat1 = particle.lat + v1*.5*particle.dt
@@ -47,7 +49,8 @@ def AdvectionRK4_3D(particle, fieldset, time):
 def AdvectionEE(particle, fieldset, time):
     """Advection of particles using Explicit Euler (aka Euler Forward) integration.
 
-    Function needs to be converted to Kernel object before execution"""
+    Function needs to be converted to Kernel object before execution.
+    """
     (u1, v1) = fieldset.UV[particle]
     particle.lon += u1 * particle.dt
     particle.lat += v1 * particle.dt
@@ -58,7 +61,8 @@ def AdvectionRK45(particle, fieldset, time):
 
     Times-step dt is halved if error is larger than tolerance, and doubled
     if error is smaller than 1/10th of tolerance, with tolerance set to
-    1e-5 * dt by default."""
+    1e-5 * dt by default.
+    """
     rk45tol = 1e-5
     c = [1./4., 3./8., 12./13., 1., 1./2.]
     A = [[1./4., 0., 0., 0., 0.],
@@ -103,11 +107,12 @@ def AdvectionRK45(particle, fieldset, time):
 
 
 def AdvectionAnalytical(particle, fieldset, time):
-    """Advection of particles using 'analytical advection' integration
+    """Advection of particles using 'analytical advection' integration.
 
     Based on Ariane/TRACMASS algorithm, as detailed in e.g. Doos et al (https://doi.org/10.5194/gmd-10-1733-2017).
     Note that the time-dependent scheme is currently implemented with 'intermediate timesteps'
-    (default 10 per model timestep) and not yet with the full analytical time integration"""
+    (default 10 per model timestep) and not yet with the full analytical time integration.
+    """
     import numpy as np
 
     import parcels.tools.interpolation_utils as i_u
diff --git a/parcels/application_kernels/advectiondiffusion.py b/parcels/application_kernels/advectiondiffusion.py
index 7a9760d08..26d57b0b9 100644
--- a/parcels/application_kernels/advectiondiffusion.py
+++ b/parcels/application_kernels/advectiondiffusion.py
@@ -1,6 +1,7 @@
-"""Collection of pre-built advection-diffusion kernels
+"""Collection of pre-built advection-diffusion kernels.
 
-See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`_ for a detailed explanation"""
+See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`_ for a detailed explanation.
+"""
 import math
 
 import parcels.rng as ParcelsRandom
@@ -9,8 +10,7 @@
 
 
 def AdvectionDiffusionM1(particle, fieldset, time):
-    """Kernel for 2D advection-diffusion, solved using the Milstein scheme
-    at first order (M1).
+    """Kernel for 2D advection-diffusion, solved using the Milstein scheme at first order (M1).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
     and variable `fieldset.dres`, setting the resolution for the central
@@ -47,8 +47,7 @@ def AdvectionDiffusionM1(particle, fieldset, time):
 
 
 def AdvectionDiffusionEM(particle, fieldset, time):
-    """Kernel for 2D advection-diffusion, solved using the Euler-Maruyama
-    scheme (EM).
+    """Kernel for 2D advection-diffusion, solved using the Euler-Maruyama scheme (EM).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
     and variable `fieldset.dres`, setting the resolution for the central
diff --git a/parcels/application_kernels/interaction.py b/parcels/application_kernels/interaction.py
index 7b6ed62aa..3cd5cc268 100644
--- a/parcels/application_kernels/interaction.py
+++ b/parcels/application_kernels/interaction.py
@@ -1,4 +1,4 @@
-"""Collection of pre-built interaction kernels"""
+"""Collection of pre-built interaction kernels."""
 import numpy as np
 
 from parcels.tools.statuscodes import OperationCode, StateCode
@@ -8,7 +8,7 @@
 
 
 def NearestNeighborWithinRange(particle, fieldset, time, neighbors, mutator):
-    """Computes the nearest neighbor within range for each particle
+    """Computes the nearest neighbor within range for each particle.
 
     Particle has to have the nearest_neighbor property. If no particle
     is in range, set nearest_neighbor property to -1.
diff --git a/parcels/collection/collectionaos.py b/parcels/collection/collectionaos.py
index 66710b152..513b748c3 100644
--- a/parcels/collection/collectionaos.py
+++ b/parcels/collection/collectionaos.py
@@ -39,7 +39,6 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, p
         :param ngrid: number of grids in the fieldset of the overarching ParticleSet - required for initialising the
         field references of the ctypes-link of particles that are allocated
         """
-
         super(ParticleCollection, self).__init__()
 
         assert pid_orig is not None, "particle IDs are None - incompatible with the collection. Invalid state."
@@ -152,9 +151,7 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, p
         self._riterator = None
 
     def __del__(self):
-        """
-        Collection - Destructor
-        """
+        """Collection - Destructor."""
         super().__del__()
 
     def iterator(self):
@@ -879,9 +876,7 @@ def clear(self):
         self._ncount = 0
 
     def cstruct(self):
-        """
-        'cstruct' returns the ctypes mapping of the particle data. This depends on the specific structure in question.
-        """
+        """Returns the ctypes mapping of the particle data. This depends on the specific structure in question."""
         cstruct = self._data_c.ctypes.data_as(c_void_p)
         return cstruct
 
@@ -959,6 +954,7 @@ class ParticleAccessorAOS(BaseParticleAccessor):
                   particle is stored in the corresponding data arrays
                   of the ParticleCollecion.
     """
+
     _index = 0
     _next_dt = None
 
diff --git a/parcels/collection/collections.py b/parcels/collection/collections.py
index 68ef700ae..eef004af9 100644
--- a/parcels/collection/collections.py
+++ b/parcels/collection/collections.py
@@ -20,7 +20,8 @@ class Collection(ABC):
     @abstractmethod
     def __init__(self):
         """
-        Collection - Constructor
+        Collection - Constructor.
+
         Initializes a collection by pre-allocating memory (where needed), initialising indexing structures
         (where needed), initialising iterators and preparing the C-JIT-glue.
         """
@@ -28,9 +29,7 @@ def __init__(self):
 
     @abstractmethod
     def __del__(self):
-        """
-        Collection - Destructor
-        """
+        """Collection - Destructor."""
         pass
 
     @property
@@ -749,9 +748,7 @@ def toArray(self):
         pass
 
     def __len__(self):
-        """
-        This function returns the length, in terms of 'number of elements, of a collection.
-        """
+        """This function returns the length, in terms of 'number of elements, of a collection."""
         return self._ncount
 
     @abstractmethod
@@ -790,8 +787,8 @@ class ParticleCollection(Collection):
     _data = None  # formerly: particle_data
 
     def __init__(self):
-        """
-        ParticleCollection - Constructor
+        """Constructor for ParticleCollection.
+
         Initializes a particle collection by pre-allocating memory (where needed), initialising indexing structures
         (where needed), initialising iterators (if maintaining a persistent iterator) and preparing the C-JIT-glue.
 
@@ -808,9 +805,7 @@ def __init__(self):
         super().__init__()
 
     def __del__(self):
-        """
-        ParticleCollection - Destructor
-        """
+        """Destructor for ParticleCollection."""
         pass
 
     @property
@@ -834,9 +829,7 @@ def pu_centers(self):
 
     @property
     def pclass(self):
-        """
-        'pclass' stores the actual class type of the particles allocated and managed in this collection
-        """
+        """Stores the actual class type of the particles allocated and managed in this collection."""
         return self._pclass
 
     @property
@@ -877,9 +870,7 @@ def particle_data(self):
 
     @abstractmethod
     def cstruct(self):
-        """
-        'cstruct' returns the ctypes mapping of the particle data. This depends on the specific structure in question.
-        """
+        """Returns the ctypes mapping of the particle data. This depends on the specific structure in question."""
         pass
 
     @abstractmethod
diff --git a/parcels/collection/collectionsoa.py b/parcels/collection/collectionsoa.py
index 733f2b948..65a234015 100644
--- a/parcels/collection/collectionsoa.py
+++ b/parcels/collection/collectionsoa.py
@@ -34,7 +34,6 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, p
         :param ngrid: number of grids in the fieldset of the overarching ParticleSet - required for initialising the
         field references of the ctypes-link of particles that are allocated
         """
-
         super(ParticleCollection, self).__init__()
 
         assert pid_orig is not None, "particle IDs are None - incompatible with the collection. Invalid state."
@@ -166,9 +165,7 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, p
         self._riterator = None
 
     def __del__(self):
-        """
-        Collection - Destructor
-        """
+        """Collection - Destructor"""
         super().__del__()
 
     def iterator(self):
@@ -781,9 +778,7 @@ def clear(self):
         raise NotImplementedError
 
     def cstruct(self):
-        """
-        'cstruct' returns the ctypes mapping of the particle data. This depends on the specific structure in question.
-        """
+        """Returns the ctypes mapping of the particle data. This depends on the specific structure in question."""
         class CParticles(Structure):
             _fields_ = [(v.name, POINTER(np.ctypeslib.as_ctypes_type(v.dtype))) for v in self._ptype.variables]
 
@@ -861,6 +856,7 @@ class ParticleAccessorSOA(BaseParticleAccessor):
                   particle is stored in the corresponding data arrays
                   of the ParticleCollecion.
     """
+
     _index = 0
     _next_dt = None
 
diff --git a/parcels/collection/iterators.py b/parcels/collection/iterators.py
index b8d2b0917..f57b6547d 100644
--- a/parcels/collection/iterators.py
+++ b/parcels/collection/iterators.py
@@ -22,6 +22,7 @@ class BaseParticleCollectionIterator(ABC):
     """Interface for the ParticleCollection iterator. Provides the
     ability to iterate over the particles in the ParticleCollection.
     """
+
     def ___init___(self):
         self._head = None
         self._tail = None
@@ -29,42 +30,33 @@ def ___init___(self):
 
     @abstractmethod
     def __next__(self):
-        """Returns a ParticleAccessor for the next particle in the
-        ParticleSet.
-        """
+        """Returns a ParticleAccessor for the next particle in the ParticleSet."""
         pass
 
     @property
     def head(self):
-        """Returns a ParticleAccessor for the first particle in the
-        ParticleSet.
-        """
+        """Returns a ParticleAccessor for the first particle in the ParticleSet."""
         return self._head
 
     @property
     def tail(self):
-        """Returns a ParticleAccessor for the last particle in the
-        ParticleSet.
-        """
+        """Returns a ParticleAccessor for the last particle in the ParticleSet."""
         return self._tail
 
     @property
     def current(self):
-        """Returns a ParticleAccessor for the particle that the iteration
-        is currently at.
-        """
+        """Returns a ParticleAccessor for the particle that the iteration currently at."""
         return self._current
 
     @abstractmethod
     def __repr__(self):
-        """Represents the current position in the iteration.
-        """
+        """Represents the current position in the iteration."""
         pass
 
 
 class BaseParticleAccessor(ABC):
-    """Interface for the ParticleAccessor. Implements a wrapper around
-    particles to provide easy access."""
+    """Interface for the ParticleAccessor. Implements a wrapper around particles to provide easy access."""
+
     _pcoll = None
 
     def __init__(self, pcoll):
@@ -82,9 +74,7 @@ def delete(self):
         self.state = OperationCode.Delete
 
     def set_state(self, state):
-        """Syntactic sugar for changing the state of the underlying
-        particle.
-        """
+        """Syntactic sugar for changing the state of the underlying particle."""
         self.state = state
 
     def succeeded(self):
diff --git a/parcels/compilation/codecompiler.py b/parcels/compilation/codecompiler.py
index b8e6f7948..4afb7feae 100644
--- a/parcels/compilation/codecompiler.py
+++ b/parcels/compilation/codecompiler.py
@@ -199,7 +199,8 @@ class CCompiler:
 
     :arg cc: C compiler executable (uses environment variable ``CC`` if not provided).
     :arg cppargs: A list of arguments to the C compiler (optional).
-    :arg ldargs: A list of arguments to the linker (optional)."""
+    :arg ldargs: A list of arguments to the linker (optional).
+    """
 
     def __init__(self, cc=None, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=os.getcwd()):
         if cppargs is None:
@@ -241,9 +242,8 @@ def _create_compile_process_(self, cmd, src, log):
 
 
 class CCompiler_SS(CCompiler):
-    """
-    single-stage C-compiler; used for a SINGLE source file
-    """
+    """Single-stage C-compiler; used for a SINGLE source file."""
+
     def __init__(self, cc=None, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=os.getcwd()):
         super().__init__(cc=cc, cppargs=cppargs, ldargs=ldargs, incdirs=incdirs, libdirs=libdirs, libs=libs, tmp_dir=tmp_dir)
 
@@ -270,7 +270,9 @@ class GNUCompiler_SS(CCompiler_SS):
 
     :arg cppargs: A list of arguments to pass to the C compiler
          (optional).
-    :arg ldargs: A list of arguments to pass to the linker (optional)."""
+    :arg ldargs: A list of arguments to pass to the linker (optional).
+    """
+
     def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=os.getcwd()):
         c_params = GNU_parameters(cppargs, ldargs, incdirs, libdirs, libs)
         super().__init__(c_params.compiler, cppargs=c_params.cppargs, ldargs=c_params.ldargs, incdirs=c_params.incdirs, libdirs=c_params.libdirs, libs=c_params.libs, tmp_dir=tmp_dir)
diff --git a/parcels/compilation/codegenerator.py b/parcels/compilation/codegenerator.py
index 83b0ad70f..daef66580 100644
--- a/parcels/compilation/codegenerator.py
+++ b/parcels/compilation/codegenerator.py
@@ -273,7 +273,8 @@ def __getattr__(self, attr):
 class IntrinsicTransformer(ast.NodeTransformer):
     """AST transformer that catches any mention of intrinsic variable
     names, such as 'particle' or 'fieldset', inserts placeholder objects
-    and propagates attribute access."""
+    and propagates attribute access.
+    """
 
     def __init__(self, fieldset=None, ptype=JITParticle):
         self.fieldset = fieldset
@@ -286,14 +287,14 @@ def __init__(self, fieldset=None, ptype=JITParticle):
         self.stmt_stack = []
 
     def get_tmp(self):
-        """Create a new temporary veriable name"""
+        """Create a new temporary veriable name."""
         tmp = "parcels_tmpvar%d" % self._tmp_counter
         self._tmp_counter += 1
         self.tmp_vars += [tmp]
         return tmp
 
     def visit_Name(self, node):
-        """Inject IntrinsicNode objects into the tree according to keyword"""
+        """Inject IntrinsicNode objects into the tree according to keyword."""
         if node.id == 'fieldset' and self.fieldset is not None:
             node = FieldSetNode(self.fieldset, ccode='fset')
         elif node.id == 'particle':
@@ -480,8 +481,7 @@ def visit_Call(self, node):
 
 
 class TupleSplitter(ast.NodeTransformer):
-    """AST transformer that detects and splits Pythonic tuple
-    assignments into multiple statements for conversion to C."""
+    """AST transformer that detects and splits Pythonic tuple assignments into multiple statements for conversion to C."""
 
     def visit_Assign(self, node):
         if isinstance(node.targets[0], ast.Tuple) \
@@ -498,9 +498,10 @@ def visit_Assign(self, node):
 
 
 class AbstractKernelGenerator(ABC, ast.NodeVisitor):
-    """Code generator class that translates simple Python kernel
-    functions into C functions by populating and accessing the `ccode`
-    attriibute on nodes in the Python AST."""
+    """Code generator class that translates simple Python kernel functions into C functions.
+
+    Works by populating and accessing the `ccode` attriibute on nodes in the Python AST.
+    """
 
     # Intrinsic variables that appear as function arguments
     kernel_vars = ['particle', 'fieldset', 'time', 'output_time', 'tol']
@@ -559,9 +560,11 @@ def visit_FunctionDef(self, node):
         pass
 
     def visit_Call(self, node):
-        """Generate C code for simple C-style function calls. Please
-        note that starred and keyword arguments are currently not
-        supported."""
+        """Generate C code for simple C-style function calls.
+
+        Please note that starred and keyword arguments are currently not
+        supported.
+        """
         pointer_args = False
         parcels_customed_Cfunc = False
         if isinstance(node.func, PrintNode):
@@ -624,8 +627,7 @@ def visit_Call(self, node):
                 raise RuntimeError("Error in converting Kernel to C. See https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_parcels_structure.ipynb#3.-Kernel-execution for hints and tips")
 
     def visit_Name(self, node):
-        """Catches any mention of intrinsic variable names, such as
-        'particle' or 'fieldset' and inserts our placeholder objects"""
+        """Catches any mention of intrinsic variable names such as 'particle' or 'fieldset' and inserts our placeholder objects."""
         if node.id == 'True':
             node.id = "1"
         if node.id == 'False':
@@ -811,30 +813,30 @@ def visit_Pass(self, node):
         node.ccode = c.Statement("")
 
     def visit_FieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         self.field_args[node.obj.ccode_name] = node.obj
 
     def visit_SummedFieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         for fld in node.obj:
             self.field_args[fld.ccode_name] = fld
 
     def visit_NestedFieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         for fld in node.obj:
             self.field_args[fld.ccode_name] = fld
 
     def visit_VectorFieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         self.vector_field_args[node.obj.ccode_name] = node.obj
 
     def visit_SummedVectorFieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         for fld in node.obj:
             self.vector_field_args[fld.ccode_name] = fld
 
     def visit_NestedVectorFieldNode(self, node):
-        """Record intrinsic fields used in kernel"""
+        """Record intrinsic fields used in kernel."""
         for fld in node.obj:
             self.vector_field_args[fld.ccode_name] = fld
 
@@ -1205,8 +1207,7 @@ def visit_NestedVectorFieldEvalNode(self, node):
 
 
 class LoopGenerator:
-    """Code generator class that adds type definitions and the outer
-    loop around kernel functions to generate compilable C code."""
+    """Code generator class that adds type definitions and the outer loop around kernel functions to generate compilable C code."""
 
     def __init__(self, fieldset, ptype=None):
         self.fieldset = fieldset
@@ -1353,8 +1354,7 @@ def generate(self, funcname, field_args, const_args, kernel_ast, c_include):
 
 
 class ParticleObjectLoopGenerator:
-    """Code generator class that adds type definitions and the outer
-    loop around kernel functions to generate compilable C code."""
+    """Code generator class that adds type definitions and the outer loop around kernel functions to generate compilable C code."""
 
     def __init__(self, fieldset=None, ptype=None):
         self.fieldset = fieldset
diff --git a/parcels/examples/example_mitgcm.py b/parcels/examples/example_mitgcm.py
index ade100bcb..372dbb8a4 100644
--- a/parcels/examples/example_mitgcm.py
+++ b/parcels/examples/example_mitgcm.py
@@ -17,8 +17,7 @@
 
 
 def run_mitgcm_zonally_reentrant(mode):
-    """Function that shows how to load MITgcm data in a zonally periodic domain"""
-
+    """Function that shows how to load MITgcm data in a zonally periodic domain."""
     data_path = path.join(path.dirname(__file__), "MITgcm_example_data/")
     filenames = {
         "U": data_path + "mitgcm_UV_surface_zonally_reentrant.nc",
diff --git a/parcels/examples/example_moving_eddies.py b/parcels/examples/example_moving_eddies.py
index fd5c629ea..2ad8ec1c6 100644
--- a/parcels/examples/example_moving_eddies.py
+++ b/parcels/examples/example_moving_eddies.py
@@ -91,11 +91,11 @@ def cosd(x):
 
 def moving_eddies_example(fieldset, outfile, npart=2, mode='jit', verbose=False,
                           method=AdvectionRK4):
-    """Configuration of a particle set that follows two moving eddies
+    """Configuration of a particle set that follows two moving eddies.
 
     :arg fieldset: :class FieldSet: that defines the flow field
-    :arg npart: Number of particles to initialise"""
-
+    :arg npart: Number of particles to initialise.
+    """
     # Determine particle class according to mode
     start = (3.3, 46.) if fieldset.U.grid.mesh == 'spherical' else (3.3e5, 1e5)
     finish = (3.3, 47.8) if fieldset.U.grid.mesh == 'spherical' else (3.3e5, 2.8e5)
@@ -163,7 +163,7 @@ def test_moving_eddies_fieldset(mode, mesh, tmpdir):
 
 
 def fieldsetfile(mesh, tmpdir):
-    """Generate fieldset files for moving_eddies test"""
+    """Generate fieldset files for moving_eddies test."""
     filename = tmpdir.join('moving_eddies')
     fieldset = moving_eddies_fieldset(200, 350, mesh=mesh)
     fieldset.write(filename)
diff --git a/parcels/examples/example_nemo_curvilinear.py b/parcels/examples/example_nemo_curvilinear.py
index a138aec92..abd5d7583 100644
--- a/parcels/examples/example_nemo_curvilinear.py
+++ b/parcels/examples/example_nemo_curvilinear.py
@@ -21,7 +21,7 @@
 
 
 def run_nemo_curvilinear(mode, outfile, advtype='RK4'):
-    """Function that shows how to read in curvilinear grids, in this case from NEMO"""
+    """Function that shows how to read in curvilinear grids, in this case from NEMO."""
     data_path = path.join(path.dirname(__file__), 'NemoCurvilinear_data/')
 
     filenames = {'U': {'lon': data_path + 'mesh_mask.nc4',
diff --git a/parcels/examples/example_peninsula.py b/parcels/examples/example_peninsula.py
index 1268a089a..732d6f1bd 100644
--- a/parcels/examples/example_peninsula.py
+++ b/parcels/examples/example_peninsula.py
@@ -23,8 +23,7 @@
 
 
 def peninsula_fieldset(xdim, ydim, mesh='flat', grid_type='A'):
-    """Construct a fieldset encapsulating the flow field around an
-    idealised peninsula.
+    """Construct a fieldset encapsulating the flow field around an idealised peninsula.
 
     :param xdim: Horizontal dimension of the generated fieldset
     :param xdim: Vertical dimension of the generated fieldset
@@ -98,9 +97,9 @@ def peninsula_example(fieldset, outfile, npart, mode='jit', degree=1,
                       verbose=False, output=True, method=AdvectionRK4):
     """Example configuration of particle flow around an idealised Peninsula
 
-    :arg filename: Basename of the input fieldset
-    :arg npart: Number of particles to intialise"""
-
+    :arg filename: Basename of the input fieldset.
+    :arg npart: Number of particles to intialise.
+    """
     # First, we define a custom Particle class to which we add a
     # custom variable, the initial stream function value p.
     # We determine the particle base class according to mode.
@@ -140,7 +139,7 @@ class MyParticle(ptype[mode]):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('mesh', ['flat', 'spherical'])
 def test_peninsula_fieldset(mode, mesh, tmpdir):
-    """Execute peninsula test from fieldset generated in memory"""
+    """Execute peninsula test from fieldset generated in memory."""
     fieldset = peninsula_fieldset(100, 50, mesh)
     outfile = tmpdir.join("Peninsula")
     pset = peninsula_example(fieldset, outfile, 5, mode=mode, degree=1)
@@ -155,7 +154,7 @@ def test_peninsula_fieldset(mode, mesh, tmpdir):
 @pytest.mark.parametrize('mode', ['scipy'])  # Analytical Advection only implemented in Scipy mode
 @pytest.mark.parametrize('mesh', ['flat', 'spherical'])
 def test_peninsula_fieldset_AnalyticalAdvection(mode, mesh, tmpdir):
-    """Execute peninsula test using Analytical Advection on C grid"""
+    """Execute peninsula test using Analytical Advection on C grid."""
     fieldset = peninsula_fieldset(101, 51, 'flat', grid_type='C')
     outfile = tmpdir.join("PeninsulaAA")
     pset = peninsula_example(fieldset, outfile, npart=10, mode=mode,
@@ -166,7 +165,7 @@ def test_peninsula_fieldset_AnalyticalAdvection(mode, mesh, tmpdir):
 
 
 def fieldsetfile(mesh, tmpdir):
-    """Generate fieldset files for peninsula test"""
+    """Generate fieldset files for peninsula test."""
     filename = tmpdir.join('peninsula')
     fieldset = peninsula_fieldset(100, 50, mesh=mesh)
     fieldset.write(filename)
@@ -176,7 +175,7 @@ def fieldsetfile(mesh, tmpdir):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('mesh', ['flat', 'spherical'])
 def test_peninsula_file(mode, mesh, tmpdir):
-    """Open fieldset files and execute"""
+    """Open fieldset files and execute."""
     gc.collect()
     fieldset = FieldSet.from_parcels(fieldsetfile(mesh, tmpdir), extra_fields={'P': 'P'}, allow_time_extrapolation=True)
     outfile = tmpdir.join("Peninsula")
diff --git a/parcels/examples/example_recursive_errorhandling.py b/parcels/examples/example_recursive_errorhandling.py
index 5c3e2f2ca..36f72eb4d 100644
--- a/parcels/examples/example_recursive_errorhandling.py
+++ b/parcels/examples/example_recursive_errorhandling.py
@@ -27,8 +27,8 @@ def test_recursive_errorhandling(mode, xdim=2, ydim=2):
 
     This behaviour can be useful if particles need to be 'pushed out'
     from e.g. land. Note however that current under-the-hood
-    implementation is not extremely efficient, so code could be slow."""
-
+    implementation is not extremely efficient, so code could be slow.
+    """
     dimensions = {'lon': np.linspace(0., 1., xdim, dtype=np.float32),
                   'lat': np.linspace(0., 1., ydim, dtype=np.float32)}
     data = {'U': np.zeros((ydim, xdim), dtype=np.float32),
@@ -44,13 +44,17 @@ def test_recursive_errorhandling(mode, xdim=2, ydim=2):
 
     def TestLon(particle, fieldset, time):
         """Kernel to check whether a longitude is larger than fieldset.minlon.
-        If not, the Kernel throws an error"""
+
+        If not, the Kernel throws an error.
+        """
         if particle.lon <= fieldset.minlon:
             return ErrorCode.Error
 
     def Error_RandomiseLon(particle, fieldset, time):
         """Error handling kernel that draws a new longitude.
-        Note that this new longitude can be smaller than fieldset.minlon"""
+
+        Note that this new longitude can be smaller than fieldset.minlon.
+        """
         particle.lon = ParcelsRandom.uniform(0., 1.)
 
     ParcelsRandom.seed(123456)
diff --git a/parcels/field.py b/parcels/field.py
index f0e991e30..88f074d52 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -93,6 +93,7 @@ class Field:
 
     * `Summed Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_SummedFields.ipynb>`_
     """
+
     def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=None, mesh='flat', timestamps=None,
                  fieldtype=None, transpose=False, vmin=None, vmax=None, cast_data_dtype='float32', time_origin=None,
                  interp_method='linear', allow_time_extrapolation=None, time_periodic=False, gridindexingtype='nemo',
@@ -267,7 +268,7 @@ def collect_timeslices(timestamps, data_filenames, _grid_fb_class, dimensions, i
     def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
                     mesh='spherical', timestamps=None, allow_time_extrapolation=None, time_periodic=False,
                     deferred_load=True, **kwargs):
-        """Create field from netCDF file
+        """Create field from netCDF file.
 
         :param filenames: list of filenames to read for the field. filenames can be a list [files] or
                a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data)
@@ -486,7 +487,7 @@ def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
     @classmethod
     def from_xarray(cls, da, name, dimensions, mesh='spherical', allow_time_extrapolation=None,
                     time_periodic=False, **kwargs):
-        """Create field from xarray Variable
+        """Create field from xarray Variable.
 
         :param da: Xarray DataArray
         :param name: Name of the Field
@@ -503,7 +504,6 @@ def from_xarray(cls, da, name, dimensions, mesh='spherical', allow_time_extrapol
         :param time_periodic: boolean whether to loop periodically over the time component of the FieldSet
                This flag overrides the allow_time_interpolation and sets it to False
         """
-
         data = da.data
         interp_method = kwargs.pop('interp_method', 'linear')
 
@@ -586,19 +586,20 @@ def set_scaling_factor(self, factor):
             self.data *= factor
 
     def set_depth_from_field(self, field):
-        """Define the depth dimensions from another (time-varying) field
+        """Define the depth dimensions from another (time-varying) field.
 
         See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timevaryingdepthdimensions.ipynb>`_
-        for a detailed explanation on how to set up time-evolving depth dimensions
-
+        for a detailed explanation on how to set up time-evolving depth dimensions.
         """
         self.grid.depth_field = field
         if self.grid != field.grid:
             field.grid.depth_field = field
 
     def calc_cell_edge_sizes(self):
-        """Method to calculate cell sizes based on numpy.gradient method
-                Currently only works for Rectilinear Grids"""
+        """Method to calculate cell sizes based on numpy.gradient method.
+
+        Currently only works for Rectilinear Grids
+        """
         if not self.grid.cell_edge_sizes:
             if self.grid.gtype in (GridCode.RectilinearZGrid, GridCode.RectilinearSGrid):
                 self.grid.cell_edge_sizes['x'] = np.zeros((self.grid.ydim, self.grid.xdim), dtype=np.float32)
@@ -618,8 +619,10 @@ def calc_cell_edge_sizes(self):
                 exit(-1)
 
     def cell_areas(self):
-        """Method to calculate cell sizes based on cell_edge_sizes
-                Currently only works for Rectilinear Grids"""
+        """Method to calculate cell sizes based on cell_edge_sizes.
+
+        Currently only works for Rectilinear Grids
+        """
         if not self.grid.cell_edge_sizes:
             self.calc_cell_edge_sizes()
         return self.grid.cell_edge_sizes['x'] * self.grid.cell_edge_sizes['y']
@@ -1038,13 +1041,13 @@ def interpolator3D(self, ti, z, y, x, time, particle=None):
             raise RuntimeError(self.interp_method+" is not implemented for 3D grids")
 
     def temporal_interpolate_fullfield(self, ti, time):
-        """Calculate the data of a field between two snapshots,
-        using linear interpolation
+        """Calculate the data of a field between two snapshots using linear interpolation.
 
         :param ti: Index in time array associated with time (via :func:`time_index`)
         :param time: Time to interpolate to
 
-        :rtype: Linearly interpolated field"""
+        :rtype: Linearly interpolated field.
+        """
         t0 = self.grid.time[ti]
         if time == t0:
             return self.data[ti, :]
@@ -1057,8 +1060,7 @@ def temporal_interpolate_fullfield(self, ti, time):
             return f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
 
     def spatial_interpolation(self, ti, z, y, x, time, particle=None):
-        """Interpolate horizontal field values using a SciPy interpolator"""
-
+        """Interpolate horizontal field values using a SciPy interpolator."""
         if self.grid.zdim == 1:
             val = self.interpolator2D(ti, z, y, x, particle=particle)
         else:
@@ -1072,7 +1074,7 @@ def spatial_interpolation(self, ti, z, y, x, time, particle=None):
             return val
 
     def time_index(self, time):
-        """Find the index in the time array associated with a given time
+        """Find the index in the time array associated with a given time.
 
         Note that we normalize to either the first or the last index
         if the sampled value is outside the time value range.
@@ -1214,9 +1216,7 @@ def chunk_data(self):
 
     @property
     def ctypes_struct(self):
-        """Returns a ctypes struct object containing all relevant
-        pointers and sizes for this field."""
-
+        """Returns a ctypes struct object containing all relevant pointers and sizes for this field."""
         # Ctypes struct corresponding to the type definition in parcels.h
         class CField(Structure):
             _fields_ = [('xdim', c_int), ('ydim', c_int), ('zdim', c_int),
@@ -1247,7 +1247,7 @@ class CField(Structure):
 
     def show(self, animation=False, show_time=None, domain=None, depth_level=0, projection='PlateCarree', land=True,
              vmin=None, vmax=None, savefile=None, **kwargs):
-        """Method to 'show' a Parcels Field
+        """Method to 'show' a Parcels Field.
 
         :param animation: Boolean whether result is a single plot, or an animation
         :param show_time: Time in seconds from start after which to show the Field (only in single-plot mode)
@@ -1266,7 +1266,9 @@ def show(self, animation=False, show_time=None, domain=None, depth_level=0, proj
             plt.show()
 
     def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
-        """Add a 'halo' to all Fields in a FieldSet, through extending the Field (and lon/lat)
+        """Add a 'halo' to all Fields in a FieldSet.
+
+        Add a 'halo' to all Fields in a FieldSet, through extending the Field (and lon/lat)
         by copying a small portion of the field on one side of the domain to the other.
         Before adding a periodic halo to the Field, it has to be added to the Grid on which the Field depends
 
@@ -1310,10 +1312,11 @@ def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
             return data
 
     def write(self, filename, varname=None):
-        """Write a :class:`Field` to a netcdf file
+        """Write a :class:`Field` to a netcdf file.
 
         :param filename: Basename of the file
-        :param varname: Name of the field, to be appended to the filename"""
+        :param varname: Name of the field, to be appended to the filename.
+        """
         filepath = str(Path(f'{filename}{self.name}.nc'))
         if varname is None:
             varname = self.name
@@ -1429,6 +1432,7 @@ class VectorField:
     :param V: field defining the meridional component
     :param W: field defining the vertical component (default: None)
     """
+
     def __init__(self, name, U, V, W=None):
         self.name = name
         self.U = U
@@ -1641,6 +1645,8 @@ def spatial_c_grid_interpolation3D_full(self, ti, z, y, x, time, particle=None):
 
     def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, applyConversion=True):
         """
+        Perform C grid interpolation in 3D.
+
         +---+---+---+
         |   |V1 |   |
         +---+---+---+
@@ -1824,7 +1830,8 @@ def ccode_eval_object(self, varU, varV, varW, U, V, W, t, z, y, x):
 
 
 class DeferredArray():
-    """Class used for throwing error when Field.data is not read in deferred loading mode"""
+    """Class used for throwing error when Field.data is not read in deferred loading mode."""
+
     data_shape = ()
 
     def __init__(self):
@@ -1847,8 +1854,9 @@ def __getitem__(self, key):
 
 
 class SummedField(list):
-    """Class SummedField is a list of Fields over which Field interpolation
-    is summed. This can e.g. be used when combining multiple flow fields,
+    """Class SummedField is a list of Fields over which Field interpolation is summed.
+
+    This can e.g. be used when combining multiple flow fields,
     where the total flow is the sum of all the individual flows.
     Note that the individual Fields can be on different Grids.
     Also note that, since SummedFields are lists, the individual Fields can
@@ -1917,8 +1925,10 @@ def __add__(self, field):
 
 
 class NestedField(list):
-    """Class NestedField is a list of Fields from which the first one to be not declared out-of-boundaries
-    at particle position is interpolated. This induces that the order of the fields in the list matters.
+    """NestedField is a class that allows for interpolation of fields on different grids of potentially varying resolution.
+
+    The NestedField class is a list of Fields where the first Field that contains the particle within the domain is then used for interpolation.
+    This induces that the order of the fields in the list matters.
     Each one it its turn, a field is interpolated: if the interpolation succeeds or if an error other
     than `ErrorOutOfBounds` is thrown, the function is stopped. Otherwise, next field is interpolated.
     NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries.
diff --git a/parcels/fieldfilebuffer.py b/parcels/fieldfilebuffer.py
index 923214cc1..85e8339d0 100644
--- a/parcels/fieldfilebuffer.py
+++ b/parcels/fieldfilebuffer.py
@@ -299,7 +299,8 @@ def __enter__(self):
         return self
 
     def __exit__(self, type, value, traceback):
-        """
+        """Function releases the file handle.
+
         This function releases the file handle. Hence access to the dataset and its header-information is lost. The
         previously executed chunking is lost. Furthermore, if the file access required file locking, the lock-handle
         is freed so other processes can now access the file again.
@@ -307,7 +308,8 @@ def __exit__(self, type, value, traceback):
         self.close()
 
     def close(self):
-        """
+        """Teardown FileBuffer object with dask.
+
         This function can be called to initialise an orderly teardown of a FileBuffer object with dask, meaning
         to release the file handle, deposing the dataset, and releasing the file lock (if required).
         """
@@ -357,9 +359,8 @@ def add_to_dimension_name_map(self, name_map):
             self._static_name_maps[pcls_dim_name].append(name_map[pcls_dim_name])
 
     def _get_available_dims_indices_by_request(self):
-        """
-        [private function - not to be called from outside the class]
-        Returns a dict mapping 'parcels_dimname' -> [None, int32_index_data_array].
+        """Returns a dict mapping 'parcels_dimname' -> [None, int32_index_data_array].
+
         This dictionary is based on the information provided by the requested dimensions.
         Example: {'time': 0, 'depth': None, 'lat': 1, 'lon': 2}
         """
@@ -382,7 +383,6 @@ def _get_available_dims_indices_by_request(self):
 
     def _get_available_dims_indices_by_namemap(self):
         """
-        [private function - not to be called from outside the class]
         Returns a dict mapping 'parcels_dimname' -> [None, int32_index_data_array].
         This dictionary is based on the information provided by the requested dimensions.
         Example: {'time': 0, 'depth': 1, 'lat': 2, 'lon': 3}
@@ -394,7 +394,6 @@ def _get_available_dims_indices_by_namemap(self):
 
     def _get_available_dims_indices_by_netcdf_file(self):
         """
-        [private function - not to be called from outside the class]
         [File needs to be open (i.e. self.dataset is not None) for this to work - otherwise generating an error]
         Returns a dict mapping 'parcels_dimname' -> [None, int32_index_data_array].
         This dictionary is based on the information provided by the requested dimensions.
@@ -419,7 +418,6 @@ def _get_available_dims_indices_by_netcdf_file(self):
 
     def _is_dimension_available(self, dimension_name):
         """
-        [private function - not to be called from outside the class]
         This function returns a boolean value indicating if a certain variable (name) is avaialble in the
         requested dimensions as well as in the actual dataset of the file. If any of the two conditions is not met,
         if returns 'False'.
@@ -430,7 +428,6 @@ def _is_dimension_available(self, dimension_name):
 
     def _is_dimension_chunked(self, dimension_name):
         """
-        [private function - not to be called from outside the class]
         This functions returns a boolean value indicating if a certain variable is available in the requested
         dimensions, the NetCDF file dataset, and is also required to be chunked according to the requested
         chunksize dictionary. If any of the two conditions is not met, if returns 'False'.
@@ -444,7 +441,6 @@ def _is_dimension_chunked(self, dimension_name):
 
     def _is_dimension_in_dataset(self, parcels_dimension_name, netcdf_dimension_name=None):
         """
-        [private function - not to be called from outside the class]
         [File needs to be open (i.e. self.dataset is not None) for this to work - otherwise generating an error]
         This function returns the index, the name and the size of a NetCDF dimension in the file (in order: index, name, size).
         It requires as input the name of the related parcels dimension (i.e. one of ['time', 'depth', 'lat', 'lon']. If
@@ -472,7 +468,6 @@ def _is_dimension_in_dataset(self, parcels_dimension_name, netcdf_dimension_name
 
     def _is_dimension_in_chunksize_request(self, parcels_dimension_name):
         """
-        [private function - not to be called from outside the class]
         This function returns the dense-array index, the NetCDF dimension name and the requested chunsize of a requested
         parcels dimension(in order: index, name, size). This only works if the chunksize is provided as a dictionary
         of tuples of parcels dimensions and their chunk mapping (i.e. dict(parcels_dim_name => (netcdf_dim_name, chunksize)).
@@ -490,16 +485,12 @@ def _is_dimension_in_chunksize_request(self, parcels_dimension_name):
         return k, dname, dvalue
 
     def _netcdf_DimNotFound_warning_message(self, dimension_name):
-        """
-        [private function - not to be called from outside the class]
-        Helper function that issues a warning message if a certain requested NetCDF dimension is not found in the file.
-        """
+        """Helper function that issues a warning message if a certain requested NetCDF dimension is not found in the file."""
         display_name = dimension_name if (dimension_name not in self.dimensions) else self.dimensions[dimension_name]
         return f"Did not find {display_name} in NetCDF dims. Please specifiy chunksize as dictionary for NetCDF dimension names, e.g.\n chunksize={{ '{display_name}': <number>, ... }}."
 
     def _chunkmap_to_chunksize(self):
         """
-        [private function - not to be called from outside the class]
         [File needs to be open via the '__enter__'-method for this to work - otherwise generating an error]
         This functions translates the array-index-to-chunksize chunk map into a proper fieldsize dictionary that
         can later be used for re-qunking, if a previously-opened file is re-opened again.
@@ -539,7 +530,6 @@ def _chunkmap_to_chunksize(self):
 
     def _get_initial_chunk_dictionary_by_dict_(self):
         """
-        [private function - not to be called from outside the class]
         [File needs to be open (i.e. self.dataset is not None) for this to work - otherwise generating an error]
         Maps and correlates the requested dictionary-style chunksize with the requested parcels dimensions, variables
         and the NetCDF-available dimensions. Thus, it takes care to remove chunksize arguments that are not in the
@@ -590,10 +580,7 @@ def _get_initial_chunk_dictionary_by_dict_(self):
         return chunk_dict, chunk_index_map
 
     def _failsafe_parse_(self):
-        """
-        [private function - not to be called from outside the class]
-        ['name' need to be initialised]
-        """
+        """['name' need to be initialised]"""
         # ==== fail - open it as a normal array and deduce the dimensions from the variable-function names ==== #
         # ==== done by parsing ALL variables in the NetCDF, and comparing their call-parameters with the   ==== #
         # ==== name map available here.                                                                    ==== #
@@ -638,7 +625,6 @@ def _failsafe_parse_(self):
 
     def _get_initial_chunk_dictionary(self):
         """
-        [private function - not to be called from outside the class]
         Super-function that maps and correlates the requested chunksize with the requested parcels dimensions, variables
         and the NetCDF-available dimensions. Thus, it takes care to remove chunksize arguments that are not in the
         Parcels- or NetCDF dimensions, or whose chunking would be omitted due to an empty chunk dimension.
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index f83e9685c..c0aea9071 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -23,12 +23,13 @@
 
 
 class FieldSet:
-    """FieldSet class that holds hydrodynamic data needed to execute particles
+    """FieldSet class that holds hydrodynamic data needed to execute particles.
 
     :param U: :class:`parcels.field.Field` object for zonal velocity component
     :param V: :class:`parcels.field.Field` object for meridional velocity component
     :param fields: Dictionary of additional :class:`parcels.field.Field` objects
     """
+
     def __init__(self, U, V, fields=None):
         self.gridset = GridSet()
         self.completed = False
@@ -55,7 +56,7 @@ def checkvaliddimensionsdict(dims):
     @classmethod
     def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
                   allow_time_extrapolation=None, time_periodic=False, **kwargs):
-        """Initialise FieldSet object from raw data
+        """Initialise FieldSet object from raw data.
 
         :param data: Dictionary mapping field names to numpy arrays.
                Note that at least a 'U' and 'V' numpy array need to be given, and that
@@ -96,7 +97,6 @@ def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
         * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_
 
         """
-
         fields = {}
         for name, datafld in data.items():
             # Use dimensions[name] if dimensions is a dict of dicts
@@ -127,7 +127,7 @@ def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
         return cls(u, v, fields=fields)
 
     def add_field(self, field, name=None):
-        """Add a :class:`parcels.field.Field` object to the FieldSet
+        """Add a :class:`parcels.field.Field` object to the FieldSet.
 
         :param field: :class:`parcels.field.Field` object to be added
         :param name: Name of the :class:`parcels.field.Field` object to be added
@@ -174,7 +174,7 @@ def add_constant_field(self, name, value, mesh='flat'):
         self.add_field(Field(name, value, lon=0, lat=0, mesh=mesh))
 
     def add_vector_field(self, vfield):
-        """Add a :class:`parcels.field.VectorField` object to the FieldSet
+        """Add a :class:`parcels.field.VectorField` object to the FieldSet.
 
         :param vfield: :class:`parcels.field.VectorField` object to be added
         """
@@ -291,7 +291,7 @@ def parse_wildcards(cls, paths, filenames, var):
     def from_netcdf(cls, filenames, variables, dimensions, indices=None, fieldtype=None,
                     mesh='spherical', timestamps=None, allow_time_extrapolation=None, time_periodic=False,
                     deferred_load=True, chunksize=None, **kwargs):
-        """Initialises FieldSet object from NetCDF files
+        """Initialises FieldSet object from NetCDF files.
 
         :param filenames: Dictionary mapping variables to file(s). The
                filepath may contain wildcards to indicate multiple files
@@ -484,7 +484,6 @@ def from_nemo(cls, filenames, variables, dimensions, indices=None, mesh='spheric
         :param chunksize: size of the chunks in dask loading. Default is None (no chunking)
 
         """
-
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_nemo'
         if kwargs.pop('gridindexingtype', 'nemo') != 'nemo':
@@ -500,20 +499,21 @@ def from_nemo(cls, filenames, variables, dimensions, indices=None, mesh='spheric
     def from_mitgcm(cls, filenames, variables, dimensions, indices=None, mesh='spherical',
                     allow_time_extrapolation=None, time_periodic=False,
                     tracer_interp_method='cgrid_tracer', chunksize=None, **kwargs):
-        """Initialises FieldSet object from NetCDF files of MITgcm fields.
-           All parameters and keywords are exactly the same as for FieldSet.from_nemo(), except that
-           gridindexing is set to 'mitgcm' for grids that have the shape
-
-           +-----------------------------+-----------------------------+-----------------------------+
-           |                             |         V[k,j+1,i]          |                             |
-           +-----------------------------+-----------------------------+-----------------------------+
-           |U[k,j,i]                     |    W[k-1:k,j,i], T[k,j,i]   |U[k,j,i+1]                   |
-           +-----------------------------+-----------------------------+-----------------------------+
-           |                             |         V[k,j,i]            +                             |
-           +-----------------------------+-----------------------------+-----------------------------+
-
-           For indexing details: https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#spatial-discretization-of-the-dynamical-equations
-           Note that vertical velocity (W) is assumed postive in the positive z direction (which is upward in MITgcm)
+        """
+        Initialises FieldSet object from NetCDF files of MITgcm fields.
+        All parameters and keywords are exactly the same as for FieldSet.from_nemo(), except that
+        gridindexing is set to 'mitgcm' for grids that have the shape
+
+        +-----------------------------+-----------------------------+-----------------------------+
+        |                             |         V[k,j+1,i]          |                             |
+        +-----------------------------+-----------------------------+-----------------------------+
+        |U[k,j,i]                     |    W[k-1:k,j,i], T[k,j,i]   |U[k,j,i+1]                   |
+        +-----------------------------+-----------------------------+-----------------------------+
+        |                             |         V[k,j,i]            +                             |
+        +-----------------------------+-----------------------------+-----------------------------+
+
+        For indexing details: https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#spatial-discretization-of-the-dynamical-equations
+        Note that vertical velocity (W) is assumed postive in the positive z direction (which is upward in MITgcm)
         """
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_mitgcm'
@@ -582,10 +582,9 @@ def from_c_grid_dataset(cls, filenames, variables, dimensions, indices=None, mes
                Note that in the case of from_nemo() and from_cgrid(), the velocity fields are default to 'cgrid_velocity'
         :param gridindexingtype: The type of gridindexing. Set to 'nemo' in FieldSet.from_nemo()
                See also the Grid indexing documentation on oceanparcels.org
-        :param chunksize: size of the chunks in dask loading
+        :param chunksize: size of the chunks in dask loading.
 
         """
-
         if 'U' in dimensions and 'V' in dimensions and dimensions['U'] != dimensions['V']:
             raise ValueError("On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. "
                              "See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb")
@@ -672,7 +671,6 @@ def from_pop(cls, filenames, variables, dimensions, indices=None, mesh='spherica
                but many POP outputs are in 'cm'
 
         """
-
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_pop'
         fieldset = cls.from_b_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
@@ -752,7 +750,6 @@ def from_mom5(cls, filenames, variables, dimensions, indices=None, mesh='spheric
 
 
         """
-
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_mom5'
         fieldset = cls.from_b_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
@@ -820,7 +817,6 @@ def from_b_grid_dataset(cls, filenames, variables, dimensions, indices=None, mes
         :param chunksize: size of the chunks in dask loading
 
         """
-
         if 'U' in dimensions and 'V' in dimensions and dimensions['U'] != dimensions['V']:
             raise ValueError("On a B-grid, the dimensions of velocities should be the (top) corners of the grid cells, so the same for U and V. "
                              "See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb")
@@ -870,7 +866,6 @@ def from_parcels(cls, basename, uvar='vozocrtx', vvar='vomecrty', indices=None,
         :param chunksize: size of the chunks in dask loading
 
         """
-
         if extra_fields is None:
             extra_fields = {}
         if 'creation_log' not in kwargs.keys():
@@ -916,7 +911,6 @@ def from_xarray_dataset(cls, ds, variables, dimensions, mesh='spherical', allow_
         :param time_periodic: To loop periodically over the time component of the Field. It is set to either False or the length of the period (either float in seconds or datetime.timedelta object). (Default: False)
                This flag overrides the allow_time_interpolation and sets it to False
         """
-
         fields = {}
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_xarray_dataset'
@@ -937,7 +931,8 @@ def from_xarray_dataset(cls, ds, variables, dimensions, mesh='spherical', allow_
 
     def get_fields(self):
         """Returns a list of all the :class:`parcels.field.Field` and :class:`parcels.field.VectorField`
-        objects associated with this FieldSet"""
+        objects associated with this FieldSet.
+        """
         fields = []
         for v in self.__dict__.values():
             if type(v) in [Field, VectorField]:
@@ -975,7 +970,6 @@ def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
         :param meridional: Create a halo in meridional direction (boolean)
         :param halosize: size of the halo (in grid points). Default is 5 grid points
         """
-
         for grid in self.gridset.grids:
             grid.add_periodic_halo(zonal, meridional, halosize)
         for attr, value in iter(self.__dict__.items()):
@@ -983,10 +977,10 @@ def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
                 value.add_periodic_halo(zonal, meridional, halosize)
 
     def write(self, filename):
-        """Write FieldSet to NetCDF file using NEMO convention
-
-        :param filename: Basename of the output fileset"""
+        """Write FieldSet to NetCDF file using NEMO convention.
 
+        :param filename: Basename of the output fileset.
+        """
         if MPI is None or MPI.COMM_WORLD.Get_rank() == 0:
             logger.info(f"Generating FieldSet output with basename: {filename}")
 
diff --git a/parcels/grid.py b/parcels/grid.py
index 3a6849418..da848ee8b 100644
--- a/parcels/grid.py
+++ b/parcels/grid.py
@@ -23,9 +23,7 @@ class CGrid(Structure):
 
 
 class Grid:
-    """Grid class that defines a (spatial and temporal) grid on which Fields are defined
-
-    """
+    """Grid class that defines a (spatial and temporal) grid on which Fields are defined."""
 
     def __init__(self, lon, lat, time, time_origin, mesh):
         self.xi = None
@@ -96,7 +94,8 @@ def ctypes_struct(self):
     @property
     def child_ctypes_struct(self):
         """Returns a ctypes struct object containing all relevant
-        pointers and sizes for this grid."""
+        pointers and sizes for this grid.
+        """
 
         class CStructuredGrid(Structure):
             # z4d is only to have same cstruct as RectilinearSGrid
@@ -240,8 +239,9 @@ def chunk_loaded(self):
 
 
 class RectilinearGrid(Grid):
-    """Rectilinear Grid
-       Private base class for RectilinearZGrid and RectilinearSGrid
+    """Rectilinear Grid class
+
+    Private base class for RectilinearZGrid and RectilinearSGrid
 
     """
 
@@ -293,7 +293,7 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
 
 
 class RectilinearZGrid(RectilinearGrid):
-    """Rectilinear Z Grid
+    """Rectilinear Z Grid.
 
     :param lon: Vector containing the longitude coordinates of the grid
     :param lat: Vector containing the latitude coordinates of the grid
diff --git a/parcels/gridset.py b/parcels/gridset.py
index 5ab5f36f8..46a431fa6 100644
--- a/parcels/gridset.py
+++ b/parcels/gridset.py
@@ -4,9 +4,7 @@
 
 
 class GridSet:
-    """GridSet class that holds the Grids on which the Fields are defined
-
-    """
+    """GridSet class that holds the Grids on which the Fields are defined."""
 
     def __init__(self):
         self.grids = []
@@ -47,10 +45,10 @@ def add_grid(self, field):
 
     def dimrange(self, dim):
         """Returns maximum value of a dimension (lon, lat, depth or time)
-           on 'left' side and minimum value on 'right' side for all grids
-           in a gridset. Useful for finding e.g. longitude range that
-           overlaps on all grids in a gridset"""
-
+        on 'left' side and minimum value on 'right' side for all grids
+        in a gridset. Useful for finding e.g. longitude range that
+        overlaps on all grids in a gridset.
+        """
         maxleft, minright = (-np.inf, np.inf)
         for g in self.grids:
             if getattr(g, dim).size == 1:
diff --git a/parcels/interaction/baseinteractionkernel.py b/parcels/interaction/baseinteractionkernel.py
index 6390aedb2..057bb7060 100644
--- a/parcels/interaction/baseinteractionkernel.py
+++ b/parcels/interaction/baseinteractionkernel.py
@@ -64,7 +64,7 @@ def check_fieldsets_in_kernels(self, pyfunc):
 
     def check_kernel_signature_on_version(self):
         """
-        returns numkernelargs
+        Returns numkernelargs.
         Adaptation of this method in the BaseKernel that works with
         lists of functions.
         """
diff --git a/parcels/interaction/interactionkernelsoa.py b/parcels/interaction/interactionkernelsoa.py
index c3683bdc4..40d6ea345 100644
--- a/parcels/interaction/interactionkernelsoa.py
+++ b/parcels/interaction/interactionkernelsoa.py
@@ -69,7 +69,7 @@ def __radd__(self, kernel):
         return kernel.merge(self, InteractionKernelSOA)
 
     def execute_python(self, pset, endtime, dt):
-        """Performs the core update loop via Python
+        """Performs the core update loop via Python.
 
         InteractionKernels do not implement ways to catch or recover from
         errors caused during execution of the kernel function(s).
@@ -122,7 +122,7 @@ def execute_python(self, pset, endtime, dt):
                 pset[particle_idx].dt = dt
 
     def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False):
-        """Execute this Kernel over a ParticleSet for several timesteps
+        """Execute this Kernel over a ParticleSet for several timesteps.
 
         InteractionKernels do not implement ways to catch or recover from
         errors caused during execution of the kernel function(s).
diff --git a/parcels/interaction/neighborsearch/base.py b/parcels/interaction/neighborsearch/base.py
index 206d42df2..4c59cda80 100644
--- a/parcels/interaction/neighborsearch/base.py
+++ b/parcels/interaction/neighborsearch/base.py
@@ -41,26 +41,26 @@ def __init__(self, inter_dist_vert, inter_dist_horiz,
 
     @abstractmethod
     def find_neighbors_by_coor(self, coor):
-        '''Get the neighbors around a certain location.
+        """Get the neighbors around a certain location.
 
         :param coor: Numpy array with [depth, lat, lon].
         :returns List of particle indices.
-        '''
+        """
         raise NotImplementedError
 
     def find_neighbors_by_idx(self, particle_idx):
-        '''Get the neighbors around a certain particle.
+        """Get the neighbors around a certain particle.
 
         Mainly useful for Structure of Array (SoA) datastructure
 
         :param particle_idx: index of the particle (SoA).
         :returns List of particle indices
-        '''
+        """
         coor = self._values[:, particle_idx].reshape(3, 1)
         return self.find_neighbors_by_coor(coor)
 
     def update_values(self, new_values, new_active_mask=None):
-        '''Update the coordinates of the particles.
+        """Update the coordinates of the particles.
 
         This is a default implementation simply rebuilds the structure.
         If the rebuilding is slow, a faster implementation can be provided.
@@ -68,7 +68,7 @@ def update_values(self, new_values, new_active_mask=None):
         :param new_values: numpy array ([depth, lat, lon], n_particles) with
                            new coordinates of the particles.
         :param new_active_mask: boolean array indicating active particles.
-        '''
+        """
         self.rebuild(new_values, new_active_mask)
 
     def rebuild(self, values, active_mask=-1):
@@ -90,7 +90,7 @@ def rebuild(self, values, active_mask=-1):
 
     @property
     def active_idx(self):
-        "Indices of the currently active mask."
+        """Indices of the currently active mask."""
         # See __init__ comments for a more detailed explanation.
         if self._active_mask is None:
             return np.arange(self._values.shape[1])
@@ -128,7 +128,8 @@ def _get_close_neighbor_dist(self, coor, subset_idx):
 
 
 class BaseFlatNeighborSearch(BaseNeighborSearch):
-    "Base class for neighbor searches with a flat mesh."
+    """Base class for neighbor searches with a flat mesh."""
+
     def _distance(self, coor, subset_idx):
         coor = coor.reshape(3, 1)
         horiz_distance = np.sqrt(np.sum((
@@ -158,7 +159,8 @@ def _distance(self, coor, subset_idx):
 
 
 class BaseSphericalNeighborSearch(BaseNeighborSearch):
-    "Base class for a neighbor search with a spherical mesh."
+    """Base class for a neighbor search with a spherical mesh."""
+
     def _distance(self, coor, subset_idx):
         vert_distances, horiz_distances = spherical_distance(
             *coor,
diff --git a/parcels/interaction/neighborsearch/basehash.py b/parcels/interaction/neighborsearch/basehash.py
index 17b9ecedd..550f94f85 100644
--- a/parcels/interaction/neighborsearch/basehash.py
+++ b/parcels/interaction/neighborsearch/basehash.py
@@ -5,23 +5,23 @@
 
 class BaseHashNeighborSearch(ABC):
     def find_neighbors_by_coor(self, coor):
-        '''Get the neighbors around a certain location.
+        """Get the neighbors around a certain location.
 
         :param coor: Numpy array with [depth, lat, lon].
         :returns List of particle indices.
-        '''
+        """
         coor = coor.reshape(3, 1)
         hash_id = self._values_to_hashes(coor)[0]
         return self._find_neighbors(hash_id, coor)
 
     def find_neighbors_by_idx(self, particle_idx):
-        '''Get the neighbors around a certain particle.
+        """Get the neighbors around a certain particle.
 
         Mainly useful for Structure of Array (SoA) datastructure
 
         :param particle_idx: index of the particle (SoA).
         :returns List of particle indices
-        '''
+        """
         hash_id = self._particle_hashes[particle_idx]
         coor = self._values[:, particle_idx].reshape(3, 1)
         return self._find_neighbors(hash_id, coor)
@@ -31,10 +31,10 @@ def _find_neighbors(self, hash_id, coor):
         raise NotImplementedError
 
     def consistency_check(self):
-        '''See if all values are in their proper place.
+        """See if all values are in their proper place.
 
         Only used for debugging purposes.
-        '''
+        """
         active_idx = self.active_idx
         if active_idx is None:
             active_idx = np.arange(self._values.shape[1])
@@ -54,13 +54,13 @@ def consistency_check(self):
         assert np.all(cur_hashes == self._particle_hashes[active_idx])
 
     def update_values(self, new_values, new_active_mask=None):
-        '''Update the locations of (some) of the particles.
+        """Update the locations of (some) of the particles.
 
         Particles that stay in the same location are computationally cheap.
         The order and number of the particles is assumed to remain the same.
 
         :param new_values: new (depth, lat, lon) values for particles.
-        '''
+        """
         if self._values is None:
             self.rebuild(new_values, new_active_mask)
             return
@@ -110,7 +110,7 @@ def _values_to_hashes(self, values, active_idx=None):
         raise NotImplementedError
 
     def _deactivate_particles(self, particle_idx):
-        "Remove particles from the hashtable."
+        """Remove particles from the hashtable."""
         # Get the hashes of the particles to be removed.
         remove_split = hash_split(self._particle_hashes[particle_idx])
         for cur_hash, remove_idx in remove_split.items():
@@ -127,7 +127,7 @@ def _deactivate_particles(self, particle_idx):
                     len(self._hashtable[cur_hash]))
 
     def _activate_particles(self, particle_idx):
-        "Add particles to the hashtable"
+        """Add particles to the hashtable"""
         # See _deactivate_particles.
         add_split = hash_split(self._particle_hashes[particle_idx])
         for cur_hash, add_idx in add_split.items():
@@ -143,7 +143,7 @@ def _activate_particles(self, particle_idx):
 
 
 def hash_split(hash_ids, active_idx=None):
-    '''Create a hashtable.
+    """Create a hashtable.
 
     Multiple particles that are found in the same cell are put in a list
     with that particular hash.
@@ -151,7 +151,7 @@ def hash_split(hash_ids, active_idx=None):
     :param hash_ids: Hash values for the particles.
     :param active_idx: Subset on which to compute the hash split.
     :returns hash_split: Dictionary with {hash: [idx_1, idx_2, ..], ..}
-    '''
+    """
     if len(hash_ids) == 0:
         return {}
     if active_idx is not None:
diff --git a/parcels/interaction/neighborsearch/bruteforce.py b/parcels/interaction/neighborsearch/bruteforce.py
index 22b2301df..f61df9005 100644
--- a/parcels/interaction/neighborsearch/bruteforce.py
+++ b/parcels/interaction/neighborsearch/bruteforce.py
@@ -5,14 +5,14 @@
 
 
 class BruteFlatNeighborSearch(BaseFlatNeighborSearch):
-    '''Brute force implementation to find the neighbors.'''
+    """Brute force implementation to find the neighbors."""
 
     def find_neighbors_by_coor(self, coor):
         return self._get_close_neighbor_dist(coor, self.active_idx)
 
 
 class BruteSphericalNeighborSearch(BaseSphericalNeighborSearch):
-    '''Brute force implementation to find the neighbors.'''
+    """Brute force implementation to find the neighbors."""
 
     def find_neighbors_by_coor(self, coor):
         return self._get_close_neighbor_dist(coor, self.active_idx)
diff --git a/parcels/interaction/neighborsearch/distanceutils.py b/parcels/interaction/neighborsearch/distanceutils.py
index 8059411c8..56a70ac7f 100644
--- a/parcels/interaction/neighborsearch/distanceutils.py
+++ b/parcels/interaction/neighborsearch/distanceutils.py
@@ -2,18 +2,18 @@
 
 
 def fast_distance(lat1, lon1, lat2, lon2):
-    '''Compute the arc distance assuming the earth is a sphere.
+    """Compute the arc distance assuming the earth is a sphere.
 
     This is not the only possible implementation. It was taken from:
     https://www.mkompf.com/gps/distcalc.html
-    '''
+    """
     g = np.sin(lat1)*np.sin(lat2)+np.cos(lat1)*np.cos(lat2)*np.cos(lon1-lon2)
     return np.arccos(np.minimum(1, g))
 
 
 def spherical_distance(depth1_m, lat1_deg, lon1_deg, depth2_m, lat2_deg,
                        lon2_deg):
-    "Compute the arc distance, uses degrees as input."
+    """Compute the arc distance, uses degrees as input."""
     R_earth = 6371000
     lat1 = np.pi*lat1_deg/180
     lon1 = np.pi*lon1_deg/180
diff --git a/parcels/interaction/neighborsearch/hashflat.py b/parcels/interaction/neighborsearch/hashflat.py
index 33d98d29a..cabecd507 100644
--- a/parcels/interaction/neighborsearch/hashflat.py
+++ b/parcels/interaction/neighborsearch/hashflat.py
@@ -8,7 +8,8 @@
 
 
 class HashFlatNeighborSearch(BaseHashNeighborSearch, BaseFlatNeighborSearch):
-    '''Neighbor search using a hashtable (similar to octtrees).'''
+    """Neighbor search using a hashtable (similar to octtrees)."""
+
     _box = None
 
     def _find_neighbors(self, hash_id, coor):
diff --git a/parcels/interaction/neighborsearch/hashspherical.py b/parcels/interaction/neighborsearch/hashspherical.py
index 8821e9ec1..df770b246 100644
--- a/parcels/interaction/neighborsearch/hashspherical.py
+++ b/parcels/interaction/neighborsearch/hashspherical.py
@@ -11,22 +11,23 @@
 
 class HashSphericalNeighborSearch(BaseHashNeighborSearch,
                                   BaseSphericalNeighborSearch):
-    '''Neighbor search using a hashtable (similar to octtrees).'''
+    """Neighbor search using a hashtable (similar to octtrees)."""
+
     def __init__(self, inter_dist_vert, inter_dist_horiz,
                  max_depth=100000):
-        '''Initialize the neighbor data structure.
+        """Initialize the neighbor data structure.
 
         :param interaction_distance: maximum horizontal interaction distance.
         :param interaction_depth: maximum depth of interaction.
         :param values: depth, lat, lon values for particles.
         :param max_depth: maximum depth of the ocean.
-        '''
+        """
         super().__init__(inter_dist_vert, inter_dist_horiz, max_depth)
 
         self._init_structure()
 
     def _find_neighbors(self, hash_id, coor):
-        '''Get neighbors from hash_id and location.'''
+        """Get neighbors from hash_id and location."""
         # Get the neighboring cells.
         neighbor_blocks = geo_hash_to_neighbors(
             hash_id, coor, self._bits, self.inter_arc_dist)
@@ -43,12 +44,12 @@ def _find_neighbors(self, hash_id, coor):
         return self._get_close_neighbor_dist(coor, potential_neighbors)
 
     def _values_to_hashes(self, values, active_idx=None):
-        '''Convert coordinates to cell ids.
+        """Convert coordinates to cell ids.
 
         :param values: array of positions of particles to convert
                        ([depth, lat, lon], # of particles to convert).
         :returns array of cell ids.
-        '''
+        """
         if active_idx is None:
             active_idx = np.arange(values.shape[1], dtype=int)
         depth = values[0, active_idx]
@@ -78,10 +79,10 @@ def _values_to_hashes(self, values, active_idx=None):
         return point_array
 
     def rebuild(self, values, active_mask=-1):
-        '''Recreate the tree with new values.
+        """Recreate the tree with new values.
 
         :param values: positions of the particles.
-        '''
+        """
         super().rebuild(values, active_mask)
         active_idx = self.active_idx
 
@@ -100,7 +101,7 @@ def rebuild(self, values, active_mask=-1):
             self._hash_idx[idx_array] = np.arange(len(idx_array))
 
     def _init_structure(self):
-        '''Initialize the basic tree properties without building'''
+        """Initialize the basic tree properties without building"""
         epsilon = 1e-12
         R_earth = 6371000
 
@@ -117,7 +118,7 @@ def _init_structure(self):
 
 
 def i_3d_to_hash(i_depth, i_lat, i_lon, lat_sign, bits):
-    '''Convert longitude and lattitude id's to hash'''
+    """Convert longitude and lattitude id's to hash"""
     point_hash = lat_sign
     point_hash = np.bitwise_or(point_hash, np.left_shift(i_depth, 1))
     point_hash = np.bitwise_or(point_hash, np.left_shift(i_lat, 1+bits[0]))
@@ -126,7 +127,7 @@ def i_3d_to_hash(i_depth, i_lat, i_lon, lat_sign, bits):
 
 
 def geo_hash_to_neighbors(hash_id, coor, bits, inter_arc_dist):
-    '''Compute the hashes of all neighboring cells in a 3x3x3 neighborhood.'''
+    """Compute the hashes of all neighboring cells in a 3x3x3 neighborhood."""
     lat_sign = hash_id & 0x1
     i_depth = (hash_id >> 1) & ((1 << bits[0])-1)
     i_lat = (hash_id >> (1+bits[0])) & ((1 << bits[1])-1)
diff --git a/parcels/kernel/basekernel.py b/parcels/kernel/basekernel.py
index ec9912174..3a8eeaebb 100644
--- a/parcels/kernel/basekernel.py
+++ b/parcels/kernel/basekernel.py
@@ -54,6 +54,7 @@ class BaseKernel:
     The py_ast argument may be derived from the code string, but for
     concatenation, the merged AST plus the new header definition is required.
     """
+
     _pyfunc = None
     _fieldset = None
     _ptype = None
@@ -134,7 +135,7 @@ def _cache_key(self):
 
     @staticmethod
     def fix_indentation(string):
-        """Fix indentation to allow in-lined kernel definitions"""
+        """Fix indentation to allow in-lined kernel definitions."""
         lines = string.split('\n')
         indent = re_indent.match(lines[0])
         if indent:
@@ -143,7 +144,8 @@ def fix_indentation(string):
 
     def check_fieldsets_in_kernels(self, pyfunc):
         """
-        function checks the integrity of the fieldset with the kernels.
+        Checks the integrity of the fieldset with the kernels.
+
         This function is to be called from the derived class when setting up the 'pyfunc'.
         """
         if self.fieldset is not None:
@@ -168,9 +170,6 @@ def check_fieldsets_in_kernels(self, pyfunc):
                     raise NotImplementedError('Analytical Advection only works with Z-grids in the vertical')
 
     def check_kernel_signature_on_version(self):
-        """
-        returns numkernelargs
-        """
         numkernelargs = 0
         if self._pyfunc is not None:
             if version_info[0] < 3:
@@ -207,9 +206,7 @@ def remove_lib(self):
                 self.src_file = src_file_or_files
 
     def get_kernel_compile_files(self):
-        """
-        Returns the correct src_file, lib_file, log_file for this kernel
-        """
+        """Returns the correct src_file, lib_file, log_file for this kernel."""
         if MPI:
             mpi_comm = MPI.COMM_WORLD
             mpi_rank = mpi_comm.Get_rank()
@@ -236,7 +233,7 @@ def get_kernel_compile_files(self):
         return src_file_or_files, lib_file, log_file
 
     def compile(self, compiler):
-        """ Writes kernel code to file and compiles it."""
+        """Writes kernel code to file and compiles it."""
         all_files_array = []
         if self.src_file is None:
             if self.dyn_srcs is not None:
@@ -315,9 +312,7 @@ def remove_deleted(self, pset, output_file, endtime):
         pset.remove_indices(indices)
 
     def load_fieldset_jit(self, pset):
-        """
-        Updates the loaded fields of pset's fieldset according to the chunk information within their grids
-        """
+        """Updates the loaded fields of pset's fieldset according to the chunk information within their grids."""
         if pset.fieldset is not None:
             for g in pset.fieldset.gridset.grids:
                 g.cstruct = None  # This force to point newly the grids from Python to C
diff --git a/parcels/kernel/kernelaos.py b/parcels/kernel/kernelaos.py
index 7e1261459..8fa7dd892 100644
--- a/parcels/kernel/kernelaos.py
+++ b/parcels/kernel/kernelaos.py
@@ -129,7 +129,7 @@ def __radd__(self, kernel):
         return kernel.merge(self, KernelAOS)
 
     def execute_jit(self, pset, endtime, dt):
-        """Invokes JIT engine to perform the core update loop"""
+        """Invokes JIT engine to perform the core update loop."""
         self.load_fieldset_jit(pset)
 
         fargs = []
@@ -145,7 +145,7 @@ def execute_jit(self, pset, endtime, dt):
             self._function(c_int(len(pset)), pdata, c_double(endtime), c_double(dt))
 
     def execute_python(self, pset, endtime, dt):
-        """Performs the core update loop via Python"""
+        """Performs the core update loop via Python."""
         # sign of dt: { [0, 1]: forward simulation; -1: backward simulation }
         sign_dt = np.sign(dt)
 
@@ -166,14 +166,14 @@ def execute_python(self, pset, endtime, dt):
             self.evaluate_particle(p, endtime, sign_dt, dt, analytical=analytical)
 
     def remove_deleted(self, pset, output_file, endtime):
-        """Utility to remove all particles that signalled deletion"""
+        """Utility to remove all particles that signalled deletion."""
         indices = [i for i, p in enumerate(pset) if p.state == OperationCode.Delete]
         if len(indices) > 0 and output_file is not None:
             output_file.write(pset, endtime, deleted_only=indices)
         pset.remove_indices(indices)
 
     def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False):
-        """Execute this Kernel over a ParticleSet for several timesteps"""
+        """Execute this Kernel over a ParticleSet for several timesteps."""
         for p in pset:
             p.reset_state()
 
diff --git a/parcels/kernel/kernelsoa.py b/parcels/kernel/kernelsoa.py
index d0818e9d8..29d2a5633 100644
--- a/parcels/kernel/kernelsoa.py
+++ b/parcels/kernel/kernelsoa.py
@@ -116,7 +116,7 @@ def __init__(self, fieldset, ptype, pyfunc=None, funcname=None,
                 self.src_file = src_file_or_files
 
     def execute_jit(self, pset, endtime, dt):
-        """Invokes JIT engine to perform the core update loop"""
+        """Invokes JIT engine to perform the core update loop."""
         self.load_fieldset_jit(pset)
 
         fargs = [byref(f.ctypes_struct) for f in self.field_args.values()]
@@ -126,7 +126,7 @@ def execute_jit(self, pset, endtime, dt):
                               c_double(endtime), c_double(dt), *fargs)
 
     def execute_python(self, pset, endtime, dt):
-        """Performs the core update loop via Python"""
+        """Performs the core update loop via Python."""
         # sign of dt: { [0, 1]: forward simulation; -1: backward simulation }
         sign_dt = np.sign(dt)
 
@@ -163,8 +163,7 @@ def __radd__(self, kernel):
         return kernel.merge(self, KernelSOA)
 
     def remove_deleted(self, pset, output_file, endtime):
-        """
-        Utility to remove all particles that signalled deletion
+        """Utility to remove all particles that signalled deletion.
 
         This deletion function is targetted to index-addressable, random-access array-collections.
         """
@@ -176,7 +175,7 @@ def remove_deleted(self, pset, output_file, endtime):
         pset.remove_indices(indices)
 
     def execute(self, pset, endtime, dt, recovery=None, output_file=None, execute_once=False):
-        """Execute this Kernel over a ParticleSet for several timesteps"""
+        """Execute this Kernel over a ParticleSet for several timesteps."""
         pset.collection.state[:] = StateCode.Evaluate
 
         if abs(dt) < 1e-6 and not execute_once:
diff --git a/parcels/particle.py b/parcels/particle.py
index 2d2aa07fb..9bf2a3798 100644
--- a/parcels/particle.py
+++ b/parcels/particle.py
@@ -13,7 +13,7 @@
 
 
 class Variable:
-    """Descriptor class that delegates data access to particle data
+    """Descriptor class that delegates data access to particle data.
 
     :param name: Variable name as used within kernels
     :param dtype: Data type (numpy.dtype) of the variable
@@ -23,6 +23,7 @@ class Variable:
              If to_write = 'once', the variable will be written as a time-independent 1D array
     :type to_write: (bool, 'once', optional)
     """
+
     def __init__(self, name, dtype=np.float32, initial=0, to_write=True):
         self.name = name
         self.dtype = dtype
@@ -47,12 +48,12 @@ def __repr__(self):
         return f"PVar<{self.name}|{self.dtype}>"
 
     def is64bit(self):
-        """Check whether variable is 64-bit"""
+        """Check whether variable is 64-bit."""
         return True if self.dtype in indicators_64bit else False
 
 
 class ParticleType:
-    """Class encapsulating the type information for custom particles
+    """Class encapsulating the type information for custom particles.
 
     :param user_vars: Optional list of (name, dtype) tuples for custom variables
     """
@@ -96,7 +97,7 @@ def _cache_key(self):
 
     @property
     def dtype(self):
-        """Numpy.dtype object that defines the C struct"""
+        """Numpy.dtype object that defines the C struct."""
         type_list = [(v.name, v.dtype) for v in self.variables]
         for v in self.variables:
             if v.dtype not in self.supported_dtypes:
@@ -108,13 +109,12 @@ def dtype(self):
 
     @property
     def size(self):
-        """Size of the underlying particle struct in bytes"""
+        """Size of the underlying particle struct in bytes."""
         return sum([8 if v.is64bit() else 4 for v in self.variables])
 
     @property
     def supported_dtypes(self):
-        """List of all supported numpy dtypes. All others are not supported"""
-
+        """List of all supported numpy dtypes. All others are not supported."""
         # Developer note: other dtypes (mostly 2-byte ones) are not supported now
         # because implementing and aligning them in cgen.GenerableStruct is a
         # major headache. Perhaps in a later stage
@@ -122,7 +122,8 @@ def supported_dtypes(self):
 
 
 class _Particle:
-    """Private base class for all particle types"""
+    """Private base class for all particle types."""
+
     lastID = 0  # class-level variable keeping track of last Particle ID used
 
     def __init__(self):
@@ -169,8 +170,7 @@ def setLastID(cls, offset):
 
 
 class ScipyParticle(_Particle):
-    """Class encapsulating the basic attributes of a particle,
-    to be executed in SciPy mode
+    """Class encapsulating the basic attributes of a particle, to be executed in SciPy mode.
 
     :param lon: Initial longitude of particle
     :param lat: Initial latitude of particle
@@ -230,7 +230,7 @@ class ScipyInteractionParticle(ScipyParticle):
 
 
 class JITParticle(ScipyParticle):
-    """Particle class for JIT-based (Just-In-Time) Particle objects
+    """Particle class for JIT-based (Just-In-Time) Particle objects.
 
     :param lon: Initial longitude of particle
     :param lat: Initial latitude of particle
diff --git a/parcels/particlefile/__init__.py b/parcels/particlefile/__init__.py
index 9e5d618fc..a994ce0b1 100644
--- a/parcels/particlefile/__init__.py
+++ b/parcels/particlefile/__init__.py
@@ -1,4 +1,4 @@
-"""Module controlling the writing of ParticleSets to NetCDF file"""
+"""Module controlling the writing of ParticleSets to NetCDF file."""
 
 from .baseparticlefile import BaseParticleFile  # noqa: F401
 from .baseparticlefile import _set_calendar  # noqa: F401
diff --git a/parcels/particlefile/baseparticlefile.py b/parcels/particlefile/baseparticlefile.py
index 424b65945..e7b676bf3 100644
--- a/parcels/particlefile/baseparticlefile.py
+++ b/parcels/particlefile/baseparticlefile.py
@@ -1,4 +1,4 @@
-"""Module controlling the writing of ParticleSets to Zarr file"""
+"""Module controlling the writing of ParticleSets to Zarr file."""
 import os
 from abc import ABC, abstractmethod
 from datetime import timedelta as delta
@@ -42,6 +42,7 @@ class BaseParticleFile(ABC):
     :param write_ondelete: Boolean to write particle data only when they are deleted. Default is False
     :param create_new_zarrfile: Boolean to create a new file. Default is True
     """
+
     write_ondelete = None
     outputdt = None
     lasttime_written = None
@@ -112,19 +113,16 @@ def __init__(self, name, particleset, outputdt=np.infty, chunks=None, write_onde
 
     @abstractmethod
     def _reserved_var_names(self):
-        """
-        returns the reserved dimension names not to be written just once.
-        """
+        """Returns the reserved dimension names not to be written just once."""
         pass
 
     def _create_variables_attribute_dict(self):
         """
-        creates the dictionary with variable attributes.
+        Creates the dictionary with variable attributes.
 
         Attention:
         For ParticleSet structures other than SoA, and structures where ID != index, this has to be overridden.
         """
-
         attrs = {'z': {"long_name": "",
                        "standard_name": "depth",
                        "units": "m",
@@ -166,7 +164,7 @@ def close(self, delete_tempfiles=True):
         pass
 
     def add_metadata(self, name, message):
-        """Add metadata to :class:`parcels.particleset.ParticleSet`
+        """Add metadata to :class:`parcels.particleset.ParticleSet`.
 
         :param name: Name of the metadata variabale
         :param message: message to be written
@@ -200,13 +198,12 @@ def _extend_zarr_dims(self, Z, store, dtype, axis):
         zarr.consolidate_metadata(store)
 
     def write(self, pset, time, deleted_only=False):
-        """Write all data from one time step to the zarr file
+        """Write all data from one time step to the zarr file.
 
         :param pset: ParticleSet object to write
         :param time: Time at which to write ParticleSet
         :param deleted_only: Flag to write only the deleted Particles
         """
-
         time = time.total_seconds() if isinstance(time, delta) else time
 
         if self.lasttime_written != time and (self.write_ondelete is False or deleted_only is not False):
diff --git a/parcels/particlefile/particlefileaos.py b/parcels/particlefile/particlefileaos.py
index 75f89c314..1d747e42a 100644
--- a/parcels/particlefile/particlefileaos.py
+++ b/parcels/particlefile/particlefileaos.py
@@ -27,7 +27,5 @@ def __del__(self):
         super().__del__()
 
     def _reserved_var_names(self):
-        """
-        returns the reserved dimension names not to be written just once.
-        """
+        """Returns the reserved dimension names not to be written just once."""
         return ['time', 'lat', 'lon', 'depth', 'id']
diff --git a/parcels/particlefile/particlefilesoa.py b/parcels/particlefile/particlefilesoa.py
index c47abf885..c5031e3f4 100644
--- a/parcels/particlefile/particlefilesoa.py
+++ b/parcels/particlefile/particlefilesoa.py
@@ -27,7 +27,5 @@ def __del__(self):
         super().__del__()
 
     def _reserved_var_names(self):
-        """
-        returns the reserved dimension names not to be written just once.
-        """
+        """Returns the reserved dimension names not to be written just once."""
         return ['time', 'lat', 'lon', 'depth', 'id']
diff --git a/parcels/particleset/baseparticleset.py b/parcels/particleset/baseparticleset.py
index fabdb5293..8dcca9fb9 100644
--- a/parcels/particleset/baseparticleset.py
+++ b/parcels/particleset/baseparticleset.py
@@ -25,6 +25,7 @@ class NDCluster(ABC):
 
 class BaseParticleSet(NDCluster):
     """Base ParticleSet."""
+
     _collection = None
     kernel = None
     interaction_kernel = None
@@ -60,9 +61,7 @@ def iterator(self):
         return self._collection.iterator()
 
     def __iter__(self):
-        """Allows for more intuitive iteration over a particleset, while
-        in reality iterating over the particles in the collection.
-        """
+        """Allows for more intuitive iteration over a particleset, while in reality iterating over the particles in the collection."""
         return self.iterator()
 
     def __getattr__(self, name):
@@ -109,11 +108,11 @@ def __create_progressbar(self, starttime, endtime):
 
     @classmethod
     def from_list(cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs):
-        """Initialise the ParticleSet from lists of lon and lat
+        """Initialise the ParticleSet from lists of lon and lat.
 
         :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity
         :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle`
-                 object that defines custom particle
+            object that defines custom particle
         :param lon: List of initial longitude values for particles
         :param lat: List of initial latitude values for particles
         :param depth: Optional list of initial depth values for particles. Default is 0m
@@ -123,12 +122,14 @@ def from_list(cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=N
                It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
                and np.float64 if the interpolation method is 'cgrid_velocity'
         Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
-       """
+        """
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt, lonlatdepth_dtype=lonlatdepth_dtype, **kwargs)
 
     @classmethod
     def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None):
-        """Initialise the ParticleSet from start/finish coordinates with equidistant spacing
+        """Create a particleset in the shape of a line (according to a cartesian grid).
+
+        Initialise the ParticleSet from start/finish coordinates with equidistant spacing
         Note that this method uses simple numpy.linspace calls and does not take into account
         great circles, so may not be a exact on a globe
 
@@ -145,7 +146,6 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None,
                It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
                and np.float64 if the interpolation method is 'cgrid_velocity'
         """
-
         lon = np.linspace(start[0], finish[0], size)
         lat = np.linspace(start[1], finish[1], size)
         if type(depth) in [int, float]:
@@ -166,7 +166,7 @@ def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'):
 
     @classmethod
     def from_field(cls, fieldset, pclass, start_field, size, mode='monte_carlo', depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None):
-        """Initialise the ParticleSet randomly drawn according to distribution from a field
+        """Initialise the ParticleSet randomly drawn according to distribution from a field.
 
         :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity
         :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle`
@@ -181,7 +181,6 @@ def from_field(cls, fieldset, pclass, start_field, size, mode='monte_carlo', dep
                It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
                and np.float64 if the interpolation method is 'cgrid_velocity'
         """
-
         lon, lat = cls.monte_carlo_sample(start_field, size, mode)
 
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, lonlatdepth_dtype=lonlatdepth_dtype, repeatdt=repeatdt)
@@ -210,8 +209,7 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
         pass
 
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
-        """Method to calculate the density of particles in a ParticleSet from their locations,
-        through a 2D histogram.
+        """Calculate 2D particle density field from ParticleSet particle locations.
 
         :param field: Optional :mod:`parcels.field.Field` object to calculate the histogram
                       on. Default is `fieldset.U`
@@ -227,8 +225,9 @@ def density(self, field_name=None, particle_val=None, relative=False, area_scale
 
     @abstractmethod
     def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
-        """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object
-        based on `fieldset` and `ptype` of the ParticleSet
+        """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object.
+
+        Conversion is based on `fieldset` and `ptype` of the ParticleSet.
         :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
         """
         pass
@@ -238,8 +237,7 @@ def InteractionKernel(self, pyfunc_inter):
 
     @abstractmethod
     def ParticleFile(self, *args, **kwargs):
-        """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile`
-        object from the ParticleSet"""
+        """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile` object from the ParticleSet."""
         pass
 
     @abstractmethod
@@ -296,8 +294,9 @@ def _impute_release_times(self, default):
     def execute(self, pyfunc=AdvectionRK4, pyfunc_inter=None, endtime=None, runtime=None, dt=1.,
                 moviedt=None, recovery=None, output_file=None, movie_background_field=None,
                 verbose_progress=None, postIterationCallbacks=None, callbackdt=None):
-        """Execute a given kernel function over the particle set for
-        multiple timesteps. Optionally also provide sub-timestepping
+        """Execute a given kernel function over the particle set for multiple timesteps.
+
+        Optionally also provide sub-timestepping
         for particle output.
 
         :param pyfunc: Kernel function to execute. This can be the name of a
@@ -514,7 +513,7 @@ def execute(self, pyfunc=AdvectionRK4, pyfunc_inter=None, endtime=None, runtime=
 
     def show(self, with_particles=True, show_time=None, field=None, domain=None, projection='PlateCarree',
              land=True, vmin=None, vmax=None, savefile=None, animation=False, **kwargs):
-        """Method to 'show' a Parcels ParticleSet
+        """Method to 'show' a Parcels ParticleSet.
 
         :param with_particles: Boolean whether to show particles
         :param show_time: Time at which to show the ParticleSet
diff --git a/parcels/particleset/particlesetaos.py b/parcels/particleset/particlesetaos.py
index cd8106cba..c2e4a2df3 100644
--- a/parcels/particleset/particlesetaos.py
+++ b/parcels/particleset/particlesetaos.py
@@ -395,7 +395,7 @@ def iterator(self):
         return super().iterator()
 
     def __getitem__(self, index):
-        """Get a single particle by index"""
+        """Get a single particle by index."""
         return self._collection.get_single_by_index(index)
 
     def __getattr__(self, name):
@@ -589,21 +589,20 @@ def add(self, particles):
         return self
 
     def remove_indices(self, indices):
-        """Method to remove particles from the ParticleSet, based on their `indices`"""
+        """Method to remove particles from the ParticleSet, based on their `indices`."""
         if type(indices) in [int, np.int32, np.intp]:
             self._collection.remove_single_by_index(indices)
         else:
             self._collection.remove_multi_by_indices(indices)
 
     def remove_booleanvector(self, indices):
-        """Method to remove particles from the ParticleSet, based on an array of booleans"""
+        """Method to remove particles from the ParticleSet, based on an array of booleans."""
         # indices = np.where(indices)[0]
         indices = np.nonzero(indices)[0]
         self.remove_indices(indices)
 
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
-        """Method to calculate the density of particles in a ParticleSet from their locations,
-        through a 2D histogram.
+        """Calculate 2D particle density field from ParticleSet particle locations.
 
         :param field: Optional :mod:`parcels.field.Field` object to calculate the histogram
                       on. Default is `fieldset.U`
@@ -615,7 +614,6 @@ def density(self, field_name=None, particle_val=None, relative=False, area_scale
         :param area_scale: Boolean to control whether the density is scaled by the area
                            (in m^2) of each grid cell. Default is False
         """
-
         field_name = field_name if field_name else "U"
         sampling_name = "UV" if field_name in ["U", "V"] else field_name
         field = getattr(self.fieldset, field_name)
@@ -658,7 +656,7 @@ def density(self, field_name=None, particle_val=None, relative=False, area_scale
 
     def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
         """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object
-        based on `fieldset` and `ptype` of the ParticleSet
+        based on `fieldset` and `ptype` of the ParticleSet.
 
         :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
         """
@@ -666,12 +664,13 @@ def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
 
     def ParticleFile(self, *args, **kwargs):
         """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile`
-        object from the ParticleSet"""
+        object from the ParticleSet.
+        """
         return ParticleFileAOS(*args, particleset=self, **kwargs)
 
     def set_variable_write_status(self, var, write_status):
         """
-        Method to set the write status of a Variable
+        Method to set the write status of a Variable.
         :param var: Name of the variable (string)
         :param write_status: Write status of the variable (True, False or 'once')
         """
diff --git a/parcels/particleset/particlesetsoa.py b/parcels/particleset/particlesetsoa.py
index 97e2d1c8c..3b18cb3bb 100644
--- a/parcels/particleset/particlesetsoa.py
+++ b/parcels/particleset/particlesetsoa.py
@@ -39,9 +39,7 @@
 
 
 def _convert_to_reltime(time):
-    """Check to determine if the value of the time parameter needs to be
-    converted to a relative value (relative to the time_origin).
-    """
+    """Check to determine if the value of the time parameter needs to be converted to a relative value (relative to the time_origin)."""
     if isinstance(time, np.datetime64) or (hasattr(time, 'calendar') and time.calendar in _get_cftime_calendars()):
         return True
     return False
@@ -266,8 +264,7 @@ def _impute_release_times(self, default):
         return np.min(self._collection.data['time']), np.max(self._collection.data['time'])
 
     def data_indices(self, variable_name, compare_values, invert=False):
-        """Get the indices of all particles where the value of
-        `variable_name` equals (one of) `compare_values`.
+        """Get the indices of all particles where the value of `variable_name` equals (one of) `compare_values`.
 
         :param variable_name: Name of the variable to check.
         :param compare_values: Value or list of values to compare to.
@@ -288,7 +285,6 @@ def populate_indices(self):
         This is only intended for curvilinear grids, where the initial index search
         may be quite expensive.
         """
-
         if self.fieldset is None:
             # we need to be attached to a fieldset to have a valid
             # gridset to search for indices
@@ -339,7 +335,7 @@ def num_error_particles(self):
             [StateCode.Success, StateCode.Evaluate], invert=True))
 
     def __getitem__(self, index):
-        """Get a single particle by index"""
+        """Get a single particle by index."""
         return self._collection.get_single_by_index(index)
 
     def cstruct(self):
@@ -421,7 +417,6 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
                It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
                and np.float64 if the interpolation method is 'cgrid_velocity'
         """
-
         if repeatdt is not None:
             logger.warning(f'Note that the `repeatdt` argument is not retained from {filename}, and that '
                            'setting a new repeatdt will start particles from the _new_ particle '
@@ -515,7 +510,9 @@ def __sizeof__(self):
         return sys.getsizeof(self._collection)
 
     def __iadd__(self, particles):
-        """Add particles to the ParticleSet. Note that this is an
+        """Add particles to the ParticleSet.
+
+        Note that this is an
         incremental add, the particles will be added to the ParticleSet
         on which this function is called.
 
@@ -548,7 +545,7 @@ def add(self, particles):
         return self
 
     def remove_indices(self, indices):
-        """Method to remove particles from the ParticleSet, based on their `indices`"""
+        """Method to remove particles from the ParticleSet, based on their `indices`."""
         # Removing particles invalidates the neighbor search structure.
         self._dirty_neighbor = True
         if type(indices) in [int, np.int32, np.intp]:
@@ -557,14 +554,13 @@ def remove_indices(self, indices):
             self._collection.remove_multi_by_indices(indices)
 
     def remove_booleanvector(self, indices):
-        """Method to remove particles from the ParticleSet, based on an array of booleans"""
+        """Method to remove particles from the ParticleSet, based on an array of booleans."""
         # Removing particles invalidates the neighbor search structure.
         self._dirty_neighbor = True
         self.remove_indices(np.where(indices)[0])
 
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
-        """Method to calculate the density of particles in a ParticleSet from their locations,
-        through a 2D histogram.
+        """Calculate 2D particle density field from ParticleSet particle locations.
 
         :param field: Optional :mod:`parcels.field.Field` object to calculate the histogram
                       on. Default is `fieldset.U`
@@ -576,7 +572,6 @@ def density(self, field_name=None, particle_val=None, relative=False, area_scale
         :param area_scale: Boolean to control whether the density is scaled by the area
                            (in m^2) of each grid cell. Default is False
         """
-
         field_name = field_name if field_name else "U"
         sampling_name = "UV" if field_name in ["U", "V"] else field_name
         field = getattr(self.fieldset, field_name)
@@ -618,9 +613,9 @@ def density(self, field_name=None, particle_val=None, relative=False, area_scale
         return density
 
     def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
-        """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object
-        based on `fieldset` and `ptype` of the ParticleSet
+        """Wrapper method to convert a `pyfunc` into a :class:`parcels.kernel.Kernel` object.
 
+        Conversion is based on `fieldset` and `ptype` of the ParticleSet.
         :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
         """
         return Kernel(self.fieldset, self.collection.ptype, pyfunc=pyfunc, c_include=c_include,
@@ -632,13 +627,13 @@ def InteractionKernel(self, pyfunc_inter):
         return InteractionKernelSOA(self.fieldset, self.collection.ptype, pyfunc=pyfunc_inter)
 
     def ParticleFile(self, *args, **kwargs):
-        """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile`
-        object from the ParticleSet"""
+        """Wrapper method to initialise a :class:`parcels.particlefile.ParticleFile` object from the ParticleSet."""
         return ParticleFile(*args, particleset=self, **kwargs)
 
     def set_variable_write_status(self, var, write_status):
         """
-        Method to set the write status of a Variable
+        Method to set the write status of a Variable.
+
         :param var: Name of the variable (string)
         :param write_status: Write status of the variable (True, False or 'once')
         """
diff --git a/parcels/plotting.py b/parcels/plotting.py
index fe51ad279..dcc6eb590 100644
--- a/parcels/plotting.py
+++ b/parcels/plotting.py
@@ -13,7 +13,7 @@
 def plotparticles(particles, with_particles=True, show_time=None, field=None, domain=None,
                   projection='PlateCarree', land=True, vmin=None, vmax=None, savefile=None,
                   animation=False, **kwargs):
-    """Function to plot a Parcels ParticleSet
+    """Function to plot a Parcels ParticleSet.
 
     :param show_time: Time in seconds from start after which to show the ParticleSet
     :param with_particles: Boolean whether particles are also plotted on Field
@@ -26,7 +26,6 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do
     :param savefile: Name of a file to save the plot to
     :param animation: Boolean whether result is a single plot, or an animation
     """
-
     show_time = particles[0].time if show_time is None else show_time
     if isinstance(show_time, datetime):
         show_time = np.datetime64(show_time)
@@ -100,7 +99,7 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do
 
 def plotfield(field, show_time=None, domain=None, depth_level=0, projection='PlateCarree', land=True,
               vmin=None, vmax=None, savefile=None, **kwargs):
-    """Function to plot a Parcels Field
+    """Function to plot a Parcels Field.
 
     :param show_time: Time in seconds from start after which to show the Field
     :param domain: dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show
@@ -112,7 +111,6 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
     :param savefile: Name of a file to save the plot to
     :param animation: Boolean whether result is a single plot, or an animation
     """
-
     if type(field) is VectorField:
         spherical = True if field.U.grid.mesh == 'spherical' else False
         field = [field.U, field.V]
diff --git a/parcels/rng.py b/parcels/rng.py
index 0f6585dca..b77405605 100644
--- a/parcels/rng.py
+++ b/parcels/rng.py
@@ -133,13 +133,13 @@ def _assign_parcels_random_ccodeconverter():
 
 
 def seed(seed):
-    """Sets the seed for parcels internal RNG"""
+    """Sets the seed for parcels internal RNG."""
     _assign_parcels_random_ccodeconverter()
     _parcels_random_ccodeconverter.lib.pcls_seed(c_int(seed))
 
 
 def random():
-    """Returns a random float between 0. and 1."""
+    """Returns a random float between 0.0 and 1.0."""
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_random
     rnd.argtype = []
@@ -148,7 +148,7 @@ def random():
 
 
 def uniform(low, high):
-    """Returns a random float between `low` and `high`"""
+    """Returns a random float between `low` and `high`."""
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_uniform
     rnd.argtype = [c_float, c_float]
@@ -157,7 +157,7 @@ def uniform(low, high):
 
 
 def randint(low, high):
-    """Returns a random int between `low` and `high`"""
+    """Returns a random int between `low` and `high`."""
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_randint
     rnd.argtype = [c_int, c_int]
@@ -166,7 +166,7 @@ def randint(low, high):
 
 
 def normalvariate(loc, scale):
-    """Returns a random float on normal distribution with mean `loc` and width `scale`"""
+    """Returns a random float on normal distribution with mean `loc` and width `scale`."""
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_normalvariate
     rnd.argtype = [c_float, c_float]
@@ -175,7 +175,7 @@ def normalvariate(loc, scale):
 
 
 def expovariate(lamb):
-    """Returns a randome float of an exponential distribution with parameter lamb"""
+    """Returns a randome float of an exponential distribution with parameter lamb."""
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_expovariate
     rnd.argtype = c_float
@@ -185,7 +185,8 @@ def expovariate(lamb):
 
 def vonmisesvariate(mu, kappa):
     """Returns a randome float of a Von Mises distribution
-    with mean angle mu and concentration parameter kappa"""
+    with mean angle mu and concentration parameter kappa.
+    """
     _assign_parcels_random_ccodeconverter()
     rnd = _parcels_random_ccodeconverter.lib.pcls_vonmisesvariate
     rnd.argtype = [c_float, c_float]
diff --git a/parcels/scripts/plottrajectoriesfile.py b/parcels/scripts/plottrajectoriesfile.py
index 0a1fa76d8..c3ce043e1 100644
--- a/parcels/scripts/plottrajectoriesfile.py
+++ b/parcels/scripts/plottrajectoriesfile.py
@@ -17,7 +17,7 @@
 def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
                          tracerlon='x', tracerlat='y', recordedvar=None, movie_forward=True,
                          bins=20, show_plt=True, central_longitude=0):
-    """Quick and simple plotting of Parcels trajectories
+    """Quick and simple plotting of Parcels trajectories.
 
     :param filename: Name of Parcels-generated NetCDF file with particle positions
     :param mode: Type of plot to show. Supported are '2d', '3d', 'hist2d',
@@ -34,7 +34,6 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
     :param show_plt: Boolean whether plot should directly be show (for py.test)
     :param central_longitude: Degrees East at which to center the plot
     """
-
     environ["HDF5_USE_FILE_LOCKING"] = "FALSE"
     try:
         pfile = xr.open_zarr(str(filename), decode_cf=True)
diff --git a/parcels/tools/loggers.py b/parcels/tools/loggers.py
index c15b180e5..04bd233b2 100644
--- a/parcels/tools/loggers.py
+++ b/parcels/tools/loggers.py
@@ -1,4 +1,4 @@
-"""Script to create a `logger` for Parcels"""
+"""Script to create a `logger` for Parcels."""
 import logging
 
 __all__ = ['logger', 'XarrayDecodedFilter']
@@ -8,8 +8,8 @@
 
 
 class DuplicateFilter:
-    """Utility class to prevent warning_once warnings from being
-    displayed more than once"""
+    """Utility class to prevent warning_once warnings from being displayed more than once."""
+
     def __init__(self):
         self.msgs = set()
 
@@ -21,13 +21,13 @@ def filter(self, record):
 
 
 def warning_once(self, message, *args, **kws):
-    """Custom logging level for warnings that need to be displayed only once"""
+    """Custom logging level for warnings that need to be displayed only once."""
     if self.isEnabledFor(warning_once_level):
         self._log(warning_once_level, message, args, **kws)
 
 
 def info_once(self, message, *args, **kws):
-    """Custom logging level for info that need to be displayed only once"""
+    """Custom logging level for info that need to be displayed only once."""
     if self.isEnabledFor(info_once_level):
         self._log(info_once_level, message, args, **kws)
 
@@ -49,6 +49,7 @@ def info_once(self, message, *args, **kws):
 
 
 class XarrayDecodedFilter(logging.Filter):
-    """Filters the warning_once from fieldfilebuffer when cf_decoding fails"""
+    """Filters the warning_once from fieldfilebuffer when cf_decoding fails."""
+
     def filter(self, record):
         return 'Filling values might be wrongly parsed' not in record.getMessage()
diff --git a/parcels/tools/statuscodes.py b/parcels/tools/statuscodes.py
index 4c366cadb..5e07d6028 100644
--- a/parcels/tools/statuscodes.py
+++ b/parcels/tools/statuscodes.py
@@ -1,4 +1,4 @@
-"""Collection of pre-built recovery kernels"""
+"""Collection of pre-built recovery kernels."""
 
 
 __all__ = ['StateCode', 'OperationCode', 'ErrorCode',
@@ -27,9 +27,7 @@ class ErrorCode:
 
 
 class DaskChunkingError(RuntimeError):
-    """
-    Error indicating to the user that something with setting up Dask and chunked fieldsets went wrong.
-    """
+    """Error indicating to the user that something with setting up Dask and chunked fieldsets went wrong."""
 
     def __init__(self, src_class_type, message):
         msg = f"[{str(src_class_type)}]: {message}"
@@ -37,7 +35,7 @@ def __init__(self, src_class_type, message):
 
 
 class FieldSamplingError(RuntimeError):
-    """Utility error class to propagate erroneous field sampling in Scipy mode"""
+    """Utility error class to propagate erroneous field sampling in Scipy mode."""
 
     def __init__(self, x, y, z, field=None):
         self.field = field
@@ -49,7 +47,7 @@ def __init__(self, x, y, z, field=None):
 
 
 class FieldOutOfBoundError(RuntimeError):
-    """Utility error class to propagate out-of-bound field sampling in Scipy mode"""
+    """Utility error class to propagate out-of-bound field sampling in Scipy mode."""
 
     def __init__(self, x, y, z, field=None):
         self.field = field
@@ -62,7 +60,9 @@ def __init__(self, x, y, z, field=None):
 
 class FieldOutOfBoundSurfaceError(RuntimeError):
     """Utility error class to propagate out-of-bound field sampling at the surface in Scipy mode.
-       Note that if ErrorThroughSurface is not part of the recovery kernel, Parcels will use ErrorOutOfBounds."""
+
+    Note that if ErrorThroughSurface is not part of the recovery kernel, Parcels will use ErrorOutOfBounds.
+    """
 
     def __init__(self, x, y, z, field=None):
         self.field = field
@@ -74,7 +74,7 @@ def __init__(self, x, y, z, field=None):
 
 
 class TimeExtrapolationError(RuntimeError):
-    """Utility error class to propagate erroneous time extrapolation sampling in Scipy mode"""
+    """Utility error class to propagate erroneous time extrapolation sampling in Scipy mode."""
 
     def __init__(self, time, field=None, msg='allow_time_extrapoltion'):
         if field is not None and field.grid.time_origin and time is not None:
@@ -90,7 +90,7 @@ def __init__(self, time, field=None, msg='allow_time_extrapoltion'):
 
 
 class KernelError(RuntimeError):
-    """General particle kernel error with optional custom message"""
+    """General particle kernel error with optional custom message."""
 
     def __init__(self, particle, fieldset=None, msg=None):
         message = (f"{particle.state}\n"
@@ -103,7 +103,8 @@ def __init__(self, particle, fieldset=None, msg=None):
 
 
 class NotTestedError(Exception):
-    "Raised when the function is not tested. Code is commented"
+    """Raised when the function is not tested. Code is commented."""
+
     def __init__(self):
         super().__init__("Method has not been tested. Code has been commented out until unit tests are available.")
 
@@ -115,13 +116,13 @@ def parse_particletime(time, fieldset):
 
 
 def recovery_kernel_error(particle, fieldset, time):
-    """Default error kernel that throws exception"""
+    """Default error kernel that throws exception."""
     msg = f"Error: {particle.exception if particle.exception else None}"
     raise KernelError(particle, fieldset=fieldset, msg=msg)
 
 
 class InterpolationError(KernelError):
-    """Particle kernel error for interpolation error"""
+    """Particle kernel error for interpolation error."""
 
     def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
         if lon and lat:
@@ -132,7 +133,7 @@ def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
 
 
 class OutOfBoundsError(KernelError):
-    """Particle kernel error for out-of-bounds field sampling"""
+    """Particle kernel error for out-of-bounds field sampling."""
 
     def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
         if lon and lat:
@@ -143,7 +144,7 @@ def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
 
 
 class ThroughSurfaceError(KernelError):
-    """Particle kernel error for field sampling at surface"""
+    """Particle kernel error for field sampling at surface."""
 
     def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
         if lon and lat:
@@ -154,7 +155,7 @@ def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
 
 
 class OutOfTimeError(KernelError):
-    """Particle kernel error for time extrapolation field sampling"""
+    """Particle kernel error for time extrapolation field sampling."""
 
     def __init__(self, particle, fieldset):
         message = (f"Field sampled outside time domain at time {parse_particletime(particle.time, fieldset)}."
@@ -163,7 +164,7 @@ def __init__(self, particle, fieldset):
 
 
 def recovery_kernel_interpolation(particle, fieldset, time):
-    """Default sampling error kernel that throws InterpolationError"""
+    """Default sampling error kernel that throws InterpolationError."""
     if particle.exception is None:
         # TODO: JIT does not yet provide the context that created
         # the exception. We need to pass that info back from C.
@@ -174,7 +175,7 @@ def recovery_kernel_interpolation(particle, fieldset, time):
 
 
 def recovery_kernel_out_of_bounds(particle, fieldset, time):
-    """Default sampling error kernel that throws OutOfBoundsError"""
+    """Default sampling error kernel that throws OutOfBoundsError."""
     if particle.exception is None:
         # TODO: JIT does not yet provide the context that created
         # the exception. We need to pass that info back from C.
@@ -185,7 +186,7 @@ def recovery_kernel_out_of_bounds(particle, fieldset, time):
 
 
 def recovery_kernel_through_surface(particle, fieldset, time):
-    """Default sampling error kernel that throws OutOfBoundsError"""
+    """Default sampling error kernel that throws OutOfBoundsError."""
     if particle.exception is None:
         # TODO: JIT does not yet provide the context that created
         # the exception. We need to pass that info back from C.
@@ -196,7 +197,7 @@ def recovery_kernel_through_surface(particle, fieldset, time):
 
 
 def recovery_kernel_time_extrapolation(particle, fieldset, time):
-    """Default sampling error kernel that throws OutOfTimeError"""
+    """Default sampling error kernel that throws OutOfTimeError."""
     raise OutOfTimeError(particle, fieldset)
 
 
diff --git a/setup.cfg b/setup.cfg
index e180298eb..db7013414 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -18,7 +18,39 @@ classifiers =
     Operating System :: OS Independent
 
 [flake8]
-ignore = E501,F403,E226,E722,W503
+ignore =
+    # missing whitespace around arithmetic operator
+    E226,
+    # do not use bare except, specify exception instead
+    E722,
+    # line break before binary operator
+    W503
+    # line too long (82 > 79 characters)
+    E501,
+    # ‘from module import *’ used; unable to detect undefined names
+    F403,
+
+    # IGNORE FOR NOW (requires more work)
+    # Missing docstring in public module
+    D100,
+    # Missing docstring in public class
+    D101,
+    # Missing docstring in public method
+    D102,
+    # Missing docstring in public function
+    D103,
+    # Missing docstring in public package
+    D104,
+    # Missing docstring in magic method
+    D105,
+    # Missing docstring in __init__
+    D107,
+    # First line should end with a period (requires writing of summaries)
+    D400,
+    # First line should be in imperative mood (requires writing of summaries)
+    D401,
+    # 1 blank line required between summary line and description (requires writing of summaries)
+    D205,
 
 [tool:pytest]
 python_files = test_*.py example_*.py *tutorial*


Test patch:
diff --git a/tests/test_advection.py b/tests/test_advection.py
index b73fbe578..35802b7e4 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -73,9 +73,7 @@ def depth_fixture(zdim=2):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_advection_zonal(lon, lat, depth, pset_mode, mode, npart=10):
-    """ Particles at high latitude move geographically faster due to
-        the pole correction in `GeographicPolar`.
-    """
+    """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`."""
     data2D = {'U': np.ones((lon.size, lat.size), dtype=np.float32),
               'V': np.zeros((lon.size, lat.size), dtype=np.float32)}
     data3D = {'U': np.ones((lon.size, lat.size, depth.size), dtype=np.float32),
@@ -103,9 +101,7 @@ def test_advection_zonal(lon, lat, depth, pset_mode, mode, npart=10):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_advection_meridional(lon, lat, pset_mode, mode, npart=10):
-    """ Particles at high latitude move geographically faster due to
-        the pole correction in `GeographicPolar`.
-    """
+    """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`."""
     data = {'U': np.zeros((lon.size, lat.size), dtype=np.float32),
             'V': np.ones((lon.size, lat.size), dtype=np.float32)}
     dimensions = {'lon': lon, 'lat': lat}
@@ -122,8 +118,7 @@ def test_advection_meridional(lon, lat, pset_mode, mode, npart=10):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['jit', 'scipy'])
 def test_advection_3D(pset_mode, mode, npart=11):
-    """ 'Flat' 2D zonal flow that increases linearly with depth from 0 m/s to 1 m/s
-    """
+    """Flat 2D zonal flow that increases linearly with depth from 0 m/s to 1 m/s."""
     xdim = ydim = zdim = 2
     dimensions = {'lon': np.linspace(0., 1e4, xdim, dtype=np.float32),
                   'lat': np.linspace(0., 1e4, ydim, dtype=np.float32),
diff --git a/tests/test_data/create_testfields.py b/tests/test_data/create_testfields.py
index 91896cbe4..1eaba9430 100644
--- a/tests/test_data/create_testfields.py
+++ b/tests/test_data/create_testfields.py
@@ -74,7 +74,8 @@ def write_simple_2Dt(field, filename, varname=None):
     """Write a :class:`Field` to a netcdf file
 
     :param filename: Basename of the file
-    :param varname: Name of the field, to be appended to the filename"""
+    :param varname: Name of the field, to be appended to the filename.
+    """
     filepath = str(f'{filename}{field.name}.nc')
     if varname is None:
         varname = field.name
diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py
index 907c7172c..ce8ab4536 100644
--- a/tests/test_diffusion.py
+++ b/tests/test_diffusion.py
@@ -30,7 +30,7 @@
 
 
 def zeros_fieldset(mesh='spherical', xdim=200, ydim=100, mesh_conversion=1):
-    """Generates a zero velocity field"""
+    """Generates a zero velocity field."""
     lon = np.linspace(-1e5*mesh_conversion, 1e5*mesh_conversion, xdim, dtype=np.float32)
     lat = np.linspace(-1e5*mesh_conversion, 1e5*mesh_conversion, ydim, dtype=np.float32)
 
@@ -76,8 +76,7 @@ def test_fieldKh_Brownian(mesh, mode, pset_mode, xdim=200, ydim=100, kh_zonal=10
 @pytest.mark.parametrize('kernel', [AdvectionDiffusionM1,
                                     AdvectionDiffusionEM])
 def test_fieldKh_SpatiallyVaryingDiffusion(mesh, mode, pset_mode, kernel, xdim=200, ydim=100):
-    """Test advection-diffusion kernels on a non-uniform diffusivity field
-    with a linear gradient in one direction"""
+    """Test advection-diffusion kernels on a non-uniform diffusivity field with a linear gradient in one direction."""
     mesh_conversion = 1/1852./60 if mesh == 'spherical' else 1
     fieldset = zeros_fieldset(mesh=mesh, xdim=xdim, ydim=ydim, mesh_conversion=mesh_conversion)
 
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index ba99cfb9a..c55c60be5 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -66,7 +66,7 @@ def generate_fieldset(xdim, ydim, zdim=1, tdim=1):
 @pytest.mark.parametrize('xdim', [100, 200])
 @pytest.mark.parametrize('ydim', [100, 200])
 def test_fieldset_from_data(xdim, ydim):
-    """ Simple test for fieldset initialisation from data. """
+    """Simple test for fieldset initialisation from data."""
     data, dimensions = generate_fieldset(xdim, ydim)
     fieldset = FieldSet.from_data(data, dimensions)
     assert len(fieldset.U.data.shape) == 3
@@ -76,7 +76,7 @@ def test_fieldset_from_data(xdim, ydim):
 
 
 def test_fieldset_extra_syntax():
-    """ Simple test for fieldset initialisation from data. """
+    """Simple test for fieldset initialisation from data."""
     data, dimensions = generate_fieldset(10, 10)
     failed = False
     try:
@@ -109,9 +109,7 @@ def test_fieldset_from_data_timedims(ttype, tdim):
 @pytest.mark.parametrize('xdim', [100, 200])
 @pytest.mark.parametrize('ydim', [100, 50])
 def test_fieldset_from_data_different_dimensions(xdim, ydim, zdim=4, tdim=2):
-    """ Test for fieldset initialisation from data using
-    dict-of-dict for dimensions. """
-
+    """Test for fieldset initialisation from data using dict-of-dict for dimensions."""
     lon = np.linspace(0., 1., xdim, dtype=np.float32)
     lat = np.linspace(0., 1., ydim, dtype=np.float32)
     depth = np.zeros(zdim, dtype=np.float32)
@@ -137,7 +135,7 @@ def test_fieldset_from_data_different_dimensions(xdim, ydim, zdim=4, tdim=2):
 @pytest.mark.parametrize('xdim', [100, 200])
 @pytest.mark.parametrize('ydim', [100, 200])
 def test_fieldset_from_parcels(xdim, ydim, tmpdir, filename='test_parcels'):
-    """ Simple test for fieldset initialisation from Parcels FieldSet file format. """
+    """Simple test for fieldset initialisation from Parcels FieldSet file format."""
     filepath = tmpdir.join(filename)
     data, dimensions = generate_fieldset(xdim, ydim)
     fieldset_out = FieldSet.from_data(data, dimensions)
@@ -262,7 +260,7 @@ def test_fieldset_float64(cast_data_dtype, mode, tmpdir, xdim=10, ydim=5):
 @pytest.mark.parametrize('indslon', [range(10, 20), [1]])
 @pytest.mark.parametrize('indslat', [range(30, 60), [22]])
 def test_fieldset_from_file_subsets(indslon, indslat, tmpdir, filename='test_subsets'):
-    """ Test for subsetting fieldset from file using indices dict. """
+    """Test for subsetting fieldset from file using indices dict."""
     data, dimensions = generate_fieldset(100, 100)
     filepath = tmpdir.join(filename)
     fieldsetfull = FieldSet.from_data(data, dimensions)
@@ -367,7 +365,7 @@ def test_add_field_after_pset(pset_mode, fieldtype):
 
 @pytest.mark.parametrize('chunksize', ['auto', None])
 def test_fieldset_samegrids_from_file(tmpdir, chunksize, filename='test_subsets'):
-    """ Test for subsetting fieldset from file using indices dict. """
+    """Test for subsetting fieldset from file using indices dict."""
     data, dimensions = generate_fieldset(100, 100)
     filepath1 = tmpdir.join(filename+'_1')
     fieldset1 = FieldSet.from_data(data, dimensions)
@@ -407,7 +405,7 @@ def test_fieldset_dimlength1_cgrid(gridtype):
 
 @pytest.mark.parametrize('chunksize', ['auto', None])
 def test_fieldset_diffgrids_from_file(tmpdir, chunksize, filename='test_subsets'):
-    """ Test for subsetting fieldset from file using indices dict. """
+    """Test for subsetting fieldset from file using indices dict."""
     data, dimensions = generate_fieldset(100, 100)
     filepath1 = tmpdir.join(filename+'_1')
     fieldset1 = FieldSet.from_data(data, dimensions)
@@ -432,7 +430,7 @@ def test_fieldset_diffgrids_from_file(tmpdir, chunksize, filename='test_subsets'
 
 @pytest.mark.parametrize('chunksize', ['auto', None])
 def test_fieldset_diffgrids_from_file_data(tmpdir, chunksize, filename='test_subsets'):
-    """ Test for subsetting fieldset from file using indices dict. """
+    """Test for subsetting fieldset from file using indices dict."""
     data, dimensions = generate_fieldset(100, 100)
     filepath = tmpdir.join(filename)
     fieldset_data = FieldSet.from_data(data, dimensions)
@@ -460,7 +458,7 @@ def test_fieldset_diffgrids_from_file_data(tmpdir, chunksize, filename='test_sub
 
 
 def test_fieldset_samegrids_from_data(tmpdir, filename='test_subsets'):
-    """ Test for subsetting fieldset from file using indices dict. """
+    """Test for subsetting fieldset from file using indices dict."""
     data, dimensions = generate_fieldset(100, 100)
     fieldset1 = FieldSet.from_data(data, dimensions)
     field_data = fieldset1.U
@@ -1020,7 +1018,7 @@ def test_fieldset_frompop(pset_mode, mode):
 
 @pytest.mark.parametrize('pset_mode', pset_modes)
 def test_fieldset_from_data_gridtypes(pset_mode, xdim=20, ydim=10, zdim=4):
-    """ Simple test for fieldset initialisation from data. """
+    """Simple test for fieldset initialisation from data."""
     lon = np.linspace(0., 10., xdim, dtype=np.float32)
     lat = np.linspace(0., 10., ydim, dtype=np.float32)
     depth = np.linspace(0., 1., zdim, dtype=np.float32)
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index d44f6aa31..5dd9083c9 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -71,9 +71,7 @@ def k_sample_P_fixture():
 
 
 def fieldset(xdim=200, ydim=100):
-    """ Standard fieldset spanning the earth's coordinates with U and V
-        equivalent to longitude and latitude in deg.
-    """
+    """Standard fieldset spanning the earth's coordinates with U and V equivalent to longitude and latitude in deg."""
     lon = np.linspace(-180, 180, xdim, dtype=np.float32)
     lat = np.linspace(-90, 90, ydim, dtype=np.float32)
     U, V = np.meshgrid(lat, lon)
@@ -89,7 +87,7 @@ def fieldset_fixture(xdim=200, ydim=100):
 
 
 def fieldset_geometric(xdim=200, ydim=100):
-    """ Standard earth fieldset with U and V equivalent to lon/lat in m. """
+    """Standard earth fieldset with U and V equivalent to lon/lat in m."""
     lon = np.linspace(-180, 180, xdim, dtype=np.float32)
     lat = np.linspace(-90, 90, ydim, dtype=np.float32)
     U, V = np.meshgrid(lat, lon)
@@ -109,8 +107,8 @@ def fieldset_geometric_fixture(xdim=200, ydim=100):
 
 
 def fieldset_geometric_polar(xdim=200, ydim=100):
-    """ Standard earth fieldset with U and V equivalent to lon/lat in m
-        and the inversion of the pole correction applied to U.
+    """Standard earth fieldset with U and V equivalent to lon/lat in m
+    and the inversion of the pole correction applied to U.
     """
     lon = np.linspace(-180, 180, xdim, dtype=np.float32)
     lat = np.linspace(-90, 90, ydim, dtype=np.float32)
@@ -131,7 +129,7 @@ def fieldset_geometric_polar_fixture(xdim=200, ydim=100):
 
 
 def test_fieldset_sample(fieldset, xdim=120, ydim=80):
-    """ Sample the fieldset using indexing notation. """
+    """Sample the fieldset using indexing notation."""
     lon = np.linspace(-170, 170, xdim, dtype=np.float32)
     lat = np.linspace(-80, 80, ydim, dtype=np.float32)
     v_s = np.array([fieldset.UV[0, 0., 70., x][1] for x in lon])
@@ -141,7 +139,7 @@ def test_fieldset_sample(fieldset, xdim=120, ydim=80):
 
 
 def test_fieldset_sample_eval(fieldset, xdim=60, ydim=60):
-    """ Sample the fieldset using the explicit eval function. """
+    """Sample the fieldset using the explicit eval function."""
     lon = np.linspace(-170, 170, xdim, dtype=np.float32)
     lat = np.linspace(-80, 80, ydim, dtype=np.float32)
     v_s = np.array([fieldset.UV.eval(0, 0., 70., x)[1] for x in lon])
@@ -431,13 +429,12 @@ def test_partialslip_nearland_vertical(pset_mode, mode, boundaryslip, npart=20):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('lat_flip', [False, True])
 def test_fieldset_sample_particle(pset_mode, mode, k_sample_uv, lat_flip, npart=120):
-    """ Sample the fieldset using an array of particles.
+    """Sample the fieldset using an array of particles.
 
     Note that the low tolerances (1.e-6) are due to the first-order
     interpolation in JIT mode and give an indication of the
     corresponding sampling error.
     """
-
     lon = np.linspace(-180, 180, 200, dtype=np.float32)
     if lat_flip:
         lat = np.linspace(90, -90, 100, dtype=np.float32)
@@ -463,7 +460,7 @@ def test_fieldset_sample_particle(pset_mode, mode, k_sample_uv, lat_flip, npart=
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_fieldset_sample_geographic(fieldset_geometric, pset_mode, mode, k_sample_uv, npart=120):
-    """ Sample a fieldset with conversion to geographic units (degrees). """
+    """Sample a fieldset with conversion to geographic units (degrees)."""
     fieldset = fieldset_geometric
     lon = np.linspace(-170, 170, npart)
     lat = np.linspace(-80, 80, npart)
@@ -480,7 +477,7 @@ def test_fieldset_sample_geographic(fieldset_geometric, pset_mode, mode, k_sampl
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_fieldset_sample_geographic_noconvert(fieldset_geometric, pset_mode, mode, k_sample_uv_noconvert, npart=120):
-    """ Sample a fieldset without conversion to geographic units. """
+    """Sample a fieldset without conversion to geographic units."""
     fieldset = fieldset_geometric
     lon = np.linspace(-170, 170, npart)
     lat = np.linspace(-80, 80, npart)
@@ -497,7 +494,7 @@ def test_fieldset_sample_geographic_noconvert(fieldset_geometric, pset_mode, mod
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_fieldset_sample_geographic_polar(fieldset_geometric_polar, pset_mode, mode, k_sample_uv, npart=120):
-    """ Sample a fieldset with conversion to geographic units and a pole correction. """
+    """Sample a fieldset with conversion to geographic units and a pole correction."""
     fieldset = fieldset_geometric_polar
     lon = np.linspace(-170, 170, npart)
     lat = np.linspace(-80, 80, npart)
@@ -516,11 +513,10 @@ def test_fieldset_sample_geographic_polar(fieldset_geometric_polar, pset_mode, m
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_meridionalflow_spherical(pset_mode, mode, xdim=100, ydim=200):
-    """ Create uniform NORTHWARD flow on spherical earth and advect particles
+    """Create uniform NORTHWARD flow on spherical earth and advect particles.
 
-    As flow is so simple, it can be directly compared to analytical solution
+    As flow is so simple, it can be directly compared to analytical solution.
     """
-
     maxvel = 1.
     dimensions = {'lon': np.linspace(-180, 180, xdim, dtype=np.float32),
                   'lat': np.linspace(-90, 90, ydim, dtype=np.float32)}
@@ -544,7 +540,7 @@ def test_meridionalflow_spherical(pset_mode, mode, xdim=100, ydim=200):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_zonalflow_spherical(pset_mode, mode, k_sample_p, xdim=100, ydim=200):
-    """ Create uniform EASTWARD flow on spherical earth and advect particles
+    """Create uniform EASTWARD flow on spherical earth and advect particles.
 
     As flow is so simple, it can be directly compared to analytical solution
     Note that in this case the cosine conversion is needed
@@ -579,9 +575,7 @@ def test_zonalflow_spherical(pset_mode, mode, k_sample_p, xdim=100, ydim=200):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_random_field(pset_mode, mode, k_sample_p, xdim=20, ydim=20, npart=100):
-    """Sampling test that tests for overshoots by sampling a field of
-    random numbers between 0 and 1.
-    """
+    """Sampling test that tests for overshoots by sampling a field of random numbers between 0 and 1."""
     np.random.seed(123456)
     dimensions = {'lon': np.linspace(0., 1., xdim, dtype=np.float32),
                   'lat': np.linspace(0., 1., ydim, dtype=np.float32)}
diff --git a/tests/test_interaction_kernel.py b/tests/test_interaction_kernel.py
index e7a998f97..ad6346fcf 100644
--- a/tests/test_interaction_kernel.py
+++ b/tests/test_interaction_kernel.py
@@ -14,9 +14,7 @@
 
 
 def DummyMoveNeighbor(particle, fieldset, time, neighbors, mutator):
-    """A particle boosts the movement of its nearest neighbor, by adding
-    0.1 to its lat position.
-    """
+    """A particle boosts the movement of its nearest neighbor, by adding 0.1 to its lat position."""
     if len(neighbors) == 0:
         return StateCode.Success
 
@@ -37,7 +35,7 @@ def DoNothing(particle, fieldset, time):
 
 
 def fieldset(xdim=20, ydim=20, mesh='spherical'):
-    """ Standard unit mesh fieldset """
+    """Standard unit mesh fieldset."""
     lon = np.linspace(0., 1., xdim, dtype=np.float32)
     lat = np.linspace(0., 1., ydim, dtype=np.float32)
     U, V = np.meshgrid(lat, lon)
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index 1aee70990..7a34105ee 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -33,7 +33,7 @@ def DoNothing(particle, fieldset, time):
 
 
 def fieldset(xdim=20, ydim=20):
-    """ Standard unit mesh fieldset """
+    """Standard unit mesh fieldset."""
     lon = np.linspace(0., 1., xdim, dtype=np.float32)
     lat = np.linspace(0., 1., ydim, dtype=np.float32)
     U, V = np.meshgrid(lat, lon)
diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py
index 98b8b55e4..c19934fbc 100644
--- a/tests/test_kernel_language.py
+++ b/tests/test_kernel_language.py
@@ -41,7 +41,7 @@ def expr_kernel(name, pset, expr, pset_mode):
 
 
 def fieldset(xdim=20, ydim=20):
-    """ Standard unit mesh fieldset """
+    """Standard unit mesh fieldset."""
     lon = np.linspace(0., 1., xdim, dtype=np.float32)
     lat = np.linspace(0., 1., ydim, dtype=np.float32)
     U, V = np.meshgrid(lat, lon)
@@ -64,7 +64,7 @@ def fieldset_fixture(xdim=20, ydim=20):
     ('Div', '24 / 4', 6),
 ])
 def test_expression_int(pset_mode, mode, name, expr, result, npart=10):
-    """ Test basic arithmetic expressions """
+    """Test basic arithmetic expressions."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = pset_type[pset_mode]['pset'](None, pclass=TestParticle,
@@ -84,7 +84,7 @@ class TestParticle(ptype[mode]):
     ('Pow', '2 ** 3', 8),
 ])
 def test_expression_float(pset_mode, mode, name, expr, result, npart=10):
-    """ Test basic arithmetic expressions """
+    """Test basic arithmetic expressions."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = pset_type[pset_mode]['pset'](None, pclass=TestParticle,
@@ -110,7 +110,7 @@ class TestParticle(ptype[mode]):
     ('CheckNaN', 'math.nan != math.nan', True),
 ])
 def test_expression_bool(pset_mode, mode, name, expr, result, npart=10):
-    """ Test basic arithmetic expressions """
+    """Test basic arithmetic expressions."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = pset_type[pset_mode]['pset'](None, pclass=TestParticle,
@@ -126,7 +126,7 @@ class TestParticle(ptype[mode]):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_while_if_break(pset_mode, mode):
-    """Test while, if and break commands"""
+    """Test while, if and break commands."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32, initial=0.)
     pset = pset_type[pset_mode]['pset'](pclass=TestParticle, lon=[0], lat=[0])
@@ -145,7 +145,7 @@ def kernel(particle, fieldset, time):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_nested_if(pset_mode, mode):
-    """Test nested if commands"""
+    """Test nested if commands."""
     class TestParticle(ptype[mode]):
         p0 = Variable('p0', dtype=np.int32, initial=0)
         p1 = Variable('p1', dtype=np.int32, initial=1)
@@ -164,7 +164,7 @@ def kernel(particle, fieldset, time):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pass(pset_mode, mode):
-    """Test pass commands"""
+    """Test pass commands."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.int32, initial=0)
     pset = pset_type[pset_mode]['pset'](pclass=TestParticle, lon=0, lat=0)
@@ -190,7 +190,7 @@ def kernel(particle, fieldset, time):
 
 @pytest.mark.parametrize('pset_mode', pset_modes)
 def test_parcels_tmpvar_in_kernel(pset_mode):
-    """Tests for error thrown if variable with 'tmp' defined in custom kernel"""
+    """Tests for error thrown if variable with 'tmp' defined in custom kernel."""
     error_thrown = False
     pset = pset_type[pset_mode]['pset'](pclass=JITParticle, lon=0, lat=0)
 
@@ -210,7 +210,7 @@ def kernel_pnum(particle, fieldset, time):
 
 @pytest.mark.parametrize('pset_mode', pset_modes)
 def test_abs(pset_mode):
-    """Tests for error thrown if using abs in kernel"""
+    """Tests for error thrown if using abs in kernel."""
     error_thrown = False
     pset = pset_type[pset_mode]['pset'](pclass=JITParticle, lon=0, lat=0)
 
@@ -227,7 +227,7 @@ def kernel_abs(particle, fieldset, time):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_if_withfield(fieldset, pset_mode, mode):
-    """Test combination of if and Field sampling commands"""
+    """Test combination of if and Field sampling commands."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32, initial=0.)
     pset = pset_type[pset_mode]['pset'](fieldset, pclass=TestParticle, lon=[0], lat=[0])
@@ -268,7 +268,7 @@ def kernel(particle, fieldset, time):
                       reason="py.test FD capturing does not work for jit on python3 or Win"))
      ])
 def test_print(fieldset, pset_mode, mode, capfd):
-    """Test print statements"""
+    """Test print statements."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32, initial=0.)
     pset = pset_type[pset_mode]['pset'](fieldset, pclass=TestParticle, lon=[0.5], lat=[0.5])
@@ -329,7 +329,7 @@ def random_series(npart, rngfunc, rngargs, mode):
     ('randint', [0, 20]),
 ])
 def test_random_float(pset_mode, mode, rngfunc, rngargs, npart=10):
-    """ Test basic random number generation """
+    """Test basic random number generation."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32 if rngfunc == 'randint' else np.float32)
     pset = pset_type[pset_mode]['pset'](pclass=TestParticle,
diff --git a/tests/test_particles.py b/tests/test_particles.py
index 48d55505a..afad5e1c4 100644
--- a/tests/test_particles.py
+++ b/tests/test_particles.py
@@ -48,7 +48,7 @@ class TestParticle(ptype[mode]):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_variable_init(fieldset, pset_mode, mode, npart=10):
-    """Test that checks correct initialisation of custom variables"""
+    """Test that checks correct initialisation of custom variables."""
     class TestParticle(ptype[mode]):
         p_float = Variable('p_float', dtype=np.float32, initial=10.)
         p_double = Variable('p_double', dtype=np.float64, initial=11.)
@@ -72,7 +72,7 @@ def addOne(particle, fieldset, time):
 @pytest.mark.parametrize('mode', ['jit'])
 @pytest.mark.parametrize('type', ['np.int8', 'mp.float', 'np.int16'])
 def test_variable_unsupported_dtypes(fieldset, pset_mode, mode, type):
-    """Test that checks errors thrown for unsupported dtypes in JIT mode"""
+    """Test that checks errors thrown for unsupported dtypes in JIT mode."""
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=type, initial=10.)
     error_thrown = False
@@ -86,7 +86,7 @@ class TestParticle(ptype[mode]):
 @pytest.mark.parametrize('pset_mode', pset_modes)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_variable_special_names(fieldset, pset_mode, mode):
-    """Test that checks errors thrown for special names"""
+    """Test that checks errors thrown for special names."""
     for vars in ['z', 'lon']:
         class TestParticle(ptype[mode]):
             tmp = Variable(vars, dtype=np.float32, initial=10.)
@@ -102,7 +102,7 @@ class TestParticle(ptype[mode]):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('coord_type', [np.float32, np.float64])
 def test_variable_init_relative(fieldset, pset_mode, mode, coord_type, npart=10):
-    """Test that checks relative initialisation of custom variables"""
+    """Test that checks relative initialisation of custom variables."""
     lonlat_type = np.float64 if coord_type == 'double' else np.float32
 
     class TestParticle(ptype[mode]):