Repository: OceanParcels/parcels
Pull request #393 (2018-06-29 13:26:24)

List of velocity fields for advection
This PR introduced a new feature to Parcels: `FieldLists`.  With these `FieldLists`, you can run a kernel over multiple `Fields` and the values of the individual `Field` interpolations are automatically added. Note that the `Fields` can be on different `Grids`

This is useful for example if particles need to be advected with a combination of different velocity fields such as ocean currents, windage, etc). See also #326, where @willirath suggested such a feature

If the velocity fields are `(U_ocean, V_ocean)` for the currents and `(U_stokes, V_stokes)` for Stokes drift, then they can be combined as
```
fset = FieldSet([U_ocean, U_stokes], [V_ocean, V_stokes])
```
And then simply use the built-in `AdvectionRK4` kernel as normal

Lists can be of arbitrary lengths. This also works in 3D and for fields other than velocity `Fields`.

I'd really like some feedback on whether this works. Note that before merging a few minor things need to be fixed
1) Confirmation that the implementation also works on `mesh=spherical` and Curvilinear `Grids` (now tested only for `mesh=flat`)
2) Some under-the-hood cleaning of the `fieldset.ugrid` attribute, which is a rather hacky implementation
3) Create a tutorial to explain how to use these `FieldLists`. This tutorial can then also explain how to work with different `Grids`



Patch:
diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py
index 843db7b0..88b9e369 100644
--- a/parcels/codegenerator.py
+++ b/parcels/codegenerator.py
@@ -1,4 +1,5 @@
 from parcels.field import Field, VectorField
+from parcels.fieldset import FieldList, VectorFieldList
 from parcels.loggers import logger
 import ast
 import cgen as c
@@ -20,6 +21,12 @@ class FieldSetNode(IntrinsicNode):
         if isinstance(getattr(self.obj, attr), Field):
             return FieldNode(getattr(self.obj, attr),
                              ccode="%s->%s" % (self.ccode, attr))
+        elif isinstance(getattr(self.obj, attr), VectorFieldList):
+            return VectorFieldListNode(getattr(self.obj, attr),
+                                       ccode="%s->%s" % (self.ccode, attr))
+        elif isinstance(getattr(self.obj, attr), FieldList) or isinstance(getattr(self.obj, attr), list):
+            return FieldListNode(getattr(self.obj, attr),
+                                 ccode="%s->%s" % (self.ccode, attr))
         elif isinstance(getattr(self.obj, attr), VectorField):
             return VectorFieldNode(getattr(self.obj, attr),
                                    ccode="%s->%s" % (self.ccode, attr))
@@ -54,6 +61,32 @@ class VectorFieldEvalNode(IntrinsicNode):
         self.var3 = var3  # third variable for UVW interpolation
 
 
+class FieldListNode(IntrinsicNode):
+    def __getitem__(self, attr):
+        return FieldListEvalNode(self.obj, attr)
+
+
+class FieldListEvalNode(IntrinsicNode):
+    def __init__(self, fields, args, var):
+        self.fields = fields
+        self.args = args
+        self.var = var  # the variable in which the interpolated field is written
+
+
+class VectorFieldListNode(IntrinsicNode):
+    def __getitem__(self, attr):
+        return VectorFieldListEvalNode(self.obj, attr)
+
+
+class VectorFieldListEvalNode(IntrinsicNode):
+    def __init__(self, field, args, var, var2, var3):
+        self.field = field
+        self.args = args
+        self.var = var  # the variable in which the interpolated field is written
+        self.var2 = var2  # second variable for UV interpolation
+        self.var3 = var3  # third variable for UVW interpolation
+
+
 class ConstNode(IntrinsicNode):
     def __getitem__(self, attr):
         return attr
@@ -178,7 +211,24 @@ class IntrinsicTransformer(ast.NodeTransformer):
 
         # If we encounter field evaluation we replace it with a
         # temporary variable and put the evaluation call on the stack.
