{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "\"Open" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Eigenvalue Problems -- Shallow Water on the Sphere" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Overview:** This notebook describes how to solve an eigenvalue problem to find the fastest growing mode in a shear flow on the sphere, and initialize a simulation with that mode.\n", "\n", "**About Dedalus:** [Dedalus](http://dedalus-project.org) is an open-source Python package for solving partial differential equations (PDEs) using global spectral methods.\n", "These methods provide highly accurate numerical solutions for PDEs with smooth solutions in simple domains like boxes and spheres.\n", "Dedalus implements modern parallel algorithms utilizing sparse polynomial bases, but all with an easy-to-use symbolic interface.\n", "The code is being used in a wide range of fields, often for problems involving fluid dynamics.\n", "\n", "**Author:** [Keaton Burns](http://keaton-burns.com)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Setup" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "This cell checks if Dedalus is installed and performs some other basic setup.\n", "\n", "If Dedalus is not installed and you are using Google Colab, it will automatically be installed.\n", "This may take a few minutes the first time you run the notebook, but subsequent sessions during the next day or so should have the installation cached.\n", "No need to worry about the details -- just execute the cell.\n", "\n", "If you are not using Google Colab, follow the installation instructions in the [Dedalus Docs](https://dedalus-project.readthedocs.io/en/latest/pages/installation.html) to install Dedalus locally on your computer.\n", "Installation using conda is typically straightforward for Mac and Linux.\n", "No promises on Windows.\n", "Execute the cell to confirm Dedalus is installed and importable." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "env: OMP_NUM_THREADS=1\n", "env: NUMEXPR_MAX_THREADS=1\n", "Dedalus already installed :)\n" ] } ], "source": [ "# Set environment variables for best performance\n", "%env OMP_NUM_THREADS=1\n", "%env NUMEXPR_MAX_THREADS=1\n", "\n", "# Minimize logging output\n", "import logging\n", "logging.disable(logging.DEBUG)\n", "\n", "# Check if running on google colab\n", "import os\n", "using_google_colab = bool(os.getenv(\"COLAB_RELEASE_TAG\"))\n", "\n", "# Check for Dedalus\n", "try:\n", " import dedalus.public as de\n", " print(\"Dedalus already installed :)\")\n", "except:\n", " print(\"Dedalus not installed yet.\")\n", " if using_google_colab:\n", " print(\"Installing for Google Colab.\")\n", " print()\n", " # Step 1: Install FFTW\n", " !apt-get install libfftw3-dev\n", " !apt-get install libfftw3-mpi-dev\n", " # Step 2: Set paths for Dedalus installation\n", " import os\n", " os.environ['MPI_INCLUDE_PATH'] = \"/usr/lib/x86_64-linux-gnu/openmpi/include\"\n", " os.environ['MPI_LIBRARY_PATH'] = \"/usr/lib/x86_64-linux-gnu\"\n", " os.environ['FFTW_INCLUDE_PATH'] = \"/usr/include\"\n", " os.environ['FFTW_LIBRARY_PATH'] = \"/usr/lib/x86_64-linux-gnu\"\n", " # Step 3: Install Dedalus using pip\n", " !pip3 install cython \"mpi4py<4.0\" numpy setuptools wheel\n", " !CC=mpicc pip3 install --no-cache --no-build-isolation http://github.com/dedalusproject/dedalus/zipball/master/\n", " !pip3 install -q ipympl\n", " # Step 4: Check installation\n", " print()\n", " try:\n", " import dedalus.public as de\n", " print(\"Dedalus successfully installed :)\")\n", " except:\n", " print(\"Error installing Dedalus :(\")\n", " raise\n", " else:\n", " print(\"See website for installation instructions:\")\n", " print(\"https://dedalus-project.readthedocs.io/en/latest/pages/installation.html\")\n", "\n", "# Setup interactive matplotlib\n", "if using_google_colab:\n", " from google.colab import output\n", " output.enable_custom_widget_manager()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Content" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "First let's import everything we need to run the rest of the notebook." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "np.seterr(over=\"raise\")\n", "import matplotlib.pyplot as plt\n", "import dedalus.public as d3\n", "import logging\n", "logger = logging.getLogger(__name__)\n", "%matplotlib widget" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Let's consider the shallow water equations on the sphere, with $u$ being the horizontal fluid velocity, and the total fluid depth being $H + h$ for constant $H$.\n", "\n", "$$\\partial_t u + u \\cdot \\nabla u = - g \\nabla h - 2 \\Omega \\times u -\\nu \\nabla^4 u$$\n", "$$\\partial_t h + \\nabla \\cdot ((H + h) u) = -\\nu \\nabla^4 h$$\n", "\n", "Here we've adding a regularizing hyperviscosity to both equations.\n", "[Gelewsky et al. (2004)](https://doi.org/10.3402/tellusa.v56i5.14436) is a classic test problem for spherical shallow-water solvers which examines the evolution of a balanced, barotropically unstable, mid-latitude jet under a prescribed perturbation.\n", "Here we'll go a bit farther and find and evolve the most unstable mode of the jet using the eigenvalue problem in Dedalus. \n", "\n", "This will give us three problems in total:\n", "* An LBVP to find the height field balancing the prescribed jet profile.\n", "* An EVP to find the most unstable mode of the jet.\n", "* An IVP to examine the nonlinear evolution of the instability." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Setup domain" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "First let's set some basic parameters and build the domain." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Simulation units\n", "meter = 1 / 6.37122e6\n", "hour = 1\n", "second = hour / 3600\n", "\n", "# Parameters\n", "Nphi = 256\n", "Ntheta = 128\n", "dealias = (3/2, 3/2)\n", "R = 6.37122e6 * meter\n", "Omega = 7.292e-5 / second\n", "g = 9.80616 * meter / second**2\n", "H = 1e4 * meter\n", "dtype = np.float64" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now let's build two bases for the sphere. First, one to just represent zonally constant fields, and second, one to represent full 2D fields on the sphere:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Bases\n", "coords = d3.S2Coordinates('phi', 'theta')\n", "dist = d3.Distributor(coords, dtype=dtype)\n", "full_basis = d3.SphereBasis(coords, (Nphi, Ntheta), radius=R, dealias=dealias, dtype=dtype)\n", "zonal_basis = d3.SphereBasis(coords, (1, Ntheta), radius=R, dealias=dealias, dtype=dtype)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Balanced zonal jet (LBVP)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We start with prescribing a mid-latitude zonal jet, and solving a zonally-symmetric LBVP to find the height field that balances this jet profile.\n", "This will be the background state of our eigenvalue problem.\n", "\n", "First we build the background fields using the zonally symmetric basis:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Backgroudn fields\n", "u0 = dist.VectorField(coords, name='u', bases=zonal_basis)\n", "h0 = dist.Field(name='h', bases=zonal_basis)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next let's setup the zonal jet (the details here are from the reference above):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Setup zonal jet\n", "phi, theta = dist.local_grids(zonal_basis)\n", "lat = np.pi / 2 - theta + 0*phi\n", "umax = 80 * meter / second\n", "lat0 = np.pi / 7\n", "lat1 = np.pi / 2 - lat0\n", "en = np.exp(-4 / (lat1 - lat0)**2)\n", "jet = (lat0 <= lat) * (lat <= lat1)\n", "u_jet = umax / en * np.exp(1 / (lat[jet] - lat0) / (lat[jet] - lat1))\n", "u0['g'][0][jet] = u_jet" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can now solve for the balancing height field (ignoring hyperdiffusivity).\n", "This comes from taking the divergence of the momentum equation, and using a gauge freedom to fix the mean of $h$ to be 0:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2024-09-27 10:17:56,462 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 9.9e+01/s\n" ] } ], "source": [ "# Substitutions\n", "zcross = lambda A: d3.MulCosine(d3.skew(A))\n", "\n", "# Find balanced height field\n", "c = dist.Field(name='c')\n", "problem = d3.LBVP([h0, c], namespace=locals())\n", "problem.add_equation(\"g*lap(h0) + c = - div(u0@grad(u0) + 2*Omega*zcross(u0))\")\n", "problem.add_equation(\"ave(h0) = 0\")\n", "solver = problem.build_solver()\n", "solver.solve()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot these backgrounds:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "91cd319f0a074f4abfdd69a23b99b23f", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABc80lEQVR4nO3deXwTdf4/8FeaNukdetALSrlEQEAQEIoiIMsloqjroqvIpeICCwgsyhdZTgFBERdRFFkOT34KuoKuHAKKUpSjKvdyFQptgV4pvdIjn98foaFpjqbDpMkkr+fjkUebTyYznwntp2/mNZ8ZlRBCgIiIiIh8hp+7O0BERERE9YsFIBEREZGPYQFIRERE5GNYABIRERH5GBaARERERD6GBSARERGRj2EBSERERORjWAASERER+RgWgEREREQ+hgUgERERkY9hAUhERETkY1gAEhEREfkYFoBEREREPoYFIBEREZGPYQFIdaJSqWp9zJkzx93dRFpaGlQqFdatW1fv23b1Z1Bz/cePH8ecOXOQlpbmsm0Sudq6deugUqlw8OBBm68/+OCDaNq0qaR1jxw5UvJ758yZA5VKhezs7FqXXbhwIb766itJ23E1V48Ttf37SbFv3z7MmTMH+fn5Vq/17t0bvXv3lm1bvsjf3R0gZUlJSbHZXlFRgWeeeQaXL1/GAw88UM+98i0pKSlo3Lix+fnx48cxd+5c9O7dW/IfOSJvNmvWLEyaNMnl21m4cCH+/Oc/Y+jQoS7fVl0pcZzYt28f5s6di5EjR6JBgwYWr73zzjvu6ZQXYQFIddK9e3eb7RMnTsT58+fx3nvv4e67767nXvkWe/8GRGRbixYt3N0FtykvL4dKpVLs+u1p27ZtvW/T2zACplv24YcfYsWKFRgzZgyef/55i9dyc3Mxbtw4NGrUCBqNBs2bN8fMmTNhMBgsllOpVJgwYQI+/PBDtGnTBsHBwbjzzjuxdetWi+XOnDmDUaNG4bbbbkNwcDAaNWqEIUOG4MiRI3Xu97Vr16DRaDBr1iyr106ePAmVSoV//etf5rasrCyMHTsWjRs3hkajQbNmzTB37lxUVFTUuq2jR4/i4YcfRkREBAIDA9GxY0esX7/earn8/HxMnToVzZs3h1arRUxMDB544AGcPHnSvEz1CHjdunV4/PHHAQB9+vQxx/Dr1q3D/Pnz4e/vj/T0dKvtjB49GlFRUSgtLa2170SeSgiBd955Bx07dkRQUBAiIiLw5z//GefOnbNYzlYEnJ+fjzFjxiAyMhKhoaEYPHgwzp07Z/cUjitXruDJJ5+ETqdDbGwsRo8eDb1eb35dpVKhqKgI69evN/8eOoooq05TWbJkCV599VU0adIEgYGB6NKlC77//nur5U+fPo2//vWviImJgVarRZs2bbBy5UqLZfbs2QOVSoUPP/wQU6dORaNGjaDVavHBBx/YHScAoGnTphg5cqTVNmvGrPbWf+bMGfMyeXl5GDVqFCIjIxESEoIhQ4ZY/Xvs2LEDDz/8MBo3bozAwEC0bNkSY8eOtYjZ58yZg3/84x8AgGbNmpn7vGfPHpt9A+T/e+PteASQbklqairGjh2Lrl27Wg1GpaWl6NOnD86ePYu5c+eiQ4cO2Lt3LxYtWoTffvsN33zzjcXy33zzDQ4cOIB58+YhNDQUS5YswSOPPIJTp06hefPmAICMjAxERUVh8eLFaNiwIXJzc7F+/Xp069YNqampuP32253ue8OGDfHggw9i/fr1mDt3Lvz8bv5/aO3atdBoNHjqqacAmIq/u+++G35+fvjnP/+JFi1aICUlBQsWLEBaWhrWrl1rdzunTp1Cjx49EBMTg3/961+IiorCRx99hJEjR+LKlSuYPn06AOD69eu49957kZaWhpdeegndunVDYWEhfvzxR2RmZqJ169ZW6x48eDAWLlyI//u//8PKlStx1113ATAd8RBC4NVXX8V7772HBQsWmN+Tm5uLzz77DBMmTEBgYKDTnxdRfaisrLT5nyohhFXb2LFjsW7dOkycOBGvvfYacnNzMW/ePPTo0QO///47YmNjbW7DaDRiyJAhOHjwIObMmYO77roLKSkpGDhwoN1+PfbYYxg2bBjGjBmDI0eOYMaMGQCAf//73wBMp2bcf//96NOnj/k/leHh4bXu79tvv42kpCQsX74cRqMRS5YswaBBg/DDDz8gOTkZgCm+7dGjB5o0aYI33ngDcXFx2LZtGyZOnIjs7GzMnj3bYp0zZsxAcnIyVq1aBT8/P3Tp0gV5eXk2xwkpaq4/JibG/NqYMWPQr18/fPLJJ0hPT8crr7yC3r17448//jDHuGfPnkVycjKeffZZ6HQ6pKWlYdmyZbj33ntx5MgRBAQE4Nlnn0Vubi5WrFiBzZs3Iz4+HoD9I3+u+Hvj9QSRRNeuXRNJSUmiYcOG4uLFi1avr1q1SgAQ/+///T+L9tdee00AENu3bze3ARCxsbGioKDA3JaVlSX8/PzEokWL7PahoqJClJWVidtuu028+OKL5vbz588LAGLt2rUO9+Hrr7+26ktFRYVISEgQjz32mLlt7NixIjQ0VFy4cMHi/a+//roAII4dO2axL7NnzzY/f+KJJ4RWq7X6jAYNGiSCg4NFfn6+EEKIefPmCQBix44dDvtcc/2ff/65ACB2795tteyIESNETEyMMBgM5rbXXntN+Pn5ifPnzzvcDlF9Wrt2rQDg8JGUlGRePiUlRQAQb7zxhsV60tPTRVBQkJg+fbq5bcSIERbv/eabbwQA8e6771q8d9GiRVa/X7NnzxYAxJIlSyyWHTdunAgMDBRGo9HcFhISIkaMGOHU/laNUQkJCaKkpMTcXlBQICIjI8Wf/vQnc9uAAQNE48aNhV6vt1jHhAkTRGBgoMjNzRVCCLF7924BQNx3331W23M0TiQlJdnsd69evUSvXr3Mzx2tv+rf75FHHrFo//nnnwUAsWDBApufg9FoFOXl5eLChQsCgPjPf/5jfm3p0qUCgM2xqmbf6uPvjbdhBEySVFZW4oknnsClS5ewceNGJCYmWi2za9cuhISE4M9//rNFe1XUUDPm6NOnD8LCwszPY2NjERMTgwsXLpjbKioqsHDhQrRt2xYajQb+/v7QaDQ4ffo0Tpw4Uef9GDRoEOLi4iyO4G3btg0ZGRkYPXq0uW3r1q3o06cPEhISUFFRYX4MGjQIAPDDDz/Y3cauXbvQt29fq89o5MiRKC4uNk+s+e9//4tWrVrhT3/6U533w55Jkybh6tWr+PzzzwGYjny8++67GDx4sGJOBCffsmHDBhw4cMDqce+991ost3XrVqhUKjz99NMWv5NxcXG48847zVGhLVW/r3/5y18s2p988km773nooYcsnnfo0AGlpaW4evVqHffQ0qOPPmpxJD4sLAxDhgzBjz/+iMrKSpSWluL777/HI488guDgYIt9feCBB1BaWor9+/dbrPOxxx67pT7VxtH6q1KTKj169EBSUhJ2795tbrt69SpeeOEFJCYmwt/fHwEBAUhKSgIASeM44Jq/N96OETBJMn36dHz//fd4/fXX0adPH5vL5OTkIC4uzuoE4ZiYGPj7+yMnJ8eiPSoqymodWq0WJSUl5udTpkzBypUr8dJLL6FXr16IiIiAn58fnn32WYvlnOXv74/hw4djxYoVyM/PR4MGDbBu3TrEx8djwIAB5uWuXLmCLVu2ICAgwOZ6HF0iIicnxxxfVJeQkGB+HTCdk9ikSZM674MjnTp1Qs+ePbFy5Uo89dRT2Lp1K9LS0vDee+/Juh0iubRp0wZdunSxatfpdBbns165cgVCCLsxr6MYLycnB/7+/oiMjLRot7cuwHp80mq1ACBp3KkuLi7OZltZWRkKCwtRWFiIiooKrFixAitWrLC5jprjj63xRk6O1m9vf6rGOaPRiP79+yMjIwOzZs1C+/btERISAqPRiO7du0v+PF3x98bbsQCkOvv000+xbNkyDBs2DFOnTrW7XFRUFH755RcIISx+Ka9evYqKigpER0fXedsfffQRnnnmGSxcuNCiPTs72+oyAc4aNWoUli5dis8++wzDhg3D119/jcmTJ0OtVpuXiY6ORocOHfDqq6/aXEdVMWdLVFQUMjMzrdozMjLM6wZM5yReunRJ0j44MnHiRDz++OM4fPgw3n77bbRq1Qr9+vWTfTtE9Sk6OhoqlQp79+41F2PV2WqrEhUVhYqKCuTm5loUgVlZWS7pqyO2tpmVlQWNRoPQ0FAEBARArVZj+PDhGD9+vM11NGvWzOJ5XWflBgYGWk2UAEzjqq1x2tH67e1Py5YtAZgmxP3+++9Yt24dRowYYV6m+kQSKVzx98bbMQKmOvnjjz/w7LPPol27dlizZo3DZfv27YvCwkKrC6Nu2LDB/HpdqVQqq4H9m2++weXLl+u8ript2rRBt27dsHbtWnzyyScwGAwYNWqUxTIPPvggjh49ihYtWqBLly5WD0cFYN++fbFr1y5zwVdlw4YNCA4ONl/WZdCgQfjf//6HXbt21an/tR2JeOSRR9CkSRNMnToVO3fuxLhx49xy2QYiOT344IMQQuDy5cs2fyfbt29v9729evUCAGzcuNGi/bPPPrulPkk5grR582aL2fjXr1/Hli1b0LNnT6jVagQHB6NPnz5ITU1Fhw4dbO6rraNZtvoG2B4nmjZtij/++MOi7X//+x9OnTpVp30BgI8//tji+b59+3DhwgXzjN2qsafmOG4rlajLUVZX/L3xdjwCSE7Ly8vD0KFDYTAY8NJLL9m99ErDhg3RokULPPPMM1i5ciVGjBiBtLQ0tG/fHj/99BMWLlyIBx54QNK5bg8++CDWrVuH1q1bo0OHDjh06BCWLl1qcWFkKUaPHo2xY8ciIyMDPXr0sJpNPG/ePOzYsQM9evTAxIkTcfvtt6O0tBRpaWn49ttvsWrVKrt9mD17tvkcwn/+85+IjIzExx9/jG+++QZLliyBTqcDAEyePBkbN27Eww8/jJdffhl33303SkpK8MMPP+DBBx+0G7W3a9cOAPD+++8jLCwMgYGBaNasmfmPglqtxvjx4/HSSy8hJCTE5uUeiJTmnnvuwfPPP49Ro0bh4MGDuO+++xASEoLMzEz89NNPaN++Pf72t7/ZfO/AgQNxzz33YOrUqSgoKEDnzp2RkpJiLhaqXxGgLtq3b489e/Zgy5YtiI+PR1hYWK1XJlCr1ejXrx+mTJkCo9GI1157DQUFBZg7d655mbfeegv33nsvevbsib/97W9o2rQprl+/jjNnzmDLli1O/afR0TgxfPhwPP300xg3bhwee+wxXLhwAUuWLEHDhg3r/BkcPHgQzz77LB5//HGkp6dj5syZaNSoEcaNGwcAaN26NVq0aIGXX34ZQghERkZiy5Yt2LFjh9W6qor4t956CyNGjEBAQABuv/12i3P3qrji743Xc+sUFFKUqhlgtT2qzybLyckRL7zwgoiPjxf+/v4iKSlJzJgxQ5SWllqsG4AYP3681TZrzk7Ly8sTY8aMETExMSI4OFjce++9Yu/evVYzwpydBVxFr9eLoKAgAUCsXr3a5jLXrl0TEydOFM2aNRMBAQEiMjJSdO7cWcycOVMUFhZa7Ev1WYRCCHHkyBExZMgQodPphEajEXfeeafNvuXl5YlJkyaJJk2aiICAABETEyMGDx4sTp486XD9y5cvF82aNRNqtdrmfqelpQkA4oUXXnDq8yCqb1WzSA8cOGDz9cGDB1vM5K3y73//W3Tr1k2EhISIoKAg0aJFC/HMM8+IgwcPmpepOQtYCCFyc3PFqFGjRIMGDURwcLDo16+f2L9/vwAg3nrrLfNyVbOAr127ZrO/1Weo/vbbb+Kee+4RwcHBAoDFmFRT1Rj12muviblz54rGjRsLjUYjOnXqJLZt22Zz+dGjR4tGjRqJgIAA0bBhQ9GjRw+L2bVVY/Tnn39uc5v2xgmj0SiWLFkimjdvLgIDA0WXLl3Erl277M4CtrX+qs9j+/btYvjw4aJBgwYiKChIPPDAA+L06dMWyx4/flz069dPhIWFiYiICPH444+Lixcv2hzbZsyYIRISEoSfn5/FLOaafRNC/r833k4lhI2LKxGRV1mxYgUmTpyIo0eP4o477nB3d4g80ieffIKnnnoKP//8M3r06OHSbaWlpaFZs2ZYunQppk2b5tJtEdnCCJjIi6WmpuL8+fOYN28eHn74YRZ/RDd8+umnuHz5Mtq3bw8/Pz/s378fS5cuxX333efy4o/IE7AAJPJijzzyCLKystCzZ0+sWrXK3d0h8hhhYWH47LPPsGDBAhQVFSE+Ph4jR460uGsOkTdjBExERETkY3gZGCIiIiIfwwKQiIiIyMewACQiIiLyMZwE4gSj0YiMjAyEhYXxDgpEXkYIgevXryMhIUHyBYA9EcctIu92q2MXC0AnZGRkIDEx0d3dICIXSk9Pv+U7yngSjltEvkHq2MUC0AlVt51JT09HeHi4m3tDRHIqKChAYmKizdtLKRnHLSLvdqtjFwtAJ1TFJ+Hh4RxIibyUt8WkHLeIfIPUsct7TnghIiIiIqewACQiIiLyMSwAiYiIiHwMC0AiIiIiH8MCkIiIiMjHsAAkIiIi8jEsAImIiIh8DAtAIiIiIh/DApCIiIjIx7AAJCIiIvIxLACJiIiIfAwLQCIiIiIfwwKQiIiIyMewAJRZpr4E+85mI1Nf4u6uEBE5T38ZOP+j6SsReT1/d3fAm3y8/wJe+c9RCAH4qYBFj7bHsK5N3N0tIiLH9q0AdswChABUfsCQt4C7nnF3r4jIhXgEUCaZ+hLMulH8AYBRAP+3+SiPBBKRZ9NfBrbfKP4AQBiBLZN5JJDIy7EAlMn57CIYhWVbpRBIyy52T4eIiJyRexZAjcFLVAK559zSHSKqHywAZdIsOgSqGm0qFdA0Otgt/SEickpkC+s2lQqIbF7/fSGiesMC0JVE7YsQEXkcjl1EXo8FoEzOZxdZjZkCYARMRJ4t96yNRsEImMjLsQCUCSNgIlIkRsBEPokFoCsxRiEiJeLYReT1WADKhBEwESkSI2Ain8QCUCYhGrXN9mANP2Ii8mABIXbaefoKkTdjdSKTorJKm+3FZcZ67gkRUR2UF9lpZ3pB5M1YAMqEk0CISJE4CYTIJ7EAdCWeSE1ESsSxi8jrsQCUCSeBEJEicRIIkU9iASgTRsBEpEiMgIl8EgtAV2KMQkRKxLGLyOuxAJQJI2AiUiRGwEQ+iQWgTBgBE5EiMQIm8kksAF2JMQoRKRHHLiKvxwJQJoyAiUiRGAET+SQWgDJhBExEisQImMgnsQB0JcYoRKREHLuIvB4LQJkwAiYiRWIETOSTWADKhBEwESkSI2Ain8QC0JUYoxCREnHsIvJ6LABlwgiYiBSJETCRT1J0Adi0aVOoVCqrx/jx4wEAI0eOtHqte/fuLukLI2AicpYnjV2MgIl8k7+7O3ArDhw4gMrKSvPzo0ePol+/fnj88cfNbQMHDsTatWvNzzUaTf11kDEKEdnAsYuI3E3RBWDDhg0tni9evBgtWrRAr169zG1arRZxcXEu74ujCDheF+Ty7RORcnjS2OUwAtY1cv32icgtFB0BV1dWVoaPPvoIo0ePhkp1M4zds2cPYmJi0KpVKzz33HO4evWqS7bPCJiIpHD32MUImMg3KfoIYHVfffUV8vPzMXLkSHPboEGD8PjjjyMpKQnnz5/HrFmzcP/99+PQoUPQarV212UwGGAwGMzPCwoKpHWKMQoR1UKusUu2cQvg2EXkA7ymAFyzZg0GDRqEhIQEc9uwYcPM37dr1w5dunRBUlISvvnmGzz66KN217Vo0SLMnTu3TttnBExEUsg1dkkZtwAwAibyUV4RAV+4cAE7d+7Es88+63C5+Ph4JCUl4fTp0w6XmzFjBvR6vfmRnp5eax8YARNRXck5dkkZtwAwAibyUV5xBHDt2rWIiYnB4MGDHS6Xk5OD9PR0xMfHO1xOq9U6jIidxhiFiByQc+ySbdwCOHYR+QDFHwE0Go1Yu3YtRowYAX//m/VsYWEhpk2bhpSUFKSlpWHPnj0YMmQIoqOj8cgjj8jeD14ImojqwlPGLl4Imsg3Kf4I4M6dO3Hx4kWMHj3aol2tVuPIkSPYsGED8vPzER8fjz59+mDjxo0ICwuTvR9VEXD1IpARMBHZ4yljFyNgIt+k+AKwf//+EMI6rwgKCsK2bdvc0KNqGKMQkR0cu4jInRQfAXsKRsBEpEiMgIl8EgtAmXAWMBEpEiNgIp/EAtCVGKMQkRJx7CLyeiwAZcIImIgUiREwkU9iASiTEI3aZnuwhh8xEXmwgBA77Tx9hcibsTqRSVFZpc324jJjPfeEiKgOyovstDO9IPJmLABlwkkgRKRInARC5JNYALoST6QmIiXi2EXk9VgAyoSTQIhIkTgJhMgnsQCUCSNgIlIkRsBEPokFoCsxRiEiJeLYReT1WADKhBEwESkSI2Ain8QCUCaMgIlIkRgBE/kkFoCuxBiFiJSIYxeR12MBKBNGwESkSIyAiXwSC0CZMAImIkViBEzkk1gAuhJjFCJSIo5dRF6PBaBMGAETkSIxAibySSwAZcIImIgUiREwkU9iAehKjFGISIk4dhF5PRaAMmEETESKxAiYyCexAJQJI2AiUiRGwEQ+iQWgKzFGISIl4thF5PVYAMqEETARKRIjYCKfxAJQJoyAiUiRGAET+SQWgK7EGIWIlIhjF5HXYwEoE0bARKRIjICJfBILQJkwAiYiRYpsAdQcvRgBE3k9FoCuxBiFiBRBOHxKRN6HBaBMGAETkSIxAibySSwAZdIsOgR+NTNgAH9cyq/3vhAROY0RMJFPYgEok3hdEF4a1NqqffF3J5GRX4xMfQn2nc1Gpr7EDb0jInKEETCRr/F3dwe8SftGOqs2IYCH3v4ZOUVlEALwUwGLHm2PYV2buKGHREQ1OIqAdY3qvTtEVD94BFBG9mLg7EJT8QcARgHM2HwEl/N4biAReYCAEDvtvIIBkTdjASijeF0QFj3aHmqVqQpUq1R4qpv1kT6jAPot+wFPf/ALlnx3EtuOZSFLX2pznYyOicilyotst+dfrN9+EFG9YgQss2Fdm+C+Vg2Rll1svgbgp79ehLHGOTXF5Ub8dCYbP53JNrfFhGlxZ2ID3NlYhzsTG+DM1ULM33ocRkbHROQqkS0AlR8gjJbtX4wCyq4Ddz3jnn4RkUuphBA83bcWBQUF0Ol00Ov1CA8Pr/P7Nx64iP/bfBSVQkCtUmH+0DvQoXED/H4pH3+k6/H7pXz878p1qyKxJpUKWPxoB9yREI54XSAiQzRQqWxkzkTktFv9/fZUddqvwxuALZOsi0CVGph8hOcCEnmgWx27eASwHtQ8KhivCwIAtGukw1PdTMsUl1XgWEYBfk/Px++X9Pj1fA6uFBgs1iME8NKmP8zPtf5+iNcFIk4XiARdEOJ0gYhvEIT48EDENwhEvC4IEcEBThWJmfoSnM8uQrPoEHP/iMhH3PUMoAk1HfWrTlQC6b8CuVGmI4UsBIm8BgvAehKvC3JYWAVr/NG1aSS6No0EYCrI7lm8y+qoYOu4MGQXliG70ABDhRFpOcVIy7E/oaSqSDRt31QYxumCkFCtbduxLPzfl0cYNRP5ssRudqLgkaavKj9gyFuMhIm8BAtAD1U1oaR6dLzw0XbmwsxQUYmrBQZk5JcgU19641Fi/pqlL0V2YZlTRWJ1RgG8vOkITmZeR5OoYESGaNAgWIPIYA0iQgIQEaxBsEbN6JnI2+gamQq8LZNNR/6ggsUFAYXRFBO36MsjgURegOcAOsGd5whl6kusomNnlZbfKBJvFITmr/mlyCooQWZ+KXKKyurcJ42/HyKDNWgQHIDIEA0iQjSICA64USRqEBF8sy0iWIPIEOlFI6NpcjWeA1iD/rLpGoBF16wjYcAUBXf8K9D2YSD6tmrvOcuYmKge3erYpegCcM6cOZg7d65FW2xsLLKysgAAQgjMnTsX77//PvLy8tCtWzesXLkSd9xxR522461/IAAgLbsI97+xxyJqVgEY2ikBZZUCeUVlyC0qQ35xOXKLy1BWYbS7Lkc0aj/zEcSqorDmc3NBeaOA3Pp7BqNpcjl3/H7Xx9h1y/ulvwwsb2cdCVfXsA0QkQSc3m5ajjExUb3x+Ukgd9xxB3bu3Gl+rlarzd8vWbIEy5Ytw7p169CqVSssWLAA/fr1w6lTpxAWFuaO7nqcptEhDqPm6oQQKCmvRG5RGfKKypFXXIa8YlOBmFdcbioWi8uQX1yG3KKbz8sqjCirNOJKgcFqYouzjAJ4adMR7Dl1DfG6IIQH+UMXFABdUADCAwOgC672fVAAAgP8GFOTR/P4satmJKxSA/0XAJoQ4MTXwLkfgGsnTI8qwgh8PRHQ6oCW9wNajrNEnkrxBaC/vz/i4uKs2oUQWL58OWbOnIlHH30UALB+/XrExsbik08+wdixY+u7qx7L3izlmlQqFYI1/gjW+KNxhHPrrl405heX3ygWy24Uh+U3isWqtnJzQWmwc6Txv0eznNquRu2H8CB/hN8oEqsXh7qgAKsC0rxccABCNf7ws3VLFycxtiZnKGLsuusZ0zl/ueeAyOY3493OI4CSfODnt4CfltXcA+DzG0cAo1oC8XfeeHQE4jsAQTUGD8bHRG6h+ALw9OnTSEhIgFarRbdu3bBw4UI0b94c58+fR1ZWFvr3729eVqvVolevXti3bx8LwBpqm6UslZSiEQDOXytC32U1omkV8EKvFhACKCgth76kHAU3Hvobj4LSClQaBcoqjTdmS9f9HEc/FRAWaLtQNLUFWBSWptdMy207loVXvjrK2JpqpZixS9fIdmEW1ADo+izw8/IaMbEKCI0FCrOAnDOmx9FNN19ukAQkdDQVhUXXgF/eY3xM5AaKLgC7deuGDRs2oFWrVrhy5QoWLFiAHj164NixY+ZzaWJjYy3eExsbiwsXLjhcr8FggMFwM6osKCiQv/PkULOGzkfT1QkhUFRWaSoIi8vNhWLNYrGgtMLcrq/Wbqgwwihgbr8VVbH1lt8zERseaFFMMrr2ba4Yu9wybtmKiYcsNxVxRdlA5u9A5m83vv4O5KUB+RdMj+P/sVxXVXycfdpUIEY0AyKbWR8xJCJZKLoAHDRokPn79u3bIzk5GS1atMD69evRvXt3ALD6YyqEqPUP7KJFi6xO0Kb652w0XZ1KpUKo1h+hWn80alD3I5ql5ZUoKK1WKJbYLhRNRWQ59CUV5sLyuqHC5jqr3+6vNlKj6/CgAIRppUfXjK3rlyvGLreNW/Zi4pBooGVf06NKSR6Q+YepGDy9A0j7scbKBLDvX5ZNgQ1MhWBk85tFYdXX0DjAz8Yt7RkrE9VK0QVgTSEhIWjfvj1Onz6NoUOHAgCysrIQHx9vXubq1atW/7OuacaMGZgyZYr5eUFBARITE13SZ3LMVdG0PYEBagQGqBETFljn96bnFqHXUsvY2k8F/GNAawCoVjR6RnRdFVvvPnkVr357grG1G8kxdrl13LIXE9cUFAE072V6tHvMxixjFdDmYVN8nHceKLwClOYDGammR03+QUBEU8uiMOcs8CtjZaLaeFUBaDAYcOLECfTs2RPNmjVDXFwcduzYgU6dOgEAysrK8MMPP+C1115zuB6tVgutVlsfXSYvkhgpLbYGnI+ua8bXckfXVbH1N39kolFEMKJCNIgKNV2mJzpUi8iq58Ea+KttHHkhSeQYuxQ3bjmKj6uUFZli49xzQO55U1GYe970XH8JqCixnolcnTACX//dFDdHtQTCE4DwRje+JgBh8YC/gj4zIhkpugCcNm0ahgwZgiZNmuDq1atYsGABCgoKMGLECKhUKkyePBkLFy7Ebbfdhttuuw0LFy5EcHAw/vrXv7q76+SlpMTWgPzRdc342l50nVtkQEm59YzrH0/XHltXXbcxOuRmYWgqGKs/1yIq1HRtR7UT8bSvRNEcu26wFx9X0YQAsXeYHjVVlgP5F6sVheeByweB9F+slz2z0/SwJaShdWFo/r6RqUjUBNvfB8bNpFCKLgAvXbqEJ598EtnZ2WjYsCG6d++O/fv3IykpCQAwffp0lJSUYNy4ceaLqW7fvp3XACSXqu/YGpAeXdu657SfCpjS/3ZUVBqRW1SGnMIy5BQZkFNoukRPbnEZhADyi8uRX1yOc9eKat2OSgXzBb+jqhWHpiOLGkSGaPH7pXx8sPecT0TRHLuqcTY+rkkdAES1MD2q2Lp4tcoP6DUDKC8CCjJuPC6bvlYaTDORi66Zzku0JyjCdoF45Tjwy7uMm0mRFH0nkPrizXcCIdp44GKdYutKo0B+cRlyqhWHuUWm8xZzbxSKptcM5ouES/V458a4M7EBbosJRcuYUESFyh/Xeevvt7fuV60Ob3AcK1cRAijOvVkMFlwGrmdaFoj6y6bCsS7iOwINmpguhRMWa5qoUv37kGjAT13raohq49O3gqsvPjuQks+4lXtO16ai0oi84nJToVhYhuyiMuQWGkxF4o1CMS27CKeuFNa6rsgQDVreKAZviwnFbTFhaBkTithwrcUM2bpEyd76++2t++WUqvsZ24qV60IIwFBgWRRWfZ91FMg4XPd1qvxMsXNoLBAWB4TGmApDi+9jTa8HOPjZZfTs81gA1gOfHkiJ6oGtKFqlAp68uwky80tw+mohLuWV2H1/mNYfLWND0bJhKIrKKvDfo1kQTkbJ3vr77a375THsxc2Dl5nOTyy8YprNfP3Kje+vmKJmR/dWrkmru1kMmgvGWNO1En/70FSgqvyAB5eb7s5CPoUFYD3gQErkerVF0cVlFTh3rQinr17HmauFOH2lEGeuFuJCbjEqjfaHMbVKhZ9e7mP3SKC3/n576355FGfj5irGStMFss2FYZapMDR/fxW4fqOtorRufQlLMB0JtCgWaxxRDIkB1Io+9Z+qudXfcf4kEJFHqG0GdbDGH+0a6dCukc6i3VBRibTsYvyaloPPD6bjj0uWd8CoFAJp2cVePauY3KS2Wcw1+alNxVhYLBDvYLmq6NmqMMwCMo8A5/dYv+d6hunhkMp0DmLNI4oW5yvGmNo1IbZXwejZa7AAJCKPIWUGtdZfjX1nszFv63HYyjPUKhWaRju4jAfRrZA6i9kRlQoI1JkeDVtZvmYveh72sanN4uhitSOKhVdNRyqrZj1fOeq4D5qwm8VgVbGovwSc2ALgRvQ8cDFw9/Om/pLisAAkIo9XUlaJS3nFSM8rRnpuCdJzb35/IbcIRYZKm++ripJ59I+8hr0LaLd+wPH7jEagOMdx7Fz1tbwYKLsO5F43He2zRRiB/04Hdsw29Sks3vQIjzfF0dW/hsaaLttDHoUFIBHVK1szdMsrjcjML71R1FUr9G58zS401Hk7swa3wQMd4ln8kfepa/QMmO6ZHNrQ9EB7x8sarlebvHLjiOKlg8CxTdbLVpQAOWdMD7tUNy64baM4DIu/eVeWQJ310URGzi7DApCIXEYIgeuGCuQVmS4i/Z/fLmP9vguoSmqbR4fAUGFEpr4EDuZxADDN9G0cGYzEiCAkVvsaGKDG8DW/WLxfrVKx+CPv5orouYo2zPSIbnmzTX8ZOP5ljehZDYzcamoryDSdg2jx9cbDWAEUXTU9HF1wOyDYFDVXFYfFucDZXbCInLuNdc0++yAWgETklJrFXF5xGfKKypFXXGa+4HPejTuF5BeXIbeoHPnFZahwUNmdy755kV2Nvx8aRwQhMSIYiZFVX4PNz3VBARbX+qvO1j2YWfwRyche9JzUw/H7jEagONt0/cTrWbaLxIIMoDTfFD3nnjM9aqqKnHfNByKaARFJQIMkIKKp6WuDJqaHo9v2kQUWgEReytHFkF1RzDkSrFEjOECN7KIyq9fmDGmLQe3j0TBUCz8n7hdsi9R7MBNRHUiOnmNMD0fKim8eMSzIBC78BBxaZ72c4TqQ9YfpYUtIzM3isEGTaoViEqBLtDwX0cfjZRaARApmNAoUllVAX1wOfUk5CkrLUVBSjh3Hr2Dz4cvmqLVdQjiCtf6yFXMRwRpEhASY7+8bEay58X0AGlRvu7FMYIDa5sWe1SoVBrSLQ2x43e5hbIs77sFM5HNcFT1rgi3v7ZzUw3SdxZqR89ObTecd5l0A8i8C+RdufH/BdOmcqqj50gHrbaj8TPFyRBJQWQFc+hXmeHnwm0CXkfLvlwdjAUjkhLrcWqyuyiuNKCgpR0FpBfQlNwq5G1+rvi8orf785nLXS8trPXcOAI5mFNhsl1rMSRGvC2JUS0TOsRc5t+hte3khgJI866Iw/+LNYrGiBCi4ZHpYvNcIbJ0E/PKu6V7OsW2BmDuAmDamCSr2LnOj8COIvBOIE3hFfd+28cBFzNh8BEYHtxYrLa+0W7zpqxVsVYVcQbXlispsX8KkLrT+fggPCoAuKAB+KuB/Nu6rO/H+lujWPEqWYu5WuPK+w1J46++3t+4X+Rg57+tceNVUFP5vO7B3qXPvC2wAxLS9URS2uVkYnvga2DLJVDyq/EzFqqO7wLgAbwVXDziQej8hBAoNNwu1qiNtF3OKsOi/J1Hzl6RdQjhKyiuhL6lAQWk5yirqcH9PO0K1/tAFBSAs0PS16hFu8X2N1wJNr1cv5OxFrY5uh+bLvPX321v3i+iW2buY9pC3TJe8uXoMuHrCdM9l4eR/0FVqYPKRej0SyFvBkSK5IlKtqDSaY1SLo3ClNeJUO0fk6nJKnK1I1U8Fi2KteoFWs3gLD7RcLizQH/5qP1k+B0atREQO2IuXax7BqzAA2f8DrhwHrt54XDluHSEDpvVsfBpo9xjQ4n7TUcLq0bEHxsU8AugE/k9aXo4i1dLyStvFW7HlOXLVi7yq8+cKDRW33DeNv9+NAs1UrGn9/ZByLtdiGZUKWDi0PZpGh1gUdqFaf7uXKXEHT4taPZW3/n57634RyUZqvHzlBLCqh+URxJpC44DmvU3FYNE1YMcs2eNiRsD1gANp3VVdZiS3sAw5Ny41kltkQFpOEVbtOWcVqUaGaFBoqJAlSg3RqM3Rqa0jcrogf+iCrY/C1YxSq2w8cNHqaFrNcwBJubz199tb94vIIxzeYHkEsdd0QBMKnNsNpP1smnBij0xxMSNgctqtxK5CCBSUVCCnyIDcIlNRl1NoKupuFnhVbaZHWaXzxVxutevDVUWp1gWa/804NdDOOXIyRqlVeI05IiKyYO+aiD0mAOWlQPovpmLw+NfW91MWlUDGb6b3uDEaZgHoI2rGrq8+0h4D74i7UcjdLOrMhdyNI3Y5N47g5RVJu25csEaNyBANokJMlxMJDFDju6NZFkcA/VTA2pFd0SImFOFBAQjV+Eu+ILCr8BpzRERkwd41EQMCgea9TI+uz1lPOAGATWOApGTg3B63zSRmBOwEpUUpRqPAtUID0nOLcSmvBMcz9Hh/73lZ1h2q9UfkjWKuqqiLCtWav48MrdYeokWQhpEqeTal/X47y1v3i0hxLOJiP9O1BfU2JpLUMRpmBOxlnIlphRDIKSozF3jpeTe+5hbjcl4JLuWXOH0uXVigv+1CLkSDqFANIkNMbVGhGtmuG8dIlYiIfEbNuDg8AfjpTeD7uZbLiUrg+FdA26H1EgezAPQgNWPa8X1aom18uEWBdymvBJfySlBS7vjaRH4qU2yZGBmEyBAN/nvEOnbdM603mkSFuHan7GCkSkREPqNmXNxhGLBrvnU0vO3/gO2v1EsczALQAxQZKrDtWBZe3nTEXKQZBbBi1xm771GpgNiwQCRGBqFxRDASI0xfG0cGITEiGHG6QARUmwxhK3Z1V/FHRETk02pei7A6YTS1t+jr0iOBLADrQc1Yt7S8Eocv5CHlXA72nc3B7+n5didYtIwJQZt4nbnAqyr4EhoEQuvvfBzL2JWIiMiDVEXDx74Ets+0fE1UmiJjFoDKVT3WVQFoFh1i8xy9OF0gsvSlFm1qlQofjukmW7HG2JWIiMiD6BoBdzxiin2rn6ilUpnOF3QhFoAulJFfbBHrCgDnsosAALHhWiQ3j0KPFtFIbhGFxMhgmzEtCzYiIiIfUw/XZ2EB6AKZ+hIcyyjA27tO2/w3fOMvd+LRTo2sbhvGmJaIiMjH5J6FdcUnGAErzcYDF/Hy5iOwd3VFtUqFHi2i7N4zljEtERGRD4lsYbo+YM0ZwRmHgWY9XbZZee+Z5eMy9SWYUaP4U8F0yRUAjHWJiIjIkq4R8Ke51u0755puFeciPAIoo/PZRag5mVcAWPFEJ0SFahnrEhERkbWEjtZtLp4JzAJQRkcu6a3a1CoVOjeNYOFHREREtkW2sG5z8UxgRsAyydSX4LXvTlq1Tx90O4s/IiIiqhsXzwRmASgTW/EvAHRo1KDe+0JEREQKknvWRuONmcAuwgJQJs2iQ1BzXq9KBTSNDnZLf4iIiEghGAF7mXq4kCMRERF5IUbAynA+u8jWZRyRll3sju4QERGRUjACVi5GwERERCQJI2AvwwiYiIiIpGAErAyMgImIiEgSRsDKxQiYiIiIJGEE7GUYARMREZEUjIDtW7RoEbp27YqwsDDExMRg6NChOHXqlMUyI0eOhEqlsnh0795d9r4wAiYiZ3nS2EVEHoARcN388MMPGD9+PPbv348dO3agoqIC/fv3R1FRkcVyAwcORGZmpvnx7bffyt4XRsBE5CxPGruIyAO4IQL2d9ma68F3331n8Xzt2rWIiYnBoUOHcN9995nbtVot4uLi6rt7jICJyCaPH7uIyP0YATtPr9cDACIjIy3a9+zZg5iYGLRq1QrPPfccrl696nA9BoMBBQUFFo/aMAImIqnkGLukjFtE5CEYAUsnhMCUKVNw7733ol27dub2QYMG4eOPP8auXbvwxhtv4MCBA7j//vthMBjsrmvRokXQ6XTmR2JiYq3bD9GobbYHa7zmIyYiF5Br7JIybhGRhwgIsdPuutPIVEIIrwgqx48fj2+++QY//fQTGjdubHe5zMxMJCUl4bPPPsOjjz5qcxmDwWAxyBYUFCAxMRF6vR7h4eE237PvbDb+uvoXq/ZPn+uO5BZRddwbIqovBQUF0Ol0Dn+/XUmusUvKuEVEHuL8j8D6IdbtI7YCzXrafMutjl2KPgewyt///nd8/fXX+PHHHx0OoAAQHx+PpKQknD592u4yWq0WWq22Tn2omgRSvZrmJBAickTOsUvKuEVEHoKTQOpGCIG///3v+PLLL7Fnzx40a9as1vfk5OQgPT0d8fHx9dBB12+CiJTH48cuInI/TgKxb/z48fjoo4/wySefICwsDFlZWcjKykJJSQkAoLCwENOmTUNKSgrS0tKwZ88eDBkyBNHR0XjkkUdk7QsngRCRszxp7CIiD+CGSSCKPgL47rvvAgB69+5t0b527VqMHDkSarUaR44cwYYNG5Cfn4/4+Hj06dMHGzduRFhYmKx9YQRMRM7ypLGLiDwAI+C6qW3+SlBQELZt21ZPvbGBETAR2eDxYxcRuZ+La4hbLgD1ej2+/PJL7N27F2lpaSguLkbDhg3RqVMnDBgwAD169JCjnx7PUQQcrwtyR5eIiIhICRxFwLpGLtmk5HMAMzMz8dxzzyE+Ph7z5s1DUVEROnbsiL59+6Jx48bYvXs3+vXrh7Zt22Ljxo1y9tkj8VZwREREJImSIuA777wTzzzzDH799VeLi5dWV1JSgq+++grLli1Deno6pk2bJrmjisQImIiIiKTw1Aj42LFjaNiwocNlgoKC8OSTT+LJJ5/EtWvXpG5KERgBExERkSRKioBrK/5udXmlYQRMREREkigpAq7u66+/ttmuUqkQGBiIli1bOnWhU6/DCJiIiIik8NQIuLqhQ4dCpVJZXdqgqk2lUuHee+/FV199hYiICDk26XEYARMREZEkSoqAq9uxYwe6du2KHTt2QK/XQ6/XY8eOHbj77ruxdetW/Pjjj8jJyfHqSSCMgImIiEgSpUbAkyZNwvvvv29xzb++ffsiMDAQzz//PI4dO4bly5dj9OjRcmxOORgBExERkRRKuBfw2bNnER4ebtUeHh6Oc+dM97G77bbbkJ2dLcfmPBLvBUxERESSuOFewLIUgJ07d8Y//vEPi0u9XLt2DdOnT0fXrl0BAKdPn0bjxo3l2JxHYgRMREREkig1Al6zZg0efvhhNG7cGImJiVCpVLh48SKaN2+O//znPwCAwsJCzJo1S47NKQcjYCIiIpJCCbOAb7/9dpw4cQLbtm3D//73Pwgh0Lp1a/Tr1w9+fqaDjEOHDpVjUx6Ls4CJiIhIEjfMApalAARMl3wZOHAgevfuDa1WC5WqZiDq3aoi4OpFICNgIiIiqpUbImBZzgE0Go2YP38+GjVqhNDQUJw/fx4AMGvWLKxZs0aOTSgTI2AiIiKSQgmzgBcsWIB169ZhyZIl0Gg05vb27dvjgw8+kGMTHo+zgImIiEgSpc4C3rBhA95//3089dRTUKvV5vYOHTrg5MmTcmzC43EWMBEREUmi1Aj48uXLaNmypVW70WhEeXm5HJtQJkbAREREJIUSIuA77rgDe/futWr//PPP0alTJzk24fEYARMREZEkboiAZZkFPHv2bAwfPhyXL1+G0WjE5s2bcerUKWzYsAFbt26VYxMej7OAiYiISBKlRsBDhgzBxo0b8e2330KlUuGf//wnTpw4gS1btqBfv35ybEKZGAETERGRFEq4EDQADBgwAAMGDJBrdYrDC0ETERGRJG64ELQsRwCJs4CJiIhIIiXdCzgiIsLpu33k5uZK3YyyMQImIiIiKTw1Al6+fLn5+5ycHCxYsAADBgxAcnIyACAlJQXbtm3DrFmzbrmTSsAImIiIiCRR0r2AR4wYYf7+sccew7x58zBhwgRz28SJE/H2229j586dePHFF2+tlwoQolHbbA/WMGUnIiIiBwJC7LS77jQyWaqTbdu2YeDAgVbtAwYMwM6dO+XYhMcrKqu02V5cZqznnhAREZGilBfZaXfdtYRlKQCjoqLw5ZdfWrV/9dVXiIqKkmMTHo+TQIiIiEgSJU0CqW7u3LkYM2YM9uzZYz4HcP/+/fjuu+/wwQcfyLEJZeIkECIiIpJCCbeCGzlyJPbt24cGDRpg8+bN2LRpE3Q6HX7++WeMHDlSjk14PN4KjoiIiCRR6q3gAKBbt274+OOP5Vqd4vBWcERERCSJkm4FV1Rk54RFmZb3CoyAiYiISApPjYBbtmyJhQsXIiMjw+4yQgjs2LEDgwYNwr/+9S+pm1IERsBEREQkiZIi4D179uCVV17B3Llz0bFjR3Tp0gUJCQkIDAxEXl4ejh8/jpSUFAQEBGDGjBl4/vnn5ey3x2EETERERJIoaRbw7bffjs8//xyXLl3C559/jh9//BH79u1DSUkJoqOj0alTJ6xevRoPPPAA/Px89GLIjICJiIhICk+9FVyVxo0b48UXX/SJu304wlvBERERkSRuuBWcjx6akx8vBE1ERESSKGkWMDmBETARERFJ4amzgMkSZwETERGRJG6YBcwCUCaMgImIiEgSRsBehhEwERERSaGUCHjv3r14+umnkZycjMuXLwMAPvzwQ/z0009ybeKWvPPOO2jWrBkCAwPRuXNn7N27V9b1MwImIldw9dhFRB5AqRHwpk2bMGDAAAQFBSE1NRUGgwEAcP36dSxcuFCOTdySjRs3YvLkyZg5cyZSU1PRs2dPDBo0CBcvXpRtG4yAiUhu9TF2EZEHUGoEvGDBAqxatQqrV69GQECAub1Hjx44fPiwHJu4JcuWLcOYMWPw7LPPok2bNli+fDkSExPx7rvvunbDjICJ6Ba4bewiIvdTQgR86tQp3HfffVbt4eHhyM/Pl2MTkpWVleHQoUPo37+/RXv//v2xb98+m+8xGAwoKCiweNSGETARyamuY5eUcYuIPIRSI+D4+HicOXPGqv2nn35C8+auO3zpjOzsbFRWViI2NtaiPTY2FllZWTbfs2jRIuh0OvMjMTGx1u0wAiYiOdV17JIybhGRh1BqBDx27FhMmjQJv/zyC1QqFTIyMvDxxx9j2rRpGDdunBybuGUqlWV5JoSwaqsyY8YM6PV68yM9PV3aRhkBE9Etcnbskm3cIiLP4On3AgaA6dOnQ6/Xo0+fPigtLcV9990HrVaLadOmYcKECXJsQrLo6Gio1Wqr/zFfvXrV6n/WVbRaLbRabZ22w3sBE5Gc6jp2SRm3iMhDKPlewK+++iqys7Px66+/Yv/+/bh27Rrmz58v1+ol02g06Ny5M3bs2GHRvmPHDvTo0UO27TACJiI51dfYRUQewA0RsCxHAKsEBwejS5cucq5SFlOmTMHw4cPRpUsXJCcn4/3338fFixfxwgsvuHbDjICJ6Ba4bewiIvfz1Aj40UcfdXrZzZs3S92MLIYNG4acnBzMmzcPmZmZaNeuHb799lskJSXJtg1GwEQkt/oYu4jIA7ghApZcAOp0OvP3Qgh8+eWX0Ol05iOAhw4dQn5+fp0KRVcaN26cSyekVEXA1YtARsBEdKtcPXYRkQdQUgS8du1a8/cvvfQS/vKXv2DVqlVQq9UAgMrKSowbNw7h4eG33kulYgRMREREUijhQtD//ve/MW3aNHPxBwBqtRpTpkzBv//9bzk24fF4IWgiIiKSRKkXgq6oqMCJEyes2k+cOAGj0SjHJjxeiEZtsz1YI9tEayIiIvJGASF22l13Gpkss4BHjRqF0aNH48yZM+jevTsAYP/+/Vi8eDFGjRolxyY8XlFZpc324jLfKICJiIhIovIiO+2uSxFlKQBff/11xMXF4c0330RmZiYA0+3hpk+fjqlTp8qxCY/HSSBEREQkiZImgVTn5+eH6dOnY/r06eYbkPv05I8qnARCREREUihhEkh14eHhPln8cRIIERERSeKGSSCyHAFs1qyZzZuTVzl3znU74CkYARMREZEkSo2AJ0+ebPG8vLwcqamp+O677/CPf/xDjk0oEyNgIiIiksJTbwVX3aRJk2y2r1y5EgcPHpRjEx6Pt4IjIiIiSdxwKziXXqRu0KBB2LRpkys34TGqIuDqGAETERFRrdwQAbu0APziiy8QGRnpyk14NkbAREREJIUSIuBOnTpZTAIRQiArKwvXrl3DO++8I8cmPB4jYCIiIpLEDRGwLAXgww8/bFEA+vn5oWHDhujduzdat24txyY8HmcBExERkSRKnQU8Z84cOVbjfRgBExERkRRKuBC0Wq3G1atXrdpzcnKgVqvl2ITH44WgiYiISBI3XAhalgJQCNtlqsFggEajkWMTHo+zgImIiEgSpUXA//rXvwAAKpUKH3zwAUJDQ82vVVZW4scff/SZcwBtYgRMREREUnjyLOA333wTgOkI4KpVqyziXo1Gg6ZNm2LVqlW31kOF4CxgIiIikkRps4DPnz8PAOjTpw82b96MiIgIWTqlRJwFTERERJIoLQKusnv3bjlW430YARMREZEUnhoBT5kyBfPnz0dISAimTJnicNlly5ZJ3YxiMAImIiIiSZQUAaempqK8vBwAcPjwYYsLQfsiRsBEREQkiZIi4Oqx7549e+Toi/dhBExERERSKOFC0KNHj8b169et2ouKijB69Gg5NuHxeCFoIiIikkSpF4Jev349SkpKrNpLSkqwYcMGOTbh8XghaCIiIpJESREwABQUFEAIASEErl+/jsDAQPNrlZWV+PbbbxETE3PLnVQsRsBEREQkhafOAgaABg0aQKVSQaVSoVWrVlavq1QqzJ0791Y2oRicBUxERESSKGkWMGCaCCKEwP33349NmzYhMjLS/JpGo0FSUhISEhJuuZNKwFnAREREJInSIuBevXoBMN0RJDExEX5+spxS6D0YARMREZEUnhwBV0lKSgIAFBcX4+LFiygrK7N4vUOHDnJsxqMxAiYiIiJJlBYBV7l27RpGjRqF//73vzZfr6yslGMzHo0RMBEREUnihghYlsx28uTJyMvLw/79+xEUFITvvvsO69evx2233Yavv/5ajk0oEyNgIiIikkIJEfCuXbvwn//8B127doWfnx+SkpLQr18/hIeHY9GiRRg8eLAcm/FojICJiIhIEjdEwLIcASwqKjJf7y8yMhLXrl0DALRv3x6HDx+WYxMeL0SjttkerOHEGCIiInIgIMROu+tOI5OlOrn99ttx6tQpAEDHjh3x3nvv4fLly1i1ahXi4+Pl2ITHKyqzfZ5jcZmxnntCREREilJeZKfddbeTlSUCnjx5MjIzMwEAs2fPxoABA/Dxxx9Do9Fg3bp1cmzC43ESCBEREUkS2QKoWUV48nUAqzz11FPm7zt16oS0tDScPHkSTZo0QXR0tBybUCZOAiEiIiKnCIdP5eaSE9SCg4Nx1113+VTx52gSCBEREZFdjiaBuIjkI4BTpkxxetlly5ZJ3YxiMAImIiIiSZQUAaempjq1nEqlkroJ5WMETERERE6p3whYcgG4e/duOftRZ2lpaZg/fz527dqFrKwsJCQk4Omnn8bMmTOh0WjMy9kqQN9991288MILsvaH1wEkImd42thFRB5AqbeCc4eTJ0/CaDTivffeQ8uWLXH06FE899xzKCoqwuuvv26x7Nq1azFw4EDzc51OJ3t/GAETkTM8bewiIg+gpAjY3QYOHGgxMDZv3hynTp3Cu+++azWINmjQAHFxcfXdRUbARGRFEWMXEbmBF8wCdhe9Xo/IyEir9gkTJiA6Ohpdu3bFqlWrYDQ6vjizwWBAQUGBxaM2nAVMRFLJMXZJGbeIyEMoaRawpzl79ixWrFiBN954w6J9/vz56Nu3L4KCgvD9999j6tSpyM7OxiuvvGJ3XYsWLcLcuXPrtH1GwEQkhVxjl5Rxi4g8hBsiYJUQwqOCyjlz5tQ6iB04cABdunQxP8/IyECvXr3Qq1cvfPDBBw7f+8Ybb2DevHnQ6/V2lzEYDDAYDObnBQUFSExMhF6vR3h4uM33ZOpL0GPRLssCEMC+GfdzEgiRBysoKIBOp3P4++0Md49dUsYtIvIQ+svAm21rNKqAF4/ZnQRyq2OXxx0BnDBhAp544gmHyzRt2tT8fUZGBvr06YPk5GS8//77ta6/e/fuKCgowJUrVxAbG2tzGa1WC61WW6d+cxYwkW9z99glZdwiIg/BWcBAdHS003cQuXz5Mvr06YPOnTtj7dq18POr/ZTG1NRUBAYGokGDBrfYU0uMgIl8m1LHLiLyAJwF7LyMjAz07t0bTZo0weuvv45r166ZX6uaNbdlyxZkZWUhOTkZQUFB2L17N2bOnInnn3++fv6n7FHhOhF5AkWMXUTkBgq5ELS7bd++HWfOnMGZM2fQuHFji9eqTmsMCAjAO++8gylTpsBoNKJ58+aYN28exo8fL3t/GAETkTM8bewiIg/ghgjY4yaBeCJnTrS0OQlEBex7mZNAiDyZXJNAPI237heRV9JfBt68A1YR8GTXTQLxqusAehyW1kREROQUXghakXghaCIiIpLEDReCZgEok6pZwNVxFjARERHVyjwLuBoXzwJmAehKjICJiIjIKYyAFYkRMBEREUnCCFi5GAETERGRJIyAvQwjYCIiInIKI2BFYgRMREREkjACVi5GwERERCQJI2AvwwiYiIiInMIIWJEYARMREZEkjICVixEwERERScII2MswAiYiIiKnMAJWJEbAREREJAkjYOUK0ahttgdr+BETERGRAwEhdtpddxoZqxOZFJVV2mwvLjPWc0+IiIhIUcqL7LS7LkVkASgTTgIhIiIiSTgJxMtwEggRERE5hZNAFImTQIiIiEgSTgJRLkbAREREJAkjYC/DCJiIiIicwghYkRgBExERkSSMgJWLETARERFJwgjYyzACJiIiIqcwAlYkRsBEREQkCSNg5WIETERERJIwAvYyjICJiIjIKYyAFYkRMBEREUnCCFi5GAETERGRJIyAvQwjYCIiInIKI2BFYgRMREREkjACVi5GwERERCQJI2AvwwiYiIiInMIIWJEYARMREZEkjICVixEwERERScII2MswAiYiIiKnMAJWJEbAREREJAkjYOViBExERESSMAL2MoyAiYiIyCmMgBWJETARERFJwghYuRgBExERkSSMgOumadOmUKlUFo+XX37ZYpmLFy9iyJAhCAkJQXR0NCZOnIiysrL66SAjYCKywePHLiJyg/qNgP1du3rXmzdvHp577jnz89DQUPP3lZWVGDx4MBo2bIiffvoJOTk5GDFiBIQQWLFihaz9cBQBx+uCZN0WESmfp4xdROQBHEXAukYu2aTiC8CwsDDExcXZfG379u04fvw40tPTkZCQAAB44403MHLkSLz66qsIDw+XrR8hGrXN9mCNog+yEpGLeMrYRUQeICDETrvrTiNTfHXy2muvISoqCh07dsSrr75qEZGkpKSgXbt25gEUAAYMGACDwYBDhw7ZXafBYEBBQYHFozZFZZU224vLjHXYGyLyFXKPXVLGLSLyEOVFdtpdN5FU0UcAJ02ahLvuugsRERH49ddfMWPGDJw/fx4ffPABACArKwuxsbEW74mIiIBGo0FWVpbd9S5atAhz586tU1+qJoFUj4E5CYSIbHHF2CVl3CIiD2GeBFKtivC1SSBz5syxOjm65uPgwYMAgBdffBG9evVChw4d8Oyzz2LVqlVYs2YNcnJyzOtTqWrOzQWEEDbbq8yYMQN6vd78SE9Pl7YznARC5DPcPXbJNm4RkZv4+CSQCRMm4IknnnC4TNOmTW22d+/eHQBw5swZREVFIS4uDr/88ovFMnl5eSgvL7f633V1Wq0WWq22Tv3mJBAi3+busUvKuEVEHoKTQIDo6GhER0dLem9qaioAID4+HgCQnJyMV199FZmZmea27du3Q6vVonPnzvJ0+AZGwES+TaljFxF5ADdEwB5XADorJSUF+/fvR58+faDT6XDgwAG8+OKLeOihh9CkSRMAQP/+/dG2bVsMHz4cS5cuRW5uLqZNm4bnnnuufmbRMQImohoUMXYRkRvwVnBO0Wq12LhxI3r37o22bdvin//8J5577jl8+umn5mXUajW++eYbBAYG4p577sFf/vIXDB06FK+//rrs/eGt4IjIGZ42dhGRB3DDreAUewTwrrvuwv79+2tdrkmTJti6davL+8MImIic4WljFxF5AM4C9jKMgImIiMgpjIAViREwERERSeKGCJgFoEyqIuDqGAETERFRrcwRcDWMgBWMETARERE5hRGwIjECJiIiIkkYASsXI2AiIiKShBGwl2EETERERE5hBKxIjICJiIhIEkbAysUImIiIiCRhBOxlGAETERGRUxgBKxIjYCIiIpKEEbByMQImIiIiSRgBexlGwEREROQURsCKxAiYiIiIJGEErFyMgImIiEgSRsBehhEwEREROYURsCIxAiYiIiJJGAErFyNgIiIikoQRsJdhBExEREROYQSsSIyAiYiISBJGwMrFCJiIiIgkYQTsZRgBExERkVMYASsSI2AiIiKShBGwcoVo1DbbgzX8iImIiMiBgBA77a47jYzViUyKyipttheXGeu5J0RERKQo5UV22l2XIrIAlAkngRAREZEknATiZTgJhIiIiJzCSSCKxEkgREREJAkngSgXI2AiIiKShBGwl2EETERERE5hBKxIjICJiIhIEkbAysUImIiIiCRhBOxlGAETERGRUxgBKxIjYCIiIpKEEbByMQImIiIiSRgBexlGwEREROQURsCKxAiYiIiIJGEErFyMgImIiEgSRsBehhEwEREROYURsCIxAiYiIiJJGAE7b8+ePVCpVDYfBw4cMC9n6/VVq1bJ3h9GwETkDE8bu4jIA7ghAvZ32ZpdrEePHsjMzLRomzVrFnbu3IkuXbpYtK9duxYDBw40P9fpdPXSR0bARFSTIsYuInKD+o2AFVsAajQaxMXFmZ+Xl5fj66+/xoQJE6BSWVbRDRo0sFjWFRxFwPG6IJdum4iUw9PGLiLyAI4iYF0jl2xSsRFwTV9//TWys7MxcuRIq9cmTJiA6OhodO3aFatWrYLRaJR9+4yAiUgKd49dROQBGAFLt2bNGgwYMACJiYkW7fPnz0ffvn0RFBSE77//HlOnTkV2djZeeeUVu+syGAwwGAzm5wUFBdI6xQiYiGoh19gl27hFRG7i47OA58yZY/cE6arHwYMHLd5z6dIlbNu2DWPGjLFa3yuvvILk5GR07NgRU6dOxbx587B06VKHfVi0aBF0Op35UXNgtoWzgIl8m7vHLinjFhF5CDfMAlYJITzqOFV2djays7MdLtO0aVMEBgaan8+fPx8rVqzA5cuXERAQ4PC9P//8M+69915kZWUhNjbW5jK2/iedmJgIvV6P8PBwm+/J1Jegx6JdFkWgSgXse/l+ngNI5MEKCgqg0+kc/n47w91jl5Rxi4g8hP4y8OYdsDjsp1IBk4/ZPQfwVscuj4uAo6OjER0d7fTyQgisXbsWzzzzTK0DKACkpqYiMDAQDRo0sLuMVquFVqt1ug/2O3frqyAiZXD32CXbuEVEbsJZwHWya9cunD9/3maEsmXLFmRlZSE5ORlBQUHYvXs3Zs6cieeff172gZKzgImoLjxl7CIiD+CGWcCKLwDXrFmDHj16oE2bNlavBQQE4J133sGUKVNgNBrRvHlzzJs3D+PHj5e9H1WzgGtGwJwFTES2eMrYRUQewDwLuEYEzFnA9n3yySd2Xxs4cKDFRVTrHSNgIrLDo8cuInIDH58FrFScBUxERESS8F7AysULQRMREZEkbrgQNAtAV2IETERERE5hBKxIjICJiIhIEkbAyhWiUdtsD9bwIyYiIiIHAkLstLvuNDJWJzIpKqu02V5cxpu3ExERkQPlRXbaXZcisgCUCSeBEBERkSScBOJlOAmEiIiInMJJIIrESSBEREQkCSeBKBcjYCIiIpKEEbCXYQRMRERETmEErEiMgImIiEgSRsDKxQiYiIiIJGEE7GUYARMREZFTGAErEiNgIiIikoQRsHIxAiYiIiJJGAF7GUbARERE5BRGwIrECJiIiIgkYQSsXIyAiYiISBJGwF6GETARERE5hRGwIjECJiIiIkkYASsXI2AiIiKShBGwl2EETERERE5hBKxIjICJiIhIEkbAytUsOgR+NY7eqlUqRsBERETkWGQLQFWjJFOpGQErQbwuCIsebQ+1ylQFqlUqLHy0HeJ1QW7uGREREXk0XSNgyFumog8wfR2y3NTuIv4uW7MPGta1Ce5r1RBp2cVoGh3M4o+IiIicc9czQIu+ptg3srlLiz+ABaDs4nVBLPyIiIio7nSNXF74VWEETERERORjWAASERER+RgWgEREREQ+hgUgERERkY9hAUhERETkY1gAEhEREfkYFoBEREREPoYFIBEREZGPYQFIRERE5GNYABIRERH5GBaARERERD6GBSARERGRj/F3dweUQAgBACgoKHBzT4hIblW/11W/596C4xaRd7vVsYsFoBOuX78OAEhMTHRzT4jIVa5fvw6dTufubsiG4xaRb5A6dqmEt/231wWMRiMyMjIQFhYGlUrlcNmCggIkJiYiPT0d4eHh9dRD+XnLfgDesy/cD9cQQuD69etISEiAn5/3nBVjb9zytM/fXfg53MTPwkRpn8Otjl08AugEPz8/NG7cuE7vCQ8PV8QPUG28ZT8A79kX7of8vOnIX5Xaxi1P+vzdiZ/DTfwsTJT0OdzK2OU9/90lIiIiIqewACQiIiLyMSwAZabVajF79mxotVp3d+WWeMt+AN6zL9wPkgM/fxN+DjfxszDxtc+Bk0CIiIiIfAyPABIRERH5GBaARERERD6GBSARERGRj2EBSERERORjWADW4p133kGzZs0QGBiIzp07Y+/evQ6X/+GHH9C5c2cEBgaiefPmWLVqldUymzZtQtu2baHVatG2bVt8+eWXruq+Bbn35dixY3jsscfQtGlTqFQqLF++3IW9v0nu/Vi9ejV69uyJiIgIRERE4E9/+hN+/fVXV+4CAPn3Y/PmzejSpQsaNGiAkJAQdOzYER9++KErd8HMFb8nVT777DOoVCoMHTpU5l57vry8PAwfPhw6nQ46nQ7Dhw9Hfn6+w/cIITBnzhwkJCQgKCgIvXv3xrFjxyyWMRgM+Pvf/47o6GiEhITgoYcewqVLl+q87UmTJqFz587QarXo2LGjVV/S0tKgUqmsHt99951PfQ4AcOTIEfTq1QtBQUFo1KgR5s2bV+d7uHr653Dx4kUMGTIEISEhiI6OxsSJE1FWVmZ+/VZ+Htz1t7i27cr1+bqFILs+++wzERAQIFavXi2OHz8uJk2aJEJCQsSFCxdsLn/u3DkRHBwsJk2aJI4fPy5Wr14tAgICxBdffGFeZt++fUKtVouFCxeKEydOiIULFwp/f3+xf/9+xe3Lr7/+KqZNmyY+/fRTERcXJ958802X7oOr9uOvf/2rWLlypUhNTRUnTpwQo0aNEjqdTly6dElR+7F7926xefNmcfz4cXHmzBmxfPlyoVarxXfffeey/XDVvlRJS0sTjRo1Ej179hQPP/ywS/fDEw0cOFC0a9dO7Nu3T+zbt0+0a9dOPPjggw7fs3jxYhEWFiY2bdokjhw5IoYNGybi4+NFQUGBeZkXXnhBNGrUSOzYsUMcPnxY9OnTR9x5552ioqKiTtv++9//Lt5++20xfPhwceedd1r15fz58wKA2Llzp8jMzDQ/DAaDT30Oer1exMbGiieeeEIcOXJEbNq0SYSFhYnXX3/daz6HiooK0a5dO9GnTx9x+PBhsWPHDpGQkCAmTJhgXkbqz4O7/hY7s125Pl93YAHowN133y1eeOEFi7bWrVuLl19+2eby06dPF61bt7ZoGzt2rOjevbv5+V/+8hcxcOBAi2UGDBggnnjiCZl6bZsr9qW6pKSkeikAXb0fQpgGsrCwMLF+/fpb77Ad9bEfQgjRqVMn8corr9xaZ2vhqn2pqKgQ99xzj/jggw/EiBEjfK4APH78uABg8QcpJSVFABAnT560+R6j0Sji4uLE4sWLzW2lpaVCp9OJVatWCSGEyM/PFwEBAeKzzz4zL3P58mXh5+dn/s9CXbc9e/ZshwVgampqnfa9Om/4HN555x2h0+lEaWmpuW3RokUiISFBGI1Gr/gcvv32W+Hn5ycuX75sXubTTz8VWq1W6PV6IYT0nwd3/S2ubbtyfb7uwgjYjrKyMhw6dAj9+/e3aO/fvz/27dtn8z0pKSlWyw8YMAAHDx5EeXm5w2XsrVMOrtqX+lZf+1FcXIzy8nJERkbK0/Ea6mM/hBD4/vvvcerUKdx3333ydb4GV+7LvHnz0LBhQ4wZM0b+jitASkoKdDodunXrZm7r3r07dDqd3c/2/PnzyMrKsvh8tVotevXqZX7PoUOHUF5ebrFMQkIC2rVrZ15GyrYdeeihhxATE4N77rkHX3zxRZ3e6w2fQ0pKCnr16mVxgeEBAwYgIyMDaWlpTq/Dkz+HlJQUtGvXDgkJCRb7aDAYcOjQIYt+1eXnwV1/i53Zrlyfr7uwALQjOzsblZWViI2NtWiPjY1FVlaWzfdkZWXZXL6iogLZ2dkOl7G3Tjm4al/qW33tx8svv4xGjRrhT3/6kzwdr8GV+6HX6xEaGgqNRoPBgwdjxYoV6Nevn/w7cYOr9uXnn3/GmjVrsHr1atd0XAGysrIQExNj1R4TE+PwswXg8N8jKysLGo0GERERDpep67ZtCQ0NxbJly/DFF1/g22+/Rd++fTFs2DB89NFHTq/DGz4Hez/z1fvqzDo8+XOwtY8RERHQaDTmZaT8PLjrb7Ez25Xr83UXf7duXQFUKpXFcyGEVVtty9dsr+s65eKKfXEHV+7HkiVL8Omnn2LPnj0IDAyUobf2uWI/wsLC8Ntvv6GwsBDff/89pkyZgubNm6N3797yddzJvkndl+vXr+Ppp5/G6tWrER0dLX9n3WzOnDmYO3euw2UOHDgAwPbPqDPjhZQxpuYyUrddXXR0NF588UXz8y5duiAvLw9LlizBmTNnfOZzsNcXAFizZg169Ojh8L1K+RxqW8bRz8PTTz8t6z7I9bdYrmVqqq+/+46wALQjOjoaarXaqkK/evWqVbVfJS4uzuby/v7+iIqKcriMvXXKwVX7Ut9cvR+vv/46Fi5ciJ07d6JDhw7ydr4aV+6Hn58fWrZsCQDo2LEjTpw4gUWLFrmsAHTFvhw7dgxpaWkYMmSI+XWj0QgA8Pf3x6lTp9CiRQuZ96T+TJgwAU888YTDZZo2bYo//vgDV65csXrt2rVrDj9bwHTUIT4+3txe/d8jLi4OZWVlyMvLszgqcfXqVXMhEhcXV+dtO6t79+744IMPfOpzsPczDwB/+9vfMG3aNIfvV8LnEBcXh19++cXi9by8PJSXlzv8rKp+Huxx199iZ7Yr1+frLoyA7dBoNOjcuTN27Nhh0b5jxw67/2jJyclWy2/fvh1dunRBQECAw2Vc+YPgqn2pb67cj6VLl2L+/Pn47rvv0KVLF/k7X019/nsIIWAwGG6903a4Yl9at26NI0eO4LfffjM/HnroIfTp0we//fYbEhMTXbY/9SE6OhqtW7d2+AgMDERycjL0er3FJYl++eUX6PV6u59ts2bNEBcXZ/H5lpWV4YcffjC/p3PnzggICLBYJjMzE0ePHjUvI2XbzkpNTUV8fLxPfQ7Jycn48ccfLS6Jsn37diQkJOCuu+7yis8hOTkZR48eRWZmpsU+arVadO7c2e5nU/XzYI+7/hY7s125Pl+3cfk0EwWrmgK+Zs0acfz4cTF58mQREhIi0tLShBBCvPzyy2L48OHm5aumnr/44ovi+PHjYs2aNVZTz3/++WehVqvF4sWLxYkTJ8TixYvr9TIwcu6LwWAQqampIjU1VcTHx4tp06aJ1NRUcfr0aUXtx2uvvSY0Go344osvLC5NcP36dUXtx8KFC8X27dvF2bNnxYkTJ8Qbb7wh/P39xerVq122H67al5p8cRawEKZLb3To0EGkpKSIlJQU0b59e6vLftx+++1i8+bN5ueLFy8WOp1ObN68WRw5ckQ8+eSTNi9L0bhxY7Fz505x+PBhcf/999u87Edt2z59+rRITU0VY8eOFa1atTKPB1WX9Vi3bp34+OOPxfHjx8XJkyfF0qVLRUBAgFi2bJlPfQ75+fkiNjZWPPnkk+LIkSNi8+bNIjw8XNJlYDz1c6i6DEzfvn3F4cOHxc6dO0Xjxo0tLgMj9efBXX+La9uunJ+vO7AArMXKlStFUlKS0Gg04q677hI//PCD+bURI0aIXr16WSy/Z88e0alTJ6HRaETTpk3Fu+++a7XOzz//XNx+++0iICBAtG7dWmzatMnVuyGEkH9fqqb013zUXI+n70dSUpLN/Zg9e7ai9mPmzJmiZcuWIjAwUERERIjk5GSLSw+4kit+T6rz1QIwJydHPPXUUyIsLEyEhYWJp556SuTl5VksA0CsXbvW/NxoNIrZs2eLuLg4odVqxX333SeOHDli8Z6SkhIxYcIEERkZKYKCgsSDDz4oLl68WOdt9+rVy+bvzvnz54UQpj/4bdq0EcHBwSIsLEx07txZfPjhhz73OQghxB9//CF69uwptFqtiIuLE3PmzHH6EjBK+RwuXLggBg8eLIKCgkRkZKSYMGGCxaVvbuXnwV1/ix1tVwj5Pl93UAlRx0uRExEREZGi8RxAIiIiIh/DApCIiIjIx7AAJCIiIvIxLACJiIiIfAwLQCIiIiIfwwKQiIiIyMewACQiIiLyMSwAiYiIiHwMC0AiIiIiH8MCkIiIiMjHsAAkIiIi8jEsAImIiIh8DAtAIiIiIh/DApCIiIjIx7AAJCIiIvIxLACJiIiIfAwLQCIiIiIf8/8BrMZhvWblHFIAAAAASUVORK5CYII=", "text/html": [ "\n", "
\n", "
\n", " Figure\n", "
\n", " \n", "
\n", " " ], "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "u0.change_scales(1)\n", "h0.change_scales(1)\n", "\n", "fig, ax = plt.subplots(1,2)\n", "ax[0].plot(u0['g'][0][0], lat[0]*180/np.pi, '.-', color='C0')\n", "ax[0].set_title('Zonal velocity')\n", "ax[0].set_ylabel('latitude (deg)')\n", "ax[1].plot(h0['g'][0], lat[0]*180/np.pi, '.-', color='C1', label='height')\n", "ax[1].set_title('Height perturbation')\n", "plt.tight_layout()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Finding the most unstable mode (EVP)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Instead of prescribing some simple pertubation, we can use Dedalus to solve eigenvalue problems to find the most unstable linear eigenmode of the balanced jet.\n", "To do this, we create new fields representing perturbations to the jet and the eigenvalue, and build an EVP for linearized the equations around the background state.\n", "\n", "We need to pass the eigenvalue field $\\sigma$ to the `EVP` class on instantiation, and can use it multiplicatively in the equations." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Perturbation fields\n", "u1 = dist.VectorField(coords, name='u1', bases=full_basis)\n", "h1 = dist.Field(name='h1', bases=full_basis)\n", "sigma = dist.Field(name='sigma') # eigenvalue\n", "\n", "# Parameters\n", "nu = 1e5 * meter**2 / second / 32**2 # Hyperdiffusion constant\n", "\n", "# Eigenvalue problem\n", "problem = d3.EVP([u1, h1], eigenvalue=sigma, namespace=locals())\n", "problem.add_equation(\"sigma*u1 + u1@grad(u0) + u0@grad(u1) + nu*lap(lap(u1)) + g*grad(h1) + 2*Omega*zcross(u1) = 0\")\n", "problem.add_equation(\"sigma*h1 + div(h0*u1) + div(h1*u0) + nu*lap(lap(h1)) + H*div(u1) = 0\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The eigenvalue solver (like all solvers in Dedalus) splits up the problem into its *linearly separable subproblems*.\n", "Each of these *subproblems* is represented by a separate matrix -- these are the things constructed when you build a solver object, and together these form the diagonal blocks of the total linear system.\n", "Here we have NCCs (the background flow) that depend on the latitude, meaning the system is linearly coupled over all $\\ell$ for each $m$.\n", "\n", "In any event, the eigenvalue solver allows you to find the eigenvalues for each subproblem independently.\n", "The subproblems are objects in the `solver.subproblems` list.\n", "Each has a `.group` attribute that describes the corresponding mode (wavenumber or spherical harmonic order/degree).\n", "You can use the `solver.subproblems_by_group` dictionary to find the object associated with a given horizontal mode, here in the form `(m, None)` to indicate the matrices that couple all $\\ell$ for a given $m$.\n", "\n", "Here lets loop over the subproblems and compute the fastest growing mode for the first 15 spherical harmonic orders." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/kburns/Software/miniforge3/envs/d3/lib/python3.11/site-packages/scipy/sparse/_index.py:137: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n", " self._set_arrayXarray_sparse(i, j, x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2024-09-27 10:17:59,429 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.4e-01/s\n", "2024-09-27 10:18:03,268 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.0e-01/s\n", "2024-09-27 10:18:06,936 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.3e-01/s\n", "2024-09-27 10:18:10,634 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.2e-01/s\n", "2024-09-27 10:18:14,256 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.0e-01/s\n", "2024-09-27 10:18:17,957 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.9e-01/s\n", "2024-09-27 10:18:21,557 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.2e-01/s\n", "2024-09-27 10:18:25,295 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.0e-01/s\n", "2024-09-27 10:18:30,580 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 3s, Remaining: 0s, Rate: 3.1e-01/s\n", "2024-09-27 10:18:36,035 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.6e-01/s\n", "2024-09-27 10:18:39,657 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.1e-01/s\n", "2024-09-27 10:18:43,815 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.2e-01/s\n", "2024-09-27 10:18:47,492 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.9e-01/s\n", "2024-09-27 10:18:50,831 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.1e-01/s\n", "2024-09-27 10:18:55,194 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.0e-01/s\n", "2024-09-27 10:18:58,282 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.3e-01/s\n" ] } ], "source": [ "# Eigenvalue solver\n", "m_max = 15\n", "growth_rate = np.zeros(m_max + 1)\n", "solver = problem.build_solver()\n", "for m in range(m_max + 1):\n", " sp = solver.subproblems_by_group[(m, None)]\n", " solver.solve_dense(sp)\n", " growth_rate[m] = np.max(solver.eigenvalues.real)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now let's plot the growth rates:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'growth rate')" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "d74317dc3361425687bfd69baf13dc48", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9k0lEQVR4nO3dfXxU5Z3///dkcjMWk5EEc4eIAWwlBquZFAwadVcJNxZFsaIVsdrVhmIhifoLN1oKrqZaVy1VwiJhW5afSrcIhW2kpN5ElHSjSUCBVG2NBHHyjYE6Qd3cMDnfP/hmyjgDBMjJzOS8no/HPB7ONdec8zkac965rnOuYzMMwxAAAAAsIyrUBQAAAKB/EQABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALCY6FAXEMm6u7v16aefKj4+XjabLdTlAACAXjAMQ4cOHVJ6erqioqw5FkYAPA2ffvqphg0bFuoyAADAKdi3b5/OOeecUJcREgTA0xAfHy/pyA9QQkJCiKsBAAC90dbWpmHDhvnO41ZEADwNPdO+CQkJBEAAACKMlS/fsubENwAAgIURAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFkMABAAAsBgWggYwIHm7DdU0HlTLoXYlxzs0NiNR9ijrLvoKAEcjAAIYcLbscmvJ5j1ye9p9bWlOhxZPzdSkrLQQVgYA4YEpYAADypZdbs1eW+cX/iSp2dOu2WvrtGWXO0SVAUD4IAACMJW321D13w7o9zv2q/pvB+TtNkzd15LNexRsDz1tSzbvMbUGAIgETAEDME1/T8XWNB4MGPk7miHJ7WlXTeNB5Y5M6vP9A0CkYAQQgClCMRXbcujY4e9U+gHAQEUABNDnQjUVmxzv6NN+ADBQEQAB9LmTmYrtS2MzEpXmdOhYi73YdGQKemxGYp/uFwAiDQEQQJ8L1VSsPcqmxVMzJSkgBPa8Xzw1k/UAAVgeARBAnwvlVOykrDSVzcxWqtN/26lOh8pmZrMOIAAoggLg8uXLlZGRIYfDIZfLpW3bth23f1VVlVwulxwOh0aMGKEVK1b4fX7VVVfJZrMFvK699lozDwOwhFBPxU7KStObJf+sF+6+VL+85WK9cPelerPknwl/APD/REQAXLdunQoLC7Vo0SLV19crLy9PkydPVlNTU9D+jY2NmjJlivLy8lRfX6+FCxdq7ty5Wr9+va/PSy+9JLfb7Xvt2rVLdrtd3/ve9/rrsIABKxymYu1RNuWOTNL1Fw9V7sgkpn0B4Cg2wzDCfkXUcePGKTs7W2VlZb620aNHa9q0aSotLQ3oX1JSok2bNqmhocHXVlBQoJ07d6q6ujroPp5++mn99Kc/ldvt1qBBg3pVV1tbm5xOpzwejxISEk7yqICBj0eyAQhHnL8jYCHozs5O1dbWav78+X7t+fn52r59e9DvVFdXKz8/369t4sSJKi8vV1dXl2JiYgK+U15erltuuaXX4Q/AiU3KStOEzFTVNB5Uy6F2JccfmfZlNA4AQivsA2Bra6u8Xq9SUlL82lNSUtTc3Bz0O83NzUH7Hz58WK2trUpL8x95qKmp0a5du1ReXn7cWjo6OtTR0eF739bWdjKHAlhSz1QsACB8RMQ1gJJks/mPGBiGEdB2ov7B2qUjo39ZWVkaO3bscWsoLS2V0+n0vYYNG9bb8gEAAMJG2AfAIUOGyG63B4z2tbS0BIzy9UhNTQ3aPzo6WklJ/iMRX331lV588UX9y7/8ywlrWbBggTwej++1b9++kzwaAACA0Av7ABgbGyuXy6XKykq/9srKSo0fPz7od3JzcwP6b926VTk5OQHX//32t79VR0eHZs6cecJa4uLilJCQ4PcCAACINGEfACWpuLhYq1at0urVq9XQ0KCioiI1NTWpoKBA0pGRuVmzZvn6FxQUaO/evSouLlZDQ4NWr16t8vJy3X///QHbLi8v17Rp0wJGBgEAAAaqsL8JRJJmzJihAwcOaOnSpXK73crKylJFRYWGDx8uSXK73X5rAmZkZKiiokJFRUV69tlnlZ6ermXLlmn69Ol+2/3ggw/05ptvauvWrf16PAAAAKEUEesAhivWEQIAIPJw/o6QKWAAAAD0HQIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFhMRKwDCACRwtttqKbxoFoOtSs53qGxGYmyRx37ueUAEAoEQADoI1t2ubVk8x65Pe2+tjSnQ4unZmpSVloIKwMAf0wBA0Af2LLLrdlr6/zCnyQ1e9o1e22dtuxyh6gyAAhEAASA0+TtNrRk8x4Fe6xST9uSzXvk7ebBSwDCAwEQAE5TTePBgJG/oxmS3J521TQe7L+iAOA4CIAAcJpaDh07/J1KPwAwGwEQAE5TcryjT/sBgNkIgABwmsZmJCrN6dCxFnux6cjdwGMzEvuzLAA4JgIgAJwme5RNi6dmSlJACOx5v3hqJusBAggbBEAA6AOTstJUNjNbqU7/ad5Up0NlM7NZBxBAWGEhaMAieEKF+SZlpWlCZir/ngGEPQIgYAE8oaL/2KNsyh2ZFOoyAOC4mAIGBjieUAEA+DoCIDCA8YQKAEAwBEBgAOMJFQCAYAiAwADGEyoAAMEQAIEBjCdUAACCIQACAxhPqAAABEMABAYwnlABAAiGAAgMcDyhAgDwdSwEDVgAT6gAAByNAAhYBE+oAAD0YAoYAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxBEAAAACLiZgAuHz5cmVkZMjhcMjlcmnbtm3H7V9VVSWXyyWHw6ERI0ZoxYoVAX0+//xzzZkzR2lpaXI4HBo9erQqKirMOgQAAICwEBEBcN26dSosLNSiRYtUX1+vvLw8TZ48WU1NTUH7NzY2asqUKcrLy1N9fb0WLlyouXPnav369b4+nZ2dmjBhgj7++GP97ne/0/vvv6/nnntOQ4cO7a/DAgAACAmbYRhGqIs4kXHjxik7O1tlZWW+ttGjR2vatGkqLS0N6F9SUqJNmzapoaHB11ZQUKCdO3equrpakrRixQr94he/0F/+8hfFxMScUl1tbW1yOp3yeDxKSEg4pW0AAID+xfk7AkYAOzs7VVtbq/z8fL/2/Px8bd++Peh3qqurA/pPnDhR77zzjrq6uiRJmzZtUm5urubMmaOUlBRlZWXp0UcfldfrNedAAAAAwkTYPwu4tbVVXq9XKSkpfu0pKSlqbm4O+p3m5uag/Q8fPqzW1lalpaXpo48+0quvvqrbbrtNFRUV+vDDDzVnzhwdPnxYP/3pT4Nut6OjQx0dHb73bW1tp3l0AAAA/S/sRwB72Gw2v/eGYQS0naj/0e3d3d1KTk7WypUr5XK5dMstt2jRokV+08xfV1paKqfT6XsNGzbsVA8HAAAgZMI+AA4ZMkR2uz1gtK+lpSVglK9Hampq0P7R0dFKSkqSJKWlpemb3/ym7Ha7r8/o0aPV3Nyszs7OoNtdsGCBPB6P77Vv377TOTQAAICQCPsAGBsbK5fLpcrKSr/2yspKjR8/Puh3cnNzA/pv3bpVOTk5vhs+LrvsMv31r39Vd3e3r88HH3ygtLQ0xcbGBt1uXFycEhIS/F4AAACRJuwDoCQVFxdr1apVWr16tRoaGlRUVKSmpiYVFBRIOjIyN2vWLF//goIC7d27V8XFxWpoaNDq1atVXl6u+++/39dn9uzZOnDggObNm6cPPvhAf/jDH/Too49qzpw5/X58AAAA/SnsbwKRpBkzZujAgQNaunSp3G63srKyVFFRoeHDh0uS3G6335qAGRkZqqioUFFRkZ599lmlp6dr2bJlmj59uq/PsGHDtHXrVhUVFemiiy7S0KFDNW/ePJWUlPT78QEAAPSniFgHMFyxjhAAAJGH83eETAEDAACg7xAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALCY61AUAAE6ft9tQTeNBtRxqV3K8Q2MzEmWPsoW6LABhigAIABFuyy63lmzeI7en3deW5nRo8dRMTcpKC2FlAMIVU8AAEMG27HJr9to6v/AnSc2eds1eW6ctu9whqgxAOCMAAkCE8nYbWrJ5j4wgn/W0Ldm8R97uYD0AWBkBEAAiVE3jwYCRv6MZktyedtU0Huy/ogBEBAIgAESolkPHDn+n0g+AdRAAASBCJcc7+rQfAOsgAAJAhBqbkag0p0PHWuzFpiN3A4/NSOzPsgBEAAIgAEQoe5RNi6dmSlJACOx5v3hqJusBAghAAASACDYpK01lM7OV6vSf5k11OlQ2M5t1AAEExULQABDhJmWlaUJmKk8CAdBrBEAAGADsUTbljkwKdRkAIgRTwAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAi4mYALh8+XJlZGTI4XDI5XJp27Ztx+1fVVUll8slh8OhESNGaMWKFX6f//rXv5bNZgt4tbfz0HQAADCwRUQAXLdunQoLC7Vo0SLV19crLy9PkydPVlNTU9D+jY2NmjJlivLy8lRfX6+FCxdq7ty5Wr9+vV+/hIQEud1uv5fDwUPTAQDAwGYzDMMIdREnMm7cOGVnZ6usrMzXNnr0aE2bNk2lpaUB/UtKSrRp0yY1NDT42goKCrRz505VV1dLOjICWFhYqM8///yU62pra5PT6ZTH41FCQsIpbwcAAPQfzt8RMALY2dmp2tpa5efn+7Xn5+dr+/btQb9TXV0d0H/ixIl655131NXV5Wv74osvNHz4cJ1zzjn67ne/q/r6+uPW0tHRoba2Nr8XcLK83Yaq/3ZAv9+xX9V/OyBvd9j/DQYAGGDC/lFwra2t8nq9SklJ8WtPSUlRc3Nz0O80NzcH7X/48GG1trYqLS1NF1xwgX79619rzJgxamtr0y9/+Utddtll2rlzp84///yg2y0tLdWSJUv65sBgSVt2ubVk8x65Pf+41jTN6dDiqZmalJUWwsoAAFYS9iOAPWw2/4eaG4YR0Hai/ke3X3rppZo5c6a+/e1vKy8vT7/97W/1zW9+U7/61a+Ouc0FCxbI4/H4Xvv27TvVw4EFbdnl1uy1dX7hT5KaPe2avbZOW3a5Q1QZAMBqwn4EcMiQIbLb7QGjfS0tLQGjfD1SU1OD9o+OjlZSUvCHpUdFRek73/mOPvzww2PWEhcXp7i4uJM8AuDItO+SzXsUbLLXkGSTtGTzHk3ITJU96th/2AAA0BfCfgQwNjZWLpdLlZWVfu2VlZUaP3580O/k5uYG9N+6datycnIUExMT9DuGYWjHjh1KS2MaDn2vpvFgwMjf0QxJbk+7ahoP9l9RAADLCvsAKEnFxcVatWqVVq9erYaGBhUVFampqUkFBQWSjkzNzpo1y9e/oKBAe/fuVXFxsRoaGrR69WqVl5fr/vvv9/VZsmSJ/vjHP+qjjz7Sjh079MMf/lA7duzwbRPoSy2Here+ZG/7AQBwOsJ+CliSZsyYoQMHDmjp0qVyu93KyspSRUWFhg8fLklyu91+awJmZGSooqJCRUVFevbZZ5Wenq5ly5Zp+vTpvj6ff/657rnnHjU3N8vpdOqSSy7RG2+8obFjx/b78WHgS47v3fqSve0HAMDpiIh1AMMV6wiht7zdhi5/7FU1e9qDXgdok5TqdOjNkn/mGkAAMBnn7wiZAgYinT3KpsVTMyUdCXtH63m/eGom4Q8A0C8IgEA/mZSVprKZ2Up1+k/zpjodKpuZzTqAAIB+ExHXAAIDxaSsNE3ITFVN40G1HGpXcrxDYzMSGfkDAPQrAiDQz+xRNuWODL4eJQAA/YEpYAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDGmBsDPP/9cq1at0oIFC3Tw4EFJUl1dnfbv32/mbgEAAHAc0WZt+N1339U111wjp9Opjz/+WHfffbcSExO1YcMG7d27V2vWrDFr1wAAADgO00YAi4uL9YMf/EAffvihHA6Hr33y5Ml64403zNotAAAATsC0APj222/rRz/6UUD70KFD1dzcbNZuAQAAcAKmBUCHw6G2traA9vfff19nn322WbsFAADACZgWAK+//notXbpUXV1dkiSbzaampibNnz9f06dPN2u3AAAAOAHTAuATTzyhzz77TMnJyfrf//1fXXnllRo1apTi4+P1yCOPmLVbAAAAnIBpdwEnJCTozTff1Kuvvqq6ujp1d3crOztb11xzjVm7BAAAQC/YDMMwzNjwmjVrNGPGDMXFxfm1d3Z26sUXX9SsWbPM2G2/amtrk9PplMfjUUJCQqjLAQAAvcD528QAaLfb5Xa7lZyc7Nd+4MABJScny+v1mrHbfsUPEAAAkYfzt4nXABqGIZvNFtD+ySefyOl0mrVbAAAAnECfXwN4ySWXyGazyWaz6eqrr1Z09D924fV61djYqEmTJvX1bgEAANBLfR4Ap02bJknasWOHJk6cqDPPPNP3WWxsrM477zyWgQEAAAihPg+AixcvliSdd955mjFjht9j4AAAABB6pi0Dc8cdd5i1aQAAAJwG024C8Xq9euKJJzR27FilpqYqMTHR73Wyli9froyMDDkcDrlcLm3btu24/auqquRyueRwODRixAitWLHimH1ffPFF2Ww23/Q1AADAQGZaAFyyZImefPJJ3XzzzfJ4PCouLtaNN96oqKgo/exnPzupba1bt06FhYVatGiR6uvrlZeXp8mTJ6upqSlo/8bGRk2ZMkV5eXmqr6/XwoULNXfuXK1fvz6g7969e3X//fcrLy/vVA4TAAAg4pi2DuDIkSO1bNkyXXvttYqPj9eOHTt8bX/+85/1/PPP93pb48aNU3Z2tsrKynxto0eP1rRp01RaWhrQv6SkRJs2bVJDQ4OvraCgQDt37lR1dbWvzev16sorr9Sdd96pbdu26fPPP9fGjRt7XRfrCAEAEHk4f5s4Atjc3KwxY8ZIks4880x5PB5J0ne/+1394Q9/6PV2Ojs7VVtbq/z8fL/2/Px8bd++Peh3qqurA/pPnDhR77zzjrq6unxtS5cu1dlnn60f/vCHva4HAAAg0pkWAM855xy53W5J0qhRo7R161ZJ0ttvvx3weLjjaW1tldfrVUpKil97SkqKmpubg36nubk5aP/Dhw+rtbVVkvTWW2+pvLxczz33XK9r6ejoUFtbm98LAAAg0pgWAG+44Qa98sorkqR58+bpoYce0vnnn69Zs2bprrvuOuntff2pIsd60sjx+ve0Hzp0SDNnztRzzz2nIUOG9LqG0tJSOZ1O32vYsGEncQQAAADhwbRlYH7+85/7/vmmm27SsGHD9NZbb2nUqFG67rrrer2dIUOGyG63B4z2tbS0BIzy9UhNTQ3aPzo6WklJSdq9e7c+/vhjTZ061fd5d3e3JCk6Olrvv/++Ro4cGbDdBQsWqLi42Pe+ra2NEAgAACKOKQGwq6tL99xzjx566CGNGDFC0pEbOcaNG3fS24qNjZXL5VJlZaVuuOEGX3tlZaWuv/76oN/Jzc3V5s2b/dq2bt2qnJwcxcTE6IILLtB7773n9/mDDz6oQ4cO6Ze//OUxQ11cXNxJTV8DAACEI1OmgGNiYrRhw4Y+215xcbFWrVql1atXq6GhQUVFRWpqalJBQYGkIyNzs2bN8vUvKCjQ3r17VVxcrIaGBq1evVrl5eW6//77JUkOh0NZWVl+r7POOkvx8fHKyspSbGxsn9UOAAAQbkybAr7hhhu0ceNGvynTUzVjxgwdOHBAS5culdvtVlZWlioqKjR8+HBJktvt9lsTMCMjQxUVFSoqKtKzzz6r9PR0LVu2jGcQAwAAyMR1AB955BE98cQTuvrqq+VyuTRo0CC/z+fOnWvGbvsV6wgBABB5OH+bGAAzMjKOvVObTR999JEZu+1X/AABABB5OH+bOAXc2Nho1qYBAABwGkxbBxAAAADhiQAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGJMuwtYkj7//HPV1NSopaXF96zdHkc/uQMAAAD9x7QAuHnzZt1222368ssvFR8fL5vN5vvMZrMRAAEAAELEtCng++67T3fddZcOHTqkzz//XH//+999r4MHD5q1WwAAAJyAaQFw//79mjt3rr7xjW+YtQsAAACcAtMC4MSJE/XOO++YtXkAAACcoj69BnDTpk2+f7722mv1wAMPaM+ePRozZoxiYmL8+l533XV9uWsAAAD0ks0wDKOvNhYV1bsBRZvNJq/X21e7DRkeJg0AQOTh/N3HI4BfX+oFAAAA4ce0awDXrFmjjo6OgPbOzk6tWbPGrN0CAADgBPp0CvhodrtdbrdbycnJfu0HDhxQcnIyU8AAACAkOH+bOAJoGIbf4s89PvnkEzmdTrN2CwAAgBPo8yeBXHLJJbLZbLLZbLr66qsVHf2PXXi9XjU2NmrSpEl9vVsAQAh4uw3VNB5Uy6F2Jcc7NDYjUfaowD/+AYSXPg+A06ZNkyTt2LFDEydO1Jlnnun7LDY2Vuedd56mT5/e17sFAPSzLbvcWrJ5j9yedl9bmtOhxVMzNSkrLYSVATgR064B/M1vfqNbbrlFcXFxZmw+LHANAQCr2rLLrdlr6/T1E0jP2F/ZzGxCIMIW528TrwH84IMP9MYbb+irr74yaxcAgBDwdhtasnlPQPiT5GtbsnmPvN2mjC8A6AOmBcDa2lpNnz5dgwcPVm5urhYsWKAtW7boiy++MGuXAIB+UNN40G/a9+sMSW5Pu2oaD/ZfUQBOimkBcMuWLfr73/+u119/Xddff73q6+s1Y8YMJSYm6tJLLzVrtwAAk7UcOnb4O5V+APpfn98EcjS73a7c3FwlJiZq8ODBio+P18aNG/W3v/3NzN0CAEyUHO/o034A+p9pI4BlZWW65ZZblJaWpry8PG3dulV5eXmqra3VZ599ZtZuAQAmG5uRqDSnQ8da7MWmI3cDj81I7M+yAJwE00YA58yZo7PPPlv33XefCgoKLHuXDQAMNPYomxZPzdTstXWySX43g/SEwsVTM1kPEAhjpo0AvvTSS7rtttv04osvKjk5WePGjVNJSYlefvllbgQBgAg3KStNZTOzler0n+ZNdTpYAgaIAKatA3g0j8ejbdu26Xe/+52ef/552Ww2dXR0mL1b07GOEACr40kgiEScv02+CeTgwYOqqqrS66+/rtdff127du1SUlKSrrzySjN3CwDoJ/Yom3JHJoW6DAAnybQAeNFFF2nPnj1KTEzUFVdcobvvvltXXXWVsrKyzNolAAAAesG0AHjPPfcQ+AAAAMKQaQHw3nvv9f1zz2WGNhvXhQAAAISaaXcBS9KaNWs0ZswYnXHGGTrjjDN00UUX6T//8z/N3CUAAABOwLQRwCeffFIPPfSQ7r33Xl122WUyDENvvfWWCgoK1NraqqKiIrN2DQAAgOMwbQTwV7/6lcrKyvTYY4/puuuu0/XXX6/HH39cy5cv17Jly056e8uXL1dGRoYcDodcLpe2bdt23P5VVVVyuVxyOBwaMWKEVqxY4ff5Sy+9pJycHJ111lkaNGiQLr74YkYnAQCAJZgWAN1ut8aPHx/QPn78eLnd7pPa1rp161RYWKhFixapvr5eeXl5mjx5spqamoL2b2xs1JQpU5SXl6f6+notXLhQc+fO1fr16319EhMTtWjRIlVXV+vdd9/VnXfeqTvvvFN//OMfT+5AAQAAIoxpC0FnZWXp+9//vhYuXOjX/q//+q9at26d3nvvvV5va9y4ccrOzlZZWZmvbfTo0Zo2bZpKS0sD+peUlGjTpk1qaGjwtRUUFGjnzp2qrq4+5n6ys7N17bXX6uGHH+5VXSwkCQBA5OH8beI1gEuWLNGMGTP0xhtv6LLLLpPNZtObb76pV155Rb/97W97vZ3Ozk7V1tZq/vz5fu35+fnavn170O9UV1crPz/fr23ixIkqLy9XV1eXYmJi/D4zDEOvvvqq3n//fT322GO9rg0AACASmRYAp0+frpqaGj355JPauHGjDMNQZmamampqdMkll/R6O62trfJ6vUpJSfFrT0lJUXNzc9DvNDc3B+1/+PBhtba2Ki3tyDMqPR6Phg4dqo6ODtntdi1fvlwTJkw4Zi0dHR1+j7Bra2vr9XEAAACEC1MCYFdXl+655x499NBDWrt2bZ9s8+trCBqGcdx1BYP1/3p7fHy8duzYoS+++EKvvPKKiouLNWLECF111VVBt1laWqolS5ac4hEAAACEB1NuAomJidGGDRv6ZFtDhgyR3W4PGO1raWkJGOXrkZqaGrR/dHS0kpL+8czKqKgojRo1ShdffLHuu+8+3XTTTUGvKeyxYMECeTwe32vfvn2ncWQAAAChYdpdwDfccIM2btx42tuJjY2Vy+VSZWWlX3tlZWXQu4wlKTc3N6D/1q1blZOTE3D939EMw/Cb4v26uLg4JSQk+L0AAAAijWnXAI4aNUoPP/ywtm/fLpfLpUGDBvl9Pnfu3F5vq7i4WLfffrtycnKUm5urlStXqqmpSQUFBZKOjMzt379fa9askXTkjt9nnnlGxcXFuvvuu1VdXa3y8nK98MILvm2WlpYqJydHI0eOVGdnpyoqKrRmzRq/O40BAAAGItMC4KpVq3TWWWeptrZWtbW1fp/ZbLaTCoAzZszQgQMHtHTpUrndbmVlZamiokLDhw+XdGTNwaPXBMzIyFBFRYWKior07LPPKj09XcuWLdP06dN9fb788kv9+Mc/1ieffKIzzjhDF1xwgdauXasZM2ac5pEDAACEN9PWAbQC1hECACDycP428RpAAAAAhCfTpoCLi4uDtttsNjkcDo0aNUrXX3+9EhMTzSoBAAAAQZg2BfxP//RPqqurk9fr1be+9S0ZhqEPP/xQdrtdF1xwgd5//33f00EyMzPNKMF0DCEDABB5OH+bOAV8/fXX65prrtGnn36q2tpa1dXVaf/+/ZowYYJuvfVW7d+/X1dccYWKiorMKgEAAABBmDYCOHToUFVWVgaM7u3evVv5+fnav3+/6urqlJ+fr9bWVjNKMB1/QQAAEHk4f5s4AujxeNTS0hLQ/tlnn/meoXvWWWeps7PTrBIAAAAQhKlTwHfddZc2bNigTz75RPv379eGDRv0wx/+UNOmTZMk1dTU6Jvf/KZZJQAAACAI06aAv/jiCxUVFWnNmjU6fPiwJCk6Olp33HGHnnrqKQ0aNEg7duyQJF188cVmlGA6hpABAIg8nL/7YSHoL774Qh999JEMw9DIkSN15plnmrm7fsUPEAAAkYfzt4nrAPY488wzddFFF5m9GwAAAPQSTwIBAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiTF8HEAhX3m5DNY0H1XKoXcnxDo3NSJQ9yhbqsgAAMB0BEJa0ZZdbSzbvkdvT7mtLczq0eGqmJmWlhbAyAADMxxQwLGfLLrdmr63zC3+S1Oxp1+y1ddqyyx2iygAA6B8EQFiKt9vQks17FOwB2D1tSzbvkbfb1EdkAwAQUgRAWEpN48GAkb+jGZLcnnbVNB7sv6IAAOhnBEBYSsuhY4e/U+kHAEAkIgDCUpLjHX3aDwCASEQAhKWMzUhUmtOhYy32YtORu4HHZiT2Z1kAAPQrAiAsxR5l0+KpmZIUEAJ73i+emsl6gACAAY0ACMuZlJWmspnZSnX6T/OmOh0qm5nNOoAAgAGPhaBhSZOy0jQhM5UngQAALIkACMuyR9mUOzIp1GUAANDvmAIGAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiIiYALl++XBkZGXI4HHK5XNq2bdtx+1dVVcnlcsnhcGjEiBFasWKF3+fPPfec8vLyNHjwYA0ePFjXXHONampqzDwEAACAsBARAXDdunUqLCzUokWLVF9fr7y8PE2ePFlNTU1B+zc2NmrKlCnKy8tTfX29Fi5cqLlz52r9+vW+Pq+//rpuvfVWvfbaa6qurta5556r/Px87d+/v78OCwAAICRshmEYoS7iRMaNG6fs7GyVlZX52kaPHq1p06aptLQ0oH9JSYk2bdqkhoYGX1tBQYF27typ6urqoPvwer0aPHiwnnnmGc2aNatXdbW1tcnpdMrj8SghIeEkjwoAAIQC5+8IGAHs7OxUbW2t8vPz/drz8/O1ffv2oN+prq4O6D9x4kS988476urqCvqdr776Sl1dXUpMTDxmLR0dHWpra/N7AQAARJqwD4Ctra3yer1KSUnxa09JSVFzc3PQ7zQ3Nwftf/jwYbW2tgb9zvz58zV06FBdc801x6yltLRUTqfT9xo2bNhJHg0AAEDohX0A7GGz2fzeG4YR0Hai/sHaJenxxx/XCy+8oJdeekkOh+OY21ywYIE8Ho/vtW/fvpM5BAAAgLAQHeoCTmTIkCGy2+0Bo30tLS0Bo3w9UlNTg/aPjo5WUlKSX/sTTzyhRx99VH/605900UUXHbeWuLg4xcXFncJRAAAAhI+wHwGMjY2Vy+VSZWWlX3tlZaXGjx8f9Du5ubkB/bdu3aqcnBzFxMT42n7xi1/o4Ycf1pYtW5STk9P3xQMAAIShsA+AklRcXKxVq1Zp9erVamhoUFFRkZqamlRQUCDpyNTs0XfuFhQUaO/evSouLlZDQ4NWr16t8vJy3X///b4+jz/+uB588EGtXr1a5513npqbm9Xc3Kwvvvii348PAACgP4X9FLAkzZgxQwcOHNDSpUvldruVlZWliooKDR8+XJLkdrv91gTMyMhQRUWFioqK9Oyzzyo9PV3Lli3T9OnTfX2WL1+uzs5O3XTTTX77Wrx4sX72s5/1y3EBAACEQkSsAxiuWEcIAIDIw/k7QqaAAQAA0HcIgAAAABYTEdcAYmDzdhuqaTyolkPtSo53aGxGouxRx17jEQAAnB4CIEJqyy63lmzeI7en3deW5nRo8dRMTcpKC2FlAAAMXEwBI2S27HJr9to6v/AnSc2eds1eW6ctu9whqgwAgIGNAIiQ8HYbWrJ5j4Ldgt7TtmTzHnm7uUkdAIC+RgBESNQ0HgwY+TuaIcntaVdN48H+KwoAAIsgACIkWg4dO/ydSj8AANB7BECERHK8o0/7AQCA3iMAIiTGZiQqzenQsRZ7senI3cBjMxL7sywAACyBAIiQsEfZtHhqpiQFhMCe94unZrIeIAAAJiAAImQmZaWpbGa2Up3+07ypTofKZmazDiAAACZhIWiE1KSsNE3ITOVJIAAA9CMCIELOHmVT7sikUJcBAIBlMAUMAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABYTMQFw+fLlysjIkMPhkMvl0rZt247bv6qqSi6XSw6HQyNGjNCKFSv8Pt+9e7emT5+u8847TzabTU8//bSJ1QMAAISPiAiA69atU2FhoRYtWqT6+nrl5eVp8uTJampqCtq/sbFRU6ZMUV5enurr67Vw4ULNnTtX69ev9/X56quvNGLECP385z9Xampqfx0KAABAyNkMwzBCXcSJjBs3TtnZ2SorK/O1jR49WtOmTVNpaWlA/5KSEm3atEkNDQ2+toKCAu3cuVPV1dUB/c877zwVFhaqsLDwpOpqa2uT0+mUx+NRQkLCSX0XAACEBufvCBgB7OzsVG1trfLz8/3a8/PztX379qDfqa6uDug/ceJEvfPOO+rq6jKtVgAAgEgQHeoCTqS1tVVer1cpKSl+7SkpKWpubg76nebm5qD9Dx8+rNbWVqWlpZ1SLR0dHero6PC9b2trO6XtAAAAhFLYjwD2sNlsfu8NwwhoO1H/YO0no7S0VE6n0/caNmzYKW8LAAAgVMI+AA4ZMkR2uz1gtK+lpSVglK9Hampq0P7R0dFKSko65VoWLFggj8fje+3bt++UtwUAABAqYR8AY2Nj5XK5VFlZ6ddeWVmp8ePHB/1Obm5uQP+tW7cqJydHMTExp1xLXFycEhIS/F4AAACRJuyvAZSk4uJi3X777crJyVFubq5WrlyppqYmFRQUSDoyMrd//36tWbNG0pE7fp955hkVFxfr7rvvVnV1tcrLy/XCCy/4ttnZ2ak9e/b4/nn//v3asWOHzjzzTI0aNar/DzIMeLsN1TQeVMuhdiXHOzQ2I1H2qFOfMgcAAOEpIpaBkY4sBP3444/L7XYrKytLTz31lK644gpJ0g9+8AN9/PHHev311339q6qqVFRUpN27dys9PV0lJSW+wChJH3/8sTIyMgL2c+WVV/pt53gG0m3kW3a5tWTzHrk97b62NKdDi6dmalLWqd00AwBAOBpI5+9TFTEBMBwNlB+gLbvcmr22Tl//QegZ+yubmU0IBAAMGAPl/H06wv4aQJjL221oyeY9AeFPkq9tyeY98nbzdwIAAAMFAdDiahoP+k37fp0hye1pV03jwf4rCgAAmIoAaHEth44d/k6lHwAACH8EQItLjnf0aT8AABD+CIAWNzYjUWlOh4612ItNR+4GHpuR2J9lAQAAExEALc4eZdPiqZmSFBACe94vnprJeoAAAAwgBEBoUlaaymZmK9XpP82b6nSwBAwAAANQRDwJBOablJWmCZmpPAkEAAALIADCxx5lU+7IpFCXAQAATMYUMAAAgMUQAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFsM6gACAiOPtNli4HjgNBEAAQETZssutJZv3yO1p97WlOR1aPDWTR1cCvcQUMAAgYmzZ5dbstXV+4U+Smj3tmr22Tlt2uUNUGRBZCIAAgIjg7Ta0ZPMeGUE+62lbsnmPvN3BegA4GgEQABARahoPBoz8Hc2Q5Pa0q6bxYP8VBUQoAiAAICK0HDp2+DuVfoCVEQABABEhOd7Rp/0AKyMAAgAiwtiMRKU5HTrWYi82HbkbeGxGYn+WBUQkAiAAICLYo2xaPDVTkgJCYM/7xVMzWQ8Q6AUCIAAgYkzKSlPZzGylOv2neVOdDpXNzGYdQKCXWAgaABBRJmWlaUJmKk8CAU4DARAAEHHsUTbljkwKdRlAxGIKGAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAi2EdwDDk7TZY4BQAAJgmYkYAly9froyMDDkcDrlcLm3btu24/auqquRyueRwODRixAitWLEioM/69euVmZmpuLg4ZWZmasOGDWaV32tbdrl1+WOv6tbn/qx5L+7Qrc/9WZc/9qq27HKHujQAsDxvt6Hqvx3Q73fsV/XfDsjbbQzo/WLgiogRwHXr1qmwsFDLly/XZZddpn//93/X5MmTtWfPHp177rkB/RsbGzVlyhTdfffdWrt2rd566y39+Mc/1tlnn63p06dLkqqrqzVjxgw9/PDDuuGGG7RhwwbdfPPNevPNNzVu3Lj+PkRJR8Lf7LV1+vr/1s2eds1eW8dzLgEghLbscmvJ5j1ye9p9bWlOhxZPzTT1d3Oo9iuFbkaKmTDz2QzDCPs/I8aNG6fs7GyVlZX52kaPHq1p06aptLQ0oH9JSYk2bdqkhoYGX1tBQYF27typ6upqSdKMGTPU1taml19+2ddn0qRJGjx4sF544YVe1dXW1ian06lPP/1UCQkJp3p4ko78sE9Y9mf9n0MdQT+3SUpJiNPWn1zK/wQA0M8qGz5T0e92B/yB3vPb+KmbLtSE0WcPmP327Lv0j3/1Oy+lxMdpwcRRpu2zv/bb1tam9PR0eTye0z5/R6qwHwHs7OxUbW2t5s+f79een5+v7du3B/1OdXW18vPz/domTpyo8vJydXV1KSYmRtXV1SoqKgro8/TTTx+zlo6ODnV0/OMHsq2tTZKUnp5+MocUVNywMUr9fmCY7WFIam7rUPKF49Wx773T3h8AoJdsURpaUC57/BDZbP5/gBuSDKNb965+XftX/FAyuiN/v5LO+Gauzp628EgZR+27ue1/Ne+/dumzjY/qfz+o7tN9hnK/VhT21wC2trbK6/UqJSXFrz0lJUXNzc1Bv9Pc3By0/+HDh9Xa2nrcPsfapiSVlpbK6XT6XsOGDTuVQwrKfubgPu0HAOgbcedcqOiEswNCWA+bLUrRCWcr7pwLB8R+ZYtS4tX3/L99+O/bZouSZBz53NbHESJU+7WosB8B7BHw149hHPN/imP1/3r7yW5zwYIFKi4u9r1va2vTsGHD+mQKuObjv+vO/9x5wn7//V//v8aeRwgEgP7yh13/R//fhoYT9lvzXxt1bVbKCfuF+35PdD7qCZ6vvre3T89H/bnfnilgKwv7ADhkyBDZ7faAkbmWlpaAEbweqampQftHR0crKSnpuH2OtU1JiouLU1xcXED7oEGDNGjQoF4dz7FcMfobSnO+r2ZPe8C1HtKR6z1SnQ5dMXoo1wACQD8aNsTZ636ney4Ih/22dX3ey362iN2v1+s9re8PBGE/jhobGyuXy6XKykq/9srKSo0fPz7od3JzcwP6b926VTk5OYqJiTlun2Nt02z2KJsWT82U9I+Le3v0vF88NZPwBwD9bGxGotKcjoDfzT1sOnJX7tiMxAGx3+R4R5/2C/f9WlXYB0BJKi4u1qpVq7R69Wo1NDSoqKhITU1NKigokHRkanbWrFm+/gUFBdq7d6+Ki4vV0NCg1atXq7y8XPfff7+vz7x587R161Y99thj+stf/qLHHntMf/rTn1RYWNjfh+czKStNZTOzler0/+FOdTpYAgYAQiRUf6CHar9WC7xWFRHLwEhHFoJ+/PHH5Xa7lZWVpaeeekpXXHGFJOkHP/iBPv74Y73++uu+/lVVVSoqKtLu3buVnp6ukpISX2Ds8bvf/U4PPvigPvroI40cOVKPPPKIbrzxxl7X1LMMTF/fRs76RwAQfqy0DmDPurSS/C5L6jkTmTUo0V/7Nev8HUkiJgCGI36AAMBarLQw8kAOvJy/CYCnhR8gAMBANlADL+fvCLgLGAAAhIY9yqbckUmW2a+VRMRNIAAAAOg7BEAAAACLIQACAABYDAEQAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxPAnkNPQ8Ra+trS3ElQAAgN7qOW9b+Wm4BMDTcOjQIUnSsGHDQlwJAAA4WYcOHZLT6Qx1GSFhM6wcf09Td3e3Pv30U8XHx8tm69uHY7e1tWnYsGHat2+fJR5UzfEObBzvwMbxDmwD8XgNw9ChQ4eUnp6uqChrXg3HCOBpiIqK0jnnnGPqPhISEgbM/3C9wfEObBzvwMbxDmwD7XitOvLXw5qxFwAAwMIIgAAAABZDAAxTcXFxWrx4seLi4kJdSr/geAc2jndg43gHNqsdr1VwEwgAAIDFMAIIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgGFo+fLlysjIkMPhkMvl0rZt20JdkilKS0v1ne98R/Hx8UpOTta0adP0/vvvh7qsflNaWiqbzabCwsJQl2Ka/fv3a+bMmUpKStI3vvENXXzxxaqtrQ11WaY5fPiwHnzwQWVkZOiMM87QiBEjtHTpUnV3d4e6tD7xxhtvaOrUqUpPT5fNZtPGjRv9PjcMQz/72c+Unp6uM844Q1dddZV2794dmmL7wPGOt6urSyUlJRozZowGDRqk9PR0zZo1S59++mnoCj5NJ/rve7Qf/ehHstlsevrpp/utPvQtAmCYWbdunQoLC7Vo0SLV19crLy9PkydPVlNTU6hL63NVVVWaM2eO/vznP6uyslKHDx9Wfn6+vvzyy1CXZrq3335bK1eu1EUXXRTqUkzz97//XZdddpliYmL08ssva8+ePfq3f/s3nXXWWaEuzTSPPfaYVqxYoWeeeUYNDQ16/PHH9Ytf/EK/+tWvQl1an/jyyy/17W9/W88880zQzx9//HE9+eSTeuaZZ/T2228rNTVVEyZM8D03PdIc73i/+uor1dXV6aGHHlJdXZ1eeuklffDBB7ruuutCUGnfONF/3x4bN27U//zP/yg9Pb2fKoMpDISVsWPHGgUFBX5tF1xwgTF//vwQVdR/WlpaDElGVVVVqEsx1aFDh4zzzz/fqKysNK688kpj3rx5oS7JFCUlJcbll18e6jL61bXXXmvcddddfm033nijMXPmzBBVZB5JxoYNG3zvu7u7jdTUVOPnP/+5r629vd1wOp3GihUrQlBh3/r68QZTU1NjSDL27t3bP0WZ6FjH+8knnxhDhw41du3aZQwfPtx46qmn+r029A1GAMNIZ2enamtrlZ+f79een5+v7du3h6iq/uPxeCRJiYmJIa7EXHPmzNG1116ra665JtSlmGrTpk3KycnR9773PSUnJ+uSSy7Rc889F+qyTHX55ZfrlVde0QcffCBJ2rlzp958801NmTIlxJWZr7GxUc3NzX6/v+Li4nTllVda4veXdOR3mM1mG7Cj3N3d3br99tv1wAMP6MILLwx1OThN0aEuAP/Q2toqr9erlJQUv/aUlBQ1NzeHqKr+YRiGiouLdfnllysrKyvU5ZjmxRdfVF1dnd5+++1Ql2K6jz76SGVlZSouLtbChQtVU1OjuXPnKi4uTrNmzQp1eaYoKSmRx+PRBRdcILvdLq/Xq0ceeUS33nprqEszXc/vqGC/v/bu3RuKkvpVe3u75s+fr+9///tKSEgIdTmmeOyxxxQdHa25c+eGuhT0AQJgGLLZbH7vDcMIaBto7r33Xr377rt68803Q12Kafbt26d58+Zp69atcjgcoS7HdN3d3crJydGjjz4qSbrkkku0e/dulZWVDdgAuG7dOq1du1bPP/+8LrzwQu3YsUOFhYVKT0/XHXfcEery+oUVf391dXXplltuUXd3t5YvXx7qckxRW1urX/7yl6qrqxvw/z2tgingMDJkyBDZ7faA0b6WlpaAv6oHkp/85CfatGmTXnvtNZ1zzjmhLsc0tbW1amlpkcvlUnR0tKKjo1VVVaVly5YpOjpaXq831CX2qbS0NGVmZvq1jR49ekDe0NTjgQce0Pz583XLLbdozJgxuv3221VUVKTS0tJQl2a61NRUSbLc76+uri7dfPPNamxsVGVl5YAd/du2bZtaWlp07rnn+n5/7d27V/fdd5/OO++8UJeHU0AADCOxsbFyuVyqrKz0a6+srNT48eNDVJV5DMPQvffeq5deekmvvvqqMjIyQl2Sqa6++mq999572rFjh++Vk5Oj2267TTt27JDdbg91iX3qsssuC1jW54MPPtDw4cNDVJH5vvrqK0VF+f9atdvtA2YZmOPJyMhQamqq3++vzs5OVVVVDcjfX9I/wt+HH36oP/3pT0pKSgp1Saa5/fbb9e677/r9/kpPT9cDDzygP/7xj6EuD6eAKeAwU1xcrNtvv105OTnKzc3VypUr1dTUpIKCglCX1ufmzJmj559/Xr///e8VHx/vGzlwOp0644wzQlxd34uPjw+4vnHQoEFKSkoakNc9FhUVafz48Xr00Ud18803q6amRitXrtTKlStDXZpppk6dqkceeUTnnnuuLrzwQtXX1+vJJ5/UXXfdFerS+sQXX3yhv/71r773jY2N2rFjhxITE3XuueeqsLBQjz76qM4//3ydf/75evTRR/WNb3xD3//+90NY9ak73vGmp6frpptuUl1dnf77v/9bXq/X9zssMTFRsbGxoSr7lJ3ov+/XA25MTIxSU1P1rW99q79LRV8I7U3ICObZZ581hg8fbsTGxhrZ2dkDdlkUSUFf//Ef/xHq0vrNQF4GxjAMY/PmzUZWVpYRFxdnXHDBBcbKlStDXZKp2trajHnz5hnnnnuu4XA4jBEjRhiLFi0yOjo6Ql1an3jttdeC/j97xx13GIZxZCmYxYsXG6mpqUZcXJxxxRVXGO+9915oiz4NxzvexsbGY/4Oe+2110Jd+ik50X/fr2MZmMhmMwzD6KesCQAAgDDANYAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAMu76qqr9JOf/ESFhYUaPHiwUlJStHLlSn355Ze68847FR8fr5EjR+rll18OdakA0CcIgAAg6Te/+Y2GDBmimpoa/eQnP9Hs2bP1ve99T+PHj1ddXZ0mTpyo22+/XV999VWoSwWA02YzDMMIdREAEEpXXXWVvF6vtm3bJknyer1yOp268cYbtWbNGklSc3Oz0tLSVF1drUsvvTSU5QLAaWMEEAAkXXTRRb5/ttvtSkpK0pgxY3xtKSkpkqSWlpZ+rw0A+hoBEAAkxcTE+L232Wx+bTabTZLU3d3dr3UBgBkIgAAAABZDAAQAALAYAiAAAIDFcBcwAACAxTACCAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBi/i/BLXKmcXSxJgAAAABJRU5ErkJggg==", "text/html": [ "\n", "
\n", "
\n", " Figure\n", "
\n", " \n", "
\n", " " ], "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.axhline(0, c='k')\n", "plt.plot(np.arange(m_max + 1), growth_rate, 'o')\n", "plt.xlabel('m')\n", "plt.ylabel('growth rate')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We see that the jet is unstable to a range of low-m perturbations, but it's most unstable for $m=7$.\n", "We can solve again for that subproblem and set the perturbation variables to the most unstable mode using the `solver.set_state` method:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/kburns/Dropbox/Git/dedalus3master/dedalus/tools/array.py:353: ComplexWarning: Casting complex values to real discards the imaginary part\n", " dest[:] = src\n" ] } ], "source": [ "sp = solver.subproblems_by_group[(7, None)]\n", "solver.solve_dense(sp)\n", "index = np.argmax(solver.eigenvalues.real)\n", "solver.set_state(index, sp.subsystems[0])" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot e.g. the height perturbation in the eigenmode:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "526649613ba040c6818246659378fb51", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACPaElEQVR4nO3deXwV1d0/8M/cPTtLICFlVVGhWKVgaXDDtgJuhdrWrUWolurjgog+WooLLoCiUn6Ka2td2lp9+ijuC/hUUSq4UNyR1hYFkcgWkpDlbnN+f8yZMzP33iQ3ydzcm9zP+/W6Jpk7d+Z875kZD/M954wmhBAgIiIiorzhyXYBiIiIiKh7sQFIRERElGfYACQiIiLKM2wAEhEREeUZNgCJiIiI8gwbgERERER5hg1AIiIiojzDBiARERFRnmEDkIiIiCjPsAFIRERElGfYACQiIiLKM2wAEhEREeUZNgCJiIiI8owv2wXoCXRdx1dffYWSkhJompbt4hAREXWKEAINDQ2oqqqCx9M994BaWloQiURc214gEEAoFHJte/mKDcA0fPXVVxgyZEi2i0FEROSKbdu2YfDgwRnfT0tLC0YMK0bNzrhr26ysrMSWLVvYCOwiNgDTUFJSAgAYcv3V8HT3ASeyeMdRE9nbdzbjBvI39mzGDeRv7NmKO9v1DeRn7Fk8zvWWFmxbeKP6/1qmRSIR1OyM44sNw1Fa0vU7jvUNOoaN+xyRSIQNwC5iAzANZtrXEwqxAdhd2ADMjmw3CPI19nxsBJnyMfZsX9+Abu/OVFyiobik6/vUkf3vrrdgA5Byk3lxzoELZV4RWm40Cqh7mOcX67x75eH1LS50xF04zOJC7/pGCABHARMRERHlHd4BJCIioozSIaCj67cA3dgGGdgAzHWJqZnuTBnY95WtFJEmspMmyWZqzL7PfI09m3Hby9Fd8r3Os13fZjm6W7ZizwIdOtxI3rqzFQKYAiYiIiLKO7wDSERERBkVFwJx0fW7rG5sgwxsAHaE0Ny7Xd/ZdIMbaYrOxOBmmqKjMbiZmuloHNmMu7OfSdTZGPK1zjv7mUTZPMdN2Tree1p9d+Vzdtk6x015kk4md7ABSERERBnFQSC5hw1AIiIiyigdAnE2AHMKB4EQERER5RneAexu+T7jfrbiZ9+Y7Mjzx21lTb4+Ui+f5fjxzhRw7mEDkIiIiDKKo4BzD1PARERERHmGdwCJiIgoo3T5cmM75I4efQcwFovh6quvxogRI1BQUIADDjgAN9xwA3TdOkSEEFi4cCGqqqpQUFCASZMm4eOPP85iqYmIiPJLXI4CduNF7ujRDcBbbrkF9957L1asWIFNmzZh6dKluPXWW3HnnXeqdZYuXYply5ZhxYoVeOedd1BZWYkTTjgBDQ0NWSw5ERERUfb06BTwunXrMG3aNJx88skAgOHDh+Mvf/kL3n33XQDG3b/ly5djwYIFOO200wAADz/8MCoqKvDoo4/i/PPPz1rZiYiI8kVcGC83tkPu6NF3AI8++mj83//9H/75z38CAN5//32sXbsWJ510EgBgy5YtqKmpweTJk9VngsEgjjvuOLz55ptZKTMRERFRtvXoO4BXXXUV6urqcOihh8Lr9SIej2PRokU466yzAAA1NTUAgIqKCsfnKioq8MUXX7S63XA4jHA4rP6ur6/PQOmJiIjyAweB5J4e3QB8/PHH8ac//QmPPvoovvnNb+K9997D3LlzUVVVhZkzZ6r1NM05QaYQImmZ3ZIlS3D99ddnrNxERET5RIeGOLo+WbXuwjbI0KNTwP/93/+NX//61zjzzDNx2GGHYcaMGbjsssuwZMkSAEBlZSUA606gaefOnUl3Be3mz5+Puro69dq2bVvmgiAiIiLqZj26AdjU1ASPxxmC1+tV08CMGDEClZWVWL16tXo/EolgzZo1mDhxYqvbDQaDKC0tdbyIiIioc3Th3ovc0aNTwKeeeioWLVqEoUOH4pvf/CY2btyIZcuW4dxzzwVgpH7nzp2LxYsXY+TIkRg5ciQWL16MwsJCnH322VkuPRERUX6Iu5QCdmMbZOjRDcA777wT11xzDS688ELs3LkTVVVVOP/883Httdeqda688ko0NzfjwgsvRG1tLSZMmIBVq1ahpKQkiyUnIiIiyh5NCD5ZuT319fUoKyvDsJsXwRMKdW1jWg583SKL/4LKVvzZjNmUj7Fn83jPdp3na+zZvsbla+xpxq23tOCLXy9AXV1dt3RvMv//+ebHg1Bc0vVeZ/sbdEz85o5uK39v1qP7ABIRERFRx/XoFDARERHlPl1o0F24O+vGNsjABiARERFlFAeB5B6mgImIiIjyDO8AEhERUUbF4UHchXtOcRfKQgY2AImIiCijhEt9AAX7ALqGKWAiIiKiPMM7gN3N/NdLtufKIurtzHOMdwyIso6DQHIP7wASERER5RneASQiIqKMigsP4sKFQSBMnrmGDUAiIiLKKB0adBeSjjrYAnQLU8BEREREeYZ3AImIiCijOAgk97ABSERERBnlXh9ApoDdwhQwERERUZ7hHUAiIiLKKGMQSNfTt25sgwxsAHaEkK9UMnRMdmQOW9fnlk53exmIPd24Mzafdr7GnqW4s3qcq0KkW4AM7DpbdZ7FmFURcj12IGvHe1bPcZfpLj0LmKOA3cMUMBEREVGe4R1AIiIiyigOAsk9vANIRERElGd4B9Atqf5R0sF+I119Zn2qz6fdh6Qr/6jKYuxdilltpHP7TvnZbqzzLsfOOu9kAVIsSzOWrJ7jaiOd2XEry7N0vHdrfbf2+W6o86wf6y7S4eGTQHIMG4BERESUUXGhId7VfwHJ7ZA7mAImIiIiyjO8A9gRbU0Dk0hLsW4b/3DJ1D9qhNZOyiDxPTfurttjz1LMQBqpEjdjN2MRKZal2nU2Yxet/N5ZibFnMW6gG2O3x5IDsbcad0eXp8Ne52nE5XbsidvL2vGei/UNpI5TtPFeN4i7NA1MnClg17ABSERERBmlCw90F0YB6xwF7BqmgImIiIjyDO8AdoCmG6/Ubzr/FG2851gv1XtaO+8nrt7W7X60kTJIkRrROpMuaSv2VtZJe1vtfK6t8jriNlc0F6ZI53c49hQpQFXeVsrdajwdrHOg7fpMWecp1m/v2Gl95+2s3pH0V1eO946Ocret3+b5kHLHyesn1XeKVHB753ir66T6WDvHe0fSgh36HjtyjtsK4kbsnT5G3bq+AamP9za6uqQTd6vr2T/S0XM8aQfJ28rYk3TawRRw7uEdQCIiIqI8wzuARERElFE63JnCpbUkHHUc7wB2gJkCTvVCwkvTjVvtmplqbG8EsQYIj3xp8uURQBov4ZUvj1CfhflqjyyXo6z2GNJ4qbjNbaWK3b5LLSH1YY/dIxyvdmN3vIxtJMbu2JetTGYZU9Vhh+LWk7eZVp3DKqtZ3+nWuSN2e52nw1ZWVcYuxK6WtxG7o87NmFPUeYfiTqzzVPtKEbNZ54mxtxt33HjZ67yt4zwle313ts7VZ1Mf76ni7urxnvL6pqeob9FKHTjqvCvnuEh5fWuzXZGqnjpwfdNs9Z20rXSkOMc7HXsr1zcVf0JdtFbf2WBOBO3GqzPuvvtujBgxAqFQCOPGjcMbb7zR5vpr1qzBuHHjEAqFcMABB+Dee+9NWueJJ57A6NGjEQwGMXr0aKxcudLx/pIlS3DkkUeipKQEAwcOxPTp07F58+ZOlT8T2AAkIiKiXuvxxx/H3LlzsWDBAmzcuBHHHHMMTjzxRGzdujXl+lu2bMFJJ52EY445Bhs3bsRvfvMbzJkzB0888YRaZ926dTjjjDMwY8YMvP/++5gxYwZOP/10vPXWW2qdNWvW4KKLLsL69euxevVqxGIxTJ48GY2NjRmPOR2aEBxT3Z76+nqUlZVhxA2L4AmFUq6T9C9Q27/SEu98GL+LpPeStyHSv6sDyH/1afbNp+70nmIghGP9xGXt7TbVwIfE2M1/amgi5XclEr6XxO20vvPkgrQZu27FnrReqvXb2nWKeB1xJ34vrcRurS+c67dbAGdBUt2V0GzHRGKcjvUzFbtHtLp+yjpP9651QkFajd1cR09clry9Th3vqc5x9V4HYu/Q3VvnyvaY0j7H7T/RfuytHbfOc9e2zJP6PHYc54nbanXnyYVJ+/omf6aKvSODIVKdz62d4/YiJK+fmdgBQNO1tM5xvaUFn1+zAHV1dSgtLW2nAF1n/v9zxYYJKCjueq+z5v0xXDzurQ6Vf8KECfj2t7+Ne+65Ry0bNWoUpk+fjiVLliStf9VVV+GZZ57Bpk2b1LILLrgA77//PtatWwcAOOOMM1BfX48XX3xRrTN16lT07dsXf/nLX1KWY9euXRg4cCDWrFmDY489Nq2yZxL7AHaAFtegxVs5Y9v4n70mGz/C9n7iBdVxcVCNJaR3lbJvTJcXINvFsNWRyzI9YBUudYOh3SKkit3jXKTKk+rC6BEptmFtq02O/9HJ2GUjT+0qxSbsF0brO0jROGxLGw16zWPbbYo6Fx7n55CQxnS8lyjp4DHjtt4zY3KMDE2sTx0p/6eZduzmTzN2W523NdrZ+keBrfGTbn2bZVXlto53TXcet9b6KX5PsX66h5oKzWM7x9WxjISVbJ/32HbSRoOx3UKYJ5P5w3a8q9MsVUMnqRGspR17qw0dj+0cT1XnSf8QTHWud+AcBwAh1FeQss4Tj3fd/nvy+m3t3lHn9usVjDpPmtauteub+bcbsad7fZPvJV7fWv1/WIbp0KB36I5G69sBjIalXTAYRDAYTFo/Eolgw4YN+PWvf+1YPnnyZLz55psp97Fu3TpMnjzZsWzKlCl44IEHEI1G4ff7sW7dOlx22WVJ6yxfvrzVstfV1QEA+vXr1+o63YkpYCIiIupRhgwZgrKyMvVKdScPAHbv3o14PI6KigrH8oqKCtTU1KT8TE1NTcr1Y7EYdu/e3eY6rW1TCIF58+bh6KOPxpgxY9KKMdN4B5CIiIgyKi48iLvwJBBzG9u2bXOkgFPd/bPTNOfdRyFE0rL21k9c3pFtXnzxxfjggw+wdu3aNsvZnXr8HcDt27fj5z//Ofr374/CwkIcccQR2LBhg3pfCIGFCxeiqqoKBQUFmDRpEj7++OMslpiIiIi6orS01PFqrQFYXl4Or9ebdGdu586dSXfwTJWVlSnX9/l86N+/f5vrpNrmJZdcgmeeeQavvvoqBg8enHaMmdaj7wDW1tbiqKOOwvHHH48XX3wRAwcOxL///W/06dNHrbN06VIsW7YMDz30EA4++GDcdNNNOOGEE7B582aUlJR0aH9aFPB4ncsS+wBZ/T401bwWspOGo1+Y+XmzCe4RVnPc3nE+oTO1/R8XQiT8IjTV31ANdtCFKqTqH2TvC5TQN8bxtJM2Okzb+/sk9n8SHg0anP3CUvVTUX1jPLaYE/vGac6Y1WdTxW6W2yyPGae9z0uqvjG22NXPxA7Wtr8ddW6vbzMWGHWu+n7at5XY/81W56kGBiTGbmzL7M2tOcumaVYfUDhjcqyn4tScMZubTtEvLLF/l6PO1R/yh1ekrG+1Ha8ZuxWno58U0qjzhL6O0DXrPEvVzysxTnvsCd+LY5mt3EnnuK4lnbPqPZGiP6CW4hy3H+8pYk8aoieEFZ/50zzXdC11lzJbfRs/Yf1srX+ofZepYvdY33Fi3aWcniXluS6S1m/v+qbKrZvXNKvO1Wlg/rT3AUyM3Xact3W8W/09bbGYv9iuzSLh/wuOY9p+jqdxfXPUubDiU2U1f40nnP+2WFRMcS3peq7FkBXuPQmkY9sIBAIYN24cVq9ejR/96Edq+erVqzFt2rSUn6mursazzz7rWLZq1SqMHz8efr9frbN69WpHP8BVq1Zh4sSJ6m8hBC655BKsXLkSr732GkaMGNGhsmdaj24A3nLLLRgyZAgefPBBtWz48OHqdyEEli9fjgULFuC0004DADz88MOoqKjAo48+ivPPP7+7i0xERJR3dKFBd2Mi6E5sY968eZgxYwbGjx+P6upq3H///di6dSsuuOACAMD8+fOxfft2PPLIIwCMEb8rVqzAvHnzMHv2bKxbtw4PPPCAY3TvpZdeimOPPRa33HILpk2bhqeffhqvvPKKI8V70UUX4dFHH8XTTz+NkpISdcewrKwMBQUFXfkaXNGjU8DPPPMMxo8fj5/+9KcYOHAgxo4di9/97nfq/S1btqCmpsYxmicYDOK4445rdfQPAITDYdTX1zteRERE1POcccYZWL58OW644QYcccQReP311/HCCy9g2LBhAIAdO3Y45gQcMWIEXnjhBbz22ms44ogjcOONN+KOO+7Aj3/8Y7XOxIkT8dhjj+HBBx/Et771LTz00EN4/PHHMWHCBLXOPffcg7q6OkyaNAmDBg1Sr8cff7z7gm9Dj74D+J///Af33HMP5s2bh9/85jd4++23MWfOHASDQZxzzjmqtZ1qpM4XX3zR6naXLFmC66+/Pmm5JwZ4os5lifPcOf5WaQGZrtBSpMbM+/fyaR4AVOpQ8wgkThuhpcpT2Kd1MFMjcfNjmtV5NXH4v7CWecz149b27KmDVlPAjjSo9VN4E9JT9n9qJHxnztgT0iaacMacWADzp26bHiGesKNUKV3dit1cP93UmGPuRrOuzOLIOte9mlHficVOiE/F7RVW7I70UIov3l7fgFXnurDqWH0sOf2t2b6npNgFkqcLSdi9o/weJKf57cdEYrrMI1QKWE2n4bXS5eZO261zc/oXRx3K7yEhF+jo6mBbPyk1Zp8uJNXxnnCca5qtvh07A1KmQb3O+gbgrPN0znHY0vwqFtvOdNvv5mYTzgv109ENInXcqgi2+ravb1zjnHXs2IQZn1nnXmGLXe0h/eubjFFd3+JWnSeV3Zbu9yTE3uoUWAm7std5UlcX+zme6nixXd8SY7fnqROPdXXltD/iw1bnKmWc+L1oAlpC1wj7ca7izVIKWHcpBdzZJ4FceOGFuPDCC1O+99BDDyUtO+644/CPf/yjzW3+5Cc/wU9+8pNW38/1aZZ7dANQ13WMHz8eixcvBgCMHTsWH3/8Me655x6cc845ar2Ojv6ZP38+5s2bp/6ur6/HkCFDXC49ERFRftCFB7oLo4Dd2AYZevQ3OWjQIIwePdqxbNSoUepWbmVlJQB0aPQPYKSJE0cYEREREfUWPfoO4FFHHZX0YOV//vOfKq8/YsQIVFZWYvXq1Rg7diwAY1bwNWvW4JZbbunw/rwt1gBGU1JKTKUAbYPX5E/dljKCLS1kfE5A85upEeNeveYV8Hjl7ylSJGZ6TW0/7lFpIXg8snwaNPNfTLb0B2CkBBJTI1rMSHWb7xsbTp0aApwpQDM24VMZCxWomeoSXlsa1PwyfQKajNPjs2I3wkidAjZj11XKW4MeU7lz57o6kkZAQtfUaDiPLXbzu0hKm9i3Z6tzM3bdPJNsGV5dVprjsVjmd2XG6TNTgroVs63OE2MXQrOOK93YmB4z00QeWzYoOf1tTwsBRj0nfQdxtJ0OTEgFCq+t3s1jX9OSHgVmpUEFYIvZ+Cng8TqPfU/io8TgrHPzONflcQ6PBiHrPfGxf5qtq4Mj9sRj3979oa0UsC1e8zj3qHNc1rn9QqGOfQFN1rt57Kvj3SNUzKmOd3Wc6x7oSee4LDNgdYOwn7u2mFv7DtJOAZt1LI93IaxMq7ox47X9ro53GZv9XLfVeXv1DVtseswDzWOOArbVeeKTXeTf9jjt34GWeD1sJwUszJjN70DYskuy/CoDb+vqAH+K63k71zczbnWOqzrXgLj90VJwPhkkoatDquMckZS7y7g4NMRdeBKIG9sgQ49uAF522WWYOHEiFi9ejNNPPx1vv/027r//ftx///0AjJNz7ty5WLx4MUaOHImRI0di8eLFKCwsxNlnn53l0hMREeUHpoBzT49uAB555JFYuXIl5s+fjxtuuAEjRozA8uXL8bOf/Uytc+WVV6K5uRkXXnghamtrMWHCBKxatarDcwASERER9RY9ugEIAKeccgpOOeWUVt/XNA0LFy7EwoULu7wvXzPgjTuXJaY6dJke0PyALnMjekCu47WNXlTpJJkG8Qt4/MbGffKn16vD7zN+N1MkHlu6QE9IkURjXsRleiAWNQqia1aeSo0alJ/X4pqVGolYPz2p0oKJWQp7asSW9gYA3Q9AxmymRDQzbWb7rlQq0B+HR6ZJzNh9MlXi88bhTZEeissNx+Je+dODuEwBxyNeuS+507hmHyhnxGZLdavYo1bciaNjHfHbUoG6mXX2y2XyZxxWJloNFNSg0kJm6lfVeSAOr4rZ+On16o76BuRcWip2Wdcx+R1EvWqEnJkWQhzQ1GhoKy1kxq1itsWOttLfiel+r6xvAHFZ55oGK0Xnc+ZIhc9Kg3oD1vFur28Abda5rnsQNWM26zxqm4k3njhbu3Usp4rdngZPWd8A7CO+dVtXB7PezXNcjRD22jah6lyHZh7nAescBwC/z6r/ts7xWNyjjnkzZnNkpRAeaHHnEGhN16yYE+s6miIN3kqdJ6b5NZ8zbuM9s0uKNSTX6upg1rkOrzrHrXO9reubuqbZjvd4VKZG1crWMF01ENsWm9d2fQNkV5eELhEpz3Hb8W6mvc3jHX6oa1niIQcPrNSv/On1x+HzWfUNAB6PnvJYB4zj3X59A4xzPC5jELJwmtkPwTbps+rOErPFLutdb065u4yLw530bbz9VShNvJdKRERElGd6/B1AIiIiym3sA5h72ADsgMB+AW9AwD7hrZUWkqkaMyUUB2Jm6s8cqOi1pSw8VloIMFKBgaBx3z4UNO7Vh3wxBHzGMr/MYfk8Vo4mJvOPUZmTaon60RI1qjQstx8VQFyOEjVHm6nJQuNWSsTbIn9GrJSBSo3FRHJqyMwqegHdJ2OX6aC47R69NUmwTA9pttFxZtovGEcgIGMPWLEDQMAbUzGb6SFdaCr2SNyItyXmQ0vE+PLD6vuRhfRqQNRKibUau5kaixgxA1bq0J4eskb+ampkoBlzXLfCTXw+qPBYI2DN1K8/ZMUdlGmhkN8oiN8Td9S3EZNH1XckZuy8WdW5HxF5SpspQ0StXKRKg8o4veHk1JgnCnhkGtGRDrXVtxk7YBzvSceGxzYK3DbROWCkvs3Ub9A8zgMxFNhiBgC/N+6obzN2wDjewzLmlqhR5y0aEDVHi5r7sj0PObGrg7clOTXmiYlWU6GOrg5+mZoM2NZL+H4cz2w1R3z6dUd9A0DIL49zXwxBr/G7vc7N2KMyFRiO+9AsYw57jO8gAuPvuNAAc6J6WxpUxS5PDDNub8QZu4rbmUWG8Djr24zd/t2Y+1LMZfJ49waMmALBKILyXDfrPOiNwS/TwfYUsFnfYXmOm8d7S9SHFo9RkKiZIo9pVnJRN2Oy4jTPcVX/UcATFY5y2+s8sauD7tOsLg7265uZEk8cNe4VaoS3X17Xg4GYOreDfuu6nur6Bhh1bsYeNmP3+NT1LW72a1F1riVN7O+NGOe5I/YmZEVceBB3ofHmxjbIwG+SiIiIKM/wDiARERFllIBmPTKxi9shd7AB2AH+/Tp8fitPYKT5rHQQYKUAoWnwmOkDW1pQPRfSb6WFACAQjKGowLhXXxIwfhb5Iyj0Gff3AzI95IFQJ1HMTP3KNMF+TxxAyHjPHDVmu8erJsOV6VBP1EoP+Mw0aAvgC8uyyRSJJ2YbaWbG7rFSQmbsMfNZrB4jBe5cX6j3zNFxvpCV8i4OGQUplrEX+418RcgbhS8hH6dDQ0vcSAE1xYyfXk9ATRxrjhC1HqppmxjVnv5sccbua5HpqgjgiZgpMTNtbh+5nRy7kCMQ9RQpIZWx8AloCanfopARZ3EwjCIZs1nnIW8UnoThqBHdp+q7yWsORzTqXNc9iEYT8s6pRoGaqcCwMbLd+F3GHhXW6Fhbnav6NlPe5vEuNCsFaJ8MO2EUJfxmKjCOUMjYQVHQiLck2GKrb5kSNQsN6zm7Zrq/KeZHoxZwvBeNeREzR1MmTgQc1ZLS/b4WwGurb+OnUKNoHSO3Zfwq/Wl71rRK8yemjL3WOW6OAvUHYyi01TdgneuFvqiK2X68x+QXaB7v+6NW7lWNDI6Zo4GtPib2iY9VzGa9N5t/C3gj5rlufme2OjfT2n5b7OpaZjseElPGXgFhjnyV6f6ATPcXh8Iq9mLb8e6TOUvzeNehWd071PXNiF0I6xyPOU40M+aEc73FXu+2czzqvM7ZN5XUpcdvC9DsApIqZWzObBDQ4QvKLh0y9pJQWF3f7Od44vXNrPOI7lP1vV8LApAjg1V9Jxzv8RTneIv92i6vwU0phnp3A6aAcw+/SSIiIqI8wzuARERElFG60NQAl65uh9zBO4BEREREeYZ3ADvA26zDG9Otvi4+Tc1Kbk6tYu8rpaZFCcqfBQJ6gewbUyj7gBQa/WDKClvQJ2R0zin1y/5B/hbVL8hvm3sgLDtiNcaMDUdkxzP79Ciqn0jEC63FI8tv/PTJaQD8jdbvvmajf4ivRTj6gxkbsXWO8Tr7gpnxAnA8DUU3Y5Y/9ZDcRlEM/iLZB6zQiLNPqBllwRYZu/GzQM7JEvTE4JF9ZMz5n5p1v/pXYIs8hOO2p0OYT0gQETNuL3yNxvq+Rhl7U+rYAcAT1VW/IMeUP+ZUGHKZ7rP6QllThFhxxwtkfyZZ51phDKEi2e9N9vcsk3XeJ9CC4oTYU9W5Jy5UHyH1pATZFywS86qnoKDZK2PX4GuSse834zW/A2H1ibL1+9RitvoGAK8G3ey7Kh9vYu8jZ06FYz/e4zJmUSSfeFFo1XlpgRFn36A83gMtKJId8YLyePdouqpvM3ZTS9ynvgOzzqNRL/Sw8bvHPN7NuO11bdZ/s7D1AZTTcMThPNZhq3O/9btzCiT5u3m8h2x1Xiin/5B1XlrQgtKQEXsfebwXy455Rb6wo74BICq8aJZ9/2K61bczrjtjj8m40eKFt8U8x60695uxNzmPc1+LsKZ/idmOd1vfP+NNLel8V3Xut9V7oTxPC+PwFRmxF8pz3KzzPsFmlAZaZMzG91LgiSad42Hdp/oDmtc31cdX96rrmy6Pd0+zF94m5/XNPN79zbZzvNmscyt2e52rOja/DvPpJvYpv2z9YM36No93yOM9UBRBceI5Hmxx1DdgHe+AUd9GvHIqp6gneWqYmBcx8xxPON59jZoVu+2n34y5RZ6T2eoDCI96ak1Xt0PuYAOQiIiIMoop4NzDpjQRERFRnuEdwA7wRnR4dR16QKbgPBriQXlrvkD+LIb6af4eK5GpneIYQsVG2qNvsXGPvl9I/gw2oV/AyE8Ve5PTA+Z0F03xAJrlHBzm1BD7wsY0ILXNhWjYL6eB2S+nR2nwwd9gfNbfYGzLb6bBGgX8jWY6SKbBwjo8MXOuB7lzDdB9MmYzBSyzELGQhmihjL1Ixl4CREtkakTG7iuW00AUN6NvgZES6V9gxN4/0GilfuV8HanSn/tlPjmqe7E/avy+t6XQ+A4aC9C8X+bh9vtkvEYh/fuNl/m7EbsOf5MzNWLGbU+BWtPdaLYnQGhW7MXJ9Q4AsRIdolSmAItlKqiwBf0LjS+/X1DWu8zPFfnCKDTnKpF0aGiSdW1+By1xP+ojRpy1zUbsDY1GnUf3B+CRMfvqPVbsZn3vFzJ2q869YRl7NKHOYaXD4poHQh7z6ni31Xm0RMZuq3OPqm8j9r5Fsq5DTegfNApUIuenKPamTn+a9b1ffgdmndeGC1DbaMTeJOtc7PfD12DGnHC874eqa79Mf3mbdXijyfVtpT9lvH6zy4OGWMiM2apzVe+lMr1WIqexKY6gpMiIr588zstDjegbMH4vlbGH5Lwd9umdWmQ/gmjMq9KB9TL2fc0FqGsqAABE9hvfi1ZvrBNo8FhpT3Wci6TYParORfL0Tj4NcRm7eezHg6nPcQCIFQvESuT3KOu8qDiMvkXGOV5eYB7vxs8+/mZ1nKdKfzbLfHKz7lfHfp28vu1tMup8//4QdHV9k+d4g5by+gYAviYdvmYrZuOnLQ1qq/O4z6pvM3Yjbk3FHlN1LhCX9e2X1/XiIuN471/YiHJ5be8rr+ulvhYUyq4O9umdzHO7PmbUq/n0j/3RgDrH9zUa77XsDwANfhmzdY6bP9U53mTFrrp3yK4O8Yj1vXcnHR7oLtxzcmMbZGADkIiIiDIqLjTr8XVd3A65g01pIiIiojzDO4Ad4IkZN7HNFHDcryFaKGdtLzXWifQxfkbLdOh95Ei4UiMdMrBkPwYV1QMAvlFQZyzzG3/38+1HicdIC9lHxDXK4YV1cSMV0KL71cjAvWEjLbB7v5GT2F8fAvYZaZNgnVGuQB0QMHaBYL186kidTMs2xuBtkqMuzbSAbXp74ZWj70I+iEJjn+ZIuFiBHIlYrCFcJmOXP2NlcXj7GKmO/qVGGmRQsVGIQYX1GBQ0Yq/wGz/LvE0qFWaKymGGDfEQdseMfJOZItkfC2J3s5GP2dNg/AzXBeHbJ1O/dZr8DiDjFSp2/34rdk+LsU8tKtOPwpxSX4Pwy1GGISNu3e9TaUEzHRYp01TMkTI5Wq+PHL1cGsaA0v0yZiM3VVWwD5WyUAN8xrJClfK20jJNss73xQsRlbl2M020LxLCrkajvuvrjO9DrzPq3L/Pi4AZs1nndToCDWbsxj68TTLucBxa3Jl6FV4vRFCOpi6QdV5gS3/KlLe9zqNlxja0MmO7ZWVNqCyRMRcaBaqScQ8M1KOPfBp9ofm4Alj13SRTgLtiJSoF3CRHu++R6f7d9cVoqTeWeeqMMgb3aVbsdUZdBOtlOqwhBp88zr3Nxk8tEnPUNwAIvxd6UI7ANLs8+KyUd6RU1rs61wViZTLNX2bE0q/EOt4HF+wDAFXn5b4GlMhHUtjrGzDOa/Mc/zpqfLkx3Yt95jneZBzn++qKEJf17asz6imwT8ZdDwRl7IEG+cSZhjh8Zn23yNjtde6RcQbkaPpCH5CQ7o8UaVbMst6jfWTKuyyKkjLz+mbWeT0Ghaz6BoB+XiMNWuQJJ434bdSD2Bs34muOGLE1xgLYI9OfOxuM472p3vgutDofgvus65v5U53jZuyN1vGuznH79c08x4PyOPdqatRvTHbpiZSY57p1jsf6yGtIWRj9ZX1Xlcjrm7yuVwbr1bW9xGN8P/ZrnJnmb9ALsFvm06Py+6iPGinvXU3F2FtvfAexOuN499Z5EZDXNzP2oDredescl7F7wrbYpZhwXmu7CweB5B7eASQiIiLKM7wDSERERBklhEfd9e3qdsgdbAB2gBbVoQldpUHjIQ0xOTJMpQL7yhFx/cMY2NdIAY4o2wsAOKhoJw4M7QQADPHvAQAM9BrrlHhi6nZsVA5La9B92BU30h/mKDldWCND97XIlGiDkTJAbQDBvcZWgsbmUbBXILRHjkzcZ6SpvHVyJuD9TUDYWCZkWkjTPEDASE9ohTLFiAJoITM1Zo2ABYBIsS3t3d9ILRT0bcE3+u4DABxYYhTkoEIj7hHBnfiGrxYA0EemAAttD0MPy9v7e3Qjpu1aXzTI31WKJBJUo+LCMjXi3+NDcI/x2ZDxdSO0Vz6MfW8UvjqZXq+XsTc1AxGjvELIOpMpbwSD0IqN1IuQqUBofsivHbFCKxVo1rfoZ6Ry+/UzUl3DymoxstiIWdV5YA8qvWZaSI4elKMBdRj1DQA7ZZ3HhQe7NSM9ZI74rmspQEODrJdaI3ZV53uB0F5jewUy9sDeMLwyZq1Bzg7bIus8ZqUhNb+s84ICiBK5fZkOFV7NUd+ALRXYP4ZAf+O7rexjxHZg6W4cVLgLADAiaPw063yAt9FR34BR5/tk2rsmbmy4UQ+qyY/NEbDm6NeWhiC8e2Xq11bnoT3GdkO1ciLtWludN8l6Dxv1JIRu1XdAjqYtLIDwFMovxPihRnwXaojIka/hfrLO+kXRp59x/g4t2wcAOKjEiPfA0E4MC+w2vhef8b308UQc9Q0AjfK83iVToACwVw4zbY77URcxjv36/Ubs8X0BBPYYnwmq49yq80CtEZ9vn63OzfqOOlN/ms8HFBjb95TI4z3gVWlQawSs7frWX6Y/ZZ1X9G3A8FKjICOLzHN8l7q+DTBTv3KUtwfW9W2fTPfXxErRIoz6NDv4N8UC6vrWVC/LWGvVuRl7wV5Z53ti8NfK61uD7foGAJGI8/oGAMGAur6ZXQDiBT7r+qZSwMYq4b4CcXl9K+prbH9wn304sMSo45EFRuzDAubxvk9d34KaNeK3SV7D9sSNfX8V64s6+bvZzcMc9byvsUClfv17fDJ2IFQr61se72ad+/e1QNsvY280YrfXueYztqEFkBVxaIij6+lbN7ZBBjaliYiIiPIM7wASERFRRunCnQEcCVNXUhewAdgRmvEyJwSOBTXIgamIFctnpZYat9z7ljVhWKmR9hpVssP4GfoKB8gUwRA5IWiZR6bxtAK1m/26kV7xIooG3TlaMKz7sD8qR8o1y0lwm6yJYM2RYWaaILQ7iuBuIy3g2Wukq0S9kZIS+xuhR5yTD2s+PzxFRjpI88uUgQhZkwKbKTEz7iIgVizTjaXGtirL6lXq95vF2wEAhwa/AgAM9+1DhUy9FXuMjXjgQVQYcdYLI20SgbGtgBZX6W/z2ccN4RDCzXIi2P1yIth6TY36Vanf3TI1sqcRWm2DjN34qTc3q7QQZFrIE5IjSwFA/q5GiPqsSb+t2AV0OdlzqRwJOaR0HwDgkJKv8c2CLwEABwe+BgBU+VpQ7jHSO0EtBLtmEYZHxtworNGxahJsmQZtaAkiJuvbv985EjJUK1CwW45K3SOPoT37gX1Gfev7jXScqnNbGtRTKOvc5wOEUTY1EbBfQzxkpkJl7MVyQuXSqBrtfGCpkQ77ZvFXOCRoHPMHylRghUwrlnkK4JGJB10mQuv0FsRl7Hvk8a4LD5rNkZJy4uumZpmqbfDBX+8cCRnaq6NA1beR/lJ1XlcPvdn4PkTMSol5gjK1WCzTr0ErN2ZOeG5N+m3FHpfHe1FpMwaXGgU4pNRIAY4uMI73kYEaDJEPHh7gNbZboBWr7YflSMxa3ajrJj0Cr0wVmnXeFA+goUVOft5kHu8e+M0R3vuc6f7grjB8e+SswHWyzhv2Q5cpYCR0dfAUFFhpcPk8XOG1Jjw30/7xQiAmJ/nW5PWtnxzdP7x0L0YXm9c3I/YD/btR5TP2VarJFKYmu5BAx37dPL6N+tqnRVXDQE2KHAmq+kajPN5tdR6qlalfebwHdjXBUyuva/XGd2Ae7446N9P9okhd3wDjhE55jqs61xEsNcpdVWbU+cEluzCqyLiuWdc3owzlHj8KZXcC83iPihhq5bW9RcgJoTUdUbO+zeubPN4jTX541PXNKEewzkr9hnYZ5fHtlTNf19ZBb5Cxy+MdQreubzLdj6Ls5IB1l/oAurENMvCbJCIiIsozvANIREREGaVDU4877Op2yB1sAHZAvNAHzedTkyDHQ8YLAOJB+VzZoJGOKQxEUCyfb2tO9hvyRBFAHKlERQwx+d5+mQ7dp/uxR44ONCdK3RspQr0cJRYJyxRti7zFHwbkY4TV8x+9ER2a+Wxf3ZwM1ZwA12OlgGxpUM1MfwaDMu4AYgXGeuboODlHL+JBQITkaLSAUe4if0Q927dIjoQLyYlv/bYRcXGZkopDR7NMiTXIiVr3xo3hd3tixdgjJ0rdFzHyMo3hAHQZu0/G7g3bY3Y++1KL6UDChMdG7PJX+R1oQStuUWCkSeIFxn6iBR4rJWbWeUjAEzDrWz771m/EW+JtURO/ms+59UJTaU8z5a3LUaENehQNMtW9R44C3hsvwl75ANJaOSFwOOyD1uxVMRtxQsVtPutUTf4ai6lRzlbo8gKq+aGZI77NOg8FoYeM2O11LrNTkIN1VZ37AzEUB4yCFMtuDYWeCEJaVMZse6A0gKiIwyO/g7D8Dhr0OPbpRnx75AjY3dFi1Ebks1/D8jiUde4Na7bj3IpdHefqeHfGDdjq2ut1xAwAoiAAXdZ3XJ7j9jqPF8jvVp7jBcEoiv3yGbBe57N9A1ocXnMksZxwOoqYqu/9MgW4TzdHxBZgl5zwfI+s870tBWgOy3Rdi6zzFg1eVd/GT+vZvnHrOI9bsav6hozdTIMGg2oUsJCTfscKvI7rG2Cc43pQnuNBeY7LOi/xtSTF7reN8jbjNY/3GOJoUNc3+TzveLGa7H1PWE5qHw4i2uKT8VnnuPnTer6tHN0bj1sxy+9bxe31quubeY7bYzcne48VeNQz3dX1TdY5QnGE1DlufPFFvrDt+mY+09liXt+i8rreJKJokJszj/e9sWJV3/vkiO+mFqN+RNgLb4uWELuwXdfMczxhInvYZjQQGjRf4jnuBxHABiARERFlGJ8FnHvYACQiIqKM4iCQ3MNvkoiIiCjP8A5gB4RLfYj7fYgUm/2DoJ4OAZ/sdyL7uMV1jxrWv1c+LmR7tK/aVr1uDNcv8ljTsJjTnezT+wAAamJ98GWkHwBgS3O5say5BPXNsu9KWE7nEDNuiWs6YN4d19X0FV545LB/rzD616inPoQjVr8ROd0JggGIQtnnqlhOR1AWQKTU2Fe00NkPTg8IaB7nxEwR3Yt9UaOPy9c+4xECZj+4FuFDH/lwdL9mTXHTaPYH0vs7vqut4f7Y1mT8vrtZPjC+xQ9EZJ8es2ufSJ66I14k+0jGCuCV8WkFRrm89qci+GR/GdlHRi8KIlZq/B7uY3xXkRJryh+zzoVXqK8tLvtymU/sqI0V4StbfQNAg96g+oN6YcRu1nm9XoZd8VIAUHX+eUt/bG8yvj+zzmMtPkd9G78YP3SfBj0o+2/JOtfixVZ9m9OdmH2lNE099cWq8xCipTLmMlnnRRriKnZzuhCrzsMxOXVHzCjjzmip6hMVlf3O9ujGtCEhR53LJ9rofdR3tTVs1P/nzf3xdZNxvJp9oiD7gmkxDWbXQvVUHr+GuHxajafENp0PjD5v3oSnYMDrUX1cdXmcx0qCiKj6lue4OQ1IyIrZPN6F0NAUk/UdNVb82mvUlxc6GoX51Bf5hArEEJelatCN2HbFjDr/KtoXn7fIY7+pDwDjqS9mPzgtIuvcdqoJcwaXkCxrkR8+XU49Yk71UlTg6A8IAJDHgygIQpffVVQe75FSHyJF1tQ3gDzezfqWx1pUPqWlIRbCzogRg1nnceHBnrhR3/brG2Cc//vkOV4TNeLcGumPLU3G9c2s84bmEPSIsQ+fecjY6tycqiYeknEWh+SRBmhyOh/NfGyNEEZ9A+qpL6LQOsejZQEZu9dR3wAg/GbcAkKY57h1vH8ddV7fGnWjzks9Ler6ZtZ5k16InbJvs3m8bwv3x5fNxu97m4ydh1vk+RrRrOubpPuAuHmOF8r15IXJ4/VYU73EbB+U9Y2Q2a8ZwFZ0Ox2aO/MAchCIa9gAJCIioowSLo0CFmwAuoYpYCIiIqI8wzuAHRDu60Us4EVUZhbiBYDwmXkJuSwqU3rNIWz3GKmRiEyX7I4U4z/+gQCsaSOCHis1pWaEN1Nj0QI1LcLuFiM9sLepEC3mDPlRmRIzi+CzUheRYvkEC82nUkRemerwho1UoBbXrc/Kf1QJvwfxgEwpyekgokUeRIudT4LQzdS3BxAy/RmJGOXf02Q92L4xZqy4PdAHANDX14RCrzULvqlJzruwX/7cG5FTYYQLVeq3dr+x81iLD4iZU5lY5THLFtblkyZk6sNf6IVXPlFBTaFge56Q+cQLPWBOf+FBtNCMXU5fUmylQXXbLAp6zFhvv9x+jc9I8cR0D/bIaUw+9xvprTJfs6O+AatDc0M8pFKo5vQnO5uLUdts/G4+FUFEPVYqzEwBytkdjDqXU7fI9Kav2A+vTE9qURm7rc6FTI2ptFKBB9GihNiLrH2Y+4Ss82jUi30txhdjdn9oiftREzCO/T5+I91vpr79HisFHJZf5P54SHUZ2NUij/fmItQ2GsuizfIyZatzsw7sdS7kw+5jsu68fYxCeyM6tLizm4LwaND9Vn0DchoQGbN5/phpf90HdayJuLF+Y0sAu7xGec3U1r6oWef9UepzTo/i0XRV302yH0G93MHeaCF2txjH+W55/uxvCkLINKjHVmdmFwSzjGad634PfCXyiSFmnUR1R9rYiEWeHwGPOsfN7yxSpFnXNzMN6hPqmItFrOsbAGzXylS3h10R49j/t38giuW8Jfb6BoxrnHmOm9/VnkgRdpnneKOxLNLigxZ13p/QzcMgZKXohcectsdru76ZU8PYznF1nZDfVdCDmKx383uMFlndPMzjS2UsdQ0tYWPh7kbr+lYnj1vz+mbGXeiNqOubWectul+d4+b1bXe4SF0v65qMbZmpby2uqVs05vU2VqghLLen+4yFvmI5PVJLATzRFNc31TXG+FzY60zLdxdduJQC5ihg17ABSERERBnFUcC5h98kERERUZ7pVXcAlyxZgt/85je49NJLsXz5cgCAEALXX3897r//ftTW1mLChAm466678M1vfrPD24+UAN6glQ6zpwK1qEwjynRVo66plKiZ1ggF+iAkh7T5PfLpEB45i78m1K3tuExhhuM+NcKyWaYfIhEf4jJFYKbhzDS0+dQCe9mihRq8fWQayUwnydSIFkdyCthjpQzMbcQDVtpJpYVUKhDQ5OjMqGamrj1olCM3dwaMfFLQZ8Qb8MVU7D6PlQKOyZjN0YURGXdL1IeWiJypX456FjGPSjfqclRqrFCzpURlSqdEk3F7YGaiPDEZu21gpIrd/Fp9miN2Y5u2Ed9mJjKuQcincjQL46D4WnYBqG0sQEimhYKyzoPemKO+ATjqPCxHF5p13hLxqWMoZo4G1TU1GtWM06T7rCe1eGQ5vBGvFbMcGOgYSWqmmGx1bqab7LGb3435WbPOYx4/6mQa3Bytu6exyIrZ7zzeU9V5XHjQEvXL2L0ydj+i8skfIpaQCvQLxGWcahRwAIjK0fleOULczLZrcaHq25H+NmP3WbGbMZtPPFFPPvHIlByg6jyMIHbJ77muSaZEA8ao0AJ/FH5v3BF7qnPcPN5bYtZxbtZ5tMUHYe7TTGEGBWLyD2Eeo/K8D5dq8JrdQuQx5LHFbrLXufBZ35/5U8Vu/t9BAJr8TuNNxsIG2R+iJezHHpkSDfn7AAACvjgCXqPe7fUNGMd70jke86n0alTGrke9VlcHs4wyPQvNOkbNOvdErGPeEzPPFbm6cF7fjNisc9x+vKvrm1xP1XmLFxEYK+6Vx2hDcwhfya4OBfI4D9iu756E3HtM96jYw1Hz+uZH2KxvdbzLEd+a9YQp81qv+6zrfMS8vkWt411d2+3XN4/zZzThyUDdhSng3NNr7gC+8847uP/++/Gtb33LsXzp0qVYtmwZVqxYgXfeeQeVlZU44YQT0NDQkKWSEhEREWVXr2gA7t+/Hz/72c/wu9/9Dn37WnOvCSGwfPlyLFiwAKeddhrGjBmDhx9+GE1NTXj00UezWGIiIqL8octpYNx4kTt6RQr4oosuwsknn4wf/OAHuOmmm9TyLVu2oKamBpMnT1bLgsEgjjvuOLz55ps4//zzO7SfWAEgbOkwAPBEzdyCvG0fM9NEHkTlqMSoTNE2+gTglaO0EiaV1TxCpTVNQmhqu+Zdez1uGwXqk6PMYI6I06HL9ICckxWemAZz7l2PmQKM2yYSNtMkKbICqq+tZovZXCbX94Y1aHJUpJBpIuH1Iuw30iUtsoyaOVraI+AxvwO1/VZihzHC2BxlbKbDNA2AnKA1bo6082uIxcyUiPwZM+O1pT9tsauYE0ZJwvrarVF4XqjUrxrY2OyBMFP/Zkp0v0zZ+gWazJhVXevJ9W6LWwjnsSR0a4S1WuYVgJn2luk13Zz7tVBTx58zdrnMnBvWng5NFbuZMvJaP83vyhM2c5HyvagPQk6kHZEPnQ/7QtB85gHijNfjFSoPmyp2Vde6pupbrSPrXNd0RGXaVis06xwpYzd+alY6MFWd279ie33DOg48MQAtcp/m8R7WEJexN8sRmWad7/MJaB7zODe/A6hliYTuUXOyq7iFpob/iqDxuZhXgy5T/+bxrtmOdytmWe645qzv1mK3pQkT0/2eqGal/mV6VTTJOvf7EZYxN5h17hEpr29A68e74/omy2h2bdFDxpvCZ07+bLum2WM34zTPcfO70JH6OE84x4XHnvqVb4WtOhfyHI83ypHHvgCaZcx1tusbAHi8uuP61lrs9uNcyG4B5udEUEdczVAgY4+mPscBQIsZXUQcy1Kc4/HsDAJmCjgH9fgG4GOPPYZ//OMfeOedd5Leq6mpAQBUVFQ4lldUVOCLL75odZvhcBjhcFj9XV9f71JpiYiIiLKvR6eAt23bhksvvRR/+tOfEAqFWl1P0xLuJAiRtMxuyZIlKCsrU68hQ4a4VmYiIqJ8Y94BdONF7ujRdwA3bNiAnTt3Yty4cWpZPB7H66+/jhUrVmDz5s0AjDuBgwYNUuvs3Lkz6a6g3fz58zFv3jz1d319PYYMGWKlCFQ6xJYqkClXETZH7WmOFJrxUzh+BwDdfMam11pmpk/gESpF4xi6aaZTZCrQTIeKEJLTpXENSEx/2lIkSekS28hgR7rM3L2Z5THfiwEejzVBL2BMsmvFaaaF5d8+YRuFaKVIk2M3g7V/B/KnT4d68GfQSpfq5kg5M4Vm+6nSn7Z4E5+pa3+ucKpUoZaQcvFEoSaRFuY/KGwpRDOFZa9ze32by4zPWWkzq4uLLXZ7Ks1MtckLoW6rc1XumFXXKi2kym+vf7nMXueJqUJbulwObIU3Yo3CVd+B4zj3Ji0DgLhXqLSw9WxdKz6rILbvQa0ndx6wpQztdZ5Q76peY1pSFwBP3FnfatcJGVr7eWGm3MxzHLZzXE+IU9jqXzeP5ZSx24/3hPrXBDTzWbRqRK5m1XdCyhNxrZVzHEnLzO8n6Ti3nevm8eKNA7DVt1EMq84T61j4AOFJuL6Z54In4fom40zq1uUVjvq2x6vrmhoVq2KPackpYNvxntT9QU+REhcpUuhmnXs06/rmtb4Lc1Jt1UVGXdet3+117qxvc6ewvgdAfVfwwnF9A4CYrY6RWNe61Q3C0f3BngoHEM/SbR+mgHNPj24Afv/738eHH37oWPaLX/wChx56KK666ioccMABqKysxOrVqzF27FgAQCQSwZo1a3DLLbe0ut1gMIigfJICERERUW/ToxuAJSUlGDNmjGNZUVER+vfvr5bPnTsXixcvxsiRIzFy5EgsXrwYhYWFOPvss7NRZCIiorzDO4C5p0c3ANNx5ZVXorm5GRdeeKGaCHrVqlUoKSnJdtGIiIiIsqLXNQBfe+01x9+apmHhwoVYuHCha/vw2PsOpZheAHBOKWBNM6Cp2fXN2e3NpxAIn7D6tdnWSewXpnls/UjMXZr7tJXHPp2GmkU+Zj4lQFM/7VMImMsSnt/eer9A82dCHz2h2fsFOfsO6T7N6htnxu4VVh8hkdAX0AtomnM6DWhWzCYhEmIGrCelxDWrX6AtdrOMajoFzdl/xv57q9PGJMYu4/R4k/vGOWOXfYXMM9AHNQ2ImjJHQ/L0GVpCfdvKY58yB7a+cZr5JA1z+gjzaQGx5CklHDG31S/Q9v0nPmlAeDXrd3Wcm39bTzJRfZ18wroSqSmWrOmC1LFv32eq49y8M2BOqxGzfib2jRIxWx9QM3bdVp2ppkdKjF2z1bctduNnwrktvwvziSuaWiaPbQh1nTDrWrP1AU55jgvr+DbKbx3n9jr3pDjHgYQ6t/eNS3WOqxXlohTTxtiP96TrW9xW12b5zelmvFAHm3NaLDjYj3er3s1HwdhjT/gZB4RtWixjfVh9o23Heap+gSruFLF7UlzfAKN+Hdc3Gbt5LKhpkjwJ9W3uC0adC/t3j9bPcRWvre+vscx2fTNjz9INNKOLbdd3nvi/XOq8XtcAJCIiotzCFHDu6dHTwBARERFRx/EOYCeY/wDRUt2LTkyV2ta3Lzdv42u6/IBupc2s9INInYYxt6Ul/60yFuYTBIQ1fQESpiCwT2kDmb7RPEKlclQMqf7BJVL8bsbugTWdRmIZbU8fUc8kt62vvh8zHSasVJfWxs1/TbPFZxbHVh71fdueboKE1KVKwXXkWemJRbJNnaJitx0TKmbh/I6Fbu3fKrhteynSoCotaP5ti1+YH9CRXO+6lQoy00/mtB6abV9tssWdlDazx67q1SyPlep0pPQSZoGBZqv3FCealnhMAxBm6j4hVad5NNvxZO5Is85jexq/tdhTHe+w1V9CXdunzlF1rmvW9534vTj+KW7tINU5bn23CfUK64kd9joXifswrw0eLfk81dqp/hTXt6RpmmwpWvV921P26rzvWDLPcbyrfdq+A/v1BLb616Gub45liYHaz7fWzmukjl0dS6quE65vchvWuWLfeTvXNccCkbxPzYpNXd9s+048zrN1A413AHMPG4BERESUUWwA5h6mgImIiIjyDO8AdoAa8WZPebR29942QtCaPd+WbkzxjxhNpUnk37otTWWWAbb0iv2pCbCtADgetI6EUWNt5QCEZksVpPjngSpPYrraVg6RIr3aVtyODdvTiAnr68J8WLpIDkHYY3ZuU+hacprSLiGF5UiN29axp6WTtpVihGCq0dFJ8du2pVKG5gg+WKlcNTpaF61uQwjN+h5TPcUlVUwqlSpX8SKJ0FqJOXEbZpxeW3UmpuoSN2yW1YzdtqmkUc5tdEWwjwK2ulckHA8J5XakQs31EuJ31Hka3x9sf7d5oyIxNt2e4pRPlxDC6iLS3jkOONKOKbun2MptfizxPBZe6zLR1qh3e53br2/qZyuxa/Z0vxqlb7u+2eoz5fVNlkckXsvso8CTDj57GtRW523UZ7vXN8Co48R6b+sSK2CbmcBa3V7fjjK2VufqGmleG+S2Uh2j9jIldHXpbrwDmHt4B5CIiIgoz/AOIBEREWWUEJp117qL2yF3sAHYAerB57ZUV1spYPU5W5rA+r311I7KwdgmYFZpKttE0CJpiJhtG8KWHlApsdTrOMqhwZq4VA09TU6NqLBbSTkkpQXtE8cmFNueurBSI1Zazj6aEwBEqjRoirgcIy1VGjxxx7bv0ZMQd0L5VV20lUa2FyFplLGt3Inb0DX71y33o6kyqZF+WorYbdtKnCTWGHXdSl7KVkYzKB2atchW5+3Vt6M4HqTs/pBYBPsIWvWAe9t37KjvFPtxlEdoSSOsYUsBJ6ZE7fWj0r4arFS0GXuadZ6USk/RBSC5ALAmI4btf2y2OhetpUEBZ/pb/q0ldfdIUV51/RJqYl7zFLCPjlXdAdqLPUV3D5E4wlcdAMJ2Dmrqv+r7NUeje0Tq65tZHvv1DTBGWNtG4AK2utNSdGtpJbXb5vUtjXPckRJOOIagW/VpXt+EbR+OczzF/tQHUl3f5PZTdfcQCed4qm4e3UGH5spE0G5sgwxMARMREVGvdvfdd2PEiBEIhUIYN24c3njjjTbXX7NmDcaNG4dQKIQDDjgA9957b9I6TzzxBEaPHo1gMIjRo0dj5cqVjvdff/11nHrqqaiqqoKmaXjqqafcDKnL2AAkIiKijDIHgbjx6qjHH38cc+fOxYIFC7Bx40Ycc8wxOPHEE7F169aU62/ZsgUnnXQSjjnmGGzcuBG/+c1vMGfOHDzxxBNqnXXr1uGMM87AjBkz8P7772PGjBk4/fTT8dZbb6l1Ghsbcfjhh2PFihUd/8K6gSZEyuQG2dTX16OsrAwjFi6GJxSynptp/+bSGR2pwRp5mCplpEYlmsts6T576jgxnWYvSMLJodlSwIkjbDVbykCzpUiTnoOaYhRte6NqW0sBQxPJ6RJ7ykitL9R7SHzPHmKK2JOf56lZaSdbatSemnHGmzyC1LG+2p8zZsdbjnS/LaZUaeHE9e2xayliT5nPlKslxm6feDdh8mnN9p6jzhPWS5VCbXdkqDpuE49p2L4DKzaRMFm1Y5Ru4vW+tTpP7NpgL39SfVrPwXakDPXUx1Bro4ATipHQ3cMZk+McT1zfdmzY6zydUdSOsiad61ry9crRPUC+Z1+WeI7bt2dyHI+O4gCe1Oe48bft+uU4zq3POraVVBDnm46YUl3fZNm1hC4gjmce68nHUFsj/R3FSKhH5zEtV/S2cY5rSP/6Jt9MeX0zY0u4ljmOablMb2rB59csQF1dHUpLS5Fp5v8/v7PyUviKgl3eXqwxjLd/9P86VP4JEybg29/+Nu655x61bNSoUZg+fTqWLFmStP5VV12FZ555Bps2bVLLLrjgArz//vtYt24dAOCMM85AfX09XnzxRbXO1KlT0bdvX/zlL39J2qamaVi5ciWmT5+ebqgZxzuARERE1KPU19c7XuFwOOV6kUgEGzZswOTJkx3LJ0+ejDfffDPlZ9atW5e0/pQpU/Duu+8iGo22uU5r28xFbAASERFRRrmdAh4yZAjKysrUK9WdPADYvXs34vE4KioqHMsrKipQU1OT8jM1NTUp14/FYti9e3eb67S2zVzEUcAd4RFGesMcraX+0w5bGidlWjiRLU2l3rc9w1SNLlOfbaNPhCM1Zv5Msb5mbT9pAmhNJH027Y4DCSms1iaJVSP2zL/NSU4hrPImxW3fQQr2uFvPoCWNENRgG4lnG52YcmRga2yp2sQ6T10IYY2KVQeHsArnKGsaMdv/TkwVqjJa23ekqxNH3Qot/fo2P5dG7PavWNW3uR+PVe/Jh3kr8Scd5+ZPLXX92+pb7Tux3tM93tM8zhOHFzvq3H6Oy22mdY6rbSMp9pQTQtuuPdaoV1udu3GOpyobnJcSR53bh+yam0z3+mb+TPP6BsB5jbPXeUdiT5WqTXW8pzjH1X5sx3mXrm/yZ3v1DSBrt33cngZm27ZtjhRwMNh2ellLGFUuhEha1t76ics7us1cwwYgERER9SilpaVp9QEsLy+H1+tNujO3c+fOpDt4psrKypTr+3w+9O/fv811WttmLmIKmIiIiDJKuJT+7ehdxEAggHHjxmH16tWO5atXr8bEiRNTfqa6ujpp/VWrVmH8+PHw+/1trtPaNnMR7wASERFRrzVv3jzMmDED48ePR3V1Ne6//35s3boVF1xwAQBg/vz52L59Ox555BEAxojfFStWYN68eZg9ezbWrVuHBx54wDG699JLL8Wxxx6LW265BdOmTcPTTz+NV155BWvXrlXr7N+/H5999pn6e8uWLXjvvffQr18/DB06tJuibx0bgB1gPeHA6jPS5kPXzc85+ook/EyxotnbzN41RvVdsfWvSmtCdMdUDimmjUhk26bqIyVsfWISOse0F3+b/d/a/KD8mL0fpH1nXY3dzvqS1aqJmxdeq/NSV+q8tX+8avauQGY/Ic2F2FNNG+HYcfLHVL+wFHXe6dhTlc3+qITEboe67YPpHu9J/R/bq/OE/mC28zkx9rSPc7V9q7xtdUlz1Ln9HE/cabp13pG7I7aukY467+w5bttuq+UVmuP6Zq6u+r2aU/OkeuJNyu0lF6Qz1zegi7En/EzVB9R+jrd7fUsoa+qdt1YY5yLVeznx+OpmKb7eTm+no8444wzs2bMHN9xwA3bs2IExY8bghRdewLBhwwAAO3bscMwJOGLECLzwwgu47LLLcNddd6Gqqgp33HEHfvzjH6t1Jk6ciMceewxXX301rrnmGhx44IF4/PHHMWHCBLXOu+++i+OPP179PW/ePADAzJkz8dBDD3UiEndxHsA0mPMYDb9pETyhkPWGCw1AxwWkrc7UjguL+w1Aa5nt/wr2Bk8rF9VONQBTdUpurQO5vdHkZgOwtdgBYx7AVPG62QBM+ikyE3sbDUDN3mBI1ZG8re+qrV23d7wnLvMkx5ZyLrSuNAC7GLsrDcCOnOOJO+1kAzD13Ha2uJPeQ+rvqq3dphu7ek8kvZe8DXcbgBmPPY3jvEPnuP1zre48uTBtxi7/UaU3t+CL33T/PICH/+/l8BZ2fR7AeFMY7//k9m4rf2/GPoBEREREeYYp4A4wU8D2jGGH75+2kRZqe+e2j5tTCHRwG+ncvbGnxuy7VksSUhhpl8Geekn8l7Gw/lDpIdt7KdPgHfze20wLpShrUqrG9iV0ts7bXS3xO7XfBBGdi72tOnekh2yFVB/RbNvoQn2397mUqTH7QdfV4z3d78t28zvxeHe9zhO7e9jq1V7nid9Lu7tNuiuUtMt26zzxe+vwtSrVHb22OM5xq7AZub6plR27txa5cH1LYvty076+OT7Xzq5TZCnaLGNH69Mlbk8DQ13HBiARERFllC40aC403jrzLGBKjSlgIiIiojzDO4AdITv2ioRFHZGyw3QrK2lI0Wk4VefedAvRTjrQsc2E7SalJVKtn04ROpsOte07RZY6jR13YKcpcoGu1nlrK5lPh7DHnqpTeiZiVzto5eNu1Xcrx3tSakztuO2Ptr3zdt5PPIEEUta3uUpHYm+1zls5Z1OlwVMVMf0CpLlea/Vte69Ldd7O9S1p/Q6Oe3FuN4117CdXYjcFx3rO1dMuglvXN1sZ0t95mjvNWgq4E10pWtkOuYN3AImIiIjyDO8AEhERUUZxEEjuYQOwM1KMIMsIW2rMtijFei7vN2mWWrSauurUrtscMefMvThSY13ecQe0UcYu7TpXY0+VGktR553edVdSY13eeZo7bafSO/3/nTTrvNtjd8w+btuPG+d4B4uh9tVGFwDXd5oYewquH+/pnuNd2nluYgMw9zAFTERERJRneAeQiIiIMorTwOQeNgA7Qo4CVtobkZnqFn4Xj902UwZua2PEXGspYcd7jm11cN8pUmOJb2VMqi+5jTS4Q0djb2ebbaYFM6G11Fg31rfth+OtbtVGGtyhtTRdZ+LPRuypunuonSN5WaKu1H2KqQcydbynnOUgUYpR9732+pYFHAWce5gCJiIiIsozvANIREREGWXcAXRjEIgLhSEAvANIRERElHd6dANwyZIlOPLII1FSUoKBAwdi+vTp2Lx5s2MdIQQWLlyIqqoqFBQUYNKkSfj444/dKYBme7X3flvrtUVoSZ1CNOF8ZUWm4zal+Bdjt8SeagfpxNTR2FOt006dZ1xHY3e7vrN1vKfaeGfq28U+XDlT561+LsWro7rhHE+5re46x9uSretbFpjTwLjxInf06AbgmjVrcNFFF2H9+vVYvXo1YrEYJk+ejMbGRrXO0qVLsWzZMqxYsQLvvPMOKisrccIJJ6ChoSGLJSciIsofwsUXuaNH9wF86aWXHH8/+OCDGDhwIDZs2IBjjz0WQggsX74cCxYswGmnnQYAePjhh1FRUYFHH30U559/fjaKTURERJRVPfoOYKK6ujoAQL9+/QAAW7ZsQU1NDSZPnqzWCQaDOO644/Dmm292fAdt3ZfPUApISZEaSyxWV1IGbX6+rQ1nMmZThmNvU7p1nglpxNzVOm93hRyLHchwnWf6HG/rs9k8zu07Srkcma3zNGPPap13evvtfD6bsXcjpoBzT4++A2gnhMC8efNw9NFHY8yYMQCAmpoaAEBFRYVj3YqKCnzxxRetbiscDiMcDqu/6+vrM1BiIiKiPOFW/raHN4RzSZcbgHV1dVi5ciXeeOMNfP7552hqasKAAQMwduxYTJkyBRMnTnSjnO26+OKL8cEHH2Dt2rVJ72ma818MQoikZXZLlizB9ddf73oZiYiIiHJBp1PAO3bswOzZszFo0CDccMMNaGxsxBFHHIHvf//7GDx4MF599VWccMIJGD16NB5//HE3y5zkkksuwTPPPINXX30VgwcPVssrKysBWHcCTTt37ky6K2g3f/581NXVqde2bducK6RzP74zo0HT0c7t78SUQbvFTDe1kImY7Z9LRxupEqDj6ZIOpVUyUd/pxN1OzGbRspoq6midZzH2rB7vGTrOXT/eOxp7OrJc5+mVMUPXNxd19BzPqbSxW+lfpoBd0+k7gIcffjjOOeccvP322yrlmqi5uRlPPfUUli1bhm3btuGKK67odEFTEULgkksuwcqVK/Haa69hxIgRjvdHjBiByspKrF69GmPHjgUARCIRrFmzBrfcckur2w0GgwgGg66WlYiIKF/xUXC5p9MNwI8//hgDBgxoc52CggKcddZZOOuss7Br167O7qpVF110ER599FE8/fTTKCkpUXf6ysrKUFBQAE3TMHfuXCxevBgjR47EyJEjsXjxYhQWFuLss892vTxEREREPUGnG4DtNf66un467rnnHgDApEmTHMsffPBBzJo1CwBw5ZVXorm5GRdeeCFqa2sxYcIErFq1CiUlJV0vgCY6djvarTvXaT3Z3LZbN//FZG4s3bizeLfe9fRHNmPPdp1n4zg3ZTt2exnaXd/FfQutQ8FkLd2XqVHB2Yg92/VtL0N7u+5hd8LcGsHLUcDucWUU8DPPPJNyuaZpCIVCOOigg5LSs24QadwL1jQNCxcuxMKFC13fPxEREVFP5EoDcPr06dA0LalBZi7TNA1HH300nnrqKfTt29eNXRIREVFP4dYADt4BdI0rE0GvXr0aRx55JFavXq1Gzq5evRrf+c538Nxzz+H111/Hnj17XB8EkveyOSIqmzOT5vMFIB/rOxdk81hnnXe/bF/feuE1zhwE4saL3OHKHcBLL70U999/v2POv+9///sIhUL41a9+hY8//hjLly/Hueee68buiIiIiKgLXGkA/vvf/0ZpaWnS8tLSUvznP/8BAIwcORK7d+92Y3dERETUk/BJIDnHlRTwuHHj8N///d+OqV527dqFK6+8EkceeSQA4F//+pdjkuZeI5/TJNnE1Fh2MDWWf1jf2dHLYuezgHOPK3cAH3jgAUybNg2DBw/GkCFDoGkatm7digMOOABPP/00AGD//v245ppr3NgdEREREXWBKw3AQw45BJs2bcLLL7+Mf/7znxBC4NBDD8UJJ5wAj8e4yTh9+nQ3dkVEREQ9UZ4mTnKVKw1AwJjyZerUqZg0aRKCwSA0jbdpiYiIiHKRK30AdV3HjTfeiG984xsoLi7Gli1bAADXXHMNHnjgATd2QW3Jdn846l753jcqW9j/Mf/kc39fl7EPYO5xpQF400034aGHHsLSpUsRCATU8sMOOwy///3v3dgFERER9VTCxRe5wpUG4COPPIL7778fP/vZz+D1etXyb33rW/j000/d2AURERERucSVPoDbt2/HQQcdlLRc13VEo1E3dpH7OvoQcXJHBx+g7ipNsL6zQWhMy3W3bF7fsnmOk4s0+XJjO+QGV+4AfvOb38Qbb7yRtPyvf/0rxo4d68YuiIiIqKdiCjjnuHIH8LrrrsOMGTOwfft26LqOJ598Eps3b8YjjzyC5557zo1dEBEREZFLXLkDeOqpp+Lxxx/HCy+8AE3TcO2112LTpk149tlnccIJJ7ixC8plHCmXHdke/c06737scpAdHP3ddbwDmHNcmwdwypQpmDJlilubIyIiIqIMca0BSERERJSSW3cze8sd0RzQ6QZg3759037ax969ezu7G+oIjo7sfhz9nX9Y59nB0cA9mhDGy43tkDs63QBcvny5+n3Pnj246aabMGXKFFRXVwMA1q1bh5dffhnXXHNNlwtJRERERO7pdANw5syZ6vcf//jHuOGGG3DxxRerZXPmzMGKFSvwyiuv4LLLLutaKYmIiKjncmsAB+8AusaVUcAvv/wypk6dmrR8ypQpeOWVV9zYRc/B0ZFEmdebRkcS5QPznHXjRa5wpQHYv39/rFy5Mmn5U089hf79+7uxCyIiIiJyiSujgK+//nqcd955eO2111QfwPXr1+Oll17C73//ezd2QURERD2UW8kxJtjc40oDcNasWRg1ahTuuOMOPPnkkxBCYPTo0fj73/+OCRMmuLELorZxBDRR5nEENHUW+wDmHNfmAZwwYQL+/Oc/u7U5IiIiIsqQTvcBbGxszOj6RERE1EtwEEjO6XQD8KCDDsLixYvx1VdftbqOEAKrV6/GiSeeiDvuuKOzuyIiIiIiF3U6Bfzaa6/h6quvxvXXX48jjjgC48ePR1VVFUKhEGpra/HJJ59g3bp18Pv9mD9/Pn71q1+5WW4iygWa4L/IiTKtNzwFhX0Ac06nG4CHHHII/vrXv+LLL7/EX//6V7z++ut488030dzcjPLycowdOxa/+93vcNJJJ8HjcWW2GSIiIuqJ2ADMOV0eBDJ48GBcdtllfNoHERERUQ/h2ihgIk4RQdQNekM6kPIP7wDmHDYAiYiIKLPcGsHLGwyuYec8IiIiojzDO4BEPR1TgkSZxy4uXcJHweUe3gEkIiIiyjOuNQDfeOMN/PznP0d1dTW2b98OAPjjH/+ItWvXurWLLrn77rsxYsQIhEIhjBs3Dm+88Ua2i0RERJQfhIuvTuhoG2DNmjUYN24cQqEQDjjgANx7771J6zzxxBMYPXo0gsEgRo8ejZUrV3Z5v93JlQbgE088gSlTpqCgoAAbN25EOBwGADQ0NGDx4sVu7KJLHn/8ccydOxcLFizAxo0bccwxx+DEE0/E1q1bs100IiIiyqCOtgG2bNmCk046Cccccww2btyI3/zmN5gzZw6eeOIJtc66detwxhlnYMaMGXj//fcxY8YMnH766Xjrrbc6vd/upgkhupxRHzt2LC677DKcc845KCkpwfvvv48DDjgA7733HqZOnYqamho3ytppEyZMwLe//W3cc889atmoUaMwffp0LFmypN3P19fXo6ysDMNuuQmeUCi9nWarn0gudJDIx9hzoV9QtuLPZuz5XOf5Gnu2r3E9PHa9pQVfXHU16urqUFpa6kKh2mb+/3PoLTfBU5Dm/z/boDe3YGsHy9/RNsBVV12FZ555Bps2bVLLLrjgArz//vtYt24dAOCMM85AfX09XnzxRbXO1KlT0bdvX/zlL3/p1H67myt3ADdv3oxjjz02aXlpaSn27dvnxi46LRKJYMOGDZg8ebJj+eTJk/Hmm2+m/Ew4HEZ9fb3jRURERJ2jwRoI0qWX3F7i/6PNzGOizrQB1q1bl7T+lClT8O677yIajba5jrnNzuy3u7nSABw0aBA+++yzpOVr167FAQcc4MYuOm337t2Ix+OoqKhwLK+oqGj1zuSSJUtQVlamXkOGDOmOohIREfVO5jyAbrwADBkyxPH/6dbuqHWmDVBTU5Ny/Vgsht27d7e5jrnNzuy3u7kyDcz555+PSy+9FH/4wx+gaRq++uorrFu3DldccQWuvfZaN3bRZZrmvG0vhEhaZpo/fz7mzZun/q6vr2cjkIiIKEds27bNkQIOBoNtrt+RNkBr6ycuT2ebHd1vd3KlAXjllVeirq4Oxx9/PFpaWnDsscciGAziiiuuwMUXX+zGLjqtvLwcXq83qcW9c+fOpJa5KRgMtnswERERUZpcfhRcaWlpWn0AO9MGqKysTLm+z+dD//7921zH3GZn9tvdXJsGZtGiRdi9ezfefvttrF+/Hrt27cKNN97o1uY7LRAIYNy4cVi9erVj+erVqzFx4sQslYqIiIgyrTNtgOrq6qT1V61ahfHjx8Pv97e5jrnNntD2cPVJIIWFhRg/frybm3TFvHnzMGPGDIwfPx7V1dW4//77sXXrVlxwwQXZLpr7+FQIou6hieyPBCbqKVy+A9gR7bUB5s+fj+3bt+ORRx4BYIz4XbFiBebNm4fZs2dj3bp1eOCBB9ToXgC49NJLceyxx+KWW27BtGnT8PTTT+OVV15xzH2c622PTjcATzvttLTXffLJJzu7G1ecccYZ2LNnD2644Qbs2LEDY8aMwQsvvIBhw4ZltVxERET5IJuPgmuvDbBjxw7H3HwjRozACy+8gMsuuwx33XUXqqqqcMcdd+DHP/6xWmfixIl47LHHcPXVV+Oaa67BgQceiMcffxwTJkxIe7/Z1ul5AH/xi1+o34UQWLlyJcrKytQdwA0bNmDfvn047bTT8OCDD7pT2izpUfMAmvJxjrB8jNmO8wB2v3yNPV/jBnp87NmaB3D4okXp//+zDXpLCz5fsKDbyt+bdfoOoL1Rd9VVV+H000/HvffeC6/XCwCIx+O48MILWUFERET5LospYErNlUEgf/jDH3DFFVeoxh8AeL1ezJs3D3/4wx/c2AURERH1VFl+FjAlc6UBGIvFHI9MMW3atAm6rruxCyIiIiJyiSujgH/xi1/g3HPPxWeffYbvfve7AID169fj5ptvdvQVJCIiovyTzUEglJorDcDbbrsNlZWV+O1vf4sdO3YAMB4Pd+WVV+Lyyy93YxdERETUU9ke49bl7ZArXGkAejweXHnllbjyyitRX18PABz8QURERJSjXJ0IGmDDj4iIiBJwFHDOcaUBOGLEiDYfbvyf//zHjd0QERERkQtcaQDOnTvX8Xc0GsXGjRvx0ksv4b//+7/d2AURERH1UBwEkntcaQBeeumlKZffddddePfdd93YBREREfVUTAHnHFfmAWzNiSeeiCeeeCKTuyAiIiKiDnJ9EIjd//7v/6Jfv36Z3AURERHlOpdSwLwD6B5XGoBjx451DAIRQqCmpga7du3C3Xff7cYuiIiIqKdiCjjnuNIAnDZtmqMB6PF4MGDAAEyaNAmHHnqoG7sgIiIiIpe40gBcuHChG5shIiKi3oh3AHOOK4NAvF4vdu7cmbR8z5498Hq9buyCiIiIiFziyh1AIVI3ycPhMAKBgBu7ICIioh6K8wDmni41AO+44w4AgKZp+P3vf4/i4mL1Xjwex+uvv84+gEREREQ5pksNwN/+9rcAjDuA9957ryPdGwgEMHz4cNx7771dKyERERERuapLDcAtW7YAAI4//ng8+eST6Nu3ryuFoi7g/XGi7iFaf/45ESXgIJCc40ofwFdffdWNzRAREVEvxD6AuafTDcB58+bhxhtvRFFREebNm9fmusuWLevsboiIiIjIZZ1uAG7cuBHRaBQA8I9//MMxETQRERGRA+/e5ZRONwDtad/XXnvNjbIQERERUTdwZSLoc889Fw0NDUnLGxsbce6557qxCyIiIuqphIsvcoUrDcCHH34Yzc3NScubm5vxyCOPuLELIiIi6qHMQSBuvMgdXRoFXF9fDyEEhBBoaGhAKBRS78XjcbzwwgsYOHBglwtJRERERO7pUgOwT58+0DQNmqbh4IMPTnpf0zRcf/31XdkFERER9XScBzDndKkB+Oqrr0IIge9973t44okn0K9fP/VeIBDAsGHDUFVV1eVCEhERUc/FeQBzT5cagMcddxwA44kgQ4YMgcfjSpdCIiIiIsogV54EMmzYMABAU1MTtm7dikgk4nj/W9/6lhu7ISIiop6IKeCc40oDcNeuXfjFL36BF198MeX78Xjcjd0QERFRT8QGYM5xJWc7d+5c1NbWYv369SgoKMBLL72Ehx9+GCNHjsQzzzzjxi6IqDWcG4Eo84RmvIh6CVfuAP7tb3/D008/jSOPPBIejwfDhg3DCSecgNLSUixZsgQnn3yyG7shIiKiHoiDQHKPK3cAGxsb1Xx//fr1w65duwAAhx12GP7xj3+4sYskn3/+Oc477zyMGDECBQUFOPDAA3Hdddcl9T/cunUrTj31VBQVFaG8vBxz5sxJWoeIiIgon7hyB/CQQw7B5s2bMXz4cBxxxBG47777MHz4cNx7770YNGiQG7tI8umnn0LXddx333046KCD8NFHH2H27NlobGzEbbfdBsDoe3jyySdjwIABWLt2Lfbs2YOZM2dCCIE777wzI+UiIiKiBOwDmHNcaQDOnTsXO3bsAABcd911mDJlCv785z8jEAjgoYcecmMXSaZOnYqpU6eqvw844ABs3rwZ99xzj2oArlq1Cp988gm2bdum5iO8/fbbMWvWLCxatAilpaUZKRsRERHZsAGYc1xpAP7sZz9Tv48dOxaff/45Pv30UwwdOhTl5eVu7CItdXV1jsmo161bhzFjxjgmo54yZQrC4TA2bNiA448/vtvKRkRERJQrXGkAJiosLMS3v/3tTGy6Vf/+979x55134vbbb1fLampqUFFR4Vivb9++CAQCqKmpaXVb4XAY4XBY/V1fX+9+gYmIiPIEB4Hknk43AOfNm5f2usuWLUt73YULF7b7/OB33nkH48ePV39/9dVXmDp1Kn7605/il7/8pWNdTUseti+ESLnctGTJEj7DmIiIyC1MAeecTjcAN27cmNZ6bTW0Urn44otx5plntrnO8OHD1e9fffUVjj/+eFRXV+P+++93rFdZWYm33nrLsay2thbRaDTpzqDd/PnzHQ3c+vp6DBkypANREBEREeWuTjcAX331VTfLoZSXl6fdb3D79u04/vjjMW7cODz44INJzyKurq7GokWLsGPHDjUaedWqVQgGgxg3blyr2w0GgwgGg50PgoiIiBSmgHNPRvoAdoevvvoKkyZNwtChQ3HbbbepuQcB484fAEyePBmjR4/GjBkzcOutt2Lv3r244oorMHv2bI4AJiIiorzVYxuAq1atwmeffYbPPvsMgwcPdrwnhPFPBK/Xi+effx4XXnghjjrqKBQUFODss89W08QQERFRN2AfwJzTYxuAs2bNwqxZs9pdb+jQoXjuuecyXyDiczKJugNzYNQTsQGYc1x5FBwRERER9Rw99g4gERER9QyafLmxHXIHG4BERESUWUwB5xymgImIiIjyDO8AEhERUUZxHsDcwzuARERERHmGdwCJiIgos9gHMOewAUhERESZx8ZbTmEKmIiIiCjP8A4gEXUen/5ClHm9YOQDB4HkHjYAiYiIKLPYBzDnMAVMRERElGd4B5CIiIgyiing3MM7gERERER5hncAiYiIKLPYBzDn8A4gERERZZSZAnbjlSm1tbWYMWMGysrKUFZWhhkzZmDfvn1tfkYIgYULF6KqqgoFBQWYNGkSPv74Y8c64XAYl1xyCcrLy1FUVIQf/vCH+PLLLx3rLFq0CBMnTkRhYSH69OnjcmSpsQFIREREee/ss8/Ge++9h5deegkvvfQS3nvvPcyYMaPNzyxduhTLli3DihUr8M4776CyshInnHACGhoa1Dpz587FypUr8dhjj2Ht2rXYv38/TjnlFMTjcbVOJBLBT3/6U/zXf/1XxuJLxBQwERERZVaOp4A3bdqEl156CevXr8eECRMAAL/73e9QXV2NzZs345BDDkkuihBYvnw5FixYgNNOOw0A8PDDD6OiogKPPvoozj//fNTV1eGBBx7AH//4R/zgBz8AAPzpT3/CkCFD8Morr2DKlCkAgOuvvx4A8NBDD2UmwBR4B5B6Bw4NI8o8oXHyb+oc4eILQH19veMVDoe7VLx169ahrKxMNf4A4Lvf/S7Kysrw5ptvpvzMli1bUFNTg8mTJ6tlwWAQxx13nPrMhg0bEI1GHetUVVVhzJgxrW63u7ABSERERD3KkCFDVF+9srIyLFmypEvbq6mpwcCBA5OWDxw4EDU1Na1+BgAqKiocyysqKtR7NTU1CAQC6Nu3b6vrZAtTwERERJRRbs8DuG3bNpSWlqrlwWAw5foLFy5U6dXWvPPOO8a2teS720KIlMsdZUp4P53PpLNOprEBSERERJnlch/A0tJSRwOwNRdffDHOPPPMNtcZPnw4PvjgA3z99ddJ7+3atSvpDp+psrISgHGXb9CgQWr5zp071WcqKysRiURQW1vruAu4c+dOTJw4sd3yZxJTwERERNQrlZeX49BDD23zFQqFUF1djbq6Orz99tvqs2+99Rbq6upabaiNGDEClZWVWL16tVoWiUSwZs0a9Zlx48bB7/c71tmxYwc++uijrDcAeQeQiIiIMkoTApro+i1AN7aRyqhRozB16lTMnj0b9913HwDgV7/6FU455RTHCOBDDz0US5YswY9+9CNomoa5c+di8eLFGDlyJEaOHInFixejsLAQZ599NgCgrKwM5513Hi6//HL0798f/fr1wxVXXIHDDjtMjQoGgK1bt2Lv3r3YunUr4vE43nvvPQDAQQcdhOLi4ozEzAYgERER5b0///nPmDNnjhqx+8Mf/hArVqxwrLN582bU1dWpv6+88ko0NzfjwgsvRG1tLSZMmIBVq1ahpKRErfPb3/4WPp8Pp59+Opqbm/H9738fDz30ELxer1rn2muvxcMPP6z+Hjt2LADg1VdfxaRJkzIRLjQhMtSc7kXq6+tRVlaGYbfcBE8olN6HsjVVQjanQ8nm9BDZngYmX2PP9pQg+Rh7PsZsl634e8k5rre04IurrkZdXV1afei6yvz/5xE/XwRvIM3/f7YhHmnBe39a0G3l7814B5CIiIgyyu1RwNR1HARCRERElGd4B9BtuZAiIerteBuAqGfJ8UfB5SM2AImIiCijmALOPUwBExEREeUZ3gEkIiKizGIKOOfwDiARERFRnuEdQCIiIsoo9gHMPWwAEhERUWYxBZxzekUKOBwO44gjjoCmaer5eaatW7fi1FNPRVFREcrLyzFnzhxEIpHsFJSIiIgoB/SKO4BXXnklqqqq8P777zuWx+NxnHzyyRgwYADWrl2LPXv2YObMmRBC4M4778xSaYmIiPIP07e5pcc3AF988UWsWrUKTzzxBF588UXHe6tWrcInn3yCbdu2oaqqCgBw++23Y9asWVi0aBGfI0hdx4m/iboHWw89mxDGy43tkCt6dAr466+/xuzZs/HHP/4RhYWFSe+vW7cOY8aMUY0/AJgyZQrC4TA2bNjQ6nbD4TDq6+sdLyIiIqLeosc2AIUQmDVrFi644AKMHz8+5To1NTWoqKhwLOvbty8CgQBqampa3faSJUtQVlamXkOGDHG17ERERPnEHAXsxovckXMNwIULF0LTtDZf7777Lu68807U19dj/vz5bW5P05JTdEKIlMtN8+fPR11dnXpt27aty3ERERER5Yqc6wN48cUX48wzz2xzneHDh+Omm27C+vXrEQwGHe+NHz8eP/vZz/Dwww+jsrISb731luP92tpaRKPRpDuDdsFgMGm7RERE1EmcBibn5FwDsLy8HOXl5e2ud8cdd+Cmm25Sf3/11VeYMmUKHn/8cUyYMAEAUF1djUWLFmHHjh0YNGgQAGNgSDAYxLhx4zITABERETlouvFyYzvkjpxrAKZr6NChjr+Li4sBAAceeCAGDx4MAJg8eTJGjx6NGTNm4NZbb8XevXtxxRVXYPbs2RwBTERERHkr5/oAusnr9eL5559HKBTCUUcdhdNPPx3Tp0/Hbbfdlu2iERER5Q/h4otc0WPvACYaPnw4RIr5gYYOHYrnnnsuCyUiIiIigM8CzkW9+g4gERERESXrNXcACfynEVF34NNfiDqOTwLJOWwAEhERUUYxBZx7mAImIiIiyjO8A0hERESZxYmgcw7vABIRERHlGd4BJCIiooxiH8DcwwYgUU/FKyFR9+DI767jKOCcwxQwERERUZ7hHUAiIiLKKKaAcw8bgEQ9lZmW4hWRKLM0wTRwV3EUcM5hCpiIiIgoz/AOIBEREWUUU8C5h3cAiYiIiPIM7wAS9VT8pzBR92D/v67ThfFyYzvkCjYAiYiIKLM4CCTnMAVMRERElGd4B5CIiIgySoNLg0C6vgmS2AAkIiKizOKj4HIOU8BEREREeYZ3AKnrOEKOqHtw5Df1UJwHMPfwDiARERFRnuEdQCIiIsosTgOTc9gApK4z78kzFUyUWUJjDox6JE0IaC4M4HBjG2RgCpiIiIgoz/AOIBEREWWWLl9ubIdcwQYgdV02U79Mh2UH0/3Zka3jnfWdHb3o+sYUcO5hCpiIiIgoz/AOIBEREWUWRwHnHDYAqeuyOQrY3GcvSpWkLZuxc+R3dmRrFHC+17cmeH3rKj4KLucwBUxERESUZ3gHkIiIiDKKj4LLPbwDSERERJRneAfQLfnaNwbI32lgsl3n+Rp7Pt8CyMdpYHKhvrMVfy7E7hb2Acw5Pf4O4PPPP48JEyagoKAA5eXlOO200xzvb926FaeeeiqKiopQXl6OOXPmIBKJZKm0RERE+UfT3XuRO3r0HcAnnngCs2fPxuLFi/G9730PQgh8+OGH6v14PI6TTz4ZAwYMwNq1a7Fnzx7MnDkTQgjceeedWSw5ERERUfb02AZgLBbDpZdeiltvvRXnnXeeWn7IIYeo31etWoVPPvkE27ZtQ1VVFQDg9ttvx6xZs7Bo0SKUlpa6V6BcmCYhH6eIyPZUKNmub7Mc3S3bdZ7N1Fg+xp7tmO1lyAZOA9N1TAHnnB6bAv7HP/6B7du3w+PxYOzYsRg0aBBOPPFEfPzxx2qddevWYcyYMarxBwBTpkxBOBzGhg0bslFsIiKi/CNcfJEremwD8D//+Q8AYOHChbj66qvx3HPPoW/fvjjuuOOwd+9eAEBNTQ0qKiocn+vbty8CgQBqampa3XY4HEZ9fb3jRURERNRb5FwDcOHChdA0rc3Xu+++C103eoIuWLAAP/7xjzFu3Dg8+OCD0DQNf/3rX9X2NC35tr0QIuVy05IlS1BWVqZeQ4YMab/gQsvuSDG3JlnqKDPubMeeDdmOORux2+s730ZGZjP2bJ7jQPbrO9txZyP2bNe5yzQhXHuRO3KuAXjxxRdj06ZNbb7GjBmDQYMGAQBGjx6tPhsMBnHAAQdg69atAIDKysqkO321tbWIRqNJdwbt5s+fj7q6OvXatm1bBiIlIiKiXFFbW4sZM2aomz8zZszAvn372vyMEAILFy5EVVUVCgoKMGnSJEdXNMDIKl5yySUoLy9HUVERfvjDH+LLL79U73/++ec477zzMGLECBQUFODAAw/Eddddl/EZS3JuEEh5eTnKy8vbXW/cuHEIBoPYvHkzjj76aABANBrF559/jmHDhgEAqqursWjRIuzYsUM1GFetWoVgMIhx48a1uu1gMIhgMOhCNERERNQTBoGcffbZ+PLLL/HSSy8BAH71q19hxowZePbZZ1v9zNKlS7Fs2TI89NBDOPjgg3HTTTfhhBNOwObNm1FSUgIAmDt3Lp599lk89thj6N+/Py6//HKccsop2LBhA7xeLz799FPouo777rsPBx10ED766CPMnj0bjY2NuO222zIWryZEz72fOnfuXPzv//4v/vCHP2DYsGG49dZb8eyzz+LTTz9F3759EY/HccQRR6CiogK33nor9u7di1mzZmH69Okdmgamvr4eZWVlGHbLTfCEQul9KN/SY3b5GHs+T5Sbr7FzIvDsyLcuH3YuxK63tOCLq65GXV2duzNhtML8/+fx354PnzfN/3+2IRZvwav/WOJ6+Tdt2oTRo0dj/fr1mDBhAgBg/fr1qK6uxqeffuqYYcQkhEBVVRXmzp2Lq666CoBxt6+iogK33HILzj//fNTV1WHAgAH44x//iDPOOAMA8NVXX2HIkCF44YUXMGXKlJTlufXWW3HPPfeo8Q6ZkHMp4I649dZbceaZZ2LGjBk48sgj8cUXX+Bvf/sb+vbtCwDwer14/vnnEQqFcNRRR+H000/H9OnTM9qiJiIiop5l3bp1KCsrU40/APjud7+LsrIyvPnmmyk/s2XLFtTU1GDy5MlqWTAYxHHHHac+s2HDBkSjUcc6VVVVGDNmTKvbBYC6ujr069evq2G1KedSwB3h9/tx2223tdmgGzp0KJ577rluLBURERHZuTWAw9xG4uwcXe26VVNTg4EDByYtHzhwYKuzhpjLE8cUVFRU4IsvvlDrBAIBdWPKvk5r2/33v/+NO++8E7fffnuH4+iIHn0HMOsSRwV294ixxNGg3ZkayWbsqeLujtizXd9A9uo8F+s728d7d8hm7Plc57l2jveG0cACVj/ALr2MzQ0ZMsQxW8eSJUtS7jbdmUWAzs0akupz6XymtXW++uorTJ06FT/96U/xy1/+ss1tdFWPvgNIRERE+Wfbtm2OPoCt3f27+OKLceaZZ7a5reHDh+ODDz7A119/nfTerl27Wp01pLKyEoBxl88caAoAO3fuVJ+prKxEJBJBbW2t4y7gzp07MXHiRMf2vvrqKxx//PGorq7G/fff32aZ3cAGIBEREWWWy6OAS0tL0xoEku7MItXV1airq8Pbb7+N73znOwCAt956C3V1dUkNNdOIESNQWVmJ1atXY+zYsQCASCSCNWvW4JZbbgFgzFji9/uxevVqnH766QCAHTt24KOPPsLSpUvVtrZv347jjz9ezWns8WQ+QcsUMBEREeW1UaNGYerUqZg9ezbWr1+P9evXY/bs2TjllFMcI4APPfRQrFy5EoCR+p07dy4WL16MlStX4qOPPsKsWbNQWFiIs88+GwBQVlaG8847D5dffjn+7//+Dxs3bsTPf/5zHHbYYfjBD34AwLjzN2nSJAwZMgS33XYbdu3ahZqamjafWOYG3gHsivb6ZGS630g6289Uv5Fsxt7etrMVM5D9Ou+NsWfzOE93+70x9lyu82zXN9B7Y88UHYAbX5nuwjZa8ec//xlz5sxRI3Z/+MMfYsWKFY51Nm/ejLq6OvX3lVdeiebmZlx44YWora3FhAkTsGrVKjUHIAD89re/hc/nw+mnn47m5mZ8//vfx0MPPQSv1wvAmJ/4s88+w2effYbBgwc79pfJmfp69DyA3aVT8wACvWbeqE7Jx3nCgOzXeb7Gnu3/KeZr7Pk456eph8aerXkAvz/mSvi8XX/AQiwexv99tLTbyt+bMQVMRERElGeYAu6IbD0EPtvy8e5Gtu/kAfkZez7ezTLla+zZvsblc+zdqQc8Ci7fsAFIREREmcUGYM5hCpiIiIgoz/AOYEd0dTb2zqQa3ExPdLbsbqQpOhuHG/F3pvzZjNmtzwM9L/ZsxdyVz6WSrXOd53j3ftbUU2PvTrwDmHPYACQiIqLM6gHTwOQbpoCJiIiI8gzvAHZEpkcB5+soQCB/Y8/XuIH8jD0X0nX5eLznY8ymxNiz9F1oQkBzIX3rxjbIwDuARERERHmGdwCJiIgoszgIJOewAdgRXR0FbJfqNnx33ZpPFUN3pSlaizFfY++OuFuLrzvibiu+3l7nuXaOp/OeW7J5nufjOQ6kf55nKyWtu/T/T50NQLcwBUxERESUZ3gHkIiIiDKLKeCcwwZgR2RiFHCujRDrTvk4CtSUr7Hn4yhQUz7Gnu3rG5CfI75zoetF8o5darzlwDHVSzAFTERERJRneAeQiIiIMosp4JzDO4BEREREeYZ3ADsiE9PAZLt/Tjb6qWQzdjPebMad+Ht3yYXYsxV3LvTHyqfYsx1zW393h0w/Nao16Zzj2ZwGxo3+e5wGxjVsABIREVFmCd14ubEdcgVTwERERER5hncAO6K3TQOTj9MjmPIxdk6Fkh35XOf5Gnu2r3G5OA0MB4HkHDYAiYiIKLPYBzDnMAVMRERElGd4B7AjeuMo4GzIhZjzOfZsyPao0FyIPRuyfZznY+xA9ruY5OIoYKaAcw7vABIRERHlGd4BJCIioswScOkOYNc3QQY2ADuCo4DdkY8xm/I19lweFZlp+Vjn2a5vID/rPJ2YOQqYJKaAiYiIiPJMj24A/vOf/8S0adNQXl6O0tJSHHXUUXj11Vcd62zduhWnnnoqioqKUF5ejjlz5iASiWSpxERERHlI1917kSt6dAr45JNPxsEHH4y//e1vKCgowPLly3HKKafg3//+NyorKxGPx3HyySdjwIABWLt2Lfbs2YOZM2dCCIE777yz4zvkKGB35GPMpnyNPdtx53Ps2ZDtrhZAftZ5OjFzFDBJPfYO4O7du/HZZ5/h17/+Nb71rW9h5MiRuPnmm9HU1ISPP/4YALBq1Sp88skn+NOf/oSxY8fiBz/4AW6//Xb87ne/Q319fZYjICIiIsqOHtsA7N+/P0aNGoVHHnkEjY2NiMViuO+++1BRUYFx48YBANatW4cxY8agqqpKfW7KlCkIh8PYsGFDtopORESUX8w7gG68yBU9NgWsaRpWr16NadOmoaSkBB6PBxUVFXjppZfQp08fAEBNTQ0qKiocn+vbty8CgQBqampa3XY4HEY4HFZ/824hERFRF/BRcDkn5xqACxcuxPXXX9/mOu+88w7GjRuHCy+8EAMHDsQbb7yBgoIC/P73v8cpp5yCd955B4MGDQJgNBQTCSFSLjctWbIkdRncnAYm233RgNzup5Ip2e6blK+x52vcQH7Gns/XNyC3+15m+3ygnJFzDcCLL74YZ555ZpvrDB8+HH/729/w3HPPoba2FqWlpQCAu+++G6tXr8bDDz+MX//616isrMRbb73l+GxtbS2i0WjSnUG7+fPnY968eerv+vp6DBkypAtRERER5S8hdAjR9RG8bmyDDDnXACwvL0d5eXm76zU1NQEAPB5nN0aPxwNdDhOvrq7GokWLsGPHDnVHcNWqVQgGg6qfYCrBYBDBYLCzIRARERHltJxrAKaruroaffv2xcyZM3HttdeioKAAv/vd77BlyxacfPLJAIDJkydj9OjRmDFjBm699Vbs3bsXV1xxBWbPnq3uGnZIJqaByYZ0HhieSbkQezZk4kkyHZHt2LMl28d5PseeDdm+vgG5nXrP5jQwbvTf4yAQ1/TYUcDl5eV46aWXsH//fnzve9/D+PHjsXbtWjz99NM4/PDDAQBerxfPP/88QqEQjjrqKJx++umYPn06brvttiyXnoiIKI9wFHDO6bF3AAFg/PjxePnll9tcZ+jQoXjuuee6qUREREREua9HNwC7nRspvGymJ7KZkrLHm6+xZzPuxN8zLTHWfKzzbNd3qr+7Q7Zjz/b1LdXfmdTRczxb6WldBzQXBnBwEIhr2AAkIiKizBIuzQPIFLBremwfQCIiIiLqHN4B7G6Jt9+7M1WQuK/uTAWk2lc+xt7dqTH7/rJd30D26ry7017ZqvNU+8rX2LN5vOdbnadB6DqECylgzgPoHt4BJCIiIsozvANIREREmcU+gDmHDcCOyMTotlyYMLW7cYLc7MiFiXFzeYLcTMl2nedz7NmS7eO8rdiz9b3oLh2LbAC6hilgIiIiojzDO4BERESUWUIAcGMeQN4BdAsbgERERJRRQhcQLqSABRuArmEKmIiIiCjP8A4gERERZZbQ4U4KmPMAuoV3AImIiIjyDO8AEhERUUaxD2DuYQOQiIiIMosp4JzDBmAazH9x6C0tGdh4Hk6UmgvPqczH2PNxYlxTvsaezbiB/I092xNBt8H8/1h330mLIerKg0BiiHZ9IwQA0ATvp7bryy+/xJAhQ7JdDCIiIlds27YNgwcPzvh+WlpaMGLECNTU1Li2zcrKSmzZsgWhUMi1beYjNgDToOs6Nm/ejNGjR2Pbtm0oLS3NdpEyrr6+HkOGDMmbeAHGnA8x51u8AGPOh5g7Eq8QAg0NDaiqqoLH0z3jQFtaWhCJRFzbXiAQYOPPBUwBp8Hj8eAb3/gGAKC0tDQvLiimfIsXYMz5IN/iBRhzPkg33rKysm4ojSUUCrHBloM4DQwRERFRnmEDkIiIiCjPsAGYpmAwiOuuuw7BYDDbRekW+RYvwJjzQb7FCzDmfJBv8ZI7OAiEiIiIKM/wDiARERFRnmEDkIiIiCjPsAFIRERElGfYAEzD3XffjREjRiAUCmHcuHF44403sl0kVyxcuBCapjlelZWV6n0hBBYuXIiqqioUFBRg0qRJ+Pjjj7NY4o57/fXXceqpp6KqqgqapuGpp55yvJ9OjOFwGJdccgnKy8tRVFSEH/7wh/jyyy+7MYqOaS/mWbNmJdX7d7/7Xcc6PSnmJUuW4Mgjj0RJSQkGDhyI6dOnY/PmzY51els9pxNzb6rne+65B9/61rfUPHfV1dV48cUX1fu9rX6B9mPuTfVL2cEGYDsef/xxzJ07FwsWLMDGjRtxzDHH4MQTT8TWrVuzXTRXfPOb38SOHTvU68MPP1TvLV26FMuWLcOKFSvwzjvvoLKyEieccAIaGhqyWOKOaWxsxOGHH44VK1akfD+dGOfOnYuVK1fisccew9q1a7F//36ccsopiMfj3RVGh7QXMwBMnTrVUe8vvPCC4/2eFPOaNWtw0UUXYf369Vi9ejVisRgmT56MxsZGtU5vq+d0YgZ6Tz0PHjwYN998M9599128++67+N73vodp06apRl5vq1+g/ZiB3lO/lCWC2vSd73xHXHDBBY5lhx56qPj1r3+dpRK557rrrhOHH354yvd0XReVlZXi5ptvVstaWlpEWVmZuPfee7uphO4CIFauXKn+TifGffv2Cb/fLx577DG1zvbt24XH4xEvvfRSt5W9sxJjFkKImTNnimnTprX6mZ4e886dOwUAsWbNGiFEftRzYsxC9P567tu3r/j973+fF/VrMmMWovfXL2Ue7wC2IRKJYMOGDZg8ebJj+eTJk/Hmm29mqVTu+te//oWqqiqMGDECZ555Jv7zn/8AALZs2YKamhpH7MFgEMcdd1yviT2dGDds2IBoNOpYp6qqCmPGjOnR38Nrr72GgQMH4uCDD8bs2bOxc+dO9V5Pj7murg4A0K9fPwD5Uc+JMZt6Yz3H43E89thjaGxsRHV1dV7Ub2LMpt5Yv9R9+CzgNuzevRvxeBwVFRWO5RUVFaipqclSqdwzYcIEPPLIIzj44IPx9ddf46abbsLEiRPx8ccfq/hSxf7FF19ko7iuSyfGmpoaBAIB9O3bN2mdnnoMnHjiifjpT3+KYcOGYcuWLbjmmmvwve99Dxs2bEAwGOzRMQshMG/ePBx99NEYM2YMgN5fz6liBnpfPX/44Yeorq5GS0sLiouLsXLlSowePVo1Znpj/bYWM9D76pe6HxuAadA0zfG3ECJpWU904oknqt8PO+wwVFdX48ADD8TDDz+sOhP31tjtOhNjT/4ezjjjDPX7mDFjMH78eAwbNgzPP/88TjvttFY/1xNivvjii/HBBx9g7dq1Se/11npuLebeVs+HHHII3nvvPezbtw9PPPEEZs6ciTVr1qj3e2P9thbz6NGje139UvdjCrgN5eXl8Hq9Sf9a2rlzZ9K/NnuDoqIiHHbYYfjXv/6lRgP35tjTibGyshKRSAS1tbWtrtPTDRo0CMOGDcO//vUvAD035ksuuQTPPPMMXn31VQwePFgt78313FrMqfT0eg4EAjjooIMwfvx4LFmyBIcffjj+3//7f726fluLOZWeXr/U/dgAbEMgEMC4ceOwevVqx/LVq1dj4sSJWSpV5oTDYWzatAmDBg3CiBEjUFlZ6Yg9EolgzZo1vSb2dGIcN24c/H6/Y50dO3bgo48+6jXfw549e7Bt2zYMGjQIQM+LWQiBiy++GE8++ST+9re/YcSIEY73e2M9txdzKj29nhMJIRAOh3tl/bbGjDmV3la/1A26fdhJD/PYY48Jv98vHnjgAfHJJ5+IuXPniqKiIvH5559nu2hddvnll4vXXntN/Oc//xHr168Xp5xyiigpKVGx3XzzzaKsrEw8+eST4sMPPxRnnXWWGDRokKivr89yydPX0NAgNm7cKDZu3CgAiGXLlomNGzeKL774QgiRXowXXHCBGDx4sHjllVfEP/7xD/G9731PHH744SIWi2UrrDa1FXNDQ4O4/PLLxZtvvim2bNkiXn31VVFdXS2+8Y1v9NiY/+u//kuUlZWJ1157TezYsUO9mpqa1Dq9rZ7bi7m31fP8+fPF66+/LrZs2SI++OAD8Zvf/EZ4PB6xatUqIUTvq18h2o65t9UvZQcbgGm46667xLBhw0QgEBDf/va3HVMt9GRnnHGGGDRokPD7/aKqqkqcdtpp4uOPP1bv67ourrvuOlFZWSmCwaA49thjxYcffpjFEnfcq6++KgAkvWbOnCmESC/G5uZmcfHFF4t+/fqJgoICccopp4itW7dmIZr0tBVzU1OTmDx5shgwYIDw+/1i6NChYubMmUnx9KSYU8UKQDz44INqnd5Wz+3F3Nvq+dxzz1XX4AEDBojvf//7qvEnRO+rXyHajrm31S9lhyaEEN13v5GIiIiIso19AImIiIjyDBuARERERHmGDUAiIiKiPMMGIBEREVGeYQOQiIiIKM+wAUhERESUZ9gAJCIiIsozbAASERER5Rk2AIl6sEmTJmHu3Lm9ar+zZs3C9OnTu7ydzZs3o7KyEg0NDa2u89BDD6FPnz5d3pfdzp07MWDAAGzfvt3V7RIRuYkNQCLqsCeffBI33nij+nv48OFYvnx59gqUwoIFC3DRRRehpKSkW/c7cOBAzJgxA9ddd1237peIqCPYACSiDuvXr1+3N6w64ssvv8QzzzyDX/ziF1nZ/y9+8Qv8+c9/Rm1tbVb2T0TUHjYAiXqR2tpanHPOOejbty8KCwtx4okn4l//+pd630x5vvzyyxg1ahSKi4sxdepU7NixQ60Ti8UwZ84c9OnTB/3798dVV12FmTNnOtKy9hTwpEmT8MUXX+Cyyy6DpmnQNA0AsHDhQhxxxBGO8i1fvhzDhw9Xf8fjccybN0/t68orr0Ti48mFEFi6dCkOOOAAFBQU4PDDD8f//u//tvk9/M///A8OP/xwDB482LH8oYcewtChQ1FYWIgf/ehH2LNnT9Jnn332WYwbNw6hUAgHHHAArr/+esRiMfX+p59+iqOPPhqhUAijR4/GK6+8Ak3T8NRTT6l1DjvsMFRWVmLlypVtlpOIKFvYACTqRWbNmoV3330XzzzzDNatWwchBE466SREo1G1TlNTE2677Tb88Y9/xOuvv46tW7fiiiuuUO/fcsst+POf/4wHH3wQf//731FfX+9o3CR68sknMXjwYNxwww3YsWOHozHZnttvvx1/+MMf8MADD2Dt2rXYu3dvUqPp6quvxoMPPoh77rkHH3/8MS677DL8/Oc/x5o1a1rd7uuvv47x48c7lr311ls499xzceGFF+K9997D8ccfj5tuusmxzssvv4yf//znmDNnDj755BPcd999eOihh7Bo0SIAgK7rmD59OgoLC/HWW2/h/vvvx4IFC1KW4Tvf+Q7eeOONtL8LIqJuJYioxzruuOPEpZdeKoQQ4p///KcAIP7+97+r93fv3i0KCgrE//zP/wghhHjwwQcFAPHZZ5+pde666y5RUVGh/q6oqBC33nqr+jsWi4mhQ4eKadOmpdyvEEIMGzZM/Pa3v3WU7brrrhOHH364Y9lvf/tbMWzYMPX3oEGDxM0336z+jkajYvDgwWpf+/fvF6FQSLz55puO7Zx33nnirLPOavV7Ofzww8UNN9zgWHbWWWeJqVOnOpadccYZoqysTP19zDHHiMWLFzvW+eMf/ygGDRokhBDixRdfFD6fT+zYsUO9v3r1agFArFy50vG5yy67TEyaNKnVMhIRZZMvy+1PInLJpk2b4PP5MGHCBLWsf//+OOSQQ7Bp0ya1rLCwEAceeKD6e9CgQdi5cycAoK6uDl9//TW+853vqPe9Xi/GjRsHXdddLW9dXR127NiB6upqtczn82H8+PEqDfzJJ5+gpaUFJ5xwguOzkUgEY8eObXXbzc3NCIVCjmWbNm3Cj370I8ey6upqvPTSS+rvDRs24J133lF3/AAjTd3S0oKmpiZs3rwZQ4YMQWVlpXrf/l3ZFRQUoKmpqdUyEhFlExuARL2ESOg7Z19u9ssDAL/f73hf07Skz9rXb2vbbfF4PEmfs6ei02E2Op9//nl84xvfcLwXDAZb/Vx5eXnSAIx0YtB1Hddffz1OO+20pPdCoVDSd9mWvXv3YsCAAWmtS0TU3dgHkKiXGD16NGKxGN566y21bM+ePfjnP/+JUaNGpbWNsrIyVFRU4O2331bL4vE4Nm7c2ObnAoEA4vG4Y9mAAQNQU1PjaHi99957jn0NGjQI69evV8tisRg2bNjgiCkYDGLr1q046KCDHK8hQ4a0Wp6xY8fik08+cSwbPXq0Y18Akv7+9re/jc2bNyft66CDDoLH48Ghhx6KrVu34uuvv1afeeedd1KW4aOPPmrzLiURUTbxDiBRLzFy5EhMmzYNs2fPxn333YeSkhL8+te/xje+8Q1MmzYt7e1ccsklWLJkCQ466CAceuihuPPOO1FbW9vmna/hw4fj9ddfx5lnnolgMIjy8nJMmjQJu3btwtKlS/GTn/wEL730El588UWUlpaqz1166aW4+eabMXLkSIwaNQrLli3Dvn371PslJSW44oorcNlll0HXdRx99NGor6/Hm2++ieLiYsycOTNleaZMmYJf/vKXiMfj8Hq9AIA5c+Zg4sSJWLp0KaZPn45Vq1Y50r8AcO211+KUU07BkCFD8NOf/hQejwcffPABPvzwQ9x000044YQTcOCBB2LmzJlYunQpGhoa1CAQ+/fT1NSEDRs2YPHixWl/70RE3Yl3AIl6kQcffBDjxo3DKaecgurqaggh8MILLySlfdty1VVX4ayzzsI555yD6upqFBcXY8qUKUl96uxuuOEGfP755zjwwANV2nPUqFG4++67cdddd+Hwww/H22+/7RhtDACXX345zjnnHMyaNQvV1dUoKSlJ6qd344034tprr8WSJUswatQoTJkyBc8++yxGjBjRanlOOukk+P1+vPLKK2rZd7/7Xfz+97/HnXfeiSOOOAKrVq3C1Vdf7fjclClT8Nxzz2H16tU48sgj8d3vfhfLli3DsGHDABj9IZ966ins378fRx55JH75y1+qbdi/n6effhpDhw7FMccc09ZXTUSUNZroTOceIsobuq5j1KhROP300x1P/8h1d999N55++mm8/PLLGd3P3//+dxx99NH47LPP1OCa73znO5g7dy7OPvvsjO6biKizmAImIocvvvgCq1atwnHHHYdwOIwVK1Zgy5YtPa4x86tf/Qq1tbVoaGhw9aklK1euRHFxMUaOHInPPvsMl156KY466ijV+Nu5cyd+8pOf4KyzznJtn0REbuMdQCJy2LZtG84880x89NFHEEJgzJgxuPnmm3Hsscdmu2g54ZFHHsGNN96Ibdu2oby8HD/4wQ9w++23o3///tkuGhFR2tgAJCIiIsozHARCRERElGfYACQiIiLKM2wAEhEREeUZNgCJiIiI8gwbgERERER5hg1AIiIiojzDBiARERFRnmEDkIiIiCjPsAFIRERElGfYACQiIiLKM2wAEhEREeUZNgCJiIiI8gwbgERERER5hg1AIiIiojzz/wESjUvUtx10wQAAAABJRU5ErkJggg==", "text/html": [ "\n", "
\n", "
\n", " Figure\n", "
\n", " \n", "
\n", " " ], "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "phi, theta = dist.local_grids(full_basis)\n", "latitude = (np.pi / 2 - theta) / np.pi * 180\n", "longitude = phi / np.pi * 180\n", "plt.figure()\n", "plt.pcolormesh(longitude.ravel(), latitude.ravel(), h1['g'].T)\n", "plt.xlabel('longitude (deg)')\n", "plt.ylabel('latitude (deg)')\n", "plt.colorbar()\n", "plt.tight_layout()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Evolving the most unstable mode (IVP)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we can run an nonlinear IVP, starting with the background jet plus a small amount of the most unstable mode, to see how it saturates.\n", "First we make new fields for the total variables and set their initial conditions using the background and fasting growing mode." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Full fields\n", "u = dist.VectorField(coords, name='u', bases=full_basis)\n", "h = dist.Field(name='h', bases=full_basis)\n", "u0.change_scales(1)\n", "u1.change_scales(1)\n", "h0.change_scales(1)\n", "h1.change_scales(1)\n", "u['g'] = u0['g'] + 1e-3*u1['g']\n", "h['g'] = h0['g'] + 1e-3*h1['g']" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we can run the full forward simulation using our original nonlinear equation set:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2024-09-27 10:19:06,357 subsystems 0/1 INFO :: Building subproblem matrices 1/127 (~1%) Elapsed: 0s, Remaining: 3s, Rate: 3.9e+01/s\n", "2024-09-27 10:19:06,521 subsystems 0/1 INFO :: Building subproblem matrices 13/127 (~10%) Elapsed: 0s, Remaining: 2s, Rate: 6.9e+01/s\n", "2024-09-27 10:19:06,693 subsystems 0/1 INFO :: Building subproblem matrices 26/127 (~20%) Elapsed: 0s, Remaining: 1s, Rate: 7.2e+01/s\n", "2024-09-27 10:19:06,878 subsystems 0/1 INFO :: Building subproblem matrices 39/127 (~31%) Elapsed: 1s, Remaining: 1s, Rate: 7.1e+01/s\n", "2024-09-27 10:19:07,081 subsystems 0/1 INFO :: Building subproblem matrices 52/127 (~41%) Elapsed: 1s, Remaining: 1s, Rate: 6.9e+01/s\n", "2024-09-27 10:19:07,281 subsystems 0/1 INFO :: Building subproblem matrices 65/127 (~51%) Elapsed: 1s, Remaining: 1s, Rate: 6.8e+01/s\n", "2024-09-27 10:19:07,480 subsystems 0/1 INFO :: Building subproblem matrices 78/127 (~61%) Elapsed: 1s, Remaining: 1s, Rate: 6.8e+01/s\n", "2024-09-27 10:19:07,678 subsystems 0/1 INFO :: Building subproblem matrices 91/127 (~72%) Elapsed: 1s, Remaining: 1s, Rate: 6.8e+01/s\n", "2024-09-27 10:19:07,877 subsystems 0/1 INFO :: Building subproblem matrices 104/127 (~82%) Elapsed: 2s, Remaining: 0s, Rate: 6.7e+01/s\n", "2024-09-27 10:19:08,070 subsystems 0/1 INFO :: Building subproblem matrices 117/127 (~92%) Elapsed: 2s, Remaining: 0s, Rate: 6.7e+01/s\n", "2024-09-27 10:19:08,225 subsystems 0/1 INFO :: Building subproblem matrices 127/127 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 6.7e+01/s\n", "2024-09-27 10:19:08,263 __main__ 0/1 INFO :: Starting main loop\n", "2024-09-27 10:19:25,957 __main__ 0/1 INFO :: Iteration=1, Time=1.666667e-01, dt=1.666667e-01\n", "2024-09-27 10:19:34,956 __main__ 0/1 INFO :: Iteration=101, Time=1.683333e+01, dt=1.666667e-01\n", "2024-09-27 10:19:44,291 __main__ 0/1 INFO :: Iteration=201, Time=3.350000e+01, dt=1.666667e-01\n", "2024-09-27 10:19:53,342 __main__ 0/1 INFO :: Iteration=301, Time=5.016667e+01, dt=1.666667e-01\n", "2024-09-27 10:20:02,421 __main__ 0/1 INFO :: Iteration=401, Time=6.683333e+01, dt=1.666667e-01\n", "2024-09-27 10:20:12,287 __main__ 0/1 INFO :: Iteration=501, Time=8.350000e+01, dt=1.666667e-01\n", "2024-09-27 10:20:22,301 __main__ 0/1 INFO :: Iteration=601, Time=1.001667e+02, dt=1.666667e-01\n", "2024-09-27 10:20:31,953 __main__ 0/1 INFO :: Iteration=701, Time=1.168333e+02, dt=1.666667e-01\n", "2024-09-27 10:20:42,143 __main__ 0/1 INFO :: Iteration=801, Time=1.335000e+02, dt=1.666667e-01\n", "2024-09-27 10:20:51,137 __main__ 0/1 INFO :: Iteration=901, Time=1.501667e+02, dt=1.666667e-01\n", "2024-09-27 10:21:01,459 __main__ 0/1 INFO :: Iteration=1001, Time=1.668333e+02, dt=1.666667e-01\n", "2024-09-27 10:21:11,282 __main__ 0/1 INFO :: Iteration=1101, Time=1.835000e+02, dt=1.666667e-01\n", "2024-09-27 10:21:20,457 __main__ 0/1 INFO :: Iteration=1201, Time=2.001667e+02, dt=1.666667e-01\n", "2024-09-27 10:21:29,778 __main__ 0/1 INFO :: Iteration=1301, Time=2.168333e+02, dt=1.666667e-01\n", "2024-09-27 10:21:38,953 __main__ 0/1 INFO :: Iteration=1401, Time=2.335000e+02, dt=1.666667e-01\n", "2024-09-27 10:21:42,588 solvers 0/1 INFO :: Simulation stop time reached.\n", "2024-09-27 10:21:42,589 solvers 0/1 INFO :: Final iteration: 1441\n", "2024-09-27 10:21:42,590 solvers 0/1 INFO :: Final sim time: 240.16666666666174\n", "2024-09-27 10:21:42,593 solvers 0/1 INFO :: Setup time (init - iter 0): 16.02 sec\n", "2024-09-27 10:21:42,597 solvers 0/1 INFO :: Warmup time (iter 0-10): 4.5 sec\n", "2024-09-27 10:21:42,598 solvers 0/1 INFO :: Run time (iter 10-end): 135.8 sec\n", "2024-09-27 10:21:42,598 solvers 0/1 INFO :: CPU time (iter 10-end): 0.03772 cpu-hr\n", "2024-09-27 10:21:42,599 solvers 0/1 INFO :: Speed: 1.028e+06 mode-stages/cpu-sec\n" ] } ], "source": [ "# Timestepping parameters\n", "timestep = 600 * second\n", "stop_sim_time = 240 * hour\n", "\n", "# Problem\n", "problem = d3.IVP([u, h], namespace=locals())\n", "problem.add_equation(\"dt(u) + nu*lap(lap(u)) + g*grad(h) + 2*Omega*zcross(u) = - u@grad(u)\")\n", "problem.add_equation(\"dt(h) + nu*lap(lap(h)) + H*div(u) = - div(h*u)\")\n", "\n", "# Solver\n", "solver = problem.build_solver(d3.RK222)\n", "solver.stop_sim_time = stop_sim_time\n", "\n", "# Analysis\n", "snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=1*hour)\n", "snapshots.add_task(h, name='height')\n", "snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity')\n", "\n", "# Main loop\n", "try:\n", " logger.info('Starting main loop')\n", " while solver.proceed:\n", " solver.step(timestep)\n", " if (solver.iteration-1) % 100 == 0:\n", " logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))\n", "except:\n", " logger.error('Exception raised, triggering end of main loop.')\n", " raise\n", "finally:\n", " solver.log_stats()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let's make a movie of the solution using some plotting tools from another script:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2024-09-27 10:21:44,506 matplotlib.animation 0/1 INFO :: Animation.save using \n", "2024-09-27 10:21:44,508 matplotlib.animation 0/1 INFO :: MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 600x600 -pix_fmt rgba -framerate 20.0 -i pipe: -vcodec h264 -pix_fmt yuv420p -y /var/folders/gl/8q1_pm2s1490lvyfvm_8yby80000gn/T/tmp8tzvqhz0/temp.m4v\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import matplotlib\n", "from matplotlib import animation\n", "from IPython.display import HTML\n", "import h5py\n", "\n", "# Plot parameters\n", "task = 'vorticity'\n", "cmap = plt.cm.RdBu_r\n", "dpi = 100\n", "figsize = (6, 6)\n", "\n", "def build_s2_coord_vertices(phi, theta):\n", " phi = phi.ravel()\n", " phi_vert = np.concatenate([phi, [2*np.pi]])\n", " phi_vert -= phi_vert[1] / 2\n", " theta = theta.ravel()\n", " theta_mid = (theta[:-1] + theta[1:]) / 2\n", " theta_vert = np.concatenate([[np.pi], theta_mid, [0]])\n", " return np.meshgrid(phi_vert, theta_vert, indexing='ij')\n", "\n", "# Create figure\n", "with h5py.File('snapshots/snapshots_s1.h5', mode='r') as file:\n", " fig = plt.figure(figsize=figsize, dpi=dpi)\n", " ax = fig.add_axes([0, 0, 1, 1], projection='3d')\n", " # Plot writes\n", " dset = file['tasks'][task]\n", " phi = dset.dims[1][0][:].ravel()\n", " theta = dset.dims[2][0][:].ravel()\n", " phi_vert, theta_vert = build_s2_coord_vertices(phi, theta)\n", " x = np.sin(theta_vert) * np.cos(phi_vert)\n", " y = np.sin(theta_vert) * np.sin(phi_vert)\n", " z = np.cos(theta_vert)\n", " data = dset[0]\n", " clim = np.max(np.abs(dset[:]))\n", " norm = matplotlib.colors.Normalize(-clim, clim)\n", " fc = cmap(norm(data))\n", " surf = ax.plot_surface(x, y, z, facecolors=fc, cstride=1, rstride=1, linewidth=0, antialiased=False, shade=False, zorder=5)\n", " ax.set_box_aspect((1,1,1))\n", " ax.set_xlim(-0.7, 0.7)\n", " ax.set_ylim(-0.7, 0.7)\n", " ax.set_zlim(-0.7, 0.7)\n", " ax.axis('off')\n", "\n", " def animate(i):\n", " data = dset[i]\n", " fc = cmap(norm(data))\n", " surf.set_facecolors(fc.reshape(fc.size//4, 4))\n", " return surf\n", "\n", " anim = animation.FuncAnimation(fig, animate, frames=dset.shape[0], interval=50, blit=True)\n", " video = HTML(anim.to_html5_video())\n", " plt.close(fig)\n", "\n", "video.data = video.data.replace('autoplay', '')\n", "video" ] } ], "metadata": { "kernelspec": { "display_name": "d3", "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.11.6" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }