PlotTransects.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cdf = \"model_results/test_do_transect.nc\"\n",
    "\n",
    "transect_shp = \"gis/transect_sjf_thru_puget_sound.shp\"\n",
    "variable = \"controlDOXG\"\n",
    "start_date = \"2014-01-01\"\n",
    "\n",
    "# TODO:\n",
    "# - transect name\n",
    "# - location labels along the transect\n",
    "\n",
    "from pathlib import Path\n",
    "from os import path\n",
    "from netCDF4 import Dataset\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.animation as animation\n",
    "#%matplotlib widget\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'netCDF4._netCDF4.Dataset'>\n",
       "root group (NETCDF4 data model, file format HDF5):\n",
       "    dimensions(sizes): time(216), node(236), siglay(10)\n",
       "    variables(dimensions): int32 node(node), float32 time(time), float32 siglay(siglay), float32 zeta(time, node), float32 controlDOXG(time, siglay, node)\n",
       "    groups: "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = Dataset(data_cdf)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04166667"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timestep = data['time'][2] - data['time'][1]\n",
    "timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_base = Path(data_cdf).stem\n",
    "transect_base = Path(transect_shp).stem\n",
    "out_file = lambda name: path.join(\"notebook_outs\", out_base + \"_\" + transect_base + \"_\" + name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>depth</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>node_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1159</th>\n",
       "      <td>257.800720</td>\n",
       "      <td>POINT (374011.682 5370035.740)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1252</th>\n",
       "      <td>191.153000</td>\n",
       "      <td>POINT (377292.015 5370100.970)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1339</th>\n",
       "      <td>191.626007</td>\n",
       "      <td>POINT (380243.734 5368997.670)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1338</th>\n",
       "      <td>191.858002</td>\n",
       "      <td>POINT (383314.237 5367358.850)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1423</th>\n",
       "      <td>185.763000</td>\n",
       "      <td>POINT (386738.325 5366096.140)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              depth                        geometry\n",
       "node_id                                            \n",
       "1159     257.800720  POINT (374011.682 5370035.740)\n",
       "1252     191.153000  POINT (377292.015 5370100.970)\n",
       "1339     191.626007  POINT (380243.734 5368997.670)\n",
       "1338     191.858002  POINT (383314.237 5367358.850)\n",
       "1423     185.763000  POINT (386738.325 5366096.140)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transect = gpd.read_file(transect_shp)\n",
    "transect.set_index('node_id', inplace=True)\n",
    "transect.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>depth</th>\n",
       "      <th>geometry</th>\n",
       "      <th>distance</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>node_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1159</th>\n",
       "      <td>257.800720</td>\n",
       "      <td>POINT (374011.682 5370035.740)</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1252</th>\n",
       "      <td>191.153000</td>\n",
       "      <td>POINT (377292.015 5370100.970)</td>\n",
       "      <td>3.280981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1339</th>\n",
       "      <td>191.626007</td>\n",
       "      <td>POINT (380243.734 5368997.670)</td>\n",
       "      <td>6.432158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1338</th>\n",
       "      <td>191.858002</td>\n",
       "      <td>POINT (383314.237 5367358.850)</td>\n",
       "      <td>9.912635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1423</th>\n",
       "      <td>185.763000</td>\n",
       "      <td>POINT (386738.325 5366096.140)</td>\n",
       "      <td>13.562131</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              depth                        geometry   distance\n",
       "node_id                                                       \n",
       "1159     257.800720  POINT (374011.682 5370035.740)   0.000000\n",
       "1252     191.153000  POINT (377292.015 5370100.970)   3.280981\n",
       "1339     191.626007  POINT (380243.734 5368997.670)   6.432158\n",
       "1338     191.858002  POINT (383314.237 5367358.850)   9.912635\n",
       "1423     185.763000  POINT (386738.325 5366096.140)  13.562131"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the distance between adjacent nodes in the transect\n",
    "transect['distance'] = transect['geometry'][1:].distance(transect['geometry'][:-1], align=False)\n",
    "# Set the first one to zero since it's currently NaN\n",
    "transect.loc[transect.index == transect.index[0], 'distance'] = 0\n",
    "# Transform into a running total in km\n",
    "transect['distance'] = transect['distance'].cumsum() / 1000\n",
    "transect.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAF3CAYAAABAE4T5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB+QUlEQVR4nO3dd5xjdb3/8dcnZZJM77OzM9sLsLssCyxdBKSIFbzqtXd/XLx4r/Va71WvXhWvXq8oCu61YAEVC0VUBEQ6LCywbO91dman95JMku/vjyQzmZlkJsmk5/N8POaxk3NOzvkmk03e+VYxxqCUUkoppcCS6QIopZRSSmULDUZKKaWUUkEajJRSSimlgjQYKaWUUkoFaTBSSimllArSYKSUUkopFaTBSCmllFJZS0R+IiIdIrIjbFu1iDwoIvuD/1bNcn+riLwoIvfFcj0NRkoppZTKZrcBV0/b9hngb8aYVcDfgrej+QiwO9aLaTBSSimlVNYyxjwG9EzbfA3ws+DvPwOujXRfEWkGXgP8KNbraTBSSimlVK5pMMa0AQT/rY9y3HeATwH+WE9sm3fRslxtba1ZunRppouhlFJKpcXzzz/fZYypS9f1rrq0xHT3+BK+/wvb3TuBsbBNm4wxm+ZbLhF5LdBhjHleRC6N9X55H4yWLl3Kli1bMl0MpZRSKi1E5Gg6r9fd4+Op+xYnfH/nkv1jxpiNcd6tXUQajTFtItIIdEQ45iLg9SLyasAJlIvIL40x75ztxNqUppRSSqlccy/wnuDv7wHumX6AMeazxphmY8xS4K3Aw3OFItBgpJRSSqksJiK/Ap4GThGRFhH5AHAjcKWI7AeuDN5GRBaKyJ/nc728b0pTSimlVO4yxrwtyq7LIxzbCrw6wvZHgEdiuZ7WGCmllFJKBWkwUkoppZQK0qY0pZRSBcuY2P+N59hM/qvmJ+eCkYhcDdwEWIEfGWNuzHCRlMIY8Psnf0K30/lvrmzz+8HnC/yEfp9t21z75rp+OvbH8kGVjn3ZUIZY9mVD6FAqmpwKRiJiBb5PoAd6C/CciNxrjNmVyXKFv+HP9Xuqjk32NcJ/Qh9A8W5L1jG5cG6VPBYLWK2Bn9Dv0/8N/R76EZn8N9rvie63WmffH7odIhL539Dv5sTDk/swYfvM1H+n7GOWfVPvF9g/275p1wkv34x9ke7HLPtiL0ukckQ7T7RrRz8+9mvGtD9se6xl8i8aBcB2whG1TL5FY4gYrMecs14z/F8iXNe3ZHTK+W/4HGoecioYAecCB4wxhwBE5NcE1kuJGox27IDly2f/Fjjbt8O5AoV++5hb+AdY+AfcbNvme4zNlrpzR9sW6UN0Pv+aLV/EEnwjtog/sF38E2/OlontZsq+0H0Eg8XiRyD4r5myP3ROCTvP1PtOnmPkQMPkcQSuVbr6ZOSyTbuvMPUaVosfi/ixWgI/ltC/4sdi0f9QCsY3DgBg31Ke4ZLMz/jGgVkfQ7IeZ+g8IRqM5ifXglETcDzsdgtw3mx3KCmBiy6K7Rtj+O3p30QjfXDF+3uqjk3FNZIZHlLF/7t4J0rNHeMbB+DSwO/p/nAY3NEMQNm6lsltg2Mzjitb2DJjWzzXMIBrXeLnUOkVy4f49A/ouYTOFe1+cwWLbJfLZS9kuRaMJMK2GV8xReQ64DqAxYsX84tfpLpYuSPdYUJbmXJfKCjNdUyZhpy8FG/YiVUsoUGDhcqEXAtGLcCisNvNQOv0g4KLz20C2LhxY8rq5kMhY/obh/5nVkolU3g4TXcAtW8pT1o40jCkckGuBaPngFUisgw4QWDtk7fHcsdINSXz/s8eofJF/1OrQqG1RCqZkvHeGUvtZoi+dlU0KewBknzGGC/wYeCvwG7gTmPMzlnv1Ls7r/uiKJVM4R8ssTahxfNhpHJTqr/w6RdKlU1yrcYIY8yfgXktEKeUik26Qo/WPmW/+YaX+dw/0utw+usl3tdqpEEGSkGO1RgppfJD6MMo3hoqlTuS1TQW7XUR2jffWstMvu7sW8q1tiwL5VyNkVLpEv6NMpOdX/ORBqLCkGjH7XS/JrTGUoXTGiOlwkT6Bjr9TTqXPsjTWVbtb6QiibdWJFOvIX3tqhANRkqRnx/qoceTiW/fie7Pt79BoQmFoESah7Lh/2Cmr6+ygzalqYIzvdNlPr0ZztYfI1NSce1CbtrUTsNqumTONaU0GKkkSdYH1Vwfosn8MMinQATZ9Xjm2xl2ynIkMdRAFWJIyJeAlE2vW6VAg5FKUKzNIdHetHP9zbBQP4zTJZGh1/n698jHYJjr//9VftNglET5POwy0TeyZL8Bhn8IZPoDIR3Xz+fXVLLlSw1KuFj//2TzY9cQpHKNBqMkyccPsGx9Q5s+YizZHwYd+5sibq9fdSJqWZJdhnhfT9n6t0q3bAwGiUrkb5pNASlXX5OZ/sKVKO1nlDwFPSotGZNr5csEXdOHqefSm1oyyxotFM21L1MjarJhJI9Krnz4myaj/LP9f8vG82aDfPgsygZaY8TkiynWtJ0PL77Zvlnm4ptyur7lhd5UI9UepbMcoWup/JKsv2kmaz3ifQyxfBkJ//821//BeM6fj/Lh8ynTNBiFCVVFRgtK+fiCy6cP1/l+GMTzBjrbm3O0UVXpHlEX6fHE+mGSq3K1GSRVUxqk+7lIZihK5LhY79Oxvynq/4VcfP2o5NJgNE14+MnHIASpDUPxfptLtkRrwhL9Vhnt8Ua6ViwLYcZirr9fvN/A802uhaN8+XISz+NIdS1OLOePFo5y7fWjkk+DkUqa8Dej6W9M6f4gnj56bTbR3kSP7ls85faS1cdmPUcijzHeN+H5hKLpx+VzOEqnRIJNOptbU32tfAl2SoVoMFLzFuu3s3jN94M70eYmmBmKpm+LFJLmE45g7g/LZNd4FULtUTpMX2Q4luPzQT4EIv2CoCLRYJRm+VRNm67q8FS9ccUTiiIdEy0cQWJlTnczIGgoSob5zNae6nCRivea+ZY5FzpD59P7dL4QkY8A/w8Q4P+MMd+JcMylwHcAO9BljLkkkWtpMEqTVM+9k0g5EpXuN7ZUfKubTyiK9fyJljldfTXyORSl6/9YNoeiWE0vRy6UGaL/H8uF4KXiIyLrCISicwEPcL+I/MkYsz/smErgB8DVxphjIlKf6PU0GKVBpt9M0jXPTyols/YoWaEoWq1RumQyFGW6D1k2iOX/VbbXPMw1SCCeNetikcz3j0J8zRWw04BnjDEjACLyKPAG4L/Djnk78AdjzDEAY0xHohfTYJRi6QpFqb5OtnwLm09NzGyPIVk1Reky379HokEz2nULra/GfP+/ZfrLUqxlyIcvVXMptNdujtoBfFVEaoBR4NXAlmnHrAbsIvIIUAbcZIz5eSIX02CUQtHeVKJ9C5vt21myv7nFKt43s7k6KCdDPG9ksZQ/VaEoVW+4yfyA0Q7YqZfppvPp186GUJYLYv27Te8bGE/TZL6QUQv27cXzOUWtiIQHnU3GmE2hG8aY3SLyDeBBYAh4CfBOO4cNOBu4HHABT4vIM8aYffEWRoNRisS6+vx875Ns8/nQnR4wUhmSYgkdmQxFs5UhG0OIBqTYJbqGWS73K5yPdNYUpbtWKh+aU7NElzFm42wHGGN+DPwYQES+Bkx/UluC5xkGhkXkMeAMQINRpmT6zSdRyXgjiXUUV7hkBKX5lj1TzWfz6Z+TyeaIbG0KiSYVH0jz+X+eq+8R85Frr5npZnsNTf97agBKLRGpN8Z0iMhi4B+AC6Ydcg9ws4jYgCLgPOB/E7lWQQWjZL5wc/1NLpk1Q8k4RyY7MWdarE1umfjmrTVIk3L9/3y6pWM6j0hrqCXjXIkqxGa0NPp9sI/ROHCDMaZXRK4HMMbcGmxuux/YBviBHxljdiRyoYIJRtNXjZ/PCzaX3yAzHYjiOXchhaVs6gA62wzmc90vWx5DMuXy/3el8oUx5uII226ddvubwDfne62CCEbJfGPLpTfJZH5jy0SzU+iaqQhIyXo8yRyyP1uwSFdtUTJGu+VTOMql/++ZlOlpPLKVNq/lpoIIRsmQrDfIdHXGTVffoXRIZkDKlseUzxJtGsy2QKWhaHbZHkqyRTJaKVR65X0w8o3aI26PNkx+tuMSNdcbSLK+ZedqZ+RYzTcgZfvjg8ivhVz8AIrlNV2/6kRWPjYNRLPLxr9ZKiWr1kcDUu7I+2AUi2gv/HTO9JpoOMr3MBRJIs1XufQ486U5KpYO3KF92fRhG++isIUim/5GIcnqfK1UOA1GQdPDUSYWSkznf+xcCgqRxFN7lOrHmsqlQRJ9TaSyf1a8Mhn0Eh0lpOFoqkIPHcmqNdLaotygwShMIi/+XHzDyPVQFC6fHku4ZI0ezPR6biGz1R6F9uVLTZlKv1x8H1bZy5LpAmSbwR3NMX9TzLX/jEf3Lc7bIKGiy6a/+fT/M3PdVunVsb9J/yZRJKMGUWshc4PWGCUokTeLTE1qmE0fjCq1ov2t07GGXbyyeUHaQv8Ay5UwlM5yJntyYG1Wy15aYxSnSN+o5hKtpkYDS/7Qv2XsYvk/lIoPvFg/iAo9FClV6DQYxSEZtUTx7p8v/cBW4bKltigWiXwJmU2omTwbFmtW2SWWGspsfG0YY3jnTYf50/P9E9tuvKstgyXKDxqMUiSe/jwaXrLTvoONmS5CXPR1ND/Z+MGXTrnSfJZJ05eWyoSeQS+d/eM8s2+IN37zEL96ooev3Dkwsf/JPcM0V1WcnbEC5gHtYxSjeN40EvmASsXoIf2gnL9QOFq9Ij++hWXLKLVMS/b0HKpwZPK10tY7TvN122Zs39Paj9uzCEdR4BiIPLGxio3WGCXRfEd9JTPIaCian3yuLcql10YqO2HH0rSmVDZ5ZOdgxO2DY2427/MDoWCk5kODUQzmqi1K5jD40LkSrXXSIfnzFykU5VpQmou+RpSalCvNiE/uGcJqkYj7HnppBJ/P0DHgS3Op8k9BNKXFsjRBNoo2e7F+qKXOfAJQNs02HYv5NqvN9/Fm43IgSmWzJ/cMsbK+hpP9Q/SPjk1sF4HdLR46Bxz4/f4MljA/FFSNUfgbcKwjXrLhTTu8JkhDUerMFYpiDU259HeaT81ktNvxyIb/X0rlgoERH9uPjbKsrpLqkmIAVtRV8c+XnU+500lbr5uTfdqMlgwFFYxgZiCa7Y1Z37QLR7KaypasPjZRg5JL4SgZIyjjfby58sVEqVSIt2/b1iMjGAOLqyuoKQ0EozMXN7K8roaqYhddg27aNBglRd4HI6+7aM5j9M23sKWi/1AoIOVa7dF8g08qHqv+/1Spli2LHM/mxcMjACwor6QmWGNUXuwEoLLYRc/IKC1d3uQXsgDlfTCKVaLrA+XKh55KjnhDVHjtUa68VrJtlvZc6xuYizR8Zk6s4ejFwyNUFDsodzmpKy8BoMIZCEiVxU76RkY50qE1RsmgwSiMvjmo+YoUIKY3r2VbQNL5s5SGz8yJdamaFw+PsrCiHIDTmxbw/pdtZGFl4HZVsQuf37Dt6BiuorlbSdTsCmJUWjxSPZGjKlyhpjVI7gi2ZLwOw8sx/Xz6Olcq8453eVjfFKghslosrFnYMLGvssQFwLaTdkqdjoyUL59oMFIFLd3zEyW7Y/Z8zxMpmIUHOFUYCr22vGN/U8b7Gc1Wc+TzGQZGfZQ4I89oXVUcCEYtbSdZWluN16dzGc2HNqUpFadkhKnw5rVsk22hqNA/tFNNn9/0ea6lh83He+K+X/+ID2PAZY9cl1EeVktUUmTHGKOTGc2DBqMEZduHR6bl28zQc5nt8ebSayNVZc2l50CpdPn0X7Zzw73b8BtD2bqWmNfr6x0O1AAVR+k/VOwowiKBGbHLXM4klrgwaVNanPQNf6rwgBD++/RFV/cdbMybhVhDMv2Yktkcl4raq/DyZWvtWL6Jt0lIa4vSxxjD3q5BLBYbY16wRAhC0ZrVeoYCw/Cd9shNaRYRypwO+kfHKCt2QXdvkktfWDQYqYQkMkt0poOEii5V4ShZMtkHJJfoc5S92gbHGPb4AB89Yz6KozSLwcyAFKoxihaMAOxWK4B2vk4CbUqLg9YWzd++g41TfjJdllRJ9WtFX4sqmlhrgbS2aFI6AuX+7qGJ34/0jUU8JrymKPz3UI2Ryx59KL4luLhsiUOH68+X1hgVoGiBINbanEwHmmRZvaIt5eEollqYaHMfpVs2hy2tCYlNKOzMtXB2IYcivzET/XFS7dnjPVQXF7GyppT9XZPB6GjfKC9bVBrxPpFGp/UOhfoYRa8xsgYfU4lNP9bnK+ueQRH5EvD/gM7gps8ZY/4c3PdZ4AOAD/hXY8xf01WubP7QiFWii6SGB6Z8CUXpMj0cxbsmWTY3b8Uj25vq8lWkgFTIoahnbJRX3ncH7z3tDG5YuxFIXeg2xvDOO5+lqaKMv3/wIvZ1DU7sO9o3GvE+0foY9QZrjGYLRqGwV2RNT+hLNxH5GPBBwADbgfcZY8bC9lcAvwQWE8g23zLG/DSRa2VdMAr6X2PMt8I3iMga4K3AWmAh8JCIrDbG6IQNMZhPoEllGMpkv6N0hbz5hOp8COQhGo4yp5DDULj/ePYR+j1ubtu3jRvWbkxpTeSx/lE6htz0jo7TM+rlqaPdNFWWc6JvgJaByMEomt5hH0U2C7ZgP6JI1i9q5ETfABV5OCpNRJqAfwXWGGNGReROAnngtrDDbgB2GWNeJyJ1wF4Rud0Y44n3etkajCK5Bvi1McYNHBaRA8C5wNPJOPlszRm5/uGktTwz6XOSGRqOUkfDz1TGGJ5uP0FFkYO11XU81nqMB1sOAzDs8WBfdAyIHjTma0tLYL6icZ+fW587yo72Ad541lr+tH0vJwejB6NITWk9Q95Za4sAXnHqCs5dtoiy/O18bQNcIjIOFAOt0/YboExEBCgFeoCEVtXN1mD0YRF5N7AF+IQxphdoAp4JO6YluC1u8TZn5LJcCAC5PFotl8uu8kehhqKnTrbw8Scf4Ovnv4LLmpZObB/1jvOvT/yVR1sDIXxxWQWdo8PUlZZwxdqV/GrzSzxxrJ/Xra5OWdmeOdaD3WLBa/zc/MQerCKsX7SQJw4c5XjvIMYYJIa+TsYYDrS5o85hFCLBIfuZ4B+zM7Z7wTzO8FKtiGwJ27DJGLMpdMMYc0JEvgUcA0aBB4wxD0w7yc3AvQQCUxnwlkQnuszIqDQReUhEdkT4uQa4BVgBbADagP8J3S3CqUyU818nIltEZEvH0PjEwp3ZuIBnKuVCKMqEZD8v+jwrlRkvdbfT4x7jnx79M13ukYntP93zEo+2HuOaM9Zw5ZpVVFcWs2HxQt574dmsrKsF4G8HO6Oddt4O9QzxixePsq6pgRV1NYyMezltYT0ljiI2LFrI/o6eiYkew0WqLbr98R6e2DPExuVLUlbeLNBljNkY9rMpfKeIVBFoNVpGoCtNiYi8c9o5XglsDe7fANwsIuWJFCYjNUbGmCtiOU5E/g+4L3izBVgUtruZmVVpofNvAjYBrC5pjhie8lkuflCnq+YlF5+bTNs1dIw/tD9B9/gARRYb/7ni3TitOiQ4m9SvOlGQtUa9Y5NNUs/1n+TqumUY4M6Du1ldX8vFpyyLeL9ltdXcveMIN165Eqct0Jx2cnCMzcd7uGbNwnmX6ztPHACE129Yw662Dg50dHPu0sDf54o1K3F7vdyx9RD9Y+N89tJVLKkspswxs6nMGMO37mmnuaaKl60onC/1EVwBHDbGdAKIyB+ACwl0tg55H3CjMcYAB0TkMHAq8Gy8F8u6eYxEJPyT6w3AjuDv9wJvFRGHiCwDVpHAA853ufzBn8tlz1eP9Wznk3s38fzwXgaK+nhh4AB/Hn0i08VSCoCusGD0kYfu5+zf/ZjPPP03WocHOX959CBx+Wkr6BsZ4/ubj09s+/wDO3jPb59jc0v/vMo07vPzxz2trG1aQJnLycalzbzvoo2ctjDQ1CQivGb9qVy9bjV/2dvKy259hJXfup+fHds641yP7x5i+7FRXrZ6WdqmGMhSx4DzRaQ42IfocmB3hGMuBxCRBuAU4FAiF8u6YAT8t4hsF5FtwGXAxwCMMTuBO4FdwP3ADToiLSBbJkxMhlx9DLla7tn8seMZ/uvQ7TRVlfNvV1/Cv1x+IXVlJfy18/kpxw14h3lx4EDM502kObsQa0PU3LrGRlhYOdlaMjTu4e4j+7h45VLWLYre5+WUBXU0VpRx2/OH2dnezy3PHOTunYEGiJuePjKvMj16uJPe0XHOXBR4T7BaLKxtapjSn0hEuGLNKv7t6pfzrvPPpKmygk//8gTtfYGPNJ/P0NLt4b/vPklZsYv1jfXzKlOuM8ZsBn4HvEBgqL4F2CQi14vI9cHDvgJcKCLbgb8BnzbGdCVyvazrfG2Medcs+74KfDWNxcl6+fiBrB2aM+/oaDs/OP5HTl1Qz7svPAt7sLnh3GWL+NO2PexddIRTbEsB2HT8z/y1+3luOeMGVtoXzXJWpZKrc3SEmrJi2voHMAauOWMN569YPPF6jUZEOH/FYu56YSdX//QJBt1erCKsrK/h4f0t9I2todI5+yiwaO7a2YrLbmNlQ+2cx9aVlVJXVkptWQn/++AT/McdvXztnZWc95ndHOkMjDJ/zYY1FM3xeAqBMeaLwBenbb41bH8rcFUyrpWNNUYqBvlSQxRNPj+2bOc3fr595A84bFb+8ZwzpnzIbFzajEWEe/sCzWkD3hH+3vMSAL/qeTAj5VWFq2tshHKng1Af5oby0jlDUchZi5uwW624fYbrXn4O/3rFRVx+2krGxr384qWTCZXH4/Nz35421jY1zDrn0HRNVRWc3ryAnz92nP++u50jnR4uPmU5py1s4MJZmgRVamRdjVG+mesDPpGakUIJDVpzlBl/6nyWXcNHedvGDZS5pnayLnM6WNvUwBOte3lvRS+Pde/AY7wsr63mqRMH6KrtpdZalfQy6ZIgcyu0Dtgen49+j3vKEPWaspKY7+8qsvPuC8+i2G5nSW3gNes3hgqXk99tb+Ffzou/9vPvBzvoHxvnrEXxd+C+ZsMa9p3s5H/+eJJSh53Xrz81puH8Kvm0xiiFYgkw8YacQglFuSgf/jb948P8qOUvrK6v46xlkd/cX7l2NV6fn88f+TH3dW5maXU1b9y4Dq/fz28HHk5ziVWh6nEHOl6XOx1cvCow+qyqxBXXOU5rrJ8IRRBYVuPMxQvZ3trJsf7IC73O5q5drRQX2VheP3cz2nSVxS7efeHZAFy4apmGogzSYJQi8XxIxtIslu9NZyo77Bk+zojfzRWnrYr6xrygooy3nbeBw/1dtLi7uHDFYhrKy1hWW8Wz3fswkacXUyqpOkcD8xaVOh28bsNp3PjGq5MycmvD4oX4jeH2bfE1p7m9Pv60p411TQuwWRP7aD1lQR1fvuZKrjh1ZUL3V8lRcMEom8NFtLJlc5lVfjk21gHAgsrZmyROb17Aa9afyqKqCk5fHBj9s6qhlhN9/XRL96z31SVBVDJ0jQWCUYmzCItIXH16ZtNUWU5lsYv7987djL+va5DR8cBIsvv3tTPo9k6MRktUsaMIi0VrizKpoIJROgLGfGt2wu+vtUTJ/5tpn6XZHR1tp8zhoNgx9wSOl526go9c+TLswQ+kFfU1GGDL2PTpRZRKvl53oKmrJIbXajxEhHVNDexs66JzeDzqcdtP9nP+Dx7mdT9/mr2dg/x4y2Gqi52sqK9LanlU+hVUMEqlZIeYQg9E4fS5SJ9jYx3Ul5cmdN8l1ZXYrRZeGIx9TqP5ONI7zB93R5z8XhWAiWA0xxpiiVjX1BDoM7erPeJ+t9fHp/6yDbvVypaWbs77wcM8driLc5cv1tqePFCQo9KSOdpJP7RzSyr/XrlWG7Vj8Aj1jkrqiyqBwPIDR8c62LAgsefIZrWyrLaaPb2tmFqDRFzeMDDBY7zNaR37m6aMTNva2scbb3+a7hEPz9xwJafWFidUZpW7+txjWERw2pP/MbasthqHzcZDBzr553Oap+zzG8Pbf/0sTx/r4e3nnkFNWQl9I6MYv2Ft03wWUlXZomBqjJLxgajNXCpfdHsG+Ld9/8fNJ38/uW18kBGfm4bSsoTPu6K+hrb+QU5KoK9Sy1gnLWPJWayzY38TY14fNz6yh6t/+jij3sDC2Q8c7EnK+VVu6fOMUVJUlJLRW1aLhWW1Vexo656x0OvtW4/xt4MdvPGstZy1tJklNVWcsWghG5Y0xTyHkspuBROMkkkDUW4q1L9bpMd9d8dTeI2P7f1H8eIF4NhYoNmgvjz2uWCmW1lfA8DzwX5GNx76DZ86sAmv8SZ8znBfuOsENz66l9MaG/jElRdjs1jYfFyDUSHqdY9RHGHh1WRZUV9Dx+Awe7omh+37jeG/Ht7N8toqzl+R16vdFzQNRhFobVD2me/fI10d73PBqM/NfZ3P4LTbGHJ72OMPrLN4bDRQy9Mwx4i02TRVVmCzWNg5fAS3f5wDo610jg1y/8jTSSn7cx2tLKut4h0XnEVFSTFNVeXsau9j0D3OTU/u55cvHuXkYPzzz6jc0+ceS0n/opAVdYGQ/8CByeW2dnUM0D7k5oIVi3WeoTxWsMEo1skXw5vOlJpLtr1OIpXn/q4tDPnGeNNZ6wF4fmQPAMfGOnHZ7FNmEo6XzWphUXUFB/raOTByAp/xYxXh961P4jf+hM8L4PX72dXbxaKqyoltS2qqON7TzzcfP8AXH9rFh+/dyuU/fgKvf37XUtkv1TVGTVXlOGxWfvjMXr775H48Pj9PHAmEpCW11Sm7rsq8gghGscwPNNsK9dn2Yafik+6/XzbXOPqMj9+3P8Gymmo2LGmkusTFtv4jQGCofkN56by/CS+pqeJ4Xz9bh/cDcOWaVbT09/H42IvzOu+B/h7cPh+LqyumXMvr9/ODp/dzWmMdbzr7dE70D/P7Xcnp16SyV6iPUapYLRZeuXY1FouVLzy0i5dveox7drVSW+qiukQ7++ezgghGs8nmDzGV2zL9uop0/cd7d9Du6eXS1SuAQJ+g/V2deMQ9r6H64ZbWVuHz+/lz5xYqXS4uPXU5NaXFfGf/PRwPdsQ+um/uhTHdPi9HB/snbu/oCdy3aUowqgTA6ze84pTlnL20ieIiOz9+/ihHeoe5Z5cO589Hxhh63WOUJnkOo+lefspyPnbVxbz3orM52D3I08d6OGdp/GuoqdxS8MFI5bdsDCeZYozhzpOPUVdSymnNgUnoVtbXMjo+zuaxHfR5h1lQlviItJAlNYG1pzrG+llUXYHNauWDF5+DH8OnDv2QYV9sfYB+suclrr7vDvb0B2bS3t7dgctuo6Z0sg9UZbGLqmIXzVUVLK2rwW61snFpM1uOnuTSTY/ynt8+x507cmsaBTW3Ye84Xr8/pTVG4dY1LeBjV17Mx6+6mMvXrErLNVXmaDBSKsUyEY4iXXP38HH2j5zgklOWTawptSI4iuy+9mcBklJjVOZ0UFMaaGpYUhUISXVlpbzl3DPoGhniaV9sTWpbOlrxGcPXtj0JwPaeDpqrKmash/X+l23kPReePdEEeP7yxRhjKLIX0VBeyqf/so0Bd/QZjFXu6QtO7lhclLo+RtPVl5eysLI8bddTmZP3wcjtTt9/HJVdsqm2JhvsG2kBYM3C+oltFS4ndWUlvNhxFIC6isRHpIVbGqw1WlRdOWPbIU9L1Pt5fD48Ph/GGF7q7qDIauWZ4y38ZN9L7OnrntLxOqSxsnzKqur15aV85IqX8a9XXMSbN55O78gYV/3kSe7cdjwpjy1bdexvynQR0iY063UsS9coFa+8D0ZKZVq2zIh9YqwLh9VGmWvqqLOVwTXO7FbrlIAxH+uaGqh0OWmumfyGXeIooszp4NhIV8SZr0+ODHHlH2/nk5sf4shgPwMeN69atxpXkZ0btzxJdbGLC1bM3TcJoLm6AqfdztLaat688XR6Rr1cd9cL/P2wznmUD/b1BZpXq0u1E7RKvoJcEkSpfBatpuzEWBd1ZSUzRp2trK/l6YPHqC8rmdFMlajTmxs5vXlmORrKS2kd7MMYM6UcI95xPvTon2kbGaLr+Ajn1C0MlG1BHeXFLo5193LV2tU4Elj+4bzli1nf3Mh/3vsQm7Yc57JlOtQ61z3TfoIyR6CpVKlk0xojlTNysWksW2qLAFrcXdSWzfyGvaIuEBTS8SGzoLyM9oEhRqbNhP3DnS+wq7eLy05ZzrjPz6ZdL+Cw2WgoK+WMRY28bsOahEJRiKvIzrqmBh47cILhcd98H4bKIGMMz3a0sryuRidZVCmhwUjllFwKR9kUisb9XtrdvdRFWAet1OngtetP5cKVS1NejoaKUjw+HweG+6Zsf6jlECvra7j69FMoLrLTPjzMouqKpK5UfvbSZoY949y+7WTSzqnS7/jQAG0jQ6wKDhxQKtk0GCmVR6IFxzZ3D34MdSWRO1dfeuoKltZWpbJoQKDGCGBX/+QEjC1DA+zv72VtYwNWi4W1CxsAWBLWcTsZVjfUUlxk54+7NRjlsq1dgTX90vF6VYVJg5HKOblUa5QtTrgDSxnUlCenc3WiGioCwWhP/2Qn6EdaAyPiTg2Ollu/KPD3XZbkZResFgsLKspo6R/GTFsxXeWOHvcoAOUuZ4ZLovKVBiOVkzQcTQqfvT1a892JscAonrqy5AzHT1RxkZ1yp4PD/b0T2/5+4ij1pSXUBst26oI6PnTZ+ZzSWJf069eWltA5OMyYdjPKWf0eNwBOu07FolJDg5FSKZLq8BYKRKtXtE0JRNNvA+zsGKHYXpQV874srqlky7ET/HzfNtqGB9ncfoLTwuZWEhFWpKhjbW1pCUNuDy2DOuFjrhrwuHHZbUntf6ZUOB2ur1QMXuoNBIozqjxx3S8UXOZrviGrw99OXXlma4tC3nT26fz62Zf4ry1P8F9bngBgbWNDWq4dqjHb2THEqiodtp+LBjxuXGmc8VoVHg1GSs0hFIoSlYxwtHpFW1zhaPqxHf42VkQYqp8JpU4H77/4HLa3nGRkzIOzyMay+vSElFBz3d7OIThFg1EuCtQYaTBSqaNNaSovZVsfpFSWZ65ze4ybPn8vdSXZMxmeRYQzFjVywaolnLmkKW3z0dSWBMLhvu6htFxPJV+/1hgVJBH5mIjsFJEdIvIrEXFO2y8i8l0ROSAi20TkrESvpcFIJdV8a1ey0RlVnomf+chUWOv0dwBQV5o9wShT7DYrlcUujvcNZ7ooKkH9HndaF49VmSciTcC/AhuNMesAK/DWaYe9ClgV/LkOuCXR62kwUkkTCkUv9RalJSBFCxrZVlsULlVlCz/vjGY0X2DenkwP1c8WdaXFnBwYybsh+/WrTmS6CGkxMO6mWJvSCpENcImIDSgGWqftvwb4uQl4BqgUkYTecLWPkUqKSEHopd6iedeyzGV6/51kBo9EO1zPJVkdsmPV6Q9MiFeb4aH62aK2rIStx1oZHjeUFunIplwzoDVGWcfrttNzcMF8TlErIlvCbm8yxmwK3TDGnBCRbwHHgFHgAWPMA9PO0QQcD7vdEtwW95utBiM1b7PVDqUrHOWaVISjaM9Du+8kZQ4nznmsNZZP6kpLGB33crTfzdo6rUXLJW6fF7fPp52v80+XMWZjtJ0iUkWgRmgZ0Af8VkTeaYz5ZfhhEe6aULVwwTWl5WMfmEyK5fnMxec8k02Bydbpb4+4eGyhqgp2wN7fM5rhkqh4hSZ31M7XBecK4LAxptMYMw78Abhw2jEtwKKw283MbG6LScEFI5U88YSHXAxHIakqe7qa0zr8J6kr1Wa0EIfNCsCgO7+mv+7Y35TpIqTcgAajQnUMOF9EiiUwhPVyYPe0Y+4F3h0cnXY+0G+MSehNVuvWC8z0D/lEm7lyOejMJZ8e26gZZdA/QF1p/n9oxspuDQSjkfH8CkaFYGI5EA1GBcUYs1lEfge8AHiBF4FNInJ9cP+twJ+BVwMHgBHgfYleT4NRnpvrQz6dISAd/Y3mK59CEUBncERaptdIyyYajHLXRI2R9jEqOMaYLwJfnLb51rD9BrghGdcqqGAUPpx8PjUls9030gdrusNANn+450I4ylbxzn4NcNwXWLm+oUKDUYjdFuhBMKzBKOdoHyOVDgUVjBIR6YM8fBh3toSQbClHLLI1HEV7DlNR1nT1LzrsPUix3aFD9cOEaoxGNRjlnFCNkQ7XV6mkwWgW4TVMs+3PlExfP9/1cJxjvEj7WOAD1CUuziu6iCJxZLhksTvsO8CS+oq0LbmRCyaDkT/DJVHxGvQEvqTo1BMqlQr21ZXOWotkXCv8HPkQiLKt1ijSc/oUP6OVnYHpxIIe9N7DW4rez1r7GUm9fngTWbJqk0b8w5z0tXJm9alJOV++CAWjMa83wyVJnkIYkQYwOO6myGrFatEB1Sp1CjYYzSXZ4SPRIBBejnwIROGyLRyF8+Kmg71cvHIZr15/CgBHe/r4/ZYd/GDo25xdfBZvsr+XcktFQuefra9QskLSUd8hABZXVyZ8jnxktwY+VLXGKPcMjY9rbZFKuYJ+hWVr0MjWcqVCNoSjSM/3Sfbixcvq+gbswXlvVtbX8IlXvoy/7T7Iw3u2cmD8s3y0+D+ot6RuosZQSEokIB3yHkCARbXlSS5VbrNaLFgtgturfYxyzdC4R4ORSjl9haXRXCPWCikQhUt3OIrleT7BTqxiYXnD1Bohm9XKK9et5vTmBfzwkc3cNPYVPub8ErWW+lQVF0hsluwj3oMsKC/HqUObZ7BbrYzlUTCqX3WiIJrTBsc9ODQYqRQrmIbabA0doZXos7V86ZKuxx/rdU6wnSU1VVHfhBdWlnPdJefi9nu4yf1levydySzmvPmNn8O+gyyuqcx0UbKS3WrFk0fBqFAMjXtw2jQYqdQqmGCksl+2hMMxBuniCCvramc9rqmqgusuOZcR7yg3ub9Cr78nTSWcW4f/JKNmhCVVNZkuSlayWy24fdrHKNcMjru1BlSlnAYjVTBiDV6t7AIMqxtmD0YAi6or+eDLz2XQO8R3PV9mxIzMs5TJcdh7AIAldWUZLkl2slutedfHqH7ViUwXIeWGxj24tClNpZgGI5VVsqFZ8QQ7cFjtMXdaXlpbxdvPP5OO0V52lT+U4tLF5rD3IC5bEXXlpZkuSlayW614fPkVjApBYFSa1hip1CqIYJTpD1oVv1RMlxCrE2xneV11XHOlrKgLNFm1W4/GXbZUOOw7wOKaCiw6sWNEdquV8TyrMcp3fmMY1lFpKg0KIhip3JSJQDtIJwO0s7IuvlFmTruNCpeTdk/mO2GPmVHafCdYUl2d6aJkrSKbFY/Pj9+YTBclafJ9VNqIdxwD2vlapZwGI5X34qst2gHAKQviDxX1ZSV0jPRjyOyH7RHvIQyGJdVVGS1HNrNbLXh9Prx51P863/sYDY0HpvTQPkYq1TQYqayW7lqjE+ygtMhJQ2X8i67WlZfSOTiEKcpsB+wj3oMALK7TiR2jsVutjPt8ePOnwijvhYKRQ/sYqRTL+2A04tM+FrluPuEonvsa/JxgO6saahJadLW+rIQxr5c+Z2vc902mk/5WqlzFuHQF8qhCwUhH7OeOUDAqslszXBKV7zISjETkzSKyU0T8IrJx2r7PisgBEdkrIq8M2362iGwP7vuu6HLhBSUdNUc9HGeMQVbWJjaLdX1ZYARYu+VIEksVv25fJ9WlroyWIdvZrRbGfX6tMcohg55AMNLO1yrVMlVjtAP4B+Cx8I0isgZ4K7AWuBr4gYiEvh7cAlwHrAr+XJ220qqcFG+YOsFOAFY3JtY3pz44NP6ktCR0/3CD/gF2jG9N6L5dpoOakvibAgtJaObrfKsxyud+RkNebUpT6ZGRYGSM2W2M2Rth1zXAr40xbmPMYeAAcK6INALlxpinjTEG+DlwbfpKrLJBquc46uIw5Q4XVQnWtpS7nBTZrLR7T867LH8d+yO3Dv0vQ46OiPsPeQ8w6B+Ysd1j3Az4B6guLp53GfKZ3WrF6/fj8WuVUa7QztcqXbKtj1ETcDzsdktwW1Pw9+nbVQGKJRwlEqB6aWFBZeITIlpEqCsroWNs/kuD7BnfiQFOuHbO2Oc1Xr47eCN/sf9ixr5ufxcANcVaYzSbIlugInrI481wSVSsQk1pDh2ur1IsZcFIRB4SkR0Rfq6Z7W4RtplZtke79nUiskVEtowxGG/RVQ6IFnwSrVXy46ePEzSUzW8kV31ZKR1DgxjreMLn6Pf30eYPNIm0ycwJI9t8JxhnnJMR+jJ1BxezrS5zJnz9QmC3Bt76hsd1ksdcMTEqTYORSrGUvcKMMVckcLcWYFHY7WagNbi9OcL2aNfeBGwCqJPlWleep0IB6Iwqz5TbiRikHR/j8w5GdWWlbD3Wiru6B+dwQ0Ln2OfdPfF7m3fmy7zFFwhLncODGNc44pvsc9HlCwUj7Xw9G7s1UGM04smvTkb5PMnj0LgHp82GxaLjblRqZVtT2r3AW0XEISLLCHSyftYY0wYMisj5wdFo7wbuyWRBVfZIRt+j3mBLbWPV/Prm1JeXYIB2+5GEz7F3fBcuWxHNVRWcHO2esb/FdwyA3pFRPM6+Kfu6/Z0UWa2UOnQZnNmEglE+1RjlcyiCYDDS/kUqDTI1XP8NItICXAD8SUT+CmCM2QncCewC7gduMMaE3rk+BPyIQIfsg8Bf0l5wlbd6gl3bFlTMb9HV0JD9k3IsofsbY9jr3cmK+moaK8poHxjE2NxTjgkFI2Og0z71Ol3+TqpLShKah6mQTAQjT/4Eo3w3ND6OQ4ORSoNMjUq7yxjTbIxxGGMajDGvDNv3VWPMCmPMKcaYv4Rt32KMWRfc9+Hg6DSlkqKXFqpcJfN+460rLUGAdpPYsOkufwc9/m5W1dXTUFHGoNvNgKttYr8xhhbfUZqrKgDonDY1QLe/k+oSbUabS6jz9Uge1Rjlu8Fxt9YYqbSIORiJSEnYnEJK5ZUejrOgYv4juew2K5UlLjo8XQndf693FwCrFlSzIDgvUpv14MT+bn8nY2aMMxYGmk3amQxgxhi6/B06Ii0GtmDn61ENRjljaNyjHa9VWkQNRiJiEZG3i8ifRKQD2AO0BWes/qaIrEpfMZVKHT9e+mmbd8frkPqyUjqG+zHE37F37/guKpwu6sqLWVBRBkArkyPTQs1oKxqqKXEU0Tk+GcCGzRBu46Zag9Gc8q2PUb73LwIYHPfg0skdVRrMVmP0d2AF8FlggTFmkTGmHrgYeAa4UUTemYYyKpVS/bTjx5fcYDQ4jM8xNOtxJ7zH6PVPdq72Gz97vbtYGVyrrcLlxGm30eabHJnW4juGIDRWllJXWkLnaN/Evi5/YDJIXQ5kbkUTo9LyIxgVgkHtfF2wROQUEdka9jMgIh+ddsw7RGRb8OcpETkj0evN9iq7whgzYzIWY0wP8Hvg9yKi8V3lvN5gx+v5jkgLqS8vYdzno9dxgjp35LA1bjx8Z+jrFNmsfLbka5RKOSd8xxk2Q6yqDVTGiggN5aWcHJmcMLLFd5T6slLsNiu1ZSXsO9mFKXMjXsfE5I46h9HcJobr50mNUSEYHh/XYFSggitlbAAIduk5Adw17bDDwCXGmF4ReRWBKXvOS+R6UWuMwkORiFSJyHoROSv0M/0YpXJVDy0IQsM8R6SF1IUWk7XOnJwxZNv4i4yYEfrGB/klNwVHo032LwpZUB4cmWYfA+C47yhNVYGwVVdawsDYGGPOQK1Tl09rjGIVmuBx1Jtf8xjlK5/fz4h3XJvSFMDlwEFjzJQ3WGPMU8aY3uDNZ5g692Fc5ozfIvIV4L0EhsiHRoIZ4BWJXlSpbNLLcaqLS7DbkjO2oL4s0MfnJMdZF+WYp92PU+ks5mWrF3Pftj08WnMPe8cPUldSRkXJZI1PQ0UZmw8fp9/Vis1dR5+/l4XlCwGoDV6n3X6UpSyky99JaZFDO6jGIFRjNDquS4LkgmFv4Du4U1/bKrDQ/K/mOOYDzGNKn1heZf8IrDDGeBK9iFLZrJcWmiqT12G5zOnAabfR7m2PuL/P38Me7w4uX7mSl5+ynAOd3dzVfg9irJzTMLUTbWhkWqvlABbfKAALK4M1RsFg1MEJlhIYsVZTrovHxiI0XN/t1aa0XDAYXA5Em9Kyk9dtp/1w43xOUSsiW8JubwquYDGFiBQBryfQ9zkiEbmMQDB6WaKFiWW4/g6gMtELKJXNfIzTz0kayiqSdk4RCXTAHu2NuH+z+0kMho1LFmER4a3nnEFxkZ1xM86q2ropx06OTDs2sRRIU00gLNWUBkJQhwl0zu7yd+gcRjGyWixYRHS4fo6YWCdNg1G+6jLGbAz7mRGKgl4FvGCMifitU0TWE5gI+hpjzMxlA2IUy6vs68CLIrIDmJiC1xjz+kQvqlS26KMVg5+G8rKknreurIQDHd2YUg/inVyewxjD057HWV5TQ21FIMSUOh2884KzeHDHflY1VU05T5nTgavIzklvG16flwqni5Lgch8Om40Kl5PO8W58xkevv4cNxSuT+jjymd1qwa19jHKCBiMV9DaiNKOJyGLgD8C7jDH75nORWF5lPwO+AWyHBCZmUSqLTayRVpncJqj68lKeP3qC0dpOiocmm8cO+fbT6W/nsiUbphy/vK6af7ps5gAKEWFBcGTamG+UhTVTA1xtWQmdI/30+Lvx46e6WJvSYmW3WrUpLUcMejQYFToRKQauBP4pbNv1AMaYW4EvADXAD4JLInmNMRsTuVYsr7IuY8x3Ezm5UtmuhxYsCPXlyZ0UcWLNNNsRljMZjJ5xP06RxcYZSxpiPldDeRlbj7Xi8fWxtnxqjVBdaQnbW07SXRKYAbtGR6TFTINR7hjSPkYFzxgzQiD4hG+7Nez3DwIfTMa1Yulj9LyIfF1ELpg+XF+pXNdLC7WlpdisyV3tJhS0TjK5yKvbuHlh/FnWL2qM65vvgopSxrxe/PhpqpzaF6q2rIRhj4dj7AegulznMIqV3WbF49NK8FwwEYxsOlxfpV4s785nBv89P2ybDtdXeaGX4yxN0vxF4WpKSrCI8HTvi1R717PadhpbPVsYM2Ocs3hRXOcK7/+0sHpqWetKAwFst2c3luBs2So2dqsFj9YY5QStMVLpNOerzBhzWToKolS6eXEzQAf1ZauTfm6b1cJVa1fxyJ5DfM/731TYyrD5HdQUl7B8QWVc5wqNTHPabVRNayoLDdk/OHKUKlcxVkvM60IXPLvViseX+8GoENZJGxr3IExOs6BUKkUNRsF10O4wxkSsaxaRFUCjMeaJVBVOqWQzGDyMMMoAHRwADA3lyVkjbbor1qziktXL2dXWwQtHWtl9sp3XrDiNYMfAmJU6iihxFFFfVopl2n2rS4oRAZ/x6YzXcSqyWnF7vfiNmfG8quwSWict3v87SiVithqjGgLD9J8Hngc6ASewErgE6AI+k/ISKpUkW7mHLfwOP1NnO15UnfymtBC7zcoZixo5Y1EjXp8fqyX+N3YR4TXrT6XCObOZzGa1UF1cTPfwiI5Ii5PdamXI7cFnIIE/i0ojr6sPpy4HotIkajAyxtwkIjcT6Et0EbAeGAV2E5gn4Fi0+yqVjQ7wFLVlLs5dsoRSZxGlziJqy5zUlqUnUNisiTdznbsser+k2rISuodHqClOXcDLR3arlXGfD68f7NoCmdUG3LqArEqfWV9pxhgf8GDwR6mc5WGEHo5zRdNqLl2zJNPFSaq60hL20kl1qdYYxcNutQSCkZn7WJU59atO0PGYm+IirTFS6aHfk1RB6Aiugby0pmrOY3NNqMarusyR4ZLkFrvVyrjXh47Yz171q07g8xt2tA9MrBGoVKpp3aQqCO0EZohfWp9/b67rmxvpGRqjqSr/Hlsq2W1Wxn3+nK8xql91Iq9Hph3sGWJ03Mfi+rq5D1YqCbTGSKVNO/t5ktvwMJqBa+9jQVlFXnbgLHc5ef2Zp+lQ/TjZrRY8Ph/jvhxPRnmqflVgNveX2voBaExTX0Cl5qwxEhEH8EZgafjxxpgvp65YKp8YDLt4kKf5OX58mLITXDT4OYT0DAUy+OlgP2fV1qfleio3OGyBt7NBj5fGkqI5jlaZ8lJbH3abjTrtQ6fSJJavmPcA1wBeYDjsR6k5efHwCLfwJD9ldUMtl65ewa7BHRwq+WPaytDLCTyMsriqZu6DVcEIdebtGhnPcEnUdKHaIoCdzsUsrCzXGlGVNrH0MWo2xlyd8pKovPA4P2I/T1As5RRTxagZpJ82rjztFK5ctwIMnOjr47HuO6mxraXSuyLlZQr1L1qWh/2LVOKKiwK1RJ3DHiC5iwinSz73LQKo+NLd7G5sZHG5TkWh0ieWCP6UiJye8pKonDdMD3t5hOaaEpY1FVFS20151Qjvu/A8Xnn6SiwiWCzC288/E4fNxt+LvsV4GvobtbOfErsjbfMVqdwwUWM07MlwSRKTr6EoVFtU8aW7GR0d5eTJk9SGrReoVKrNtiTIdgKLxdqA94nIIcANCGCMMevTU0SVK/bwMH58vGXjBuoqoi9PUeZ08I7zN7Dp0c08V/4dLhz4bErL1c4+ltZW6nICaopiR7DGaCT3glG+hqLpjh0LzCNco/2LVBrN1pT22rSVQuU8P1528zdWN9TPGopCVjXUctmpK3h4zzZOLX6B6pGzUlKuMQbop40Lq9em5Pwqd5UEa4y6czAY5avw2iKAw4cPA1Du1M7xKn2iNqUZY44aY44C/xX6PXxb+oqocsERtjBCHxcuWxbzfS4/bSUljiK22H+BITVDptvZD8DSmoqUnF/lrlAfo97R3ApG+V5bFApFMBmMqiOsE6hUqsTSx2jKV20RsQJnp6Y4Klft5AGqnCWsaY595JfDbuPKNSs50n+SjrLHUlKudvZhEWFxrXa8VlPZrBYcNiv9ozoqLVuEhyKAI0eOYLNaKXPprO4qfaIGIxH5rIgMAutFZEBEBoO3OwgM4VcKgB6O08ZuLli5GEucy5Sfv3wxVcUunuNODL6kl62dfTRVVGK3WZN+bpX7XEVFDLgzX2PUsb8p72uCEnH48GGqS4uxaP9AlUazNaV93RhTBnzTGFNujCkL/tQYY1LbW1bllF08gFWsnDPLCvDR2KxWrl63mtbBHk5U3J/Ucvnx0skhltTm3/poKjlKiuwMjmVPjVEoIEULSfkcnlbd/v0pt9vb23nhhReodGkzmkqvWOYx+pyI/APwMgKj1B43xtyd0lKpnOFhhP08wYZFjZS5Eltu48zFTfx97yE2j9/NQq7CQnKW7ejmKF48LKmuTsr5VP4pdhQx7PFgjIlr1OLhnmEayhwU2+e/3ORcISh8ssNc9FJXOw+0HGJPbxftI8MsKa/ga+ddRkXRZOCZHooGBwfZsGED3V1dvO38M9NdZFXgYulj9H3gemA7sAO4XkS+P/tdVL7q4TibuYOn+DlP8BMe4ibGGeOC5bF3up7OYhFete4UukeGOFr1h6SVdaLjtfYvUlGUFNkZcY/j8cd+n9FxHxfd+nc+/+DeeV8/lhqgXK8l+vDj9/OTPVs5PDaAtdjGwyeOcu1ff8udB3by0z0v8ZrnHuDVV7+S3bt3T9znl7/8JSdPnuS6S89n/UJdykelVyxfdy4B1hljDICI/IxASFIF6Gl+Tis7KbLZsIpgtVhYX9s471Xr1yysZ0lNJc8O/5Ul/ENSao1Oso9KZzFVpXNPH6AKU3FRESMeD24fOGLshvZCay8j4z7+drAz7pqm+apfdSKngtLJkSHaR4d5w4a1XLR6KQCHOnv45dMv8O/PPgrAsroa/v7ww6w//XTu//OfeMWVV3HLLbewqLaapVX6pUalXyzBaC+wGDgavL0I2JayEqmsNUwPJ9jJFaeu4ur1q5J6bhHhFaeu5KdPbqGt5i80db9+3ufsYB8r6yrnXziVt4qL7IyMjzPiNZQXxRZwnjnWA0BL7wAnBr00lyce4mMJOtOb0nIpHG3v7gBgUfXkdBnL66r599ddTvfQCF6fj8bKcgbH3Hz/4ad50z+8kY/+27+xfft23nL+WTopq8qIWJrSaoDdIvKIiDwC7ALqROReEbk3paVTWeUATwKGMxen5k35tIX11JYW85L3AQxxtG1EMEIfQ3TTXKH9i1R0xY4ijIlvWZBnjnVjEcFvDH871D2v69+zq5W99pci7qtfdSLn+xdt7+nAIsKCyqk1PxYR6spKaAxuL3M6eM+FZzEyNsaXvvQlVjTUaROayphYaoy+kPJSqJywn8dZXFVNQ2Vqpue3iHDxqmXc9eJO+sqfo2rgvITP1c0RAJqqdI0lFV1o9uv2YQ/r6uYe/eQ3hs3HezhjUSMvHW/jsaM9vGfDgoSv/8WHdlLmKubx/1cbc2frXKktgkCNUWNFGUUxTJfRWFnOhy+/EAM0VWoTmsqcOWuMjDGPAkcAe/D3Z4EXjDGPBm+rAtDNUXo4nrLaopCNy5px2e3stN89r/N0BYNRc40GIxXd5HppsQ3Z390xyIDby9rGOhZWlvP0kZP87PkjnBwci/vaxhjaBsdoGxhl3G9iqiHKpVBkjGF7TyfNVbHPOr+wslxDkcq4OYORiPw/4HfAD4ObmoG7U1gmlYX28wQWsXDm4sS/HcfCYbNx/orF7O4+wqjrcMLn6eII1cWluIqSM/Rf5afi0Hppw+6Yjt98PNB0trimmtUNtbT2DfKR+17io3/eFfe1+91e3F4/fSOjDHjmXhInl0IRwLGhAQY8bpZUV2a6KCoPiEiliPxORPaIyG4RuSDKceeIiE9E3pTotWLpY3QDcBEwAGCM2Q9o428B8ePnAE9y6oJ6Sl2pX8zxopVLEBH2lv024XN0c5jmaq0tUrMLrZfWFawxuv7uF/jKw3uiHr+3cwiX3UZViYur1q7iE1ddzIZFjTx6sI0hT3wzt7cHa5m8fj/H+mevccq1UASTHa+bqnWdQpUUNwH3G2NOBc4Adk8/ILhk2TeAv87nQrEEI7cxZqJnoojYIEUrfqqs1MoORujlzEXxz2ydiMpiFxsWNfJSz3bGbX1RjxulnwE6Zmx3M8wAHSwsr0xdIVVeCPUx6h7x4PX7uWvnCf6wu4Pg7CQzDHrGcdrtiAg2q5XGynLOXtrM6LiX3+7q4HDPMH87MPM1GUl489uB7tGox+ViKALY3deF1SIsqNAvKGp+RKQceDnwYwBjjMcY0xfh0H8Bfg8RPhjiEEswelREPge4RORK4LfAH+dzUZVb9vMETmsR6xbFvkDsfL189XLcXi+Hq34Xcf8xXuBOPskfrZ/Db5n6bbs7OLNEUxx9G1RhchbZEYHeUQ/7uoZwe/30DI/ijlL5M+j24rRP7Ui8qr4Wl93GHVuP8cbbn+btv9nMcAwzRp4cmnzdHu4bmdfjSJbHWo+xq7crKec6MTRIVbELqyWWjxmlZrUc6AR+KiIvisiPRKQk/AARaQLeANw634vFMirtM8AHCEzq+E/An4EfzffCKnt5GGGMQdwM4WaYwzzLWYsWYLembyHW5uoKVtXX8GjH3+iyeTnb+y4clODHxxbuZCv3UuZwMugeo6vu79R3vmrivl0E+iY115SmrbwqN1lEqCkp5ljvIC+19QEwMDpG75ifxtKZr/cht5ci29S3TZvVwrqmBTx3pGVi25a2IS5ZMnsn4vAao6O9mQ9Go95x/uWJ+6l0OnnoNe+Y9//31uFBKot1ctVC4HHbOXpkXjWbtSKyJez2JmPMprDbNuAs4F+MMZtF5CYC2eQ/wo75DvBpY4xvvvNfzRmMjDF+EbkbuNsY0zmvq6msNkgnz3EnB3hixr6NixenvTzvuehs7t++nycPPMoh63Oc7X8rB8xTtLGH85Yu4bVnnMKN9z/CAXmMeqYGo3KHizKnI+1lVrlnRV0N206c5PkTgWkoDHCob5TG0pnBesjtxWGb+ba5YfFCnjvSwprGena1dfBi20BMwchhs2IRoXUgelNaujx84gijXi+jQ0P84shO3r9i/bzOd2JkkGUNOo+YikmXMWbjLPtbgBZjzObg7d8RCEbhNgK/DoaiWuDVIuJNZG3XqMFIAmf/IvBhQIKbfMD3jDFfjvdCKnu5GeZF7mYn9yNiuHjFchZWVOIqsuJy2KhwOagtK5n7REnmtNu59qw1nLOsid9t2cHjvT/BbrHx1rPOZOPyhQBsWNTI5kPHOLeoG5sn0NTXzRGaq3XIr4rNivoaNh8+zu93nJiYuPFQ7wgXNc8MRgNuL84IgXt1Qy03XHYBCyvL+fe7/sr29oE5r3tyaIwKlxOb1UrHUOaD0Z+OHqDC6aTUWcSPdr7Au5auTbjWyOPz0Tk6wsaKpcktpCpIxpiTInJcRE4xxuwFLicw2XT4MRMLdorIbcB9iS54P1uN0UcJjEY7xxhzOHix5cAtIvIxY8z/JnJBlV2O8xIPczNuhjh7cTOvXHsq1WXZVdPSVFXBv1x+IdtOtLGgrJIFYRNMnrWkiScPHKW18q8s7ng744zRRytnVazOYIlVLllRFwjUvaMeVjfUsq+9iyO9kYPKoMdLRenMt00RYVldoHakprSEQ92Dc1735OAYZU4HTrud7uER/MZgSfESGNGuMeBx82jrUS5YsYSFVeX8+tmXeLb/JBdVJ9Y80jocePzVpen/QqXy1r8At4tIEXAIeJ+IXA9gjJl3v6JwswWjdwNXGmMmeuIZYw6JyDuBBwANRnngWX5FaTF86PxLaK7N3j45FouwYdHCGdsXV1dSW1rCPu9mFvN2ejiGwdBUoR2vVWwqip3UlpbQNTTM+uYG9rV3cSxKZ+ghtxenffYeCAsqSmnpG8TnN1gt0YPOyUE3laUVlDiKONDRxZDHUO5ITTA60N/Dd7c9x0Mth/jvl13BaxdNXevwweOHGPf7OXPxQqzWQGfpXf1dCQejEyOBYFThSP30HqowGGO2EmguCxcxEBlj3jufa802XMAeHorCLtgJSVj6XGVcF0fo5igvW7ksq0PRbESEs5Y0cbinA3fJ4YkZrxfl6ONRmbGiPlDbs6i6ijKng5ODM2uMjDEMe+YORo0V5XQNDXNyOPq8RsaYQFNasYPq0mLcXh/HBmZOMpmsofq37Hiev7UewWq1cseRqZNRenw+frjrBRrKSmmurqCutAQB9g/0JHw939uvAaDCqcFI5Z7ZgtFsqyrGvuJinths7810EZJuH49iFQtnLU3tbNapdvaSJgxwtOwvdHGEYruDyuK5171SKuTcZYs4o7mRBRVlVBY76YzQ52fU68NnDI45a4zKMAZePBm9n9Ggx8vouI8Kl4Oa0kDT8J6u1I1M63WPsbCyjNObF7CzvYNRn3di3492v8iRwX6u2bAGEcFus1JdUszxgf6I5xr3+/jx7q2Meb0R96+6/fscPXoUEdH/hyonzRaMzhCRgQg/g8Dp6SpgJm2290785BsfXg7wBGsXLqAkx7/V1ZQWs7S2il3DW+niME1VZcx3uKYqLEtqqnjXhWdhtVioLHbRMzyKzz91kschdyAIOOZYELUxOKHh1rbo/YxCQ/XLnU5qSgLBaH/3cMLln0ufZwxXkZ3TGusZ9YzzePdxAI4N9XPLzuc5o7mR1Y11E8fXl5fSOjCIP8JEl891tPGNF5/irhP7ol7v6NGjVJboHEYqN0V91RpjrMaY8gg/ZcaYvG1Ky+cwFO4ozzPGEGdnYBh+Kpy1uImOwUG6OUJTZWWmi6NyWFWxi76RMYbGpwUjTygYzV5jVFNajM1iYVdH9BqjUDAqdTomgtGhnlQGIzfFdjurG2qxiPBQ6xGMMXxly+OICNdsOG3K8fXlpXQODjPon9k40BbsWH1odGaN0qrbvw/AkSNHqHLpHEYqNxVUnI8WeGIJQ+eNV6WyaCnlZ+YsvPt4lDKHi9Oa8mOekQ2LGie+nS6sqMxsYVROqyp24fH5aJnW52cgVGM0R1Oa1WKhvryUw1FGpt23p43fbg9MBlnmdGK3WSl3OWjpT11TWr97jJKiIlxFdpbUVPL08ePcsX8Hj7Ye46q1qymfNhFjQ1kpXr+fA0Mz3xPbRoYA6BiZGuRCoQgCwUgnd1S5KiPBSETeLCI7RcQvIhvDti8VkVER2Rr8uTVs39kisl1EDojId2UebSXhQSiWmqFcCkV+vHRwkB38lYf5Pr/hY/yEd7HF/jNMMCAN08NxtrJxSTOWWUbN5JJiRxGnBZsCmqt1iLBKXOgD/eC0IfuhprQi69wLBiwoL+XkwBAe39Rap5ODY7zzN8/yixeP4bLbqAj2wakpKaF9cCTqGm1zOTrYz+v//BsO9M/sMO3z+xkc91AcXBfughVL6Bga4j+3PM7CinIuWrVkxn3qywODF/YMdM/YFwpG3aOTQS48FHm9Xk6cOEFtuQ6AULkpliVBUmEH8A/ADyPsO2iM2RBh+y3AdcAzBJYluRr4S6oKmKv+zI20shOAMoeLxTUVCA280Ho/g85WLh77BPt5AoNh49LmDJc2ua5Ys4pKl4vacg1GKnFVJYFgdLh3BJisUQ01pRXNUWMEgaA+Nj7OsBeKwroktQfXR3vrOes5fVHjRLNcTWkx+9o7cfvAmcC78s3bn2NPXzd/6TrMv1RMrQUeHA80hxUXBfoSnrWkieV11Ww91srahQ0R+wGFgtG+gZlfHE8Gg1Hv6NiMfQBPPPEEPp9P5zBSOSsjwcgYsxuIuYOsiDQC5caYp4O3fw5cSxqCUS7VFvXTRis7edmK5VxyylIqS5yICMYYHt1bwX3btjHs+DxD7nGWVNfQEDZRYj5orqqgWReOVfNUFawxOtoXucbIOUcfIwh00HZ7fcEao8n3ue6RQEipLi2Z0leppqSYgVE3HSM+FpcHtsc6VP/IQB9/PLofgPbRmf2U+j2BJsFix2TX0MpiF5eeuiLqOYuL7JQ5HRzt75uxr3U4EIz6R0fxGT+n3nHLxL7Ozk7e+c53Ul9ZwbqwztxK5ZJs7GO0LLh67qMicnFwWxOBtVJCWoLbUiqXQhHAYZ4F4OWnLKOq1DURPEWES09dwbsvPIuO8VYGaGfjkkWZLKpSWau4yI7dauHEtPXLBmPsYwSBDtp+YxgZnzqXUddwIKSUTpv4MDRkf29X/B2wb9u7DavFgt1qoWt0Zj+lPnegZsdVFN+YmbqyEk4ODc1o3gvVGPWPulnw7S9M2XfrrbfS2trKuy86G+cco/eUylYpC0Yi8pCI7Ijwc80sd2sDFhtjzgQ+DtwhIuWEf+WaFLUxXkSuE5EtIrJljLmn5k8XH15a2IaJXvR5OcRmFlVWUV0aee6Q9c2N/PNlF3DBsqWctaQhJWVQKteJCFXFLtqnTfIY66g0gKLgMX1jU+f6CdUYlUQJRi+09gGx1xb5jeGB44dY01hHVbGLntGZ8y+FaoziDUblTgcDY268Ye9Xgx43w95xShxFuL1eTh4/NuU+W7Zsob6ijAUl2vFa5a6UBSNjzBXGmHURfu6Z5T5uY0x38PfngYPAagI1ROEdYpqB1lnOs8kYs9EYs9FJWXIeUBIcZjN/5ut0l26e++A4DdJJF4c5valx1uMW11TyxnPW4ijKVPcypbJfZbGL7uHRKbUloRqjohhqQkJzHQ24p9YYdY94sIjMCCkLK8upLy/l249u5+4nYv/itL27g66xEU5vWkC5y0nf6NiMGp5+T7DGyB5fMCpzOhgcc+Nm8jGEOl4vqa4E4MEHH+QTn/gE3uBkjy++8AILy7PnPVepRGRVU5qI1ImINfj7cmAVcMgY0wYMisj5wdFo7waiBqxs1UcbAC1lf0v6uQ8TCFvrm7UmSKn5qioJzGU0FpZrBt3jOGzWmBZ6DdUYhYb4h3QNuylx2Gecw261cv0l51FZ7OK6R/5EazCAzOXBlkNYRVjdWE+5y8nA2Bjj06bn6HMH+xjFWWNU6nTg9nrp905OWxAKRotrA90Mvvadm/j2t7/NT753E93d3RxvaWFRXX5MAaIKV6aG679BRFqAC4A/ichfg7teDmwTkZeA3wHXG2NC408/BPwIOECgJinlHa+TPcnjAO0A7O07ADKe1HMf5lkWlldSW5FfHaqVyoSqYheDY266RyeT0VAM66SFhI6bHoy6Rz2UFEWeab7c5eQNZ63D4/PxRF9LxGOme6jlMCvqayguslPudNA/6mbUTL3mRI1RnMGozOkA4OToZEibCEbBGqOTnYHlNG/+7nfZunUrAE2V5XFdR6lsk5FgZIy5yxjTbIxxGGMajDGvDG7/vTFmrTHmDGPMWcaYP4bdZ0uwKW6FMebDJtEJPzJokA4sIvSPjjFQ/WzSzjtEN+3s5/Tm3F7zTKlsEZrL6EDPZGfmIY83pv5FMNncNhSxxij6EjyhEXEnhufuG9nrHuPQQB+nNgRGf1W4nPj8fjrcUztg93vcOO22uJfnCAWj9rHJDuGeqy5GRFgUNvrz1MZ6th85xo9/GJh9ZUGpfjlTuS2rmtLS4bzxqhk/s0lmrdEg7axpbMAqwtGix5N23iM8B8DpzbP3L1JKxSYUUA6HTfI45I49GDmiNqV5Zg1GFa5AGGmLoSnt2GBgSY664JxD5a7AoIsTo1Pv2+cem5jDKB6TNUaBoLXq9u9z/PhxKopduIrsOO02qopdvHnj6ditVn71299SXVoy6+NTKhcUVA/caCEoHcPyvbgZoZ+mioV4/T729u9nHT6E+Q9pPcxmGkrLWZBn8xIplSmhGqMjvZO1LwNub0yTO8JkMBryTBuuP+JhTWX04GCzWilzOiIOu5/uaDAYVZcGyloeClXTglG/xx13/yKYDEYDC2snZrZua2ujPLj9vGWLWFBZToXLySeuuphjfQPUa8drlQcKKhhl0gAdAFSXFFNesoDfbtnOcPULlPacM6/zjtBHG3u5onlVMoqplCLQLCXA8f6pNUZOW2zD0ENNaYNhNUY+v6Fv1EOpwzHntbtGYghGQ/0IgfcUgHJnoMZoem1Tv8cdd/8imJxrqb2jY2JbV2cnJcFzvW7DmonttWUl1JbpTNcqPxRcU1qmhDpe15Y7Wde0AIsIxxyPzfu8gWY0o81oSiWRzWqh3OWkLWwuo8E4Ol+HaoyGxyeDUe+oBwOUOmcPKRXBYfe+ObpRHh3sp6rYhd0aCGGhGqOOabNf93vGEqoxsloslDiK6OqZXH+to6ODUlfkedKUyhcajNJkMFhjVFtWTImjiOV11ewd3APin+OeszvMs9SWlLGwSr+tKZVMlcVOuobCgpHbG/NsznarBREY8UwGo66RwLD3kjlqjCqLnfSPjDGGd9bjjg32Ux3W0dlutVJcZJ/RDNfndlMc5xxGIWVOB30DAxO3u7u7KdNgpPKcBqM0GaADh9U+8c3t9OYFdA4NMVKxPeFzjjFIK7s4vWlBzOvOKaViU1XsomdklHF/oOYmnuH6IoLDZmMkrMaoa3jqYq7RVLhcjI6P0zMeeZHWkKND/dROW6i13OWkd2zyfsYY+j1jlJYk1v+wzFFE39AIxhiGh4cZHRujxDl7sFMq12kwSpNBOqguKZ4IMKc3LUCAY65Hot5njAG8uKPuP8jTGPyc3qTD9JVKtsriwCSPR/bX07JnAeM+f8zBCAL9jEbD1koLLQdSPMeorcriQI3M8ZHoQ/YHPG563WPUTw9GTsfE7Nerbv8+jT/8Bj5j5rxmNKHZr8eGh+jqCsxZNH2dN6XyjQajNBmgnZqyyY6b5S4nS2qr2D+8i0jLvp1gB7/mo/zR/lm8zPzmOEA7z/IrllbXsKhWR4IolWxVxS58fj+to4MMewMTshZZY2+ScthsjIXVGHVPNKXNUWMUDEYtIwNRjzk2FBiRNr3Dc7nLyeCYm6YffBWAnmD/oET6GEFg9usht5vh3l46OzsD5U+wWU6pXKHBKA0M/kCNkat0yvb1TQtoGxjgZOnUOY328Rh/4UZKXVa6xtt4qvRbUxae9ePlYW7GYvXztnPP1GY0pVKgKrgQ6pHh/okOzaXO2GtLHDYboxGa0qLNfB1S6Zp7kscjA4FgVDVtMsUKl4OB0TEGugIhprU1sKRkaQLzGEGgxsjj9dHRemIiGMXaz0qpXKXBKA1G6MOHl+qSqd/uzlzSRGWxi3uHbuF+23/SyUGe5/c8wi0sr63mI1dexNWnn8KeoZ3sLvnVxP228Ds6OMCbz15PTbm29yuVCqH+OwcHeycmU5weRGZTZLPi9vomFnXtHvHgstuwWWd/251rksfdb3klX9+/lWJHEbXT+g6VO534jeHYoYOBsh8M/FsVR6Cber5AWY4dOTTRlFaSYO2TUrlC5zFKg9AcRrWlU0dzlDkdfPrqS3jywFEe3nOQu/h3AM5etIg3nbMOu83CK05dwYneAZ5s+SOVjuUYdylbuZdzlixmw1JdMFapVAn0CYTDA324nIFakulBZDYOm42BsTFO7m/CKsKJ9l1zzmEEgUkeSx1FE6PLQpMrAvj9fs4sL6fSUcT7r7wY+7Tam9Ds10cOHeQ84MCBA4hAdUls8y9NVxoMRieOHqPfFvhd+xipfKfBKA0m5zCa+eZkt1m59NTlnL9iEY/vO4pdbFxy2pKJ5jER4S3nrqdzcIiHhn+AhRJqS8q49sw1M86llEoem9VCdUkxJ4YGKfPaKbbb4+rE7LBZ8Yz78OHHipUe9ygljthqWyqKnfS43VNCEQSGyw8PD3P12lXUFc8cNh+ay+j4kSNAoMaoqqQEmzWx5q/Q7Nctx4/hLqvCarHE1QFdqVykTWlpMEgHgkwsMxCJ027nyrUruXTN0hl9hhw2G++9aCMWC7gtA7zzvDNxFGk7v1KpVltaQvvgEMcG+6mJc3FUh82G2+vFG+wf2DM2GnOwqnS56Bkaxhfs9B0S6jNUFmXIfEWwxujY0aNAIBglWlsUfp22YB+jEqdD+zSqvKfRPw0GaKfC5Zqzb8FsakqLueEVFzA85qW5tnTuOyil5q2urITDXT1YvFBdGd8kqkU2G26vDy+BSVzbR4c5tbo+pvtWFDs51NXDUE8PFfWTTeYTnantkb8YhYJMa1sbAPv372dlReKTv5YUFWERofVkO1a/dc6O40rlA60xSoNBOqiN89tmJA3lZSyvT/2Ct0qpgLrSEjxeH8eHBuL+P+ywWfF4vYwbP0PjnsC8Q9WVMd230uVi1DPOiaOHp2wPBaPyKDVGVouFUkcRnV3dDA4O0tXVRV1leVzlDmexCI2VZew5fCRQY1Sk36VVZojIERHZLiJbRWRLlGMuDe7fKSKPJnotfZWnwQDtLCmpzHQxlFJxqgubJ2j6LNNzcdhtGKDiIx/E7XTBb380seDrXBZWBuYme/LRR1hzzvkT29uCNUHRghEEOmB39/dPjEib75eyFXU1PH3wKG5bEVUxdB5XKoUuM8Z0RdohIpXAD4CrjTHHRCS26tkItMYoRpvtvWy298Z9v3HGGGWA6mJt/lIq14RPoFhdGl9fnaLgiLG+3m4OHToEQGWMna8XB2uWnnhk6pfe1tZWSp2OWTtTl7sc9A+PsG/fvuA159f8tbyumnGfn2PHW3SdNJXN3g78wRhzDMAY05HoibTGKAaJBKKQ0OKxsX5TVEplj0qXC5vFgtfvpyrO/8MOW+Dttb+vj0MnTgJQPcsAjHDFjiLqykrYumPHlO2tra0TQ/KjqXA6ae0dYM/OnQDUzKPzNcCy2uqJ30s1GKkIxjw29h+tmc8paqc1j20yxmyadowBHhARA/wwwv7VgF1EHgHKgJuMMT9PpDAajOYwn1AEYXMYlWkVtFK5xmIRakqL6RoanpiROlaOYI1Rf18vhw4dosThwBXH5IhLaqrY09bB2PAwzuDksK2trZTNUQMUWhbkuWefpbKkGOc8l/AocRTRWFFGW//gxLxGSiVZlzFm4xzHXGSMaQ02kT0oInuMMY+F7bcBZwOXAy7gaRF5xhizL97CaFNaBKFms/mGIgibw6hMa4yUykULK8upLyvFYolvmPqUGqNDh6iOs6/P4upKhtwetjz2yMS2EydaqJij1qnc5Qh8tX74b6ysq5712FgtC55HZ71WmWKMaQ3+2wHcBZw77ZAW4H5jzHCwH9JjwBmJXKtgaow223s5b7xqyu1knGcuA3TgtNnj+qaolMoe1565lnGfL+77FQWD0eDAQGA+oQgTMs5mSU0lAH/761/o8xnOOuss2ts7WLdm9az3CzW1eTzjnLZoYdzljmRFXQ1PHTiqC8iqjBCREsBijBkM/n4V8OVph90D3CwiNqAIOA/430SuVzDBCObfLDYbg0GY+Y1ykI7g0gI6KZpSuagkwc7Loaa03t4ejh49ystXL4/r/gsqyrBbrdz0k5/Re9P3WbtqJT6fj4qSufsYhSyvqoi/4BGc3rSAd55/JktrknM+peLUANwV/By1AXcYY+4XkesBjDG3GmN2i8j9wDbAD/zIGLMj6hlnUVDBKJkGaKeDg3RxmC4O0cVhltjXctn4J6YcN0g7i5Mwh5FSKrc4gktnHD52nPHxcWrL4xuZarVYWFRdwaHOHurKS9m5/wAAZXMMmQ8tC9JUVU5ZgovHTmexCBsWJ6f2Sal4GWMOEaFZzBhz67Tb3wS+Od/raTCK03njVRzjBe7nW4DBKhYaKypYaHWxv3sLG0t2Uja8FgA/fgbopKp4aUbLrJRKv1Afoz2HjgBQO0dNTySvWX8q3SNjrKyt4qt/+js+vz/qrNchpcFO3uuaG+O+nlKqQIJRspvQ9vAoZUVOPnDxRhZUlmKzWhh2e/iv+x5mp+u3nB8MRiP04sdLdXHiU/IrpXJTaB6jvUePAdAQ5wSREBiZtiQ4CnrjkiY2Hz4+paksEotF+NTVl1Cs/RqVSkhBBKNkGmeMFrZybvNCmmsmp9ovcRRx3vJFPHVgH2c4juNyLwobkTa/eUSUUrnHZrFgEWF4zE1NWem8g8qr15/KivoaKmLoxB1tkVml1Nx0uH4czhuv4jgv4cXDuoUz29tfvnoZAHvLfw1MTu5YW67BSKlCIyITHbAbK+Y/832Jo4izljTN+zxKqdlpMIrTYZ6l2O5gxYKZozOqS4rZsHghW3tewmPrZoAOLAiVcQ7TVUrlh9CQ/UU1yZlPSCmVehqM4uBjnGO8wNqmBqyWyE/dZacsx+PzcbDqTgZop8JVHPVYpVR+C9UYNSVp2LxSKvW0j1GMQqPRxhnj9Mbow1YbK8s5tbGOF7qfwUGD9i9SqoCFhuw3zHO9MqVU+mhVRhwO8SxOaxGrF84+8/Vlp65g2OOhh+M6Ik2pAlZks+EqsmtzulI5RINRDM4br8KPl6NsYU1THTbr7E/b8tpqFldXAmgwUqqAndJQxzlLm3Xme6VySN4Ho2GJf42jSFrZjZth1jXOPSpERHjFaSsAqCub/2gUpVRuesVpK3j9hjWZLoZSKg7axyhGh3kWu8XGaU2xjS5Zu7CBf77sApbWxL7grFJKKaUyS4NRDJ6xdXPE+xynNdZht80+HX+IiLC8TofoKqWUUrkk75vSkmHYHGKUftYt1MnVlFJKqXymwSgG/WYrVrFyWlNNpouilFJKqRTSYDQHY3z0+19kdUMtriJteVRKKaXymQajOXT6H2acXjYuWZzpoiillFIqxTQYzcJt2jnpv481jY2sX1yX6eIopZRSKsU0GEVhjJ/jvl9SZBXecOZanaBNKaWUKgAajKLo9j/OsDnI6zasoarUkeniKKWUUioNNBhF4DHdtPnvYlV9Pectj75grFJKKaXyS0EFI2P8+MzIHMcYjvtux2ox/MNZp2sTmlJKKVVACmr8eYvvV/SYJ7GJE4fUYjP1FEk1Voqx4sIqxbhNO0NmD9esX0ddua6IrZRSShWSggpGHrqoKnaxuqGBnuFheob30DM6hs/vn3Lc0poaLlqtw/OVUkqpQlNQwcjHKM3lpbz5nLVTto/7fIx6xhkb9zLm9bKwohyLNqEppZRSBaeggpFhlCKbfcZ2u9WK3WWl3JWBQimllFIqaxRU52sfbhwRgpFSSimlFBRaMDJjOGwFVUmmlFJKqTgUTDAyxo8fjwYjpZRSKgeJiFVEXhSR+yLsqxCRP4rISyKyU0Tel+h1CiYY+RkD0KY0pZRSKjd9BNgdZd8NwC5jzBnApcD/iEhRIhfJSDASkW+KyB4R2SYid4lIZdi+z4rIARHZKyKvDNt+tohsD+77rsQ586IvGIycdq0xUkoppXKJiDQDrwF+FOUQA5QFs0Ep0AN4E7lWplLCg8BnjTFeEfkG8Fng0yKyBngrsBZYCDwkIquNMT7gFuA64Bngz8DVwF9ivaBfg5FSSimVdKMGdrjndYpaEdkSdnuTMWbTtGO+A3wKKItyjpuBe4HW4DFvMcb4oxw7q4zUGBljHjDGhJLcM0Bz8PdrgF8bY9zGmMPAAeBcEWkEyo0xTxtjDPBz4Np4rukzoWBUMK2HSimlVC7oMsZsDPuZEopE5LVAhzHm+VnO8UpgK4FKlQ3AzSJSnkhhsiElvJ/Jmp8m4HjYvpbgtqbg79O3RyQi14nIFhHZ4jVDwGRTmqvImrSCK6WUUirlLgJeLyJHgF8DrxCRX0475n3AH0zAAeAwcGoiF0tZMBKRh0RkR4Sfa8KO+TyBNsDbQ5sinMrMsj0iY8ymUPK0SSkQ1pRWpJ2vlVJKqVxhjPmsMabZGLOUQHebh40x75x22DHgcgARaQBOAQ4lcr2Udbgxxlwx234ReQ/wWuDyYPMYBGqCFoUd1kygvbCFyea28O0xm+h8rcP1lVJKqZwnItcDGGNuBb4C3CYi2wlUpnzaGNOVyHkzkhJE5Grg08AlxpiRsF33AneIyLcJtBOuAp41xvhEZFBEzgc2A+8GvhfPNf3BPkYO7XytlFJK5SRjzCPAI8Hfbw3b3gpclYxrZCol3Aw4gAeDo+6fMcZcb4zZKSJ3ArsINLHdEByRBvAh4DbARaBPUswj0mCyxshh0z5GSimllIosI8HIGLNyln1fBb4aYfsWYF2i1/Qzht1qxWrJhv7mSimllMpGBZMSfIzprNdKKaWUmlXBBCO/GdNmNKWUUkrNqmCCkY9RnfVaKaWUUrMqmGDkZ0yDkVJKKaVmpcFIKaWUUiqogILRKEVWDUZKKaWUiq5ggpGOSlNKKaXUXAonGBk3DrsGI6WUUkpFVxDByG/GMfhw6DppSimllJpFYQSjieVANBgppZRSKrqCCEahddJ0VJpSSimlZlNgwUhnvlZKKaVUdAURjPwmEIxcRRqMlFJKKRVdYQSjUI1RkTalKaWUUiq6gghGE01p2vlaKaWUUrMojGBkRgHtfK2UUkqp2RVEMPLjBsChwUgppZRSsyiIYORjDEEosmrna6WUUkpFVxDByM8oRTYbIpLpoiillFIqixVEMPKZMe14rZRSSqk5FUQw8uPW/kVKKaWUmlNBBCMfozoiTSmllFJzKohg5GdMlwNRSiml1JwKJBiN4rDZM10MpZRSSmW5AglGYxqMlFJKKTWngghGPsYo0lFpSimlVM4SEauIvCgi90XYJyLyXRE5ICLbROSsRK+T98HIGIPPeHBoMFJKKaVy2UeA3VH2vQpYFfy5Drgl0YvkfTACAxic2pSmlFJK5SQRaQZeA/woyiHXAD83Ac8AlSLSmMi1CiAY+QFdJ00ppZTKYd8BPkXoQ32mJuB42O2W4La45X1aMBgAHDpcXymllEq6UfGxzdab+AnGqRWRLWFbNhljNoVuiMhrgQ5jzPMicmmUs0Ra88skUpy8D0ahcOnSGiOllFIqG3UZYzbOsv8i4PUi8mrACZSLyC+NMe8MO6YFWBR2uxloTaQwed+UFqoxchVpjZFSSimVa4wxnzXGNBtjlgJvBR6eFooA7gXeHRyddj7Qb4xpS+R6BVCNEqwxKiqAh6qUUkoVCBG5HsAYcyvwZ+DVwAFgBHhfoufN+7RgQp2vdbi+UkopldOMMY8AjwR/vzVsuwFuSMY18r4pLdT3SheRVUoppdRc8j4YTdYYaR8jpZRSSs0u74MRGKwWCzarBiOllFJKzS7vg5HBr/2LlFJKKRWTvA9G4Mehy4EopZRSKgZ5H4wMBqfOeq2UUkqpGOR9MCouteocRkoppZSKSd4nBp/Pp32MlFJqmut/+AtKKqsyXQw1T8N9vdz6T+/KdDHySt7XGPl9PlxFRZkuhlJKZRUNRflB/47Jl/fByOf34yzSztdKKaWUmlv+ByOtMVJKKaVUjPK+840xBqd2vlZKqVm1t7fzhS98ge7ubiwWC294wxt429veBsDAwACf+cxnaGtro7GxkRtvvJHy8nL6+/v51Kc+xa5du3jd617Hpz71qRnn/fjHP05LSwt33nlnxOvu3r2bL33pS7jdbi666CI++clPIiK88MIL/M///A/79+/n61//OpdffnnE+3s8Hr74xS+ye/duKioq+PrXv87ChQsB+O53v8sTTzwBwAc/+EGuuuqqGfe/6aabeOyxx7Db7TQ3N/PFL36RsrIydu7cyVe/+lUg8Dly3XXXcdlll824/6ZNm7jrrruoqgo0ad1www1cdNFFjI+P87WvfY1du3ZhsVj45Cc/ydlnnx3z9QF++tOfcs8992CxWPi3f/s3LrjggojPgUquvK8xAnQeI6WUmoPVauVjH/sYv/vd77jtttv47W9/y6FDhwC47bbbOPfcc7nrrrs499xzue222wAoKiriQx/6EB/96EcjnvPhhx/G5XLNet2vf/3rfP7zn+euu+7i+PHjPP300wAsWLCAL33pS1x99dWz3v+ee+6hrKyMu+++m3e84x1873vfA+CJJ55gz549/OpXv+JnP/sZv/jFLxgeHp5x//POO48777yTX//61yxevJif/vSnAKxYsYJf/OIX3HHHHXzve9/ja1/7Gj6fL2IZ3v72t3PHHXdwxx13cNFFFwFw1113AfCb3/yG73//+/zv//4vfr8/5usfOnSIBx54gDvvvJObb76ZG2+8MeL9VfIVRjDSBWSVUmpWtbW1nHrqqQAUFxezdOlSOjs7AXjkkUd47WtfC8BrX/taHnnkEQBcLhcbNmygKEJ3hZGREW6//XY+8IEPRL1mV1cXw8PDrF+/HhHhNa95zcS5Fy5cyKpVq7BYZv+YevTRRyfKdvnll/Pss89ijOHQoUOcffbZWK1WXC4Xq1at4qmnnppx//PPPx9rcMmo008/nY6ODgCcTufEdo/Hg4jMWo7pDh8+zDnnnANAdXU1ZWVl7N69O+brP/roo1x11VUUFRWxcOFCFi1axI4dO+Iqg0pMQQQjpy4gq5RSMWttbWXv3r2sW7cOgJ6eHmpra4FAgOrt7Z3zHLfeeivvfOc7cTqdUY/p6OigoaFh4nZ9ff1EMIhVR0cHCxYsAAK1XqWlpfT397N69WqefPJJxsbG6OvrY8uWLbS3t896rnvvvZcLL7xw4vaOHTv4x3/8R97ylrfw2c9+diLAfOUrX2HXrl0Tx91555289a1v5ctf/jIDAwMArFq1ikcffRSfz0drayu7d++euP70+0e6fmdn54znJhRUVWoVRFVKkSW+pK+UUoVqZGSET33qU3ziE5+gpKQkoXPs27eP48eP8/GPf5zW1ta47htvzUy0c5x//vns2rWL97///VRWVrJ+/fqJYBPJT37yE6xWK6961asmtq1bt44777yTw4cP88UvfpGLLrqIoqIi/uM//mPimDe96U188IMfBAJh8Dvf+Q5f+MIXuOaaazhy5Ajvete7WLBgwZTrh98/2vWNMREfl0q9wghG1oKoGFNKqXnxer186lOf4lWvehWveMUrJrZXV1fT1dVFbW0tXV1dEx2No9m2bRu7d+/mda97HT6fj56eHq677rqJWiSAl7/85bzpTW+aUovT0dFBXV3drOf+wQ9+MNGh+o477qC+vp6TJ09SX1+Pz+djaGiI8vJyAN7//vfz/ve/H4DPf/7zLF68OOI577vvPh5//HFuueWWiOFj2bJluFwuDhw4wJo1a6bsq66unvj92muvnehvZbVa+fjHPz6x7/3vfz+LFi2K+fr19fUznptQrZ1KrYJIDDrztVJKzc4Yw5e//GWWLVvGO97xjin7LrnkEu677z4g8CF+6aWXznquN73pTdx///388Y9/5Ec/+hGLFy9m06ZNWCyWiU7K119/PbW1tZSUlLB9+3aMMfzpT3/ikksumfXc//zP/zxxDggErFDZ/va3v3HOOecgIvj9fvr7+wHYv38/+/fv5/zzz59xvqeeeoqf/exnfPvb357S7Nfa2jrR2bqtrY2jR49OjHYL19XVNfH7I488wsqVKwEYGxtjdHQUgM2bN2O1Wlm+fHnM17/kkkt44IEH8Hg8tLa2cuzYsYmmTZVaGUkMIvJN4HWABzgIvM8Y0yciS4HdwN7goc8YY64P3uds4DbABfwZ+IiJVNcYgVM7Xyul1Kxeeukl/vznP7Ny5Ure/va3A5NDz9/73vfymc98hnvuuYcFCxbwjW98Y+J+r3vd6xgeHmZ8fJxHHnmEm2++OWIAiOYzn/nMlOH6oT42u3bt4pOf/CQDAwM89thj/PCHP4w45P/aa6/lP/7jP7j22mspLy/n61//OhCo/Qo1cZWUlPBf//VfEZvS/vu//xuPx8MNN9wABJrPPve5z7F161Zuu+02bDYbIsJnPvMZKisrgUAfoTe+8Y2sWbOG7373u+zduxcRYeHChXzuc58DAv2yPvzhD2OxWKivr+fLX/7yxDXD7x/t+suXL+fKK6/kzW9+M1arlc985jNzdkRXySExZovkXlTkKuBhY4xXRL4BYIz5dDAY3WeMmRGLReRZ4CPAMwSC0XeNMX+J4VrmK9dehUtnv1ZKqQmf+M19mS6CSpL/ectrp9z+3wce97f09qdt1FGxZYlZZft0wvffNn7D88aYjUks0rxkJH4aYx4wxniDN58Bmmc7XkQagXJjzNPBWqKfA9fGej1tSlNKKaVULLKhXu79QHjNzzIReVFEHhWRi4PbmoCWsGNagtsiEpHrRGSLiGwRwKKj0pRSSikVg5RVpYjIQ8CCCLs+b4y5J3jM5wEvcHtwXxuw2BjTHexTdLeIrAUiJZuobYDGmE3AJgCb1ZL+tkKllFJK5aSUBSNjzBWz7ReR9wCvBS4PdaI2xrgBd/D350XkILCaQA1ReHNbMxDT5BgSMVMppZRSSs2UkaY0Ebka+DTwemPMSNj2OhGxBn9fDqwCDhlj2oBBETlfApM8vBu4J8ZrJb38SimllMpPmepjdDNQBjwoIltF5Nbg9pcD20TkJeB3wPXGmJ7gvg8BPwIOEBjiP+eINADtXqSUUjMN9829rIfKfvp3TL6MDNcyxqyMsv33wO+j7NsCxD27ldYYKaXUTLf+07syXQSlslI2jEpLKYsGI6WUUipniYhTRJ4VkZdEZKeI/GeEY94hItuCP0+JyBmJXi/vJ/jRmUKVUkqpnOYGXmGMGRIRO/CEiPzFGPNM2DGHgUuMMb0i8ioCI9PPS+RieR+MtClNKaWUyl3BketDwZv24I+ZdsxTYTfnnDh6NnlfnZL3D1AppZQKl4fz1IiIVUS2Ah3Ag8aYzbMc/gFiHKAV8VqZWCstnYqL7L7KYtd4psuRaQasAr5MlyMb6HMRoM/DJH0uAvR5CMj156F3ZNQ66hlP2wKhInI/UDuPUziBsbDbm4ITNUe6ViVwF/AvxpgdEfZfBvwAeJkxpjuRwuR9MBKRLdm0OF2m6PMwSZ+LAH0eJulzEaDPQ4A+D9lNRL4IDBtjvjVt+3oCoelVxph9iZ5fW5qUUkoplbWCkz9XBn93AVcAe6Ydsxj4A/Cu+YQiKIDO10oppZTKaY3Az4IrY1iAO40x94nI9QDGmFuBLwA1wA+Cg668idb6FUIwithOWYD0eZikz0WAPg+T9LkI0OchQJ+HLGKM2QacGWH7rWG/fxD4YDKul/d9jJRSSimlYqV9jJRSSimlgvI2GInI1SKyV0QOiMhnMl2edBORIyKyPbhI75bgtmoReVBE9gf/rcp0OZNNRH4iIh0isiNsW9THLSKfDb5G9orIKzNT6tSI8lx8SUROBF8XW0Xk1WH78vK5EJFFIvJ3EdkdXE7gI8HtBfW6mOV5KMTXRMQlJgrtNaEiy8umtGAHrX3AlUAL8BzwNmPMrowWLI1E5Aiw0RjTFbbtv4EeY8yNwbBYZYz5dKbKmAoi8nICM6T+3BizLrgt4uMWkTXAr4BzgYXAQ8BqY0zOzl8SLspz8SVgKMIw17x9LkSkEWg0xrwgImXA88C1wHspoNfFLM/DP1J4rwkBSsKXmAA+AvwDBfSaUJHla43RucABY8whY4wH+DVwTYbLlA2uAX4W/P1nBN4U84ox5jGgZ9rmaI/7GuDXxhi3MeYwcIDAaycvRHkuosnb58IY02aMeSH4+yCwG2iiwF4XszwP0eTl8wCBJSaMMZGWmCio14SKLF+DURNwPOx2C7O/AeQjAzwgIs+LyHXBbQ3GmDYIvEkC9RkrXXpFe9yF+jr5sARWoP5JWFNBQTwXIrKUwOiWzRTw62La8wAF+JqQyEtMFOxrQk3K12AUaZ2Y/GsznN1FxpizgFcBNwSbVdRUhfg6uQVYAWwA2oD/CW7P++dCREqB3wMfNcYMzHZohG1581xEeB4K8jVhjPEZYzYQWGz0XBFZN8vhef1cqKnyNRi1AIvCbjcDrRkqS0YYY1qD/3YQmCL9XKA92M8g1N+gI3MlTKtoj7vgXifGmPbgB4If+D8mmwPy+rkI9iP5PXC7MeYPwc0F97qI9DwU6msixBjTBzwCXE0BvibUTPkajJ4DVonIMhEpAt4K3JvhMqWNiJQEO1ciIiXAVcAOAs/Be4KHvQe4JzMlTLtoj/te4K0i4hCRZcAq4NkMlC9tQm/6QW8g8LqAPH4ugh1tfwzsNsZ8O2xXQb0uoj0PBfqaiLbEREG9JlRkeTnztTHGKyIfBv4KWIGfGGN2ZrhY6dQA3BV4H8QG3GGMuV9EngPuFJEPAMeAN2ewjCkhIr8CLgVqRaQF+CJwIxEetzFmp4jcCewCvMAN+TTKJMpzcamIbCDQDHAE+CfI++fiIuBdwPZgnxKAz1F4r4toz8PbCvA1EW2JiacprNeEiiAvh+srpZRSSiUiX5vSlFJKKaXipsFIKaWUUipIg5FSSimlVJAGI6WUUkqpIA1GSimllFJBGoyUSpCI+IKrke8MrtL9cRGxBPdtFJHvznLfpSLy9vSVNjYickREalN07tNlcgX3HhE5HPz9oVRcL86yzfr3EJFGEbkv+Pt7ReTmeVzrobBlN5RSWUaDkVKJGzXGbDDGrAWuBF5NYK4gjDFbjDH/Ost9lwJZF4xSyRizPfh8bSAwYd6/BW9fETpGRDI1t9pSZv97fJzArNDJ8Avgn5N0LqVUkmkwUioJgkuvXEdgMU4RkUvDahguCaspeTE4K/mNwMXBbR8L1lg8LiIvBH8uDN73UhF5RER+JyJ7ROT24AzGiMg5IvJUsLbqWREpCy6M+U0ReU4Ci4L+U6TyisjdElhgeKdMLjI8/ZiPi8iO4M9Hg9uWishuEfm/4H0fCM4cHCrPNhF5OliGHZHOG+E6j4jI10TkUeAjIvI6EdkcfK4eEpGG4HFfksAip4+IyCER+dfg9hIR+VPwedghIm8Jbj9bRB4NPs6/yuRSDyuD530p+FyvmP73iFDMNwL3Ryj7a4KPt1ZEbhORW0Tk78HyXRIs724RuS3sbvcCb4vluVFKZYAxRn/0R38S+AGGImzrJTDz+KXAfcFtfySwqC9AKYHZyCf2B7cXA87g76uALcHfLwX6CazNZAGeBl4GFAGHgHOCx5UHz3sd8O/BbQ5gC7AsQjmrg/+6CCwBURO8fQSoBc4GtgMlwTLvJLAa+1ICM/9uCB5/J/DO4O87gAuDv98I7JjlubsNeFPw90eAH4Ttq2Jy8tkPAv8T/P1LwFPBx1ULdAN2AqHl/8LuXxHc/hRQF9z2FgIz4ENgRfk3BH93Bp/7KX+PaWVdBjwfdvu9wM0Els94HKgKe0y/JrDg6DXAAHB68O/2fOg5Cx67P/Sc64/+6E92/eTlkiBKZVCkVbifBL4tIrcDfzDGtAQrfcLZgZslsDSDD1gdtu9ZY0wLgASWclhKICy1GWOeAzDB1eJF5CpgvYi8KXjfCgJB6/C06/2riLwh+Pui4DHdYftfBtxljBkOnvcPwMUEajsOG2O2Bo97HlgqgXWnyowxTwW33wG8NsJzEc1vwn5vBn4TrOEpmlb2Pxlj3IBbRDoIhNDtwLdE5BsEws3jElgpfR3wYPC5tgJtwdq6JmPMXQDGmLHg45utbI1A57RtlwEbgatCz33QH40xRkS2A+3GmO3B8+8k8HfbGjyuA1jI1OdcKZUFtClNqSQRkeUEQk1H+HZjzI0Eaj5cwDMicmqEu38MaAfOIPCBWxS2zx32u49AzZAQWNtqRjGAfzHBvjzGmGXGmAemlfNSAotmXmCMOQN4kUDNyfTzRBOtPPMxHPb794CbjTGnE1i3K7xsM65tjNnHZA3X10XkC8Hy7Ax7Hk43xlyVYDlHmfn8HALKmBpgw8vnn1ZWP1PXpnQGz6uUyjIajJRKAhGpA24l8IFupu1bYQIdj79BoGnrVGCQwAdrSAWBGiA/gYU+rXNccg+wUETOCV6jLNhx+a/Ah0TEHty+WkRKpt23Aug1xowEQ9r5Ec7/GHCtiBQH7x9qNorIGNMLDIpI6FxvnaP8s6kATgR/f89sBwKIyEJgxBjzS+BbwFnAXqBORC4IHmMXkbXB2p0WEbk2uN0hIsXM/HuE20egtifcUeAfgJ+LyNo4HltolfsFBJotlVJZRoORUolzBTvr7gQeAh4A/jPCcR8Ndgp+iUAtwV+AbYA32AH4Y8APgPeIyDMEaiGGI5xngjHGQ6DfzPeC532QQC3EjwisAP5CsPPzD2FGk/n9gE1EtgFfAZ6JcP4XCPSZeZZAn5wfGWNenOP5+ACwSQIrlAuB5r5EfAn4rYg8DnTFcPzpwLPBZsbPA/8VfH7eBHwj+PxsBS4MHv8uAk2J2wj0Q1rAzL/HhGBz4kERWTlt+17gHcGyrojj8Z0NPGOM8cZxH6VUmsi0L7dKKZUQESk1xgwFf/8M0GiM+UiGi5UUwf5YZxtj/j0J57oJuNcY87f5l0wplWza+VoplSyvEZHPEnhfOUpg9FZeMMbcJSI1STrdDg1FSmUvrTFSSimllArSPkZKKaWUUkEajJRSSimlgjQYKaWUUkoFaTBSSimllArSYKSUUkopFaTBSCmllFIq6P8Dwzpn3l9flI8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mycmap = 'plasma'\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "# Draw the submerged earth\n",
    "ax.fill_between(transect['distance'], -1 * transect['depth'].max(),\n",
    "                -1 * transect['depth'], color='tab:brown', zorder=2)\n",
    "# Give a nice line for the bottom\n",
    "ax.plot(transect['distance'], -1 * transect['depth'], 'k', zorder=3)\n",
    "ax.set(xbound=(0, transect['distance'].max()), ybound=(-1 * transect['depth'].max(), 10),\n",
    "       xlabel=\"Distance along Transect (km)\", ylabel=\"Depth (m)\")\n",
    "node_indices = np.searchsorted(data['node'][:], transect.index.to_numpy())\n",
    "dates = pd.Timestamp(start_date) + pd.to_timedelta(data['time'][:], 'd')\n",
    "date_text = ax.text(transect['distance'].max() - 0.02 * (transect['distance'].max() -\n",
    "                                                        transect['distance'].min()),\n",
    "                    -1 * transect['depth'].max() + 10, dates[0].strftime(\"%Y-%m-%d %X\"),\n",
    "                    backgroundcolor='w', alpha=0.8, ha='right', va='baseline')\n",
    "\n",
    "all_vals = data[variable][:,:, node_indices]\n",
    "norm = colors.Normalize(all_vals.min(), all_vals.max())\n",
    "\n",
    "sfc = None\n",
    "tcf = None\n",
    "colorbar = None\n",
    "\n",
    "def animate(t):\n",
    "    global sfc, tcf, cbar_ax, fig, colorbar\n",
    "    zetas = data['zeta'][t, node_indices]\n",
    "    if sfc is None:\n",
    "        sfc, = ax.plot(transect['distance'], zetas, 'b')\n",
    "    else:\n",
    "        sfc.set_data(transect['distance'], zetas)\n",
    "    artists = [sfc]\n",
    "    # Take the difference of the water surface elevations (zetas) and the depths\n",
    "    # (transect['depth']) at each point to get the thickness of the water column.\n",
    "    # Multiply this by the sigma layer depth fractions to get the z value that\n",
    "    # corresponds to each layer at each node, with some linear algebra needed to\n",
    "    # make the multiplication and final array dimensions work out how we want\n",
    "    zs = (transect['depth'].to_numpy() - zetas)[np.newaxis].T * data['siglay'][:]\n",
    "    # Offset z by the tide height (the adjusted zero point)\n",
    "    zs += zetas[np.newaxis].T\n",
    "\n",
    "    # Repeat the distances to get matching x coordinates for the z's\n",
    "    xs = np.repeat(transect['distance'].to_numpy(), zs.shape[1])\n",
    "    xs.shape # (2360,)\n",
    "\n",
    "    vals = data[variable][t, :, node_indices].T.flatten()\n",
    "    if tcf is not None:\n",
    "        for c in tcf.collections:\n",
    "            c.remove()\n",
    "    tcf = ax.tricontourf(xs, zs.data.flatten(), vals, cmap=mycmap)\n",
    "    artists += tcf.collections\n",
    "    if colorbar is None:\n",
    "        colorbar = fig.colorbar(tcf)\n",
    "        artists += colorbar.lines\n",
    "    ax.set(xbound=(0, transect['distance'].max()), ybound=(-1 * transect['depth'].max(), 10))\n",
    "    date_text.set_text(dates[t].strftime(\"%Y-%m-%d %X\"))\n",
    "    return artists\n",
    "\n",
    "anim = animation.FuncAnimation(fig, animate, data['time'].size, blit=True)\n",
    "anim.save(out_file(\"{0}.mp4\".format(variable)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ssm-analysis]",
   "language": "python",
   "name": "conda-env-ssm-analysis-py"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}