{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial on the Analytical Advection kernel in Parcels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While [Lagrangian Ocean Analysis](https://www.sciencedirect.com/science/article/pii/S1463500317301853) has been around since at least the 1980s, the [Blanke and Raynaud (1997)](https://journals.ametsoc.org/doi/full/10.1175/1520-0485%281997%29027%3C1038%3AKOTPEU%3E2.0.CO%3B2) paper has really spurred the use of Lagrangian particles for large-scale simulations. In their 1997 paper, Blanke and Raynaud introduce the so-called *Analytical Advection* scheme for pathway integration. This scheme has been the base for the [Ariane](http://stockage.univ-brest.fr/~grima/Ariane/) and [TRACMASS](http://www.tracmass.org/) tools. We have also implemented it in Parcels, particularly to facilitate comparison with for example the Runge-Kutta integration scheme.\n", "\n", "In this tutorial, we will briefly explain what the scheme is and how it can be used in Parcels. For more information, see for example [Döös et al (2017)](https://www.geosci-model-dev.net/10/1733/2017/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most advection schemes, including for example Runge-Kutta schemes, calculate particle trajectories by integrating the velocity field through time-stepping. The Analytical Advection scheme, however, does not use time-stepping. Instead, the trajectory within a grid cell is analytically computed assuming that the velocities change linearly between grid cells. This yields Ordinary Differential Equations for the time is takes to cross a grid cell in each direction. By solving these equations, we can compute the trajectory of a particle within a grid cell, from one face to another. See [Figure 2 of Van Sebille et al (2018)](https://www.sciencedirect.com/science/article/pii/S1463500317301853#fig0002) for a schematic comparing the Analytical Advection scheme to the fourth order Runge-Kutta scheme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the Analytical scheme works with a few limitations:\n", "1. The velocity field should be defined on a C-grid (see also the [Parcels NEMO tutorial](https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb)).\n", "\n", "And specifically for the implementation in Parcels\n", "2. The `AdvectionAnalytical` kernel only works for `Scipy Particles`.\n", "3. Since Analytical Advection does not use timestepping, the `dt` parameter in `pset.execute()` should be set to `np.inf`. For backward-in-time simulations, it should be set to `-np.inf`.\n", "4. For time-varying fields, only the 'intermediate timesteps' scheme ([section 2.3 of Döös et al 2017](https://www.geosci-model-dev.net/10/1733/2017/gmd-10-1733-2017.pdf)) is implemented. While there is also a way to also analytically solve the time-evolving fields ([section 2.4 of Döös et al 2017](https://www.geosci-model-dev.net/10/1733/2017/gmd-10-1733-2017.pdf)), this is not yet implemented in Parcels. \n", "\n", "We welcome contributions to the further development of this algorithm and in particular the analytical time-varying case. See [here](https://github.com/OceanParcels/parcels/blob/master/parcels/kernels/advection.py) for the code of the `AdvectionAnalytical` kernel." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we will show how this `AdvectionAnalytical` kernel performs on one idealised time-constant flow and two idealised time-varying flows: a radial rotation, the time-varying double-gyre as implemented in e.g. [Froyland and Padberg (2009)](https://www.sciencedirect.com/science/article/abs/pii/S0167278909000803) and the Bickley Jet as implemented in e.g. [Hadjighasem et al (2017)](https://aip.scitation.org/doi/10.1063/1.4982720).\n", "\n", "First import the relevant modules." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "from parcels import FieldSet, ParticleSet, ScipyParticle, JITParticle, Variable\n", "from parcels import AdvectionAnalytical, AdvectionRK4, plotTrajectoriesFile\n", "import numpy as np\n", "from datetime import timedelta as delta\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Radial rotation example\n", "\n", "As in [Figure 4a of Lange and Van Sebille (2017)](https://doi.org/10.5194/gmd-10-4175-2017), we define a circular flow with period 24 hours, on a C-grid" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def radialrotation_fieldset(xdim=201, ydim=201):\n", " # Coordinates of the test fieldset (on C-grid in m)\n", " a = b = 20000 # domain size\n", " lon = np.linspace(-a/2, a/2, xdim, dtype=np.float32)\n", " lat = np.linspace(-b/2, b/2, ydim, dtype=np.float32)\n", " dx, dy = lon[2]-lon[1], lat[2]-lat[1]\n", "\n", " # Define arrays R (radius), U (zonal velocity) and V (meridional velocity)\n", " U = np.zeros((lat.size, lon.size), dtype=np.float32)\n", " V = np.zeros((lat.size, lon.size), dtype=np.float32)\n", " R = np.zeros((lat.size, lon.size), dtype=np.float32)\n", "\n", " def calc_r_phi(ln, lt):\n", " return np.sqrt(ln**2 + lt**2), np.arctan2(ln, lt)\n", "\n", " omega = 2 * np.pi / delta(days=1).total_seconds()\n", " for i in range(lon.size):\n", " for j in range(lat.size):\n", " r, phi = calc_r_phi(lon[i], lat[j])\n", " R[j, i] = r\n", " r, phi = calc_r_phi(lon[i]-dx/2, lat[j])\n", " V[j, i] = -omega * r * np.sin(phi)\n", " r, phi = calc_r_phi(lon[i], lat[j]-dy/2)\n", " U[j, i] = omega * r * np.cos(phi)\n", "\n", " data = {'U': U, 'V': V, 'R': R}\n", " dimensions = {'lon': lon, 'lat': lat}\n", " fieldset = FieldSet.from_data(data, dimensions, mesh='flat')\n", " fieldset.U.interp_method = 'cgrid_velocity'\n", " fieldset.V.interp_method = 'cgrid_velocity'\n", " return fieldset\n", "\n", "fieldsetRR = radialrotation_fieldset()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now simulate a set of particles on this fieldset, using the `AdvectionAnalytical` kernel. Keep track of how the radius of the Particle trajectory changes during the run." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: Particle initialisation from field can be very slow as it is computed in scipy mode.\n" ] } ], "source": [ "def UpdateR(particle, fieldset, time):\n", " particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]\n", "\n", "class MyParticle(ScipyParticle):\n", " radius = Variable('radius', dtype=np.float32, initial=0.)\n", " radius_start = Variable('radius_start', dtype=np.float32, initial=fieldsetRR.R)\n", "\n", "pset = ParticleSet(fieldsetRR, pclass=MyParticle, lon=0, lat=4e3, time=0)\n", "\n", "output = pset.ParticleFile(name='radialAnalytical.nc', outputdt=delta(hours=1))\n", "pset.execute(pset.Kernel(UpdateR) + AdvectionAnalytical,\n", " runtime=delta(hours=24),\n", " dt=np.inf, # needs to be set to np.inf for Analytical Advection\n", " output_file=output)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the trajectory and calculate how much the radius has changed during the run." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV5fn48c91MskgE8IIBEKCylTCkqFB0apFcdZVi3611FXt+tb1tfqrWu23/bbVOlocVasVrYpSZwGJLIMaZK+EQCAJAhmMAJnn+v3xPEkjhpFxck6S6/16nVfOuZ91nZPkuc59389z36KqGGOMMQAefwdgjDEmcFhSMMYY08CSgjHGmAaWFIwxxjSwpGCMMaaBJQVjjDENLCmYTkdE7hWR505gvRdF5OH2iOkox/9QRGb46/iN4jihz8t0DWL3KZj2JiLbgCSgDjgIfAD8WFUrWrCvTOAVVU1uwbYvAoWq+j8t2FaBdFXNa+62bUlEHgTSVPX7/ozDdB5WUzD+cqGqRgGjgDFAS07MwW0eVRsJ5Nga6yhxmvZjScH4laoWAR8CwwBE5AYR2SAiB0QkX0R+VL+uiGSKSKGI3CUiXwOvudv2EZEK99FHRB4UkVcabTdJRJaJyF4R2SEi1zcVi4hME5GV7nrLRGTEUdZb5D5d5R7zyiZi+5uIxInIeyKyR0TK3efJjfaTJSI3NXr9X+57LxeRj0UkpdGyoSIyT0TKRGSX2+RzHnAvcKUbxyp33T4iMtddN09EfthoPw+KyJsi8oqI7Aeub+LzGt/o81rl1sbql13v/l4OiMhWEbn26L9d0xFZUjB+JSL9gAuAr9yi3cA0oDtwA/BHERnVaJNeQDyQAvwAOB8oVtUo91F8xP774ySOPwM9gFOBlU3EMQp4AfgRkAD8FZgrImFHrquqZ7hPR7rHfL2J2Gbi/H/9zX3dHzgMPHmUz+FinBP8pW6ci3GSHiISDcwHPgL6AGnAAlX9CPgN8Lobx0h3d68Bhe66lwO/EZGzGx1uOvAmEAu8ekQcfYH3gYfd9/IL4C0R6SEikcATwPmqGg1MaOqzNB2bJQXjL++IyF5gCfApzskNVX1fVbeo41Pg38DkRtt5gQdUtUpVD5/Aca4F5qvqa6pao6qlqtrUieyHwF9Vdbmq1qnqS0AVML4Z7+kbsbnHektVD6nqAeAR4MyjbPsj4FFV3aCqtTifx6lubWEa8LWq/p+qVqrqAVVd3tRO3CQ7CbjLXXcl8BxwXaPVPlPVd1TV28Rn+H3gA1X9wF0+D/gSJ3HXv8dhItJNVXeq6rpmfD6mA7CkYPzlYlWNVdUUVb21/uQkIueLSLbb9LEX52SU2Gi7Papa2Yzj9AO2nMB6KcDP3SaTve6x++F82z5R34hNRCJE5K8iUuA21SwCYkUk6CjHf7zRscsAAfo24z3gxlvmJqF6Be5+6u04xvYpwBVHfA6TgN6qehC4ErgZ2Cki74vIyScYl+kgLCmYgOE21bwF/B5IUtVYnCuTpNFqR14ud7zL53YAg07g8DuAR9xEVf+IUNXXTjD8pmL5OXASME5VuwP1zU7Ct+0AfnTE8bup6rLjvIcjj1kMxLtNTvX6A0XH2ObIOP5+RByRqvoYgKp+rKrnAL2BjcCzx9iX6YAsKZhAEgqEAXuAWhE5Hzj3ONvsAhJEJOYoy18FporI90QkWEQSROTUJtZ7FrhZRMaJI1JEvnvEyfXI46YeJ7ZonH6EvSISDzxwjHX/AtwjIkMBRCRGRK5wl70H9BKRn4hImIhEi8i4RnEMEBEPgKruAJYBj4pIuNtZfiNH9B0cwyvAhSLyHREJcveRKSLJIpIkIhe5fQtVQAXOZcWmE7GkYAKG2+RxB/AGUA5cA8w9zjYbcTpW893mjj5HLN+O0wT1c5wmmZXAyCb28yVOv8KT7rHzgOuPcegHgZfcY37vKOv8CegGlADZOB3FR3sfc4DfArPdpqa1OJ3o9Z/LOcCFwNdALjDF3fSf7s9SEVnhPr8aGIBTa5iD088x7xjvpXEcO3A6ou/FSc47gP/GOVd4cD7HYpzP8kzg1hPZr+k47OY1Y/zEvbT1OVV92d+xGFPPagrG+IGIROA0P231dyzGNGZJwZh2JiI9cZqBPsW5JNeYgGHNR8YYYxpYTcEYY0yDDj8YVmJiog4YMKBF2x48eJDIyMi2DagNWFzNY3E1j8XVPJ01rpycnBJV7fGtBaraoR8ZGRnaUgsXLmzxtr5kcTWPxdU8FlfzdNa4gC+1iXOqNR8ZY4xpYEnBGGNMA0sKxhhjGlhSMMYY08CSgjHGmAZtlhTcERW/EpH33NcDRWS5iOSKyOsiEuqWh7mv89zlAxrt4x63fJOIfKetYjPGGHNi2rKmcCewodHr3wJ/VNV0nFEnb3TLbwTKVTUN+KO7HiIyBLgKGAqcBzx9lMlIjAloOQXl3DdnDS+trSSnoNzf4RjTLG1y85o4k5F/F2e6wZ+JiABn4Qx9DPASzlDDz+AMy/ugW/4m8KS7/nRgtqpWAVtFJA8YC3zWFjEa4wt1XqXsYDWlB6soOVDNF9vKeHJhHnVeZ/iYRX9dxv3ThjBhUCIJkaHERoQS5Glqjh1jAkObjH0kIm8Cj+JMKvILnHHos93aQP28sR+q6jARWQucp6qF7rItwDicRJGtqq+45c+727zZxPFm4kyMTlJSUsbs2bNbFHdFRQVRUVEt2taXLK7maeu4ar3KgWplX5Wy3/15oFrZV63sb1S2v1o5UH38qd8aEyAqFLqHCt1DhWj3Z/cwITrE+dm4vFswON+Z2k5X+T22lc4a15QpU3JUdfSR5a2uKYjINGC3quaISGZ9cROr6nGWHWubbxaqzgJmAYwePVozMzObWu24srKyaOm2vmRxNc+JxFVZU0dJRRUlFdWUHKhyn7uvj3i+91BNk/sID/GQGBVGYlQYJ/cMo0d0aMPrxKgwEqJC2b2/kl/8czXVdV4AQoKE//nuEBKiQimtqKb0YDWlFVXuc+eYG/dWsb+y6WOGBnlIiAolPjKUhKgwEiNDSYhynidEOsePd8sSo8IIDzl+i2tH/j36Q1eLqy2ajyYCF4nIBUA40B1nxqlYEQlW1VogGWe2JoBCnInIC0UkGIjBmcWpvrxe422M+YacbWUszi0hPSmK4vI6Ktd+3ejk7jTllFRUUXrQSQIHqmqb3E90WDCJ0WEkRoWS1iOK8anx3zjRNz7xR4QGndC39r5xEby9opDioiJuv3AcGSlxx92mutbb0AxVnzBKK6opqXCTyEEnoeTvqaCkoorKGm+T+4kMDXISRlQoCZFO4qhPIolu2dKiGr78eBNTTu55QrGZrqXVSUFV7wHuAXBrCr9Q1WtF5J/A5cBsYAbwrrvJXPf1Z+7yT1RVRWQu8A8R+QPQB0gHPm9tfKbzKKmoYmleCe98VcTCTXu+uXB5TsPT2IgQ90QeytA+3d2T+39OivVJ4ES/WTdXRkocGSlxZGWVnvBJNzTYQ6+YcHrFhJ/Q+oeqa92k0SiJHKx2nrtJpGjvYVYX7qXsYDW13iMr3Xk8tTCP747ozWUZyYwbGE9EaIcfH9O0AV/+FdyFM9/sw8BXwPNu+fPA392O5DKcK45Q1XUi8gawHqgFblNVmxS8C6uqrSNnWzmLcktYnLuHdcX7AacZp54AE/sEcfdlp9Mj2mlKCQnq/LffRIQGExEfTL/4iOOu6/Uq+ytrKD1YzXOLtzL78+0oTtvsh2u/5r3VOwkJEkb1j2NyeiIT0xIZkRxrHeJdVJsmBVXNArLc5/k4Vw8duU4lcMVRtn8E5wom0wWpKrm7K1i0eQ+Lc0tYvrWUyhovwR4hIyWOX5w7mMnpPaiu83Ld88upqfUSEuwhs18Iw/rG+Dv8gOXxCLERzpVPl2ck89aX26lTCAn28Lfrx1DrVZbklbAkt4Tf/3szv//3ZrqHB3P6oAQmpfdgUloiAxIi2rzD2wQmqy8avyqtqGJJXgmL3drArv1VAKT2iOSqMf2ZnJ7IuNQEosK++af66k3jyc4vZXxqAge2rvJH6B1SRkocvxwTTlVsCuNTExqatyan94Dznd/Hsi2lLMktYUleCR+v2wVA39huTEpLZJJbk4iPDPXn2zA+ZEnBtKujNQnFRoQwMS2RM9ITmZTeg76x3Y65n/p2e4CsrT4Pu1NJiwsiMzOtyWUJUWFcOLIPF47sg6qyrfSQW4vYwwdrd/L6lzsAGNqnO5PSE5mUlsiYAfE+6Zsx/mFJwfiUqpK3u6IhCSzPL+NwTR3BHmFUoyahYX1jrA07wIgIAxMjGZgYyXXjU6it87KmaF9DLeKFJVv566f5hAZ7GDsgnolpiUxOT2RI7+547HfZYVlSMG2ucZPQktwSvt5fCThNQleO6XfUJiET2IKDPJzWP47T+sfx47PTOVhVy+fbypwkkVvCbz/ayG8/griIECakObWISWmJJ9QZbgKH/VeaVquqrSOnoLyhX2BtUcuahEzHEhkWzJSTejLlpJ4A7D5QydK8EpbklrIkbw/vr94JQEpCBJPcWsTpqYnERIT4M2xzHJYUTLNZk5BpSs/ocC45LZlLTktGVdmyp4LFuSUN95a8unw7HoHhybFMSktgUloPRqXEEhYcRE5BOe9tqSZ6YLndUOdnlhTMceUUlPPJxl0U7ajm/T2rWGxNQuY4RIS0ntGk9YzmhokDqanzsnLH3ob+iL98ms9TC7fQLSSIk3pFs7ZoH3Ve5b1t2bx603hLDH5k/8HmmOasKOLn/1xJ/Q2xkaE7yTypJ5PTncsTk+OsvdgcX0iQhzED4hkzIJ6fnjOYA5U1ZOeXsSR3D/9aXdxwx3VljZcP1uy0pOBHlhRMk3J3HeBPC3Ib2oXBuXv4lsxB3H5Wuv8CM51CdHgI5wxJ4pwhSVx0al+ueTabqlpnPKcXlmylpKKKH5+VTlrPwBudtLOzpGC+IW/3AR5fkMd7q4uJCAnislF9eW/1TmrrvAQJnD4o0d8hmk4mIyWOf/xwPK/N/4ILJo5k+dYy/v5ZAXNXFXPhiD7ccXYaaT2j/R1ml2FJwQCQt7uCJxbk8i83Gdxy5iB+ODmVuMhQrhmXQnZ+KWF7C6xab3wiIyWOA4NCyTw5ibNOTmLm5FSeXbyVlz/bxr9WF/Pd4b254+x0BidZcvA1SwpdXP4eJxnMXVVMeEgQPzpjEDPPSP3GMAb/GfWz0I+Rmq4kISqMu88/mZlnpPLc4nxeWraN99fs5ILhvbnjrHRO6mXJwVcsKXRRW0sO8ucFubyzsoiw4CB+ODmVmWekkhAV5u/QjGkQHxnKL887mR9OTuW5Jfm8tKyA91fv5ILhvbjj7HRO7tXd3yF2OpYUupiC0oM8sSCPd1YWERIk3DhpID86cxCJlgxMAIuLDOW/v+Mkh+eXbOVvS7fxwZqvOX+YkxxO6W3Joa1YUugitpce4s+f5PL2V0UEe4QbJgzgR2cOoke0JQPTccRGhPLzc0/ixkkDecFNDh+u/ZrvDE3ijrPTGdrHhlBvLUsKndyOskM8+Ukeb60oJMgjzDh9ADdnptIz+sRm+DImEMVGhPKzc0/ixkmpvLB0Ky8s3crH63Zx7hAnOdj8Gi3X6qQgIuHAIiDM3d+bqvqAiAzEmYozHlgBXKeq1SISBrwMZAClwJWqus3d1z3AjUAdcIeqftza+LqqHWWHeGphHm/mFOLxCN8fn8KtmYPo2d2Sgek8YiJC+Ok5g/mvSQP529KtPL9kK/9ev4uppyTxk6mWHFqiLWoKVcBZqlohIiHAEhH5EPgZ8EdVnS0if8E52T/j/ixX1TQRuQr4LXCliAzBmZpzKM4czfNFZLBNydk8RXsP8+QnebyZswNBuHZcf27JTDvhuX+N6YhiuoXwk6mDuWHiQF5cuo3nl+Qz7c+7mHpKT+48ezDDky05nKhWJwVVVaDCfRniPhQ4C7jGLX8JeBAnKUx3nwO8CTwpzjx/04HZqloFbHXncB4LfNbaGLuC4r2HeWphHm986SSDq8b059Ypg+gdYyOTmq4jplsId05N54ZJA3hp6TaeW7KVC59cwlkn9+TOs9MZ2S/W3yEGPHHO6a3ciUgQkAOkAU8BvwOyVTXNXd4P+FBVh4nIWuA8VS10l20BxuEkimxVfcUtf97d5s0mjjcTmAmQlJSUMXv27BbFXVFRQVRU4N1G35y4yiq9vJdfw6IdtShwRnIw01JDSOjW9pPXd4bPqz1ZXM3ji7gO1yrzCmr4eFsNB2tgRI8gLh4UQmrsic8U11k/rylTpuSo6ugjy9uko9lt4jlVRGKBOcApTa3m/mxqLGU9RnlTx5sFzAIYPXq0ZmZmNjdkALKysmjptr50InHt2l/J0wvzeO3zHXhVuWJMf26bMsinA9R15M/LHyyu5vFVXOcDFVW1vLRsG88tzufX2ZWcObgHd05NZ1T/49+h39U+rza9+khV94pIFjAeiBWRYFWtBZKBYne1QqAfUCgiwUAMUNaovF7jbYxr9/5Kns7awj8+347Xq1yekcxtU9JsditjjiEqLJjbpqQxY8IAXv5sG88uyufSp5dxxuAe3Hm2M8Bjdn4p41MTuvxQLm1x9VEPoMZNCN2AqTidxwuBy3GuQJoBvOtuMtd9/Zm7/BNVVRGZC/xDRP6A09GcDnze2vg6iwUbdvF01hZWF+7Fq3DZqL78+Kx0SwbGNENUWDC3ZqYx4/QB/D27gFmL8rnsmWXUzwUVGuzp8vM5tEVNoTfwktuv4AHeUNX3RGQ9MFtEHga+Ap53138e+LvbkVyGc8URqrpORN4A1gO1wG125ZHjH8sLuHfOWgCCRHji6lP57og+fo7KmI4rMiyYm88cxHXjU7j11RV8unkPANW1XrLzSy0ptIaqrgZOa6I8H+fqoSPLK4ErjrKvR4BHWhtTZ/LuyiLuf3ddoxJlW+khv8VjTGcSGRbMHWenk51fSlWtF6/C4Zqu/V3U7mgOUKrKk5/k8X/zNnNKr2jySw5SW+clJNjD+NQEf4dnTKdRP5/Dp5t2s3DTbp78JA+vV/nFuSfh6YJzjFtSCEC1XuUX/1zNWysKueS0vjx22XDWFu23jjBjfKR+ePgfn53Or95dx9NZW8jfc5A/XDnS36G1O0sKAWbvoWp+90Ulm8oL+enUwdxxdhoi0vBHa4zxnZAgD7+5ZBiDekTyyAcbuPKvh/mvdK+/w2pXlhQCyLaSg9zw4hfs2Ovl8atOZfqpff0dkjFdjohw0+RUBiREcsfsr/h1iZf04fu6zDhKbX/bq2mRL7aVccnTS9l7qJpfjg23hGCMn00dksSbN09ABK74y2d8vO5rf4fULiwpBIB3viri2meXExcRypxbJzI47sRvwTfG+M6QPt351enhDO4Vzc2v5PCXT7fQFkMDBTJLCn6kqvxp/mZ+8vpKRqXE8vatExiQGOnvsIwxjcSGeXh95nguGN6bxz7cyC/fXE11beftZ7A+BT+pqq3j7rfWMOerIi4blcyjlw4nNNhytDGBKDwkiD9fdRqDekTxxIJctpcd4i/fzyAuMtTfobU5Owv5QfnBaq577nPmfFXEL84dzO+vGGEJwZgA5/EIPztnMH+68lS+2r6Xi59eypY9FcffsIOxM1E7y99TwSVPL2Vl4V6euPo0bj8rHWc6CWNMR3DxaX15beY4KiprueSppSzNK/F3SG3KkkI7Wp5fyqXPLGN/ZS2v/XAcF4208YuM6YgyUuJ557aJ9IoJ5wcvfM4/lm/3d0htxpJCO3l7RSHff345CZGhvHPrRDJS4v0dkjGmFfrFR/DWLROYlJbIvXPW8NB766nzdvwrkywp+Jiq8od/b+Jnb6xizIB43r5lIv0TbLhrYzqD6PAQnp8xmusnDOD5JVuZ+fKXVFTV+jusVrGk4EOVNXXcOXslT3ySx/dGJ/PiDWOJiQjxd1jGmDYUHOThwYuG8tD0oWRt3sPlzyyjaO9hf4fVYpYUfKS0oorvP7ecuauK+eV5J/Hby+wKI2M6s+tOH8AL14+hqPww059cylfby/0dUovYWcoHtuyp4JKnl7GmaB9PXTOKWzPT7AojY7qAMwf34O1bJ9At1MNVs7L516qON6Nwq5OCiPQTkYUiskFE1onInW55vIjME5Fc92ecWy4i8oSI5InIahEZ1WhfM9z1c0VkRmtj84fPtpRy6dPLOFRdy2szx/PdEb39HZIxph2lJ0Xzzq0TGZEcw49f+4rH5+d2qKEx2uKO5lrg56q6QkSigRwRmQdcDyxQ1cdE5G7gbuAu4Hyc+ZfTgXHAM8A4EYkHHgBGA+ruZ66qdog6WE5BOc8tzuff678mNTGKF64fY/MnG9NFJUSF8cpN47jn7TX8cf5mviwoY3RKHJPSewT8EPhtMR3nTmCn+/yAiGwA+gLTgUx3tZeALJykMB14WZ3UmS0isSLS2113nqqWAbiJ5TzgtdbG6Gs5BeVc+dfPqPUqHoH7pw2xhGBMFxcWHMT/XTGS8JAg/rF8O4tzS3jm0y28etP4gE4M0pbVGhEZACwChgHbVTW20bJyVY0TkfeAx1R1iVu+ACdZZALhqvqwW34/cFhVf9/EcWYCMwGSkpIyZs+e3aJ4KyoqiIqKatG2jb24toqsQucyNA9waXoI0wa1fEyUtoqrrVlczWNxNU9njeu9LdW8mVsDgACXtfL80FZxTZkyJUdVRx9Z3mYD4olIFPAW8BNV3X+MjtWmFugxyr9dqDoLmAUwevRozczMbHa8AFlZWbR023o1dV5+tfxToJYggZBgD1dPHdOqbwJtEZcvWFzNY3E1T2eNK3pgOe9ty6ayxhlZ9ZIzRzExLdHvcR1Nm1x9JCIhOAnhVVV92y3e5TYL4f7c7ZYXAv0abZ4MFB+jPKDNWpTP9vJD3HP+yfzs3JMCvmpojGlfGSlxvHrTeK4Z2x8FPtm4+7jb+FOrawriVAmeBzao6h8aLZoLzAAec3++26j8dhGZjdPRvE9Vd4rIx8Bv6q9SAs4F7mltfL5UUHqQJxbkct7QXvzozEH+DscYE6Dq51gXgb8t3cpFI/swsl/s8Tf0g7aoKUwErgPOEpGV7uMCnGRwjojkAue4rwE+APKBPOBZ4FYAt4P5IeAL9/Hr+k7nQKSq3DdnLSHu3YzGGHM8d51/Mj2iw7jrrdXU1AXmRD1tcfXREpruDwA4u4n1FbjtKPt6AXihtTG1h3dWFrEkr4RfTx9Kr5hwf4djjOkAuoeH8ND0Ycz8ew6zFuVz25Q0f4f0LXZHcwuUH6zmofc2cGq/WK4dl+LvcIwxHci5Q3txwfBePL4gl60lB/0dzrdYUmiBRz/cwP7DNTx66XCCPDZ8hTGmeR68cChhwR7ufms13gAbbtuSQjNl55fyxpeF3DQ5lVN6d/d3OMaYDqhn93Duu+AUlm8t440vd/g7nG+wpNAMVbV13DtnDf3iu3Hn2en+DscY04FdOaYf41Pj+c0HG9i9v9Lf4TSwpNAMTy/cQv6egzx88XC6hQb5OxxjTAcmIjx66Qgqa708+K91/g6ngSWFE5S3u4JnsrZw0cg+nDm4h7/DMcZ0AgMTI7nz7HQ+WPM1H6/72t/hAJYUTojXq9w7Zw3hIR7unzbE3+EYYzqRmWekcnKvaH717lr2V9b4OxxLCifizZxCPt9axr0XnEKP6DB/h2OM6URCgjz89rIR7DlQxf9+tNHf4VhSOJ6Siioe+WADYwfE873R/Y6/gTHGNNPIfrHcMHEgr2Rv54tt/h3IwZLCcTz83noOVdfym0uH4bF7EowxPvLzcweTHNeNu99aTVVtnd/isKRwDItz9/DOymJuyUwjrWe0v8MxxnRiEaHBPHLJcLbsOchTC7f4LQ5LCkdxuLqO++asJTUxklszbQRUY4zvnTm4B5ec1pdnsvLY9PUBv8RgSeEonvgkl+1lh3j4kmGEh9g9CcaY9nH/tCFEh4dw99urqfPDEBiWFJqw8ev9PLson8szkpkwqPUzJBljzImKjwzlV9OG8NX2vbySXdDux7ekcASvV7n37TV07xbCfRec4u9wjDFd0PRTnZtk//ejjRTtPdyux26r6ThfEJHdIrK2UVm8iMwTkVz3Z5xbLiLyhIjkichqERnVaJsZ7vq5IjKjLWJrrlc/386K7Xv5n++eQlxk6yfXNsaY5hIRHr54GF6FH/9jBU8tzCWnoLxdjt1WNYUXgfOOKLsbWKCq6cAC9zXA+UC6+5gJPANOEgEewJmicyzwQKOpOdvF/PW7eOhf6xnetzuXnNa3PQ9tjDHf0C8+gqvG9mPF9r38/uPNXPtcdrskhjZJCqq6CDjyjovpwEvu85eAixuVv6yObCBWRHoD3wHmqWqZqpYD8/h2ovGZnIJyZv79S6rrvGzeVcGK7Xvb69DGGNOkBLe1QoGaWi/Z+aU+P2arp+M8hiRV3QmgqjtFpKdb3hdoPIB4oVt2tPJvEZGZOLUMkpKSyMrKalGAFRUVDdv+a0s19R39NbVeXpv/BQcG+af5qHFcgcTiah6Lq3ksrm8L31eHR8Cr4BEI21tAVlahT+PyZVI4mqZuC9ZjlH+7UHUWMAtg9OjRmpmZ2aJAsrKyqN92f2wRb+WuRIDQEA9XTx1DRkq7tl41GVcgsbiax+JqHovr2zKBuvgtPPrhRh64aBjfH/+f6X99FZcvrz7a5TYL4f7c7ZYXAo0HEUoGio9R3i7qs8+MCQN49abxfksIxhjT2CWjnAaTypr2GfrCl0lhLlB/BdEM4N1G5T9wr0IaD+xzm5k+Bs4VkTi3g/lct6xdrC3aR2iwh/u+e4olBGNMwOgZHU7vmHDWFO1rl+O1SfORiLyGU9NJFJFCnKuIHgPeEJEbge3AFe7qHwAXAHnAIeAGAFUtE5GHgC/c9X6tqu02XODaov2c0iuakCC7dcMYE1iG941hdWEHSgqqevVRFp3dxLoK3HaU/bwAvNAWMTWHqrKueB/TRvZp70MbY8xxjewXy7/X72Lf4RpiuoX49Fj2tRgoLD/M/spahvWJ8XcoxhjzLSOSnXPTmnaoLVhSwOlPABjWt7ufIzHGmG8b0TcWgFWFvr9/ypICsLZ4H8EeYXCSzYW1IVkAACAASURBVJlgjAk8MREhpCREWE2hvawt2k96UrQNkW2MCVgjkmNZbTUF31NV1hbtY1gfazoyxgSukckxFO+rZM+BKp8ep8snhV37qyg9WM2wvtbJbIwJXCOSnX4FX9cWunxSsE5mY0xHMLRPdzyCz+9XsKRQvA8ROKW3JQVjTOCKDAsmrWeU1RR8bW3RflITI4kI9cfYgMYYc+KczuZ9OPcA+0aXTwrri/dZf4IxpkMYmRxD6cFqn07R2aWTwv5qpXhfpd3JbIzpEIa7nc2+vF+hSyeF7fudoWiHWiezMaYDOKV3NCFBwipLCr6xbb8XgKFWUzDGdABhwUGc3Ku7Tzubu3RSKNjvpX98hM9HHTTGmLYyIjmGNYX78Pqos7nLJwW7P8EY05GMSI7hQFUtuw9ZUmhT+w7XsPuQWtORMaZDqb+zOX+f1yf7D7ikICLnicgmEckTkbt9dZz1xfsB7HJUY0yHkt4zivAQD1v3+WbO5oBKCiISBDwFnA8MAa4WkSG+ONa6Yqf3fqgNhGeM6UCCgzwMSIgkZ1cdOQXlbb7/gEoKwFggT1XzVbUamA1M98WBFuXuITwICkoP+WL3xhjjEzkF5eTurqCsUrn2uew2Twziy9ulm0tELgfOU9Wb3NfXAeNU9fYj1psJzARISkrKmD17drOOk1dexyPLK1GUUI/wyzHhpMUFzlwKFRUVREVF+TuMb7G4msfiah6L68S8t6WaN3NrAOdb/aXpIUwbFNrs/UyZMiVHVUcfWR5oA/5IE2XfylqqOguYBTB69GjNzMxs1kHWLcxD2QQIdQpVsSlkZqa1IFzfyMrKornvqT1YXM1jcTWPxXViogeW81buMhQIDfFw9dQxZKTEtdn+A635qBDo1+h1MlDc1gcZn5pAsMfJP8FBHsanJrT1IYwxxieG9e2OxwOD4zy8etP4Nk0IEHhJ4QsgXUQGikgocBUwt60PkpESx4MXDQXgZ+cMbvMP1RhjfCV3VwV1XpjaP8Qn566ASgqqWgvcDnwMbADeUNV1vjjWpaP64hGoqKr1xe6NMcYn6icGS+num9N3oPUpoKofAB/4+jgRocH0jfKwcofvJ8I2xpi2srZ4H9FhwfSIaKoLtvUCqqbQ3lJjPKzasdenE1YYY0xbWlu0nyF9uuMRSwptLjXGw/7KWrbZvQrGmA6gts7Lhp37fToSQ9dOCrHOvQmrrAnJGNMBbNlzkKpa3w7k2aWTQp9IoVtIEKt8PBG2Mca0hfpOZl/OFtmlk0KQRxjeN8ZqCsaYDmFt8T7CQzyk9vDdHdZdOikAjOwXw9ri/dTU+WYYWmOMaSvrivczpHd3gjy+6WQGSwqM7BdLda2XTV8f8HcoxhhzVF6vsr7Yt53MYEmBke6EFXa/gjEmkBWUHaKiqtan/QlgSYHkuG4kRIZav4IxJqDVdzIP9fEUwl0+KYgII/vF2hVIxpiAtrZ4H6FBHtJ7Rvv0OF0+KYDThJS7u8LGQTLGBKx1RfsZ3CuK0GDfnrYtKeBcgaQKawr3+TsUY4z5FlVlbfE+n/cngCUFAEa4nc3WhGSMCURFew+z91ANQ3185RFYUgAgPjKU/vER1tlsjAlIa4v2AzCsj287mcGSQoOR/WItKRhjAtL64n0EeYRTeltSaDcjk2Mo3lfJ7gOV/g7FGGO+YW3xftJ6RBEeEuTzY7UqKYjIFSKyTkS8IjL6iGX3iEieiGwSke80Kj/PLcsTkbsblQ8UkeUikisir7vTcbabU/s5/Qqrd1hnszEmsKwt2ufz+xPqtbamsBa4FFjUuFBEhuDMrzwUOA94WkSCRCQIeAo4HxgCXO2uC/Bb4I+qmg6UAze2MrZmGdonhiCPWGezMSagzF+/i90HqojpFtIux2tVUlDVDaq6qYlF04HZqlqlqluBPGCs+8hT1XxVrQZmA9NFRICzgDfd7V8CLm5NbM3VLTSIfnHd+NeqYnIKytvz0MYY06ScgnJufTUHgFezt7fLuclXczT3BbIbvS50ywB2HFE+DkgA9qpqbRPrf4uIzARmAiQlJZGVldWiICsqKhq2zSuvY3tZJV6Fq/66jLvGhJMW5/v2u+PFFUgsruaxuJrH4vq2f22pprrOmS64ts7La/O/4MCgUJ/GddykICLzgV5NLLpPVd892mZNlClN10z0GOs3SVVnAbMARo8erZmZmUdb9ZiysrKo33bdwjzqKz21XqiKTSEzM61F+22txnEFEoureSyu5rG4vm295kLuZgQIDfFw9dQxZKTE+TSu4yYFVZ3agv0WAv0avU4Git3nTZWXALEiEuzWFhqv3y7GpyYQFuKhssaLAgMSItrz8MYY8w37K2t4cVkBAxIiuCwjmQmDEhsSgi/56pLUucBVIhImIgOBdOBz4Asg3b3SKBSnM3quqiqwELjc3X4GcLRaiE9kpMTx6k3jueXMQUSEBPH37AKcsIwxpv397qNNlFRU8fhVp/Hjs9LbJSFA6y9JvURECoHTgfdF5GMAVV0HvAGsBz4CblPVOrcWcDvwMbABeMNdF+Au4GcikofTx/B8a2JriYyUOO46/2Tuv3AI2fllzP5ix/E3MsaYNpZTUM4rywu4fsJARrqXy7eXVnU0q+ocYM5Rlj0CPNJE+QfAB02U5+NcneR3V43px7sri/jN+xuYclJPesWE+zskY0wXUV3r5d6319C7ezg/P3dwux/f7mhugojw2KUjqK7zcv+7a60ZyRjTbp5dnM+mXQf49fRhRIb56gLRo7OkcBQDEiP52TmDmbd+Fx+s+drf4RhjuoCtJQd5fEEuFwzvxdQhSX6JwZLCMdw4aSDD+8bwwNy17D1U7e9wjDGdmKpy35w1hAV7ePDCoX6Lw5LCMQQHeXjssuGUH6rh4fc3+DscY0wn9taKIpZtKeXu80+mZ3f/9WNaUjiOoX1iuPnMVN7MKWTR5j3+DscY0wmVVlTx8PvrGZ0Sx9Vj+vs1FksKJ+DHZ6WT2iOSe+es4aDN42yMaWOPvL+Bg1W1/ObS4Xg8TQ3w0H4sKZyA8JAgfnvZCArLD/N//97s73CMMZ3I4tw9vP1VETefOYjBSdH+DseSwokaMyCe68an8LdlW1mx3UZRNca03uHqOu6bs5aBiZHcNsU/Y60dyZJCM/zyvJPo1T2cu99aTXWt19/hGGM6uCc+yWV72SEeuWRYu8yqdiIsKTRDdHgIj1wyjM27Kng6K8/f4RhjOrANO/cza1E+V7iD3QUKSwrNdNbJSUw/tQ9PLcxj864D/g7HGNMB1XmVe95eQ2y3EO694BR/h/MNlhRa4FfThhAVFswv31xNndeGwDDGNM8r2QWs3LGX+6cNIS6yXaejPy5LCi2QEBXGgxcNZeWOvby0bJu/wzHGdCA79x3mdx9vYnJ6ItNP7ePvcL7FkkILXTSyD1NO6sHvPt7EjrJD/g7HGNNBPPDuOmq9Xh65eDjO9PSBxZJCC4kID18yHI/AvXPW2Eiqxpjj+mjt1/x7/S5+MnUw/QN0dsfWTrLzOxHZKCKrRWSOiMQ2WnaPiOSJyCYR+U6j8vPcsjwRubtR+UARWS4iuSLyujszW0DrG9uNu88/mcW5Jby1osjf4RhjAtiByhoemLuWU3p358ZJA/0dzlG1tqYwDximqiOAzcA9ACIyBGeqzaHAecDTIhIkIkHAU8D5wBDganddgN8Cf1TVdKAcuLGVsbWLa8elMGZAHA+9t549B6r8HY4xJgDlFJTzgxc+Z9f+Kh69dDghQYHbSNOqyFT13+4UmwDZQLL7fDowW1WrVHUrkIczq9pYIE9V81W1GpgNTBenYe0s4E13+5eAi1sTW3vxeIRHLx3B4eo67nztK55amEdOgd3xbIxx5BSUc/WsbL7avpcgjwT8FYvSVm3hIvIv4HVVfUVEngSyVfUVd9nzwIfuquep6k1u+XXAOOBBd/00t7wf8KGqDjvKsWYCMwGSkpIyZs+e3aKYKyoqiIqKatG2R3pxbSVZhXUAhHrgl2PCSYtr2R2KbRlXW7K4msfiap7OGtdfVlWSvdM5N3iAS9NDmDao9a3jrY1rypQpOao6+sjy4871JiLzgV5NLLpPVd9117kPqAVerd+sifWVpmsmeoz1m6Sqs4BZAKNHj9bMzMyjrXpMWVlZtHTbI63x5pJV6AyWV+OFwzH9ycxM93tcbcniah6Lq3k6W1yqyl8X5ZO9cyMizskvJNjD1VPHkJES57e4jue4SUFVpx5ruYjMAKYBZ+t/qh2FQL9GqyUDxe7zpspLgFgRCXaboxqv3yFMGJTIUyF5VNV4UeCTjbuZMWEA3cND/B2aMaadVdd6uf+dtbz+5Q6mjejNteP6s2L7XsanJrRJQvClVs0KLSLnAXcBZ6pq44v15wL/EJE/AH2AdOBznBpBuogMBIpwOqOvUVUVkYXA5Tj9DDOAd1sTW3vLSInj1ZvGk51fyoHDNTy3ZCuXPb2MF64fQ7/4wLz0zBjT9vYdquGWV3NYtqWUH5+Vxk+nDsbjEU4PoPGNjqVVSQF4EggD5rk3YWSr6s2quk5E3gDW4zQr3aaqdQAicjvwMRAEvKCq69x93QXMFpGHga+A51sZW7vLSIlr+BZwxuAe3PLqCqY/tZS/XpfBmAHxfo7OGONrBaUHueHFL9hRdoj/u2Ikl2UkH3+jANOqpFDfMXyUZY8AjzRR/gHwQRPl+ThXJ3UKE9ISeee2idz44hdc82w2j146gss74B+IMebEfLmtjJl/z8Gryis3jmNcaoK/Q2qRwL1YthMYmBjJnFsnMnZgPL/45yoe+3Aj3gC/HM0Y03zvrizimmeXE9MthDm3TuywCQEsKfhcTEQIL94wlmvH9ecvn27h5ldybJ5nYzoJVeVP8zdz5+yVnNY/ljm3TmBgYqS/w2oVSwrtICTIw8MXD+PBC4cwf8MuLv/LZxTvPezvsIwxrVBZU8dPX1/Jn+bnctmoZP5+4zhiIwJ+dJ7jsqTQTkSE6ycO5IXrx1BYdojpTy1l5Y69/g7LGNMCpRVVfP+55byzspj//s5J/P6KEYQGd47Taed4Fx1I5kk9efvWCXQLCeLKv37G3FUd6nYMY7q8vN0VXPL0MlYX7ePJa07jtilpATkEdktZUvCD9KRo3rltIiOTY7njta/4w7zNNvS2MR3AsrwSLn16KYeqa5k9czzTRgTeJDmtZUnBT+IjQ/n7TWO5PCOZJxbkcvtrX1FZU+fvsIwxR/HGFzv4wQufk9Q9nDm3TmRU/8C+M7mlWnvzmmmFsOAgfnf5CNJ7RvHYRxspLDvEsz/41vhUxhg/8qry2Icb+cunW5icnshT147q1MPXWFLwMxHhR2cOYmBiJD95fSUXPbmUW5ocG9YY094OV9fx9Moqvty1hWvG9ef/XTQ0oOdCaAud+911IOcO7cWbN0/AI/DI8ko+Wvu1v0MypkvbfaCSq2Z9Rs6uOv7nu6fwyMXDOn1CAEsKAWVIn+68c/tEkqM83PxKDk9n5VkHtDF+sPHr/Vzy1DI276rgx6eFcdPk1E51hdGxWFIIMD2jw7l7bDgXjezD/360iZ//cxVVtdYBbUx7WbhpN5c/8xm1Xi//vPl0RiV1rVZ2SwoBKDRIePyqU/nZOYN5e0UR1z67nNIKm//ZGF97+bNt3PjiF/SPj+Cd2yYyrG+Mv0Nqd10rBXYgIsIdZ6eT2iOSn7+xiulPLeXn55xE8b7DHWKiDmM6ki+2lfHoBxtYsX0vU0/pyeNXnUZkWNc8PXbNd92BTBvRh35xEVz/t8/56RsrESAsxMOrN423xGBMK6kqsxbl89hHG1GFII9zNWBXTQhgzUcdwsh+sVw5xpnFVIHKGi8vL9tmw3Ab00KqyoINu7joyaU8+qGTENwFfL61zK+x+VurkoKIPCQiq0VkpYj8W0T6uOUiIk+ISJ67fFSjbWaISK77mNGoPENE1rjbPCFdpav/BJ0zpBfhIR5EnDlN311VzPmPL+a91cXUWXIw5oSoKp9s3MX0p5Zy40tfsvdwNbdlDiI8xEOQQEiwh/EdeC6EttDaOtLvVPV+ABG5A/gVcDNwPs68zOnAOOAZYJyIxAMPAKNxvvTmiMhcVS1315kJZOPMzHYe8GEr4+s0Gs8BPXZgPMV7D/PnT/K4/R9fkdYzlx+flca0EX0I8lguNeZIqkrWpj38af5mVhXuIzmuG7+9bDiXjkomJMjDWackkZ1fav11tH46zv2NXkbinOgBpgMvq3ORfbaIxIpIbyATmKeqZQAiMg84T0SygO6q+plb/jJwMZYUvqHxHNDg9Dd8uHYnf16Qx52zV/L4/FxuPyuNi0b2IbgL3GRjzPGoKp9u3sOf5ueycsde+sZ247FLnWTQeKjrI/+3ujJp7c1RIvII8ANgHzBFVfeIyHvAY6q6xF1nAXAXTlIIV9WH3fL7gcNAlrv+VLd8MnCXqk47yjFn4tQqSEpKypg9e3aLYq+oqCAqKqpF2/pSc+PyqpKzq465W2rYccBLUoQwLTWE0/sEE9yGNYfO8nm1F4uredoyLlVlbUkd7+TVsGWfl4Rw4cJBIUzq2/z/ic76eU2ZMiVHVb812NpxawoiMh/o1cSi+1T1XVW9D7hPRO4BbsdpHmrqU9cWlDdJVWcBswBGjx6tmZmZx3wPR5OVlUVLt/WllsR1FvBzrzJvwy6eWJDL82v3M684mNumpDZUkf0RV3uwuJqnM8elqizOLeFP8zezYrtTM/jNJWlcnpHc4klwOvPn1ZTjJoX6b+8n4B/A+zhJoRDo12hZMlDslmceUZ7llic3sb5pBo9H+M7QXpw7JIlPNu7m8QW53PXWGp5YkMdtU1r3j2FMIFNVluSV8Kf5ueQUlNMnJpxHLhnGFRn97G++mVrVpyAi6aqa6768CNjoPp8L3C4is3E6mvep6k4R+Rj4jYjUN96dC9yjqmUickBExgPLcZqj/tya2LoyEeHsU5I46+SeZG3ew+Pzc7l3zhqe/CSXWzIH8b0x/QgLDvJ3mMa0mqqyNK+UP83fzJcF5fSOCefhi4dxxehk+xtvodZeffSYiJwEeIECnCuPwLl66AIgDzgE3ADgnvwfAr5w1/t1faczcAvwItANp4PZOplbSUSYclJPMgf3YHFuCY8vyOX+d9fx1MIt3HxmKleN7U94iP3jmI5HVVm2xUkGX2xzksFDFw/je5YMWq21Vx9ddpRyBW47yrIXgBeaKP8SsJkEfEBEOGNwDyanJ/LZllL+tCCXB/+1nqeztvCjMwdxzdj+dAu1fyQT+FTV+Ruen8vn28ro1T2ch6YPtdpvG+q693J3QSLChLREJqQlkp1fyuPzc3novfU8k7WFH52RyrXj+xMRan8SJjAt2+L0GXy+tYyk7mH8evpQvje6n9V225idAbqo8akJjJ+ZwOdby3hiQS6PfLCBZz7dwg8np/KD01O69NgvJrB85jYTLXeTwf+7aChXjrFk4Cv2n9/FjR0Yzys3jSOnoIwnFuTx2482MmvRFm5yk0N0J56L1gS27HwnGWTnl9EzOowHLxxi/WDtwJKCASAjJZ6X/mssK3fs5YkFufzu403MWpTPjZMGMmPCAGK6WXIwvpVXXse6hXl0Dw/mgzVf81l+KT2iw3jgwiFcbcmg3VhSMN9war9YXrh+DGsK9/H4glz+MG8zzy7O578mDiR4Xy3rFubZ+DCmzX2+tZTHPq+kVjcBEBsRwq+mDeGacZYM2pslBdOk4ckxPDdjNGuL9vHkJ3k8vqD+dpRNhAQJf776NM4b1tuvMZqOS1UpKD3E4tw9LM4tIWvTHmrdMQwEuGHCAP5r0kC/xthVWVIwxzSsbwx/uS6DB+eu48Vl2wCoqVNufmUFAxIimJiWyKS0RE4flEBsRKh/gzUBbd+hGpZuKWFxbgmLc/dQWH4YgL6x3ThzcCILN+5GcYavnpTew7/BdmGWFMwJuXBkH/6RvY06heAgD9eO68/2skO8u7KYV5dvRwSG941pSBIZKXFW7e/iqmu9fLW9nCV5JSzKLWFN4V68ClFhwZw+KIGZZ6QyOb0HAxIiEBGem7OAqtgUa570M0sK5oRkpMTxyzHh3/qnranzsrpwL0tyS1maV8Kzi/J5JmsLYcEeRg+Ia0gSQ/vE2FwPnZyqsmXPQZa4TULZ+aUcrK7DI05f1e1npXNGeiIj+8U2OUBjWlwQmZlpfojcNGZJwZywpv5pQ4I8ZKTEk5ESz51T0zlYVcvnW8tYklfC0rwS/vejTfwvm4jpFsKEQQkNSSLF/XZoOrayg9UszXOag5bkllC8rxKAlIQILj6tL5PTe3D6oAS7eq0DsaRg2lRkWDBTTu7JlJN7ArD7QCWfbSllSa6TJD5c+zXgtCNPSktkYnoiEwYlkBgV5s+wzQmqqq0jp6CcJblO38Da4n2oQvfwYCYMSuS2sxKZnNaD/gkR/g7VtJAlBeNTPaPDmX5qX6af2hdVZWvJQZbmlbAkr4QP1+7k9S93AHBK7+5MSnNqEmMHxttwGwFCVcndXdHQObw8v4zDNXUEe4TT+sfy06mDmZSeyIi+MTbbXydh/3mm3YgIqT2iSO0RxXWnD6DOq6wt2tfQ1PTSsgKeXbyVkCBhVP+4hpqEnXDaV0lFFUvzSli0uYQleXvYtb8KgNTESL43OplJ6T0Ynxpvd7t3UpYUjN8EeYSR/WIZ2S+W26akcbi6ji8L/tMf8Yf5m/m/eZuJDgtm/KAEJg5KYFJ6IoN6RFl/RCvlFJQ3TFQ/tE93vtxW3nDPwPqdztTrsREhTExLZHJaIpPSE0mOsyahrsCSggkY3UKDmJzeg8nuNeplB6ud/gg3ScxbvwuApO5hTExLJDm2G1V1XuIP135jOj/TNK9XKT9UzbKiGl6cl01NnRcRJznX1CkhQUJGShz//Z2TmJxuV4x1VZYUTMCKjwzluyN6890Rzp3T20sPsXSL0x8xb93XHKiqa1j3z6s+omf3cBIiQ4mPDCUhKqzheXxkKAmRYSREhZIQGUpcZGibzFkdCKpq6yg7WE3JgWpKKqrcRzWljZ7X/yw7WIX3iJnPVeG05Bhun5LO2IHxNjquaZukICK/AH4H9FDVEnHq9o/jzL52CLheVVe4684A/sfd9GFVfcktz+A/M699ANzpTtZjDAD9EyLon9Cfq8f258lPnHGZ6k9yaT2j6BPbjdKKavL3HOTLbeWUH6r+1kmwXvfwYBKjwv6TNKK+mTwaJ5P4yNB2m+dXVTlYXUfJgW+e1Evrfx6s+kYC2F9Z2+R+uoUEkRgdSmJUGMlxEZzWP5aEyDASo0JZtSGX97fWUef1EhLs4d4LhtjNYqZBq5OCiPQDzgG2Nyo+H0h3H+OAZ4BxIhIPPACMBhTIEZG5qlrurjMTyMZJCudhU3Kaozh9UCKhC/OoqfUSJHD/tKHfOrHVeZW9h6opO1hN6cH//CytqPpPWUU1BaWHWLG9nLKDR08i0eHBjWoezsm1cUKpTx71CSUsOKhh1M9xA+MZmBhJ6cFqSg5UsafxSb7im9/wSyqqqKr1NhlDbEQICZHOif6UPt1JdJ8nRDnxJEaHkRgZRmJ06DGv3sqqKeD754xs6FOwhGAaa4uawh+BXwLvNiqbDrzsftPPFpFYEekNZALz6udlFpF5wHkikgV0V9XP3PKXgYuxpGCOIiMljldvGk92filhewuaPLEFecRpRooKI/0E9un1KvsO1/wngVRUNTxvnFAKyw+xqnAv5QerqT1KFukWEsThmjpg01GPF+SRhpN8QlQog3pEkRgd9o2yxKgwekSHERfRtrWVjJQ4SwamSa1KCiJyEVCkqquOuBqkL7Cj0etCt+xY5YVNlB/tuDNxahUkJSWRlZXVovgrKipavK0vWVwnbqhARchhn8TVDUgGkoOA7u6jQTCqQRyqhQPVyv5q5YD72F+trN5TR95exRnzE0b2CGJCn2C6hwrdw4SYUCEiBDwiQB1w2H3g1KEPOI8SnEdbCsTfI1hczeWruI6bFERkPtCriUX3AfcC5za1WRNl2oLyJqnqLGAWwOjRozUzM/Noqx5TVlYWLd3Wlyyu5gnEuHIKyrn6r8uoU2fUz19dPi5gvpkH4ucFFldz+Squ4yYFVZ3aVLmIDAcGAvW1hGRghYiMxfmm36/R6slAsVueeUR5llue3MT6xnRIRxtA0JhA1+JGSlVdo6o9VXWAqg7AObGPUtWvgbnAD8QxHtinqjuBj4FzRSROROJwahkfu8sOiMh498qlH/DNPgpjOpy0uCBum5JmCcF0KL66KPkDnMtR83AuSb0BQFXLROQh4At3vV/XdzoDt/CfS1I/xDqZjTGm3bVZUnBrC/XPFbjtKOu9ALzQRPmXwLC2iscYY0zzdY7bOo0xxrQJSwrGGGMaWFIwxhjTwJKCMcaYBtLRx5wTkT1AQQs3T6TtbxhtCxZX81hczWNxNU9njStFVXscWdjhk0JriMiXqjra33EcyeJqHoureSyu5ulqcVnzkTHGmAaWFIwxxjTo6klhlr8DOAqLq3ksruaxuJqnS8XVpfsUjDHGfFNXrykYY4xpxJKCMcaYBp0+KYjIL0RERSTRfS0i8oSI5InIahEZ1WjdGSKS6z5mNCrPEJE17jZPyBHTzDUznofc464UkX+LSJ8Aiet3IrLRPfYcEYlttOwe9xibROQ7jcrPc8vyROTuRuUDRWS5G+/rIhLairiuEJH/3975xEZVRXH4O6EVjFhpFbWhJNjEGCEaNQ3RaEwjRrEQqjsSV+JK4r+4MOgkGBcuwIUsNOnCmEgCIhSNEWOwGhtdCDUCxWoDtsWEakNjENENihwX93Tmzfhm2nkzb2Zozi+5mfvOfXfuN++e9t53z5u5P4jIJRHpKiirG9ccuGMZUmzvHRGZFpGRiK1NRAbs8w7Yz9Un8rWETMtF5EsRGbU+fK5BuBaJyJCIDBvXo+0iAwAABIVJREFUq2aP9Q8RWWjHY1a+IvJesT5YId8CETkqIgfqwqWq8zYRNvo5SPhy23Vm6yH8LLcAdwOHzd4GTNhrq+VbrWwIuMfqfAo8UgFTSyT/LNDXIFwPAU2W3wZss/xKYBhYSNhUaRxYYGkc6ASusHNWWp29wEbL9wFPVcB1K3ALYTOmroi9rlyzMBdlSNHX7wfuAkYitu3AFstvifRp2b6WkKmdsMcKwNXASeu3enMJsNjyzcBhay/WP4DN5P5ONwLvl/LBKvTlC8Bu4EApv02La77fKbwBvEj+1p69wE4NOgQsEZF24GFgQFXPqurvwACw1spaVPUbDVd8J/BoUiBVPR85vCrCVm+uz1T1oh0eIrcTXi+wR1UvqOopwh4Zqy2NqeqEqv4N7AF67W7lAaDf6r9bIdeoqp6IKaor1yyKZUipLQBU9SvgbIG5l/A5If/zluVrFTBNqeoRy/8JjBL2Xq83l6rqX3bYbEkp7h9R3n5gjflTMR9MLBHpANYBb9txKb9NhWveDgoisgH4RVWHC4qWAacjx5NmK2WfjLFXwvaaiJwGHge2NgpXRJvIbXJULte1wLnIAFNNrqgalasUW611g4ZdDbHX681e7rWrWLa0cSdhVl53LluiOQZMEwaZcYr7R7Z9K/+D4E9pXK8dhInsJTsu5bepcKW181pNJCKfAzfGFGWAlwlLIv+rFmPTBPZEXKr6kapmgIyIvAQ8DbzSCFx2Tga4COyaqVaknbgJRWpccdXS5qpAtWwriarmU3NqTGQxsB94XlXPS/HQV824VPVf4A4JsbMPCcuUxdqoCZeIrAemVfU7Eemepe3UuC7rQUFVH4yzi8hthLW0YXPADuCIiKwmjJrLI6d3AL+avbvAPmj2jpjzy+aK0W7gE8KgUHcuC+CtB9bYkhQluChi/41w299ks5dqXq+oUueqQKXYaqkzItKuqlO2DDNt9nJ9LbFEpJkwIOxS1Q8ahWtGqnpORAYJMYVi/jHDNSkiTcA1hKW6avfzvcAGEekBFgEthDuH2nJVGhS5HBLwM7lA8zryg1lDmgtmnSIEslot32Zl39q5MwHdngpYbo7knwH6G4RrLfAjsLTAvor8oNUEIZDaZPmbyAVTV1mdfeQHxjZXoQ8HyQ80NwRXEdaiDCn7+QryA82vkx/Q3Z7U1xLyCCHWtaPAXm+upcASy18JfE2YDMX6B2Fr4WhAd28pH6xSX3aTCzTXlCtVJ22URP6gIMBbhDXE78n/R7OJEJQZA56I2LuAEavzJvZN8IQs++29jgMfA8sahGuMsA55zFJfpCxjbZwg8oQT4WmRk1aWidg7CU9GjZlDL6yA6zHCzOcCcAY42Ahcc+COZUixvfeAKeAfu15PEtaXvwB+steZyUTZvpaQ6T7CssXxiF/1NADX7cBR4xoBtpbyD8KsfZ/Zh4DO2XywCv3ZTW5QqCmX/8yFy+VyubKat08fuVwul6t8+aDgcrlcrqx8UHC5XC5XVj4ouFwulysrHxRcLpfLlZUPCi6Xy+XKygcFl8vlcmX1HyOBz9P/J+BfAAAAAElFTkSuQmCC\n", "text/plain": [ "