-        if isinstance(node.value, FieldNode):
+        if isinstance(node.value, FieldListNode):
+            tmp = [self.get_tmp() for _ in node.value.obj]
+            # Insert placeholder node for field eval ...
+            self.stmt_stack += [FieldListEvalNode(node.value, node.slice, tmp)]
+            # .. and return the name of the temporary that will be populated
+            return ast.Name(id='+'.join(tmp))
+        elif isinstance(node.value, VectorFieldListNode):
+            tmp = [self.get_tmp() for _ in node.value.obj.U]
+            tmp2 = [self.get_tmp() for _ in node.value.obj.U]
+            tmp3 = [self.get_tmp() if node.value.obj.W else None for _ in node.value.obj.U]
+            # Insert placeholder node for field eval ...
+            self.stmt_stack += [VectorFieldListEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
+            # .. and return the name of the temporary that will be populated
+            if all(tmp3):
+                return ast.Tuple([ast.Name(id='+'.join(tmp)), ast.Name(id='+'.join(tmp2)), ast.Name(id='+'.join(tmp3))], ast.Load())
+            else:
+                return ast.Tuple([ast.Name(id='+'.join(tmp)), ast.Name(id='+'.join(tmp2))], ast.Load())
+        elif isinstance(node.value, FieldNode):
             tmp = self.get_tmp()
             # Insert placeholder node for field eval ...
             self.stmt_stack += [FieldEvalNode(node.value, node.slice, tmp)]
@@ -261,8 +311,6 @@ class KernelGenerator(ast.NodeVisitor):
         self.fieldset = fieldset
         self.ptype = ptype
         self.field_args = OrderedDict()
-        # Hack alert: JIT requires U field to update fieldset indexes
-        self.field_args['U'] = fieldset.U
         self.vector_field_args = OrderedDict()
         self.const_args = OrderedDict()
 
@@ -567,10 +615,25 @@ class KernelGenerator(ast.NodeVisitor):
         """Record intrinsic fields used in kernel"""
         self.field_args[node.obj.name] = node.obj
 
+    def visit_FieldListNode(self, node):
+        """Record intrinsic fields used in kernel"""
+        for fld in node.obj:
+            self.field_args[fld.name] = fld
+
     def visit_VectorFieldNode(self, node):
         """Record intrinsic fields used in kernel"""
         self.vector_field_args[node.obj.name] = node.obj
 
+    def visit_VectorFieldListNode(self, node):
+        """Record intrinsic fields used in kernel"""
+        for fld in node.obj.U:
+            self.field_args[fld.name] = fld
+        for fld in node.obj.V:
+            self.field_args[fld.name] = fld
+        if hasattr(node.obj, 'W') and node.obj.W:
+            for fld in node.obj.W:
+                self.field_args[fld.name] = fld
+
     def visit_ConstNode(self, node):
         self.const_args[node.ccode] = node.obj
 
@@ -600,6 +663,39 @@ class KernelGenerator(ast.NodeVisitor):
         node.ccode = c.Block([c.Assign("err", ccode_eval),
                               conv_stat, c.Statement("CHECKERROR(err)")])
 
+    def visit_FieldListEvalNode(self, node):
+        self.visit(node.fields)
+        self.visit(node.args)
+        cstat = []
+        for fld, var in zip(node.fields.obj, node.var):
+            ccode_eval = fld.ccode_eval(var, *node.args.ccode)
+            ccode_conv = fld.ccode_convert(*node.args.ccode)
+            conv_stat = c.Statement("%s *= %s" % (var, ccode_conv))
+            cstat += [c.Assign("err", ccode_eval), conv_stat, c.Statement("CHECKERROR(err)")]
+        node.ccode = c.Block(cstat)
+
+    def visit_VectorFieldListEvalNode(self, node):
+        self.visit(node.field)
+        self.visit(node.args)
+        cstat = []
+        if node.field.obj.W:
+            Wlist = node.field.obj.W
+        else:
+            Wlist = [None] * len(node.field.obj.U)
+        for U, V, W, var, var2, var3 in zip(node.field.obj.U, node.field.obj.V, Wlist, node.var, node.var2, node.var3):
+            vfld = VectorField(node.field.obj.name, U, V, W)
+            ccode_eval = vfld.ccode_eval(var, var2, var3, U, V, W, *node.args.ccode)
+            ccode_conv1 = U.ccode_convert(*node.args.ccode)
+            ccode_conv2 = V.ccode_convert(*node.args.ccode)
+            statements = [c.Statement("%s *= %s" % (var, ccode_conv1)),
+                          c.Statement("%s *= %s" % (var2, ccode_conv2))]
+            if var3:
+                ccode_conv3 = W.ccode_convert(*node.args.ccode)
+                statements.append(c.Statement("%s *= %s" % (var3, ccode_conv3)))
+            conv_stat = c.Block(statements)
+            cstat += [c.Assign("err", ccode_eval), conv_stat, c.Statement("CHECKERROR(err)")]
+        node.ccode = c.Block(cstat)
+
     def visit_Return(self, node):
         self.visit(node.value)
         node.ccode = c.Statement('return %s' % node.value.ccode)
diff --git a/parcels/examples/example_moving_eddies.py b/parcels/examples/example_moving_eddies.py
index b1f7c028..95fc3969 100644
--- a/parcels/examples/example_moving_eddies.py
+++ b/parcels/examples/example_moving_eddies.py
@@ -152,6 +152,10 @@ def test_moving_eddies_file(fieldsetfile, mode):
 def test_periodic_and_computeTimeChunk_eddies(mode):
     filename = path.join(path.dirname(__file__), 'MovingEddies_data', 'moving_eddies')
     fieldset = FieldSet.from_parcels(filename)
+    fieldset.add_constant('halo_west', fieldset.U.grid.lon[0])
+    fieldset.add_constant('halo_east', fieldset.U.grid.lon[-1])
+    fieldset.add_constant('halo_south', fieldset.U.grid.lat[0])
+    fieldset.add_constant('halo_north', fieldset.U.grid.lat[-1])
     fieldset.add_periodic_halo(zonal=True, meridional=True)
     pset = ParticleSet.from_list(fieldset=fieldset,
                                  pclass=ptype[mode],
diff --git a/parcels/examples/tutorial_FieldLists.ipynb b/parcels/examples/tutorial_FieldLists.ipynb
new file mode 100644
index 00000000..cb06d1d6
--- /dev/null
+++ b/parcels/examples/tutorial_FieldLists.ipynb
@@ -0,0 +1,231 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Tutorial on how to combine different Fields for advection into a `FieldList` object"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In some oceanographic applications, you may want to advect particles using a combination of different velocity data sets. For example, particles at the surface are transported by a combination of geostrophic, Ekman and Stokes flow. And often, these flows are not even on the same grid.\n",
+    "\n",
+    "One option would be to write a  `Kernel` that computes the movement of particles due to each of these flows. However, in Parcels it is possible to directly combine different flows (without interpolation) and feed them into the built-in `AdvectionRK4` kernel. For that, we use so-called `FieldList` objects.\n",
+    "\n",
+    "This tutorial shows how to use these `FieldLists` with a very idealised example. We start by importing the relevant modules."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from parcels import Field, FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile, AdvectionRK4\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, let's first define a zonal and meridional velocity field on a 1kmx1km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is zero everywhere."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xdim, ydim = (10, 20)\n",
+    "U = Field('U', np.ones((ydim, xdim), dtype=np.float32),\n",
+    "          lon=np.linspace(0., 1e3, xdim, dtype=np.float32),\n",
+    "          lat=np.linspace(0., 1e3, ydim, dtype=np.float32))\n",
+    "V = Field('V', np.zeros((ydim, xdim), dtype=np.float32), grid=U.grid)\n",
+    "fieldset = FieldSet(U, V)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We then run a particle and plot its trajectory"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gn/T/parcels-501/2fccd08be708344bcd942ed95f41d66c.so\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFAJJREFUeJzt3Xu0XnV95/H3B4LcFNCQKCZgQk1HLCOVHjFKS2fEWqSM\n8dYODowUL8x00YJYx0HbNTO90KVLl61UZZqCipVSHUCHxcIKQx1aq6GeBGGCUJsFApEoUTJcBmyT\n5jt/PDvLh/DLOQ9w9nmSc96vtc569rNvv+9OsvI5+/fbl1QVkiTtbK9xFyBJ2j0ZEJKkJgNCktRk\nQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1LRh3AU/HoYceWsuWLRt3GZK0R1m7du0PqmrRdOvt\n0QGxbNkyJicnx12GJO1Rktw9ynp2MUmSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLU\nZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1NRrQCQ5N8n6JLcleddOy96TpJIc\n2n1PkguTbEhya5Jj+6xNkjS13gIiydHAO4HjgGOAU5Ks6JYdDvwCcM/QJq8FVnQ/ZwEX9VWbJGl6\nfZ5BHAWsqapHq2obcCPwhm7ZHwLvBWpo/VXAZ2pgDXBIksN6rE+SNIU+A2I9cEKShUkOAE4GDk/y\nOuC7VXXLTusvAe4d+r6xm/c4Sc5KMplkcvPmzX3VLknz3oK+dlxVtyf5IHA98AhwC7AN+C3gNY1N\n0tpNY7+rgdUAExMTT1guSZoZvQ5SV9UlVXVsVZ0APAB8B1gO3JLkO8BSYF2S5zE4Yzh8aPOlwH19\n1idJ2rW+r2Ja3H0eAbyRwRjD4qpaVlXLGITCsVX1PeBq4K3d1UwrgQeralOf9UmSdq23LqbOlUkW\nAluBs6tqyxTrXstgnGID8ChwZs+1SZKm0GtAVNXPTbN82dB0AWf3WY8kaXTeSS1JajIgJElNBoQk\nqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKa\nDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmA\nkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTb0GRJJzk6xPcluSd3XzPpTkjiS3JvlCkkOG\n1n9fkg1J/j7JL/ZZmyRpar0FRJKjgXcCxwHHAKckWQFcDxxdVS8Bvg28r1v/xcCpwE8BJwGfSLJ3\nX/VJkqbW5xnEUcCaqnq0qrYBNwJvqKrruu8Aa4Cl3fQq4C+q6h+r6i5gA4NwkSSNQZ8BsR44IcnC\nJAcAJwOH77TO24AvddNLgHuHlm3s5kmSxmBBXzuuqtuTfJBBl9IjwC3AjjMHkvxW9/2yHbNau9l5\nRpKzgLMAjjjiiBmuWpK0Q6+D1FV1SVUdW1UnAA8A/wCQ5AzgFOC0qtoRAht5/BnGUuC+xj5XV9VE\nVU0sWrSoz/IlaV7r+yqmxd3nEcAbgcuTnAT8Z+B1VfXo0OpXA6cm2TfJcmAF8Hd91idJ2rXeupg6\nVyZZCGwFzq6qLUk+BuwLXJ8EBgPZ/7GqbkvyeeBbDLqezq6qf+65PknSLvQaEFX1c415L5xi/QuA\nC/qsSZI0Gu+kliQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS\n1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVLTSAGRgdOT/Jfu\n+xFJjuu3NEnSOI16BvEJ4BXAW7rvDwMf76UiSdJuYcGI6728qo5NcjNAVW1J8owe65IkjdmoZxBb\nk+wNFECSRcD23qqSJI3dqAFxIfAFYHGSC4CvAn/QW1WSpLEbqYupqi5LshY4EQjw+qq6vdfKJElj\nNWVAJHnO0Nf7gcuHl1XVA30VJkkar+nOINYyGHcIcASwpZs+BLgHWN5rdZKksZlyDKKqllfVkcCX\ngX9TVYdW1ULgFOCq2ShQkjQeow5Sv6yqrt3xpaq+BPx8PyVJknYHo94H8YMkvw18lkGX0+nAD3ur\nSpI0dqOeQbwFWMTgUtcvAov58V3VkqQ5aNTLXB8Azu25FknSbmSkgEjyFbq7qIdV1atmvCJJ0m5h\n1DGI9wxN7we8Cdg23UZJzgXeyeDS2D+tqj/q7q34HLAM+A7wK92znQJ8FDgZeBT41apaN2J9kqQZ\nNtIYRFWtHfr526p6N/DyqbZJcjSDcDgOOAY4JckK4HzghqpaAdzQfQd4LbCi+zkLuOipHJAkaWaM\n+j6I5wz9HJrkF4HnTbPZUcCaqnq0qrYBNwJvAFYBl3brXAq8vpteBXymBtYAhyQ57Mke0CjW3r2F\nj39lA2vv3tLH7nfLtj3m2TXfjtk/69k1W22P2sU0fEf1NuAu4O3TbLMeuCDJQuAxBl1Hk8Bzq2oT\nQFVtSrK4W38JcO/Q9hu7eZtGrHEka+/ewq/8ydf55+3FXoEXPe9ZPGu/fWayiV16+EdbueN7D7O9\nmNW2x9XuONv2mOd+u+Nse3c45irYd5+9uOwdK/mZFzy7l7ZGvcz1qKo6sruzekVVvQb4xlQbdA/z\n+yBwPfCXwC1MPW6R1m6esFJyVpLJJJObN28esfwfW3PnD9m+fbDb7QUP/WjaoZQZ89CPttE1Patt\nj6vdcbbtMc/9dsfZ9u5wzAVs3badNXf2d0vaqGcQXwOO3Wne1xvzHqeqLgEuAUjyBwzOCr6f5LDu\n7OEwBg8BpFt2+NDmS4H7GvtcDawGmJiYeEKATGflkQvZd5+92LptO/ss2IuPnvrS3tJ3Z2vv3sJp\nF6+Z9bbH1e442/aY536742x7dzrmlUcu7K2tVO36/9gkz2PQzfNZ4N/x49/yDwL+e1W9aMqdJ4ur\n6v4kRwDXMXht6fuBH1bVB5KcDzynqt6b5JeAX2fQFfVy4MKqmvK91xMTEzU5OTnKcT7O2ru3sObO\nH7LyyIWz9pc67rY9Zo95LrY7zrb35GNOsraqJqZdb5qAOAP4VWCCwfjBDg8Dn66qKR/Yl+RvgIXA\nVuDdVXVDNybxeQZPh70H+OWqeqC7zPVjwEkMLnM9s6qm/N//qQaEJM1nMxIQQzt7U1VdOSOVzSAD\nQpKevFEDYroXBp1eVZ8FliV5987Lq+ojT6NGSdJubLpB6gO7z2c2lj3pAWJJ0p5jyoCoqj/pJv9X\nVf3t8LIkx/dWlSRp7Ea9D+KPR5wnSZojphuDeAXwSmDRTmMQBwF791mYJGm8phuDeAaD8YcFwLOG\n5j8EvLmvoiRJ4zfdGMSNwI1JPl1Vd89STZKk3cCoj9p4NMmHgJ9i8D4IwBcGSdJcNuog9WXAHcBy\n4HcYvOhnyof1SZL2bKMGxMLuwXtbq+rGqnobsLLHuiRJYzZqF9PW7nNT91C9+xg8bVWSNEeNGhC/\nn+Rg4DcZ3P9wEPCu3qqSJI3dSAFRVdd0kw8C/xogiQEhSXPYqGMQLU94eJ8kae54OgHRekWoJGmO\neDoB4dNcJWkOm+5ZTA/TDoIA+/dSkSRptzDdozaeNdVySdLc9XS6mCRJc5gBIUlqMiAkSU0GhCSp\nyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoM\nCElSU68BkeS8JLclWZ/k8iT7JTkxybok30zy1SQv7NbdN8nnkmxIclOSZX3WJkmaWm8BkWQJcA4w\nUVVHA3sDpwIXAadV1U8Dfw78drfJ24EtVfVC4A+BD/ZVmyRpen13MS0A9k+yADgAuA8o4KBu+cHd\nPIBVwKXd9BXAiUnSc32SpF1Y0NeOq+q7ST4M3AM8BlxXVdcleQdwbZLHgIeAld0mS4B7u223JXkQ\nWAj8oK8aJUm71mcX07MZnBUsB54PHJjkdOA84OSqWgp8CvjIjk0au6nGfs9KMplkcvPmzf0UL0nq\ntYvp1cBdVbW5qrYCVwHHA8dU1U3dOp8DXtlNbwQOB+i6pA4GHth5p1W1uqomqmpi0aJFPZYvSfNb\nnwFxD7AyyQHdWMKJwLeAg5P8ZLfOLwC3d9NXA2d0028G/qqqnnAGIUmaHX2OQdyU5ApgHbANuBlY\nzeBM4cok24EtwNu6TS4B/izJBgZnDqf2VZskaXrZk39Jn5iYqMnJyXGXIUl7lCRrq2piuvW8k1qS\n1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElN\nBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRA\nSJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVJTrwGR5LwktyVZn+Ty\nJPtl4IIk305ye5JzunWT5MIkG5LcmuTYPmuTJE1tQV87TrIEOAd4cVU9luTzwKlAgMOBF1XV9iSL\nu01eC6zofl4OXNR9SpLGoO8upgXA/kkWAAcA9wG/BvxuVW0HqKr7u3VXAZ+pgTXAIUkO67k+SdIu\n9BYQVfVd4MPAPcAm4MGqug74CeDfJplM8qUkK7pNlgD3Du1iYzdPkjQGvQVEkmczOCtYDjwfODDJ\n6cC+wI+qagL4U+CTOzZp7KYa+z2rC5fJzZs391O8JKnXLqZXA3dV1eaq2gpcBbySwZnBld06XwBe\n0k1vZDA2scNSBl1Sj1NVq6tqoqomFi1a1FvxkjTf9RkQ9wArkxyQJMCJwO3AF4FXdev8PPDtbvpq\n4K3d1UwrGXRJbeqxPknSFHq7iqmqbkpyBbAO2AbcDKwG9gcuS3Ie8Ajwjm6Ta4GTgQ3Ao8CZfdUm\nSZpeqp7Qzb/HmJiYqMnJyXGXIUl7lCRru3HgKXkntSSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKT\nASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTXv0C4OSbAbufoqb\nHwr8YAbL2RN4zPODxzw/PJ1jfkFVLZpupT06IJ6OJJOjvFFpLvGY5wePeX6YjWO2i0mS1GRASJKa\n5nNArB53AWPgMc8PHvP80Psxz9sxCEnS1ObzGYQkaQrzMiCSnJTk75NsSHL+uOvpW5LDk3wlye1J\nbkty7rhrmg1J9k5yc5Jrxl3LbElySJIrktzR/X2/Ytw19SnJed2/6fVJLk+y37hr6kOSTya5P8n6\noXnPSXJ9kn/oPp890+3Ou4BIsjfwceC1wIuBtyR58Xir6t024Der6ihgJXD2PDhmgHOB28ddxCz7\nKPCXVfUi4Bjm8PEnWQKcA0xU1dHA3sCp462qN58GTtpp3vnADVW1Arih+z6j5l1AAMcBG6rqzqr6\nJ+AvgFVjrqlXVbWpqtZ10w8z+E9jyXir6leSpcAvARePu5bZkuQg4ATgEoCq+qeq+r/jrap3C4D9\nkywADgDuG3M9vaiqvwYe2Gn2KuDSbvpS4PUz3e58DIglwL1D3zcyx/+zHJZkGfBS4KbxVtK7PwLe\nC2wfdyGz6EhgM/Cprmvt4iQHjruovlTVd4EPA/cAm4AHq+q68VY1q55bVZtg8EsgsHimG5iPAZHG\nvHlxKVeSZwJXAu+qqofGXU9fkpwC3F9Va8ddyyxbABwLXFRVLwX+Hz10O+wuuj73VcBy4PnAgUlO\nH29Vc8t8DIiNwOFD35cyR09LhyXZh0E4XFZVV427np4dD7wuyXcYdCG+Kslnx1vSrNgIbKyqHWeH\nVzAIjLnq1cBdVbW5qrYCVwGvHHNNs+n7SQ4D6D7vn+kG5mNAfANYkWR5kmcwGNS6esw19SpJGPRL\n315VHxl3PX2rqvdV1dKqWsbg7/evqmrO/2ZZVd8D7k3yL7pZJwLfGmNJfbsHWJnkgO7f+InM4UH5\nhquBM7rpM4D/OdMNLJjpHe7uqmpbkl8HvszgqodPVtVtYy6rb8cD/x74P0m+2c17f1VdO8aa1I/f\nAC7rfvm5EzhzzPX0pqpuSnIFsI7BlXo3M0fvqE5yOfCvgEOTbAT+K/AB4PNJ3s4gLH95xtv1TmpJ\nUst87GKSJI3AgJAkNRkQkqQmA0KS1GRASJKaDAjNW0ke6Xn/F+94KGKS9z+F7ZcNP71Tmm1e5qp5\nK8kjVfXM3bWt7rlZ13RPKpVmnWcQ0pAkL0hyQ5Jbu88juvmfTnJhkq8luTPJm7v5eyX5RPdOgmuS\nXDu07H8nmUjyAQZPHP1mkst2PjNI8p4k/62b/pkktyT5OnD20Dp7J/lQkm90tf2HWfxj0TxlQEiP\n9zHgM1X1EuAy4MKhZYcBPwucwuAuVoA3AsuAfwm8A3jCC3qq6nzgsar66ao6bZr2PwWcU1U77+ft\nDJ5W+jLgZcA7kyx/MgcmPVkGhPR4rwD+vJv+MwaBsMMXq2p7VX0LeG4372eB/9HN/x7wlafacJKD\ngUOq6sah9nd4DfDW7lEpNwELgRVPtS1pFPPuWUzSkzQ8SPePQ9PZ6fPJ2Mbjfznb8ZrMsOtHzwf4\njar68lNoT3pKPIOQHu9r/Pi1lacBX51m/a8Cb+rGIp7L4IFqLVu7R64DfB9YnGRhkn0ZdFnRvf3t\nwSQ7zlqGu6O+DPzajn0k+cm5/DIg7R48g9B8dkD3ZMwdPsLgHcefTPKfGLydbbqnoV7J4DHT64Fv\nM+j+ebCx3mrg1iTrquq0JL/brXsXcMfQemd27T/KIBR2uJjBWMe67tHWm+nhFZPSMC9zlZ6mJM+s\nqkeSLAT+Dji+G4+Q9mieQUhP3zVJDgGeAfye4aC5wjMISVKTg9SSpCYDQpLUZEBIkpoMCElSkwEh\nSWoyICRJTf8f8JPPHrVCUx0AAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x1237671d0>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<module 'matplotlib.pyplot' from '/Users/erik/miniconda2/envs/py2_parcels/lib/python2.7/site-packages/matplotlib/pyplot.pyc'>"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[900])\n",
+    "output_file = pset.ParticleFile(name='FieldListParticle_adv.nc', outputdt=1)\n",
+    "pset.execute(AdvectionRK4, runtime=10, dt=1, output_file=output_file)\n",
+    "plotTrajectoriesFile('FieldListParticle_adv.nc')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The trajectory plot shows a particle moving eastward on the 1 m/s flow, as expected"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, let's define another set of velocities (`Ustokes, Vstokes`) on a different, higher-resolution grid. This flow is southwestward at -0.2 m/s in each direction.\n",
+    "\n",
+    "Note that it is **very important to specify the `fieldtype` of the fields, as Parcels will otherwise not perform the unit conversion correctly**."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gf = 10  # factor by which the resolution of this grid is higher than of the original one.\n",
+    "Ustokes = Field('Ustokes', -0.2*np.ones((ydim*gf, xdim*gf), dtype=np.float32),\n",
+    "                lon=np.linspace(0., 1e3, xdim*gf, dtype=np.float32),\n",
+    "                lat=np.linspace(0., 1e3, ydim*gf, dtype=np.float32),\n",
+    "                fieldtype='U')\n",
+    "Vstokes = Field('Vstokes', -0.2*np.ones((ydim*gf, xdim*gf), dtype=np.float32), grid=Ustokes.grid, \n",
+    "                fieldtype='V')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now comes the trick of the `FieldLists`. We can simply define a new `FieldSet` with a list of different `Fields`, as in `U=[U, Ustokes]`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fieldset = FieldSet(U=[U, Ustokes], V=[V, Vstokes])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And if we then run the particle again and plot its trajectory, we see that it moves slightly southward too (and less far eastward), because of the Stokes drift."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gn/T/parcels-501/2c8d709ac096936d044d226a1fbf84ad.so\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VWXW9/HvL4UqJUBAelGkSFMCRhBwBMWOAirKWFDE\ncVQQnnFQn3eaOjOM4wOKOChWdJBRAR1wQEB0CCDF0HuXroYivWe9f5wdjUgJmMM+Sdbnus6VfXa5\n99pemJW77PuWmeGcc87ltriwA3DOOZc/eYJxzjkXFZ5gnHPORYUnGOecc1HhCcY551xUeIJxzjkX\nFZ5gnHPORYUnGOecc1HhCcY551xUJIQdQJjKlStnNWrUCDsM55zLU2bPnr3VzJJPdV6BTjA1atQg\nPT097DCccy5PkbQuJ+d5E5lzzrmo8ATjnHMuKjzBOOeciwpPMM4556LCE4xzzrmoiGqCkdRL0iJJ\niyU9GuwrI2mipJXBz6RgvyQNlLRK0gJJF5+gzKaSFgbnDZSkk5XrnHMuHFFLMJIaAPcDzYHGwPWS\nagOPA5PMrDYwKfgOcA1QO/j0AAafoOjBwfGsc68O9p+o3Fw3e90OXvp8FbPX7YjWLZxzLs+LZg2m\nHjDDzPaZ2RFgMnAz0AEYGpwzFLgp2O4AvG0RM4DSkipmLzD4XtLMpltkree3j7n+eOXmqtnrdnDH\nqzN4bvxyur46w5OMc86dQDQTzCKgtaSykooB1wJVgQpmtgUg+Fk+OL8ysCHb9RuDfdlVDvYf75wT\nlfsjknpISpeUnpGRcdoPNWPNNg4dycSAA0cyGTYjR+8bOedcgRO1BGNmS4G/AROBT4D5wJGTXKLj\nFXMG55wqriFmlmJmKcnJp5zp4CdSa5WlcGIccYoEM2ruJnoOn8u2PQdPuyznnMvPojpVjJm9DrwO\nIOkvRGoc30iqaGZbgiavb4PTNxKp4WSpAmw+psiNwf7jnXOicnNV0+pJDOueyow120ipnsSMNdsZ\n9PlKpqzM4I83XsiNjSsRjDtwzrkCLdqjyMoHP6sBHYHhwGjg7uCUu4F/B9ujgbuC0WSpwM6sJq8s\nwffdklKD0WN3HXP98crNdU2rJ/HQL87nklpl6dWuNv/p2YrqZYvT61/z6D40nS0790fr1s45l2co\n0lcepcKlKUBZ4DDQx8wmSSoLvA9UA9YDt5jZ9iBhDCIyKmwf0M3M0oNy5plZk2A7BXgLKAqMAx4x\nMztRuSeLLyUlxXJrssujmcab09by3ITlJMTF8cS1dbm9WTXi4rw245zLXyTNNrOUU54XzQQT63Iz\nwWRZv20fj49awBert3FJzTL069SImuWK5+o9nHMuTDlNMP4mfy6rVrYYw7pfwt86NWTJll1c/Xwa\nQ9JWc+RoZtihOefcWeUJJgokcVuzanzapw2tL0jmL2OX0XHwFyzdsivs0Jxz7qzxBBNFFUoWYcid\nTRl0x0Vs2rGfG16cSv+JKzh45GjYoTnnXNR5gokySVzfqBKf9mnDDY0rMXDSSq4fOJU5630GAOdc\n/uYJ5ixJKl6IAbc14c17mrHn4BE6Df6Cpz9ewr5DJ3v31Dnn8i5PMGfZL+qWZ0Lv1nS9pBqvT11L\n++fTmLZqa9hhOedcrvMEE4ISRRJ55qaGvNcjlYS4OLq+NpO+Ixawc//hsENzzrlc4wkmRJfUKsu4\nXq34VZvzGDFnI1f2n8yExV+HHZZzzuUKTzAhK5IYz+PX1OWjX7ek7DmF6fHObB56dw4Zu33yTOdc\n3uYJJkY0rFKK0Q+35DdXXcDExd9w5YDJfDh3IwV5pgXnXN7mCSaGJMbH8fAVtRnb6zJqlStO7/fm\n0+2tL9n0nU+e6ZzLezzBxKDzy5fgg1+14A831Gfmmu1c1X8y78xYR2am12acc3mHT3aZy5Nd5rYN\n2/fx5IcLmbJyK81rlOHOS6uxfvt+UmuVpWn1pLDDc84VQD6bcg7khQQDYGaMmL2RP4xezL5DRxFQ\nODGOYd1TPck45846n005H5HELSlVuadFDSCyRvSBw5mMnnfsgp/OORc7or2iZW9JiyUtkjRcUhFJ\nV0iaE+wbKikhODdJ0oeSFkiaJanBCcqcImle8Nks6aNg/+WSdmY79vtoPlsY2tarQJHEOLLWMPvn\njK94bvxyDhz2yTOdc7EnIVoFS6oM9ATqm9l+Se8DdwB/Atqa2QpJTxFZ3vh14ElgnpndLKku8BLQ\n9thyzaxVtnuM5MdLI08xs+uj9Uxha1o9iWHdU5mxZhsNKpVk9PwtDPp8FeMWbeHZzo1oWr1M2CE6\n59z3ot1ElgAUDWopxYC9wEEzWxEcnwh0CrbrA5MAzGwZUENShRMVLKkEcAXwUZRij0lNqyfx0C/O\np02d8vzfrY0Zem9zDhzOpPPL0/nj6MXsPeiTZzrnYkPUEoyZbQKeA9YDW4CdwPtAoqSszqHOQNVg\nez7QEUBSc6A6UOUkt7gZmGRm2VfxulTSfEnjJF2Yaw8Tw9pckMz43q25K7U6Q6d/xVUD0khbkRF2\nWM45F70EIykJ6ADUBCoBxYGuQBdggKRZwG4g60/ufkCSpHnAI8DcbMeO53ZgeLbvc4DqZtYYeJET\n1Gwk9ZCULik9IyN//CI+p3ACf+rQgPcfuJTCiXHc9cYsfvPBfHbu88kznXPhidowZUm3AFeb2X3B\n97uAVDP7dbZzrgK6m9mtx1wrYC3Q6JgaStbxssAKoLKZHTjB/b8CUszshHPh55VhyqfjwOGjDJy0\nklfS1lCmeCGe7nAhVzeoGHZYzrl8JBaGKa8HUiUVCxJGW2CppPJBgIWBvsDLwffSkgoF13YH0o6X\nXAK3AB9nTy6Szg3uk9XEFgdsi8JzxbQiifH89uq6/PuhlpQvUZhf/XMOD/5zNt/uPm4eds65qIlm\nH8xMYASRpquFwb2GAI9JWgosAMaY2WfBJfWAxZKWAdcAvbLKkjRWUqVsxXfhx81jEOnPWSRpPjAQ\n6GIF+C3SBpVL8dFDLfnt1XWYtOxbruyfxgfpG3zyTOfcWeNv8uezJrLjWZ2xh8dHLuDLr3bQqnY5\n/nJzQ6qWKRZ2WM65PCoWmshcjDgv+Rze63EpT3W4kDnrdtD++TTemrbWJ890zkWVJ5gCIi5O3HVp\nDcb3bk2zGmX445gl3PrKdFZ9uyfs0Jxz+ZQnmAKmSlIx3urWjP63NmZVxh6ufWEKL32+isNHM8MO\nzTmXz3iCKYAk0fHiKkzs3YYr61fg7+OX02HQNBZt2hl2aM65fMQTTAGWXKIwL3W9mJd/2ZSMPQfp\n8NI0/vbJMp880zmXKzzBOK5ucC6f9m5Dp4srM/i/q7n2hSl8+dX2sMNyzuVxnmAcAKWKJfJs58b8\n875LOHQ0k1tens7v/72IPT55pnPuDHmCcT9yWe1yTOjdmntb1uSdGeu4qv9kPl/+bdhhOefyIH/R\nsgC8aHmmZq/bQd+RC1j17R7aXFCOhpVL84u65X2ZZucKOH/R0v1sTasn8Z+el3FL08pMXrGVQZ+v\nosuQ6cz2/hnnXA54gnEnVTghnhrlzvl+mebDR42+oxby7S6fPNM5d3KeYNwppdYqS6GEOOIFCXFi\n3ba9tO0/mfe/9MkznXMnlhB2AC72Na2exLDuqcxYs43UWmUpU7wQfUcu4LcjFzB6/mb+cnNDqpX1\nyTOdcz/mnfzeyX9GMjONd2etp9+4ZRzNNH7Tvg73tKhBfFZbmnMu3/JOfhdVcXHil6nVmdC7Nam1\nyvD0x0vo/PIXrPxmd9ihOedihCcY97NUKl2UN+5pxvO3NeGrrXu5buBUBk5ayaEjPnmmcwVdVBOM\npN6SFktaJGm4pCKSrpA0J9g3VFJCcG6SpA8lLZA0S1KDE5T5lqS1kuYFnybBfkkaKGlVUMbF0Xw2\n9wNJ3HRRZSb2aUP7BufSf+IKbhw0lQUbvws7NOdciKKWYCRVBnoCKWbWAIgH7gCGElnOuAGwDrg7\nuORJYJ6ZNQLuAl44SfGPmVmT4DMv2HcNUDv49AAG5/YzuZMrd05hXrz9Il69K4Ud+w5x00vT+OvY\npT55pnMFVLSbyBKAokEtpRiwFzhoZiuC4xOBTsF2fWASgJktA2pIqnAa9+oAvG0RM4DSkirmxkO4\n03Nl/QpM7NOG25pV5ZW0NVz9fBoz1mwLOyzn3FkWtQRjZpuA54D1wBZgJ/A+kCgpa/RBZ6BqsD0f\n6AggqTlQHahyguL/HDSDDZBUONhXGdiQ7ZyNwT4XgpJFEvlrx0a82/0SMg26DJnB/364kN0HDocd\nmnPuLIlmE1kSkVpFTaASUBzoCnQBBkiaBewGsqbr7QckSZoHPALMzXYsuyeAukAzoAzQN+uWxzn3\nJ2OwJfWQlC4pPSMj4wyfzuVUi/PLMf7R1nS/rCbDZ63nqgFpfLbsm7DDcs6dBdFsImsHrDWzDDM7\nDIwCWpjZdDNrZWbNgTRgJYCZ7TKzbmbWhEgfTDKw9thCzWxL0Ax2EHgTaB4c2sgPtSGI1H42H+f6\nIWaWYmYpycnJufe07oSKForn/11fn5EPtqBEkQTufSudR/81l+17D4UdmnMuiqKZYNYDqZKKSRLQ\nFlgqqTxA0LTVF3g5+F5aUqHg2u5AmpntOrbQrH6VoMybgEXBodHAXcFoslRgp5ltid7judN1UbUk\nPn6kFb3a1uY/C7fQrv9kRs/f7NPNOJdPRbMPZiYwApgDLAzuNQR4TNJSYAEwxsw+Cy6pByyWtIzI\niLBeWWVJGiupUvB1mKSFQZnlgGeC/WOBNcAq4FXg19F6NnfmCiXE0fvKCxjzyGVUTSpKz+Fzuf/t\n2Xy90yfPdC6/8alifKqY0BzNNN6ctpbnJiwnMS6OJ6+rR5dmVYlUTp1zscqninExLz5OdG9Vi/GP\ntqZB5VI8MWohd7w6k3Xb9oYdmnMuF3iCcaGrXrY4795/CX/t2JBFm3bS/vk0XpuyhqOZBbd27Vx+\n4AnGxQRJ3N68GhP7tOGy88vxzH+W0nHwFyz/2ifPdC6v8j4Y74OJOWbGmAVb+OPoxew+cJibL6pM\nlaRitDy/HE2rJ4UdnnMFXk77YDzBeIKJWdv3HqLXv+YyZeVWIDICbfj9qZ5knAuZd/K7PK9M8UKk\n1ir7/RQNh45k0m/cUvYf8skzncsLPMG4mJZaqyyFE+OIV2TU2Zdf7aD982l8sXpr2KE5507Bm8i8\niSzmzV63gxlrtpFaqyyHj2by+MgFfLVtH7c3r8YT19alZJHEsEN0rkDxPpgc8ASTN+0/dJTnP13B\nq1PWkFyiMH++qSHt6p/Oyg7OuZ/D+2BcvlW0UDxPXFuPjx5qSVKxQnR/O51Hhs9l256DYYfmnMvG\nE4zLsxpVKc3ohy+jz5UX8MmiyOSZ/563ySfPdC5GeIJxeVqhhDh6tq3Nf3q2onrZ4vT61zzuG5rO\n5u/2hx2acwWeJxiXL1xQoQQjH2zB766vz/TV27hqQBrDZq4j06ebcS40nmBcvhEfJ+67rCbjH21N\n46ql+N8PF3H7qzNYu9Unz3QuDJ5gXL5TrWwx/nnfJTzbqRFLtuzi6ufTeGXyao4czQw7NOcKFE8w\nLl+SxK3NqvJpnza0viCZv45bRsfBX7B0y08WSXXORUlUE4yk3pIWS1okabikIpKukDQn2DdUUkJw\nbpKkDyUtkDRLUoMTlDlM0vLg+jckJQb7L5e0U9K84PP7aD6byxsqlCzCkDub8tIdF7P5u/3c8OJU\n+k9YzsEjPt2Mc9EWtQQjqTLQE0gxswZAPHAHMBToEuxbB9wdXPIkMM/MGgF3AS+coOhhQF2gIVAU\n6J7t2BQzaxJ8nsrtZ3J5kySua1SRib3bcGPjSgz8bBXXD5zKnPU7wg7NuXwt2k1kCUDRoJZSDNgL\nHDSzFcHxiUCnYLs+MAnAzJYBNST95PVsMxtrAWAWUCXKz+DyiaTiheh/WxPe7NaMvQeP0GnwFzw1\nZgn7Dh0JOzTn8qWoJRgz2wQ8B6wHtgA7gfeBRElZUwx0BqoG2/OBjgCSmgPVOUnyCJrG7gQ+ybb7\nUknzJY2TdOEJrushKV1SekZGxhk/n8u7flGnPBP6tOHO1Oq8MW0t7Z9PY+pKnzzTudwWzSayJKAD\nUBOoBBQHugJdgAGSZgG7gaw/H/sBSZLmAY8Ac7MdO55/AGlmNiX4PgeobmaNgReBj453kZkNMbMU\nM0tJTk7+OY/o8rBzCifwVIcGvP/ApSTExfHL12fy2xHz2bn/cNihOZdvRLOJrB2w1swyzOwwMApo\nYWbTzayVmTUH0oCVAGa2y8y6mVkTIn0wycDa4xUs6Q/B8T5Z+4Lr9wTbY4nUlMpF8flcPtC8ZhnG\n9WrFg5efx8g5m7iy/2TGL/467LCcyxeimWDWA6mSikkS0BZYKqk8gKTCQF/g5eB7aUmFgmu7E6md\n/GRMqaTuQHvgdjPLzLb/3OA+WU1sccC2qD2dyzeKJMbT9+q6fPTrlpQ9pzAPvDObh4bNYdLSb3jp\n81XMXueDAZw7EwnRKtjMZkoaQaTp6giRJq8hwDOSrieSAAab2WfBJfWAtyUdBZYA92WVJWks0N3M\nNhNJSOuA6UE+GRWMGOsMPCjpCLCfyEg1nyfE5VjDKqUY/XBLhqStYcDEFfxn4RYEFE6MY1h3X6rZ\nudPl68H4ejDuOJ76eDFvTP3q++8PtKnFE9fUCy8g52KIrwfj3M9wXcNKFEmMQ8H3odO+4p3pX/nk\nmc6dhhzVYIK+ja5ALTN7SlI14FwzmxXtAKPJazDuZLKWaq5VrjjvzlrPlJVbaVYjiX6dGnFe8jlh\nh+dcaHJ1yWRJg4FM4AozqxcMQZ5gZs1+fqjh8QTjcsrMGDF7I09/vIQDRzJ5tF1terSqRUK8NwK4\ngie3m8guMbOHgAMAZrYDKHTyS5zLPyRxS0pVPv2fNlxRpzzPfrKcm/4xjcWbd4YdmnMxK6cJ5rCk\neMAAJCUTqdE4V6CUL1GEl+9syuCuF/P1zoPcOGgafx+/jAOHffJM546V0wQzEPgQKC/pz8BU4C9R\ni8q5GHdNw4p82qc1NzWpzEufr+bagVNI/2p72GE5F1NyPExZUl0iL0sKmGRmS6MZ2NngfTAuN0xe\nkcGToxayeed+7r60Bo+1r0PxwlF7xcy50OVKJ7+kMie72Mzy9J9snmBcbtl78Ah/H7+codO/olKp\novy1Y0NaX+Bz3bn8Kbc6+WcD6cHPDGAFkbnDMoJ9zjmgeOEE/njjhXzwwKUUTozjrjdm8ZsP5vPd\nvkNhh+ZcaE6aYMysppnVAsYDN5hZOTMrC1xPZPJK51w2KTXKMLZnKx76xXl8OHcT7fqnMW7hlrDD\nci4UOe3kbxbMUAyAmY0D2kQnJOfytiKJ8TzWvi6jH25JhZKFeXDYHH71zmy+3XUg7NCcO6tymmC2\nSvp/kmpIqi7pf/GZip07qQsrleKjh1ry26vr8Nnyb2nXfzIfpG+gIM//5wqWnCaY24msv/IhkYW8\nygf7nHMnkRgfx68vP59xvVpR59wSPDZiAXe9MYsN2/eFHZpzUeezKfsoMneWZGYaw2auo9+4ZRjw\n2/Z1uPPSGsTH6ZTXOhdLcnsuss8J3uLPzsyuOLPwYoMnGBeGjTv28b8fLmLyigyaVk/ib50acn75\nEmGH5VyO5fZcZL8BHgs+vwPmERm+fKogektaLGmRpOGSiki6QtKcYN9QSQnBuUmSPpS0QNIsSQ1O\nUGZNSTMlrZT0XtYqmJIKB99XBcdr5PDZnDurqiQV461uzeh/a2NWZ+zh2hemMuizlRw+6rMvufwl\nRwnGzGZn+0wzsz7AJSe7RlJloCeQYmYNgHjgDmAokdUmGxBZmfLu4JIngXlm1gi4C3jhBEX/DRhg\nZrWBHfyw8uV9wA4zOx8YEJznXEySRMeLqzCxdxuurF+B5yas4MZB01i0ySfPdPlHjhKMpDLZPuUk\ntQfOzcGlCUDRoJZSDNgLHDSzFcHxiUCnYLs+MAnAzJYBNSRVOCYOAVcAI4JdQ4Gbgu0OwXeC422D\n852LWcklCvNS14t5+ZdN2brnIB1emka/ccuYvnorL32+itnrdoQdonNnLKcTJs0m0gcj4Aiwlh9q\nDsdlZpskPQesB/YDE4D3gWclpZhZOtAZqBpcMh/oCEyV1ByoDlQBvslWbFngOzM7EnzfCFQOtisD\nG4J7H5G0Mzh/aw6f0bnQXN3gXC6tVZY/j13Cy5NX88rk1UhQKCGOYd1TaVo9KewQnTttOe2DqWdm\ntYI3+2ub2VXAlye7IFiUrANQE6gEFCeyKmYXYICkWcBuIgkLoB+QJGke8AgwN9ux74s9zq0sB8ey\nx9VDUrqk9IyMjJM9gnNnValiiTzbuTG3plTBgEyDg4czSVvxbdihOXdGcppgvjjOvumnuKYdsNbM\nMszsMJGpZVqY2XQza2VmzYE0InObYWa7zKybmTUh0geTTKSmlN1WoHTWwAAiNZzNwfZGgtpQcLwU\n8JPJOM1siJmlmFlKcrJPRuhiz23NqlEkIQ4R+Qtp2Iz1fL7ck4zLe06aYCSdK6kpkX6UiyRdHHwu\nJ9KncjLrgVRJxYK+kLbAUknlg7ILA32Bl4PvpbNGhAHdgTQz25W9QIuMqf6cSNMaRAYI/DvYHs0P\nAwY6A59ZQX7Jx+VZTasnMez+VH7Tvg5/ubkhpYsXotubX9LnvXns2OuTZ7q841TT9d8N3AOk8ONh\nybuBt8zspBNeSvoTcBuRpq65RBLHM0Qmy4wDBpvZ88G5lwJvA0eBJcB9wdLMSBoLdDezzZJqAf8C\nygRl/tLMDkoqArwDXESk5tLFzNacLD5/D8blBQePHOWlz1bxj/+uplTRRP7U4UKua1gRH8PiwpLb\nL1p2MrORuRJZDPEE4/KSJZt30XfkAhZu2slV9Svw9E0NqFCySNhhuQIotxYc+6WZ/VPS/3D8N/n7\n/7www+UJxuU1R45m8vrUtfSfuIJCCXH8v+vqcWtKVa/NuLMqt97kLx78PAcoccznnJ8VoXPutCXE\nx/FAm/P45NHW1KtYkr4jF/LL12eyfptPnuliT06byFqa2bRT7ctrvAbj8rLMTOPdWevpN24ZRzON\n37Svwz0tfPJMF325PRfZiznc55w7S+LixC9TqzOhd2tSa5Xh6Y+X0GnwF6z4ZnfYoTkHnOJN/mBk\nVwsgWVKfbIdKEplbzDkXskqli/LGPc3497zN/GnMYq4bOIVHrqjNr9qcR6GEnP4N6VzuO9W/vkJE\n+loS+HH/yy5+eBfFORcySdx0UWUm9mnD1Q0q0n/iCm4cNJX5G74LOzRXgOW0D6a6ma07C/GcVd4H\n4/KriUu+4f99tJCM3Qe5v1UtHm13AUULeaODyx057YPJ6WSX+yT9HbgQ+H7gfV5fcMy5/OrK+hW4\npFYZ/jp2Ka+krWH84q/p16kRqbXKhh2aK0By2kA7DFhGZOLKPwFfcYrJLp1z4SpZJJG/dmzEu90v\nIdOgy5AZPPnhQnYdOBx2aK6AyGmCKWtmrwOHzWyymd0LpEYxLudcLmlxfjnGP9qa+1vV5F+z1nNV\n/zQ+W/bNqS907mfKaYLJ+pNni6TrJF1EZCZj51weULRQPP97XX1G/bolpYomcu9b6fT611y27TkY\ndmguH8tpgnlGUingf4DfAK8Bj0YtKudcVDSpWpoxj1zGo+1qM3bhFq4ckMbo+ZvxicddNOQowZjZ\nx2a208wWmdkvzKwpcF6UY3PORUGhhDgebXcBHz/SiqplitFz+Fzufzudr3ceCDs0l8/kaJjycS+U\n1ptZtVyO56zyYcquoDuaabw5bS3PTVhOYlwcXVOrUaJIAqm1yvkyze6EcnuY8nHv8TOudc7FgPg4\n0b1VLa6sX4GHhs3h5cmRJZQKJ6zi3ftTPcm4n+XnzCPhjbbO5RPVyxbnmobnfv9X48Ejmbw4aSVH\nM/1/c3fmTrVk8m5Ju47z2Q1UOlXhknpLWixpkaThkopIukLSnGDfUEkJwbmlJI2RND+4pttxyish\naV62z1ZJWSti3iMpI9ux7mf438S5Aim1VjkKJ8YRL4gT/HdFBh3/MY3lX/vkme7MnHEfzCkLlioD\nU4H6ZrZf0vvAJ0Re1GxrZiskPQWsM7PXJT0JlDKzvpKSgeXAuWZ2wkXIJc0GeptZmqR7gBQzezin\nMXofjHM/NnvdDmas2UZqzTJs3nmAP45ezK4Dh/n15efz0C/O98kzHXB2+mByIgEoKukwUAzYCxw0\nsxXB8YnAE8DrRJrcSiiyNN85wHbgyIkKllQbKA9MiV74zhUsTasnfd/v0hRoeX45nhqzmBcmrWTc\noi0827kxTaqWDjdIl2dE7c8RM9sEPAesB7YAO4H3gURJWZmvM1A12B4E1AM2AwuBXmaWeZJb3A68\nZz+ugnWStEDSCElVT3Shcy5nyhQvxPNdLuKNe1LYfeAIHf8xjWc+XsK+Qyf828+570UtwUhKAjoQ\nmb+sEpHll7sCXYABkmYBu/mhltIemBec2wQYJKnkSW7RBRie7fsYoIaZNQI+BYaeIK4ektIlpWdk\nZJzp4zlXoFxRtwITerfm9ubVeG3qWq5+fgpfrNoadlguxkWzQbUdsNbMMszsMDAKaGFm082slZk1\nB9KAlcH53YBRFrEKWAvUPV7BkhoDCWY2O2ufmW0zs6x5L14lUsP/CTMbYmYpZpaSnJycG8/pXIFQ\nokgif765If/qkUqc4I7XZvL4yAXs3O+TZ7rji2aCWQ+kSioW9Ku0BZZKKg8gqTDQF3g52/ltg2MV\ngDrAmhOUfTs/rr0gqWK2rzcCS3PpOZxz2aTWKssnj7bmgTa1eD99A1cNmMzEJT55pvupaPbBzARG\nAHOI9KnEAUOAxyQtBRYAY8zss+CSp4EWkhYCk4C+ZrYVQNK8Y4q/lWMSDNAzGN48H+gJ3JP7T+Wc\nAyiSGM8T19Tjo4daklSsEPe/nc7D785hq0+e6bKJ2jDlvMCHKTv38x06kskrk1fz4merKF44nj/c\ncCEdmlQi0nDh8qOcDlP2Qe3OuZ+lUEIcj7StzX96XkaNcsV59L153PvWl2z+bn/YobmQeYJxzuWK\n2hVKMOKZ95RhAAAVvElEQVRXLfj99fWZsWY7Vw1I450Z68j06WYKLE8wzrlcEx8n7r2sJhN6t6ZJ\n1dL87qNFdHl1Bmu37g07NBcCTzDOuVxXtUwx3rmvOc92asTSLbu4+vk0Xp68miNHT/butMtvPME4\n56JCErc2q8qnfdrQ5oJk+o1bxk3/mMaSzbvCDs2dJZ5gnHNRVaFkEV65sykv3XExX+88wI2DpvJ/\nE5Zz8MjRsENzUeYJxjkXdZK4rlFFJvZuw41NKvHiZ6u4buBUZq/bEXZoLoo8wTjnzpqk4oXof2sT\n3urWjP2HjtL55S/405jF7D3ok2fmR/6ipb9o6Vwo9hw8wrOfLOPt6etILlGIy+uUp0uzar5Mcx7g\nL1o652LaOYUTeKpDA57p0ICtuw/xQfpGbn15OmnLfZbz/MITjHMuVDsPHCZrVpmjZjw4bDafLPo6\n3KBcrvAE45wLVWqtshRKiCNekWlnkksU5lf/nM2vh83m290Hwg7P/QzRXjLZOedOqmn1JIZ1T2XG\nmm2k1ipLoyqlGJK2hhc+Xcm0Vdv4/fX16XhxZZ88Mw/yTn7v5HcuJq36dg99Ry5g9rodtL4gmb/c\n3IAqScXCDsvhnfzOuTzu/PLn8MEDl/KnGy8k/avttB+QxtvTv/LJM/MQTzDOuZgVFyfublGD8Y+2\n5uLqSfz+34u5bch0VmfsCTs0lwNRTTCSegerTC6SNFxSEUlXSJoT7BsqKSE4t5SkMZLmB9d0O0GZ\n/5W0XNK84PP9EsyS3pO0StJMSTWi+WzOubOnaplivH1vc567pTErvtnDNS9M4R//XcVhnzwzpkUt\nwUiqTGTp4hQzawDEA3cAQ4Euwb51wN3BJQ8BS8ysMXA58H+SCp2g+K5m1iT4fBvsuw/YYWbnAwOA\nv0XjuZxz4ZBE56ZVmNinNW3rlufZT5Zz00vTWLRpZ9ihuROIdhNZAlA0qKUUA/YCB81sRXB8ItAp\n2DaghCJDRc4BtgOnM39EByLJC2AE0FY+7MS5fKd8iSIM/mVTBne9mG92HaTDS9P4+/hlHDjsk2fG\nmqglGDPbBDwHrAe2ADuB94FESVmjDzoDVYPtQUA9YDOwEOhlZieq/74ZNI/9LlsSqQxsCO59JLhf\n2WMvlNRDUrqk9IwMf2PYubzqmoYV+bRPa26+qDIvfb6aawdOIf2r7WGH5bKJZhNZEpFaRU2gElAc\n6Ap0AQZImgXs5odaSntgXnBuE2CQpJLHKbqrmTUEWgWfO7NueZxzfzLcxMyGmFmKmaUkJyef6eM5\n52JA6WKFeO6Wxrx9b3MOHs7kllem84d/L2KPT54ZE6LZRNYOWGtmGWZ2GBgFtDCz6WbWysyaA2nA\nyuD8bsAoi1gFrAXqHltoUDPCzHYD7wLNg0MbCWpDQZNcKSLNbM65fK71BclM6N2auy+twdsz1tF+\nQBqTV3gLRdiimWDWA6mSigXNWG2BpdlHfQF9gZeznd82OFYBqAOsyV6gpARJ5YLtROB6YFFweDQ/\nDBjoDHxmBfktUucKmOKFE/jjjRfywQOXUiQxjrvfmMX/vD+f7/YdCju0AiuafTAziXS2zyHSpxIH\nDAEek7QUWACMMbPPgkueBlpIWghMAvqa2VYASfOCcwoD4yUtINKctgl4NTj2OlBW0iqgD/B4tJ7N\nORe7UmqU4T89W/HwL87no3mbaNc/jXELt4QdVoHkU8X4VDHO5VuLN+/ktyMWsHjzLq6+8Fye6nAh\n5UsWCTusPM+ninHOFXgXVirFvx9qSd+r6/LZ8m9p138yH6RvoCD/YX02eYJxzuVrCfFxPHj5eYzr\n1Yo655bgsRELuOuNWWzYvi/s0PI9byLzJjLnCozMTGPYzHX0G7cMA37bvg4NKpdi5trtpNYq68s1\n51BOm8h8PRjnXIERFyfuvLQGV9SrwJOjFvLHMUuQIi/RFUqIY1j3VE8yucibyJxzBU7l0kV5q1sz\nrm14LmaQaXDwcCZfrN4admj5iicY51yBJIn7LqtF4YTIr0EDRsze6JNn5iLvg/E+GOcKtNnrdjBj\nzTYkeHPaV2zfe4j7W9Xi0Xa1KZIYH3Z4MSmnfTCeYDzBOOcCO/cd5i9jl/Je+gZqlStOv06NaF6z\nTNhhxRx/D8Y5505TqWKJ/K1zI4Z1v4TDmZnc+sp0fvfRInYfOBx2aHmSJxjnnDtGy/PLMf7R1tzb\nsib/nBmZPPPz5d+e+kL3I55gnHPuOIoVSuD3N9Rn5IMtKF44gW5vfkmf9+axY69PnplTnmCcc+4k\nLq6WxMc9L6PnFeczev5m2vWfzMcLNvt0MzngCcY5506hcEI8fa6qw5hHLqNyUlEefncuD7wzm292\nHQg7tJjmCcY553KoXsWSjHqwBU9eW5fJKzJo138y73253mszJ+AJxjnnTkNCfBw9Wp/H+EdbU79i\nSfqOXEjX12ayfptPnnksTzDOOXcGapQrzvD7U/nzzQ1YsHEn7Z9P4/Wpazma6bWZLFFNMJJ6S1os\naZGk4ZKKSLpC0pxg31BJCcG5pSSNkTQ/uKbbccorJuk/kpYF5/TLduweSRmS5gWf7tF8Nueci4sT\nXS+pzsQ+rbn0vLI8/fESOg3+ghXf7A47tJgQtQQjqTLQE0gxswZAPHAHMBToEuxbB9wdXPIQsMTM\nGgOXA/8nqdBxin7OzOoCFwEtJV2T7dh7ZtYk+LwWlQdzzrljVCxVlNfvTuGFLk1Yv30f1w2cwguf\nruTQkcywQwtVtJvIEoCiQS2lGLAXOGhmK4LjE4FOwbYBJSQJOAfYDhzJXpiZ7TOzz4PtQ8AcoEqU\nn8E5505JEh2aVGZi79Zc06AiAz5dwY2DpjJ/w3dhhxaaqCUYM9sEPAesB7YAO4H3gURJWXPYdAaq\nBtuDgHrAZmAh0MvMTpj+JZUGbgAmZdvdSdICSSMkVT3BdT0kpUtKz8jIOPMHdM654yh7TmEG3n4R\nr92Vwnf7DnPzP6bxl7FL2X/oaNihnXXRbCJLAjoANYFKQHGgK9AFGCBpFrCbH2op7YF5wblNgEGS\nSp6g7ARgODDQzNYEu8cANcysEfApkaa4nzCzIWaWYmYpycnJP/9BnXPuONrVr8CEPq25rVk1hqSt\n4ZoX0pi+elvYYZ1V0WwiawesNbMMMzsMjAJamNl0M2tlZs2BNGBlcH43YJRFrALWAnVPUPYQYKWZ\nPZ+1w8y2mdnB4OurQNMoPJNzzuVYySKJ/LVjQ969/xIMuP3VGTz54UJ2FZDJM6O5ZPJ6IFVSMWA/\n0BZIl1TezL6VVBjoC/w52/ltgSmSKgB1gDXHFirpGaAU0P2Y/RXNbEvw9UZgaRSeyTnnTluL88rx\nSa/WDPh0Ba9NWcNnS7+lW8saHMk0UmuVzbfLNEd1PRhJfwJuI9IMNpdIUngGuJ5I7WlwVi1EUiXg\nLaAikSWy+5nZP4Nj88ysiaQqwAZgGZBVWxlkZq9J+iuRxHKEyACBB81s2cni8/VgnHNn27wN39Hz\n3bms3xF5MbNwQhzv3p+ap5KMLziWA55gnHNhGDhpJQMmriDrt+/1jSry4u0XERlEG/t8wTHnnItR\nLc8vR+HEOOIEEny8YAvdh6azZef+sEPLVV6D8RqMcy4Es9ftYMaabTSvWYb5G77juQnLSYyL44lr\n69GlWVXi4mK3NuNNZDngCcY5FyvWb9vH46MW8MXqbaTWKkO/jo2oUa542GEdlzeROedcHlKtbDGG\ndb+Efh0bsnjTLto/n8aQtNUcOZp3p5vxBOOcczFCEl2aV2Ninza0qp3MX8Yuo9PgL1j29a6wQzsj\nnmCccy7GnFuqCK/e1ZQXb7+IjTv2c/3AqfSfuIKDR/LWdDOeYJxzLgZJ4obGlZjYpw03NK7EwEkr\nueHFqcxdvyPs0HLME4xzzsWwMsULMeC2Jrx5TzN2HzhCx8Ff8PTHS9h36MipLw6ZJxjnnMsDflG3\nPBN6t6brJdV4fepa2j+fxrRVW8MO66Q8wTjnXB5Rokgiz9zUkPd6pJIQF0fX12by+MgF7Nwfm5Nn\neoJxzrk85pJaZRnXqxUPtKnF++kbuLL/ZCYs/jrssH7CE4xzzuVBRRLjeeKaenz0UEvKFC9Ej3dm\n8/C7c9i65+CpLz5LPME451we1qhKacY8chm/ueoCJiz+hnb9J/Ph3I3EwiwtnmCccy6PS4yP4+Er\najO212XUKlec3u/N5963vmTzd+FOnukJxjnn8onzy5fgg1+14A831GfGmu1c2X8y78xYR2ZmOLWZ\nqCYYSb0lLZa0SNJwSUUkXSFpTrBvqKSE4NxSksZImh9c0+0EZTaVtFDSKkkDFSygIKmMpImSVgY/\n887qPc45l0vi40S3ljWZ0Ls1F1VL4ncfLaLLkBmsydhz1mOJWoKRVBnoCaSYWQMgHrgDGAp0Cfat\nA+4OLnkIWGJmjYHLgf+TVOg4RQ8GegC1g8/Vwf7HgUlmVhuYFHx3zrkCqWqZYrxzX3Oe7dyIZV/v\n4poXpvDy5LM7eWa0m8gSgKJBLaUYsBc4aGYrguMTgU7BtgElghrJOUSWPf7Rq6qSKgIlzWy6RXqw\n3gZuCg53IJK8CH7ehHPOFWCSuDWlKp/2acPldZLpN24ZN/1jGiNnb+Slz1cxe110p52JWoIxs03A\nc8B6YAuwE3gfSJSUtY5AZ6BqsD0IqAdsBhYCvczs2FRbGdiY7fvGYB9ABTPbEtx7C1A+Vx/IOefy\nqPIli/DKnSkM7noxG7bv538+mM9z45fT9bUZUU0y0WwiSyJSq6gJVAKKA12BLsAASbOA3fxQS2kP\nzAvObQIMklTy2GKPc6vT6r2S1ENSuqT0jIyM07nUOefytGsaVuSuS6sDkV+ch49kMmPNtqjdL5pN\nZO2AtWaWYWaHgVFAi6B5q5WZNQfSgJXB+d2AURaxClgL1D2mzI1AlWzfqxCp8QB8EzShZTWlfXu8\noMxsiJmlmFlKcnJyLjymc87lHZfXKU+RxDjiBYkJcaTWKhu1e0UzwawHUiUVC/pV2gJLJZUHkFQY\n6Au8nO38tsGxCkAdYE32AoOmr92SUoMy7wL+HRwezQ8DBu7Ott8551ygafUkhnVPpc9VdRjWPZWm\n1aM34DYhWgWb2UxJI4A5RJrB5gJDgGckXU8kuQ02s8+CS54G3pK0kEhTWF8z2wogaZ6ZNQnOexB4\nCygKjAs+AP2A9yXdRyRZ3RKtZ3POubysafWkqCaWLIqF6QTCkpKSYunp6WGH4ZxzeYqk2WaWcqrz\n/E1+55xzUeEJxjnnXFR4gnHOORcVnmCcc85FhScY55xzUVGgR5FJyiAy4eaZKAdszcVwckusxgWx\nG5vHdXo8rtOTH+OqbmanfFO9QCeYn0NSek6G6Z1tsRoXxG5sHtfp8bhOT0GOy5vInHPORYUnGOec\nc1HhCebMDQk7gBOI1bggdmPzuE6Px3V6Cmxc3gfjnHMuKrwG45xzLio8wZwBSVdLWi5plaTHw44H\nQNIbkr6VtCjsWLKTVFXS55KWSlosqVfYMQFIKiJplqT5QVx/Cjum7CTFS5or6eOwY8ki6StJCyXN\nkxQzs8RKKi1phKRlwb+zS2MgpjrBf6eszy5Jj4YdF4Ck3sG/+UWShksqErV7eRPZ6ZEUD6wAriSy\nANqXwO1mtiTkuFoDe4C3zaxBmLFkFyz+VtHM5kgqAcwGboqB/14CipvZHkmJwFQiy3TPCDOuLJL6\nAClASTO7Pux4IJJggJSsZTRihaShwBQze01SIaCYmX0XdlxZgt8Zm4BLzOxM37vLrVgqE/m3Xt/M\n9kt6HxhrZm9F435egzl9zYFVZrbGzA4B/yKyNHSozCwN2B52HMcysy1mNifY3g0sBSqHGxUEK6fu\nCb4mBp+Y+GtLUhXgOuC1sGOJdcGy6q2B1wHM7FAsJZdAW2B12MklmwSgqKQEoBg/rAqc6zzBnL7K\nwIZs3zcSA78w8wJJNYCLgJnhRhIRNEPNI7K89kQzi4m4gOeB3wKZYQdyDAMmSJotqUfYwQRqARnA\nm0GT4muSiocd1DG6AMPDDgLAzDYBzxFZlHELsNPMJkTrfp5gTp+Osy8m/vKNZZLOAUYCj5rZrrDj\nATCzo8FKqVWA5pJCb1oMVnv91sxmhx3LcbQ0s4uBa4CHgmbZsCUAFxNZHfciYC8QE/2iAEGT3Y3A\nB2HHAiApiUiLS02gElBc0i+jdT9PMKdvI1A12/cqRLGKmR8EfRwjgWFmNirseI4VNKn8F7g65FAA\nWgI3Bv0d/wKukPTPcEOKMLPNwc9vgQ+JNBeHbSOwMVvtcwSRhBMrrgHmmNk3YQcSaAesNbMMMzsM\njAJaROtmnmBO35dAbUk1g79OugCjQ44pZgWd6a8DS82sf9jxZJGULKl0sF2UyP94y8KNCszsCTOr\nYmY1iPzb+szMovYXZk5JKh4M0iBogroKCH3Eopl9DWyQVCfY1RYIdQDJMW4nRprHAuuBVEnFgv83\n2xLpF42KhGgVnF+Z2RFJDwPjgXjgDTNbHHJYSBoOXA6Uk7QR+IOZvR5uVEDkL/I7gYVBfwfAk2Y2\nNsSYACoCQ4MRPnHA+2YWM0OCY1AF4MPI7yQSgHfN7JNwQ/reI8Cw4A++NUC3kOMBQFIxIqNNHwg7\nlixmNlPSCGAOcASYSxTf6Pdhys4556LCm8icc85FhScY55xzUeEJxjnnXFR4gnHOORcVnmCcc85F\nhScY586QpD2nPutnlf+apPrB9pNncH2NWJtd2xUsPkzZuTMkaY+ZnROr9wrmfvs4lmbXdgWL12Cc\ny0WSqkuaJGlB8LNasP8tSQMlfSFpjaTOwf44Sf8I1uf4WNLYbMf+KylFUj8is9/OkzTs2JqJpN9I\n+mOw3TRY42Y68FC2c+Il/V3Sl0FsMfPyn8u/PME4l7sGEVmTpxEwDBiY7VhF4DLgeqBfsK8jUANo\nCHQHfrJYlpk9Duw3syZm1vUU938T6Glmx5ZzH5GZc5sBzYD7JdU8nQdz7nR5gnEud10KvBtsv0Mk\noWT5yMwyg8XWKgT7LgM+CPZ/DXx+pjeWVAoobWaTs90/y1XAXcF0PTOBskDtM72Xcznhc5E5F13Z\nOzkPZtvWMT9PxxF+/Mdh1pK34sRLRwh4xMzGn8H9nDsjXoNxLnd9QWQWZICuRJanPZmpQKegL6YC\nkQlLj+dwsOwBwDdAeUllJRUm0uSWtezATklZtabszWnjgQezypB0QQwuzOXyGa/BOHfmigUzV2fp\nD/QE3pD0GJGVFk81s+9IIlOmLwJWEGm+2nmc84YACyTNMbOukp4Kzl3Lj5cZ6Bbcfx+RpJLlNSJ9\nPXOCadozgJty9JTOnSEfpuxcyCSdY2Z7JJUFZhFZOfLrsONy7ufyGoxz4fs4WPysEPC0JxeXX3gN\nxjnnXFR4J79zzrmo8ATjnHMuKjzBOOeciwpPMM4556LCE4xzzrmo8ATjnHMuKv4/I8pZwrO8fZ0A\nAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x123733b90>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<module 'matplotlib.pyplot' from '/Users/erik/miniconda2/envs/py2_parcels/lib/python2.7/site-packages/matplotlib/pyplot.pyc'>"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[900])\n",
+    "output_file = pset.ParticleFile(name='FieldListParticle_adv_stokes.nc', outputdt=1)\n",
+    "pset.execute(AdvectionRK4, runtime=10, dt=1, output_file=output_file)\n",
+    "plotTrajectoriesFile('FieldListParticle_adv_stokes.nc')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "What happens under the hood is that each `Field` in the `FieldList` is interpolated separately to the particle location, and that the different velocities are added in each step of the RK4 advection. So `FieldLists` are effortless to users\n",
+    "\n",
+    "Note that `FieldLists` work for any type of `Field`, not only for velocities. Any call to a `Field` interpolation (`fieldset.fld[time, lon, lat, depth]`) will return the sum of all `Fields` in the list if `fld` is a `FieldList`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "language": "python",
+   "name": "python2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/parcels/examples/tutorial_periodic_boundaries.ipynb b/parcels/examples/tutorial_periodic_boundaries.ipynb
index bbf4c993..aa22149e 100644
--- a/parcels/examples/tutorial_periodic_boundaries.ipynb
+++ b/parcels/examples/tutorial_periodic_boundaries.ipynb
@@ -88,7 +88,11 @@
     "editable": true
    },
    "source": [
-    "Extending the fieldset with a halo is very simply done using the `add_periodic_halo()` method. Halos can be added either in the zonal direction, the meridional direction, or both, by setting `zonal` and/or `meridional` to `True`"
+    "Extending the fieldset with a halo is very simply done using the `add_periodic_halo()` method. Halos can be added either in the zonal direction, the meridional direction, or both, by setting `zonal` and/or `meridional` to `True`.\n",
+    "\n",
+    "But before we apply the halo, we first define two new fieldset constants `halo_east` and `halo_west`. They store the original zonal extend of the grid (so *before* adding the halo) and will be used later in the `periodicBC` kernel.\n",
+    "\n",
+    "***Note that some hydrodynamic data, such as the global ORCA grid used in NEMO, already has a halo.*** In these cases, **do not** extent the fieldset with the halo but only add the periodic boundary kernel, where you use the explicit values for halo_east and halo_west"
    ]
   },
   {
@@ -101,6 +105,9 @@
    },
    "outputs": [],
    "source": [
+    "fieldset.add_constant('halo_west', fieldset.U.grid.lon[0])\n",
+    "fieldset.add_constant('halo_east', fieldset.U.grid.lon[-1])\n",
+    "\n",
     "fieldset.add_periodic_halo(zonal=True)"
    ]
   },
