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 79c46106f..fca003d17 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -10,7 +10,10 @@ repos:
     rev: '6.0.0'
     hooks:
     -   id: flake8
-        additional_dependencies: ["flake8-docstrings"]
+-   repo: https://github.com/pycqa/pydocstyle
+    rev: 4.0.0  # pick a git hash / tag to point to
+    hooks:
+    -   id: pydocstyle
 -   repo: https://github.com/pycqa/isort
     rev: 5.12.0
     hooks:
diff --git a/docs/Makefile b/docs/Makefile
index 8db75bac0..580347a22 100755
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -5,7 +5,7 @@
 SPHINXOPTS    =
 SPHINXBUILD   = sphinx-build
 PAPER         =
-BUILDDIR      = ../../oceanparcels_website/gh-pages
+BUILDDIR      = _build
 
 # User-friendly check for sphinx-build
 ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1)
diff --git a/docs/conf.py b/docs/conf.py
index abbdd8f86..0b6013f07 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -33,6 +33,8 @@
     'sphinx.ext.autodoc',
     'sphinx.ext.todo',
     "sphinx.ext.linkcode",
+    "sphinx.ext.napoleon",
+    "numpydoc",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -98,7 +100,7 @@
 # show_authors = False
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = 'sphinx'
+pygments_style = 'monokai'
 
 # A list of ignored prefixes for module index sorting.
 # modindex_common_prefix = []
@@ -109,7 +111,6 @@
 # If true, `todo` and `todoList` produce output, else they produce nothing.
 todo_include_todos = True
 
-
 # -- Options for HTML output ----------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
@@ -137,6 +138,21 @@
 # pixels large.
 # html_favicon = None
 
+# numpydoc support
+# ----------------
+numpydoc_class_members_toctree = False  # https://stackoverflow.com/a/73294408
+
+# full list of numpydoc error codes: https://numpydoc.readthedocs.io/en/latest/validation.html
+numpydoc_validation_checks = {
+    "GL05",
+    "GL06",
+    "GL07",
+    "GL10",
+    "PR05",
+    "PR10",
+    "RT02",
+}
+
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
diff --git a/docs/index.rst b/docs/index.rst
index 8a3ec22ad..e185fa71d 100755
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -199,14 +199,14 @@ parcels.plotting module
     :undoc-members:
 
 parcels.scripts.plottrajectoriesfile module
------------------------------------
+-------------------------------------------
 
 .. automodule:: parcels.scripts.plottrajectoriesfile
     :members:
     :undoc-members:
 
 parcels.scripts.get_examples module
----------------------------
+-----------------------------------
 
 .. automodule:: parcels.scripts.get_examples
     :members:
diff --git a/docs/make.bat b/docs/make.bat
index 5076bd6f3..9a7cbd8fd 100755
--- a/docs/make.bat
+++ b/docs/make.bat
@@ -5,7 +5,7 @@ REM Command file for Sphinx documentation
 if "%SPHINXBUILD%" == "" (
 	set SPHINXBUILD=sphinx-build
 )
-set BUILDDIR=../../oceanparcels_website/gh-pages
+set BUILDDIR=_build
 set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% .
 set I18NSPHINXOPTS=%SPHINXOPTS% .
 if NOT "%PAPER%" == "" (
diff --git a/environment_py3_linux.yml b/environment_py3_linux.yml
index 6b178a6a4..986206f15 100644
--- a/environment_py3_linux.yml
+++ b/environment_py3_linux.yml
@@ -38,10 +38,11 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
-  - flake8-docstrings
   - pre_commit
+  - pydocstyle
 
   # Docs
+  - numpydoc
   - sphinx<6
   - pydata-sphinx-theme
   - sphinx-autobuild
diff --git a/environment_py3_osx.yml b/environment_py3_osx.yml
index 09287a2f6..3a1c656d3 100644
--- a/environment_py3_osx.yml
+++ b/environment_py3_osx.yml
@@ -38,10 +38,11 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
-  - flake8-docstrings
   - pre_commit
+  - pydocstyle
 
   # Docs
+  - numpydoc
   - sphinx<6
   - pydata-sphinx-theme
   - sphinx-autobuild
diff --git a/environment_py3_win.yml b/environment_py3_win.yml
index 795f9c25b..ce9e46f06 100644
--- a/environment_py3_win.yml
+++ b/environment_py3_win.yml
@@ -35,10 +35,11 @@ dependencies:
 
   # Linting
   - flake8>=2.1.0
-  - flake8-docstrings
   - pre_commit
+  - pydocstyle
 
   # Docs
+  - numpydoc
   - sphinx<6
   - pydata-sphinx-theme
   - sphinx-autobuild
diff --git a/parcels/application_kernels/EOSseawaterproperties.py b/parcels/application_kernels/EOSseawaterproperties.py
index 86e8142aa..50f6cf156 100644
--- a/parcels/application_kernels/EOSseawaterproperties.py
+++ b/parcels/application_kernels/EOSseawaterproperties.py
@@ -11,14 +11,13 @@ def PressureFromLatDepth(particle, fieldset, time):
     Returns
     -------
     p : array_like
-           pressure [db]
+        pressure [db]
 
     References
     ----------
-    .. [1] Saunders, Peter M., 1981: Practical Conversion of Pressure to Depth.
+    1. Saunders, Peter M., 1981: Practical Conversion of Pressure to Depth.
        J. Phys. Oceanogr., 11, 573-574.
        doi: 10.1175/1520-0485(1981)011<0573:PCOPTD>2.0.CO;2
-
     """
     # Angle conversions.
     deg2rad = math.pi / 180.0
@@ -29,31 +28,36 @@ def PressureFromLatDepth(particle, fieldset, time):
 
 
 def AdiabticTemperatureGradient(particle, fieldset, time):
-    """
-    Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
+    """Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
+
 
     Parameters
     ----------
-    :param particle.S: salinity [psu (PSS-78)]
-    :type particle.S: array_like
-    :param particle.T: temperature [℃ (ITS-90)]
-    :type particle.T: array_like
-    :param particle.pressure: pressure [db]
-    :type particle.pressure: array_like
+    particle :
+        The particle object with the following attributes:
+            - S : array_like
+                Salinity in psu (PSS-78).
+            - T : array_like
+                Temperature in ℃ (ITS-90).
+            - pressure : array_like
+                Pressure in db.
+    fieldset :
+        The fieldset object
 
     Returns
     -------
-    :return: adiabatic temperature gradient [℃ db :sup:`-1`]
-    :type: array_like
+    array_like
+        Adiabatic temperature gradient in ℃ db⁻¹.
+
 
     References
     ----------
-    .. [1] Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
+    1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
        computation of fundamental properties of seawater. UNESCO Tech. Pap. in
        Mar. Sci., No. 44, 53 pp.
        http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf
 
-    .. [2] Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
+    2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
        temperature gradient and potential temperature of sea water. Deep-Sea
        Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6
 
@@ -80,28 +84,35 @@ def PtempFromTemp(particle, fieldset, time):
 
     Parameters
     ----------
-    :param particle.S: salinity [psu (PSS-78)]
-    :type particle.S: array_like
-    :param particle.T: temperature [℃ (ITS-90)]
-    :type particle.T: array_like
-    :param particle.pressure: pressure [db]
-    :type particle.pressure: array_like
-    :param fieldset.refpressure: reference pressure [db], default = 0
-    :type fieldset.refpressure: array_like
+    particle :
+        The particle object with the following attributes:
+            - S : array_like
+                Salinity in psu (PSS-78).
+            - T : array_like
+                Temperature in ℃ (ITS-90).
+            - pressure : array_like
+                Pressure in db.
+    fieldset :
+        The fieldset object with the following attributes:
+            - refpressure : array_like, optional
+                Reference pressure in db (default is 0).
+    time : float
+        Simulation time (not used in this function but required for consistency with other kernels).
 
     Returns
     -------
-    :return: potential temperature relative to PR [℃ (ITS-90)]
-    :type: array_like
+    array_like
+        Potential temperature relative to reference pressure in ℃ (ITS-90).
+
 
     References
     ----------
-    .. [1] Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
+    1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
        computation of fundamental properties of seawater. UNESCO Tech. Pap. in
        Mar. Sci., No. 44, 53 pp.  Eqn.(31) p.39.
        http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf
 
-    .. [2] Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
+    2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
        temperature gradient and potential temperature of sea water. Deep-Sea
        Res. Vol20,401-408. doi:10.1016/0011-7471(73)90063-6
 
@@ -173,28 +184,34 @@ def TempFromPtemp(particle, fieldset, time):
 
     Parameters
     ----------
-    :param particle.S: salinity [psu (PSS-78)]
-    :type particle.S: array_like
-    :param particle.T: potential temperature [℃ (ITS-90)]
-    :type particle.T: array_like
-    :param particle.pressure: pressure [db]
-    :type particle.pressure: array_like
-    :param fieldset.refpressure: reference pressure [db], default = 0
-    :type fieldset.refpressure: array_like
+    particle :
+        The particle object with the following attributes:
+            - S : array_like
+                Salinity in psu (PSS-78).
+            - T : array_like
+                Potential temperature in ℃ (ITS-90).
+            - pressure : array_like
+                Pressure in db.
+    fieldset :
+        The fieldset object with the following attributes:
+            - refpressure : array_like, optional
+                Reference pressure in db (default is 0).
+    time : float
+        Simulation time (not used in this function but required for consistency with other kernels).
 
     Returns
     -------
-    :return: temperature [℃ (ITS-90)]
-    :type: array_like
+    array_like
+        Temperature in ℃ (ITS-90).
 
     References
     ----------
-    .. [1] Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
+    1. Fofonoff, P. and Millard, R.C. Jr UNESCO 1983. Algorithms for
        computation of fundamental properties of seawater. UNESCO Tech. Pap. in
        Mar. Sci., No. 44, 53 pp.  Eqn.(31) p.39.
        http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf
 
-    .. [2] Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
+    2. Bryden, H. 1973. New Polynomials for thermal expansion, adiabatic
        temperature gradient and potential temperature of sea water. Deep-Sea
        Res.  Vol20,401-408. doi:10.1016/0011-7471(73)90063-6
 
diff --git a/parcels/application_kernels/TEOSseawaterdensity.py b/parcels/application_kernels/TEOSseawaterdensity.py
index fa34abc31..18bff5880 100644
--- a/parcels/application_kernels/TEOSseawaterdensity.py
+++ b/parcels/application_kernels/TEOSseawaterdensity.py
@@ -5,20 +5,22 @@
 
 
 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
 
-    References:
-    Roquet, F., Madec, G., McDougall, T. J., Barker, P. M., 2014: Accurate
-    polynomial expressions for the density and specific volume of
-    seawater using the TEOS-10 standard. Ocean Modelling.
-    McDougall, T. J., D. R. Jackett, D. G. Wright and R. Feistel, 2003:
-    Accurate and computationally efficient algorithms for potential
-    temperature and density of seawater.  Journal of Atmospheric and
-    Oceanic Technology, 20, 730-741.
+    References
+    ----------
+    1. Roquet, F., Madec, G., McDougall, T. J., Barker, P. M., 2014: Accurate
+       polynomial expressions for the density and specific volume of
+       seawater using the TEOS-10 standard. Ocean Modelling.
+
+    2. McDougall, T. J., D. R. Jackett, D. G. Wright and R. Feistel, 2003:
+       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]
diff --git a/parcels/application_kernels/advectiondiffusion.py b/parcels/application_kernels/advectiondiffusion.py
index 26d57b0b9..c41919bdd 100644
--- a/parcels/application_kernels/advectiondiffusion.py
+++ b/parcels/application_kernels/advectiondiffusion.py
@@ -1,6 +1,6 @@
 """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
 
@@ -87,9 +87,9 @@ def DiffusionUniformKh(particle, fieldset, time):
 
     Assumes that fieldset has constant fields `Kh_zonal` and `Kh_meridional`.
     These can be added via e.g.
-        fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh)
-
-        fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh)
+    `fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh)`
+    or
+    `fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh)`
     where mesh is either 'flat' or 'spherical'
 
     This kernel assumes diffusivity gradients are zero and is therefore more efficient.
diff --git a/parcels/collection/collectionaos.py b/parcels/collection/collectionaos.py
index 513b748c3..e33dfd2cd 100644
--- a/parcels/collection/collectionaos.py
+++ b/parcels/collection/collectionaos.py
@@ -36,8 +36,11 @@ class ParticleCollectionAOS(ParticleCollection):
 
     def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, partitions=None, ngrid=1, **kwargs):
         """
-        :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
+        Parameters
+        ----------
+        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__()
 
@@ -180,7 +183,10 @@ def __getitem__(self, index):
         Access a particle in this collection using the fastest access
         method for this collection - by its index.
 
-        :param index: int or np.int32 index of a particle in this collection
+        Parameters
+        ----------
+        index : int
+            Index of the particle to access
         """
         return self.get_single_by_index(index)
 
@@ -189,7 +195,10 @@ def __getattr__(self, name):
         Access a single property of all particles.
         CAUTION: this function is not(!) in-place and is quite slow
 
-        :param name: name of the property
+        Parameters
+        ----------
+        name : str
+            Name of the property to access
         """
         pdtype = None
         for var in self._ptype.variables:
@@ -909,8 +918,7 @@ def setallvardata(self, var, val):
             setattr(self._data[i], var, val)
 
     def toArray(self):
-        """
-        This function converts (or: transforms; reformats; translates) this collection into an array-like structure
+        """This function converts (or: transforms; reformats; translates) this collection into an array-like structure
         (e.g. Python list or numpy nD array) that can be addressed by index. In the common case of 'no ID recovery',
         the global ID and the index match exactly.
 
@@ -924,10 +932,15 @@ def toArray(self):
         return self._data.tolist()
 
     def set_variable_write_status(self, var, write_status):
-        """
-        Method to set the write status of a Variable
-        :param var: Name of the variable (string)
-        :param status: Write status of the variable (True, False or 'once')
+        """Method to set the write status of a Variable
+
+        Parameters
+        ----------
+        var :
+            Name of the variable (string)
+        write_status :
+            Write status of the variable (True, False or 'once')
+
         """
         var_changed = False
         for v in self._ptype.variables:
@@ -948,11 +961,15 @@ class ParticleAccessorAOS(BaseParticleAccessor):
     """Wrapper that provides access to particle data in the collection,
     as if interacting with the particle itself.
 
-    :param pcoll: ParticleCollection that the represented particle
-                  belongs to.
-    :param index: The index at which the data for the represented
-                  particle is stored in the corresponding data arrays
-                  of the ParticleCollecion.
+    Parameters
+    ----------
+    pcoll :
+        ParticleCollection that the represented particle
+        belongs to.
+    index :
+        The index at which the data for the represented
+        particle is stored in the corresponding data arrays
+        of the ParticleCollecion.
     """
 
     _index = 0
@@ -967,11 +984,18 @@ def __init__(self, pcoll, index):
         self._next_dt = None
 
     def __getattr__(self, name):
-        """Get the value of an attribute of the particle.
+        """
+        Get the value of an attribute of the particle.
+
+        Parameters
+        ----------
+        name : str
+            Name of the requested particle attribute.
 
-        :param name: Name of the requested particle attribute.
-        :return: The value of the particle attribute in the underlying
-                 collection data array.
+        Returns
+        -------
+        any
+            The value of the particle attribute in the underlying collection data array.
         """
         if name in BaseParticleAccessor.__dict__.keys():
             result = super().__getattr__(name)
@@ -982,11 +1006,15 @@ def __getattr__(self, name):
         return result
 
     def __setattr__(self, name, value):
-        """Set the value of an attribute of the particle.
+        """
+        Set the value of an attribute of the particle.
 
-        :param name: Name of the particle attribute.
-        :param value: Value that will be assigned to the particle
-                      attribute in the underlying collection data array.
+        Parameters
+        ----------
+        name : str
+            Name of the particle attribute.
+        value : any
+            Value that will be assigned to the particle attribute in the underlying collection data array.
         """
         if name in BaseParticleAccessor.__dict__.keys():
             super().__setattr__(name, value)
@@ -1022,12 +1050,17 @@ def __iter__(self):
 class ParticleCollectionIteratorAOS(BaseParticleCollectionIterator):
     """Iterator for looping over the particles in the ParticleCollection.
 
-    :param pcoll: ParticleCollection that stores the particles.
-    :param reverse: Flag to indicate reverse iteration (i.e. starting at
-                    the largest index, instead of the smallest).
-    :param subset: Subset of indices to iterate over, this allows the
-                   creation of an iterator that represents part of the
-                   collection.
+    Parameters
+    ----------
+    pcoll :
+        ParticleCollection that stores the particles.
+    reverse :
+        Flag to indicate reverse iteration (i.e. starting at
+        the largest index, instead of the smallest).
+    subset :
+        Subset of indices to iterate over, this allows the
+        creation of an iterator that represents part of the
+        collection.
     """
 
     def __init__(self, pcoll, reverse=False, subset=None):
diff --git a/parcels/collection/collections.py b/parcels/collection/collections.py
index eef004af9..68cfde945 100644
--- a/parcels/collection/collections.py
+++ b/parcels/collection/collections.py
@@ -877,10 +877,17 @@ def cstruct(self):
     def __getattr__(self, name):
         """
         Access a single property of all particles.
-        NOTE: This is a fallback implementation, and it is NOT efficient.
-        Specific datastructures may implement a more efficient variant.
 
-        :param name: name of the property
+        Parameters
+        ----------
+        name : str
+            Name of the property to access
+
+
+        Notes
+        -----
+        This is a fallback implementation, and it is NOT efficient.
+        Specific datastructures may implement a more efficient variant.
         """
         for v in self.ptype.variables:
             if v.name == name:
@@ -908,10 +915,15 @@ def setvardata(self, var, index, val):
     def set_variable_write_status(self, var, write_status):
         """
         Method to set the write status of a Variable
-        :param var: Name of the variable (string)
-        :param status: Write status of the variable (True, False or 'once')
 
-         This function depends on the specific collection in question and thus needs to be specified in specific
+        This function depends on the specific collection in question and thus needs to be specified in specific
          derivatives classes.
+
+        Parameters
+        ----------
+        var : str
+            Name of the variable
+        write_status : bool, str
+            Write status of the variable (True, False or 'once')
         """
         pass
diff --git a/parcels/collection/collectionsoa.py b/parcels/collection/collectionsoa.py
index 65a234015..8665cdb39 100644
--- a/parcels/collection/collectionsoa.py
+++ b/parcels/collection/collectionsoa.py
@@ -31,8 +31,11 @@ class ParticleCollectionSOA(ParticleCollection):
 
     def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, partitions=None, ngrid=1, **kwargs):
         """
-        :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
+        Parameters
+        ----------
+        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__()
 
@@ -194,7 +197,10 @@ def __getitem__(self, index):
         Access a particle in this collection using the fastest access
         method for this collection - by its index.
 
-        :param index: int or np.int32 index of a particle in this collection
+        Parameters
+        ----------
+        index : int
+            Index of the particle to access
         """
         return self.get_single_by_index(index)
 
@@ -202,7 +208,10 @@ def __getattr__(self, name):
         """
         Access a single property of all particles.
 
-        :param name: name of the property
+        Parameters
+        ----------
+        name : str
+            Name of the property to access
         """
         for v in self.ptype.variables:
             if v.name == name and name in self._data:
@@ -341,9 +350,14 @@ def _recursive_ID_lookup(self, low, high, sublist):
         """Identify the middle element of the sublist and perform binary
         search on it.
 
-        :param low: Lowerbound on the indices to search for IDs.
-        :param high: Upperbound on the indices to search for IDs.
-        :param sublist: (Sub)list of IDs to look for.
+        Parameters
+        ----------
+        low :
+            Lowerbound on the indices to search for IDs.
+        high :
+            Upperbound on the indices to search for IDs.
+        sublist : list
+            Sublist of IDs to look for.
         """
         raise NotTestedError
         # median = floor(len(sublist) / 2)
@@ -832,10 +846,15 @@ def toArray(self):
         raise NotImplementedError
 
     def set_variable_write_status(self, var, write_status):
-        """
-        Method to set the write status of a Variable
-        :param var: Name of the variable (string)
-        :param status: Write status of the variable (True, False or 'once')
+        """Method to set the write status of a Variable
+
+        Parameters
+        ----------
+        var :
+            Name of the variable (string)
+        status :
+            Write status of the variable (True, False or 'once')
+        write_status :
         """
         var_changed = False
         for v in self._ptype.variables:
@@ -850,11 +869,15 @@ class ParticleAccessorSOA(BaseParticleAccessor):
     """Wrapper that provides access to particle data in the collection,
     as if interacting with the particle itself.
 
-    :param pcoll: ParticleCollection that the represented particle
-                  belongs to.
-    :param index: The index at which the data for the represented
-                  particle is stored in the corresponding data arrays
-                  of the ParticleCollecion.
+    Parameters
+    ----------
+    pcoll :
+        ParticleCollection that the represented particle
+        belongs to.
+    index :
+        The index at which the data for the represented
+        particle is stored in the corresponding data arrays
+        of the ParticleCollecion.
     """
 
     _index = 0
@@ -869,11 +892,18 @@ def __init__(self, pcoll, index):
         self._next_dt = None
 
     def __getattr__(self, name):
-        """Get the value of an attribute of the particle.
+        """
+        Get the value of an attribute of the particle.
 
-        :param name: Name of the requested particle attribute.
-        :return: The value of the particle attribute in the underlying
-                 collection data array.
+        Parameters
+        ----------
+        name : str
+            Name of the requested particle attribute.
+
+        Returns
+        -------
+        any
+            The value of the particle attribute in the underlying collection data array.
         """
         if name in BaseParticleAccessor.__dict__.keys():
             result = super().__getattr__(name)
@@ -884,11 +914,15 @@ def __getattr__(self, name):
         return result
 
     def __setattr__(self, name, value):
-        """Set the value of an attribute of the particle.
+        """
+        Set the value of an attribute of the particle.
 
-        :param name: Name of the particle attribute.
-        :param value: Value that will be assigned to the particle
-                      attribute in the underlying collection data array.
+        Parameters
+        ----------
+        name : str
+            Name of the particle attribute.
+        value : any
+            Value that will be assigned to the particle attribute in the underlying collection data array.
         """
         if name in BaseParticleAccessor.__dict__.keys():
             super().__setattr__(name, value)
@@ -937,12 +971,17 @@ def __getitem__(self, items):
 class ParticleCollectionIteratorSOA(BaseParticleCollectionIterator):
     """Iterator for looping over the particles in the ParticleCollection.
 
-    :param pcoll: ParticleCollection that stores the particles.
-    :param reverse: Flag to indicate reverse iteration (i.e. starting at
-                    the largest index, instead of the smallest).
-    :param subset: Subset of indices to iterate over, this allows the
-                   creation of an iterator that represents part of the
-                   collection.
+    Parameters
+    ----------
+    pcoll :
+        ParticleCollection that stores the particles.
+    reverse :
+        Flag to indicate reverse iteration (i.e. starting at
+        the largest index, instead of the smallest).
+    subset :
+        Subset of indices to iterate over, this allows the
+        creation of an iterator that represents part of the
+        collection.
     """
 
     def __init__(self, pcoll, reverse=False, subset=None):
diff --git a/parcels/compilation/codecompiler.py b/parcels/compilation/codecompiler.py
index 4afb7feae..b8f57eb13 100644
--- a/parcels/compilation/codecompiler.py
+++ b/parcels/compilation/codecompiler.py
@@ -197,9 +197,14 @@ def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=N
 class CCompiler:
     """A compiler object for creating and loading shared libraries.
 
-    :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).
+    Parameters
+    ----------
+    cc :
+        C compiler executable (uses environment variable ``CC`` if not provided).
+    cppargs :
+        A list of arguments to the C compiler (optional).
+    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()):
@@ -268,9 +273,13 @@ def compile(self, src, obj, log):
 class GNUCompiler_SS(CCompiler_SS):
     """A compiler object for the GNU Linux toolchain.
 
-    :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).
+    Parameters
+    ----------
+    cppargs :
+        A list of arguments to pass to the C compiler
+        (optional).
+    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()):
diff --git a/parcels/examples/example_moving_eddies.py b/parcels/examples/example_moving_eddies.py
index 2ad8ec1c6..fd874de36 100644
--- a/parcels/examples/example_moving_eddies.py
+++ b/parcels/examples/example_moving_eddies.py
@@ -25,17 +25,28 @@ def moving_eddies_fieldset(xdim=200, ydim=350, mesh='flat'):
     """Generate a fieldset encapsulating the flow field consisting of two
     moving eddies, one moving westward and the other moving northwestward.
 
-    :param xdim: Horizontal dimension of the generated fieldset
-    :param xdim: Vertical dimension of the generated fieldset
-    :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical: Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat  (default): No conversion, lat/lon are assumed to be in m.
-
+    Parameters
+    ----------
+    xdim :
+        Horizontal dimension of the generated fieldset (Default value = 200)
+    xdim :
+        Vertical dimension of the generated fieldset (Default value = 200)
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical: Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat  (default): No conversion, lat/lon are assumed to be in m.
+    ydim :
+         (Default value = 350)
+
+
+    Notes
+    -----
     Note that this is not a proper geophysical flow. Rather, a Gaussian eddy is moved
     artificially with uniform velocities. Velocities are calculated from geostrophy.
+
     """
     # Set Parcels FieldSet variables
     time = np.arange(0., 8. * 86400., 86400., dtype=np.float64)
@@ -93,8 +104,21 @@ def moving_eddies_example(fieldset, outfile, npart=2, mode='jit', verbose=False,
                           method=AdvectionRK4):
     """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.
+
+    Parameters
+    ----------
+    fieldset :
+        :class FieldSet: that defines the flow field
+    outfile :
+
+    npart :
+         Number of particles to initialise. (Default value = 2)
+    mode :
+         (Default value = 'jit')
+    verbose :
+         (Default value = False)
+    method :
+         (Default value = AdvectionRK4)
     """
     # Determine particle class according to mode
     start = (3.3, 46.) if fieldset.U.grid.mesh == 'spherical' else (3.3e5, 1e5)
diff --git a/parcels/examples/example_peninsula.py b/parcels/examples/example_peninsula.py
index 732d6f1bd..8e86bd3c9 100644
--- a/parcels/examples/example_peninsula.py
+++ b/parcels/examples/example_peninsula.py
@@ -25,22 +25,30 @@
 def peninsula_fieldset(xdim, ydim, mesh='flat', grid_type='A'):
     """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
-    :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical: Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat  (default): No conversion, lat/lon are assumed to be in m.
-    :param grid_type: Option whether grid is either Arakawa A (default) or C
-
-    The original test description can be found in Fig. 2.2.3 in:
-    North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of
-    recommended practices for modelling physical - biological
-    interactions during fish early life.
-    ICES Cooperative Research Report No. 295. 111 pp.
-    http://archimer.ifremer.fr/doc/00157/26792/24888.pdf
+    Parameters
+    ----------
+    xdim :
+        Horizontal dimension of the generated fieldset
+    xdim :
+        Vertical dimension of the generated fieldset
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical: Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat  (default): No conversion, lat/lon are assumed to be in m.
+    grid_type :
+        Option whether grid is either Arakawa A (default) or C
+
+        The original test description can be found in Fig. 2.2.3 in:
+        North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of
+        recommended practices for modelling physical - biological
+        interactions during fish early life.
+        ICES Cooperative Research Report No. 295. 111 pp.
+        http://archimer.ifremer.fr/doc/00157/26792/24888.pdf
+    ydim :
+
 
     """
     # Set Parcels FieldSet variables
@@ -97,8 +105,25 @@ 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.
+    Parameters
+    ----------
+    fieldset :
+
+    outfile : str
+        Basename of the input fieldset.
+    npart : int
+        Number of particles to intialise.
+    mode :
+         (Default value = 'jit')
+    degree :
+         (Default value = 1)
+    verbose :
+         (Default value = False)
+    output :
+         (Default value = True)
+    method :
+         (Default value = AdvectionRK4)
+
     """
     # First, we define a custom Particle class to which we add a
     # custom variable, the initial stream function value p.
diff --git a/parcels/field.py b/parcels/field.py
index 88f074d52..60087c755 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -45,53 +45,77 @@ def _isParticle(key):
 class Field:
     """Class that encapsulates access to field data.
 
-    :param name: Name of the field
-    :param data: 2D, 3D or 4D numpy array of field data.
-
-           1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
-              whichever is relevant for the dataset, use the flag transpose=True
-           2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
-              use the flag transpose=False
-           3. If data has any other shape, you first need to reorder it
-    :param lon: Longitude coordinates (numpy vector or array) of the field (only if grid is None)
-    :param lat: Latitude coordinates (numpy vector or array) of the field (only if grid is None)
-    :param depth: Depth coordinates (numpy vector or array) of the field (only if grid is None)
-    :param time: Time coordinates (numpy vector) of the field (only if grid is None)
-    :param mesh: String indicating the type of mesh coordinates and
-           units used during velocity interpolation: (only if grid is None)
-
-           1. spherical: Lat and lon in degree, with a
-              correction for zonal velocity U near the poles.
-           2. flat (default): No conversion, lat/lon are assumed to be in m.
-    :param timestamps: A numpy array containing the timestamps for each of the files in filenames, for loading
-           from netCDF files only. Default is None if the netCDF dimensions dictionary includes time.
-    :param grid: :class:`parcels.grid.Grid` object containing all the lon, lat depth, time
-           mesh and time_origin information. Can be constructed from any of the Grid objects
-    :param fieldtype: Type of Field to be used for UnitConverter when using SummedFields
-           (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-    :param transpose: Transpose data to required (lon, lat) layout
-    :param vmin: Minimum allowed value on the field. Data below this value are set to zero
-    :param vmax: Maximum allowed value on the field. Data above this value are set to zero
-    :param cast_data_dtype: Cast Field data to dtype. Supported dtypes are np.float32 (default) and np.float64.
-           Note that dtype can only be float32 in JIT mode
-    :param time_origin: Time origin (TimeConverter object) of the time axis (only if grid is None)
-    :param interp_method: Method for interpolation. Options are 'linear' (default), 'nearest',
-           'linear_invdist_land_tracer', 'cgrid_velocity', 'cgrid_tracer' and 'bgrid_velocity'
-    :param allow_time_extrapolation: boolean whether to allow for extrapolation in time
-           (i.e. beyond the last available time snapshot)
-    :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).
-           The last value of the time series can be provided (which is the same as the initial one) or not (Default: False)
-           This flag overrides the allow_time_interpolation and sets it to False
-    :param chunkdims_name_map (opt.): gives a name map to the FieldFileBuffer that declared a mapping between chunksize name, NetCDF dimension and Parcels dimension;
-           required only if currently incompatible OCM field is loaded and chunking is used by 'chunksize' (which is the default)
-    :param to_write: Write the Field in NetCDF format at the same frequency as the ParticleFile outputdt,
-           using a filenaming scheme based on the ParticleFile name
-
+    Parameters
+    ----------
+    name : str
+        Name of the field
+    data : np.ndarray
+        2D, 3D or 4D numpy array of field data.
+
+        1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
+        whichever is relevant for the dataset, use the flag transpose=True
+        2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
+        use the flag transpose=False
+        3. If data has any other shape, you first need to reorder it
+    lon : np.ndarray or list
+        Longitude coordinates (numpy vector or array) of the field (only if grid is None)
+    lat : np.ndarray or list
+        Latitude coordinates (numpy vector or array) of the field (only if grid is None)
+    depth : np.ndarray or list
+        Depth coordinates (numpy vector or array) of the field (only if grid is None)
+    time : np.ndarray
+        Time coordinates (numpy vector) of the field (only if grid is None)
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation: (only if grid is None)
+
+        1. spherical: Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat (default): No conversion, lat/lon are assumed to be in m.
+    timestamps : np.ndarray
+        A numpy array containing the timestamps for each of the files in filenames, for loading
+        from netCDF files only. Default is None if the netCDF dimensions dictionary includes time.
+    grid : parcels.grid.Grid
+        :class:`parcels.grid.Grid` object containing all the lon, lat depth, time
+        mesh and time_origin information. Can be constructed from any of the Grid objects
+    fieldtype : str
+        Type of Field to be used for UnitConverter when using SummedFields
+        (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+    transpose : bool
+        Transpose data to required (lon, lat) layout
+    vmin : float
+        Minimum allowed value on the field. Data below this value are set to zero
+    vmax : float
+        Maximum allowed value on the field. Data above this value are set to zero
+    cast_data_dtype : str
+        Cast Field data to dtype. Supported dtypes are "float32" (np.float32 (default)) and "float64 (np.float64).
+        Note that dtype can only be "float32" in JIT mode
+    time_origin : parcels.tools.converters.TimeConverter
+        Time origin of the time axis (only if grid is None)
+    interp_method : str
+        Method for interpolation. Options are 'linear' (default), 'nearest',
+        'linear_invdist_land_tracer', 'cgrid_velocity', 'cgrid_tracer' and 'bgrid_velocity'
+    allow_time_extrapolation : bool
+        boolean whether to allow for extrapolation in time
+        (i.e. beyond the last available time snapshot)
+    time_periodic : bool, float or datetime.timedelta
+        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).
+        The last value of the time series can be provided (which is the same as the initial one) or not (Default: False)
+        This flag overrides the allow_time_extrapolation and sets it to False
+    chunkdims_name_map : str, optional
+        Gives a name map to the FieldFileBuffer that declared a mapping between chunksize name, NetCDF dimension and Parcels dimension;
+        required only if currently incompatible OCM field is loaded and chunking is used by 'chunksize' (which is the default)
+    to_write : bool
+        Write the Field in NetCDF format at the same frequency as the ParticleFile outputdt,
+        using a filenaming scheme based on the ParticleFile name
+
+    Examples
+    --------
     For usage examples see the following tutorials:
 
-    * `Nested Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`_
+    * `Nested Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`__
 
-    * `Summed Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_SummedFields.ipynb>`_
+    * `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,
@@ -270,41 +294,62 @@ def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
                     deferred_load=True, **kwargs):
         """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)
-               In the latter case, time values are in filenames[data]
-        :param variable: Tuple mapping field name to variable name in the NetCDF file.
-        :param dimensions: Dictionary mapping variable names for the relevant dimensions in the NetCDF file
-        :param indices: dictionary mapping indices for each dimension to read from file.
-               This can be used for reading in only a subregion of the NetCDF file.
-               Note that negative indices are not allowed.
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param timestamps: A numpy array of datetime64 objects containing the timestamps for each of the files in filenames.
-               Default is None if dimensions includes time.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation in time
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param deferred_load: boolean whether to only pre-load data (in deferred mode) or
-               fully load them (default: True). It is advised to deferred load the data, since in
-               that case Parcels deals with a better memory management during particle set execution.
-               deferred_load=False is however sometimes necessary for plotting the fields.
-        :param gridindexingtype: The type of gridindexing. Either 'nemo' (default) or 'mitgcm' are supported.
-               See also the Grid indexing documentation on oceanparcels.org
-        :param chunksize: size of the chunks in dask loading
-        :param netcdf_decodewarning: boolean whether to show a warning id there is a problem decoding the netcdf files.
-               Default is True, but in some cases where these warnings are expected, it may be useful to silence them
-               by setting netcdf_decodewarning=False.
-
+        Parameters
+        ----------
+        filenames : list of str or dict
+            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)
+            In the latter case, time values are in filenames[data]
+        variable : tuple of str or str
+            Tuple mapping field name to variable name in the NetCDF file.
+        dimensions : dict
+            Dictionary mapping variable names for the relevant dimensions in the NetCDF file
+        indices :
+            dictionary mapping indices for each dimension to read from file.
+            This can be used for reading in only a subregion of the NetCDF file.
+            Note that negative indices are not allowed. (Default value = None)
+        mesh :
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        timestamps :
+            A numpy array of datetime64 objects containing the timestamps for each of the files in filenames.
+            Default is None if dimensions includes time.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation in time
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            boolean whether to loop periodically over the time component of the FieldSet
+            This flag overrides the allow_time_extrapolation and sets it to False (Default value = False)
+        deferred_load : bool
+            boolean whether to only pre-load data (in deferred mode) or
+            fully load them (default: True). It is advised to deferred load the data, since in
+            that case Parcels deals with a better memory management during particle set execution.
+            deferred_load=False is however sometimes necessary for plotting the fields.
+        gridindexingtype : str
+            The type of gridindexing. Either 'nemo' (default) or 'mitgcm' are supported.
+            See also the Grid indexing documentation on oceanparcels.org
+        chunksize :
+            size of the chunks in dask loading
+        netcdf_decodewarning : bool
+            Whether to show a warning id there is a problem decoding the netcdf files.
+            Default is True, but in some cases where these warnings are expected, it may be useful to silence them
+            by setting netcdf_decodewarning=False.
+        grid :
+             (Default value = None)
+        **kwargs :
+            Keyword arguments passed to the :class:`Field` constructor.
+
+        Examples
+        --------
         For usage examples see the following tutorial:
 
-        * `Timestamps <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timestamps.ipynb>`_
+        * `Timestamps <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timestamps.ipynb>`__
+
         """
         # Ensure the timestamps array is compatible with the user-provided datafiles.
         if timestamps is not None:
@@ -489,20 +534,30 @@ def from_xarray(cls, da, name, dimensions, mesh='spherical', allow_time_extrapol
                     time_periodic=False, **kwargs):
         """Create field from xarray Variable.
 
-        :param da: Xarray DataArray
-        :param name: Name of the Field
-        :param dimensions: Dictionary mapping variable names for the relevant dimensions in the DataArray
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation in time
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
+        Parameters
+        ----------
+        da : xr.DataArray
+            Xarray DataArray
+        name : str
+            Name of the Field
+        dimensions : dict
+            Dictionary mapping variable names for the relevant dimensions in the DataArray
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation in time
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            boolean whether to loop periodically over the time component of the FieldSet
+            This flag overrides the allow_time_extrapolation and sets it to False (Default value = False)
+        **kwargs :
+            Keyword arguments passed to the :class:`Field` constructor.
         """
         data = da.data
         interp_method = kwargs.pop('interp_method', 'linear')
@@ -573,11 +628,17 @@ def reshape(self, data, transpose=False):
     def set_scaling_factor(self, factor):
         """Scales the field data by some constant factor.
 
-        :param factor: scaling factor
+        Parameters
+        ----------
+        factor :
+            scaling factor
 
+
+        Examples
+        --------
         For usage examples see the following tutorial:
 
-        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_
+        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__
         """
         if self._scaling_factor:
             raise NotImplementedError(f'Scaling factor for field {self.name} already defined.')
@@ -588,8 +649,11 @@ def set_scaling_factor(self, factor):
     def set_depth_from_field(self, 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>`_
+        Notes
+        -----
+        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.
+
         """
         self.grid.depth_field = field
         if self.grid != field.grid:
@@ -1043,10 +1107,12 @@ def interpolator3D(self, ti, z, y, x, time, particle=None):
     def temporal_interpolate_fullfield(self, ti, time):
         """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.
+        Parameters
+        ----------
+        ti :
+            Index in time array associated with time (via :func:`time_index`)
+        time :
+            Time to interpolate to
         """
         t0 = self.grid.time[ti]
         if time == t0:
@@ -1249,15 +1315,28 @@ def show(self, animation=False, show_time=None, domain=None, depth_level=0, proj
              vmin=None, vmax=None, savefile=None, **kwargs):
         """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)
-        :param domain: dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show
-        :param depth_level: depth level to be plotted (default 0)
-        :param projection: type of cartopy projection to use (default PlateCarree)
-        :param land: Boolean whether to show land. This is ignored for flat meshes
-        :param vmin: minimum colour scale (only in single-plot mode)
-        :param vmax: maximum colour scale (only in single-plot mode)
-        :param savefile: Name of a file to save the plot to
+        Parameters
+        ----------
+        animation : bool
+            Whether result is a single plot, or an animation (Default value = False)
+        show_time : float
+            Time in seconds from start after which to show the Field (only in single-plot mode) (Default value = None)
+        domain : dict
+            dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show (Default value = None)
+        depth_level :
+            depth level to be plotted (default 0)
+        projection :
+            type of cartopy projection to use (default PlateCarree)
+        land : bool
+            Whether to show land. This is ignored for flat meshes (Default value = True)
+        vmin : float
+            minimum colour scale (only in single-plot mode) (Default value = None)
+        vmax : float
+            maximum colour scale (only in single-plot mode) (Default value = None)
+        savefile : str, opptional
+            Name of a file to save the plot to (Default value = None)
+        **kwargs :
+            Additional keyword arguments to pass to :func:`parcels.plotting.plotfield`
         """
         from parcels.plotting import plotfield
         plt, _, _, _ = plotfield(self, animation=animation, show_time=show_time, domain=domain, depth_level=depth_level,
@@ -1272,13 +1351,19 @@ def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
         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
 
-        See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb>`_
+        See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb>`__
         for a detailed explanation on how to set up periodic boundaries
 
-        :param zonal: Create a halo in zonal direction (boolean)
-        :param meridional: Create a halo in meridional direction (boolean)
-        :param halosize: size of the halo (in grid points). Default is 5 grid points
-        :param data: if data is not None, the periodic halo will be achieved on data instead of self.data and data will be returned
+        Parameters
+        ----------
+        zonal : bool
+            Create a halo in zonal direction.
+        meridional : bool
+            Create a halo in meridional direction.
+        halosize : int
+            Size of the halo (in grid points). Default is 5 grid points
+        data :
+            if data is not None, the periodic halo will be achieved on data instead of self.data and data will be returned (Default value = None)
         """
         dataNone = not isinstance(data, (np.ndarray, da.core.Array))
         if self.grid.defer_load and dataNone:
@@ -1314,8 +1399,12 @@ def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
     def write(self, 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.
+        Parameters
+        ----------
+        filename : str
+            Basename of the file (i.e. '{filename}{Field.name}.nc')
+        varname : str
+            Name of the field, to be appended to the filename. (Default value = None)
         """
         filepath = str(Path(f'{filename}{self.name}.nc'))
         if varname is None:
@@ -1427,10 +1516,16 @@ class VectorField:
     """Class VectorField stores 2 or 3 fields which defines together a vector field.
     This enables to interpolate them as one single vector field in the kernels.
 
-    :param name: Name of the vector field
-    :param U: field defining the zonal component
-    :param V: field defining the meridional component
-    :param W: field defining the vertical component (default: None)
+    Parameters
+    ----------
+    name : str
+        Name of the vector field
+    U : parcels.field.Field
+        field defining the zonal component
+    V : parcels.field.Field
+        field defining the meridional component
+    W : parcels.field.Field
+        field defining the vertical component (default: None)
     """
 
     def __init__(self, name, U, V, W=None):
@@ -1644,8 +1739,7 @@ def spatial_c_grid_interpolation3D_full(self, ti, z, y, x, time, particle=None):
         return (u, v, w)
 
     def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, applyConversion=True):
-        """
-        Perform C grid interpolation in 3D.
+        """Perform C grid interpolation in 3D.
 
         +---+---+---+
         |   |V1 |   |
@@ -1863,13 +1957,24 @@ class SummedField(list):
     still be queried through their list index (e.g. SummedField[1]).
     SummedField is composed of either Fields or VectorFields.
 
-    See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_SummedFields.ipynb>`_
+
+    Parameters
+    ----------
+    name : str
+        Name of the SummedField
+    F : list of Field
+        List of fields. F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
+    V : list of Field
+        List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
+    W : list of Field
+        List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None)
+
+
+    Examples
+    --------
+    See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_SummedFields.ipynb>`__
     for a detailed tutorial
 
-    :param name: Name of the SummedField
-    :param F: List of fields. F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
-    :param V: List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
-    :param W: List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None)
     """
 
     def __init__(self, name, F, V=None, W=None):
@@ -1934,13 +2039,23 @@ class NestedField(list):
     NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries.
     NestedField is composed of either Fields or VectorFields.
 
-    See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`_
+    Parameters
+    ----------
+    name : str
+        Name of the NestedField
+    F : list of Field
+        List of fields (order matters). F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
+    V : list of Field
+        List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
+    W : list of Field
+        List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None)
+
+
+    Examples
+    --------
+    See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`__
     for a detailed tutorial
 
-    :param name: Name of the NestedField
-    :param F: List of fields (order matters). F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField
-    :param V: List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None)
-    :param W: List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None)
     """
 
     def __init__(self, name, F, V=None, W=None):
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index c0aea9071..a4cdc4481 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -25,9 +25,15 @@
 class FieldSet:
     """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
+    Parameters
+    ----------
+    U : parcels.field.Field
+        Field object for zonal velocity component
+    V : parcels.field.Field
+        Field object for meridional velocity component
+    fields : dict mapping str to Field
+        Additional fields to include in the FieldSet. These fields can be used
+        in custom kernels.
     """
 
     def __init__(self, U, V, fields=None):
@@ -58,44 +64,54 @@ def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
                   allow_time_extrapolation=None, time_periodic=False, **kwargs):
         """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
-               the built-in Advection kernels assume that U and V are in m/s
-
-               1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
-                  whichever is relevant for the dataset, use the flag transpose=True
-               2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
-                  use the flag transpose=False (default value)
-               3. If data has any other shape, you first need to reorder it
-        :param dimensions: Dictionary mapping field dimensions (lon,
-               lat, depth, time) to numpy arrays.
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable
-               (e.g. dimensions['U'], dimensions['V'], etc).
-        :param transpose: Boolean whether to transpose data on read-in
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-
-        Usage examples
-        ==============
-
-        * `Analytical advection <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_analyticaladvection.ipynb>`_
+        Parameters
+        ----------
+        data :
+            Dictionary mapping field names to numpy arrays.
+            Note that at least a 'U' and 'V' numpy array need to be given, and that
+            the built-in Advection kernels assume that U and V are in m/s
+
+            1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
+            whichever is relevant for the dataset, use the flag transpose=True
+            2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
+            use the flag transpose=False (default value)
+            3. If data has any other shape, you first need to reorder it
+        dimensions : dict
+            Dictionary mapping field dimensions (lon,
+            lat, depth, time) to numpy arrays.
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable
+            (e.g. dimensions['U'], dimensions['V'], etc).
+        transpose : bool
+            Whether to transpose data on read-in (Default value = False)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        **kwargs :
+            Keyword arguments passed to the :class:`Field` constructor.
+
+        Examples
+        --------
+        For usage examples see the following tutorials:
 
-        * `Diffusion <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`_
+        * `Analytical advection <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_analyticaladvection.ipynb>`__
 
-        * `Interpolation <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_interpolation.ipynb>`_
+        * `Diffusion <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`__
 
-        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_
+        * `Interpolation <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_interpolation.ipynb>`__
 
+        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__
         """
         fields = {}
         for name, datafld in data.items():
@@ -129,14 +145,21 @@ def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
     def add_field(self, field, name=None):
         """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
+        Parameters
+        ----------
+        field : parcels.field.Field
+            Field object to be added
+        name : str
+            Name of the :class:`parcels.field.Field` object to be added
+
 
+        Examples
+        --------
         For usage examples see the following tutorials:
 
-        * `Nested Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`_
+        * `Nested Fields <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_NestedFields.ipynb>`__
 
-        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_
+        * `Unit converters <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__ (Default value = None)
 
         """
         if self.completed:
@@ -166,17 +189,29 @@ def add_constant_field(self, name, value, mesh='flat'):
         """Wrapper function to add a Field that is constant in space,
            useful e.g. when using constant horizontal diffusivity
 
-        :param name: Name of the :class:`parcels.field.Field` object to be added
-        :param value: Value of the constant field (stored as 32-bit float)
-        :param units: Optional UnitConverter object, to convert units
-                      (e.g. for Horizontal diffusivity from m2/s to degree2/s)
+        Parameters
+        ----------
+        name : str
+            Name of the :class:`parcels.field.Field` object to be added
+        value : float
+            Value of the constant field (stored as 32-bit float)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
         """
         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.
 
-        :param vfield: :class:`parcels.field.VectorField` object to be added
+        Parameters
+        ----------
+        vfield : parcels.VectorField
+            class:`parcels.field.VectorField` object to be added
         """
         setattr(self, vfield.name, vfield)
         for v in vfield.__dict__.values():
@@ -293,64 +328,85 @@ def from_netcdf(cls, filenames, variables, dimensions, indices=None, fieldtype=N
                     deferred_load=True, chunksize=None, **kwargs):
         """Initialises FieldSet object from NetCDF files.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}.
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable names in the netCDF file(s).
-               Note that the built-in Advection kernels assume that U and V are in m/s
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable
-               (e.g. dimensions['U'], dimensions['V'], etc).
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also `this tuturial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param timestamps: list of lists or array of arrays containing the timestamps for
-               each of the files in filenames. Outer list/array corresponds to files, inner
-               array corresponds to indices within files.
-               Default is None if dimensions includes time.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param deferred_load: boolean whether to only pre-load data (in deferred mode) or
-               fully load them (default: True). It is advised to deferred load the data, since in
-               that case Parcels deals with a better memory management during particle set execution.
-               deferred_load=False is however sometimes necessary for plotting the fields.
-        :param interp_method: Method for interpolation. Options are 'linear' (default), 'nearest',
-               'linear_invdist_land_tracer', 'cgrid_velocity', 'cgrid_tracer' and 'bgrid_velocity'
-        :param gridindexingtype: The type of gridindexing. Either 'nemo' (default) or 'mitgcm' are supported.
-               See also the Grid indexing documentation on oceanparcels.org
-        :param chunksize: size of the chunks in dask loading. Default is None (no chunking). Can be None or False (no chunking),
-               'auto' (chunking is done in the background, but results in one grid per field individually), or a dict in the format
-               '{parcels_varname: {netcdf_dimname : (parcels_dimname, chunksize_as_int)}, ...}', where 'parcels_dimname' is one of ('time', 'depth', 'lat', 'lon')
-        :param netcdf_engine: engine to use for netcdf reading in xarray. Default is 'netcdf',
-               but in cases where this doesn't work, setting netcdf_engine='scipy' could help
-
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}.
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable names in the netCDF file(s).
+            Note that the built-in Advection kernels assume that U and V are in m/s
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable
+            (e.g. dimensions['U'], dimensions['V'], etc).
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None) (Default value = None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tuturial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        timestamps :
+            list of lists or array of arrays containing the timestamps for
+            each of the files in filenames. Outer list/array corresponds to files, inner
+            array corresponds to indices within files.
+            Default is None if dimensions includes time.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        deferred_load : bool
+            boolean whether to only pre-load data (in deferred mode) or
+            fully load them (default: True). It is advised to deferred load the data, since in
+            that case Parcels deals with a better memory management during particle set execution.
+            deferred_load=False is however sometimes necessary for plotting the fields.
+        interp_method : str
+            Method for interpolation. Options are 'linear' (default), 'nearest',
+            'linear_invdist_land_tracer', 'cgrid_velocity', 'cgrid_tracer' and 'bgrid_velocity'
+        gridindexingtype : str
+            The type of gridindexing. Either 'nemo' (default) or 'mitgcm' are supported.
+            See also the Grid indexing documentation on oceanparcels.org
+        chunksize :
+            size of the chunks in dask loading. Default is None (no chunking). Can be None or False (no chunking),
+            'auto' (chunking is done in the background, but results in one grid per field individually), or a dict in the format
+            '{parcels_varname: {netcdf_dimname : (parcels_dimname, chunksize_as_int)}, ...}', where 'parcels_dimname' is one of ('time', 'depth', 'lat', 'lon')
+        netcdf_engine :
+            engine to use for netcdf reading in xarray. Default is 'netcdf',
+            but in cases where this doesn't work, setting netcdf_engine='scipy' could help
+        **kwargs :
+            Keyword arguments passed to the :class:`parcels.Field` constructor.
+
+
+        Examples
+        --------
         For usage examples see the following tutorials:
 
-        * `Basic Parcels setup <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/parcels_tutorial.ipynb>`_
+        * `Basic Parcels setup <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/parcels_tutorial.ipynb>`__
 
-        * `Argo floats <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_Argofloats.ipynb>`_
+        * `Argo floats <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_Argofloats.ipynb>`__
 
-        * `Timestamps <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timestamps.ipynb>`_
+        * `Timestamps <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timestamps.ipynb>`__
 
-        * `Time-evolving depth dimensions <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timevaryingdepthdimensions.ipynb>`_
+        * `Time-evolving depth dimensions <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_timevaryingdepthdimensions.ipynb>`__
 
         """
         # Ensure that times are not provided both in netcdf file and in 'timestamps'.
@@ -426,62 +482,76 @@ def from_nemo(cls, filenames, variables, dimensions, indices=None, mesh='spheric
                   tracer_interp_method='cgrid_tracer', chunksize=None, **kwargs):
         """Initialises FieldSet object from NetCDF files of Curvilinear NEMO fields.
 
-        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_curvilinear.ipynb>`_
-        for a detailed tutorial on the setup for 2D NEMO fields and `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb>`_
+        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_curvilinear.ipynb>`__
+        for a detailed tutorial on the setup for 2D NEMO fields and `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb>`__
         for the tutorial on the setup for 3D NEMO fields.
 
-        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb>`_
+        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb>`__
         for a more detailed explanation of the different methods that can be used for c-grid datasets.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files,
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable names in the netCDF file(s).
-               Note that the built-in Advection kernels assume that U and V are in m/s
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable.
-               Watch out: NEMO is discretised on a C-grid:
-               U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).
-
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |         V[k,j+1,i+1]        |                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j+1,i]                   |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|U[k,j+1,i+1]                 |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |         V[k,j,i+1]          +                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-
-               To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes,
-               which are located on the corners of the cells.
-               (for indexing details: https://www.nemo-ocean.eu/doc/img360.png )
-               In 3D, the depth is the one corresponding to W nodes
-               The gridindexingtype is set to 'nemo'. See also the Grid indexing documentation on oceanparcels.org
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param tracer_interp_method: Method for interpolation of tracer fields. It is recommended to use 'cgrid_tracer' (default)
-               Note that in the case of from_nemo() and from_cgrid(), the velocity fields are default to 'cgrid_velocity'
-        :param chunksize: size of the chunks in dask loading. Default is None (no chunking)
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files,
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable names in the netCDF file(s).
+            Note that the built-in Advection kernels assume that U and V are in m/s
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable.
+            Watch out: NEMO is discretised on a C-grid:
+            U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).
+
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |         V[k,j+1,i+1]        |                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j+1,i]                   |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|U[k,j+1,i+1]                 |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |         V[k,j,i+1]          +                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+
+            To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes,
+            which are located on the corners of the cells.
+            (for indexing details: https://www.nemo-ocean.eu/doc/img360.png )
+            In 3D, the depth is the one corresponding to W nodes
+            The gridindexingtype is set to 'nemo'. See also the Grid indexing documentation on oceanparcels.org
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        tracer_interp_method : str
+            Method for interpolation of tracer fields. It is recommended to use 'cgrid_tracer' (default)
+            Note that in the case of from_nemo() and from_cgrid(), the velocity fields are default to 'cgrid_velocity'
+        chunksize :
+            size of the chunks in dask loading. Default is None (no chunking)
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_c_grid_dataset` constructor.
 
         """
         if 'creation_log' not in kwargs.keys():
@@ -499,8 +569,7 @@ 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.
+        """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
 
@@ -530,60 +599,74 @@ def from_c_grid_dataset(cls, filenames, variables, dimensions, indices=None, mes
                             tracer_interp_method='cgrid_tracer', gridindexingtype='nemo', chunksize=None, **kwargs):
         """Initialises FieldSet object from NetCDF files of Curvilinear NEMO fields.
 
-        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb>`_
+        See `here <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb>`__
         for a more detailed explanation of the different methods that can be used for c-grid datasets.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files,
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable
-               names in the netCDF file(s).
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable.
-               Watch out: NEMO is discretised on a C-grid:
-               U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).
-
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |         V[k,j+1,i+1]        |                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j+1,i]                   |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|U[k,j+1,i+1]                 |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |         V[k,j,i+1]          +                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-
-               To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes,
-               which are located on the corners of the cells.
-               (for indexing details: https://www.nemo-ocean.eu/doc/img360.png )
-               In 3D, the depth is the one corresponding to W nodes.
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param tracer_interp_method: Method for interpolation of tracer fields. It is recommended to use 'cgrid_tracer' (default)
-               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.
-
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files,
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable
+            names in the netCDF file(s).
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable.
+            Watch out: NEMO is discretised on a C-grid:
+            U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).
+
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |         V[k,j+1,i+1]        |                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j+1,i]                   |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|U[k,j+1,i+1]                 |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |         V[k,j,i+1]          +                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+
+            To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes,
+            which are located on the corners of the cells.
+            (for indexing details: https://www.nemo-ocean.eu/doc/img360.png )
+            In 3D, the depth is the one corresponding to W nodes.
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        tracer_interp_method : str
+            Method for interpolation of tracer fields. It is recommended to use 'cgrid_tracer' (default)
+            Note that in the case of from_nemo() and from_cgrid(), the velocity fields are default to 'cgrid_velocity'
+        gridindexingtype : str
+            The type of gridindexing. Set to 'nemo' in FieldSet.from_nemo()
+            See also the Grid indexing documentation on oceanparcels.org (Default value = 'nemo')
+        chunksize :
+            size of the chunks in dask loading. (Default value = None)
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_netcdf` constructor.
         """
         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. "
