Repository: OceanParcels/parcels
Pull request #488 (2018-11-12 13:20:52)

Nested fields


Patch:
diff --git a/environment_py2_win.yml b/environment_py2_win.yml
index b0fc4c11..bb8408a7 100644
--- a/environment_py2_win.yml
+++ b/environment_py2_win.yml
@@ -22,6 +22,7 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - dask
   - cftime
   - pip:
       - pytest>=2.7.0
diff --git a/environment_py3_win.yml b/environment_py3_win.yml
index de2f7e0b..473c7e75 100644
--- a/environment_py3_win.yml
+++ b/environment_py3_win.yml
@@ -21,6 +21,7 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - dask
   - cftime
   - ipykernel<5.0
   - pip:
diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py
index fe215f36..ec812c9e 100644
--- a/parcels/codegenerator.py
+++ b/parcels/codegenerator.py
@@ -1,4 +1,4 @@
-from parcels.field import Field, VectorField, SummedField, SummedVectorField
+from parcels.field import Field, VectorField, SummedField, SummedVectorField, NestedField
 from parcels.tools.loggers import logger
 import ast
 import cgen as c
@@ -19,6 +19,13 @@ 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), NestedField):
+            if isinstance(getattr(self.obj, attr)[0], VectorField):
+                return NestedVectorFieldNode(getattr(self.obj, attr),
+                                             ccode="%s->%s" % (self.ccode, attr))
+            else:
+                return NestedFieldNode(getattr(self.obj, attr),
+                                       ccode="%s->%s" % (self.ccode, attr))
         elif isinstance(getattr(self.obj, attr), SummedVectorField):
             return SummedVectorFieldNode(getattr(self.obj, attr),
                                          ccode="%s->%s" % (self.ccode, attr))
@@ -85,6 +92,32 @@ class SummedVectorFieldEvalNode(IntrinsicNode):
         self.var3 = var3  # third variable for UVW interpolation
 
 