@@ -133,22 +140,7 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
-   "source": [
-    "The fieldset constants `halo_east` and `halo_west` in the kernel above were set during the call to `add_periodic_halo()`, and store the original zonal extent of the fieldset. If we would have also used `meridional=True` in `add_periodic_halo()`, `fieldset.halo_north` and `fieldset.halo_south` would also have been added.\n",
-    "\n",
-    "***Note that some hydrodynamic data, such as the global ORCA grid used in NEMO, already has a halo.*** In these cases, **do not** extent the fieldset with the halo but only add the periodic boundary kernel, where you use the explicit values for halo_east and halo_west"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "deletable": true,
-    "editable": true
-   },
+   "metadata": {},
    "source": [
     "Now define a particle set and execute it as usual"
    ]
diff --git a/parcels/field.py b/parcels/field.py
index 98374ec4..c23e9075 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -160,6 +160,8 @@ class Field(object):
            2. flat: No conversion, lat/lon are assumed to be in m.
     :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 FieldLists
+           (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
@@ -176,7 +178,7 @@ class Field(object):
                       'Kh_meridional': GeographicSquare()}
 
     def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=None, mesh='flat',
-                 transpose=False, vmin=None, vmax=None, time_origin=None,
+                 fieldtype=None, transpose=False, vmin=None, vmax=None, time_origin=None,
                  interp_method='linear', allow_time_extrapolation=None, time_periodic=False, **kwargs):
         self.name = name
         self.data = data