@@ -614,61 +697,76 @@ def from_pop(cls, filenames, variables, dimensions, indices=None, mesh='spherica
         """Initialises FieldSet object from NetCDF files of POP fields.
             It is assumed that the velocities in the POP fields is in cm/s.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files,
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable names in the netCDF file(s).
-               Note that the built-in Advection kernels assume that U and V are in m/s
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable.
-               Watch out: POP is discretised on a B-grid:
-               U and V velocity nodes are not located as W velocity and T tracer nodes (see http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf ).
-
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j+1,i],V[k,j+1,i]        |                             |U[k,j+1,i+1],V[k,j+1,i+1]    |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j,i],V[k,j,i]            |                             +U[k,j,i+1],V[k,j,i+1]        |
-               +-----------------------------+-----------------------------+-----------------------------+
-
-               In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
-                      T node is at the cell centre and interpolated constant per cell as a C-grid.
-               In 3D: U and V nodes are at the middle of the cell vertical edges,
-                      They are interpolated bilinearly (independently of z) in the cell.
-                      W nodes are at the centre of the horizontal interfaces.
-                      They are interpolated linearly (as a function of z) in the cell.
-                      T node is at the cell centre, and constant per cell.
-                      Note that Parcels assumes that the length of the depth dimension (at the W-points)
-                      is one larger than the size of the velocity and tracer fields in the depth dimension.
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param tracer_interp_method: Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
-               Note that in the case of from_pop() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
-        :param chunksize: size of the chunks in dask loading
-        :param depth_units: The units of the vertical dimension. Default in Parcels is 'm',
-               but many POP outputs are in 'cm'
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files,
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable names in the netCDF file(s).
+            Note that the built-in Advection kernels assume that U and V are in m/s
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable.
+            Watch out: POP is discretised on a B-grid:
+            U and V velocity nodes are not located as W velocity and T tracer nodes (see http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf ).
+
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j+1,i],V[k,j+1,i]        |                             |U[k,j+1,i+1],V[k,j+1,i+1]    |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j,i],V[k,j,i]            |                             +U[k,j,i+1],V[k,j,i+1]        |
+            +-----------------------------+-----------------------------+-----------------------------+
+
+            In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
+            T node is at the cell centre and interpolated constant per cell as a C-grid.
+            In 3D: U and V nodes are at the middle of the cell vertical edges,
+            They are interpolated bilinearly (independently of z) in the cell.
+            W nodes are at the centre of the horizontal interfaces.
+            They are interpolated linearly (as a function of z) in the cell.
+            T node is at the cell centre, and constant per cell.
+            Note that Parcels assumes that the length of the depth dimension (at the W-points)
+            is one larger than the size of the velocity and tracer fields in the depth dimension.
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        tracer_interp_method : str
+            Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
+            Note that in the case of from_pop() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
+        chunksize :
+            size of the chunks in dask loading (Default value = None)
+        depth_units :
+            The units of the vertical dimension. Default in Parcels is 'm',
+            but many POP outputs are in 'cm'
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_b_grid_dataset` constructor.
 
         """
         if 'creation_log' not in kwargs.keys():
@@ -696,59 +794,71 @@ def from_mom5(cls, filenames, variables, dimensions, indices=None, mesh='spheric
                   tracer_interp_method='bgrid_tracer', chunksize=None, **kwargs):
         """Initialises FieldSet object from NetCDF files of MOM5 fields.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files,
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable names in the netCDF file(s).
-               Note that the built-in Advection kernels assume that U and V are in m/s
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable.
-
-               +-------------------------------+-------------------------------+-------------------------------+
-               |U[k,j+1,i],V[k,j+1,i]          |                               |U[k,j+1,i+1],V[k,j+1,i+1]      |
-               +-------------------------------+-------------------------------+-------------------------------+
-               |                               |W[k-1:k+1,j+1,i+1],T[k,j+1,i+1]|                               |
-               +-------------------------------+-------------------------------+-------------------------------+
-               |U[k,j,i],V[k,j,i]              |                               +U[k,j,i+1],V[k,j,i+1]          |
-               +-------------------------------+-------------------------------+-------------------------------+
-
-               In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
-                      T node is at the cell centre and interpolated constant per cell as a C-grid.
-               In 3D: U and V nodes are at the midlle of the cell vertical edges,
-                      They are interpolated bilinearly (independently of z) in the cell.
-                      W nodes are at the centre of the horizontal interfaces, but below the U and V.
-                      They are interpolated linearly (as a function of z) in the cell.
-                      Note that W is normally directed upward in MOM5, but Parcels requires W
-                      in the positive z-direction (downward) so W is multiplied by -1.
-                      T node is at the cell centre, and constant per cell.
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param tracer_interp_method: Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
-               Note that in the case of from_mom5() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
-        :param chunksize: size of the chunks in dask loading
-
-
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files,
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable names in the netCDF file(s).
+            Note that the built-in Advection kernels assume that U and V are in m/s
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable.
+
+            +-------------------------------+-------------------------------+-------------------------------+
+            |U[k,j+1,i],V[k,j+1,i]          |                               |U[k,j+1,i+1],V[k,j+1,i+1]      |
+            +-------------------------------+-------------------------------+-------------------------------+
+            |                               |W[k-1:k+1,j+1,i+1],T[k,j+1,i+1]|                               |
+            +-------------------------------+-------------------------------+-------------------------------+
+            |U[k,j,i],V[k,j,i]              |                               +U[k,j,i+1],V[k,j,i+1]          |
+            +-------------------------------+-------------------------------+-------------------------------+
+
+            In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
+            T node is at the cell centre and interpolated constant per cell as a C-grid.
+            In 3D: U and V nodes are at the midlle of the cell vertical edges,
+            They are interpolated bilinearly (independently of z) in the cell.
+            W nodes are at the centre of the horizontal interfaces, but below the U and V.
+            They are interpolated linearly (as a function of z) in the cell.
+            Note that W is normally directed upward in MOM5, but Parcels requires W
+            in the positive z-direction (downward) so W is multiplied by -1.
+            T node is at the cell centre, and constant per cell.
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        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
+        tracer_interp_method : str
+            Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
+            Note that in the case of from_mom5() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
+        chunksize :
+            size of the chunks in dask loading (Default value = None)
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_b_grid_dataset` constructor.
         """
         if 'creation_log' not in kwargs.keys():
             kwargs['creation_log'] = 'from_mom5'
@@ -765,57 +875,70 @@ def from_b_grid_dataset(cls, filenames, variables, dimensions, indices=None, mes
                             tracer_interp_method='bgrid_tracer', chunksize=None, **kwargs):
         """Initialises FieldSet object from NetCDF files of Bgrid fields.
 
-        :param filenames: Dictionary mapping variables to file(s). The
-               filepath may contain wildcards to indicate multiple files,
-               or be a list of file.
-               filenames can be a list [files], a dictionary {var:[files]},
-               a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
-               or a dictionary of dictionaries {var:{dim:[files]}}
-               time values are in filenames[data]
-        :param variables: Dictionary mapping variables to variable
-               names in the netCDF file(s).
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the netCF file(s).
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable.
-               U and V velocity nodes are not located as W velocity and T tracer nodes (see http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf ).
-
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j+1,i],V[k,j+1,i]        |                             |U[k,j+1,i+1],V[k,j+1,i+1]    |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |                             |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|                             |
-               +-----------------------------+-----------------------------+-----------------------------+
-               |U[k,j,i],V[k,j,i]            |                             +U[k,j,i+1],V[k,j,i+1]        |
-               +-----------------------------+-----------------------------+-----------------------------+
-
-               In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
-                      T node is at the cell centre and interpolated constant per cell as a C-grid.
-               In 3D: U and V nodes are at the midlle of the cell vertical edges,
-                      They are interpolated bilinearly (independently of z) in the cell.
-                      W nodes are at the centre of the horizontal interfaces.
-                      They are interpolated linearly (as a function of z) in the cell.
-                      T node is at the cell centre, and constant per cell.
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param tracer_interp_method: Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
-               Note that in the case of from_pop() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
-        :param chunksize: size of the chunks in dask loading
-
+        Parameters
+        ----------
+        filenames :
+            Dictionary mapping variables to file(s). The
+            filepath may contain wildcards to indicate multiple files,
+            or be a list of file.
+            filenames can be a list [files], a dictionary {var:[files]},
+            a dictionary {dim:[files]} (if lon, lat, depth and/or data not stored in same files as data),
+            or a dictionary of dictionaries {var:{dim:[files]}}
+            time values are in filenames[data]
+        variables : dict
+            Dictionary mapping variables to variable
+            names in the netCDF file(s).
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the netCF file(s).
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable.
+            U and V velocity nodes are not located as W velocity and T tracer nodes (see http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf ).
+
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j+1,i],V[k,j+1,i]        |                             |U[k,j+1,i+1],V[k,j+1,i+1]    |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |                             |W[k:k+2,j+1,i+1],T[k,j+1,i+1]|                             |
+            +-----------------------------+-----------------------------+-----------------------------+
+            |U[k,j,i],V[k,j,i]            |                             +U[k,j,i+1],V[k,j,i+1]        |
+            +-----------------------------+-----------------------------+-----------------------------+
+
+            In 2D: U and V nodes are on the cell vertices and interpolated bilinearly as a A-grid.
+            T node is at the cell centre and interpolated constant per cell as a C-grid.
+            In 3D: U and V nodes are at the midlle of the cell vertical edges,
+            They are interpolated bilinearly (independently of z) in the cell.
+            W nodes are at the centre of the horizontal interfaces.
+            They are interpolated linearly (as a function of z) in the cell.
+            T node is at the cell centre, and constant per cell.
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        tracer_interp_method : str
+            Method for interpolation of tracer fields. It is recommended to use 'bgrid_tracer' (default)
+            Note that in the case of from_pop() and from_bgrid(), the velocity fields are default to 'bgrid_velocity'
+        chunksize :
+            size of the chunks in dask loading (Default value = None)
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_netcdf` constructor.
         """
         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. "
@@ -845,26 +968,41 @@ def from_parcels(cls, basename, uvar='vozocrtx', vvar='vomecrty', indices=None,
                      chunksize=None, **kwargs):
         """Initialises FieldSet data from NetCDF files using the Parcels FieldSet.write() conventions.
 
-        :param basename: Base name of the file(s); may contain
-               wildcards to indicate multiple files.
-        :param indices: Optional dictionary of indices for each dimension
-               to read from file(s), to allow for reading of subset of data.
-               Default is to read the full extent of each dimension.
-               Note that negative indices are not allowed.
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param extra_fields: Extra fields to read beyond U and V
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
-        :param deferred_load: boolean whether to only pre-load data (in deferred mode) or
-               fully load them (default: True). It is advised to deferred load the data, since in
-               that case Parcels deals with a better memory management during particle set execution.
-               deferred_load=False is however sometimes necessary for plotting the fields.
-        :param chunksize: size of the chunks in dask loading
-
+        Parameters
+        ----------
+        basename : str
+            Base name of the file(s); may contain
+            wildcards to indicate multiple files.
+        indices :
+            Optional dictionary of indices for each dimension
+            to read from file(s), to allow for reading of subset of data.
+            Default is to read the full extent of each dimension.
+            Note that negative indices are not allowed.
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        extra_fields :
+            Extra fields to read beyond U and V (Default value = None)
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        deferred_load : bool
+            boolean whether to only pre-load data (in deferred mode) or
+            fully load them (default: True). It is advised to deferred load the data, since in
+            that case Parcels deals with a better memory management during particle set execution.
+            deferred_load=False is however sometimes necessary for plotting the fields.
+        chunksize :
+            size of the chunks in dask loading (Default value = None)
+        uvar :
+             (Default value = 'vozocrtx')
+        vvar :
+             (Default value = 'vomecrty')
+        **kwargs :
+            Keyword arguments passed to the :func:`Fieldset.from_netcdf` constructor.
         """
         if extra_fields is None:
             extra_fields = {}
@@ -889,27 +1027,38 @@ def from_xarray_dataset(cls, ds, variables, dimensions, mesh='spherical', allow_
                             time_periodic=False, **kwargs):
         """Initialises FieldSet data from xarray Datasets.
 
-        :param ds: xarray Dataset.
-               Note that the built-in Advection kernels assume that U and V are in m/s
-        :param variables: Dictionary mapping parcels variable names to data variables in the xarray Dataset.
-        :param dimensions: Dictionary mapping data dimensions (lon,
-               lat, depth, time, data) to dimensions in the xarray Dataset.
-               Note that dimensions can also be a dictionary of dictionaries if
-               dimension names are different for each variable
-               (e.g. dimensions['U'], dimensions['V'], etc).
-        :param fieldtype: Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
-               (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-        :param mesh: String indicating the type of mesh coordinates and
-               units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`_:
-
-               1. spherical (default): Lat and lon in degree, with a
-                  correction for zonal velocity U near the poles.
-               2. flat: No conversion, lat/lon are assumed to be in m.
-        :param allow_time_extrapolation: boolean whether to allow for extrapolation
-               (i.e. beyond the last available time snapshot)
-               Default is False if dimensions includes time, else True
-        :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
+        Parameters
+        ----------
+        ds : xr.Dataset
+            xarray Dataset.
+            Note that the built-in Advection kernels assume that U and V are in m/s
+        variables : dict
+            Dictionary mapping parcels variable names to data variables in the xarray Dataset.
+        dimensions : dict
+            Dictionary mapping data dimensions (lon,
+            lat, depth, time, data) to dimensions in the xarray Dataset.
+            Note that dimensions can also be a dictionary of dictionaries if
+            dimension names are different for each variable
+            (e.g. dimensions['U'], dimensions['V'], etc).
+        fieldtype :
+            Optional dictionary mapping fields to fieldtypes to be used for UnitConverter.
+            (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
+        mesh : str
+            String indicating the type of mesh coordinates and
+            units used during velocity interpolation, see also `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_unitconverters.ipynb>`__:
+
+            1. spherical (default): Lat and lon in degree, with a
+            correction for zonal velocity U near the poles.
+            2. flat: No conversion, lat/lon are assumed to be in m.
+        allow_time_extrapolation : bool
+            boolean whether to allow for extrapolation
+            (i.e. beyond the last available time snapshot)
+            Default is False if dimensions includes time, else True
+        time_periodic : bool, float or datetime.timedelta
+            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
+        **kwargs :
+            Keyword arguments passed to the :func:`Field.from_xarray` constructor.
         """
         fields = {}
         if 'creation_log' not in kwargs.keys():
@@ -951,13 +1100,20 @@ def add_constant(self, name, value):
         stored as 32-bit floats. While constants can be updated during
         execution in SciPy mode, they can not be updated in JIT mode.
 
-        Tutorials using fieldset.add_constant:
-        `Analytical advection <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_analyticaladvection.ipynb>`_
-        `Diffusion <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`_
-        `Periodic boundaries <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb>`_
+        Parameters
+        ----------
+        name : str
+            Name of the constant
+        value :
+            Value of the constant (stored as 32-bit float)
 
-        :param name: Name of the constant
-        :param value: Value of the constant (stored as 32-bit float)
+
+        Examples
+        --------
+        Tutorials using fieldset.add_constant:
+        `Analytical advection <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_analyticaladvection.ipynb>`__
+        `Diffusion <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`__
+        `Periodic boundaries <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb>`__
         """
         setattr(self, name, value)
 
@@ -966,9 +1122,14 @@ def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
         through extending the Field (and lon/lat) by copying a small portion
         of the field on one side of the domain to the other.
 
-        :param zonal: Create a halo in zonal direction (boolean)
-        :param meridional: Create a halo in meridional direction (boolean)
-        :param halosize: size of the halo (in grid points). Default is 5 grid points
+        Parameters
+        ----------
+        zonal : bool
+            Create a halo in zonal direction (Default value = False)
+        meridional : bool
+            Create a halo in meridional direction (Default value = False)
+        halosize : int
+            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)
@@ -979,7 +1140,10 @@ def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
     def write(self, filename):
         """Write FieldSet to NetCDF file using NEMO convention.
 
-        :param filename: Basename of the output fileset.
+        Parameters
+        ----------
+        filename : str
+            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}")
@@ -999,8 +1163,12 @@ def computeTimeChunk(self, time, dt):
         with default option deferred_load. The loaded time steps are at or immediatly before time
         and the two time steps immediately following time if dt is positive (and inversely for negative dt)
 
-        :param time: Time around which the FieldSet chunks are to be loaded. Time is provided as a double, relatively to Fieldset.time_origin
-        :param dt: time step of the integration scheme
+        Parameters
+        ----------
+        time :
+            Time around which the FieldSet chunks are to be loaded. Time is provided as a double, relatively to Fieldset.time_origin
+        dt :
+            time step of the integration scheme
         """
         signdt = np.sign(dt)
         nextTime = np.infty if dt > 0 else -np.infty
diff --git a/parcels/grid.py b/parcels/grid.py
index da848ee8b..42346b2b8 100644
--- a/parcels/grid.py
+++ b/parcels/grid.py
@@ -96,7 +96,6 @@ def child_ctypes_struct(self):
         """Returns a ctypes struct object containing all relevant
         pointers and sizes for this grid.
         """
-
         class CStructuredGrid(Structure):
             # z4d is only to have same cstruct as RectilinearSGrid
             _fields_ = [('xdim', c_int), ('ydim', c_int), ('zdim', c_int),
@@ -266,9 +265,14 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
         """Add a 'halo' to the Grid, through extending the Grid (and lon/lat)
         similarly to the halo created for the Fields
 
-        :param zonal: Create a halo in zonal direction (boolean)
-        :param meridional: Create a halo in meridional direction (boolean)
-        :param halosize: size of the halo (in grid points). Default is 5 grid points
+        Parameters
+        ----------
+        zonal : bool
+            Create a halo in zonal direction
+        meridional : bool
+            Create a halo in meridional direction
+        halosize : int
+            size of the halo (in grid points). Default is 5 grid points
         """
         if zonal:
             lonshift = (self.lon[-1] - 2 * self.lon[0] + self.lon[1])
@@ -295,18 +299,26 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
 class RectilinearZGrid(RectilinearGrid):
     """Rectilinear Z Grid.
 
-    :param lon: Vector containing the longitude coordinates of the grid
-    :param lat: Vector containing the latitude coordinates of the grid
-    :param depth: Vector containing the vertical coordinates of the grid, which are z-coordinates.
-           The depth of the different layers is thus constant.
-    :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (TimeConverter object) of the time axis
-    :param mesh: String indicating the type of mesh coordinates and
-           units used during velocity interpolation:
-
-           1. spherical (default): Lat and lon in degree, with a
-              correction for zonal velocity U near the poles.
-           2. flat: No conversion, lat/lon are assumed to be in m.
+    Parameters
+    ----------
+    lon :
+        Vector containing the longitude coordinates of the grid
+    lat :
+        Vector containing the latitude coordinates of the grid
+    depth :
+        Vector containing the vertical coordinates of the grid, which are z-coordinates.
+        The depth of the different layers is thus constant.
+    time :
+        Vector containing the time coordinates of the grid
+    time_origin : parcels.tools.converters.TimeConverter
+        Time origin of the time axis
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical (default): Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat: No conversion, lat/lon are assumed to be in m.
     """
 
     def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat'):
@@ -328,23 +340,31 @@ class RectilinearSGrid(RectilinearGrid):
     """Rectilinear S Grid. Same horizontal discretisation as a rectilinear z grid,
        but with s vertical coordinates
 
