6_calculate_assemblage_apparent_density.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "52753700-bcad-43bc-a1cd-9e8904c04d08",
   "metadata": {},
   "source": [
    "# Conversion of the volume to mass concentrations\n",
    "This notebook executes the following steps:\n",
    "* Calculate a theoretical apparent densities for the LISST-200X size classes using the formula from Fall et al. 2021. \n",
    "* Average the densities of the different size classes to match the grain size populations\n",
    "* Convert the predicted TVC grids to Suspended Particulate Matter Concentration (SPMC) grids.\n",
    "* Saves the SPMC grids to netCDF files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "73904840-b128-40a9-b6b1-caca1b109431",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "rc(\"pdf\", fonttype=42)\n",
    "from timbers_code.grainsize import get_assemblage_apparent_density\n",
    "from glob import glob\n",
    "import os\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3147a937",
   "metadata": {},
   "source": [
    "## Set the parameters\n",
    "* lisst_size_classes: path to the file that defines the size classes. From LISST-200X Manual\n",
    "* d_p: the primary particle size in µm, below which rho_a = 2500 kg/m3 (Fettweis et al., 2008)\n",
    "* d_p_li: lower interval of primary particle size in µm for Kwinte study area (Fettweis et al., 2008)\n",
    "* d_p_ui: upper interval of primary particle size in µm for Kwinte study area (Fettweis et al., 2008) \n",
    "* F: Fractal Dimension (Fettweis et al., 2021) \n",
    "* F_li: Fractal Dimension, lower interval for Kwinte study area (Fettweis personal communication) \n",
    "* F_ui: Fractal Dimension, upper interval for Kwinte study area (Fettweis personal communication) \n",
    "* rho_p: primary particle density in kg/m3 (Fettweis et al., 2008)\n",
    "* raster_path: path where output netCDF files will be stored\n",
    "* raster_files: netCDF files of predicted TVC data\n",
    "* bin_boundaries: list of grain size boundaries used to calculate TVC for different grain size populations\n",
    "* variable_schema: mapping of which which TVC variables (in each of the netCDF files for the different grain size populations/assemblages) should be used to calculate the SPMC variables\n",
    "* filepath_apparent_densities_plot: file to save the plot of the apparent densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d01c5fdd-bd80-4d7d-a1f5-7c85ee72d09f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# input files\n",
    "lisst_size_classes = 'sensor_data/lisst_size_classes.csv'\n",
    "# input parameters\n",
    "d_p = 2 #µm, the primary particle size, below which rho_a = 2500 kg/m3 (Fettweis et al., 2008)\n",
    "d_p_li = 1 #µm,lower interval of primary particle size for Kwinte study area (Fettweis et al., 2008)\n",
    "d_p_ui = 3 #µm, upper interval of primary particle size for Kwinte study area (Fettweis et al., 2008) \n",
    "F = 2 # Fractal Dimension (Fettweis et al., 2021) \n",
    "F_li = 1.9 # Fractal Dimension, lower interval for Kwinte study area (Fettweis personal communication) \n",
    "F_ui = 2.1 # Fractal Dimension, upper interval for Kwinte study area (Fettweis personal communication) \n",
    "# F_li = 2.04 # Fractal Dimension, lower interval for Kwinte study area (Fettweis et al., 2008) \n",
    "# F_ui = 2.12 # Fractal Dimension, upper interval for Kwinte study area (Fettweis et al., 2008) \n",
    "rho_p = 2500 #kg/m3 (Fettweis et al., 2008)\n",
    "# get predicted raster files\n",
    "raster_path = 'G:/Algorithm2/predictions_final'\n",
    "raster_files = glob(f'{raster_path}/*/*TVC*.nc') + glob(f'{raster_path}/*/*Total_Volume_Concentration*.nc')\n",
    "## For size classes based on Fetweiss et al 2014\n",
    "# bin_boundaries = np.array([15.,200.]) #assemblage boundaries\n",
    "# variable_schema = pd.DataFrame({'SPMC_mean': ['TVC_(0.0,_15.0]_mean','TVC_(15.0,_200.0]_mean','TVC_(200.0,_459.0]_mean','Total_Volume_Concentration_[ppm]_mean'],\n",
    "#                                  'SPMC_ci_lower':['TVC_(0.0,_15.0]_obs_ci_lower','TVC_(15.0,_200.0]_obs_ci_lower','TVC_(200.0,_459.0]_obs_ci_lower','Total_Volume_Concentration_[ppm]_obs_ci_lower'],\n",
    "#                                  'SPMC_ci_upper':['TVC_(0.0,_15.0]_obs_ci_upper','TVC_(15.0,_200.0]_obs_ci_upper','TVC_(200.0,_459.0]_obs_ci_upper','Total_Volume_Concentration_[ppm]_obs_ci_upper']},\n",
    "#                                 index = pd.Series(['(0.0, 15.0]', '(15.0, 200.0]', '(200.0, 459.0]','(0.0, 459.0]']).rename('Assemblage'))\n",
    "## For size classes based on Lee et al 2012\n",
    "bin_boundaries = np.array([3.,20.,200.]) #assemblage boundaries\n",
    "variable_schema = pd.DataFrame({'SPMC_mean': ['TVC_(0.0,_3.0]_mean','TVC_(3.0,_20.0]_mean','TVC_(20.0,_200.0]_mean','TVC_(200.0,_420.0]_mean','Total_Volume_Concentration_[ppm]_mean'],\n",
    "                                 'SPMC_ci_lower':['TVC_(0.0,_3.0]_obs_ci_lower','TVC_(3.0,_20.0]_obs_ci_lower','TVC_(20.0,_200.0]_obs_ci_lower','TVC_(200.0,_420.0]_obs_ci_lower','Total_Volume_Concentration_[ppm]_obs_ci_lower'],\n",
    "                                 'SPMC_ci_upper':['TVC_(0.0,_3.0]_obs_ci_upper','TVC_(3.0,_20.0]_obs_ci_upper','TVC_(20.0,_200.0]_obs_ci_upper','TVC_(200.0,_420.0]_obs_ci_upper','Total_Volume_Concentration_[ppm]_obs_ci_upper']},\n",
    "                                index = pd.Series(['(0.0, 3.0]','(3.0, 20.0]','(20.0, 200.0]', '(200.0, 459.0]','(0.0, 459.0]']).rename('Assemblage'))\n",
    "# output files\n",
    "filepath_apparent_densities_plot = 'publication_output/apparent_densities.pdf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8b1352-7210-4019-ae0c-1adabe732ff9",
   "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>SPMC_mean</th>\n",
       "      <th>SPMC_ci_lower</th>\n",
       "      <th>SPMC_ci_upper</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Assemblage</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(0.0, 3.0]</th>\n",
       "      <td>TVC_(0.0,_3.0]_mean</td>\n",
       "      <td>TVC_(0.0,_3.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(0.0,_3.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(3.0, 20.0]</th>\n",
       "      <td>TVC_(3.0,_20.0]_mean</td>\n",
       "      <td>TVC_(3.0,_20.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(3.0,_20.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(20.0, 200.0]</th>\n",
       "      <td>TVC_(20.0,_200.0]_mean</td>\n",
       "      <td>TVC_(20.0,_200.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(20.0,_200.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(200.0, 459.0]</th>\n",
       "      <td>TVC_(200.0,_420.0]_mean</td>\n",
       "      <td>TVC_(200.0,_420.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(200.0,_420.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(0.0, 459.0]</th>\n",
       "      <td>Total_Volume_Concentration_[ppm]_mean</td>\n",
       "      <td>Total_Volume_Concentration_[ppm]_obs_ci_lower</td>\n",
       "      <td>Total_Volume_Concentration_[ppm]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            SPMC_mean  \\\n",
       "Assemblage                                              \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_mean   \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_mean   \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_mean   \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_mean   \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_mean   \n",
       "\n",
       "                                                SPMC_ci_lower  \\\n",
       "Assemblage                                                      \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_obs_ci_lower   \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_obs_ci_lower   \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_obs_ci_lower   \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_obs_ci_lower   \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_obs_ci_lower   \n",
       "\n",
       "                                                SPMC_ci_upper  \n",
       "Assemblage                                                     \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_obs_ci_upper  \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_obs_ci_upper  \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_obs_ci_upper  \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_obs_ci_upper  \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_obs_ci_upper  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variable_schema"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3e2d2ca-c06a-4ff5-90f0-16299a745e06",
   "metadata": {},
   "source": [
    "## Read in size classes, exported from the LISST manual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a3c12ace-5fd5-4f43-870a-7845b0c4fe55",
   "metadata": {},
   "outputs": [],
   "source": [
    "size_classes = pd.read_csv(lisst_size_classes)\n",
    "size_classes = size_classes.set_index('Size_bin')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1796f2dc",
   "metadata": {},
   "source": [
    "## Calculate the apparent densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0cbabad9-1748-4484-8fc4-cc04eab0d5d0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAJQCAYAAAAOv1mVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAB0QklEQVR4nO3deXxU5aH/8c8zk33fEwj7vgRkCYKAC6ICakVtXVotam29VWytdrV6q1btra1tf/W61Vat29W61H0XtYiiEnbCvgQIewLZ95nn98cMIUCABJKcM8n3/XrNa86cOTN8wzXle885z/MYay0iIiIiEjo8TgcQERERkdZRgRMREREJMSpwIiIiIiFGBU5EREQkxKjAiYiIiISYMKcDdLS0tDTbp08fp2OIiIiIHNPChQuLrLXph+7vcgWuT58+5OXlOR1DRERE5JiMMZub269LqCIiIiIhRgVOREREJMSowImIiIiEmC53D5yIiIi0r/r6egoLC6mpqXE6SsiIioqiR48ehIeHt+h4FTgRERFpU4WFhcTHx9OnTx+MMU7HcT1rLcXFxRQWFtK3b98WfUaXUEVERKRN1dTUkJqaqvLWQsYYUlNTW3XGUgVORERE2pzKW+u09u9LBU5EREQkxKjAiYiIiIQYFTgRERHp1AoKCsjJyXE6RptSgRMREZFOzVqL3+93Okab0jQiIiIi0m7uejOfldvL2vQ7h3VP4I5vDD/qMQUFBcyYMYMpU6bw/PPPk5SUxA9+8AO++OILsrOzef3114mOjmbJkiX88Ic/pKqqiv79+/PEE0+QnJzc7Hf+/e9/57HHHqOuro4BAwbwzDPPEBMT06Y/W0vpDJyIiIh0SmvWrGHWrFksXryYrVu3Mnv2bPLz80lKSuKVV14BYNasWdx3330sW7aMESNGcNdddx3x+y6++GIWLFjA0qVLGTp0KI8//nhH/SiH0Rk4ERERaTfHOlPWnnr37s2ECRMoKCigb9++jBo1CoCxY8dSUFBAaWkpJSUlnH766QBcddVVXHLJJUf8vhUrVnD77bdTUlJCRUUF06ZN64gfo1k6AyciIiKdUmxsbON2ZGRk47bX66WhoaHV33f11Vfz4IMPsnz5cu644w5HlwpTgRMREZEuKTExkeTkZD777DMAnnnmmcazcc0pLy+nW7du1NfX89xzz3VUzGbpEqqIiIh0WU899VTjIIZ+/frx5JNPHvHYu+++m/Hjx9O7d29GjBhBeXl5ByY9mLHWOvaHOyE3N9fm5eU5HUNERKTTWrVqFUOHDnU6Rshp7u/NGLPQWpt76LG6hCoiIiISYnQJVURERKSJ2bNn8/nnnx+076abbuKaa65xKNHhVOBEREREmnjooYecjnBMuoQqIiIiEmJ0Bq6N/c9X/0NpXanTMdrUuMxxfHPQN52OISIiIkEqcG1s1d5VFFUXOR2jzZTUljC3cC4XD7wYY4zTcURERAQVuDb39IynnY7Qpl5e+zJ3zb+LreVb6ZXQy+k4IiIigu6Bk2PIScsBYEXRCoeTiIiIHJ+CggJycnI69M+88847uf/++wH4zW9+w0cffdSm368zcHJU/ZP6E+mNZEXxCs7td67TcURERFrNWovf73fsz//tb3/b5t+pAidHFe4JZ0jKEPKL8p2OIiIioejdX8HO5W37nVkjYMbvj3pIQUEBM2bMYMqUKTz//PMkJSXxgx/8gC+++ILs7Gxef/11oqOjWbJkSeNSWv379+eJJ54gOTm52e/8+9//zmOPPUZdXR0DBgzgmWeeISYm5phxr776as4//3y+9a1vHdeP2xxdQpVjGp46nFV7V9Hgb3A6ioiISIutWbOGWbNmsXjxYrZu3crs2bPJz88nKSmJV155BYBZs2Zx3333sWzZMkaMGMFdd911xO+7+OKLWbBgAUuXLmXo0KE8/vjjHfWjHEZn4OSYctJy+L/V/8fG0o0MSh7kdBwREQklxzhT1p569+7NhAkTKCgooG/fvowaNQqAsWPHUlBQQGlpKSUlJZx++ukAXHXVVVxyySVH/L4VK1Zw++23U1JSQkVFBdOmTeuIH6NZOgMnxzQ8bTiALqOKiEhIiY2NbdyOjIxs3PZ6vTQ0tP6q0tVXX82DDz7I8uXLueOOO6ipqWmTnMdDBU6OqU9CH+LC48gvVoETEZHOIzExkeTkZD777DMAnnnmmcazcc0pLy+nW7du1NfX89xzz3VUzGbpEqock8d4GJY6TFOJiIhIp/PUU081DmLo168fTz755BGPvfvuuxk/fjy9e/dmxIgRlJeXd2DSgxlrrWN/uBNyc3NtXl6e0zFCzp8X/plnVj7DV9/5ighvhNNxRETExVatWsXQoUOdjhFymvt7M8YstNbmHnqsLqFKi+Sk5tDgb2DtvrVORxEREenydAlVWmT/QIYVRSsaV2cQERHpjGbPns3nn39+0L6bbrqJa6655rBj7733Xl566aWD9l1yySXcdttt7ZpRBU5apHtsd5Ijk3UfnIiIdHoPPfRQi4+97bbb2r2sNUeXUKVFjDEMTxuukagiIiIuoAInLZaTlsPG0o1U1Vc5HUVERKRLU4GTFstJzcFv/azau8rpKCIiIl2aCpy0WNOBDCIiIuIcFThpsbToNLJis7SkloiIhJSCggJyctwxg8L27dv51re+dcLfowInrTI8dTgrinUGTkREQoe1Fr/f73QMALp3787LL798wt+jaUSkVXLScpizZQ6ltaUkRiY6HUdERFzuvq/vY/Xe1W36nUNShvDLk3951GMKCgqYMWMGU6ZM4fnnnycpKYkf/OAHfPHFF2RnZ/P6668THR3NkiVLGpfS6t+/P0888QTJycnNfuff//53HnvsMerq6hgwYADPPPMMMTExzR579dVXc/755zeebYuLi6OiooKCggLOP/98Vqw4sZMhOgMnrTI8NXAfnC6jioiI261Zs4ZZs2axePFitm7dyuzZs8nPzycpKYlXXnkFgFmzZnHfffexbNkyRowYwV133XXE77v44otZsGABS5cuZejQoTz++OMd9aMcRmfgpFX2D2TIL85nYvZEh9OIiIjbHetMWXvq3bs3EyZMoKCggL59+zJq1CgAxo4dS0FBAaWlpZSUlHD66acDcNVVV3HJJZcc8ftWrFjB7bffTklJCRUVFUybNq0jfoxm6QyctEpCRAK9E3prJKqIiLhebGxs43ZkZGTjttfrpaGhodXfd/XVV/Pggw+yfPly7rjjDmpqao54bFhYWON9d9Za6urqWv3nHY0KnLSaBjKIiEhnkJiYSHJyMp999hkAzzzzTOPZuOaUl5fTrVs36uvree6554763X369GHhwoUAvP7669TX17ddcHQJVY5DTloO72x6hz1Ve0iPSXc6joiIyHF76qmnGgcx9OvXjyeffPKIx959992MHz+e3r17M2LECMrLy4947A9+8ANmzpzJySefzNSpUw86G9gWjLW2Tb/Q7XJzc21eXp7TMULa4t2LmfXuLB6Y8gBTek1xOo6IiLjMqlWrGDp0qNMxQk5zf2/GmIXW2txDj9UlVGm1wcmD8RiPLqOKiIg4RJdQpdViwmPon9RfU4mIiEinNHv2bD7//POD9t10001cc801hx1777338tJLLx2075JLLuG2225r14wqcHJcclJz+GTrJ1hrMcY4HUdERFwmlP99eOihh1p87G233dYmZa21t7TpEqocl5y0HEpqS9hWsc3pKCIi4jJRUVEUFxe3upR0VdZaiouLiYqKavFndAZOjsv+CX1XFK+gR3wPh9OIiIib9OjRg8LCQvbs2eN0lJARFRVFjx4t//dUBU6Oy6CkQYR7wskvymd6n+lOxxERERcJDw+nb9++Tsfo1HQJVY5LuDecISlDtCKDiIiIA1Tg5LgNSx3GyuKV+Pw+p6OIiIh0KSpwctxy0nKoaqiioKzA6SgiIiJdigqcHLec1BwA8os1H5yIiEhH6vACZ4yJMsZ8bYxZaozJN8bcFdyfYoz50BizLvic3OQztxpj1htj1hhjpjXZP9YYszz43gMmVCecCVF9E/sSHRat++BEREQ6mBNn4GqBM621JwGjgOnGmAnAr4A51tqBwJzga4wxw4DLgeHAdOBhY4w3+F2PANcBA4MPDYfsQF6Pl2Gpw7Qig4iISAfr8AJnAyqCL8ODDwvMBJ4K7n8KuDC4PRN4wVpba63dBKwHTjbGdAMSrLXzbWCmwKebfEY6SE5qDqv3rqbeV+90FBERkS7DkXvgjDFeY8wSYDfwobX2KyDTWrsDIPicETw8G9ja5OOFwX3Zwe1D9zf3511njMkzxuRpUsG2lZOWQ52/jnUl65yOIiIi0mU4UuCstT5r7SigB4GzaTlHOby5+9rsUfY39+c9Zq3NtdbmpqentzqvHNnw1OCKDLoPTkREpMM4OgrVWlsCfErg3rVdwcuiBJ93Bw8rBHo2+VgPYHtwf49m9ksH6hHfg8TIRI1EFRER6UBOjEJNN8YkBbejgbOA1cAbwFXBw64CXg9uvwFcboyJNMb0JTBY4evgZdZyY8yE4OjTWU0+Ix3EGMPw1OEayCAiItKBnDgD1w34xBizDFhA4B64t4DfA2cbY9YBZwdfY63NB14EVgLvAbOttfun/r8e+AeBgQ0bgHc78geRgOGpw1lfsp7qhmqno4iIiHQJHb6YvbV2GTC6mf3FwNQjfOZe4N5m9ucBR7t/TjpATloOPutjzd41jMoY5XQcERGRTk8rMcgJy0kLdGgNZBAREekYKnBywjJiMsiIzmBFsQqciIhIR1CBkzYxPE0DGURERDqKCpy0ieGpwykoK6CsrszpKCIiIp2eCpy0if33wa0qXuVwEhERkc5PBU7ahFZkEBER6TgqcNImkqKS6BHXQysyiIiIdAAVOGkzOWk5OgMnIiLSAVTgpM3kpOWwo3IHxdXFTkcRERHp1FTgpM3svw9Ol1FFRETalwqctJmhqUMxGF1GFRERaWcqcNJmYsNj6ZfYT2fgRERE2pkKnLSp4WnDWVG0Amut01FEREQ6LRU4aVM5aTnsrdnLzsqdTkcRERHptFTgpE3lpAZWZNDC9iIiIu1HBU7a1OCUwYR5wjSQQUREpB2pwEmbivBGMCh5EPlFGsggIiLSXlTgpM3lpOaQX5yP3/qdjiIiItIpqcBJmxueNpyK+gq2lG1xOoqIiEinpAInbW7/igwayCAiItI+VOCkzfVP6k+UN0r3wYmIiLQTFThpc2GeMIamDtVIVBERkXaiAiftIicth1V7V1Hvq3c6ioiISKejAiftYkzGGGp9tVoXVUREpB2owEm7GJM5BoCFuxY6nERERKTzUYGTdpESlULfxL4s2r3I6SgiIiKdjgqctJsxGWNYvGsxPr/P6SgiIiKdigqctJuxmWMpry9nfcl6p6OIiIh0Kipw0m7GZo4FIG9XnsNJREREOhcVOGk33eO60y22G4t26T44ERGRtqQCJ+1qTOYYFu1ehLXW6SgiIiKdhgqctKsxGWMoqi5iS7kWthcREWkrKnDSrnIzcwF0GVVERKQNqcBJu+qb2JfkyGQNZBAREWlDKnDSrowxgfvgdAZORESkzajASbsbkzGGwopCdlXucjqKiIhIp6ACJ+1u/3xwWlZLRESkbajASbsbnDKYmLAYLWwvIiLSRlTgpN2FecIYlTFKZ+BERETaiAqcdIixmWNZt28dpbWlTkcREREJeSpw0iHGZIwBYPHuxQ4nERERCX0qcNIhRqSPINwTrvvgRERE2oAKnHSISG8kI9JGaD44ERGRNqACJx1mTOYYVhavpKq+yukoIiIiIU0FTjrM2MyxNNgGlhUtczqKiIhISFOBkw4zKn0UHuPRZVQREZETpAInHSYuIo7ByYM1kEFEROQEqcBJhxqbOZZle5ZR76t3OoqIiEjIUoGTDjUmcww1vhpW7l3pdBQREZGQpQInHWr/hL66jCoiInL8VOCkQ6VGp9InoY8GMoiIiJwAFTjpcGMzx7Jo9yL81u90FBERkZCkAicdbmzmWMrrylm3b53TUUREREKSCpx0uDGZgfvgFu3WZVQREZHjoQInHa57bHeyYrM0kEFEROQ4qcBJhzPGMCZjDIt2LcJa63QcERGRkKMCJ44YmzmWPdV7KCwvdDqKiIhIyFGBE0eMzRwLQN6uPIeTiIiIhB4VOHFEv8R+JEUmaSCDiIjIcVCBE0fsvw9OAxlERERaTwVOHDMmcwxby7eyp2qP01FERERCigqcOGb/fXALd+ssnIiISGuowIljhqQMITosmoU7VeBERERaQwVOHBPmCWNU+igNZBAREWklFThx1NjMsazbt47S2lKno4iIiIQMFThx1JjMMVgsS3YvcTqKiIhIyFCBE0eNSBtBmCdMAxlERERaocMLnDGmpzHmE2PMKmNMvjHmpuD+O40x24wxS4KPc5t85lZjzHpjzBpjzLQm+8caY5YH33vAGGM6+ueRExMVFsWItBGaD05ERKQVnDgD1wD81Fo7FJgAzDbGDAu+9xdr7ajg4x2A4HuXA8OB6cDDxhhv8PhHgOuAgcHH9A78OaSNjMkYw8qilVQ3VDsdRUREJCSEdfQfaK3dAewIbpcbY1YB2Uf5yEzgBWttLbDJGLMeONkYUwAkWGvnAxhjngYuBN5tx/jHdP2zCymurHMywhFdOCqb74zv5XSMw4zNHMvjKx5n+Z7lnNztZKfjiIiIuJ6j98AZY/oAo4GvgrtuNMYsM8Y8YYxJDu7LBrY2+VhhcF92cPvQ/c39OdcZY/KMMXl79rTvrP8eY/AYXPfYtq+aP76/mnqfv11//uMxKmMUBqPLqCIiIi3U4Wfg9jPGxAGvAD+x1pYZYx4B7gZs8PlPwPeA5u5rs0fZf/hOax8DHgPIzc1t9pi28tAVY9rz64/bnFW7uPapPP6zZg9nDct0Os5B4iPiGZwyWAMZREREWsiRM3DGmHAC5e05a+2/Aay1u6y1PmutH/g7sP9aWiHQs8nHewDbg/t7NLNfmnHaoHRSYiN4dfE2p6M0a2zmWJbtWUa9v97pKCIiIq7nxChUAzwOrLLW/rnJ/m5NDrsIWBHcfgO43BgTaYzpS2CwwtfBe+nKjTETgt85C3i9Q36IEBTu9fCNkd34cNUuymrcV5LGZIyhuqGaVcWrnI4iIiLiek6cgZsEfBc485ApQ/4QnBJkGTAFuBnAWpsPvAisBN4DZltrfcHvuh74B7Ae2IDDAxjc7qIxPahr8PPu8h1ORznMmMzApedFu7SsloiIyLE4MQp1Hs3fv/bOUT5zL3BvM/vzgJy2S9e5ndQjkb5psby6eBuXjXPXaNS06DT6JPRh4a6FXJ1ztdNxREREXE0rMXQhxhguGp3Nlxv3sq3EfXOujckcw6Ldi/Bb942UFRERcRMVuC7mwlGBmVZec+FghrGZYymrK2N9yXqno4iIiLiaClwX0ys1htzeyby6eBvWtuuMKq02JkP3wYmIiLSEClwXdNGYbNbvriB/e5nTUQ6SHZdNZkymCpyIiMgxqMB1QeeN6EaE1+O6OeGMMYzJHMPCXQtdd3ZQRETETVTguqCkmAimDEnn9SXbaXDZ0lpjM8ayu3o3heWFxz5YRESki1KB66IuGt2Doopa5q0vcjrKQcZljQPgy51fOpxERETEvVTguqgpQ9JJjA533WjUvol9yYrN4ottXzgdRURExLVU4LqoyDAv543sxvv5u6isbXA6TiNjDJO6T+LLHV9qXVQREZEjUIHrwi4enU11vY/383c6HeUgk7InUVFfwfI9y52OIiIi4koqcF3Y2N7J9EyJdt1o1PHdxuM1Xj7f/rnTUURERFxJBa4LM8Zw0ahsPl9fxK6yGqfjNEqISGBk+kjdByciInIEKnBd3IWjs/FbeGPJdqejHGRi94nkF+ezr2af01FERERcRwWui+uXHsdJPZP4t8suo07qPgmLZf72+U5HERERcR0VOOHi0dms2lHG6p3uWVprWOowkiKTdB+ciIhIM1TghPNHdiPMY1w1mMHr8XJKt1P4YvsXWlZLRETkECpwQmpcJKcPSuf1xdvx+91TliZmT6Souoi1+9Y6HUVERMRVVOAECAxm2FlWw5cbi52O0mhS90kAuowqIiJyCBU4AeDsYZnER4a5ajBDekw6g5IH8fk2FTgREZGmVOAEgKhwLzNGZPHeip1U1/mcjtNoUvdJLNq9iKr6KqejiIiIuIYKnDS6cHQ2FbUNfLhql9NRGk3KnkSDv4EFOxc4HUVERMQ1VOCk0YS+qXRLjOI1F11GHZ0xmuiwaOZtm+d0FBEREddQgZNGHo9h5qhs/rN2D0UVtU7HASDCG8G4rHF8sV3LaomIiOynAicHuXhMNj6/5c2l7llaa2L3iWwp38LWsq1ORxEREXEFFTg5yKDMeIZ3T3DVZdTJ2ZMBTSciIiKynwqcHOai0dksLSxlw54Kp6MA0Cu+F9lx2SpwIiIiQSpwcpgLTuqOx+Cas3DGGCZ1n8TXO76m3lfvdBwRERHHqcDJYTISopg0II1XF29zzdJak7InUdVQxZI9S5yOIiIi4jgVOGnWxWOyKdxXTd7mfU5HAeDkrJMJM2FalUFERAQVODmCacOziInw8qpLLqPGRcRxUsZJmk5EREQEFTg5gpiIMKYNz+LtZdupqXfH0lqTsyezau8qiqqLnI4iIiLiKBU4OaKLRmdTVtPAp2t2Ox0FCMwHBzB/+3yHk4iIiDhLBU6OaGL/VNLjI/n3IndcRh2SMoSUqBRNJyIiIl2eCpwcUZjXw7k5Wcxdt4faBucvo3qMh4ndJ/LFti/wW7/TcURERByjAidHddqgdGrq/SwscMdo1IndJ7Kvdh+r9q5yOoqIiIhjVODkqCb0SyXca5i7zh0DB/bfB/fFNo1GFRGRrksFTo4qNjKMMb2S+WzdHqejAJAancrQlKHM2zbP6SgiIiKOUYGTYzptUDr528soqqh1OgoQWJVh2Z5lVNS5Y61WERGRjqYCJ8d06sA0AD5f757LqA22ga92fuV0FBEREUeowMkxDe+eSFJMOJ+55D64UemjiA2P1bJaIiLSZanAyTF5PYZJA9L4bN0erHV+cftwbzgnZ53MF9u/cEUeERGRjqYCJy1y2sA0dpXVsm63O+47m9R9EtsqtrG5bLPTUURERDqcCpy0yOSB6QDMXeuO0agTswPTiWhVBhER6YpU4KRFspOi6Z8e65r74HrG96R3Qm/dByciIl2SCpy02KkD0/lqU7ErltWCwGjUvF151PrcMb2JiIhIR1GBkxY7dWCaq5bVmtR9EtUN1SzatcjpKCIiIh1KBU5azG3Lao3LGke4J5wvtmtZLRER6VpU4KTF3LasVkx4DGMyxmggg4iIdDkqcNIqbltWa2L2RNbtW8euyl1ORxEREekwKnDSKm5bVmtS90kAuowqIiJdigqctMrw7okku2hZrUHJg0iPTleBExGRLkUFTlrF6zFMdNGyWsYYTul+Cl9s/wKf3x3Tm4iIiLQ3FThpNbctqzU5ezJldWXkF+c7HUVERKRDqMBJq7ltWa1Tup2CwWg0qoiIdBkqcNJqbltWKykqieGpw7WsloiIdBkqcHJc9i+rVVPvjvvOJmVPYnnRckprS52OIiIi0u5U4OS4nDYouKzWZpcsq5U9Cb/1M3/HfKejiIiItDsVODku4/sGltVyy2XUEWkjSI5M5uMtHzsdRUREpN2pwMlxcduyWmGeMKb0msLcwrnU+eqcjiMiItKuVODkuLltWa2zep1FZX0l87frMqqIiHRuKnBy3Ny2rNaEbhOIC4/joy0fOR1FRESkXanAyXHbv6zW3LXuKHDh3nBO73k6n2z9hHp/vdNxRERE2o0KnBw3r8cwyUXLagGc3etsSmtLyduZ53QUERGRdqMCJyfktIHp7C53z7JaE7MnEh0WzUebdRlVREQ6LxU4OSGTg/fBuWVZreiwaCZnT2bOljla3F5ERDotFTg5Id1dtqwWwNm9z6a4ppile5Y6HUVERKRdqMDJCXPbslqnZp9KuCecDzd/6HQUERGRdtHhBc4Y09MY84kxZpUxJt8Yc1Nwf4ox5kNjzLrgc3KTz9xqjFlvjFljjJnWZP9YY8zy4HsPGGNMR/884r5lteIi4pjYfSJztsxxzeAKERGRtuTEGbgG4KfW2qHABGC2MWYY8CtgjrV2IDAn+Jrge5cDw4HpwMPGGG/wux4BrgMGBh/TO/IHkYD9y2rNdcmqDABn9T6LHZU7yC/OdzqKiIhIm+vwAmet3WGtXRTcLgdWAdnATOCp4GFPARcGt2cCL1hra621m4D1wMnGmG5AgrV2vg2cZnm6yWekA8VGhjG2dzLzXHQf3JSeU/Aar0ajiohIp+ToPXDGmD7AaOArINNauwMCJQ/ICB6WDWxt8rHC4L7s4Pah+5v7c64zxuQZY/L27HHPWaLO5NSB7lpWKzEykXFZ4/hoy0e6jCoiIp2OYwXOGBMHvAL8xFpbdrRDm9lnj7L/8J3WPmatzbXW5qanp7c+rByT25bVgsBo1M1lm1lXss7pKCIiIm3KkQJnjAknUN6es9b+O7h7V/CyKMHn3cH9hUDPJh/vAWwP7u/RzH5xgNuW1QI4s9eZGIwuo4qISKfjxChUAzwOrLLW/rnJW28AVwW3rwJeb7L/cmNMpDGmL4HBCl8HL7OWG2MmBL9zVpPPSAdz47JaadFpjM4YrelERESk03HiDNwk4LvAmcaYJcHHucDvgbONMeuAs4OvsdbmAy8CK4H3gNnW2v0Tjl0P/IPAwIYNwLsd+pPIQfYvq7V2lzuW1YLAZdT1JespKC1wOoqIiEibcWIU6jxrrbHWjrTWjgo+3rHWFltrp1prBwaf9zb5zL3W2v7W2sHW2neb7M+z1uYE37vRuuXUTxe1f1mtz1w0ncjUXlMB+GiLLqOKiEjnoZUYpM24cVmtbnHdyEnN0X1wIiLSqajASZty27JaEJjUN784n+0VGuMiIiKdgwqctCm3LasFgQIHMGfLHIeTiIiItA0VOGlTblxWq3dCbwYmD9RlVBER6TRU4KRN7V9W6zMXzQcHcHavs1m8ezFF1e7KJSIicjxU4KTNnTownZU7ythT7o5ltSBwGdVimbNZl1FFRCT0qcBJm9u/rNYXG9xztmtA0gD6JPThwy2a1FdEREKfCpy0OTcuq2WM4azeZ5G3M4+SmhKn44iIiJwQFThpc/uX1Zq7bg9+v3vmVj6r11n4rI9Ptn7idBQREZETogIn7eKsoZnsKa9l8dYSp6M0GpY6jO6x3bUqg4iIhDwVOGkXU4ZkEO41vJ+/0+kojYwxTO09lfnb51NR5571WkVERFpLBU7aRWJ0OJMGpPHeip24aYnas3ufTb2/nrmFc52OIiIictxU4KTdTB+exZa9VazcUeZ0lEYnpZ9EWnSaLqOKiEhIU4GTdnP2sEw8Bt5f4Z7LqB7jYWqvqczbNo/qhmqn44iIiBwXFThpN6lxkZzcN4V3XVTgIDCpb3VDNZ9v+9zpKCIiIsdFBU7a1YycbqzbXcH63e4ZNJCbmUtSZBIfbtakviIiEppU4KRdTRueBeCq0ahhnjCm9JzC3MK51PnqnI4jIiLSaipw0q6yEqMY3SuJ91x4GbWivoIvd3zpdBQREZFWU4GTdjd9eBbLt5WydW+V01EaTeg2gbjwOD7arNGoIiISelTgpN1Nz3HfZdQIbwSn9TiNT7Z+QoO/wek4IiIiraICJ+2ud2osQ7sluO4y6tm9z6aktoS8XXlORxEREWkVFTjpENOHZ7Fwyz52l9U4HaXRxO4TifJG6TKqiIiEHBU46RAzRmRhLby/cpfTURrFhMcwOXsyH2/5GL/1Ox1HRESkxVTgpEMMzIijX1qsq1ZlgMBo1D3Ve1i6Z6nTUURERFpMBU46hDGG6TlZzN9YzL5K98y9dnqP0wn3hPN+wftORxEREWkxFTjpMNNzsvD5LR+tcs9l1LiIOKb0nMI7G9+h3lfvdBwREZEWaVGBM8aktOCR1M5ZJcSNyE4kOynaVdOJAMwcMJN9tfuYu22u01FERERaJKyFx20PPsxRjvECvU44kXRaxhimDc/i2a82U1HbQFxkS//za18Tu08kLTqN19a/xtReU52OIyIickwtvYS6ylrbz1rb90gPoLg9g0rnMGNEFnUNfj5ZvdvpKI3CPGF8o983+KzwM4qqi5yOIyIickwtLXCntNEx0sWN6ZVMWlyk6yb1vaD/Bfisj3c2vuN0FBERkWNqUYGz1h5z9tWWHCPi9RjOGZ7JJ2t2U1PvczpOowHJA8hJzeG1Da9hrXU6joiIyFEds8AZY35ojPm7MeZyY8xbxpjrOyKYdF4zcrKoqvMxd+0ep6McZOaAmazbt47Ve1c7HUVEROSoWnIG7kzgOuBGa+35wEntG0k6uwn9UkmMDuc9l41GndF3BuGecF7f8LrTUURERI6qJQWu2AauKd0XfF3bjnmkCwj3ejhraCYfrdxFXYN7lrBKjExkSs8pvL3xbc0JJyIirtaSAvdXAGvtm8HX/26/ONJVTM/JoqymgS83umvw8swBMympLWFuoeaEExER9zpmgbPWHnpD0BftlEW6kFMHphET4eVdl41Gndh9IunR6by24TWno4iIiBxRq5bSMsb8A9hljNlqjPkqOLjhR+2UTTqxqHAvU4Zk8OHKnfj87hn1GeYJ4/z+52tOOBERcbXWroV6KpBpre0JXAy8CsS2eSrpEmbkZFFUUUdewV6noxxkZv+ZmhNORERcrbUF7ksgGcBau81a+4619vdtH0u6gjMGZxAR5nHdaNT+Sf0ZkTZCc8KJiIhrtbbAPQb8xxjzM2PMqcaYxPYIJV1DXGQYpw1M4/0VO11XlGb2D8wJt2rvKqejiIiIHKa1Be5Z4EUgDLgB+MIYs6HNU0mXMT2nG9tLa1hWWOp0lINM7zs9MCfces0JJyIi7tPaAldorb3DWvt7a+23rbXDgWHtEUy6hrOGZhDmMa67jJoYmciZvc7k7U1vU+erczqOiIjIQVpb4JYYY25qusNaq4l95bglxURwSv9U3nPpZdTS2lLNCSciIq7T2gKXCfzQGLM9uC7qvcaYS9ojmHQd04ZnsamokrW7KpyOcpBTup9CenS6LqOKiIjrtKrAWWsvtdYOBfoCvwHWAie3RzDpOs4Znokx8O6KHU5HOUjjnHDbNCeciIi4S4sKnDHmFGOM2f/aWltrrV1krX3KWvvz9osnXUFGfBS5vZN5z2WrMgBc2P9CfNbH2xvfdjqKiIhIo5aegbsKWGiMecEYc7UxJqs9Q0nXMz2nG6t3llNQVOl0lIP0S+rHyLSRvLZec8KJiIh7tKjAWWt/aK0dA9xJYCLffxpj5htjfmeMOc0Y423PkNL5TRueCeC60agQWOB+fcl6zQknIiKu0dp74FZba/9irZ0OnAnMAy4BvmqPcNJ19EiOYUR2ousWtweY1mcaEZ4IDWYQERHXaO1i9uuMMf82xtwBnAOstNb+yFqb2z7xpCuZnpPF0q0lbC+pdjrKQTQnnIiIuE1rpxH5G7ATKAZmACuMMcuNMb81xoS3eTrpUqbnBG6tfN+ll1FLa0v5T+F/nI4iIiLS6gJ3pbX2Bmvtg9baHwKTgY+BMuDPbZ5OupT+6XEMyoxz5WjUU7qdQkZ0hi6jioiIK7S2wJUaY0buf2GtXQJMsNbeD0xqy2DSNU0fnsWCgr0UVbhrgQ+vx8v5/c9n3rZ5mhNOREQc19oC90PgSWPM48aYHxljHgT8wfci2jaadEXnjuyG38KbS7c7HeUwM/vP1JxwIiLiCq0dhbqKwMoL7wEZwHrgfGNMLPBC28eTrmZIVgI52Qm8lFfodJTDaE44ERFxi9aOQr3HWuuz1r5krf1va+3/A0qstZXW2nvaJ6J0NZfm9mTljjJWbCt1Osph9s8Jt3LvSqejiIhIF9baS6jZxpjv7H9hjEkHPmrbSNLVXXBSdyK8Hl5e6L6zcNP7TteccCIi4rjWFrj/An5gjDnZGDMO+AS4v+1jSVeWFBPBOcMzeW3JNmobfE7HOUhCRAJTe03lnU3vaE44ERFxTEsXs3/aGPMT4BRgNvAY8AhwobVWd3RLm7sktyclVfV8tHK301EOoznhRETEaS09A/dU8NjvAf8H9AH2AVcaY77VPtGkK5s8II1uiVG8tHCr01EOM6HbBM0JJyIijmrpYvZzrLV/ttZeZa0dBaQBPwU2ABPaMZ90UV6P4ZtjejB37R52ltY4HecgXo+Xb/T/huaEExERx7T2HjgArLUN1tpl1tpnrLU/a+tQIgDfGtsDv4VXFrlvMMMFAy7AZ328teEtp6OIiEgX1NJ74Ba1xTEirdEnLZaT+6bwUt5W18271i+xH6PSR/HS2pfwW/+xPyAiItKGWnoGbqgxZtlRHssJXFYVaVOXjO1BQXEVeZv3OR3lMN8Z+h22lG/hs8LPnI4iIiJdTFgLjxvSgmPcNd+DdArnjujGnW/k81LeVsb1SXE6zkHO6n0WGTEZPLPqGU7vebrTcUREpAtp6SCGzc09gGuAR4H/Bma2Z1DpmmIjwzhvZDfeWraDytoGp+McJNwTzreHfJuvdnzFun3rnI4jIiJdyHENYmgiCfgSuBcY3JIPGGOeMMbsNsasaLLvTmPMNmPMkuDj3Cbv3WqMWW+MWWOMmdZk/1hjzPLgew8YY8wJ/iziUpfm9qSqzsc7y3c4HeUw3xr4LSK9kTy36jmno4iISBdyogVuL+AFdge3W+KfwPRm9v/FWjsq+HgHwBgzDLgcGB78zMPGGG/w+EeA64CBwUdz3ymdwNjeyfRLi3XlAvdJUUmc3+983tr4Fvtq3HefnoiIdE6tLnDGGI8xJgzAWvtbApdQHwBatPK4tXYuLS97M4EXrLW11tpNwHrgZGNMNyDBWjvfBoYnPg1c2LqfREKFMYZvju3B1wV72VRU6XScw1w59EpqfbW8vPZlp6OIiEgX0aoCZ4y5EdgFbDbGLDXGXGut3W6tvdZa+5cTzHJjcETrE8aY5OC+bKDpVPyFwX3Zwe1D9x8p93XGmDxjTN6ePXtOMKY44ZtjeuAx8LILV2YYkDyACd0m8MLqF6j31zsdR0REuoDWnoH7KTDCWptN4JLlZGPMnW2Q4xGgPzAK2AH8Kbi/ufva7FH2N8ta+5i1Ntdam5uenn6CUcUJWYlRnDYonVcWbsPnd9eccADfHfZddlfv5sOCD52OIiIiXUBrC1wFgfvdsNbuAK4FLj7RENbaXdZan7XWD/wdODn4ViHQs8mhPYDtwf09mtkvndiluT3ZWVbDvPXuW75qcvZkeif01mAGERHpEK0tcI8ALxljBgRf9wKqTjRE8J62/S4C9o9QfQO43BgTaYzpS2CwwtfB8lhujJkQHH06C9DK4p3c1KEZJMWE82Ke+y6jeoyH7wz5DsuKlrF0z1Kn44iISCfXqgJnrX0YeA74hzFmL4FBBWuMMZcYYwa25DuMMc8D84HBxphCY8y1wB+CU4IsA6YANwf/vHzgRWAl8B4w21q7f8Lg64F/BDNsAN5tzc8ioScyzMuFo7L5MH8XJVV1Tsc5zMwBM4kLj+O5lToLJyIi7csc7xqTwZGow4DRBO5dO8lae2bbRWsfubm5Ni8vz+kYcpzyt5dy3gPzuOuC4Vw1sY/TcQ7zxwV/5LlVz/HeN98jKzbL6TgiIhLijDELrbW5h+4/7nngrLUN1tpl1tqnrLU3h0J5k9A3vHsiw7ol8JILR6MCfHvIt7FY/rXmX05HERGRTuxEJ/IV6XCX5vZgxbYyVm4vczrKYXrE92BKzym8tPYlqhuqnY4jIiKdlAqchJyZo7KJ8HpcexbuiqFXUFpbytsb33Y6ioiIdFIqcBJykmMjOHtYJq8t3kZdg9/pOIfJzcxlSMoQnlv1HMd7j6mIiMjRqMBJSPpWbg/2VdUzZ9Uup6McxhjDFUOvYH3Jer7c8aXTcUREpBNSgZOQdNrAdLISonhpofsWuAeY0XcGKVEpPLvqWaejiIhIJ6QCJyHJ6zFcPCabT9fsZldZjdNxDhPpjeTSwZcyt3Aum8s2Ox1HREQ6GRU4CVmX5PbEb+Hfi7Y5HaVZlw2+jDBPmJbXEhGRNqcCJyGrb1os4/ok81LeVlcOFkiLTmNGnxm8tv41yurcN+WJiIiELhU4CWmXjO3JxqJKFm3Z53SUZl0x7AqqG6p5dd2rTkcREZFORAVOQtq5I7sRE+HlpTx3DmYYnjqcMRljeH718/j8vmN/QEREpAVU4CSkxUWGce6Ibry5dDtVdQ1Ox2nWlcOuZFvFNj7d+qnTUUREpJNQgZOQd2luTyrrfLy7fKfTUZo1pecUusd215QiIiLSZlTgJOSN65NMn9QYXsxz59JaYZ4wvj3k2+TtymP13tVOxxERkU5ABU5CnjGGS3J78tWmvWwurnQ6TrMuGngR0WHRPLtSZ+FEROTEqcBJp3DxmGy8HsNzX21xOkqzEiMTuaD/Bbyz6R2KqoucjiMiIiFOBU46hW6J0czIyeL5r7ZQXlPvdJxmXTH0Cur99by09iWno4iISIhTgZNO4/un9qO8toEXXTqlSN/EvkzOnsy/Vv+LOl+d03FERCSEqcBJpzGqZxLj+iTzxLxNNPj8Tsdp1neHfpfimmLeL3jf6SgiIhLCVOCkU/n+qf3YVlLNe/nunFLklO6n0C+xH0/lP+XK5b9ERCQ0qMBJp3LW0Ez6pMbw9882ubIgGWP4Xs73WLNvDR9v+djpOCIiEqJU4KRT8XoM107uy9KtJeRtduf6qOf1O4/eCb15eOnD+K07L/WKiIi7qcBJp/OtsT1Jignn73M3Oh2lWWGeMP5r5H+xdt9aPtr8kdNxREQkBKnASacTHeHlyvG9+XDVLjYVuXNi33P7nkufhD48svQRnYUTEZFWU4GTTmnWxN6Eezw8MW+T01Ga5fV4uf6k61lfsp4PNn/gdBwREQkxKnDSKWXERzFzVHdeWriVfZXunHNtWp9p9EvsxyNLHsHn9zkdR0REQogKnHRa3z+1HzX1fp77arPTUZrl9Xi5ftT1bCzdqHnhRESkVVTgpNManBXPaYPSeWr+Zmob3HmG65ze5zAgaQCPLNVZOBERaTkVOOnUfnBqX/aU1/L6ku1OR2mWx3i4YdQNFJQV8M6md5yOIyIiIUIFTjq1yQPSGJIVz+MundgXYGqvqQxKHsTflv2NBn+D03FERCQEqMBJp2aM4fun9mPNrnLmrityOk6zPMbDDSfdwOayzby98W2n44iISAhQgZNO74KTupMRH8k/PnPnxL4AZ/Y6k6EpQ3UWTkREWkQFTjq9iDAPV03sw2frili9s8zpOM0yxnD9SdeztXwrb2540+k4IiLicipw0iVcMb4X0eFe/vGZOyf2BTij5xkMSx3G35b9jXp/vdNxRETExVTgpEtIiongktwevL5kG7vLapyO0yxjDLNHzWZbxTbeWP+G03FERMTFVOCky/jepL40+C1PzS9wOsoRnZp9KiPSRvDYsseo9+ksnIiINE8FTrqMPmmxnDMsk2e/3EJVnTsHChhjuGHUDWyv3M6r6191Oo6IiLiUCpx0KT84tR+l1fW8vLDQ6ShHNKn7JEamj+Tvy/9Onc+d67iKiIizVOCkSxnbO5lRPZN4fN4mfH53Tuy7/164nZU7eXWdzsKJiMjhVOCkSzHG8INT+7G5uIoPV+5yOs4RndLtFEZnjOax5Y9R66t1Oo6IiLiMCpx0OdOGZ9IjOdrVE/vuvxdud9VuXln7itNxRETEZVTgpMsJ83r43qS+5G3ex+It+5yOc0Tjs8YzNnMs/1j+D2oa3Dn1iYiIOEMFTrqkS8f1JD4qzNUT++6/F25P9R5eXvuy03FERMRFVOCkS4qLDOM743vx7oodbN1b5XScIxqXNY6Ts07m8RWPU91Q7XQcERFxCRU46bKuntgHjzE88bl7z8IB3DDqBoqqi3hxzYtORxEREZdQgZMuq1tiNN84qTsvLthKabV7Vz0YmzmW8d3G88SKJ6iqd+/ZQhER6TgqcNKlff/UvlTW+Xj2y81ORzmq2aNms7dmr87CiYgIoAInXdzw7olMGZzOY3M3uvos3OiM0UzsPpHHVzxOaW2p03FERMRhKnDS5f1s2mBKq+t5bO4Gp6Mc1S1jb6GsroyHljzkdBQREXGYCpx0ecO7J3LBSd15Yl4Bu8vdO9/a4JTBXDroUv615l+s3bfW6TgiIuIgFTgR4JazB1Hv8/Pgx+udjnJUN46+kfiIeH7/9e+x1p1ruYqISPtTgRMB+qTFctm4nvzfV1vYUuzekZ6JkYn8ePSPWbBzAR9s/sDpOCIi4hAVOJGgH08dSJjX8OcP1zgd5ai+OfCbDE4ezP1592tyXxGRLkoFTiQoMyGKayb15fWl21m1o8zpOEfk9Xi5dfyt7KzcyRMrnnA6joiIOEAFTqSJH57Wn/jIMO5/391n4cZmjmVG3xk8sfwJCssLnY4jIiIdTAVOpInEmHB+eEZ/5qzeTV7BXqfjHNUtY2/B6/Hyp7w/OR1FREQ6mAqcyCGumdiX9PhI7ntvtatHembFZvGDET/goy0fMX/7fKfjiIhIB1KBEzlEdISXH08dyIKCfXy6Zo/TcY5q1vBZ9Ijrwe+//j31fveuJCEiIm1LBU6kGZeP60nv1Bj+8P4a/H73noWL9Ebyy5N/ycbSjbyw+gWn44iISAdRgRNpRrjXwy1nD2LVjjLeXLbd6ThHdXqP05mUPYmHlzxMcXWx03FERKQDqMCJHME3RnZnSFY8f/5wLfU+v9NxjsgYwy/H/ZKahhoeWPyA03FERKQDqMCJHIHHY/jl9CFsLq7iXwu2Oh3nqPom9uXKYVfy6rpXWVG0wuk4IiLSzlTgRI7ijMHpjOuTzF/nrKO6zud0nKP6r5H/RUpUCv/z9f/gt+49YygiIidOBU7kKIwx/GL6EPaU1/LkF5ucjnNUcRFx3Dz2ZpbtWcZbG99yOo6IiLSjMKcDdDr/uhIqi5xO0bYGTIXTfu50CseM65PC1CEZPPrpBq44uTeJMeFORzqib/T/Bi+ueZG/LPwLZ/Y8k7iIOKcjiYhIO9AZuLbmCQdvJ3rUlsPH98CWr5z+m3XUz6YNpry2gUfnbnA6ylF5jIdbx99KUXURf1v2N6fjiIhIO9EZuLZ2yZNOJ2hbtRXwwGj46A645l0wxulEjhjaLYGZJ3Xnyc83cc3EPmQkRDkd6Yhy0nK4aMBFPLvyWS4aeBH9Evs5HUlERNpYh5+BM8Y8YYzZbYxZ0WRfijHmQ2PMuuBzcpP3bjXGrDfGrDHGTGuyf6wxZnnwvQeM6aLNor1FxsEZv4Qt82HNu06ncdQtZw+mwWd54ON1Tkc5ppvG3ERUWBR/+PoPrl4OTEREjo8Tl1D/CUw/ZN+vgDnW2oHAnOBrjDHDgMuB4cHPPGyM8QY/8whwHTAw+Dj0O6WtjLkKUgfAR3eCr8HpNI7plRrDt0/uxQtfb2VzcaXTcY4qNTqVG0bdwOfbP+c/hf9xOo6IiLSxDi9w1tq5wN5Dds8EngpuPwVc2GT/C9baWmvtJmA9cLIxphuQYK2dbwOnF55u8hlpa95wmPobKFoDS//P6TSO+tGZAwj3evjzh2udjnJMlw+5nH6J/bjv6/uo9dU6HUdERNqQWwYxZFprdwAEnzOC+7OBpjOoFgb3ZQe3D93fLGPMdcaYPGNM3p497l6c3LWGXgDZufDJ76Cuyuk0jslIiOJ7k/vw+pLt5G8vdTrOUYV7wvnlyb+ksKKQZ1Y+43QcERFpQ24pcEfS3H1t9ij7m2Wtfcxam2utzU1PT2+zcF2KMXD2b6F8B3z1iNNpHHXdaf1JjA7n/vfXOB3lmCZ2n8jUXlN5bNljbC1z92oSIiLScm4pcLuCl0UJPu8O7i8EejY5rgewPbi/RzP7pT31mQSDpsO8/wdVh14F7zoSo8O5/oz+fLJmD19vcv/fw69O/hVe4+X2z2/H53f3ahIiItIybilwbwBXBbevAl5vsv9yY0ykMaYvgcEKXwcvs5YbYyYER5/OavIZaU9T74C6Cph7v9NJHHXVKX3ITIjknrdX4vO7e5RnVmwWvzr5VyzavYhnVz3rdBwREWkDTkwj8jwwHxhsjCk0xlwL/B442xizDjg7+BprbT7wIrASeA+Yba3dfwrheuAfBAY2bAC69hwXHSVzGJz0HVjwd9i32ek0jomO8HL7ecNYVljKM/MLnI5zTBf0v4ApPafwwKIH2Fiy0ek4IiJygkxXmyMqNzfX5uXlOR0jtJVug/8dA8NmwsWPOZ3GMdZarn5yAQs37+PDW06jW2K005GOqqi6iItev4jsuGyeOfcZwj3uXRJMREQCjDELrbW5h+53yyVUCSWJ2TD+h7DsRdixzOk0jjHGcM+FOTT4/dz1xkqn4xxTWnQa/z3hv8kvzufx5Y87HUdERE6ACpwcn8k/gajEwBJbXVjPlBhumjqI9/J38uHKXU7HOaZz+pzDjL4z+NvSv7GqeJXTcURE5DipwMnxiU6G034GGz6GDZ84ncZR3z+1L4Mz47nj9RVU1rp/pYrbxt9GclQyv573a+p8dU7HERGR46ACJ8dv3A8gsWfgLJzf73Qax4R7Pfzu4hFsL60JiRUaEiMTuXPinawvWc9DSx5yOo6IiBwHFTg5fuFRMOU22LEU8v/tdBpHje2dzBXje/Hk55tYsc3dKzQAnNbjNL458Jv8M/+fLNm9xOk4IiLSSipwcmJGXgqZOTDnt9DQtS/H/WL6EFJiI7n138tdPzccwM9yf0ZWTBa3zbuNqvquuzyaiEgoUoGTE+Pxwll3QclmyHvC6TSOSowO545vDGP5tlKeDoG54eIi4rhn8j1sKd/C/1v0/5yOIyIiraACJyduwFTocyrM/QPUlDmdxlHnj+zG6YPSuf/9NeworXY6zjGNyxrHlUOv5PnVz/Plji+djiMiIi2kAicnbv9C91XF8PlfnU7jqP1zw/ms5c438p2O0yI/HvNj+iT04Tef/4byunKn44iISAuowEnbyB4Dwy+G+Q9B+U6n0zhq/9xw7+fv4oN89/9dRIdFc+/ke9lVtYs/Lvij03FERKQFVOCk7Zx5O/jr4dP/cTqJ475/al+GZMVzxxv5VITA3HAj00dybc61vLr+VT7d+qnTcURE5BhU4KTtpPaH3O/Bomdgj/vnQ2tP4V4P9140gp1lNfwlBOaGA/jhST9kUPIg7vziTkpqSpyOIyIiR6ECJ23rtF9AeDTMucvpJI5rOjfc8kL3zw0X4Y3gd5N/R2ldKfd+da/TcURE5ChU4KRtxaXDpJtg9Vuw5Sun0zju59OGkBoXya2vLqPB5/7VKganDOaGk27gvYL3eG/Te07HERGRI1CBk7Y34QaIzYAP/7tLL7EFB+aGW7GtjKfnb3Y6Totck3MNI9JGcM9X97Cnao/TcUREpBkqcNL2IuNg6n/D1q9g0T+dTuO480Z044zB6fzpgzVsL3H/3HBhnjDumXwPNQ013Dn/Tqx1/6oSIiJdjQqctI/R3w1M7vvBb6B0m9NpHGWM4e6ZoTU3XL/Eftw05ibmFs7l2VXPOh1HREQOoQIn7cMYuOAB8DfAWzdDFz+L0zMlhp+cNYgPVu7i/RCYGw7gyqFXMrXXVP6U9ycW7FzgdBwREWlCBU7aT0q/wKXUde/D8pedTuO4aycH5oa7M0TmhjPGcM+ke+gZ35Of/edn7KwMjeIpItIVqMBJ+xr/Q8jOhXd/ARVd+4b4pnPD3f/+GqfjtEhcRBx/nfJXahpquOXTW6jz1TkdSUREUIGT9ubxwswHobYc3vul02kcN7Z3MrMm9OafXxTwyZrdTsdpkX5J/fjd5N+xvGg5//O1VtkQEXEDFThpfxlD4bSfw4pXYPU7Tqdx3K3nDmVIVjw/fXEpu8pqnI7TIlN7T+X7I77Py2tf5pW1rzgdR0Sky1OBk44x+WbIGA5v3wLVJU6ncVRUuJcHvzOa6jofP3lhCT5/aAzwuHHUjUzsPpF7v7qX5XuWOx1HRKRLU4GTjhEWATP/Fyp2BSb47eIGZMRz18zhzN9YzEOfrHc6Tot4PV7uO/U+MmIyuPnTmymuLnY6kohIl6UCJx0neyycciMseho2fup0GsddMrYHF47qzv/7aC1fbQyNMpQUlcRfzvgLJbUl/Hzuz2nwu380rYhIZ6QCJx1ryq8D04u88WOoq3Q6jaOMMdxz0Qh6p8Zy0wtL2FsZGiM8h6YO5Y5T7mDBzgX8ZeFfnI4jItIlqcBJxwqPhgsehJLN8PG9TqdxXFxkGP/77dHsrazj5y8tDZllq77R/xt8Z8h3eHrl07y76V2n44iIdDkqcNLx+kyC3Gvhy4dhq2b4z8lO5NZzhzBn9W6e+LzA6Tgt9rNxP2NMxhju+OIO1uwNjXntREQ6CxU4ccZZd0JCNrxxIzTUOp3GcVdP7MPZwzL5/burWFZY4nScFgn3hPOnM/5EXHgcN396M6W1pU5HEhHpMlTgxBlRCXD+X2DPaph7v9NpHGeM4Y/fGkl6XCQ/en4x5TX1TkdqkbToNP58xp/ZUbmDWz+7Fb/1Ox1JRKRLUIET5ww6B0ZeBvP+DDtXOJ3GcUkxETzw7dEU7qvm16+uCJn74UZljOJX437FZ9s+45GljzgdR0SkS1CBE2dN/z1EJcHrs8GnKSly+6Rwy9mDeHPpdv61YKvTcVrs0sGXcuGAC3l06aN8uvVTp+OIiHR6KnDirJgUOPePsGMJfPmQ02lc4frT+zN5QBp3vpnP2l3lTsdpEWMMt42/jWGpw7j1s1vZXLbZ6UgiIp2aCpw4b/hFMOR8+OR3UBQaqxK0J4/H8OfLTiIuMowb/28R1XU+pyO1SFRYFH854y+Ee8L5ySc/oaq+yulIIiKdlgqcOM8YOPd+8EbCmz8Gv26Ez4iP4s+XjmLtrgp++1a+03FarHtcd/5w+h/YWLqRn/7np9T7Q2MwhohIqFGBE3dI6AbT7oXNn8PCJ5xO4wqnDUrn+jP68/zXW3lz6Xan47TYhG4T+M2E3zBv2zxun3e7RqaKiLQDFThxj9FXQr8z4MM7YF+B02lc4ZazBzGmVxK3/ns5m4tDZ+mxbw76JjeNuYl3Nr3DHxb8IWRG1IqIhAoVOHEPY+AbD4DxwouzoL7G6USOC/d6eODbo/EY+NHzi6lrCJ2zWdfmXMt3h32X51Y9x2PLHnM6johIp6ICJ+6S3Bsu/hvsWArv/MzpNK7QIzmGP3zrJJYVlvL7d1c7HafFjDH8LPdnfKPfN3hwyYO8uOZFpyOJiHQaKnDiPoNnwGk/h8XPwMKnnE7jCtNzsrh6Yh+e+HwTz3+9xek4LeYxHu6adBen9TiNe768h/cL3nc6kohIp6ACJ+50xq3Q/8zAWbhtC51O4wq3nzeU0welc/trK/h0zW6n47RYuCec+0+/n9EZo/nVZ79i/vb5TkcSEQl5KnDiTh4vfPNxiMuCF6+CymKnEzkuzOvhoSvGMDgzntnPLSJ/e+gsHh8dFs3/Tv1f+ib25aZPbmJFkZZOExE5ESpw4l4xKXDpU1CxG165FvyhMaFte4qLDOPJa8aREB3O9/65gB2l1U5HarGEiAT+dtbfSIlK4fqPrmdj6UanI4mIhCwVOHG37DFw3v2w8ZPASg1CZkIUT14zjspaH9c8uYDymtCZLDc9Jp3Hzn4Mr/Fy3QfXsbNyp9ORRERCkgqcuN+YWYHHZ/fD6necTuMKQ7ISeOTKMazfXcENzy2i3hc604v0SujFo2c/SmV9Jdd9eB37avY5HUlEJOSowElomPFH6DYKXv0vKN7gdBpXOHVgOr+7eASfrSvi9ldXhNRkuUNShvDAmQ+wrXwbs+fM1rqpIiKtpAInoSE8Ci57JjC44V/fhbrQWZWgPV2a25MfnzmAf+Vt5aFP1jsdp1XGZY3jj6f/kfzifH7yyU+o89U5HUlEJGSowEnoSOoVGJm6eyW8eROE0Bmn9nTz2YO4aHQ293+wltcWb3M6Tquc2etM7jzlTubvmM+v5/0anwaqiIi0iAqchJYBU+HM22D5S/D1351O4wrGGO775kgm9Evh5y8v5cuNoTXlykUDL+KWsbfwfsH7/M/X/xNSl4JFRJyiAiehZ/JPYdAMeP9W2PKV02lcISLMw9+uzKV3aizXPZ3H+t3lTkdqlWtyruGa4dfwrzX/4i8L/6ISJyJyDCpwEno8HrjoUUjsCS9dBeW7nE7kCokx4Tx59Tgiwrxc/eQC9pTXOh2pVW4eezOXDb6MJ/Of5J4v78FvQ2dkrYhIR1OBk9AUnQSXPQvVJfDyNeALnbnQ2lPPlBieuDqX4oo6rn1qAVV1DU5HajFjDLeNv43v5XyPF9e+yK2f3Uq9X/93FRFpjgqchK6sHPjGX2Hz5/DRnU6ncY2RPZL432+PZsW2Un78/BJ8/tC5HGmM4eaxN3PTmJt4Z9M73PzJzdQ01DgdS0TEdVTgJLSddBmcfB3MfxDyX3U6jWucNSyTO74xnI9W7eLut1aG3D1l3x/xfW4ffztzC+cye85sKus1bYyISFMqcBL6zrkXepwMr82GXSudTuMaV03sw/cn9+WfXxTw+LxNTsdptcuGXMbvTv0dC3ct5Pvvf5+SmhKnI4mIuIYKnIS+sIjAoveR8fDMRbA39MpKe/n1uUOZkZPFPW+v4tkvNzsdp9XO73c+/2/K/2PtvrVc8/417K7a7XQkERFXUIGTziGhO8x6DXy18PRMKNvhdCJX8HgMf7lsFFOHZHD7ayt48vPQK7dn9DyDh896mG0V27jq3asoLC90OpKIiONU4KTzyBgKV74CVcXwzIVQtdfpRK4QFe7lkSvHMm14Jne9uZLH5obeWrLju43nH+f8g7K6Mq569yo2lITezyAi0pZU4KRzyR4L334+cBn12W9CbWhNaNteIsI8PPidMZw/shu/e2c1D368zulIrTYyfST/nP5P/Pi5+r2ryS/KdzqSiIhjVOCk8+l7WuCeuB1L4flvQ32104lcIdzr4f9dNoqLg+um/vnDtSE3OnVg8kCenv40seGxXPvBtSzYucDpSCIijlCBk85p8Ay46G9QMA9e0kS/+4V5PfzxkpO4NLcHD8xZx33vrQm5EtczoSdPTX+KzJhMrv/oeuYWznU6kohIh1OBk85r5CVw3v2w9l147Xrwa2kmAK/H8PuLR3LlhF48+p8N3P3WqpArcZmxmTw5/Un6Jfbjpo9v4r1N7zkdSUSkQ4U5HUCkXY37PtSUwZy7IDIBzvsTGON0Ksd5PIa7Z+YQ7vXwxOebqPf5ueuC4Xg8ofN3kxKVwuPTHufGOTfyi7m/oLS2lMuGXOZ0LBGRDqECJ53fqbdATSl8/v8gKhHOusPpRK5gjOE35w8jIszD3/6zkXqfn99dNCKkSlx8RDyPnv0oP/vPz7jnq3tYV7KOX477JeHecKejiYi0KxU46RrOujNQ4ub9GaISYPLNTidyBWMMv5o+hEivhwc+Xk+dz88fv3US3hAqcdFh0Tww5QH+uuivPJn/JOtL1vOn0/9EanSq09FERNqN7oGTrsGYwOXTnG8FFr7Pe8LpRK5hjOGWcwZzy9mD+Peibdz8ryU0+ELrfkGvx8stubfw+1N/z4qiFVz+9uWsLNayaiLSeanASdfh8cJFj8LAafDWLbD8ZacTucqPpw7kVzOG8MbS7fzo+cXUNYRWiQM4r995PDXjKQBmvTuLdza+43AiEZH24aoCZ4wpMMYsN8YsMcbkBfelGGM+NMasCz4nNzn+VmPMemPMGmPMNOeSS8jwhgfmiOs9CV79L1ij0YtN/fD0/vz3+cN4d8VObnhuIbUNPqcjtdrw1OG8cN4LDE8dzi8/+yV/zvszPn/o/RwiIkfjqgIXNMVaO8pamxt8/StgjrV2IDAn+BpjzDDgcmA4MB142BjjdSKwhJjw6MBqDVkj4KWrYNNnTidylWsn9+XuC3P4aNVurnt6IVV1DU5HarXU6FT+cc4/uGzwZTyZ/yQ3zLmB0tpSp2OJiLQZNxa4Q80EngpuPwVc2GT/C9baWmvtJmA9cHLHx5OQFJUAV7wCyX3g+cuhMM/pRK7y3Qm9ue+bI/hs3R6++ch8CvdVOR2p1cK94dw+4XbuOOUOvt75Nd9++9us37fe6VgiIm3CbQXOAh8YYxYaY64L7su01u4ACD5nBPdnA1ubfLYwuO8wxpjrjDF5xpi8PXv2tFN0CTmxqfDdVyE2DZ76Bqx93+lErnLZuF48cfU4CvdVceFDn5NXsNfpSMflW4O+xRPTnqCqvoor3rmCOVvmOB1JROSEua3ATbLWjgFmALONMacd5djm5jlodjp5a+1j1tpca21uenp6W+SUziKhO3zvA0gbGFg3ddHTTidylTMGZ/Da7EnER4Xz7b9/yYt5W4/9IRcanTGaF85/gX6J/fjJJz/hkSWP4LehN0hDRGQ/VxU4a+324PNu4FUCl0R3GWO6AQSfdwcPLwR6Nvl4D2B7x6WVTiM+E65+G/qdAW/8CD79PYTY0lLtqX96HK/dMInxfVP5xcvLuPutlSE3zQhAVmwW/5zxTy7ofwEPL32Ymz+5mcr6SqdjiYgcF9cUOGNMrDEmfv82cA6wAngDuCp42FXA68HtN4DLjTGRxpi+wEDg645NLZ1GZDx8519w0nfg0/+BN38MvtC7eb+9JMaE889rxnH1xD48Pm8T33sqj9LqeqdjtVqkN5J7Jt3DL8b9gv8U/ocr3r6CLWVbnI4lItJqrilwQCYwzxizlEARe9ta+x7we+BsY8w64Ozga6y1+cCLwErgPWC2tVZzBcjx84bDhQ/DaT8PXEp94TtQpzM0+4V5Pdx5wXD+5+IRzN9QxEUPf86motD7+zHG8N1h3+XRsx+lqKaIS968hFfXvYrVWVcRCSGmq/2PVm5urs3L04hDOYa8J+Dtn0K3UfCdFyFO90429dXGYq5/bhENPj8PXTGGUweG5t/P9ort3DbvNvJ25TG111R+c8pvSIlKcTqWiEgjY8zCJlOrNXLTGTgR98j9Hlz2HOxeBY+fDcUbnE7kKuP7pfL67El0T4rm6icX8OTnm0LyDFb3uO7845x/8NOxP2Vu4Vwufv1i5hbOdTqWiMgxqcCJHMmQc+GqN6CmFB4/B7YtdDqRq/RMieHl6ydy5pAM7npzJb9+dXlILr/l9Xi5Oudqnj/veVKiU5g9ZzZ3z7+bqvrQm/tORLoOFTiRo+l5Mlz7AUTEwD/P11xxh4iLDONvV47lxikDeP7rrVz5j68orqh1OtZxGZwymOfPe56rh1/NS2tf4tK3LmXZnmVOxxIRaZYKnMixpA2Eaz+C1AGaK64ZHo/hZ9MG89fLR7G0sIQLHvycVTvKnI51XCK9kfw096c8Pu1x6nx1zHp3Fg8veZh6f+iNuBWRzk0FTqQl4jPhmneg3+nBueLu01xxh5g5KpsX/+sUGvx+vvnIF7y+ZJvTkY7buKxxvHLBK5zb91weWfoIs96ZRUFpgdOxREQaqcCJtFRkfGBE6knfhk9/B2/epLniDnFSzyTeuHEyQ7slcNMLS/jR84sprQrNs1fxEfH87tTfcf/p97O1YiuXvHkJ/1r9r5AcrCEinY8KnEhreMPhwkfg1J/BoqfguW9Cxe5jf64LyUyI4l/XTeBn5wzi3eU7mP7XuXy+vsjpWMdtWp9p/PuCfzMmcwz3fHUPN8y5gaLq0P15RKRzUIETaS1jYOp/wwX/C1u+hEcnw8ZPnU7lKmFeDzeeOZB/3zCR6AgvV/zjK+5+ayU19aE513ZGTAaPnvUot558Kwt2LuCi1y/io80fOR1LRLowFTiR4zVmFnx/DkQlwtMXwsf36pLqIUb2SOLtH53KrFN68/i8TVzw4DxWbg/NAQ7GGL4z9Du8+I0X6R7XnZs/vZkfzfkRW8u3Oh1NRLogrcQgcqJqK+Cdn8PS/4Pek+Cb/4CE7k6ncp1P1+zm5y8vo6Sqjp+eM5gfnNoPr8c4Heu41Pvr+b9V/8fDSx6mwd/A90d8n2tyriEqLMrpaCLSyRxpJQYVOJG2suT5wPJbYZFw0d9g0DlOJ3KdvZV1/Prfy3kvfycn903hz5eeRI/kGKdjHbddlbv4U96feLfgXXrE9eDW8bdyWo/TnI4lIp2IClyQCpy0qz1r4eVrYNcKmPgjmHpHYOCDNLLW8vLCQu56cyUGuPOC4Vw8JhtjQvNsHMBXO77id1/9jo2lGzmjxxn88uRf0iO+h9OxRKQTUIELUoGTdldfDe/fBnmPQ3YufOsJSO7tdCrX2bq3ilteXMKCgn2cOyKLey8cQXJshNOxjlu9r55nVz3LI0sfwW/9XDviWr6X8z0ivZFORxOREKYCF6QCJx0m/1V448eAgZkPwrALnE7kOj6/5W9zN/CXD9eSHBPB/ZecxGmD0p2OdUJ2Vu7kT3l/4r2C93RZVUROmApckAqcdKi9m+Dl78H2RXDydXD23RCuG90PtWJbKT/51xLW767gstye/GL6YFLjQvvM1Zc7vuR3X/2OTaWbOKPnGfxynC6rikjrqcAFqcBJh2uogzl3wfwHIWsEfOufkDbA6VSuU1Pv408frOGJzwuIjfDys2mD+c7JvQjzhu5sR/W+ep5Z9QyPLn0Uv/U3jlbVZVURaSkVuCAVOHHMmnfhtevBVw/n/jGwJFcI37jfXtbtKufON/P5fH0xQ7Li+e3MHE7um+J0rBOys3Inf1zwRz7Y/AHZcdlcf9L1nN/vfLwer9PRRMTlVOCCVODEUaWF8Mr3Yct86Hs6nPdnnY1rhrWWd1fs5J63VrK9tIaZo7rz63OHkpkQ2pef52+fz18W/oVVe1fRN7EvN4y6gXN6n4PHhO5ZRhFpXypwQSpw4ji/DxY+CR/9FhqqYfItMPlm3RvXjKq6Bh75dAN/m7uRcI/hx1MHcs2kvkSEhW7hsdYyZ8scHlz8IBtKNzA4eTA3jr6R03ucHtJTqYhI+1CBC1KBE9co3wUf3AbLX4KU/nDen6D/FKdTudLm4krufmslH63aTb/0WO78xvCQH63q8/t4Z9M7PLL0EbaWb2Vk2khuHH0jE7pNUJETkUYqcEEqcOI6Gz4OrOCwdyOMuATOuRfiM51O5UqfrN7NXW/mU1BcxTnDMvnv84fRMyV0V3KAwLJcr69/nUeXPsquql2MyxrHj0b/iNEZo52OJiIuoAIXpAInrlRfA/P+AvP+DGHRcNZvYOz3wBO6lwrbS22Dj8fnbeJ/56zHby3Xn9GfH57en6jw0B4QUOur5eW1L/P3ZX+nuKaYydmTuXH0jQxPHe50NBFxkApckAqcuFrRenj7Ftj0n8AqDuf/BbqNdDqVK+0oreZ376zmzaXb6ZEcza0zhjIjJwuPJ7QvP1bVV/H86ud5YsUTlNWVcVavs5g9ajYDkjXYRaQrUoELUoET17M2cF/c+7+GqmIYfz1MuRUi451O5krzNxRz5xv5rNlVzpCseG6aOpBpw0O/yJXXlfPsymd5auVTVNVXcW6/c7k251oGJg90OpqIdCAVuCAVOAkZ1ftgzm8h70lI6A4z7oMh52vuuGb4/Ja3lm3nr3PWsXFPJUOy4vnJWQM5Z1joF7mSmhKezH+S51c/T3VDNRO7T2TWsFlM7D5Rgx1EugAVuCAVOAk5WxfAWzfDruUw4GyY+t/Q7SSnU7mSz295c+l2Hpizjo1FlQztlhA8I5cZ8mWnpKaEl9a+xP+t/j+KqosYkDSAWcNmcW6/c7Wyg0gnpgIXpAInIcnXAF89CnP/ADWlMGwmnPFryBjidDJXavD5eXPZdh6Ys55NRZUM65bATWcN5JxhoV/k6nx1vFfwHk/nP82afWtIiUrh8iGXc9ngy0iJCu0VK0TkcCpwQSpwEtKqS+DLh2H+w1BXASMvhdN/Can9nU7mSg0+P28s3c7/fhwocsO7B87Ind0Jipy1lq93fs3TK59mbuFcIr2RnN/vfGYNm0W/pH5OxxORNqICF6QCJ51CZTF88Vf46jHw1cGo78Dpv4CkXk4nc6UGn5/Xl2znfz9eR0FxFcO7J/CTswZx1tCMkC9yABtLN/LMymd4c8Ob1PpqOTX7VGYNn8X4rPGd4ucT6cpU4IJU4KRTKd8VmDsu74nA6NWxV8OpP4WEbk4nc6UGn5/XgkVuc3EVOdkJ3DhlAGcNzSTMG/pz7u2t2cuLa17k+dXPs7dmL4OSBzFr2Cym9ZlGVJiWahMJRSpwQSpw0imVFsLcP8LiZ8ETBuO+H1hfNTbN6WSudGiRy06K5soJvbl8XE+SYyOcjnfCan21vLPxHZ5e+TTrS9YTHxHPuX3P5aKBFzEsZZjOyomEEBW4IBU46dT2boL/3AfL/hVY0WHC9TDxRohOdjqZK/n8lo9W7eKfnxcwf2MxkWEeZo7qzqxT+pCTneh0vBNmrWXBzgW8uv5VPtz8IbW+WgYnD+aigRdxXt/zSIpKcjqiiByDClyQCpx0CXvWwKf/A/mvQmQinDIbxl2rM3JHsXZXOU99UcC/F22jut5Hbu9krprYh+k5WYR3gsurZXVlvLfpPf697t/kF+cT7gnnzF5nctGAi5jQbQJeT2gvRSbSWanABanASZeyczl88jtY8w54I2HEJTD+Os0jdxSl1fW8lLeVZ77czObiKjITIrlifG++fXIv0uM7x3xra/au4bX1r/HWxrcoqS0hKzaLmf1ncuGAC+kR38PpeCLShApckAqcdEm7V8PXj8HS56G+CnpNhPH/FVjZwRvmdDpX8vstn67dzT+/2MzctXuI8Ho4b2Q3Zp3Sm9G9Oscl6TpfHZ9s/YRX17/KF9u+wGIZnzWeiwZexNReUzXwQcQFVOCCVOCkS6veB4ufC5S5ks2Q0CNwaXXMVRCb6nQ619qwp4Jn5m/m5YWFVNQ2cFKPRGad0ocZI7KIiegcBXhn5U5eW/8ar61/jW0V24gPj2dq76mc0/scJnSbQLg33OmIIl2SClyQCpwI4PfB2vcDqzts+g+ERQUvr/4XZI1wOp1rVdQ28O9FhTz1RQEb9lQSE+Fl2vAsLhqdzcT+qZ1iKhK/9ZO3M4/X1r/GJ1s/oaK+gviIeKb0nMK0PtOY0G0CEd7QH6krEipU4IJU4EQOsWtl8PLqC9BQDb0nB4rc4HN1efUIrLUsKNjHq4u38fay7ZTVNJAeH8kFJ3XnotHZDO+e0Cmm6qjz1TF/+3w+2PwBn2z5hPL6cuLD4zmj5xmc0+ccJnafqDIn0s5U4IJU4ESOoHpfYB65rx+Dki2Q2DNweXXk5ZoY+Chq6n18umY3ry7exserd1PvswzMiOPC0dnMHNWdHskxTkdsE3W+Or7c8SUfFHzAx1s/pryunLjwOE7veTrn9D6HSdmTiPR2jkEeIm6iAhekAidyDH4frH0veHl1LmCg72kw8jIY+g2ISnA6oWuVVNXxzvKdvLq4kAUF+wAY3zeFi0ZnM2NENxKjO8d9ZPW+er7a+VVjmSutLSUmLIbTe57OtN7TmJQ9SQMgRNqIClyQCpxIKxSth+UvBSYG3rcpcK/c4BmBMtd/KoTp8tmRbN1bxetLtvHvxdvYuKeSiDAPU4dkMHNUd04blN5pBj/U++tZsGMBH2z+gDlb5lBSW0KkN5LcrFxOzT6VydmT6Z3Q2+mYIiFLBS5IBU7kOFgLhXmw/EVY8QpUFQdWdxh+MYy8FHqOh05wz1d7sNayfFspry7exptLt1NUUUdEmIfJA9I4a2gmU4dmkJnQOc5W1fvryduZx9zCuczbNo+CsgIAesb3ZHL2ZCZnT2Zc1jiiw6KdDSoSQlTgglTgRE6Qrx42fAzLXoTVbwcGPiT1DoxiHXkZpA9yOqFrNfj8fF2wl49W7ubDVTvZurcagJN6JHL2sEzOGpbJ4Mz4TjEAAmBr+VbmbZvHvG3z+HrH19T4agJn5zJzGwtd74TenebnFWkPKnBBKnAibai2PFDilv0LNn4K1h9Y5WH//XJJvZxO6FrWWtbuquCjVbv4cOUulmwtAaBHcjRnDc3k7GGZnNw3pVMs4wVQ66tl4c6FfLbts4POzvWI68Hk7Mmc2uNUcjNziQnvHIM+RNqKClyQCpxIOynfCSv+HShzO5YE9mUMh8HTYdB0yB4LWm/ziHaX1TBn9W4+WrmLeeuLqG3wEx8VxpTBGZw1LJPTB6V3mkEQAIXlhQfOzu38muqGasJMGMPThpObmUtuVi6jM0YTGx7rdFQRR6nABanAiXSAovWw9t3AZMGbvwDrg5g0GHgODJoG/c/UaNajqKprYN66Ij5cuYuPV++muLKOMI/hpJ5JTOqfyin90xjdK4mo8M5RiGt9tSzctZCvd3xN3q488ovyabANeI2XoSlDGZs5ltysXMZkjiEhQv/dSNeiAhekAifSwar3wfo5gTK37gOoKQFPOPSZFDgzN2gapPRzOqVr+fyWJVv38dGq3XyxoZjlhSX4LUSGeRjXJ4VT+qcysX8qI7ITO8VKEABV9VUs3bOUvF155O3MY3nRcur99RgMg1MGB87QZQYKXXJU51iXVuRIVOCCVOBEHORrgMKvA/PMrXkPitYE9qcNDhS5QdOh58mgdTePqKymnq827uWLDUXM31DM6p3lAMRHhjG+Xwqn9E9j0oBUBmXE4/F0jsEBNQ01LC9aTt7OPPJ25bF0z1JqfbUADEgawJiMMYxIH8GItBH0TeyLx3SOIisCKnCNVOBEXGTvRlj7QeBya8Hn4K+H8FjoNSFwhq7PqdB9tArdURRV1DJ/QzFfbChm/oYiCoqrAEiNjWBC/1RO6ZdKbp9kBmbE4+0kha7OV0d+cf5Bha6yvhKA2PBYclJzyEnLYUTaCEakjyAjJsPhxCLHTwUuSAVOxKVqymDjJ7DpMyiYB3tWBfaHxwTmmeszOfDoPkYTCB/FtpJqvlgfODv3+YYidpUFzlTFRYZxUs9ERvdMZnSvJEb1TCI1rnMsfeW3fgpKC1hWtIwVRStYXrSctXvX0mAbAMiIyWBE2ghy0nIYmTaSYanDiIuIczi1SMuowAWpwImEiMoi2Px5oMwVfA678wP7w6IDl1n7nBo4S5c9FsI6RxFpa9ZaCoqrWLxlH4u3lLBoyz5W7yzH5w/8736f1BhG9woUutE9kxnSLb5TTVuyeu/qxkK3fM9ytpRvAcBg6JfYj+FpwxmcPJjBKYMZlDxI99OJK6nABanAiYSoymLY8kWw0M2DXSsC+8OioMc46HUKZI8JnKGLz3Q2q4tV1TWwvLCUxVtLWLR5H4u2lFBUEThLFxXuYWR2UqDQ9UpiePdEeiRHd5qJdktrSxsL3YqiFeQX51NUXdT4fkZ0BoNSBjEoeRCDkwOlrk9iH8I8nWPZMwlNKnBBKnAinUTV3sAUJQXzYPM82JUfmEgYICE7cO/c/kLXfTREJzka162stWwrqW48Q7d4Swn520up9wX+bYiPCmNYtwSGdktgWPcEhnVLYGBmHJFhnWMKk+LqYtbuW8vafWtZs3cNa/etZUPpBhr8gcuvEZ4I+if1D5S64Jm6wcmDSYpKcja4dBkqcEEqcCKdVG0F7FwG2xbB9kWB532bDryf0v9AocseA1kjIUKz/jenpt7Hyh1lrNxexsodZazaUcbqHeVU1/sACPMYBmTEBUpdsNgN7ZZASmznuDex3lfPxtKNjcVuf7krriluPCY1KpW+iX0PevRL7EdWbJZGwUqbUoELUoET6UKq9sL2xcFCF3wu3xF4z3ghY2hg6a+MYYHtjGEQnwWd5JJhW/L5LZuLKxuL3aodgXK3f5AEQFZCFEO7xTM4K4H+6bH0z4ijf3pcp1lBoqi6KFDo9q5lY+lGNpVuYmPpRsrqyhqPifJG0SexD30TDi53vRN6ExUW5WB6CVUqcEEqcCJdXNmOA2foti+CnSugcveB96OSAkUus0mpyxgK0brBvTnFFbWs2lHOyh2lgeftZWwsqmi8BAuQHh/JgPQ4+mfE0j89UOoGZMTRLTEq5O+vs9ayr3Yfm0o3NRa6/dvbK7ZjCfw9GAzd47rTK74XPeN70iO+x0HPWjJMjkQFLkgFTkQOU1kEu1fC7lVNnldB7YEzK8R3O/hMXcaQwGVZ3Vt3mAafny17q9iwp5INeypYv7ui8bm8pqHxuJgIL/3SYwPlLj2Ovumx9EqJoXdKLIkxoX/Wrqahhs1lm9lUtqmx1G0t28rWiq2U1pYedGxKVAo94nrQI75HY6nrGd+THnE9SI9J12XZLkwFLkgFTkRaxFoo2wa7Vh5c7vasAd+By4bEpAaKXGr/4HO/wHNKP633eghrLUUVdY2Fbn+p27inkm0l1QcdmxAVRq/UGHqlxNArJTb4HHh0T4oK+WXDyurKKCwvpLC8kK3lW9lavpXCisDrHZU78O8fkANEeiPpFtuNrNgsusV2O7AdF9jOjMnU5dlOTAUuSAVORE6IrwH2FQQmGt67EYo3HHgu337wsbHpTcpdv8Bzcl9I6hW4JBvilw/bUlVdA5uLq9iyt4qte6sO2t66r+qgS7JejyE7KZpeKTH0TImhZ0o03ROj6ZYYRfekaDIToogIC92CV++rZ0fljkCpCxa8HZU72Fm5kx2VO9hTveewz6REpRxe8GK7kRGTQUZMBmnRaUR4O8cgk65GBS5IBU5E2k1dJezdBHs3BIvdhsDr4g1QsfPgY8NjIaknJPaAxJ7B7V6B10k9A5dsPZ1jqo4T5fNbdpbVsKU4UOi2BB+bgwVvb2XdQccbA2lxkXRPjCIrMYpuidF0Tzr4OSM+MmTP4tX56thVtaux0O1/3lG5g50VO9leuZ3qhurDPpcUmUR6TDoZ0YFClxGTceB1TFrj/nAtXecqKnBBKnAi4ojaisCZun0FULoVSrYGnvdvV+89+HhPGCR0D5S7/QUvPitQ7OKzIC4L4jK0TixQWdvAjtIadpRWs6Okhu1Nn0tr2FFSTWWd76DPeD2G9LhIMhIiyYiPJD0+kvT4KNLjD7zOiI8kLS6SqPDQKtLWWsrqythZuZPdVbvZU70n8Fy1hz3Ve9hTtYfd1bspri7GZ32HfT45Mpm0mDRSolJIiUohNSqV1OjUxteN+6NTiQ6LduAn7FpU4IJU4ETEleoqobTw8GJXWhjYLtsOh/1jayA27UChi2/66Bbclxm4lNuFlxuz1lJW03BYwdtZVsPu8lr2BB/FlbU0909iYnT4YcUuNS6SlNgIUmIiSImLIDU2gpTYCOIiw0JmZK3P72Nf7b7GYte07BVXF7O3Zm/jo7K+stnviA6Lbix5KVEpJEclkxSZRGJkIkmRSSRFJQWeg/sSIxMJ9+j/6WgNFbggFTgRCUl+H1TuCcxjV74r+LwzcGm2fOeB/ZW7D6xI0VRkQmDARWxaoNAdtJ0W3E47sN0FC1+Dz8/eyrrGUre7vCb4XHvIcw019c38HQMRXg/JseGkxEaSEnzeX+6SYyNIjgknKTqCxOhwkmLCSYwJJy4iDI/H3aWvpqGGfTX7KK4JFLtDC17jo3ovJbUl1PnrjvhdceFxBwpek7KXGJlIfEQ88RHxJEQkHPYcEx7TJUfjqsAFqcCJSKfmawgUvabFrrIYqooC06VU7oGq4sB2VRH4G5r/nsiEwECL6OTAVCmN24c8og55L7zzj4a01lJV52NvZV3jo7iyjn3B572VtYe913T6lEN5TOAsX2J0OIkxESQFt5Niwhv3J0SFEx8VRnzj84Ftt13itdZS3VBNaW0pJbUllNSWHLTd9FFaU9r4fnl9+VG/12M8xIXHNRa6/eUuPiKeuIg4YsNjiQs/5PmQ/bHhsSG3tu2RClxo/RQiInJ03jBI6BZ4HIu1UFMSKHiVe5qUvGC5q9534FG67cB2M/dNNQqLhqhEiIwPTKMSmXDguel24779xwU/ExEbeLj4MqQxhtjIMGIjw+iZ0rLl2Ooa/OyrqqOkqp7S6npKquoora5vfJRU1VNSfeC9zcWVlFTXU1Zdj/8Y51kivJ7DSt3+7bjIMGIjvYG8EYHMcZFeYhq3w4iJ8AaPC2uT0bvGGGLCY4gJj6FbXAv+Owzy+X1U1FdQXldOWV0Z5XXlB22X1ZVRVltGeX1wf20ZRaVFlNeVU1FfQVVDVYv+nOiw6INKXUx4DDFhwUd4DNFh0Qf2HeF5/zFp0WnH+9d0wlTgRES6KmMOnDlLG9Cyz1gLteUHylxNycFFr3of1JRCTVlgIuTa8sD9e7VlgX1HuJfqkGAHilxEXOC5abmLiAs8IuMOvA6PDj5ig88xgeeImAPb4TGODfqICPOQmRBFZkLrzlD6/Zby2gbKquspr2mgvKaeitqGxu2ymgPbTZ8Liqoaj62s8+E7Vgvcn9PrISbSS2xEGNERXqLDvURHeIkJPqLC92+HNdk+cFx0eOCYqHBP8Dn4CDvw2nuEy8Vej7fxPrnj4fP7qGqoorK+koq6CirqKwLb9RVU1VdRUR/cV1fZ+F55fTnV9dXsrtpNVUMVVfVVjc/7V9E4kozoDOZcOue4srYFFTgREWk5YwJnzKISILl36z/vazhQ7PaXuv3PdeWBwRy1FYHnuorgozLwqNh9YHv/e83d73c0nvAmhS4awqICl33DogP3/YUHn8Oig/ubPJoeFxYFYRHgjWzy3Ny+SPBGHHhu5ZlFj8c0XkY9XtZaahv8VNY2UFnro7KugcraBipqG6iq8wVK3iHblbU+aup9VNU1NF4u3rbPR1Wdj+rg/iPdB3gs4V5DVJiXyCZFLzrcS2SYh4gwD5FhHiLDvI3bEYe8jgz3EOH1EBnuJdIbeB3uDewLD/MQ4Q0nIiyFCG86qWGGbhGB9yPDgscFn8O95ogDTqy11PhqDip01Q3VVDVUUV0feHZ6sErIFzhjzHTgr4AX+Ie19vcORxIRkSPxhkFMSuBxoqyFhppAoauvgvrq4HZ18LF/uwrqqg5sNz6CxzXUBr+nInD5uCG4r/G96iPfK9ha3ojgI7yZ7Sb7PGGHvB8eKJ/esOBzc6/Dmuw/8Np4wojyhBHl8ZLqCQvs94RBmBciwiFh/z7vgfc8EcHXXjDeg983HvCE4cdDjd9Q1QDV9bax3NU0PvzUNhzYbnxusq+23hd8HTi2rsFPRW0DtfV+6nyB92sb/NQ1+APPvuMrjUcS4fUQ5jWNhS58/2uP58B28L0wj4fwsDDCPQmEe5NIi4/g/H5tGqdVQrrAGWO8wEPA2UAhsMAY84a1dqWzyUREpN0Zc+BMWnvzNQRKXkNNoNj56gLlzlcLDXWB/YfuO+i5JrhdB7568Ncf2PbVNb9dX9Vk//7PNASf6wOlcv/+1p6JbAMeICb4ANOk7O1/9gSejeeQ95p5vf94jwcivRDlafLegYc1XvyY4MODH4Mv+Oy3HnwY/BZ8NnBM0+fAI/CeD0OD/5B9TZ4b9r/2G3w+Q4Pf4rPBz2Bp8JvAPZyM6PC/9/1CusABJwPrrbUbAYwxLwAzARU4ERFpO94w8Abvu3Mjv79JsWtS9PwNwYevyXZD8P2GI7xfH3ztCwxYafq+9R/y2nfg2Kavrf/A55tu+/1Njmn6XpNn6z/wnrXB/fVgfRjrx2v9eP2+wHv7vwvb5LPBz9lD9x36CH7e2kM+32T7aPfBRWQD93TM/32bEeoFLhvY2uR1ITD+0IOMMdcB1wH06tWrY5KJiIh0FI8HPJFdcv6+dtdY6A4pdg5PwxbqM+I1dwfhYX+j1trHrLW51trc9PT0DoglIiIinYIJXh72hgUGp4RHHRjh7KBQL3CFQM8mr3sA2x3KIiIiItIhQr3ALQAGGmP6GmMigMuBNxzOJCIiItKuQvoeOGttgzHmRuB9AtOIPGGtzXc4loiIiEi7CukCB2CtfQd4x+kcIiIiIh0l1C+hioiIiHQ5KnAiIiIiIUYFTkRERCTEqMCJiIiIhBgVOBEREZEQowInIiIiEmJU4ERERERCjAqciIiISIhRgRMREREJMSpwIiIiIiFGBU5EREQkxKjAiYiIiIQYFTgRERGREKMCJyIiIhJiVOBEREREQowKnIiIiEiIUYETERERCTEqcCIiIiIhxlhrnc7QoYwxpcC6YxyWBhR1QBw3SARKnQ5Bx+Royz/jRL/reD7fms+09NiWHKffB2d0pd+J9v59aM3x+p04WFf6nXDL78Ohn+9trU0/7AhrbZd6AI+14Jg8p3O66e+js+Royz/jRL/reD7fms+09Fj9PrTffyOhkMUtvxPt/fvQmuP1O9F+/424PYtbfh9a+vmueAn1TacDuIxb/j46Ikdb/hkn+l3H8/nWfKalx7rl//5u4aa/j670O9Hevw+tOd5N/w24gZv+Pto7i1t+H1r0+S53CbUljDF51tpcp3OIuIF+H0QOpt8JcYOueAauJR5zOoCIi+j3QeRg+p0Qx+kMnIiIiEiI0Rk4ERERkRCjAiciIiISYlTgREREREKMCpyIiIhIiFGBOwZjTKwx5iljzN+NMVc4nUfEacaYfsaYx40xLzudRcQNjDEXBv+NeN0Yc47TeaRr6JIFzhjzhDFmtzFmxSH7pxtj1hhj1htjfhXcfTHwsrX2B8AFHR5WpAO05nfCWrvRWnutM0lFOkYrfydeC/4bcTVwmQNxpQvqkgUO+CcwvekOY4wXeAiYAQwDvm2MGQb0ALYGD/N1YEaRjvRPWv47IdIV/JPW/07cHnxfpN11yQJnrZ0L7D1k98nA+uDZhTrgBWAmUEigxEEX/fuSzq+VvxMinV5rfidMwH3Au9baRR2dVbomFZIDsjlwpg0CxS0b+DfwTWPMI7hrTTiR9tbs74QxJtUY8ygw2hhzqzPRRBxxpH8nfgScBXzLGPNDJ4JJ1xPmdAAXMc3ss9baSuCajg4j4gJH+p0oBvSPlHRFR/qdeAB4oKPDSNemM3AHFAI9m7zuAWx3KIuIG+h3QuRg+p0Q11CBO2ABMNAY09cYEwFcDrzhcCYRJ+l3QuRg+p0Q1+iSBc4Y8zwwHxhsjCk0xlxrrW0AbgTeB1YBL1pr853MKdJR9DshcjD9TojbGWut0xlEREREpBW65Bk4ERERkVCmAiciIiISYlTgREREREKMCpyIiIhIiFGBExEREQkxKnAiIiIiIUYFTkS6NGOMNcY80+R1mDFmjzHmrVZ+z6fGmNzg9jvGmKQ2jioi0khroYpIV1cJ5Bhjoq211cDZwLYT+UJr7bltkkxE5Ah0Bk5EBN4Fzgtufxt4fv8bxphYY8wTxpgFxpjFxpiZwf3RxpgXjDHLjDH/AqKbfKbAGJMW3H7NGLPQGJNvjLmuyTEVxph7jTFLjTFfGmMyO+IHFZHOQQVORAReAC43xkQBI4Gvmrx3G/CxtXYcMAX4ozEmFrgeqLLWjgTuBcYe4bu/Z60dC+QCPzbGpAb3xwJfWmtPAuYCP2jrH0pEOi8VOBHp8qy1y4A+BM6+vXPI2+cAvzLGLAE+BaKAXsBpwLNNPr/sCF//Y2PMUuBLoCcwMLi/Dth/n93C4J8vItIiugdORCTgDeB+4Awgtcl+A3zTWrum6cHGGICjLiZtjDkDOAs4xVpbZYz5lEABBKi3Bxaj9qH/PRaRVtAZOBGRgCeA31prlx+y/33gRybY2Iwxo4P75wJXBPflELj0eqhEYF+wvA0BJrRLchHpclTgREQAa22htfavzbx1NxAOLDPGrAi+BngEiDPGLAN+AXzdzGffA8KCx9xN4DKqiMgJMwfO4IuIiIhIKNAZOBEREZEQowInIiIiEmJU4ERERERCjAqciIiISIhRgRMREREJMSpwIiIiIiFGBU5EREQkxPx/JZud8qiroDQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "apparent_densities = get_assemblage_apparent_density(size_classes = size_classes.Median, bin_boundaries = bin_boundaries,\n",
    "                                rho_p=rho_p, d_p=d_p, d_p_li=d_p_li,d_p_ui=d_p_ui,\n",
    "                                F=F, F_li=F_li, F_ui=F_ui,plot=True)\n",
    "plt.savefig(filepath_apparent_densities_plot, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "074e8ee4-0b66-420d-9c6d-5c147eccb9dc",
   "metadata": {},
   "source": [
    "# Merge the average apparent densities for each grain size population on the variable schema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "14654d46-0527-44fe-8a1e-666bb2214713",
   "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>rho_a</th>\n",
       "      <th>rho_a_li</th>\n",
       "      <th>rho_a_ui</th>\n",
       "      <th>SPMC_mean</th>\n",
       "      <th>SPMC_ci_lower</th>\n",
       "      <th>SPMC_ci_upper</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Assemblage</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(0.0, 3.0]</th>\n",
       "      <td>2328.658630</td>\n",
       "      <td>1000.332554</td>\n",
       "      <td>3174.833405</td>\n",
       "      <td>TVC_(0.0,_3.0]_mean</td>\n",
       "      <td>TVC_(0.0,_3.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(0.0,_3.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(3.0, 20.0]</th>\n",
       "      <td>757.798133</td>\n",
       "      <td>318.871040</td>\n",
       "      <td>1213.433144</td>\n",
       "      <td>TVC_(3.0,_20.0]_mean</td>\n",
       "      <td>TVC_(3.0,_20.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(3.0,_20.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(20.0, 200.0]</th>\n",
       "      <td>93.086999</td>\n",
       "      <td>31.924500</td>\n",
       "      <td>183.024290</td>\n",
       "      <td>TVC_(20.0,_200.0]_mean</td>\n",
       "      <td>TVC_(20.0,_200.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(20.0,_200.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(200.0, 459.0]</th>\n",
       "      <td>15.866572</td>\n",
       "      <td>4.477126</td>\n",
       "      <td>37.806383</td>\n",
       "      <td>TVC_(200.0,_420.0]_mean</td>\n",
       "      <td>TVC_(200.0,_420.0]_obs_ci_lower</td>\n",
       "      <td>TVC_(200.0,_420.0]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(0.0, 459.0]</th>\n",
       "      <td>614.428377</td>\n",
       "      <td>258.262330</td>\n",
       "      <td>921.853798</td>\n",
       "      <td>Total_Volume_Concentration_[ppm]_mean</td>\n",
       "      <td>Total_Volume_Concentration_[ppm]_obs_ci_lower</td>\n",
       "      <td>Total_Volume_Concentration_[ppm]_obs_ci_upper</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      rho_a     rho_a_li     rho_a_ui  \\\n",
       "Assemblage                                              \n",
       "(0.0, 3.0]      2328.658630  1000.332554  3174.833405   \n",
       "(3.0, 20.0]      757.798133   318.871040  1213.433144   \n",
       "(20.0, 200.0]     93.086999    31.924500   183.024290   \n",
       "(200.0, 459.0]    15.866572     4.477126    37.806383   \n",
       "(0.0, 459.0]     614.428377   258.262330   921.853798   \n",
       "\n",
       "                                            SPMC_mean  \\\n",
       "Assemblage                                              \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_mean   \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_mean   \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_mean   \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_mean   \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_mean   \n",
       "\n",
       "                                                SPMC_ci_lower  \\\n",
       "Assemblage                                                      \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_obs_ci_lower   \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_obs_ci_lower   \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_obs_ci_lower   \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_obs_ci_lower   \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_obs_ci_lower   \n",
       "\n",
       "                                                SPMC_ci_upper  \n",
       "Assemblage                                                     \n",
       "(0.0, 3.0]                        TVC_(0.0,_3.0]_obs_ci_upper  \n",
       "(3.0, 20.0]                      TVC_(3.0,_20.0]_obs_ci_upper  \n",
       "(20.0, 200.0]                  TVC_(20.0,_200.0]_obs_ci_upper  \n",
       "(200.0, 459.0]                TVC_(200.0,_420.0]_obs_ci_upper  \n",
       "(0.0, 459.0]    Total_Volume_Concentration_[ppm]_obs_ci_upper  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "apparent_densities = apparent_densities.merge(variable_schema, left_index=True, right_index=True, how = 'outer')\n",
    "apparent_densities"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bbfe098-b990-4247-bc6c-659ee91ccfcf",
   "metadata": {},
   "source": [
    "# Calculate the SPMC for the different grain size populations\n",
    "This is done by multiplying apparent density (in g/cm3) with TVC (in ul/l). Dividing by 1000 gives the SPMC in kg/m3.\n",
    "\n",
    "The SPMC grids are then written to a netCDF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a213ef85-a149-4b22-903a-d697621b6857",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 4min 11s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for raster_file in raster_files:\n",
    "    # read the raster and it's data variables\n",
    "    raster = xr.open_dataset(raster_file)\n",
    "    data_vars = list(raster.data_vars)\n",
    "    data_vars.remove('spatial_ref')\n",
    "    for var in data_vars:\n",
    "        # if the var matches a SPMC_mean colun, multiply with rho_a\n",
    "        match = apparent_densities[apparent_densities.SPMC_mean == var]\n",
    "        if len(match) != 0:\n",
    "            assemblage = match.index.str.replace(' ','_')[0]\n",
    "            raster[f'SPMC_{assemblage}_mean'] = raster[var] *match['rho_a'].values[0]/1000\n",
    "        # if the var matches a SPMC_ci_lower column, multiply with rho_a_li\n",
    "        match = apparent_densities[apparent_densities.SPMC_ci_lower == var]\n",
    "        if len(match) != 0:\n",
    "            assemblage = match.index.str.replace(' ','_')[0]\n",
    "            raster[f'SPMC_{assemblage}_ci_lower'] = raster[var] * match['rho_a_li'].values[0]/1000\n",
    "        # if the var matches a SPMC_ci_upper column, multiply with rho_a_ui\n",
    "        match = apparent_densities[apparent_densities.SPMC_ci_upper == var]\n",
    "        if len(match) != 0:\n",
    "            assemblage = match.index.str.replace(' ','_')[0]\n",
    "            raster[f'SPMC_{assemblage}_ci_upper'] = raster[var] *match['rho_a_ui'].values[0]/1000\n",
    "        # remove the original variable\n",
    "        raster = raster.drop_vars(var)\n",
    "    # write the resulting raster to a file\n",
    "    file_path = os.path.dirname(raster_file)\n",
    "    file_name = f'{os.path.split(file_path)[-1]}_SPMC_{assemblage}.nc'\n",
    "    raster.to_netcdf(os.path.join(file_path, file_name))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pcp-env2",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "bfbba0520bba5db8d0b67600201a76128af681101e25bae349f2447023d45a53"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}