@@ -192,10 +194,11 @@ class Field(object):
         self.lat = self.grid.lat
         self.depth = self.grid.depth
         self.time = self.grid.time
-        if self.grid.mesh is 'flat' or (name not in self.unitconverters.keys()):
+        fieldtype = self.name if fieldtype is None else fieldtype
+        if self.grid.mesh is 'flat' or (fieldtype not in self.unitconverters.keys()):
             self.units = UnitConverter()
         elif self.grid.mesh is 'spherical':
-            self.units = self.unitconverters[name]
+            self.units = self.unitconverters[fieldtype]
         else:
             raise ValueError("Unsupported mesh type. Choose either: 'spherical' or 'flat'")
         self.interp_method = interp_method
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index d9593cbc..7904c05a 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -9,7 +9,7 @@ from glob import glob
 from copy import deepcopy
 
 
-__all__ = ['FieldSet']
+__all__ = ['FieldSet', 'FieldList']
 
 
 class FieldSet(object):
@@ -22,14 +22,15 @@ class FieldSet(object):
     def __init__(self, U, V, fields=None):
         self.gridset = GridSet()
         if U:
-            self.add_field(U)
+            self.add_field(U, 'U')
+            self.time_origin = self.U.grid.time_origin if isinstance(self.U, Field) else self.U[0].grid.time_origin
         if V:
