{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial on Parallel Parcels with MPI "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parcels can be run in Parallel with MPI. To do this, first follow the installation instructions at http://oceanparcels.org/#parallel_install."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Splitting the ParticleSet with MPI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have installed a parallel version of Parcels, you can simply run your script with\n",
    "\n",
    "```\n",
    "mpirun -np <np> python <yourscript.py>\n",
    "```\n",
    "Where `<np>` is the number of processors you want to use\n",
    "\n",
    "Parcels will then split the `ParticleSet` into `<np>` smaller ParticleSets, based on a `sklearn.cluster.KMeans` clustering. Each of those smaller `ParticleSets` will be executed by one of the `<np>` MPI processors.\n",
    "\n",
    "Note that in principle this means that all MPI processors need access to the full `FieldSet`, which can be Gigabytes in size for large global datasets. Therefore, efficient parallelisation only works if at the same time we also chunk the `FieldSet` into smaller domains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chunking the FieldSet with dask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic idea of field-chunking in Parcels is that we use the `dask` library to load only these regions of the `Field` that are occupied by `Particles`. The advantage is that if different MPI procesors take care of Particles in different parts of the domain, each only needs to load a small section of the full `FieldSet` (although note that this load-balancing functionality is still in [development](#Future-developments)). Furthermore, the field-chunking in principle makes the `indices` keyword superfluous, as Parcels will determine which part of the domain to load itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default behaviour is for `dask` to control the chunking, via a call to `da.from_array(data, chunks='auto')`. The chunksizes are then determined by the layout of the NetCDF files. \n",
    "\n",
    "However, in tests we have experienced that this may not necessarily be the most efficient chunking. Therefore, Parcels provides control over the chunk-size via the `field_chunksize` keyword in `Field` creation, which requires a tuple that sets the typical size of chunks for each dimension.\n",
    "\n",
    "It is strongly encouraged to explore what the best value for chunksize is for your experiment, which will depend  on the `FieldSet`, the `ParticleSet` and the type of simulation (2D versus 3D). As a guidance, we have found that chunksizes in the zonal and meridional direction of approximately `500` are typically most effective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below shows an example script to explore the scaling of the time taken for `pset.execute()` as a function of zonal and meridional `field_chunksize` for a dataset from the [CMEMS](http://marine.copernicus.eu/) portal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "%pylab inline\n",
    "from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4\n",
    "from datetime import timedelta as delta\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def set_cmems_fieldset(cs):\n",
    "    files = sorted(glob('/Users/erik/Desktop/CMEMSdata_tmp/mercatorglorys12v1_gl12_mean_201607*.nc'))\n",
    "    variables = {'U': 'uo', 'V': 'vo'}\n",
    "    dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time'}\n",
    "\n",
    "    if cs not in ['auto', False]:\n",
    "        cs = (1, cs, cs)\n",
    "    return FieldSet.from_netcdf(files, variables, dimensions, field_chunksize=cs)\n",
    "\n",
    "func_time = []\n",
    "chunksize = [50, 100, 200, 400, 800, 1000, 1500, 2000, 2500, 4000, 'auto', False]\n",
    "for cs in chunksize:\n",
    "    \n",
    "    fieldset = set_cmems_fieldset(cs)\n",
    "    pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=[0], lat=[0], repeatdt=delta(hours=1))\n",
    "\n",
    "    tic = time.time()\n",
    "    pset.execute(AdvectionRK4, dt=delta(hours=1))\n",
    "    func_time.append(time.time()-tic)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAAGuCAYAAADClqRVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl43GW99/HPPUtmkjTJpEmXpEmaLlAoJU2ge1E2saCIBUQWlUUf1IMeVM6pwvERPC6PnIMrej2P4gIIKCBCRRE5bog2BSkkpdAFuqRZmm5ptqaZJJO5nz9mEpI2TabNTCaZ3/t1Xbky8/vN8o3XSPPJ9/59b2OtFQAAAABg4nMluwAAAAAAQHwQ8AAAAAAgRRDwAAAAACBFEPAAAAAAIEUQ8AAAAAAgRRDwAAAAACBFEPAAAAAAIEUQ8AAAAAAgRRDwAAAAACBFeJJdQCzy8/NtaWlpsssAAAAAgKR45ZVXDlprp4z0uAkR8EpLS7Vhw4ZklwEAAAAASWGM2R3L41iiCQAAAAApgoAHAAAAACmCgAcAAAAAKWJCXIMHAAAAwHl6enpUX1+vYDCY7FLGjN/vV1FRkbxe70k9n4AHAAAAYFyqr69XVlaWSktLZYxJdjkJZ61VU1OT6uvrNWvWrJN6DZZoAgAAABiXgsGg8vLyHBHuJMkYo7y8vFF1LAl4AAAAAMYtp4S7PqP9eQl4AAAAAJAiCHgAAAAAEAdr167V5s2bk1pDwgKeMabYGPNXY8wWY8wbxpjPRI9/2RjTYIypjn69J1E1AAAAAMBYSemAJykk6d+stadLWibpU8aY+dFz37HWlke/fp/AGgAAAADgpK1evVpnn322zjjjDN13332SpEmTJvWff+KJJ3TjjTeqsrJSTz/9tNasWaPy8nLt2LFD1dXVWrZsmcrKynT55Zerubk54fUmbJsEa22jpMbo7XZjzBZJMxL1fgAAAABS3P3vPfbYGaulJTdL3UekR6469nz5dVLFh6SOJunx6wefu+mZEd/yZz/7mSZPnqzOzk4tXrxYV1555ZCPW7FihS677DJdeuml+sAHPiBJKisr0/e//32de+65uvPOO/Wf//mf+u53vzvie47GmFyDZ4wplVQh6aXooU8bY14zxvzMGJM7FjUAAAAAwIm69957tXDhQi1btkx1dXV66623Ynpea2urWlpadO6550qSbrjhBr3wwguJLFXSGGx0boyZJOnXkj5rrW0zxvw/SV+VZKPfvyXpo0M87+OSPi5JJSUlw77H2qoG3fPcNu1p6VRhIF1rVs3T6gqahQAAAEBKGa7jlpYx/PnMvJg6dgM9//zz+tOf/qT169crIyND5513noLB4KCtDEazZ10iJLSDZ4zxKhLuHrHWPilJ1tp91tpea21Y0o8lLRnqudba+6y1i6y1i6ZMmXLc91hb1aA7nnxNDS2dspIaWjp1x5ObtLaqIf4/EAAAAADHaG1tVW5urjIyMrR161a9+OKLkqRp06Zpy5YtCofDeuqpp/ofn5WVpfb2dklSTk6OcnNz9fe//12S9NBDD/V38xIpkVM0jaSfStpirf32gOMFAx52uaTXR/M+9zy3TZ094UHHOnt6dc9z20bzsgAAAAAc7uKLL1YoFFJZWZm+9KUvadmyZZKku+++W5deeqkuuOACFRS8HW+uueYa3XPPPaqoqNCOHTv04IMPas2aNSorK1N1dbXuvPPOhNecyCWaKyV9RNImY0x19Nh/SLrWGFOuyBLNGkmfGM2b7GnpPKHjAAAAABALn8+nZ599dshzfYNUBlq5cuUx2yT0df3GSiKnaP5DkhniVFy3RSgMpKthiDBXGEiP59sAAAAAwLg3JlM0E2nNqnlK97oHHUv3urVm1bwkVQQAAAAAyZHwKZqJ1jct83OPVctKmsEUTQAAAAAONeE7eFIk5OVNStN1S0u07vYLCHcAAAAAHCklAp4k+TxuBXt6k10GAAAAACRN6gQ8r0tdR22XAAAAAABOkjIBz+9xqytEBw8AAABA/LjdbpWXl/d/1dTUHPexNTU1WrBgwdgVN4QJP2Slj9/rUpAOHgAAAIA4Sk9PV3V19cgPHCdSpoPno4MHAAAAYAzU1NToHe94h8466yydddZZqqysPOYxb7zxhpYsWaLy8nKVlZXprbfekiQ9/PDD/cc/8YlPqLc3vhkmpTp4Bw+Hkl0GAAAAgAS56Q83HXNsVekqXXPaNeoMdeqWP91yzPn3z32/Vs9dreZgs257/rZB5+6/+P4R37Ozs1Pl5eWSpFmzZumpp57S1KlT9cc//lF+v19vvfWWrr32Wm3YsGHQ8374wx/qM5/5jD70oQ+pu7tbvb292rJlix577DGtW7dOXq9Xt9xyix555BFdf/31J/I/w7BSKOAxRRMAAABAfA21RLOnp0ef/vSnVV1dLbfbrTfffPOY5y1fvlxf//rXVV9fryuuuEKnnHKK/vznP+uVV17R4sWLJUXC49SpU+Nab8oEPJ/Hpa4Q1+ABAAAAqWq4jlu6J33Y87n+3Jg6drH4zne+o2nTpmnjxo0Kh8Py+/3HPOa6667T0qVL9cwzz2jVqlX6yU9+ImutbrjhBn3jG9+ISx1DSZlr8OjgAQAAABgLra2tKigokMvl0kMPPTTkdXQ7d+7U7Nmzdeutt+qyyy7Ta6+9pgsvvFBPPPGE9u/fL0k6dOiQdu/eHdfaCHgAAAAAcAJuueUWPfjgg1q2bJnefPNNZWZmHvOYxx57TAsWLFB5ebm2bt2q66+/XvPnz9fXvvY1vfvd71ZZWZkuuugiNTY2xrU2Y62N6wsmwqJFi+zRFy0e7Ru/36IHKmu07WuXjFFVAAAAABJpy5YtOv3005Ndxpgb6uc2xrxirV000nNTpoPn87rVFQprIgRWAAAAAEiElAl4fm/kR2HQCgAAAACnSpmA5/O4JUldPQQ8AAAAIFU4bYXeaH/elAl4fR28YIhBKwAAAEAq8Pv9ampqckzIs9aqqalpyG0XYpVC++BFOnhM0gQAAABSQ1FRkerr63XgwIFklzJm/H6/ioqKTvr5KRPwuAYPAAAASC1er1ezZs1KdhkTSuos0aSDBwAAAMDhUibg+ejgAQAAAHC4lAl4fi8dPAAAAADOljoBr3+JJh08AAAAAM6UMgHv7SWadPAAAAAAOFPKBDw6eAAAAACcLnUCXt9G51yDBwAAAMChUibg9W10zhRNAAAAAE6VOgGPDh4AAAAAh0udgOeJDlkh4AEAAABwqJQJeMYY+TwulmgCAAAAcKyUCXhSZLNzlmgCAAAAcKqUCng+j4ttEgAAAAA4VkoFPL/XzUbnAAAAABwrxQIeHTwAAAAAzpVSAc/ncStIBw8AAACAQ6VUwPN7XeqigwcAAADAoVIs4NHBAwAAAOBcKRXwfB46eAAAAACcK7UCHh08AAAAAA6WUgHP73HTwQMAAADgWCkV8HxeF/vgAQAAAHCslAp4fo+bffAAAAAAOFZqBTyvS8EeOngAAAAAnCmlAp7P41YobBXqpYsHAAAAwHlSKuD5vZEfpytEwAMAAADgPCkV8HyeyI/DMk0AAAAATpRSAc/vdUuigwcAAADAmVIy4NHBAwAAAOBEKRXw3l6iSQcPAAAAgPOkVMB7e4kmHTwAAAAAzpNSAc/npYMHAAAAwLlSK+B56OABAAAAcK6UCnh+OngAAAAAHCzFAh4dPAAAAADOlVIBr2+KZhcdPAAAAAAOlFIBr38fPDp4AAAAABwoNQMeG50DAAAAcKCUCngs0QQAAADgZCkV8Lxul9wuwxJNAAAAAI6UUgFPinTx2CYBAAAAgBOlXMDze91skwAAAADAkVIv4NHBAwAAAOBQKRfwfF43UzQBAAAAOFLqBTyPS10hOngAAAAAnCflAp6fDh4AAAAAh0q5gEcHDwAAAIBTpVzA83vd6qKDBwAAAMCBUjDgMUUTAAAAgDOlXMDzedgHDwAAAIAzpVzAo4MHAAAAwKkSFvCMMcXGmL8aY7YYY94wxnwmenyyMeaPxpi3ot9z4/m+fq9bQTp4AAAAABwokR28kKR/s9aeLmmZpE8ZY+ZLul3Sn621p0j6c/R+3Pg8LnXRwQMAAADgQAkLeNbaRmvtq9Hb7ZK2SJoh6f2SHow+7EFJq+P5vn0dPGttPF8WAAAAAMa9MbkGzxhTKqlC0kuSpllrG6VICJQ0NZ7v5fO4ZK3U3UsXDwAAAICzJDzgGWMmSfq1pM9aa9tO4HkfN8ZsMMZsOHDgQMzv5/e6JYnNzgEAAAA4TkIDnjHGq0i4e8Ra+2T08D5jTEH0fIGk/UM911p7n7V2kbV20ZQpU2J+T1804AXZ7BwAAACAw3gS9cLGGCPpp5K2WGu/PeDU05JukHR39PtvRnyxg29J97938LEzVktLbpa6j0iPXNV/eFV7UKekdci9cb90zk1SR5P0+PXHvubij0oLrpRa66UnP3Hs+RWfluZdEnnv33722PPv/HdpzvlS42vSH+449vyFd0olS6Xal6Q/f+XY8xd/Qyook3b8VXrhm8eef993pfxTpG3PSpU/OPb8FT+Scoqk138tvfyzY89/8OdSZp5U9YhU/Ytjz3/oV1JahvTPH0tvrD32/E3PRL6vu1d687nB57x+6cO/jtz+239LO/82+HxGrnT1w5Hbf/qyVPfy4PPZhdKVP47cfvZ2ae+mwefz5kiX3Ru5/fStUtOOweennyldcnfk9q9vltr2DD5fvFh615cjtx/7sHSkefD52edK534+cvvhK6We4ODzp66SVt4auX3050467mevX/l1UsWH+Ozx2eOzdzQ+e5HbfPaOPc9nL/Kdz56OwWePz57EZ2+oz94wEhbwJK2U9BFJm4wx1dFj/6FIsHvcGPMxSbWShvhf7OS5XEaSFAqzRBMAAACAsyQs4Flr/yHJHOf0hSfyWjVpXt1UMHgWy6rsLF0jqdNldMuAc80d3drmytCn0zP0KUnNbpduKzh2jsvVmZm6WNJej0d3DHH+hox0nSdpl9ejrwxx/uPpfi2XtNWXpv8a4vxn/D6VS6r2+/S9Ic5/wZem0yStT/frviHO3+n1aJak5zPS9eAQ57/h8Wi6pD9kZuqxIc5/2+1SrqS1WZn6zRDn/6/LKF3So9lZem6I8/dHvz+Qk62/HXXe5/Hph9HbPwxk66Wjzgd8AX0nevu7uQFtDA0+Py0zT9G/x+i/8nK11Qw+PzN7sr4cvf3l/MnandY+6Pxpk3P1hejt26fmaV9maND5hbkB9f0d5HNTp6ilyzvo/NJAtj4Zvf3J6VPUFeoadP7cnGzdGL199OdOOv5nr8/7szK1Wnz2+Ozx2Tsan70IPnt89o7GZy+Czx6fvYH47A3/2RvOmEzRHEsm2sHrYYomAAAAAIcxE2G/uEWLFtkNGzbE9NiXdjbp6vte1CP/a6lWzs1PcGUAAAAAkHjGmFestYtGelzKdfD8TNEEAAAA4FApF/B83siPxD54AAAAAJwm5QKe30MHDwAAAIAzpV7A61+iSQcPAAAAgLOkXMDzefqWaNLBAwAAAOAsKRfw6OABAAAAcKqUC3h9HTyuwQMAAADgNCkX8FwuozS3iymaAAAAABwn5QKeFNkqgQ4eAAAAAKdJzYDncTNkBQAAAIDjeI53whgzOYbnh621LXGsJy78Xpe6GLICAAAAwGGOG/Ak7Yl+mWEe45ZUEteK4sDvdStIBw8AAACAwwwX8LZYayuGe7IxpirO9cSFz0MHDwAAAIDzDHcN3vIYnh/LY8YcHTwAAAAATnTcgGetDUqSMWaOMcYXvX2eMeZWY0xg4GPGG7/XxUbnAAAAABwnlimav5bUa4yZK+mnkmZJ+kVCqxolpmgCAAAAcKJYAl7YWhuSdLmk71prPyepILFljQ4dPAAAAABOFEvA6zHGXCvpBkm/ix7zJq6k0fN73Gx0DgAAAMBxYgl4NykyTOXr1tpdxphZkh5ObFmj4/O61BWigwcAAADAWYbbJkGSZK3dLOnWAfd3Sbo7kUWNlo8OHgAAAAAHOm4Hzxhz30hPjuUxyeDzsg8eAAAAAOcZroO32hgz3DYIRtL5ca4nLvwet7p7wwqHrVwuk+xyAAAAAGBMDBfw1sTw/L/Hq5B48nvdkqSuUFjpae4kVwMAAAAAY+O4Ac9a++BYFhJPPk9k5Wmwp5eABwAAAMAxYpmiOeEM7OABAAAAgFOkaMB7u4MHAAAAAE6RkgHP56GDBwAAAMB5ht0HzxhTJOkaSe+QVCipU9Lrkp6R9Ky1dlwmKDp4AAAAAJzouAHPGHO/pBmSfifpvyTtl+SXdKqkiyV90Rhzu7X2hbEo9ET0XYNHwAMAAADgJMN18L5lrX19iOOvS3rSGJMmqSQxZY1O3xRNlmgCAAAAcJLhtkkYKtwNPN8taXvcK4oDOngAAAAAnOi4Q1aMMb81xrzPGOMd4txsY8xXjDEfTWx5J6f/Gjw6eAAAAAAcZLglmjdLuk3Sd40xhyQdUOQavFJJOyT9wFr7m4RXeBL6p2jSwQMAAADgIMMt0dwr6fOSPm+MKZVUoMgUzTettUfGpLqT5KODBwAAAMCBYt0Hz0rKtNZWS7LGmKwE1jRqdPAAAAAAONGIAc8Yc7OkJyT9KHqoSNLaRBY1Wn3X4DFFEwAAAICTxNLB+5SklZLaJMla+5akqYksarTS3C4ZwxRNAAAAAM4SS8Drim6JIEkyxngUWbI5bhlj5PO4CHgAAAAAHCWWgPc3Y8x/SEo3xlwk6VeSfpvYskbP73WzRBMAAACAo8QS8G5XZIuETZI+Ien31tovJrSqOPB73HTwAAAAADjKcPvg9flXa+33JP2474Ax5jPRY+OWz+uigwcAAADAUWLp4N0wxLEb41xH3NHBAwAAAOA0x+3gGWOulXSdpFnGmKcHnMqS1JTowkbL73Up2EMHDwAAAIBzDLdEs1JSo6R8Sd8acLxd0muJLCoefB63ukJ08AAAAAA4x3EDnrV2t6TdkpaPXTnx4/O61B4MJbsMAAAAABgzIw5ZMca06+1979IkeSV1WGuzE1nYaPm9bh1o70p2GQAAAAAwZkYMeNbarIH3jTGrJS1JWEVx4vO41M0UTQAAAAAOEssUzUGstWslXZCAWuLK72WKJgAAAABniWWJ5hUD7rokLdLbSzbHLZ/HpSAdPAAAAAAOEstG5+8bcDskqUbS+xNSTRz5vW510cEDAAAA4CCxXIN301gUEm9+Lx08AAAAAM4y4jV4xpgHjTGBAfdzjTE/S2xZo+fzuNUbturpJeQBAAAAcIZYhqyUWWtb+u5Ya5slVSSupPjweyM/WhddPAAAAAAOEUvAcxljcvvuGGMmK7Zr95LK73VLEpM0AQAAADhGLEHtW5IqjTFPKDI984OSvp7QquLA56GDBwAAAMBZYhmy8nNjzAZF9r4zkq6w1m5OeGWjRAcPAAAAgNPEutH5ZEkd1trvSzpgjJmVwJriwuch4AEAAABwllimaN4l6QuS7oge8kp6OJFFxYOPISsAAAAAHCaWDt7lki6T1CFJ1to9krISWVQ8+OngAQAAAHCYWAJet7XWKjJgRcaYzMSWFB/92yT00MEDAAAA4AyxBLzHjTE/khQwxtws6U+SfpLYskav7xq8rhAdPAAAAADOEMsUzW8aYy6S1CZpnqQ7rbV/THhlo9TXwQvSwQMAAADgECMGPGPM/Gig++OAY+dZa59PZGGj5WObBAAAAAAOE+sSzc+biHRjzPclfSPRhY2Wn43OAQAAADhMLAFvqaQSSZWSXpa0R9LKRBYVD2x0DgAAAMBpYgl4PZI6JaVL8kvaZa0d920xn4dr8AAAAAA4SywB72VFAt5iSedIutYY80RCq4oDj9slj8swRRMAAACAY4w4ZEXSx6y1G6K390p6vzHmIwmsKW78XjcdPAAAAACOMWIHz1q7wRhzjjHmJkkyxuRL+kfCK4sDn8dFBw8AAACAY4wY8Iwxd0n6gqQ7oofSJD2cyKLihQ4eAAAAACeJ5Rq8yyVdJqlDkqy1eyRljfQkY8zPjDH7jTGvDzj2ZWNMgzGmOvr1npMtPBY+r0tBOngAAAAAHCKWgNdtrbWSrCQZYzJjfO0HJF08xPHvWGvLo1+/j/G1TorP41YXHTwAAAAADhHrRuc/khQwxtws6U+SfjzSk6y1L0g6NMr6RsXv5Ro8AAAAAM4x4hRNa+03jTEXSWqTNE/SndbaP47iPT9tjLle0gZJ/2atbR7Faw3L73Gz0TkAAAAAx4ilgydr7R+ttWustf8+ynD3/yTNkVQuqVHSt473QGPMx40xG4wxGw4cOHBSb+bzutQVYokmAAAAAGeIKeAdzRhz38k8z1q7z1rba60NK7LMc8kwj73PWrvIWrtoypQpJ/N2dPAAAAAAOMpJBTxJPzqZJxljCgbcvVzS68d7bDz4vC62SQAAAADgGLHsg3fVEIdnx/C8X0paL2meMabeGPMxSf9tjNlkjHlN0vmSPneiBZ8Iv8fNkBUAAAAAjjHikBVFNjj/VQzHBrHWXjvE4Z/GWFdc+OngAQAAAHCQ4wY8Y8wlkt4jaYYx5t4Bp7IlhRJdWDz4vFyDBwAAAMA5huvg7VFkK4PLJL0y4Hi7Ery0Ml78nsgUTWutjDHJLgcAAAAAEuq4Ac9au1HSRmPML6KPK7HWbhuzyuLA53VLkrpCYfmjtwEAAAAgVcUyRfNiSdWS/iBJxphyY8zTCa0qTnyeyI/XxXV4AAAAABwgloD3ZUX2q2uRJGtttaTSxJUUP/7+Dh7X4QEAAABIfbEEvJC1tjXhlSRAX8BjkiYAAAAAJ4hlm4TXjTHXSXIbY06RdKukysSWFR/9SzTp4AEAAABwgFg6eP8q6QxJXZJ+KalN0mcTWVS80MEDAAAA4CQjdvCstUckfVHSF40xbkmZ1tpgwiuLA783kl+DdPAAAAAAOMCIHTxjzC+MMdnGmExJb0jaZoxZk/jSRs/niQ5ZoYMHAAAAwAFiWaI531rbJmm1pN9LKpH0kYRWFSf9HbweOngAAABAKllb1aCVd/9Fs25/Rivv/ovWVjUku6RxIZYhK15jjFeRgPcDa22PMcYmuK646L8GjyWaAAAAQMpYW9WgO57cpM5oI6ehpVN3PLlJkrS6YkYyS0u6WDp4P5JUIylT0gvGmJmKDFoZ99joHAAAAEg9//3c1v5w16ezp1f3PLctSRWNH7EMWblX0r0DDu02xpyfuJLihw4eAAAAMPHtbw+qurZF1XUtqqpt0Z6WoWc+7mnpHOPKxp8RA54xJk/SXZLOkWQl/UPSVyQ1Jba00evr4LFNAgAAADAxBHt69XpDayTM1bWourZFDdHg5nEZnVaQpcw0tzq6j23iFAbSx7rccSeWa/AelfSCpCuj9z8k6TFJ70pUUfHS18Fjo3MAAABg/AmHrXY1dfR356rrWrSlsU2hcGTkx4xAusqLA7ppZanKiwNaMCNHfq/7mGvwJCnd69aaVfOS9aOMG7EEvMnW2q8OuP81Y8zqRBUUT3TwAAAAgPGjuaP77c5cXYs21rWotbNHkpSZ5lZZUUA3v3O2KooDKi8JaGqWf8jX6Rukcs9z27SnpVOFgXStWTXP8QNWpNgC3l+NMddIejx6/wOSnklcSfFjjFGax6UutkkAAAAAxlR3KKzNjW2qrm3u787VNB2RJLmMdOq0LF2yYLrKiwOqKMnV3KmT5HaZmF9/dcUMAt0QYgl4n5B0m6SHovfdkjqMMbdJstba7EQVFw9+j0tdITp4AAAAQKJYa1Xf3KlXB4S5N/a0qTv6e/jULJ/KiwP64OJiVRTn6syiHE3yxRJFcKJimaKZNRaFJIrf62ajcwAAACCO2oI9eq2uVVUDAl1TR7ckye916cwZObph+UyVF+eqoiSgghy/jIm9O4eTl/Kx2eelgwcAAACcrFBvWNv2tatqwCCUHQcOy0bmoGjOlEydN2+qyksCqigOaN70LHndsWy3jURI+YDn99DBAwAAAGLV2Nqp6tq3tyjY1NDaP61ycmaayosDumxhocqLA1pYHFBOujfJFWOg1A94LNEEAAAAhtTRFdKmvj3nosst97V1SZLS3C7NL8zW1YuLVVESUHlxQCWTM1hqOc7FFPCMMW5J0wY+3lpbm6ii4snHkBUAAABA4bDV9gOH+7tzVbXNenNfu6JbzmlmXoaWzc5TeXEkzM0vzJbP405u0ThhIwY8Y8y/SrpL0j5JfUnJSipLYF1x4/e6daQ7lOwyAAAAgDF1oL0res1cs6pqW/RafasOd0V+L87ye1ReHNC7509TeUlAC4sCypvkS3LFiIdYOnifkTTPWtuU6GISwe916VAHHTwAAACkrmBPr97Y09o/CKWqtkUNLZ2SJLfL6LTpWVpdUajy4lyVFwc0Oz9TrhPYcw4TRywBr05Sa6ILSRSfx62uENfgAQAAIDVYa7XrYEf/RMvquhZt3tOmUHStZWGOXxUlubpxRanKSwJaUJij9DSWWjpFLAFvp6TnjTHPSOrqO2it/XbCqoojn9elYA8dPAAAAExMzR3dqq6PTLTsC3StnT2SpIw0t8qKcnTzO2ervDiyTcHUbH+SK0YyxRLwaqNfadGvCYUOHgAAACaK7lBYWxrbBnXndh3skCQZI506NUsXnzE9MtWyJKBTpmbJzVJLDDBiwLPW/udYFJIofq9LXXTwAAAAMM5Ya1Xf3Nm/31x1XbNe39Om7ugE+ClZPpUXB/SBs4tUURJQWVFAk3wpv8sZRum4nxBjzHettZ81xvxWkamZg1hrL0toZXHi97oVpIMHAACAJGsL9ui1ulZV1zX3d+cOHu6WFNna68wZObp+2UyVlwRUUZKrwhw/e87hhA33J4CHot+/ORaFJIrP41JPr1Vv2NK+BgAAwJgI9Yb15r7Dqqpr7r92bvuBw7LRtsnsKZl656lTVFEcCXPzpmfJ63Ylt2ikhOMGPGvtK9Hvfxu7cuLP741MDOoK9SojjZY2AAAA4m9va7B/v7mquhZtqm9VZ09kFVluhlflxQFdWlao8pKAyosCysnwJrlipKqUTzx+T+QvIcGesDIm3IgYAAAAjDdHukPaVN+cWrRMAAAgAElEQVTav99cdV2L9rYFJUlet9H8whxdvbg4MtWyJKCSyRkstcSYSfmA5xvQwQMAAABORDhstePAYVUNCHNv7mtXb3TPuZLJGVoya7LKiyNTLecXZPevIAOSYcSAZ4y5ylr7q5GOjVd+79sdPAAAAGA4Bw939V8zV1XXrNfqWtXeFZIkZfk9Ki8O6F2nz4kEuuKA8ib5klwxMFgsHbw7JB0d5oY6Ni75PZG/oAR76OABAADgbcGeXr2x5+0956pqm1Xf3ClJcruMTpuepcvKC/uXWs7OnyQXQ/swzg23TcIlkt4jaYYx5t4Bp7IlhRJdWLz4oh28rhAdPAAAAKey1qqm6Uhki4Joh25zY5t6eiNLLQty/KooCej65TNVXpyrM2fkKD2NpZaYeIbr4O2RtEHSZZJeGXC8XdLnEllUPNHBAwAAcJ6WI939nbnquhZtrGtR85EeSVJGmltnzsjRx86Z3d+dm5btT3LFQHwMt03CRkkbjTG/sNb2jGFNcdU3ZIWABwAAkJq6Q2Ft3Rtdahntzu082CFJMkY6ZeokvXv+9MgWBcUBnTJ1kjzsOYcUFcs1eEuMMV+WNDP6eCPJWmtnJ7KwePF5WKIJAACQKqy1qm/uHNSde72htf93vfxJPpUXB3Tl2UWqKA7ozKIcZfnZcw7OEUvA+6kiSzJfkTTh2mB+OngAAAATVnuwR5vqWwdtU3DwcJekyB/yF8zI0YeXzVRFtDs3I5DOnnNwtFgCXqu19tmEV5Ig/R08tkkAAAAY13rDVm/ua48GuWZV17Xorf2HZSNzUDQ7P1PvPCVf5SUBVRTn6rSCLHlZagkMEkvA+6sx5h5JT0rq6jtorX01YVXFkZ+NzgEAAMalfW1BVdVG9purrm3RpoZWHemO/M4WyPCqvDig95xZoIqSXC0sylEgIy3JFQPjXywBb2n0+6IBx6ykC+JfTvyx0TkAAEDydXb3alNDq6rrmvuXWja2BiVJXrfR/IJsXXV2UX93bmZeBkstgZMwYsCz1p4/FoUkio9tEgAAAMZUOGy18+DhaHcuMtly27529YYjay2LJ6drUelklRdHrps7ozC7f9UVgNEZMeAZY6ZJ+j+SCq21lxhj5ktabq39acKriwOv28hlmKIJAACQKE2Hu/onWlbVtmhjfYvagyFJUpbPo4XFAf3LuXMiga4koPxJviRXDKSuWJZoPiDpfklfjN5/U9JjikzXHPeMMfJ73XTwAAAA4qAr1Ks39rT17zdXVdesukOdkiSXkU6bnq33LSyMbCBeHNCcKZPkcrHUEhgrsQS8fGvt48aYOyTJWhsyxkyotOTzuOjgAQAAnCBrrXY3HXm7O1fXos17WtXTG1lqWZDjV3lxQB9eOlPl0T3nMtJi+fUSQKLE8v/ADmNMniKDVWSMWSapNaFVxRkdPAAAgJG1HulRdX1LtDsX2aag+UiPJCnd61ZZUY4+es4sVRQHVF6cq+k5/iRXDOBosQS82yQ9LWmOMWadpCmSPpDQquLM73UrSAcPAACgX09vWFsb2yNTLaMdup0HOiRJxkhzp0zSRfOnqbw4V+XFAZ06bZI87DkHjHuxTNF81RhzrqR5koykbdbanoRXFkc+j0tddPAAAIBDWWu1pzWoqtrm/mvnNjW09l/Ckj8pTeXFAV15VpHKiwMqK8pRlt+b5KoBnIxYpmj6Jd0i6RxFlmn+3RjzQ2ttMNHFxYuPDh4AAHCQw10hvRa9Zq7v+rkD7V2SpDSPSwsKs/XhZTP7tykoyk1nzzkgRcSyRPPnktolfT96/1pJD0m6KlFFxZvf4+IaPAAAkJJ6w1Zv7W+PbB4e7c69ub9dNjIHRbPyM3XO3HxVlETC3GnTs5XmYaklkKpiCXjzrLULB9z/qzFmY6IKSgSf163Wzgm1qhQAAGBI+9uCqoruN1dd16xN9a3q6I78ITsn3avy4oAuXjC9P9AFMtKSXDGAsRRLwKsyxiyz1r4oScaYpZLWJbas+PJ7XNpPBw8AAEwwnd29en1Pq6prI/vNVde2aE9r5CoZj8tofmG2rjy7KBrmclWal8FSS8DhYgl4SyVdb4ypjd4vkbTFGLNJkrXWliWsujjxsU0CAAAY58Jhq50HO6LXzDWrqrZFW/e2qzccWWtZlJuus2bm6qPFAVWU5OqMwmz5ve4kVw1gvIkl4F2c8CoSzM9G5wAAYJw51NEd2WuuNjIMZWNdi9qCIUnSJJ9HC4tz9MlzZ/dvUzAly5fkigFMBLEEPI+kemttlzHmPEllkn5urW1JaGVxxEbnAAAgmbpCvdq8p61/omVVbYtqDx2RJLmMNG96tt5bVhjZQLwkoDlTJsntYqklgBMXS8D7taRFxpi5kn6qyKbnv5D0nkQWFk8+j0vBHjp4AAAg8ay1qj10pD/IVde1aPOeNnX3Rn4XmZ7tV3lxQNctLVF5cUBnzshRpi+WX8kAYGSx/NckbK0NGWOukPRda+33jTFViS4snvxet7pCvbLWcuExAACIq9bOHm0csN9cdV2LDnV0S5LSvW6dWZSjm1aWRvacKwmoICc9yRUDSGWxBLweY8y1kq6X9L7oMW/iSoo/v9elsJV6eq3SPAQ8AABwcnp6w9q2tz2ygXh0m4IdBzokScZIc6dM0oWnTVV5dIuCedOy5HGz5xyAsRNLwLtJ0iclfd1au8sYM0vSw4ktK758nsiEqa5QLxt7AgCAmFhr1dga7N9vrrquRZsaWvsv+8jLTFNFSUCXV8xQeXGuyopzlO2fUH8DB5CCRgx41trNkm4dcH+XpLsTWVS8+b2RUBfsCSvLn+RiAADAuHS4K6TX6qPLLKPXzu1v75IkpXlcOqMwW9ctmanykoAqigMqyk3n0g8A444jruj1RfeIYZImAACQpN6w1fb9h1VV29x/3dyb+9oV3XJOpXkZWjk3P3LdXHFApxdkswoIwITgjIAX/Q8ye+EBAOBM+9uD/fvNVde26LX6FnV0R/7wm5Pu1cLigFadMT1y7VxRQLmZaUmuGABOTswBzxiTaa3tSGQxieKngwcAgGMEe3r1ekProG0KGlo6JUkel9HpBdm64qwiVUQHoczKz2SpJYCUMWLAM8askPQTSZMklRhjFkr6hLX2lkQXFy99Aa8rRMADACCVhMNWu5o6+q+Zq6pr1tbGdoWiay1nBNJVXhLQTStLVVES0BmFOf2/FwBAKoqlg/cdSasU2eBc1tqNxph3JrSqOOtfoslm5wAATGjNHd3RINeiqtpmbaxrUVswJEma5POorChHH3/nbFWU5GphcY6mMl0NgMPEtETTWlt31NKFCdUK61+iSQcPAIAJozsU1ubGNlVHB6FU1bVod9MRSZLLSKdOy9J7ywpUXhxQRUmu5kyZJLeLpZYAnC2WgFcXXaZpjTFpimyZsCWxZcVXXwcvSAcPAIBxyVqrukOdqqp7e6rlGw1t6u6N/Ns9NcunipKArllcovLigMqKcpTpc8SsOAA4IbH8l/GTkr4naYakekn/I+lTIz3JGPMzSZdK2m+tXRA9NlnSY5JKJdVI+qC1tvlkCj8RXIMHAMD40hbs0cYB+81V17WoqaNbUmT/2rIZAd24srR/m4KCHD+DUAAgBrFsdH5Q0odO4rUfkPQDST8fcOx2SX+21t5tjLk9ev8LJ/HaJ2TgRucAAGBshXrD2rq3vT/IVde1aPv+w/3n506dpPNPm9of5uZNz5LXzZ5zAHAyYpmiOUvSvyrSdet/vLX2suGeZ619wRhTetTh90s6L3r7QUnPawwCns/DNgkAAMTD2qoG3fPcNu1p6VRhIF1rVs3T6ooZgx7T2NrZvz1BdW2LNjW0qjP6b/DkzDRVFAf0/oWFKi8JqKwooJx0bzJ+FABISbEs0Vwr6aeSfitptC2wadbaRkmy1jYaY6aO8vVi0tfBY6NzAABO3tqqBt3x5Kb+sNbQ0qnbn3xN2/e3K9PnVXX0+rl9bV2SpDS3S2fMyNY1S4ojg1CKc1U8OZ2llgCQQLEEvKC19t6EV3IUY8zHJX1ckkpKSkb1WnTwAAAYnZ7esL7x7Jb+cNcn2BPWD/66Q5I0My9Dy2fnRZZaluTq9IKs/n+DAQBjI5aA9z1jzF2KDFfp6jtorX31JN5vnzGmINq9K5C0/3gPtNbeJ+k+SVq0aJE9iffq53YZed2GDh4AAMfR0xvW3tag6pqPqL65M/oVud3Q3KnG1k6Fh/nX+NUvXaTJmWljVzAAYEixBLwzJX1E0gV6e4mmjd4/UU9LukHS3dHvvzmJ1zgpfo+bDh4AwLFONMAZIxVk+1WUm6GlsyarKDddP1+/Wy2dPce89oxAOuEOAMaJWALe5ZJmW2u7T+SFjTG/VGSgSr4xpl7SXYoEu8eNMR+TVCvpqhMr9+T5vG6maAIAUlY8AlxRbkb/9+k5fqV5Bk+ynD1l0qBr8CQp3evWmlXzxurHBACMIJaAt1FSQMMspxyKtfba45y68EReJ158Hhf74AEAJqyxCHAj6ZuWOdIUTQBA8sQS8KZJ2mqMeVmDr8EbdpuE8cbvdamLDh4AYJwK9YbV2BocFNwG3t7bFlTvgASXiAAXi9UVMwh0ADCOxRLw7kp4FWPA7+UaPABA8pxMgJue7VdRbrqW9Ae4t0NcQU56QgIcAGBiGzHgWWv/NhaFJFpkiSYdPABAYhDgAADjwXEDnjHmH9bac4wx7YpMzew/Jclaa7MTXl0c0cEDAIwGAQ4AMBEM18HLlCRrbdYY1ZJQPo9LbcFjRzsDACAR4AAAqWG4gDeqzcXHG7/XzZAVAHAwAhwAwAmGC3hTjTG3He+ktfbbCagnYfxet4JskwAAE9baqoZhx/MT4AAAGD7guSVNUuSauwnP53Gx0TkATFBrqxoGbbDd0NKpf//VRv3ipd1yuYzqmzvV2HpsgJuWFQlwi0tzB20hUJSbroKAXz6PO1k/EgAACTFcwGu01n5lzCpJsMgSTTp4ADAR3f3slv5w1ycUttqwu1lnleRq0UwCHAAA0vABLyU6d318XpeCbJMAABOGtVav7G7W/ZU12tvWdZzHSE/8y4oxrgwAgPFruIB34ZhVMQZ8Hre6Q2FZa2VMSmVXAEgpwZ5e/e61Rj1QuUuvN7Qp2+/RJJ9Hh7tCxzy2MJCehAoBABi/jhvwrLWHxrKQRPN7IxfKd4XC8ntZsgMA482+tqAefnG3fvFSrZo6unXqtEn6P5efqdUVhfqfN/YNugZPktK9bq1ZNS+JFQMAMP4M18FLKf7odRjBnl4CHgCME9ZavVrbogcqa/Tspkb1Wqt3nT5NN60o1fI5ef0rLvqmZQ43RRMAADgo4PkGdPAAAMnVFerV7zc16oF1NdpY36osv0c3rijV9ctLVZKXMeRzVlfMINABADACxwS8gR08AEBy7G8L6pGXavXIS7U6eLhLc6Zk6qurF+iKihnK9DnmnyQAABLGMf+a9i3LZC88ABh71XUtemDdLj2zqVGhsNUF86bqxpWlOmduPoOvAACII8cEPJ+nb4kmHTwAGAvdobCefb1R96+rUXVdi7J8Hn1kWamuXz5TpfmZyS4PAICU5JiARwcPAMbGgfYu/eKlWj380m4daO/S7PxMfeX9Z+iKs4o0iWWYAAAklGP+pe0bssI1eACQGJvqW3V/5S79bmOjunvDOm/eFN20cpbeMTdfLhfLMAEAGAuOCXh9Q1aYogkA8dPTG9YfXt+rBypr9MruZmWmuXXd0hJdv3ymZk+ZlOzyAABwHOcEPDp4ABA3TYe79Mt/1uqhF3drX1uXSvMydNf75usDZxcpy+9NdnkAADiWYwKej20SAGDUXm9o1QOVNXp64x51h8J656lTdPcVpTr31CkswwQAYBxwTMDzs9E5AJyUUG9Yz72xTw9U7tLLNc3KSHPrmsXFun55qeZOZRkmAADjiWMCns9LBw8ATsShjm798p+1evjF3WpsDapkcoa+dGlkGWZOOsswAQAYj5wT8Dx08AAgFpv3tOnByhqtrW5QVyisc+bm66vvX6DzT5sqN8swAQAY1xwV8IyRuujgAcAxQr1h/WnLPt2/rkYv7TqkdK9bHzi7SDeuKNUp07KSXR4AAIiRYwKeMUY+j0tBOngA0K/lSLcefblOD63frYaWThXlpuuL7zldH1xUrJwMlmECADDROCbgSZFJmnTwAEDaujeyDPOpqgYFe8JaMSdPd71vvi48fRrLMAEAmMAcFfD8XpeCPXTwADhTb9jqT1v26YF1NVq/s0l+r0uXV8zQDStKddr07GSXBwAA4sBhAc+tYIgOHgBnaT3So8c31OnB9TWqb+7UjEC6br/kNF29qFi5mWnJLg8AAMSRowKez+NSFx08AA7x1r52PVBZoydfbVBnT6+Wzpqs//3e0/Wu06fJ43YluzwAAJAAjgp4dPAApLresNVft+7XA5U1+sf2g/J5XFpdHlmGOb+QZZgAAKQ6RwU8n8fFRucAUlJrZ49+taFOP1+/W7WHjqggx6/PXzxP1ywu0WSWYQIA4BiOCnh+r1uHu0LJLgMA4mb7/sN6sLJGv361Xke6e7WkdLJuv+Q0vXs+yzABAHAiRwU8n8etg4e7k10GAIxKOGz1tzcP6P7KGr3w5gGluV26rLxQN64o1YIZOckuDwAAJJGzAp7XxT54ACas9mCPnnilXg9W1qim6YimZfv07+8+VdcsKVH+JF+yywMAAOOAowKe3+NWV4gpmgAmlp0HDuvn63frVxvq1NHdq7Nn5urf3j1PFy+YLi/LMAEAwADOCnhehqwAmBjCYasX3jqgBypr9Py2yDLMSxcW6MYVpSorCiS7PAAAME45KuD56OABGOcOd4X06+gyzJ0HOzQly6fbLjpV1y4p0ZQslmECAIDhOSrg0cEDMF7VHOzoX4bZ3hVSeXFA37umXJcsKFCah2WYAAAgNg4LeG6Fwlah3jDjwwEknbVW/9h+UA+sq9Fftu2Xx2X03jMLdMOKUlWU5Ca7PAAAMAE5KuD5on8F7woR8AAkT0dXSE9WNejByhpt339Y+ZN8uvWCU/ShpSWamu1PdnkAAGACc0zAW1vVoP/7/A5J0oXf/ptuv/g0ra6YkeSqADhJbdMR/Xx9jR7bUKf2YEhlRTn6ztUL9Z4zC+TzuJNdHgAASAGOCHhrqxp0x5Ob1Bm9/m5va1B3PLlJkgh5ABLKWqvKHU26f12N/rx1n9zG6JIzC3TTylJVFAdkjEl2iQAAIIU4IuDd89y2/nDXp7OnV/c8t42AByAhjnSH9FR0Geab+w4rLzNNnz5/rj60dKam57AMEwAAJIYjAt6els4TOg4AJ6vu0BE9/OJuPfpynVo7e3RGYba+edVCXVpWIL+XZZgAACCxHBHwCgPpahgizBUG0pNQDYBUY63VizsP6YHKXfrj5n0yxujiBdN104pSnT0zl2WYAABgzDgi4K1ZNW/QNXhSZKLmmlXzklgVgIlkbVWD7nlum/a0dKowkK41q+Zp1RnT9ZvqBj1QWaOte9uVm+HVv5w3Rx9eNlMFOfwBCQAAjD1jrU12DSNatGiR3bBhw6heY+AvZ5K0dNZkPfqJ5fEoD0CKO3pQkyR5XEZet1FnT1inF2TrppWlumxhIcswAQBAQhhjXrHWLhrpcY7o4EmRaZl9A1XW/Gqjfr+pUR1dIWX6HPM/AYCTNNSgplDYyuM2euzjy7Rk1mSWYQIAgHHBkbt9X7OkWB3dvXrmtcZklwJgnOrpDWtDzSHd++e3hryGV5K6esJaOjuPcAcAAMYNR7avzirJ1ZwpmXr05Vp9cHFxsssBMA6Ew1Zb9rapcnuTKncc1D93HVJHd6+Mkbxuo57eY5ezM6gJAACMN44MeMYYXbO4RF///Ra9ua9dp07LSnZJAMaYtVY1TUe0bvtBrd/RpPU7m3Soo1uSNDs/U5efNUMr5+Rr2ew8/e3NA8dcg5fudTOoCQAAjDuODHiSdPlZM/Tfz23VYy/X6UuXzk92OQDGwL62oCp3HNS67U2q3H5Qe1qDkqTp2X6dN2+KVs7J14q5ecdMwOy7fvfoKZp9xwEAAMYLxwa8/Ek+XTR/mp6qatDnL54nn4fJd0CqaT3So/U7m6Kh7qB2HOiQJAUyvFo+O0//cn6+Vs7J06z8zBGvoxs4qAkAAGC8cmzAk6QPLirW7zft1Z8279d7ywqSXQ6AUTrSHdKGmmat23FQldub9PqeVlkbWU65ZNZkXb24WCvm5Gt+QbZcLgajAACA1OPogPeOU6aoMMevR1+uJeABE1BPb1jVdS2q3N6kdTsOqqq2WT29Vl63UUVxrj5z4SlaOTdfC4sCSvM4cmgwAABwGEcHPLfL6KpFxbr3L2+pvvmIinIzkl0SgGGEw1abG9u0fkck0P1z1yEdiU66PKMwWx9dOUsr5uZrcWmuMtIc/Z83AADgUI7/DeiqRUW69y9v6Vcb6vW5i05NdjkABhg46bJyR2TaZfORHknS7CmZuvKsIq2cm6dls/MUyEhLcrUAAADJ5/iAV5SboXPm5utXG+p064WnyM11OUBS7W19e9Ll+h1vT7osyPHrgtOmaeXcPK2Yk6/pOf4kVwoAADD+OD7gSdI1i0v0qV+8qn9sP6hzT52S7HIAR2k50q0XdzZpXfQ6up3RSZe5GV4tn5OnW+bka+XcfJXmZYw46RIAAMDpCHiS3jV/qnIzvHrs5VoCHpBgR7pDermmWZXbD2rdjoN6Y0+brJUy0iKTLq9dXKIVc/N0+nQmXQIAAJwoAp4kn8etK84q0s/X1+jg4S7lT/IluyQgZXSHwtpY3xK5jm57k6rqBky6LMnVZy88VSvn5qmMSZcAAACjRsCLunpxsX76j1166tUG3fzO2ckuB5iw+iZd9l1H93LN25MuFxTm6KPnzNLKOflaxKRLAACAuOO3q6hTp2XprJKAHttQp//1jllc6wPEyFqrXQc7tG5Hkyq3H9T6nU1qiU66nDMlUx84u0gr5uRr+ew85WR4k1wtAABAaiPgDXD14mJ94deb9Gpts86eOTnZ5QDj1t7WoNZFr6Fbv6NJjdFJl4U5fr3r9LcnXU7LZtIlAADAWCLgDXBpWaG+8tvNevSfdQQ8YICWI939m4tX7mgaNOlyxZx8rYgGOiZdAgAAJBcBb4BMn0fvW1io31Tv0Z3vm68sP8vJ4ExHukP6565DqtzRpHXbD2pzY2TSZWZ00uV1S0q0fA6TLgEAAMYbAt5RPri4WI++XKffvdaoa5eUJLscYEx0h8KqrotOutxxUNV1LerptUpzu1RREtDn3nWqVszJ08LigLxuJl0CAACMVwS8o1QUB3TqtEl69OU6Ah5SVt+ky8h1dE16edchdfZEJl2eOSNHHztntlbOzdOimZOVnuZOdrkAAACIEQHvKMYYXb24RF/93WZt3dum06ZnJ7skYNSstdp5sEOV2yPX0A2cdDl36iR9cFGRVszN17JZTLoEAACYyAh4Q7i8Yob+69mteuzlOt31vjOSXQ5wUhpbO7Vue5Mqd0Q2GN/bFpl0OSOQrotOn6aVc/O1fE4eky4BAABSSFICnjGmRlK7pF5JIWvtomTUcTyTM9N00RnT9FRVg75w8Wnye1mihvGvuaNb63dGhqKs39GknQcjky4nZ6Zp+Zw8rZyTrxVz8jSTSZcAAAApK5kdvPOttQeT+P7DumZxsZ55rVH/s3mfLltYmOxy4HBrqxp0z3PbtKelU4WBdK1ZNU8XzZ+mf9Yc6l92OXDS5dLZebpuaYlWzMnXadOzmHQJAADgECzRPI6Vc/I1I5Cux1+uI+AhqdZWNeiOJzeps6dXktTQ0qnbHq+WtZKVlOZ26ayZAd32rv/f3p1H11nXeRx/f7O12Zo0adpma9pa2tJCaWu6aI8I6FBApWUErYMj4yyI4nEZKcKcM0cYN0RnHOfMGRlExPEUSwWqtcqiojIy0gXoXlsq3ZJAF9KkS9I0y3f+uE/Sm3CzNrnPzc3ndc49z3Of+8u93/t8+3vSb57f83um885phcwp00yXIiIiIiNVWAWeA8+amQP/7e4PhhRHt1JSjA9VlvPtX+/lcG0D5QVZYYckI9TXfrm7o7hr1+aQMyqN7350vma6FBEREZEOYf2Zf4m7zweuBW43s8u7NjCzW81ss5ltPnbsWPwjBG6qLMMM1mw+HMrny8h1trmVJ16q4i//6wWOnmqK2eZMUwvvuqhIxZ2IiIiIdAilwHP3mmB5FFgLLIzR5kF3r3T3yqKioniHCEBJfibvnl7ETzZX0drmocQgI8v+42f4yvpdLP76b/jCT7ZS19hMXmbsE+0l+Zlxjk5EREREEl3ch2iaWTaQ4u6ngvWrgX+Jdxx99eHKcj656mWe33uMK2eODzscSULNrW38etcRVm04xB/2HSctxVg6eyI3L57EO6YW8rMtNZ2uwQPITE9l5dIZIUYtIiIiIokojGvwJgBrg2na04BH3f3pEOLok/dcPIHC7AxWbzqkAk8GVU1dI6s3HmL1psMcPdVEaX4md1w9nQ9VljM+6t50y+eVArxlFs327SIiIiIi7eJe4Ln7a8Bl8f7cgcpIS+GDby/j4T/s59ipJopyR4UdkgxjbW3O868eY9WGQ/xm9xEcuGJ6EV9bVMGVM8eT2s3tDJbPK1VBJyIiIiK90m0S+uBDleU8+PxrPPlyFZ9499vCDkeGoTdPN7FmcxWPbjzI4dpGxuVkcNu738ZHFk7SDK0iIiIiMmhU4PXBtPE5VFaM5bFNh7n18qkEw0tFeuTubDpwglUbDvLU9jc419rGoikF3Ll0JktnTyQjTfeqExEREZHBpQKvjz68oJyVj29j04ETLJxSEHY4ksBOnm1m7cvVrNpwkL1HTpM7Oo2/WjSJmxdN4qIJuWGHJyIiIiJJTAVeH71vTjH3/nwXj206rAJPYtpRXc+qDQf52ZYaGs61Ms9qkgoAABADSURBVKcsj/s/OIf3X1ZMVoa6moiIiIgMPf2vs4+yMtL4wGUlrH2lii9dP4sxo9PDDkkSQOO5Vn6+rYZVGw6x9XAdo9NTWHZZKTcvnsScsvywwxMRERGREUYFXj+sWFDOjzceYt2WGj66uCLscCRE+46eZtWGgzzxUhUnz7YwbXwO93xgFjfMLyMvU8W/iIiIiIRDBV4/zCnLY+bEXNZsPqwCbwQ619LGs7veYNWLh/jja2+Snmpcc0kxH100iYVTCjT5joiIiIiETgVeP5gZKxaUc8/Pd7Gzpp7ZJXlhhyRxUHWigR9vPMRjm6o4frqJsrGZ3HnNDG56e7nuiygiIiIiCUUFXj8tn1fK1576E2s2HebeZSrwkslPX6nmm8/soaaukeL80bzv0mJeO3aG3+45CsBVM8dz8+IKLr+oqNsbkouIiIiIhEkFXj/lZ2VwzeyJrH2lmruvu5jR6alhhySD4KevVHP3k9tpbG4FoKbuLN/73/3kjkrl9iunsWLhJErzM0OOUkRERESkZyrwBmDFgnLWba3hmZ1vsGxuadjhSD+4O3UNzRx48wwH32zoWK7fVkNzq7+lfW5mOl+4ekYIkYqIiIiI9J8KvAFYPLWQgux0Vj6+jc+t3kJJfiYrl85g+TwVe4nA3Tl2uilSwB3vXMgdfPMMJ8+2dLQ1g5K8zJjFHcDrdWfjFbaIiIiIyAVTgTcA67bWcLKxhZa2SFFQXdfI3U9uB1CRFydtbc4bJ89GFW6R4u1AsGw419rRNjXFKB+bSUVhNvMm5VNRmM3kwiwqCrMpL8hkVFoqS+57juq6xrd8TomGZYqIiIjIMKICbwC++cyejuKuXWNzK19/ajfL5pZouvxB0tLaRk3dWQ7WBoXb8fMF3MHaBs61tHW0zUhNobwgk8mF2bxjaiGTx2V1FHIl+Zmkp6b0+Fkrl87odA0eQGZ6KiuXanimiIiIiAwfKvAGoCbGmR6AIyebmHPvs1xSkselZXnMLhnDpaV5TC7MJkWzLsZ0rqWNqhMNnYZRti8P1zZ0KqRHp6cwuTCbqUXZXDVzfEcBN6kwi+K8zAua2bL9zGv7LJoadisiIiIiw5EKvAEoyc+MOZwvPzOd6+YUs7O6nkf+70DHGaacUWnMKhkTFH6R5dSinBEz1f7Z5lYO1ca4Hq72DNUnGok+GZozKo3J47KYVTyGay+ZyOTCbCoKs5g8LpvxuaOG9Ozo8nmlKuhEREREZFhTgTcA3Q3nu+f62R0FQnNrG68eOc2O6np21NSzvbqeRzce5OwLbR3tZwVn+GaXjOHSsjymFeWQ1stQwkR1uqklMnSy0/Vwkeev13eeqCQ/K52KwmzmTxrLDfPKOq6Hm1yYRUF2hoa4ioiIiIgMkLnHnj0wkVRWVvrmzZvDDqOT6Jti93U4X0trG38+doYd1ZGCb2dNPTtrTnZMCDIqLYWLi8dwSWl74ZfH9Am5ZKQlRtFX39Ac83q4A282cPx0U6e243JGdRRuFYVZkbNwwXp+VkZI30BEREREZHgys5fcvbLXdirwwtXa5uw/Hin62gu/XTUnOdUUmco/IzWFGRNzuaQ0r6Pwmz4hd0husO7u1J4516lwi17WNTR3aj9xzOjzhdu48wVcRWE2OaN0clhEREREZLCowBvG2tqcg7UNHUXfjpp6dlSfpL4xUmClpRjTJ+SeP9NXmses4jEdRV9PZxfdnaOnmmJeD3fweENHYQmQYpHrDTuug4u6Hq58bBaZGYNfZIqIiIiIyFupwEsy7k7ViUa2R53p21Fdz4ngrFpqijGtKIcxo9PYUlXX6cbdaSnGxcW5NLc6B99s6HTtYFqKUV4QGUJZURBcCxfcYqBsbOQecSIiIiIiEq6+FngaRzdMmEUKsfKCLK67tBiIFH019WfPn+mrruf3e4/R5RZ9tLQ5u18/xRUzilgybVzUpCbZlOSPHrYTu4iIiIiISGcq8IYxM6M0P5PS/EyWzp4IwJS7fhGzbWub89AtC+IZnoiIiIiIxJlO3SSZkvzMfm0XEREREZHkoQIvyaxcOoPMLjNsZqansnLpjJAiEhERERGReNEQzSTTPltmf+/RJyIiIiIiw58KvCS0fF6pCjoRERERkRFIQzRFRERERESShAo8ERERERGRJKECT0REREREJEmowBMREREREUkSKvBERERERESShAo8ERERERGRJKECT0REREREJEmowBMREREREUkSKvBERERERESShAo8ERERERGRJKECT0REREREJEmowBMREREREUkS5u5hx9ArMzsF7Ak7jhFuHHA87CBEeUgAykH4lIPEoDyETzkIn3KQGEZKHircvai3RmnxiGQQ7HH3yrCDGMnMbLNyED7lIXzKQfiUg8SgPIRPOQifcpAYlIfONERTREREREQkSajAExERERERSRLDpcB7MOwARDlIEMpD+JSD8CkHiUF5CJ9yED7lIDEoD1GGxSQrIiIiIiIi0rvhcgZPREREREREepHQBZ6ZXWNme8xsn5ndFXY8yc7MDpjZdjPbYmabg20FZvYrM3s1WI4NtpuZ/UeQm21mNj/c6IcnM3vYzI6a2Y6obf3e52Z2S9D+VTO7JYzvMlx1k4N7zKw66AtbzOy6qNfuDnKwx8yWRm3X8WqAzKzczH5rZrvNbKeZfTbYrr4QRz3kQf0hTsxstJltNLOtQQ7uDbZPMbMNwb/rx8wsI9g+Kni+L3h9ctR7xcyN9KyHHDxiZvuj+sHcYLuOR0PEzFLN7BUzWx88Vz/oK3dPyAeQCvwZmApkAFuBWWHHlcwP4AAwrsu2+4G7gvW7gG8E69cBTwEGLAY2hB3/cHwAlwPzgR0D3edAAfBasBwbrI8N+7sNl0c3ObgHuCNG21nBsWgUMCU4RqXqeHXBOSgG5gfrucDeYF+rLyRGHtQf4pcDA3KC9XRgQ/BvfA2wItj+APDJYP1TwAPB+grgsZ5yE/b3Gw6PHnLwCHBjjPY6Hg1dLv4ReBRYHzxXP+jjI5HP4C0E9rn7a+5+DlgNLAs5ppFoGfDDYP2HwPKo7f/jES8C+WZWHEaAw5m7Pw/Udtnc332+FPiVu9e6+wngV8A1Qx99cugmB91ZBqx29yZ33w/sI3Ks0vHqArj76+7+crB+CtgNlKK+EFc95KE76g+DLPg3fTp4mh48HLgKeDzY3rUvtPeRx4H3mJnRfW6kFz3koDs6Hg0BMysD3gc8FDw31A/6LJELvFLgcNTzKnr+RSMXzoFnzewlM7s12DbB3V+HyC9/YHywXfkZOv3d58rF0Ph0MNzm4fahgSgHQy4YWjOPyF/N1RdC0iUPoP4QN8GwtC3AUSJFwZ+BOndvCZpE78+OfR28Xg8UohxckK45cPf2fvDVoB9828xGBdvUD4bGvwN3Am3B80LUD/oskQs8i7FNU34OrSXuPh+4FrjdzC7voa3yE3/d7XPlYvB9F3gbMBd4HfjXYLtyMITMLAd4Avicu5/sqWmMbcrDIImRB/WHOHL3VnefC5QROdtwcaxmwVI5GAJdc2BmlwB3AzOBBUSGXX4xaK4cDDIzez9w1N1fit4co6n6QTcSucCrAsqjnpcBNSHFMiK4e02wPAqsJfKL5Uj70MtgeTRorvwMnf7uc+VikLn7keAXfBvwPc4P6VAOhoiZpRMpKla5+5PBZvWFOIuVB/WHcLh7HfA7Itd15ZtZWvBS9P7s2NfB63lEhpwrB4MgKgfXBEOY3d2bgB+gfjCUlgDXm9kBIkO8ryJyRk/9oI8SucDbBFwUzJiTQeSiyXUhx5S0zCzbzHLb14GrgR1E9nn7zE+3AD8L1tcBHwtmj1oM1LcPpZIL1t99/gxwtZmNDYZOXR1skwHqcj3pDUT6AkRysCKYsWsKcBGwER2vLkhwrcT3gd3u/m9RL6kvxFF3eVB/iB8zKzKz/GA9E3gvkWshfwvcGDTr2hfa+8iNwHPu7nSfG+lFNzn4U9Qfm4zItV/R/UDHo0Hk7ne7e5m7TyZy/HjO3W9G/aDP0npvEg53bzGzTxPpDKnAw+6+M+SwktkEYG3kuEUa8Ki7P21mm4A1ZvZ3wCHgpqD9L4nMHLUPaAA+Hv+Qhz8z+zFwBTDOzKqALwH30Y997u61ZvZlIv+pAvgXd+/rpCEjXjc5uMIiU2A7kdllPwHg7jvNbA2wC2gBbnf31uB9dLwauCXAXwPbg+teAP4J9YV46y4PH1F/iJti4Idmlkrkj/Br3H29me0CVpvZV4BXiBTiBMsfmdk+ImcsVkDPuZFedZeD58ysiMiwvy3AbUF7HY/i54uoH/SJRQpcERERERERGe4SeYimiIiIiIiI9IMKPBERERERkSShAk9ERERERCRJqMATERERERFJEirwREREREREkoQKPBERERERkSShAk9ERBKKmX3GzHab2Qkzu6uXtn9jZv/ZzWunB/DZvzOzyv7+XDfv9YiZ3dh7SzCz28zsY4PxuSIiMrIl7I3ORURkxPoUcK277w87kHhx9wfCjkFERJKDzuCJiEjCMLMHgKnAOjP7fPvZOTMrMrMnzGxT8FgS42enmNkfg9e/3IfPutPMtpvZVjO7L+qlm8xso5ntNbN3BW07nSk0s/VmdkWwftrMvhq8z4tmNiHGZ305OKOXYmb3mdkuM9tmZt8KXr/HzO4wsxIz2xL1aDWzir58fxEREVCBJyIiCcTdbwNqgCuBE1EvfQf4trsvAD4IPBTjx78DfDdo80ZPn2Nm1wLLgUXufhlwf9TLae6+EPgc8KU+hJ0NvBi8z/PAP3T5rPuB8cDHgXzgBmC2u88BvhLd1t1r3H2uu88Fvgc84e4H+/j9RURENERTRESGhfcCs8ys/fkYM8vt0mYJkeIH4EfAN3p5vx+4ewOAu9dGvfZksHwJmNyH2M4B66N+5i+iXvtnYIO73wpgZieBs8BDZvaLqJ/rJDhD9/fAu6Lifcv3d/dTfYhPRERGEBV4IiIyHKQA73D3xuiNUQVPO+/j+1kPbZuCZSvnf0+20HnUy+io9WZ39xg/A7AJeLuZFbh7rbu3mNlC4D3ACuDTwFWdAjMrBr4PXO/u7RPFxPz+IiIiXWmIpoiIDAfPEimGADCzuTHavECkaAK4uQ/v97dmlhW8X0Ev7Q8Ac4Nr6MqBhX0JGngauA/4hZnlmlkOkOfuvyQyBLTT9zCzdGAN8EV339sl3t6+v4iIiAo8EREZFj4DVAYTk+wCbovR5rPA7Wa2Ccjr6c3c/WlgHbDZzLYAd/Ty+S8A+4HtwLeAl/sauLv/hMj1dOuAXGC9mW0Dfg98vkvzdwILgHujJlopoW/fX0REBDs/qkRERERERESGM53BExERERERSRKaZEVERJKWmV1KZEbNaE3uviiMeERERIaahmiKiIiIiIgkCQ3RFBERERERSRIq8ERERERERJKECjwREREREZEkoQJPREREREQkSajAExERERERSRL/D0YszPBO9xk/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 7)) \n",
    "\n",
    "ax.plot(chunksize[:-2], func_time[:-2], 'o-')\n",
    "ax.plot([0, 4300], [func_time[-2], func_time[-2]], '--', label=chunksize[-2])\n",
    "ax.plot([0, 4300], [func_time[-1], func_time[-1]], '--', label=chunksize[-1])\n",
    "plt.xlim([0, 4300])\n",
    "plt.legend()\n",
    "ax.set_xlabel('field_chunksize')\n",
    "ax.set_ylabel('Time spent in pset.execute() [s]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows that in this case, `field_chunksize='auto'` and `field_chunksize=False` (the two dashed lines) are roughly the same speed, but that the fastest run is for `field_chunksize=(1, 400, 400)`. \n",
    "\n",
    "Important to note is that for small chunksize (`field_chunksize=(1, 50, 50)` in this case), the execution time is actually _longer_ than for no chunking at all, due to `dask` overhead. But this again will depend on the `FieldSet` used and the specifics of your experiment.\n",
    "\n",
    "Also note that this is for a 2D application. For 3D applications, the `field_chunksize=False` will almost always be slower than `field_chunksize='auto'` or any tuple. That can be seen in the plot below, which shows the same analysis but then for a set of simulations using the full 3D CMEMS code. In this case, the `field_chunksize='auto'` is about a factor 100(!!) faster than running without chunking. But choosing too small chunksizes can make the code even slower, again highlighting that it is wise to explore which chunksize is best for your experiment before you perform it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4UAAAGuCAYAAAAnJKN8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcXGWd9/3Pr6srXQ0EEiCJWViCMihgSDBAFBVcE5eBuM2AOiAzj848uPvcjDBz36LO+IjDjAvO7YIjioIj3ggBFWUARR0VJEiQNexLEiARSNi6k16u+486nVTvlaSrT3edz/v16lfVuc6p6l/3qao+376uc51IKSFJkiRJKqaWvAuQJEmSJOXHUChJkiRJBWYolCRJkqQCMxRKkiRJUoEZCiVJkiSpwAyFkiRJklRghkJJkiRJKjBDoSRJkiQVmKFQkiRJkgqsNe8CGmHvvfdO+++/f95lSJIkSVIubrzxxj+llGbUs21ThsL999+flStX5l2GJEmSJOUiIh6sd1uHj0qSJElSgRkKJUmSJKnADIWSJEmSVGBNeU6hJEmSpOLp6upizZo1dHZ25l3KuKlUKsybN49yubzDz2EolCRJktQU1qxZw9SpU9l///2JiLzLabiUEo8//jhr1qxh/vz5O/w8Dh+VJEmS1BQ6OzvZa6+9ChEIASKCvfbaa6d7Rg2FkiRJkppGUQJhn7H4eQ2FkiRJklRghkJJkiRJysmKFSu4/fbbc63BUChJkiRJOTEUSpIkSVKTWb58OS95yUs45JBDOPfccwHYbbfdtq6/+OKLec973sNvf/tbLr/8ck477TQWLlzIvffey6pVq1iyZAkLFizgLW95C08++WTD6/WSFJIkSZKa07feNLjtkOVw5Hthy3Nw4TsGr1/4Tlj0Lnj2cfjBSf3XnfKTur7teeedx5577klHRwdHHHEEb3vb24bc7mUvexnHHXccb37zm3n7298OwIIFC/jyl7/MMcccwyc+8Qk+9alP8cUvfrGu77uj7CmUJEmSpDF0zjnncNhhh7FkyRIefvhh7r777roet2nTJjZu3MgxxxwDwMknn8yvfvWrRpYKFLSncMVNazn7ytWs29jBnGntnLb0IJYvmpt3WZIkSZLG0kg9e1N2GXn9rnvV3TNY69prr+Xqq6/md7/7HbvssgvHHnssnZ2d/S4dsbPXFRxrhespXHHTWs645BbWbuwgAWs3dnDGJbew4qa1eZcmSZIkaZLbtGkT06dPZ5ddduHOO+/kuuuuA2DWrFnccccd9Pb2cumll27dfurUqTz99NMA7LHHHkyfPp1f//rXAHz3u9/d2mvYSIULhWdfuZqOrp5+bR1dPZx95eqcKpIkSZLULJYtW0Z3dzcLFizgf/2v/8WSJUsAOOuss3jzm9/Mq1/9ambPnr11+xNOOIGzzz6bRYsWce+993L++edz2mmnsWDBAlatWsUnPvGJhtdcuOGj6zZ2bFe7JEmSJNWrra2Nn/70p0Ou65tMptbRRx896JIUfb2L46VwPYVzprVvV7skSZIkNbPChcLTlh5Ee7nUr629XOK0pQflVJEkSZIk5adww0f7Zhk98/Jb2dTRzew9Knx82QudfVSSJElSIRUuFEI1GHZ09XDGJbdwyakvY/YeDh2VJEmSVEyFGz7ap1Ku/uidXb05VyJJkiRJ+SluKGytnlfYsaVnlC0lSZIkqXkVNxROqYbCzm5DoSRJkqSxUSqVWLhw4davBx54YNhtH3jgAQ499NDxK24YhTynELb1FHZ2GQolSZIkjY329nZWrVqVdxnbpbg9hVvPKTQUSpIkSWqcBx54gFe84hUcfvjhHH744fz2t78dtM1tt93GkUceycKFC1mwYAF33303ABdccMHW9r/927+lp2fs80thewrb+4aPOtGMJEmS1JRO+dkpg9qW7r+UE154Ah3dHZx69amD1h//guNZ/oLlPNn5JB+79mP91n1r2bdG/Z4dHR0sXLgQgPnz53PppZcyc+ZMrrrqKiqVCnfffTcnnngiK1eu7Pe4r33ta3z4wx/mXe96F1u2bKGnp4c77riDiy66iN/85jeUy2VOPfVULrzwQk466aTt+TWMqrCh0OGjkiRJksbaUMNHu7q6+MAHPsCqVasolUrcddddgx730pe+lM985jOsWbOGt771rRx44IFcc8013HjjjRxxxBFANXDOnDlzzGsubigsZ7OPGgolSZKkpjRSz157a/uI66dXptfVM1iPL3zhC8yaNYubb76Z3t5eKpXKoG3e+c53ctRRR/GTn/yEpUuX8h//8R+klDj55JP57Gc/OyZ1DKew5xS2lx0+KkmSJKnxNm3axOzZs2lpaeG73/3ukOcF3nfffRxwwAF86EMf4rjjjuOPf/wjr3nNa7j44otZv349AE888QQPPvjgmNdX2FDY5kQzkiRJksbBqaeeyvnnn8+SJUu466672HXXXQdtc9FFF3HooYeycOFC7rzzTk466SQOPvhg/vmf/5nXv/71LFiwgNe97nU88sgjY15fpJTG/Enztnjx4jTwxM2BUkoc8A9X8MFXvYCPvf6gcapMkiRJUqPccccdvOhFL8q7jHE31M8dETemlBbX8/jC9hRGBJXWkucUSpIkSSq0woZCqF6r0HMKJUmSJBVZwUNhyXMKJUmSpCbSjKfHjWQsft5Ch8L2ssNHJUmSpGZRqVR4/PHHCxMMU0o8/vjjQ17iYnsU9jqFAG3lksNHJUmSpCYxb9481qxZw4YNG/IuZdxUKhXmzZu3U89R6FBYKbewudueQkmSJKkZlMtl5s+fn3cZk06hh49WWkt0bDEUSpIkSSquQofC9iklOu0plCRJklRghQ6FXpJCkiRJUtEVOxQ6fFSSJElSwRU7FE4pOdGMJEmSpEIrdihs9ZIUkiRJkoqt2KGw3OLF6yVJkiQVWqFDYXu5RE9voqvH3kJJkiRJxVToUFgplwDotLdQkiRJUkEVPBRWf3yHkEqSJEkqqoKHwmpP4WYnm5EkSZJUUIZCHD4qSZIkqbgMheBlKSRJkiQVVqFDYXsWCj2nUJIkSVJRFToU9k004/BRSZIkSUVV8FDoOYWSJEmSiq3godBLUkiSJEkqtoKHQi9JIUmSJKnYDIVAZ7c9hZIkSZKKyVAIdGwxFEqSJEkqpmKHwta+2UcdPipJkiSpmBoWCiNin4j4RUTcERG3RcSHs/Y9I+KqiLg7u52etUdEnBMR90TEHyPi8JrnOjnb/u6IOHmsamwttVAuhcNHJUmSJBVWI3sKu4H/L6X0ImAJ8P6IOBg4HbgmpXQgcE22DPAG4MDs633AV6EaIoEzgaOAI4Ez+4LkWKi0lhw+KkmSJKmwGhYKU0qPpJT+kN1/GrgDmAscD5yfbXY+sDy7fzzwnVR1HTAtImYDS4GrUkpPpJSeBK4Clo1VnZUpJTbbUyhJkiSpoMblnMKI2B9YBFwPzEopPQLV4AjMzDabCzxc87A1Wdtw7WOiUm7xnEJJkiRJhdXwUBgRuwE/BD6SUnpqpE2HaEsjtA/8Pu+LiJURsXLDhg111+fwUUmSJElF1tBQGBFlqoHwwpTSJVnzY9mwULLb9Vn7GmCfmofPA9aN0N5PSunclNLilNLiGTNm1F1j+5SSE81IkiRJKqxGzj4awDeBO1JKn69ZdTnQN4PoycBlNe0nZbOQLgE2ZcNLrwReHxHTswlmXp+1jYlKa4nOLkOhJEmSpGJqbeBzHw38FXBLRKzK2v4BOAv4QUT8DfAQ8I5s3RXAG4F7gOeAUwBSSk9ExD8BN2TbfTql9MRYFdlWbuHpzu6xejpJkiRJmlQaFgpTSv/N0OcDArxmiO0T8P5hnus84Lyxq26b9nKJDU9vbsRTS5IkSdKENy6zj05klbLDRyVJkiQVl6HQS1JIkiRJKrDCh8L2cokOewolSZIkFVThQ6HDRyVJkiQVWeFDYVu5xObuXnp7U96lSJIkSdK4K3worJSrv4LN3Z5XKEmSJKl4Ch8K28slAIeQSpIkSSqkwofCSl8o7DYUSpIkSSoeQ2E2fLRji6FQkiRJUvEUPhRuGz7qOYWSJEmSiqfwobDN4aOSJEmSCqzwobDSmoVCh49KkiRJKqDCh8L2KfYUSpIkSSquwofCvolmPKdQkiRJUhEZClu9TqEkSZKk4ip8KOwbPtphKJQkSZJUQIUPhdt6Ch0+KkmSJKl4Ch8K27aeU2hPoSRJkqTiMRS2thBhKJQkSZJUTIUPhRFBpbVkKJQkSZJUSIUPhVC9LIXnFEqSJEkqIkMhUCmXnH1UkiRJUiEZCoH2ssNHJUmSJBWToRBoK5ccPipJkiSpkAyF9J1TaE+hJEmSpOIxFOLwUUmSJEnFZSikOtFMZ7ehUJIkSVLxGAqpDh/t2GIolCRJklQ8hkKynkInmpEkSZJUQIZCqqFws8NHJUmSJBWQoRCotNpTKEmSJKmYDIVA+5QWOpx9VJIkSVIBGQqp9hT29Ca6euwtlCRJklQshkKq5xQCXqtQkiRJUuEYCoHKlGoodAipJEmSpKIxFAKV1uqvYbOTzUiSJEkqGEMhDh+VJEmSVFyGQraFQoePSpIkSSoaQyHQvrWn0OGjkiRJkorFUAhUytVfg8NHJUmSJBWNoRCHj0qSJEkqLkMhTjQjSZIkqbhah1sREXvW8fjelNLGMawnF33DR70khSRJkqSiGTYUAuuyrxhhmxKw75hWlAOHj0qSJEkqqpFC4R0ppUUjPTgibhrjenLR7vBRSZIkSQU10jmFL63j8fVsM+FVvCSFJEmSpIIaNhSmlDoBIuL5EdGW3T82Ij4UEdNqt5nsSi1BuRQOH5UkSZJUOPXMPvpDoCciXgB8E5gPfK+hVeWgUi45fFSSJElS4dQTCntTSt3AW4AvppQ+CsxubFnjr1IusbnbUChJkiSpWOoJhV0RcSJwMvDjrK3cuJLyUSm3eE6hJEmSpMKpJxSeQnVCmc+klO6PiPnABY0ta/y1l0t0bLGnUJIkSVKxjHRJCgBSSrcDH6pZvh84q5FF5aFSLtHp8FFJkiRJBTNsT2FEnDvag+vZZrKotDrRjCRJkqTiGamncHlEjHTJiQBeNcb15Kat3MJTnd15lyFJkiRJ42qkUHhaHY//9VgVkrf2cokNT2/OuwxJkiRJGlfDhsKU0vnjWUjevE6hJEmSpCKqZ/bRQqiUW+gwFEqSJEkqGENhpr1c8jqFkiRJkgrHUJhx+KgkSZKkIhrxOoURMQ84AXgFMAfoAG4FfgL8NKXUNF1rbeUSm7t76e1NtLRE3uVIkiRJ0rgY6TqF3wLOA7YAnwNOBE4FrgaWAf8dEa8cjyLHQ3u5BMDm7qbJuZIkSZI0qpF6Cv8tpXTrEO23ApdExBRg38aUNf4q5Wo+7uzqoX1KKedqJEmSJGl8DNtTOEwgrF2/JaV0z9iXlI9K1lPoDKSSJEmSimSk4aM/iog/j4jyEOsOiIhPR8RfN7a88dM3fNTJZiRJkiQVyUjDR98LfAz4YkQ8AWwAKsD+wL3Av6eULmt4heNk2/BRzymUJEmSVBzDhsKU0qPA3wN/HxH7A7Opzj56V0rpuXGpbhy19fUUdttTKEmSJKk46r1OYQJ2TSmtAlJETG1gTbnYOnx0i6FQkiRJUnGMGgoj4r3AxcDXs6Z5wIo6HndeRKyPiFtr2j4ZEWsjYlX29caadWdExD0RsToilta0L8va7omI07fnh9seFXsKJUmSJBVQPT2F7weOBp4CSCndDcys43Hfpno9w4G+kFJamH1dARARBwMnAIdkj/lKRJQiogT8b+ANwMHAidm2Y85zCiVJkiQVUT2hcHNKaUvfQkS0Uh1OOqKU0q+AJ+qs43jg+ymlzSml+4F7gCOzr3tSSvdlNXw/23bMVVqzS1I4fFSSJElSgdQTCn8ZEf8AtEfE64D/A/xoJ77nByLij9nw0ulZ21zg4Zpt1mRtw7WPub4L1jt8VJIkSVKR1BMKT6d6OYpbgL8Frkgp/eMOfr+vAs8HFgKPAP+WtccQ26YR2geJiPdFxMqIWLlhw4btLqyvp9Dho5IkSZKKpJ5Q+MGU0jdSSu9IKb09pfSNiPjwjnyzlNJjKaWelFIv8A2qw0Oh2gO4T82m84B1I7QP9dznppQWp5QWz5gxY7tra9t6TqE9hZIkSZKKo55QePIQbe/ZkW8WEbNrFt8C9M1MejlwQkS0RcR84EDg98ANwIERMT8iplCdjObyHfneo2lrbSHCUChJkiSpWIa9eH1EnAi8E5gfEbVBbCrw+GhPHBH/CRwL7B0Ra4AzgWMjYiHVIaAPUB2OSkrptoj4AXA70A28P6XUkz3PB4ArgRJwXkrptu38GesSEVRaS4ZCSZIkSYUybCgEfkv1vL+92XbuH8DTwB9He+KU0olDNH9zhO0/A3xmiPYrgCtG+35joVJuocNQKEmSJKlAhg2FKaUHgQeBl45fOflqL5ecaEaSJElSoYzUUwhARDzNthk/pwBl4NmU0u6NLCwPlbLDRyVJkiQVy6ihMKU0tXY5IpazbdbQptJmKJQkSZJUMPXMPtpPSmkF8OoG1JK79nKLw0clSZIkFUo9w0ffWrPYAixmmAvIT3YOH5UkSZJUNKOGQuDPa+53U72UxPENqSZnlXKJpzq78i5DkiRJksZNPecUnjIehUwE7eUSHVvsKZQkSZJUHKOeUxgR50fEtJrl6RFxXmPLykeb5xRKkiRJKph6JppZkFLa2LeQUnoSWNS4kvJTKZfY3G1PoSRJkqTiqCcUtkTE9L6FiNiT+s5FnHQqrQ4flSRJklQs9YS7fwN+GxEXU5119C+AzzS0qpy0T2mhs9vho5IkSZKKo56JZr4TESupXpswgLemlG5veGU5qLSW6OlNdPX0Ui5t9yUcJUmSJGnSqTf57Ak8m1L6MrAhIuY3sKbcVMolADq8VqEkSZKkgqhn9tEzgY8DZ2RNZeCCRhaVl8qUaij0AvaSJEmSiqKensK3AMcBzwKklNYBUxtZVF4qrdVfx2YvSyFJkiSpIOoJhVtSSonqJDNExK6NLSk/Dh+VJEmSVDT1hMIfRMTXgWkR8V7gauA/GltWPtrLDh+VJEmSVCz1zD76rxHxOuAp4CDgEymlqxpeWQ4qW0Ohw0clSZIkFcOooTAiDs5C4FU1bcemlK5tZGF5qJSrHacOH5UkSZJUFPUOH/37qGqPiC8Dn210YXmoOHxUkiRJUsHUEwqPAvYFfgvcAKwDjm5kUXkxFEqSJEkqmnpCYRfQAbQDFeD+lFJTnnTXN3zUS1JIkiRJKop6QuENVEPhEcDLgRMj4uKGVpWTdi9JIUmSJKlgRp1oBviblNLK7P6jwPER8VcNrCk3Dh+VJEmSVDSj9hSmlFZGxMsj4hSAiNgb+O+GV5YDL0khSZIkqWhGDYURcSbwceCMrGkKcEEji8pLqSUol8Lho5IkSZIKo55zCt8CHAc8C5BSWgdMbWRReaqUSw4flSRJklQY9YTCLSmlBCSAiNi1sSXlq1IusbnbUChJkiSpGOq9eP3XgWkR8V7gauAbjS0rP5VyCx1bDIWSJEmSimHU2UdTSv8aEa8DngIOAj6RUrqq4ZXlpL1ccqIZSZIkSYVRzyUpyEJg0wbBWpVyiU6Hj0qSJEkqiHqGjw4SEeeOdSETRaW15PBRSZIkSYWxQ6EQ+PqYVjGBVKaU6Ox2+KgkSZKkYqjnOoXvGKL5gAbUMiFUWlvY7CUpJEmSJBVEPT2FZ9TZ1hQq5ZIXr5ckSZJUGMNONBMRbwDeCMyNiHNqVu0OdDe6sLy0e/F6SZIkSQUy0uyj64CVwHHAjTXtTwMfbWRReaqUW7wkhSRJkqTCGDYUppRuBm6OiO9l2+2bUlo9bpXlxOGjkiRJkoqknnMKlwGrgJ8BRMTCiLi8oVXlqFIusaW7l97elHcpkiRJktRw9YTCTwJHAhsBUkqrgP0bV1K+KuUSAJu9LIUkSZKkAqgnFHanlDY1vJIJolKu/kqcbEaSJElSEdQTCm+NiHcCpYg4MCK+DPy2wXXlpq+n0PMKJUmSJBVBPaHwg8AhwGbgP4GngI80sqg8tWeh0J5CSZIkSUUw0iUpAEgpPQf8I/CPEVECdk0pdTa8spxsGz7qOYWSJEmSmt+oPYUR8b2I2D0idgVuA1ZHxGmNLy0fbQ4flSRJklQg9QwfPTil9BSwHLgC2Bf4q4ZWlaO+4aObDYWSJEmSCqCeUFiOiDLVUHhZSqkLaNqL+PVNNNPZbSiUJEmS1PzqCYVfBx4AdgV+FRH7UZ1spin1nVPYscVzCiVJkiQ1v3ommjkHOKem6cGIeFXjSsqXs49KkiRJKpJ6JprZKyLOiYg/RMSNEfElYI9xqC0XDh+VJEmSVCT1DB/9PrABeBvw9uz+RY0sKk+V1mz20S2GQkmSJEnNb9Tho8CeKaV/qln+54hY3qiC8laZUs3Jm7s9p1CSJElS86unp/AXEXFCRLRkX38B/KTRheVlSqmFCM8plCRJklQM9YTCvwW+B2zOvr4PfCwino6IppuFNCKotJYcPipJkiSpEOqZfXTqeBQykbRPKTnRjCRJkqRCqKensHAqrS10dnlOoSRJkqTmZygcQqVc8pxCSZIkSYVgKByCoVCSJElSUdRzSQoiogTMqt0+pfRQo4rKW6Xs8FFJkiRJxTBqKIyIDwJnAo8BfUkpAQsaWFeu7CmUJEmSVBT19BR+GDgopfR4o4uZKCrlEps6uvIuQ5IkSZIarp5zCh8GNjW6kImk3Z5CSZIkSQVRT0/hfcC1EfETqhevByCl9PmGVZWzNs8plCRJklQQ9YTCh7KvKdlX0/OcQkmSJElFMWooTCl9ajwKmUgcPipJkiSpKIYNhRHxxZTSRyLiR1RnG+0npXRcQyvLUaXcQmd3LyklIiLvciRJkiSpYUbqKfxudvuv41HIRFJpLdHTm+jqSUxpNRRKkiRJal7Dzj6aUroxu/3lUF+jPXFEnBcR6yPi1pq2PSPiqoi4O7udnrVHRJwTEfdExB8j4vCax5ycbX93RJy8cz9ufdqnlADo7HYIqSRJkqTmVs8lKXbUt4FlA9pOB65JKR0IXJMtA7wBODD7eh/wVaiGSOBM4CjgSODMviDZSG3lLBR6XqEkSZKkJtewUJhS+hXwxIDm44Hzs/vnA8tr2r+Tqq4DpkXEbGApcFVK6YmU0pPAVQwOmmOu0lr9tXRu8bIUkiRJkprbqKEwIt5RT1udZqWUHgHIbmdm7XOBh2u2W5O1Ddc+VJ3vi4iVEbFyw4YNO1helcNHJUmSJBVFPT2FZ9TZtjOGms0ljdA+uDGlc1NKi1NKi2fMmLFTxVRaHT4qSZIkqRhGuiTFG4A3AnMj4pyaVbsD3Tv4/R6LiNkppUey4aHrs/Y1wD41280D1mXtxw5ov3YHv3fdKlvPKXT4qCRJkqTmNlJP4TpgJdAJ3FjzdTnVc/12xOVA3wyiJwOX1bSflM1CugTYlA0vvRJ4fURMzyaYeX3W1lDtU6q/lg57CiVJkiQ1uWF7ClNKNwM3R8T3Ukpd2/vEEfGfVHv59o6INVRnET0L+EFE/A3wENB3buIVVHsl7wGeA07JangiIv4JuCHb7tMppYGT14y5NoePSpIkSSqIkS5e3+fIiPgksF+2fQAppXTASA9KKZ04zKrXDLFtAt4/zPOcB5xXR51jpuIlKSRJkiQVRD2h8JvAR6kOHS1ESqqUs0tSGAolSZIkNbl6QuGmlNJPG17JBNLuRDOSJEmSCqKeUPiLiDgbuATY3NeYUvpDw6rKmcNHJUmSJBVFPaHwqOx2cU1bAl499uVMDH2h0NlHJUmSJDW7UUNhSulV41HIRFJqCaaUWhw+KkmSJKnpjXSdQgAiYlZEfDMifpotH5xdUqKptZVbHD4qSZIkqemNGgqBb1O9YPycbPku4CONKmiiqJRLhkJJkiRJTa+eULh3SukHQC9ASqmbAlyaot1QKEmSJKkA6gmFz0bEXlQnlyEilgCbGlrVBFApe06hJEmSpOZXz+yjHwMuB54fEb8BZgBvb2hVE0ClXHL2UUmSJElNr57ZR/8QEccABwEBrE4pdTW8spx5TqEkSZKkIhg1FEZEBTgVeDnVIaS/joivpZQ6G11cnirlEps6mj77SpIkSSq4es4p/A5wCPBl4N+Bg4HvNrKoiaDS2sJmewolSZIkNbl6zik8KKV0WM3yLyLi5kYVNFG0T/GcQkmSJEnNr56ewpuyGUcBiIijgN80rqSJodLqOYWSJEmSml89PYVHASdFxEPZ8r7AHRFxC5BSSgsaVl2OvCSFJEmSpCKoJxQua3gVE5CXpJAkSZJUBPUMH20FHk0pPQjMB44HNqWUHszamlKlXGJLdy+9vSnvUiRJkiSpYeoJhT8EeiLiBcA3qQbD7zW0qgmgUi4BsLnbIaSSJEmSmlc9obA3pdQNvBX4Ykrpo8DsxpaVv0q5+qtxCKkkSZKkZlZPKOyKiBOBk4AfZ23lxpU0MbRnPYXOQCpJkiSpmdUTCk8BXgp8JqV0f0TMBy5obFn5qxgKJUmSJBXAqLOPppRuBz5Us3w/cFYji5oIHD4qSZIkqQjq6SkspG09hU40I0mSJKl5GQqHsXX2UXsKJUmSJDWxukNhROzayEImmr5Q6PBRSZIkSc1s1FAYES+LiNuBO7LlwyLiKw2vLGftDh+VJEmSVACjTjQDfAFYClwOkFK6OSJe2dCqdtIDTz3AKT87pV/b0v2XcsILT6Cju4NTrz510GOOf8HxLH/Bcp7sfJKPXfsxOrt6aN93I19dvRsXr2vjLw/6S5bNX8ajzz7KGb8+Y9DjTz7kZI7d51ju33Q/n/7dpwetf9+C9/HSOS/lzifu5HO//9yg9R8+/MMsnLmQVetX8aU/fGnQ+o8f+XFeuOcL+d2633HuH88dtP4TL/0E8/eYz7UPX8v5t50/aP1nX/FZnrfr8/jZ/T/jotUXDVr/+WM/z/TKdFbcs4LL7rls0PqvvPYrtLe28/07v8+VD1w5aP23ln0LgG/f+m1+ueaX/da1tbbxtdd+DYCv3fw1rn/k+n7rp7VN4wuv+gIAX7zxi9y8MtHoAAAgAElEQVS84eZ+62ftOouzXlGd2+hzv/8cdz5xZ7/1++2+H5982ScB+ORvP8mDTz3Yb/0L93whHz/y4wCc/uvTeezZx/qtP2zGYXzkJR8B4KO/+CgbN2/st/6o2Ufxd4f9HQB/d/Xfsbl7c7/1x8w7hvcc+h6AQa872P7X3kC+9nztga89X3u+9mr52vO1B772fO352qs12mtvNHUNH00pPTygqenHVLa0BAC9KeVciSRJkiQ1TqRRQk9EXAx8Hvh3YAnVy1MsTimd0PjydszixYvTypUrd+o5NnV0cdin/ov/+aYX8f+84oAxqkySJEmSGi8ibkwpLa5n23p6Cv8OeD8wF1gDLMyWm1rfdQo3d3tOoSRJkqTmVc/F6/8EvGscaplQppRaiIBOZx+VJEmS1MRGDYURMR/4ILB/7fYppeMaV1b+IoJKa4mOLYZCSZIkSc2rntlHVwDfBH4EFGosZfuUEp3dhkJJkiRJzaueUNiZUjqn4ZVMQJXWFq9TKEmSJKmp1RMKvxQRZwL/BWy9YEdK6Q8Nq2qCqJRLdHhOoSRJkqQmVk8ofDHwV8Cr2TZ8NGXLTa1SLrHZUChJkiSpidUTCt8CHJBS2tLoYiaaStnho5IkSZKaWz3XKbwZmNboQiYih49KkiRJanb19BTOAu6MiBvof05hU1+SAqC9XGJTR1feZUiSJElSw9QTCs9seBUTVKVc8uL1kiRJkpraqKEwpfTL8ShkImrznEJJkiRJTW7YUBgR/51SenlEPE11ttGtq4CUUtq94dXlrN2eQkmSJElNbqSewl0BUkpTx6mWCcfho5IkSZKa3Uizj6YR1hVCpdxCR1cPKRX+VyFJkiSpSY3UUzgzIj423MqU0ucbUM+E0l4u0ZugqycxpTXyLkeSJEmSxtxIobAE7Eb1HMJCqpRLAHR29zCltZ5LOkqSJEnS5DJSKHwkpfTpcatkAmrrC4VdPexeKedcjSRJkiSNvZG6vwrbQ9inkvUOdm7xshSSJEmSmtNIofA141bFBNU+ZdvwUUmSJElqRsOGwpTSE+NZyERUad02fFSSJEmSmpGzp4ygb6KZji2GQkmSJEnNyVA4gvYp2TmF3Z5TKEmSJKk5GQpH0ObwUUmSJElNzlA4gkrZUChJkiSpuRkKR7B19lFDoSRJkqQmZSgcwdbrFHZ5TqEkSZKk5mQoHMHW2UftKZQkSZLUpAyFI/CcQkmSJEnNzlA4glJLMKXU4vBRSZIkSU3LUDiCFTetpaunl6/98l6OPuvnrLhpbd4lSZIkSdKYMhQOY8VNaznjkltI2fLajR2cccktBkNJkiRJTcVQOIyzr1w9aIKZjq4ezr5ydU4VSZIkSdLYMxQOY93Gju1qlyRJkqTJyFA4jDnT2rerXZIkSZImI0PhME5behDt2SUp+rSXS5y29KCcKpIkSZKksdeadwET1fJFc4HquYVrN3bQEvCZ5YdubZckSZKkZpBLT2FEPBARt0TEqohYmbXtGRFXRcTd2e30rD0i4pyIuCci/hgRh49XncsXzeU3p7+a//3Ow+lNMHP3ynh9a0mSJEkaF3kOH31VSmlhSmlxtnw6cE1K6UDgmmwZ4A3AgdnX+4Cvjnehr3nRTHZra2XFKi9HIUmSJKm5TKRzCo8Hzs/unw8sr2n/Tqq6DpgWEbPHs7BKucQbDn0eP7v1UToHXKZCkiRJkiazvEJhAv4rIm6MiPdlbbNSSo8AZLczs/a5wMM1j12TtfUTEe+LiJURsXLDhg1jXvDyRXN5ZnM3V9/x2Jg/tyRJkiTlJa9QeHRK6XCqQ0PfHxGvHGHbGKItDWpI6dyU0uKU0uIZM2aMVZ1bLTlgL2ZObWPFTevG/LklSZIkKS+5hMKU0rrsdj1wKXAk8FjfsNDsdn22+Rpgn5qHzwPGPZmVWoLjF87h2tXrefLZLeP97SVJkiSpIcY9FEbErhExte8+8HrgVuBy4ORss5OBy7L7lwMnZbOQLgE29Q0zHW/HL5xLd2/iJ7fk8u0lSZIkaczlcZ3CWcClEdH3/b+XUvpZRNwA/CAi/gZ4CHhHtv0VwBuBe4DngFPGv+SqQ+bszoEzd+OyVWt595L98ipDkiRJksbMuIfClNJ9wGFDtD8OvGaI9gS8fxxKG1VEsHzRXM6+cjUPP/Ec++y5S94lSZIkSdJOmUiXpJgUjjtsDgCX3+yEM5IkSZImP0Phdtpnz104Yv/pXHrTWqqdmJIkSZI0eRkKd8DxC+dyz/pnuG3dU3mXIkmSJEk7xVC4A9704tmUS8Flq9bmXYokSZIk7RRD4Q6YvusUjvmzmVy2ah09vQ4hlSRJkjR5GQp30PJFc1j/9Gauu+/xvEuRJEmSpB1mKNxBr33RLHZra2XFTQ4hlSRJkjR5GQp3UKVcYtmhz+Ontz5KZ1dP3uVIkiRJ0g4xFO6EtyyayzObu7nmjvV5lyJJkiRJO8RQuBOWHLAXM6e2scJZSCVJkiRNUobCnVBqCY47bA7Xrl7Pxue25F2OJEmSJG03Q+FOWr5oLl09iZ/c8kjepUiSJEnSdjMU7qRD5uzOC2bu5iykkiRJkiYlQ+FOigjesmguNzzwJGuefC7vciRJkiRpuxgKx8Bxh80B4LJV63KuRJIkSZK2j6FwDOyz5y4s3m86K25aS0op73IkSZIkqW6GwjGyfNFc7l7/DLc/8lTepUiSJElS3QyFY+RNL55Na0s4hFSSJEnSpGIoHCPTd53CsQfN5LJVa+npdQipJEmSpMnBUDiGli+aw2NPbeb6+x7PuxRJkiRJqouhcAy99kWz2K2tlUu9ZqEkSZKkScJQOIYq5RLLDn0eP7v1UTq7evIuR5IkSZJGZSgcY8sXzuXpzd38/M71eZciSZIkSaMyFI6xlz5/L2ZObXMIqSRJkqRJwVA4xkotwXGHzeHa1evZ+NyWvMuRJEmSpBEZChtg+aK5dPUkrrjl0bxLkSRJkqQRGQob4JA5u/OCmbuxwiGkkiRJkiY4Q2EDRATLF87h9w88wZonn8u7HEmSJEkalqGwQY5fOBeAy1aty7kSSZIkSRqeobBB9tlzFxbvN53LVq0lpZR3OZIkSZI0JENhAx2/aC53PfYMdzzydN6lSJIkSdKQDIUN9OYXz6a1JVixyglnJEmSJE1MhsIGmr7rFI49aAaXr1pHT69DSCVJkiRNPIbCBlu+aC6PPtXJ9fc/nncpkiRJkjSIobDBXvuiWezW1uo1CyVJkiRNSIbCBquUSyw79Hn89JZH6ezqybscSZIkSerHUDgOli+cy9Obu/n5nevzLkWSJEmS+jEUjoOXPn8vZk5tcwipJEmSpAnHUDgOSi3BcYfN4drVG9j43Ja8y5EkSZKkrQyF42SPXcps6ell4aev4uizfm6voSRJkqQJwVA4DlbctJav/OKerctrN3ZwxiW3GAwlSZIk5c5QOA7OvnI1HV29/do6unr4lyvvzKkiSZIkSaoyFI6DdRs7hmnv5Nxf3cuTz3qeoSRJkqR8GArHwZxp7UO2Tym18P9fcSdHffYaPnrRKm588AlSSuNcnSRJkqQiMxSOg9OWHkR7udSvrb1c4l/evoCffeQV/OXifbjq9sd421d/xxu+9Gu+e92DPLO5O6dqJUmSJBVJNGPP1OLFi9PKlSvzLqOfFTet5ewrV7NuYwdzprVz2tKDWL5o7tb1z27u5rJV67jguge5/ZGn2HVKieMXzeXdR+3HwXN2z7FySZIkSZNNRNyYUlpc17aGwoklpcSqhzdy4fUP8aOb17G5u5dF+07j3Uftx5sWzKYyoMdRkiRJkgYyFE7iUFhr43Nb+OEf1nLh9Q9y34ZnmbZLmbcfPo93HrUvB8zYLe/yJEmSJE1QhsImCYV9Ukr87t7HufD6h7jytkfp7k0c/YK9eNdR+/G6g2dRLnlqqCRJkqRtDIVNFgprrX+qkx+sfJj//P3DrN3YwcypbZxwxD6ccOS+w85yKkmSJKlYDIVNHAr79PQmfnHnei68/kGuvWsDAbz6hbN415J9OebAGbS0RN4lSpIkScrJ9oTC1kYXo8YotQSvPXgWrz14Fg8/8Rz/+fuH+MHKh7n6jsfYZ892TjxyX/5i8T7svVtb3qVKkiRJmsDsKWwiW7p7ufK2R7nguge5/v4nKJeCZYfO5t1H7cuR8/ckwt5DSZIkqQgcPlrQUFjrnvVPc8F1D/HDP6zh6c5uDpy5G+86al/e+pJ57F4p512eJEmSpAYyFBoKt+rY0sOPbl7Hhdc/yM1rNtFeLnHcYXN495L9ePG8PfIuT5IkSVIDGAoNhUO6Zc0mLrjuQS67eS2dXb0smLcH7z5qPwC+dM3drNvYwZxp7Zy29CCWL5qbc7WSJEmSdpSh0FA4ok0dXVz6hzVceP1D3L3+mUHrK+UWznrrAoOhJEmSNEkZCg2FdUkpccRnruZPz2wZcv2MqW3s0V4e9LV7e5ndK63923fZdr+9XHJSG0mSJE0oK25ay9lXri7M6DgvSaG6RASPDxMIAV7zwpls6ujiqc4u1j/dyd3rn2bTc108vbmbkf6XUC7F1vC4R3uZ3SuDg2Xt+r5QuXulld3aWg2UksZM0Q4AJElDW3HTWs645BY6unoAWLuxgzMuuQXAvwsYCgtvzrR21m7sGNQ+d1o7Z71twZCP6elNPNPZzaaOrkFfT3UOWO7o4snntvDA489uXe4dIVCWWqJfL+TuQ4TJodp3by8zta2VlhYDJXggLIEHABOFn0fS6FJKpAS9KdGbIDFgueZ2tO36lmvboXrbmxK9vf0fN+R2vYOfv/Y2UX2eoZ5/uO22LteuTyM8boTttv4+araj5ucYartL/rB269+DPh1dPZx95Wo/kzAUFt5pSw/qd9AE0F4ucdrSg4Z9TKklqsNFd9n+S1v09iae2dLNpue2BcinhgiXmzq2hc41T3Zs3aZ7hETZEjB1iF7J3dtbhw2XfV9TK2VK2xkoJ+qBjgfCE8NEfX1MFiklenoT3b3bbrt7ege19fT2Zutqlnuq6//px7cPeQDwTz++nWm7lGmJoCWCCIigeh9oaQlaojqaImDo7aKvHaB6W9se2eNr22PAdn3ra7eDbc9bu36yasbPo9qD020HntlBKNltb/+D1W0HrDXb9Q518Nt3YDzaQfLg7YY6SK5nu1Tzcwx7cN33c/f2/YyDw0lvNoSoNkwMGWJqfv7hQsdo4We425HCyfaGjqG2Gy10DAxvQz1uuO20faqfqds+nwm2fo72fZZHVD/P+7Z7bkvPkM+1bojOkSLynEJNmoPXlBLPbekZuodymPt9AfOpji629PSO+PxTK61DD3XdZXDv5KqHn+Qrv7iXzd3bnrNSbuETbz6YpYc8j5607T9jPb2p321vqva29m+v2bY3VR+f/aeup285267vfl979eB5W/vnr1rNpo7uQT/fHu2tfPS1f1b9gOw72BziYLalpW+5Zh39D1xbWvqWhznobak9qB36gLn/ge8QdQw4CB94oDzaNnkaeCAM1X+2fPatL97u91bfAUN3Fna2BaIsIPVsW+7ut1zT3jM4RPU9tnvAcu12XcM9rjfR05Po6u3ttzwwpPWrc0C46x/stgW+rr7XtkdJ/Qx8z2w7GMqWtyOsDgqvDHw/9T3vts+C6jbDB97htvvvu/9EZ/fgz95KawtHv2DvmgPlgeEkO1CuCVVDbzcwPPQdyA8d2gY9rm+73r7gM0rY82W53XbktVvPdv3+kTPcP2IG/B0b6m/GUNsNfA9sCx0jb7etzprH9fs5+v99Hfze2dZe+7dxyPdYS20Ngx9X/QfXwH82Db/dwH+CtQQ1P/O2n2ukf5Zt73YD92FLBNEy1GfS4LC3I3/njz7r58OOjvvN6a/e7uebDJxoxlCoAVJKdHb1jhgmh++17OoX/jQ5DAyOw4XSgQF44AHIcGF16HBbXf7jmk1D/hOiXAoOet7UweGuZrlriN6wPLUEtJZaaG0JSi2R3W5bLpf62luqt6Vt27W2tPRbHrjdwOdqbQlKpaDct022XLvdsI9rCcqlln7L7//eH4acSGvGbm187a9eQu0wqdqD/sEBo2b4E7X/4a/pBaJ2WNZwvTQDe0iG6aXpt83O9bYwoPZhh2IN+zMO6AXa+jMO3dvRP2jB7Y88Nexr69C5uw9xcD1EKBh0YDv8diOFia3/CKtnu9r3ec12jPBZMVw4Gfi47QkxQ/Um17PdiL3QLdv+gTfsdi2jP//Az8jh9o80UYzlP2wnCyeakQaICNqnlGifUuJ5e1S2+/GdXT081bktNL7tq78bdttPHXcILS1BKYJSdvBRvR9b2/v+y1fb3hJU12UHtH1/XPvul7a2s3W59ralpfr4N3/5v3lkU+egumbvUeGKD72i30Flv6FGvf3/2z7wgHDwAefgbQYOkRp4gNp3YDzwgHq4bYY6gN32Pbd/m4EH/4N+jprzIwaep7G1rt7+32Nbj8S27zlcr3RXT2LW1MqQAaa1JlwNFYa2Bq1BAa3vsS2UByyP+riWlq3fp3XAcqnvtdoyeQ/q/uebDh7yAOAf3/QiXrLf9BwrK46R/jP/4w++IoeKJBVVX/CbDKPj8mAolOpQKZeolEvMnFoNlHNHmKDn5JftP87V9ffxZS8c8kD448teyPRdp+RYWXGMdCD8zfcckUNFxeQBQP525Lx1SWqU5Yvm+jdgGM05fHT/PdLKM1/ev/GQ5XDke2HLc3DhOwY/aOE7YdG74NnH4QcnDV5/xF/DoW+DTWvgkr8dvP5lH4CD3gB/uht+9JHB61/5P+D5r4JH/gg/O2Pw+td8AvY9Ch66Hq759OD1yz4LsxfAvb+AX/3r4PV//kXY+0BY/VP47b8PXv/Wr8Me8+DWH8IN5w1e/xffgV33gpsuhFXfG7z+Xf8HpuwCv/8G3LZi8PpTflK9/c05cNeV/deVK/DuH1bv//Jf4L5f9l+/y3T4ywuq96/+JDx8Q//1u8+Bt32jev+np8Ojt/Rfv9fz4bhzqvcv/xA8fm//9c97MbzhrOr9H74XnlrXf/0+R8BrP1m9f9G74bkn+68/4Bg45u+r9y94G3R18qdnNnPfn56lNyWu6VnEN3reTHu5xK9n/it779bW//E5vPb+9MxmHnriObb09PKfbX/Jq974Fyx/3uO+9sbptdf3+ri1Z18+3V3dp+dM+Qovm7G5/+tjB157/fzZUjj6Q9X733oTg/i5V7jX3lYN+NzrZztee3/6+nFbP4+mlFrYd89d2Pvlp/ja87Xn556vvf7rfe2N+Wsv/vqK5hs+GhHLgC8BJeA/Ukpn5VySCqzvwP6hJ54jeqo9QKctPYi9V7WN8sjxsfdubVtrXPLKhfD8ufDI4zlXVRx9v/sHNrYS3dVLvyzaexp7lzbmXJk0/mo/jyRJE9Ok6CmMiBJwF/A6YA1wA3BiSun2obZ3ohlJkiRJRbY9E820NLqYMXIkcE9K6b6U0hbg+8DxOdckSZIkSZPeZAmFc4GHa5bXZG2SJEmSpJ0wWULhUHOi9xv3GhHvi4iVEbFyw4YN41SWJEmSJE1ukyUUrgH2qVmeB/SbUiildG5KaXFKafGMGTPGtThJkiRJmqwmSyi8ATgwIuZHxBTgBODynGuSJEmSpElvUlySIqXUHREfAK6kekmK81JKt+VcliRJkiRNepMiFAKklK4Arsi7DkmSJElqJpNl+KgkSZIkqQEMhZIkSZJUYIZCSZIkSSowQ6EkSZIkFZihUJIkSZIKzFAoSZIkSQVmKJQkSZKkAjMUSpIkSVKBRUop7xrGXEQ8DazOuw6xN/CnvIsoOPdB/twHE4P7IX/ug/y5DyYG90P+irIP9kspzahnw9ZGV5KT1SmlxXkXUXQRsdL9kC/3Qf7cBxOD+yF/7oP8uQ8mBvdD/twHgzl8VJIkSZIKzFAoSZIkSQXWrKHw3LwLEOB+mAjcB/lzH0wM7of8uQ/y5z6YGNwP+XMfDNCUE81IkiRJkurTrD2FkiRJkqQ6NF0ojIhlEbE6Iu6JiNPzrqeZRcQDEXFLRKyKiJVZ254RcVVE3J3dTs/aIyLOyfbLHyPi8Hyrn7wi4ryIWB8Rt9a0bffvPSJOzra/OyJOzuNnmayG2QefjIi12fthVUS8sWbdGdk+WB0RS2va/bzaQRGxT0T8IiLuiIjbIuLDWbvvhXEywj7wvTCOIqISEb+PiJuz/fCprH1+RFyfva4viogpWXtbtnxPtn7/mucacv9oZCPsg29HxP0174WFWbufRw0SEaWIuCkifpwt+z6oV0qpab6AEnAvcAAwBbgZODjvupr1C3gA2HtA278Ap2f3Twc+l91/I/BTIIAlwPV51z9Zv4BXAocDt+7o7x3YE7gvu52e3Z+e9882Wb6G2QefBP7HENsenH0WtQHzs8+okp9XO70PZgOHZ/enAndlv2vfC/nvA98L47sfAtgtu18Grs9e4z8ATsjavwb8v9n9U4GvZfdPAC4aaf/k/fNNhq8R9sG3gbcPsb2fR43bFx8Dvgf8OFv2fVDnV7P1FB4J3JNSui+ltAX4PnB8zjUVzfHA+dn984HlNe3fSVXXAdMiYnYeBU52KaVfAU8MaN7e3/tS4KqU0hMppSeBq4Blja++OQyzD4ZzPPD9lNLmlNL9wD1UP6v8vNoJKaVHUkp/yO4/DdwBzMX3wrgZYR8Mx/dCA2Sv6WeyxXL2lYBXAxdn7QPfC33vkYuB10REMPz+0ShG2AfD8fOoASJiHvAm4D+y5cD3Qd2aLRTOBR6uWV7DyH+gtHMS8F8RcWNEvC9rm5VSegSqBwzAzKzdfdNY2/t7d380xgeyoUDn9Q1bxH3QcNmwn0VU/zvveyEHA/YB+F4YV9mQuVXAeqpB4l5gY0qpO9uk9ne69fedrd8E7IX7YacM3Acppb73wmey98IXIqIta/O90BhfBP4e6M2W98L3Qd2aLRTGEG1Or9o4R6eUDgfeALw/Il45wrbum3wM93t3f4y9rwLPBxYCjwD/lrW7DxooInYDfgh8JKX01EibDtHmfhgDQ+wD3wvjLKXUk1JaCMyj2qvxoqE2y27dDw0wcB9ExKHAGcALgSOoDgn9eLa5+2CMRcSbgfUppRtrm4fY1PfBMJotFK4B9qlZngesy6mWppdSWpfdrgcupfqH6LG+YaHZ7fpsc/dNY23v7939McZSSo9lBwW9wDfYNtzEfdAgEVGmGkYuTCldkjX7XhhHQ+0D3wv5SSltBK6lep7atIhozVbV/k63/r6z9XtQHQ7vfhgDNftgWTbEOqWUNgPfwvdCIx0NHBcRD1Adgv5qqj2Hvg/q1Gyh8AbgwGymoSlUTxy9POeamlJE7BoRU/vuA68HbqX6++6bLetk4LLs/uXASdmMW0uATX1DvDQmtvf3fiXw+oiYng3ten3Wph004BzZt1B9P0B1H5yQzXQ2HzgQ+D1+Xu2U7NyPbwJ3pJQ+X7PK98I4GW4f+F4YXxExIyKmZffbgddSPb/zF8Dbs80Gvhf63iNvB36eUkoMv380imH2wZ01/6AKquey1b4X/DwaQymlM1JK81JK+1P9DPl5Suld+D6oW+vom0weKaXuiPgA1TdQCTgvpXRbzmU1q1nApdXPOVqB76WUfhYRNwA/iIi/AR4C3pFtfwXV2bbuAZ7j/7Z3byFWVXEcx78/K4jKjMhCIbo8FpWQF0oEuz34EkUFRhQZXSTDEiQjCCsNJhFCCBIyKnwoDHsQDekhKpAsu5jSBBJYENKTQkllKf8ezhbPTDPOUWrEs78fOJwzs/57XWazz5n/rLXXwILx73J/SPIOMBe4KMnPwHJggBP4uVfV/iQr6PwyBvBiVfW6cUrrjXIO5qaz3XjR2Zn3MYCq+i7JBmAQOAwsqqojTT2+X5282cD9wO7mPh6AZ/FaGE+jnYN7vRbG1RTg7SRn0Plj/4aq2pxkEHg3yUrgGzoJPM3z+iQ/0JkZmQ/HPz8a02jn4KMkk+ksSdwJLGzifT8aP8vwOuhJOkmxJEmSJKmN+m35qCRJkiTpBJgUSpIkSVKLmRRKkiRJUouZFEqSJElSi5kUSpIkSVKLmRRKkiRJUouZFEqSTntJFif5PsmBJM+MEftgkldHKTt4Em1/nGT6iR43Sl1vJbl77EhIsjDJA/9Fu5Kkduurf14vSWqtx4F5VbX3VHdkvFTV2lPdB0lSf3CmUJJ0WkuyFrgS2JRkydFZwCSTk2xMsqN5zB7h2CuSfNaUr+ihraeT7E7ybZKBrqJ7knyRZE+SOU3skBnJJJuTzG1eH0zyUlPP9iSXjNDWimbmcEKSgSSDSXYlWd2UP59kaZKpSXZ2PY4kuayX8UuSBCaFkqTTXFUtBPYBNwEHuorWAK9U1QzgLmDdCIevAV5rYn45XjtJ5gF3ALOq6jpgVVfxmVU1E3gKWN5Dt88Ftjf1fAo8MqytVcDFwALgAuBO4OqquhZY2R1bVfuqalpVTQNeBzZW1U89jl+SJJePSpL61q3AVUmOfn1+konDYmbTSZgA1gMvj1Hfm1X1O0BV7e8qe795/gq4vIe+/QVs7jrmtq6y54DPq+pRgCS/An8C65Js6TpuiGYm8GFgTld//zX+qvqth/5JklrEpFCS1K8mADdU1R/d3+xKko6qHuvLcWIPNc9HOPbZepihK3LO7nr9d1XVCMcA7ACuT3JhVe2vqsNJZgK3APOBJ4Cbh3QsmQK8AdxeVUc3yxlx/JIkDefyUUlSv/qQTgIFQJJpI8Rso5NoAdzXQ30PJTmnqe/CMeJ/BKY19wReCszspdPAVmAA2JJkYpLzgElV9QGd5alDxpHkLGADsKyq9gzr71jjlyTJpFCS1LcWA9ObzVkGgYUjxDwJLEqyA5h0vMqqaiuwCfgyyU5g6RjtbwP2AruB1cDXvXa8qt6jc3/gJmAisDnJLuATYMmw8BuBGf+1PCoAAABsSURBVMALXZvNTKW38UuSRI6tXpEkSZIktY0zhZIkSZLUYm40I0lSlyTX0NmJtNuhqpp1KvojSdL/zeWjkiRJktRiLh+VJEmSpBYzKZQkSZKkFjMplCRJkqQWMymUJEmSpBYzKZQkSZKkFvsHju0Dwu7QXawAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 7)) \n",
    "\n",
    "ax.plot(chunksize[:-2], func_time3D[:-2], 'o-')\n",
    "ax.plot([0, 4300], [func_time3D[-2], func_time3D[-2]], '--', label=chunksize[-2])\n",
    "ax.plot([0, 4300], [func_time3D[-1], func_time3D[-1]], '--', label=chunksize[-1])\n",
    "plt.xlim([0, 4300])\n",
    "plt.legend()\n",
    "ax.set_xlabel('field_chunksize')\n",
    "ax.set_ylabel('Time spent in pset.execute() [s]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Future developments: load-balancing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The current implementation of MPI parallelisations in Parcels is still fairly rudimentary. In particular, we will continue to develop the load-balancing of the `ParticleSet`.\n",
    "\n",
    "With load-balancing we mean that `Particles` that are close together are ideally on the same MPI processor. Practically, it means that we need to take care how `Particles` are spread over chunks and processors. See for example the two figures below:\n",
    "\n",
    "![](http://oceanparcels.org/images/parcelsParallel.png)\n",
    "*Example of load-balancing for Particles. The domain is chunked along the thick lines, and the orange and blue particles are on separate MPI processors. Before load-balancing (left panel), two chuncks in the centre of the domain have both orange and blue particles. After the load-balancing (right panel), the Particles are redistributed over the processors so that the number of chunks and particles per processor is optimised.*\n",
    "\n",
    "The difficulty is that since we don't know how the `ParticleSet` will disperse over time, we need to do this load-balancing 'on the fly'. If you to contribute to the optimisation of the load-balancing, please leave a message on [github](https://github.com/OceanParcels/parcels/issues)!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}