+class NestedFieldNode(IntrinsicNode):
+    def __getitem__(self, attr):
+        return NestedFieldEvalNode(self.obj, attr)
+
+
+class NestedFieldEvalNode(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 NestedVectorFieldNode(IntrinsicNode):
+    def __getitem__(self, attr):
+        return NestedVectorFieldEvalNode(self.obj, attr)
+
+
+class NestedVectorFieldEvalNode(IntrinsicNode):
+    def __init__(self, fields, args, var, var2, var3):
+        self.fields = fields
+        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
@@ -248,6 +281,19 @@ class IntrinsicTransformer(ast.NodeTransformer):
                 return ast.Tuple([ast.Name(id=tmp), ast.Name(id=tmp2), ast.Name(id=tmp3)], ast.Load())
             else:
                 return ast.Tuple([ast.Name(id=tmp), ast.Name(id=tmp2)], ast.Load())
+        elif isinstance(node.value, NestedFieldNode):
+            tmp = self.get_tmp()
+            self.stmt_stack += [NestedFieldEvalNode(node.value, node.slice, tmp)]
+            return ast.Name(id=tmp)
+        elif isinstance(node.value, NestedVectorFieldNode):
+            tmp = self.get_tmp()
+            tmp2 = self.get_tmp()
+            tmp3 = self.get_tmp() if list.__getitem__(node.value.obj, 0).W else None
+            self.stmt_stack += [NestedVectorFieldEvalNode(node.value, node.slice, tmp, tmp2, tmp3)]
+            if tmp3:
+                return ast.Tuple([ast.Name(id=tmp), ast.Name(id=tmp2), ast.Name(id=tmp3)], ast.Load())
+            else:
+                return ast.Tuple([ast.Name(id=tmp), ast.Name(id=tmp2)], ast.Load())
         else:
             return node
 
@@ -633,6 +679,11 @@ class KernelGenerator(ast.NodeVisitor):
         for fld in node.obj:
             self.field_args[fld.name] = fld
 
+    def visit_NestedFieldNode(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
@@ -647,6 +698,11 @@ class KernelGenerator(ast.NodeVisitor):
             for fld in node.obj.W:
                 self.field_args[fld.name] = fld
 
+    def visit_NestedVectorFieldNode(self, node):
+        """Record intrinsic fields used in kernel"""
+        for fld in node.obj:
+            self.vector_field_args[fld.name] = fld
+
     def visit_ConstNode(self, node):
         self.const_args[node.ccode] = node.obj
 
@@ -712,6 +768,44 @@ class KernelGenerator(ast.NodeVisitor):
             cstat += [c.Assign("err", ccode_eval), conv_stat, c.Statement("CHECKERROR(err)")]
         node.ccode = c.Block(cstat)
 
+    def visit_NestedFieldEvalNode(self, node):
+        self.visit(node.fields)
+        self.visit(node.args)
+        cstat = []
+        for fld in node.fields.obj:
+            ccode_eval = fld.ccode_eval(node.var, *node.args.ccode)
+            ccode_conv = fld.ccode_convert(*node.args.ccode)
+            conv_stat = c.Statement("%s *= %s" % (node.var, ccode_conv))
+            cstat += [c.Assign("err", ccode_eval),
+                      conv_stat,
+                      c.If("err != ERROR_OUT_OF_BOUNDS ", c.Block([c.Statement("CHECKERROR(err)"), c.Statement("break")]))]
+        cstat += [c.Statement("CHECKERROR(err)"), c.Statement("break")]
+        node.ccode = c.While("1==1", c.Block(cstat))
+
+    def visit_NestedVectorFieldEvalNode(self, node):
+        self.visit(node.fields)
+        self.visit(node.args)
+        cstat = []
+        for fld in node.fields.obj:
+            ccode_eval = fld.ccode_eval(node.var, node.var2, node.var3,
+                                        fld.U, fld.V, fld.W,
+                                        *node.args.ccode)
+            if fld.U.interp_method != 'cgrid_linear':
+                ccode_conv1 = fld.U.ccode_convert(*node.args.ccode)
+                ccode_conv2 = fld.V.ccode_convert(*node.args.ccode)
+                statements = [c.Statement("%s *= %s" % (node.var, ccode_conv1)),
+                              c.Statement("%s *= %s" % (node.var2, ccode_conv2))]
+            else:
+                statements = []
+            if node.var3:
+                ccode_conv3 = fld.W.ccode_convert(*node.args.ccode)
+                statements.append(c.Statement("%s *= %s" % (node.var3, ccode_conv3)))
+            cstat += [c.Assign("err", ccode_eval),
+                      c.Block(statements),
+                      c.If("err != ERROR_OUT_OF_BOUNDS ", c.Block([c.Statement("CHECKERROR(err)"), c.Statement("break")]))]
+        cstat += [c.Statement("CHECKERROR(err)"), c.Statement("break")]
+        node.ccode = c.While("1==1", 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/tutorial_NestedFields.ipynb b/parcels/examples/tutorial_NestedFields.ipynb
new file mode 100644
index 00000000..40dcb7d7
--- /dev/null
+++ b/parcels/examples/tutorial_NestedFields.ipynb
@@ -0,0 +1,262 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Tutorial on how to combine different Fields into a `NestedField` object"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In some applications, you may have access to different fields that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a field covering the entire region and another one only covering part of it, but with a higher resolution. The set of those fields form what we call nested fields.\n",
+    "\n",
+    "It is possible to combine all those fields with kernels, either with different if/else statements depending on particle position, or using recovery kernels (if only two levels of nested fields).\n",
+    "\n",
+    "However, an easier way to work with nested fields in Parcels is to combine all those fields into one `NestedField` object. The Parcels code will then try to successively interpolate the different fields. \n",
+    "\n",
+    "For each Particle, the algorithm is the following:\n",
+    "\n",
+    "1. Interpolate the particle onto the first `Field` in the `NestedFields` list.\n",
+    "\n",
+    "2. If the interpolation succeeds or if an error other than `ErrorOutOfBounds` is thrown, the function is stopped.\n",
+    "\n",
+    "3. If an `ErrorOutOfBounds` is thrown, try step 1) again with the next `Field` in the `NestedFields` list \n",
+    "\n",
+    "4. If interpolation on the last `Field` in the `NestedFields` list also returns an `ErrorOutOfBounds`, then the Particle is flagged as OutOfBounds.\n",
+    "\n",
+    "This algorithm means that **the order of the fields in the `NestedField` matters**. In particular, the smallest/finest resolution fields have to be listed _before_ the larger/coarser resolution fields.\n",
+    "\n",
+    "This tutorial shows how to use these `NestedField` with a very idealised example."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from parcels import Field, NestedField, FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile, AdvectionRK4\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon / 200 * pi / 2) m/s."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dim = 21\n",
+    "lon = np.linspace(0., 2e3, dim, dtype=np.float32)\n",
+    "lat = np.linspace(0., 2e3, dim, dtype=np.float32)\n",
+    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
+    "V1_data = np.cos(lon_g / 200 * np.pi/2)\n",
+    "U1 = Field('U1', np.ones((dim, dim), dtype=np.float32), lon=lon, lat=lat)\n",
+    "V1 = Field('V1', V1_data, grid=U1.grid)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now define the same velocity field on a low resolution (dx = 2km) 20kmx4km grid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xdim = 11\n",
+    "ydim = 3\n",
+    "lon = np.linspace(-2e3, 18e3, xdim, dtype=np.float32)\n",
+    "lat = np.linspace(-1e3, 3e3, ydim, dtype=np.float32)\n",
+    "lon_g, lat_g = np.meshgrid(lon, lat)\n",
+    "V2_data = np.cos(lon_g / 200 * np.pi/2)\n",
+    "U2 = Field('U2', np.ones((ydim, xdim), dtype=np.float32), lon=lon, lat=lat)\n",
+    "V2 = Field('V2', V2_data, grid=U2.grid)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now combine those fields into a `NestedField` and create the fieldset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "U = NestedField('U', [U1, U2])\n",
+    "V = NestedField('V', [V1, V2])\n",
+    "fieldset = FieldSet(U, V)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/ff29215bba0bd5244da08955369801d5.so\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4HNWZ7/Hva3lhsUFiE8Y74CwGBscS2BOyKENYbzIGwuqZxNnGZAbuHZ7JvQ8kZOFCSMjcySSTSQIxgWASjHFiFg9xAoZBEBZvchy8sQhj2cLGBtwCK4CxpPf+UdVSqd3VUktd3S3593kePeo6VdX1nmqp3q5Tp+qYuyMiIpLNkFIHICIi5UtJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoQMOmb2dTP7RS+Wu8PMvlOMmGK2/3szm12q7Ufi6NX+kv2T6T4JKTYz2wxUA+3AX4AlwP9099Y+vFcd8Gt3H9uHde8Amt39G31Y14HJ7t6Y77qFZGbXAce7+9+XMg4ZvHQmIaXyaXcfCUwDTgH6cqAeWvCoCqScY4saKHFK6ShJSEm5+yvA74ETAczsC2a20cx2m9kmM7s8vayZ1ZlZs5ldbWavAneH6x5jZq3hzzFmdp2Z/Tqy3kfM7GkzazGzrWb2+WyxmNmnzGxNuNzTZvZXMcs9Eb78c7jNS7LE9kszqzKzB83sNTNLha/HRt6n3sy+HJn+Ylj3lJk9ZGYTIvNOMLOlZrbLzHaETURnA18HLgnj+HO47DFmtjhcttHM/iHyPteZ2W/N7Ndm9hbw+Sz7a0Zkf/05PFtLz/t8+LnsNrOXzezv4j9dGQyUJKSkzGwccC7wp7BoJ/Ap4BDgC8APzWxaZJWjgcOACcDngHOAbe4+MvzZlvH+4wkSyX8CRwJTgTVZ4pgG3A5cDhwO/BxYbGYjMpd194+FL08Ot3lPltjmEPx//TKcHg+8A/wkZj+cR3DAvyCM848ESRAzGwU8AvwBOAY4HnjU3f8AfBe4J4zj5PDt7gaaw2UvBL5rZqdHNjcT+C1QCdyVEccY4HfAd8K6/G9gkZkdaWYHAz8GznH3UcCHs+1LGVyUJKRU7jezFuBJ4HGCgx3u/jt3f8kDjwMPAx+NrNcBfNvd97j7O73Yzt8Bj7j73e6+193fcPdsB7Z/AH7u7svdvd3d5wF7gBl51KlbbOG2Frn72+6+G7gR+HjMupcD33P3je7eRrA/poZnE58CXnX3H7j7u+6+292XZ3uTMOl+BLg6XHYN8Avgs5HFnnH3+929I8s+/HtgibsvCecvBVYRJPJ0HU80swPdfbu7r89j/8gApCQhpXKeu1e6+wR3/6f0wcrMzjGzZWFTSQvBwemIyHqvufu7eWxnHPBSL5abAHw1bGJpCbc9juDbeG91i83MDjKzn5tZU9i08wRQaWYVMdv/j8i2dwEGjMmjDoTx7gqTUlpT+D5pW3OsPwG4KGM/fAQY7e5/AS4BvgJsN7PfmdkHehmXDFBKElI2wqadRcC/AdXuXknQ88kii2V2x+upe95W4LhebH4rcGOYuNI/B7n73b0MP1ssXwXeD0x390OAdDOVsa+twOUZ2z/Q3Z/uoQ6Z29wGHBY2UaWNB17JsU5mHL/KiONgd78JwN0fcvczgNHAc8CtOd5LBgElCSknw4ERwGtAm5mdA5zZwzo7gMPN7NCY+XcBnzSzi81sqJkdbmZTsyx3K/AVM5tugYPN7H9kHGwzt3tsD7GNIrgO0WJmhwHfzrHsLcDXzOwEADM71MwuCuc9CBxtZleZ2QgzG2Vm0yNxTDSzIQDuvhV4GviemR0QXnz/EhnXHnL4NfBpMzvLzCrC96gzs7FmVm1mfxtem9gDtBJ0Y5ZBTElCykbYRPK/gIVACpgFLO5hnecILtRuCptHjsmYv4WgyeqrBE04a4CTs7zPKoLrEj8Jt90IfD7Hpq8D5oXbvDhmmR8BBwKvA8sILjzH1eM+4PvAgrBpah3BRfn0fjkD+DTwKvAi8Ilw1d+Ev98ws9Xh68uAiQRnFfcRXCdZmqMu0Ti2ElzY/jpBst4K/B+CY8UQgv24jWBffhz4p968rwxcuplOpETCrrS/cPc7Sx2LSBydSYiUgJkdRNBc9XKpYxHJRUlCpMjM7CiCZqPHCboAi5QtNTeJiEgsnUmIiEisAf9wryOOOMInTpyYdd5f/vIXDj744OIG1AuKKz+KK3/lGpviyk+ScTU0NLzu7kf2uKC7D+ifmpoaj/PYY4/FzislxZUfxZW/co1NceUnybiAVd6LY6yam0REJJaShIiIxFKSEBGRWEoSIiISS0lCRERi9TtJhE+JXBEOc7jezP5vWD7JzJab2Ytmdo+ZDQ/LR4TTjeH8iZH3+lpY/ryZndXf2EREpH8KcSaxB/gbD4ZOnAqcbWYzCJ5o+UN3n0zwVM0vhct/CUi5+/HAD8PlMLMpwKXACcDZwM9iBmcREZEi6ffNdGF/29Zwclj448DfEDzqGWAewaOVbyZ4DPF1YflvgZ+YmYXlC9x9D/CymTUCpwLP9DfGbOY2zGX+2vlJvHWPWlpaqNxcWZJt56K48lOucUH5xqa48pMrrlknzWJOzZzEYyjIHdfhN/4GggHaf0ow1GKLB2P1QjAoe3r4xDGEwye6e5uZvUkw8PwYgmfuk2WdzO3NIRhonurqaurr67PG1draGjvvZ2t+RmNrI8ePPL5XdSyk9vZ2Wlpair7dniiu/JRrXFC+sSmu/MTF1djaSEtLC+/b/b7EYyhIknD3doJB2ysJBjn5YLbFwt/Zhm70HOXZtjcXmAtQW1vrdXV1WeOqr68nbl7l5kpqK2up/3x91vlJyhVXKSmu/JRrXFC+sSmu/MTFVXdHUFaMmAvau8ndW4B6YAbBgO/pJDSWYDQrCM4QxgGE8w8lGOWqszzLOiIiUgKF6N10ZHgGgZkdCHwS2Ag8BlwYLjYbeCB8vTicJpz/3+F1jcXApWHvp0nAZGBFf+MTEZG+K0Rz02iCsX4rCJLOQnd/0Mw2EIzX+x3gT8Bt4fK3Ab8KL0zvIujRhLuvN7OFwAagDbgibMYSEZESKUTvpmeBD2Up30TQOymz/F3gopj3uhG4sb8xiYhIYeiOaxERiaUkISIisZQkREQklpKEiIjEUpIQEZFYShIiIhJLSUJERGIpSYiISCwlCRERiaUkISIisZQkREQklpKEiIjEUpIQEZFYShIiIhJLSUJERGIpSYiISKxCDF86zsweM7ONZrbezP45LL/OzF4xszXhz7mRdb5mZo1m9ryZnRUpPzssazSza/obm4iI9E8hhi9tA77q7qvNbBTQYGZLw3k/dPd/iy5sZlMIhiw9ATgGeMTM3hfO/ilwBtAMrDSzxe6+oQAxiohIHxRi+NLtwPbw9W4z2wiMybHKTGCBu+8BXg7Huk4Pc9oYDnuKmS0Il1WSEBEpkUKcSXQys4kE410vB04DrjSzzwGrCM42UgQJZFlktWa6ksrWjPLpMduZA8wBqK6upr6+Pms8ra2tsfNaWloAYucnKVdcpaS48lOucUH5xqa48hMXVzGPXwVLEmY2ElgEXOXub5nZzcANgIe/fwB8EbAsqzvZr494tm25+1xgLkBtba3X1dVljam+vp64eZWbKwFi5ycpV1ylpLjyU65xQfnGprjyExdXMY9fBUkSZjaMIEHc5e73Arj7jsj8W4EHw8lmYFxk9bHAtvB1XLmIiJRAIXo3GXAbsNHd/z1SPjqy2PnAuvD1YuBSMxthZpOAycAKYCUw2cwmmdlwgovbi/sbn4iI9F0hziROAz4LrDWzNWHZ14HLzGwqQZPRZuByAHdfb2YLCS5ItwFXuHs7gJldCTwEVAC3u/v6AsQnIiJ9VIjeTU+S/TrDkhzr3AjcmKV8Sa71RESkuHTHtYiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJFYhhi8dZ2aPmdlGM1tvZv8clh9mZkvN7MXwd1VYbmb2YzNrNLNnzWxa5L1mh8u/aGaz+xubiIj0TyHOJNqAr7r7B4EZwBVmNgW4BnjU3ScDj4bTAOcQjGs9GZgD3AxBUgG+DUwHTgW+nU4sIiJSGoUYvnQ7sD18vdvMNgJjgJlAXbjYPKAeuDosv9PdHVhmZpVmNjpcdqm77wIws6XA2cDd/Y2xrDTOZerrP4NHKkuz/Ymz4Pg5pdm2iAw4/U4SUWY2EfgQsByoDhMI7r7dzI4KFxsDbI2s1hyWxZVn284cgrMQqqurqa+vzxpPa2tr7LyWlhaA2PlJmfr6zzj4vUZaWo4v6nYBRu5tpLWlhTXN78s6P9f+KiXFlb9yjU1x5ScurmIevwqWJMxsJLAIuMrd3zKz2EWzlHmO8n0L3ecCcwFqa2u9rq4u64bq6+uJm1e5OfgmHzc/MY9U0tJyPJUXrinudgEeqaOS+Drn2l+lpLjyV66xKa78xMVVzONXQXo3mdkwggRxl7vfGxbvCJuRCH/vDMubgXGR1ccC23KUi4hIiRSid5MBtwEb3f3fI7MWA+keSrOBByLlnwt7Oc0A3gybpR4CzjSzqvCC9ZlhmYiIlEghmptOAz4LrDWzdBvK14GbgIVm9iVgC3BROG8JcC7QCLwNfAHA3XeZ2Q3AynC569MXsUVEpDQK0bvpSbJfTwA4PcvyDlwR8163A7f3NyYRESkM3XEtIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhIrEINX3q7me00s3WRsuvM7BUzWxP+nBuZ9zUzazSz583srEj52WFZo5ldU4jYRESk7woxMh3AHcBPgDszyn/o7v8WLTCzKcClwAnAMcAjZva+cPZPgTMIxrteaWaL3X1DgWKUAWj+8i384um3mb9lFZd//DhqJlSVOqTEzV++hXtWbqH6kAP2izo3NKVYtLqZ13fv4chRI7hg2thBX+eBpCBJwt2fMLOJvVx8JrDA3fcAL5tZI3BqOK/R3TcBmNmCcFklif1M+qDR0LSL519tBWDThh08vGEHxx81ki+eNolZ08eXOMrCStd51eY3eGHHX8LSNwd9nW95/CUe2bADj5TftXwLp06s4upzPqhkUQYsGE20AG8UJIkH3f3EcPo64PPAW8Aq4KvunjKznwDL3P3X4XK3Ab8P3+Zsd/9yWP5ZYLq7X5llW3OAOQDV1dU1CxYsyBpTa2srI0eOzDrvqjVXAfCjqT/Ku679MfX1q2hvb2dt9X8WdbvpbQOsOSJ7nXPtr2JpTLXzveXv0t7DcudOHMrFHxhRlJjiFGp/JVHncvgss0nHVb9lL/M2vEeuo48Bs6cMp278sKLFVW7i4irE8esTn/hEg7vX9rRcoZqbsrkZuAHw8PcPgC+SfTxsJ/v1kax/Q+4+F5gLUFtb63V1dVkDqK+vJ25e5eZKgNj5iXmkkpaWluJvN9w2xNc51/4qhoamFL9asabHgyXAks1ttB10eEmbYwqxvxqaUtyx/E8Fr3OpP8s49fX1bDvwWOZtWJszQUDwzz9vw3tsoyrxz7mc91e2uIp5/Eqsd5O773D3dnfvAG6lq0mpGRgXWXQssC1HuewH5i/fwkW3PE3Trrd7vc7DG3Zw8c+fYf7yLQlGlpx0nZtT7/R6nYc37OCSuc/Q0JRKMLLk1G/Zy7X39Zwg0pyBX+eBLrEkYWajI5PnA+meT4uBS81shJlNAiYDK4CVwGQzm2Rmwwkubi9OKj4pHw1NKb55/1o6shw5zpxSzTemH8Cif/wwp0zc95tke4fzzfvXDrgDyPzlW7g2R50X/eOHY+vc1u7c8vhLRYiysOYv35K1iWlM1YH83fTxLPrHD/Pd80/i+CMP3qe5YaDWeTAoSHOTmd0N1AFHmFkz8G2gzsymEnwZ2AxcDuDu681sIcEF6TbgCndvD9/nSuAhoAK43d3XFyI+KW/3NjTTnnHkqBhi3DDzRGZNH099fT01E6r4zVc+zPzlW/jm/Wu7Ld/usGh184C5yNnQlOJbD6wj83JgtM5p6Tpnfvv+7+d20tCUGnh1zigfWmH8+NIPddajZkIVs6aPHxR1HiwK1bvpsizFt+VY/kbgxizlS4AlhYhJBo5nX+l+FnD8USP5/mf+KuvBYNb08bz/6FFcvehZGne2dpa/vntP4nEWyrJNb9AWOYUw4Iwp1bHt7umkce39azsTS0eHD6jE2Oc6RxLFQKvzYKE7rqWkbvvjJta+srtzemiFxSaItJoJVXz/M3/F0CFdjRKPbtw5YK5NbN31l27Tl3/sWOZ+rjZnnWdNH8+N551ERVhlB+5ZsXXA1Hnj9je7Tfe6zuefxJBonVcOnDoPFkoSUjLzl2/hht9t7Jw24OLacb36plgzoYpLTunq59DuzjfuX1v2B5D5y7ewYGVz57QBow7sXRfPWdPHc+mpXU1R7T4wrsf8elkTDz77aud0vnW+LFrnAXoNaiBTkpCSaGhK8c0H1nUrqxhifGba2F6/xwXTxnY7m+hw+NYD68r2ABJcoN+3zjOOPbzX73HBtLGdZxPQdT2mXKWvRUQN9joPNkoSUhLLNr1Be6SNeojB9TNPzKu9uWZCFdfPPBGLHkA6nGWb3ihkqAVz7+pm2r3/dT79g9Xdysr5esxvG7Z268FlDP46DzZKElISr6S63w8x56PH9umxE7Omj+fyjx7bOe3A7nf29je8gmtoSrFgRVdTWMUQ4zvnndSnOl/+8eMGxPWYhqYUC1du7ZyuGGLMnjK8MHV+rjzrPBgpSUjRNTSluDty8MinjTqbzHV/8eTLZdfkFJxFdE2f/oGj+vwspmzXY8qxmW1Rw7517usjNvapc0d51nkwUpKQort7RVO3ewTybaPONOPYw7t9y2wPu0qWk82vd+/RdMSo/j13KvN6TDnWeeOrb3Wb3h/qPBgpSUhRNTSluG/1K53TFUMs7zbqTOlrE+k/Zgd+29BcNt8yG5pSPBO5TjK0Ir8L9Nl0Xo8Jp8uxzmu2tHROF7LOaeVW58FKSUKKKnp3tQGXnDKuII/Azuweureto2y+ZS5cuaXz4m0+3Xx7Mmv6eC6q7TrwtrV1lM1F+189s7nzJrhC13nmycd0TpdTnQcrJQkpmoamFPes6roWMWzokH5/u4z6TM3YbjeblcO3zIamFL+NJKtC13nquK4DbwdQddDwgr13XzU0pfivP3c9m7PQdZ4eaZoslzoPZkoSUjTRRzMYcGFNYUcgq5lQxd+ePKZzuhy+Zd67upn2juB1EnVOvf1et4fhrdv2ZuyyxbIo42xxf6jzYKYkUUINTSl++lhjXt92812nL9tIykHDKzpfO3DiMYcWfBunTDqs83Wpv2U2NKVYmOCZEwQX7YdF7jQr9dnT/ljnwS7JQYckh4amFJf+/Bn2djhDjF71mf+PR1/gR0tfxIHhFcbdc/465ze0uU+8xE2/f44Oh6FDjKdq91B9SOlGc3t4fdejGYYQfCMstPS3zHR7eCm/ZS7b9AZ725M7c4Lg7Omi2nHcFd4zkL4Wc2aJnoF37+rmRM8WIb7OevBfMnQmUSL/+oeN7A3/mTo8eNplrpuD5j7xEj8MEwTAez08X//Xz2zmu0ue67xg2tbhvPx6K7vfbStQDfLT0JRi2aZdndNDhw7pV7fXOOX0LbMqcv9GUmdOEHQNTdc5fS2mMdWbse4Kq6EpxW8SPotIi3aHLZfrT4OVkkQJvNPmLH+5+x+0E//coYamFN/7/XP7lKefr59t+W8tzj4Ux2utpXmcwaMbd3Tr7ZLEN0zo+paZ1t5euusSG1/terptUmdOENT54ow6P7er+EmiGGdOadnqXOrrT4OVkkQJ7Ho3++CNcTcHLVi5ZZ8BaqDr+fqZ7m1ozjriGcDe9FXUItvWEgzRaQYjhiX3DRP2/WZdqusSTeEjwYcYDB+WzJlTWrcbzcwYOSzbUPLJmnHs4Z3P0RpWkexnDGFvtkid1cspGQVJEmZ2u5ntNLN1kbLDzGypmb0Y/q4Ky83MfmxmjWb2rJlNi6wzO1z+RTObXYjYys3ud9vY/V73wVfSsp02NzSluHd190dL57qBqqEpxT0NXaf8QwwiN6nS8vbeop+WNzSlWBx2iaww41ufOiHR9uOaCVV861NTgKAp7/oH1xe/zpt38ccXXgdgSJHq/H/Oej8QfNmY/9x7Ra/ze23tkS8nvR3Fuu9qJlTxlY9NAoIvTKX4nPcHhTqTuAM4O6PsGuBRd58MPBpOA5xDMK71ZGAOcDMESYVg2NPpwKnAt9OJZTB5LfL0SgMumz6ev4vcBJbZbTN4Wmr35aMXuDNvGlu26Q3aIqf8l57a/Xn87sV/lMGCFV03k7l7Ys0uUW9Frr3s2Vv8G+tuf6rrZrJi1Tk68tvejuI/TvuWxzd1vi7W03gPGhH0vXHgvb1qckpCQZKEuz8B7MoongnMC1/PA86LlN/pgWVApZmNBs4Clrr7LndPAUvZN/EMaA1NqW7XBNIX9k4Y03VBM7PbZvQibPri5wXTxlIxpKssejZRmeVi6QXTxjK0RBdzG5pS3PunroNVRUWyzS5pM449PHYfJa2hKcUf1m3vnC5qnUt0M2FDU4o/vvBa53Sx6lx1UFdvvVJ3eR6skuwCW+3u2wHcfbuZHRWWjwG2RpZrDsviyvdhZnMIzkKorq6mvr4+awCtra2x81pagufKxM1PwoMvvceF3vUt/7Sjjd0v/5nVL3Xvtrl05QaOeSf4VvbbZe90rm/A6nXPccw7w/nr0RU8+UpwcXJvWwd3P7KS3ccNZ9Hqd7Mu/5HRXd8H9rZ18JP/Ws7sEw/oFl+u/dVX89bt6XYmlK5z/cu9f4++xnXa6AqeyLKPCiUurnnr9nS7mawvde6rU4+u4JntydU5zi/XvUv6aleuOhf6b2z1S11naF1/75viV4iRxN9+IcTFVczjVynuk8h2Rc1zlO9b6D4XmAtQW1vrdXV1WTdUX19P3LzKzZUAsfOTMGpSivaHjQ53RgwbwpWfnk7NhCpGTUqxeNMzvBceWZ561bny0ycD8OIfnu5cf9jQIVz2yVOomVDFtgO38OR9a4FgJ0078QOMOnoUax7KvvyoSSn2PkTn8ultRNvJc+2vvmhoSvHU0me6xZOucz76Gte2A7fwRMY+qivAc6JyxdXQlOLJAtS5r7YduIVnEqxzNg1NKZ56uHd1LvTf2KhJKR7cvIy9bR3d/t7zVei4CiUurmIev5JMEjvMbHR4FjEa2BmWNwPjIsuNBbaF5XUZ5fUJxld0NROq2D36EHamdnPXl2d0/jHH3Rz0zp622G6j2W4aW/tKS7cHyUWXr5lQxYsHj+CNvwTNXekug0kevIrZJTKb6D4ykuuCGpV5TagUdU4bYkWsc8I30MWpmVDFXV+ewbJNbzDj2MN1Q10CkkwSi4HZwE3h7wci5Vea2QKCi9RvhonkIeC7kYvVZwJfSzC+khh1wFDaDzSOy/hjvmDaWBau2sredseBhau20hEZsSXzxqT0TWPps4+Fq7biHfHLAxx96AGdSaK/Yzj0RrpLpHtxukRm2/6IoUN4t60DrDjt1dF9muTNZLm2P2LoEPa0BY0/xajzoQd0HUaS7t6cTc2EKiWHBBWqC+zdwDPA+82s2cy+RJAczjCzF4EzwmmAJcAmoBG4FfgnAHffBdwArAx/rg/L9guZNwe1tXu3Nt7Mb2eZN421tXsvH6rWdf9A0t5+r62oXSIz1Uyo4lufPiE4myhSV9jtb77TVdNsN7ckrGZCFd/+9AlAcbr/NjSluP7BjUBw5pJ0V18pvkL1brrM3Ue7+zB3H+vut7n7G+5+urtPDn/vCpd1d7/C3Y9z95PcfVXkfW539+PDn18WIraBJNprKSruG2n0prGouAFe3np3L+mD9d725LvC3lzf2Pm6WF0iM0WbW4rRFXZuCbqBZipmnZdteoP3IjdoFqN5S4pLd1yXkZoJVVmf7xN3VpB59pE2ZfQhWZc/5IBhmBWnK2xDU4pnXuo6ESxWl8hM0Wc5Jd0ttGHzLp59peuBgqWsc7G6wo46oOvJvh2uLqiDkZJEmbnklO49UYb1MOxj8DiG3O+RNuqAoRwZGWc4ydHb7lrWVJRnNfUks1kuyTrf+seus4hS1/ljY7quEyRZ54fX7+x8neTzqaR0lCTKzKzp4/nu+Sdx8thDOXNKNQt6eBx4zYQq7rn8w5wxpZqTxx7Kd8/P/cjxI0eOSPybdUNTigfWdI1jXYoLuFHFeGJoQ1OKhzfs6JwudZ1PGzM08ZsJG5pSPNX4eud0Uk/2ldLSeBJlaFbGozd6UjOhils/V9urZUcdMJSLE34Wf+Y41qX6Rp2WbpabvyKoc/rRJ4WMKXMc61LX+fiqCi6YNpbfrArOIJKoc+Y41qWusyRDZxL7oVyP9eivpMex7qsTczz6pL+SHse6rz6U4PjXSY9jLeVDSWI/VDOhiotquv6hCzkWdClvrMolyXGRl216PdFxrPsqyTr/ZtXWsjpblOQoSeyn/mps92+Zu9/ZW5D3fSX1dufrJEdjy1fmiHX3rNiacyTAfLywo7Xz9f5Q52KOQCelpySxn8rshfKLJ1/u95CXDU0p7l7ZdfAo1qMwemOfEevcY0cCzEdDU4rFa7qaXfaHOkevOQHUve9InUUMYkoS+6kZxx7eNZIZwY1fT73Sv7OJ+Suaut1kXIxHf+Sj2+htxI8EmI9fPvVyt3vJ94c6b9zxVrfpIyLdqmXwUZLYT9VMqOL6mSd2jlrnwOOvtPe5OWL+8i3c29DV7bViiHH9zBPL6htmZ53DaQfuWdn3Jpj6LXt58NnIuBFlXOfoaIYLV23t89nE/GVNrG5q6ZyOu7tfBg8lif3YrOndR63rcPrUHNHQlOKbD6zr1h3yklPG5dWNt1gy69ze0bcmmIamFHdu7D6WQTnX+fQPHtk53dbHR7I0NKX4xgOdIxRjwMW148oqKUomL35NAAAQzElEQVThKUns54LusN2bI/Lt6RQMsdrV6FIxpLy/XV5Q0/8639vQHHl4YfnXufqQA7tNvx4ZRre3Fq7aMqDqLIWhJLGfq5lQxT98ZFLntJN/T6fnX+3eRv3lj0wq62+X/a1zQ1OKBSu7mqjKsZkpU+a1iUc37sirma2hKcVvVkaGoR0AdZbCUJIQRh04rFt/+rl/3NTrA8gdT73M4j93tctb+H7lLjPGfOr840df7Na75/QPHFWWzUxRNROquOSUaE8n+Mb9a3td5x889DwdkemBUGcpDCUJCZ4aGvmW2dHLA8j8ZU1c918bupWVW++eOJm9u3pb518/s5nHX3itW9lA6d2TeTbR4fDN+9f2eD3m9ic38XRGc9xAqbP0n5KE7NMDBoIDyLX3xR805y/fwtfvX9etbIgxYJogMnt3Qc91/t6SjXzjgfXdyiqMAdMun+1zbne4etGzsYniOw+u7xxUKG0g1Vn6L/EkYWabzWytma0xs1Vh2WFmttTMXgx/V4XlZmY/NrNGM3vWzKYlHZ8EZk0fz+wpw7sdQBz4+n1rmXPnqs6DSENTin+Yt5Kv37e22/pm8J3zcj+BttzMmj6e75x3Uq/qfOHNT/HzJzZ1W9+AG847aUAkxbRZ08dzxpTqbmWNO1u58JanuWlJVzJoaErxmZ89yS+e3Nxt2SE28Oos/VOsp8B+wt1fj0xfAzzq7jeZ2TXh9NXAOcDk8Gc6cHP4W4qgbvwwtlHV7ZHXAA9v2MHSDTt4/9GjeO7V3VnX/eQHqwdUgkibNX089c/vjK3z+MMOomnX2/usZ8DsKcMHZJ0v//hxPPrczm490tzhlic28fgLr9GB8/yrrfusZwy8LwLSf6VqbpoJzAtfzwPOi5TfGQ5xugyoNLPRpQhwf3X5x49jaJYhUR1iE8TQCuMrHz8u4ciSk6vO2RIEwOUfO5a68eV/gT6bmglV3DDzRLJUmY2v7s6aICCosxLE/sc84cHazexlIEXwP/dzd59rZi3uXhlZJuXuVWb2IHCTuz8Zlj8KXB0dBzssnwPMAaiurq5ZsGBB1m23trYycuTIrPOuWnMVAD+a+qN+1jA/U1+/ivb2dtZW/2dRt5veNsCaI7LXOb2/GlPtLNn0Hqtf68i6XNS0oyo4d9Iwjq+q6HHZvsr1ORZKb+tswDkTh3LxB0YUJa6+6k1sjal2Fj6/hxdach8DonUuRlylMNDiKsTx6xOf+ESDu/c4EE0xmptOc/dtZnYUsNTMnsuxbJbvNuzzF+zuc4G5ALW1tV5XV5f1zerr64mbV7k5yFFx8xPzSCUtLS3F3264bYivc3p/1QFfJrg4/Y3713a7gSptSBGvQeT6HAuljp7rfOrEKq4+54Od7fHFiKuvehNbHUGdb1qykZ8/sWnffzT2rXMx4iqFgRZXMY9fiScJd98W/t5pZvcBpwI7zGy0u28Pm5PSA+U2A+Miq48FtiElMWv6eN5/9CgWrW7udofukaNGcMG0wTl+QGadW95+jz1tHVxySn6jBQ4k15z7Qc444ej9qs7Se4kmCTM7GBji7rvD12cC1wOLgdnATeHvB8JVFgNXmtkCggvWb7r79n3fWYqlZkLVoEwGuajOIl2SPpOoBu4zs/S25rv7H8xsJbDQzL4EbAEuCpdfApwLNAJvA19IOD4REckh0STh7puAk7OUvwGcnqXcgSuSjElERHpPd1yLiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISq+yShJmdbWbPm1mjmV1T6nhERPZnZZUkzKwC+ClwDjAFuMzMppQ2KhGR/VfSw5fm61SgMRzRjnCs65nAhpJGVUg7H6cS4JG64m87tQaqphZ/uyIyYJVbkhgDbI1MNwPTMxcysznAHIDq6mrq6+uzvllra2vsvJaWFoDY+Umpy9h+UdlEduypYXsf9lcpKa78lWtsiis/cXEV8/hVbknCspT5PgXuc4G5ALW1tV5XV5f1zerr64mbV7m5EiB2fnI8Z1xJqwTeHzOvlHHlorjyV66xKa78xMVVzONXWV2TIDhzGBeZHgtsK1EsIiL7vXJLEiuByWY2ycyGA5cCi0sck4jIfqusmpvcvc3MrgQeAiqA2919fYnDEhHZb5VVkgBw9yXAklLHISIi5dfcJCIiZURJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkViJJQkzu87MXjGzNeHPuZF5XzOzRjN73szOipSfHZY1mtk1ScUmIiK9k/TIdD9093+LFpjZFIKxq08AjgEeMbP3hbN/CpwBNAMrzWyxu29IOEYREYlRiuFLZwIL3H0P8LKZNQKnhvMa3X0TgJktCJdVkhARKZGkk8SVZvY5YBXwVXdPAWOAZZFlmsMygK0Z5dOzvamZzQHmAFRXV1NfX591462trbHzWlpaAGLnJylXXKWkuPJTrnFB+camuPITF1cxj1/9ShJm9ghwdJZZ1wI3AzcAHv7+AfBFwLIs72S/PuLZtuvuc4G5ALW1tV5XV5c1vvr6euLmVW6uBIidn6RccZWS4spPucYF5Rub4spPXFzFPH71K0m4+yd7s5yZ3Qo8GE42A+Mis8cC28LXceUiIlICSfZuGh2ZPB9YF75eDFxqZiPMbBIwGVgBrAQmm9kkMxtOcHF7cVLxiYhIz5K8JvGvZjaVoMloM3A5gLuvN7OFBBek24Ar3L0dwMyuBB4CKoDb3X19gvGJiEgPEksS7v7ZHPNuBG7MUr4EWJJUTCIikh/dcS0iIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYnVryRhZheZ2Xoz6zCz2ox5XzOzRjN73szOipSfHZY1mtk1kfJJZrbczF40s3vCIUxFRKSE+nsmsQ64AHgiWmhmUwjGqD4BOBv4mZlVmFkF8FPgHGAKcFm4LMD3gR+6+2QgBXypn7GJiEg/9Wv4UnffCGBmmbNmAgvcfQ/wspk1AqeG8xrdfVO43gJgppltBP4GmBUuMw+4Dri5P/Hl8njT4wDU3VGX1CZitbS0ULm5sujb7Yniyk+5xgXlG5viyk9cXGteXcPUo6cWJYakxrgeAyyLTDeHZQBbM8qnA4cDLe7elmX5fZjZHGAOQHV1NfX19VmXa21tjZ2X1tLSknN+Etrb20uy3Z4orvyUa1xQvrEprvzExTXxgInUDK/p8fhWCD0mCTN7BDg6y6xr3f2BuNWylDnZm7c8x/JZuftcYC5AbW2t19XVZV2uvr6euHleF/v2icsVVykprvyUa1xQvrEprvyUQ1w9Jgl3/2Qf3rcZGBeZHgtsC19nK38dqDSzoeHZRHR5EREpkaS6wC4GLjWzEWY2CZgMrABWApPDnkzDCS5uL3Z3Bx4DLgzXnw3EnaWIiEiR9LcL7Plm1gz8NfA7M3sIwN3XAwuBDcAfgCvcvT08S7gSeAjYCCwMlwW4GviX8CL34cBt/YlNRET6r7+9m+4D7ouZdyNwY5byJcCSLOWb6OoBJSIiZUB3XIuISCwlCRERiaUkISIisZQkREQklgW9TwcuM3sNaIqZfQTBPRjlRnHlR3Hlr1xjU1z5STKuCe5+ZE8LDfgkkYuZrXL32p6XLC7FlR/Flb9yjU1x5acc4lJzk4iIxFKSEBGRWIM9ScwtdQAxFFd+FFf+yjU2xZWfksc1qK9JiIhI/wz2MwkREekHJQkREYk1oJOEmf0/M3vOzJ41s/vMrDIsn2hm75jZmvDnlsg6NWa21swazezHFo69amaHmdlSM3sx/F2VUMxnm9nz4favSWIbGdsbZ2aPmdlGM1tvZv8cll9nZq9E9tG5kXW+Fsb3vJmdlVTsZrY5/CzWmNmqsCzr52CBH4fbftbMpkXeZ3a4/ItmNrufMb0/sk/WmNlbZnZVKfaXmd1uZjvNbF2krGD7J+5/oY9xlfx/MSaugn1uFgxxsDyM6x4Lhjvoa1z3RGLabGZrir2/es3dB+wPcCYwNHz9feD74euJwLqYdVYQPNrcgN8D54Tl/wpcE76+Jv1eBY63AngJOBYYDvwZmJLwPhoNTAtfjwJeAKYQjCH+v7MsPyWMawQwKYy3IonYgc3AERllWT8H4Nzw8zJgBrA8LD8M2BT+rgpfVxXw83oVmFCK/QV8DJgW/Vsu5P6J+1/oY1wl/1+MiatgnxvB8AeXhq9vAf6xr3FlzP8B8K1i76/e/gzoMwl3f9i7xsVeRjCiXSwzGw0c4u7PeLBH7wTOC2fPBOaFr+dFygvpVKDR3Te5+3vAgnC7iXH37e6+Ony9m2Acj9jxw8N4Frj7Hnd/GWgM4y5W7HGfw0zgTg8sIxjJcDRwFrDU3Xe5ewpYCpxdoFhOB15y97g7+tNxJbK/3P0JYFeW7fV7//Twv5B3XOXwvxizv+Lk9bmF39r/BvhtIeMK3/di4O5c71HKY9eAThIZvkiQXdMmmdmfzOxxM/toWDaGYGjVtGa6DpjV7r4dggMrcFQCMY4BtsZsP3FmNhH4ELA8LLoybB64PXKKGhdjErE78LCZNZjZnLAs7nMoZlxpl9L9n7fU+wsKt39y/S/0V7n9LxbiczscaIkkwkLtr48CO9z9xUhZqfdXN2WfJMzsETNbl+VnZmSZa4E24K6waDsw3t0/BPwLMN/MDiE4TctUzD7AJdu+mY0EFgFXuftbwM3AccBUgv31gx5iTCL209x9GnAOcIWZfSzHssWMi7C9+W+B34RF5bC/csk3jqT2W7n9Lxbqc0sq3svo/kWk1PtrH/0ama4Y3P2TueZbcCHuU8Dp4WkY7r4H2BO+bjCzl4D3EWTf6GnwWGBb+HqHmY129+3hqd3OwtYEwu2Pi9l+YsxsGEGCuMvd7wVw9x2R+bcCD/YixoLG7u7bwt87zew+glP9uM8hLq5moC6jvL4/cYXOAVan91M57K9QofZPrv+FPinH/8UCfm6vEzThDQ3PJgqxv4YCFwA1kXjL7thV9mcSuZjZ2QRjY/+tu78dKT/SzCrC18cCk4FN4anYbjObEbYFfg54IFxtMZDu+TE7Ul5IK4HJYS+J4QTNGYsT2E6nsJ63ARvd/d8j5aMji50PpHteLAYuNbMRZjaJYN+tKHTsZnawmY1Kvya48LmO+M9hMfA5C8wA3gw/z4eAM82sKmxKODMs669u3/BKvb8iCrJ/evhfyFu5/i8W6nMLk95jwIWFiCv0SeA5d+9sRir1/sqqkFfBi/1DcLFpK7Am/LklLP8MsJ6gZ8Jq4NORdWoJ/lBeAn5C113nhwOPAi+Gvw9LKOZzCXoYvQRcW4R99BGC09JnI/vpXOBXwNqwfDEwOrLOtWF8zxPp8VLI2Al6j/w5/Fmffr+4z4HgdPun4bbXArWR9/pi+LfQCHyhAPvsIOAN4NBIWdH3F0GS2g7sJfgm+aVC7p+4/4U+xlXy/8WYuAr2uYV/syvCuv4GGNHXuMLyO4CvZCxbdscuPZZDRERiDejmJhERSZaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYn1/wFIbEZYMCCTCwAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x10c63afd0>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[1000])\n",
+    "output_file = pset.ParticleFile(name='NestedFieldParticle.nc', outputdt=50)\n",
+    "pset.execute(AdvectionRK4, runtime=14000, dt=10, output_file=output_file)\n",
+    "plt = plotTrajectoriesFile('NestedFieldParticle.nc', show_plt=False)\n",
+    "plt.plot([0,2e3,2e3,0,0],[0,0,2e3,2e3,0], c='orange')\n",
+    "plt.plot([-2e3,18e3,18e3,-2e3,-2e3],[-1e3,-1e3,3e3,3e3,-1e3], c='green');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As we observe, there is a change of dynamic at lon=2000, which corresponds to the change of grid.\n",
+    "\n",
+    "The analytical solution to the problem:\n",
+    "\\begin{align}\n",
+    "dx/dt &= 1;\\\\\n",
+    "dy/dt &= \\cos(x \\pi/400);\\\\\n",
+    "\\text{with  } x(0) &= 0, y(0) = 1000\n",
+    "\\end{align}\n",
+    "\n",
+    "is\n",
+    "\\begin{align}\n",
+    "x(t) &= t;\\\\\n",
+    "y(t) &= 1000 + 400/\\pi \\sin(t \\pi / 400)\n",
+    "\\end{align}\n",
+    "which is captured by the High Resolution field (orange area) but not the Low Resolution one (green area)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Keep track of the field interpolated\n",
+    "\n",
+    "For different reasons, you may want to keep track of the field you have interpolated. You can do that easily by creating another field that share the grid with original fields.\n",
+    "Watch out that this operation has a cost of a full interpolation operation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "F1 = Field('F1', np.ones((U1.grid.ydim, U1.grid.xdim), dtype=np.float32), grid=U1.grid)\n",
+    "F2 = Field('F2', 2*np.ones((U2.grid.ydim, U2.grid.xdim), dtype=np.float32), grid=U2.grid)\n",
+    "F = NestedField('F', [F1, F2])\n",
+    "fieldset.add_field(F)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: Compiled SampleParticleSampleNestedFieldIndex ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/57e7f8f8b9c55074a13384b46fceb0eb.so\n",
+      "WARNING: dt or runtime are zero, or endtime is equal to Particle.time. The kernels will be executed once, without incrementing time\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Particle (1000, 500) interpolates Field #1\n",
+      "Particle (10000, 500) interpolates Field #2\n"
+     ]
+    }
+   ],
+   "source": [
+    "from parcels import Variable\n",
+    "\n",
+    "def SampleNestedFieldIndex(particle, fieldset, time, dt):\n",
+    "    particle.f = fieldset.F[time, particle.lon, particle.lat, particle.depth]\n",
+    "\n",
+    "class SampleParticle(JITParticle):\n",
+    "    f = Variable('f', dtype=np.int32)\n",
+    "    \n",
+    "pset = ParticleSet(fieldset, pclass= SampleParticle, lon=[1000], lat=[500])\n",
+    "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n",
+    "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n",
+    "\n",
+    "pset[0].lon = 10000\n",
+    "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n",
+    "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "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.15"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/parcels/field.py b/parcels/field.py
index b122ed94..5969b279 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -14,7 +14,7 @@ from .grid import (RectilinearZGrid, RectilinearSGrid, CurvilinearZGrid,
                    CurvilinearSGrid, CGrid, GridCode)
 
 
-__all__ = ['Field', 'VectorField', 'SummedField', 'SummedVectorField']
+__all__ = ['Field', 'VectorField', 'SummedField', 'SummedVectorField', 'NestedField']
 
 
 class Field(object):
@@ -1156,6 +1156,48 @@ class SummedVectorField(list):
             return self.eval(*key)
 
 
+class NestedField(list):
+    """Class NestedField is a list of Fields from which the first one to be not declared out-of-boundaries
+    at particle position is interpolated. This induces that the order of the fields in the list matters.
+    Each one it its turn, a field is interpolated: if the interpolation succeeds or if an error other
+    than `ErrorOutOfBounds` is thrown, the function is stopped. Otherwise, next field is interpolated.
+    NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries.
+    NestedField is composed of either Fields or VectorFields.
+
+    :param name: Name of the Nested field
+    :param U: List of fields (order matters). U can be a scalar Field, a VectorField, or the zonal component of the VectorField
+    :param V: List of fields defining the meridional component (default: None)
+    :param W: List of fields defining the vertical component (default: None)
+    """
+
+    def __init__(self, name, U, V=None, W=None):
+        if V is None:
+            for Ui in U:
+                self.append(Ui)
+        elif W is None:
+            for (i, Ui, Vi) in zip(range(len(U)), U, V):
+                self.append(VectorField(name+'_%d' % i, Ui, Vi))
+        else:
+            for (i, Ui, Vi, Wi) in zip(range(len(U)), U, V, W):
+                self.append(VectorField(name+'_%d' % i, Ui, Vi, Wi))
+        self.name = name
+
+    def __getitem__(self, key):
+        if isinstance(key, int):
+            return list.__getitem__(self, key)
+        else:
+            for iField in range(len(self)):
+                try:
+                    val = list.__getitem__(self, iField).eval(*key)
+                    break
+                except FieldSamplingError:
+                    if iField == len(self)-1:
+                        raise
+                    else:
+                        pass
+            return val
+
+
 class NetcdfFileBuffer(object):
     """ Class that encapsulates and manages deferred access to file data. """
 
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index 6fee225a..040d4f66 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -1,4 +1,4 @@
-from parcels.field import Field, VectorField, SummedField, SummedVectorField
+from parcels.field import Field, VectorField, SummedField, SummedVectorField, NestedField
 from parcels.gridset import GridSet
 from parcels.grid import RectilinearZGrid
 from parcels.tools.loggers import logger
@@ -98,6 +98,12 @@ class FieldSet(object):
             for fld in field:
                 self.gridset.add_grid(fld)
                 fld.fieldset = self
+        elif isinstance(field, NestedField):
+            setattr(self, name, field)
+            for fld in field:
+                setattr(self, fld.name, fld)
+                self.gridset.add_grid(fld)
+                fld.fieldset = self
         elif isinstance(field, list):
             raise NotImplementedError('FieldLists have been replaced by SummedFields. Use the + operator instead of []')
         else:
@@ -112,6 +118,10 @@ class FieldSet(object):
         """
         setattr(self, vfield.name, vfield)
         vfield.fieldset = self
+        if isinstance(vfield, NestedField):
+            for f in vfield:
+                setattr(self, f.name, f)
+                f.fieldset = self
 
     def check_complete(self):
         assert self.U, 'FieldSet does not have a Field named "U"'
@@ -132,11 +142,15 @@ class FieldSet(object):
         if not hasattr(self, 'UV'):
             if isinstance(self.U, SummedField):
                 self.add_vector_field(SummedVectorField('UV', self.U, self.V))
+            elif isinstance(self.U, NestedField):
+                self.add_vector_field(NestedField('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'):
             if isinstance(self.U, SummedField):
                 self.add_vector_field(SummedVectorField('UVW', self.U, self.V, self.W))
+            elif isinstance(self.U, NestedField):
+                self.add_vector_field(NestedField('UVW', self.U, self.V, self.W))
             else:
                 self.add_vector_field(VectorField('UVW', self.U, self.V, self.W))
 


Test patch:
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index 48eae470..0f7ce300 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -1,5 +1,5 @@
-from parcels import (FieldSet, Field, ParticleSet, ScipyParticle, JITParticle, Geographic,
-                     AdvectionRK4, AdvectionRK4_3D, Variable)
+from parcels import (FieldSet, Field, NestedField, ParticleSet, ScipyParticle, JITParticle, Geographic,
+                     AdvectionRK4, AdvectionRK4_3D, Variable, ErrorCode)
 import numpy as np
 import pytest
 from math import cos, pi
@@ -453,3 +453,53 @@ def test_summedfields(mode, with_W, k_sample_p, mesh):
     assert np.isclose(pset[0].p, 60)
     assert np.isclose(pset[0].lon*conv, 0.6, atol=1e-3)
     assert np.isclose(pset[0].lat, 0.9)
+
+
+@pytest.mark.parametrize('mode', ['jit', 'scipy'])
+def test_nestedfields(mode, k_sample_p):
+    xdim = 10
+    ydim = 20
+
+    U1 = Field('U1', 0.1*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 1., xdim, dtype=np.float32),
+               lat=np.linspace(0., 1., ydim, dtype=np.float32))
+    V1 = Field('V1', 0.2*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 1., xdim, dtype=np.float32),
+               lat=np.linspace(0., 1., ydim, dtype=np.float32))
+    U2 = Field('U2', 0.3*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 2., xdim, dtype=np.float32),
+               lat=np.linspace(0., 2., ydim, dtype=np.float32))
+    V2 = Field('V2', 0.4*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 2., xdim, dtype=np.float32),
+               lat=np.linspace(0., 2., ydim, dtype=np.float32))
+    U = NestedField('U', [U1, U2])
+    V = NestedField('V', [V1, V2])
+    fieldset = FieldSet(U, V)
+
+    P1 = Field('P1', 0.1*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 1., xdim, dtype=np.float32),
+               lat=np.linspace(0., 1., ydim, dtype=np.float32))
+    P2 = Field('P2', 0.2*np.ones((ydim, xdim), dtype=np.float32),
+               lon=np.linspace(0., 2., xdim, dtype=np.float32),
+               lat=np.linspace(0., 2., ydim, dtype=np.float32))
+    P = NestedField('P', [P1, P2])
+    fieldset.add_field(P)
+
+    def Recover(particle, fieldset, time, dt):
+        particle.lon = -1
+        particle.lat = -1
+        particle.p = 999
+        particle.time = particle.time + dt
+
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0], lat=[.3])
+    pset.execute(AdvectionRK4+pset.Kernel(k_sample_p), runtime=1, dt=1)
+    assert np.isclose(pset[0].lat, .5)
+    assert np.isclose(pset[0].p, .1)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0], lat=[1.3])
+    pset.execute(AdvectionRK4+pset.Kernel(k_sample_p), runtime=1, dt=1)
+    assert np.isclose(pset[0].lat, 1.7)
+    assert np.isclose(pset[0].p, .2)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=[0], lat=[2.3])
+    pset.execute(AdvectionRK4+pset.Kernel(k_sample_p), runtime=1, dt=1, recovery={ErrorCode.ErrorOutOfBounds: Recover})
+    assert np.isclose(pset[0].lat, -1)
+    assert np.isclose(pset[0].p, 999)