-            self.add_field(V)
+            self.add_field(V, 'V')
 
         # Add additional fields as attributes
         if fields:
             for name, field in fields.items():
-                self.add_field(field)
+                self.add_field(field, name)
 
     @classmethod
     def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
@@ -83,14 +84,22 @@ class FieldSet(object):
         v = fields.pop('V', None)
         return cls(u, v, fields=fields)
 
-    def add_field(self, field):
+    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
         """
-        setattr(self, field.name, field)
-        self.gridset.add_grid(field)
-        field.fieldset = self
+        name = field.name if name is None else name
+        if isinstance(field, list):
+            setattr(self, name, FieldList(field))
+            for fld in field:
+                self.gridset.add_grid(fld)
+                fld.fieldset = self
+        else:
+            setattr(self, name, field)
+            self.gridset.add_grid(field)
+            field.fieldset = self
 
     def add_vector_field(self, vfield):
         """Add a :class:`parcels.field.VectorField` object to the FieldSet
@@ -101,23 +110,32 @@ class FieldSet(object):
         vfield.fieldset = self
 
     def check_complete(self):
-        assert(self.U), ('U field is not defined')
-        assert(self.V), ('V field is not defined')
-        ugrid = self.U.grid
+        assert self.U, 'FieldSet does not have a Field named "U"'
+        assert self.V, 'FieldSet does not have a Field named "V"'
+        for attr, value in vars(self).items():
+            if type(value) is Field:
+                assert value.name == attr, 'Field %s.name (%s) is not consistent' % (value.name, attr)
+
         for g in self.gridset.grids:
             g.check_zonal_periodic()
