{ "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 numpy setuptools wheel\n", " !CC=mpicc pip3 install --no-cache --no-build-isolation dedalus\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", " # Step 5: Install ipympl and restart kernel\n", " !pip3 install -q ipympl\n", " get_ipython().kernel.do_shutdown(restart=True)\n", " else:\n", " print(\"See website for installation instructions:\")\n", " print(\"https://dedalus-project.readthedocs.io/en/latest/pages/installation.html\")" ] }, { "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": [ { "name": "stdout", "output_type": "stream", "text": [ "env: OMP_NUM_THREADS=1\n", "env: NUMEXPR_MAX_THREADS=1\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", "logger = logging.getLogger(__name__)\n", "\n", "# Check if running on google colab\n", "import os\n", "using_google_colab = bool(os.getenv(\"COLAB_RELEASE_TAG\"))\n", "\n", "# Setup interactive matplotlib\n", "if using_google_colab:\n", " from google.colab import output\n", " output.enable_custom_widget_manager()\n", "\n", "# Import libraries\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import dedalus.public as d3\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", "R = 1 # Earth radius\n", "hour = 1\n", "\n", "# Physical parameters\n", "meter = R / 6.37122e6\n", "second = hour / 3600\n", "Omega = 7.292e-5 / second\n", "g = 9.80616 * meter / second**2\n", "H = 1e4 * meter\n", "\n", "# Numerical parameters\n", "Nphi = 256\n", "Ntheta = 128\n", "dealias = 3/2\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": [ "2025-10-07 22:22:22,265 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 8.2e+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": "f8befa0f73f84fd6889421d90408cde7", "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 dozen 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:201: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil and dok are more efficient.\n", " self._set_arrayXarray_sparse(i, j, x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2025-10-07 22:22:26,043 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 4.6e-01/s\n", "2025-10-07 22:22:31,499 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.5e-01/s\n", "2025-10-07 22:22:37,100 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.6e-01/s\n", "2025-10-07 22:22:42,834 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.6e-01/s\n", "2025-10-07 22:22:48,222 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.6e-01/s\n", "2025-10-07 22:22:53,695 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.6e-01/s\n", "2025-10-07 22:22:59,062 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.4e-01/s\n", "2025-10-07 22:23:04,285 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.3e-01/s\n", "2025-10-07 22:23:09,877 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.5e-01/s\n", "2025-10-07 22:23:15,203 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.3e-01/s\n", "2025-10-07 22:23:20,020 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.5e-01/s\n", "2025-10-07 22:23:24,998 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.5e-01/s\n", "2025-10-07 22:23:29,644 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 5.4e-01/s\n" ] } ], "source": [ "# Eigenvalue solver\n", "m_max = 12\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": "c2c3800d449446f186454343822cb69c", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8KklEQVR4nO3df3RU9Z3/8ddk8mMsJiMJ5pdGGMAupNFqJgWDRu1Wwg8XpWJFW8Rqqw2rhSSl3whoI3g01bbatUo4/Mi2LFvFr4CFbaRJtSCWdKMkoQVStSUSxJlvGqgT1M0PJvf7B5up40wgQCYzk/t8nDPnOJ987r3vO3JyX/l87nyuxTAMQwAAADCNmHAXAAAAgKFFAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGRiw11ANOvt7dUHH3ygxMREWSyWcJcDAAAGwDAMHT9+XJmZmYqJMedYGAHwHHzwwQfKysoKdxkAAOAsHD58WBdffHG4ywgLAuA5SExMlHTyH1BSUlKYqwEAAAPR0dGhrKws33XcjAiA56Bv2jcpKYkACABAlDHz7VvmnPgGAAAwMQIgAACAyRAAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmw0LQABAlvL2G6luOqe14p1ITbZrkSJY1xrwL2QI4ewRAAIgC2/e5tHzbAbk8nb62DLtN5bOyNT0nI4yVAYhGTAEDQITbvs+lBRsa/MKfJLk9nVqwoUHb97nCVBmAaEUABIDP8PYaqvvrUf2q6Yjq/npU3l4jrLUs33ZAwSroa1u+7UBYawQQfZgCBoBPibSp1vqWYwEjf59mSHJ5OlXfckz541KGrjAAUY0RQAD4X5E41dp2vP/wdzb9AEAiAAKApMidak1NtA1qPwCQCIAAIOnMplqH0iRHsjLsNvW32ItFJ6eoJzmSh7IsAFGOAAgAitypVmuMReWzsiUpIAT2vS+flc16gADOCAEQABTZU63TczJUOS9X6Xb/Y6fbbaqcl8s6gADOWNQEwJUrV8rhcMhms8npdGrXrl2n7L9z5045nU7ZbDaNHTtWq1at8vv59ddfL4vFEvC68cYbQ3kaACJUpE+1Ts/J0Btl/6zn771K/3b7FXr+3qv0Rtk/E/4AnJWoCIAbN25UcXGxli1bpsbGRhUUFGjGjBlqbW0N2r+lpUUzZ85UQUGBGhsbtXTpUi1cuFCbNm3y9dm8ebNcLpfvtW/fPlmtVn3ta18bqtMCEEGiYarVGmNR/rgU3XzFRcofl8K0L4CzZjEMI+JXD508ebJyc3NVWVnpa5s4caJmz56tioqKgP5lZWXaunWrmpubfW1FRUXau3ev6urqgh7jpz/9qX7wgx/I5XJpxIgRA6qro6NDdrtdHo9HSUlJZ3hWACJRpK0DCGDwcf2OgoWgu7u7tWfPHj344IN+7YWFhdq9e3fQberq6lRYWOjXNm3aNK1bt049PT2Ki4sL2GbdunW6/fbbBxz+AAxP03MyNDU7XfUtx9R2vFOpiSenfRltAzCcRHwAbG9vl9frVVpaml97Wlqa3G530G3cbnfQ/idOnFB7e7syMvz/iq+vr9e+ffu0bt26U9bS1dWlrq4u3/uOjo4zORUAUaJvqhUAhquouAdQkiwW/7++DcMIaDtd/2Dt0snRv5ycHE2aNOmUNVRUVMhut/teWVlZAy0fAAAgYkR8ABw1apSsVmvAaF9bW1vAKF+f9PT0oP1jY2OVkuL/V/0nn3yiF154Qd/+9rdPW8uSJUvk8Xh8r8OHD5/h2QAAAIRfxAfA+Ph4OZ1O1dbW+rXX1tZqypQpQbfJz88P6F9TU6O8vLyA+/9efPFFdXV1ad68eaetJSEhQUlJSX4vAACAaBPxAVCSSktLtXbtWlVVVam5uVklJSVqbW1VUVGRpJMjc/Pnz/f1Lyoq0qFDh1RaWqrm5mZVVVVp3bp1Wrx4ccC+161bp9mzZweMDAIAAAxXEf8lEEmaO3eujh49qhUrVsjlciknJ0fV1dUaPXq0JMnlcvmtCehwOFRdXa2SkhI999xzyszM1DPPPKM5c+b47fedd97RG2+8oZqamiE9HwAAgHCKinUAIxXrCAEAEH24fkfJFDAAAAAGDwEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEwmKtYBBABEJm+vofqWY2o73qnURJsmOZJljen/Oe0AIgMBEABwVrbvc2n5tgNyeTp9bRl2m8pnZWt6TkYYKwNwOkwBAwDO2PZ9Li3Y0OAX/iTJ7enUgg0N2r7PFabKAAwEARAAcEa8vYaWbzugYI+R6mtbvu2AvL08aAqIVARAAMAZqW85FjDy92mGJJenU/Utx4auKABnhAAIADgjbcf7D39n0w/A0CMAAgDOSGqibVD7ARh6BEAAwBmZ5EhWht2m/hZ7sejkt4EnOZKHsiwAZ4AACAA4I9YYi8pnZUtSQAjse18+K5v1AIEIRgAEAJyx6TkZqpyXq3S7/zRvut2mynm5rAMIRDgWggYQNjxFIrpNz8nQ1Ox0/h8CUYgACCAseIrE8GCNsSh/XEq4ywBwhpgCBjDkeIoEAIQXARDAkOIpEgAQfgRAAEOKp0gAQPgRAAEMKZ4iAQDhRwAEMKR4igQAhB8BEMCQ4ikSABB+BEAAQ4qnSABA+BEAAQw5niIBAOHFQtAAwoKnSABA+BAAAYQNT5EAgPBgChgAAMBkCIAAAAAmQwAEAAAwGQIgAACAyRAAAQAATIYACAAAYDIEQAAAAJMhAAIAAJhM1ATAlStXyuFwyGazyel0ateuXafsv3PnTjmdTtlsNo0dO1arVq0K6PPhhx/q/vvvV0ZGhmw2myZOnKjq6upQnQIAAEBEiIoAuHHjRhUXF2vZsmVqbGxUQUGBZsyYodbW1qD9W1paNHPmTBUUFKixsVFLly7VwoULtWnTJl+f7u5uTZ06Ve+9955eeuklvf3221qzZo0uuuiioTotAACAsLAYhmGEu4jTmTx5snJzc1VZWelrmzhxombPnq2KioqA/mVlZdq6dauam5t9bUVFRdq7d6/q6uokSatWrdKPfvQj/fnPf1ZcXNxZ1dXR0SG73S6Px6OkpKSz2gcAABhaXL+jYASwu7tbe/bsUWFhoV97YWGhdu/eHXSburq6gP7Tpk3TW2+9pZ6eHknS1q1blZ+fr/vvv19paWnKycnR448/Lq/X228tXV1d6ujo8HsBAABEm4gPgO3t7fJ6vUpLS/NrT0tLk9vtDrqN2+0O2v/EiRNqb2+XJB08eFAvvfSSvF6vqqur9dBDD+knP/mJHnvssX5rqaiokN1u972ysrLO8ewAAACGXsQHwD4Wi8XvvWEYAW2n6//p9t7eXqWmpmr16tVyOp26/fbbtWzZMr9p5s9asmSJPB6P73X48OGzPR0AAICwiQ13AaczatQoWa3WgNG+tra2gFG+Punp6UH7x8bGKiUlRZKUkZGhuLg4Wa1WX5+JEyfK7Xaru7tb8fHxAftNSEhQQkLCuZ4SAABAWEX8CGB8fLycTqdqa2v92mtrazVlypSg2+Tn5wf0r6mpUV5enu8LH1dffbX+8pe/qLe319fnnXfeUUZGRtDwBwAAMFxEfACUpNLSUq1du1ZVVVVqbm5WSUmJWltbVVRUJOnk1Oz8+fN9/YuKinTo0CGVlpaqublZVVVVWrdunRYvXuzrs2DBAh09elSLFi3SO++8o1//+td6/PHHdf/99w/5+QEAAAyliJ8ClqS5c+fq6NGjWrFihVwul3JyclRdXa3Ro0dLklwul9+agA6HQ9XV1SopKdFzzz2nzMxMPfPMM5ozZ46vT1ZWlmpqalRSUqLLL79cF110kRYtWqSysrIhPz8AAIChFBXrAEYq1hECACD6cP2OkilgAAAADB4CIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMJnYcBcAAMBg8/Yaqm85prbjnUpNtGmSI1nWGEu4ywIiBgEQADCsbN/n0vJtB+TydPraMuw2lc/K1vScjDBWBkQOpoABAMPG9n0uLdjQ4Bf+JMnt6dSCDQ3avs8VpsqAyEIABAAMC95eQ8u3HZAR5Gd9bcu3HZC3N1gPwFwIgACAYaG+5VjAyN+nGZJcnk7VtxwbuqKACEUABAAMC23H+w9/Z9MPGM4IgACAYSE10Tao/YDhjAAIABgWJjmSlWG3qb/FXiw6+W3gSY7koSwLiEgEQADAsGCNsah8VrYkBYTAvvfls7JZDxAQARAAMIxMz8lQ5bxcpdv9p3nT7TZVzstlHUDgf7EQNABgWJmek6Gp2ek8CQQ4BQIgAGDYscZYlD8uJdxlABGLKWAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMlETQBcuXKlHA6HbDabnE6ndu3adcr+O3fulNPplM1m09ixY7Vq1Sq/n//85z+XxWIJeHV28pBwAAAwvEVFANy4caOKi4u1bNkyNTY2qqCgQDNmzFBra2vQ/i0tLZo5c6YKCgrU2NiopUuXauHChdq0aZNfv6SkJLlcLr+XzcZDwgEAwPBmMQzDCHcRpzN58mTl5uaqsrLS1zZx4kTNnj1bFRUVAf3Lysq0detWNTc3+9qKioq0d+9e1dXVSTo5AlhcXKwPP/zwrOvq6OiQ3W6Xx+NRUlLSWe8HAAAMHa7fUTAC2N3drT179qiwsNCvvbCwULt37w66TV1dXUD/adOm6a233lJPT4+v7aOPPtLo0aN18cUX61/+5V/U2Ng4+CcARAhvr6G6vx7Vr5qOqO6vR+Xtjfi//QAAIRLxj4Jrb2+X1+tVWlqaX3taWprcbnfQbdxud9D+J06cUHt7uzIyMjRhwgT9/Oc/12WXXaaOjg7927/9m66++mrt3btXl156adD9dnV1qaury/e+o6PjHM8OGBrb97m0fNsBuTz/uMc1w25T+axsTc/JCGNlAIBwiPgRwD4Wi/9DvA3DCGg7Xf9Pt1911VWaN2+evvjFL6qgoEAvvviiPv/5z+tnP/tZv/usqKiQ3W73vbKyss72dIAhs32fSws2NPiFP0lyezq1YEODtu9zhakyAEC4RHwAHDVqlKxWa8BoX1tbW8AoX5/09PSg/WNjY5WSEvzh4DExMfrSl76kd999t99alixZIo/H43sdPnz4DM8GGFreXkPLtx1QsMnevrbl2w4wHQwAJhPxATA+Pl5Op1O1tbV+7bW1tZoyZUrQbfLz8wP619TUKC8vT3FxcUG3MQxDTU1NysjofzosISFBSUlJfi8gktW3HAsY+fs0Q5LL06n6lmNDVxQAIOwiPgBKUmlpqdauXauqqio1NzerpKREra2tKioqknRyZG7+/Pm+/kVFRTp06JBKS0vV3NysqqoqrVu3TosXL/b1Wb58uX7zm9/o4MGDampq0re+9S01NTX59gkMB23HB7au5UD7AQCGh4j/EogkzZ07V0ePHtWKFSvkcrmUk5Oj6upqjR49WpLkcrn81gR0OByqrq5WSUmJnnvuOWVmZuqZZ57RnDlzfH0+/PBD3XfffXK73bLb7bryyiv1+uuva9KkSUN+fkCopCYObF3LgfYDAAwPUbEOYKRiHSFEOm+voWueeE1uT2fQ+wAtktLtNr1R9s+yxvT/pSoAGE64fkfJFDCAs2ONsah8Vrakk2Hv0/rel8/KJvwBgMkQAIFhbnpOhirn5Srd7j/Nm263qXJeLusAAoAJRcU9gADOzfScDE3NTld9yzG1He9UaqJNkxzJjPwBgEkRAAGTsMZYlD8u+DqYAABzYQoYAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYTEgD4Icffqi1a9dqyZIlOnbsmCSpoaFBR44cCeVhAQAAcAqxodrxH//4R91www2y2+167733dO+99yo5OVlbtmzRoUOHtH79+lAdGgAAAKcQshHA0tJSffOb39S7774rm83ma58xY4Zef/31UB0WAAAApxGyAPjmm2/qO9/5TkD7RRddJLfbHarDAgAA4DRCFgBtNps6OjoC2t9++21deOGFoTosAAAATiNkAfDmm2/WihUr1NPTI0myWCxqbW3Vgw8+qDlz5oTqsAAAADiNkAXAH//4x/rb3/6m1NRU/c///I+uu+46jR8/XomJiXrsscdCdVgAAACcRsi+BZyUlKQ33nhDr732mhoaGtTb26vc3FzdcMMNoTokAAAABsBiGIYRih2vX79ec+fOVUJCgl97d3e3XnjhBc2fPz8Uhx1SHR0dstvt8ng8SkpKCnc5AABgALh+hzAAWq1WuVwupaam+rUfPXpUqamp8nq9oTjskOIfEAAA0YfrdwjvATQMQxaLJaD9/fffl91uD9VhAQAAcBqDfg/glVdeKYvFIovFoq985SuKjf3HIbxer1paWjR9+vTBPiwAAAAGaNAD4OzZsyVJTU1NmjZtms4//3zfz+Lj4zVmzBiWgQEAAAijQQ+A5eXlkqQxY8Zo7ty5fo+BAwAAQPiFbBmYu+66K1S7BgAAwDkIWQD0er16+umn9eKLL6q1tVXd3d1+Pz927FioDg0AAIBTCNm3gJcvX66nnnpKt912mzwej0pLS3XLLbcoJiZGjzzyyBnvb+XKlXI4HLLZbHI6ndq1a9cp++/cuVNOp1M2m01jx47VqlWr+u37wgsvyGKx+O5fBAAAGM5CFgD/8z//U2vWrNHixYsVGxurO+64Q2vXrtUPfvAD/eEPfzijfW3cuFHFxcVatmyZGhsbVVBQoBkzZqi1tTVo/5aWFs2cOVMFBQVqbGzU0qVLtXDhQm3atCmg76FDh7R48WIVFBSc1XkCAABEm5AtBD1ixAg1NzfrkksuUUZGhn79618rNzdXBw8e1JVXXimPxzPgfU2ePFm5ubmqrKz0tU2cOFGzZ89WRUVFQP+ysjJt3bpVzc3NvraioiLt3btXdXV1vjav16vrrrtOd999t3bt2qUPP/xQL7/88oDrYiFJAACiD9fvEI4AXnzxxXK5XJKk8ePHq6amRpL05ptvBjwe7lS6u7u1Z88eFRYW+rUXFhZq9+7dQbepq6sL6D9t2jS99dZb6unp8bWtWLFCF154ob71rW8NqJauri51dHT4vQAAAKJNyALgV7/6Vb366quSpEWLFunhhx/WpZdeqvnz5+uee+4Z8H7a29vl9XqVlpbm156Wlia32x10G7fbHbT/iRMn1N7eLkn6/e9/r3Xr1mnNmjUDrqWiokJ2u933ysrKGvC2AAAAkSJk3wL+4Q9/6PvvW2+9VVlZWfr973+v8ePH66abbjrj/X32sXL9PWruVP372o8fP6558+ZpzZo1GjVq1IBrWLJkiUpLS33vOzo6CIEAACDqhCQA9vT06L777tPDDz+ssWPHSjp5H9/kyZPPeF+jRo2S1WoNGO1ra2sLGOXrk56eHrR/bGysUlJStH//fr333nuaNWuW7+e9vb2SpNjYWL399tsaN25cwH4TEhLOaPoaAAAgEoVkCjguLk5btmwZlH3Fx8fL6XSqtrbWr722tlZTpkwJuk1+fn5A/5qaGuXl5SkuLk4TJkzQn/70JzU1NfleN910k7785S+rqamJUT0AADCshfQewDP5Ru2plJaWau3ataqqqlJzc7NKSkrU2tqqoqIiSSenZufPn+/rX1RUpEOHDqm0tFTNzc2qqqrSunXrtHjxYkmSzWZTTk6O3+uCCy5QYmKicnJyFB8fPyh1AwAARKKQ3QM4fvx4Pfroo9q9e7ecTqdGjBjh9/OFCxcOeF9z587V0aNHtWLFCrlcLuXk5Ki6ulqjR4+WJLlcLr81AR0Oh6qrq1VSUqLnnntOmZmZeuaZZzRnzpzBOTkAAIAoFrJ1AB0OR/8HtVh08ODBUBx2SLGOEAAA0YfrdwhHAFtaWkK1awAAAJyDkN0DCAAAgMhEAAQAADAZAiAAAIDJEAABAABMhgAIAABgMiH7FrAkffjhh6qvr1dbW5vvUWt9Pr1wMwAAAIZOyALgtm3b9I1vfEMff/yxEhMTZbFYfD+zWCwEQAAAgDAJ2RTw9773Pd1zzz06fvy4PvzwQ/3973/3vY4dOxaqwwIAAOA0QhYAjxw5ooULF+pzn/tcqA4BAACAsxCyADht2jS99dZbodo9AAAAztKg3gO4detW33/feOON+v73v68DBw7osssuU1xcnF/fm266aTAPDQAAgAGyGIZhDNbOYmIGNqBosVjk9XoH67Bhw8OkAQCIPly/B3kE8LNLvQAAACDyhOwewPXr16urqyugvbu7W+vXrw/VYQEAAHAagzoF/GlWq1Uul0upqal+7UePHlVqaipTwAAAICy4fodwBNAwDL/Fn/u8//77stvtoTosAAAATmPQnwRy5ZVXymKxyGKx6Ctf+YpiY/9xCK/Xq5aWFk2fPn2wDwsAAIABGvQAOHv2bElSU1OTpk2bpvPPP9/3s/j4eI0ZM0Zz5swZ7MMCABDRvL2G6luOqe14p1ITbZrkSJY1JnCmDBgKgx4Ay8vLJUljxozR7bffroSEhME+BAAAUWX7PpeWbzsgl6fT15Zht6l8Vram52SEsTKYVcjuAXznnXf0+uuv65NPPgnVIQAAiHjb97m0YEODX/iTJLenUws2NGj7PleYKoOZhSwA7tmzR3PmzNHIkSOVn5+vJUuWaPv27froo49CdUgAACKKt9fQ8m0HFGy5jb625dsOyNsbkgU5gH6FLABu375df//737Vjxw7dfPPNamxs1Ny5c5WcnKyrrroqVIcFACBi1LccCxj5+zRDksvTqfqWY0NXFKAQ3AP4aVarVfn5+UpOTtbIkSOVmJiol19+WX/9619DeVgAACJC2/H+w9/Z9AMGS8hGACsrK3X77bcrIyNDBQUFqqmpUUFBgfbs2aO//e1voTosAAARIzXRNqj9gMESshHA+++/XxdeeKG+973vqaioyLQrbQMAzGuSI1kZdpvcns6g9wFaJKXbTy4JAwylkI0Abt68Wd/4xjf0wgsvKDU1VZMnT1ZZWZleeeUVvggCADAFa4xF5bOyJZ0Me5/W9758VjbrAWLIhexZwJ/m8Xi0a9cuvfTSS/rlL38pi8Wirq6uUB825HiWIABgIFgHMLJw/Q7xl0COHTumnTt3aseOHdqxY4f27dunlJQUXXfddaE8LAAAEWV6ToamZqfzJBBEjJAFwMsvv1wHDhxQcnKyrr32Wt177726/vrrlZOTE6pDAgAQsawxFuWPSwl3GYCkEAbA++67j8AHAAAQgUIWAB944AHff/fdZmixMNQNAAAQbiH7FrAkrV+/XpdddpnOO+88nXfeebr88sv1H//xH6E8JAAAAE4jZCOATz31lB5++GE98MADuvrqq2UYhn7/+9+rqKhI7e3tKikpCdWhAQAAcAohWwbG4XBo+fLlmj9/vl/7L37xCz3yyCNqaWkJxWGHFF8jBwAg+nD9DuEUsMvl0pQpUwLap0yZIpfLdcb7W7lypRwOh2w2m5xOp3bt2nXK/jt37pTT6ZTNZtPYsWO1atUqv59v3rxZeXl5uuCCCzRixAhdccUVTE8DAABTCFkAHD9+vF588cWA9o0bN+rSSy89o31t3LhRxcXFWrZsmRobG1VQUKAZM2aotbU1aP+WlhbNnDlTBQUFamxs1NKlS7Vw4UJt2rTJ1yc5OVnLli1TXV2d/vjHP+ruu+/W3Xffrd/85jdndqIAAABRJmRTwJs2bdLcuXN1ww036Oqrr5bFYtEbb7yhV199VS+++KK++tWvDnhfkydPVm5uriorK31tEydO1OzZs1VRURHQv6ysTFu3blVzc7OvraioSHv37lVdXV2/x8nNzdWNN96oRx99dEB1MYQMAED04fodwhHAOXPmqL6+XqNGjdLLL7+szZs3a9SoUaqvrz+j8Nfd3a09e/aosLDQr72wsFC7d+8Ouk1dXV1A/2nTpumtt95ST09PQH/DMPTqq6/q7bff1rXXXjvg2gAAAKJRSL4F3NPTo/vuu08PP/ywNmzYcE77am9vl9frVVpaml97Wlqa3G530G3cbnfQ/idOnFB7e7syMk4+d9Hj8eiiiy5SV1eXrFarVq5cqalTp/ZbS1dXl98zjDs6Os72tAAAAMImJCOAcXFx2rJly6Du87OLSBuGccqFpYP1/2x7YmKimpqa9Oabb+qxxx5TaWmpduzY0e8+KyoqZLfbfa+srKyzOBMAAIDwCtkU8Fe/+lW9/PLL57yfUaNGyWq1Boz2tbW1BYzy9UlPTw/aPzY2Vikp/3gOY0xMjMaPH68rrrhC3/ve93TrrbcGvaewz5IlS+TxeHyvw4cPn8OZAQAAhEfIFoIeP368Hn30Ue3evVtOp1MjRozw+/nChQsHtJ/4+Hg5nU7V1tb63TtYW1urm2++Oeg2+fn52rZtm19bTU2N8vLyFBcX1++xDMPwm+L9rISEBCUkJAyobgAAgEgV0oWg+z2oxaKDBw8OeF8bN27UnXfeqVWrVik/P1+rV6/WmjVrtH//fo0ePVpLlizRkSNHtH79ekknl4HJycnRd77zHd17772qq6tTUVGRnn/+ec2ZM0fSyencvLw8jRs3Tt3d3aqurlZZWZkqKyv17W9/e0B18S0iAACiD9fvEI4ADuaTPubOnaujR49qxYoVcrlcysnJUXV1tUaPHi3p5KLTn14T0OFwqLq6WiUlJXruueeUmZmpZ555xhf+JOnjjz/Wv/7rv+r999/XeeedpwkTJmjDhg2aO3fuoNUNAAAQiUI2AmgG/AUBAED04fodwhHA0tLSoO0Wi0U2m03jx4/XzTffrOTk5FCVAAAAgCBCNgL45S9/WQ0NDfJ6vfqnf/onGYahd999V1arVRMmTNDbb7/tezpIdnZ2KEoIOf6CAAAg+nD9DuEyMDfffLNuuOEGffDBB9qzZ48aGhp05MgRTZ06VXfccYeOHDmia6+9ViUlJaEqAQAAAEGEbATwoosuUm1tbcDo3v79+1VYWKgjR46ooaFBhYWFam9vD0UJIcdfEAAARB+u3yEcAfR4PGprawto/9vf/uZ7hNoFF1yg7u7uUJUAAACAIEI6BXzPPfdoy5Ytev/993XkyBFt2bJF3/rWtzR79mxJUn19vT7/+c+HqgQAAAAEEbIp4I8++kglJSVav369Tpw4IUmKjY3VXXfdpaefflojRoxQU1OTJOmKK64IRQkhxxAyAADRh+v3EKwD+NFHH+ngwYMyDEPjxo3T+eefH8rDDSn+AQEAEH24fodwHcA+559/vi6//PJQHwYAAAADFLJ7AAEAABCZCIAAAAAmQwAEAAAwGQIgAACAyRAAAQAATIYACAAAYDIEQAAAAJMhAAIAAJhMyBeCBszG22uovuWY2o53KjXRpkmOZFljLOEuCwAAHwIgMIi273Np+bYDcnk6fW0ZdpvKZ2Vrek5GGCsDAOAfmAIGBsn2fS4t2NDgF/4kye3p1IINDdq+zxWmygAA8EcABAaBt9fQ8m0HZAT5WV/b8m0H5O0N1gMAgKFFAAQGQX3LsYCRv08zJLk8napvOTZ0RQEA0A8CIDAI2o73H/7Oph8AAKFEAAQGQWqibVD7AQAQSgRAYBBMciQrw25Tf4u9WHTy28CTHMlDWRYAAEERAIFBYI2xqHxWtiQFhMC+9+WzslkPEAAQEQiAwCCZnpOhynm5Srf7T/Om222qnJfLOoAAgIjBQtDAIJqek6Gp2ek8CQQAENEIgMAgs8ZYlD8uJdxlAADQL6aAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmQwAEAAAwmagJgCtXrpTD4ZDNZpPT6dSuXbtO2X/nzp1yOp2y2WwaO3asVq1a5ffzNWvWqKCgQCNHjtTIkSN1ww03qL6+PpSnAAAAEBGiIgBu3LhRxcXFWrZsmRobG1VQUKAZM2aotbU1aP+WlhbNnDlTBQUFamxs1NKlS7Vw4UJt2rTJ12fHjh2644479Lvf/U51dXW65JJLVFhYqCNHjgzVaQEAAISFxTAMI9xFnM7kyZOVm5uryspKX9vEiRM1e/ZsVVRUBPQvKyvT1q1b1dzc7GsrKirS3r17VVdXF/QYXq9XI0eO1LPPPqv58+cPqK6Ojg7Z7XZ5PB4lJSWd4VkBAIBw4PodBSOA3d3d2rNnjwoLC/3aCwsLtXv37qDb1NXVBfSfNm2a3nrrLfX09ATd5pNPPlFPT4+Sk5P7raWrq0sdHR1+LwAAgGgT8QGwvb1dXq9XaWlpfu1paWlyu91Bt3G73UH7nzhxQu3t7UG3efDBB3XRRRfphhtu6LeWiooK2e123ysrK+sMzwYAACD8Ij4A9rFYLH7vDcMIaDtd/2DtkvTkk0/q+eef1+bNm2Wz2frd55IlS+TxeHyvw4cPn8kpAAAARITYcBdwOqNGjZLVag0Y7WtrawsY5euTnp4etH9sbKxSUlL82n/84x/r8ccf129/+1tdfvnlp6wlISFBCQkJZ3EWAAAAkSPiRwDj4+PldDpVW1vr115bW6spU6YE3SY/Pz+gf01NjfLy8hQXF+dr+9GPfqRHH31U27dvV15e3uAXDwAAEIEiPgBKUmlpqdauXauqqio1NzerpKREra2tKioqknRyavbT39wtKirSoUOHVFpaqubmZlVVVWndunVavHixr8+TTz6phx56SFVVVRozZozcbrfcbrc++uijIT8/AACAoRTxU8CSNHfuXB09elQrVqyQy+VSTk6OqqurNXr0aEmSy+XyWxPQ4XCourpaJSUleu6555SZmalnnnlGc+bM8fVZuXKluru7deutt/odq7y8XI888siQnBcAAEA4RMU6gJGKdYQAAIg+XL+jZAoYAAAAg4cACAAAYDJRcQ8g0B9vr6H6lmNqO96p1ESbJjmSZY3pf31IAABAAEQU277PpeXbDsjl6fS1ZdhtKp+Vrek5GWGsDACAyMYUMKLS9n0uLdjQ4Bf+JMnt6dSCDQ3avs8VpsoAAIh8BEBEHW+voeXbDijY19f72pZvOyBvL19wBwAgGAIgok59y7GAkb9PMyS5PJ2qbzk2dEUBABBFCICIOm3H+w9/Z9MPAACzIQAi6qQm2ga1HwAAZkMARNSZ5EhWht2m/hZ7sejkt4EnOZKHsiwAAKIGARBRxxpjUfmsbEkKCIF978tnZbMeIAAA/SAAIipNz8lQ5bxcpdv9p3nT7TZVzstlHUAAAE6BhaARtabnZGhqdjpPAgEA4AwRABHVrDEW5Y9LCXcZAABEFaaAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmQwAEAAAwGQIgAACAyRAAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmQwAEAAAwGQIgAACAyRAAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkoiYArly5Ug6HQzabTU6nU7t27Tpl/507d8rpdMpms2ns2LFatWqV38/379+vOXPmaMyYMbJYLPrpT38awuoBAAAiR1QEwI0bN6q4uFjLli1TY2OjCgoKNGPGDLW2tgbt39LSopkzZ6qgoECNjY1aunSpFi5cqE2bNvn6fPLJJxo7dqx++MMfKj09fahOBQAAIOwshmEY4S7idCZPnqzc3FxVVlb62iZOnKjZs2eroqIioH9ZWZm2bt2q5uZmX1tRUZH27t2rurq6gP5jxoxRcXGxiouLz6iujo4O2e12eTweJSUlndG2AAAgPLh+R8EIYHd3t/bs2aPCwkK/9sLCQu3evTvoNnV1dQH9p02bprfeeks9PT0hqxUAACAaxIa7gNNpb2+X1+tVWlqaX3taWprcbnfQbdxud9D+J06cUHt7uzIyMs6qlq6uLnV1dfned3R0nNV+AAAAwiniRwD7WCwWv/eGYQS0na5/sPYzUVFRIbvd7ntlZWWd9b4AAADCJeID4KhRo2S1WgNG+9ra2gJG+fqkp6cH7R8bG6uUlJSzrmXJkiXyeDy+1+HDh896XwAAAOES8QEwPj5eTqdTtbW1fu21tbWaMmVK0G3y8/MD+tfU1CgvL09xcXFnXUtCQoKSkpL8XgAAANEm4gOgJJWWlmrt2rWqqqpSc3OzSkpK1NraqqKiIkknR+bmz5/v619UVKRDhw6ptLRUzc3Nqqqq0rp167R48WJfn+7ubjU1NampqUnd3d06cuSImpqa9Je//GXIzw8AAGAoRcUyMNLJhaCffPJJuVwu5eTk6Omnn9a1114rSfrmN7+p9957Tzt27PD137lzp0pKSrR//35lZmaqrKzMFxgl6b333pPD4Qg4znXXXee3n1Mx09fIvb2G6luOqe14p1ITbZrkSJY15uzvpwQAIFzMdP3uT9QEwEhkln9A2/e5tHzbAbk8nb62DLtN5bOyNT3n7L5RDQBAuJjl+n0qUTEFjPDZvs+lBRsa/MKfJLk9nVqwoUHb97nCVBkAADhbBED0y9traPm2Awo2RNzXtnzbAXl7GUQGACCaEADRr/qWYwEjf59mSHJ5OlXfcmzoigIAAOeMAIh+tR3vP/ydTT8AABAZCIDoV2qibVD7AQCAyEAARL8mOZKVYbepv8VeLDr5beBJjuShLAsAAJwjAiD6ZY2xqHxWtiQFhMC+9+WzslkPEACAKEMAxClNz8lQ5bxcpdv9p3nT7TZVzstlHUAAAKJQbLgLQOSbnpOhqdnpPAkEAIBhggCIAbHGWJQ/LiXcZQAAgEHAFDAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGRYBgYAABPz9hqs82pCBEAAAExq+z6Xlm87IJen09eWYbepfFY2T3oa5pgCBgDAhLbvc2nBhga/8CdJbk+nFmxo0PZ9rjBVhqFAAAQAwGS8vYaWbzsgI8jP+tqWbzsgb2+wHhgOCIAAAJhMfcuxgJG/TzMkuTydqm85NnRFYUgRAAEAMJm24/2Hv7Pph+hDAAQAwGRSE22D2g/RhwAIAIDJTHIkK8NuU3+LvVh08tvAkxzJQ1kWhhABEAAAk7HGWFQ+K1uSAkJg3/vyWdmsBziMEQABADCh6TkZqpyXq3S7/zRvut2mynm5rAM4zLEQNAAAJjU9J0NTs9N5EogJEQABADAxa4xF+eNSwl0GhhhTwAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAzrAEYgb6/BopwAACBkomYEcOXKlXI4HLLZbHI6ndq1a9cp++/cuVNOp1M2m01jx47VqlWrAvps2rRJ2dnZSkhIUHZ2trZs2RKq8gds+z6XrnniNd2x5g9a9EKT7ljzB13zxGvavs8V7tIAABgS3l5DdX89ql81HVHdX4/K22uEu6RhJyoC4MaNG1VcXKxly5apsbFRBQUFmjFjhlpbW4P2b2lp0cyZM1VQUKDGxkYtXbpUCxcu1KZNm3x96urqNHfuXN15553au3ev7rzzTt1222367//+76E6rQDb97m0YEODXJ5Ov3a3p1MLNjQQAgEAwx4DIUPDYhhGxMfqyZMnKzc3V5WVlb62iRMnavbs2aqoqAjoX1ZWpq1bt6q5udnXVlRUpL1796qurk6SNHfuXHV0dOiVV17x9Zk+fbpGjhyp559/fkB1dXR0yG6364MPPlBSUtLZnp6kk3/tTH3mD/p/x7uC/twiKS0pQTXfvYrpYADAsFTb/DeVvLRfnw0mfVe9p2/9gqZOvPCcj9PR0aHMzEx5PJ5zvn5Hq4i/B7C7u1t79uzRgw8+6NdeWFio3bt3B92mrq5OhYWFfm3Tpk3TunXr1NPTo7i4ONXV1amkpCSgz09/+tN+a+nq6lJX1z8CWkdHhyQpMzPzTE4pqISsy5T+9cAw28eQ5O7oUuoXpqjr8J/O+XgAAEQUS4wuKlona+IoWSz+Ax2GJMPo1QNVO3Rk1bckozc8NQ4jET8F3N7eLq/Xq7S0NL/2tLQ0ud3uoNu43e6g/U+cOKH29vZT9ulvn5JUUVEhu93ue2VlZZ3NKQVlPX/koPYDACCaJFz8BcUmXRgQ/vpYLDGKTbpQCRd/YYgrG54ifgSwT8BfA4bR7z+S/vp/tv1M97lkyRKVlpb63nd0dCgrK2tQpoDr3/u77v6Pvaft91//9z81aQwhEAAwvPx63//T/9nSfNp+6//vy7oxJ+20/U6lbwrYzCI+AI4aNUpWqzVgZK6trS1gBK9Penp60P6xsbFKSUk5ZZ/+9ilJCQkJSkhICGgfMWKERowYMaDz6c+1Ez+nDPvbcns6A+59kE7e/5But+naiRdxDyAAYNjJGmUfcL9zveZ6vd5z2n44iPgp4Pj4eDmdTtXW1vq119bWasqUKUG3yc/PD+hfU1OjvLw8xcXFnbJPf/sMNWuMReWzsiX942bXPn3vy2dlE/4AAMPSJEeyMuy2gGtgH4ukDPvJtXFx7iI+AEpSaWmp1q5dq6qqKjU3N6ukpEStra0qKiqSdHJqdv78+b7+RUVFOnTokEpLS9Xc3KyqqiqtW7dOixcv9vVZtGiRampq9MQTT+jPf/6znnjiCf32t79VcXHxUJ+ez/ScDFXOy1W63ebXnm63qXJerqbnZISpMgAAQouBkKEVFcvASCcXgn7yySflcrmUk5Ojp59+Wtdee60k6Zvf/Kbee+897dixw9d/586dKikp0f79+5WZmamysjJfYOzz0ksv6aGHHtLBgwc1btw4PfbYY7rlllsGXFPfMjCD/TVyngQCADCr7ftcWr7tgN+auBl2m8pnZQ/aQEiort/RJGoCYCTiHxAAAIMv1AMhXL+j4EsgAADAXKwxFuWPSwl3GcNaVNwDCAAAgMFDAAQAADAZAiAAAIDJEAABAABMhgAIAABgMgRAAAAAkyEAAgAAmAwBEAAAwGQIgAAAACbDk0DOQd9T9Do6OsJcCQAAGKi+67aZn4ZLADwHx48flyRlZWWFuRIAAHCmjh8/LrvdHu4ywsJimDn+nqPe3l598MEHSkxMlMUyeA+plk7+dZKVlaXDhw+b9kHVA8VnNXB8VgPHZzVwfFYDx2c1cKH8rAzD0PHjx5WZmamYGHPeDccI4DmIiYnRxRdfHNJjJCUl8UtigPisBo7PauD4rAaOz2rg+KwGLlSflVlH/vqYM/YCAACYGAEQAADAZAiAESohIUHl5eVKSEgIdykRj89q4PisBo7PauD4rAaOz2rg+KxCiy+BAAAAmAwjgAAAACZDAAQAADAZAiAAAIDJEAABAABMhgAYgVauXCmHwyGbzSan06ldu3aFu6SIU1FRoS996UtKTExUamqqZs+erbfffjvcZUWFiooKWSwWFRcXh7uUiHXkyBHNmzdPKSkp+tznPqcrrrhCe/bsCXdZEefEiRN66KGH5HA4dN5552ns2LFasWKFent7w11a2L3++uuaNWuWMjMzZbFY9PLLL/v93DAMPfLII8rMzNR5552n66+/Xvv37w9PsWF2qs+qp6dHZWVluuyyyzRixAhlZmZq/vz5+uCDD8JX8DBBAIwwGzduVHFxsZYtW6bGxkYVFBRoxowZam1tDXdpEWXnzp26//779Yc//EG1tbU6ceKECgsL9fHHH4e7tIj25ptvavXq1br88svDXUrE+vvf/66rr75acXFxeuWVV3TgwAH95Cc/0QUXXBDu0iLOE088oVWrVunZZ59Vc3OznnzySf3oRz/Sz372s3CXFnYff/yxvvjFL+rZZ58N+vMnn3xSTz31lJ599lm9+eabSk9P19SpU33PmDeTU31Wn3zyiRoaGvTwww+roaFBmzdv1jvvvKObbropDJUOMwYiyqRJk4yioiK/tgkTJhgPPvhgmCqKDm1tbYYkY+fOneEuJWIdP37cuPTSS43a2lrjuuuuMxYtWhTukiJSWVmZcc0114S7jKhw4403Gvfcc49f2y233GLMmzcvTBVFJknGli1bfO97e3uN9PR044c//KGvrbOz07Db7caqVavCUGHk+OxnFUx9fb0hyTh06NDQFDVMMQIYQbq7u7Vnzx4VFhb6tRcWFmr37t1hqio6eDweSVJycnKYK4lc999/v2688UbdcMMN4S4lom3dulV5eXn62te+ptTUVF155ZVas2ZNuMuKSNdcc41effVVvfPOO5KkvXv36o033tDMmTPDXFlka2lpkdvt9vtdn5CQoOuuu47f9QPg8XhksVgYlT9HseEuAP/Q3t4ur9ertLQ0v/a0tDS53e4wVRX5DMNQaWmprrnmGuXk5IS7nIj0wgsvqKGhQW+++Wa4S4l4Bw8eVGVlpUpLS7V06VLV19dr4cKFSkhI0Pz588NdXkQpKyuTx+PRhAkTZLVa5fV69dhjj+mOO+4Id2kRre/3ebDf9YcOHQpHSVGjs7NTDz74oL7+9a8rKSkp3OVENQJgBLJYLH7vDcMIaMM/PPDAA/rjH/+oN954I9ylRKTDhw9r0aJFqqmpkc1mC3c5Ea+3t1d5eXl6/PHHJUlXXnml9u/fr8rKSgLgZ2zcuFEbNmzQL3/5S33hC19QU1OTiouLlZmZqbvuuivc5UU8ftefmZ6eHt1+++3q7e3VypUrw11O1CMARpBRo0bJarUGjPa1tbUF/KWIk7773e9q69atev3113XxxReHu5yItGfPHrW1tcnpdPravF6vXn/9dT377LPq6uqS1WoNY4WRJSMjQ9nZ2X5tEydO1KZNm8JUUeT6/ve/rwcffFC33367JOmyyy7ToUOHVFFRQQA8hfT0dEknRwIzMjJ87fyu719PT49uu+02tbS06LXXXmP0bxBwD2AEiY+Pl9PpVG1trV97bW2tpkyZEqaqIpNhGHrggQe0efNmvfbaa3I4HOEuKWJ95Stf0Z/+9Cc1NTX5Xnl5efrGN76hpqYmwt9nXH311QFLCr3zzjsaPXp0mCqKXJ988oliYvwvI1arlWVgTsPhcCg9Pd3vd313d7d27tzJ7/og+sLfu+++q9/+9rdKSUkJd0nDAiOAEaa0tFR33nmn8vLylJ+fr9WrV6u1tVVFRUXhLi2i3H///frlL3+pX/3qV0pMTPSNmtrtdp133nlhri6yJCYmBtwbOWLECKWkpHDPZBAlJSWaMmWKHn/8cd12222qr6/X6tWrtXr16nCXFnFmzZqlxx57TJdccom+8IUvqLGxUU899ZTuueeecJcWdh999JH+8pe/+N63tLSoqalJycnJuuSSS1RcXKzHH39cl156qS699FI9/vjj+tznPqevf/3rYaw6PE71WWVmZurWW29VQ0OD/uu//kter9f3+z45OVnx8fHhKjv6hfdLyAjmueeeM0aPHm3Ex8cbubm5LG0ShKSgr3//938Pd2lRgWVgTm3btm1GTk6OkZCQYEyYMMFYvXp1uEuKSB0dHcaiRYuMSy65xLDZbMbYsWONZcuWGV1dXeEuLex+97vfBf0ddddddxmGcXIpmPLyciM9Pd1ISEgwrr32WuNPf/pTeIsOk1N9Vi0tLf3+vv/d734X7tKjmsUwDGMoAycAAADCi3sAAQAATIYACAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmQwAEAAAwGQIgAACAyRAAAQAATIYACAAAYDIEQACmd/311+u73/2uiouLNXLkSKWlpWn16tX6+OOPdffddysxMVHjxo3TK6+8Eu5SAWBQEAABQNIvfvELjRo1SvX19frud7+rBQsW6Gtf+5qmTJmihoYGTZs2TXfeeac++eSTcJcKAOfMYhiGEe4iACCcrr/+enm9Xu3atUuS5PV6Zbfbdcstt2j9+vWSJLfbrYyMDNXV1emqq64KZ7kAcM4YAQQASZdffrnvv61Wq1JSUnTZZZf52tLS0iRJbW1tQ14bAAw2AiAASIqLi/N7b7FY/NosFoskqbe3d0jrAoBQIAACAACYDAEQAADAZAiAAAAAJsO3gAEAAEyGEUAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJvP/AQ72iS10fZbfAAAAAElFTkSuQmCC", "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": "804f113c9f0740098ebf3385a8a253dd", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACPQElEQVR4nO3deXwV1d0/8M/cPTtLICFlValCsUrB0uCGbQXcCrWtW4tQLdXHBRF9tBQXXABFpfwU19a6tLX69FHcF/CpohRwobgjrRUFgcgWkkCSu835/THnzHLvTXKTzM29yf28X6/7usncuTPznTMzOZnvOWc0IYQAEREREeUNT7Y3gIiIiIi6FiuARERERHmGFUAiIiKiPMMKIBEREVGeYQWQiIiIKM+wAkhERESUZ1gBJCIiIsozrAASERER5RlWAImIiIjyDCuARERERHmGFUAiIiKiPMMKIBEREVGeYQWQiIiIKM/4sr0B3YGu69i+fTtKSkqgaVq2N4eIiKhDhBBoaGhAVVUVPJ6uuQfU3NyMSCTi2vICgQBCoZBry8tXrACmYfv27Rg0aFC2N4OIiMgVW7duxcCBAzO+nubmZgwbUoyanXHXlllZWYnNmzezEthJrACmoaSkBAAw6IZr4OmKA05k+S6jJrK37mzGns24gfyNPVtxZ7u8gfws82xf34D8jF3GrDc3Y+v1N5t/1zItEomgZmccX64fitKSzt9xrG/QMWTMF4hEIqwAdhIrgGlQaV9PKMQKYKblwAUya/I1dlYAsyMfK0FKPsaeEHNXN2cqLtFQXNL5derIgeOnh2AFMBepEzUXLpREPZU6v3KhIkjUw8WFjrgLp1pc6J1fCAFgL2AiIiKivMM7gERERJRROgR0dP4WoBvLIAMrgLlME9lJA2czNZa4zq6M376ubMaezTK3b0dXsa8vX4/3fCrzbJ7j9nXm6/UtS02LdOhwI3nrzlIIYAqYiIiIKO/wDiARERFlVFwIxEXn77q6sQwysALYHkLr/O3z9qYd3ExTtHfb3UwVdDQON+LvSBzZjD2bZd7R7yTKZnkr3el4dyvujsbQE8q8O5V3Z79nx9EiqINYASQiIqKMYieQ3MMKIBEREWWUDoE4K4A5hZ1AiIiIiPIM7wB2pXx8/FA+y8fHq+WCHHrcFnWRfH6cYDfBFHDuYQWQiIiIMoq9gHMPU8BEREREeYZ3AImIiCijdPlyYznkjm59BzAWi+Gaa67BsGHDUFBQgIMOOgg33ngjdN06RIQQmD9/PqqqqlBQUIAJEybg448/zuJWExER5Ze47AXsxovc0a0rgLfeeivuu+8+LFu2DBs3bsTixYtx22234a677jLnWbx4MZYsWYJly5bhnXfeQWVlJU488UQ0NDRkccuJiIiIsqdbp4DXrl2LKVOm4JRTTgEADB06FH/961/x7rvvAjDu/i1duhTz5s3D6aefDgB45JFHUFFRgcceewwXXHBB1radiIgoX8SF8XJjOeSObn0H8JhjjsH//d//4V//+hcA4P3338fq1atx8sknAwA2b96MmpoaTJw40fxOMBjE8ccfjzVr1mRlm4mIiIiyrVvfAbz66qtRV1eHww47DF6vF/F4HAsWLMDZZ58NAKipqQEAVFRUOL5XUVGBL7/8ssXlhsNhhMNh8/f6+voMbD0REVF+YCeQ3NOtK4BPPPEE/vznP+Oxxx7Dt771Lbz33nuYPXs2qqqqMH36dHM+TXMOEiqESJpmt2jRItxwww0Z224iIqJ8okNDHJ0fsFt3YRlk6NYp4P/+7//Gb37zG5x11lk4/PDDMW3aNFx++eVYtGgRAKCyshKAdSdQ2blzZ9JdQbu5c+eirq7OfG3dujVzQRARERF1sW5dAWxsbITH4wzB6/Waw8AMGzYMlZWVWLlypfl5JBLBqlWrMH78+BaXGwwGUVpa6ngRERFRx+jCvRe5o1ungE877TQsWLAAgwcPxre+9S1s2LABS5YswXnnnQfASP3Onj0bCxcuxPDhwzF8+HAsXLgQhYWFOOecc7K89URERPkh7lIK2I1lkKFbVwDvuusuXHvttbjooouwc+dOVFVV4YILLsB1111nznPVVVehqakJF110EWprazFu3DisWLECJSUlWdxyIiIiouzRhOCTldtSX1+PsrIyDLllATyhUMcXpGVpV4sc+I8pW7ED2Ys/H2NW8jX2bMYN5G/s+XiOK+2MXW9uxpe/mYe6urouad6k/n6u+XgAiks63+psf4OO8d/a0WXb35N16zaARERERNR+3ToFTERERLlPFxp0F+7UurEMMrACSERERBnFTiC5hylgIiIiojzDO4BERESUUXF4EHfhnlPchW0hAyuARERElFHCpTaAgm0AXcMUMBEREVGe4R1Aop5IjU/G/5aJMkedX7kwHmCOYyeQ3MM7gERERER5hncAiYiIKKPiwoO4cKETCG+2uoYVQCIiIsooHRp0F5KOOlgDdAtTwERERER5hncAiYiIKKPYCST3sAJIREREGeVeG0CmgN3CFDARERFRnuEdQCIiIsoooxNI59O3biyDDKwAtoeQL7sMHIvpjt2bsbFH27Ncl+PPauzpLjOfy7wnxZ7F49zchDw83vMxZnP17Vhuxo73LGVQdZeeBcxewO5hCpiIiIgoz/AOIBEREWUUO4HkHt4BJCIiIsozvAPYWS39M5JmW4/2tAlJ57vtbjfS2X+m7N9vZyxZjb0zcaf6bhfFnvUy72TsWStzN49zc+XtXIRLsWf1HDc3Is2vulzeQBfFn+Xybun7XXa8Z4AOD58EkmNYASQiIqKMigsN8c7WiuVyyB1MARMRERHlGd4BbI9Uw8AkUv+cCLSaMsjUPzFquWmlCkQLP7eXPebEaalWm83Y3YpZSYw9i3ED7Yg91e/tYY8lB2JvV9wtTUtXmue4uaoMD6OSi8d7Jm/SCK2Lyxsw4mzHEGCZPN47fY5nKYMad2kYmDhTwK5hBZCIiIgyShce6C70AtbZC9g1TAETERER5RneAWwHTTdezonOX1P+b5JuOiBxWW18LykVYPvdkRpTM6qJKUaE19qbKrJtW9Ls7U2NaG18nuor6cae4nP1e4uplHTT/K3N2oEyTyf2dpeTfb6E+du1rMRtsx9SraQF2yrvFudJ9bVWtre9qTFXj3f757aDz43YWzvOW5Uivg4vK1XsraTB04m7xfnsX2nvOZ60gvSWlXrlrcye4etbypjSib2tMlfTEv+GdRGmgHMP7wASERER5RneASQiIqKM0uHOEC5ZuoHZI/EOYDuoFLD9BfWSPYTNz1q4Sy20hDSATKEJDyA8wvFCG6/E+YXHWp59fdYvzpcmUsSRIsZW49adMWsJ62h9h8rYNaQdc1Lsan+me12xxW5uYzvjt5d5OrE7ytxR3mhXeSeXOZLKPOn4SojZUea22Nss87jxSjre0y1rRbO2sVNlniL2lso6ZewuHe+J55S9zO0xZ/ocb/Hvqi32xOM9nbi1FOXdrjK3n+MaOha7+V1bmaeKPdX1TSSXX7vOcfv87bm+2c9xW+xpHede+Uq8vqVzjUtVTilizgY1ELQbr4645557MGzYMIRCIYwZMwZvvvlmq/OvWrUKY8aMQSgUwkEHHYT77rsvaZ4nn3wSI0eORDAYxMiRI7F8+XLH54sWLcJRRx2FkpIS9O/fH1OnTsWmTZs6tP2ZwAogERER9VhPPPEEZs+ejXnz5mHDhg049thjcdJJJ2HLli0p59+8eTNOPvlkHHvssdiwYQN++9vfYtasWXjyySfNedauXYszzzwT06ZNw/vvv49p06bhjDPOwFtvvWXOs2rVKlx88cVYt24dVq5ciVgshokTJ+LAgQMZjzkdmhDsU92W+vp6lJWVYdiNC+AJhRyfOe7s2Mn/AB2fpWoYbrtzk3TbsK3/+BJLTmjWImyfmdN0ZycQzf5fbKr5W1t14h0O27vjjpwt9pa+54i9PXfybAvUUvxHrtnnSYizpdgd32tp1anKXEvxmaeNmBNXlu5/+Akb0mqZC83MmaSar8UOAi2tPqlcU09LWd5J83euzO2LcHTqSdHZKSl2l473lLFD3u1JmD+T5zjUotsRe7vKPNU5rn5P9xxPXGm6sSccRCnPdb2Fc1z93tFz3L6t7Yk91V1ZTaR/rANyu9M8x835Ez9LnqY3N2PzdfNQV1eH0tLSdmxQx6i/n8vWj0NBcedbnTXtj+GSMW+1a/vHjRuH73znO7j33nvNaSNGjMDUqVOxaNGipPmvvvpqPPvss9i4caM57cILL8T777+PtWvXAgDOPPNM1NfX46WXXjLnmTx5Mnr37o2//vWvKbdj165d6N+/P1atWoXjjjsurW3PJLYBbActrkGLJ5zBrVz41a32pIpgqvlbqTS1KOliLyBkJc9cVYpFOP5oqnSAbVlp/WFo5Y+95rGtNvGiCfv+kHMlpHWcG5mCY2Fy3bo13dzvantSXRj1Fi6gbazamMH2rmL3OD9qqTegOd1e7gnLapWjMivsb84/gub8KX5OMX+6h5oZlq1yq5nHMhJmsn3fXuatVBjT2gghzEqtgHW8qyUk/rHUVBoMCWXenopQqj/2nrbL2/GZJ1XsHTjH7eVn/8wxn5oHzvJOmD+d00yzba+6XiWVuW0epJjfWdZprNy+ASpg9aZraV/fALUPOnGOyy/Yr2/molJc38zf7dc3NW97/9vQnec4hNZyCtdxPTc3Ozn2xL9hXUSHBr1dNeCWlwMYFUu7YDCIYDCYNH8kEsH69evxm9/8xjF94sSJWLNmTcp1rF27FhMnTnRMmzRpEh588EFEo1H4/X6sXbsWl19+edI8S5cubXHb6+rqAAB9+vRpcZ6uxBQwERERdSuDBg1CWVmZ+Up1Jw8Adu/ejXg8joqKCsf0iooK1NTUpPxOTU1NyvljsRh2797d6jwtLVMIgTlz5uCYY47BqFGj0oox03gHkIiIiDIqLjyIu/AkELWMrVu3OlLAqe7+2Wma8+6jECJpWlvzJ05vzzIvueQSfPDBB1i9enWr29mVuv0dwG3btuEXv/gF+vbti8LCQhx55JFYv369+bkQAvPnz0dVVRUKCgowYcIEfPzxx1ncYiIiIuqM0tJSx6ulCmB5eTm8Xm/SnbmdO3cm3cFTKisrU87v8/nQt2/fVudJtcxLL70Uzz77LF577TUMHDgw7RgzrVvfAaytrcXRRx+NE044AS+99BL69++P//znP+jVq5c5z+LFi7FkyRI8/PDD+OY3v4mbb74ZJ554IjZt2oSSkpJ2rU+LAh6vs4lGUlsY9aEasgFI2U7FbBuj5rHPb2scn+qfCSESfrC3fVPbo9o82dt7pGobo9u+q95baTBttoOxN3b2qJjV4oWznYxjYbDav9ljT2wXlLgeW7iAsLbbbPilWe1l4IzJMZ8Zp+aM2b6xKdpBtljm5gY62wKmLHMNgFfF7oxT2NqHtVnmCW0d1b4QEKnbeSXGaSvzxP0CpI7dXt4AzLam8KSIRaRoD2hvD+VJ2Fda6tiTuqep/8B12wapMte1lptX6Vb7t5THeSvtwhzHu1lW1j5O2cYt4buO492MXSTNn1juQgApz3FbeRvTWm7bl/Ict8edIvbk9p62WGzXN3N+r3O7jWtC8jkOqHMmnXMczvJOiE8knv/2zxPaujnbQVrztFnetm0UHs1xfQPgbOusvm+PNzF2+7VPLT5V7LYyV9dRswOfbp3/5qXPXq4J7T0dw77Iz+KxpLC7hHtPAmnfMgKBAMaMGYOVK1fixz/+sTl95cqVmDJlSsrvVFdX47nnnnNMW7FiBcaOHQu/32/Os3LlSkc7wBUrVmD8+PHm70IIXHrppVi+fDlef/11DBs2rF3bnmndugJ46623YtCgQXjooYfMaUOHDjV/FkJg6dKlmDdvHk4//XQAwCOPPIKKigo89thjuOCCC7p6k4mIiPKOLjTobgwE3YFlzJkzB9OmTcPYsWNRXV2NBx54AFu2bMGFF14IAJg7dy62bduGRx99FIDR43fZsmWYM2cOZs6cibVr1+LBBx909O697LLLcNxxx+HWW2/FlClT8Mwzz+DVV191pHgvvvhiPPbYY3jmmWdQUlJi3jEsKytDQUFBZ3aDK7p1CvjZZ5/F2LFj8bOf/Qz9+/fH6NGj8fvf/978fPPmzaipqXH05gkGgzj++ONb7P0DAOFwGPX19Y4XERERdT9nnnkmli5dihtvvBFHHnkk3njjDbz44osYMmQIAGDHjh2OMQGHDRuGF198Ea+//jqOPPJI3HTTTbjzzjvxk5/8xJxn/PjxePzxx/HQQw/h29/+Nh5++GE88cQTGDdunDnPvffei7q6OkyYMAEDBgwwX0888UTXBd+Kbn0H8PPPP8e9996LOXPm4Le//S3efvttzJo1C8FgEOeee65Z207VU+fLL79scbmLFi3CDTfckDTdEwM8UVs6xJEWSnj3akhKC9p/Vu8yJSi8Alpi2kQT0FrLUwhbWgAyDRZXM6nPgMQ0j2ZLD6v5002NJY57Jjy2lKhMBeleDUJ9KTGdaEuNCBm7Zo89KT2UYuX2oQ1UOkwXVrrb/Ip9fhWfTI2oJ1skxi7nb223w1bWSal/n7X9iZsBj7DK22PFDsihJTTVyLiNMlfDv8StWNSKhHP3OWNPmD/d1JhIcZzby9wc2iFx2AvbZsOr5he2mG1l3lZ526PShS12290A3fYz7OVqO87Vu7CnxJ2bb2cvc7Os1fe8tm1SaT77pqr47Oe4vbzt39BSlLtjyB/reBdmLPLYN3OBKVL6ugZPQuwph4hJ9xx3XN9gnOfOLLVMdSIpdkAe762c42am2f6ID3V9U2WuAVBpzBRjQlrXN+v35GYQaZ7j9t/NVLct9kSaVeYioaw1T+uxJw/5o1nXNxULNKszQuJwLsKa5ijzhP0ispQC1l1KAXf0SSAXXXQRLrroopSfPfzww0nTjj/+ePzzn/9sdZk//elP8dOf/rTFz3N9mOVuXQHUdR1jx47FwoULAQCjR4/Gxx9/jHvvvRfnnnuuOV97e//MnTsXc+bMMX+vr6/HoEGDXN56IiKi/KALD3QXegG7sQwydOs9OWDAAIwcOdIxbcSIEeat3MrKSgBoV+8fwEgTJ/YwIiIiIuopuvUdwKOPPjrpwcr/+te/zLz+sGHDUFlZiZUrV2L06NEAjFHBV61ahVtvvbXd6/M2W504AWcKUKU/hdyjRtZCpixUesgLR4oAgJky1Lw6PD5d/mxM83hSp4CFTBXoZnpA/h7z2HNLxme2HpBJqcCYLVUQU5+lSJfYl5ciDaqro8iWAdTlHVaR0OsNHkCoOG2xmzF75WcpUiQqbiEAPe5xxA7dY3sShC39raYlpIU8MSvmpBRZS+mhpDS/rdxV+tset0hIFXmFo7yNdxW3gOYxpnkS9xmcZW6Wt0duiEetx2Nut+ol6bGlheyxq98TU6NtpoAT4xVWptOM3Z72MI99GbdPN49983j3CDPmVMe7eZzrRrx6XANk7PYeqyKht6t1vDtjVvvAXt5AGilgFbM83m0PJLHSw17b91T6U5W5T5jl7jHTg62XuRW7dY5rqrzlSZj0qEdY57w9zsR94EgPtpYCTnWcq1k0LTn9bWvqAH/ica47rm/GMlKXuaO8AfM4R9zqfmuWuWadb2YatJXj3NE71ibp+mY2a7FS7cI83p3pdwDOpg5+ZxlrXpH29c2I22Nd38zjXYOmDrYUKf3Ea5kWs+0HFW9zctxdIQ4NcReeBOLGMsjQrSuAl19+OcaPH4+FCxfijDPOwNtvv40HHngADzzwAADj4jR79mwsXLgQw4cPx/Dhw7Fw4UIUFhbinHPOyfLWExER5QemgHNPt64AHnXUUVi+fDnmzp2LG2+8EcOGDcPSpUvx85//3JznqquuQlNTEy666CLU1tZi3LhxWLFiRbvHACQiIiLqKbp1BRAATj31VJx66qktfq5pGubPn4/58+d3el2+JsAbh6N3nD1FAAC6X84csDIymkoLAFarS5UW8hv36j1+HT75s0+mCXzeOLwpUkNxueBY3CvfjYXGYx7EI165LrnSuGalyWwpMUD2ao7In6PyPUXvWEdP3sQevx5AkzEL+R6HredbQg9BeIWZ+vWoeANxeM2YjXf1u8eWItFtaVAVcywm90HUa/YOU6kUlSLRYPWmU+kQT8QWsy12YwUtpIcS0/1eq7zjAbkuW5k7UmIwUoGaTHt6AzJ2W5n7vMbPrZW5rnsQVTHHZLlHbSPxxhNG5RVWXPbYVdxJacEUA4EbA/saP+oJTR00P6DL2M0ewvbY7alfAJpfh0/GrsrY74u3Wd6AdZzH4l4z5rhZ5h5ocWc+zp4GTSrraOq0YKJU6X5NNfMI2OZTvWLNbtjCmfYG4A3o8Jrl7TzXPR7hiF3FHzdjto73eFSmRs0Zra7Ziee4Fge8qrzlu9n0QXem/k1JPX6ttLd5fbPt6sRDDh5YqV/5bsbt0+H3yWueTI22dLwnXt9iqswjRnkb8cnvxm0Du8ec717buW7fL62VNxKOd80P6HJ+s9y9th755nVRXdeFdX0zY487yhto/XiPxrxW+cvYdc2b3CsaKiart7vHVuaehKYunqbkuLtCHO6kb+Ntz0Jp4r1UIiIiojzT7e8AEhERUW5jG8DcwwpgOwQaBLwBYd5z172A7pO37WVaIG67P504SLDQbL3jZFrIFzS+EAjEEAoYeYqQz7hnH/DG4PM4U2O60BCTKZ9I3Ci+5ph8j/gRlquMqfyzVwOizpSYPU3glT3CvFFrmicm0xP29JB9kFfY4ralRuK2dG/i80HN3sA+KzXiDxlxhgJRBGVaKOQ3NsQvV67iN2Iy4o7qXkRkzE1R4z3s8SMiD2eVQkHUyr2b6Q8ZpzecnBqzUsEiOf2t2dL8ttiT0kgqdaRr1odmKlDAF5RlG1CxG+8F/qgZs1+mBz2aMNNC9tjDMubmqJGPa1bh6pqVgrI9DzkxLWSWuSMNLqz5E2JyNHXwG8s1U96250/rXud3ACsl5pGpQH/Idpz75b7wxRD0Gj/by1vFHpWpwHBclbkfYY/xcwTGPogLDYiqlctts/eEDFsxJ8Ues2KHM4tspb59WlK63/65lpiX8li93L0BI6ZAMIqgrbwBmHH7vXHHOa7eVcz2Mg/LHHRUzRezpQJ1FZMVpypvs/xtZZ4q/Z3Y1EH3aVZ5J8QpvCl6jXuF2cPbL493FXfIH0VQlrv9HE+MPRr3WrGr65ss8zBkeQO2MteSBkH22srcYyt3Iw7hvL4BCU0drHMcAPQ4EFPXN9X53mtLw6sezTJujz+OgIw9FLSu6yGf8bNXHuctXd8Ao6ybzeubsfyoAOL28gYcIxykOse9Cal/0ZidwYnjwoO4C5U3N5ZBBu5JIiIiojzDO4BERESUUQKa9djITi6H3MEKYDv4D+jwRXQIj5UeUKmRmHo+pJkCtL4n7IPh+lXq15keKA6FURwwchbFfuOefcgbhS8hH6dDQ3PcyEs0xox3TQsCAOK2HqLWQzVtA6Pa0p8A4Gu2UgW+ZpmuigCeiEqJyXdrjFWzt6OZDotbsZudEVOkyOy9nlXqtyhkxFkcDKPEb2yUSn+HZJ7KY+ueGJFdEZvjPuz32FYC2Ts2mpB3TtULVKWFwkavbuNnGbtMg3qiVuxWHJrZE1KVeVxoVgowYTBsezoJfpUKjCMoy70kZGxIka3MQzIdGFD5O1jP2VXp/saYHwfkDlafmb2CPcI2KKwqcy0pLeRT6aFmYUuJyvKJi6Se28KjOcrbDBAyTZgiZWyVtzMVWBiKoDhoxFwiYy/0Rc2Y7cd7TO5Adbzvj1plbvaUNMvcyseZvVxVGjQM+FS5yzL3hYUZs9lD1FbmZlMHv3W8x0TCOe5DcspYDf7s182e3kF5nBcFrdjVOV4oU4I+T9w81u1l3izL3X68m71ik3KvgCeWcK6nOsfDav8IK/VvW1Ri+tM43tWAxHKeVCljVeYB3Wzaoq5v6ngvDoTNmNU5nqrMI7rPLO/95vVNNoeIeRE3R4I23rR46nMckKM3qHPcVuaeeMI5rlnXN/Mct0a8hkelxHVrfnOwb7+zqUMgGENRgfM4L/JHzNgD8lz3QJjlHVOpX7PM4wBCxmeqzG05O6vXs1Xm9ms7YJS9L6xilsdXE1PAZOCeJCIiIsozvANIREREGaULzezk09nlkDt4B5CIiIgoz/AOYDt4IwJeIayHwXts/4kkjByvB4B4UP4ckm0uCuPwFxltQIoKjcYavUJGo6SyYDNK/UbDjSLZYMmv6fDIhh5q7KMm3W/+B9Qsi089+SIe9yCmngQSMeb3NnnhO2B87jtgbIa/0Xj3NQI+2R5EtQ/yRHWzrYijPZsa+kQ9qNy2D6whQoz3eBCIF8j2JgXyiQCFsm1fUQQlsm1MmYy9V6AZxTL2AtkuyG8bbyIsV6ba7MRsbUBUu6BIzIt4WG5Ik1fGLuNu1ODbL2NuUvtAWG2iVOxqOJCYrY2Mag/lFRCa8/8l4bGeiGEOAyTLPF6gQxTJJ14UWmVeWmDE2TtobEhpQJa5N4KgbLRmL3MVu9Ic95nxq7Z/qh2cHvbC0yzLvVHFbrwAW/k3WXF7I3LYCrW7be2ihG24G/Vz4pAvug/QVcwhW5nL8g4WGW3dVNyloWb0Cho/F8vjvMgXdpQ3AESFF02y7V9Md7bttLd1VU++QbMXXhm7T5W7LHN/I+BrdB7nvmZb+7e47XhPaPunHu2i+2Br5yfnCdjKvVCWWaEc2qQohkJ5jqvYewWbrPL2GfulQDZc82i6eY6bxzsEIjJ21e4pqnvNp9/oMnaPOt4bPVZZq9ibbOd4U0I7uJhwlDcgy1ntDq91fXO0e4TtOA8JxOU5Dnm8B4oiKE48x2WZl/qbUSAbngZtbV2jcqeqdr561H6OO4/3WMQLyHNdHe++A5oVe8K7v0nA2yyfOhK1neOJ7Td9mvmUCc1s72mVf/I5LszrmzrHQ4VGbGWFzea1vVS2by7xN5ttXVNd3w7EgnIfyPK1DfkVsx3vmjrHmzzOOA/YYrdd1xPbOMeas/MsjTg85pN7OrsccgcrgERERJRRTAHnHlaliYiIiPIM7wB2gMpIxUIaooVyOIoiY1q0VL4XC8RLZDqoWA71UtyE3gVGWqBvgXGvvm/AyMuV+ptRIMfrSJUe2C/zDlHdi/1R4+e9zYUAgH0HCgAATfuDwH5jfn+DsZH+/cZL/QwYw9kARhpUpUY8MZmqtaU/reFuNNsTIOSQBQUqbg3RYhWz3C8lOkSpTAEWy1RQoZEC6lt4AH2CRux9ZC66yBdGoRqrRFJDIzTGA+Y+UMOB1EeCqG0yYm84YAyTEN0fgEfG7Kv3OOPdD/gPCBm7So3o8IZl7FE1roO1fpXyjMu0rwh4EA8mxF4MREtU7MaXVZl7iqMolrH3LpJlHWpE36BR3iVynIZir0r3W2Wu0mH740Hsl+NRqDKvDReg9oARe+N+Y5rYb+wXX4MH/gZj2/wNMvYDRjkbP8t0lZkO053pbsB44onfiDmuhkDxaYiFrPJWsZv7oNRYniiRw9gUR1BSZMTXRx7n5SEj7t6BRpTK2EMea6gfVd7Nsh1BNOY104H1MvZ9TcZxXtdYgEiDsV+0BmOeQIPHSnuax7mMu1HA3yjL2ixzkTzUj09DXMaujn1V5tFCzTrHZZnHigViJfK8ked4kVnmTSgvMGLuI8u8l7/JPM7t6U9AprxljrFJ7oPGeAC1YSPmfY3G+/79IeiyvL3qHLeVuV+l+dVx3qjD12TFbLzbxnAxU96yzH2a+aQbM/YiK/aYKvdSdbzH4C+Waf5i2byhoBHlIaPce6vrmyzzQm/EMbQTYFzj6mNGfOrpH/ujAdQ1G9Nq5fUtLI93NPjgb0hxju+3ylvFDhhDwHhkUwfHUD/q+hYwlqV7NCvmghTHuXl9k+dqSQwFsrx7FRnX9T4h63jvJa9v6hwPphjeqTEeQJM8x9X1bV/YuKbVNhWiYb8cBsYsc1/Kc9x4F47rG2Ac7+rabu72eHZSwDo80F245+TGMsjACiARERFlVFxo1iP8Orkccger0kRERER5hncA20F4jFesQPZKK9YQLjM+i8j3WC/Z87Msgr6lRgpgQHG98V5YjwHBOgBAhd94L/Ma86h0GABEZdfShngIu2NGvkmlSPbHgtjdZORj9jQY7+E6IzXi2+eDv87470iuBoE6gWC9TPnulynp/TL1Fo5Bi8p0gFBD6msQAWP98QLZI8/vM1NEKuUdKdXMuCNlsvdnL9kTrjSM8tL9MmYjT1FVsA8AUBmsQz+fMa3QTHlbqZFG2aV0X9xIc0Z1r5kC3hcx0iE79xejod7YH/o+I33ir/MioGI2djeCdXJU/gYd/gZjHd4mYxu1cBxaQipEeGUP6qAX8QIj5SJkKigW0hAtNn62l3m0TKZ8e8nef6VGKqiypAFVhcYGVcnC6B+oRy9Z3oXqcQVSVPjQKFOAu8wyD6FR9gzcI9P9u+uL0VxvTPPs88s4je0K1BnlDQDBepkOa4jB16hiN961iHoqvDB7uQq/7HkYtMpapQKjhZpZ3ir2qCzzWFkMwTIjlr6lRi6qqrgeVSEj5koZe7ks8xJvs6O8ASPtWyfL++uosYKY7sU+mf7c3Wgc5/vqjPd4XQC+fV4zZsAo86CMPdAQl7HL470xCq1Zxm4vc4+MUx3vhT5ApgNVKjBSpI536xyP9pIp77IoSsqM8u5fIo/zQnmuh+owIOA8x4s84aRe/Q26cUzvjRWbKeADMeN9T1MhdstzvFEe71qdD8F9xncD+6zYASBYr5sx+w/IMm+MWue4/fEdMMpcD8rj3KuZ76qJQ1iWebQUCPeSKV95vPtUmZc0oqpExlxgHe/lMj9Z4jH2j/36ptL8Dbos32gJonJ/1EeN/bGrsRh7641jIiavb946o8yD+zT4zXNcHe86/PtlzDJ2T1ie69HkMhd+L+KFxnbEg1aTh2ihfBKJbMoT6aXedYgy+TSbUlXm+zGgyNiQb8jY+8sN6+PbjxKPkfa2l/kBeX1Tx7uRAja2Y6863vcbueaGugJodeocl2VeZx3z6roeqJdlciAGrzzXPeoct5W5eX3zOc+/rsJOILmHdwCJiIiI8gzvABIREVFGCeEx73x3djnkDlYA20GXPeRUiiRSbKUIon2NdENBb+O2/zd678PBJXsAAIcU7gQADAvuxDd8tQCAXjIFWGh7EHpY3treI9NC27TeZorITJFEgmavX5X69e8xijG4R0Nor7Gs0F75MPa9UfjqZCqiXo6C3CjfI1EIIXvHyfQAgkFoxUZ6QvjkiaZp5sPRo7aUGACE++hAbyP92bevkfYdWlaLg4t2AQAODhmxDwoY+6LSW48SNRiq7JamA2iQad6dcSP9oQa+3YlSMw2qegXu318AvdaYFtzrke9AwR5jeaFaOSDtXmMfe+uboDXIEVKbjWkiZqVBNL+RZtEKjOWLkgLoQTnIttdKAUdkL0AzFdg3hkBfY99W9jJSPweX7gYAHFK4C8OCxj5QZd7Pe8BR3oBV5vv0IGrixoJVmkgXHrMHbJ3sBdpcH4R3r0wL7TG+a5b5Hh2hWjmQdq2tzFV5h41ycpR5QPamLZSxewrNbTN7fBdqiMier5Hessz6GMvq1ecABpftM2Iuscp8SMDYD5U+Y7/0Uul+CKg9cED2dt4VLzLXuVd2M22K+80UcP1+4z0u0/2BPV4EzePc2J6CvXEEao11+PYZ8drLXEStFCQAaD4fUGCcW54SebwHvGZ5Wz1gZdxlQKSvTLXJMq/o3YChpcaGDC9S57ixDwb596Cf10iJF8ke3h4AUajyNmKpiRknUrPuNxu3N8oUcG1TARrr5TaqMt+rIbRHxa7Ld3k+1YbhbZBlvV/GHolAyLS3pgYyD9rKXDYBiBfJNLhfc/RyB4Bwb4G4vL4V9TaWP7i3cUwPK96D4QVG7EMC6njfZ17fgpqzx2+j8GBP3CjP7bHeAIC6eIHZzKNO9oDdd6DATP1a1zcZd61AwR6Z/pRl7t/XDG2/jP2AEbu9zDWf/FMXkj2JiwqhyWYPqid0PKQhpsq7l3zvLc+VvmH0621c34aVGWV+SNFO6/rmNzauv9eYp8QTM9NrqswbdB92yeub6ukPGGlgANinrm8Nxj7Q9vkR3GNd3wCgYK9AaI/sbb9PXt/qbGUeltc3e5kH5PVNnuNaQXYqUHFoiKPz6Vs3lkEGVqWJiIiI8gzvABIREVFG6cKdDhwJw3dSJ7AC2A7xoAbNryFm3KFHrAiIFct0Y6mRiqgsM1Jew0t3YUTRDgDAYcHtAIChvn2okKnWYo9xO94jb8JGRQz1wrh9H4GxrIAWN1MF6jmRDeEQwk1yUND9ciDYeqvnr5n63S1TI3sOQKs1euSJeuNdbzJSBiIeB2RayCNTIx7ASpNIut8aIDVuxm6chaIkhlLZE3Jw6T4j9uKd+FbBVwCAbwa+BgBUyYFgyz0hBLWQY/lNIgyPjLled6bqosKLxpjsNdhsbFfsgA/+/VavOMBIC6nUSHCPsS7vHjlK7L566PuNdJwekQNOCx2az1iuR6YCzTSRCFkDAasBsEMaYjI7GiuWA8yWRtFP9nZWqd9vFRtlfWhwBw6WaaEKmVYs8xSY5a3LRGidbmxrHBHs0WNmzIAxIHBDxIi5sUmm7fb7zPI2Y5epwILdEfj3GOkvs8zr6qE3GesQMee+9QRD8BTLnJdMCwKArtKgKgUcgpkai8vBnotKjWUOLK3DoaVGGmxkwTYAwPBADQbJBw/38xrLLdCKzeWHhbEdtbpxvDfqEXhlqtAc+DwWxP6w8d1oozre5eC/9UBwn0z97jGO9+DuMHyqvOuMc1BvMH7Xm8NAQlMHT0GB1exBPg9XeK0Bz9XA13FV5iUCWqmx3X1k7/6hpXsxstg4x0eEjNgP9hvHQZVPR6kmU5ia7E0PHft11fvbOA73aVaZNNkGOgeApuYAcEAN7G6VeahWpn5VKnCXsT2e2nqIehmzPN7tZe5R6X5hFKYWsMrcGhRZgxxwwCrzYh3BUmO7q8qMg+6QYiPOEUXbbdc3Y7+Xe/wolE0J1PGuyrxOD6NZGLGr3rFR3Wc281DHe6TRD495fTO2Q41sULBHR2iXsT2+vXIU5No6q7zl8a7KHJrHPMc9MuWtFYQA20D3ABAL2mOXx4Qs895ljRhSaqS9R5SoMt+Og2Tae5B8vnGZRzZN0QrMfbtfnuNeRNGgO3vghnUf9kdl7+8mObC7LPNAvcfs4R2qlc1bdkcR3G1cbz17jXhFvTGT2H/Aur5Jms8PT5FRFprfWK4WCiAbdJfaALqxDDJwTxIRERHlGd4BJCIioozSoZmPwOvscsgdrAC2QzyoAQEN8rG8iIcAEZK90QIyNeY3bsEX+yIokj3hQnLgW7+tR1xcpifiMhXYJKJokIN27o0bXS73xIqxRz54dF/ESCkcCAegh41i8zUbN3DloybhDQt4I/J5n+rZlzE9+dmPMu2reWGmQbWgDCoYhCgwUgRqIOhYyGOlxEIqdrmeQByFAfmsY5+xISXeZjN29YxbrzxpdeiIipj82VhGgx5Fg0x77tML5T4w8k97o0Xm81DDMm6t2WvGrB4h7I0I81mn5sCvsqev6vVqhC4vHpofmuodp1LeMk2kFwYQKzC2R/WIjAWt2FWZ+wMxFAeMDSmWKSA1uHVIi8JrPnxT9qoUcXhkeYflPmjQ4zLuAuyRPWBrZdfT2kgh9oeNbYvL2L1hzVbeVuzGvtAB9dzPFM/7VClP8z0UNNP9qsz1Ah/ispegij0esspbC8h0c1CWuT+CYq+R4lJlHtDi8KqelXKA8SisMt8vU4D7dNUjtsAc/HqPjH1fOIQmmQJGk7G93mZNxmu8zJgBeKJxK+a4s6e1UeYyZpUGDQbN8hZy0O9Ygdcc5N08zuWhoQd1BILyHJdlXuJrNmNXAx37bb281fGtjvewiKFB/rxP9vTeK3uF7o6VoDZixK7KPNrsgyeceI4b57kRu+zpqXq0x3VzQHfzOPd6rfNdnuPmuR4KQA+p2GVTlAL79U0ev6E4QgGrvAHj+d2AUeYhTT3T2aKub1EY29goU8ANAtgjz+1dsgf0nmiROch7Y7Ms82avVd6265sRtw4tps7xhIHsYRvRQLY303x+6xyX73oogLh5jltlbpW7sTxv0Fh+YSCCYr96frdV5qGEQc2VqIghJmPfb5a534zdvL5FiszRDSLq+ibL3BNJjt0b0W0De6t3NYi/x4rd1rTHjF2Wu7quE/FIICIioozis4BzDyuARERElFHsBJJ7uCeJiIiI8gzvALZDtFCDHtDMdiK6X0DzOgcliuhGG4z6aAhf+4wnOwRlG5lm4UMv+XB0f0LbkQN6EHv1vgCAbVFjhPwvmsuxtdH4eVejfEJCsx+QbUQ01cxLboLu1cyhO1Q7D080ZBv6wGhr4rU/FcEn24yotjFFQcRKjZ/DvYz2QZESa3gE9UQQIeP2eIC4bMulRrSvjRVha7SvI74G3RiuodATgRfO4U7q9TLsihvtgb6K9AEAfN5UbuyLxjLUNxk7PNYsY4ppMJtayWyA7tOgq4e6F8l2XrJ9leb3Q1PDnaj2YZpmjpAvCmXbmGL51JWyACJlxrapJ5/EC4zytscOAJG4jEGODbQzasQR0qKIynZne/RGOc0q8wPySRD79F4AgO3R3tgSNvbZF03G+9eNxebwL2aZxzSzvNU/wmqomlihD1pMDedjbIcWDDrLGwC86okQQejFcmidEjn8Ri8/IiWyDaAcAiUesmLWPHIfyDRMc9yH2qgx43Zvbys+YQxNUSLbBaoyj8ODBt0oF9UGbHu0Nz5v6gcA2NZo7I+65gJEZXlrUXn82pr2qQcpqDKPFfnhkxvsUUO9FKmD1vZF+dQXURCELmOOlsnYS32IFFtD3wBWG0B4hXmsReU53hALYbdso6vaP6on2OyJN6LI4xySo1n4sE83Yq6JGnFuiRhlvbmxHF/Lc7xBHu96xAufOmRkCMJjDc8TD8k4S+R5DUBTT/lQj60Rwipv2f5RHe+x0iAivQIydtkertDW1tWv2pYJR3kD1vH+dbTMbOd7QDfKvNTTbF7f4vIeQ6Ns27szXoLttusbAHzV1Bt7G43Pw81+uZM16/omqTKPBz2IF8r2u/LC5PFaQ72Y7QIVvx+QQ5+oczxWFkC4zIglUqzavFrXN/jk8S7bbcd1jzlUjSrzQlv57vMYQ7LYyzxqtmvuBQCoifUyr2+b5fWtpqkE9XL4Fz0syzNmHe8q26nL3RILeeGR56x6loj5JKNwxGoLKa/5CAZs1zd5nBc528h2FR2aO+MAshOIa1gBJCIioowSLvUCFqwAuoYpYCIiIqI8wzuA7RAtBPQgoKs0gQcQceO/kUjE2JV7mwrN+RtkyqAmYKSCS31NKPQ6R8FXGuNB7Jf5pr1yOIi94ULsbpLDYhwwUh2xZh8QV0OZGG9qe2KFQFg36vTCK4d3KPbC2yxHwVdDw9iepWM9AUANheBBtFAOCVGkHgqvQT6/3UxFKPGoB/vlEzpqfEZqJKZ7sCdi7Icv/Eaqo8xnpL6DHisdqRrzNsRDZkppT1gOi9FsfH9PY5GZBhVR+f+KsKWDZOxG6s4rp8kUSbF8gkS0EFpUxq7SpxogZGosrlLHaiiMIo8Vu8ocB611Qqa8o1EvamV5e+SCm+XTHGoCpejlN2JWqSK/x0oBh+WO3C/zbfuiBdjVLGOXZV57wEqDImaVuZkOKlTLkmXu8yEmy87bS+6DiA4t7mymYJa534N4SKW/rDJPTIPqPpjHmogb8x2Qw3V87S1BTK5/n0wFbwn0QbEcv0INj6KOd114zKYC9TJ9tzdaiN3NRsy7G433hgMhiIhMiSWmvAPWsWmVuQe+IjlUjhwyKbHMjVjkk1gCHjNmtc8iRZpV3ioNqtL9AojJ7VFNErZpZYjIlOjXYSO128tvpPuLvWFHeQPGEy/UOa721R55ru9qKkLtAWNaxEx9W/+fWylAmCl64THmU8OZeMsC8Ibl0DC2MjfTiAFjPpU2jxZ4HOc4YKRB1brMbJ2uoTliTNx9oMgRU0M0hG2BXmbMAFDojTjKGwCa5ULrYyHz+rY7bLzvaSxCXaNRZrrcx1pcM29PqOubWebCA91nTPQVyzJvLoAnmnx9A4wnvMRt1zdjWR5Ei+TPMlseDwHCZ5U3AMSjVplv8xhlrJr57I4U43N/fxm7MTSM/foWlU+1aVTNPaIFSde3vY2FCKuhb2R5m9con3UcWmWumTH4So3veeT+1OK64/oGAMLvMWNXx3uz39k0oavowqUUMHsBu4YVQCIiIsoo9gLOPdyTRERERHmmR90BXLRoEX7729/isssuw9KlSwEAQgjccMMNeOCBB1BbW4tx48bh7rvvxre+9a12Lz8eAuBIBVqjtkc11QPW+H1/UxA7A8bt/qDPSMsEfDH4PcbPPo8zBRzTPWbvwkjMKJbmqM9MvcRkDzER85g90/SA8R4rNG6JC698WgmAaInx7ol4oDJRnpjs2WbvTalSBapXpU8zU0AqvRoP2nrHqUykTEOLJi+ahJHW+lqmS2oPFCAk00JB2Y0x6DXevR7dTJeqW/lx3YOwTKWFVewypR6J+Mzev5qu4hRmnCiR2+23nlzhlT14PVHVq06YPQrt6UD1j6TukylReTboAVvsZirQ+q4q85jHjzpZ3uoJBrvMMo8h6JdPgElR5iptqnqNNkf9CMe8MnajAKJhH0TM+T+a7heIyzjtKVEAiBZ74I3IJwDITJQWF2Z5J6WHPLbYVbx+o5kDYL0Lj7O8ASAM48NdUS/qGmVKVDZ1KPBH4ffGHbHbyzwuY1fHe3PMOs5VU4pos89sXmGmMOXTGWJCg1CZcbkvwqUavCqFpo73eIrjXZW511qGWdYBW+zqyqjKPOJBvNGY2CDbQzSH/dgjU6Ihv7HDA+pc98aSznFdaMnnuHoP+xGVsetR9SQLWNsom2BAs6VEi60nRgDG8e6JqXNFzi6c5W3EppoA2Mrddryr+cwyb/YiIq9ve+W2qZ7KO4KlCMlzPOCzjneP5kzDquM9Evchoo7zqFHm4YjPbOogVA9YzXoah2pyocokFtIQKZXXt6h1vKcqbxW3GbvsMR+3x67K3Na8RfU815uMlR7QNfPYVKn6UKCXGbs6zr2yzD2acFzfACAc96HJFjNgHO/xiLr4yv0t09Dq6Uv2bYsWavD2ktc11URCxR1PfY4Lr+ZYRiRLKVSmgHNPj7kD+M477+CBBx7At7/9bcf0xYsXY8mSJVi2bBneeecdVFZW4sQTT0RDQ0OWtpSIiIgou3pEBXD//v34+c9/jt///vfo3dsai0wIgaVLl2LevHk4/fTTMWrUKDzyyCNobGzEY489lsUtJiIiyh+6HAbGjRe5o0ekgC+++GKccsop+OEPf4ibb77ZnL5582bU1NRg4sSJ5rRgMIjjjz8ea9aswQUXXNC+FWkwqswyxeANa9Bkr0ghU2/igEyR+QNo9sleaapnmUfA45XTVNVbcw44CliD7Apdg1BpgbiVGoEcoDWuetqpgYBjGjwydaHSvlrcGjBapXTM9JAOM8VljzExLSw8MFO/ZsfGJhlvVIOuUqL7jcMp5hdoVDGbAwjLuL3CHEw4VcwwY4e5D8xpqkdmQCAml6ersV+jmjFIMmzpT11ts23gaHs6NEXsKl4zdvVsdR3whFUuUn4W9UHIgbQjPpna8RkbpPl0W+xy0GyvSFneKn6zrBPKHLAG5dU1HVGZwtMKE+KNaY5yN941Z3kjIW7nbgc8gO51TvPEADTLdarjXe6LuM+LJtkjU5X5Pp+wyluVtUorJqRFjXg95vi1ZsxCM7v/iqDxnZhKZQU1xGSPUM12vFsxy+2OW+Vl7u5UsXusd7O81SmrBqEWgCbTq6JRlrnfj7CMucGnDjYZr/04T+N41+NW+ZsDffsE9JAuf1aDP2tWijvxXNdt53jcmtbica6uZwn7wNx/YavMhTzH4z7Z81iWeZNfNwfD17xWmad7fQOMMhcyTaq+J4I64uYIBTJ2VRapjvOY1UQkqblHC2Vuxuy1Plb7VJWLuqaIJg+iMvaoLPMDPgF4reuait3cBynO8cTrmx732MpbXtOgevzq0EPqOJfbF7OVf6rreuK5bl+/ajKSnU7ATAHnoG5fAXz88cfxz3/+E++8807SZzU1NQCAiooKx/SKigp8+eWXLS4zHA4jHA6bv9fX17u0tURERETZ161TwFu3bsVll12GP//5zwiFQi3Op2nO/xiEEEnT7BYtWoSysjLzNWjQINe2mYiIKN+oO4BuvMgd3foO4Pr167Fz506MGTPGnBaPx/HGG29g2bJl2LRpEwDjTuCAAQPMeXbu3Jl0V9Bu7ty5mDNnjvl7fX29UQkUAHTAzGLFAI/HOSizGmRXeK2Bhs10ok/YeiE6U6TCa3uusFktt54/ar77dOshkMGEdIKuQVdpL9u7ShWoaVb6xJYetD1XOFWqUEtIuai0o/BoEKoybcZm9WSzUqnCijthGjy22M1z2/Z7QioNPushmbotjWRut+pJaKZGNNv2W2mixFSRI24Vuy2l4jVjt8rcXt72mITXmzQt7hWO8jZigrkPrA2B9W7OJ1ce0BzlDcBZ5olp/piW1ATAY0+VWo97dcRrfi7f1TGjUr/qWaPCa0sZe60yV+Wvq/hsx3vic4WN4yYxdgFNPYvW7JFrlbkwj1srdqtsE1KB8dSp0aTj3Fbu6nhRZY6IZqUMbbEnlbvcVuER0BOaQQgPbOd4wo63/02zl3kg+RyH7owTtmfHmnHajnfz/E9IDzqaQajNsJ0XZpl7NOv6JtPwVo9ij5VCNc9xOMrbvn80j3Be31Tsaj+Y+wxW7HJSzF6+iWWtW+napCYv8eSydqTGbWWedH2Tx7vQtKRmIcIrUl/f5PabMdvL3HZ8m9T1TY7qoJr4iFBycxAj9sT47Ne0FMd+wvUtYZzsLsMUcO7p1hXAH/zgB/jwww8d0375y1/isMMOw9VXX42DDjoIlZWVWLlyJUaPHg0AiEQiWLVqFW699dYWlxsMBhEMBlv8nIiIiKg769YVwJKSEowaNcoxraioCH379jWnz549GwsXLsTw4cMxfPhwLFy4EIWFhTjnnHOysclERER5h3cAc0+3rgCm46qrrkJTUxMuuugicyDoFStWoKSkJNubRkRERJQVPa4C+Prrrzt+1zQN8+fPx/z58zu9bE9ctp+ytx1JaKNnH0Ilub2MZrWN8yUMbeITECKhzYgX0DTncBrQzOZXJnMIDduwMbC1HdFjznaB5mj7MVt7QM3ZdsT+c8r2UuZMtthlnB5vcts49fQB4RNW+zB19PmsGMwhc1RTM/twCrbmZyKh/U6q2HWzPZTHbCdljpofsw0bo4a2scXdartA2/4XicNo2MpctQdTcQqf5ihv5zustp1Qw8bojpjNT20x298htKQyFrY4VXkK88kwVttIs32YLWbHUBKJsdvK3JMYu9dZ3sa7XI1XM8tYDXuhwWoXpspas7WXUrHby9wcSiRulbmjvAFbmdvaxKUoc3u8KdsFmjNaMat3a6gkFa/aF9bwTI4yF9bnAKCpMtdSDRtjrdoRu1nualwP63g3y9pW5sI2hIgxP8z3xOM8VbvAVOe4da5rSHzCiOYTKa9vxv4RxhBJgKPM7dc3Fbuw73skHO9J57iW3PbXPM5TlLtmxWe2kbTHDtt8Ml7zEbTmPtBs57ZzHwif1dYZtnkS2zobQ+Y4V+o43lOd6yr2mHryjbUP7MPiqGmJ17dsMZrQd/7uXZaaMPZIPa4CSERERLmFKeDc062HgSEiIiKi9uMdwI4QLfwMW3ZQt/8iP7ON1K6GdTCr4Lb0k/kwdluqS2vlxreZKrKlEsyUqsdKT5npLXP+5LRGqhHkW5VqX9hjT0gxCR1mntwxTW2HcKZI7Gl2e0rMTJOo3+2xqy+Y+9hKMZrpUs1KN6phPczFt/UPpi3mpKcN2PexmVpU22MNG5GYZhO2UWDMtJPQkp4mAKQubwAQcdt3zSdvaLbjSe1bNbSFbVmJ744Ft/Cz+j3FPtASj3M1dIkmbOnMhEJM+CWxqYO9zM3zQZUrNGv/eRI+0zXr+FJPF/FozrKS760WfWKzgITvArDOb48GTU+xHfZjsh0cx7u5TufyYTuPhP18VunpxHPcHqz9fEtVxmq2xOuDZhWj/Wk7Sdc3te8ErC+kWoF90Ymz2ZqCmOvUrNg0+36Q8wHO41ykKGtzWak2I0WZ2zffOs7VNU1+QbfOO+v8EEg6V+zhpTjezdnsx1ni+W8vc/McUOcb0Ob1rIvwDmDuYQWQiIiIMooVwNzDFDARERFRnuEdwHbQvUavPS3x1j7gSDcAcKRXzVRAa/+4CNv9/hQpGl2oh6WLFjMowr4M2wPXU6Y2ErbbkTpM0RPOnpa2rzPpu+o9cbmpUoy2ZZmpFNWDz/xIs3pH61YKKHEZ9getm/vPvq2pYrKVFQDHQ+HNxdvTQ+ksw2N9zyp3W1rTvr32bU3RZMDRC7CV1Ky9F7CVikqVXrW229hmYTYPMOez7QNHmbfUO1KD9W+k/dhv7VhP6GXsSNHJhQkhUpe33FZz/9n2o/l0kDSOd6Eln8fCa2aqU/d6T3Es28vb8VmqVQtbCwezl7616MTe3ZomUseuJ8Zu+z3FhljpUlv6u5XY0r2+GcuyPkuaJ2G71bbaYzY+8ljnuNksI/m7jjJXx449RdvCMeq4ftliM893+/a3dOzYyzrlOZ7iK2Z5qm3UrLKwheYo78RlpYw9Rbm3QKQ4x0WWbvvwDmDu4R1AIiIiojzDO4BERESUUUJo1p3MTi6H3MEKYDsYD7m3ZQlaSa8JW+8rx8CxiT29bGlKM1UA65a9vTcnAAg9RVrIvj1m6ld+T7enDBK3VZg9YK1eifYubtZ2W6nLFOtNTF3a0oJWL2NbeiNxGbpmLsNMjQhru8yefq2kSAAkDQSt2XugJu14JPU8VoOUemDNrlJGImFdjuXYNydFWadKD9rTnwCMwWsTUjSaRzPKO8V6rI2Cs8yT0oLJKVFHGi8xDRZPSHuj7TJPmRb0tLbdmrkuNYu9vAGjd2fK8jaXId/scSalxlJsry39rcpbnQKaLQ2f1Fvbzh57YlrQfLcfmGqhwnYOao5FCWFtiOrRLpK6wcJ5jtvSqub2J6TBHen1NFK7mqeF61uqcxxoualLwjFkDkyvCfP6Zl9Pq+d44vYIzXl9UxvUQg9bo1mLPMe9Vpkn9oYW3uTvmlKUuT0tjNauTWobbQNvm8eXrrVd3io+FUbS9c05j2M7NFiDcav9nqKpS1fQobkyELQbyyADU8BERETUo91zzz0YNmwYQqEQxowZgzfffLPV+VetWoUxY8YgFArhoIMOwn333Zc0z5NPPomRI0ciGAxi5MiRWL58uePzN954A6eddhqqqqqgaRqefvppN0PqNFYAiYiIKKNUJxA3Xu31xBNPYPbs2Zg3bx42bNiAY489FieddBK2bNmScv7Nmzfj5JNPxrHHHosNGzbgt7/9LWbNmoUnn3zSnGft2rU488wzMW3aNLz//vuYNm0azjjjDLz11lvmPAcOHMARRxyBZcuWtX+HdQFNiJQJHrKpr69HWVkZDrpuITyhUOqesEorKWBoIjldYk8Zabaf1WcJqTTHsZ8it5e0bbYec4mpA62FdEJiD1JNT9Gzso30UGLvOHtMiWlhx/wJg5vC1hPS6n2ZIr1mn5yYekkx+LRmG5Q36dm3QkvqBeqIv42eoebv5qDDsN4TU4T2/ZPUmzZFj+8UsTviTtETOrk85ffiSEoZQk8+hlrrBWzbjITUmLMc7eWbNL8jdmu/pCzvhJU6tjUx9Ws7thNTkY50qX1aqvK2L1Num7kJieVvK/O0z3G1TNv5YK2rpfx9ckzyYauOmJypcfWZ7V1P3o+t9fRP2gzbuWs/ps3Yvc6YUjYPsO+XNGJ3Pq/YFluq6xuM6YnpUk23UsZJ8wNIdby3eo630gxCpCj/jlzfzI+SjvNUMVn7JfH6pjeGsXn+b1FXV4fS0lJkmvr7+d3ll8FXFOz08mIHwnj7x/+vXds/btw4fOc738G9995rThsxYgSmTp2KRYsWJc1/9dVX49lnn8XGjRvNaRdeeCHef/99rF27FgBw5plnor6+Hi+99JI5z+TJk9G7d2/89a9/TVqmpmlYvnw5pk6dmm6oGcc7gERERNSt1NfXO17hcDjlfJFIBOvXr8fEiRMd0ydOnIg1a9ak/M7atWuT5p80aRLeffddRKPRVudpaZm5iBVAIiIiyii3U8CDBg1CWVmZ+Up1Jw8Adu/ejXg8joqKCsf0iooK1NTUpPxOTU1NyvljsRh2797d6jwtLTMXsRdwOwivfKWbNE+4pW9PGSXNKqw7+uYgsRC27nnyzfH9VtpCONKCLcySooegBltPPHODbE8iTid2Dc6YW9pUW3dTs1esmQ9RuQwtRewtxN1S2k6kGAxbE+YCkwZztfe0Tuxx2RZbmbcWe9Iu1jVrHSpFJjQkh9xGmadqAtBS+XtsYZobIpCYXmsz9nYc58ZBZ8znKPOE5xRDS6O8HctFUuwpB4S2pdusHt62lGQnylstN+W2oYVz3LHfNccijR/SPMfVe6oN0Jz71jEQcOKBKLT0Yk+Vqk11vCd0KddszTGE7Tjv1PVNvrc4ALgtfWuWuf2ao/aLfXmtaeUcb7X8bQNBd+j6ppaV6hxP2kZr+YkDP1u9gruW28PAbN261ZECDgZbTy9rCb2shRBJ09qaP3F6e5eZa1gBJCIiom6ltLQ0rTaA5eXl8Hq9SXfmdu7cmXQHT6msrEw5v8/nQ9++fVudp6Vl5iKmgImIiCijhEvp3/beRQwEAhgzZgxWrlzpmL5y5UqMHz8+5Xeqq6uT5l+xYgXGjh0Lv9/f6jwtLTMX8Q4gERER9Vhz5szBtGnTMHbsWFRXV+OBBx7Ali1bcOGFFwIA5s6di23btuHRRx8FYPT4XbZsGebMmYOZM2di7dq1ePDBBx29ey+77DIcd9xxuPXWWzFlyhQ888wzePXVV7F69Wpznv379+Ozzz4zf9+8eTPee+899OnTB4MHD+6i6FvGCmA7GA/WttpPtPrQebTR/q3VL8qv2duM2FeWzvIcQxokDPVgZ07UzFkTFy+8VuOltmK2rc65ra20kbE1f7KNzm9tT9JK24o/oW1MymFMUizH3BMeW9sskX5521bZ8naKhA9tbbDM5kG67YutPRmhxWUj7TJPfDqEJpJjT/s4N5ePNsvb/pnxBA7NOX9Hjvf23BnQnOUNyHATGsB1KHb7e8LMqjWtSJhNQDOfjJJ2mSe1eU0xDFSqbVOzpHrCT0eO91TneMK6HGWeMFRJp85xx4aknqTB1t7XtnirbZztHOzk9S3VjPYyt069FLF39hxvo8wThyXqailOrw4vp73OPPNM7NmzBzfeeCN27NiBUaNG4cUXX8SQIUMAADt27HCMCThs2DC8+OKLuPzyy3H33XejqqoKd955J37yk5+Y84wfPx6PP/44rrnmGlx77bU4+OCD8cQTT2DcuHHmPO+++y5OOOEE8/c5c+YAAKZPn46HH364A5G4i+MApkGNYzT0pgXwhELm9A5VABP/MNrHCEzxxyQjFUD7ohIbEustXFDcrAAm7pcWYrfmd78CqNkrDKk6D7S2r1pbdTtjdzwyLFWlyY0/DqliV/OkEbsrFUD7e6oOE0nzd74CmHpsO1uZJ8bX2r5qabWtVADTir0zZZ5mBTBl7PI9edzQFN9rbRPaqgAmlGeLZe5CBbClilCqsfNS/lPoQgXQUa7tjT1T1/WEz0VTM76YN6/LxwE84n+vgLew8+MAxhvDeP+nd3TZ9vdkbANIRERElGeYAm4P9SSHhBRWm+y34RP/M7alBM1Uge2zlCmDdt6zbTVFkGJbRcK2acLa3nbfL05zHyWmBe0jc5jDRiTts/SW2fYMzo0Utkn2fdeutsepyryVTXCkxhx5qoTY0119e8pcrUvOriVOc7vME1Nj9tht83S6vG2/O9KBCRtpOwWT9lu7R61IdUevNSnOsU6XeZszJq3eWeaduL61yFbmqa8vLpzjrX0nYRtTnuNax4/1dh8nrV3f0l11WtdzxxndofW4xe1hYKjzWAEkIiKijNKFBs2FyltHngVMqTEFTERERJRneAewHYRKjyamDtrx/XQkpUOB1KmitFfcjpU6ukcmf73dq043PaSeFNBCKjhxE9u3Ea2u0raCVr6aav62VttSB4EUMzlSY0krb3/I6Zd5QiFrKb7aztg7ksrKdHmnXnErX3XjHG+jvJPmb2c/iORlt/F5UvfrFr7q5vUtRRm2eX1LmD+9lbfxeSsLdvX61plzvL2xd7RpRtZSwB1IsbewHHIH7wASERER5RneASQiIqKMYieQ3MMKYHvYx3ZC++/Am8toiXDmXhypMcd8HVlxO7SyjZ1edUvLThG7fbK7G5G4TSlSYynSMR1ebWdS/51eeZorbaPQO3zNTZXySujtCk10feyO0cdTbFsGV524Gea6WkkPZmSlQsvMud7WNS4b53git8/xdNgGArdNSjGfK6sCoFp5tJz670qsAOYepoCJiIiI8gzvABIREVFGcRiY3MMKYHuoFHBrWbNUt+872rXMlhpL/ChjUuWd24pZaW/saSyzxTR4JrSWGmup93VL6ZqOlHlCakxN7nKJPXLdLPOUy8hS7K2lxjIde0KTh1Sb45YUq0qWohd2V1zfErerS8s83XPc/nnS8jq/OV1+fcsS9gLOPUwBExEREeUZ3gEkIiKijDLuALrRCcSFjSEAvANIRERElHe6dQVw0aJFOOqoo1BSUoL+/ftj6tSp2LRpk2MeIQTmz5+PqqoqFBQUYMKECfj44487t+KE4WBSfqa1MV86Uvy3pAnrlTGpVpBOTO2NPdU85uNWkjcn43HbV5g0Dam3N1XMHS3zbMbeVpm3+L0UrwxuUkZ01fGeShtl7oYWl9We4zzV525c3zIce4uyFXsuXN+yRA0D48aL3NGtK4CrVq3CxRdfjHXr1mHlypWIxWKYOHEiDhw4YM6zePFiLFmyBMuWLcM777yDyspKnHjiiWhoaMjilhMREeUP4eKL3NGt2wC+/PLLjt8feugh9O/fH+vXr8dxxx0HIQSWLl2KefPm4fTTTwcAPPLII6ioqMBjjz2GCy64IBubTURERJRV3foOYKK6ujoAQJ8+fQAAmzdvRk1NDSZOnGjOEwwGcfzxx2PNmjWdX2EG0l5JUqQMzNVnOnXQ2oIzHXuacXdUm9/twbG3Kp2YOxN3S9/N5nFuX1GLnyHzZZ5qtd39HE/n+7keeybKvAuuby1+P0t5ZqaAc0+3vgNoJ4TAnDlzcMwxx2DUqFEAgJqaGgBARUWFY96Kigp8+eWXLS4rHA4jHA6bv9fX12dgi4mIiPKEW/lb5oBd0+kKYF1dHZYvX44333wTX3zxBRobG9GvXz+MHj0akyZNwvjx493YzjZdcskl+OCDD7B69eqkzzTN+R+DECJpmt2iRYtwww03uL6NRERERLmgwyngHTt2YObMmRgwYABuvPFGHDhwAEceeSR+8IMfYODAgXjttddw4oknYuTIkXjiiSfc3OYkl156KZ599lm89tprGDhwoDm9srISgHUnUNm5c2fSXUG7uXPnoq6uznxt3bq17Y3oaA+xdOdtJWVgLqqdqYO0MwFtzehG79+WtCPmrPSia2+ZZzH2tPdPWgtL8Upn/ra08zhPN6aMxZ6OLMbernM83RR4VzR9SbUJHbi+ZfV4T0e2zvFscCv9yxSwazp8B/CII47Aueeei7fffttMuSZqamrC008/jSVLlmDr1q248sorO7yhqQghcOmll2L58uV4/fXXMWzYMMfnw4YNQ2VlJVauXInRo0cDACKRCFatWoVbb721xeUGg0EEg0FXt5WIiChf8VFwuafDFcCPP/4Y/fr1a3WegoICnH322Tj77LOxa9eujq6qRRdffDEee+wxPPPMMygpKTHv9JWVlaGgoACapmH27NlYuHAhhg8fjuHDh2PhwoUoLCzEOeec4/r2EBEREXUHHa4AtlX56+z86bj33nsBABMmTHBMf+ihhzBjxgwAwFVXXYWmpiZcdNFFqK2txbhx47BixQqUlJS0f4Xq/nq6t6CzeKfa9TSAWmA2Yk/rafa2VbsZezbjtq8332IXWruCyVraK1M9RIHslXl7UmxuxZ/NmO0LzdfYu4BbPXjZC9g9rvQCfvbZZ1NO1zQNoVAIhxxySFJ61g0ijXvBmqZh/vz5mD9/vuvrJyIiIuqOXKkATp06FZqmJVXI1DRN03DMMcfg6aefRu/evd1YJREREXUXbnXg4B1A17gyEPTKlStx1FFHYeXKlWbP2ZUrV+K73/0unn/+ebzxxhvYs2eP651AsiabXa3y+eDPZg+w7pp3cUM2j/Vslne+lnm2r2/5eo3r4bGrTiBuvMgdrtwBvOyyy/DAAw84xvz7wQ9+gFAohF//+tf4+OOPsXTpUpx33nlurI6IiIiIOsGVCuB//vMflJaWJk0vLS3F559/DgAYPnw4du/e7cbqiIiIqDvhk0Byjisp4DFjxuC///u/HUO97Nq1C1dddRWOOuooAMC///1vxyDN1AlMjWUHU2P5J9tlnq94fetx+Czg3OPKHcAHH3wQU6ZMwcCBAzFo0CBomoYtW7bgoIMOwjPPPAMA2L9/P6699lo3VkdEREREneBKBfDQQw/Fxo0b8corr+Bf//oXhBA47LDDcOKJJ8LjMW4yTp061Y1VERERUXfEm6s5xZUKIGAM+TJ58mRMmDABwWAQmsbbtERERES5yJU2gLqu46abbsI3vvENFBcXY/PmzQCAa6+9Fg8++KAbq8hNbCuSHWwD0vWy3RaOZd712B4uO3ro8c42gLnHlQrgzTffjIcffhiLFy9GIBAwpx9++OH4wx/+4MYqiIiIqLsSLr7IFa5UAB999FE88MAD+PnPfw6v12tO//a3v41PP/3UjVUQERERkUtcaQO4bds2HHLIIUnTdV1HNBp1YxWUSjsfIu4qtU7eju9aQmNaLhs0kZ1jPZvneD7j9S0DNPlyYznkBlfuAH7rW9/Cm2++mTT9b3/7G0aPHu3GKoiIiKi7Ygo457hyB/D666/HtGnTsG3bNui6jqeeegqbNm3Co48+iueff96NVRARERGRS1y5A3jaaafhiSeewIsvvghN03Dddddh48aNeO6553DiiSe6sYrcxZ5y2cGnoeSfHto7Muex9ze5gXcAc45r4wBOmjQJkyZNcmtxRERERJQhrlUAiYiIiFJy644u7wq7psMVwN69e6f9tI+9e/d2dDVElEvYO5Koa/SwHv9CGC83lkPu6HAFcOnSpebPe/bswc0334xJkyahuroaALB27Vq88soruPbaazu9kURERETkng5XAKdPn27+/JOf/AQ33ngjLrnkEnParFmzsGzZMrz66qu4/PLLO7eVRERE1H251YGDdwBd40ov4FdeeQWTJ09Omj5p0iS8+uqrbqyCclUPSlF0G+wdmX9Y5tnBHv/uUcewGy9yhSsVwL59+2L58uVJ059++mn07dvXjVUQERERkUtc6QV8ww034Pzzz8frr79utgFct24dXn75ZfzhD39wYxVERETUTbl1M5U3ZN3jSgVwxowZGDFiBO6880489dRTEEJg5MiR+Mc//oFx48a5sQpqDZ8XSpR57AFN1HFsA5hzXBsHcNy4cfjLX/7i1uKIiIiIKEM63AbwwIEDGZ2fiIiIegh2Ask5Ha4AHnLIIVi4cCG2b9/e4jxCCKxcuRInnXQS7rzzzo6uioiIiIhc1OEU8Ouvv45rrrkGN9xwA4488kiMHTsWVVVVCIVCqK2txSeffIK1a9fC7/dj7ty5+PWvf+3mdhNlnyb43yhRprGNc8/ANoA5p8MVwEMPPRR/+9vf8NVXX+Fvf/sb3njjDaxZswZNTU0oLy/H6NGj8fvf/x4nn3wyPB5XRpshIiKi7ogVwJzT6U4gAwcOxOWXX86nfRARERF1E671As57HCKCqGsIjelAou6GdwBzDiuARERElFlu9eDlTRbXsHEeERERUZ7hHUAiIiLKKD4KLvfwDiARERFRnnGtAvjmm2/iF7/4Baqrq7Ft2zYAwJ/+9CesXr3arVV0yj333INhw4YhFAphzJgxePPNN7O9SURERPlBuPjqgPbWAVatWoUxY8YgFArhoIMOwn333Zc0z5NPPomRI0ciGAxi5MiRWL58eafX25VcqQA++eSTmDRpEgoKCrBhwwaEw2EAQENDAxYuXOjGKjrliSeewOzZszFv3jxs2LABxx57LE466SRs2bIl25tGncVHAxFlnlv5O6IsaG8dYPPmzTj55JNx7LHHYsOGDfjtb3+LWbNm4cknnzTnWbt2Lc4880xMmzYN77//PqZNm4YzzjgDb731VofX29U0IUSnz+rRo0fj8ssvx7nnnouSkhK8//77OOigg/Dee+9h8uTJqKmpcWNbO2zcuHH4zne+g3vvvdecNmLECEydOhWLFi1q8/v19fUoKyvDkFtvhicUan3mbFZGsnWBzoUKWD7Gns0/yNku83yNPduVsHyNPVtxZyBmvbkZX159Derq6lBaWur68hOpv5+Db70ZnoI2/n6mQW9qxpZ2bn976wBXX301nn32WWzcuNGcduGFF+L999/H2rVrAQBnnnkm6uvr8dJLL5nzTJ48Gb1798Zf//rXDq23q7lyB3DTpk047rjjkqaXlpZi3759bqyiwyKRCNavX4+JEyc6pk+cOBFr1qxJ+Z1wOIz6+nrHi4iIiDpGg3UjuVMvubzEv9Eq85ioI3WAtWvXJs0/adIkvPvuu4hGo63Oo5bZkfV2NVcqgAMGDMBnn32WNH316tU46KCD3FhFh+3evRvxeBwVFRWO6RUVFS3emVy0aBHKysrM16BBg7piU4mIiHom1VzHjReAQYMGOf5Ot3RHrSN1gJqampTzx2Ix7N69u9V51DI7st6u5sowMBdccAEuu+wy/PGPf4Smadi+fTvWrl2LK6+8Etddd50bq+g0TXPewhdCJE1T5s6dizlz5pi/19fXsxJIRESUI7Zu3epIAQeDwVbnb08doKX5E6ens8z2rrcruVIBvOqqq1BXV4cTTjgBzc3NOO644xAMBnHllVfikksucWMVHVZeXg6v15tU4965c2dSzVwJBoNtHkxERESUJpcfBVdaWppWG8CO1AEqKytTzu/z+dC3b99W51HL7Mh6u5prw8AsWLAAu3fvxttvv41169Zh165duOmmm9xafIcFAgGMGTMGK1eudExfuXIlxo8fn6WtIiIiokzrSB2guro6af4VK1Zg7Nix8Pv9rc6jltkd6h6uPgmksLAQY8eOdXORrpgzZw6mTZuGsWPHorq6Gg888AC2bNmCCy+8MNubRkRE1PO5fAewPdqqA8ydOxfbtm3Do48+CsDo8bts2TLMmTMHM2fOxNq1a/Hggw+avXsB4LLLLsNxxx2HW2+9FVOmTMEzzzyDV1991TH2ca7XPTpcATz99NPTnvepp57q6GpcceaZZ2LPnj248cYbsWPHDowaNQovvvgihgwZktXtIiIiygfZfBRcW3WAHTt2OMbmGzZsGF588UVcfvnluPvuu1FVVYU777wTP/nJT8x5xo8fj8cffxzXXHMNrr32Whx88MF44oknMG7cuLTXm20dHgfwl7/8pfmzEALLly9HWVmZeQdw/fr12LdvH04//XQ89NBD7mxtlnAcwDZke0w4ID9jz8dx0ZR8jZ3jAGYHxwHsMPX3c+iCBW3//UyD3tyML+bN67Lt78k6fAfQXqm7+uqrccYZZ+C+++6D1+sFAMTjcVx00UUsIKKeSP1hynZFkKgnU+dXtiv+bshiCphSc6UTyB//+EdceeWVZuUPALxeL+bMmYM//vGPbqyCiIiIuqssPwuYkrlSAYzFYo5HpigbN26EruturIKIiIiIXOJKL+Bf/vKXOO+88/DZZ5/he9/7HgBg3bp1uOWWWxxtBYmIiCj/ZLMTCKXmSgXw9ttvR2VlJX73u99hx44dAIzHw1111VW44oor3FgFERERdVe2x7h1ejnkClcqgB6PB1dddRWuuuoq1NfXAwA7fxARERHlKFcHggZY8SMiIqIE7AWcc1ypAA4bNqzVhxt//vnnbqyGiIiIiFzgSgVw9uzZjt+j0Sg2bNiAl19+Gf/93//txiqIiIiom2InkNzjSgXwsssuSzn97rvvxrvvvuvGKoiIiKi7Ygo457gyDmBLTjrpJDz55JOZXAURERERtZPrnUDs/vd//xd9+vTJ5CqIiIgo17mUAuYdQPe4UgEcPXq0oxOIEAI1NTXYtWsX7rnnHjdWQURERN0VU8A5x5UK4JQpUxwVQI/Hg379+mHChAk47LDD3FgFEREREbnElQrg/Pnz3VgMERER9US8A5hzXOkE4vV6sXPnzqTpe/bsgdfrdWMVREREROQSV+4ACpG6Sh4OhxEIBNxYBRHlEj6PkyjzetCgdxwHMPd0qgJ45513AgA0TcMf/vAHFBcXm5/F43G88cYbbANIRERElGM6VQH83e9+B8C4A3jfffc50r2BQABDhw7Ffffd17ktJCIiIiJXdaoCuHnzZgDACSecgKeeegq9e/d2ZaOIiIioB2EnkJzjShvA1157zY3FEBERUQ/ENoC5p8MVwDlz5uCmm25CUVER5syZ0+q8S5Ys6ehqiIiIiMhlHa4AbtiwAdFoFADwz3/+0zEQNBEREZED797llA5XAO1p39dff92NbSEiIiKiLuDKQNDnnXceGhoakqYfOHAA5513nhurICIiou5KuPgiV7hSAXzkkUfQ1NSUNL2pqQmPPvqoG6sgIiKibkp1AnHjRe7oVC/g+vp6CCEghEBDQwNCoZD5WTwex4svvoj+/ft3eiOJiIiIyD2dqgD26tULmqZB0zR885vfTPpc0zTccMMNnVkFERERdXccBzDndKoC+Nprr0EIge9///t48skn0adPH/OzQCCAIUOGoKqqqtMbSURERN0XxwHMPZ2qAB5//PEAjCeCDBo0CB6PK00KiYiIiCiDXHkSyJAhQwAAjY2N2LJlCyKRiOPzb3/7226shoiIiLojpoBzjisVwF27duGXv/wlXnrppZSfx+NxN1ZDRERE3RErgDnHlZzt7NmzUVtbi3Xr1qGgoAAvv/wyHnnkEQwfPhzPPvusG6sgSo3jAhBlntCMFxH1GK7cAfz73/+OZ555BkcddRQ8Hg+GDBmCE088EaWlpVi0aBFOOeUUN1ZDRERE3RA7geQeV+4AHjhwwBzvr0+fPti1axcA4PDDD8c///lPN1aR5IsvvsD555+PYcOGoaCgAAcffDCuv/76pPaHW7ZswWmnnYaioiKUl5dj1qxZSfMQERER5RNX7gAeeuih2LRpE4YOHYojjzwS999/P4YOHYr77rsPAwYMcGMVST799FPouo77778fhxxyCD766CPMnDkTBw4cwO233w7AaHt4yimnoF+/fli9ejX27NmD6dOnQwiBu+66KyPbRURERAnYBjDnuFIBnD17Nnbs2AEAuP766zFp0iT85S9/QSAQwMMPP+zGKpJMnjwZkydPNn8/6KCDsGnTJtx7771mBXDFihX45JNPsHXrVnM8wjvuuAMzZszAggULUFpampFtIyIiIhtWAHOOKxXAn//85+bPo0ePxhdffIFPP/0UgwcPRnl5uRurSEtdXZ1jMOq1a9di1KhRjsGoJ02ahHA4jPXr1+OEE07osm0jIiIiyhWuVAATFRYW4jvf+U4mFt2i//znP7jrrrtwxx13mNNqampQUVHhmK93794IBAKoqalpcVnhcBjhcNj8vb6+3v0NJiIiyhPsBJJ7OlwBnDNnTtrzLlmyJO1558+f3+bzg9955x2MHTvW/H379u2YPHkyfvazn+FXv/qVY15NSx66QAiRcrqyaNEiPsOYiIjILUwB55wOVwA3bNiQ1nytVbRSueSSS3DWWWe1Os/QoUPNn7dv344TTjgB1dXVeOCBBxzzVVZW4q233nJMq62tRTQaTbozaDd37lxHBbe+vh6DBg1qRxREREREuavDFcDXXnvNze0wlZeXp91ucNu2bTjhhBMwZswYPPTQQ0nPIq6ursaCBQuwY8cOszfyihUrEAwGMWbMmBaXGwwGEQwGOx4EERERmZgCzj0ZaQPYFbZv344JEyZg8ODBuP32282xBwHjzh8ATJw4ESNHjsS0adNw2223Ye/evbjyyisxc+ZM9gAmIiKivNVtK4ArVqzAZ599hs8++wwDBw50fCaE8S+C1+vFCy+8gIsuughHH300CgoKcM4555jDxBAREVEXYBvAnNNtK4AzZszAjBkz2pxv8ODBeP755zO/QXxOJlHXYA6IqPthBTDnuPIoOCIiIiLqPrrtHUAiIiLqHjT5cmM55A5WAImIiCizmALOOUwBExEREeUZ3gEkIiKijOI4gLmHdwCJiIiI8gzvABIREVFmsQ1gzmEFkIiIiDKPlbecwhQwERERUZ7hHUCijuLTX4gyj63+ewR2Ask9rAASERFRZrENYM5hCpiIiIgoz/AOIBEREWUUU8C5h3cAiYiIiPIM7wASERFRZrENYM7hHUAiIiLKKJUCduOVKbW1tZg2bRrKyspQVlaGadOmYd++fa1+RwiB+fPno6qqCgUFBZgwYQI+/vhjxzzhcBiXXnopysvLUVRUhB/96Ef46quvHPMsWLAA48ePR2FhIXr16uVyZKmxAkhERER575xzzsF7772Hl19+GS+//DLee+89TJs2rdXvLF68GEuWLMGyZcvwzjvvoLKyEieeeCIaGhrMeWbPno3ly5fj8ccfx+rVq7F//36ceuqpiMfj5jyRSAQ/+9nP8F//9V8Ziy8RU8BERESUWTmeAt64cSNefvllrFu3DuPGjQMA/P73v0d1dTU2bdqEQw89NHlThMDSpUsxb948nH766QCARx55BBUVFXjsscdwwQUXoK6uDg8++CD+9Kc/4Yc//CEA4M9//jMGDRqEV199FZMmTQIA3HDDDQCAhx9+ODMBpsA7gD1Bpu+LE5Ex8DcH/ybqGOHiC0B9fb3jFQ6HO7V5a9euRVlZmVn5A4Dvfe97KCsrw5o1a1J+Z/PmzaipqcHEiRPNacFgEMcff7z5nfXr1yMajTrmqaqqwqhRo1pcbldhBZCIiIi6lUGDBplt9crKyrBo0aJOLa+mpgb9+/dPmt6/f3/U1NS0+B0AqKiocEyvqKgwP6upqUEgEEDv3r1bnCdbmAImIiKijHJ7HMCtW7eitLTUnB4MBlPOP3/+fDO92pJ33nnHWLaWfIdfCJFyumObEj5P5zvpzJNprAASERFRZrncBrC0tNRRAWzJJZdcgrPOOqvVeYYOHYoPPvgAX3/9ddJnu3btSrrDp1RWVgIw7vINGDDAnL5z507zO5WVlYhEIqitrXXcBdy5cyfGjx/f5vZnElPARERE1COVl5fjsMMOa/UVCoVQXV2Nuro6vP322+Z333rrLdTV1bVYURs2bBgqKyuxcuVKc1okEsGqVavM74wZMwZ+v98xz44dO/DRRx9lvQLIO4BERESUUZoQ0ETnbwG6sYxURowYgcmTJ2PmzJm4//77AQC//vWvceqppzp6AB922GFYtGgRfvzjH0PTNMyePRsLFy7E8OHDMXz4cCxcuBCFhYU455xzAABlZWU4//zzccUVV6Bv377o06cPrrzyShx++OFmr2AA2LJlC/bu3YstW7YgHo/jvffeAwAccsghKC4uzkjMrAASERFR3vvLX/6CWbNmmT12f/SjH2HZsmWOeTZt2oS6ujrz96uuugpNTU246KKLUFtbi3HjxmHFihUoKSkx5/nd734Hn8+HM844A01NTfjBD36Ahx9+GF6v15znuuuuwyOPPGL+Pnr0aADAa6+9hgkTJmQiXGhCZKg63YPU19ejrKwMQ269GZ5QKPVM2RweIttDwORr7NkeEiRfY89W3NkubyA/Y8/29Q3IXvwZiF1vbsaXV1+Durq6tNrQdZb6+3nkLxbAG2jh72c7xCPNeO/P87ps+3sy3gEkIiKijHK7FzB1HjuBEBEREeUZ3gEk6m74LzBR18iFtH9PkeOPgstHrAASERFRRjEFnHuYAiYiIiLKM7wDSERERJnFFHDO4R1AIiIiojzDO4BERESUUWwDmHtYASQiIqLMYgo45/SIFHA4HMaRRx4JTdPM5+cpW7ZswWmnnYaioiKUl5dj1qxZiEQi2dlQIiIiohzQI+4AXnXVVaiqqsL777/vmB6Px3HKKaegX79+WL16Nfbs2YPp06dDCIG77rorS1tLRESUf5i+zS3dvgL40ksvYcWKFXjyySfx0ksvOT5bsWIFPvnkE2zduhVVVVUAgDvuuAMzZszAggUL+BxBN3CgVKLM419O6u6EMF5uLIdc0a1TwF9//TVmzpyJP/3pTygsLEz6fO3atRg1apRZ+QOASZMmIRwOY/369S0uNxwOo76+3vEiIiIi6im6bQVQCIEZM2bgwgsvxNixY1POU1NTg4qKCse03r17IxAIoKampsVlL1q0CGVlZeZr0KBBrm47ERFRPlG9gN14kTtyrgI4f/58aJrW6uvdd9/FXXfdhfr6esydO7fV5WlacopSCJFyujJ37lzU1dWZr61bt3Y6LiIiIqJckXNtAC+55BKcddZZrc4zdOhQ3HzzzVi3bh2CwaDjs7Fjx+LnP/85HnnkEVRWVuKtt95yfF5bW4toNJp0Z9AuGAwmLZeIiIg6iMPA5JycqwCWl5ejvLy8zfnuvPNO3Hzzzebv27dvx6RJk/DEE09g3LhxAIDq6mosWLAAO3bswIABAwAYHUOCwSDGjBmTmQCIiIjIQdONlxvLIXfkXAUwXYMHD3b8XlxcDAA4+OCDMXDgQADAxIkTMXLkSEybNg233XYb9u7diyuvvBIzZ85kD2AiIiLKWznXBtBNXq8XL7zwAkKhEI4++micccYZmDp1Km6//fZsbxoREVH+EC6+yBXd9g5goqFDh0KkGB9o8ODBeP7557OwRURERATwWcC5qEffASQiIiKiZD3mDiARdQE++YWoa/S0W118EkjOYQWQiIiIMoop4NzDFDARERFRnuEdQCIiIsosDgSdc3gHkIiIiCjP8A4gERERZRTbAOYeVgCpe+JVgKhrsOc3uYG9gHMOU8BEREREeYZ3AImIiCijmALOPawAUvek0lK8GhBlliaYBqbOYy/gnMMUMBEREVGe4R1AIiIiyiimgHMP7wASERER5RneAaTuif8GEnUNtv8jN+jCeLmxHHIFK4BERESUWewEknOYAiYiIiLKM7wDSERERBmlwaVOIJ1fBEmsABIREVFm8VFwOYcpYCIiIqI8wzuAncUeckRdgz2/ibotjgOYe3gHkIiIiCjP8A4gERERZRaHgck5rAB2lrofzVQwUWapc4w5IKJuRxMCmgsdONxYBhmYAiYiIiLKM7wDSERERJmly5cbyyFXsALYWdlM/TIVlh1M92dHNo93lnnX4/WtR2EKOPcwBUxERESUZ3gHkIiIiDKLvYBzDiuAnZXNXsD53Csym7Gz53d2ZLvM87W8sxV7Pl/feiI+Ci7nMAVMRERElGd4B5CIiIgyio+Cyz28A0hERESUZ3gHsLPydRiYbLeHytfY8/nf33ws81wo73yMPdvXt56IbQBzTre/A/jCCy9g3LhxKCgoQHl5OU4//XTH51u2bMFpp52GoqIilJeXY9asWYhEIlnaWiIiovyj6e69yB3d+g7gk08+iZkzZ2LhwoX4/ve/DyEEPvzwQ/PzeDyOU045Bf369cPq1auxZ88eTJ8+HUII3HXXXVncciIiIqLs6bYVwFgshssuuwy33XYbzj//fHP6oYceav68YsUKfPLJJ9i6dSuqqqoAAHfccQdmzJiBBQsWoLS0tPMbkq/DwGR7KJR8jV1o2U2NZTt2+zZ0pXweCiUfY8/29U2tOxeaALiFKeCc021TwP/85z+xbds2eDwejB49GgMGDMBJJ52Ejz/+2Jxn7dq1GDVqlFn5A4BJkyYhHA5j/fr12dhsIiKi/CNcfJErum0F8PPPPwcAzJ8/H9dccw2ef/559O7dG8cffzz27t0LAKipqUFFRYXje71790YgEEBNTU2Lyw6Hw6ivr3e8iIiIiHqKnKsAzp8/H5qmtfp69913oetGS9B58+bhJz/5CcaMGYOHHnoImqbhb3/7m7k8TUu+hS+ESDldWbRoEcrKyszXoEGDUs8oNOuVDW4NrNQRuRB3V8duL+986xmZzdjt5Z3N+LMh2zFn61jPhfLOdpn3pPQvAE0I117kjpyrAF5yySXYuHFjq69Ro0ZhwIABAICRI0ea3w0GgzjooIOwZcsWAEBlZWXSnb7a2lpEo9GkO4N2c+fORV1dnfnaunVrBiIlIiKiXFFbW4tp06aZN3+mTZuGffv2tfodIQTmz5+PqqoqFBQUYMKECY6maICRVbz00ktRXl6OoqIi/OhHP8JXX31lfv7FF1/g/PPPx7Bhw1BQUICDDz4Y119/fcZHLMm5TiDl5eUoLy9vc74xY8YgGAxi06ZNOOaYYwAA0WgUX3zxBYYMGQIAqK6uxoIFC7Bjxw6zwrhixQoEg0GMGTOmxWUHg0EEg0EXoiEiIqLu0AnknHPOwVdffYWXX34ZAPDrX/8a06ZNw3PPPdfidxYvXowlS5bg4Ycfxje/+U3cfPPNOPHEE7Fp0yaUlJQAAGbPno3nnnsOjz/+OPr27YsrrrgCp556KtavXw+v14tPP/0Uuq7j/vvvxyGHHIKPPvoIM2fOxIEDB3D77bdnLF5NiO57P3X27Nn43//9X/zxj3/EkCFDcNttt+G5557Dp59+it69eyMej+PII49ERUUFbrvtNuzduxczZszA1KlT2zUMTH19PcrKyjDk1pvhCYVanzkfBwnOx5iVfI092wPl5mvs+XiOK/lY5hmIWW9uxpdXX4O6ujp3RsJog/r7ecJ35sLnbePvZxpi8Wa89s9Frm//xo0bMXLkSKxbtw7jxo0DAKxbtw7V1dX49NNPHSOMKEIIVFVVYfbs2bj66qsBGHf7KioqcOutt+KCCy5AXV0d+vXrhz/96U8488wzAQDbt2/HoEGD8OKLL2LSpEkpt+e2227Dvffea/Z3yIScSwG3x2233YazzjoL06ZNw1FHHYUvv/wSf//739G7d28AgNfrxQsvvIBQKISjjz4aZ5xxBqZOnZrRGjURERF1L2vXrkVZWZlZ+QOA733veygrK8OaNWtSfmfz5s2oqanBxIkTzWnBYBDHH3+8+Z3169cjGo065qmqqsKoUaNaXC4A1NXVoU+fPp0Nq1U5lwJuD7/fj9tvv73VCt3gwYPx/PPPd+FWERERkZ1bHTjUMhJH5+hs062amhr0798/aXr//v1bHDVETU/sU1BRUYEvv/zSnCcQCJg3puzztLTc//znP7jrrrtwxx13tDuO9ujWdwCzJrFHZFf2lkvsEdlVPcayGbPS1THb5Vt522WrzLMZez6e40quneP5UOZAds/xriBgtQPs1MtY3KBBgxyjdSxatCjlatMdWQTo2Kghqb6Xzndammf79u2YPHkyfvazn+FXv/pVq8vorG59B5CIiIjyz9atWx1tAFu6+3fJJZfgrLPOanVZQ4cOxQcffICvv/466bNdu3a1OGpIZWUlAOMun+poCgA7d+40v1NZWYlIJILa2lrHXcCdO3di/PjxjuVt374dJ5xwAqqrq/HAAw+0us1uYAWQiIiIMsvlXsClpaVpdQJJd2SR6upq1NXV4e2338Z3v/tdAMBbb72Furq6pIqaMmzYMFRWVmLlypUYPXo0ACASiWDVqlW49dZbARgjlvj9fqxcuRJnnHEGAGDHjh346KOPsHjxYnNZ27ZtwwknnGCOaezxZD5ByxQwERER5bURI0Zg8uTJmDlzJtatW4d169Zh5syZOPXUUx09gA877DAsX74cgJH6nT17NhYuXIjly5fjo48+wowZM1BYWIhzzjkHAFBWVobzzz8fV1xxBf7v//4PGzZswC9+8Qscfvjh+OEPfwjAuPM3YcIEDBo0CLfffjt27dqFmpqaVp9Y5gbeAeyIdNpmZKrdSFvLzVS7kWzGnO7yM9lmpq1lZ6u8gcy3FcrH2PPxHG/P8vOxzLN9fQO6d7tAHYAbu1B3YRkt+Mtf/oJZs2aZPXZ/9KMfYdmyZY55Nm3ahLq6OvP3q666Ck1NTbjoootQW1uLcePGYcWKFeYYgADwu9/9Dj6fD2eccQaamprwgx/8AA8//DC8Xi8AY3zizz77DJ999hkGDhzoWF8mR+rr1uMAdpV2jQOo9KBxo9KW7bHC8jX2bP9RyNfY8/EcV/KxzLN9fQNciT1b4wD+YNRV8Hk7/4CFWDyM//tocZdtf0/GFDARERFRnmEKuD2yMSRCNvSQ/3Q7LB/v7GS7zPM19mzfzcvX2PPxHFeyFXs3eBRcvmEFkIiIiDKLFcCcwxQwERERUZ7hHcD26Ojo7B295e7GrfqObK8baYrObnu2Yu/M9+w6sv3dPWYlH2Pvbue4G99VumPs3fkcV3LhXG8P3gHMOawAEhERUWZ1g2Fg8g1TwERERER5hncA2yOTvYDzscdvPvYCVPI19mz3gszH2PMxZrt87PHbWsxZ2h+aENBcSN+6sQwy8A4gERERUZ7hHUAiIiLKLHYCyTmsALZHR3sBK9m+LZ9q27sqTZEqvq5KRWQr7myXN5B7sWfrOE/nMzdl63jnOd6+z9ySa7Fne3+konfy76d9OeQKpoCJiIiI8gzvABIREVFmMQWcc1gBbA83ewHnag+xrpCvsedjL1AlH8ucPWCzIx/PcSWd2LO2f1yqACIHzqsegilgIiIiojzDO4BERESUWUwB5xzeASQiIiLKM7wD2B6dHQYGsNpfZKMdhtr2bLRTscebr7FnK+5caI+Vj2WezZhb+r0rZDv2fDvHlXTP9WwOA+NG+z0OA+MaVgCJiIgos4RuvNxYDrmCKWAiIiKiPMM7gO3hRmoj28MEANlPU2RDPg+Fko/lDWS/zPM19mxf4/I19nTjztb+YSeQnMMKIBEREWUW2wDmHKaAiYiIiPIM7wC2h5u9gLMhm73klHx+MkE+9wLOhmzHnc+xZwN7Abs7n9uYAs45vANIRERElGd4B5CIiIgyS8ClO4CdXwQZWAFsj+7eCzjbvSIB9oTNhnyNPV/jBvIz9ny+vgHsBUztxhQwERERUZ7p1hXAf/3rX5gyZQrKy8tRWlqKo48+Gq+99ppjni1btuC0005DUVERysvLMWvWLEQikSxtMRERUR7Sdfde5IpunQI+5ZRT8M1vfhN///vfUVBQgKVLl+LUU0/Ff/7zH1RWViIej+OUU05Bv379sHr1auzZswfTp0+HEAJ33XVX+1fYU3oBZ1O2ewhmQ7bLPJ9jz5ZspyPzMfZciLk7HO/sBUxSt70DuHv3bnz22Wf4zW9+g29/+9sYPnw4brnlFjQ2NuLjjz8GAKxYsQKffPIJ/vznP2P06NH44Q9/iDvuuAO///3vUV9fn+UIiIiIiLKj21YA+/btixEjRuDRRx/FgQMHEIvFcP/996OiogJjxowBAKxduxajRo1CVVWV+b1JkyYhHA5j/fr12dp0IiKi/KLuALrxIld02xSwpmlYuXIlpkyZgpKSEng8HlRUVODll19Gr169AAA1NTWoqKhwfK93794IBAKoqalpcdnhcBjhcNj8nXcLiYiIOoGPgss5OVcBnD9/Pm644YZW53nnnXcwZswYXHTRRejfvz/efPNNFBQU4A9/+ANOPfVUvPPOOxgwYAAAo6KYSAiRcrqyaNGi1NvAYWA6Lx+HiMjnMs/X2LPd3jYfY8/HmO1yfRgYyjk5VwG85JJLcNZZZ7U6z9ChQ/H3v/8dzz//PGpra1FaWgoAuOeee7By5Uo88sgj+M1vfoPKykq89dZbju/W1tYiGo0m3Rm0mzt3LubMmWP+Xl9fj0GDBnUiKiIiovwlhA4hOt+D141lkCHnKoDl5eUoLy9vc77GxkYAgMfjbMbo8Xigy27i1dXVWLBgAXbs2GHeEVyxYgWCwaDZTjCVYDCIYDDY0RCIiIiIclrOVQDTVV1djd69e2P69Om47rrrUFBQgN///vfYvHkzTjnlFADAxIkTMXLkSEybNg233XYb9u7diyuvvBIzZ8407xq2C4eB6bhsD5OQC7FnQ7bTcfkYO5D94zyfY8+GXDjOu8MwMG6032MnENd0217A5eXlePnll7F//358//vfx9ixY7F69Wo888wzOOKIIwAAXq8XL7zwAkKhEI4++micccYZmDp1Km6//fYsbz0REVEeYS/gnNNt7wACwNixY/HKK6+0Os/gwYPx/PPPd9EWEREREeW+bl0B7HId7QVsv+Wejdvv2UzHZTMlZY+3q2PP9zLPdnmn+r0rZPt4z2bMLf2eSdk8x4Hsl7f9vb3f62q6DmgudOBgJxDXsAJIREREmSVcGgeQKWDXdNs2gERERETUMbwD2BVS3XLvypRB4rq6MgWQuK5sxd3VaY9cKvNsxp7N41zJ1vHe1WnBbJV5rpzjSr7GnuMDPAtdh3AhBcxxAN3DO4BEREREeYZ3AImIiCiz2AYw57AC2B5uDgSdCwOHZkM+DxKbz7FnQ7YHRc6F2LMh28d5PsYO5P5A0LpLowOwAugapoCJiIiI8gzvABIREVFmCQHAjXEAeQfQLawAEhERUUYJXUC4kAIWrAC6hilgIiIiojzDO4BERESUWUKHOylgjgPoFt4BJCIiIsozvANIREREGcU2gLmHFUAiIiLKLKaAcw4rgGlQ/3Hozc0uLCwHnteYjwOl5mPMSr7Gns24gfyNPR/PcaUblLn6O9bVd9JiiLryIJAYop1fCAEANMH7qW366quvMGjQoGxvBhERkSu2bt2KgQMHZnw9zc3NGDZsGGpqalxbZmVlJTZv3oxQKOTaMvMRK4Bp0HUdmzZtwsiRI7F161aUlpZme5Myrr6+HoMGDcqbeAHGnA8x51u8AGPOh5jbE68QAg0NDaiqqoLH0zX9QJubmxGJRFxbXiAQYOXPBUwBp8Hj8eAb3/gGAKC0tDQvLihKvsULMOZ8kG/xAow5H6Qbb1lZWRdsjSUUCrHCloM4DAwRERFRnmEFkIiIiCjPsAKYpmAwiOuvvx7BYDDbm9Il8i1egDHng3yLF2DM+SDf4iV3sBMIERERUZ7hHUAiIiKiPMMKIBEREVGeYQWQiIiIKM+wApiGe+65B8OGDUMoFMKYMWPw5ptvZnuTXDF//nxomuZ4VVZWmp8LITB//nxUVVWhoKAAEyZMwMcff5zFLW6/N954A6eddhqqqqqgaRqefvppx+fpxBgOh3HppZeivLwcRUVF+NGPfoSvvvqqC6Non7ZinjFjRlK5f+9733PM051iXrRoEY466iiUlJSgf//+mDp1KjZt2uSYp6eVczox96Ryvvfee/Htb3/bHOeuuroaL730kvl5TytfoO2Ye1L5UnawAtiGJ554ArNnz8a8efOwYcMGHHvssTjppJOwZcuWbG+aK771rW9hx44d5uvDDz80P1u8eDGWLFmCZcuW4Z133kFlZSVOPPFENDQ0ZHGL2+fAgQM44ogjsGzZspSfpxPj7NmzsXz5cjz++ONYvXo19u/fj1NPPRXxeLyrwmiXtmIGgMmTJzvK/cUXX3R83p1iXrVqFS6++GKsW7cOK1euRCwWw8SJE3HgwAFznp5WzunEDPScch44cCBuueUWvPvuu3j33Xfx/e9/H1OmTDEreT2tfIG2YwZ6TvlSlghq1Xe/+11x4YUXOqYddthh4je/+U2Wtsg9119/vTjiiCNSfqbruqisrBS33HKLOa25uVmUlZWJ++67r4u20F0AxPLly83f04lx3759wu/3i8cff9ycZ9u2bcLj8YiXX365y7a9oxJjFkKI6dOniylTprT4ne4e886dOwUAsWrVKiFEfpRzYsxC9Pxy7t27t/jDH/6QF+WrqJiF6PnlS5nHO4CtiEQiWL9+PSZOnOiYPnHiRKxZsyZLW+Wuf//736iqqsKwYcNw1lln4fPPPwcAbN68GTU1NY7Yg8Egjj/++B4Tezoxrl+/HtFo1DFPVVUVRo0a1a33w+uvv47+/fvjm9/8JmbOnImdO3ean3X3mOvq6gAAffr0AZAf5ZwYs9ITyzkej+Pxxx/HgQMHUF1dnRflmxiz0hPLl7oOnwXcit27dyMej6OiosIxvaKiAjU1NVnaKveMGzcOjz76KL75zW/i66+/xs0334zx48fj448/NuNLFfuXX36Zjc11XTox1tTUIBAIoHfv3knzdNdj4KSTTsLPfvYzDBkyBJs3b8a1116L73//+1i/fj2CwWC3jlkIgTlz5uCYY47BqFGjAPT8ck4VM9DzyvnDDz9EdXU1mpubUVxcjOXLl2PkyJFmZaYnlm9LMQM9r3yp67ECmAZN0xy/CyGSpnVHJ510kvnz4Ycfjurqahx88MF45JFHzMbEPTV2u47E2J33w5lnnmn+PGrUKIwdOxZDhgzBCy+8gNNPP73F73WHmC+55BJ88MEHWL16ddJnPbWcW4q5p5XzoYceivfeew/79u3Dk08+ienTp2PVqlXm5z2xfFuKeeTIkT2ufKnrMQXcivLycni93qT/lnbu3Jn032ZPUFRUhMMPPxz//ve/zd7APTn2dGKsrKxEJBJBbW1ti/N0dwMGDMCQIUPw73//G0D3jfnSSy/Fs88+i9deew0DBw40p/fkcm4p5lS6ezkHAgEccsghGDt2LBYtWoQjjjgC/+///b8eXb4txZxKdy9f6nqsALYiEAhgzJgxWLlypWP6ypUrMX78+CxtVeaEw2Fs3LgRAwYMwLBhw1BZWemIPRKJYNWqVT0m9nRiHDNmDPx+v2OeHTt24KOPPuox+2HPnj3YunUrBgwYAKD7xSyEwCWXXIKnnnoKf//73zFs2DDH5z2xnNuKOZXuXs6JhBAIh8M9snxbomJOpaeVL3WBLu920s08/vjjwu/3iwcffFB88sknYvbs2aKoqEh88cUX2d60TrviiivE66+/Lj7//HOxbt06ceqpp4qSkhIztltuuUWUlZWJp556Snz44Yfi7LPPFgMGDBD19fVZ3vL0NTQ0iA0bNogNGzYIAGLJkiViw4YN4ssvvxRCpBfjhRdeKAYOHCheffVV8c9//lN8//vfF0cccYSIxWLZCqtVrcXc0NAgrrjiCrFmzRqxefNm8dprr4nq6mrxjW98o9vG/F//9V+irKxMvP7662LHjh3mq7Gx0Zynp5VzWzH3tHKeO3eueOONN8TmzZvFBx98IH77298Kj8cjVqxYIYToeeUrROsx97TypexgBTANd999txgyZIgIBALiO9/5jmOohe7szDPPFAMGDBB+v19UVVWJ008/XXz88cfm57qui+uvv15UVlaKYDAojjvuOPHhhx9mcYvb77XXXhMAkl7Tp08XQqQXY1NTk7jkkktEnz59REFBgTj11FPFli1bshBNelqLubGxUUycOFH069dP+P1+MXjwYDF9+vSkeLpTzKliBSAeeughc56eVs5txdzTyvm8884zr8H9+vUTP/jBD8zKnxA9r3yFaD3mnla+lB2aEEJ03f1GIiIiIso2tgEkIiIiyjOsABIRERHlGVYAiYiIiPIMK4BEREREeYYVQCIiIqI8wwogERERUZ5hBZCIiIgoz7ACSERERJRnWAEk6sYmTJiA2bNn96j1zpgxA1OnTu30cjZt2oTKyko0NDS0OM/DDz+MXr16dXpddjt37kS/fv2wbds2V5dLROQmVgCJqN2eeuop3HTTTebvQ4cOxdKlS7O3QSnMmzcPF198MUpKSrp0vf3798e0adNw/fXXd+l6iYjagxVAImq3Pn36dHnFqj2++uorPPvss/jlL3+ZlfX/8pe/xF/+8hfU1tZmZf1ERG1hBZCoB6mtrcW5556L3r17o7CwECeddBL+/e9/m5+rlOcrr7yCESNGoLi4GJMnT8aOHTvMeWKxGGbNmoVevXqhb9++uPrqqzF9+nRHWtaeAp4wYQK+/PJLXH755dA0DZqmAQDmz5+PI4880rF9S5cuxdChQ83f4/E45syZY67rqquuQuLjyYUQWLx4MQ466CAUFBTgiCOOwP/+7/+2uh/+53/+B0cccQQGDhzomP7www9j8ODBKCwsxI9//GPs2bMn6bvPPfccxowZg1AohIMOOgg33HADYrGY+fmnn36KY445BqFQCCNHjsSrr74KTdPw9NNPm/McfvjhqKysxPLly1vdTiKibGEFkKgHmTFjBt599108++yzWLt2LYQQOPnkkxGNRs15Ghsbcfvtt+NPf/oT3njjDWzZsgVXXnml+fmtt96Kv/zlL3jooYfwj3/8A/X19Y7KTaKnnnoKAwcOxI033ogdO3Y4KpNtueOOO/DHP/4RDz74IFavXo29e/cmVZquueYaPPTQQ7j33nvx8ccf4/LLL8cvfvELrFq1qsXlvvHGGxg7dqxj2ltvvYXzzjsPF110Ed577z2ccMIJuPnmmx3zvPLKK/jFL36BWbNm4ZNPPsH999+Phx9+GAsWLAAA6LqOqVOnorCwEG+99RYeeOABzJs3L+U2fPe738Wbb76Z9r4gIupSgoi6reOPP15cdtllQggh/vWvfwkA4h//+If5+e7du0VBQYH4n//5HyGEEA899JAAID777DNznrvvvltUVFSYv1dUVIjbbrvN/D0Wi4nBgweLKVOmpFyvEEIMGTJE/O53v3Ns2/XXXy+OOOIIx7Tf/e53YsiQIebvAwYMELfccov5ezQaFQMHDjTXtX//fhEKhcSaNWscyzn//PPF2Wef3eJ+OeKII8SNN97omHb22WeLyZMnO6adeeaZoqyszPz92GOPFQsXLnTM86c//UkMGDBACCHESy+9JHw+n9ixY4f5+cqVKwUAsXz5csf3Lr/8cjFhwoQWt5GIKJt8Wa5/EpFLNm7cCJ/Ph3HjxpnT+vbti0MPPRQbN240pxUWFuLggw82fx8wYAB27twJAKirq8PXX3+N7373u+bnXq8XY8aMga7rrm5vXV0dduzYgerqanOaz+fD2LFjzTTwJ598gubmZpx44omO70YiEYwePbrFZTc1NSEUCjmmbdy4ET/+8Y8d06qrq/Hyyy+bv69fvx7vvPOOeccPMNLUzc3NaGxsxKZNmzBo0CBUVlaan9v3lV1BQQEaGxtb3EYiomxiBZCohxAJbefs01W7PADw+/2OzzVNS/quff7Wlt0aj8eT9D17KjodqtL5wgsv4Bvf+Ibjs2Aw2OL3ysvLkzpgpBODruu44YYbcPrppyd9FgqFkvZla/bu3Yt+/fqlNS8RUVdjG0CiHmLkyJGIxWJ46623zGl79uzBv/71L4wYMSKtZZSVlaGiogJvv/22OS0ej2PDhg2tfi8QCCAejzum9evXDzU1NY6K13vvvedY14ABA7Bu3TpzWiwWw/r16x0xBYNBbNmyBYcccojjNWjQoBa3Z/To0fjkk08c00aOHOlYF4Ck37/zne9g06ZNSes65JBD4PF4cNhhh2HLli34+uuvze+88847Kbfho48+avUuJRFRNvEOIFEPMXz4cEyZMgUzZ87E/fffj5KSEvzmN7/BN77xDUyZMiXt5Vx66aVYtGgRDjnkEBx22GG46667UFtb2+qdr6FDh+KNN97AWWedhWAwiPLyckyYMAG7du3C4sWL8dOf/hQvv/wyXnrpJZSWlprfu+yyy3DLLbdg+PDhGDFiBJYsWYJ9+/aZn5eUlODKK6/E5ZdfDl3Xccwxx6C+vh5r1qxBcXExpk+fnnJ7Jk2ahF/96leIx+Pwer0AgFmzZmH8+PFYvHgxpk6dihUrVjjSvwBw3XXX4dRTT8WgQYPws5/9DB6PBx988AE+/PBD3HzzzTjxxBNx8MEHY/r06Vi8eDEaGhrMTiD2/dPY2Ij169dj4cKFae93IqKuxDuARD3IQw89hDFjxuDUU09FdXU1hBB48cUXk9K+rbn66qtx9tln49xzz0V1dTWKi4sxadKkpDZ1djfeeCO++OILHHzwwWbac8SIEbjnnntw991344gjjsDbb7/t6G0MAFdccQXOPfdczJgxA9XV1SgpKUlqp3fTTTfhuuuuw6JFizBixAhMmjQJzz33HIYNG9bi9px88snw+/149dVXzWnf+9738Ic//AF33XUXjjzySKxYsQLXXHON43uTJk3C888/j5UrV+Koo47C9773PSxZsgRDhgwBYLSHfPrpp7F//34cddRR+NWvfmUuw75/nnnmGQwePBjHHntsa7uaiChrNNGRxj1ElDd0XceIESNwxhlnOJ7+kevuuecePPPMM3jllVcyup5//OMfOOaYY/DZZ5+ZnWu++93vYvbs2TjnnHMyum4ioo5iCpiIHL788kusWLECxx9/PMLhMJYtW4bNmzd3u8rMr3/9a9TW1qKhocHVp5YsX74cxcXFGD58OD777DNcdtllOProo83K386dO/HTn/4UZ599tmvrJCJyG+8AEpHD1q1bcdZZZ+Gjjz6CEAKjRo3CLbfcguOOOy7bm5YTHn30Udx0003YunUrysvL8cMf/hB33HEH+vbtm+1NIyJKGyuARERERHmGnUCIiIiI8gwrgERERER5hhVAIiIiojzDCiARERFRnmEFkIiIiCjPsAJIRERElGdYASQiIiLKM6wAEhEREeUZVgCJiIiI8gwrgERERER5hhVAIiIiojzDCiARERFRnmEFkIiIiCjPsAJIRERElGf+P13zOSIUHZw4AAAAAElFTkSuQmCC", "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 a nonlinear IVP, starting with the background jet plus a small amount of the most unstable mode, to see how it saturates.\n", "\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": [ "2025-10-07 22:23:42,771 subsystems 0/1 INFO :: Building subproblem matrices 1/127 (~1%) Elapsed: 0s, Remaining: 4s, Rate: 2.9e+01/s\n", "2025-10-07 22:23:42,996 subsystems 0/1 INFO :: Building subproblem matrices 13/127 (~10%) Elapsed: 0s, Remaining: 2s, Rate: 5.0e+01/s\n", "2025-10-07 22:23:43,237 subsystems 0/1 INFO :: Building subproblem matrices 26/127 (~20%) Elapsed: 1s, Remaining: 2s, Rate: 5.2e+01/s\n", "2025-10-07 22:23:43,497 subsystems 0/1 INFO :: Building subproblem matrices 39/127 (~31%) Elapsed: 1s, Remaining: 2s, Rate: 5.1e+01/s\n", "2025-10-07 22:23:43,753 subsystems 0/1 INFO :: Building subproblem matrices 52/127 (~41%) Elapsed: 1s, Remaining: 1s, Rate: 5.1e+01/s\n", "2025-10-07 22:23:44,016 subsystems 0/1 INFO :: Building subproblem matrices 65/127 (~51%) Elapsed: 1s, Remaining: 1s, Rate: 5.1e+01/s\n", "2025-10-07 22:23:44,279 subsystems 0/1 INFO :: Building subproblem matrices 78/127 (~61%) Elapsed: 2s, Remaining: 1s, Rate: 5.1e+01/s\n", "2025-10-07 22:23:44,541 subsystems 0/1 INFO :: Building subproblem matrices 91/127 (~72%) Elapsed: 2s, Remaining: 1s, Rate: 5.0e+01/s\n", "2025-10-07 22:23:44,852 subsystems 0/1 INFO :: Building subproblem matrices 104/127 (~82%) Elapsed: 2s, Remaining: 0s, Rate: 4.9e+01/s\n", "2025-10-07 22:23:45,114 subsystems 0/1 INFO :: Building subproblem matrices 117/127 (~92%) Elapsed: 2s, Remaining: 0s, Rate: 4.9e+01/s\n", "2025-10-07 22:23:45,315 subsystems 0/1 INFO :: Building subproblem matrices 127/127 (~100%) Elapsed: 3s, Remaining: 0s, Rate: 4.9e+01/s\n", "2025-10-07 22:23:45,347 __main__ 0/1 INFO :: Starting main loop\n", "2025-10-07 22:24:08,242 __main__ 0/1 INFO :: Iteration=1, Time=1.666667e-01, dt=1.666667e-01\n", "2025-10-07 22:24:19,463 __main__ 0/1 INFO :: Iteration=101, Time=1.683333e+01, dt=1.666667e-01\n", "2025-10-07 22:24:31,055 __main__ 0/1 INFO :: Iteration=201, Time=3.350000e+01, dt=1.666667e-01\n", "2025-10-07 22:24:42,221 __main__ 0/1 INFO :: Iteration=301, Time=5.016667e+01, dt=1.666667e-01\n", "2025-10-07 22:24:53,116 __main__ 0/1 INFO :: Iteration=401, Time=6.683333e+01, dt=1.666667e-01\n", "2025-10-07 22:25:04,944 __main__ 0/1 INFO :: Iteration=501, Time=8.350000e+01, dt=1.666667e-01\n", "2025-10-07 22:25:16,124 __main__ 0/1 INFO :: Iteration=601, Time=1.001667e+02, dt=1.666667e-01\n", "2025-10-07 22:25:27,921 __main__ 0/1 INFO :: Iteration=701, Time=1.168333e+02, dt=1.666667e-01\n", "2025-10-07 22:25:39,409 __main__ 0/1 INFO :: Iteration=801, Time=1.335000e+02, dt=1.666667e-01\n", "2025-10-07 22:25:50,760 __main__ 0/1 INFO :: Iteration=901, Time=1.501667e+02, dt=1.666667e-01\n", "2025-10-07 22:26:02,349 __main__ 0/1 INFO :: Iteration=1001, Time=1.668333e+02, dt=1.666667e-01\n", "2025-10-07 22:26:13,837 __main__ 0/1 INFO :: Iteration=1101, Time=1.835000e+02, dt=1.666667e-01\n", "2025-10-07 22:26:25,327 __main__ 0/1 INFO :: Iteration=1201, Time=2.001667e+02, dt=1.666667e-01\n", "2025-10-07 22:26:36,481 __main__ 0/1 INFO :: Iteration=1301, Time=2.168333e+02, dt=1.666667e-01\n", "2025-10-07 22:26:47,621 __main__ 0/1 INFO :: Iteration=1401, Time=2.335000e+02, dt=1.666667e-01\n", "2025-10-07 22:26:52,078 solvers 0/1 INFO :: Simulation stop time reached.\n", "2025-10-07 22:26:52,079 solvers 0/1 INFO :: Final iteration: 1441\n", "2025-10-07 22:26:52,080 solvers 0/1 INFO :: Final sim time: 240.16666666666174\n", "2025-10-07 22:26:52,091 solvers 0/1 INFO :: Setup time (init - iter 0): 2.681 sec\n", "2025-10-07 22:26:52,091 solvers 0/1 INFO :: Warmup time (iter 0-10): 23.89 sec\n", "2025-10-07 22:26:52,092 solvers 0/1 INFO :: Run time (iter 10-end): 162.8 sec\n", "2025-10-07 22:26:52,092 solvers 0/1 INFO :: CPU time (iter 10-end): 0.04523 cpu-hr\n", "2025-10-07 22:26:52,093 solvers 0/1 INFO :: Speed: 8.571e+05 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": [ "2025-10-07 22:26:53,952 matplotlib.animation 0/1 INFO :: Animation.save using \n", "2025-10-07 22:26:53,954 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/tmpubfrdw_f/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 }