{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial on how to combine different Fields for advection into a `SummedField` 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 `SummedField` objects.\n", "\n", "This tutorial shows how to use these `SummedField` with a very idealised example. We start by importing the relevant modules." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "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_f7p00000gr/T/parcels-504/95f68bdd3e2c61236f4287b88f438342.so\n", "100% (10.0 of 10.0) |####################| Elapsed Time: 0:00:00 Time: 0:00:00\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFctJREFUeJzt3XuU3WV97/H3B4LcFQSJQICgpt5YYiGVVEVD0QqIoq6ieLAiXqKntKDLtl7as3CdVayu9tTWYw8HitdKgxbUspCqiEaPy8IxQfAE4yWIkAgIQlACCAn5nj/2b84Zh8nMzmT27Mwz79das2bv336e3/P97iSf/cuz9ySpKiRJ7dpp2AVIkgbLoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBrx1WkvcmuaiPcZ9I8lczUdNW1v/3JGcMa/1RdfT1fGnuiZ+j11Ql+SkwH3gEuB+4EviTqto4hXMtBT5dVQumMPcTwPqq+sspzC1gUVWt3da50ynJ+4CnVNXrhlmH2uQVvbbXy6pqL+Ao4HeAqYTtvGmvaprsyLWNNlvq1HAY9JoWVfUz4N+BIwCSnJlkTZL7kvwkyVtHxiZZmmR9kncluQNY3s09KMnG7uugJO9L8ulR856f5NtJ7k2yLskbxqslyclJru/GfTvJs7Yy7pvdzRu6NV8zTm0fT7JvkiuS3JVkQ3d7wajzrEjy5lH339j1viHJl5McNuqxZya5Ksk9SX7ebbecALwXeE1Xxw3d2IOSXN6NXZvkLaPO874klyb5dJJfAW8Y5/laMur5uqH7W9PIY2/ofl3uS3JzktO3/qur2c6g17RIcghwEvDd7tCdwMnAY4EzgQ8lOWrUlCcCjwcOA14PnAjcVlV7dV+3jTn/ofReDP478ATg2cD149RxFPAx4K3AfsAFwOVJdh07tqpe0N08slvzM+PUtozen5OPd/cPBR4EPrKV5+EV9EL7VV2d/4veCxlJ9ga+CnwJOAh4CnB1VX0JeD/wma6OI7vTLQfWd2P/AHh/kuNHLXcKcCmwD3DxmDoOBr4I/FXXy58ClyV5QpI9gQ8DJ1bV3sBzx3su1Q6DXtvrC0nuBb4FfINeYFFVX6yqm6rnG8BXgGNHzdsCnFtVD1XVg32sczrw1apaXlWbquruqhovnN4CXFBV11bVI1X1SeAhYMk29PQbtXVrXVZVD1TVfcB5wAu3MvetwF9X1Zqq2kzv+Xh2d1V/MnBHVf23qvp1Vd1XVdeOd5LuhfP5wLu6sdcDFwF/OGrYf1TVF6pqyzjP4euAK6vqyu7xq4CV9F6MR3o8IsnuVXV7Vd24Dc+PZhmDXtvrFVW1T1UdVlV/NBI4SU5Mck237XAvvYDZf9S8u6rq19uwziHATX2MOwx4Z7ddcW+39iH0ror79Ru1JdkjyQVJbum2Sb4J7JNk562s/w+j1r4HCHDwNvRAV+893QvLiFu684xYN8H8w4BTxzwPzwcOrKr7gdcAbwNuT/LFJE/rsy7NQga9pl23TXIZ8LfA/Krah94ncjJq2NiPe0328a91wJP7WH4dcF734jPytUdVLe+z/PFqeSfwVOCYqnosMLLlEx5tHfDWMevvXlXfnqSHsWveBjy+2+4ZcSjwswnmjK3jn8fUsWdVfQCgqr5cVS8GDgR+APzTBOfSLGfQaxAeA+wK3AVsTnIi8PuTzPk5sF+Sx23l8YuBFyV5dZJ5SfZL8uxxxv0T8LYkx6RnzyQvHROYY9d90iS17U1vX/7eJI8Hzp1g7P8E3pPkmQBJHpfk1O6xK4AnJnl7kl2T7J3kmFF1LEyyE0BVrQO+Dfx1kt26N5TfxJi9+Al8GnhZkpck2bk7x9IkC5LMT/Lybq/+IWAjvY/IqlEGvaZdt91wNvBZYAPwn4DLJ5nzA3pvPv6k22o4aMzjt9Lb/nknve2Q64EjxznPSnr79B/p1l4LvGGCpd8HfLJb89VbGfP3wO7AL4Br6L2ZurU+Pg98ELik2+ZZTe+N5pHn5cXAy4A7gB8Dx3VT/7X7fneS67rbrwUW0ru6/zy99w2umqCX0XWso/dm7XvpveCuA/6M3p/5neg9j7fRey5fCPxRP+fV7OQPTEnbqfuY5kVV9alh1yKNxyt6aTsk2YPe1s/Nw65F2hqDXpqiJAfQ24L5Br2Pl0o7JLduJKlxXtFLUuN2iH8Iaf/996+FCxdOae7999/PnnvuOb0F7eDseW6w57lhe3petWrVL6rqCZON2yGCfuHChaxcuXJKc1esWMHSpUunt6AdnD3PDfY8N2xPz0lu6WecWzeS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb1FfRJzkmyOsmNSd4+5rE/TVJJ9u/uJ8mHk6xN8r0kRw2icElSfyYN+iRHAG8BngMcCZycZFH32CHAi4FbR005EVjUfS0Dzp/mmiVJ26CfK/qnA9dU1QNVtRn4BvDK7rEPAX8O1KjxpwCfqp5rgH2SHDidRUuS+tfP/xm7GjgvyX7Ag8BJwMokLwd+VlU3JBk9/mBg3aj767tjt48elGQZvSt+5s+fz4oVK6bUwMaNG6c8d7ay57nBnueGmeh50qCvqjVJPghcBWwEbgA2A38B/P44UzLOsXrUgaoLgQsBFi9eXFP9z3H9z4TnBnueG+x5MPp6M7aqPlpVR1XVC4B7gJ8ChwM3JPkpsAC4LskT6V3BHzJq+gLgtuksWpLUv34/dXNA9/1Q4FX09uAPqKqFVbWQXrgfVVV3AJcDr+8+fbME+GVV3b61c0uSBqufPXqAy7o9+k3AWVW1YYKxV9Lbx18LPACcuX0lSpK2R19BX1XHTvL4wlG3Czhr+8qSJE0XfzJWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqXF9Bn+ScJKuT3Jjk7d2xv0nygyTfS/L5JPuMGv+eJGuT/DDJSwZVvCRpcpMGfZIjgLcAzwGOBE5Osgi4Cjiiqp4F/Ah4Tzf+GcBpwDOBE4D/kWTnwZQvSZpMP1f0TweuqaoHqmoz8A3glVX1le4+wDXAgu72KcAlVfVQVd0MrKX3IiFJGoJ5fYxZDZyXZD/gQeAkYOWYMW8EPtPdPphe8I9Y3x37DUmWAcsA5s+fz4oVK7ap8BEbN26c8tzZyp7nBnueG2ai50mDvqrWJPkgva2ajcANwMiVPEn+ort/8cih8U4zznkvBC4EWLx4cS1dunRbawdgxYoVTHXubGXPc4M9zw0z0XNfb8ZW1Uer6qiqegFwD/BjgCRnACcDp1fVSJivBw4ZNX0BcNv0lSxJ2hb9furmgO77ocCrgOVJTgDeBby8qh4YNfxy4LQkuyY5HFgE/O/pLVuS1K9+9ugBLuv26DcBZ1XVhiQfAXYFrkoCvTds31ZVNyb5LPB9els6Z1XVI4MoXpI0ub6CvqqOHefYUyYYfx5w3nbUJUmaJv5krCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Li+gj7JOUlWJ7kxydu7Y49PclWSH3ff9+2OJ8mHk6xN8r0kRw2yAUnSxCYN+iRHAG8BngMcCZycZBHwbuDqqloEXN3dBzgRWNR9LQPOH0DdAKy6ZQNX3PQwq27ZMKglJlz7H7++dsbXtueZNdd6Hla/I2vb82CkqiYekJwKvKSq3tzd/y/AQ8CbgKVVdXuSA4EVVfXUJBd0t5d34384Mm5rayxevLhWrly5TYWvumUDr77gP3hkS7FT4GlP3Ju9d9tlm84xVff9ehM/uOM+thQzuvaw1h3m2vbc/rrDXHtH6Xm3XXbi4jcv4ejD9t2mcyRZVVWLJxs3r49zrQbOS7If8CBwErASmD8S3l3YH9CNPxhYN2r++u7YbwR9kmX0rviZP38+K1as6KOU/++Kmx7mkS29F6ktBXdu2Mgju2ebzjFVdz9YdEvP6NrDWneYa9tz++sOc+0dpeeHN21h+Ve/w31PfsxA1po06KtqTZIPAlcBG4EbgM0TTBnvWXrUXxuq6kLgQuhd0S9durSfev+fvQ/fwBU/vYaHN23hMbvsxAVnbvur4VStumUDp190DZs2b2GXeTO39si69mzPg1x3pvsdvfZc7vm1L/qdga096dbNoyYk76d3lX4OQ9y6gd4Ttfyr3xnoEzTR2tf85G6WPGm/GV3bnu150OsOo9+Rte1520zn1g1JDqiqO5McCrwK+F3gcOAM4APd93/rhl8O/HGSS4BjgF9OFPLb4+jD9uW+Jz9mxn9xRtYe1rr2PLNrz6Weh9XvyNr2PBh9BT1wWbdHvwk4q6o2JPkA8NkkbwJuBU7txl5Jbx9/LfAAcOY01yxJ2gZ9BX1VHTvOsbuB48c5XsBZ21+aJGk6+JOxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4/oK+iTvSHJjktVJlifZLcnxSa5Lcn2SbyV5Sjd21ySfSbI2ybVJFg6yAUnSxCYN+iQHA2cDi6vqCGBn4DTgfOD0qno28C/AX3ZT3gRsqKqnAB8CPjiIwiVJ/el362YesHuSecAewG1AAY/tHn9cdwzgFOCT3e1LgeOTZHrKlSRtq1TV5IOSc4DzgAeBr1TV6UmOBb7QHfsVsKSqfpVkNXBCVa3v5t4EHFNVvxhzzmXAMoD58+cffckll0ypgY0bN7LXXntNae5sZc9zgz3PDdvT83HHHbeqqhZPOrCqJvwC9gW+BjwB2IVeuL8O+By9AAf4M+Ci7vaNwIJR828C9ptojaOPPrqm6utf//qU585W9jw32PPcsD09Aytrkgyvqr62bl4E3FxVd1XVpi7gnwccWVXXdmM+Azy3u70eOASg2+p5HHBPH+tIkgagn6C/FViSZI9ur/144PvA45L8VjfmxcCa7vblwBnd7T8Avta98kiShmDeZAOq6toklwLXAZuB7wIX0rtyvyzJFmAD8MZuykeBf06ylt6V/GmDKFyS1J9Jgx6gqs4Fzh1z+PPd19ixvwZO3f7SJEnTwZ+MlaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGtdX0Cd5R5Ibk6xOsjzJbuk5L8mPkqxJcnY3Nkk+nGRtku8lOWqwLUiSJjJvsgFJDgbOBp5RVQ8m+SxwGhDgEOBpVbUlyQHdlBOBRd3XMcD53XdJ0hD0u3UzD9g9yTxgD+A24D8D/7WqtgBU1Z3d2FOAT1XPNcA+SQ6c5rolSX1KVU0+KDkHOA94EPhKVZ2e5G7g74BXAncBZ1fVj5NcAXygqr7Vzb0aeFdVrRxzzmXAMoD58+cffckll0ypgY0bN7LXXntNae5sZc9zgz3PDdvT83HHHbeqqhZPNq6frZt96V2lHw7cC/xrktcBuwK/rqrFSV4FfAw4lt6WzliPejWpqguBCwEWL15cS5cunayUca1YsYKpzp2t7HlusOe5YSZ67mfr5kXAzVV1V1VtAj4HPBdYD1zWjfk88Kzu9np6e/cjFtDb6pEkDUE/QX8rsCTJHkkCHA+sAb4A/F435oXAj7rblwOv7z59swT4ZVXdPs11S5L6NOnWTVVdm+RS4DpgM/BdelsuuwMXJ3kHsBF4czflSuAkYC3wAHDmAOqWJPVp0qAHqKpzgXPHHH4IeOk4Yws4a/tLkyRNB38yVpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb19Z+DD7yI5C7glilO3x/4xTSWMxvY89xgz3PD9vR8WFU9YbJBO0TQb48kK/v5X9BbYs9zgz3PDTPRs1s3ktQ4g16SGtdC0F847AKGwJ7nBnueGwbe86zfo5ckTayFK3pJ0gQMeklq3KwO+iQnJPlhkrVJ3j3segYtySFJvp5kTZIbk5wz7JpmQpKdk3w3yRXDrmWmJNknyaVJftD9ev/usGsapCTv6H5Pr06yPMluw65pEJJ8LMmdSVaPOvb4JFcl+XH3fd/pXnfWBn2SnYF/BE4EngG8NskzhlvVwG0G3llVTweWAGfNgZ4BzgHWDLuIGfYPwJeq6mnAkTTcf5KDgbOBxVV1BLAzcNpwqxqYTwAnjDn2buDqqloEXN3dn1azNuiB5wBrq+onVfUwcAlwypBrGqiqur2qrutu30fvD//Bw61qsJIsAF4KXDTsWmZKkscCLwA+ClBVD1fVvcOtauDmAbsnmQfsAdw25HoGoqq+Cdwz5vApwCe7258EXjHd687moD8YWDfq/noaD73RkiwEfhu4driVDNzfA38ObBl2ITPoScBdwMe7LauLkuw57KIGpap+BvwtcCtwO/DLqvrKcKuaUfOr6nboXcwBB0z3ArM56DPOsTnxWdEkewGXAW+vql8Nu55BSXIycGdVrRp2LTNsHnAUcH5V/TZwPwP46/yOotuTPgU4HDgI2DPJ64ZbVVtmc9CvBw4ZdX8Bjf51b7Qku9AL+Yur6nPDrmfAnge8PMlP6W3N/V6STw+3pBmxHlhfVSN/W7uUXvC36kXAzVV1V1VtAj4HPHfINc2knyc5EKD7fud0LzCbg/47wKIkhyd5DL03by4fck0DlST09m3XVNXfDbueQauq91TVgqpaSO/X92tV1fyVXlXdAaxL8tTu0PHA94dY0qDdCixJskf3e/x4Gn7zeRyXA2d0t88A/m26F5g33SecKVW1OckfA1+m9y79x6rqxiGXNWjPA/4Q+D9Jru+OvbeqrhxiTRqMPwEu7i5ifgKcOeR6Bqaqrk1yKXAdvU+WfZdG/ymEJMuBpcD+SdYD5wIfAD6b5E30XvROnfZ1/ScQJKlts3nrRpLUB4Nekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNe7/AuHu7lXMa6xGAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "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 `SummedFields`. We can simply define a new `FieldSet` with a summation 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_f7p00000gr/T/parcels-504/8ba0de8db83dec9865083d5c33acfdb2.so\n", "100% (10.0 of 10.0) |####################| Elapsed Time: 0:00:00 Time: 0:00:00\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEICAYAAABBBrPDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FdX9//HXOwtrWMIWZRFEEERANAgBZF/cUKm477SItVgQWkv1+/uWtmq1/baAW1VcURFUQAsICgoSQEDZ932RTQUEMYisn98fd2JjzHITEu5N8nk+HvPI3JlzznzmQu4nc2buOTIznHPOuezERDoA55xz0c0ThXPOuRx5onDOOZcjTxTOOedy5InCOedcjjxROOecy5EnCldsSXpI0othlHtV0iOnI6Zsjj9V0p2ROn6GOMJ6v1zJI/8ehYsUSVuBJOAEcAiYAvzWzNLy0VYn4A0zq52Puq8CO8zs/+WjrgENzWxjXusWJEl/BhqY2W2RjMMVT35F4SLtKjNLAC4CLgby82EdV+BRFZBoji2johKniwxPFC4qmNlOYCrQFEBSH0lrJH0nabOke9LLSuokaYekIZK+BMYEdWtKSguWmpL+LOmNDPUukfSppAOStku6K6tYJPWUtDQo96mk5tmUSw1WlwXHvDGL2F6RlChpsqQ9kvYH67UztPOJpL4ZXv8yOPf9kj6UVDfDvvMlTZf0jaSvgu6iy4CHgBuDOJYFZWtKmhiU3Sjp7gzt/FnSOElvSDoI3JXF+5WS4f1aFly1pe+7K/h3+U7SFkm3Zv+v64o6TxQuKkiqA1wBLAk2fQ30BCoCfYDhki7KUOUMoApQF7gDuBzYZWYJwbIrU/tnEUomTwHVgRbA0iziuAh4GbgHqAo8D0yUVDpzWTPrEKxeEBzzrSxi60fo9+yV4PVZwGHg6Wzeh16EPvSvDeKcTSgRIqkC8BHwAVATaAB8bGYfAH8D3griuCBobgywIyh7HfA3SV0zHO4aYBxQGRidKY5awPvAI8G5/B4YL6m6pPLAk8DlZlYBaJvVe+mKD08ULtLek3QAmAPMIvSBh5m9b2abLGQWMA1on6HeSWComR0xs8NhHOdW4CMzG2Nmx8xsn5ll9eF2N/C8mS0wsxNmNgo4AqTk4Zx+EltwrPFm9r2ZfQc8CnTMpu49wGNmtsbMjhN6P1oEVxU9gS/N7F9m9oOZfWdmC7JqJEi8lwBDgrJLgReB2zMUm2dm75nZySzew9uAKWY2Jdg/HVhIKJmnn2NTSWXNbLeZrcrD++OKGE8ULtJ6mVllM6trZr9J/8CSdLmk+UG3yQFCH1DVMtTbY2Y/5OE4dYBNYZSrC/wu6G45EBy7DqG/ysP1k9gklZP0vKRtQTdPKlBZUmw2x38iw7G/AQTUysM5EMT7TZCY0m0L2km3PYf6dYHrM70PlwBnmtkh4Ebg18BuSe9LahxmXK4I8kThok7QzTMe+CeQZGaVCT0RpQzFMj+ul9vje9uBc8I4/Hbg0SB5pS/lzGxMmOFnFcvvgEZAazOrCKR3WYmf2w7ck+n4Zc3s01zOIfMxdwFVgu6qdGcBO3OokzmO1zPFUd7MHgcwsw/NrDtwJrAWeCGHtlwR54nCRaNSQGlgD3Bc0uVAj1zqfAVUlVQpm/2jgW6SbpAUJ6mqpBZZlHsB+LWk1gopL+nKTB+4mY9bP5fYKhC6L3FAUhVgaA5lnwMelHQ+gKRKkq4P9k0GzpB0v6TSkipIap0hjnqSYgDMbDvwKfCYpDLBDflfkeleRA7eAK6SdKmk2KCNTpJqS0qSdHVwr+IIkEboEWdXTHmicFEn6C4ZALwN7AduASbmUmctoZu3m4OukpqZ9n9BqPvqd4S6c5YCF2TRzkJC9ymeDo69Ebgrh0P/GRgVHPOGbMqMAMoCe4H5hG5GZ3ce7wJ/B8YG3VQrCd2oT39fugNXAV8CG4DOQdV3gp/7JC0O1m8G6hG6uniX0H2T6TmcS8Y4thO62f0QoYS9HXiA0GdGDKH3cReh97Ij8Jtw2nVFk3/hzrkICx6zfdHMXot0LM5lxa8onIsgSeUIdV1tiXQszmXHE4VzESKpBqEupFmEHg92Lip515Nzzrkc+RWFc865HBWLgcCqVatm9erVy1fdQ4cOUb58+YINqABEa1wQvbF5XHnjceVNcYxr0aJFe82seq4FzazIL8nJyZZfM2fOzHfdwhStcZlFb2weV954XHlTHOMCFloYn7He9eSccy5Hniicc87lyBOFc865HHmicM45lyNPFM4553IUVqKQNFDSSkmrJN0fbKsSTMm4IfiZGGyXpCeDqReXZ5qVLGObyZJWBOWelKSc2nXOORcZuSYKSU0JjabZitBomz0lNQT+SGgaxobAx8FrCI102TBY+gHPZtP0s8H+9LKXBduza7fALdq2n8mbjrJo2/7COoRzzhV54VxRnAfMt9A0jscJjUvzC0JDEI8KyowCegXr1wCvBY/pzic0k9eZGRsMXlc0s3nBs7yvZaqfVbsFatG2/dzywnzGbTjGrS/M92ThnHPZyHWsJ0nnAf8B2hCafOVjQnPn3m6hmcfSy+03s0RJk4HHzWxOsP1jQvP2LsxQtmVQplvwun1QpqekA1m1m0Vc/QhdkZCUlJQ8duzYPJ345E1HGb/h2I9TfLWtGUu/5mXy1EZhSktLIyEhIdJhZClaY/O48sbjypviGFfnzp0XmVnL3MrlOoSHma2R9HdgOqGZrJYBx3OoktX0jpmzUThlcotrJDASoGXLltapU6e8VKfC2fuZvHU+R46dBODTXSeoVr0SQ69qQtWE0nlqqzB88skn5PWcTpdojc3jyhuPK29Kclxh3cw2s5fM7CIz60BoRqsNwFfpXUrBz6+D4jsITQKfrjahmbAy2hFsz6pMdu0WqOS6iYzum0LvhvGM7ZfCoG7nMnXlbroNm8V/lu4ktyst55wrKcJ96qlG8PMs4FpCU05OBO4MitxJqHuKYPsdwdNPKcC3ZrY7Y3vB6+8kpQRPO92RqX5W7Ra45LqJ9DynFK3rV2Vgt4a8P6A9dauWZ+DYpfQdtZDd3x4urEM751yREe73KMZLWg1MAvqb2X7gcaC7pA2E5vF9PCg7BdhMaK7hF8gwl66kpRnavBd4MSi3CZgabM+u3UJ3blIFxt/blv935XnM3bSX7sNSGb1gGydP+tWFc67kCmuYcTNrn8W2fUDXLLYb0D+bdlpkWF8INA233dMlNkb0bV+fHk3O4I8TlvM/765k4tJdPN67OWdXi74hhp1zrrD5N7OzcVbVcozu25q/927G6t0HuWxEKiNTN3H8xMlIh+acc6eVJ4ocSOLGi8/io8Ed6XBudf42ZS3XPvspa3YfjHRozjl32niiCENSxTKMvD2Zp2+5kJ37D3PVU3MYNn09R46fiHRozjlX6DxRhEkSPZvX5KPBHbnqgpo8+fEGej45h8Vf+De6nXPFmyeKPEosX4rhN7bglbsuJu3IcXo/+ykPT17N90dz+g6ic84VXZ4o8qlz4xpMG9SBW1ufxUtztnDpiFTmbtwb6bCcc67AeaI4BRXKxPNIr2a81S+FuJgYbn1xAUPGLefbw8ciHZpzzhUYTxQFoHX9qkwd2J5fdzyHcYt30H3YLKat+jLSYTnnXIHwRFFAysTH8sfLG/Peb9pRNaE0/V5fRP83F7PnuyORDs05506JJ4oC1qx2JSbe147f9ziX6au+ovvwWby7ZIcPMuicK7I8URSC+NgY7uvSkCkDL6F+tfIMemsZfV79nJ0HfJBB51zR44miEDWoUYF3ft2WoVc1YcHmb+gxbBavz/dBBp1zRUtYgwK6/IuNEX3anU2385J46N0V/O97K5m0dBe3tzmLL745TEr9qiTX/dkEfs45FzU8UZwmdaqU47VftmLcoh0MnbiK345ZioDS8TGM7pviycI5F7W86+k0ksT1LetwV9t6QGju1x+OnWTi0swTADrnXPTwRBEBXc9Lokx8DDHBzOFvzN/KPz9cxw/HfJBB51z0CXcq1EGSVklaKWmMpDKSukhaHGwbJSkuKJso6V1JyyV9JulnkxMF5WZLWhosuyS9F2zvJOnbDPv+VHCnGx3S5+v+XY9GjOpzMb0urM3TMzdy5ZOzWbTtm0iH55xzP5HrPQpJtYABQBMzOyzpbeAW4C9AVzNbL+mvhOa3fgl4CFhqZr+Q1Bh4hqxnwmuf4Rjj+enc2LPNrOcpnFfUS66b+ON9iY6NanB1i5o8NGEF1z03jzvb1COlnD8Z5ZyLDuF2PcUBZYOrhnLAIeCIma0P9k8HegfrTYCPAcxsLVBPUlJ2DUuqAHQB3st7+MVHx3Or8+GgDtyRUpdR87byP3MOk7p+T6TDcs45FM43hiUNBB4FDgPTgNuArUBvM1so6Qmgi5k1k/Q3oIyZDZbUCvgUaG1mi7Jp+w7gajO7LnjdCRgP7AB2Ab83s1VZ1OsH9ANISkpKHjt2bJ5OPF1aWhoJCQn5qltY1u8/wUvLD/PVYXFJrThublyK8vGKdFg/isb3DDyuvPK48qY4xtW5c+dFZtYy14JmluMCJAIzgOpAPKG//G8D2gCzgc+AR4AlQfmKwCvAUuB14HPgghzan0oo4ZChfkKwfgWwIbcYk5OTLb9mzpyZ77qF6cOPZtjfp66x+g++by0fmW5TV+yKdEg/itb3zOPKG48rb4pjXMBCy+Xz1czC6nrqBmwxsz1mdgyYALQ1s3lm1t7MWgGpwIYg8Rw0sz5m1gK4I0gwW7JqWFJVoBXwfobEddDM0oL1KUC8pGphxFmslIoVf7isMf/p344aFUrz6zcWc+8bi/j6ux8iHZpzroQJJ1F8AaRIKidJhG5Mr5FUA0BSaWAI8FzwurKkUkHdvkCqmR3Mpu3rgclm9uOnn6QzguMQdF3FAPvyfmrFQ9NalXivfzv+cFkjPl77Nd2HpfLOwu0+yKBz7rTJNVGY2QJgHLAYWBHUGQk8IGkNsByYZGYzgirnAaskrQUuBwamtyVpiqSaGZq/CRiT6ZDXASslLQOeBG6yEv6pGB8bw286NWDqwPacm5TAA+OWc8fLn7H9m+8jHZpzrgQIawgPMxsKDM20+YFgyVx2HtAwm3auyPS6UxZlngaeDieukuac6gm81a8NbyzYxt+nruXSEan84dJG3NGmHjEx0XOz2zlXvPg3s4uYmBhxR5t6fDioAxfXq8KfJ63mhufnsfHrtEiH5pwrpjxRFFG1E8vxap+LGXbDBWzck8YVT8zmmZkbOXbiZKRDc84VM54oijBJXHtRbaYP6kj3Jkn834fruObpuazc+W2kQ3POFSOeKIqB6hVK88ytF/HcbcnsSTvCNc/M5e8frPVBBp1zBcITRTFyWdMz+GhQR3pfVItnP9nEFU/M5vOtPsigc+7UeKIoZiqVi+cf113AG79qzdETJ7n+uXn86T8rSTtyPNKhOeeKKE8UxdQlDasxbVAHftnubF6fv40ew2Yxc93XkQ7LOVcE+VSoxVi5UnH86aomXNn8TIaMX06fVz6n47nVaFarMp0b1/DpV51zYfErihIguW4i7w+4hOuTazFr/V6enrmRm0bOY5Hfv3DOhcETRQlROi6WetUSfpx+9dgJY8iEFXx90AcZdM7lzBNFCZJSvyql4mKIFcTFiG37DtF12Cze/twHGXTOZc/vUZQg6XN1z9+8j5T6ValSvhRDxi/nD+OXM3HZLv72i2acVbVcpMN0zkUZv6IoYZLrJtK/cwOS6yZydrXyjL07hUd6NWXp9gNcOiKVl+Zs4cRJv7pwzv2XJ4oSLiZG3JZSl2mDOpBSvwoPT17Ndc99yoavvot0aM65KOGJwgFQs3JZXr7rYkbc2IKtew9x5ZNzePLjDRw97oMMOlfShZUoJA2StErSSkljJJWR1EXS4mDbKElxQdlESe9KWi7pM0lNs2nzVUlbJC0NlhbBdkl6UtLGoI2LCu50XU4k0evCWkwf3JFLm57BsOnrufrpOSzfcSDSoTnnIijXRCGpFjAAaGlmTYFY4BZgFKHZ55oC24A7gyoPAUvNrDmhObOfyKH5B8ysRbAsDbZdTmjio4ZAP+DZvJ+WOxXVEkrz1M0X8sIdLdn//VF6PTOXx6as8UEGnSuhwu16igPKBlcN5YBDwBEzWx/snw70DtabAB8DmNlaoJ6kpDzEdA3wmoXMBypLOjMP9V0B6d4kiemDO3LjxXV4PnUzl41IZf7mEjt9uXMlVjhzZu8E/gl8AewGvgXeBuIltQyKXQfUCdaXAdcCSGoF1AVqZ9P8o0H30nBJpYNttYDtGcrsCLa5CKhYJp7Hrm3Om31bc9LgppHzGbXqCN/9cCzSoTnnThPl9kUrSYnAeOBG4ADwDjAO2AT8AygNTAOuNLMLJVUk1N10IbACaAz0NbNlmdo9E/gSKAWMBDaZ2V8lvQ88ZmZzgnIfA38ws0WZ6vcj1DVFUlJS8tixY/P1BqSlpZGQkJCvuoUpGuM6csKYsOEo07YeI7FMDHc0KUWLGtHzVZxofM/A48orjytvTiWuzp07LzKzlrkWNLMcF+B64KUMr+8A/p2pTA/g7SzqCtgKVMzlGJ2AycH688DNGfatA87MqX5ycrLl18yZM/NdtzBFa1xmZi+++5F1H/aJ1R0y2QaOWWz70o5EOiQzi973zOPKG48rb04lLmCh5ZIDzCysexRfACmSykkS0BVYI6kGQNBlNAR4LnhdWVKpoG5fINXMDmZuNP2+Q9BmL2BlsGsicEfw9FMK8K2Z7Q4jTneanFM5lsm/bc/Arg15f8Vuug2bxcRlu3wYEOeKqXDuUSwg1NW0mFBXUgyhrqIHJK0BlgOTzGxGUOU8YJWktYSeYBqY3pakKZJqBi9HS1oRtFkNeCTYPgXYDGwEXgB+c0pn6ApFqbgYBnU/l0m/vYQ6iWUZMGYJd7+2iC+/9UEGnStuwupgNrOhwNBMmx8Ilsxl5xF6tDWrdq7IsN4lmzIG9A8nLhd5jc+oyITftOOVuVv457R1dB82i4euPI+bLq5D6GLROVfU+Tez3SmLjRF929fnw/s70LRWJR6csIJbXljAtn2HIh2ac64AeKJwBaZu1fK8eXdrHru2GSt3fsulI1J5cfZmH2TQuSLOE4UrUJK4udVZTB/ckUsaVOOR99dw7bOfsu5LH2TQuaIqeh6Cd8XKGZXK8MIdLZm0fDd/nriKnk/N5hcX1qJ2YjnaNajm83U7V4R4onCFRhJXX1CTSxpUY+DYJby9cAcAT8/cyJi7UzxZOFdEeNeTK3RVypcipX5V0p+BOnr8JI9PXcPhoz7IoHNFgScKd1qk1K9K6fjQfN2xMeLzrfu5dEQqn27aG+nQnHO58K4nd1pknq/72ImT/HH8cm55YQE3tzqLB69oTMUy8ZEO0zmXBU8U7rRJrpv4k/sSUwd2YMRH63lh9mZmrP2KR3s1o1uTvIxI75w7HbzryUVM2VKxPHjFebzXvx2J5UrR97WF/HbMEvalHYl0aM65DDxRuIhrXrsyE++7hMHdz+WDlaFBBv+zdKcPMuhclPBE4aJCqbgYBnRtyPsD2lO3ankGjl3Kr0YtZNeBw5EOzbkSzxOFiyrnJlVg/L1t+d+eTZi3aR89hqcyesE2TvowIM5FjCcKF3ViY8SvLjmbD+/vwAV1KvE/767k5hfms2WvDzLoXCR4onBR66yq5XjjV635R+/mrN59kMtGpPL8rE0cP3Ey0qE5V6J4onBRTRI3XFyHjwZ3pMO51Xls6lquffZT1uz+2aSJzrlCElaikDRI0ipJKyWNkVRGUhdJi4NtoyTFBWUTJb0rabmkzyQ1zabN0ZLWBfVflhQfbO8k6VtJS4PlTwV3uq6oSqpYhpG3J/PMLRex68BhrnpqDsOmrePIcR8GxLnClmuikFQLGAC0NLOmQCxwCzAKuCnYtg24M6jyELDUzJoDdwBPZNP0aKAx0AwoS2h+7XSzzaxFsPw176fliiNJXNn8TKYP6sjVF9TkyRkb6fnkHBZ/sT/SoTlXrIXb9RQHlA2uGsoBh4AjZrY+2D8d6B2sNwE+BjCztUA9ST/7uq2ZTbEA8BlQO/+n4UqSxPKlGHZjC17pczGHjhyn97Of8tdJqzly3J+Mcq4wKJwvNUkaCDwKHAamAbcBW4HeZrZQ0hNAFzNrJulvQBkzGyypFfAp0NrMFmXTdjywABhoZrMldQLGAzuAXcDvzWxVFvX6Af0AkpKSkseOHZu3Mw+kpaWRkJCQr7qFKVrjguiK7fBxY9z6o3z8xXGqljZ+2aws51eLjXRYPxFN71dGHlfeFMe4OnfuvMjMWuZa0MxyXIBEYAZQHYgH3iOUKNoAswldDTwCLAnKVwReAZYCrwOfAxfk0P4LwIgMrysCCcH6FcCG3GJMTk62/Jo5c2a+6xamaI3LLDpjW7B5n7X6y/tWd8hke+CdpXbg+6ORDulH0fh+mXlceVUc4wIWWi6fr2YW1qCA3YAtZrYHQNIEoK2ZvQG0D7b1AM4NEs9BoE+wXcCWYPkZSUODBHRPhsR1MMP6FEn/llTNzHw8apetVmdX4eF2ZVly7ExGpm7mk3V7eLhXUy49/4xIh+ZckRfOPYovgBRJ5YIP/q7AGkk1ACSVBoYAzwWvK0sqFdTtC6Rm/PBPJ6kvcClws5mdzLD9jOA4BF1XMcC+/J6gKzlKxYohlzXmvd+0o2pCae55fRH9Ry/m4zVf8czMjSza5je9ncuPXK8ozGyBpHHAYuA4sAQYCTwiqSehD/JnzWxGUOU84DVJJ4DVwK/S25I0BehrZrsIJZZtwLwgL0yw0BNO1wH3SjpO6J7ITcElknNhaVa7EhPva8fI1M0Mn76e91fsRkDp+BhG9/UpWJ3Lq7DmozCzocDQTJsfCJbMZecBDbNp54oM61ke28yeBp4OJy7nshMfG0P/zg3Yd+gIL8/ZigE/HDvJtNVfeqJwLo/8m9muWLuyWU3KxMf8OF/3qLlbeX3eVh9k0Lk88BnuXLGWcQrW+tXK8+ZnX/C//1nFxGW7eLx3c86pHn2POzoXbTxRuGIv4xSslzU9g3GLdvDw5NVc/sRs7u/WkH7t6xMX6xfXzmXHfztciSKJ61vW4aPfdaRLoxr844N19Pr3XFbt+jbSoTkXtTxRuBKpRoUyPHd7Ms/eehFffnuEq5+ey/99uJYfjvkgg85l5onClWiXNzuTjwZ3oFeLWjwzcxNXPDmbhVu/iXRYzkUVTxSuxKtcrhT/uuECRv2yFUeOneT65+fx54mrOHTkeKRDcy4qeKJwLtDx3OpMG9SBO9vUY9S8rfQYnkrq+j2RDsu5iPNE4VwG5UvH8eerz+ede9pQOj6GO17+jN+/s4wD3x+NdGjORYwnCuey0LJeFaYMaE//zufw7pKddBuWytQVuyMdlnMR4YnCuWyUiY/lgUsbM/G+diRVLM29oxfz69cX8fXBHyIdmnOnlScK53Jxfs1KvNe/HX+4rBEz1n1Nt2GzeGfhdnysSldSeKJwLgzxsTH8plMDpg5sT6MzKvDAuOXc8fJnbP/m+0iH5lyh80ThXB6cUz2Bt/q14eFrzmfxtv1cOiKVV+du4YQPMuiKMU8UzuVRTIy4vU09PhzUgYvrVeHPk1Zzw/Pz2Pj1d5EOzblC4YnCuXyqnViOV/tczLAbLmDTnjSueGIOT8/YwLETJ3Ov7FwRElaikDRI0ipJKyWNkVRGUhdJi4NtoyTFBWUTJb0rabmkzyQ1zabNsyUtkLRB0lvp06dKKh283hjsr1dQJ+tcQZPEtRfVZvqgjnRvksQ/p63n6qfnsnKnDzLoio9cE4WkWsAAoKWZNQVigVuAUYSmKW1KaErTO4MqDwFLzaw5cAfwRDZN/x0YbmYNgf38d8rUXwH7zawBMDwo51xUq16hNM/cehHP3ZbM3rQjXPPMXB6fupZ5m/YyedNRn6/bFWnhdj3FAWWDq4ZywCHgiJmtD/ZPB3oH602AjwHMbC1QT1JSxsYUmiS7CzAu2DQK6BWsXxO8JtjfNSjvXNS7rOkZfDSoI70vqsVzszZxywsLGL/hGLe+ON+ThSuyFM6z4JIGAo8Ch4FpwG3AVqC3mS2U9ATQxcyaSfobUMbMBktqBXwKtDazRRnaqwbMD64akFQHmGpmTSWtBC4zsx3Bvk1B/b2ZYuoH9ANISkpKHjt2bL7egLS0NBISom+Ws2iNC6I3tmiL66UVR5i9878DC15zTjy/aFgqghH9VLS9X+k8rrw5lbg6d+68yMxa5lYu1xnuJCUS+iv/bOAA8A5wK3ATMFxSaULJI/034nHgCUlLgRXAkgz7fmw2i0NZGPv+u8FsJDASoGXLltapU6fcTiVLn3zyCfmtW5iiNS6I3tiiLa4KZ+/n8xfm88Px0M3tuV+KXh2a0LlRjQhHFhJt71c6jytvTkdc4XQ9dQO2mNkeMzsGTADamtk8M2tvZq2AVGADgJkdNLM+ZtaC0D2K6sCWTG3uBSqn3wAHagO7gvUdQB2AYH8lwCcIcEVOct1ERt+dwnUN4/nbL5pRuXwp+rzyOYPfWsr+Qz7IoCs6wkkUXwApksoF9wq6Amsk1YDQU0rAEOC54HXl9CeYgL5AqpkdzNighfq7ZgLXBZvuBP4TrE/kvzfGrwNmmI+V4Iqo5LqJ9DynFLe0Pov3B1zCgC4NmLhsF92GzWLy8l0+DIgrEnJNFGa2gNBN5cWEupJiCHX5PCBpDbAcmGRmM4Iq5wGrJK0FLgcGprclaYqkmsHLIcBgSRuBqsBLwfaXgKrB9sHAH0/tFJ2LDqXjYhncoxET77uEmpXLct+bS7jn9UV85YMMuiiX6z0KADMbCgzNtPmBYMlcdh7QMJt2rsiwvhlolUWZH4Drw4nLuaKoSc2KvPubtrw0ZwvDpq+n27BZ/L8rz+OGlnXwB/xcNPJvZjsXAXGxMdzT8Rw+uL8D551ZkSHjV3DbSwv4Yp8PMuiijycK5yLo7GrlGXt3Co/0asqy7d9y6YhUXprjgwy66OKJwrkIi4kRt6XUZdqgDqTUr8LDk1fT+9lPWf+VDzLoooMnCueiRM3KZXn5rosZcWMLtu07xJUOy4AhAAAZ5ElEQVRPzubJjzdw9LgPMugiyxOFc1FEEr0urMX0wR25rOmZDJu+nqufnsOy7QciHZorwTxROBeFqiWU5qmbL+SFO1qy//uj/OLfc3lsyhoOHz0R6dBcCeSJwrko1r1JEtMHd+TGi+vwfOpmLn8ilfmb90U6LFfCeKJwLspVLBPPY9c2582+rTlpcNPI+Tz07goO/nAs0qG5EsIThXNFRNsG1fjw/g7c3f5sxn72BT2GpTJj7VeRDsuVAJ4onCtCypaK5X+ubMKE37SjUtl4fvnqQgaOXcK+tCORDs0VY54onCuCWtSpzKTfXsL93RoyZcVuug9PZeIyH2TQFQ5PFM4VUaXiYri/27lM/m176lQpx4AxS7j7tYV8+a0PMugKVliDAjrnolejMyow4d62vDJ3C/+cto7uw2Zxa8pZVCgTR0r9aiTXTYx0iK6I80ThXDEQGyP6tq9P9yZJ9B+9mOdmbQagdNxG3rw7xZOFOyXe9eRcMVK3ankub3bGj/MJHzl+kqc+3uCDDLpT4onCuWImpX41SsfHECuIEXyyfg/X/nsu6770QQZd/oSVKCQNkrRK0kpJYySVkdRF0uJg26j0+a8lVZI0SdKyoE6fLNqrIGlphmWvpBHBvrsk7cmwr2/BnrJzxVty3URG901hcI9GvHNPG566+UJ27D9Mz6dmM3z6eh9k0OVZrvcoJNUCBgBNzOywpLeBW4C/AF3NbL2kvxKa5/oloD+w2syuklQdWCdptJn9OJu8mX0HtMhwjEXAhAyHfcvM7iuA83OuREqum/jjfYlkoF2Davx10iqe+HgDU1fu5h/XXRDZAF2REm7XUxxQNrhqKAccAo6Y2fpg/3Sgd7BuQAWF5nRMAL4BjmfXsKSGQA1gdt7Dd86Fo0r5Uoy46UJevqsl3/1wnGv/PZcxa4/w/dFsfzWd+5HC+YKOpIHAo8BhYBpwG7AV6G1mCyU9AXQxs2aSKgATgcZABeBGM3s/h7b/BFQ0s98Hr+8CHgP2AOuBQWa2PYt6/YB+AElJScljx44N95x/Ii0tjYSEhHzVLUzRGhdEb2weV3gOHzfeXneUmduPU72s6NO0NE2qxkY6rB9F2/uVrjjG1blz50Vm1jLXgmaW4wIkAjOA6kA88B6hRNGG0FXAZ8AjwJKg/HXAcEBAA2ALoUSQXfurgeQMr6sCpYP1XwMzcosxOTnZ8mvmzJn5rluYojUus+iNzePKm2fHf2Qd/zHD6g6ZbEPGLbMD3x+NdEhmFr3vV3GMC1houXy+mllYXU/dgC1mtsfMjhG6l9DWzOaZWXszawWkAhuC8n2ACUEcG4NE0TirhiVdAMSZ2aIMiWufmaUPXPMCoS5W51wBa1wllg/u78A9Hevz9sLt9Bg+i+mrfZBB93PhJIovgBRJ5YL7Dl2BNZJqAEgqDQwBnstQvmuwLwloBGzOpu2bgTEZN0g6M8PLq4E14Z2Kcy6vysTH8uDl5/Fe/3YklivF3a8t5L43F7PXBxl0GeSaKMxsATAOWAysCOqMBB6QtAZYDkwysxlBlYeBtpJWAB8DQ8xsL4CkpZmav4FMiQIYEDxWu4zQ01Z35efEnHPha167MhPvu4TfdT+Xaau+ovuwWby3ZKcPMuiAMIfwMLOhwNBMmx8IlsxldwE9smmnRabX9bMo8yDwYDhxOecKTqm4GH7btSGXNT2DP4xfzv1vLeU/S3fy6C+aUbNy2UiH5yLIv5ntnPuJhkkVGPfrtvypZxPmb/6GHsNTeX3+Nk76MCAllicK59zPxMaIX15yNtMGdaBFncr873sruemF+WzZeyjSobkI8EThnMtWnSrleP1XrfhH7+as2X2Qy0ak8tysTRw/4cOAlCSeKJxzOZLEDRfX4aPBHel4bnUen7qWXv+ey+pdByMdmjtNPFE458KSVLEMz9+ezDO3XMSX3/7A1U/P4V/T1nHk+IlIh+YKmScK51zYJHFl8zOZPqgjV7eoyVMzNnLlk3NYtG1/pENzhcgThXMuzxLLl2LYDS14tc/FHD56guue+5S/TFrFoSM+yGBx5InCOZdvnRrV4MNBHbg9pS6vzN1Kp3/O5IFxy/wKo5jxROGcOyUJpeP46zVNeeSapuz97ijvLNzBDc/NI3XdnkiH5gqIJwrnXIH49odjKJis+4QZ945exAcrv4xsUK5AeKJwzhWIlPpVKRUXmqu7VFwM1SuU5tdvLOI3oxfx9Xc/RDo8dwrCGuvJOedykz5X9/zN+0ipX5XmtSsxMnUzT3y0gbkb9/Gnnk249qJaKP2ywxUZfkXhnCswyXUT6d+5Acl1E4mPjaF/5wZMGdieBjUS+N07y7jzlc/Zsf/7SIfp8sgThXOuUDWokcA797ThL1efz8Kt33Dp8FRem7fVBxksQjxROOcKXUyMuLNtPT68vwMX1U3kT/9ZxY0j57FpT1qkQ3NhCCtRSBoUTCa0UtIYSWUkdZG0ONg2SlJcULaSpEmSlgV1+mTT5ieS1klaGiw/zpgn6S1JGyUtkFSvoE7WORdZdaqU47VftuKf11/A+q/SuPyJ2fz7k40c80EGo1quiUJSLUIzzbU0s6ZALHALMAq4Kdi2DbgzqNIfWG1mFwCdgH9JKpVN87eaWYtg+TrY9itgv5k1AIYDf8/fqTnnopEkrkuuzfTBHejauAb/+GAdvZ6Zy8qd30Y6NJeNcLue4oCywVVDOeAQcMTM1gf7pwO9g3UDKgTzaycA3wB5+V7/NYSSEISmYO0qf0zCuWKnRoUyPHtbMs/eehFfHTzCNc/M5f8+XMvRE37vItoonDlxJQ0EHgUOA9OA24CtQG8zWyjpCaCLmTWTVAGYCDQGKgA3mtn7WbT5CVAVOAGMBx4xM5O0ErjMzHYE5TYBrdPn3c5Qvx/QDyApKSl57Nix+Th9SEtLIyEhIV91C1O0xgXRG5vHlTfRFFfaUWPsuqPM2XmcpLJG3+ZlaZgYG+mwfiKa3q+MTiWuzp07LzKzlrkWNLMcFyARmAFUB+KB9wglijbAbOAz4BFgSVD+OkJdRgIaAFuAilm0Wyv4WYFQ8rkjeL0KqJ2h3Cagak4xJicnW37NnDkz33ULU7TGZRa9sXlceRONcc1a97VdNPR9q/fHyfan91bYdz8ci3RIP4rG98vs1OICFlouOcDMwup66gZsMbM9ZnYMmAC0NbN5ZtbezFoBqcCGoHwfYEIQx8YgUTTOIkHtDH5+B7wJtAp27QDqAARdXZUIdV8554q5DudW59FLynJnm3q8Nn8blw5PZdZ6HzMq0sJJFF8AKZLKBfcKugJrMj6lBAwBnstQvmuwLwloBGzO2KCkOEnVgvV4oCewMtg9kf/eGL8OmBFkPudcCVAmTvz56vN55542lImP4c6XP+N3by/jwPdHIx1aiZVrojCzBYRuKi8GVgR1RgIPSFoDLAcmmdmMoMrDQFtJK4CPgSEW3F+QtDQoUxr4UNJyYCmwE3gh2PcSUFXSRmAw8MdTPkvnXJHTsl4V3h/Qnvs6N+C9pTvpNiyVqSt2RzqsEimssZ7MbCgwNNPmB4Ilc9ldQI9s2mkR/DwEJGdT5gfg+nDics4Vb2XiY/n9pY24vNkZ/GHccu4dvZjLzj+Dv15zPjUqlol0eCWGfzPbORf1zq9Zif/0b8eQyxozY93XdBs2i3cWbsd7pU8PTxTOuSIhLjaGezudw9SB7Wl0RgUeGLecO17+jO3f+CCDhc0ThXOuSDmnegJv9WvDw9ecz+Jt+7l0RCqvzt3Cwq3f8MzMjT4NayHw+Sicc0VOTIy4vU09upyXxEMTVvDnSauRQl/eKhUXw+i+KSTXTYx0mMWGX1E454qsWpXL8mqfi7mi2RmYwUmDI8dO8ummvblXdmHzROGcK9Ik8atL6lM6LvRxZsC4RTt8kMEC5InCOVfkJddN5M27U3jg0kb84bJGfH/0BNc8M5fHp67lh2MnIh1ekef3KJxzxUJy3cQf70vc2qouf5uyhudmbWLaqi95vHdzWp1dJcIRFl1+ReGcK3YqlYvn79c1Z3Tf1hw7eZIbnp/H/763ku9+OBbp0IokTxTOuWKrXYNqfHh/B37Z7mzeWBAaZHDmuq9zr+h+whOFc65YK1cqjj9d1YTx97alfOk4+rzyOYPfWsr+Qz7IYLg8UTjnSoSLzkpk8oBLGNClAROX7aLbsFlMXr7LhwEJgycK51yJUToulsE9GjHpt5dQK7Es9725hHteX8RXB3+IdGhRzROFc67EOe/Miky4ty0PXdGYWev30G3YLN76/Au/usiGJwrnXIkUFxtDvw7n8OH9HWhyZkWGjF/BrS8u4It9PshgZp4onHMlWr1q5RlzdwqP/qIpy3d8y6UjUnlpzhZOnPSri3RhJQpJgyStkrRS0hhJZSR1kbQ42DYqmN8aSZUkTZK0LKjTJ4v2ykl6X9LaoMzjGfbdJWmPpKXB0rfgTtc5534uJkbc2rou0wd3oM05VXl48mp6P/sp67/6LtKhRYVcE4WkWsAAoKWZNQVigVuAUcBNwbZt/Hee6/7AajO7AOgE/EtSqSya/qeZNQYuBNpJujzDvrfMrEWwvJjPc3POuTw5s1JZXrqzJU/c1IIvvvmeK5+czRMfbeDo8ZORDi2iwu16igPKBlcN5YBDwBEzWx/snw70DtYNqCBJQALwDXA8Y2Nm9r2ZzQzWjxKaj7v2qZyIc84VBElc06IW0wd14PKmZzL8o/Vc/fQcNn9bcseMUjh3+SUNBB4FDgPTgNuArUBvM1so6Qmgi5k1k1QBmAg0BioAN5rZ+zm0XZlQouhmZpsl3QU8BuwB1gODzGx7FvX6Af0AkpKSkseOHRv2SWeUlpZGQkJCvuoWpmiNC6I3No8rbzyu8Cz5+jivrTrKgSMnubReKX7RMJ7SsYp0WD86lferc+fOi8ysZa4FzSzHBUgEZgDVgXjgPUKJog0wG/gMeARYEpS/DhhOaA6RBsAWoGI2bccBU4H7M2yrCpQO1n8NzMgtxuTkZMuvmTNn5rtuYYrWuMyiNzaPK288rvB9e/io3fnUB1Z3yGTr+I8Z9unGvZEO6Uen8n4BCy2Xz1czC6vrqRuwxcz2mNkxYALQ1szmmVl7M2sFpAIbgvJ9gAlBHBuDRNE4m7ZHAhvMbESGxLXPzI4EL18AksOI0TnnCk3FMvHc1bQ0b97dGgNufmE+D727goMlZJDBcIYZ/wJIkVSOUNdTV2ChpBpm9rWk0sAQQl1T6eW7ArMlJQGNgM2ZG5X0CFAJ6Jtp+5lmtjt4eTWwJu+n5ZxzBa/tOdX4YGAHhn+0nhdnb2bGmq/p064ex08aKfWrFtvpV3NNFGa2QNI4QvcRjgNLCF0JPCKpJ6Eb4s+a2YygysPAq5JWEOp+GmJmewEkLTWzFpJqA/8DrAUWh+5787SFnnAaIOnq4FjfAHcV2Nk659wpKlsqloeuOI8rmp3JgDeX8NjUtQCUjovhzbuL51zdYU1cZGZDgaGZNj8QLJnL7gJ6ZNNOi+DnDkJJJKsyDwIPhhOXc85FSos6lbmuZW2GT1+PAUeOn+SVuVu46KzKBH/8Fhv+zWznnMundg2qUTo+hhiBBJOX76bvqIXs/vZwpEMrUD4VqnPO5VNy3URG901h/uZ9tDq7Csu2H+Cf09bRY1gqD15xHjddXIeYmKJ/deGJwjnnTkHGubovrleFHk3O4I8TlvPQuyuYuGwnj1/bnHrVykc4ylPjXU/OOVeAzqpajtF9W/P4tc1YtfMgl45IZWTqJo6fKLrDgHiicM65AiaJm1qdxfTBHWnfsDp/m7KW3s9+ytovD0Y6tHzxROGcc4XkjEpleOGOZJ66+UJ27D9MzyfnMGz6eo4cL1rjRnmicM65QiSJqy6oyfTBHbnqgpo8+fEGrnpqDku+2B/p0MLmicI5506DKuVLMfzGFrxy18V898Nxrn32Ux6evJrvjx7PvXKEeaJwzrnTqHPjGkwb1IFbW5/FS3O2cOmIVOZu3BvpsHLkicI5506zCmXieaRXM97ql0JcTAy3vriAP45fzreHo3OQQU8UzjkXIa3rV2XqwPbc07E+by/cTvdhs5i26stIh/Uzniiccy6CysTH8uDl5/Fe/3ZUKV+Kfq8v4r43F7M37UjulU8TTxTOORcFmteuzKTfXsLve5zLtFVf0W3YLN5dsiN9QreI8kThnHNRIj42hvu6NGTKwEuoX608g95axi9f/ZxdByI7yKAnCuecizINalTgnV+3ZehVTZi/+Ru6D5vF6/O3cfJkZK4uPFE451wUio0RfdqdzbRBHbjwrET+972V3DRyPpv3pJ32WMJKFJIGSVolaaWkMZLKSOoiaXGwbZSkuKBsJUmTJC0L6vTJps1kSSskbZT0pIKZPiRVkTRd0obgZ/GbLso558JUp0o5Xv9VK/5xXXPWfnmQy5+YzXOzTu8gg7kmCkm1gAFASzNrCsQCtwCjgJuCbduAO4Mq/YHVZnYB0An4l6RSWTT9LNAPaBgslwXb/wh8bGYNgY+D1845V2JJ4oaWdfhocEc6NarO41PX0uvfcxm/aAeTNx1l0bbCHQ4k3K6nOKBscNVQDjgEHDGz9cH+6UDvYN2ACsEVQgKhea9/8h11SWcCFc1snoVu6b8G9Ap2X0MoCRH87IVzzjlqVCzD87e35NlbL2L7N4f53TvLGLfhGLe+OL9Qk4XCefRK0kDgUeAwMA24DdgK9DazhZKeALqYWTNJFYCJQGOgAnCjmb2fqb2WwONm1i143R4YYmY9JR0ws8oZyu43s591P0nqR+iKhKSkpOSxY8fm/eyBtLQ0EhIS8lW3MEVrXBC9sXlceeNx5U20xTV+/VEmbQ59kzsGuLZhPD3PyarzJnudO3deZGYtcyuX6wx3wT2Ca4CzgQPAO8CtwE3AcEmlCSWP9KuGS4GlQBfgHGC6pNlmlnEg9qzmBszT7XwzGwmMBGjZsqV16tQpL9V/9Mknn5DfuoUpWuOC6I3N48objytvoi2uCmfvZ/qL8zl67CSl4mO4udvFP860V9DC6XrqBmwxsz1mdgyYALQNuo3am1krIBXYEJTvA0ywkI3AFkJXFxntAGpneF0b2BWsfxV0TaV3UX2dnxNzzrniLH2+7msbxjO6b0qhJQkIL1F8AaRIKhfcd+gKrJFUAyC4ohgCPJehfNdgXxLQCNicsUEz2w18JyklaPMO4D/B7on898b4nRm2O+ecyyC5biI9zylVqEkCwkgUZrYAGAcsBlYEdUYCD0haAywHJpnZjKDKw0BbSSsIPbU0xMz2AkhamqHpe4EXgY3AJmBqsP1xoLukDUD34LVzzrkIyfUeBYCZDQWGZtr8QLBkLrsL6JFNOy0yrC8EmmZRZh/BFYlzzrnI829mO+ecy5EnCueccznyROGccy5Hniicc87lKKxvZkc7SXsIjTeVH9WAaJzZPFrjguiNzePKG48rb4pjXHXNrHpuhYpFojgVkhaG8xX20y1a44Lojc3jyhuPK29Kclze9eSccy5Hniicc87lyBNFMLBgFIrWuCB6Y/O48sbjypsSG1eJv0fhnHMuZ35F4ZxzLkeeKJxzzuWoRCcKSZdJWidpo6SomJtb0suSvpa0MtKxZCSpjqSZktZIWhXMehhxkspI+kzSsiCuv0Q6powkxUpaImlypGNJJ2mrpBWSlkpaGOl40kmqLGmcpLXB/7M2URBTo+B9Sl8OSro/0nEBSBoU/J9fKWmMpDKFdqySeo9CUiywntBQ5juAz4GbzWx1hOPqAKQBr5nZz0bXjZRgEqkzzWxxMN3tIqBXFLxfAsqbWZqkeGAOMNDM5kcyrnSSBgMtCc0R3zPS8UAoUQAt04f/jxaSRgGzzexFSaWAcmZ2INJxpQs+M3YCrc0sv1/wLahYahH6v97EzA5LehuYYmavFsbxSvIVRStgo5ltNrOjwFhCU75GlJmlAt9EOo7MzGy3mS0O1r8D1gC1IhsVBDMppgUv44MlKv76kVQbuJLQvCsuB5IqAh2AlwDM7Gg0JYlAV2BTpJNEBnFAWUlxQDn+O0togSvJiaIWsD3D6x1EwQdfUSCpHnAhsCCykYQE3TtLCU2bOz2YbCsajAD+AJyMdCCZGDBN0iJJ/SIdTKA+sAd4Jeiqe1FS+UgHlclNwJhIBwFgZjuBfxKaUXQ38K2ZTSus45XkRKEstkXFX6LRTFICMB6438wORjoeADM7EUyKVRtoJSniXXaSegJfm9miSMeShXZmdhFwOdA/6O6MtDjgIuBZM7sQOARExX1DgKAr7GrgnUjHAiApkVAPyNlATaC8pNsK63glOVHsAOpkeF2bQrx0Kw6CewDjgdFmNiHS8WQWdFV8AlwW4VAA2gFXB/cDxgJdJL0R2ZBCglkoMbOvgXcJdcNG2g5gR4arwXGEEke0uBxYbGZfRTqQQDdgi5ntMbNjwASgbWEdrCQnis+BhpLODv5auAmYGOGYolZw0/glYI2ZDYt0POkkVZdUOVgvS+gXaG1kowIze9DMaptZPUL/t2aYWaH9xRcuSeWDhxEIunZ6ABF/ws7MvgS2S2oUbOoKRPRBiUxuJkq6nQJfACmSygW/m10J3TcsFGHNmV0cmdlxSfcBHwKxwMtmtirCYSFpDNAJqCZpBzDUzF6KbFRA6C/k24EVwf0AgIfMbEoEYwI4ExgVPJESA7xtZlHzKGoUSgLeDX22EAe8aWYfRDakH/0WGB384bYZ6BPheACQVI7Q05H3RDqWdGa2QNI4YDFwHFhCIQ7lUWIfj3XOOReektz15JxzLgyeKJxzzuXIE4VzzrkceaJwzjmXI08UzjnncuSJwjnnXI48UTjnnMvR/wePyZuDIFy9MAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "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 `SummedField` is interpolated separately to the particle location, and that the different velocities are added in each step of the RK4 advection. So `SummedFields` are effortless to users\n", "\n", "Note that `SummedFields` 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 `SummedField`." ] }, { "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 }