{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial: advection-diffusion kernels in Parcels\n", "In Eulerian ocean models, sub-grid scale dispersion of tracers such as heat, salt, or nutrients is often parameterized as a diffusive process. In Lagrangian particle simulations, sub-grid scale effects can be parameterized as a stochastic process, randomly displacing a particle position in proportion to the local eddy diffusivity ([Van Sebille et al. 2018](https://doi.org/10.1016/j.ocemod.2017.11.008)). Parameterizing sub-grid scale dispersion may be especially important when coarse velocity fields are used that do not resolve mesoscale eddies ([Shah et al., 2017](https://doi.org/10.1175/JPO-D-16-0098.1)). This tutorial explains how to use a sub-grid scale parameterization in _Parcels_ that is consistent with the advection-diffusion equation used in Eulerian models.\n", "\n", "## Stochastic differential equations (SDE) consistent with advection-diffusion\n", "The time-evolution of a stochastic process is described by a stochastic differential equation. The time-evolution of the conditional probability density of a stochastic process is described by a Fokker-Planck equation (FPE). The advection-diffusion equation, describing the evolution of a tracer, can be written as a Fokker-Planck equation. Therefore, we can formulate a stochastic differential equation for a particle in the Lagrangian frame undergoing advection with stochastic noise proportional to the local diffusivity in a way that is consistent with advection-diffusion in the Eulerian frame. For details, see [Shah et al., 2011](https://doi.org/10.1016/j.ocemod.2011.05.008) and [van Sebille et al., 2018](https://doi.org/10.1016/j.ocemod.2017.11.008).\n", "\n", "The stochastic differential equation for a particle trajectory including diffusion is\n", "\n", "\\begin{aligned}\n", " d\\mathbf{X}(t) &\\overset{\\text{Îto}}{=} (\\mathbf{u} + \\nabla \\cdot \\mathbf{K}) dt + \\mathbf{V}(t, \\mathbf{X})\\cdot d\\mathbf{W}(t), \\\\\n", " \\mathbf{X}(t_0) &= \\mathbf{x}_0,\n", "\\end{aligned}\n", "\n", "where $\\mathbf{X}$ is the particle position vector ($\\mathbf{x}_0$ being the initial position vector), $\\mathbf{u}$ the velocity vector, $\\mathbf{K} = \\frac{1}{2} \\mathbf{V} \\cdot \\mathbf{V}^T$ the diffusivity tensor, and $d\\mathbf{W}(t)$ a Wiener increment (normally distributed with zero mean and variance $dt$). Particle distributions obtained by solving the above equation are therefore consistent with Eulerian concentrations found by solving the advection-diffusion equation. \n", "\n", "In three-dimensional ocean models diffusion operates along slopes of neutral buoyancy. To account for these slopes, the 3D diffusivity tensor $\\mathbf{K}$ (and therefore $\\mathbf{V}$) contains off-diagonal components. Three-dimensional advection-diffusion is not yet implemented in _Parcels_, but it is currently under development. Here we instead focus on the simpler case of diffusion in a horizontal plane, where diffusivity is specified only in the zonal and meridional direction, i.e. \n", "$$\\mathbf{K}(x,y)=\\begin{bmatrix}\n", "K_x(x,y) & 0\\\\\n", "0 & K_y(x,y)\n", "\\end{bmatrix}.$$ \n", "\n", "The above stochastic differential equation then becomes \n", "\n", "\n", "\\begin{align}\n", " dX(t) &= a_x dt + b_x dW_x(t), \\quad &X(t_0) = x_0,\\\\\n", " dY(t) &= a_y dt + b_y dW_y(t), \\quad &Y(t_0) = y_0,\n", "\\end{align}\n", "\n", "where $a_i = v_i + \\partial_i K_i(x, y)$ is the deterministic drift term and $b_i = \\sqrt{2K_i(x, y)}$ a stochastic noise term ($\\partial_i$ denotes the partial derivative with respect to $i$).\n", "\n", "## Numerical Approximations of SDEs\n", "The simplest numerical approximation of the above SDEs is obtained by replacing $dt$ by a finite time discrete step $\\Delta t$ and $dW$ by a discrete increment $\\Delta W$, yielding the **Euler-Maruyama (EM) scheme** ([Maruyama, 1955](https://link.springer.com/article/10.1007/BF02846028)):\n", "$$\n", "\$$\n", " X_{n+1} = X_n + a_x \\Delta t + b_x \\Delta W_{n, x},\n", "\$$\n", "$$\n", "with a similar expression for $Y$.\n", "\n", "A higher-order scheme is found by including extra terms from a Taylor expansion on our SDE, yielding the **Milstein scheme of order 1 (M1)**:\n", "$$\n", "\$$\n", " X_{n+1} = X_n + a_x \\Delta t + b_x \\Delta W_x + \\frac{1}{2}b_x \\partial_x b_x(\\Delta W_{n, x}^2 - \\Delta t),\n", "\$$\n", "$$\n", "which can be rewritten by explicitly writing $b_x\\partial_x b_x$ as $\\partial_x K_x(z)$: \n", "$$\n", "\$$\n", " X_{n+1} = X_n + v_x \\Delta t + \\frac{1}{2}\\partial_x K_x(\\Delta W_{n, x}^2 + \\Delta t) + b\\Delta W_n.\n", "\$$\n", "$$\n", "The extra term in the M1 scheme provides extra accuracy at negligible computational cost.\n", "\n", "The spatial derivatives in the EM and M1 schemes can be approximated by a central difference. Higher order numerical schemes (see [Gräwe et al., 2012](https://doi.org/10.1007/s10236-012-0523-y)) include higher order derivatives. Since Parcels uses bilinear interpolation, these higher order derivatives cannot be computed, meaning that higher order numerical schemes cannot be used.\n", "\n", "An overview of numerical approximations for SDEs in a particle tracking setting can be found in [Gräwe (2011)](https://doi.org/10.1016/j.ocemod.2010.10.002).\n", "\n", "## Using Advection-Diffusion Kernels in Parcels\n", "The EM and M1 advection-diffusion approximations are available as AdvectionDiffusionEM and AdvectionDiffusionM1, respectively. The AdvectionDiffusionM1 kernel should be the default choice, as the increased accuracy comes at negligible computational cost. \n", "\n", "The advection component of these kernels is similar to that of the Explicit Euler advection kernel (AdvectionEE). In the special case where diffusivity is constant over the entire domain, the diffusion-only kernel DiffusionUniformKh can be used in combination with an advection kernel of choice. Since the diffusivity here is space-independent, gradients are not calculated, increasing efficiency. The diffusion-step can in this case be computed after or before advection, thus allowing you to chain kernels using the + operator.\n", "\n", "Just like velocities, diffusivities are passed to Parcels in the form of Field objects. When using DiffusionUniformKh, they should be added to the FieldSet object as constant fields, e.g. fieldset.add_constant_field(\"Kh_zonal\", 1, mesh=\"flat\").\n", "\n", "To make a central difference approximation for computing the gradient in diffusivity, a resolution for this approximation dres is needed: _Parcels_ approximates the gradients in diffusivities by using their values at the particle's location ± dres (in both $x$ and $y$). A value of dres must be specified and added to the FieldSet by the user (e.g. fieldset.add_constant(\"dres\", 0.01)). Currently, it is unclear what the best value of dres is. From experience, its size of dres should be smaller than the spatial resolution of the data, but within reasonable limits of machine precision to avoid numerical errors. We are working on a method to compute gradients differently so that specifying dres is not necessary anymore.\n", "\n", "## Example: Impermeable Diffusivity Profile\n", "\n", "Let's see the AdvectionDiffusionM1 in action and see why it's preferable over the AdvectionDiffusionEM kernel. To do so, we create an idealized profile with diffusivities $K_\\text{zonal}$ uniform everywhere ($K_\\text{zonal} = \\bar{K}=0.5$) and $K_\\text{meridional}$ constant in the zonal direction, while having the following profile in the meridional direction:\n", "\n", "$$K_\\text{meridional}(y) = \\bar{K}\\frac{2(1+\\alpha)(1+2\\alpha)}{\\alpha^2H^{1+1/\\alpha}} \\begin{cases}\n", "y(L-2y)^{1/\\alpha},\\quad 0 \\leq y \\leq L/2,\\\\\n", "(L-y)(2y-1)^{1/a},\\quad H/2 \\leq y \\leq L,\n", "\\end{cases}$$\n", "with $L$ being the basin length scale, $\\alpha$ as a parameter determining the steepness in the gradient in the profile. This profile is similar to that used by [Gräwe (2011)](https://doi.org/10.1016/j.ocemod.2010.10.002), now used in the meridional direction for illustrative purposes.\n", "\n", "Let's plot $K_\\text{meridional}(y)$:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "from datetime import timedelta\n", "from parcels import ParcelsRandom\n", "from parcels import (FieldSet, ParticleSet, JITParticle,\n", " DiffusionUniformKh, AdvectionDiffusionM1, AdvectionDiffusionEM)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "K_bar = 0.5 # Average diffusivity\n", "alpha = 1. # Profile steepness\n", "L = 1. # Basin scale\n", "Ny = 103 # Number of grid cells in y_direction (101 +2, one level above and one below, where fields are set to zero)\n", "dy = 1.03/Ny # Spatial resolution\n", "y = np.linspace(-0.01, 1.01, 103) # y-coordinates for grid\n", "y_K = np.linspace(0., 1., 101) # y-coordinates used for setting diffusivity\n", "beta = np.zeros(y_K.shape) # Placeholder for fraction term in K(y) formula\n", "\n", "for yi in range(len(y_K)):\n", " if y_K[yi] < L/2:\n", " beta[yi] = y_K[yi]*np.power(L - 2*y_K[yi], 1/alpha)\n", " elif y_K[yi] >= L/2:\n", " beta[yi] = (L - y_K[yi])*np.power(2*y_K[yi] - L, 1/alpha)\n", "Kh_meridional = 0.1*(2*(1+alpha)*(1+2*alpha))/(alpha**2*np.power(L, 1+1/alpha))*beta\n", "Kh_meridional = np.concatenate((np.array([0]), Kh_meridional, np.array([0])))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1sElEQVR4nO3deZwkZ33n+c+TWZX3fdR9dPXdakndkkq3LEBIRgKMDJ4ZMDaLsY3Aa3lmmbUHr+2X7V3WO3jW6wHGjBmNkG2wGRljZGQwaCUwBgkkunW0Wn0f1d11V2VlVd5HVeYzf0RUdtbR3dXqyqqsyt/79cpXREY8mfl0duXzjXgi4gmltUYIIUTjsqx3BYQQQqwvCQIhhGhwEgRCCNHgJAiEEKLBSRAIIUSDkyAQQogG17TeFbhakUhEb9myZb2rIYQQG8rLL78c01pHl1u34YJgy5YtHDx4cL2rIYQQG4pS6vyl1knXkBBCNDgJAiGEaHASBEII0eAkCIQQosFJEAghRIOTIBBCiAZXsyBQSj2hlJpQSr1xifVKKfU5pdRppdTrSqmba1UXgPFknmePjpPKz9byY4QQYsOp5R7BXwIPXmb9Q8AO8/EI8Oc1rAsHz03z0S8dZDSRr+XHCCHEhlOzC8q01j9QSm25TJGHgS9p4844LyqlAkqpdq31aK3qBEYg5GdLhNw2Qm4bLtuGu6ZOCCFW1Xq2gp3AYNXzIXPZkiBQSj2CsddAT0/Pm/qwkNsGwO88dXjBckezhbDbTshtI+i2ETYDovpRvcznaMZiUW+qDkIIUY/WMwiWa02XvW+m1vox4DGA/v7+N3VvzTu3hXnht+9jLJFnOlMknikylSkSzxSYyhQry85OpolnimSLpWXfx2pRBF1GOATdzZUQCblthD3zwWEn4rER9tgJOCU4hBD1bT2DYAjornreBYzU8gM7A046A84Vlc3PlogvCox4ZtacFplKG+uOjSWJZ4rMZJc/CG1RVMJhPigiHruxl+G5GBpGkNjxOZpQSoJDCLF21jMIngYeVUo9CdwOJGp9fOBqOJqtdAScdKwwOOZKZaazs2ZIFIhlisTTxt5GLG3ueaSLHBlJMpUukMzPLfs+zVZV2csIe4zQqJ6PVKZGsNibrKv5zxZCNKCaBYFS6n8AbwUiSqkh4A+AZgCt9ReAfwLeCZwGssBHalWXtdBktRD12ol67YD3iuULcyWmM7PE0uYehhkUU2aQzM+fm8owlb50V5XX0UTUDIXqgFguNDx22dsQQixVy7OGfv4K6zXw67X6/Hpnb7LS5rfS5nesqHyuWCKWLpgPc6/DnJ9ffmoizYtnp5i+RDeVvcmyJCAiXmM+6rUT9diJmGHmldAQomHIuZMbhNNmpTvkojvkumLZ2VKZ6UyRyUuGRpHRRJ43RhLE0kVK5aXH3+dDY34vJ+q1LwiMqNdG1OMg4pVTcIXY6OQXvAk1Wy20+By0+K68t1Eua2Zys0ymjLCYTBUWzqcLDMazvHphmqlMEb3MOVtum3VJWLR47bR4HUR9xnzUayfstmOVM6iEqDsSBA3OYlGV0193XeHYxlypTDxbrISFERjFSmDEUgVOjqd44XRs2YPhVosi7LbR4jNCYj4gjKnDXG4sk4PgQqwdCQKxYk1Wi9mAX3lPIz9bYjJVYCJVYDKVZyJVYCJZYCKVZzJVYDyZ5/Bwgql0gWV6pvA7m429Cp+dVq+xd9Pqs9NqTlvM4JDAEOLaSRCImnA0r+yYRqmsmUrPB4YRFEZgmPOpAi8NxJlI5ZktLU2MkNtGi/diQLSaXWKt5rI2v4Ow20aTVQbaFeJSJAjEurJa1IqOZ8wfyxhP5qsehcp0IpXn+FiSydTSPQyLgojnYli0+R20+520+hy0+x2VqdsuPwfRmOQvX2wI1ccy9rT7Lllufg9jPFlgzAyMCTMsxpJ5hqZzHDw/veyV4F5HE23mXkTbfEj4jWmbz0mb30HQ1Syn1YpNR4JAbCrVexg34L9kufxsibFEnrFkfuE0kWc0mefk+OSyexe2JkslLDr8DtoDTjr8DjoCTtr9TjoCDvxOCQuxsUgQiIbkaLayJeJmS8R9yTJzpTKT6UIlIObDYtR8fvD8NGOvjzK3KC2czVbaAw46/E7aq8KieuqRbihRR+SvUYhLaLJaaPcbW/qXUi5rYukCI4k8IzM5RmZyjCbyjCZyjMzk+cGpSSZShSXXX/gcTXSYgyB2BpdOI267jFor1owEgRDXwFLVFbW/O7BsmdlSmfGksSdhhMV8UOQYms7xk3NxUouuu7A1WejwOy6GQ8BVNW8cr7A1yZlQYnVIEAhRY81WC11BF13BS59Km8zPMjKTY3g6x7A5HTKn3z9h7FVUUwrafA66gy66Qk66gi66g87KKbttPodcxS1WTIJAiDrgczTja2tmd9vyZ0QV5kqMzuQXhMTQdJaheI4Xz0wxmhxe0P3UbFV0BJx0B110zwdFyAiLrqCLiMcmB7RFhQSBEBuAvenyB7eLc2VGZnIMTmcZjM9PswxN53j26DixdHFBeZfNSk/IRU/IRW/YRU/YTa853xlwygV4DUaCQIhNwNZkuWxQZItzDE3nGIwbAXEhnuNCPMNALMO/nJykMFeulG2yKDqDzkpI9Ibc9IQvzjttMqzHZiNBIEQDcNma2NnqZWfr0oEFy2XNRKrA+akM5+NZYzqV5UI8yz8eGiWRW3jxXZvPQV/ETV/UzdaImy1hY7476JID2BuUBIEQDc5iUcbV1H4Ht28NL1mfyM5yPm6Ew7lYhoGpDOdiGb59eHTBTZCsFkVX0GmExKJHh98pp8PWMQkCIcRl+V3N3OgKcGNXYMm66UyRgakMA5MZzk1lOBsz5n8yEF9we1Vns5WtUTfbWzxsi3rY3mI8toTdshdRByQIhBBvWtBtI+i2cXNPcMFyrY3uprOTGc7G0pyZyHB6Ms3Bc9N847WRSjmrRdEbcrG1KhzmH3L19dqRb1oIseqUUuZorw7u3LawuylbnOPsZIbTE2lOT6Q5M2lM/+XkxIKhxruCTna3Gcc1dpnTbVGP7EHUgASBEGJNuWxNXN/p5/rOhYMCzpbKXIhnOT2R5tR4iuNjKU6Op/j+icnKeE5NFkVfxM3ONi+7zIDY0+ajO+SU6yKugQSBEKIuNFstbIsaxxDesbetsrw4V2YgluH4WJKT4ylOjKU5PJTgW6+PVsp47U3s6fCxt8PHde0+9nb42d4iew8rJUEghKhrtiYLu9qMrf9qmcIcJ8dTHBtNcXQ0wZGRJE/+ZJDcrHGQ2ma1sKPVYwaDj+s7/ezt8Mt1EMuQIBBCbEhuexM39QS5qepAdamsGYhlODqa5MhIgqMjSb57fIK/e3kIMA5O72z1sr/bz41dAfZ1BdjZ6mn4K6mVXjw+bp3r7+/XBw8eXO9qCCE2CK0148kCh4cTvD40w2uDMxwanCFpjvjqaLZwfYeffd0B9nUH6O8N0hG49NDjG5VS6mWtdf+y6yQIhBCNRmvN+aksh8xgeH0owRvDicpQG50BJ/1bgvRvCXHrliA7W7wb/oK4ywWBdA0JIRqOUqoyNtPD+zsB46yl46MpDp6Pc/DcND86M1W55sHnaOKWXiMY7twW5sZO/6bqTpI9AiGEWIbWmsF4jgPn4hw8H+fAuWlOT6QB4yylO7aFuWd7hLu3R9gWddf96auyRyCEEFdJKUVP2EVP2MXP3dIFwFS6wI/PTvHC6RjPn47x7NFxwBiI7+7tEd6yK8pbdkbxO5vXs+pXTfYIhBDiTbowleWFM0Yo/Oh0jOnsLE0WxR1bw9y/p4X7r2u97J3p1tK6HSxWSj0IfBawAo9rrT+9aL0f+GugB2Pv5E+01n9xufeUIBBC1KNSWfPa4AzPHh3nuWPjlW6kPe0+fvq6Vn72pk76LnG/iLWwLkGglLICJ4EHgCHgAPDzWuujVWV+B/BrrT+plIoCJ4A2rXVxufcECQIhxMYwEMvw3NFxnj06zoHzcbSG/d0B3ntTJz+zr4OQ27am9VmvYwS3Aae11mfNSjwJPAwcrSqjAa8yjrJ4gDgwV8M6CSHEmuiLuPnovVv56L1bGUvkefrQME+9OsIfPH2ET33zKG/dFeUX7ujlLTui635qai2DoBMYrHo+BNy+qMyfAU8DI4AXeL/WuryoDEqpR4BHAHp6empSWSGEqJU2v4NH7t3GI/du4/hYkqdeHeapV4Z57i8OsC3q5pfv6eN9N3Wt2/AXtTwRdrmIW9wP9Q7gNaAD2A/8mVLKt+RFWj+mte7XWvdHo9HVrqcQQqyZ3W0+/o+H9vD8J+/jM+/fj8vWxO8+9QZ3ffq7fO67p8hV3dBnrdQyCIaA7qrnXRhb/tU+AnxdG04DA8DuGtZJCCHqgq3Jws/e1MnTj97NVz92J7f0hvjTZ0/ytj/5Pk+9OkS5vHZndNYyCA4AO5RSfUopG/ABjG6gaheAtwMopVqBXcDZGtZJCCHqilKK2/pCPP7hfr76sTuJeu184m8P8b4//xHnpzJrUoeaBYHWeg54FHgGOAZ8VWt9RCn1caXUx81inwLuUkodBr4LfFJrHatVnYQQop7d1hfiG79+N3/yr/cxEMvwM//leb5/YqLmnysXlAkhRB26MJXlY3/9MsfHkvw/772Bn7/t2k6Uudzpo5tn1CQhhNhEesIuvv5rd3Hvjii//403eG1wpmafJUEghBB1ymmz8rkP3ESL18GjX3mF2dKSs+tXhQSBEELUMb+rmd9+aDdD0zkODydq8hkSBEIIUefu3BYG4MBAvCbvL0EghBB1biZrDL8WcNVmeGsJAiGEqHPfen0MgDu2hmvy/hIEQghRx165MM1/+d4pHtzbRm+4NsNYSxAIIUSdeuXCNB/78su0Bxz88b+6sWafI0EghBB16CsvXeD9/+3HOJutPPHhW2t6+0u5Z7EQQtSR0xMp/u9vHeP7Jya5d2eUz31gPwFXbW9iI0EghBB1IJ4p8tnnTvLXL13AZbPye+/aw0fu7sO6BjetkSAQQoh1dHYyzV+8cI6vvTxEYa7EB2/v4RP37yTssa9ZHSQIhBBijZXLmhfPTvHF5wf47vEJbFYL79nfwcfu3cqOVu+a10eCQAgh1sjpiTRPvTrEP7w6wvBMjrDbxr99+w5+8Y4eWryOdauXBIEQQtTQWCLPPx0e5alXhzk8nMCi4J4dUX7zHTt56Pp2HM3rc5/iahIEQgixirTWHB9L8ezRcZ47Ns7rQ8ZAcXs7fPzeu/bwnn0dtPjWb+t/ORIEQghxjfKzJX4yEOd7xyd47tg4Q9M5lIL93QF+6x27+OnrWtel73+lJAiEEOIqlcqaw8MJXjgd4/lTMV6+ME1xroy9ycJP7Yjw6Nu2c9+elnXt978aEgRCCHEFpbLm5HiKA+fivHA6xo/PTJHMzwGwp93Hh+/s5e7tEW7vC+O0rX+f/9WSIBBCiEXysyUODc5w8Pw0B87Fefn8NCmz4e8MOHno+nbu3hHhrm1hImt4vn+tSBAIIRqa1prRRJ7Xh2Z49cIMB87FOTycYLakAdjR4uHdN3Zw65Yg/b0hukNOlKr91b5rSYJACNFQZrJFXh9KcGhwhkNDCQ4NzTCZKgBgs1q4ocvPL9/Tx629IW7pDRJ013acn3ogQSCE2LQS2VmOjCY4OpLk8LDR+J+bylbWb4u6+akdEfZ1BdjXHWB3m7cuzutfaxIEQogNT2vNSCLPkeEER0eTHBlJcnQkyfBMrlKmzedgX7eff3NrN/u7Alzf5cfnqN3QzhuJBIEQYkPJFuc4NZ7mxHiKE2Mpjo0mOTqaZCY7C4BSsDXi5ubeIB+6s5fr2n1c1+HbFAd1a0WCQAhRl2ZLZQZiGY6PpTg5luLEeIqT4ykuxLNo4zgu9iYLu9u8PHR9G9d1+Nnb4WN3mxeXTZq2qyHflhBiXRXmSpyLZTk9kTYek2lOjqU4G0tXztyxWhR9ETfXd/j5uZu72NnqZVebl56Qa03G69/sJAiEEGsikZvlzKTR2J+ZSFfmL8SzlPXFcp0BJ7vavNy3p4VdrV52tnrZ1uLG3tR4B3HXigSBEGLVFOfKDE5nORfLMBDLcDaWYWAyw+nJdOUUTTBO0+yLuLmuw8d79nWwrcXDtqjx2IhX5m50EgRCiKtSLmtGEjnOxbIMxNKcjWUqDf/gdI5S1eZ9wNXMlrCbt+yMst1s7Le3eOgOOmmyWtbxXyGq1TQIlFIPAp8FrMDjWutPL1PmrcBngGYgprV+Sy3rJIS4srlSmeGZHOenspyPZ7kwlTHmp7Kcm8pQmCtXyjqbrfRF3Ozt9PPuGzvoi7jpi7rpC7sb4mKszaBmQaCUsgKfBx4AhoADSqmntdZHq8oEgP8KPKi1vqCUaqlVfYQQC2WLc1yIG437haks5+MXG/vhmYVb9vYmCz0hF71hF/fujNAX8dAXcbM16qbFa990Qy40mlruEdwGnNZanwVQSj0JPAwcrSrzQeDrWusLAFrriRrWR4iGMlcqM5rIMzidZSieY3A6y2A8y+B0jgvx7II+ewC/s5nesIsbu/z8zL52esNuekMuesNGY2+Rs3M2rVoGQScwWPV8CLh9UZmdQLNS6vuAF/is1vpLi99IKfUI8AhAT09PTSorxEZTLmsm0wWGprMMxnNmI2/OT2cZTeQXbNVbFLT7nXSHnLx1Z5TesIueSmPvIuCSbpxGVcsgWG7zQS963gTcArwdcAI/Vkq9qLU+ueBFWj8GPAbQ39+/+D2E2JTmt+hHZnIMz+QYnjanVfPVffUAUa+d7qCTW3qDdAdddAWddIdcdAddtAccNMsBWrGMWgbBENBd9bwLGFmmTExrnQEySqkfAPuAkwixyeWKJYZnsgzP5M2GPXuxsZ/OMZbMLzi/HiDisdMZcLC73cv917XSHXTSFXTRHTKmjThgmrh2tQyCA8AOpVQfMAx8AOOYQLVvAH+mlGoCbBhdR/+5hnUSYk3MlsqMJfKMJvKMJnKMzFycjszkGE3kmDbHxpnXZFG0+R10BpzcsTVMZ9BJZ8BZmXYEnNLQi5qoWRBoreeUUo8Cz2CcPvqE1vqIUurj5vovaK2PKaW+A7wOlDFOMX2jVnUSYjXM980bDXq+Mh1N5BieyTM6k2MyXaiMhzPP52iiw2zQb+oJ0BFY2NC3+hwyXIJYF0ov/mutc/39/frgwYPrXQ2xSRXmSkwkC4wm8owl84wlcowlCowlc4wl8owl8kykCswt6rNxNltpDzjo8Dtp9zvMBt9Bu//i1G2X6zfF+lFKvay17l9unfxlioagtSZdmGM8aXTXzDfqRmN/cTqVKS55rctmpc3voM3n4I5tYdp8Dtr984280dD7nc1yLr3YsCQIxIaXK5aYSOUZTxYYT+arHsbziZQxzRZLS14bcttoNRv2G7sCtJsNfpv/4sNrb5JGXmxqEgSibhXmSkymCkykCkwkqxv6hQ1+Mj+35LX2JgutPqNR39vh477dLbT67JVl7X4nLT67HHwVAgkCsQ7ShTkmzC31+UZ+MlW42OinjHUzi86qAePMmhavnRafg21RD3dtC9Pic9Dqc1Qa+lavA59TtuKFWCkJArEqSmXNdLZYacwn5xv0ZNW8uXy5Lhqb1ULUa6fFZ6cv4ua2vhAtXofZ6JsNvM9ByGWToQ6EWGUSBOKStNYkcrPE0hcb91i6WNl6n0wXiJnTqXRhycVPAF57E1GfnajHzo1dAaNhNxv3Fq/DaPy9djnYKsQ6kiBoMFprkvk5YukCU+kisXSh0rBX5s0GPpYuUiyVl7xHs1UR9diJeu3mQVY/Ua/xPGIubzGfy71jhah/8ivdBEplTTxTXNC4G4/5ZdXzyzfuFgVhj7HlHvHa2dHirWrcbZXGPeKRrXchNhsJgjqVLc4xlS4ylSkyNd/AZwrEUmaDXjUfzxaXXMUKxpZ72G0n4rURdtvZ2eol4rURMZdFPHbCbqOxD7ltclWrEA1KgmCNFOZKxDPFyhZ7ZT5TIF7d4JvLc7NLD6iCcXFTxGNspfeEXdzcGyTqsRH22CvL57fs5cwZIcRKSBC8SfMNe/VjKm1OM0ZXTLxqaz5VWHquO1zcag97bITcNrZGPYTdNkIeY8t9fnnEY8xLn7sQYrVJq4JxADVTLJlb5ubWeqbI9HwDv6jBj2eKpC/RsFsUhNzGlnnIbeOGrgBht814mI159bxctSqEWG9XDAJzBNG/0VpPr0F9aiZXLPE3L51naDrHdHbhFnw8W6Q4t/QAKhjnt4fcRqMe9tjoDbuM5y5jqz3sthFy2ytlAs5mOc9dCLGhrGSPoA3jxvOvAE8Az+iNNmQpcD6e4T89c+KSDf48paAr6GRHi5ftLR66g85KQx/22Ai6bARdzTTJnZ6EEJvEioahVkbfxU8DHwH6ga8CX9Ran6lt9Za6lmGo87MlYukC05nZShdQdTfQ1IL+/sKyY9iAERZ+Z7OxV1C1txAy9w4iHlul3z/sthF02+QWgUKIdXXNw1BrrbVSagwYA+aAIPA1pdSzWuv/sHpVrS1Hs5WuoIuu4MrKz5bKlW6keNroQlrQpWQ+zk9leeXCDNPZ4oKbhVcLuJqNg77zAeFZPjTCHrt0Lwkh1tRKjhH8W+DDQAx4HPgtrfWsUsoCnAI2TBBcrWarxRzvxrGi8uWyJpmfJZY29ijimSKxzMX5+VNHT02keWmgyPQlzv83DjjbFpxNVH1q6Px8xJx32mQETSHEm7eSPYII8D6t9fnqhVrrslLq3bWp1sZksSgCLhsBl43tLZ4rlp8rlZnOzla6omKZInHzWgLj9FPj1NMjI0li6QKpS3RVuW1WMyDmTzO1y7UFQogVu2IQaK1//zLrjq1udRpLkzniZtRrB7xXLF+YK1X2KqbSRWNMoEXDSlyIZ3nlwjTxTHHZQeCqrzaOVILCvmA4iRavnahHhnIWolHIdQQbiL3JWrn5+ZXMDwu9eHC5qUzRHFDOmD8xliKWLjBbWpoaNqulMs5Q9YByC+bNsYncNquEhhAblATBJmW1qMrW/pVcabjpWLrAyEyeQ0OJSw437Wy2LhqgzrHscNNhj13GNBKizkgQCJSqPrZx+S6q6hvQLBi62pyfSBU4O5nhpYH4sncYmx/ltKVqqOoWr8MMDDtRM0CiXrmNpBBrRYJAXJXqPY097Zcvu/CewwUmU/mq21Eady2bPxC+3F6Gz9Fk3obSTqvXcXG+ahr12rE3SWAIcS0kCETN2Jvmr9twXbZcqayZyly8Z/Fk8uKtLSfM+ZcG4kyk8sseywi5jQPc1QHR4nPQ6rXT5jducRl22+RqcCEuQYJArDurRVWu19h7mXLlsmYmN8t4Ms9YMs9EMs94ssC4OZ1I5Tk+lmQytXQPw6Ig4jFCos3voN0MiHa/gzZzWZvfIaO7ioYkf/Viw7BYVGVIjz3tvkuWK5U1U+kCY1VBMWGGx3iywIWpLD8ZiJPILT2G4XM00e530up30O5zGFMzJNrM4JA7tInNRoJAbDpWi6LF7B66nFyxxFgyz2gix3gyz2giz9j8I5nn+GiSyXRhydXf9iYL7X4H7X4n7QEHHX7jlN75+faAA5+juYb/QiFWlwSBaFhOm5W+iJu+iPuSZWZLZSZTBUYT+aqwyDGaMOZfPDPFeKqwZIwpj73JCIuAk46q0OgMOCshIkODiHohQSDEZTRbLVe8iG+uVGbSvNZiZCbHaCLHyIyxpzGayHPUPDNqsZDbRocZDp0BF51BJ50BJ11B4/OCLumCEmujpkGglHoQ+CxgBR7XWn/6EuVuBV4E3q+1/lot6yTEamuyWowtfr+TW3qXH9q2MFdiPFFgJJEzwyLP8EyO4ekcZycz/PBUjGxx4X2qXTbjSvLOgHNJSHQGnLT6HHJxnlgVNQsCpZQV+DzwADCEcXObp7XWR5cp98fAM7WqixDrzd5kpSfsoie8/Km0WmtmsrMMz+QYms5VQmJ4JsvITJ7DwwnimeKC1zRbFR0BJ91BF11BJ92hi9PuoIuIxyZ7FGJFarlHcBtwWmt9FkAp9STwMHB0UbnfAP4euLWGdRGirimlCJo3Mbq+079smWxxjpGqoBiazjEYzzI0neO5Y+PE0guDwtlsrQqG+aAwwqIn7JID2qKilkHQCQxWPR8Cbq8uoJTqBN4L3IcEgRCX5bI1sb3Fe8lhQLLFuUo4DMazDFYFxYFz8SXDmAddzfSE3fSGXPSGXfSEXGyJGM+jXrvsTTSQWgbBcn9Fiy8L/QzwSa116XJ/dEqpR4BHAHp6elarfkJsKi5bEztbvexsXT4oEtlZBqeNkLgQz3I+nuXCVJZXB6f55usjCy7CczZb6QkZXVmVoAi72RJ20RlwylXam0wtg2AI6K563gWMLCrTDzxphkAEeKdSak5r/Q/VhbTWjwGPgXHP4lpVWIjNzO9qxu/yL9v1NFsqMzyd43w8y/mpDOensuYjww9PTZKfLVfKNlsV3SEXW81Tb7eY060RD60+2ZPYiGoZBAeAHUqpPmAY+ADwweoCWuu++Xml1F8C31wcAkKI2mu2WthiNuoQXbBOa81EqsD5qSznYhkGpjLGNGac7VSYuxgSzmYrWyJutkbcbIm46It46Iu42R714HfJMYl6VbMg0FrPKaUexTgbyAo8obU+opT6uLn+C7X6bCHE6lFKmQP6ObitL7RgXbmsGU3mORfLcDaWYWAyw0AszdHRJN85MrbgQruIx872FjfbWzxsi3rY3mI82nwO2YtYZ0ovd/f0Otbf368PHjy43tUQQlzBbKnMYDzL2ckMZybTnJ5Ic9qcVh+4dtusbGvxsD3qMaYtHna2eukJueQ6iVWklHpZa92/3Dq5slgIURPNVgtbox62Rj3cT2tludaayXSB0xNpzkwYwXBmMsOPzkzx9VeHK+UczRZ2tBgHv3e1edjV5mNXq1eOQ9SABIEQYk0pdXHY8bu2RRasS+VnOTOZ4eR4ihNjKU6Op/jhqUn+/pWhShmfo4ldbUZA7G7zsrvdx552Hx67NGdvlnxzQoi64XU0s787wP7uwILl05kiJ8ZTCwLi6UMj/M1LF7uYtoRd7O3wc12Hj+s6fOxt911xBFphkCAQQtS9oNvGHVvD3LE1XFmmtWY0kefYaJKjI0mOjCQ5PJzgW4dHK2UiHrsRCh0+ru/ws6/bT2fAKV1Li0gQCCE2JKVUZWTYt++5eAwimZ/lmBkMR0eN6X//wVnmzDOYwm4b+7oD3NjlZ193gH1dAUJu23r9M+qCBIEQYlPxOZq5fWuY26v2HgpzJU6MpTg0lODQ4AyHBmf45xMTlZsOdYec7OsyuqRu6Q2yt8OPralxrp6W00eFEA0pXZjj8FCC14dmODQ0w6HBBMMzOcA4Y2l/d4Bbt4To3xLi5p4A3g0+SJ+cPiqEEIt47E3cuS3Mndsu7jlMJPMcPD/NgXNxDp6b5vP/fJqyBouC3W0+bt0SNF6zNbKprpSWPQIhhLiEdGGO1y7MGMFwPs4r52fIzZawKLih08/d2yPcsz3Czb1BHM31fevRy+0RSBAIIcQKFefKHBqa4flTMV44HePVwRlKZY29ycJtfSHesjPK/XtazTGb6osEgRBC1EC6MMdPBqZ4/tQUz5+e5OR4GoAdLR7uv66VB65rZX9XAEsdDJUhQSCEEGtgMJ7luWPjPHdsnJfOxpkrayIeOw9c18LD+zu5bUto3UJBgkAIIdZYIjvL909O8OzRcb53fIJssURnwMnD+zt4382dl7zTXK1IEAghxDrKFud49ug4T706zA9PxSiVNTd0+vnFO3p4eH/nmhxoliAQQog6MZkq8I+HRvjqwUGOj6UIu238wh29fOiOXqJee80+V4JACCHqjNaaH5+d4onnB3ju2AQ2q4V/3d/FJx7YScSz+oEgF5QJIUSdUUpx17YId22LcHYyzRefH+BvDwzyjddGePS+7Xzk7i3Ym9bm2oTGGUxDCCHq1Naohz967w0884l7ub0vxKe/fZwH/vQHvD40syafL0EghBB1YlvUwxd/6Va+/Cu3USpr/tUXfszfHRys+edKEAghRJ35qR1Rnn70bvp7g/zW117nz753qqafJ0EghBB1KOyx86Vfvo2f3d/B//fsSX5wcrJmnyVBIIQQdarJauE/vu9GdrR4+PdffY38bKkmnyNBIIQQdcxps/Jb79hNLF3ktcGZmnyGBIEQQtS5W7cEAXj5/HRN3l+CQAgh6tz8ndM6Ao6avL8EgRBC1LmnXhkG4Pa+8BVKvjkSBEIIUce+f2KCx58f4AO3dtMRcNbkMyQIhBCiTj13dJzf+Mqr7G7z8ofv2Vuzz5GxhoQQos6Uy5rPfe8Un3nuFDd0+vlvH7qlpkNVSxAIIUQdefHsFJ/65lGOjCT5uZu7+KP3Xl/z+xXUNAiUUg8CnwWswONa608vWv8LwCfNp2ng17TWh2pZJyGEqEcDsQx//O3jfOfIGB1+B5/7+Zv4mRvbUar2t7asWRAopazA54EHgCHggFLqaa310apiA8BbtNbTSqmHgMeA22tVJyGEqCdaa34yEOfx5wd47tg4zmYr//sDO/novVvX5K5l82q5R3AbcFprfRZAKfUk8DBQCQKt9Y+qyr8IdNWwPkIIURdyxRLfOTLKF58f4I3hJAFXM//rW7fx4bu20OKtzbUCl1PLIOgEqsdPHeLyW/u/Anx7uRVKqUeARwB6enpWq35CCLFmSmXNS2en+Pqrw3znjTHShTm2Rd380Xuv5303deG0rd0ewGK1DILlOraWvS+mUuptGEFwz3LrtdaPYXQb0d/fv7HurSmEaFhaaw4PJ/jW4VG+8eoIY8k8HnsTD13fxntv6uSOrWEsltofA7iSWgbBENBd9bwLGFlcSCl1I/A48JDWeqqG9RFCiJorzJX48Zkpnj06znePTTCWzGO1KN6yM8rvvmsP9+9pXdet/+XUMggOADuUUn3AMPAB4IPVBZRSPcDXgQ9prU/WsC5CCFEzo4kcz5+K8c8nJviXE5NkiiVcNiv37ohy/3Wt3Le7hZDbtt7VvKSaBYHWek4p9SjwDMbpo09orY8opT5urv8C8PtAGPiv5ilSc1rr/lrVSQghVkMiN8uLZ6d44XSM50/HODuZAaDFa+fhmzp5YE8rd24Lr+mZP9dCab2xutz7+/v1wYMH17saQogGkszP8sr5aQ6ci/PC6SleH5qhrMHZbOX2rSHu2R7h7u0RdrV666LPfzlKqZcvtaEtVxYLIcQio4kcB85Nc/BcnAPnpjk+lkRrsFoU+7r8PPq27dy9PcJNPUFsTRt/yDYJAiFEQ8vPljg2muT1oQSvXpjmwLnpyvj/LpuVm3uC/Lu37+DWLSH2dwdw2zdfs7n5/kVCCHEJpbLm7GSa1wZnODQ0w+tDCY6NJpktGV3kUa+dW7cE+ZV7+rh1S4g97V6arBt/i/9KJAiEEJvSXKnM2ViGoyNJjowkODyc4PBQgkzRuAG8x97EjV1+fvWntrKvy8++7gBtPseajO1TbyQIhBAbXrY4x/GxFEdGkhwdSXJ0JMHxsRSFuTIAtiYLe9q8/NwtXezrCrCvO8DWiLtuD+yuNQkCIcSGUS5rBqeznBhLcXI8xfGxFMdGk5yNZZg/AdLvbGZvh48P3dHL3k4f17X72RZ1N0QXz5slQSCEqDtaayZTBU6Mpzgxlqo0/CfH0+RmS5VyXUEne9p9vPvGDvZ2+Njb6afD35jdO9dCgkAIsW601owk8pyZSHN6Is3pSWN6ajzFdHa2Ui7isbOrzcMHbutmd5uXna1edrR68WzCM3jWg3yLQoiaK86VOT+V4YzZ0J+eSHNm0nieLV7cwvc7m9kWdfPg9W3savWys83LrlYvYY99HWu/+UkQCCFWhdaa8WSBs7E052JZBmJpBmIZzsYynJ/KUipfHMWgw+9gW4uHf9PfzfYWD9tbPGyLeoh4bNKtsw4kCIQQV2U6U+RsLMO5WIYB82E09pkFW/e2Jgt9YTc7Wjw8dH1bpbHfFvVsyouyNjL53xBCLKC1ZiJV4PxUlvNTGS7Es8Z83Hg+U9V3b7UouoNO+iJu7twapi/ioi/ioS/qpt3nkNMzNwgJAiEa0GypzPB0jvPxLBemMpWG/sJUlvPxDPnZcqWsRUFn0ElvyM07b2hna8RNn/noDrloltMyNzwJAiE2oXJZM5kuMDSdZTCeYzCeZXB+fjrLaCK/oM/e0WyhJ+SiJ+Tmnh0ResMuekIutoTddAad0thvchIEQmxAWmtmsrMMTefMBn5hQz88natcVTsv6rXTHXRyS2+Q7qCLnrCL3pCL3rCbFq9dunEamASBEHWoXDb66YdnsgzP5Bmezhnz0zmGZ3IMT+cqY+bM8zub6Q452dXq5f49rXQHnXQFXXSHjOlGuUmKWHsSBEKsg8JcidGZvNGomw179XQ0kauMiDkv4GqmM+CkN+zmrm0RuoJOukOuytTnaF6nf43Y6CQIhFhlc6UyE6kCo4kcIzP5BdPRRJ6RmRyxdHHBa5SCVq+DzqCT/d0B3nlDO51BJ10BJ51BJx0Bp1xFK2pG/rKEuArlsmYqU1y2cZ+fTqQKCw7EgjHkcbvfQXvAyXXtPjoCTtr9DrOxd9Hmd2yKO12JjUmCQAhTca7MRCrPWCLPaCLPeNKYjiWNZWOJPBOp/JIuG1uThQ6/g3a/kzu3henwO2kPOOgIOCvz0m0j6pkEgdj0tNakC3OMJwuVxt2Y5hhLFBhLGtNYurDktY5mC+1+J20+B7f1hWjzO2jzOWjzO+jwO+kIOAi5ZVgEsbFJEIgNLT9bYjyZrzTy48k8E6kCY4mL8+PJ/IKhD+YFXM2VRv2GTj+tPgftfoc5NRp/n7NJGnmx6UkQiLqUny0xmSowkSowmapu6AtMpIxGfiyRJ5mfW/Jae5OFVp+DVp+dvR0+7tvdQqvPTqvPQYvXaOzb/A45nVIIkwSBWDNaa1KFOaOBNxv0+cZ+IplnMj2/vEAiN7vk9U0WRYvXTovPURnbpsXnqDT6rT4HrV7ZihfiakkQiGs2VyoTzxSZTBeqtuLn5/OVxn0ilV8whs08W5PFaOC9drZFPdy5LUzUY6fFZ6fF6yDqtdPmdxBy2eTqVyFqQIJALKtc1szkZplMGQdR5xv2ynz64vOpTLFyv9hqXkeT2cA72N8dMLfmjefz81GPbMELsd4kCBpIuaxJ5GaNxjxdYCpdXNrIm8+n0kXmyktbd1uThajHboxbE3JxU0+QqNd4biy3EfUYW/FOm/TBC7ERSBBscLPz3TLzW+dm4z6VKRJLLWzw45nlG/cmiyLssVUa8z1tvouNu9dOxHNx3muXrXchNhsJgjozf0B1Kl0knikQSxeZSheZMhv3yXSBWMps6NOFBTcJqTa/5R7x2Gj3O7i+00fEYzTqYY/NWOe1E3bbCErfuxANTYJgDWSLRsM+lTEb9Or5zMLl8UyRYmnpAVUw+twjZuO+o8XDHVtDlcY94rEtaOg9suUuhFihmgaBUupB4LOAFXhca/3pReuVuf6dQBb4Ja31K7Ws07XSWpMtloibDfi0OY1njEY9vriRTxfJzS69mAmMq1bDbqMRb/U52NPuI+yxEXHbCbltxrzHmA+5bXLeuxCiJmoWBEopK/B54AFgCDiglHpaa320qthDwA7zcTvw5+Z0zcwfQI1nja3x+a3yeKZAPDN7sYGvavQX3/Bjns1qqTTgYY+drVEPYbcxH65aPj/vsskOmRBi/dWyJboNOK21PguglHoSeBioDoKHgS9prTXwolIqoJRq11qPrnZlZrJFPvfd04wlc2ZDbzbu2dklI0XOc9ushDw2Qi4bLV47u9uMLfaQ21gWctsIeWyEzS126Y4RQmxEtQyCTmCw6vkQS7f2lyvTCSwIAqXUI8AjAD09PW+qMj86M8UTLwzQHTLGkOmLuLmlN2h2u9grjXn1Q7pihBCNoJZBsNym8eJN75WUQWv9GPAYQH9///Kb71cwf8HTFz98KztbvW/mLYQQYlOq5Z0whoDuquddwMibKLMqXDYr3SEnTXKapBBCLFDLPYIDwA6lVB8wDHwA+OCiMk8Dj5rHD24HErU4PgDwtt0t/HD3fbV4ayGE2NBqFgRa6zml1KPAMxinjz6htT6ilPq4uf4LwD9hnDp6GuP00Y/Uqj5CCCGWV9PzF7XW/4TR2Fcv+0LVvAZ+vZZ1EEIIcXlyt2whhGhwEgRCCNHgJAiEEKLBSRAIIUSDkyAQQogGJ0EghBANTunlbjZbx5RSk8D5N/nyCBBbxerUgtTx2tV7/aD+61jv9YP6r2O91a9Xax1dbsWGC4JroZQ6qLXuX+96XI7U8drVe/2g/utY7/WD+q9jvdevmnQNCSFEg5MgEEKIBtdoQfDYeldgBaSO167e6wf1X8d6rx/Ufx3rvX4VDXWMQAghxFKNtkcghBBiEQkCIYRocJsmCJRSDyqlTiilTiulfnuZ9Uop9Tlz/etKqZtX+tr1rJ9Sqlsp9c9KqWNKqSNKqX9Xi/pdSx2r1luVUq8qpb5Zj3VUSgWUUl9TSh03v88766x+nzD/j99QSv0PpZRjteu3wjruVkr9WClVUEr95tW8dj3rV2e/lUt+h+b6mv9WrorWesM/MG58cwbYCtiAQ8B1i8q8E/g2xn2S7wBeWulr17l+7cDN5rwXOLna9bvWOlat//fAV4Bv1tv/s7nur4BfNedtQKBe6gd0AgOA03z+VeCX1uk7bAFuBf4I+M2ree0616+efivL1nGtfitX+9gsewS3Aae11me11kXgSeDhRWUeBr6kDS8CAaVU+wpfu27101qPaq1fAdBap4BjGI3GaruW7xClVBfwLuDxGtTtmuuolPIB9wJfBNBaF7XWM/VSP3NdE+BUSjUBLmpz/+4r1lFrPaG1PgDMXu1r17N+9fRbucx3uFa/lauyWYKgExisej7E0j+AS5VZyWvXs34VSqktwE3AS6tcvxV9/hXKfAb4D0C5BnVbyedfqcxWYBL4C3OX/HGllLte6qe1Hgb+BLgAjGLcv/v/X+X6rbSOtXjtSq3KZ9TBb+VyPkPtfytXZbMEgVpm2eLzYi9VZiWvvVbXUj9jpVIe4O+B/01rnVzFuq3o8y9XRin1bmBCa/3y6lfryp+/wjJNwM3An2utbwIywGr3cV/LdxjE2KrsAzoAt1LqF1e5fpf8/DV47Upd82fUyW9l+Reu3W/lqmyWIBgCuqued7F0t/pSZVby2vWsH0qpZow/7L/RWn99leu2GnW8G3iPUuocxm7yfUqpv66zOg4BQ1rr+S3Er2EEQ73U735gQGs9qbWeBb4O3LXK9VtpHWvx2pW6ps+oo9/KpazVb+XqrPdBitV4YGztncXYmpo/eLN3UZl3sfAg3U9W+tp1rp8CvgR8pl6/w0Vl3krtDhZfUx2BHwK7zPk/BP7feqkfcDtwBOPYgMI4sP0b6/EdVpX9QxYejK2L38pl6lc3v5VL1XHRupr9Vq7637TeFVjF/5x3YpwlcAb4XXPZx4GPV/2RfN5cfxjov9xr66V+wD0Yu52vA6+Zj3fWUx0XvUdN/7iv8f95P3DQ/C7/AQjWWf3+T+A48AbwZcC+Tt9hG8ZWbxKYMed9dfRbWbZ+dfZbueR3uFa/lat5yBATQgjR4DbLMQIhhBBvkgSBEEI0OAkCIYRocBIEQgjR4CQIhBCiwUkQCCFEg5MgEA1LKfUxpdQXzPlmpdSXlVJ/ZV6dKkTDkCAQjexG4HVzZNJvAxe01h/WxhAPQjQMCQLRyG4A4sD3gb/TWv/u+lZHiPUhVxaLhqWUmgbmgF/WWv/jetdHiPUiewSiISmluoE0xng/7VcovpL3U0qpTy1a9g6l1IcWL7/K9/3La62bEFfStN4VEGKd3IgxauRHgReVUge01q8qpdqAvwW+BewFfgQ8gDGK5AmMgeFcGKNO/l8Yw0U/jTFaaJNSygb8KcZgY7cDvwtcp5RyYty20A5Ma61/z7xT1ZfN19+htX6/UuoPgBDGQGV/jHHfBCFqSvYIRKO6ATistR4FfhX4W6WUH+OuVl/XWv8nwA/8d+DvgF7gEcCJ0Uh7zLJPaq0/jTGy6WvArwF/pbX+HYxGP2gu/w3gK1rrXwd2m3XYB/yD1vo/A3NKqU6g2Xz/OzDul/BKrb4AIeZJEIhGdQNGtxBa62cxbhb/BEaD/ox5CumU1roMXG+WvQn4ba31H2qt/xez7LPm++3HaPBvAg4rpbxArGr5XnO5Dciar9kHPGPOa+BTGHsBfwUMY9z8/MAq/7uFWEK6hkRD0lr/wqLnvweglPoixjjzN2Lc/Bxgi9b6glLqG8BfKqUGge8BOzC6iwC2A6cwGvYvYDT2J6uWfxV4zFz+H6tec1IpFQHGMO5V/JtAGHgVIyg+t6r/cCGWIWcNCSFEg5OuISGEaHASBEII0eAkCIQQosFJEAghRIOTIBBCiAYnQSCEEA1OgkAIIRqcBIEQQjQ4CQIhhGhw/xM2WZvZKfQ5rgAAAABJRU5ErkJggg==\n", "text/plain": [ "