CMIP6_diff_osa_Taylor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import pandas as pd\n",
    "import math\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates\n",
    "from matplotlib.patches import Ellipse\n",
    "import matplotlib.transforms as transforms\n",
    "import seaborn as sns\n",
    "from datetime import datetime\n",
    "import calendar\n",
    "from typing import List, Tuple, Dict\n",
    "import warnings\n",
    "from matplotlib import ticker, cm\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import geopandas as gp\n",
    "import regionmask\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('default')\n",
    "import cartopy.crs as ccrs\n",
    "import os\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.facecolor'] = 'white'\n",
    "from collections import namedtuple\n",
    "\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
    "import matplotlib.colors as colors\n",
    "\n",
    "from statsmodels.tools.sm_exceptions import ConvergenceWarning\n",
    "warnings.simplefilter('ignore', ConvergenceWarning)\n",
    "\n",
    "import geopandas as gpd\n",
    "import rioxarray as rio\n",
    "from shapely.geometry import mapping\n",
    "from texttable import Texttable\n",
    "import gcsfs\n",
    "import geopandas as gpd\n",
    "from CMIP6_IO import CMIP6_IO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def list_all_model_combinations(prefix, period):\n",
    "    io = CMIP6_IO()\n",
    "    files_on_gcs = []\n",
    "    counter=0\n",
    "    for blob in io.bucket.list_blobs(prefix=prefix):\n",
    "        if period in blob.name and \"srf\" not in blob.name:\n",
    "            counter+=1\n",
    "            files_on_gcs.append(blob.name) \n",
    "        \n",
    "    # Variable names\n",
    "    var_names = np.unique([file.split('/')[-1].split('_')[0] for file in files_on_gcs]).tolist()\n",
    "    var_names = list(map(lambda x: x.replace(\"sw\",\"sw_srf\"), var_names))\n",
    " \n",
    "    print(var_names)\n",
    "        \n",
    "    # Model names\n",
    "    model_names = np.unique([file.split('/')[-1].split('_')[1] for file in files_on_gcs])\n",
    "    # Factor names\n",
    "    factor_names = [\"snow_sensitivity\", \"no_chl\", \"no_wind\", \"no_osa\", \"no_meltpond\", \"no_clouds\"]\n",
    "    #factor_names = [\"osa\", \"no_chl\", \"no_wind\", \"no_osa\", \"no_meltpond\", \"snow_sensitivity\", \"no_ice\"]\n",
    "    # Member IDS\n",
    "    file_names = np.unique([file.split('/')[-1].split('_')[2] for file in files_on_gcs])\n",
    "    file_names = [file for file in file_names if file.startswith('r')]\n",
    "    member_names = np.unique(file_names)\n",
    "\n",
    "    s=\"_\"\n",
    "    combined_names=[]\n",
    "    for var_name in var_names:\n",
    "        for model_name in model_names:\n",
    "            for member_name in member_names:\n",
    "                combined_names.append(s.join([var_name,model_name,member_name]))\n",
    "\n",
    "    print(f\"Total number of files in bucket {counter} and unique variables {len(var_names)} and models {len(model_names)} and combined {len(combined_names)}\")\n",
    "    combos = namedtuple('combos', ['combined_names', 'factor_names', 'model_names'])\n",
    "\n",
    "    return combos(np.unique(combined_names), factor_names, model_names), files_on_gcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_LME_records():\n",
    "    lme_file = \"gs://actea-shared/Shapefiles/LME66/LMEs66.shp\"\n",
    "    return gpd.read_file(lme_file)\n",
    "\n",
    "def convert_to_180(ds):\n",
    "    return (ds.assign_coords(lon=(((ds.lon + 180) % 360) - 180))).sortby('lon')\n",
    "\n",
    "def get_data_within_LME(ds, LME):\n",
    "    # Extract the polygon defining the boundaries of the LME\n",
    "    shdf = get_LME_records()\n",
    "    shdf_sel = shdf[shdf[\"LME_NAME\"] == LME]\n",
    "    \n",
    "   # ds = ds.rio.write_crs(4326)\n",
    "    ds.rio.set_spatial_dims(y_dim=\"lat\", x_dim=\"lon\", inplace=True)\n",
    "    ds.rio.write_crs(4326,  inplace=True)\n",
    "    \n",
    "    clipped = ds.rio.clip(geometries=shdf_sel.geometry.apply(mapping), crs=ds.rio.crs)\n",
    "    if \"x\" in ds.dims:\n",
    "        clipped = clipped.rename({\"x\": \"lon\", \"y\": \"lat\"})\n",
    "   \n",
    "    return clipped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wrap it into a simple function\n",
    "def season_mean(ds, calendar=\"standard\"):\n",
    "    # Make a DataArray with the number of days in each month, size = len(time)\n",
    "    month_length = ds.time.dt.days_in_month\n",
    "\n",
    "    # Calculate the weights by grouping by 'time.season'\n",
    "    weights = (\n",
    "        month_length.groupby(\"time.season\") / month_length.groupby(\"time.season\").sum()\n",
    "    )\n",
    "\n",
    "    # Test that the sum of the weights for each season is 1.0\n",
    "    np.testing.assert_allclose(weights.groupby(\"time.season\").sum().values, np.ones(4))\n",
    "\n",
    "    # Calculate the weighted average\n",
    "    return (ds * weights).groupby(\"time.season\").sum(dim=\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_seasonal_maps(ds1_seasonal, ds2_seasonal, ds_diff, var_name):\n",
    "    \n",
    "    land_10m = cfeature.NaturalEarthFeature('physical', 'land', '10m',\n",
    "                                        edgecolor='face',\n",
    "                                        facecolor=\"lightgray\")\n",
    "    \n",
    "    # Quick plot to show the results\n",
    "    notnull = pd.notnull(ds1_seasonal[var_name][0])\n",
    "\n",
    "    fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(14, 12), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
    "    for i, season in enumerate((\"DJF\", \"MAM\", \"JJA\", \"SON\")):\n",
    "        ds1_seasonal[var_name].sel(season=season).where(notnull).plot.pcolormesh(\n",
    "            ax=axes[i, 0],\n",
    "            vmin=0.001,\n",
    "            vmax=1,\n",
    "            levels=[0.01,0.02,0.03,0.03,0.04,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.5,0.75,1.0],\n",
    "            cmap=\"RdBu_r\",\n",
    "            add_colorbar=True,\n",
    "            transform=ccrs.PlateCarree(),\n",
    "            extend=\"max\",\n",
    "           \n",
    "        )\n",
    "\n",
    "        ds2_seasonal[var_name].sel(season=season).where(notnull).plot.pcolormesh(\n",
    "            ax=axes[i, 1],\n",
    "            vmin=0.001,\n",
    "            vmax=1,\n",
    "            levels=[0.01,0.02,0.03,0.03,0.04,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.5,0.75,1.0],\n",
    "            cmap=\"RdBu_r\",\n",
    "            add_colorbar=True,\n",
    "            transform=ccrs.PlateCarree(),\n",
    "            extend=\"max\",\n",
    "        )\n",
    "\n",
    "        ds_diff[var_name].sel(season=season).where(notnull).plot.pcolormesh(\n",
    "            ax=axes[i, 2],\n",
    "            vmin=-0.01,\n",
    "            vmax=0.01,\n",
    "            cmap=\"RdBu_r\",\n",
    "                levels=[-0.008,-0.007,-0.006,-0.005,-0.004,-0.003,-0.002,-0.001,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008],\n",
    "            add_colorbar=True,\n",
    "            transform=ccrs.PlateCarree(),\n",
    "            extend=\"both\",\n",
    "        )\n",
    "\n",
    "        axes[i, 0].set_ylabel(season)\n",
    "        axes[i, 1].set_ylabel(\"\")\n",
    "        axes[i, 2].set_ylabel(\"\")\n",
    "        axes[i, 0].add_feature(land_10m)\n",
    "        axes[i, 1].add_feature(land_10m)\n",
    "        axes[i, 2].add_feature(land_10m)\n",
    "        axes[i, 0].coastlines()\n",
    "        axes[i, 1].coastlines()\n",
    "        axes[i, 2].coastlines()\n",
    "        \n",
    "    for ax in axes.flat:\n",
    "        ax.axes.get_xaxis().set_ticklabels([])\n",
    "        ax.axes.get_yaxis().set_ticklabels([])\n",
    "        ax.axes.axis(\"tight\")\n",
    "        ax.set_xlabel(\"\")\n",
    "\n",
    "    axes[0, 0].set_title(\"Taylor\")\n",
    "    axes[0, 1].set_title(f\"{var_name}\")\n",
    "    axes[0, 2].set_title(\"Difference\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "\n",
    "    fig.suptitle(f\"Seasonal {var_name}\", fontsize=16, y=1.02)\n",
    "    plt.show()\n",
    "    \n",
    "def create_distributions(ds1, ds2, name1, name2):\n",
    "    ds_avg1 = ds1.mean({\"lat\",\"lon\"})\n",
    "    ds_avg2 = ds2.mean({\"lat\",\"lon\"})\n",
    "    df1 = ds_avg1.to_dataframe()\n",
    "    df2 = ds_avg2.to_dataframe()\n",
    "    df1[\"Status\"]=name1\n",
    "    df2[\"Status\"]=name2\n",
    "    \n",
    "    df = pd.concat([df1,df2], axis=0)\n",
    "    print(df)\n",
    "    sns.kdeplot(df, x=\"osa\", hue=\"Status\", fill=False, alpha=0.5)\n",
    "    plt.title(f\"{name1} vs {name2}\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ghi', 'osa', 'par', 'uv', 'uva', 'uvb', 'uvi']\n",
      "Total number of files in bucket 896 and unique variables 7 and models 4 and combined 392\n",
      "[CMIP6_IO] Opening file actea-shared/light_sensitivity/ssp245/osa_CanESM5_r10i1p2f1_1979-01-01-1989-01-16_scenario_no_osa_ssp245.nc\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tk/.conda/envs/opendrift/lib/python3.11/site-packages/h5py/__init__.py:36: UserWarning: h5py is running against HDF5 1.14.3 when it was built against 1.14.2, this may cause problems\n",
      "  _warn((\"h5py is running against HDF5 {0} when it was built against {1}, \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[CMIP6_IO] Opening file actea-shared/light_sensitivity/ssp245/osa_CanESM5_r10i1p2f1_1979-01-01-1989-01-16_scenario_osa_ssp245.nc\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tk/.conda/envs/opendrift/lib/python3.11/site-packages/xarray/groupers.py:326: FutureWarning: 'M' is deprecated and will be removed in a future version, please use 'ME' instead.\n",
      "  self.index_grouper = pd.Grouper(\n",
      "/home/tk/.conda/envs/opendrift/lib/python3.11/site-packages/xarray/groupers.py:326: FutureWarning: 'M' is deprecated and will be removed in a future version, please use 'ME' instead.\n",
      "  self.index_grouper = pd.Grouper(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 osa  spatial_ref  Status\n",
      "time                                     \n",
      "1979-01-31  0.460825            0  Taylor\n",
      "1979-02-28  0.488180            0  Taylor\n",
      "1979-03-31  0.483477            0  Taylor\n",
      "1979-04-30  0.476071            0  Taylor\n",
      "1979-05-31  0.474527            0  Taylor\n",
      "...              ...          ...     ...\n",
      "1988-09-30  0.314324            0     OSA\n",
      "1988-10-31  0.369926            0     OSA\n",
      "1988-11-30  0.401416            0     OSA\n",
      "1988-12-31  0.448855            0     OSA\n",
      "1989-01-31  0.478211            0     OSA\n",
      "\n",
      "[242 rows x 3 columns]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2yElEQVR4nO3deXxcdb3/8deZfSbJZF/bpPtGW7rSUgoUBFpBoSx69aIColy4iIgVlV5+4n69ICp6xfWyiCggi8gqrdBCS8tS2kKhe5umbfZ9z6zn98dJAqELbZrkTJL38/GYB+acMzOfOabJO9/VME3TRERERGSIcNhdgIiIiEhfUrgRERGRIUXhRkRERIYUhRsREREZUhRuREREZEhRuBEREZEhReFGREREhhSFGxERERlSFG5ERERkSFG4ERlmDMM4psfq1av77D3POusszjrrrD57Pbu1trbyP//zP8yaNYvk5GSSkpKYOXMm//3f/01ra+thr7/99tuZMWMGwWCQlJQUxo0bx7/927/x8ssvH/Y9ampq8Hq9GIbBhg0b+vsjiQwpLrsLEJGBtX79+h5f//CHP2TVqlW89NJLPY6fdNJJA1nWoFFZWcm5557Lnj17uPHGG7njjjsAeOmll/jRj37EQw89xL/+9S9yc3MBiMViLF68mC1btvDNb36TefPmAbBr1y6efvpp1qxZw6JFiw55nz//+c+Ew2EA7rnnHubOnTtAn1Bk8DO0t5TI8HbVVVfx2GOP0dLS0m/v0dVq0xetQaZp0tHRgd/vP+HX6o0lS5bw0ksvsWrVKk4//fQe59auXcvZZ5/NOeecwz//+U8AVq1axcc+9jHuvfdevvjFLx7yevF4HIfj0Eb06dOnU1VVxahRo9i5cyfl5eW2fWaRwUbdUiJyiLvvvpszzzyTnJwckpKSmD59OnfccQeRSKT7mh/+8Ie4XC4OHDhwyPOvvvpqMjMz6ejoOOJ71NXVcf311zNixAg8Hg9jx47l1ltvJRQK9bjOMAxuuOEGfve73zFlyhS8Xi9/+tOfDvuaF198MaNGjSIejx9ybv78+cyePbv760cffZT58+eTmppKIBBg7NixXH311Ue9Lxs2bGDFihV86UtfOiTYAJx++ulcffXVvPDCC7z11lsA1NbWApCfn3/Y1zxcsHn99dd59913+cIXvsA111xDY2Mjjz/++FFrE5H3KdyIyCH27NnD5Zdfzp///GeeeeYZvvSlL/HTn/6Ua6+9tvuaa6+9FpfLxe9///sez62rq+Phhx/mS1/6Ej6f77Cv39HRwdlnn80DDzzAsmXLePbZZ/n85z/PHXfcwaWXXnrI9U8++SS//e1vue2223jhhRc444wzDvu6V199Nfv37z+ki2379u288cYb3S0n69ev5zOf+Qxjx47l4Ycf5tlnn+W2224jGo0e9b6sXLkSsELUkXSd67p27ty5uN1uvva1r/GXv/yF8vLyo74HWN1QXZ/ns5/9LIFAoPuYiBwDU0SGtSuvvNJMSko64vlYLGZGIhHzgQceMJ1Op1lXV9fjuTk5OWYoFOo+dvvtt5sOh8MsLi7uPrZo0SJz0aJF3V//7ne/MwHzb3/7W4/3uv32203AXLFiRfcxwExNTe3xvkcSiUTM3Nxc8/LLL+9x/Fvf+pbp8XjMmpoa0zRN88477zQBs6Gh4SNf84Ouu+46EzC3b99+xGu2bdtmAuZ//ud/dh+75557zOTkZBMwATM/P9+84oorzFdeeeWQ57e2tprBYNA89dRTu49deeWVpmEY5u7du4+rXpHhSi03InKITZs2cdFFF5GZmYnT6cTtdnPFFVcQi8XYuXNn93Vf+9rXqKqq4tFHHwWs8SO//e1v+cQnPsHo0aOP+PovvfQSSUlJfOpTn+px/KqrrgLgxRdf7HH8Yx/7GOnp6R9Zt8vl4vOf/zxPPPEEjY2NgDWg989//jNLly4lMzMTgFNOOQWAf/u3f+Nvf/sbpaWlH/nax8rsHMZoGEb3sauvvpqDBw/y17/+lRtvvJHCwkIefPBBFi1axE9/+tMez//b3/5GU1NTjy6yq6++GtM0ue+++/qsTpGhTOFGRHrYv38/Z5xxBqWlpfzyl79kzZo1vPnmm9x9990AtLe3d187a9YszjjjjO5zzzzzDPv27eOGG2446nvU1taSl5fXIwAA5OTk4HK5usepdDnSeJXDufrqq+no6ODhhx8G4IUXXqC8vLzHYN4zzzyTJ598kmg0yhVXXMHIkSOZNm0aDz300FFfu6ioCIDi4uIjXrNv3z4ACgsLexxPTU3l3//93/nlL3/J66+/zjvvvENubi633norDQ0N3dfdc889+Hw+Pv7xj9PQ0EBDQwMnn3wyo0eP5v777ycWix3zvRAZrhRuRKSHJ598ktbWVp544gk+//nPc/rppzN37lw8Hs9hr7/xxhtZv349Gzdu5Ne//jUTJ07kvPPOO+p7ZGZmUllZ2d3K0aWqqopoNEpWVlaP4x8OQUdz0kknMW/evO5Wjvvuu4+CggIWL17c47qlS5fy4osv0tjYyOrVqxk5ciSXX375IVPlP6jrcz355JNHvKbr3Efdg6lTp/LZz36WSCTS3Rq2c+dO1q5dS0dHB0VFRaSnp3c/9u3bR2lpKS+88MJH3QKRYU/hRkR66AoSXq+3+5hpmvzxj3887PWXXHIJRUVFfOMb3+Bf//oX119//UeGkXPOOYeWlpZDQsIDDzzQff5EfPGLX+T1119n7dq1PP3001x55ZU4nc7DXuv1elm0aBG33347YHXJHcncuXNZvHgx99xzD6+++uoh59euXcu9997Lxz/+cebMmQNYrVRd69V82Pbt2wEoKCgA3h9I/Mc//pFVq1b1eDz33HO43W7uvffeY7wLIsOYvUN+RMRuHx5QvG3bNtPj8ZhnnXWW+dxzz5lPPPGEed5555kTJkwwAXPVqlWHvEbXQOCkpKTDDtL98IDi9vZ28+STTzZTUlLMn//85+bKlSvN7373u6bb7TYvuOCCHs8FzK985SvH9ZkaGhpMv99vjhw50gTMHTt29Dj/ne98x/ziF79oPvjgg+bq1avNJ5980jz77LNNt9ttvvvuu0d97YqKCnPatGlmIBAwb7nlFnPlypXmypUrzeXLl5uBQMCcNm2aWVFR0X39o48+ahYUFJjf/OY3zb///e/mK6+8Yj722GPmZZddZgLmFVdcYZqmNRg6Ly/PnDJlyhHf+9JLLzXdbrdZVVV1XPdDZLhRuBEZ5g43W+rpp582Z8yYYfp8PnPEiBHmN7/5TfP5558/YrjZt2+fCZjXXXfdYd/jw+HGNE2ztrbWvO6668z8/HzT5XKZo0aNMpcvX252dHT0uK434cY0TfPyyy83AXPhwoWHnHvmmWfM888/3xwxYoTp8XjMnJwc84ILLjDXrFlzTK/d0tJi/vd//7c5c+ZMMxAImIFAwDz55JPNH/3oR2ZLS0uPaw8cOGD+v//3/8yFCxeaeXl5psvlMlNSUsz58+eb//u//2tGo1HTNE3zySefNAHzrrvuOuL7/vOf/zQB82c/+9lx3AmR4UcrFIvICfvf//1fbrzxRt59912mTp1qdzkiMswp3IhIr23atIni4mKuvfZaFi5ceNSBtiIiA0XhRkR6bfTo0VRUVHDGGWfw5z//mby8PLtLEhFRuBEREZGhRVPBRUREZEhRuBEREZEhReFGREREhhSX3QUMtHg8TllZGSkpKce1pLuIiIjYxzRNmpubKSgowOE4etvMsAs3ZWVlh2xoJyIiIoPDgQMHGDly5FGvGXbhJiUlBbBuTjAYtLkaERERORZNTU0UFhZ2/x4/mmEXbrq6ooLBoMKNiIjIIHMsQ0o0oFhERESGFIUbERERGVIUbkRERGRIGXZjbkRERHorFosRiUTsLmPI8ng8HznN+1go3IiIiHwE0zSpqKigoaHB7lKGNIfDwZgxY/B4PCf0Ogo3IiIiH6Er2OTk5BAIBLQIbD/oWmS3vLycoqKiE7rHCjciIiJHEYvFuoNNZmam3eUMadnZ2ZSVlRGNRnG73b1+HQ0oFhEROYquMTaBQMDmSoa+ru6oWCx2Qq+jcCMiInIM1BXV//rqHivciIiIyJCicCMiIiJDisKNiIhIgqiqquLaa6+lqKgIr9dLXl4eS5YsYf369YDVbfPkk08e9+uOHj2au+66q2+LTWCaLSUiIpIgLrvsMiKRCH/6058YO3YslZWVvPjii9TV1dld2qCilhsRETlm8bjJgbo29la3UNrQTmsoandJQ0ZDQwNr167l9ttv5+yzz2bUqFHMmzeP5cuX84lPfILRo0cDcMkll2AYRvfXe/bsYenSpeTm5pKcnMwpp5zCv/71r+7XPeussygpKeHrX/86hmF0D9r93ve+x8yZM3vUcNddd3W/LsDq1auZN28eSUlJpKWlsXDhQkpKSvrzNvQJtdyIiMhHag/H2FBSx7byJlpDndN0ox04iHPSiHTmjc8nNXBiq8oOd8nJySQnJ/Pkk09y6qmn4vV6e5x/8803ycnJ4b777uPjH/84TqcTgJaWFi644AJ+9KMf4fP5+NOf/sSFF17Ijh07KCoq4oknnmDGjBn8x3/8B9dcc80x1xONRrn44ou55ppreOihhwiHw7zxxhuDYtaYwo2IiBxVRWMHz7xTRnNHFOIx/O3lBEMVdLS30Bhx8e5B2PaWi3Nnjuekk08BZ+8XXxvOXC4X999/P9dccw2/+93vmD17NosWLeKzn/0sJ598MtnZ2QCkpaWRl5fX/bwZM2YwY8aM7q9/9KMf8fe//52nnnqKG264gYyMDJxOJykpKT2e91GamppobGzkk5/8JOPGjQNgypQpffRp+5e6pURE5Ii2lTfxtw0HaO6Iku5o50LnOq4Jvs7lOSVcPbqez4xupcgfIhaNsmLDNt554V5oOGB32YPWZZddRllZGU899RRLlixh9erVzJ49m/vvv/+Iz2ltbeVb3/oWJ510EmlpaSQnJ7N9+3b2799/QrVkZGRw1VVXsWTJEi688EJ++ctfUl5efkKvOVAUbkRE5LD21bSy4r1KYnGT8Z46PmusYLyzCqcvBSadD6d9lYLzvsql/34NM6eehOny8WKpmy1rnoT6fXaXP2j5fD7OO+88brvtNtatW8dVV13Fd7/73SNe/81vfpPHH3+cH//4x6xZs4bNmzczffp0wuHwUd/H4XBgmmaPYx/e8fy+++5j/fr1nHbaaTzyyCNMnDiR1157rfcfboAo3IiIyCGqmjt4dks5cdNkSlIzn4yvxmdEIXsinPIlKJgJHms7AsPp4qwF85l72jngT2N1ZYCaN5+AumJ7P8QQcdJJJ9Ha2gqA2+0+ZGuCNWvWcNVVV3HJJZcwffp08vLy2LdvX49rPB7PIc/Lzs6moqKiR8DZvHnzIe8/a9Ysli9fzrp165g2bRp//etf++aD9SOFGxER6aEjEuOpzWWEo3EKvW2cF34RAxPyZ8DUS8HtP+Q5hmFw+sRcRk+ZS9SfyT8rkoi+9xSEWmz4BINTbW0tH/vYx3jwwQd55513KC4u5tFHH+WOO+5g6dKlgLVezYsvvkhFRQX19fUAjB8/nieeeILNmzfz9ttvc/nllxOPx3u89ujRo3nllVcoLS2lpqYGsGZRVVdXc8cdd7Bnzx7uvvtunn/++e7nFBcXs3z5ctavX09JSQkrVqxg586dg2LcjcKNiIj08PLOapo7oqS5Y3zSXI3TjFktNhM/DkeZKWMYBounFeAvmEq1mcr6CgdsfxY+1PUhh5ecnMz8+fP5xS9+wZlnnsm0adP4zne+wzXXXMOvf/1rAH72s5+xcuVKCgsLmTVrFgC/+MUvSE9P57TTTuPCCy9kyZIlzJ49u8dr/+AHP2Dfvn2MGzeue2DylClT+M1vfsPdd9/NjBkzeOONN7j55pu7nxMIBNi+fTuXXXYZEydO5D/+4z+44YYbuPbaawfojvSeYX64w22Ia2pqIjU1lcbGRoLBoN3liIgklL3VLfxjcxkG8OnARkaEiyElF2ZdAc5jm2C7u6qFp9/chaN8M1cUVpM+9WMwcm7/Ft6POjo6KC4uZsyYMfh8PrvLGdKOdq+P5/e3Wm5ERASwuqNe3FYFwKxAlRVsnC6YsvSYgw3A+Jxkxo7IIZ4+hvV1ybB3NYRb+6lqkUMp3IiICABv7qujJRQlzRXhtLC1lxHjzoGkzON+rdPGZWEEC9gRzqKqzYSSdX1crciRKdyIiAiN7RE2728AYJFnO26ikDEGCmb16vWyU7xMzk+BjNG8WpsMZZugvb4PKxY5MoUbERFh/Z4aonGTQm8LY0I7wXDA+HOPOoD4oywYm4XDn84+M4/yNgcUv9KHFYscmcKNiMgwV9nUwbbyZjDjnGFusvLMyDmQlHVCr5sacFutN+mj2Nzgh8qt0FrbN0WLHIXCjYjIMPfaXitwTPHWkhuvshbnG3V6n7z2zMI08CSzyxxBS9QBpRv65HVFjkbhRkRkGKtq7mBvdSuGGWe++bZ1cNRCcPfNlOfcoI+CNB+xlAK2NPmhYgtEOvrktUWOROFGRGQYe7PYGuQ70VtLutkI3hTIn9mn7zGzMB18qWxpzyIWjUD52336+iIfpnAjIjJM1bWG2VVljbU5JdbVarPguNa0ORbjc5JJ9rppDYxkV4sPSt+CD20PINKXFG5ERIapDfvqME0Y66oh22gEXxDyZvT5+zgdBlNHBCE5m21tQehohLo9ff4+MrC+973vMXPmTLvLOCyFGxGRYag1FGV7RTOYJvN4zzpYdGqft9p0mZIXBMNJiWMkrVGHNfZG+pVhGEd9XHXVVXaX2G/657tYREQS2pbSRmJxk3xnA/lGnbXTd97J/fZ+6Uke8lN9lIez2dHiY3btboi0H3aHcekb5eXl3f/7kUce4bbbbmPHjh3dx/x+++99JBLB7Xb3+euq5UZEZJiJxU22HGwETGYaO62DI+aAs+9/yXzQ5PwgeJLZHs6CeAyqt/fr+w13eXl53Y/U1FQMw+j+2u12c9111zFy5EgCgQDTp0/noYce6n7uAw88QGZmJqFQqMdrXnbZZVxxxRWHfb94PM4PfvADRo4cidfrZebMmfzzn//sPr9v3z4Mw+Bvf/sbZ511Fj6fjwcffLBfPrvCjYjIMLOrqpmWUJTkeAsTOAgOF4yY3e/vOzE3GYdhUOnMpy7shIp3+/09+4tpmoSj8QF/mKbZJ/V3dHQwZ84cnnnmGd59913+4z/+gy984Qu8/vrrAHz6058mFovx1FNPdT+npqaGZ555hi9+8YuHfc1f/vKX/OxnP+POO+/knXfeYcmSJVx00UXs2rWrx3Xf/va3ufHGG9m2bRtLlizpk8/zYeqWEhEZZrr2kJruLMFpAHnTwZPU7+8b8LgYlRmgOJrN9hY/p3kOQlsdBDL6/b37WiRmcveq3QP+vl85ezweV++3xOgyYsQIbr755u6vv/rVr/LPf/6TRx99lPnz5+P3+7n88su57777+PSnPw3AX/7yF0aOHMlZZ5112Ne88847+fa3v81nP/tZAG6//XZWrVrFXXfdxd1339193U033cSll156wp/haNRyIyIyjFQ1dVDe2IEz1sF0On85jzxlwN5/cn4KOL3sihd0FrR1wN5b3heLxfjxj3/MySefTGZmJsnJyaxYsYL9+/d3X3PNNdewYsUKSktLAbjvvvu46qqrMA6z31hTUxNlZWUsXLiwx/GFCxeybdu2Hsfmzp3bD5+oJ7XciIgMI1tKGwEY76wkyRmzdv5Oyhyw9x+dmYTTYVDnyqEuvIeM6h0wum+2ehhIbqfBV84eb8v79oWf/exn/OIXv+Cuu+5i+vTpJCUlcdNNNxEOh7uvmTVrFjNmzOCBBx5gyZIlbNmyhaeffvqor/vh4GOa5iHHkpL6v5VQ4UZEZJgIR+Od07/jTKdz1kxB/4+1+SCf28nIdD8lkQz2tPnJ8FRBez340we0jhNlGEafdA/ZZc2aNSxdupTPf/7zgDUYeNeuXUyZMqXHdV/+8pf5xS9+QWlpKeeeey6FhYWHfb1gMEhBQQFr167lzDPP7D6+bt065s2b138f5Ahs7Zb6yU9+wimnnEJKSgo5OTlcfPHFPaapHcnLL7/MnDlz8Pl8jB07lt/97ncDUK2IyOC2s7KZcDROWqyWka5ma9G+zIFvfRiXnQxON3vjedaBml1Hf4L0ufHjx7Ny5UrWrVvHtm3buPbaa6moqDjkus997nOUlpbyxz/+kauvvvqor/nNb36T22+/nUceeYQdO3Zwyy23sHnzZr72ta/118c4IlvDzcsvv8xXvvIVXnvtNVauXEk0GmXx4sW0trYe8TnFxcVccMEFnHHGGWzatIn/+q//4sYbb+Txxx8fwMpFRAafdzu7pKY7ijEMoGAWOAb+18DYbKtbotzIthb0q9k54DUMd9/5zneYPXs2S5Ys4ayzziIvL4+LL774kOuCwSCXXXYZycnJhz3/QTfeeCPf+MY3+MY3vsH06dP55z//yVNPPcWECRP650MchWH21byyPlBdXU1OTg4vv/xyj2atD/r2t7/NU0891WOA0nXXXcfbb7/N+vXrP/I9mpqaSE1NpbGxkWAw2Ge1i4gksurmEA++VoIj0sKX3StJ8hiw4CsDMkvqcB56Yz8VtQ2cG3qJ6WkdcNpXbavlo3R0dFBcXMyYMWPw+fpmt/TB5LzzzmPKlCn86le/6vf3Otq9Pp7f3wk1W6qx0fqrIiPjyNMC169fz+LFi3scW7JkCRs2bCASiRxyfSgUoqmpqcdDRGS4ea/M+vk61lFJkisOWRNtDRPjspPB5WNPPAdMU11TCaiuro6HH36Yl156ia985St2l3NcEibcmKbJsmXLOP3005k2bdoRr6uoqCA3N7fHsdzcXKLRKDU1NYdc/5Of/ITU1NTux5EGQ4mIDFWxuMmOimaIx5hqdoaI/L7fIPN4dHVN7TdzicRRuElAs2fP5tprr+X2229n0qRJdpdzXBJmttQNN9zAO++8w9q1az/y2sNNNTvccYDly5ezbNmy7q+bmpoUcERkWNlX20pbOEZSpJbRgVbwp0H6aFtrykzykOJz0RzO4GC7hzEN+yAW7beNO+X47du3z+4Sei0hvou++tWv8tRTT/HKK68wcuTIo16bl5d3yIjuqqoqXC4XmZmHrtXg9Xrxer19Wq+IyGCytczqjp/sOIDDwGq1OcwfgwPJMAxGZyaxpSNCSWsqY2LV0LgfMsbaWpcMDbZ2S5mmyQ033MATTzzBSy+9xJgxYz7yOQsWLGDlypU9jq1YsYK5c+f2y86iIiKDWXs4RnFNK0TamOI4CIbD2m4hAYzKDAAGJWa+daBur631yNBha7j5yle+woMPPshf//pXUlJSqKiooKKigvb29u5rli9f3mMH0uuuu46SkhKWLVvGtm3buPfee7nnnnt67JEhIiKW7RVNxOImOfEqsr1RyBwH3hS7ywKgMCOAYUCdkUZTxAF1xXaXJEOEreHmt7/9LY2NjZx11lnk5+d3Px555JHua8rLy3vsdTFmzBiee+45Vq9ezcyZM/nhD3/Ir371Ky677DI7PoKISEKzViQ2OcnYZx3IPfKEjYHmczvJT/WBP4397V5orYGORrvLkiHA1jE3x7LEzv3333/IsUWLFrFx48Z+qEhEZOhoaAtT0diBEWpgoq8e3D5bViQ+mqKMJMoaOthHHtMotlpvCmbaXZYMcgkzFVxERPrWjopmAIroXNsm56SEm41kjbuB/fEs4iYadyN9QuFGRGQIMk2THZXNEI8yKQG7pLrkBX143Q5C7jQqQ26o3wfxuN1lySCncCMiMgRVt4SobQnj6qhlnK8NApkQLLC7rEM4HAaF6QHwpHAwnATRELQcuoGj9N6BAwf40pe+REFBAR6Ph1GjRvG1r32N2tra7mv27t3Lv//7v1NQUIDP52PkyJEsXbqUnTsP3ffrr3/9K06nk+uuu24gP8ZxUbgRERmCurqkRlOOz2lC3jTb17Y5kpHpfjAMDhqdu4Q37D/6E+SY7d27l7lz57Jz504eeughdu/eze9+9ztefPFFFixYQF1dHeFwmPPOO4+mpiaeeOIJduzYwSOPPMK0adO6t0X6oHvvvZdvfetbPPzww7S1tdnwqT5aYnW+iojICTPNzu0Wou1Mdhy0Qk3uVLvLOqKR6da4m7JYOjETnPUlUHSqzVUNDV/5ylfweDysWLECv98PQFFREbNmzWLcuHHceuutXHvttezdu5eXXnqJUaNGATBq1CgWLlx4yOvt27ePdevW8fjjj7Nq1Soee+yxHsu1JAq13IiIDDFljR00d0TxtFczOhCCtFHgS7W7rCPKSvbgczsJe4JUhdzQeADiMbvLOjrThGh44B/HMMu4S11dHS+88ALXX399d7DpkpeXx+c+9zkeeeQRsrOzcTgcPPbYY8RiR7/v9957L5/4xCdITU3l85//PPfcc0+vbl9/U8uNiMgQs6OiCTAZTyluB1aXVAIzDIOR6X52R6IcbE8mP1YPzeWQevTteGwVi8Canw38+57xDXB5junSXbt2YZomU6ZMOez5KVOmUF9fj9vt5le/+hXf+ta3+P73v8/cuXM5++yz+dznPsfYse9vhxGPx7n//vv53//9XwA++9nPsmzZMnbv3s348Ym1xIBabkREhpBY3GRnZQuEmpnkqQGnG7ISf0fnkel+wOAgudYBjbvpdx/cdPorX/kKFRUVPPjggyxYsIBHH32UqVOn9tjuaMWKFbS2tnL++ecDkJWVxeLFi7n33nttqf9o1HIjIjKEHKhroz0cI9BRRVFGGLKnHfNf+nbqHncTz7DG3TTsh1Gn2VzVUTjdViuKHe97jMaPH49hGGzdupWLL774kPPbt28nPT2drKwsAFJSUrjooou46KKL+NGPfsSSJUv40Y9+xHnnnQdYXVJ1dXUEAoHu14jH42zatIkf/vCHOJ3OE/tsfUgtNyIiQ4i13UKMiUbnDuAJsknmR+ked+MeJONuDMMKjQP9OI4Zb5mZmZx33nn85je/6bFnI0BFRQV/+ctf+MxnPoNxmNc0DIPJkyfT2toKQG1tLf/4xz94+OGH2bx5c49HS0sLzz///Indzz6mcCMiMkREYnH2VLdAWx2T/M3gC1qDiQeBrnE3eAIcjKRALApNZXaXNej9+te/JhQKsWTJEl555RUOHDjAP//5T8477zxGjBjBj3/8YzZv3szSpUt57LHH2Lp1K7t37+aee+7h3nvvZenSpQD8+c9/JjMzk09/+tNMmzat+3HyySfzyU9+MuEGFqtbSkRkiCiuaSUcjRMMV5IfjEDu3IRd2+ZwCtL87K5qoczIAeqg8SCkFdpd1qA2YcIENmzYwPe+9z0+85nPUFtbS15eHhdffDHf/e53ycjIIB6PM3r0aL7//e+zb98+DMPo/vrrX/86YHVJXXLJJTgch7aJXHbZZXzmM5+hsrKS3Nzcgf6Ih6VwIyIyROyoaIZYmEmOMivT5A6OLqkuI9Ks6crl8XRME4ymUpsrGhpGjRrFfffdd8TzWVlZ/PKXvzzqa7zzzjtHPHfppZcSiUR6XV9/ULeUiMgQ0BGJUVzTCq3VTEput7ZaSMq0u6zjkp3ixe00aHemUB9xWuNujmNdF5EuCjciIkPA7qoWYnGTrEg5WZ5owq9tczhOh0FO0AeeJMrCfoh0QFvtRz9R5EMUbkREhoAdFc0QbmGiuwbD6YSck+wuqVdGpPnBcHSOu8FqvRE5Tgo3IiKDXGsoyoH6NmipYnJKB2SOB7f/o5+YgPJTfYC13g1gDSoWOU4KNyIig9zOymbMuEl+rJxUd2zQrG1zOAWdg4rrSaE9ZkBj4gwqNjX+p9/11T1WuBERGeR2VDRDez2T/I3gCUDG2I9+UoLyuZ1kJnvAG6SswwPt9RBqsbUmt9taFbitrc3WOoaDcDgMcMKrHWsquIjIINbYFqG8sQOjtZKJ2SHImQOOxFkGvzfyU/3UtoQpI4txlFpdUzmTbavH6XSSlpZGVVUVAIFA4LCr+sqJicfjVFdXEwgEcLlOLJ4o3IiIDGI7KpshHqWQSpJc8UE5S+rDCtJ8vFvaSDnZQCk02RtuAPLy8gC6A470D4fDQVFR0QmHR4UbEZFBbEdFk7W2TVIrJGVBcmKsEHsi8lOtcTeVsRRrE80E2IbBMAzy8/PJyclJuAXrhhKPx3PYVZCPl8KNiMggVd0coqYljLO1kvH5IchfOKi2WziS9IAbr9tByJ1MbdhFTnOltYlmAnS3OZ3OhNr9Wg5PA4pFRAapHRXNEGljjLMGn8sYtGvbfJhhGOQFfeD2UxFNgngUWtQdJMdO4UZEZBAyTdMab9NSyaSUDsgcB95ku8vqM3lBH2BQ3rWYXwJ0TcngoXAjIjIIlTV20NQWxtNexZhAaFCvbXM4uZ2L+VXG060DzQo3cuwUbkREBqEdFU3QUc94bxNur99alXgI6VqpuM5MpiNmqOVGjovCjYjIIBOLm+ysbHl/u4XcqQkx2LYvBTwugn43pieZqpAb2uog0m53WTJIKNyIiAwyJbWttHd0kBSqptAfHnJdUl3yU33gdFNO5z5Tar2RY6RwIyIyyOyoaIbWaiYmteFIyR4Sa9scTm7Q6pqqIMs6oHAjx0jhRkRkEAlH4+yp/kCXVN7JQ2Jtm8PpGndTEUvFNFG4kWOmcCMiMojsqW4h0t5CWryBXF/cGm8zRGWneHEYBm2OAM1RB7RUgHbmlmOgcCMiMojsqOhc2ya5AyNrHHiS7C6p37idDjKSPeBJoirsgXAbhJrsLksGAYUbEZFBoi0cpaSmuWeX1BCXm+IFw0mVkWkdaK60tyAZFBRuREQGiZ2VLcRba8l1t5OREhhya9scTk7noOIqusJNuY3VyGChcCMiMkjsqGiC5gpru4X8GdAHuycnupwULwBVsRTrQItabuSjDf1/GSIiQ0BjW4Sy6jqMjgYmpYSscDMMZKd4MQxodQRoiTqgWYOK5aMp3IiIDAJby61Wm0J/mOTcseAL2l3SgHA7HWQmWYOKK0MeCLdCqNnusiTBKdyIiCQ40zTZWloPLVVMDXZAwSy7SxpQ2Sk+a1Cxo3Pcjbqm5CMo3IiIJLiD9e001ZbjIcy4TC+kj7G7pAGVE+wcd2N2bsOgQcXyERRuREQS3HtlnQOJkztwjxgeA4k/qGsbhup4Z1ecpoPLRxhe/0JERAaZUDTG7oMVEGpmamoI8of+2jYflp1sDSpuNgK0Rh1quZGPpHAjIpLAdlW2EGkoJ8MdJW/EGPCm2F3SgPO4HKQHulYqdmtQsXwkhRsRkQS29WAdtFZxUrAdY8TwGkj8QbnBD69UXGFvQZLQFG5ERBJUfWuY0tISjHiUKdk+SB9td0m2yU758ErFCjdyZAo3IiIJamtZIzSVMzoQJnnMXDAMu0uyTddKxZXxzm45hRs5CoUbEZEEFI+bbNuzDyLtnJQeg7zpdpdkq67p4M0k0R4zoEXhRo5M4UZEJAHtr2ujufoAPkecseMng8tjd0m28rqcpAfcnYOKPRBq0aBiOSKFGxGRBLS1+CB0NDA5GMJVONfuchJCTtBaqbiye9yN1ruRw1O4ERFJMG3hKLv37ARg6qh88KfZW1CCyO1aqdjQSsVydAo3IiIJ5r3iMmLNNeR5I+RMPtXuchJGTteMqXiqdUB7TMkRKNyIiCQQ0zTZsnUrYDK9MA2CBXaXlDCyO2dMNcb9dMQMtdzIESnciIgkkJKKWhprK/A64kyadord5SQUn9tJqt8NnmRrpeJQi/UQ+RCFGxGRBPLOu2+DGWdKjh931vDa/ftY5AZ94HBSZWRZB9Q1JYehcCMikiCampvZu/8AACdPnzmsF+07kq71bqrNNOtAS5V9xUjCUrgREUkQb296AzMepzDdS+aok+wuJyFlJXeFm6B1oFXhRg6lcCMikgAi7S28u2c/ADOnzVCrzRF0DSqujweIxlHLjRyWwo2ISALY/vbrdERNUpP9jJ2gVpsjSfI48XucxN0BasMuaKuFWMTusiTBKNyIiNjMbG9k885iAGZMnYrDqR/NR2IYBtnJXnB6qI4ng2lCa43dZUmC0b8gERGbHdiyhpqQA48/hamT1WrzUbJSvIBBtWZMyREo3IiI2Km5kjd3HQTgpJOm4fO4bC4o8WV3DiquoWtQcbWN1UgiUrgREbGLaVL13mr2t3lwJGcxe8p4uysaFLJSrB3Sq2NJmCZquZFDKNyIiNilZicbimvBMJg06SRr9V35SJlJXpwOg5AjQHPUYc2YMk27y5IEonAjImKHaJiGrS+xs8ULwRHMmVBod0WDhtNhkJ7kAXeA6ogXoiHoaLS7LEkgCjciInYoeZW3KmKYLh9jxk3uXr9Fjk12shcMBzVGhnVA427kAxRuREQGWks1jXs38F6zHzLGMndslt0VDTrZXeNujHTrgMbdyAco3IiIDKR4DLY/zZt1fmL+TIpGFjIyPWB3VYNO1zYMNfHOGVMKN/IBCjciIgNp/3oa66p5rzUFMscxf2yG3RUNSl3deA1xP+G4AS3qlpL3KdyIiAyU5grY9ypv1icRzxhHUXa6Wm16KeBxkeR1YrqTrG0Y2uutgcUiKNyIiAyMaAi2/oO6kMF7sZGQlM2p4zLtrmpQy07xgtNNtZliHdCgYulka7h55ZVXuPDCCykoKMAwDJ588smjXr969WoMwzjksX379oEpWESkN0wTdjwPbXWsbcohnjGBsdnJjEjz213ZoNY97kbbMMiH2LrOd2trKzNmzOCLX/wil1122TE/b8eOHQSDwe6vs7Oz+6M8EZG+UboRqrZR2uFlj38ahsvN6eM1Q+pEdY27qe4eVKyWG7HYGm7OP/98zj///ON+Xk5ODmlpaX1fkIhIX6vdA7v/hWnCWscccAWZVpBKZrLWtTlR78+YsrZhMNRyI50G5ZibWbNmkZ+fzznnnMOqVauOem0oFKKpqanHQ0RkQLRUw9YnwYyz3TuNMkcebqehsTZ9JCPgweUwCDsDNEac1pibeNzusiQBDKpwk5+fzx/+8Acef/xxnnjiCSZNmsQ555zDK6+8csTn/OQnPyE1NbX7UVioJc5FZAC0N8CWv0E0TEdKIWsikwGDeWMySfZq5+++4HAYZCR7wOWnJuqFWAQ6GuwuSxLAoPoXNmnSJCZNmtT99YIFCzhw4AB33nknZ5555mGfs3z5cpYtW9b9dVNTkwKOiPSvjkZ4+yHoaIJAJq/5F9Fa3056wM3sojS7qxtSspO9VDWFqHJkMZ791qDigNYOGu4GVcvN4Zx66qns2rXriOe9Xi/BYLDHQ0Sk37Q3wOaHrP/606kadxlvl3cAcPbkHFzOQf9jN6FkdQ4qrjHTrAMtVfYVIwljULXcHM6mTZvIz8+3uwwREWiutLqiQi3gSyV28mdZuaWRuGkyITeZUZlJdlc45GQnd82Y0lo38j5bw01LSwu7d+/u/rq4uJjNmzeTkZFBUVERy5cvp7S0lAceeACAu+66i9GjRzN16lTC4TAPPvggjz/+OI8//rhdH0FExFK7xxo8HA1DUhac/Bk2lEWoagrhczs5e1KO3RUOSV3TwZtMHx0xA59abgSbw82GDRs4++yzu7/uGhtz5ZVXcv/991NeXs7+/fu7z4fDYW6++WZKS0vx+/1MnTqVZ599lgsuuGDAaxcRAawF+vathZJXrf+dVgjTLqO6w8HrxdbU5LMnZ5OkQcT9wud2kuJz0Ry3tmEY4Wy0VoN2aar9cGaYpmnaXcRAampqIjU1lcbGRo2/EZET094AO56D+hLr6xGzYdw5RHDw8Bv7qWkJMy4nmQtPzscwDFtLHcqe3FRKcU0rHwuvZkagFmZ/AVJH2l2W9LHj+f2tPyVERI5XPA5lm2DvKmv6sdMFEz8OedMBeHlrJTUtYZK8Ts6ZnKNg08+ykr0U17RSY2QAtda4G4WbYU3hRkTkeNSXwO6V7y/1n1YIky7onn68o6KZLaWNGAZ8fGq+uqMGQFaKB4Carg00tQ3DsKd/dSIix6LhgDW2pn6f9bXbB6PPtLqiOltmqpo7+Nc2a5zNvNEZFGUGbCp2eOnehiHWuQ2DZkwNewo3IiJHEotA1TYo2whN5dYxhxPyZ8DoM8DzfnhpC0d5anMZ4WicoowAp47VFgsDJT3gwdm5DUNT1EFqa5U1uFvdgcOWwo2IyIe11kD521DxDkSsBfhwOK0xNUULwJ/W4/JILM4zb5fT3BElPeDmEyfn43DoF+tAcToM0pM81MT81ITdpLo7INwC3hS7SxObKNyIiIAVaKq2QfV263938aVCwSzIPxk8hy7CF4ubPLelnNKGdrxuBxfNHIHP7RzAwgUgO9lDTXOIGiODcZRbKxUr3AxbCjciMjxFQ9Y4mvp9ULcX2mrfP+dwQvoYK9RkjAXH4bdMiMdNVrxXwd7qVlwOg4tmFJCR5BmY+qUHa9xNM7WkAeVWQM0cZ3NVYheFGxEZHuIxaCq1wkz9PmsMjRl//3xXoMmZDJkTrAHDRxHrDDbbK5pxGAafnFHAyHQNILZL96DirhlTrVqpeDhTuBGRock0rfVOusJMw35rgPAH+dMhffT7j48INF2isTjPbilnb3UrDsPggul5jMnSvlF2yky2Wszq4wGicXBpxtSwpnAjIkNHewM0lLwfaMJtPc97ApA26v0w86GBwceiNRTlmXfKKGvowOWwWmwUbOyX7HXhczvpiASoi7jIaa21Fls8QpeiDG0KNyIyeEXarUX1usJMe33P805XzzCTlH1C04Ormjp46u0ymjuieN0OLjy5gMIMdUUlAsMwyEr2cDDiozrkJSfean0/JGlK/nCkcCMig0c8Bo0HrQHA9fugpdLqfupiOCCY/36YCY6wxtKcINM02XSggbW7aojFTTKSPFw0o4B0DR5OKFkpXg7Wt1NrZAKtVrekws2wpHAjIomtvd4KM3XFVqD58LiZpCxrIHD6KEgr6vPdoJs6Iry4rZJ9NVYX19jsJJZMzdN07wSUldQ5qJg0YH/noOLJdpYkNlG4EZHEEo9D4wGo2QV1e6Ctrud5TxJkjHk/0PTTWiaxuMnmA/W8treOcDSOy2FwxsRsZoxM1UaYCap7j6lYsnVAg4qHLYUbEbFfPGa1ylTvgNpdPQcCGw5rh+eMMdaaM8m5/bqsvmma7Kpq4dXdNTS0Wa1EI9L8nDMlh8zkvm0Vkr6VmeTFMKAVH21Rg4A20By2FG5ExD4t1VDxNlS+1zPQuH3WWjNZE6yxM33c1XQ4sbjJzspmNpTUU9McAiDgcbJwfBZTC4JqrRkEPC4HqX43DdEkasJuijoaIBoGl8ZGDTcKNyIysGJRqNwCZZuhueL9454AZE+GrInW2Jk+GAh8LJo6IrxX2sR7ZY00d0StUlwOZhelM3tUGl6XxtYMJlnJXhraItREkygy66GtBoIFdpclA0zhRkQGRqTD2l374AYIt1rHDIe1RH7+jM5tDgYmSDS2RSiubWVnRTOlDe3dxwMeJzML05hRmKYBw4NUZrKH3VVQQyZQb427UbgZdhRuRKR/hVvhwOtQtsnqIgDwBWHEXMibdtjNKPuSaZq0hmOU1rezv66NA3VtNLb3nHFVmBFg+ohUxmUn4XJq0bfBLLt7G4agdUDjboYlhRsR6R+xCBx8E/avfz/UJGVB0amQc1Kft9LE4iYtHVGaOiI0dUSoaQlT0xyiuiVEezjW41qnwyAv1ce47GQm5iaT4nP3aS1in649pmrjAeImODRjalhSuBGRvle7B3a+AB2N1tcpuTD6TKsLyjAwTZNINE5HNEZHJEYkZhKNxa3/xuNEYyaRWJxovPO/MZNY/NBj0bj1nPZwjNZwtMd6fh9kGJCZ7KUoI0BRRoARaX48LrXQDEWpfjdup0HEGaAx4iRdG2gOSwo3ItJ3wm2w6wWo2o5pQp0jjdrsU6n3F9FQGaGh+ACN7RHaw3HiR0oiJ8DlMEj2uUjxuclM9pCd7CU7xUtGkge3upuGBYfDICPJS2UkQE3ETXq4zeoa7efuT0ksCjci0jdqdtP87vMUN8QobkunzDuGjuQiKHMCtYd9itNh4HU58LgcuJ0O3E4Dp8P6r8vhwOU0PvS/Hbgcnf91Grgc1jm/x0mKz4Xf7dSUbSEr2UNlUwfVZhoTqICWKmudJBk2FG5E5IREIhF2blzN1l17ONieBO4AZE8ETzJup0FWspe0gIf0gJv0JA9pfjd+jxOvy4nbaSiMSJ/LSvngNgwV0FqjcDPMKNyISK+0h2Ns3lvO2xvX097aDHgwggXkjZrM2NwgRRkBslO8OB0KLzKw3p8x1bk1hwYVDzsKNyJyXCKxOJv2N/DmzgOEy96FaIigF6ZNn82UKVMJauaR2Cwz2VqRuDHuJxw38GhQ8bCjcCMix2xvdQsvba+iub4aqreT7e7glJEuJpx2MY6kTLvLEwEg4HGR5HXSGkmiNuwiv7UGTLNf9ySTxNKrcFNcXMyYMeq/FBkuOiIxVu+oZlt5E7TWkNKwjYXZTUwemYMx/VJw++0uUaSHrGQvrR1Rajq85MdaoL0eAhl2lyUDpFdzI8ePH8/ZZ5/Ngw8+SEdHR1/XJCIJpKKxg7+8vp9t5U0YLZXMjmzkysIqpowbjTHzswo2kpCykr1gGNQY6daB1hp7C5IB1atw8/bbbzNr1iy+8Y1vkJeXx7XXXssbb7zR17WJiM3ePtDA3zYcoKk9Qmq0mn/zvcGirGbcI2bAlKUDtheUyPHqWqm4hlTrgMbdDCu9CjfTpk3j5z//OaWlpdx3331UVFRw+umnM3XqVH7+859TXa2R6SKDWTxusmpHFS9tryIWN5ngbeByz1oKfGEYMRsmnQ8OLYoniSsrxRpUXBNLtlau1oypYeWEfjq5XC4uueQS/va3v3H77bezZ88ebr75ZkaOHMkVV1xBeXl5X9UpIgMkEovzzJZyNu9vAOD0nBCfiK3C54hbu3dPWKyBmZLwMgIeHIZBh8NPS8yhDTSHmRMKNxs2bOD6668nPz+fn//859x8883s2bOHl156idLSUpYuXdpXdYrIAAhH4/xjcxl7qlpwOQw+MdbJKY3/xCAOuVOtFhsFGxkEXE4H6Ulu8CRRE3JZA4pjUbvLkgHSq9lSP//5z7nvvvvYsWMHF1xwAQ888AAXXHABjs5m6jFjxvD73/+eyZMn92mxItJ/QtEY/9hURmlDOx6Xg6WTkhhZ/DfrF0LGGJj8CQUbGVSykr3UtoSoCQcYYzZAWw2k5NldlgyAXoWb3/72t1x99dV88YtfJC/v8N8oRUVF3HPPPSdUnIgMjEgszlObrWDjdTu4ZFoW+XsftTbCTMmFqZdo8LAMOlnJXnZgUGtkAg3WuBuFm2GhV+Fm5cqVFBUVdbfUdDFNkwMHDlBUVITH4+HKK6/skyJFpP/E4ibPvlPOwXqrxeaymQXkljxtTZ31psD0T4PLa3eZIseta6Xi6njnNgwtmjE1XPRqzM24ceOoqTl0zYC6ujot7icyiJimycqtFRTXtOJ2GiydWUBu7WtQtxecLph2mRVwRAahrung9WYSMROtdTOM9CrcmKZ52OMtLS34fL4TKkhEBs76vbVsK2/GYRh84uQCRsZKYf/r1snJn4Rgvr0FipyAoM+Fx+Ug5gpQF3ZprZth5Li6pZYtWwaAYRjcdtttBAKB7nOxWIzXX3+dmTNn9mmBItI/3itr5PW9dQCcMyWHMUkReOsZ6+TIUyBnio3ViZw4wzDITvZSGg5QE3KRHWqBSLtW1R4GjivcbNq0CbBabrZs2YLH4+k+5/F4mDFjBjfffHPfVigifa6soZ0Xt1l/xc4bk8G0/GTY9GeIdFitNePOtrlCkb6RleKhtMFFrZkKdFjjbtJH2V2W9LPjCjerVq0C4Itf/CK//OUvCQaD/VKUiPSfllCUZ98pJxY3GZ+TzGnjMmH3i9BUDm4fnHSxZkbJkJGZ1LkNg5EOVFrjbhRuhrxezZa67777+roOERkA1syoMlpCUTKTPSyemotRvQMOvmldMPmT4E+ztUaRvpSV0hlu4snWAY27GRaOOdxceuml3H///QSDQS699NKjXvvEE0+ccGEi0vfW7KqmrKEDr9vBhScX4I22wI7nrJNF8yFrgr0FivSxzCRr+EQzATpiBj5NBx8WjjncpKamYnSuTpqamtpvBYlI/9hd1cKmzv2ilkzNIz3ghrcfh2gIggUwZpG9BYr0A5/bSdDvpikcoCbsYmRrNZimVtse4o453HywK0rdUiKDS2N7hBVbKwCYMyqdcdnJcPAtqN9nrWcz5UKNs5EhKyvZQ1Obn5p2LyNjLdDRAP50u8uSftSrdW7a29tpa2vr/rqkpIS77rqLFStW9FlhItI34nGTF96tIBSJk5/qY+H4LGirg70vWReM/RgEMuwtUqQfZSV7wXBQQ2eg0WJ+Q16vws3SpUt54IEHAGhoaGDevHn87Gc/Y+nSpfz2t7/t0wJF5MRs3F/fvRnm+dPycWLCtqetDTHTR8GI2XaXKNKvulYqrqFzSIXG3Qx5vQo3Gzdu5IwzzgDgscceIy8vj5KSEh544AF+9atf9WmBItJ7Vc0drNtTC8CiidmkBtxw4DVoKgOXRzt9y7CQ1bnHVG08CdPE2kBThrRehZu2tjZSUqz9ZlasWMGll16Kw+Hg1FNPpaSkpE8LFJHeicbivPBuBbG4ybicZKYWBK2/WPettS4Yfx74NDlAhr70gAenwyDsDNAUdSrcDAO9Cjfjx4/nySef5MCBA7zwwgssXrwYgKqqKi3sJ5Ig1u2ppaYlTMDj5NwpORimCdufhXjMmvKdN93uEkUGhMNhkJHkAXcS1SGXNeYsFrW7LOlHvQo3t912GzfffDOjR49m/vz5LFiwALBacWbNmtWnBYrI8TtQ18bG/fUAnHtSLgGPC0rfguYKcHlh4hJ1R8mwkpXsBZeHmpgfzDi01dpdkvSjXq1Q/KlPfYrTTz+d8vJyZsyY0X38nHPO4ZJLLumz4kTk+IWiMVZsrcQ0YfqIVGvad3sDFL9sXTD2LPCm2FmiyIDLTvGwrdygxsgAGqyuqZRcu8uSftKrcAOQl5dHXl5ej2Pz5s074YJE5MSs21NLU3uEoN/NmROzrQXLdq2EWARSR0KBWldl+OmaMVVrdg6d0DYMQ1qvwk1rayv/8z//w4svvkhVVRXxeLzH+b179/ZJcSJyfMoa2nn7QAMA507JweNyQOVWqN1tLdI36Xx1R8mwlNkZburjASJxcLdoUPFQ1qtw8+Uvf5mXX36ZL3zhC+Tn53dvyyAi9onFTf61zeqOOqkgyKjMJIi0w+6V1gVFCyApy94iRWyS5HHi9zhpdweoC7vI1YypIa1X4eb555/n2WefZeHChX1dj4j00pv76qjtnB115oRs6+CeVRBus0JN0QJ7CxSxkWEYZCV7OdCRRHWHi9xQsxX+3X67S5N+0KvZUunp6WRkaLl2kURR2xLijeI6AM6alIPf44TGg1D+tnXBxI9be0iJDGNZyR5wuKgxO9d3UuvNkNWrcPPDH/6Q2267rcf+UiJiD9O0uqNicZOx2UlMzE2GeBx2vmBdkH8ypBXaW6RIAugeVEyadUDjboasXv0p97Of/Yw9e/aQm5vL6NGjcbvdPc5v3LixT4oTkY/2zsFGyho68LgcnD05xxoDV7rRWo3Y5bWmfovI+3tMmcnWAbXcDFm9CjcXX3xxH5chIr3RFo7y6h5rh+PTxmUS9Lkh3PqBNW0WgSfJxgpFEkdmsgfDgDYjQGvUQZKmgw9ZvQo33/3ud/u6DhHphbW7aghF4uQEvcwYmWYd3LsaoiFrgbJ8rWkj0sXtdJDmd1MfSqIm7CKptdpaB0ozfoecXo25AWhoaOD//u//WL58OXV11kDGjRs3Ulpa2mfFiciRlTW0815ZEwBnT8rB4TA6BxG/Y10wYTE4ev1PXGRIykrxgttPTcQD0TB0NNpdkvSDXrXcvPPOO5x77rmkpqayb98+rrnmGjIyMvj73/9OSUkJDzzwQF/XKSIfEI+bvLTdalKfWhCkIM1vDSLetcK6IG+6tRqxiPSQlexll+GghjSgBVprwJ9mc1XS13r1Z92yZcu46qqr2LVrFz6fr/v4+eefzyuvvNJnxYnI4b1T2kh1cwif28npEzoX5ivfBM2V1iDicWfbW6BIgspK9gBQQ9d0cI27GYp6FW7efPNNrr322kOOjxgxgoqKihMuSkSOrDUUZV3nIOKF4zOtHb/DrbC3cxDxGA0iFjmSrhlTdfFk4ibWrEIZcnoVbnw+H01NTYcc37FjB9nZ2SdclIgc2drd1iDi3KCPaQWdf33ufdkaRJyco40xRY4i1e/G43IQdQaojzg1HXyI6lW4Wbp0KT/4wQ+IRCKAtaz1/v37ueWWW7jsssv6tEAReV9pQztby5owDPjY5K5BxKXvr0SsQcQiR2UYBplJHvAEqA27oK0O4jG7y5I+1qufgnfeeSfV1dXk5OTQ3t7OokWLGD9+PCkpKfz4xz/u6xpFhJ6DiKcVpJKX6rOmsX5wELFWIhb5SJnJXnB5qYkGwIxDW63dJUkf69VsqWAwyNq1a1m1ahVvvfUW8Xic2bNnc+655/Z1fSLSafPBBmo6BxEvHN85iLhiCzRXgMujlYhFjpE1qNig2kgH6q1xN8k5dpclfei4W27i8Tj33nsvn/zkJ/nqV7/Kn/70J9auXUtZWRmmaR7Xa73yyitceOGFFBQUYBgGTz755Ec+5+WXX2bOnDn4fD7Gjh3L7373u+P9CCKDTmsoyvo91l+Xp4/PsjbGjIasBfsARi0Eb7J9BYoMIu9vwxC0DmjczZBzXOHGNE0uuugivvzlL1NaWsr06dOZOnUqJSUlXHXVVVxyySXH9eatra3MmDGDX//618d0fXFxMRdccAFnnHEGmzZt4r/+67+48cYbefzxx4/rfUUGm1d31xCOxslL9TG1oPMHcsk6a5aUPx1GzLW3QJFBpCvcNJl+QnFD4WYIOq5uqfvvv59XXnmFF198kbPP7rmOxksvvcTFF1/MAw88wBVXXHFMr3f++edz/vnnH/P7/+53v6OoqIi77roLgClTprBhwwbuvPNODWSWIauyqYOt5dbsxLMmZVuDiNvr4eCb1gXjzwFnr3qYRYYlv8dJis9Fc0cStSEXBQo3Q85xtdw89NBD/Nd//dchwQbgYx/7GLfccgt/+ctf+qy4D1u/fj2LFy/ucWzJkiVs2LChe+bWh4VCIZqamno8RAYL0zR5eUc1pglT8lPIT/VbJ/a8ZM3wSB8NmeNtrVFkMMpM9oDH2mOKjiaItNtdkvSh4wo377zzDh//+MePeP7888/n7bffPuGijqSiooLc3Nwex3Jzc4lGo9TU1Bz2OT/5yU9ITU3tfhQWajaJDB47K1sobWjH7TTeH0Rcvw+qd4LhgPHnatM/kV7ISvaCw0WN2blWlBbzG1KOK9zU1dUdEi4+KDc3l/r6+hMu6miMD/0g7xrE/OHjXZYvX05jY2P348CBA/1an0hficTirNllNZfPHZ1Bis9t7R+1+1/WBQWzIFmLZor0RvegYiPNOqCuqSHluDrqY7EYLteRn+J0OolGoydc1JHk5eUdsr1DVVUVLpeLzMzMwz7H6/Xi9Xr7rSaR/vJWST3NHVFSfC7mjEq3DpZvhpZqcPtg9Om21icymHWHm3gypglGS6XNFUlfOq5wY5omV1111RHDQigU6pOijmTBggU8/fTTPY6tWLGCuXPn4na7+/W9RQZSc0eEDfvqADhjQjZup8MaE1DcuTHt6DPAE7CxQpHBLSPJg8MwCDmTaI46CCrcDCnHFW6uvPLKj7zmWGdKAbS0tLB79+7ur4uLi9m8eTMZGRkUFRWxfPlySktLeeCBBwC47rrr+PWvf82yZcu45pprWL9+Pffccw8PPfTQ8XwMkYT36u4aIjGTEWl+JuZ2rl+z71Ur4CRlaf8okRPkdBhkJLmpCVuDioOttVa3r7YvGRKOK9zcd999ffrmGzZs6DHzatmyZYAVou6//37Ky8vZv39/9/kxY8bw3HPP8fWvf527776bgoICfvWrX2kauAwpZQ3tbCtvxjCsqd+GYUBrLZS+ZV0w/hxwOO0tUmQIyEr2UtPio6bVx9h4k7UNg8axDQm2Lo5x1llnHXVV4/vvv/+QY4sWLWLjxo39WJWIfUzT5OWd1sDGk/KD5AR91ok9L1p74GRNgIyxNlYoMnRkpXihwqDGyASaoKVS4WaIUPubSALZVt5MRWMHHpfj/anftXush8MJ4z5mb4EiQ0jXoOJaOqeDt2o6+FChcCOSIMLROK/uttZrmjcmgySvy1qob/eL1gUj5kAgw8YKRYaWzGQPAHXxJKJxtNbNEKJwI5Ig3iqppyUUJdXvZlZhmnWwdKM1DsATsDbHFJE+k+J14XU7iLsD1EVcVreUDAkKNyIJoCUU5a2SrqnfWbicDgi3wb411gVjFllr24hInzEMw+qa8iRRE3Zb/+ZCLXaXJX1A4UYkAazrnPpdkOZjfE7n1O+SVyEasgY45p1sb4EiQ1R2shcMJzV0LpSp1pshQeFGxGbVzaHuXb/PmPDBqd+dswLHnaO1N0T6Sde4m9qubRg07mZI0E9MEZut2WXt+j0xN4WCtM5dv/eusqZ+Z46HjDH2FigyhL2/DUOKdUAzpoYEhRsRG+2raaWktg2nw+D0rqnfdcVQs8va9VtTv0X6VVayF8OAFiNAa9ShlpshQuFGxCbxuNm96/eMwjRSA527fu95ybpgxGxIOvyGsCLSNzwuB+kBD7iTqA65oK0OYhG7y5ITpHAjYpOt5U3UtITxuZ3MH9O5fk3lFusvR5dXU79FBkhOihdcHqpiSVZ3cGuN3SXJCVK4EbFBOBpn3R7rB+j8sRn43E5rZtTel60LRp+uXb9FBkh2ihcwqDI6W0o1Y2rQU7gRscGGkjpaQzHSAm5mjEyzDu5/DcKt4E+Hgtm21icynOSkWGtIVXUPKq62sRrpCwo3IgOsuSPCxpJ6AE4fn4XTYUBHIxx4w7pg3NngtHVPW5FhJSdozZhqNJPoiBlquRkCFG5EBtj6PbWHLti392WIRyGtELIm2lugyDDjczsJ+t3gSaY67LLGvZmm3WXJCVC4ERlAVc0d3Qv2nTmxc8G+pjKofM+6YNw5YBg2VigyPOWkeMHtpyrstca/dTTaXZKcAIUbkQFimiZrdtZgmjApL4X8VL/112HXrt950yCYb2+RIsNUTooXDAfV3dswaL2bwUzhRmSA7KttY3+dtWDfwnGdC/ZVb4fGg9YYmzGL7C1QZBizZkxBlZlmHdC4m0FN4UZkAMTjJms7F+yb2bVgXywKe1dbFxSeCr6gfQWKDHM5QWvGVF08mXBcg4oHO4UbkQHwwQX75nUt2Fe6AdobwJsMhfNtrU9kuEv2ukjyOjE9SdSEXNBcYXdJcgIUbkT6WSQW57W9tQDMG9O5YF+4FUpetS4YswhcHhsrFBHoXO/Gk0x1xA2hZgi12F2S9JLCjUg/e/tAA80dUVJ8LmaMTLUO7lsL0TCk5ELedHsLFBGgc1Cxw0mV2TWoWF1Tg5XCjUg/6ojEeHOftWDfgnGZuJwOaKmGsk3WBZr6LZIwuhbzq+qaMaWuqUFL4UakH23YV09HJEZWsocpeZ0DhveusqaAZ02A9FH2Figi3bKTrUHFtfFkYibQonAzWCnciPST5o4Im/ZbrTanjc/C4TCgrhhq94DhgHEfs7lCEfmgoN+F1+0g5k6iNqxBxYOZwo1IP3l9bx3RuMmIND9js5Ks1pq9q6yTBbMgkGFvgSLSg2EY3YOKq0Ju6GiCcJvdZUkvKNyI9IO61jDvllnLty+ckGVts1C1FZorrZlRoxfaXKGIHI41qNj1gZWK1XozGCnciPSDV3db2yyMzU5iRJq/c8G+l62ThaeCJ8neAkXksLoHFXetVKyuqUFJ4Uakj5U3trO7qgXDgIXjO7dZKNtobcTnTYbCefYWKCJHlJNiDSqujqcQN1G4GaQUbkT6kGmarN1VA8BJ+UGykr0QaX9/wb7RZ4DTbWOFInI0aX43bqdBxJ1MfcSptW4GKYUbkT60r7aNg/XtuBwGp47LtA7uXw+RDkjKgryT7S1QRI7K4TCsTTS7BhW3N1h/oMigonAj0kdM02TtbqvVZkZhGkFf5w/Gg29ZF4w9Gxz6JyeS6HJSfOBwUalxN4OWftKK9JHtFc3UNIfwuBycMrpzmve+NRCPQloRZI6zt0AROSa5nTuEV9H571hdU4OOwo1IH4jG4qzbY22OecroDPwep/XXXuV71gXjztY2CyKDRF5qZ7iJpVgrFTeX21uQHDeFG5E+sKW0kab2CMleF7OK0qyDe1dbC/flTIFggZ3lichxSA+48bodRFxdKxWr5WawUbgROUGhaIzXi+sAmD82A7fTAXV7ra0WHE4Yu8jmCkXkeBiGQW7nSsWVHW5or7cmBcigoXAjcoI27W+gPRwjPeBmakEqxOOwp3ObhRGzwZ9ub4EictzyUn3gdFMRT7UOaNzNoKJwI3ICOiIx3iqxNsdcMC4Lp8OAqvegpQpcXig6zeYKRaQ3ugYVV3Rtw6AZU4OKwo3ICdiwr55wNE52ipeJuckQi0DxK9bJUaeBJ2BvgSLSK12DimvjKYTjhvaYGmQUbkR6qTUUZfOBrlabTGtzzNK3rJ2EfUEYMcfmCkWkt5K9LlJ8LkxPEpUhDSoebBRuRHrpjX11RGIm+ak+xmYlWQMOS9ZZJ7XNgsiglxv8wKDitloNKh5EFG5EeqGpI8KWg40AnDYuy2q1OfAaREPWNgu502yuUEROlDWo2POBQcXqmhosFG5EeuH1vXXE4iaFGQGKMgMQaoGDb1onxyzSNgsiQ0Be16Bis3NQcZMW8xss9BNY5DjVt4bZWtYEwGkf3BwzFoVgPmRNsLE6EekrOUEvhgHNjiCtUQc0l9ldkhwjhRuR4/Ta3lripsnY7CQK0vzW5phlm6yTYxZpmwWRIcLrcpKR5AFPMhUhNzQp3AwWCjcix6G6OcSOymbAmiEFwL61EI9B+mjIGGNfcSLS57oHFYc8VvdzR5PdJckxULgROQ7r99ZimjAxN4WcFB+01kDlu9bJMWfaW5yI9Lm8oA8czvfH3WgTzUFB4UbkGFU0drCnqgXD+ECrTfEr1uaYWRMgdYS9BYpIn8tPfX9QsWmirqlBQuFG5Bit21MDwJT8oNUP31QO1TusMTZjtDmmyFCUmezF5TAIOZNpiDjVcjNIKNyIHIMDdW2U1LbhdBicOqar1eZl67+5UyE5277iRKTfOB0GOUEveFOsQcXN5VZrrSQ0hRuRj2CaJuv31AIwbUSQ1IAb6kugrhgMB4w+3eYKRaQ/5QZ94A5QEfZBNGytViwJTeFG5CPsq22jtKEdl8Ng3phM66+2vautkwUzwZ9uZ3ki0s/yUn1gGFTS2WrbVGpvQfKRFG5EjsI0ze6xNjMK00j2uqB2tzWo0Omydv4WkSGta6XiajONmAk0KtwkOoUbkaPYU91CVVMIj8vB3NHpVqtN11ibEXPBm2JvgSLS71L9bvweJ1FPClUht1puBgGFG5EjiMffH2szqzCNgMcFVVuhpRpcXig61eYKRWQgGIZhTQn3plDW4bbWt9IO4QlN4UbkCHZWNVPTEsbrdjB7VLq1CnHxGutk4Xxw++0tUEQGTEGaH5weymJp1gG13iQ0hRuRw/hgq83cURn43E4ofxva68ETgJGn2FyhiAykrsX8ys0MLeY3CCjciBzG1vImGtoi+D1OZhamQSwCJa9aJ0ctBJfH1vpEZGDlBn04HQatjhSaok613CQ4hRuRD4nG4ry212q1OWV0Bh6XA0rfsjbN8wUhf6a9BYrIgHM7HeSkeMEXtMbdNJVpMb8EpnAj8iHvlTXR3BEl2evi5JGp1sDB/eutk6NPt6aAi8iwk5/mB3cSZSEfREPWwGJJSAo3Ih8QicV5o7gOgHljMnA7HXDwDSvgBDIhd7rNFYqIXQo6F/MrI8s6oK6phKVwI/IB7xxspCUUJeh3M7UgCOFWOPCGdXLMmeDQPxmR4So/zZohWUsqHTEDGg/aXJEciX5Si3QKR+O8uc9qtZk/JgOX0wEl663BxCl5kD3J5gpFxE7JXhepfjemN0hFh1vhJoEp3Ih02nyggfZwjLSAm5Pyg9DRCGUbrZNjF4Fh2FugiNiuIM0P3iBlIY+1NESo2e6S5DAUbkSAjkiMDSVWq82CcZk4HAbse9VauC+tCNLH2FyhiCSCkel+cLg4GO/cRFOtNwlJ4UYE2Li/nlAkTlayh4k5KdBaCxXvWCfVaiMinUZ0jrupMNOJxFG4SVAKNzLstYdjbNrfAHyw1eYVaw2LzPGQOtLeAkUkYaQF3CR5ncQ8QSpCbmg8YHdJchgKNzLsbSipIxyNkxP0Mi47GZoroGq71VozdpHd5YlIAjEMgxFpAfAFKW33QEuVteaNJBSFGxnWWkJR3j7QAMBp47IwDAOKX7FO5kyB5Bz7ihORhDQi3Q9OL6WxNKuFV11TCUfhRoa1N/fVEYmZFKT5GJ0ZgIb9ULsHDAeMPsPu8kQkAXWNuyk304mZKNwkIIUbGbaaOiJsOdgIdLbaAOx92TqZPwMCGbbVJiKJKyvZg8/tJOJOpUrjbhKS7eHmN7/5DWPGjMHn8zFnzhzWrFlzxGtXr16NYRiHPLZv3z6AFctQ8cbeOmJxk8KMAIUZAajba/0F5nDBqNPsLk9EEpRhGFbXlC+V0vbOTTRjEbvLkg+wNdw88sgj3HTTTdx6661s2rSJM844g/PPP5/9+/cf9Xk7duygvLy8+zFhwoQBqliGioa2MO+VNQHWDClME/autk6OmG3t/i0icgQj0vzg9nEgErTWw1LXVEKxNdz8/Oc/50tf+hJf/vKXmTJlCnfddReFhYX89re/PerzcnJyyMvL6344nc4BqliGitf21hE3TcZkJVk/pKq2WbMeXB4oWmB3eSKS4Aoz/IC1iWbMxBqvJwnDtnATDod56623WLx4cY/jixcvZt26dUd97qxZs8jPz+ecc85h1apVR702FArR1NTU4yHDW21LiO0VH2i1icfenyFVOB88ARurE5HBIDvZi9/jJOxOtfaZUrhJKLaFm5qaGmKxGLm5uT2O5+bmUlFRcdjn5Ofn84c//IHHH3+cJ554gkmTJnHOOefwyiuvHPF9fvKTn5Camtr9KCws7NPPIYPP+r21mCaMz0kmN+iD8retPWI8ARh5it3licggYBgGhekB8KdyoN0DzeUQDdtdlnRy2V2A8aFl7U3TPORYl0mTJjFp0vs7My9YsIADBw5w5513cuaZZx72OcuXL2fZsmXdXzc1NSngDGMVjR3sqmzBMDpbbaJh2LfWOjlqIbi89hYoIoNGYYafnZU+DkRTOTXeCk0HIWOs3WUJNrbcZGVl4XQ6D2mlqaqqOqQ152hOPfVUdu3adcTzXq+XYDDY4yHDk2marNlVDcCU/CBZyV4o3QDhVvClQv5MewsUkUGlMD0AGJSbGdY+U+qaShi2hRuPx8OcOXNYuXJlj+MrV67ktNOOfRrupk2byM/P7+vyZAgqqW3jYH07LodhtdpE2mH/a9bJMWeC0/aGTBEZRNICblJ8LmLeVMo7PAo3CcTWn+bLli3jC1/4AnPnzmXBggX84Q9/YP/+/Vx33XWA1aVUWlrKAw88AMBdd93F6NGjmTp1KuFwmAcffJDHH3+cxx9/3M6PIYOAaZqs3V0DwMmFaQR9btizxtoTJjkbcqfaXKGIDDaGYTAyPcC2ljQOtHooaiq3fqaoe9t2toabz3zmM9TW1vKDH/yA8vJypk2bxnPPPceoUaMAKC8v77HmTTgc5uabb6a0tBS/38/UqVN59tlnueCCC+z6CDJI7Khspro5hMflYN7oDOhogoNvWSfHnGVtkikicpwKM/xsK/exP5rGQrPFar3J0tprdjNM0zTtLmIgNTU1kZqaSmNjo8bfDBOxuMmf1u2jsT3CaeMymT82E7Y/Z82SSiuEmZ9TuBGRXmnuiPB/a4ox6vZwbebb+Itmw8TFH/1EOW7H8/vb9u0XRPrbltJGGtsjJHmdzCpKh9ZaqHjHOjn2LAUbEem1FJ+brBQvpi+V/W0eqC+2uyRB4UaGuHA0zut7awGYPyYTj8sBxaut7RayJkDqSHsLFJFBb1RGAHxp7Gv3QVsddDTaXdKwp3AjQ9rG/fW0hWOkBdxMG5FqbXBXvdNqrRmzyO7yRGQIGJ2ZBA4X++NZmCZQp9YbuyncyJDVFo7yVkk9AKeNy8Jp8P7mmLnTrFlSIiInqCDNh9tp0OJKpybsUtdUAlC4kSHrjeI6wtE4OUEvE3OToW4v1JeAwwljzrC7PBEZIlxOByPTra6pkjaP9XNmeM3VSTgKNzIkNbZFeOeg1e99+vgsDNOEPS9ZJ0fMtlYkFhHpI6MyA+BNoSQUsBYIbS63u6RhTeFGhqRX99QQi5sUZQQYlZlkzY5qrQG3z9pDSkSkD43OTALDQamZTThuQO0eu0sa1hRuZMgpa2hnR0UzhgFnTMiyNscs7tw5ftTp4PbbW6CIDDlpATepfjcxf4Y1JbxO4cZOCjcypJimySs7rc0xT8oPkhP0wYHXrM0x/elWl5SISB8zDIOx2UngT2dvmxeayiHUYndZw5bCjQwpOytbKG/swONycNr4LGubhQOvWyfHnmUNJhYR6Qdjs5LB6aU4mtk5JXyv3SUNWwo3MmREY/HuzTHnjEon2euCfWsgFoXUEZA9yeYKRWQoG5Hux+Ny0ObJpCLkVteUjRRuZMjYdKCBpvYIKT4Xc0alQ3MlVGyxTo47R9ssiEi/cjoMa2CxP529rV6r5SYes7usYUnhRoaE1lCUN4rrAFg4Pgu3w7Cmfpsm5Ey2Wm5ERPrZ2Owk8KSwtyPZmszQeNDukoYlhRsZEl7bW0s4Gic36GNyXgrU7IT6feBwWWNtREQGwOjMJAyHQY0zm8aIA2p22V3SsKRwI4NeVVMHW0qtBfsWTcrGiEdh94vWycJ51iwpEZEB4Pc4KUjzQyCDPa0+6w8trVY84BRuZFAzTZNVO6owTZicl8KINL81O6qjEbwpULTA7hJFZJgZn5MM/nR2t/mtn0UtlXaXNOwo3MigtrW8ibIGa+r36ROyrB8k+9dbJ8d9DFweewsUkWFnfE4yGE7KyKEl6rBab2RAKdzIoNURibF2lzX1e/6YDFJ8btizypr6nVYIOVNsrlBEhqOgz01+qg8zkMmeVi9U77C7pGFH4UYGrfV7a2kLx8hI8jCrKB0a9kPVNmvK9/jzNPVbRGwzITcZAhnsavVb+9q11dld0rCicCODUnVziLcPNABw9qQcnJiwa6V1Mn8mpOTaVpuIyPjsFHC4OWhm0RY11DU1wBRuZNAxTZNV261BxBNzUyjKDEDZRmipApcXxpxpd4kiMsylBtzkBru6pnxWq7IMGIUbGXS2lTdT2tCO22lwxsTO/aP2rrZOjj0LPAE7yxMRAbq6prLY0eqH5gp1TQ0ghRsZVNrCUV7ZZe36PX9sJkGfG3avhFjEWoW4YJbNFYqIWCbmpoDTw0Ezm+aoQ603A0jhRgaVV3ZW0x6OkZXiZXZRurX6Z/VOMBww8eMaRCwiCSPV72ZEmh8zKZudzT6o2mp3ScOGwo0MGsU1rWwrb8Yw4LwpuTjjEdi1wjpZeAok59hboIjIh0zOT4FAJttaA9asqZYqu0saFhRuZFAIR+O8uM1a5XNWUTp5qT7Yt8Yab+NLhVGn21yhiMihJuam4HS5qXZkUxNyqvVmgCjcyKCwbk8NzR1Rgn43C8ZmQnMlHNxgnZy4RCsRi0hC8rmdjM5KgqRstrf4oXKr9poaAAo3kvDKG9vZ3LmmzblTcvA4gB3PgRmH7EmQOc7W+kREjmZKXgoEMtjemky8vREaSuwuachTuJGEFo3F+dfWSkwTpuQHGZWZZO0d1VwBbh9MOM/uEkVEjmpMVhI+j4dmby4lbR4of8fukoY8hRtJaOv31lLTEibgcbJoYrbVHVXyqnVy/HnWzt8iIgnM5XQwJT8FknN4t8kPNTsg0mF3WUOawo0krAN1bbxVUg/AuSfl4ncB25+BeAyyJkDuVHsLFBE5RtNGpII3hb2RdFpDcajWmjf9SeFGElJHJMYL71VgmtYPhXHZyVCyzppG6fZrTRsRGVSykr0UpPmJJ+Wytdmnrql+pnAjCWn1jmqaO6Kk+t2cOTHLGmNTss46OXEJeJPtLVBE5DhNLUjt7JoKYDaWac2bfqRwIwlnV2Uz28qbMAz4+LQ8vIZpdUeZcciZDDlT7C5RROS4TcxNweP10eDOsQYWl260u6QhS+FGEkpLKMq/tll/zcwbnUFBmh/2vAQt1daGmBMW21yhiEjveFwOTioIQjCfzY0BqHxXA4v7icKNJIx43OSf71bQEYmRE/Qyf2wmVO+A0resCyZ/EjxJ9hYpInICZhWmYfjSKI6kU9cetwKO9DmFG0kY6/fWcqCuDY/LwfnT8nGGGmH7s9bJwnlarE9EBr20gIcx2UkQzLNab0rf0orF/UDhRhLC3uoW3iiuA+DcKblk+J2w7SmIhiCYD2PPsrdAEZE+MqswHZJy2NaaTEdzPdTttbukIUfhRmzX2B7hhfesTTFnFqUxKS/F2hSzsdTaM+qkpeBw2lyliEjfKMzwkxUMEA7ksaXJDwdet7ukIUfhRmwVjcV59p1yOiIx8lN9nDkh2/orZv9r1gWTLgB/ur1Fioj0IcMwmDMqA4IFbGxMJlJbAk1ldpc1pCjciK1e2VVNZVMHPreTC07Ox9lRD1v/YfVBF8zUtG8RGZIm5aWQmpJCmy/X2pJh/3q7SxpSFG7ENm8faODtA40YBpw/LY+gKw7vPm5NjQzmW3tHiYgMQU6HwSmjMyA4grcakohW7YLWWrvLGjIUbsQW+2paWb2jGoAFYzMZnRmAbU9Da421+vC0y8DpsrlKEZH+MyU/hZRgKs2ebLY2+WD/OrtLGjIUbmTA1bSEeHZLOXHT5KSCIPPGZEDxK1CzCxwumHqpdvsWkSHP5XQwZ1Q6pBbyRn0SkfL3rD/w5IQp3MiAag1FeXJTKeFonJHpfs6dkotRvf39faMmfRxSR9hbpIjIAJk+IpWU1HSaPdm83eC3/tCTE6ZwIwMmEovz1NtlNHdESQ+4uXBGAc6mg9a+UQCFp0DedHuLFBEZQC6ng9PGZUHaKN5oSKajYqe1UbCcEIUbGRCxuMlzW8qpaLRmRi2dOQJfqBa2PAqxqLX68NiP2V2miMiAm5yXQlZmBiFfDm/UJ8He1Vq1+AQp3Ei/i8dNnn+3nL3VrbgcBhfOyCfd0QrvPGKtQJw6AqZeAg59O4rI8ONwGJwxPgvSi9jclExdxX6o3W13WYOafptIv4rHTVZsrWBXZQtOh8GFMwoYGYjD249AqAWSsmD6p8HptrtUERHbjMoMMCYvi1jKCFZVp2Du+pfVqi29onAj/cY0TV7cXsW28mYchsEF0/MZneq0Wmza68GXCjM+C26/3aWKiNjKMAzOnpSDK72I/eEUdla3w8E37C5r0FK4kX5hmiard1bzbqm1SN/Hp+UxPt0J7zwMLVXgCVjBRlO+RUQASA24OWVcNmSM5pXaFDr2roO2OrvLGpQUbqTPxeMmK7dWsnl/AwDnnZTLpAwHbP4LNJVbLTUnfxYCGfYWKiKSYOaOSic9ewQtrgxervTDjuc1uLgXFG6kT0VicZ5+p4z3ypowDCvYTM0wYNNfoKUaPEkw83OQkmt3qSIiCcfldHDe1DyMrAlsbUli98FKKNtod1mDjsKN9JmOSIy/byztnhX1yZMLmJYRt1ps2mrBF4RZn4fkbLtLFRFJWCPS/MwdXwDpo/lXdZDWHavVPXWcFG6kTzR3RHh0wwFKG9rxuh1cMnsE433NsOlBaG8Af7rVYqOuKBGRj7RgXCbZI0bT7k7n+bIk4u/+XbOnjoPCjZyw0oZ2HnpjPzUtYZK9Lj49p5CR4X1WsOma7j3rc+BPs7tUEZFBwekwOH9aPp7cSRwIJ/FqSRvs/pfdZQ0a2nZZes00TTYfaOCVnTXETZOsZA8XzSggtWL9+3tFZY6Hky4Cl9feYkVEBpnMZC/nzRjFsx2tbKh8j7ydW5kQLID8k+0uLeEp3EivRGJxXtxWybbyZgAm5aVw7oR0PLufhuqd1kWF82Ds2Vp5WESklybmplAxeSxvhZr4Z6VJ4O0XGeFLhfRRdpeW0PRbR45bTUuIR9480L0436JJ2ZxfGMOz+X4r2DicMPkTMP4cBRsRkRN0+vgsxo6fQjSQzVOlKdS+9SS01tpdVkLTbx45ZvG4yRvFdfz19f1UN4cIeJxcOiuf2bF3MTb/pXPV4SDMvFzNpiIifcThMDh/egH5Y6fT4U7l7wcC1L/5iGZQHYW6peSY1LaEeOG9SiqbOgAYm53EOaM8JBc/AY2l1kW5J8GEJeD22VipiMjQ43E5WDq7iEcjEWr3bOSx4mYuNR8mc74WRD0chRs5qnA0zlsl9by5r45Y3MTrdnD2hAwmh9/DeOc1iEfB5YGJH4fcqXaXKyIyZPk9Tj41fyyPAzV7NvJYcQsXxR8mf95lWhj1QxRu5LBM02RbeTPr9tTQ3GGtrTAmM8C5ea0k73/EWrsGrEFtky7QNG8RkQEQ8Lj41PxxPGFA1e6NPFbcynntf2Pyggsgc5zd5SUMhRs5xIG6Ntbsqunuggr63ZyeDxNb1mHs2GNd5E2xBgxnTwbDsLFaEZHhxe9x8qlTx/FPv4e9W9/i+XKDypeeZ+GcmbjGLNTPZBRupJNpmhTXtLKhpJ7S+nbA6uOdnxNnZvQtXPt3WRc6nNYU76LTrO4oEREZcF6XkwtnjeLVZD8bNm5gY0MlB9a+x/nVpWTOON+a3DGMGaY5vLYbbWpqIjU1lcbGRoLB4f1/PkAsbrKjopm3SuqoaQkD4DRgWlqYU433CDTttS40DMiZAqPP0OA1EZEEsqe6hZWvbaa9cjdOYszJDDNvzim4C+cMqVac4/n9rXAzTFU1d7C1rIkdFc20hWMAeIwo0301zGI7KZHONRQMhzULqug0SMq0sWIRETmSllCUFzftZu/OLRBqJsUV49SRXk6afTqOzLF2l9cnFG6OYjiHm8b2CHuqW9ha1kR1c8g6GI+SFGtklq+C6ezGZ3RuzOZ0Q+40qwtKLTUiIgnPNE32VDazesMmmiv3QTxGujvK7MIgU6bPxZ09blC35CjcHMVwCjfxuEl5UwfF1a0U17R0dzsRDeEM1TPWWc1JzgOM9rXj6Pp+D+ZD/kzIOUljakREBqFILM47xRW8+fY7tNeXg2nic8SZku3hpAkTyB4zDcOXYneZx21QhZvf/OY3/PSnP6W8vJypU6dy1113ccYZZxzx+pdffplly5bx3nvvUVBQwLe+9S2uu+66Y36/oRxuQtEYlY0hShvaKW9sp7yxg3A0BpEOCDXhCDWRb9QywVPH5JQO/M7O/+sDGZA1AXKmaq0EEZEhIhSN8V5xGZvf20pjbQXErSEI6Z4Y43JTGV00ivzCcbiSMwdFi87x/P62dbbUI488wk033cRvfvMbFi5cyO9//3vOP/98tm7dSlFR0SHXFxcXc8EFF3DNNdfw4IMP8uqrr3L99deTnZ3NZZddZsMnsEckFqexPUJ9a5ialjA1LSFqmkM0trRiRtoh0g6RVgi34o02M8bXxpikMKOzQvicpvVNnJIHmRMgayIkZQ2Kb2wRETl2XpeT2RMKmTluJPsq69m2czt79x+kvqOFDQda2HDgPZzGu+QmGRRkpZOfl092dj7BzFyMQb7SvK0tN/Pnz2f27Nn89re/7T42ZcoULr74Yn7yk58ccv23v/1tnnrqKbZt29Z97LrrruPtt99m/fr1x/SeidpyY5omkZhJeyRGezhGWzhKWzhGeyhCa0eItvY2mlraaWptp7WjA2JhiEWs/0ZDVqAx4wRdMQp8YQr8EQp8ETI9URxOJyTnQlqR9UgdCS6v3R9ZREQGWCgao6S0gj1793CgoorW1hYw4z2ucRkm6QEXGSlJpAVTSAokkZScRHJSCkkpQfyBFFwe34D/UTwoWm7C4TBvvfUWt9xyS4/jixcvZt26dYd9zvr161m8eHGPY0uWLOGee+4hEongdrv7rd6PEorG2HKwkWjcJB43iZnm+/+7uYp4ay3RWIxYLEY0ZhKOxojETMKxGKGoSSQWx4ybYEatpsN4zNra4EPfdF28jjhp7hhZ3iiZgSjZniiZ3hhJSUngz4DkbCvQJOdCINNan0ZERIY1r8vJxFEjmDhqBKZp0tTaQWnpAcoqyiivrqW+uY1oNEJ1a4zq1iaoaDrs6zgN8LodeJxOvG4nXrcLj9uN2+XE6XSSkexlzmnnDvCne59t4aampoZYLEZubs8xHrm5uVRUVBz2ORUVFYe9PhqNUlNTQ35+/iHPCYVChEKh7q8bGxsBKwH2pZZQlJWb9x3+ZEMJNJUd0+s4DRO/M47fGSfgiONzxgk44/jdBkGfl2CSl5SAH18gCdxJ4PGDNwi+dGK+IE3OD/1fGgdaWk/os4mIyNBkACPz8xiZnwdYE1Gamluor6umoaGOpuYmWts7aG0P0RYK0xaKEI9bHT5H+82Sn+xgwrR5fVpr1+/tY+lwsn2FYuNDzVqmaR5y7KOuP9zxLj/5yU/4/ve/f8jxwsLC4y1VREREjtEXv3R9v7xuc3MzqampR73GtnCTlZWF0+k8pJWmqqrqkNaZLnl5eYe93uVykZl5+AXmli9fzrJly7q/jsfj1NXVkZmZedQQ1VeampooLCzkwIEDCTXGZzDQves93bsTo/vXe7p3vad7d3SmadLc3ExBQcFHXmtbuPF4PMyZM4eVK1dyySWXdB9fuXIlS5cuPexzFixYwNNPP93j2IoVK5g7d+4Rx9t4vV683p6DZ9PS0k6s+F4IBoP6Zu0l3bve0707Mbp/vad713u6d0f2US02XRz9XMdRLVu2jP/7v//j3nvvZdu2bXz9619n//793evWLF++nCuuuKL7+uuuu46SkhKWLVvGtm3buPfee7nnnnu4+eab7foIIiIikmBsHXPzmc98htraWn7wgx9QXl7OtGnTeO655xg1ahQA5eXl7N+/v/v6MWPG8Nxzz/H1r3+du+++m4KCAn71q18NqzVuRERE5OhsH1B8/fXXc/31hx90dP/99x9ybNGiRWzcuLGfq+o7Xq+X7373u4d0jclH073rPd27E6P713u6d72ne9d3bN9+QURERKQv2TrmRkRERKSvKdyIiIjIkKJwIyIiIkOKwo2IiIgMKQo3feA3v/kNY8aMwefzMWfOHNasWXPEa8vLy7n88suZNGkSDoeDm266aeAKTUDHc++eeOIJzjvvPLKzswkGgyxYsIAXXnhhAKtNLMdz79auXcvChQvJzMzE7/czefJkfvGLXwxgtYnleO7dB7366qu4XC5mzpzZvwUmuOO5f6tXr8YwjEMe27dvH8CKE8fxfu+FQiFuvfVWRo0ahdfrZdy4cdx7770DVO0gZsoJefjhh023223+8Y9/NLdu3Wp+7WtfM5OSksySkpLDXl9cXGzeeOON5p/+9Cdz5syZ5te+9rWBLTiBHO+9+9rXvmbefvvt5htvvGHu3LnTXL58uel2u82NGzcOcOX2O957t3HjRvOvf/2r+e6775rFxcXmn//8ZzMQCJi///3vB7hy+x3vvevS0NBgjh071ly8eLE5Y8aMgSk2AR3v/Vu1apUJmDt27DDLy8u7H9FodIArt19vvvcuuugic/78+ebKlSvN4uJi8/XXXzdfffXVAax6cFK4OUHz5s0zr7vuuh7HJk+ebN5yyy0f+dxFixYN63BzIveuy0knnWR+//vf7+vSEl5f3LtLLrnE/PznP9/XpSW83t67z3zmM+b/+3//z/zud787rMPN8d6/rnBTX18/ANUltuO9d88//7yZmppq1tbWDkR5Q4q6pU5AOBzmrbfeYvHixT2OL168mHXr1tlU1eDQF/cuHo/T3NxMRkZGf5SYsPri3m3atIl169axaNGi/igxYfX23t13333s2bOH7373u/1dYkI7ke+9WbNmkZ+fzznnnMOqVav6s8yE1Jt799RTTzF37lzuuOMORowYwcSJE7n55ptpb28fiJIHNdtXKB7MampqiMVih+xinpube8ju5dJTX9y7n/3sZ7S2tvJv//Zv/VFiwjqRezdy5Eiqq6uJRqN873vf48tf/nJ/lppwenPvdu3axS233MKaNWtwuYb3j8ze3L/8/Hz+8Ic/MGfOHEKhEH/+858555xzWL16NWeeeeZAlJ0QenPv9u7dy9q1a/H5fPz973+npqaG66+/nrq6Oo27+QjD+19qHzEMo8fXpmkeckwOr7f37qGHHuJ73/se//jHP8jJyemv8hJab+7dmjVraGlp4bXXXuOWW25h/Pjx/Pu//3t/lpmQjvXexWIxLr/8cr7//e8zceLEgSov4R3P996kSZOYNGlS99cLFizgwIED3HnnncMq3HQ5nnsXj8cxDIO//OUv3bth//znP+dTn/oUd999N36/v9/rHawUbk5AVlYWTqfzkNRdVVV1SDqXnk7k3j3yyCN86Utf4tFHH+Xcc8/tzzIT0oncuzFjxgAwffp0Kisr+d73vjesws3x3rvm5mY2bNjApk2buOGGGwDrF45pmrhcLlasWMHHPvaxAak9EfTVz7xTTz2VBx98sK/LS2i9uXf5+fmMGDGiO9gATJkyBdM0OXjwIBMmTOjXmgczjbk5AR6Phzlz5rBy5coex1euXMlpp51mU1WDQ2/v3UMPPcRVV13FX//6Vz7xiU/0d5kJqa++70zTJBQK9XV5Ce14710wGGTLli1s3ry5+3HdddcxadIkNm/ezPz58weq9ITQV997mzZtIj8/v6/LS2i9uXcLFy6krKyMlpaW7mM7d+7E4XAwcuTIfq130LNtKPMQ0TW175577jG3bt1q3nTTTWZSUpK5b98+0zRN85ZbbjG/8IUv9HjOpk2bzE2bNplz5swxL7/8cnPTpk3me++9Z0f5tjree/fXv/7VdLlc5t13391jSmlDQ4NdH8E2x3vvfv3rX5tPPfWUuXPnTnPnzp3mvffeawaDQfPWW2+16yPYpjf/Zj9ouM+WOt7794tf/ML8+9//bu7cudN89913zVtuucUEzMcff9yuj2Cb4713zc3N5siRI81PfepT5nvvvWe+/PLL5oQJE8wvf/nLdn2EQUPhpg/cfffd5qhRo0yPx2POnj3bfPnll7vPXXnlleaiRYt6XA8c8hg1atTAFp0gjufeLVq06LD37sorrxz4whPA8dy7X/3qV+bUqVPNQCBgBoNBc9asWeZvfvMbMxaL2VC5/Y733+wHDfdwY5rHd/9uv/12c9y4cabP5zPT09PN008/3Xz22WdtqDoxHO/33rZt28xzzz3X9Pv95siRI81ly5aZbW1tA1z14GOYpmna1GgkIiIi0uc05kZERESGFIUbERERGVIUbkRERGRIUbgRERGRIUXhRkRERIYUhRsREREZUhRuREREZEhRuBEREZEhReFGREREhhSFGxERERlSFG5EZFAIhULceOON5OTk4PP5OP3003nzzTcBqK+v53Of+xzZ2dn4/X4mTJjAfffd1/3cb3/720ycOJFAIMDYsWP5zne+QyQSseujiEg/c9ldgIjIsfjWt77F448/zp/+9CdGjRrFHXfcwZIlS9i9eze33XYbW7du5fnnnycrK4vdu3fT3t7e/dyUlBTuv/9+CgoK2LJlC9dccw0pKSl861vfsvETiUh/0caZIpLwWltbSU9P5/777+fyyy8HIBKJMHr0aG666SbWrFlDVlYW99577zG93k9/+lMeeeQRNmzY0J9li4hN1HIjIglvz549RCIRFi5c2H3M7XYzb948tm3bxn/+539y2WWXsXHjRhYvXszFF1/Maaed1n3tY489xl133cXu3btpaWkhGo0SDAbt+CgiMgA05kZEEl5XA7NhGIccNwyD888/n5KSEm666SbKyso455xzuPnmmwF47bXX+OxnP8v555/PM888w6ZNm7j11lsJh8MD/jlEZGAo3IhIwhs/fjwej4e1a9d2H4tEImzYsIEpU6YAkJ2dzVVXXcWDDz7IXXfdxR/+8AcAXn31VUaNGsWtt97K3LlzmTBhAiUlJbZ8DhEZGOqWEpGEl5SUxH/+53/yzW9+k4yMDIqKirjjjjtoa2vjS1/6Erfddhtz5sxh6tSphEIhnnnmme7QM378ePbv38/DDz/MKaecwrPPPsvf//53mz+RiPQnhRsRGRT+53/+h3g8zhe+8AWam5uZO3cuL7zwAunp6Xg8HpYvX86+ffvw+/2cccYZPPzwwwAsXbqUr3/969xwww2EQiE+8YlP8J3vfIfvfe979n4gEek3mi0lIiIiQ4rG3IiIiMiQonAjIiIiQ4rCjYiIiAwpCjciIiIypCjciIiIyJCicCMiIiJDisKNiIiIDCkKNyIiIjKkKNyIiIjIkKJwIyIiIkOKwo2IiIgMKQo3IiIiMqT8f73aEOswnTc8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "io = CMIP6_IO()\n",
    "model_var=\"osa\"\n",
    "LMES = ['Barents Sea', 'Northern Bering - Chukchi Seas']\n",
    "scenario = \"ssp245\"\n",
    "lightpath = \"light_sensitivity\"\n",
    "prefix=f'{lightpath}/{scenario}'\n",
    "models = ['CanESM5']\n",
    "stats  = {}\n",
    "factors = [\"snow_sensitivity\", \"no_chl\", \"no_wind\", \"no_osa\", \"no_meltpond\", \"no_clouds\"]\n",
    "factors = [\"no_osa\"]\n",
    "var_name=\"osa\"\n",
    "\n",
    "period=\"1979-01-01-1989-01-16\"\n",
    "df = None\n",
    "create_csv=True\n",
    "combos, files_on_gcs = list_all_model_combinations(prefix, period)\n",
    "max_plots=1\n",
    "counter=0\n",
    "for LME in LMES:\n",
    "    for factor in factors:\n",
    "            for combined_name in combos.combined_names:\n",
    "                if model_var in combined_name:\n",
    "                    \n",
    "                    model = combined_name.split(\"_\")[1]\n",
    "                    if model in models:\n",
    "                        file1 = f\"light_sensitivity/ssp245/{combined_name}_{period}_scenario_{factor}_ssp245.nc\"\n",
    "                        file2 = f\"light_sensitivity/ssp245/{combined_name}_{period}_scenario_osa_ssp245.nc\"\n",
    "\n",
    "                        ds1 = io.open_dataset_on_gs(file1)\n",
    "                        ds2 = io.open_dataset_on_gs(file2)\n",
    "                    \n",
    "                        if ds1 and ds2:\n",
    "                            ds1=ds1.resample(time='1M').mean() #.sel(time=slice(\"1979-01-01\", \"1983-01-01\"))\n",
    "                            ds2=ds2.resample(time='1M').mean()\n",
    "\n",
    "                            ds1 = get_data_within_LME(ds1, LME).dropna(dim = 'time', how = 'all')\n",
    "                            ds2 = get_data_within_LME(ds2, LME).dropna(dim = 'time', how = 'all')\n",
    "                            \n",
    "                            ds1 = ds1.where(ds1[var_name] != 0, np.nan)\n",
    "                            ds2 = ds2.where(ds2[var_name] != 0, np.nan)\n",
    "                            \n",
    "                            ds1_seasonal = season_mean(ds1)\n",
    "                            ds2_seasonal = season_mean(ds2)\n",
    "                            ds_diff = ds1_seasonal - ds2_seasonal\n",
    "                           # create_seasonal_maps(ds1_seasonal, ds2_seasonal, ds_diff, var_name)\n",
    "                            create_distributions(ds1, ds2, \"Taylor\", \"OSA\")     \n",
    "                            \n",
    "                            counter+=1\n",
    "                            \n",
    "                        if counter==1:\n",
    "                            break"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}