-            if g is ugrid or len(g.time) == 1:
+            if len(g.time) == 1:
                 continue
-            assert isinstance(g.time_origin, type(ugrid.time_origin)), 'time origins of different grids must be have the same type'
+            assert isinstance(g.time_origin, type(self.time_origin)), 'time origins of different grids must be have the same type'
             if g.time_origin:
-                g.time = g.time + (g.time_origin - ugrid.time_origin) / np.timedelta64(1, 's')
+                g.time = g.time + (g.time_origin - self.time_origin) / np.timedelta64(1, 's')
                 if g.defer_load:
-                    g.time_full = g.time_full + (g.time_origin - ugrid.time_origin) / np.timedelta64(1, 's')
-                g.time_origin = ugrid.time_origin
+                    g.time_full = g.time_full + (g.time_origin - self.time_origin) / np.timedelta64(1, 's')
+                g.time_origin = self.time_origin
         if not hasattr(self, 'UV'):
-            self.add_vector_field(VectorField('UV', self.U, self.V))
+            if isinstance(self.U, FieldList):
+                self.add_vector_field(VectorFieldList('UV', self.U, self.V))
+            else:
+                self.add_vector_field(VectorField('UV', self.U, self.V))
         if not hasattr(self, 'UVW') and hasattr(self, 'W'):
