extreme_precip_wap_profiles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Omega Profiles\n",
    "\n",
    "We plot the OMEGA profiles (`wap`) averaged over tropical ocean areas. \n",
    "\n",
    "- The overall mean is plotted.\n",
    "- The convective mean is plotted by selecting only points with OMEGA < 0.\n",
    "- Extreme precipitation points are plotted by selecting points where pr >= the 99th percentile.\n",
    "- Long-duration events are defined as events longer than 5 days as determined by our detection/tracking method.\n",
    "- Large events are defined as events with an average size of 12 gridpoints or more as determined by our detection/tracking method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sys\n",
    "from pathlib import Path\n",
    "from importlib import import_module\n",
    "import esmlab\n",
    "\n",
    "to_import = \"/Users/brianpm/Code/yellowstone/py\"\n",
    "current_path = sys.path\n",
    "if to_import not in current_path:\n",
    "    sys.path.append(to_import)\n",
    "import extreme_precipitation_functions as prf  # get_duration_and_size(events)\n",
    "import cmip_file_utils as cmip  # quick_load_var(loc, var, table, model, experiment, member)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how to reload:\n",
    "from importlib import reload\n",
    "reload(sys.modules['cmip_file_utils'])\n",
    "import cmip_file_utils as cmip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tropics_quantiles(da, threshold=None, q=0.99):\n",
    "    da_tropics = da.sel(lat=slice(-30, 30))  # SHOULD BE ADJUSTABLE\n",
    "    da_tropics_val = np.nanquantile(da_tropics, q)\n",
    "    if threshold is not None:\n",
    "        da_climo = da.mean(dim=\"time\")\n",
    "        da_process = np.where(\n",
    "            (da_climo >= threshold).sel(lat=slice(-30, 30)).broadcast_like(da_tropics),\n",
    "            da_tropics,\n",
    "            np.nan,\n",
    "        )\n",
    "        da_threshold_val = np.nanquantile(da_process, q)\n",
    "        print(\n",
    "            f\"Compare the tropics region ({da_tropics_val}) to the rainy tropics ({da_threshold_val})\"\n",
    "        )\n",
    "        return da_threshold_val\n",
    "    else:\n",
    "        print(f\"Quantile value: {da_tropics_val}\")\n",
    "        return da_tropics_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wap_profs(wap, pr, land):\n",
    "    assert 'lat' in wap.coords\n",
    "    weight = np.cos(np.radians(wap['lat']))\n",
    "    if land is not None:\n",
    "        wap_ocn = wap.where(land <= 0)\n",
    "        pr_ocn = pr.where(land <= 0)\n",
    "    else:\n",
    "        wap_ocn = wap\n",
    "        pr_ocn = pr\n",
    "    wap_mean = esmlab.weighted_mean(wap_ocn.mean(dim=('time','lon')), dim=['lat'], weights=weight)\n",
    "    wap_conv = esmlab.weighted_mean(wap_ocn.where(wap_ocn < 0).mean(dim=('time','lon')), dim=['lat'], weights=weight)\n",
    "    pr_crit = get_tropics_quantiles(pr_ocn, threshold=None, q=0.99)\n",
    "    wap_extreme = esmlab.weighted_mean(wap_ocn.where(pr_ocn >= pr_crit).mean(dim=('time','lon')), dim=['lat'], weights=weight)\n",
    "    return wap_mean, wap_conv, wap_extreme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "dloc = \"/Volumes/Glyph6TB/CMIP6\"\n",
    "# tup1 = (\"MRI-ESM2-0\", \"amip\", \"r1i1p1f1\")\n",
    "# tup2 = (\"MRI-ESM2-0\", \"amip-lwoff\", \"r1i1p1f1\")\n",
    "\n",
    "# tup1 = ('IPSL-CM6A-LR', 'amip', 'r1i1p1f1')\n",
    "# tup2 = ('IPSL-CM6A-LR', 'amip-lwoff', 'r1i1p1f1')\n",
    "tup1 = ('IPSL-CM6A-LR', 'aqua-control', 'r1i1p1f1')\n",
    "tup2 = ('IPSL-CM6A-LR', 'aqua-control-lwoff', 'r1i1p1f1')\n",
    "# tup1 = ('CNRM-CM6-1', 'amip', 'r1i1p1f2')\n",
    "# tup2 = ('CNRM-CM6-1', 'amip-lwoff', 'r1i1p1f2')\n",
    "# tup1 = ('CNRM-CM6-1', 'aqua-control', 'r1i1p1f2')\n",
    "# tup2 = ('CNRM-CM6-1', 'aqua-control-lwoff', 'r1i1p1f2')\n",
    "# tup1 = ('CESM2', 'amip', 'r1i1p1f1')\n",
    "# tup2 = ('CESM2', 'amip-lwoff', 'r1i1p1f1')\n",
    "# tup1 = ('CESM2', 'aqua-control', 'r1i1p1f1')\n",
    "# tup2 = ('CESM2', 'aqua-control-lwoff', 'r1i1p1f1')\n",
    "\n",
    "# load the events\n",
    "ev1 = xr.open_dataset(Path(dloc) / (\"pr_events_\"+\"_\".join(tup1)+\".nc\"))['precip_events']\n",
    "ev2 = xr.open_dataset(Path(dloc) / (\"pr_events_\"+\"_\".join(tup2)+\".nc\"))['precip_events']\n",
    "\n",
    "# What if we only consider events that are longer than 1-day\n",
    "dursiz1 = prf.get_duration_and_size(ev1)\n",
    "dursiz2 = prf.get_duration_and_size(ev2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile value: 0.00044758025964256433\n",
      "CPU times: user 14.7 s, sys: 2.74 s, total: 17.5 s\n",
      "Wall time: 23.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "wap1 = cmip.quick_load_var(dloc, 'wap', 'day', *tup1, compute=False)\n",
    "wap1_trop = wap1.sel(lat=slice(-30,30)).compute()\n",
    "if 'aqua' not in tup1[1]:\n",
    "    land1 = cmip.quick_load_var(dloc, \"sftlf\", \"fx\", *tup1, compute=False)\n",
    "    land1_trop = land1.sel(lat=slice(-30,30)).compute()\n",
    "else:\n",
    "    land1 = None\n",
    "    land1_trop = None\n",
    "pr1 = cmip.quick_load_var(dloc, 'pr', 'day', *tup1, compute=False)\n",
    "pr1_trop = pr1.sel(lat=slice(-30,30)).compute()\n",
    "wap1_mean, wap1_conv, wap1_extr = wap_profs(wap1_trop, pr1_trop, land1_trop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantile value: 0.0001830892966245307\n",
      "CPU times: user 14.3 s, sys: 2.3 s, total: 16.6 s\n",
      "Wall time: 22.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# now we do it for the second case\n",
    "wap2 = cmip.quick_load_var(dloc, 'wap', 'day', *tup2, compute=False)\n",
    "pr2 = cmip.quick_load_var(dloc, 'pr', 'day', *tup2)\n",
    "\n",
    "wap2_trop = wap2.sel(lat=slice(-30,30)).compute()\n",
    "pr2_trop = pr2.sel(lat=slice(-30,30)).compute()\n",
    "\n",
    "if 'aqua' not in tup2[1]:\n",
    "    land2 = cmip.quick_load_var(dloc, \"sftlf\", \"fx\", *tup2, compute=False)\n",
    "    land2_trop = land2.sel(lat=slice(-30,30)).compute()\n",
    "else:\n",
    "    land2 = None\n",
    "    land2_trop = None\n",
    "\n",
    "wap2_mean, wap2_conv, wap2_extr = wap_profs(wap2_trop, pr2_trop, land2_trop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "duration_crit = 5  # critical value for duration\n",
    "evlong1 = dursiz1.where(dursiz1['duration'] > duration_crit).dropna()\n",
    "# evlong1 only has event labels for longer events.\n",
    "# now we have to go and get the precipitation within those events\n",
    "evmask1 = np.in1d(ev1, evlong1.index.values).reshape(ev1.shape)\n",
    "evmask1 = xr.DataArray(evmask1, dims=ev1.dims, coords=ev1.coords)\n",
    "pr1_long_events = pr1_trop.where(evmask1)\n",
    "# now what to do with them?\n",
    "# pr1_long_events  # DataArray same size as pr1_trop & ev1\n",
    "wap1_long_events = wap1_trop.where(evmask1)\n",
    "wap1_long_events_mean = esmlab.weighted_mean(wap1_long_events.mean(dim=('time','lon')), dim=[\"lat\"], weights=np.cos(np.radians(wap1_trop['lat'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "duration_crit = 5  # critical value for duration\n",
    "evlong2 = dursiz2.where(dursiz2['duration'] > duration_crit).dropna()\n",
    "# evlong1 only has event labels for longer events.\n",
    "# now we have to go and get the precipitation within those events\n",
    "evmask2 = np.in1d(ev2, evlong2.index.values).reshape(ev2.shape)\n",
    "evmask2 = xr.DataArray(evmask2, dims=ev2.dims, coords=ev2.coords)\n",
    "pr2_long_events = pr2_trop.where(evmask2)\n",
    "# now what to do with them?\n",
    "# pr1_long_events  # DataArray same size as pr1_trop & ev1\n",
    "wap2_long_events = wap2_trop.where(evmask2)\n",
    "wap2_long_events_mean = esmlab.weighted_mean(wap2_long_events.mean(dim=('time','lon')), dim=[\"lat\"], weights=np.cos(np.radians(wap2_trop['lat'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_crit = 12  # critical value for duration\n",
    "evhuge1 = dursiz1.where(dursiz1['AvgSiz'] > duration_crit).dropna()\n",
    "evmask1 = np.in1d(ev1, evhuge1.index.values).reshape(ev1.shape)\n",
    "evmask1 = xr.DataArray(evmask1, dims=ev1.dims, coords=ev1.coords)\n",
    "pr1_huge_events = pr1_trop.where(evmask1)\n",
    "wap1_huge_events = wap1_trop.where(evmask1)\n",
    "wap1_huge_events_mean = esmlab.weighted_mean(wap1_huge_events.mean(dim=('time','lon')), dim=[\"lat\"], weights=np.cos(np.radians(wap1_trop['lat'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_crit = 12  # critical value for duration\n",
    "evhuge2 = dursiz2.where(dursiz2['AvgSiz'] > duration_crit).dropna()\n",
    "evmask2 = np.in1d(ev2, evhuge2.index.values).reshape(ev2.shape)\n",
    "evmask2 = xr.DataArray(evmask2, dims=ev1.dims, coords=ev2.coords)\n",
    "pr2_huge_events = pr2_trop.where(evmask2)\n",
    "wap2_huge_events = wap2_trop.where(evmask2)\n",
    "wap2_huge_events_mean = esmlab.weighted_mean(wap2_huge_events.mean(dim=('time','lon')), dim=[\"lat\"], weights=np.cos(np.radians(wap2_trop['lat'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAEcCAYAAAAIpvc+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAB8LElEQVR4nO2dd3hURdfAfye999ASIAm9EwihBkFRURBEQVTsvfeCXV/r+9n1tWEBG2LDhh0VE1pCgJBA6CShBQLpvezO98e9gQApS3rC/J5nn+zemTtz7t3Nnp0zp4hSCo1Go9FoNK0Tu5YWQKPRaDQaTc1oRa3RaDQaTStGK2qNRqPRaFoxWlFrNBqNRtOK0Ypao9FoNJpWjFbUGo1Go9G0YrSi1miaGBEJERElIg4tLYtGo2l7tGtFLSKpIjLJfH6ViFhEpEBE8kQkQUSmVun7sIikmO17ReTLKm3LROQ6G+cUEblDRDaKSKE51tciMshsX2B+aU877rzXzONXVTkWKCILRSRHRLJF5PNq5lsgIhUi0qUOua4SkeU1tC0TkRLz2g+LyGIR6WzL9Z4KVP0caTQaTXPTrhV1NaxSSnkAPsCHwFci4iciVwKXA5PM9gjgr3rO8TpwJ3AH4Af0Br4HplTpsw24svKFudKaBew8bqzFwAGgO9ABeKlqo4i4AxcCucCcespbyW3mtfcEPI6fS1MzeqWs0WiaklNNUQOglLICHwGuQBgwAvhdKbXTbD+glJp3suOKSC/gVuASpdTfSqlSpVSRUupzpdQLVbr+BIwVEV/z9WQgEUMpV451FtAVuF8plauUKldKrT9uyguBHOA/VFH8DUEplYPxw2KoreeIyFwR2Ski+SKSLCIzqrTZi8hL5kp9l4jcWtUMfPxqVUSeFJHPqrz+WkQOiEiuiESLyIA6ZJluWkvyTJkmm8e7iMiPIpIlIjtE5Prj5vxKRD4xr2GTiESYbZ8C3YCfTIvDA1VM2deKyG7gbxGxE5FHRSRNRDLMsbxtvYcajUZTE6ekojaVxHVAAbAdWA1cISL3i0iEiNjXc+gzgL1Kqbg6+pUAPwIXm6+vAD45rs8oYCvwsYhkisgaETntuD5XAl8Ai4C+IjKsnnIfQUT8gQuAHSdx2k4gCvAGngI+q2I6vx6YCoRjWCpmnqRIvwK9MCwK64ATzP9VZI/EuI/3Y1hNxgOpZvMXwF6giynDcyJyRpXTp2HcRx+M9+Z/AEqpy4HdwHlKKQ+l1P9VOec0oB9wNnCV+ZiI8ePPo3IMjUajaQinmqIeJSI5GCvXS4AZ5mr1M+B2jC/cf4EMEZlbj/H9gXQb+36C8ePAG+ML//vj2oOBs4B/gE7Ay8APIhIAICLdMJTCQqXUQQxTfUNW1W+ISC5wGAjAuB82oZT6Wim1XyllVUp9ifHjJ9Jsvgh4TSm1RymVBTx/MkIppT5SSuUrpUqBJ4EhtaxUrwU+Ukr9acqyTym1RUS6AuOAB5VSJUqpBOADjO2OSpYrpX5RSlmAT4EhNoj3pFKqUClVjLH18IpSapdSqgB4CLhYm8U1Gk1DOdUU9WqllI9SKkApNUoptbSywTRPT8JYUd0E/EdEzq5tMNNEWmA+ooBMwCYnLKXUciAQeBRYYn7ZV6UYSFVKfWiavRcBe4CxZvvlwGZT6YCx0rxURBxFJKqKXJtskQe4QynlDQwGfDF+KNiEiFxhmptzzB9CAzGUPRgr2D1VuqedxLj2IvKCacLO4+jqOEBEulW5xgLzeFdO3OevlCFLKZV/nBxBVV4fqPK8CHCxQclWva4uHHttaYAD0LGOMTQajaZWTjVFXSemUvwaY894YB19B5jmUA+lVAzGqja4cn/TBj4D7uVEszfm/LWVNrsCCDP3bw8Ar2Aox3OUUjFV5Kp1T/d4lFJJwDPAWyIidfUXke7A+8BtgL9SygfYCFSem46hQCvpdtwQhYBbldedqjy/FJgOTMIwq4dUTquU2l3lGj3M43uAHtWIuR/wExHP4+TYV9f1mdT0PlQ9vh/D6a/q+BXAQRvn0Gg0mmrRipojoUtTRMTTdAo6BxgAxFbp5iAiLlUejsePo5TaDrwNfCEiE0TEyex7cQ2m9DeAM4Hoatq+A3xF5EpzZTkTYwW4QkRGYyikSAynr6EYPyoWUrv5W467Bpca+n2MsSc8zTxpgojUpKzcMRTWIbPv1Rz7A+cr4A4RCRbDee74+5CAYSJ2NH/gVN3D9gRKMSwVbsBztVwbGJ78V4vIGeb7GCQifZVSe4CVwPPmdQ/GMJPXuN99HAcx9p1r4wvgbhEJFREPU9YvlVIVNs6h0Wg01aIVtUEe8DCG01AO8H/AzaZ5upJ3MMzRlY/5NYx1B4YT0VvmWDuBGRie3seglMpSSv2lqikKbu7nTgPuwwi/mgtMV0odxlDGPyilkkwP9QNKqQMYoWFTRcSvBtnGHHcNxdWZd5VSZRg/Ih4zD3UFVlU3oFIqGWP/fBWGQhsErKjS5X3gd2ADhjPY4uOGeAzjR0c2hiPawiptn2CYkPcByRhOfzViOvFdDbyKcc/+5egq9xKMFfl+jB9BTyil/qxtvCo8Dzxqmvbvq6HPRxh729FACobDoM37/BqNRlMTUo2O0GiOQUQ+AL5WSv3eCGOFYCgyR73a1Gg0mrrRilrTrGhFrdFoNCeHNn1rNBqNRtOK0StqjUaj0WhaMXpFrdFoNBpNK0Yrao1Go9FoWjHtNr3h5MmT1W+//dbSYmg0Gk1bo85ER5rmpd2uqA8fPtzSImg0Go1G02DaraLWaDQajaY9oBW1RqPRaDStmDajqEVksohsFZEd9SxBqdFoNBpNm6NNKGoRscfInX0O0B+4RET6t6xUGo1Go9E0PW1CUWNUidqhlNplFoxYhFH+UKPRaDSadk1bUdRBGLWGK9lrHjsGEblBROJFJP7QoUPNJpxGo9FoNE1FW1HU1cX1VVcacp5SKkIpFREYGNgMYmk0Gk3Lk7X7EH9MfIo1/9O5I9ojbUVR78WoiVxJMEZdYY1Goznl2b1qG4NjZ2J9ZXtLi6JpAtqKol4D9BKRUBFxAi4GfmxhmTQajaZVsGPfPlICOvCnn04q1h5pEylElVIVInIb8DtgD3yklNrUwmJpNBpNq8DNVYjxmUJPt/CWFkXTBLQJRQ2glPoF+KWl5dBoNJrWRnl2NvcnxZPWpVtLi6JpAtqMotZoNBpN9ew5VEKFrGWD92ZCWloYTaPTVvaoNRqNRlMDeeWZnDfmW3Z0W9/SomiaAK2oNRqNpo3jUVrGLyufZ3pKakuLomkCtOlbo9Fo2jjZuVas6nZKLL1aWhRNE6BX1BqNRtPGiU/fwv8ir+eP7gUtLYqmCdArao1Go2njHMpM57S8ydj5RLa0KJomQCtqjUajacOUlpaStDmBwfiQUWhpaXE0TYA2fWs0Gk0bZtOmTZQDe712Eh2iiwq2R7Si1mg0mjbMunXrwBGeG345nfiipcXRNAHa9K3RaDRtmPXr1+Ps7Mzby1ZipRylFCI653d7Qq+oNRqNpg2zbt06fL19KZdviQ9+jJzcopYWSdPIaEWt0Wg0bRSLxUJiYiIWZeGWSVvY0GM7WSl7W1osTSOjFbVGo9G0UbZt20ZRURFz587lqU3zuXb5B+Sl7mlpsTSNjN6j1mg0mjbKunXrACgoKCDg4AwOuXSgaPf+FpZK09hoRa3RaDRtlCWf/MQjTuez+OGH6dU/ms2dg5hgLW1psTSNjDZ9azQaTVvCCnxzkKywx/ngr3Rucp7GGsCh28uEjDrM2DuubGkJNY2MXlFrNBpNWyClDB79Bet3H2FX/At+WPi3+60EFp3Pqpsg0SWFKdPHYmev11/tDf2OajQaTWulDHg1CRV6NyosGBbOoKI0lrigv8l18MXTZyWD8rvhd+0YlL09BQW6KEd7RK+oNRqNprWxKgse/QIVPR+pWAs4stPvdL4bPRavqDDC3xzAF4PvY8Yn03lj2TK6desGwP792pGsPaIVtUaj0bQGci3w5J/w6XzI/B4oo8B5ID8OuJclF5zGVYUR3LT4ESxLnyXY1ZON367BPziYWwcMQCkFQEBAQItegqZp0KZvjUajaSkU8OU2GPow+HaH187Bmv0X8cGzuXT2O9z+zoeEfP0gC8/pRHm8hZU+VzCmtJS5997JypUrCQ4OPpIy1N7eHqvV2tJXpGkC9Ipao9Fomptd+TD3K1gyH4pXoLAj2+8MPhx4MwtGDeK8yBCeGheK/w4nrOf/gCTPJuqfGEbechsVPXpw3333sWbNGq699tojeb2tVit79+qsZO0Rrag1Go2mOSi3wsvR8M582P0NUITFuQ+rh9zPw5FDye8ayGWjurNsVHcCPZ1h1y6W/liCZ2Ekbt//zvtrVrF18yaWLFmCi4sLUVFRREVFHRneyckJZ2fnlrs+TZOhFbVGo9E0JTFp8OjHsGIBWFJQ4kVW2CzeiBjHxyGdCAv04IaoUC4cFoyLoz0VeVb2jD9I8Kp+THjkIUq3PkpOphNPXXYh5513HlOmTGHevHmsWrWK+fPnH5nG3d0dV1fXlrtOTZOhFbVGo9E0NtnF8Ohi+GI+ZP8NKCx+pxM/9i7u7BfCAbFnVJgfH4wL4/S+HbCzM8tSZmWRFr2XnJ1dKXpsEX3uOgMHZwduuP9+ysvLefXVVwHIy8vD29v7mCnLy8vJyspq5gvVNAdaUWs0Gk1joBR8Egv/Nx+SFwF54BBK/phHeO/MsbxTZjh6TR3cmevGhTEo+FhFm3ZzIR1/m0mPTsXkbo/G220GANHR0SxcuJDHHnuMHj16AHDfffedMH1ZWRmFhYVNe42aFkErao1Go2kIm9PhoU/h9wVQshlwRYXMZNvMmTzR3Z/Ve3PwVMK140K5akwIXXyOM09XVGDdtJkd0R1IGzWX8U+H4e3maDZVcNttt9G9e3fmzp0LwHvvvYebmxuXX375McN06tRJe323U7Si1mg0mpOlpAz++xPMmw/7fwMs4DaGitnv8cNFY3lzZwapmUUEFZTy2NT+XBQRjKeL4wnDpH9YjOuLb+KzbS6Ru/fh3uVMqDSDA2+//TZJSUksXrwYNzc3ABYvXkx4ePgJY4kIJSUlTXbJmpZDK2qNRqOxlT8S4Mn5EPs5WDNBukD4/eTccykfBDjxWVwaOXGpDOnqw//O7sPkAZ1wqCn39vbtrF1qh6PnGZy1bBmewZ2PaT548CCPPfYYZ599Nueff/6R47/88suRkKyqHD58mLKyska8WE1rQStqjUajqY0DmfDo5/D1fMhLAJwgcDpceTXbbhzN++vT+CFhL+VWK2f268j148OI6O5brTIFKNhYTtHFm+iwKZyzFnyCmnMp4mB/Qr+5c+dSXFzM66+/fmSsL774goEDBzJo0KAT+nfu3Jn09PTGvHJNK0Erao1Gozmeigr44Hd4ZT5s/xEoB8dhcNabqKcvZrmv4v2YFKI/WIGLox2zR3TlmnGhhAa41z7uoUNs/PsQ5YWdsHt7EQEXnw/VKOlVq1axYMECHnzwQfr06XPk+C233MKjjz5araL28vLSub7bKc2uqEWkK/AJ0Amjsuo8pdTrIuIHfAmEAKnARUqpbPOch4BrAQtwh1Lq9+aWW6PRnAKs2wqPzoeln0B5OhAAPW6Fu66m9LoB/LhxPx8uT2bLgXwCPZ25/+w+XBrZDV93p9rHrYB9kwvpuGEkkWeNoXDzfDxdZlfb1WKxcNtttxEUFMSjjz56TFtCQgLBwcHVnnfw4EGKi4vrc9WaVk5LrKgrgHuVUutExBNYKyJ/AlcBfymlXhCRucBc4EER6Q9cDAwAugBLRaS3UsrSArJrNJr2Rm4ePP8lzJ8PGasAe3A/Fy65Gp6ZQo4/fB67m49f/oeM/FL6dvLkxZmDmTa0C87VrIZPoKSE7IRUduzy48CVLzH8zhHVOpZV8v7777Nu3ToWLVqEh4fHkeP//vsv4eHh2NtXP6eLi4v2+m6nNLuiVkqlA+nm83wR2QwEAdOBCWa3j4FlwIPm8UVKqVIgRUR2AJHAquaVXKPRtBusVvh+GTw3H9Z9C6oYpD9EvAiPXQbndSI1s5CPVmzj6/i9FJdbiOoVwEuzhhDVK6DG/efj2XNfEb5L7sX38NcMTd2Ll/sMqOXczMxMHnnkESZOnMhFF1105PjBgweZMGEC33zzDRdeeGG153bt2pUdO3ZgsVhqVOaatkmL7lGLSAgQDsQCHU0ljlIqXUQ6mN2CgNVVTttrHqtuvBuAG4Aj9Vk1Go3mCDtS4ImP4fuPoSgV8IYOV8C1V8PcSJQnxKdl88Fn8fyRfBBHOzumD+3CdVFh9Onkafs8SsGOHaxf44Zjn5mc88nVeHu41Hnaww8/TG5uLm+++eYxPwYCAwP56aefmDx5co3nVpq9KyoqtKJuZ7SYohYRD+Bb4C6lVF4tv1Cra1DVdVRKzQPmAURERFTbR6PRnGIUFcG738Ib8yHtH0DA8Qw451l4egYMd6XCYuW3Tem8H5PChj05+Lg5cuuEnlwxpjsdPOtWsFU59HUJ9g//gd+O6ZwZvRynsacfExtdE/Hx8bz//vvcddddDBgw4Ji2ffv2MWXKlFpX8pXpQ8vLy3VxjnZGiyhqEXHEUNKfK6UWm4cPikhnczXdGcgwj+8FulY5PRjQro0ajaZmlIJlq+A/8yHmS7DkA2HQ8z9wz5VwbTdwgoLSCr5cnsL8FSnszS4mxN+Np6cP4MLhwbg51ePr8eBB1q7MpcK9N2d9+TWuY0fbpKStViu33XYbHTp04MknnzymLS4ujpEjR7Jx48YTFHhVevfuzbZt2ygvLz95uTWtmpbw+hbgQ2CzUuqVKk0/AlcCL5h/f6hyfKGIvILhTNYLiGs+iTUaTZth/3544RP4dAHkbAXcwH0WXHQ1PBEF3Y3kI/tzilmwNJUvYneTX1pBZIgfj0/tzxn9OmJvg2I9nrLdVrKnZNNhRw8m3Xc3FWuexMmxr83nL1iwgNjYWD755BO8vLyOaQsLC+Oll16id+/etY7h5GR4nmtF3f5oiRX1WOByIElEEsxjD2Mo6K9E5FpgNzALQCm1SUS+ApIxPMZv1R7fGo3mCKWl8MWP8OJ8SP4dsIKMg8gH4NFZMMUTzORgG/fl8n7MLn5OTEcB5wzsxPVRYQzp6lP/+YuK2LlmL/uLvAh74h1CrzsHB0fb94izs7OZO3cu48aN47LLLjuh3dfXl3vvvbfOcXbv3g2gs5O1Q1rC63s51e87A5xRwznPAs82mVAajaZtoRSsXQ/PzodfFkJZFhAEHebCDVfB3b3Az+hqtSr+2ZzB+zG7WL0rCw9nB64cE8LVY0MI9nVrgAywf0oxfjsupK/3YYI3rcTT5fK6zzuOxx9/nMzMzBMcyACWLFnCVVddRXp6Oo6ONYd0geH1HR8fj8Wi1zHtDZ2ZTKPRtB0OHYK3P4d358OBRMAZHM+HKVfD45NghP2RZUBJuYXF6/bxwfJd7DpUSBdvFx45tx+zI7viVUscs01YrZTv3M36vU64nX4jE28JqzU2uiY2bNjA22+/zc0338zQoUNPaA8MDGTOnDk4ONT9Ve3r6wto03d7RCtqjUbTuqmogJ9+hf/Oh7gloMqBCOj1FtxzMVzhB1UWxocLSvlkVRqfrU4jq7CMQUHevH7xUM4d1BnHmgpknATpTxXjuug9fLbczYTd+3ALnl5rbHRNKKW47bbb8PPz4+mnn662z8iRIxk5cqRN4+Xm5gJQWlp60rJoWjdaUWs0mtZJcjK8PB8WfQpFB4FAcL8dLr4K5g6Cnsd235GRzwcxKSxev49yi5Uz+nbkuqhQRob62ZygpE4OHGD1xiIsXSOY+X8/4t61S72H+vzzz1m+fDkffPDBkdVwVRYtWsRPP/3E559/flLj6j3q9odW1BqNpnWgFGzbDl8uhQWfQEosYA8yBSKvhoemwFTHY761covKiU3J5Iu43fyz9RDODnbMGh7MteNCCQv0qHGqkyVvaTkV123Hd88gprz/IXZXXgENWJ3n5eVx//33ExkZydVXX11tn8OHDx/x5LaFrl2NKFaldAqJ9oZW1BqNpmWwWGB9InwdA3/EwJYYKDloNg6AwJfgpsvg1o7Q0TiakV/CmpRs4lIyiU3JYuvBfJSCAA8n7jmzN5eN6o5fXQUyToZSICGOjSvtyPDyZ+Lr7+E9e1aDlDTADz/8wIEDB/j666+xs6t+rNtuu61eY+s96vaHVtQajaZ5KC2FFWvgqxj4JwZ2rYCKPLOxG7ieCRFRcE4UzOwLg4S9OUXEpewlbnkWcSlZ7DpcCICbkz3Du/syZVBnIkP9GNrNx7YCGSdB8T8VOJ9RijiMZ9R991C29mlcHK9rlLFTUlIAGDFiRLXtubm5vPTSSzzxxBM2OZIBR0pcFhUVNYqMmtaDVtQajaZpyM+Hpavgq2hYEQN7Y0FVOjr1A59LIDwKzouCGd1Q3RS7MguJS8kibvMG4n7JYl+Okb/ay8WByFA/Lo7sSmSoPwO6eDWKY9gJKOBPYM2rZDl48+W5Y5l653J6jx+Ey0nERtdFamoqXbp0qTHVZ3x8PC+88AI33XQTQUHVljY4gcr6Bq6uro0mp6Z1oBW1RqNpHA4dgiXL4dtoWBMDGesxSs7bA+HQ8RaIHA8XjoNzA7D4K7YeyCcuJZO4FWuJ+yyLwwWGI1SAhzMjQ/24YXwYkaF+9OnoiV09MoadLJb1CvuzFYXBnxB09+Xcs6RPk8yTmppKSEhIje0TJ05k9+7ddO7c2eYxK1feeo+6/aEVtUajqR+pafBtDPwUAwnRkLvFbHABu5HQ/WEYNx4uGgWne1LuamXjvlxjxfxzCmtSs8grqQAgyMeV8b0CiQz1IzLUj9AA98bz1K4LK/A4sPxCrMFuPPLU/Uy79E/G9AxosinT0tJqDbvKyso66euvNH2XlJQ0SDZN60Mrao1GUzdKwcbN8GUM/B4Dm6KheI/Z6A0OY6HvlXB6FMyOgFHOlIiFhD05xKUcIG5RMmvTsikuN7JmhQW6M2Wwsb88IsSvYRnCGkJGBiolj+x3u5I6pD/DLhvLC5MHN+mUFouF3bt3M3v27Br73HLLLRQVFbFkyRKbx/Xx8QHQJS7bIVpRazSaE6mogNXrYVEM/B0DO5ZD+WGzsSO4jofh98PZUTB7EAy0p6C8grVp2cSlpBD3URYb9uRSZrEiAn07eTF7RNcjijnQs4XLMJaDOltB6izoWMaCLz4mvNvt0KdDk0+9f/9+Kioq6N69e419brnllpNeUXt6GvWym80SoWk2tKLWaDRQXAx/x8GiaFgeA3tWgaXAbAwD7ykwZLzh+DWzJ4QIWYVlrEnNIi5lK3H/ZrFpfy5WBQ52wsAgb64eG0JkqB8R3f3wdmtgys7GIgtYvRx2JbHceSYxY27n4WdHcE8tSrOxSUtLA6h1j3rChAknPW5BgfF+6Vzf7Q+tqDWaU5GcHFiyEr6JhtgYOLjGTM0pwEDocAWMGA8XjINpQRAAB/NKiE3JIi5hI3HfZbHtoKEYnB3sCO/mw22n9yIyxI/wbj64O7fCrxYFKkJRUfE2jj3Ssb57IePkNOge2KxipKamAjUr6oKCAiZNmsRbb73F8OHDbR43Pz8f0HHU7ZFW+N+k0WganfQDRmKRH2NgXTRkJ2LEIjmAXQR0vQvGjoeLxsKZvig3xZ6sYmJTMon7ZwNxqVmkZRrxue5O9kSE+DF9aBAjQ/0YFOzd6DHMjcoK4NsPIf5jvrv+Y17fdz2/vBzFaS30Y6JSUVeGUx2PiFBRUYGLi8tJjRscHAxQZ5UtTdtDK2qNpr2hFGzfBQtj4DfT8atgh9noBo6joc8TMCEKLh4JY91RDoodGQXGivnHVOJSsjiQZ3gP+7g5Ehnix+WjujMy1J9+nT1xaIoY5sbGYoH0dKx3BJFfDB4XnMbImzvxSm7LrvjT0tLo2LFjjfHO7u7uxMfHn/S4lWFZOjyr/aEVtUbT1rFaYe1GQzH/FQPbYqB0v9noC65REH4jnBUFlwyDQY5YUGxOzyM2JYO4RZmsSc0mq9CIYe7g6czIMH8iQ/0YGepHz0CPZolhblTmA+89B+ueZulvm7j5t878dMss+vu4EuTTsglB6oqh3r59OytXruTKK688qXF37doFQHFxcUPE07RCtKLWaNoaZWXw7zr4IhqiYyBtBVRkm41B4DXe2F+eEgUX94cQO8oqrCTtyyE2JY24j7NYm5pNfqkRw9zNz43T+3Y4opi7+bm1Xc/hQ4fgq69QKTewLWAqPl+GM2F8D/4ZEkQ3/xYKATuO1NRUhg0bVmP7l19+ybx5805aUYeEhCAi2vTdDtGKWqNp7RQWwpLVxh7z6mjYvxpU5aqpNwReABHjYUYUzAiBAKG4zML63dnEbt9B3B9ZrNudTWmFFYBeHTyYNrTLkeQinb3bScrJJ4DXFkHFXPbETuacz/fzXp8RnOlg12qUtNVqZffu3VxwwQU19rn77ru56qqrTnpsBwcHHUPdTtGKWqNpbWRmwVfL4YcYWBsNh9cBFYAd2A2B4OthVBRcFAXndAR3yCspZ21qNrFrthKXkkni3lwqrAo7gf5dvJgzsrsZw+yLv0cLxzA3JgXAD5vh4zvhinn8dtMFyKypnD0wlK3PhDZNPvAGcODAAcrKymqNoS4sLCQg4OSzoh04cICKigrt9d0O0Ypao2lpUvfC5zHwawwkRUPeJrPBCRwioff9MD4KLh0D47zBETILSlmTmkXs35uIS8kiOT0PpcDRXhgc7MP1Zo7s4d198XJpp6bQ4mL4Pxd4PgNOd6TsLA9ez9zM1YQAtDolDbbFUI8ePZqrr76aRx999KTGdnJyartbFppa0Ypao2lOlIKk7fBpNCyNga0xUJxiNnqAy1gYcglMioLLImGICwjszyk2FPOS3cSlZLEjw4hhdnG0Y1g3X+48oxeRoX6Ed/XF1amdmz/3Ao+shs/Gwpe/M//dHuzp8CKPdwjg1zujWlq6WqkrhhrggQceICrq5K/Dz88PJyenGutba9ouWlFrNE2JxQLLE+HzaPg3BlJioDzDbAwArygIvwPOjYI5QyDEAaUUqZlFxKVkEPd1NnGpmezJMvakPZ0diAjx5cJhwUSG+jEoyBsnh1Pki1kp+PFHSBsLf/SCp1+HqePY8+cuQgPdW1o6m6hU1LWZvm+88cZ6jV1cXIzFYtGZydohdSpqEQkGLgaigC5AMbAR+Bn4VSllbVIJNZq2RGkp/LIGvoyGVTGwdyVY88zG7hBwFgwbD+dHwUV9wF+wWhXbMvKJS9lL7Mos4lKyOJRv1G32d3diRIgfV48JJTLUj36dvbBva6FSjcEmYOYG2HI+fPUVby0M5/v1g/jdyZnHz+vf0tLZTFpaGgEBAbi7V//DYuvWrVx//fX8+uuvNfapifT0dCoqKigrK2sMUTWtiFoVtYjMB4KAJcB/gQzABegNTAYeEZG5SqnophZUo2mV5OXDVyvhuxiIj4ZDcaAMJYv0hy6XGI5fs6LgvG7gBhUWK5v25xG3KYXYlCzWpGaRW2w4AHX2dmFsD38iQ4045h6BzVjusTVyoBz+cxeEX4gaNoH8j7biNaoXYRsPcNXYkDYX311XDHVeXh7Z2dn1WhWHhobi6uqqw7PaIXWtqF9WSm2s5vhGYLGIOAHV58HTaNoj+zPgs+XwcwwkRkNOAkZBY3twGAY9boWoKLh0HEwIAAcoKbeQuDeXuNXbiU3JYl1aNoVlxhdxaIA7kwd0OhIqFezremor5qokAOFZ0C0GLpjBl0P28tB324kd2I1zBnVuaenqRWpqKgMHDqyxfcSIESQlJdV7fBHRmcnaIbUq6hqUdNX2MmBHbX00mjbNljT4OBr+jIHNMVC0xWxwAZdRMPgROCMKLhsN4R4gUFhawbrd2cT9vZXYlCwS9uRQZsYw9+3kyYXDjf3lyBA/OnidXD7nU4Ivi+DhcLjsctSHj5A8PobeoZ6MzC7mtdlD6eDZNu+ZUoq0tDSmTp1aY5/Vq1djsVgYO3bsSY+/a9cuioqKtOm7HWKTM5mI9AKeB/pjmL4BUEqFNZFcGk3zY7XC6s3wWQwsi4GdMVC2x2z0Bs+xMOoqOCcKLhsOYUY8cm5RuVHu8dc9xKZksXFfLharwt5OGNjFiytHdycy1J+I7r74uju12OW1erZsgSw3uDgIps+Biy9klVMml74fy6fXRhLVK5DQgLbhNFYdGRkZlJSU1OpI9uyzz+Lo6FgvRR0UFISzs7NOetIOsdXrez5G3p9XgYnA1Rj18DSatktFBfy23kjFuTIG9iwHS6bZ2An8oyD8AZgWBZcMhADjCzAjv4Q1KVnE/ZBJbEoWWw/moxQ42dsxtKsPN5/Wg8hQP4Z198WjNZZ7bG0o4F4LvD8W5syCfe/yd+6N9O7kyXBPZxZeN5JRYf4tLWWDsSU065133qm3onVxcTnpiluatoGt3yKuSqm/RESUUmnAkyISg6G8NZq2QVExfB0Li2MgLhoOrgJVaLRJD+h8HkRGwczxMKMHuBm/RfdmFxGXkk7cv4ZH9q7DxjluTvYM7+7LlEGdiQz1Y0hXH1wc9WrmpPj0U6hwgOhL4Ml/4PoQDkkp17wRz1PTBnDlmBDG9Dz5LF2tEVuSnbi7u+Pj41Ov8Q8fPkxhYaEOz2qH2KqoS0TEDtguIrcB+4AOTSeWRtMIHMqBT1fAkhjYEA1Z8UA5IOAwCEKvgnFRcEkUTOoCDsY+4q7DhcRt3ENciqGY9+UYMcxeLg5EhvpxcWRXIkP9GdDFq1Vmv2oTWIHTgB1vweXjIf4Sfkjwxyu9hIl9vFg593Q6e7ev1WFdMdQHDhwgJCSEX375hdNPP/2kxy8vL8dqtWK16ojZ9oativouwA24A3gaw/x9cqVdNJqmZmc6fBwDv8dAcgwUJGLYVR3AOQIG3g0To+DysRDhCwIWq2LrgXziYlOISzUU8+ECwxknwMOZkaF+3GCm4+zT0bPNhQO1OtKBh2+HksNwzhcwcTWMNn4gvfj7VqYN6cLEPh3o0sKlKJuCtLQ0fH198fLyqrY9ICCAW2+9lSFDhtRr/M6dO+Pr66szk7VD6oqj7gA8DPQEkoDnlVJXN8bEImIPxAP7lFJTRcQP+BIIAVKBi5RS2Wbfh4BrAQtwh1Lq98aQQdOGUQrW7oJPouGfGNgeA6WVAQhu4DEaIp+Es6PgypHQw6ieVG6xsnFfLnHRO4kzY5jzSoxyj0E+rozvFXgkVCo04BSPYW5M8vMhOxumdoPyEHh0FMyBxev2suv3Qu47uw8xD0xs1/e7rhhqBwcHXn755XqPX1ZWRmlpqTZ9t0PqWlF/AqwF3gSmAm8AVzXS3HcCm4HKn5dzgb+UUi+IyFzz9YMi0h8jM9oAjMxoS0Wkt1JKfxpPJaxW+HMjLIyG5TGwOwYq0s1GP/AbB2NuhPPGw2XhEGgkfSgpt7B+dw5xS/cRl5rJurQcisuNj05YoDtTBnc2q0r5EezbOkohtjuigXuvhPxk+GILBN4LwUZTfFo2Lg7Gvn57VtJgKOrevXvX2P7rr7+yZMkS3nrrrXqPX1BQQEVFRX1F1LRS6lLUnZRSj5jPfxeRdY0xqZmWdArwLHCPeXg6MMF8/jGwDHjQPL5IKVUKpIjIDiASWNUYsmhaKSVl8O1a+DYGYqMhfQWoHKNNgqDTBBgRBReMh1n9wM0w9+WXlLM2LZu4tYYZe8PeHMotChHo28mL2SO6HlHMgZ7tqNxjayQpCVatgu9uAKen4Fs7GABLEvfz8ZJUPrtuJM/NGNTSUjYLlTHUZ511Vo19UlJSSE5OrvccPXv2xNfXV4dntUPqUtQiIr4cDcWyr/paKZVVz3lfAx4APKsc66iUSjfHTTfN7mCkMF1dpd9e81h1wt4A3ADQrZtOmNamyCmET1fBjzGwPhoyYzHSygP2fSBkJoyJgoujYHIIOBgfyazCMtakZBxx/Nq0PxerAgc7YWCQN9eMNXJkR3T3w9tNp1ZsNhYAT70JLjGQdK3hvGdSYVEMCfbB6RRyxMvMzKSwsLDWGOpbbrmFW265pd5z2NnZ4ejoqDOTtUPqUtTeGKbvqjapylW1Ak464YmITAUylFJrRWSCLadUc6zaT6JSah4wDyAiIkJ/WlszuzNhwXL4LQY2xUDeWgwXBDtwHgL9b4DTouDycTCq45FPwcG8EmI3pROXkklcShbbDhrlHp0d7Ajv5sNtp/ciMsSP8G4+uOsY5uZFAd/8CO+8BpOXwKx34QXAzo6Y7Ye46dO1/HT7OM4PD+L88Gp/a7dbbImhXrp0KX369KFr1671niMjIwNPT8+6O2vaFHWlEA1pgjnHAtNE5FyMLGdeIvIZcFBEOpur6c4YBUDAWEFX/eQGA/ubQC5NU7Jhj+GR/Zfp+FW8yWxwAvdIiHgQzoqCK8dAb8NtQSnFnqxiYteaoVKpWaRlFgHg7mRPRIgf04cGMTLUj0HB3jg7aJNfi6AUZGXBe/7wVAFc3hvucwE7OzILSvF2daSztws3jO9BN79T0w/AlhjqSy65hHvuuYeHHnqoXnP4+fnh5ubW7vf6T0VsXnKISBDQveo59amapZR6CHjIHHMCcJ9S6jIReREj5OsF8+8P5ik/AgtF5BUMZ7JeQNzJzqtpRpSCf7bBwhiIjobUGChPNRs9wWcMjLgEpo6HK0ZARxfzNMWOjAJiV6cdMWUfyCsBwMfNkcgQPy4f1Z2Rof706+yJwylkOm21FAAPzIN3boLf9sNnl8KsSwHYdjCfs1+L5o2LwzlvSBfunHTqrvRsqUP966+/0qtXr3rP4eXlha+vrzZ9t0NszfX9X2A2kIxhnwTD0NWY5S1fAL4SkWuB3cAsAKXUJhH5ypy7ArhVe3y3Msot8N0G+DoGVkfD/uVgNQ0iEggdoiDiTpgxHmYPBg/jY2exKjan5xG73DBlr0nNJqvQiGHu4OnMyDCj1OPIUD96BnroGObWREEBfPMNHLoE3jsT3v0KzjYqWm09kI+vuyPd/Nz474WDmdSvYwsL2/Kkpqbi5eVVa9axgQMHNigFaG5uLtnZ2XTsqO93e8PWFfX5QB/T87rRUEotw/DuRimVCZxRQ79nMTzENa2BvBJYuAZ+iIG10XBoJZBvtNl1h65nw+gouCgKzutzxPGrrMJK0r4cYuON1fLa1GzyS41Qkm5+bpzet8MRxdzNT5vwWi37gcu/g2XXQvxoSO0DXQ13laKyCs7733JuiArjvrP7cFFE/fZb2xtpaWmEhITU+Jn++++/Oeecczh48GC9U4hmZ2dTXFys46jbIbYq6l2AI9CoilrTRtifBwtWwi8xsDEGcuM48lFw6g9958D4KLgsCsZ1PeL4VVxmYX2qUbgiLiWLdbuzKTXLPfbq4MG0oV2OJBfp7N3+MlG1O/bvhzvugEnPwJ6LYf0UGOwHwNq0LCxWiAz14+sbR9O746lr5q6OupKdhISEcNNNN+Hh4VHvOUJCQggODtaZydohdWUmexPDxF0EJIjIX1RR1kqpO5pWPE2LkJwBC2JgaQxsjYGiBIzkzPbgNgzCb4Mzo+CKsTDgaMGEvJJy1m7NMBVzJol7c6mwKuwE+nfxYs7I7mYMsy/+HjqGuU2xDThjD/ikw8uucJMj4Hekee63SfTp5HmkOInmKEopUlNTmTBhQo19wsLCeP311xs0j8Vioby8XCvqdkhdK+p48+9aDKcuTXvDqmBFGnwaDdExkBIDZVvNRhfwHgXjHoFzx8OVo6DL0V/8hwtKWZOUfiRHdnJ6HkqBo70wONiH680c2cO7++LlomOY2yT/bIUp4fDKpzD1Qnh6BZi/zWJ3GT/Grh8fxjc3j8HNSXvdV0dOTg75+fm1rqjfffddnJycuOaaa+o9z/bt2zlw4ABBQadW6NupQF3hWR83lyCaZqLCCj9thq9iYGU07IsBy16z0RsCx8Hwq2F6FMyJAE+nI6fuzykmbv2+IyvmnYeMco8ujnYM6+bLnWf0IjLUj/CuvrjqL+22TXQ0uPSD0zvBjEfh4jPgpmO7LF63j305xVwXFYq3q/4hVhO2eHwvW7YMf3//BinqsLAwOnTooH072iF1mb5/wkgg8ptSqvy4tjCMvN+pSqmPmkxCTcMoLIcv1sP3MRAfAxkxUJlQzq4zBEXBqCiYGQXnDwRTwSqlSM0sIm7LgSN7zHuzjUxhns4ORIT4MnO4kY5zUJA3Tg7a3NZueCMN7jwN3noHFt8Ekx8G04Vg3e5s5v27i5cvGsJ/Zw5uWTnbCLbEUC9atKjBYVVOTk54enpSUlLSoHE0rY+6TN/XY+Tifk1EsoBDGElKQoEdwP+UUj/Ucr6muTlYBB/Hws8xkBgDOasAY+WLYw/oOd1w/Lo0Cib0ADPkyWpVbMvIJy4l64hiPpRvuCP4uzsRGerHteOMdJx9O3lhr0Ol2hdKwRNPwNAouPdMeDAJbuwHxxlGUg4VUmax6jrcJ4EtK+rKGtQNCc/au3cvu3fvxt/fv95jaFondZm+D2Dk5H5AREKAzhgJmLcppYqaXjxNnWzLhgUr4I8Y2BIDhfFAOSDgOgiGXA1nRMEV42BIlyOnVVisbNqXe0Qxr0nNIrfYMJp09nZhbA9/IkONOOYegbrcY7tGATfnw/dfwqv9IAdwG3jEe3/bwXxu+Xwdb1wczoXDg7lweHALCtv2SE1Nxd3dvUYFmpKSwpQpU/jzzz+ZNGlSvedxdnbG0VFvQbRHbM5MppRKxagTrWkpFLB6P3waA8tiYGcMlCWZDY7gFQGj74ZzxhupOLv5Hjm1pNxCorm3HJuSxbq0bArLjHjL0AB3Jg/odCRUKtjXVSvmU4GKCpg4EcaOhz3PwutbjbRGJlVNseFdfQjy1SF09aGuGOqQkBB++OEHJk6c2KB5AgMD6dq1K7m5uQ0aR9P60FULWjMWBb/uhEUxsCIa9sSAZafZ6A4Bo2HokzAtCi4bCVXqKReWVrBu+6EjK+aEPTmUmTHMfTt5cuHwYEMxh/jRwav+5jZNG2Sv6Tw4OwiyJsD5U4y8gFXIyC/hjJf/5eFz+3FJZDdenDWk2cVsL6SmptZq9gY477zzGvzjuLCwkMzMTB2e1Q7Riro1UWyFL5Ngsen4dSAa1AGjzc4fOo+DyJsNx68Lw8H5qJkrt6icNckHiUs1FPPGfblYrAp7O2FgFy+uHN2dyFB/Irr74uvuVIMAmnZPmjJW0GMi4NqvoN/TMOpoc3ZhGXZ2gpeLIzdEhXFGvw41j6WxidTUVMaMGVNj+6uvvsoXX3zBmjVrGjTPvn37OHz4sN6jboecTFEOV6CbUmprnZ01tnG4DD6Nh59iYEMMZC0HTLOVQzCEnQ7jouCSKJjUD6o48GTkl7Bm6+EjpuytB/NRCpzs7Rja1YebT+tBZKgfw7r74qHLPWp+/hny8uCdS8D/a3gzGI5LCV1hsXLmq9FMHtiRZ84fxO1n1L9AhMYgNzeXnJycWj2++/bty5QpUxo8V69evejVqxfZ2dkNHkvTurC1KMd5wEuAExAqIkOB/yilpjWhbO2PXQUwf5Xh+JUcAwWrATOUwqUPDJwFp4+HK6JgWHeoYgrbm110pKJUXEoWuw4bntxuTvYM7+7LlEGdj2SFcnHUMcyaKvwL3PFfGNMdPr4E3IYfo6T3ZBWRX1JB/y5evDhrMP06ebWYqO2NytCs2kzf5557Lueee26jzamrZ7U/bF1qPQlEcrSARoLpBa6pCQXEH4ZPl8PfMbAjBkrXYRQfswOPoRB5I0weD1eNg9CjJkalFLsOFx6jmPflGDHMXi4ORIb6cXFkVyJD/RnQxUuHymiq59lnYf168PkGwqLhE454clflnq8SEISvbhrNxD7a1N2YVIZm1baifvLJJznnnHMYOXJkg+batm0b27dvr3dRD03rxVZFXaGUytWewDbwRiK8/jbsjoGKZPOgM/hFwtgHYWoUXD4GAk5ctSxJ3M8vSenEpWRxuMAo9xjg4czIUD9uMNNx9unoqcs9aurmE+B9V7g9Cu49sXlPVhGrd2UyK6Irr84eipuT3h5pCupaUSulmDdvHkFBQQ1W1F27dqVTp0464Uk7xNb/zo0icilgLyK9gDuAlU0nVhtlN3DnQ8Df0HkCjJgDF0TBrBHgVrtn9codh7lt4Xq6eLswvlfgkVCp0AAdw6w5CQ4cgPHj4eqvYM49Rrqiavh23V4+W72bGeFBBFeJFtA0Lvv27cPR0ZEOHaq3VIgIqampODg0/IeSm5sbHTp0YPfu3Q0eS9O6sPXTcTvwCEblrIXA78AzTSVUm+UjgA0wYxYs/sTm00rKLTz0XRIh/m78dtd4vcesqT+OjhASCpcHQC15SW46rQdzRnbHQW+bNClZWVn4+fnV+mN77dq1DB8+HCenhkVjZGRksGvXLh2e1Q6p8x0VEXvgR6XUI0qpEebjUaWUtq9UxQK8nwXsgzEnlwP5jb+2k5ZZxHMzBmklrWkY/v7g8zu8WLOWLim3sDk9j0BPXWq0qcnOzsbX17fG9sOHDzNmzBh++eWXBs9VVlZGeXl53R01bY46FbVSygIUiYh3M8jTdvkN2J9kPB9su6LenJ7HvOhdzBwezJieAXWfoNHUhFLw4IPQaT2Mq7nbgpWpXPFhnPYObgbqUtQBAQF89tlnDUodWklwcDD9+vVr8Dia1oetpu8SIElE/uRIhQdQSt3RJFK1ReYBnomQj82K2mJVPLQ4CW9XRx45V/+DaRpIRQW88w58PApmhNfYbergzgzs4q19H5qB7OxsOnXqVGufOXPmNMpcpaWlZGVlNcpYmtaFrZsZPwOPAdHA2ioPDcA+YAkQlgiBgdCxY11nAPDpqlQS9uTw+Hn9dbYwTcNxdIQPDsO6GbV2C/BwZlwvbb1pDupaUefl5dGnTx/i4+MbPFdKSgq7d+/GYrE0eCxN68KmFbVS6uOmFqRN8xFgBUg0VtM2rFT25xTz4u9bGd87kGlVqlppNPVm1Sr4JwAO1JxR7Lv1e3nih00kPH6WDvNrBupS1A4ODnTt2hV3d/cGz9WrVy969+5Nenp6g8fStC5szUyWgpHC4xiUUmGNLlFbwwJ8CJxhgVUbYeINdZ6ilOLxHzZiVfDs+QO1CVLTOFx+OZx/Pix+qcYuoQEeXDsuTCvpZsBqtZKbm1uronZzc2Pp0qWNMp+dnR2urq7a96AdYusedUSV5y7ALMCv8cVpg/wJpAH37IK/imzan/514wGWbs7gkXP70dVPx7BqGol3F8G2HtVmH6tkaFcfhnb1aTaRTmVyc3NRStWqqK1WKy+88AI333xzrf1sITU1laSkJFxddTnS9oZNe9RKqcwqj31KqdeA05tWtDbCPCAQ6JhovK5DUecWl/PEj5sY0MWLq8eGNLV0mlOFwkLY0B8e8DWsPNWwNi2L2e+twmrVK67moLI4Rm0KOCMjg8cee6zBlbMA/Pz88PPz0yvqdoitpu9hVV7aYaywPZtEorZEOvAjRvan5ESws4P+/Ws95b+/bSGzoJSPrhyhk01oGo9rr4WCElj7PdQQip9bXI5INXtYmibBFkXdqVMnduzYUWsucFvx9vame/fubNmypcFjaVoXtpq+X67yvAJIBS5qdGnaGvMxVi/XAQ8mQu/eUIvZaU1qFgtjd3PduFAGBeuwdE0jMnMm5LtD75q7nN63I6f3tS0iQdNwbFHUAPb2jZPkKDc3l7S0NKxWa6OMp2k92Gr6nljlcaZS6vpTvi61FfgAmIjx5ZiYWKvZu7TCwtxvEwnyceWes2r5NtVo6kPQTLjmHNheffOB3BJe+HWLNos2I7Yo6v379xMaGso///zTKPNlZWXp97gdYpOiFpE7RcRLDD4QkXUiclZTC9eq+QtIAW4A8vNh165aFfU7y3ay81Ahz8wYqCsVaRqXTz6B/15kVIyvIQ5jxY7DLIxNo7RCr7aaC1sUdceOHXniiScYNmxYjX1sJSQkhGHDhukoknaIrRrjGqXU6yJyNtABuBrD8PtHk0nW2pkH+AMzgHUbjWM1KOodGfm8/c9Opg3pouv9ahqfkhLo4FNtOctKLhgWxJkDOupc8s2ILYra3t6exx9/vFHms1gsFBQU6BV1O8RWb6bKn2jnAvOVUhuoNQiknXMQ+B64EnDGMHtDtYraalU8vHgjrk72PDa1dkczjaZeXHIDxM6DmOqbyyqsLN2cgZeLY/PKdYqTnZ2No6Mjbm41h2Dm5+czadIkdu3a1eD5duzYwZYtenujPWKrol4rIn9gKOrfRcQTMxfXKcnHGC5115uvExPBywu6dTuh66I1e4hLzeKRKf10tSJN45OWBs8+Zlh3Aqvv8uvGdK7/JJ7cYl1ZqTmpzEpWmym6rKyMnTt3UlhYWGMfWwkLC6NHjx4NHkfT+rBVUV8LzAVGKKWKAEcM83e9EBEfEflGRLaIyGYRGS0ifiLyp4hsN//6Vun/kIjsEJGtpvm95bAC7wPjgb7msaSkalOHZuSV8PyvmxkV5ses4bUUB9Zo6ssff8AH78Dv5Uc/j8cxsW8HvrpxNN6uekXdnNSVPhTA39+flJQUBg0a1OD5HB0dCQjQOdzbI7Yq6tHAVqVUjohcBjwK5DZg3teB35RSfYEhwGaMHwJ/KaV6YbhqzQUQkf7AxcAAYDLwtlkju2VYBuzAcCIDo7RgYiJU84/21E/JlFZYef6CwdrBQ9M0XHcd3LETVteshCssishQnUiwubFFUSul+OyzzygtLW3wfPv37ychIUGHZ7VDbFXU72DUpB4CPICRNPOT+kwoIl4Y69EPAZRSZUqpHGA6hlEZ8+/55vPpwCKlVKlSKgVDTUbWZ+5GYR7gC1xovt6zB3JzT9ifXpp8kJ+T0rnzjF6EBjQ84b5GcwJKmStqD9hQfZeY7YcY9vSfZBY0XBFoTg5bFHVKSgqXX345K1asaPB8zs7OuLi4NHgcTevDVkVdoQwPhenA60qp16l/ZrIw4BAwX0TWm+Fe7kBHpVQ6gPm30j06CNhT5fy95rETEJEbRCReROIPHTpUT/Fq4RDwHXAFRsZzqNaRrKC0gsd+2Eifjp5cH6XrlmiaiH/+gcmT4eddcGv1XcICPbj/7D74uOkyqs2NLYo6LCyMf/75h9NOO63B8/n7+9Ovn65r3x6xVVHni8hDwOXAz6bpub4bXg7AMOAdpVQ4UIhp5q6B6mzG1bo1KqXmKaUilFIRgYE1eNY0hE+AMo46kcFRRT1w4JFDL/2+lQN5JTx3wSCcHHSaUE0TMXIkvPUT0LPGGIwgH1dundgTe10tq9mxRVEDjBs3rlGykxUXF7Nv3z7t9d0OsVWLzAZKMeKpD2CsaF+s55x7gb1KqVjz9TcYivugiHQGMP9mVOnftcr5wcD+es5dfxSG2Xssxm55JYmJEBpqeH0DCXty+HhVKpeP6s7w7g2rhqPR1IqdHSyeCk9Xr4S3Hsjn9JeXkVVY1syCaaxWKzk5OXUq6r179+Lm5sb69esbPOf+/fvZs2ePVtTtEFtTiB4AvsWIGgY4jGEEPmnMsfaISB/z0BlAMkZ5iyvNY1cCP5jPfwQuFhFnEQkFegFx9Zm7QUQD2zjqRFZJldSh5RYrc79NpKOnC/ef3ef4ETSaxuPgQejYEWYvrfEns1Upgn3dtFWnBcjLy6uzxCUYFa8uvfRSgoMbHhUSFhbG8OHDteNqO8TW6lnXY6goP6AHxor6XQwlWx9uBz4XESdgF0aolx3wlYhcC+zGqHmNUmqTiHyFocwrgFuVUjUU8mtC5gHewMwqx0pKYOtWoyAC8H7MLrYcyGfe5cPx1MklNE2JtzdcMgfOHAzdq+/Sr7MXn1zTcn6XpzK2FuRwc3NjwYIFjTavUkqvqNshtv7UvhXD6JsHoJTazlFnr5NGKZVg7iUPVkqdr5TKNmtdn6GU6mX+zarS/1mlVA+lVB+l1K/1nbfeZGLYEy4HqiYZSk4GqxUGDyb1cCGvL93O5AGdOGtAp2YXUXOK4eICnd+B86r/N8wqLOPBbxLJL9FJTloCWxV1UVER1157LVlZWbX2s4Vdu3axbt06rajbIbYq6lKl1JGNLhFx4FQqa/spxg59dWZvQA0axCPfJ+Fkb8dT0wccf7ZG07goBZdeCs6/wR3Vd0k5XMAvG9MpKK1oXtk0gO2K+tChQ3zyySfs2bOn1n62EBQURLdu3bSibofYWpTjXxF5GHAVkTOBW4Cfmk6sVkSlE9ko4PicJomJ4OrKt7kurNiRyTPnD6Sjl45j1DQDW7fCZVYjqW81DO/ux7rHzsTRXu9PtwS2Kuru3btTUFCAs3PD0wu7uLjQrVs39u7d2+CxNK0LW/+LH8SIIk4CbgR+wchO1v5ZgZE37fjVNEBiIhX9B/DMb1sZ3t2XSyNPzPWt0TQ6IvDOWnjxXKjGW6PCYuXr+D1YrHpl1VLYqqiVUsTGxjbKKjgzM5OkpCSdmawdUqeiFhE7IEkp9b5SapZSaqb5/NT4Fngf8AIuqqYtMZF4zyAKSyt44YJB2OlYVU1z8N13EL3ccK2s5jt5/Z4c7v8mkT1ZRc0umsbAVkWdnJzMaaedRmxsbK39bKG8vJyyMh2K1x6pU1ErpazABhE59ZaL2cBXwBzg+CygBw/CoUP8bt+Bmyf0pFfH+iZq02hOkv/8B/YvNspaVhNcMLSrDz/dNk5/JluQ7OxsHBwccHevPX1wv379+Oijj4iIiGjwnJ06dWqUcTStD1v3qDsDm0QkDiOTGABKqWlNIlVr4TOghGrN3iVr1+MC5Pbqx4MTdGk5TTPy/e/wgSfkU20i333ZxQzo4tXsYmmOYkuJSwA7OzuuvrrehQiPoby8nIyMjLo7atoctu5RPwVMBf4DvFzl0X6pdCIbAQw9sXnF4r8BuPTaKbg4tlwxL80pRkYGrCuA51xh34nNm9PzmPDSMjbub0hxO01DsTV96L59+wgLCyMlJaXBc6alpbF169YGj6NpfdS6ohYRF+AmoCeGI9mHSqlTI94jFtiIsUd9HBv35ZIbu45cv0AiIno3t2SaU5k77zSU9Z6/oMuJzV393Hh8an/6dtIr6pbEVkXt5ORE165dG6XqVVhYGEOHDiUhIQGllM5Q1o6oy/T9MVCOsRt2DtAfuLOphWoVzAM8MCphV6HCYuWhxUm8dDgVt+HDWkIyzanMXXfBDjsjS141uDvZc8240GYVSXMi2dnZBAQE1NkvMDCQf//9t1HmtLOzw9NT+yW0R+oyffdXSl2mlHoPI3lmVDPI1PLkAouASzGUdRUWrExl8+5MembuwTF8SAsIpzllUQr6jIArR1SbaT89t5iRz/3FJm32bnFsXVGXlpbyn//8h5KSkgbPuWfPHpYvXw6gk560M+pS1EfyD54yJm+Az4FiTnAi25NVxMt/bONin2Lsy8uOqUGt0TQ5774LUUPgI2DGic1O9nYM7+6rk+60AmxV1Nu2beOJJ55g27ZtDZ7T19eXTp10+uL2SF2KeoiI5JmPfGBw5XMRyWsOAZudSieyYcDwKoeV4rEfNiIC93QsNQ5qRa1pTkJCYPoMuIITwwUBfw9n3rlsOAEeDc9ypak/tpa4BBg0aBA7d+5kcCN8l3h4eNC3b1+g/a2oRSRERLaIyAcislFEPheRSSKyQkS2i0ikiLiLyEciskZE1ovI9CrnxojIOvMxxjw+QUSWicg35tifSyvd2K91j1opdeq5M8cDGzBqg1Xhp8R0lm09xONT++P387/g6Ah9dClLTTNyzjnwzDnwPPDQsU2FpRU8+G0i953Vh5CA2mN3NU1Lfn4+VqvVJkUNYG/fOF+zBQUFbN68uVHGqo2nnnrqNaqNhWkQCU888cRddfTpiVFV8QZgDcbm5DhgGvAwRoXFv5VS14iIDxAnIkuBDOBMpVSJiPQCvgAqA87DgQHAfow8lGOB5Y14XY2CTgR8PPPMvwM5kvUpr6Sc//y0iY5ezkwb2sXI8R0QAFu2GNWzNJqmZssWQ1GPscDoE5vzSspZvSuLkormrwCrORZbs5IBbNiwgZCQEDZt2tTgeXNzc4/EUbe3FbVJilIqyUzCtQn4y8yQmQSEAGcBc0UkAVgGuADdMNICvS8iScDXGE7RlcQppfaaYyaY47Q6bE14cupQmYFvHOALjAXLcEWvPZ6s8csi4pmlvJ9RwZnp6TBkCFZvH+zGjYWoKBg/HoYPByenFrwATbskJQU2boSFebDMEx4th6ddwTTUdfZ2Zc0jZ5C0L5eVOw4zpmfdHseapqG83HDtcbLhe6B37948/vjj9OzZs8HzBgUFMX78eJYtW9bgsWrDhpVvU1Fa5bm1ymsrhi6zABcqpY4JJheRJ4GDwBCMxWlVz72qY1popTqxVQrVoizASO8SbT5iwHeJE18wCourYn+/YtYED+G+GWuxK45haPpGxsYl0f3nnwGwuroiI0cilYp71Cjw8Kh5Po3GFs45BypLIc57EP78EO5KAx93I+e3C4gI36zdyy9J6cQ9PEnnnm8h/Pz8gKMr69pwdXXlqaeeapR5K/fGgVM1hvp34HYRuV0ppUQkXCm1HiOYca9SyioiVwJtbktXK+rjEQzjRwiG0w7AAWA52EcLXaPd6PqTGxeoYKyO0zjUu5R1kdm867KJkpJ/GXgwiTHbkunz77PYqadR9vYwbBgyfryx6h43Dvz9W+jiNO2CebfBX/3A3w1mR8OvsZBxH7gKj0/tz71n9cHOTnh72Q7sRbjxNJ3itjnx8fFBRMjMzKyzb1ZWFjNnzuTTTz8lKCioQfPu2rWLhIQEwIipPgV5GngNSDSdwlIxMmq+DXwrIrOAf6iSBrutoBW1LXTCiCKfab7OAVaAXYzQMdqFc37vzDkVnVFyBtk9y0kalMOXY3ZyuHQZfQ5vYPT+ZIa8/iaOLxtZV1X//kcVd1QUdO3aMtelaZt07QpXXWU8H7YNlr4O5TfCGg8ctpXjfZ3h9W0vwqb97TM4ozVjb2+Pj48PWVlZdfYtLS1l165dFBU1vNJZSEgIffv2ZcuWLQ0eq7WhlErF8ByqfH1VDW03VnPudqCqW/1D5vFlGHvZlf1uazSBGxlpp04HREREqPj4+OaZrAgj5ahpKmclRhw2kN+tnM098/jHbw9p5f8Skr2OUfuTidi3Gbdi44ed6t79WMXdp49Rc1ijORkmPQV/vwb7d0GnYx2ZLnh7Bd383Hjt4vCWke0Uo2fPnowaNYrPPvusWef9+++/iYmJ4YknnmjIMPrLp5WhV9SNgRsw0XyAkSZmHRANntGORC73JzLHHxhKcUcL2/rk8+GQA2yqWE6ngjWM3JfM6O+X4PPppwBYAwOxq1TaUVEwZAg46LdKUwcfXwdLOkFHH7gyFpZGw677wFm4/fReFJUZHuGJe3Po4uOq462bED8/P5tM30opFi5cyIUXXtjgfN/79+8nLi7uVN2fbtfob/+mwBEYaT7ux/BJ3AjEgGu0PUOifRgS7QP0pcznGlL6FbK47yHiJRavotWM2LeRsf+uotPixQBYPT2xGzPmqOKOjIRGSOKvaWcEBcGNpuWvx1b45v8g9yo4GMDEdF84yxGlFDd/to7h3X154xK9um4q/P39bTJ979u3j8suuwx/f38mT57coDmdnZ1xdXU94nWuaT9o03dLoICdHDWVRwO7jKYKdyt7+xazOjiT5Q4JOJasZNj+JMalbyb0gFEKTzk5IZGRRxX3mDHgXUOVBo3mtBch+hnYsx2CO7A/pxg7ETp4OvPq0m2cPaATA4P056cxmTNnDrGxsezYsaPOvtHR0YwZMwaHRrCa/fnnn8TFxfHII480ZBi9JG9l6BV1SyAYOXZ6AteYx/YBMeAQbUdIjDshP7hzMd2wOp3Hgb4lrB2ezesumykvj2HwgSTGpW2m76oXsX/+eZSdHQwZcjQkbNw46NixxS5P08pYeDl84gKd/eHOBLr8/Sesu5cyq+LfbYcQ0Iq6kfHz87NpRQ0wZsyYRslOVlhYyJ7KED5Nu0Ir6tZCEEZJzcqympkYnuXRQpdoV7r85sp5li4ou9PJ7FXGhv45fDZyF1kVMfQ7lMCY9M0Mffc9nN54AwDVu/dRxR0VZeSJ1ntXpyZBneCh243n7lth27OQfglOBPHjxNEwwJ6cojIuem8V95/dlzP76x95DcXf35+cnBwsFkutSjgvL4+goCC+/fZbzjrrrAbNuW/fPvbu3dsoK3NN60K/o60Vf4wMttPM1wXAKpAYISDamTP+7cgZJR2B0eSGlLOxRy7/Dd9DmnU5YdnrGbM/mRGLvsb9ww8BUEFBxyru/v3h1Iy1PLV5brbxABj/DsQ8ACmb8eoWRFSvQFwc7VBKsTe7mGBfV+2YVE/8/PxQSpGTk4N/LXkT3N3dmTNnDj16NDzWvVevXgwePJitW7fW3VlTb0TkKuAPpdR+8/UHwCtKqeSmmlMr6raCB3Cm+QAj8V08EAPe0Y6MXR7A2PwAIJzCLhVs7ZnPO4MOsFmtonN+PKP3JTP6t7/wW7QIAOXnh4wde1RxDxtmFBrRnDrMvxDeLYYuHbB7dguPLf0Z/rqHHYcKmPRKNK/NHsr54Q1LwnGqUqmcs7KyalXU9vb2vPvuuzW2nwxKKZRS+sdV03MVhnvwfgCl1HVNPaFW1G0VZ4w6L2OBuRhZahOBaHCPcWBYtC/Don2BfpT6XcnOPgV81fsQa+zX4lsUS+T+TYyLTaDLTz8BoNzckNGjjzqojRoFbm4tdXWa5qBHB3jxHuP54W2w4knYMYMenXvw8bQIBvTyIbe4nB8T9jEroisujm0u82KLUZlGtK59aqUUd999NzfccAP9+/evtW9dbNu2jaSkpAaHebVWROQK4D4Md9xE4FGM6uyBwCHgaqXUbhFZAORhVMjqBDyglPpGRL4EPlZK/WKOtwD4CfgeeAGYgPHN+pZS6j2zzwPA5RixO79iLI8igM9FpBijRM6vplwjgFCl1APmuVcBw5VSt4vIZcAdgBNG1o1blFI2V9DRirq9YI9RsC0cuBPjo7wViAHnaHv6x3jTf5U3N9GTCo9ZpPYr4q/Rh4l2TMKldBXD9m9k/JYthP39FKIUysEBiYg4qrjHjgXzy0fTDnl9OrxuJOCR0z7htJhbYGsiG1z8eOLHTQT7uTGxT4cWFrLtUKmo64qlLisr47PPPmPs2LENVtRdu3alS5cuNuUYbwgi8hpNUOZSKXVXLXMOAB4BxiqlDouIH/Ax8IlS6mMRuQZ4AzjfPKUzRmmlvsCPwDfAImA28IuIOAFnADcD1wK5SqkRIuIMrBCRP8xzzwdGKqWKRMRPKZUlIrcB9yml4k3ZKsX8BlgFPGC+ng08KyL9zOdjlVLlIvI2MAf4xNaboxV1e0UwPmZ9gevNY2kYnuUxdvSM9qDnDx5cTggW56ns71tMXHgWL4/eRkXFCoYcSGL8gc30e/U1HF580Th/0KCjijsqyojb1bQ//jcFXnsMunRkyCcH2bZ5E/ZhZ7E2LYt3/93FSzOH4O2mt0lqo6rpuzacnZ3Zv38/jo2w7eTu7k6HDh3Izc1t8FitkNOBb5RShwFMhTkauMBs/xT4vyr9vzdLVyaLSKV35K/AG6YyngxEK6WKReQsYLCIVCaJ9gZ6AZOA+Uqposo5axNQKXVIRHaJyChgO9AHo8b1rcBwYI2p1F0xamTbjFbUpxLdzcdl5usMjGIjMWaxkV/cuNAajHKYSEavUtb1zeajEankWFfSP2MD4w9sZshHC3B5+20AVFiY4aBW6aTWs6f2LG8PDPKHDx80nsdtw+HHx2HjZLw6h5JXWE52URlerg4oha7QVQO2rqgBNm7cSM+ePfHy8mrQnOnp6ezYsaPJ96hrW/k2IYJhJ6yNqu1Vy1cKgFKqRESWAWdjrHC/qNJ+u1Lq92MmFJlsw5zH8yVwEbAF+M6s4iUYJveHTnKsI2i331OZDhi/R18F1gLZwK8g9wsd/V045+/OvPr1aOZ/ey/X7fuAku4/8/TUGK6Y+Rb/Of16/nHuTP7iH+C666B3b1TnzjBrFrzxBiQkgMXmLRhNa2X+ZLCWQPhges35gy9vPZOQjD18ujqNIU/9waH80rrHOAXx9vZGRGzaoz7ttNP46KOPGjxnWVkZ5eXl7dWZ7C/gIhHxBzBN3ys5GtA6B1huwziLgKuBKIyymJh/bxYRR3Ps3iLiDvwBXCMiblXmBMgHPGsYfzGGufwSDKVdKftMEelQOY6IdLdB1iO0yIpaRO4GrsP4tZKEcePcMC4sBKM82UVKqWyz/0MY+wgW4I7jf/loGgkvDINQZSbDYiAOoyZ3tBOnrezAaYUdgGHkB5ezuUceb/Q7yGaHNXTNiWfMgWRG/7OCgG++AUB5eRme5ZUr7ogIcNb5pdss/zkdXnkYunRi2rpMRqxZiY/DaezPKeZgXgnh3XzrHuMUwd7eHl9f3zoVtYiwcOFCIiIiGjxn9+7d6d27N7t3727wWK0NpdQmEXkW+FdELMB6DOesj0TkfkxnMhuG+gNjb/hHpVSZeewDDL2zzlz9HgLOV0r9JiJDgXgRKQN+AR4GFgDvVnEmqypntogkA/2VUnHmsWQReRT4Q0TsMKpB3IqxGWkTzZ5CVESCMH759Df3B77CuAH9gSyl1AsiMhfwVUo9KCL9MUwUkUAXYCnQuy6PuVadQrStUoHx71GZ9jQGML+HigOMYiP/dsggznEDAYWxjEpPZlz6FoLTzdSnLi5HU5+OHw+jR4NnTT9MNa2aC/+BxWfDimheKfTj3WU7Wf/4mbg76920Snr16sWIESNYuHBhs8xXXFzMF198QUFBAXfccUdDhmqXS/K2TEv9VzkAriJSjrGS3o9RI3SC2f4xRp3QB4HpwCKlVCmQIiI7MJT2qmaWWeOAEYAwArgHI2BhMxANrjFmsZEVPkBvyrwuJKVvAb9HHuZf5824la5mRPomxqdtocfyF7B79lmUvT0ydOhRxT1uHAQGttjlaU6CbyZCRSk4CndO/5W7f56NzIzhF4cOHMwr4eqxoS0tYYtjaxrRM844gyuvvJIrrriiQfNt27aNPXv24OurLRvtjWZX1EqpfSLyErAbw7j6h1LqDxHpqJRKN/ukV9rzMZJrrq4yxF7z2AmIyA3ADQDdunVrqkvQVGIHDDAfN2NsZKQAMeAUbUefGC/6xHlxLWFYXM9lT58iVg/K5PmRO1CWVYSnb+S0jK30f/sdHF97zRizb9+jSViioqD7SW3laJoLARyNhZf9dWOh4AEI6ozly314L/4My5KnsXc/tePw/f39OXz4cJ39unXrhncjFNXp378/GzZsIC8vr8FjaVoXza6oRcQXY5UcCuQAX5vB4DWeUs2xau31Sql5wDwwTN8Nk1Rz0ggQZj6uNI+lY3iWRwsh0e6E/OzOxaobVseJHOhdwtqe2cwLTyVXYhl4MJHTMrYw9LMvcJ03zzi/a9djFXe/ftqzvLVxnhec96jx9NedsOx1WHUuD+YEUlJh4fWLT81ymn5+fjal85w/f36jzOfg4ICLiwsFBQWNMp6m9dASpu9JQIpS6hCAiCwGxgAHRaSzuZruzNE4s71A1yrnB2OmbtO0AToDs8wHGJ7lK8Auxiw2stSV8yq6oGQ0mT3KSAjN4bEBe9nnsIZehxMYf3ALI37+He/PPwdA+fsfDQmLioLwcNBFCFoP346GQ6XQzZ7bL/iNzj/NhtA/Se87CDcnB7xdT534a1trUn/22Wf079+fYcOGNWi+zZs3s3Xr1lpTlmraJi3xDbcbGGW6vBdjZIeJBwox1mEvmH9/MPv/CCwUkVcwnMl6Yfgia9oivsBU8wHGux4LEi0ExDgzaXlHJhV3BIaT262cjWG5vNLrAFuc19MtZy3jDmxh9Mp4Onz/PQDK3R0ZM+ao4h45ElxdW+TSNBipHLoZqUaDp0ZB2j3QpQvP/N8KJm78g5kLnwB395aVsZnw8/MjJyeHioqKWitavfDCC5x33nkNVtSBgYFHwsI07YuW2KOOFZFvgHUc9SOeh1F24isRuRZDmc8y+28yPcOTzf63nkyOVE0rxx0j59Dp5usyjE9GdJViI7kBwEAKO17C1p75fDz2IPHOmwgoWsPY9M2M25xM16VLjdSnjo7IiBHHpj718WmhizvFucYdrnkCgBdX7Mf131fgrygWB4fTu6Nnu6+BXZn0JCcnh4CAgBr7/fXXX3h4eDR4vsDAQAIDA/UedTuk2cOzmgsdntVOsGDUqakMCYsGDhpNpT4WdvQpIKbTIf5x3YZHWRwj05OZcHAzPXZvxb6iHCWCDB58bOrTzp1b7HJOWQ4CqeWoEfZ8ddb3nB9zHc5LfzDej3bKwoULmTNnDlu3bqV379419jt8+DAuLi4NVtYHDx7kiy++wNXVlRtvvLEhQ+kleStDb+5pWjf2wBDzcRuGG+EOIBqcY+wZEO3NgFiz2IjbZFL7FrGq/2GeHrEL1BqGpycxIWMr/T/4EKf//c8Ys2fPoyFhUVEQFqYd1JqajkBHRwSYNXEyknorZR0688DrMdy4bzn9Hrut3cXU25pGdNq0aQwfPpw333yzQfMdPnyYgoIC3E+RrYVTCa2oNW0LwfBS6IWRqw4Md8MYcIi2o2eMBz2XGMVGrE4T2de7mLjQLN4asoc8uzgGH0xk4qFtDPnmO9wrvW07dz5WcQ8cCHY6u25TYfeIGzzyNI5KceusTHomPQ7hoZTMnNWuSmnaWury7rvvpkuXLg2er3///qxbt46ysrK6O7dRzKxkSRi6azNwZWXRjPaMVtSatk8wRmbdS8zXhzE8y6PNYiN/unGhJRhlN4qMnqWs657NQ333sc95PX0PJTDx0FaG/xWN71dfAaB8fIzUp5WKe/hwcHJqoYtrv4gIvb4aCZvzyT5TeOPcxTy06hacflgEZ5zR0uI1GFsraM2aNavWdlspKyujpKQEe/v282OnGoqVUkMBRORz4CbglRaVqBnQilrT/gjAiNSfbr7OB1aBxAgdo104J7oz55R2BiLIDikjMSSXt0cdYo/DFnxLVjEyM5lR6zfR6eefAVCursjIkTBmDAwdajx69NCr7sagL9DXBV9gjt/p2LndQEWHTixZtp2z437F9brL2mwd9MoV9aFDh2rt99RTT+Hn58ftt9/eoPkSExNJT08nODi4QeO0IWKAwS0tRHOgFbWm/eMJnGU+wCiAtwaIBt8YJ05bGchpeYEY6eYvICugjGT/XH4ZvgurxOBbHEv/bYn0ivkvDmZFMKubOzJkCDJ0yFHlPXAguJ3a2bgaQs8v/YFn2ZtdhNcFG3DZ8SB4uZBzxVV4ZaRjFxzUpmLmvb29CQsL49VXX2XOnDl07Nix2n4HDhygtLThVcj69+9PXFxck6+oQ+b+/BowtJGHTUh9YcpdtnYWEQfgHOC3RpajVdJ2PvUaTWPhDIwzH2A4qO3G2PlKAr9EJ8YlBTI2IQCxjASgwsFKWodsdrvFU+gQj5MlkQ5pW+i1/hPcS94xhrGzw9q7N/bh4YbiHmIq8Rq+oDXVE+zrhs/iCKybcyk93443pq7msT9Pg0ceh2eeaGnxbMbOzo5vvvmGsWPHMnPmTP766y+cqtlCeeeddxplPnd3d5ydndu76dtVRBLM5zHAh3WdICJXAYeVUkuqHLNTSlmbRMImQIdnaTQ1UYpR/j2RI0pcJSpk/1EP8Xy3MlJ9N5HhEociAc/STXTP30an3INH+lR06Ih9eDgSPvSoAu/VC9r3F2qjkfBZDr0+WIfb+z144qsc7v/fTXjeMxPuv7elRbOJL7/8kosvvpgbb7yRd999t9o+27dv5+6772bBggW1xlzXRnJyMt999x2hoaFceumlDRG51YZAiEiBUsrjuGPdgXsx5N6JUT6yK7AJGImR1NgNWIixMTYJI8nWd8edlwNMxEjElQ44AgOBizByLB7pq5R6remu8kT0ilqjqQlnjoaGmQgCmRxR3J5JTgxKDMe6cSh2hUe/39J997PXfTVFDmtxtCTRce0Wui1diqOlAgCLiytq0CAchoUfNZ0PGnTKZO06GYZe5gOXGRlxIiyOuGZPA7su/BS7j8EffEn3qT1h+rSWFbIWZs+eTUJCAi+88AJDhw7lpptuOqGPvb09KSkpZGVl1VtRu7q6Ym9v395X1NVxC4ZyLQYGKaVeE5H/AVHAxcDlmCtqc3X9q1LqcxH5b9XzMFbovyulForIX0qpM0TkYYyyQxcf17dZ0YpaozlZ/DEKsk44esjOKkblsCQgETondaFj0gxk+wzEaijwMocSdgbGk+kUh9UuAc+0ZLolLcTzvfcAUCJUhPXEYViV1ffQodCpk47zNpn2eBd4/CEAUudu57wP4mDzSg5OPAvL9h108XI2rBWtjGeeeYbExERuv/12+vfvz/jx449pDwsLY9OmTQ2aIzQ0FHd391rTlbZT7IBPlVKJAGLkUPXBKMTrYP6tSm4N510FVKZ1q/QALMP4yX5M3+ZGm741mqakGCP5rWk+V0kK6wawP1SpeBU5bjvZ77GKYod1OFkSCSzaQqf8o3Vnyv0DsBs6FPuqq+/evduUY1WToIC1wCB46q1NPHL/I9gHJCIHdrbKHza5ubmMHDmSrKws4uPjqy3F+/zzzyMizJ0796THP3z4MB9++CF9+vTh/PPPb4iore/m1YJp+n4Ow1ydj6GItwMHgNnA58AjGKZvH46uro8/L61K2yKl1MUich+w3OxzpK9S6qlmvEStqDWaFuEgR8znJII1UcEmsCs1viOt5HDAK5ZslzUoMfa+u+TvxNFSDoDF2QVLv/44Dg9HKp3XBg9ud9m9bCUrt4yKZw7TYVIRbynF1Ctfp6vjH9jtSAQXl5YW7whbt24lMjKSHj16sHz5ctyOixJ48cUXSUxM5NNPPz3psVevXs0ff/xBeHg45513XkPEbFOK+lRAK2qNprVgwUiPaipvlaSwbFA4pFTGa5dT4pDMAY/VFDuuxcmSRIfCzXiW5h4ZojQkDIfw41bfQUGtcoXZVGxNzcP1/C0E94ph3f9dRf6DvzK+4h/sF77VKpT2L7/8wtSpU5k9ezYLFy5stGpXFRUVvPzyywwaNIhzzz23IUOdOh+WNsIpbjvTaFoR9kAf8zHTcFxzQKAAw3810RGXpCF03TAYlahwyLHDsP/uI9t1DdkucajcBLx+X0Pgd4uPDFvm44caMgSnYUOPrr779gXH9lkbuk+IFyREApGUbDvMwD+8wD6ZvEPZrN+QSdSBDdhdc0mLJaw599xzee6553jooYcYOnQoDz744DHtOTk5zJw5kzfeeIP+/fvbPK5SCovFcio6k7V7tKLWaFo7HhhBJkZIN/YIKDF2yxIFkoLx2RCE+4bpOGwV7MoFyKdCEsh0W0Mxa3Fel4T/8rdxshiJNSyOTpT17Yfj8GE4VA0b825fpSfH9Q6AA9PAYRrrdx2i6xUbscu+DPqHUjw8EtfCvBbJfPbggw+SkJDAQw89dMIK2MHBgYKCgjozmh1PXFwc5eXlWlG3Q7TpW6NpT5QD2zhiPq/YYMWSoHDeb49Rzn0bRQ7x5LrEoyQBr9JNeJQdzUVdHNwNGToU54hhSKXpvFu3dmM6z4wrxTe7kM0D7dh25p/M2HwR/PAjTGvQnm69KCoqYty4cezatYvY2Fj69OnToPEyMjJ45513GD9+PBMnTmzIUO3jzW5H6BW1RtOecMSI+jQjPx2wM/7Jc4GNDpDYH5fEftitm4N9suBYJhjOsQnkOsdRnL0O56UbcVryE4LxI77M05vygYNxGTHs6N53v35tslCJf6Qz4Ew/qyL7+tMo+etTKsKH894tK7j1u7txWfAMnH1WneM0Bm5ubnz33XdEREQwffp0YmNj8a5i0fj444+JiYnhgw8+sGm8yiIgekXd/tCKWqM5FfAGxhoPOwQX7I3t7T1AYmdI6ozrurNx3GDFZZc9dhQBSVTIOgrK41EJCTjFvoe9tQQAi4MjxT174zAsHJcRw4+azn19W+wSTwY7O2Hc3R3g7ssorbDQvawQ+/xBUCJ889tuTl+0EL/7phr525uQ7t2788033zBp0iTmzJnDDz/8cETRenh4kJqaavNYycnJAKdiHHW7R5u+NRrNsVSmTk0Ca4KieL0FuyTB9RAYbukJFNuvpcRxHc4VibhVHN1LLeoUhGXIENxGDD+6+g4JaTOm8wqLlTdmbeKe72bBNdPY//JTeKak4tm9Y5PuZb/zzjvccsstPPzwwzz77LP1GiMuLo5ff/2VyZMnM3LkyIaI0zberFMIrag1Go1tZAIbgUQoWWehfL0Vl232OBYfAjYACRQ4xoMk4Fa2EzszIVSZuwfF/QbiHDHs6Op7wABwdm65a6kFVaGQnQK94aEbNvDcBzOQ8UPg3++abk6luOmmm5g3bx6LFi1i9uzZR9rOPPNMHnzwQSZNmlTrGAUFBbz88stMmTKFiIiIhojT6hS1iFgwPC8cgM3AlUqpogaO+QtwqVIqp4b2u4B59Z1HRCYAZUqplfUU8QjaRqLRaGzDHzjNeLhgb5jPrUBKR0g6C8uGSah4i6HM08pAbQI2YFe0HocNa7FftwDefRsAi70DhWG9YMgQPEYOxy483DCd1zPPdWMiDmKEyAE3PNaDLK+f8L/Ylaff28SdDzyMV6QL/Pll484pwptvvsmmTZu4+uqr6dOnD0OHDgWMqli2eIBnZmYC7XaPulgpNRRARD4HbgJeqeskEXFQSlVU16aUqivY/C7gM6C+PwgmYARXNlhR6xW1RqNpfCpTpyZBUXwFpeutOCeDW04akABsoNR+HZCAsyX9yGmFgZ0pGzQI16qr77CwFot5rsq/qzPoc92/dDormz+umUa3e76nz9Bs5P8earQ5Dh48SEREBHZ2dsTHxxMYGGjzuUuWLGHt2rXMmDGDwYMHN0SM1riiPlI1S0RuAgYD9wNvYhTJcACeVEr9YObsngK4AO7AHOBLwMvsd7NSKkZEUoEIjE/rV0AwRjaDp4GOwEvAVoy0ohNF5CzgKYzc3zuBq5VSBeY4HwPnYbhzzgJKgNUYaYwOAbcrpWLqff1aUWs0mmbDTJ1akWAlP74ClQieO7JxLE8CErCygXL7dThatmKHBYAyF3cK+vbHYVg4HpHDsRsWbjh5ubq22GV8/kcqM6Z9hWuHJWQk/kbm7zvo75UN55zW4LHj4+OJiopi5MiR/Pnnnzg6OvLnn3/y4osvsmTJkmprWgOkpaWxYMECZs6cyYABAxoiQu2KWngNGNqQCaohAcVdNU5pKmoRcQC+BX7DKGWZrJT6TER8gDggHENRPgMMVkplici9gItS6lkRsQfclFL5VRT1acBkpdT15lzeSqncynal1GERCQAWA+copQpF5EHAWSn1H7Pfy0qpN0XkFmCYUuo6EXkSKFBKvdTQm6NN3xqNpvnoaDwcJtnhi6lwLB1hR0dU4hkUrbVQvNaC06ZSvNO3Awk4lSTgvWE9bFiI3UfzALCKHflde1ARPgSPkRE4Dzcd1zp0aJbLmHNWCBQ9AHYP8P0/O7ns8m9QFc8hWYexeHphb6mod/haREQEH3zwAZdddhl33303//vf//Dx8aGgoICKiooaFXWlt3c7NX27ikiC+TwG+BDDpDzNLJwBxgq6stLJn0qpygQBa4CPRMQR+F4pVTlOJUnAS2bZyyU1rHxHAf2BFWbKVydgVZX2ylSAa4ELTv7yakcrao1G07KYqVOlj+AxywEPHABnKBgOm4ZTvt5Kblw51g0WPLfuxbVwI3YqAe/dCVj2rMb+h2+ODFXk3YHCgQNxGjEMr5HDjZSpPXtCUygv0xp/7fhQMv68D3f76/khtZA+Z31An6zHkK0boUePeg09Z84cEhISeOmllxgyZAjXX389K1fWvtW5du1aoBkUdS0r3ybkyB51JWY5ywuVUluPOz4SKKx8rZSKFpHxGObwT0XkRaXUJ1Xat4nIcOBc4HkR+UMp9Z/j5hcM5X9JDfKVmn8tNIFe1Ypao9G0TszUqY4j7Qi4yfQQV70hvTdqwwzy1lRQtLYCx8QcfHYn42BNxC03AdcVCbDiNQTDh6jC0ZW8sH5YIobgPToCp+HDYNAgcHdvFDEd7O3ocpoX4MXEknISrj+PsBRn9jp6sWH6EqZsegOnP9+D0NCTGveFF14gMTGRW2+9lf79+zN27FjmzJnDxIkTue66607oX7mf3U5X1NXxO3C7iNyulFIiEq6UWn98J7Oc5T6l1Psi4g4MAz6p0t4FyDJN6AXAVWZTPuAJHMbYb35LRHoqpXaIiBsQrJTaVot8+Rj74g1GK2qNRtN2EKALSBfB+xxHvHEEXKG8M2w7g9L1FrJXl2FZX4pX8jY8c5JxKE/Ab2sCauti5PP5ACiE/A5hFA8eiMuY4XiNijBW3506NUg8LxdHxj9rVFZxyy3BqcwfyXakPCePXz/bzpnRP+L6wtU2xWTb29uzaNEiIiMjufDCC4mPjycoKIjs7Oxq+3fs2PHIeacITwOvAYnm6joVmFpNvwnA/SJSjuGFfcVx7YOAF0XEipGE92bz+DzgVxFJN53JrgK+EJHKuMJHMRL21sRPwDciMh3tTFY92plMo9GQa5QLzV5dTsGacpzW78E3bRPOZRsxvM8TML7fDUrcAsjvOQA1egg+YyNxihgGvXs32HR+ILeEpZNWMid+KvK/N9l/yaUEpO3HaWC3OquYJScnM2rUKPr06UN0dDSuNTjRLVu2jH///ZfrrruOoKCghojb6ry+T3W0otZoNKcWZurU4rUWslaVUrEmG89Nyfgc3oSd2gBsQLERoRyACnsX8jv3pXToIFxPG47X6AhkyBDw8DipacvLrDgetlIcILw5OZb7/4lC7noCXn28znN//PFHpk+fzmWXXcYDDzzAtddey48//kinKhaAr7/+muTkZG688cZjjtcDrahbGdr0rdFoTi0E6Aau3ewJmuEGuAFBUHom1mRF5upS8lcV4bR2G76pybgXbcR3bwLsXQJLPgUM03mhdwiFvQZA1FD8JozAcfhw6NKlxnSpjk520MUOV+DM1/qRt+BPvG7tyVMPbuDhN2/A6faZ8N/7qz132rRpPP300zz22GOEhITg7u5OSUnJMX369etHcnLyqWT6PmXQilqj0WgAnMEuXAgMdyHwZheMiJxRkAkFays4vLwE6+o9eG7ciE9GMh65iXjEJ0D8EnjVGKLMyY/84IGUhw/C7czheI0dAX36nGDeDh/sC6+cDkDfkVD4w3ScOndn8T97GfvYl3S8tgdcff4x5zzyyCNs2LCB5557jl9++YWQkJBj2gsLDUdnrajbH9r0rdFoNCeLFSw7FRkrSshfkYlz7CZ8Ujfhmb8JOzZg5FE1Inas4kyBfx+Keg/CbvwQ/M6OxCF8KFQpaVnJ828nc+9d/8FpOOz5+SMcvttO52EuEG7kNC0oKGDs2LHs3r2bmTNn4ufnx3//+18AFixYQFpaGnfdddcx5TLrgTZ9tzKaTFGLyEcYHngZSqmB5jE/jFRuIRgeHBcppbLNtoeAazHi0O5QSv1uHh8OLABcgV+AO5UNQmtFrdFomp1iyFtbzsGYAvh3Kx6bNuJ7MBkXM/OaEeljdnXrTkHXgVQMG4zHucPxHD8CunY9Yjp/+MMknrnuMcQnEcnaicLQoCmpqYwYMQJHR0dmzJjB228b+dP/+ecfoqOjuffee/E4yf3z49CKupXRlKbvBcD/qBKvBswF/lJKvSAic83XD4pIf+BijHL3XYClItJbKWUB3gFuwIhj+wWYDPzahHJrNBpN/XAFr3GOeI3zhYdM0zlQvt/KweXF5P+VhnNcEj6pm/HO3UTg1g2w9Wf4wji9wt6HgoCBFPcdyPWnDWPv2/fTbXgAr/6xneumv4WHyxJCDybz9ddfc+aZZ7Jv3z6sVit2dnZHvMG16bv90WSK2swGE3Lc4ekYMW1gJDFfBjxoHl+klCoFUkRkBxBp5lD1UkqtAhCRT4Dz0Ypao9G0IRy72BF8kTtc1B8jEyVggewNZWT8dQj1TxKeyRvxObgZn4Mb8Tm4AP59FwCFI2cHjWXj0AsY1fc6Vu3Jx+vRdHYHn0kP0/P7999/Z/v27UD7VNRVi3KcijS3M1lHpVQ6gFIqXUQqE/MGYayYK9lrHis3nx9/vFpE5AaM1TfdunWrqZtGo9G0PPbgO8wJ32FBcH8QhrEQynKs7Ps3n6JfN+OyJhHv1GRCMzbhYnkOWbodx8NllJQKHUsLuPqii1j877+UlJTg5uaGs7Nzu1TU9aW2Mpdtidbi9V3dnoiq5Xi1KKXmYWSTISIion16yWk0mnaNk48dodO9YfpR03lVIjyA+EuAS3itrIw37e1PSeUsIudhZAdzAjKBOUqpg2bVqi4YvlCHReROYCFGRfU1GL+IhptVsS4D7jDHiAVuMbdcWxXNragPikhnczXdGcgwj+/FKFlWSTCw3zweXM1xjUajOeWpqZJWkyLyGk1S5lLddZLnLAdGmXm+rwMeAO4124YD45RSxSLyP+BvpdTzIjIZ0+oqIv2A2cBYpVS5iLyNUbv6kxNmamGauxr7j8CV5vMrgR+qHL9YRJxFJBToBcSZZvJ8ERll5nK9oso5Go1Gozl1CQZ+F5Ek4H4MZ+RKflRKFZvPxwGLAJRSvwGVydLPwFDoa8wSmmcAYc0g90nTZCtqEfkCw3EsQET2Ak8ALwBfici1wG6MAt8opTaJyFdAMlAB3FrF/HAzR8OzfkU7kmk0Gk3LcfIr36biTeAVpdSPIjIBeLJKW2GV5zWFmwnwsVLqoSaRrhFpSq/vmup2nlFD/2eBZ6s5Hg8MbETRNBqNRtP28Qb2mc+vrKXfcuAi4L8ichbgax7/C/hBRF5VSmWYeT48lVJpTSZxPWlu07dGo9FoNCeLm4jsrfK4B2MF/bWIxFA1k8yJPAWcJSLrgHOAdCBfKZWM4Yz2h4gkAn8CnZv0KuqJTiGq0Wg0mqq0q8xkZv1oi1KqQkRGA+8opYa2sFgnRWsJz9JoNBqNpinohuEbZQeUAde3sDwnjVbUGo1Go2m3KKW2A+EtLUdDaLembxE5BDSlU0AAte+LtBa0nI2LlrNx0XI2Lo0h52Gl1OTGEEbTOLRbRd3UiEi8UiqipeWoCy1n46LlbFy0nI1LW5FTc3Jor2+NRqPRaFoxWlFrNBqNRtOK0Yq6/sxraQFsRMvZuGg5GxctZ+PSVuTUnAR6j1qj0Wg0mlaMXlFrNBqNRtOK0Ypao9FoNJpWjFbUNiIifiLyp4hsN//6VtOnj4gkVHnkichdrU1Os5+PiHwjIltEZLOZWq81ypkqIknm/Wz2nLC2ymn2tReR9SKypDllNOe25fPpIiJxIrJBRDaJyFOtVM6uIvKP+bncJCJ3tkY5zX4fiUiGiGxsRtkmi8hWEdkhInOraRcRecNsTxSRYc0lm6Zp0IraduYCfymlemFUXTnhH0QptVUpNdTMIzscKAK+a1YpbZDT5HXgN6VUX2AIsLmZ5KvEVjkBJpr3tSXiQ09Gzjtp/vtYiS1ylgKnK6WGAEOBySIyqvlEBGyTswK4VynVDxgF3Coi/ZtRRrD9fV8ANFtyEBGxB97CKC7RH7ikmntzDtDLfNwAvNNc8mmaCKWUftjwALYCnc3nnYGtdfQ/C1jRGuUEvIAUTGfC1nw/gVQgoA3IGYzxhX46sKS1ylmlvxuwDhjZmuU0+/0AnNla5QRCgI3NJNdo4Pcqrx8CHjquz3vAJdVdi360zYdeUdtOR6VUOoD5t0Md/S8GvmhyqU7EFjnDgEPAfNNU+4GIuDenkNh+PxVGGbq1InJDs0l3FFvlfA14ALA2k1zHY5Ocpnk+AcgA/lRKxTafiMBJ/h+JSAhGnuZWLWczEgTsqfJ6r3nsZPto2hC6KEcVRGQp0KmapkdOchwnYBrGr91GpxHkdACGAbcrpWJF5HUM095jjSQi0Gj3c6xSar+IdAD+FJEtSqnoxpHQoKFyishUIEMptVZEJjSiaMfP0+D7qZSyAENFxAf4TkQGKqUadX+1Ef+PPIBvgbuUUnmNIdtx4zeKnM1MdSUoj4+xtaWPpg2hFXUVlFKTamoTkYMi0lkplS4inTFWJDVxDrBOKXWw0YWkUeTcC+ytspr6htr3XltKTpRS+82/GSLyHRAJNKqibgQ5xwLTRORcwAXwEpHPlFKXtTI5q46VIyLLMPZXG1VRN4acIuKIoaQ/V0otbkz5GlPOFmAv0LXK62Bgfz36aNoQ2vRtOz8CV5rPr8TYN6uJS2gZszfYIKdS6gCwR0T6mIfOAJKbR7wj1CmniLiLiGflc4x9/2bzrjWx5X4+pJQKVkqFYGx5/N3YStoGbLmfgeZKGhFxBSYBW5pLQBNb5BTgQ2CzUuqVZpStKifz/96crAF6iUioabm7GEPWqvwIXGF6f48CcivN+Jo2SktvkreVB+CP4Sy03fzrZx7vAvxSpZ8bkAl4t3I5hwLxQCLwPeDb2uTE2EvfYD42AY+01vtZpf8EWsaZzJb7ORhYb77nG4HHW6mc4zBMtYlAgvk4t7XJab7+AkgHyjFWstc2g2znAtuAnZX/E8BNwE3mc8HwDN8JJAERzf0+60fjPnQKUY1Go9FoWjHa9K3RaDQaTStGK2qNRqPRaFoxWlFrNBqNRtOK0Ypao9FoNJpWjFbUGo1Go9G0YrSi1mg0Go2mFaMVtUajadWISJiIfCgi37S0LBpNS6AVtaZWRGSZiJx93LG7ROTtkxjDR0RuOe7YynrKU1DDcYsYNas3isjXIuJWn/Gbk/re25rugY1zrjT/nvCe2Hi+q4j8axb3aJZ7rpTapZS69jg5nEQkWkR0GmRNu0crak1dfIGRprAqNlUGM1MY2gE+wDFKQSk1prEENClWRs3qgUAZRqam6mRpck5irnrf2/pS5b77cNx7YiPXAIuVUdyj1nt+sojIIBFZctyj2qpVSqkyjIxhsxsyp0bTFtCKWlMX3wBTRcQZjpQd7AIsF5HLRCTOXFW9Z66yQkRks7kqXIdRHOAFoIfZ70VznALz7xUikigiG0Tk08pJReR7MUpbbpKTL28ZA/SsTpYaZHYXkZ9NGTaKyOzqjlVev4gcyTcuIveJyJO2zmXrvTVf13U+InKPKd9GEbmryvGa7mvlavyY90RE7qzS51kRuaOGezuH6vNexwA9zfPrfO+qu79KqSSl1NTjHrUVw/jelEejad+0dA5T/Wj9D+BnYLr5fC7wItAP+AlwNI+/DVwBhGDUZB5V5fwQYONxYxYAAzCK2geYx/yqtFfmVnbFyEvtX3leDTIWmH8dMBTJzcfLUovMFwLvVxnLu7pj1V0LcB/wpK1z2XJv6zq/yrUOx8jl7A54YORDD6/jvhYcfx3m83XmczuMHNH+1cjqBByo7Z7X9t4dN1a197eG99YfeNeU66Eqx+2BQy39/6Ef+tHUD72/o7GFShPtD+bfazAqbg0H1ogIGF/KGRglKNOUUqttGPd04Bul1GEApVRWlbY7RGSG+bwr0Auj2ElNuIpIgvk8BqP6UpfjZKlJ5oXASyLyX4yCGjEiknT8MRuux5a5jqe6e2vr+eOA75RShQAishiIwihoUdN9PQGlVKqIZIpIONARWK+Uqu5eBwA5VV5Xd8/BtvfO5vtrynKCWV0pZRGRMhHxVErl13aNGk1bRitqjS18D7wiIsMAV6XUOhEZC3yslHqoakfTfFto47hCNQXtRWQCRgnG0UqpIjHqJrvUMVaxUmroceNwnCxSncxm3+EYVYmeF5E/lFL/qe4YUMGxW0ZV5bJpruP4nuPu7UmcL7UcP9lqOx8AVwGdgI9q6FPMsddb3T2fgA3vnVJqWw3392RxBkrqcZ5G02bQe9SaOlFKFQDLML7AKx2d/gJmVjr7iIifiHSvYYh8wLOa438BF4mIf+UY5nFvINv8ou8LjGqUC6lBZhHpAhQppT4DXgKGVXfMHOMg0EFE/M295aknM9fxnWq4t7aeHw2cLyJuYtTrnoGxsq3pvlbl+PfkO2AyMAL4vboLUkplA/YiUtuPJpveu1rur82Y13dIKVV+sudqNG0JvaLW2MoXwGJML2WlVLKIPAr8IYaHczlwK3Dg+BOVUpkissJ0wvpVKXW/eXyTiDwL/CsiFoxayVcBvwE3iUgixl6rLWb0OqlFZm/gRRGxmsduBgZVcwylVLmI/AeIBVKALSc5V1o13Y+5t7aeb1o2FgBx5qEPlFLrwXAI48T7WlW+E94TEfkHyFGGR3dN/IFhcl9aQ7ut71219/ckmQj8Uo/zNJo2ha5HrdFoMH8MrANmKaW219IvHLhHKXV5swlXsyyLMZzLtra0LBpNU6JN3xrNKY6I9Ad2AH/VpqQBzBX7P9WFijUnIuIEfK+VtOZUQK+oNRqNRqNpxegVtUaj0Wg0rRitqDUajUajacVoRa3RaDQaTStGK2qNRqPRaFoxWlFrNBqNRtOK0Ypao9FoNJpWjFbUGo1Go9G0YrSi1mg0Go2mFaMVtUaj0Wg0rZj/B4wNepH6IXkSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plev = wap1_trop['plev'] * 0.01 # always go back to orig to re-scale\n",
    "color0 = 'gray'\n",
    "color1 = 'black'\n",
    "color2 = 'C0'\n",
    "ax.plot(wap1_mean, plev, color=color0, label='mean')\n",
    "ax.plot(wap1_conv, plev, color=color1, label='convective')\n",
    "ax.plot(wap1_extr, plev, color=color2, label='P$_{\\mathrm{extreme}}$')\n",
    "ax.plot(wap1_long_events_mean, plev, color='magenta', label='Persistent')\n",
    "ax.plot(wap1_huge_events_mean, plev, color='red', label='Large')\n",
    "\n",
    "linsty = (0, (.25,2))\n",
    "ax.plot(wap2_mean, plev, color=color0, linestyle=linsty, dash_capstyle = \"round\")\n",
    "ax.plot(wap2_conv, plev, color=color1, linestyle=linsty, dash_capstyle = \"round\")\n",
    "ax.plot(wap2_extr, plev, color=color2, linestyle=linsty, dash_capstyle = \"round\")\n",
    "ax.plot(wap2_long_events_mean, plev, color='magenta', linestyle=linsty, dash_capstyle = \"round\")\n",
    "ax.plot(wap2_huge_events_mean, plev, color='red', linestyle=linsty, dash_capstyle = \"round\")\n",
    "\n",
    "\n",
    "ax.invert_yaxis()\n",
    "ax.set_ylabel(\"Pressure (hPa)\")\n",
    "ax.set_xlabel(\"Vertical Pressure Velocity (Pa s$^{-1}$)\")\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.set_title(f\"{tup1[0]}, {tup1[1]}\", loc='left')\n",
    "fig.legend(bbox_to_anchor=(0.9, .1), loc='lower left', frameon=False)\n",
    "figname = f\"/Users/brianpm/Dropbox/Manuscripts/M_cre_pex/wap_profiles_{tup1[0]}_{tup1[1]}_{tup2[1]}.pdf\"\n",
    "fig.savefig(figname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}