-    :param lon: Vector containing the longitude coordinates of the grid
-    :param lat: Vector containing the latitude coordinates of the grid
-    :param depth: 4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid,
-           which are s-coordinates.
-           s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
-           or any generalised vertical discretisation.
-           The depth of each node depends then on the horizontal position (lon, lat),
-           the number of the layer and the time is depth is a 4D array.
-           depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
-    :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (TimeConverter object) of the time axis
-    :param mesh: String indicating the type of mesh coordinates and
-           units used during velocity interpolation:
-
-           1. spherical (default): Lat and lon in degree, with a
-              correction for zonal velocity U near the poles.
-           2. flat: No conversion, lat/lon are assumed to be in m.
+    Parameters
+    ----------
+    lon :
+        Vector containing the longitude coordinates of the grid
+    lat :
+        Vector containing the latitude coordinates of the grid
+    depth :
+        4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid,
+        which are s-coordinates.
+        s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
+        or any generalised vertical discretisation.
+        The depth of each node depends then on the horizontal position (lon, lat),
+        the number of the layer and the time is depth is a 4D array.
+        depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
+    time :
+        Vector containing the time coordinates of the grid
+    time_origin : parcels.tools.converters.TimeConverter
+        Time origin of the time axis
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical (default): Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat: No conversion, lat/lon are assumed to be in m.
     """
 
     def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'):
@@ -391,9 +411,14 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
         """Add a 'halo' to the Grid, through extending the Grid (and lon/lat)
         similarly to the halo created for the Fields
 
-        :param zonal: Create a halo in zonal direction (boolean)
-        :param meridional: Create a halo in meridional direction (boolean)
-        :param halosize: size of the halo (in grid points). Default is 5 grid points
+        Parameters
+        ----------
+        zonal : bool
+            Create a halo in zonal direction
+        meridional : bool
+            Create a halo in meridional direction
+        halosize : int
+            size of the halo (in grid points). Default is 5 grid points
         """
         if zonal:
             lonshift = self.lon[:, -1] - 2 * self.lon[:, 0] + self.lon[:, 1]
@@ -429,18 +454,26 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
 class CurvilinearZGrid(CurvilinearGrid):
     """Curvilinear Z Grid.
 
-    :param lon: 2D array containing the longitude coordinates of the grid
-    :param lat: 2D array containing the latitude coordinates of the grid
-    :param depth: Vector containing the vertical coordinates of the grid, which are z-coordinates.
-           The depth of the different layers is thus constant.
-    :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (TimeConverter object) of the time axis
-    :param mesh: String indicating the type of mesh coordinates and
-           units used during velocity interpolation:
-
-           1. spherical (default): Lat and lon in degree, with a
-              correction for zonal velocity U near the poles.
-           2. flat: No conversion, lat/lon are assumed to be in m.
+    Parameters
+    ----------
+    lon :
+        2D array containing the longitude coordinates of the grid
+    lat :
+        2D array containing the latitude coordinates of the grid
+    depth :
+        Vector containing the vertical coordinates of the grid, which are z-coordinates.
+        The depth of the different layers is thus constant.
+    time :
+        Vector containing the time coordinates of the grid
+    time_origin : parcels.tools.converters.TimeConverter
+        Time origin of the time axis
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical (default): Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat: No conversion, lat/lon are assumed to be in m.
     """
 
     def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat'):
@@ -461,23 +494,31 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat
 class CurvilinearSGrid(CurvilinearGrid):
     """Curvilinear S Grid.
 
-    :param lon: 2D array containing the longitude coordinates of the grid
-    :param lat: 2D array containing the latitude coordinates of the grid
-    :param depth: 4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid,
-           which are s-coordinates.
-           s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
-           or any generalised vertical discretisation.
-           The depth of each node depends then on the horizontal position (lon, lat),
-           the number of the layer and the time is depth is a 4D array.
-           depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
-    :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (TimeConverter object) of the time axis
-    :param mesh: String indicating the type of mesh coordinates and
-           units used during velocity interpolation:
-
-           1. spherical (default): Lat and lon in degree, with a
-              correction for zonal velocity U near the poles.
-           2. flat: No conversion, lat/lon are assumed to be in m.
+    Parameters
+    ----------
+    lon :
+        2D array containing the longitude coordinates of the grid
+    lat :
+        2D array containing the latitude coordinates of the grid
+    depth :
+        4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid,
+        which are s-coordinates.
+        s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
+        or any generalised vertical discretisation.
+        The depth of each node depends then on the horizontal position (lon, lat),
+        the number of the layer and the time is depth is a 4D array.
+        depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
+    time :
+        Vector containing the time coordinates of the grid
+    time_origin : parcels.tools.converters.TimeConverter
+        Time origin of the time axis
+    mesh : str
+        String indicating the type of mesh coordinates and
+        units used during velocity interpolation:
+
+        1. spherical (default): Lat and lon in degree, with a
+        correction for zonal velocity U near the poles.
+        2. flat: No conversion, lat/lon are assumed to be in m.
     """
 
     def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'):
diff --git a/parcels/interaction/neighborsearch/base.py b/parcels/interaction/neighborsearch/base.py
index 4c59cda80..a5d44fc1a 100644
--- a/parcels/interaction/neighborsearch/base.py
+++ b/parcels/interaction/neighborsearch/base.py
@@ -18,10 +18,15 @@ def __init__(self, inter_dist_vert, inter_dist_horiz,
                  max_depth=100000, periodic_domain_zonal=None):
         """Initialize neighbor search
 