-            self.add_vector_field(VectorField('UVW', self.U, self.V, self.W))
+            if isinstance(self.U, FieldList):
+                self.add_vector_field(VectorFieldList('UVW', self.U, self.V, self.W))
+            else:
+                self.add_vector_field(VectorField('UVW', self.U, self.V, self.W))
 
     @classmethod
     def from_netcdf(cls, filenames, variables, dimensions, indices=None,
@@ -314,14 +332,6 @@ class FieldSet(object):
         :param halosize: size of the halo (in grid points). Default is 5 grid points
         """
 
-        # setting FieldSet constants for use in PeriodicBC kernel. Note using U-Field values
-        if zonal:
-            self.add_constant('halo_west', self.U.grid.lon[0])
-            self.add_constant('halo_east', self.U.grid.lon[-1])
-        if meridional:
-            self.add_constant('halo_south', self.U.grid.lat[0])
-            self.add_constant('halo_north', self.U.grid.lat[-1])
-
         for grid in self.gridset.grids:
             grid.add_periodic_halo(zonal, meridional, halosize)
         for attr, value in iter(self.__dict__.items()):
@@ -426,3 +436,60 @@ class FieldSet(object):
                 return nextTime
             else:
                 return time + nSteps * dt
+
+
+class FieldList(list):
+    def eval(self, time, x, y, z, applyConversion=True):
+        tmp = 0
+        for fld in self:
+            tmp += fld.eval(time, x, y, z, applyConversion)
+        return tmp
+
+    def __getitem__(self, key):
+        if isinstance(key, int):
+            return list.__getitem__(self, key)
+        else:
+            return self.eval(*key)
+
+
+class VectorFieldList(list):
+    """Class VectorFieldList stores 2 or 3 FieldLists which defines together a vector field.
+    This enables to interpolate them as one single vector FieldList in the kernels.
+
+    :param name: Name of the vector field
+    :param U: FieldList defining the zonal component
+    :param V: FieldList defining the meridional component
+    :param W: FieldList defining the vertical component (default: None)
+    """
+
+    def __init__(self, name, U, V, W=None):
+        self.name = name
+        self.U = U
+        self.V = V
+        self.W = W
+
+    def eval(self, time, x, y, z):
+        zonal = meridional = vertical = 0
+        if self.W is not None:
+            for (U, V, W) in zip(self.U, self.V, self.W):
+                vfld = VectorField(self.name, U, V, W)
+                vfld.fieldset = self.fieldset
+                (tmp1, tmp2, tmp3) = vfld.eval(time, x, y, z)
+                zonal += tmp1
+                meridional += tmp2
+                vertical += tmp3
+            return (zonal, meridional, vertical)
+        else:
+            for (U, V) in zip(self.U, self.V):
+                vfld = VectorField(self.name, U, V)
+                vfld.fieldset = self.fieldset
+                (tmp1, tmp2) = vfld.eval(time, x, y, z)
+                zonal += tmp1
+                meridional += tmp2
+            return (zonal, meridional)
+
+    def __getitem__(self, key):
+        if isinstance(key, int):
+            return list.__getitem__(self, key)
+        else:
+            return self.eval(*key)
diff --git a/parcels/gridset.py b/parcels/gridset.py
index ce530940..79ca9fbe 100644
--- a/parcels/gridset.py
+++ b/parcels/gridset.py
@@ -35,3 +35,15 @@ class GridSet(object):
             self.grids.append(grid)
             self.size += 1
         field.igrid = self.grids.index(field.grid)
+
+    def dimrange(self, dim):
+        """Returns maximum value of a dimension (lon, lat, depth or time)
+           on 'left' side and minimum value on 'right' side for all grids
+           in a gridset. Useful for finding e.g. longitude range that
+           overlaps on all grids in a gridset"""
+
+        maxleft, minright = (0, np.infty)
+        for g in self.grids:
+            maxleft = max(maxleft, getattr(g, dim)[0])
+            minright = min(minright, getattr(g, dim)[-1])
+        return maxleft, minright
diff --git a/parcels/kernels/error.py b/parcels/kernels/error.py
index d7b36395..7f545f88 100644
--- a/parcels/kernels/error.py
+++ b/parcels/kernels/error.py
@@ -29,8 +29,8 @@ class KernelError(RuntimeError):
 
 
 def parse_particletime(time, fieldset):
-    if fieldset is not None and fieldset.U.grid.time_origin:
-        time = fieldset.U.grid.time_origin + np.timedelta64(int(time), 's')
+    if fieldset is not None and fieldset.time_origin:
+        time = fieldset.time_origin + np.timedelta64(int(time), 's')
     return time
 
 
diff --git a/parcels/particleset.py b/parcels/particleset.py
index 89702501..882039d0 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -47,16 +47,18 @@ class ParticleSet(object):
 
         lon = convert_to_list(lon)
         lat = convert_to_list(lat)
-        depth = np.ones(len(lon)) * fieldset.U.grid.depth[0] if depth is None else depth
+        if depth is None:
+            mindepth, _ = self.fieldset.gridset.dimrange('depth')
+            depth = np.ones(len(lon)) * mindepth
         depth = convert_to_list(depth)
         assert len(lon) == len(lat) and len(lon) == len(depth)
 
         time = time.tolist() if isinstance(time, np.ndarray) else time
         time = [time] * len(lat) if not isinstance(time, list) else time
         time = [np.datetime64(t) if isinstance(t, datetime) else t for t in time]
-        self.time_origin = fieldset.U.grid.time_origin
+        self.time_origin = fieldset.time_origin
         if len(time) > 0 and isinstance(time[0], np.timedelta64) and not self.time_origin:
-            raise NotImplementedError('If fieldset.U.grid.time_origin is not a date, time of a particle must be a double')
+            raise NotImplementedError('If fieldset.time_origin is not a date, time of a particle must be a double')
         time = [((t - self.time_origin) / np.timedelta64(1, 's')) if isinstance(t, np.datetime64) else t for t in time]
 
         assert len(lon) == len(time)
@@ -159,8 +161,8 @@ class ParticleSet(object):
             p = np.reshape(start_field.data, (1, start_field.data.size))
             inds = np.random.choice(start_field.data.size, size, replace=True, p=p[0] / np.sum(p))
             lat, lon = np.unravel_index(inds, start_field.data[0, :, :].shape)
-            lon = fieldset.U.grid.lon[lon]
-            lat = fieldset.U.grid.lat[lat]
+            lon = start_field.grid.lon[lon]
+            lat = start_field.grid.lat[lat]
             for i in range(lon.size):
                 lon[i] = add_jitter(lon[i], lonwidth, start_field.grid.lon[0], start_field.grid.lon[-1])
                 lat[i] = add_jitter(lat[i], latwidth, start_field.grid.lat[0], start_field.grid.lat[-1])
@@ -266,7 +268,7 @@ class ParticleSet(object):
             endtime = np.datetime64(endtime)
         if isinstance(endtime, np.datetime64):
             if not self.time_origin:
-                raise NotImplementedError('If fieldset.U.grid.time_origin is not a date, execution endtime must be a double')
+                raise NotImplementedError('If fieldset.time_origin is not a date, execution endtime must be a double')
             endtime = (endtime - self.time_origin) / np.timedelta64(1, 's')
         if isinstance(runtime, delta):
             runtime = runtime.total_seconds()
@@ -285,7 +287,8 @@ class ParticleSet(object):
         # Set particle.time defaults based on sign of dt, if not set at ParticleSet construction
         for p in self:
             if np.isnan(p.time):
-                p.time = self.fieldset.U.grid.time[0] if dt >= 0 else self.fieldset.U.grid.time[-1]
+                mintime, maxtime = self.fieldset.gridset.dimrange('time')
+                p.time = mintime if dt >= 0 else maxtime
 
         # Derive _starttime and endtime from arguments or fieldset defaults
         if runtime is not None and endtime is not None:
@@ -296,7 +299,8 @@ class ParticleSet(object):
         if runtime is not None:
             endtime = _starttime + runtime * np.sign(dt)
         elif endtime is None:
-            endtime = self.fieldset.U.grid.time[-1] if dt >= 0 else self.fieldset.U.grid.time[0]
+            mintime, maxtime = self.fieldset.gridset.dimrange('time')
+            endtime = maxtime if dt >= 0 else mintime
 
         if abs(endtime-_starttime) < 1e-5 or dt == 0 or runtime == 0:
             dt = 0
@@ -385,12 +389,12 @@ class ParticleSet(object):
             show_time = np.datetime64(show_time)
         if isinstance(show_time, np.datetime64):
             if not self.time_origin:
-                raise NotImplementedError('If fieldset.U.grid.time_origin is not a date, showtime cannot be a date in particleset.show()')
+                raise NotImplementedError('If fieldset.time_origin is not a date, showtime cannot be a date in particleset.show()')
             show_time = (show_time - self.time_origin) / np.timedelta64(1, 's')
         if isinstance(show_time, delta):
             show_time = show_time.total_seconds()
         if np.isnan(show_time):
-            show_time = self.fieldset.U.grid.time[0]
+            show_time, _ = self.fieldset.gridset.dimrange('time')
         if domain is not None:
             def nearest_index(array, value):
                 """returns index of the nearest value in array using O(log n) bisection method"""


Test patch:
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index afb3dde9..bca477fb 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -1,4 +1,5 @@
-from parcels import FieldSet, Field, ParticleSet, ScipyParticle, JITParticle, Geographic, AdvectionRK4, Variable
+from parcels import (FieldSet, Field, ParticleSet, ScipyParticle, JITParticle, Geographic,
+                     AdvectionRK4, AdvectionRK4_3D, Variable)
 import numpy as np
 import pytest
 from math import cos, pi
@@ -403,3 +404,49 @@ def test_sampling_multiple_grid_sizes(mode):
 
     pset.execute(AdvectionRK4, runtime=10, dt=1)
     assert np.isclose(pset[0].lon, 0.8)
+
+
+@pytest.mark.parametrize('mode', ['jit', 'scipy'])
+@pytest.mark.parametrize('with_W', [True, False])
+@pytest.mark.parametrize('mesh', ['flat', 'spherical'])
+def test_list_of_fields(mode, with_W, k_sample_p, mesh):
+    xdim = 10
+    ydim = 20
+    zdim = 4
+    gf = 10  # factor by which the resolution of grid1 is higher than of grid2
+    U1 = Field('U1', 0.2*np.ones((zdim*gf, ydim*gf, xdim*gf), dtype=np.float32),
+               lon=np.linspace(0., 1., xdim*gf, dtype=np.float32),
+               lat=np.linspace(0., 1., ydim*gf, dtype=np.float32),
+               depth=np.linspace(0., 20., zdim*gf, dtype=np.float32),
+               mesh=mesh, fieldtype='U')
+    U2 = Field('U2', 0.1*np.ones((zdim, ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 1., xdim, dtype=np.float32),
+               lat=np.linspace(0., 1., ydim, dtype=np.float32),
+               depth=np.linspace(0., 20., zdim, dtype=np.float32),
+               mesh=mesh, fieldtype='U')
+    V1 = Field('V1', np.zeros((zdim*gf, ydim*gf, xdim*gf), dtype=np.float32), grid=U1.grid, fieldtype='V')
+    V2 = Field('V2', np.zeros((zdim, ydim, xdim), dtype=np.float32), grid=U2.grid, fieldtype='V')
+    fieldset = FieldSet([U1, U2], [V1, V2])
+
+    conv = 1852*60 if mesh == 'spherical' else 1.
+    assert np.allclose(fieldset.U.eval(0, 0, 0, 0)*conv, 0.3)
+    assert np.allclose(fieldset.U[0, 0, 0, 0]*conv, 0.3)
+
+    P1 = Field('P1', 30*np.ones((zdim*gf, ydim*gf, xdim*gf), dtype=np.float32), grid=U1.grid)
+    P2 = Field('P2', 20*np.ones((zdim, ydim, xdim), dtype=np.float32), grid=U2.grid)
+    fieldset.add_field([P1, P2], name='P')
+    assert np.allclose(fieldset.P[0, 0, 0, 0], 50)
+
+    if with_W:
+        W1 = Field('W1', 2*np.ones((zdim * gf, ydim * gf, xdim * gf), dtype=np.float32), grid=U1.grid)
+        W2 = Field('W2', np.ones((zdim, ydim, xdim), dtype=np.float32), grid=U2.grid)
+        fieldset.add_field([W1, W2], name='W')
+        pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0], lat=[0.9])
+        pset.execute(AdvectionRK4_3D+pset.Kernel(k_sample_p), runtime=2, dt=1)
+        assert np.isclose(pset[0].depth, 6)
+    else:
+        pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0], lat=[0.9])
+        pset.execute(AdvectionRK4+pset.Kernel(k_sample_p), runtime=2, dt=1)
+    assert np.isclose(pset[0].p, 50)
+    assert np.isclose(pset[0].lon*conv, 0.6, atol=1e-3)
+    assert np.isclose(pset[0].lat, 0.9)
diff --git a/tests/test_grids.py b/tests/test_grids.py
index 9ee857ae..a839f446 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -101,7 +101,7 @@ def test_avoid_repeated_grids():
     temp0_field = Field('temp', u_data, lon=lon_g0, lat=lat_g0, time=time_g0, transpose=True)
 
     other_fields = {}
-    other_fields['temp0'] = temp0_field
+    other_fields['temp'] = temp0_field
 
     field_set = FieldSet(u_field, v_field, fields=other_fields)
     assert field_set.gridset.size == 2