-        :param inter_dist_vert: interaction distance (vertical) in m
-        :param inter_dist_horiz: interaction distance (horizontal in m
-        :param max_depth: maximum depth of the particles (i.e. 100km)
-        :param zperiodic_bc_domain: zonal domain if zonal periodic boundary
+
+        Parameters
+        ----------
+        inter_dist_vert : float
+            Interaction distance (vertical) in m.
+        inter_dist_horiz : float
+            interaction distance (horizontal) in m
+        max_depth : float, optional
+            Maximum depth of the particles (default is 100000m).
         """
         self.inter_dist_vert = inter_dist_vert
         self.inter_dist_horiz = inter_dist_horiz
@@ -43,8 +48,16 @@ def __init__(self, inter_dist_vert, inter_dist_horiz,
     def find_neighbors_by_coor(self, coor):
         """Get the neighbors around a certain location.
 
-        :param coor: Numpy array with [depth, lat, lon].
-        :returns List of particle indices.
+        Parameters
+        ----------
+        coor :
+            Numpy array with [depth, lat, lon].
+
+        Returns
+        -------
+        type
+            List of particle indices.
+
         """
         raise NotImplementedError
 
@@ -53,8 +66,16 @@ def find_neighbors_by_idx(self, particle_idx):
 
         Mainly useful for Structure of Array (SoA) datastructure
 
-        :param particle_idx: index of the particle (SoA).
-        :returns List of particle indices
+        Parameters
+        ----------
+        particle_idx :
+            index of the particle (SoA).
+
+        Returns
+        -------
+        type
+            List of particle indices
+
         """
         coor = self._values[:, particle_idx].reshape(3, 1)
         return self.find_neighbors_by_coor(coor)
@@ -65,18 +86,26 @@ def update_values(self, new_values, new_active_mask=None):
         This is a default implementation simply rebuilds the structure.
         If the rebuilding is slow, a faster implementation can be provided.
 
-        :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.
+        Parameters
+        ----------
+        new_values :
+            numpy array ([depth, lat, lon], n_particles) with
+            new coordinates of the particles.
+        new_active_mask :
+            boolean array indicating active particles. (Default value = None)
         """
         self.rebuild(new_values, new_active_mask)
 
     def rebuild(self, values, active_mask=-1):
         """Rebuild the neighbor structure from scratch.
 
-        :param values: numpy array with coordinates of particles
-                       (same as update).
-        :param active_mask: boolean array indicating active particles.
+        Parameters
+        ----------
+        values :
+            numpy array with coordinates of particles
+            (same as update).
+        active_mask :
+            boolean array indicating active particles. (Default value = -1)
         """
         if values is not None:
             self._values = values
@@ -102,20 +131,36 @@ def _distance(self, coor, subset_idx):
 
         Distance depends on the mesh (spherical/flat).
 
-        :param coor: Numpy array with 3D coordinates ([depth, lat, lon]).
-        :param subset_idx: Indices of the particles to compute the distance to.
-        :returns horiz_dist: distance in the horizontal direction
-        :returns vert_dist: distance in the vertical direction.
+        Parameters
+        ----------
+        coor :
+            Numpy array with 3D coordinates ([depth, lat, lon]).
+        subset_idx :
+            Indices of the particles to compute the distance to.
+
+        Returns
+        -------
+        type
+            horiz_dist: distance in the horizontal direction
+
         """
         raise NotImplementedError
 
     def _get_close_neighbor_dist(self, coor, subset_idx):
         """Compute distances and remove non-neighbors.
 
-        :param coor: Numpy array with 3D coordinates ([depth, lat, lon]).
-        :param subset_idx: Indices of the particles to compute the distance to.
-        :returns neighbor_idx: Indices within the interaction distance.
-        :returns distances: Distance between coor and the neighbor particles.
+        Parameters
+        ----------
+        coor :
+            Numpy array with 3D coordinates ([depth, lat, lon]).
+        subset_idx :
+            Indices of the particles to compute the distance to.
+
+        Returns
+        -------
+        type
+            neighbor_idx: Indices within the interaction distance.
+
         """
         vert_distance, horiz_distance = self._distance(coor, subset_idx)
         rel_distances = np.sqrt((horiz_distance/self.inter_dist_horiz)**2
diff --git a/parcels/interaction/neighborsearch/basehash.py b/parcels/interaction/neighborsearch/basehash.py
index 550f94f85..9150241a3 100644
--- a/parcels/interaction/neighborsearch/basehash.py
+++ b/parcels/interaction/neighborsearch/basehash.py
@@ -7,8 +7,16 @@ class BaseHashNeighborSearch(ABC):
     def find_neighbors_by_coor(self, coor):
         """Get the neighbors around a certain location.
 
-        :param coor: Numpy array with [depth, lat, lon].
-        :returns List of particle indices.
+        Parameters
+        ----------
+        coor :
+            Numpy array with [depth, lat, lon].
+
+        Returns
+        -------
+        type
+            List of particle indices.
+
         """
         coor = coor.reshape(3, 1)
         hash_id = self._values_to_hashes(coor)[0]
@@ -19,8 +27,16 @@ def find_neighbors_by_idx(self, particle_idx):
 
         Mainly useful for Structure of Array (SoA) datastructure
 
-        :param particle_idx: index of the particle (SoA).
-        :returns List of particle indices
+        Parameters
+        ----------
+        particle_idx :
+            index of the particle (SoA).
+
+        Returns
+        -------
+        type
+            List of particle indices
+
         """
         hash_id = self._particle_hashes[particle_idx]
         coor = self._values[:, particle_idx].reshape(3, 1)
@@ -59,7 +75,12 @@ def update_values(self, new_values, new_active_mask=None):
         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.
+        Parameters
+        ----------
+        new_values :
+            new (depth, lat, lon) values for particles.
+        new_active_mask :
+             (Default value = None)
         """
         if self._values is None:
             self.rebuild(new_values, new_active_mask)
@@ -103,9 +124,18 @@ def _values_to_hashes(self, values, active_idx=None):
 
         The hashes correspond to the cells that particles reside in.
 
-        :param values: 3D coordinates to be hashed.
-        :param active_idx: Active particle indices (relative to values).
-        :returns all_hashes: An array of length len(values) with hashes.
+        Parameters
+        ----------
+        values :
+            3D coordinates to be hashed.
+        active_idx :
+            Active particle indices (relative to values). (Default value = None)
+
+        Returns
+        -------
+        type
+            all_hashes: An array of length len(values) with hashes.
+
         """
         raise NotImplementedError
 
@@ -148,9 +178,18 @@ def hash_split(hash_ids, active_idx=None):
     Multiple particles that are found in the same cell are put in a list
     with that particular hash.
 
-    :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, ..], ..}
+    Parameters
+    ----------
+    hash_ids :
+        Hash values for the particles.
+    active_idx :
+        Subset on which to compute the hash split. (Default value = None)
+
+    Returns
+    -------
+    type
+        hash_split: Dictionary with {hash: [idx_1, idx_2, ..], ..}
+
     """
     if len(hash_ids) == 0:
         return {}
diff --git a/parcels/interaction/neighborsearch/hashflat.py b/parcels/interaction/neighborsearch/hashflat.py
index cabecd507..66c8a7ffc 100644
--- a/parcels/interaction/neighborsearch/hashflat.py
+++ b/parcels/interaction/neighborsearch/hashflat.py
@@ -33,9 +33,18 @@ def update_values(self, new_values, new_active_mask=None):
     def _check_box(self, new_values, new_active_mask):
         """Check whether particles have moved out of the overall box.
 
-        :param new_values: New particle coordinates (depth, lat, lon) to be checked.
-        :param new_active_mask: New active mask for the particles.
-        :returns True if box is still big enough, False if not.
+        Parameters
+        ----------
+        new_values :
+            New particle coordinates (depth, lat, lon) to be checked.
+        new_active_mask :
+            New active mask for the particles.
+
+        Returns
+        -------
+        type
+            True if box is still big enough, False if not.
+
         """
         if self._box is None:
             return False
@@ -111,9 +120,18 @@ def _values_to_hashes(self, values, active_idx=None):
 def hash_to_neighbors(hash_id, bits):
     """Compute neighboring cells from a hash.
 
-    :param hash_id: hash value of the current cell.
-    :param bits: key to compute the hashesh.
-    :returns neighbors: List of cells neighboring hash_id.
+    Parameters
+    ----------
+    hash_id :
+        hash value of the current cell.
+    bits :
+        key to compute the hashesh.
+
+    Returns
+    -------
+    type
+        neighbors: List of cells neighboring hash_id.
+
     """
     coor = np.zeros((len(bits),), dtype=np.int32)
     new_coor = np.zeros((len(bits),), dtype=np.int32)
diff --git a/parcels/interaction/neighborsearch/hashspherical.py b/parcels/interaction/neighborsearch/hashspherical.py
index df770b246..67d293689 100644
--- a/parcels/interaction/neighborsearch/hashspherical.py
+++ b/parcels/interaction/neighborsearch/hashspherical.py
@@ -11,17 +11,21 @@
 
 class HashSphericalNeighborSearch(BaseHashNeighborSearch,
                                   BaseSphericalNeighborSearch):
-    """Neighbor search using a hashtable (similar to octtrees)."""
+    """Neighbor search using a hashtable (similar to octtrees).
+
+
+    Parameters
+    ----------
+    inter_dist_vert : float
+        Interaction distance (vertical) in m.
+    inter_dist_horiz : float
+        interaction distance (horizontal) in m
+    max_depth : float, optional
+        Maximum depth of the particles (default is 100000m).
+    """
 
     def __init__(self, inter_dist_vert, inter_dist_horiz,
                  max_depth=100000):
-        """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()
@@ -46,9 +50,19 @@ def _find_neighbors(self, hash_id, coor):
     def _values_to_hashes(self, values, active_idx=None):
         """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.
+        Parameters
+        ----------
+        values :
+            array of positions of particles to convert
+            ([depth, lat, lon], # of particles to convert).
+        active_idx :
+             (Default value = None)
+
+        Returns
+        -------
+        type
+            array of cell ids.
+
         """
         if active_idx is None:
             active_idx = np.arange(values.shape[1], dtype=int)
@@ -81,7 +95,12 @@ def _values_to_hashes(self, values, active_idx=None):
     def rebuild(self, values, active_mask=-1):
         """Recreate the tree with new values.
 
-        :param values: positions of the particles.
+        Parameters
+        ----------
+        values :
+            positions of the particles.
+        active_mask :
+             (Default value = -1)
         """
         super().rebuild(values, active_mask)
         active_idx = self.active_idx
diff --git a/parcels/kernel/basekernel.py b/parcels/kernel/basekernel.py
index 3a8eeaebb..b5240e14e 100644
--- a/parcels/kernel/basekernel.py
+++ b/parcels/kernel/basekernel.py
@@ -43,13 +43,22 @@
 class BaseKernel:
     """Base super class for base Kernel objects that encapsulates auto-generated code.
 
-    :arg fieldset: FieldSet object providing the field information (possibly None)
-    :arg ptype: PType object for the kernel particle
-    :arg pyfunc: (aggregated) Kernel function
-    :arg funcname: function name
-    :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
-
-    Note: A Kernel is either created from a compiled <function ...> object
+    Parameters
+    ----------
+    fieldset : parcels.Fieldset
+        FieldSet object providing the field information (possibly None)
+    ptype :
+        PType object for the kernel particle
+    pyfunc :
+        (aggregated) Kernel function
+    funcname : str
+        function name
+    delete_cfiles : bool
+        Whether to delete the C-files after compilation in JIT mode (default is True)
+
+    Notes
+    -----
+    A Kernel is either created from a compiled <function ...> object
     or the necessary information (funcname, funccode, funcvars) is provided.
     The py_ast argument may be derived from the code string, but for
     concatenation, the merged AST plus the new header definition is required.
@@ -344,13 +353,22 @@ def load_fieldset_jit(self, pset):
                     g.lat = np.array(g.lat, order='C')
 
     def evaluate_particle(self, p, endtime, sign_dt, dt, analytical=False):
-        """
-        Execute the kernel evaluation of for an individual particle.
-        :arg p: object of (sub-)type (ScipyParticle, JITParticle) or (sub-)type of BaseParticleAccessor
-        :arg fieldset: fieldset of the containing ParticleSet (e.g. pset.fieldset)
-        :arg analytical: flag indicating the analytical advector or an iterative advection
-        :arg endtime: endtime of this overall kernel evaluation step
-        :arg dt: computational integration timestep
+        """Execute the kernel evaluation of for an individual particle.
+
+        Parameters
+        ----------
+        p :
+            object of (sub-)type (ScipyParticle, JITParticle) or (sub-)type of BaseParticleAccessor
+        fieldset :
+            fieldset of the containing ParticleSet (e.g. pset.fieldset)
+        analytical :
+            flag indicating the analytical advector or an iterative advection (Default value = False)
+        endtime :
+            endtime of this overall kernel evaluation step
+        dt :
+            computational integration timestep
+        sign_dt :
+
         """
         variables = self._ptype.variables
         # back up variables in case of OperationCode.Repeat
diff --git a/parcels/kernel/kernelaos.py b/parcels/kernel/kernelaos.py
index 8fa7dd892..6e8db57f4 100644
--- a/parcels/kernel/kernelaos.py
+++ b/parcels/kernel/kernelaos.py
@@ -28,11 +28,18 @@
 class KernelAOS(BaseKernel):
     """Kernel object that encapsulates auto-generated code.
 
-    :arg fieldset: FieldSet object providing the field information
-    :arg ptype: PType object for the kernel particle
-    :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
-
-    Note: A Kernel is either created from a compiled <function ...> object
+    Parameters
+    ----------
+    fieldset :
+        FieldSet object providing the field information
+    ptype :
+        PType object for the kernel particle
+    delete_cfiles : bool
+        Whether to delete the C-files after compilation in JIT mode (default is True)
+
+    Notes
+    -----
+    A Kernel is either created from a compiled <function ...> object
     or the necessary information (funcname, funccode, funcvars) is provided.
     The py_ast argument may be derived from the code string, but for
     concatenation, the merged AST plus the new header definition is required.
diff --git a/parcels/kernel/kernelsoa.py b/parcels/kernel/kernelsoa.py
index 29d2a5633..bfcddf096 100644
--- a/parcels/kernel/kernelsoa.py
+++ b/parcels/kernel/kernelsoa.py
@@ -28,11 +28,18 @@
 class KernelSOA(BaseKernel):
     """Kernel object that encapsulates auto-generated code.
 
-    :arg fieldset: FieldSet object providing the field information
-    :arg ptype: PType object for the kernel particle
-    :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
-
-    Note: A Kernel is either created from a compiled <function ...> object
+    Parameters
+    ----------
+    fieldset :
+        FieldSet object providing the field information
+    ptype :
+        PType object for the kernel particle
+    delete_cfiles : bool
+        Whether to delete the C-files after compilation in JIT mode (default is True)
+
+    Notes
+    -----
+    A Kernel is either created from a compiled <function ...> object
     or the necessary information (funcname, funccode, funcvars) is provided.
     The py_ast argument may be derived from the code string, but for
     concatenation, the merged AST plus the new header definition is required.
diff --git a/parcels/particle.py b/parcels/particle.py
index 9bf2a3798..cc4af961e 100644
--- a/parcels/particle.py
+++ b/parcels/particle.py
@@ -15,13 +15,18 @@
 class Variable:
     """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
-    :param initial: Initial value of the variable. Note that this can also be a Field object,
-             which will then be sampled at the location of the particle
-    :param to_write: Boolean or 'once'. Controls whether Variable is written to NetCDF file.
-             If to_write = 'once', the variable will be written as a time-independent 1D array
-    :type to_write: (bool, 'once', optional)
+    Parameters
+    ----------
+    name : str
+        Variable name as used within kernels
+    dtype :
+        Data type (numpy.dtype) of the variable
+    initial :
+        Initial value of the variable. Note that this can also be a Field object,
+        which will then be sampled at the location of the particle
+    to_write : bool, 'once', optional
+        Boolean or 'once'. Controls whether Variable is written to NetCDF file.
+        If to_write = 'once', the variable will be written as a time-independent 1D array
     """
 
     def __init__(self, name, dtype=np.float32, initial=0, to_write=True):
@@ -55,7 +60,10 @@ def is64bit(self):
 class ParticleType:
     """Class encapsulating the type information for custom particles.
 
-    :param user_vars: Optional list of (name, dtype) tuples for custom variables
+    Parameters
+    ----------
+    user_vars :
+        Optional list of (name, dtype) tuples for custom variables
     """
 
     def __init__(self, pclass):
@@ -172,12 +180,22 @@ def setLastID(cls, offset):
 class ScipyParticle(_Particle):
     """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
-    :param depth: Initial depth of particle
-    :param fieldset: :mod:`parcels.fieldset.FieldSet` object to track this particle on
-    :param time: Current time of the particle
-
+    Parameters
+    ----------
+    lon : float
+        Initial longitude of particle
+    lat : float
+        Initial latitude of particle
+    depth : float
+        Initial depth of particle
+    fieldset : parcels.fieldset.FieldSet
+        mod:`parcels.fieldset.FieldSet` object to track this particle on
+    time : float
+        Current time of the particle
+
+
+    Notes
+    -----
     Additional Variables can be added via the :Class Variable: objects
     """
 
@@ -232,16 +250,25 @@ class ScipyInteractionParticle(ScipyParticle):
 class JITParticle(ScipyParticle):
     """Particle class for JIT-based (Just-In-Time) Particle objects.
 
-    :param lon: Initial longitude of particle
-    :param lat: Initial latitude of particle
-    :param fieldset: :mod:`parcels.fieldset.FieldSet` object to track this particle on
-    :param dt: Execution timestep for this particle
-    :param time: Current time of the particle
-
+    Parameters
+    ----------
+    lon : float
+        Initial longitude of particle
+    lat : float
+        Initial latitude of particle
+    fieldset : parcels.fieldset.FieldSet
+        mod:`parcels.fieldset.FieldSet` object to track this particle on
+    dt :
+        Execution timestep for this particle
+    time :
+        Current time of the particle
+
+
+    Notes
+    -----
     Additional Variables can be added via the :Class Variable: objects
 
     Users should use JITParticles for faster advection computation.
-
     """
 
     def __init__(self, *args, **kwargs):
diff --git a/parcels/particlefile/baseparticlefile.py b/parcels/particlefile/baseparticlefile.py
index e7b676bf3..755d2b1ce 100644
--- a/parcels/particlefile/baseparticlefile.py
+++ b/parcels/particlefile/baseparticlefile.py
@@ -33,14 +33,27 @@ def _set_calendar(origin_calendar):
 class BaseParticleFile(ABC):
     """Initialise trajectory output.
 
-    :param name: Basename of the output file. This can also be a Zarr store object.
-    :param particleset: ParticleSet to output
-    :param outputdt: Interval which dictates the update frequency of file output
-                     while ParticleFile is given as an argument of ParticleSet.execute()
-                     It is either a timedelta object or a positive double.
-    :param chunks: Tuple (trajs, obs) to control the size of chunks in the zarr output.
-    :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
+    Parameters
+    ----------
+    name : str
+        Basename of the output file. This can also be a Zarr store object.
+    particleset :
+        ParticleSet to output
+    outputdt :
+        Interval which dictates the update frequency of file output
+        while ParticleFile is given as an argument of ParticleSet.execute()
+        It is either a timedelta object or a positive double.
+    chunks :
+        Tuple (trajs, obs) to control the size of chunks in the zarr output.
+    write_ondelete : bool
+        Whether to write particle data only when they are deleted. Default is False
+    create_new_zarrfile : bool
+        Whether to create a new file. Default is True
+
+    Returns
+    -------
+    BaseParticleFile
+        ParticleFile object that can be used to write particle data to file
     """
 
     write_ondelete = None
@@ -117,10 +130,10 @@ def _reserved_var_names(self):
         pass
 
     def _create_variables_attribute_dict(self):
-        """
-        Creates the dictionary with variable attributes.
+        """Creates the dictionary with variable attributes.
 
-        Attention:
+        Notes
+        -----
         For ParticleSet structures other than SoA, and structures where ID != index, this has to be overridden.
         """
         attrs = {'z': {"long_name": "",
@@ -166,8 +179,12 @@ def close(self, delete_tempfiles=True):
     def add_metadata(self, name, message):
         """Add metadata to :class:`parcels.particleset.ParticleSet`.
 
-        :param name: Name of the metadata variabale
-        :param message: message to be written
+        Parameters
+        ----------
+        name : str
+            Name of the metadata variabale
+        message : str
+            message to be written
         """
         self.metadata[name] = message
 
@@ -200,9 +217,14 @@ def _extend_zarr_dims(self, Z, store, dtype, axis):
     def write(self, pset, time, deleted_only=False):
         """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
+        Parameters
+        ----------
+        pset :
+            ParticleSet object to write
+        time :
+            Time at which to write ParticleSet
+        deleted_only :
+            Flag to write only the deleted Particles (Default value = False)
         """
         time = time.total_seconds() if isinstance(time, delta) else time
 
diff --git a/parcels/particlefile/particlefileaos.py b/parcels/particlefile/particlefileaos.py
index 1d747e42a..14d911020 100644
--- a/parcels/particlefile/particlefileaos.py
+++ b/parcels/particlefile/particlefileaos.py
@@ -9,13 +9,25 @@
 class ParticleFileAOS(BaseParticleFile):
     """Initialise trajectory output.
 
-    :param name: Basename of the output file. This can also be a Zarr store.
-    :param particleset: ParticleSet to output
-    :param outputdt: Interval which dictates the update frequency of file output
-                     while ParticleFile is given as an argument of ParticleSet.execute()
-                     It is either a timedelta object or a positive double.
-    :param chunks: Tuple (trajs, obs) to control the size of chunks in the zarr output.
-    :param write_ondelete: Boolean to write particle data only when they are deleted. Default is False
+    Parameters
+    ----------
+    name : str
+        Basename of the output file. This can also be a Zarr store.
+    particleset :
+        ParticleSet to output
+    outputdt :
+        Interval which dictates the update frequency of file output
+        while ParticleFile is given as an argument of ParticleSet.execute()
+        It is either a timedelta object or a positive double.
+    chunks :
+        Tuple (trajs, obs) to control the size of chunks in the zarr output.
+    write_ondelete :
+        Boolean to write particle data only when they are deleted. Default is False
+
+    Returns
+    -------
+    ParticleFileAOS
+        ParticleFile object that can be used to write particle data to file
     """
 
     def __init__(self, name, particleset, outputdt=np.infty, chunks=None, write_ondelete=False):
diff --git a/parcels/particlefile/particlefilesoa.py b/parcels/particlefile/particlefilesoa.py
index c5031e3f4..aeb9ffc48 100644
--- a/parcels/particlefile/particlefilesoa.py
+++ b/parcels/particlefile/particlefilesoa.py
@@ -9,13 +9,26 @@
 class ParticleFileSOA(BaseParticleFile):
     """Initialise trajectory output.
 
-    :param name: Basename of the output file.  This can also be a Zarr store.
-    :param particleset: ParticleSet to output
-    :param outputdt: Interval which dictates the update frequency of file output
-                     while ParticleFile is given as an argument of ParticleSet.execute()
-                     It is either a timedelta object or a positive double.
-    :param chunks: Tuple (trajs, obs) to control the size of chunks in the zarr output.
-    :param write_ondelete: Boolean to write particle data only when they are deleted. Default is False
+    Parameters
+    ----------
+    name : str
+        Basename of the output file.  This can also be a Zarr store.
+    particleset :
+        ParticleSet to output
+    outputdt :
+        Interval which dictates the update frequency of file output
+        while ParticleFile is given as an argument of ParticleSet.execute()
+        It is either a timedelta object or a positive double.
+    chunks :
+        Tuple (trajs, obs) to control the size of chunks in the zarr output.
+    write_ondelete :
+        Whether to write particle data only when they are deleted. Default is False
+
+    Returns
+    -------
+    ParticleFileSOA
+        ParticleFile object that can be used to write particle data to file
+
     """
 
     def __init__(self, name, particleset, outputdt=np.infty, chunks=None, write_ondelete=False):
diff --git a/parcels/particleset/baseparticleset.py b/parcels/particleset/baseparticleset.py
index 8dcca9fb9..e40378608 100644
--- a/parcels/particleset/baseparticleset.py
+++ b/parcels/particleset/baseparticleset.py
@@ -68,7 +68,10 @@ def __getattr__(self, name):
         """
         Access a single property of all particles.
 
-        :param name: name of the property
+        Parameters
+        ----------
+        name : str
+            Name of the property
         """
         for v in self._collection.ptype.variables:
             if v.name == name:
@@ -110,18 +113,29 @@ def __create_progressbar(self, starttime, endtime):
     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.
 
-        :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
-        :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
-        :param time: Optional list of start time values for particles. Default is fieldset.U.time[0]
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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')
+        Parameters
+        ----------
+        fieldset :
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        lon :
+            List of initial longitude values for particles
+        lat :
+            List of initial latitude values for particles
+        depth :
+            Optional list of initial depth values for particles. Default is 0m
+        time :
+            Optional list of start time values for particles. Default is fieldset.U.time[0]
+        repeatdt :
+            Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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')
+        **kwargs :
+            Keyword arguments passed to the particleset constructor.
         """
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt, lonlatdepth_dtype=lonlatdepth_dtype, **kwargs)
 
@@ -133,18 +147,28 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None,
         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
 
-        :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
-        :param start: Starting point for initialisation of particles on a straight line.
-        :param finish: End point for initialisation of particles on a straight line.
-        :param size: Initial size of particle set
-        :param depth: Optional list of initial depth values for particles. Default is 0m
-        :param time: Optional start time value for particles. Default is fieldset.U.time[0]
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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'
+        Parameters
+        ----------
+        fieldset :
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        start :
+            Starting point for initialisation of particles on a straight line.
+        finish :
+            End point for initialisation of particles on a straight line.
+        size :
+            Initial size of particle set
+        depth :
+            Optional list of initial depth values for particles. Default is 0m
+        time :
+            Optional start time value for particles. Default is fieldset.U.time[0]
+        repeatdt :
+            Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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)
@@ -155,12 +179,23 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None,
     @classmethod
     @abstractmethod
     def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'):
-        """
-        Converts a starting field into a monte-carlo sample of lons and lats.
-
-        :param start_field: :mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
-
-        returns list(lon), list(lat)
+        """Converts a starting field into a monte-carlo sample of lons and lats.
+
+        Parameters
+        ----------
+        start_field : parcels.field.Field
+            mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
+        size :
+
+        mode :
+             (Default value = 'monte_carlo')
+
+        Returns
+        -------
+        list of float
+            A list of longitude values.
+        list of float
+            A list of latitude values.
         """
         pass
 
@@ -168,18 +203,28 @@ def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'):
     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.
 
-        :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
-        :param start_field: Field for initialising particles stochastically (horizontally)  according to the presented density field.
-        :param size: Initial size of particle set
-        :param mode: Type of random sampling. Currently only 'monte_carlo' is implemented
-        :param depth: Optional list of initial depth values for particles. Default is 0m
-        :param time: Optional start time value for particles. Default is fieldset.U.time[0]
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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'
+        Parameters
+        ----------
+        fieldset : parcels.fieldset.FieldSet
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        start_field : parcels.field.Field
+            Field for initialising particles stochastically (horizontally)  according to the presented density field.
+        size :
+            Initial size of particle set
+        mode :
+            Type of random sampling. Currently only 'monte_carlo' is implemented (Default value = 'monte_carlo')
+        depth :
+            Optional list of initial depth values for particles. Default is 0m
+        time :
+            Optional start time value for particles. Default is fieldset.U.time[0]
+        repeatdt :
+            Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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)
 
@@ -192,34 +237,50 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
         This creates a new ParticleSet based on locations of all particles written
         in a netcdf ParticleFile at a certain time. Particle IDs are preserved if restart=True
 
-        :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
-        :param filename: Name of the particlefile from which to read initial conditions
-        :param restart: Boolean to signal if pset is used for a restart (default is True).
-               In that case, Particle IDs are preserved.
-        :param restarttime: time at which the Particles will be restarted. Default is the last time written.
-               Alternatively, restarttime could be a time value (including np.datetime64) or
-               a callable function such as np.nanmin. The last is useful when running with dt < 0.
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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'
+        Parameters
+        ----------
+        fieldset :
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        filename : str
+            Name of the particlefile from which to read initial conditions
+        restart : bool
+            Boolean to signal if pset is used for a restart (default is True).
+            In that case, Particle IDs are preserved.
+        restarttime :
+            time at which the Particles will be restarted. Default is the last time written.
+            Alternatively, restarttime could be a time value (including np.datetime64) or
+            a callable function such as np.nanmin. The last is useful when running with dt < 0.
+        repeatdt :
+            Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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'
+        **kwargs :
+            Keyword arguments passed to the particleset constructor.
         """
         pass
 
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
         """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`
-        :param particle_val: Optional numpy-array of values to weigh each particle with,
-                             or string name of particle variable to use weigh particles with.
-                             Default is None, resulting in a value of 1 for each particle
-        :param relative: Boolean to control whether the density is scaled by the total
-                         weight of all particles. Default is False
-        :param area_scale: Boolean to control whether the density is scaled by the area
-                           (in m^2) of each grid cell. Default is False
+        Parameters
+        ----------
+        field_name : str, optional
+            Name of the field from the fieldset to calculate the histogram on.
+            Defaults to using "U".
+        particle_val :
+            Optional numpy-array of values to weigh each particle with,
+            or string name of particle variable to use weigh particles with.
+            Default is None, resulting in a value of 1 for each particle
+        relative : bool
+            Whether the density is scaled by the total weight of all particles.
+            Default is False
+        area_scale : bool
+            Whether the density is scaled by the area (in m^2) of each grid cell.
+            Default is False
         """
         pass
 
@@ -228,7 +289,15 @@ def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
         """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)
+
+        Parameters
+        ----------
+        delete_cfiles : bool
+            Whether to delete the C-files after compilation in JIT mode (default is True)
+        pyfunc :
+
+        c_include :
+             (Default value = "")
         """
         pass
 
@@ -246,8 +315,12 @@ def _set_particle_vector(self, name, value):
 
         This is a fallback implementation, it might be slow.
 
-        :param name: Name of the attribute (str).
-        :param value: New value to set the attribute of the particles to.
+        Parameters
+        ----------
+        name : str
+            Name of the attribute (str).
+        value :
+            New value to set the attribute of the particles to.
         """
         for p in self:
             setattr(p, name, value)
@@ -259,7 +332,12 @@ def error_particles(self):
 
         This is a fallback implementation, it might be slow.
 
-        :return: Collection iterator over error particles.
+
+        Returns
+        -------
+        iterator
+            Collection iterator over error particles.
+
         """
         error_indices = [
             i for i, p in enumerate(self)
@@ -277,8 +355,16 @@ def _impute_release_times(self, default):
 
         This is a fallback implementation, it might be slow.
 
-        :param default: Default release time.
-        :return: Minimum and maximum release times.
+        Parameters
+        ----------
+        default :
+            Default release time.
+
+        Returns
+        -------
+        type
+            Minimum and maximum release times.
+
         """
         max_rt = None
         min_rt = None
@@ -299,29 +385,44 @@ def execute(self, pyfunc=AdvectionRK4, pyfunc_inter=None, endtime=None, runtime=
         Optionally also provide sub-timestepping
         for particle output.
 
-        :param pyfunc: Kernel function to execute. This can be the name of a
-                       defined Python function or a :class:`parcels.kernel.Kernel` object.
-                       Kernels can be concatenated using the + operator
-        :param endtime: End time for the timestepping loop.
-                        It is either a datetime object or a positive double.
-        :param runtime: Length of the timestepping loop. Use instead of endtime.
-                        It is either a timedelta object or a positive double.
-        :param dt: Timestep interval to be passed to the kernel.
-                   It is either a timedelta object or a double.
-                   Use a negative value for a backward-in-time simulation.
-        :param moviedt:  Interval for inner sub-timestepping (leap), which dictates
-                         the update frequency of animation.
-                         It is either a timedelta object or a positive double.
-                         None value means no animation.
-        :param output_file: :mod:`parcels.particlefile.ParticleFile` object for particle output
-        :param recovery: Dictionary with additional `:mod:parcels.tools.error`
-                         recovery kernels to allow custom recovery behaviour in case of
-                         kernel errors.
-        :param movie_background_field: field plotted as background in the movie if moviedt is set.
-                                       'vector' shows the velocity as a vector field.
-        :param verbose_progress: Boolean for providing a progress bar for the kernel execution loop.
-        :param postIterationCallbacks: (Optional) Array of functions that are to be called after each iteration (post-process, non-Kernel)
-        :param callbackdt: (Optional, in conjecture with 'postIterationCallbacks) timestep inverval to (latestly) interrupt the running kernel and invoke post-iteration callbacks from 'postIterationCallbacks'
+        Parameters
+        ----------
+        pyfunc :
+            Kernel function to execute. This can be the name of a
+            defined Python function or a :class:`parcels.kernel.Kernel` object.
+            Kernels can be concatenated using the + operator (Default value = AdvectionRK4)
+        endtime :
+            End time for the timestepping loop.
+            It is either a datetime object or a positive double. (Default value = None)
+        runtime :
+            Length of the timestepping loop. Use instead of endtime.
+            It is either a timedelta object or a positive double. (Default value = None)
+        dt :
+            Timestep interval to be passed to the kernel.
+            It is either a timedelta object or a double.
+            Use a negative value for a backward-in-time simulation. (Default value = 1.)
+        moviedt :
+            Interval for inner sub-timestepping (leap), which dictates
+            the update frequency of animation.
+            It is either a timedelta object or a positive double.
+            None value means no animation. (Default value = None)
+        output_file :
+            mod:`parcels.particlefile.ParticleFile` object for particle output (Default value = None)
+        recovery :
+            Dictionary with additional `:mod:parcels.tools.error`
+            recovery kernels to allow custom recovery behaviour in case of
+            kernel errors. (Default value = None)
+        movie_background_field :
+            field plotted as background in the movie if moviedt is set.
+            'vector' shows the velocity as a vector field. (Default value = None)
+        verbose_progress : bool
+            Boolean for providing a progress bar for the kernel execution loop. (Default value = None)
+        postIterationCallbacks :
+            Optional) Array of functions that are to be called after each iteration (post-process, non-Kernel) (Default value = None)
+        callbackdt :
+            Optional, in conjecture with 'postIterationCallbacks) timestep inverval to (latestly) interrupt the running kernel and invoke post-iteration callbacks from 'postIterationCallbacks' (Default value = None)
+        pyfunc_inter :
+             (Default value = None)
         """
         # check if pyfunc has changed since last compile. If so, recompile
         if self.kernel is None or (self.kernel.pyfunc is not pyfunc and self.kernel is not pyfunc):
@@ -515,16 +616,31 @@ def show(self, with_particles=True, show_time=None, field=None, domain=None, pro
              land=True, vmin=None, vmax=None, savefile=None, animation=False, **kwargs):
         """Method to 'show' a Parcels ParticleSet.
 
-        :param with_particles: Boolean whether to show particles
-        :param show_time: Time at which to show the ParticleSet
-        :param field: Field to plot under particles (either None, a Field object, or 'vector')
-        :param domain: dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show
-        :param projection: type of cartopy projection to use (default PlateCarree)
-        :param land: Boolean whether to show land. This is ignored for flat meshes
-        :param vmin: minimum colour scale (only in single-plot mode)
-        :param vmax: maximum colour scale (only in single-plot mode)
-        :param savefile: Name of a file to save the plot to
-        :param animation: Boolean whether result is a single plot, or an animation
+        Parameters
+        ----------
+        with_particles : bool
+            Whether to show particles (Default value = True)
+        show_time :
+            Time at which to show the ParticleSet (Default value = None)
+        field :
+            Field to plot under particles (either None, a Field object, or 'vector') (Default value = None)
+        domain :
+            dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show (Default value = None)
+        projection :
+            type of cartopy projection to use (default PlateCarree)
+        land : bool
+            Whether to show land. This is ignored for flat meshes (Default value = True)
+        vmin : float
+            minimum colour scale (only in single-plot mode) (Default value = None)
+        vmax : float
+            maximum colour scale (only in single-plot mode) (Default value = None)
+        savefile : str
+            Name of a file to save the plot to (Default value = None)
+        animation : bool
+            Whether result is a single plot, or an animation (Default value = False)
+        **kwargs :
+            Keyword arguments passed to the :func:`parcels.plotting.plotparticles`.
+
         """
         from parcels.plotting import plotparticles
         plotparticles(particles=self, with_particles=with_particles, show_time=show_time, field=field, domain=domain,
diff --git a/parcels/particleset/particlesetaos.py b/parcels/particleset/particlesetaos.py
index c2e4a2df3..802ce1f46 100644
--- a/parcels/particleset/particlesetaos.py
+++ b/parcels/particleset/particlesetaos.py
@@ -42,23 +42,34 @@ def _convert_to_reltime(time):
 class ParticleSetAOS(BaseParticleSet):
     """Container class for storing particle and executing kernel over them.
 
-    :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity.
-           While fieldset=None is supported, this will throw a warning as it breaks most Parcels functionality
-    :param pclass: Optional :mod:`parcels.particle.JITParticle` or
-                 :mod:`parcels.particle.ScipyParticle` 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
-    :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0]
-    :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-    :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-           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'
-    :param pid_orig: Optional list of (offsets for) the particle IDs
-    :param partitions: List of cores on which to distribute the particles for MPI runs. Default: None, in which case particles
-           are distributed automatically on the processors
-
-    Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
+    Parameters
+    ----------
+    fieldset :
+        mod:`parcels.fieldset.FieldSet` object from which to sample velocity.
+        While fieldset=None is supported, this will throw a warning as it breaks most Parcels functionality
+    pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+        Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+    lon :
+        List of initial longitude values for particles
+    lat :
+        List of initial latitude values for particles
+    depth :
+        Optional list of initial depth values for particles. Default is 0m
+    time :
+        Optional list of initial time values for particles. Default is fieldset.U.grid.time[0]
+    repeatdt : datetime.timedelta or float, optional
+        Optional interval on which to repeat the release of the ParticleSet. Either timedelta object, or float in seconds.
+    lonlatdepth_dtype :
+        Floating precision for lon, lat, depth particle coordinates.
+        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'
+    pid_orig :
+        Optional list of (offsets for) the particle IDs
+    partitions :
+        List of cores on which to distribute the particles for MPI runs. Default: None, in which case particles
+        are distributed automatically on the processors
+
+        Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
     """
 
     def __init__(self, fieldset=None, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, pid_orig=None, **kwargs):
@@ -273,8 +284,7 @@ def __del__(self):
         super().__del__()
 
     def delete(self, key):
-        """
-        This is the generic super-method to indicate obejct deletion of a specific object from this collection.
+        """This is the generic super-method to indicate obejct deletion of a specific object from this collection.
 
         Comment/Annotation:
         Functions for deleting multiple objects are more specialised than just a for-each loop of single-item deletion,
@@ -290,8 +300,7 @@ def delete(self, key):
             self.delete_by_ID(key)
 
     def delete_by_index(self, index):
-        """
-        This method deletes a particle from the  the collection based on its index. It does not return the deleted item.
+        """This method deletes a particle from the  the collection based on its index. It does not return the deleted item.
         Semantically, the function appears similar to the 'remove' operation. That said, the function in OceanParcels -
         instead of directly deleting the particle - just raises the 'deleted' status flag for the indexed particle.
         In result, the particle still remains in the collection. The functional interpretation of the 'deleted' status
@@ -300,8 +309,7 @@ def delete_by_index(self, index):
         self._collection[index].state = OperationCode.Delete
 
     def delete_by_ID(self, id):
-        """
-        This method deletes a particle from the  the collection based on its ID. It does not return the deleted item.
+        """This method deletes a particle from the  the collection based on its ID. It does not return the deleted item.
         Semantically, the function appears similar to the 'remove' operation. That said, the function in OceanParcels -
         instead of directly deleting the particle - just raises the 'deleted' status flag for the indexed particle.
         In result, the particle still remains in the collection. The functional interpretation of the 'deleted' status
@@ -313,16 +321,28 @@ def delete_by_ID(self, id):
     def _set_particle_vector(self, name, value):
         """Set attributes of all particles to new values.
 
-        :param name: Name of the attribute (str).
-        :param value: New value to set the attribute of the particles to.
+        Parameters
+        ----------
+        name : str
+            Name of the attribute.
+        value :
+            New value to set the attribute of the particles to.
         """
         [setattr(p, name, value) for p in self._collection.data]
 
     def _impute_release_times(self, default):
         """Set attribute 'time' to default if encountering NaN values.
 
-        :param default: Default release time.
-        :return: Minimum and maximum release times.
+        Parameters
+        ----------
+        default :
+            Default release time.
+
+        Returns
+        -------
+        type
+            Minimum and maximum release times.
+
         """
         max_rt = None
         min_rt = None
@@ -339,12 +359,21 @@ 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`.
 
-        :param variable_name: Name of the variable to check.
-        :param compare_values: Value or list of values to compare to.
-        :param invert: Whether to invert the selection. I.e., when True,
-                       return all indices that do not equal (one of)
-                       `compare_values`.
-        :return: Numpy array of indices that satisfy the test.
+        Parameters
+        ----------
+        variable_name : str
+            Name of the variable to check.
+        compare_values :
+            Value or list of values to compare to.
+        invert :
+            Whether to invert the selection. I.e., when True,
+            return all indices that do not equal (one of)
+            `compare_values`. (Default value = False)
+
+        Returns
+        -------
+        np.ndarray
+            Numpy array of indices that satisfy the test.
         """
         compare_values = np.array([compare_values, ]) if type(compare_values) not in [list, dict, np.ndarray] else compare_values
         result = []
@@ -373,7 +402,10 @@ def particle_data(self):
     def error_particles(self):
         """Get an iterator over all particles that are in an error state.
 
-        :return: Collection iterator over error particles.
+        Returns
+        -------
+        iterator
+            Collection iterator over error particles.
         """
         error_indices = [
             i for i, p in enumerate(self)
@@ -384,7 +416,10 @@ def error_particles(self):
     def num_error_particles(self):
         """Get the number of particles that are in an error state.
 
-        :return: The number of error particles.
+        Returns
+        -------
+        int
+            Number of error particles.
         """
         return np.sum([True for p in self._collection if p.state not in [StateCode.Success, StateCode.Evaluate]])
 
@@ -402,7 +437,10 @@ def __getattr__(self, name):
         """
         Access a single property of all particles.
 
-        :param name: name of the property
+        Parameters
+        ----------
+        name : str
+            Name of the property to access.
         """
         for v in self._collection.ptype.variables:
             if v.name == name:
@@ -449,12 +487,25 @@ def lonlatdepth_dtype_from_field_interp_method(field):
 
     @classmethod
     def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'):
-        """
-        Converts a starting field into a monte-carlo sample of lons and lats.
+        """Converts a starting field into a monte-carlo sample of lons and lats.
+
+        Parameters
+        ----------
+        start_field :
+            mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
+
+            returns list(lon), list(lat)
+        size :
 
-        :param start_field: :mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
+        mode :
+             (Default value = 'monte_carlo')
 
-        returns list(lon), list(lat)
+        Returns
+        -------
+        list of float
+            A list of longitude values.
+        list of float
+            A list of latitude values.
         """
         if mode == 'monte_carlo':
             data = start_field.data if isinstance(start_field.data, np.ndarray) else np.array(start_field.data)
@@ -500,19 +551,29 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
         This creates a new ParticleSet based on locations of all particles written
         in a zarr ParticleFile at a certain time. Particle IDs are preserved if restart=True
 
-        :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
-        :param filename: Name of the particlefile from which to read initial conditions
-        :param restart: Boolean to signal if pset is used for a restart (default is True).
-               In that case, Particle IDs are preserved.
-        :param restarttime: time at which the Particles will be restarted. Default is the last time written.
-               Alternatively, restarttime could be a time value (including np.datetime64) or
-               a callable function such as np.nanmin. The last is useful when running with dt < 0.
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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'
+        Parameters
+        ----------
+        fieldset :
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        filename : str
+            Name of the particlefile from which to read initial conditions
+        restart : bool
+            BSignal if pset is used for a restart (default is True).
+            In that case, Particle IDs are preserved.
+        restarttime :
+            time at which the Particles will be restarted. Default is the last time written.
+            Alternatively, restarttime could be a time value (including np.datetime64) or
+            a callable function such as np.nanmin. The last is useful when running with dt < 0.
+        repeatdt : datetime.timedelta or float, optional
+            Optional interval on which to repeat the release of the ParticleSet. Either timedelta object, or float in seconds.
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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'
+        **kwargs :
+            Keyword arguments passed to the particleset constructor.
         """
         if repeatdt is not None:
             logger.warning(f'Note that the `repeatdt` argument is not retained from {filename}, and that '
@@ -563,13 +624,21 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
                    lonlatdepth_dtype=lonlatdepth_dtype, repeatdt=repeatdt, **kwargs)
 
     def __iadd__(self, particles):
-        """Add particles to the ParticleSet. Note that this is an
-        incremental add, the particles will be added to the ParticleSet
+        """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.
 
-        :param particles: Another ParticleSet containing particles to add
-                          to this one.
-        :return: The current ParticleSet
+        Parameters
+        ----------
+        particles :
+            Another ParticleSet containing particles to add
+            to this one.
+
+        Returns
+        -------
+        type
+            The current ParticleSet
         """
         self.add(particles)
         return self
@@ -579,9 +648,17 @@ def add(self, particles):
         incremental add, the particles will be added to the ParticleSet
         on which this function is called.
 
-        :param particles: Another ParticleSet containing particles to add
-                          to this one.
-        :return: The current ParticleSet
+        Parameters
+        ----------
+        particles :
+            Another ParticleSet containing particles to add
+            to this one.
+
+        Returns
+        -------
+        type
+            The current ParticleSet
+
         """
         if isinstance(particles, BaseParticleSet):
             particles = particles.collection
@@ -604,15 +681,21 @@ def remove_booleanvector(self, indices):
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
         """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`
-        :param particle_val: Optional numpy-array of values to weigh each particle with,
-                             or string name of particle variable to use weigh particles with.
-                             Default is None, resulting in a value of 1 for each particle
-        :param relative: Boolean to control whether the density is scaled by the total
-                         weight of all particles. Default is False
-        :param area_scale: Boolean to control whether the density is scaled by the area
-                           (in m^2) of each grid cell. Default is False
+        Parameters
+        ----------
+        field_name : str, optional
+            Name of the field from the fieldset to calculate the histogram on.
+            Defaults to using "U".
+        particle_val :
+            Optional numpy-array of values to weigh each particle with,
+            or string name of particle variable to use weigh particles with.
+            Default is None, resulting in a value of 1 for each particle
+        relative : bool
+            Whether the density is scaled by the total weight of all particles.
+            Default is False
+        area_scale : bool
+            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
@@ -658,7 +741,10 @@ 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.
 
-        :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
+        Parameters
+        ----------
+        delete_cfiles : bool
+            Whether to delete the C-files after compilation in JIT mode (default is True).
         """
         return KernelAOS(self.fieldset, self.collection.ptype, pyfunc=pyfunc, c_include=c_include, delete_cfiles=delete_cfiles)
 
@@ -669,9 +755,13 @@ def ParticleFile(self, *args, **kwargs):
         return ParticleFileAOS(*args, particleset=self, **kwargs)
 
     def set_variable_write_status(self, var, write_status):
-        """
-        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')
+        """Method to set the write status of a Variable.
+
+        Parameters
+        ----------
+        var :
+            Name of the variable (string)
+        write_status :
+            Write status of the variable (True, False or 'once')
         """
         self._collection.set_variable_write_status(var, write_status)
diff --git a/parcels/particleset/particlesetsoa.py b/parcels/particleset/particlesetsoa.py
index 3b18cb3bb..2dab5120e 100644
--- a/parcels/particleset/particlesetsoa.py
+++ b/parcels/particleset/particlesetsoa.py
@@ -53,25 +53,38 @@ class ParticleSetSOA(BaseParticleSet):
     and individual particles can only be deleted as a set procedually (i.e. by 'particle.delete()'-call during
     kernel execution).
 
-    :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity.
-           While fieldset=None is supported, this will throw a warning as it breaks most Parcels functionality
-    :param pclass: Optional :mod:`parcels.particle.JITParticle` or
-                 :mod:`parcels.particle.ScipyParticle` 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
-    :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0]
-    :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-    :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-           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'
-    :param pid_orig: Optional list of (offsets for) the particle IDs
-    :param partitions: List of cores on which to distribute the particles for MPI runs. Default: None, in which case particles
-           are distributed automatically on the processors
-    :param periodic_domain_zonal: Zonal domain size, used to apply zonally periodic boundaries for particle-particle
-           interaction. If None, no zonally periodic boundaries are applied
-
-    Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
+    Parameters
+    ----------
+    fieldset :
+        mod:`parcels.fieldset.FieldSet` object from which to sample velocity.
+        While fieldset=None is supported, this will throw a warning as it breaks most Parcels functionality
+    pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+        Optional :mod:`parcels.particle.JITParticle` or
+        :mod:`parcels.particle.ScipyParticle` object that defines custom particle
+    lon :
+        List of initial longitude values for particles
+    lat :
+        List of initial latitude values for particles
+    depth :
+        Optional list of initial depth values for particles. Default is 0m
+    time :
+        Optional list of initial time values for particles. Default is fieldset.U.grid.time[0]
+    repeatdt : datetime.timedelta or float, optional
+        Optional interval on which to repeat the release of the ParticleSet. Either timedelta object, or float in seconds.
+    lonlatdepth_dtype :
+        Floating precision for lon, lat, depth particle coordinates.
+        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'
+    pid_orig :
+        Optional list of (offsets for) the particle IDs
+    partitions :
+        List of cores on which to distribute the particles for MPI runs. Default: None, in which case particles
+        are distributed automatically on the processors
+    periodic_domain_zonal :
+        Zonal domain size, used to apply zonally periodic boundaries for particle-particle
+        interaction. If None, no zonally periodic boundaries are applied
+
+        Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
     """
 
     def __init__(self, fieldset=None, pclass=JITParticle, lon=None, lat=None,
@@ -248,16 +261,28 @@ def __del__(self):
     def _set_particle_vector(self, name, value):
         """Set attributes of all particles to new values.
 
-        :param name: Name of the attribute (str).
-        :param value: New value to set the attribute of the particles to.
+        Parameters
+        ----------
+        name : str
+            Name of the attribute (str).
+        value : any
+            New value to set the attribute of the particles to.
         """
         self.collection._data[name][:] = value
 
     def _impute_release_times(self, default):
         """Set attribute 'time' to default if encountering NaN values.
 
-        :param default: Default release time.
-        :return: Minimum and maximum release times.
+        Parameters
+        ----------
+        default :
+            Default release time.
+
+        Returns
+        -------
+        type
+            Minimum and maximum release times.
+
         """
         if np.any(np.isnan(self._collection.data['time'])):
             self._collection.data['time'][np.isnan(self._collection.data['time'])] = default
@@ -266,12 +291,22 @@ def _impute_release_times(self, default):
     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`.
 
-        :param variable_name: Name of the variable to check.
-        :param compare_values: Value or list of values to compare to.
-        :param invert: Whether to invert the selection. I.e., when True,
-                       return all indices that do not equal (one of)
-                       `compare_values`.
-        :return: Numpy array of indices that satisfy the test.
+        Parameters
+        ----------
+        variable_name : str
+            Name of the variable to check.
+        compare_values :
+            Value or list of values to compare to.
+        invert :
+            Whether to invert the selection. I.e., when True,
+            return all indices that do not equal (one of)
+            `compare_values`. (Default value = False)
+
+        Returns
+        -------
+        np.ndarray
+            Numpy array of indices that satisfy the test.
+
         """
         compare_values = np.array([compare_values, ]) if type(compare_values) not in [list, dict, np.ndarray] else compare_values
         return np.where(np.isin(self._collection.data[variable_name], compare_values, invert=invert))[0]
@@ -312,7 +347,10 @@ def populate_indices(self):
     def error_particles(self):
         """Get an iterator over all particles that are in an error state.
 
-        :return: Collection iterator over error particles.
+        Returns
+        -------
+        iterator
+            Collection iterator over error particles.
         """
         error_indices = self.data_indices('state', [StateCode.Success, StateCode.Evaluate], invert=True)
         return ParticleCollectionIterableSOA(self._collection, subset=error_indices)
@@ -328,7 +366,10 @@ def active_particles_mask(self, time, dt):
     def num_error_particles(self):
         """Get the number of particles that are in an error state.
 
-        :return: The number of error particles.
+        Returns
+        -------
+        int
+            Number of error particles.
         """
         return np.sum(np.isin(
             self._collection.data['state'],
@@ -339,8 +380,7 @@ def __getitem__(self, index):
         return self._collection.get_single_by_index(index)
 
     def cstruct(self):
-        """
-        'cstruct' returns the ctypes mapping of the combined collections cstruct and the fieldset cstruct.
+        """'cstruct' returns the ctypes mapping of the combined collections cstruct and the fieldset cstruct.
         This depends on the specific structure in question.
         """
         cstruct = self._collection.cstruct()
@@ -352,12 +392,24 @@ def ctypes_struct(self):
 
     @classmethod
     def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'):
-        """
-        Converts a starting field into a monte-carlo sample of lons and lats.
+        """Converts a starting field into a monte-carlo sample of lons and lats.
+
+        Parameters
+        ----------
+        start_field :
+            mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
+        size :
 
-        :param start_field: :mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally)  according to the presented density field.
+        mode :
+             (Default value = 'monte_carlo')
 
-        returns list(lon), list(lat)
+
+        Returns
+        -------
+        list of float
+            A list of longitude values.
+        list of float
+            A list of latitude values.
         """
         if mode == 'monte_carlo':
             data = start_field.data if isinstance(start_field.data, np.ndarray) else np.array(start_field.data)
@@ -403,19 +455,29 @@ def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime
         This creates a new ParticleSet based on locations of all particles written
         in a zarr ParticleFile at a certain time. Particle IDs are preserved if restart=True
 
-        :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
-        :param filename: Name of the particlefile from which to read initial conditions
-        :param restart: Boolean to signal if pset is used for a restart (default is True).
-               In that case, Particle IDs are preserved.
-        :param restarttime: time at which the Particles will be restarted. Default is the last time written.
-               Alternatively, restarttime could be a time value (including np.datetime64) or
-               a callable function such as np.nanmin. The last is useful when running with dt < 0.
-        :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
-        :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
-               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'
+        Parameters
+        ----------
+        fieldset : parcels.fieldset.FieldSet
+            mod:`parcels.fieldset.FieldSet` object from which to sample velocity
+        pclass : parcels.particle.JITParticle or parcels.particle.ScipyParticle
+            Particle class. May be a particle class as defined in parcels, or a subclass defining a custom particle.
+        filename : str
+            Name of the particlefile from which to read initial conditions
+        restart : bool
+            BSignal if pset is used for a restart (default is True).
+            In that case, Particle IDs are preserved.
+        restarttime :
+            time at which the Particles will be restarted. Default is the last time written.
+            Alternatively, restarttime could be a time value (including np.datetime64) or
+            a callable function such as np.nanmin. The last is useful when running with dt < 0.
+        repeatdt : datetime.timedelta or float, optional
+            Optional interval on which to repeat the release of the ParticleSet. Either timedelta object, or float in seconds.
+        lonlatdepth_dtype :
+            Floating precision for lon, lat, depth particle coordinates.
+            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'
+        **kwargs :
+            Keyword arguments passed to the particleset constructor.
         """
         if repeatdt is not None:
             logger.warning(f'Note that the `repeatdt` argument is not retained from {filename}, and that '
@@ -512,13 +574,19 @@ def __sizeof__(self):
     def __iadd__(self, particles):
         """Add particles to the ParticleSet.
 
-        Note that this is an
-        incremental add, the particles will be added to the ParticleSet
+        Note that this is an incremental add, the particles will be added to the ParticleSet
         on which this function is called.
 
-        :param particles: Another ParticleSet containing particles to add
-                          to this one.
-        :return: The current ParticleSet
+        Parameters
+        ----------
+        particles :
+            Another ParticleSet containing particles to add
+            to this one.
+
+        Returns
+        -------
+        type
+            The current ParticleSet
         """
         self.add(particles)
         return self
@@ -534,8 +602,16 @@ def add(self, particles):
         incremental add, the particles will be added to the ParticleSet
         on which this function is called.
 
-        :param particles: Another ParticleSet containing particles to add to this one.
-        :return: The current ParticleSet
+        Parameters
+        ----------
+        particles :
+            Another ParticleSet containing particles to add to this one.
+
+        Returns
+        -------
+        type
+            The current ParticleSet
+
         """
         if isinstance(particles, BaseParticleSet):
             particles = particles.collection
@@ -562,15 +638,21 @@ def remove_booleanvector(self, indices):
     def density(self, field_name=None, particle_val=None, relative=False, area_scale=False):
         """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`
-        :param particle_val: Optional numpy-array of values to weigh each particle with,
-                             or string name of particle variable to use weigh particles with.
-                             Default is None, resulting in a value of 1 for each particle
-        :param relative: Boolean to control whether the density is scaled by the total
-                         weight of all particles. Default is False
-        :param area_scale: Boolean to control whether the density is scaled by the area
-                           (in m^2) of each grid cell. Default is False
+        Parameters
+        ----------
+        field_name : str, optional
+            Name of the field from the fieldset to calculate the histogram on.
+            Defaults to using "U".
+        particle_val :
+            Optional numpy-array of values to weigh each particle with,
+            or string name of particle variable to use weigh particles with.
+            Default is None, resulting in a value of 1 for each particle
+        relative : bool
+            Whether the density is scaled by the total weight of all particles.
+            Default is False
+        area_scale : bool
+            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
@@ -616,7 +698,15 @@ def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
         """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)
+
+        Parameters
+        ----------
+        delete_cfiles : bool
+            Whether to delete the C-files after compilation in JIT mode (default is True)
+        pyfunc :
+
+        c_include :
+             (Default value = "")
         """
         return Kernel(self.fieldset, self.collection.ptype, pyfunc=pyfunc, c_include=c_include,
                       delete_cfiles=delete_cfiles)
@@ -631,10 +721,13 @@ def ParticleFile(self, *args, **kwargs):
         return ParticleFile(*args, particleset=self, **kwargs)
 
     def set_variable_write_status(self, var, write_status):
-        """
-        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')
+        """Method to set the write status of a Variable.
+
+        Parameters
+        ----------
+        var :
+            Name of the variable (string)
+        write_status :
+            Write status of the variable (True, False or 'once')
         """
         self._collection.set_variable_write_status(var, write_status)
diff --git a/parcels/plotting.py b/parcels/plotting.py
index dcc6eb590..3fba20c2a 100644
--- a/parcels/plotting.py
+++ b/parcels/plotting.py
@@ -15,16 +15,32 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do
                   animation=False, **kwargs):
     """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
-    :param field: Field to plot under particles (either None, a Field object, or 'vector')
-    :param domain: dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show
-    :param projection: type of cartopy projection to use (default PlateCarree)
-    :param land: Boolean whether to show land. This is ignored for flat meshes
-    :param vmin: minimum colour scale (only in single-plot mode)
-    :param vmax: maximum colour scale (only in single-plot mode)
-    :param savefile: Name of a file to save the plot to
-    :param animation: Boolean whether result is a single plot, or an animation
+    Parameters
+    ----------
+    show_time :
+        Time in seconds from start after which to show the ParticleSet (Default value = None)
+    with_particles : bool
+        Whether particles are also plotted on Field (Default value = True)
+    field : parcels.Field or str, optional
+        Field to plot under particles (either None, a Field object, or 'vector') (Default value = None)
+    domain :
+        dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show (Default value = None)
+    projection :
+        type of cartopy projection to use (default PlateCarree)
+    land : bool
+        Whether to show land. This is ignored for flat meshes (Default value = True)
+    vmin : float
+        minimum colour scale (only in single-plot mode) (Default value = None)
+    vmax : float
+        maximum colour scale (only in single-plot mode) (Default value = None)
+    savefile :
+        Name of a file to save the plot to (Default value = None)
+    animation : bool
+        Whether result is a single plot, or an animation (Default value = False)
+    particles :
+
+    **kwargs :
+        Keyword arguments passed to the :func:`plotfield` function.
     """
     show_time = particles[0].time if show_time is None else show_time
     if isinstance(show_time, datetime):
@@ -101,15 +117,30 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
               vmin=None, vmax=None, savefile=None, **kwargs):
     """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
-    :param depth_level: depth level to be plotted (default 0)
-    :param projection: type of cartopy projection to use (default PlateCarree)
-    :param land: Boolean whether to show land. This is ignored for flat meshes
-    :param vmin: minimum colour scale (only in single-plot mode)
-    :param vmax: maximum colour scale (only in single-plot mode)
-    :param savefile: Name of a file to save the plot to
-    :param animation: Boolean whether result is a single plot, or an animation
+    Parameters
+    ----------
+    field : parcels.field.Field
+        Field to plot
+    show_time : float
+        Time in seconds from start after which to show the Field (Default value = None)
+    domain :
+        dictionary (with keys 'N', 'S', 'E', 'W') defining domain to show (Default value = None)
+    depth_level :
+        depth level to be plotted (default 0)
+    projection :
+        type of cartopy projection to use (default PlateCarree)
+    land : bool
+        Whether to show land. This is ignored for flat meshes (Default value = True)
+    vmin : float, optional
+        minimum colour scale (only in single-plot mode) (Default value = None)
+    vmax : float, optional
+        maximum colour scale (only in single-plot mode) (Default value = None)
+    savefile : str, optional
+        Name of a file to save the plot to (Default value = None)
+    animation : bool
+        Whether result is a single plot, or an animation
+    **kwargs :
+        Provide "titlestr" as the title of the plot, or "cartopy_features" to be used in :func:`create_parcelsfig_axis`.
     """
     if type(field) is VectorField:
         spherical = True if field.U.grid.mesh == 'spherical' else False
diff --git a/parcels/scripts/plottrajectoriesfile.py b/parcels/scripts/plottrajectoriesfile.py
index c3ce043e1..f79235721 100644
--- a/parcels/scripts/plottrajectoriesfile.py
+++ b/parcels/scripts/plottrajectoriesfile.py
@@ -19,20 +19,33 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
                          bins=20, show_plt=True, central_longitude=0):
     """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',
-                'movie2d' and 'movie2d_notebook'. The latter two give animations,
-                with 'movie2d_notebook' specifically designed for jupyter notebooks
-    :param tracerfile: Name of NetCDF file to show as background
-    :param tracerfield: Name of variable to show as background
-    :param tracerlon: Name of longitude dimension of variable to show as background
-    :param tracerlat: Name of latitude dimension of variable to show as background
-    :param recordedvar: Name of variable used to color particles in scatter-plot.
-                Only works in 'movie2d' or 'movie2d_notebook' mode.
-    :param movie_forward: Boolean whether to show movie in forward or backward mode (default True)
-    :param bins: Number of bins to use in `hist2d` mode. See also https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist2d.html
-    :param show_plt: Boolean whether plot should directly be show (for py.test)
-    :param central_longitude: Degrees East at which to center the plot
+    Parameters
+    ----------
+    filename : str
+        Name of Parcels-generated NetCDF file with particle positions
+    mode :
+        Type of plot to show. Supported are '2d', '3d', 'hist2d',
+        'movie2d' and 'movie2d_notebook'. The latter two give animations,
+        with 'movie2d_notebook' specifically designed for jupyter notebooks (Default value = '2d')
+    tracerfile :
+        Name of NetCDF file to show as background (Default value = None)
+    tracerfield :
+        Name of variable to show as background (Default value = 'P')
+    tracerlon :
+        Name of longitude dimension of variable to show as background (Default value = 'x')
+    tracerlat :
+        Name of latitude dimension of variable to show as background (Default value = 'y')
+    recordedvar :
+        Name of variable used to color particles in scatter-plot.
+        Only works in 'movie2d' or 'movie2d_notebook' mode. (Default value = None)
+    movie_forward : bool
+        Whether to show movie in forward or backward mode (default True)
+    bins :
+        Number of bins to use in `hist2d` mode. See also https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist2d.html (Default value = 20)
+    show_plt : bool
+        Whether plot should directly be shown (for py.test) (Default value = True)
+    central_longitude :
+        Degrees East at which to center the plot (Default value = 0)
     """
     environ["HDF5_USE_FILE_LOCKING"] = "FALSE"
     try:
diff --git a/parcels/tools/converters.py b/parcels/tools/converters.py
index c23b1b679..b9cee1581 100644
--- a/parcels/tools/converters.py
+++ b/parcels/tools/converters.py
@@ -15,7 +15,10 @@
 def convert_to_flat_array(var):
     """Convert lists and single integers/floats to one-dimensional numpy arrays
 
-    :param var: list or numeric to convert to a one-dimensional numpy array
+    Parameters
+    ----------
+    var : np.ndarray, float or array_like
+        list or numeric to convert to a one-dimensional numpy array
     """
     if isinstance(var, np.ndarray):
         return var.flatten()
@@ -37,10 +40,12 @@ def _get_cftime_calendars():
 
 
 class TimeConverter:
-    """ Converter class for dates with different calendars in FieldSets
+    """Converter class for dates with different calendars in FieldSets
 
-    :param: time_origin: time origin of the class. Currently supported formats are
-            float, integer, numpy.datetime64, and netcdftime.DatetimeNoLeap
+    Parameters
+    ----------
+    time_origin : float, integer, numpy.datetime64 or netcdftime.DatetimeNoLeap
+        time origin of the class.
     """
 
     def __init__(self, time_origin=0):
@@ -58,8 +63,16 @@ def reltime(self, time):
         """Method to compute the difference, in seconds, between a time and the time_origin
         of the TimeConverter
 
-        :param: time: input time
-        :return: time - self.time_origin
+        Parameters
+        ----------
+        time :
+
+
+        Returns
+        -------
+        type
+            time - self.time_origin
+
         """
         time = time.time_origin if isinstance(time, TimeConverter) else time
         if self.calendar in ['np_datetime64', 'np_timedelta64']:
@@ -85,8 +98,16 @@ def reltime(self, time):
     def fulltime(self, time):
         """Method to convert a time difference in seconds to a date, based on the time_origin
 
-        :param: time: input time
-        :return: self.time_origin + time
+        Parameters
+        ----------
+        time :
+
+
+        Returns
+        -------
+        type
+            self.time_origin + time
+
         """
         time = time.time_origin if isinstance(time, TimeConverter) else time
         if self.calendar in ['np_datetime64', 'np_timedelta64']:
@@ -132,9 +153,8 @@ def __le__(self, other):
 
 
 class UnitConverter:
-    """ Interface class for spatial unit conversion during field sampling
-        that performs no conversion.
-    """
+    """Interface class for spatial unit conversion during field sampling that performs no conversion."""
+
     source_unit = None
     target_unit = None
 
@@ -152,7 +172,8 @@ def ccode_to_source(self, x, y, z):
 
 
 class Geographic(UnitConverter):
-    """ Unit converter from geometric to geographic coordinates (m to degree) """
+    """Unit converter from geometric to geographic coordinates (m to degree)"""
+
     source_unit = 'm'
     target_unit = 'degree'
 
@@ -170,9 +191,10 @@ def ccode_to_source(self, x, y, z):
 
 
 class GeographicPolar(UnitConverter):
-    """ Unit converter from geometric to geographic coordinates (m to degree)
-        with a correction to account for narrower grid cells closer to the poles.
+    """Unit converter from geometric to geographic coordinates (m to degree)
+    with a correction to account for narrower grid cells closer to the poles.
     """
+
     source_unit = 'm'
     target_unit = 'degree'
 
@@ -190,7 +212,8 @@ def ccode_to_source(self, x, y, z):
 
 
 class GeographicSquare(UnitConverter):
-    """ Square distance converter from geometric to geographic coordinates (m2 to degree2) """
+    """Square distance converter from geometric to geographic coordinates (m2 to degree2)"""
+
     source_unit = 'm2'
     target_unit = 'degree2'
 
@@ -208,9 +231,10 @@ def ccode_to_source(self, x, y, z):
 
 
 class GeographicPolarSquare(UnitConverter):
-    """ Square distance converter from geometric to geographic coordinates (m2 to degree2)
-        with a correction to account for narrower grid cells closer to the poles.
+    """Square distance converter from geometric to geographic coordinates (m2 to degree2)
+    with a correction to account for narrower grid cells closer to the poles.
     """
+
     source_unit = 'm2'
     target_unit = 'degree2'
 
@@ -233,8 +257,7 @@ def ccode_to_source(self, x, y, z):
 
 
 def convert_xarray_time_units(ds, time):
-    """ Fixes DataArrays that have time.Unit instead of expected time.units
-    """
+    """Fixes DataArrays that have time.Unit instead of expected time.units"""
     da = ds[time] if isinstance(ds, xr.Dataset) else ds
     if 'units' not in da.attrs and 'Unit' in da.attrs:
         da.attrs['units'] = da.attrs['Unit']
diff --git a/setup.cfg b/setup.cfg
index db7013414..3c6c9a449 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,6 +30,12 @@ ignore =
     # ‘from module import *’ used; unable to detect undefined names
     F403,
 
+[pydocstyle]
+ignore =
+    # Numpy docstring format http://www.pydocstyle.org/en/latest/error_codes.html#default-conventions
+    D107, D203, D212, D213, D402, D413, D415, D416, D417,
+
+    # TODO: once vvv is fixed, replace `ignore = ` section with `convention = numpy`
     # IGNORE FOR NOW (requires more work)
     # Missing docstring in public module
     D100,
@@ -44,17 +50,18 @@ ignore =
     # 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,
+    # First word of the docstring should not be `This`
+    D404,
     # 1 blank line required between summary line and description (requires writing of summaries)
-    D205,
+    D205
 
 [tool:pytest]
 python_files = test_*.py example_*.py *tutorial*
 
+
 [isort]
 profile = black
 skip_gitignore = True


Test patch:
diff --git a/tests/test_data/create_testfields.py b/tests/test_data/create_testfields.py
index 1eaba9430..f52c8d773 100644
--- a/tests/test_data/create_testfields.py
+++ b/tests/test_data/create_testfields.py
@@ -73,8 +73,14 @@ def generate_perlin_testfield():
 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.
+    Parameters
+    ----------
+    field : parcels.field.Field
+        Field to write to file
+    filename : str
+        Base name of the file to write to
+    varname : str, optional
+        Name of the variable to write to file. If None, defaults to field.name
     """
     filepath = str(f'{filename}{field.name}.nc')
     if varname is None: