{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab 7-3 - Predicting snowmelt rates with the Temperature-Index method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lab, we will use SWE and air temperature measurements from the East River Valley (SNOTEL sites, plus the Kettle Ponds measurements)." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "from metloom.pointdata import SnotelPointData\n", "import altair as alt\n", "import numpy as np\n", "from metpy.units import units\n", "import numpy as np\n", "import altair as alt\n", "alt.data_transformers.disable_max_rows()\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "start_date = '1990-01-01'\n", "end_date = '2024-01-10'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Download SNOTEL SWE data" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
datetimesitegeometrySWESWE_unitsMAX AIR TEMPMAX AIR TEMP_unitsPRECIPITATIONPRECIPITATION_unitsdatasource
01990-01-01 08:00:00+00:00380:CO:SNTLPOINT Z (-106.95327 38.89435 10190.00000)2.3in26.42degF0.0inNRCS
11990-01-02 08:00:00+00:00380:CO:SNTLPOINT Z (-106.95327 38.89435 10190.00000)2.3in23.54degF0.0inNRCS
21990-01-03 08:00:00+00:00380:CO:SNTLPOINT Z (-106.95327 38.89435 10190.00000)2.3in18.50degF0.0inNRCS
31990-01-04 08:00:00+00:00380:CO:SNTLPOINT Z (-106.95327 38.89435 10190.00000)2.3in12.02degF0.1inNRCS
41990-01-05 08:00:00+00:00380:CO:SNTLPOINT Z (-106.95327 38.89435 10190.00000)2.4in19.58degF0.0inNRCS
\n", "
" ], "text/plain": [ " datetime site \\\n", "0 1990-01-01 08:00:00+00:00 380:CO:SNTL \n", "1 1990-01-02 08:00:00+00:00 380:CO:SNTL \n", "2 1990-01-03 08:00:00+00:00 380:CO:SNTL \n", "3 1990-01-04 08:00:00+00:00 380:CO:SNTL \n", "4 1990-01-05 08:00:00+00:00 380:CO:SNTL \n", "\n", " geometry SWE SWE_units MAX AIR TEMP \\\n", "0 POINT Z (-106.95327 38.89435 10190.00000) 2.3 in 26.42 \n", "1 POINT Z (-106.95327 38.89435 10190.00000) 2.3 in 23.54 \n", "2 POINT Z (-106.95327 38.89435 10190.00000) 2.3 in 18.50 \n", "3 POINT Z (-106.95327 38.89435 10190.00000) 2.3 in 12.02 \n", "4 POINT Z (-106.95327 38.89435 10190.00000) 2.4 in 19.58 \n", "\n", " MAX AIR TEMP_units PRECIPITATION PRECIPITATION_units datasource \n", "0 degF 0.0 in NRCS \n", "1 degF 0.0 in NRCS \n", "2 degF 0.0 in NRCS \n", "3 degF 0.1 in NRCS \n", "4 degF 0.0 in NRCS " ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "snotel_point_butte = SnotelPointData(\"380:CO:SNTL\", \"Butte\")\n", "SNOTEL_VARS = [\n", " snotel_point_butte.ALLOWED_VARIABLES.SWE,\n", " snotel_point_butte.ALLOWED_VARIABLES.TEMPMAX,\n", " snotel_point_butte.ALLOWED_VARIABLES.PRECIPITATION\n", "]\n", "df_butte_longterm = snotel_point_butte.get_daily_data(\n", " pd.to_datetime(start_date), pd.to_datetime(end_date),\n", " SNOTEL_VARS\n", ")\n", "df_butte_longterm = df_butte_longterm.reset_index()\n", "df_butte_longterm.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add columns for year, and then a date object without a year specified (or, actually, where all the years have been replaced with the year 2000). The year column is useful for grouping the datasets into separate years. The \"date_no_year\" column is useful for plotting SWE evolution from each year on the same chart." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "df_butte_longterm['year'] = df_butte_longterm['datetime'].dt.year\n", "df_butte_longterm['date_no_year'] = df_butte_longterm['datetime'].apply(\n", " lambda dt: dt.replace(year=1999) if dt.month in [10,11,12] else dt.replace(year=2000)\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot annual SWE pattern for all years" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ], "text/plain": [ "alt.Chart(...)" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "alt.Chart(df_butte_longterm).mark_line().encode(\n", " alt.X('date_no_year'),\n", " alt.Y('SWE:Q'),\n", " alt.Color('year:Q')\n", ").properties(width=600)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Calculate melt rate from daily ∆ SWE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We calculate \n", "$$Melt Rate = - ∆SWE$$\n", "where ∆SWE is calculated as SWE on the second day minus SWE on the first day.\n", "\n", "To do the calculation with our dataframe, we first group by year, then calculate the difference in SWE on consecutive days, using the `diff` function." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/x_/2h52bcjx2px15bhmdpdd748h0000gn/T/ipykernel_45676/3218391937.py:1: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.\n", " delta_swe = df_butte_longterm.groupby('year').apply(lambda df: df.set_index('datetime')[['SWE']].diff())\n" ] } ], "source": [ "delta_swe = df_butte_longterm.groupby('year').apply(lambda df: df.set_index('datetime')[['SWE']].diff())\n", "delta_swe = delta_swe.rename(columns={'SWE': 'Delta SWE'})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot above shows us that the ablation period usually doesn't usually start until April 1st, so we will focus on data from after April 1st when we calculate melt rates. Also, we see the snow is always gone before the end of June. So we remove data from before April and after June." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Delta SWE
yeardatetime
19901990-04-01 08:00:00+00:000.0
1990-04-02 08:00:00+00:000.0
1990-04-03 08:00:00+00:000.0
1990-04-04 08:00:00+00:000.0
1990-04-05 08:00:00+00:000.0
\n", "
" ], "text/plain": [ " Delta SWE\n", "year datetime \n", "1990 1990-04-01 08:00:00+00:00 0.0\n", " 1990-04-02 08:00:00+00:00 0.0\n", " 1990-04-03 08:00:00+00:00 0.0\n", " 1990-04-04 08:00:00+00:00 0.0\n", " 1990-04-05 08:00:00+00:00 0.0" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "delta_swe = delta_swe[delta_swe.index.get_level_values(1).month > 3] \n", "delta_swe = delta_swe[delta_swe.index.get_level_values(1).month < 7]\n", "delta_swe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also remove days where SWE increased. We don't want to estimate a melt rate relationship on days where we did not observe melting. Also, let's calculate Melt Rate by simply negative the ∆SWE." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "delta_swe = delta_swe[delta_swe['Delta SWE'] < 0]\n", "\n", "delta_swe['MELT RATE'] = - delta_swe['Delta SWE']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok we now have a dataset of melt rates during snowpack ablation, in units of inches/day, for 20 years of data, in a nice table." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Delta SWEMELT RATE
yeardatetime
19901990-04-12 08:00:00+00:00-0.30.3
1990-04-13 08:00:00+00:00-0.10.1
1990-04-15 08:00:00+00:00-0.30.3
1990-04-16 08:00:00+00:00-0.40.4
1990-04-17 08:00:00+00:00-0.40.4
\n", "
" ], "text/plain": [ " Delta SWE MELT RATE\n", "year datetime \n", "1990 1990-04-12 08:00:00+00:00 -0.3 0.3\n", " 1990-04-13 08:00:00+00:00 -0.1 0.1\n", " 1990-04-15 08:00:00+00:00 -0.3 0.3\n", " 1990-04-16 08:00:00+00:00 -0.4 0.4\n", " 1990-04-17 08:00:00+00:00 -0.4 0.4" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "delta_swe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's add columns of maximum daily air temperature and precipitation, from the Snotel dataset that we initially downloaded. We can join our two tables on the datetime index." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
yearDelta SWEMELT RATEMAX AIR TEMPPRECIPITATION
datetime
1990-04-12 08:00:00+00:001990-0.30.347.300.1
1990-04-13 08:00:00+00:001990-0.10.140.460.0
1990-04-15 08:00:00+00:001990-0.30.351.980.0
1990-04-16 08:00:00+00:001990-0.40.454.500.1
1990-04-17 08:00:00+00:001990-0.40.445.140.1
\n", "
" ], "text/plain": [ " year Delta SWE MELT RATE MAX AIR TEMP \\\n", "datetime \n", "1990-04-12 08:00:00+00:00 1990 -0.3 0.3 47.30 \n", "1990-04-13 08:00:00+00:00 1990 -0.1 0.1 40.46 \n", "1990-04-15 08:00:00+00:00 1990 -0.3 0.3 51.98 \n", "1990-04-16 08:00:00+00:00 1990 -0.4 0.4 54.50 \n", "1990-04-17 08:00:00+00:00 1990 -0.4 0.4 45.14 \n", "\n", " PRECIPITATION \n", "datetime \n", "1990-04-12 08:00:00+00:00 0.1 \n", "1990-04-13 08:00:00+00:00 0.0 \n", "1990-04-15 08:00:00+00:00 0.0 \n", "1990-04-16 08:00:00+00:00 0.1 \n", "1990-04-17 08:00:00+00:00 0.1 " ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "delta_swe = delta_swe.reset_index().set_index('datetime').join(\n", " df_butte_longterm.set_index('datetime')[['MAX AIR TEMP', 'PRECIPITATION']]\n", ")\n", "delta_swe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's convert our units to metric. Currently, the SWE column is mm/day and the MAX AIR TEMP column is ˚F." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "delta_swe['Delta SWE'] = (delta_swe['Delta SWE'].values * units('inches')).to(units('millimeters')).magnitude\n", "delta_swe['MAX AIR TEMP'] = (delta_swe['MAX AIR TEMP'].values * units('degF')).to(units('degC')).magnitude\n", "delta_swe['PRECIPITATION'] = (delta_swe['PRECIPITATION'].values * units('inches')).to(units('millimeters')).magnitude\n", "delta_swe['MELT RATE'] = (delta_swe['MELT RATE'].values * units('inches')).to(units('millimeters')).magnitude" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's plot the relationship between ∆SWE and max daily air temp. There are lots of data points here, so let's plot a scatterplot as well as a \"2d histogram\" representing point density." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ], "text/plain": [ "alt.HConcatChart(...)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "chart_scatter = alt.Chart(delta_swe).mark_point(size=5).encode(\n", " alt.X('MAX AIR TEMP:Q').title('Max daily air temp. (˚C)'),\n", " alt.Y('MELT RATE:Q').title('Melt rate (mm/day)')\n", ")\n", "chart_2d_histogram = alt.Chart(delta_swe).mark_rect().encode(\n", " alt.X('MAX AIR TEMP:Q', bin=alt.Bin(maxbins=30)).title('Max daily air temp. (˚C)'),\n", " alt.Y('MELT RATE:Q', bin=alt.Bin(maxbins=30)).title('Melt rate (mm/day)'),\n", " alt.Color('count()').title('Point Density')\n", ")\n", "(chart_scatter | chart_2d_histogram).resolve_scale(x='shared', y='shared')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Fit a Temperature-Index equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can fit a line to our data,\n", "\n", "$$ M = M_f (T_a - T_0), $$\n", "\n", "and find the find the degree day factor $M_f$. Note that because we are assuming $T_0 = 0$, we can just ignore it and fit a line assuming a zero intercet,\n", "\n", "$$ M = M_f (T_a). $$\n", "\n", "Considering the data above, we will remove outliers where the max daily air temp is greater than 25˚C, because these do not seem to follow the major pattern." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# filter data points\n", "src = delta_swe[delta_swe['MAX AIR TEMP'] < 25]\n", "\n", "# Perform linear regression with zero intercept\n", "x = src['MAX AIR TEMP']\n", "y = src['MELT RATE']\n", "DEGREE_DAY_FACTOR = np.sum(x * y) / np.sum(x * x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the line on top of the charts above!" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ], "text/plain": [ "alt.HConcatChart(...)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create the plot for the line of best fit\n", "# Line of best fit using the slope (DEGREE_DAY_FACTOR) and zero intercept\n", "line = alt.Chart(src).mark_line(color='red').encode(\n", " alt.X('MAX AIR TEMP:Q'),\n", " alt.Y('MELT RATE:Q')\n", ").transform_calculate(\n", " \"MELT RATE\", f\"{DEGREE_DAY_FACTOR} * datum['MAX AIR TEMP']\"\n", ")\n", "\n", "# Add it to the charts above.\n", "# Add the estimated degree day factor (the DEGREE_DAY_FACTOR of the line of best fit)\n", "# to the charts.\n", "(\n", " (chart_scatter + line).properties(title=f'M_f = {round(DEGREE_DAY_FACTOR,2)} mm/˚C/day') \n", " | \n", " (chart_2d_histogram + line).properties(title=f'M_f = {round(DEGREE_DAY_FACTOR,2)} mm/˚C/day')\n", ").resolve_scale(x='shared', y='shared')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Predict snowpack ablation for one year" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using this Temperature-Index model, we can predict snowmelt for a given year, even if we only have an air temperature measurement.\n", "Let's try to predict snow melt at the Butte Snotel site using the fitted Temperature-Index model we created above. \n", "We will do this for the 2022-2023 winter season (the same year as the SOS campaign).\n", "\n", "Using the Temperature-Index model/equation, we can simply estimate snowmelt for a given day. For this modeling exercise, we will \"initialize the model\" with the *measured* April 1 SWE. We will then apply the model (i.e. estimate daily snow melt) using the measured air temps. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's grab the data we want for the single snow season." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# Grab Snotel-measured data from the 2022-23 water year\n", "df_butte_2023 = df_butte_longterm.set_index('datetime').loc['20230401': '20230630']\n", "\n", "## convert SWE from inches to mm and air temp from F to C\n", "df_butte_2023['SWE'] = (df_butte_2023['SWE'].values * units(\"inches\")).to(units(\"mm\"))\n", "df_butte_2023['MAX AIR TEMP'] = (df_butte_2023['MAX AIR TEMP'].values * units(\"degF\")).to(units(\"degC\"))\n", "df_butte_2023['PRECIPITATION'] = (df_butte_2023['PRECIPITATION'].values * units(\"inches\")).to(units(\"mm\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can calculate \"predicted melt rate\" simply from air temperature, and the DEGREE_DAY_FACTOR we calculated above, from 20 years of historical data." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# calculate melt rate from max air temp. Use the DEGREE_DAY_FACTOR calculated above\n", "df_butte_2023['melt rate'] = df_butte_2023['MAX AIR TEMP'] * DEGREE_DAY_FACTOR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we will use daily melt rates to model snowpack ablation. Recall, we need to \"initialize\" the modeled snowpack - we do this with measured SWE on April 1. Then, we will loop over each day, applying the estimated melt rate to the modeled snowpack." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# get the initial SWE value from measurements\n", "swe_measured_initial_value = df_butte_2023['SWE'].iloc[0]\n", "\n", "# calculate modeled SWE, modeling ablation using the estimated melt rates\n", "# First, create a list that will hold the modeled, daily SWE\n", "modeled_swe = [swe_measured_initial_value]\n", "\n", "# Second, loop over melt rates from each day, calculate the new\n", "# SWE using the melt rate from the previous day.\n", "for melt_rate in df_butte_2023['melt rate']:\n", " updated_swe = modeled_swe[-1] - melt_rate\n", " if updated_swe >= 0:\n", " modeled_swe.append(updated_swe)\n", " else:\n", " modeled_swe.append(0)\n", "\n", "# Add the modeled SWE back into our dataframe. We remove the last predicted SWE value because we have one day of modeled SWE beyond our measured SWE.\n", "df_butte_2023['SWE predicted'] = modeled_swe[:-1]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAHWCAYAAAAciQ/OAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcp5JREFUeJzt3Xd0FOXbxvHvppJAEhJKQiBA6L0jVYo0ERBsgAiIgqAiEFBRXkVFFBQVUREU9QcoUhQFUUFBBaQqBELvvSSEElIgpD7vHyurkR6zu9nk+pyzx+zszM69d8blypRnLMYYg4iIiIi4JDdnFyAiIiIi2acwJyIiIuLCFOZEREREXJjCnIiIiIgLU5gTERERcWEKcyIiIiIuTGFORERExIUpzImIiIi4MA9nF+BMmZmZnDx5Ej8/PywWi7PLEREREbExxpCYmEhoaChubtfe/5avw9zJkycJCwtzdhkiIiIi13Ts2DFKlSp1zdfzdZjz8/MDrE3y9/d3cjUiIiIif0tISCAsLMyWV64lX4e5y4dW/f39FeZEREQkV7rRqWC6AEJERETEhSnMiYiIiLgwhTkRERERF5Zrw9zvv/9Oly5dCA0NxWKxsHDhwiyvG2N45ZVXCA0NxcfHh1atWrFjxw7nFCsiIiLiJLk2zF24cIHatWszefLkq74+YcIEJk6cyOTJk9mwYQMhISG0a9eOxMREB1cqIiIi4jy59mrWjh070rFjx6u+Zoxh0qRJvPDCC9x7770AzJw5k+DgYGbPns2gQYMcWapcRdnnf7T9fPiNTvlm3eI4+j2LiFjl2j1z13Po0CFiYmJo3769bZq3tzctW7Zk7dq111wuJSWFhISELA8RERERV+aSYS4mJgaA4ODgLNODg4Ntr13N+PHjCQgIsD109wcRERFxdS4Z5i779yB6xpjrDqw3atQo4uPjbY9jx47Zu0QRERERu8q158xdT0hICGDdQ1eiRAnb9NjY2Cv21v2Tt7c33t7edq9PRERExFFccs9ceHg4ISEhLFu2zDYtNTWVlStX0rRpUydWJiIiIuJYuXbPXFJSEvv377c9P3ToEFFRUQQFBVG6dGkiIiIYN24cFStWpGLFiowbNw5fX1969erlxKpFREREHCvXhrmNGzfSunVr2/MRI0YA8PDDDzNjxgxGjhxJcnIyTz75JHFxcTRq1IilS5fi5+fnrJJFREREHC7XhrlWrVphjLnm6xaLhVdeeYVXXnnFcUWJiIiI5DIuec6ciIiIiFgpzImIiIi4MIU5ERERERemMCciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmMKciIiIiAtTmBMRERFxYQpzIiIiIi5MYU5ERETEhSnMiYiIiLgwhTkRERERF6YwJyIiIuLCXDbMpaen8+KLLxIeHo6Pjw/lypXj1VdfJTMz09mliYiIiDiMh7MLyK4333yTjz76iJkzZ1K9enU2btzII488QkBAAMOGDXN2eSIiIiIO4bJhbt26dXTt2pVOnToBULZsWebMmcPGjRudXJmIiIiI47jsYdbmzZvz66+/snfvXgC2bNnC6tWrueuuu5xcmYiIiIjjuOyeueeee474+HiqVKmCu7s7GRkZvP766zz44IPXXCYlJYWUlBTb84SEBEeUKiIiImI3Lrtnbt68ecyaNYvZs2ezadMmZs6cydtvv83MmTOvucz48eMJCAiwPcLCwhxYsYiIiEjOc9kw9+yzz/L888/Ts2dPatasSZ8+fRg+fDjjx4+/5jKjRo0iPj7e9jh27JgDKxYRERHJeS57mPXixYu4uWXNou7u7tcdmsTb2xtvb297lyYiIiLiMC4b5rp06cLrr79O6dKlqV69Ops3b2bixIk8+uijzi5NRERExGFcNsx98MEHjB49mieffJLY2FhCQ0MZNGgQL730krNLExEREXEYlw1zfn5+TJo0iUmTJjm7FBERERGncdkwJ7lD2ed/tP18+I1OTqzkxuxdqyv1Ii/T70HEMfT/Wu7hslezioiIiIjCnIiIiIhLU5gTERERcWEKcyIiIiIuTGFORERExIUpzImIiIi4MIU5ERERERemMCciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmMKciIiIiAtTmBMRERFxYQpzIiIiIi7MpcPciRMn6N27N0WKFMHX15c6deoQGRnp7LJEREREHMbD2QVkV1xcHM2aNaN169YsWbKE4sWLc+DAAQoXLuzs0kREREQcxmXD3JtvvklYWBjTp0+3TStbtqzzChIRERFxApc9zLpo0SIaNGjAAw88QPHixalbty6ffPLJdZdJSUkhISEhy0NERETElblsmDt48CBTp06lYsWK/Pzzzzz++OMMHTqUzz///JrLjB8/noCAANsjLCzMgRWLiIiI5DyXDXOZmZnUq1ePcePGUbduXQYNGsRjjz3G1KlTr7nMqFGjiI+Ptz2OHTvmwIpFREREcp7LhrkSJUpQrVq1LNOqVq3K0aNHr7mMt7c3/v7+WR4iIiIirsxlw1yzZs3Ys2dPlml79+6lTJkyTqpIRERExPFcNswNHz6c9evXM27cOPbv38/s2bOZNm0agwcPdnZpIiIiIg7jsmGuYcOGLFiwgDlz5lCjRg3Gjh3LpEmTeOihh5xdmoiIiIjDuOw4cwCdO3emc+fOzi5DRERExGlcds+ciIiIiCjMiYiIiLg0hTkRERERF6YwJyIiIuLCFOZEREREXJjCnIiIiIgLU5gTERERcWEuPc5cTqnx8s+4efsCcPiNTk6uRkREclrZ53+0/azveclrtGdORERExIUpzImIiIi4MIU5ERERERemMCciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sLyTJgbP348FouFiIgIZ5ciIiIi4jB5Isxt2LCBadOmUatWLWeXIiIiIuJQLh/mkpKSeOihh/jkk08IDAx0djkiIiIiDuXyYW7w4MF06tSJtm3b3nDelJQUEhISsjxEREREXJmHswv4L+bOncumTZvYsGHDTc0/fvx4xowZY+eqRERERBzHZffMHTt2jGHDhjFr1iwKFChwU8uMGjWK+Ph42+PYsWN2rlJERETEvlx2z1xkZCSxsbHUr1/fNi0jI4Pff/+dyZMnk5KSgru7e5ZlvL298fb2dnSpIiIiInbjsmGuTZs2bNu2Lcu0Rx55hCpVqvDcc89dEeRERERE8iKXDXN+fn7UqFEjy7SCBQtSpEiRK6aLiIiI5FUue86ciIiIiLjwnrmrWbFihbNLEBEREXEo7ZkTERERcWEKcyIiIiIuTGFORERExIUpzImIiIi4MIU5ERERERemMCciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmMKciIiIiAvzcHYBkr+Uff5H28+H3+jkxEpu7Fq12vsz/PP97bWOW+VKv7dblZc/261SL+RqbnW70HbkeNozJyIiIuLCFOZEREREXJjCnIiIiIgLc9kwN378eBo2bIifnx/FixenW7du7Nmzx9lliYiIiDiUy4a5lStXMnjwYNavX8+yZctIT0+nffv2XLhwwdmliYiIiDiMy17N+tNPP2V5Pn36dIoXL05kZCQtWrRwUlUiIiIijuWye+b+LT4+HoCgoCAnVyIiIiLiOC67Z+6fjDGMGDGC5s2bU6NGjWvOl5KSQkpKiu15QkKCI8oTERERsZs8sWfuqaeeYuvWrcyZM+e6840fP56AgADbIywszEEVioiIiNiHy4e5IUOGsGjRIpYvX06pUqWuO++oUaOIj4+3PY4dO+agKkVERETsw2UPsxpjGDJkCAsWLGDFihWEh4ffcBlvb2+8vb0dUJ2IiIiIY7hsmBs8eDCzZ8/mu+++w8/Pj5iYGAACAgLw8fFxcnUiIiIijuGyh1mnTp1KfHw8rVq1okSJErbHvHnznF2aiIiIiMO47J45Y4yzSxARERFxOpfdMyciIiIiCnMiIiIiLk1hTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmMKciIiIiAtTmBMRERFxYQpzIiIiIi5MYU5ERETEhSnMiYiIiLgwhTkRERERF6YwJyIiIuLCFOZEREREXJjCnIiIiIgLU5gTERERcWEKcyIiIiIuTGFORERExIUpzImIiIi4MJcPc1OmTCE8PJwCBQpQv359Vq1a5eySRERERBzGpcPcvHnziIiI4IUXXmDz5s3cfvvtdOzYkaNHjzq7NBERERGHcOkwN3HiRPr378+AAQOoWrUqkyZNIiwsjKlTpzq7NBERERGHcNkwl5qaSmRkJO3bt88yvX379qxdu9ZJVYmIiIg4loezC8iuM2fOkJGRQXBwcJbpwcHBxMTEXHWZlJQUUlJSbM/j4+MByEy5aJuWkJBgh2rzrmv17lanO7OmnHqfnPps/3yf//peOSU3/j9i799zfpSXe5GXP5u95cbv+fzicv+MMdef0bioEydOGMCsXbs2y/TXXnvNVK5c+arLvPzyywbQQw899NBDDz30cJnHsWPHrpuJXHbPXNGiRXF3d79iL1xsbOwVe+suGzVqFCNGjLA9z8zM5Ny5cxQpUgSLxWLXekVERERuhTGGxMREQkNDrzufy4Y5Ly8v6tevz7Jly7jnnnts05ctW0bXrl2vuoy3tzfe3t5ZphUuXNieZYqIiIhkW0BAwA3ncdkwBzBixAj69OlDgwYNaNKkCdOmTePo0aM8/vjjzi5NRERExCFcOsz16NGDs2fP8uqrrxIdHU2NGjVYvHgxZcqUcXZpIiIiIg5hMeZGl0iIiOQNrVq1ok6dOkyaNClfrVtE8jaXHWdORMSeVqxYgcVi4fz58zmy3LfffsvYsWNzrkARkb+49GFWERFXERQU5OwSRCSP0p45EcmTLly4QN++fSlUqBAlSpTgnXfeyfL6rFmzaNCgAX5+foSEhNCrVy9iY2MBOHz4MK1btwYgMDAQi8VCv379AOtQARMmTKBcuXL4+PhQu3Zt5s+ff8PlWrVqRUREhG39ZcuW5bXXXrPVWKZMGb777jtOnz5N165dKVSoEDVr1mTjxo1Z6l67di0tWrTAx8eHsLAwhg4dyoULF3K6fSLiQhTmRCRPevbZZ1m+fDkLFixg6dKlrFixgsjISNvrqampjB07li1btrBw4UIOHTpkC15hYWF88803AOzZs4fo6Gjee+89AF588UWmT5/O1KlT2bFjB8OHD6d3796sXLnyustdzbvvvkuzZs3YvHkznTp1ok+fPvTt25fevXuzadMmKlSoQN++fW2jv2/bto0OHTpw7733snXrVubNm8fq1at56qmn7NFCEXEVOXE3BhGR3CQxMdF4eXmZuXPn2qadPXvW+Pj4mGHDhl11mT///NMAJjEx0RhjzPLlyw1g4uLibPMkJSWZAgUKXHHnmf79+5sHH3zwmssZY0zLli2zrLtMmTKmd+/etufR0dEGMKNHj7ZNW7dunQFMdHS0McaYPn36mIEDB2Z531WrVhk3NzeTnJx8/aaISJ6lc+ZEJM85cOAAqampNGnSxDYtKCiIypUr255v3ryZV155haioKM6dO0dmZiYAR48epVq1ald93507d3Lp0iXatWuXZXpqaip169a95Tpr1apl+/nynWtq1qx5xbTY2FhCQkKIjIxk//79fPnll7Z5jDFkZmZy6NAhqlatess1iIjrU5gTkTzH3GDEpQsXLtC+fXvat2/PrFmzKFasGEePHqVDhw6kpqZec7nLge/HH3+kZMmSWV77991lboanp6ft58u3FLzatMvrzczMZNCgQQwdOvSK9ypduvQtr19E8gaFORHJcypUqICnpyfr16+3hZy4uDj27t1Ly5Yt2b17N2fOnOGNN94gLCwM4IoLDby8vADIyMiwTatWrRre3t4cPXqUli1bXnXdV1sup9SrV48dO3ZQoUKFHH9vEXFdugBCRPKcQoUK0b9/f5599ll+/fVXtm/fTr9+/XBzs37llS5dGi8vLz744AMOHjzIokWLrhgDrkyZMlgsFn744QdOnz5NUlISfn5+PPPMMwwfPpyZM2dy4MABNm/ezIcffsjMmTOvuVxOee6551i3bh2DBw8mKiqKffv2sWjRIoYMGZJj6xAR16MwJyJ50ltvvUWLFi24++67adu2Lc2bN6d+/foAFCtWjBkzZvD1119TrVo13njjDd5+++0sy5csWZIxY8bw/PPPExwcbLtidOzYsbz00kuMHz+eqlWr0qFDB77//nvCw8Ovu1xOqFWrFitXrmTfvn3cfvvt1K1bl9GjR1OiRIkcW4eIuB7dzktERETEhWnPnIiIiIgLU5gTERERcWEKcyIiIiIuTGFORERExIUpzImIiIi4MIU5ERERERemMCciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhTkRERMSFKcyJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmIezC3CmzMxMTp48iZ+fHxaLxdnliIiIiNgYY0hMTCQ0NBQ3t2vvf8vXYe7kyZOEhYU5uwwRERGRazp27BilSpW65uv5Osz5+fkB1ib5+/s7uRoRERGRvyUkJBAWFmbLK9eSr8Pc5UOr/v7+CnMiIiKSK93oVDBdACEiIiLiwhTmRERERFyYwpyIiIiIC8vX58w5izGGhEvpGGOcXcpVWSwW/At4aLgWEZFcJiMjg7S0NGeXITnE09MTd3f3//w+CnMOEhN/iXUHz7B2/1nWHjjLifPJzi7punw83QkL8iEs0JdSgT6EBflSzM/7ivksFgvlixWkaog/bm4KfyIi9mCMISYmhvPnzzu7FMlhhQsXJiQk5D/tQFGYs6M9x6LxnPcgy1Jr8VViDQ6YUMA1Ak9yWgZ7TyWx91TSTc1f2NeTJuWK0LR8EZqUL0r5YgW1Z09EJIdcDnLFixfH19dX3695gDGGixcvEhsbC0CJEiWy/V4Kc3aUue83yiVtZhCbGeQNJ91COV68BV7VO1OhflsKeF+5pys3SM80RMdf4ti5ixyLu8ixc8kci7tI3IXUK+fNMOw4Gc/5i2ks2R7Dku0xAJQIKMCA28vxUKPSFPD877uQRUTyq4yMDFuQK1KkiLPLkRzk4+MDQGxsLMWLF8/2IVeLya0nbjlAQkICAQEBxMfH22WcuUtxJ/n1209plPoHRc78iSXjH2HIJxDqPAQN+0NQuRxftyOlZWSy9Xg86w6cYe2Bs2w8EkdqeiYAIf4FeOqOCnRvEIaXR85db5OZaYhNTPkrbF7kUlombasWp7h/gRxbh4hIbnDp0iUOHTpE2bJlbf/4S96RnJzM4cOHCQ8Pp0CBrP+G3WxOUZizY5jLIiURDvwGe36CvT9B8rm/XrBAxXZw2yAofwdc595rruJSWgYLN5/g/V/3cTL+EgBhQT4Ma1OJbnVC8XDP3mf889A5Plt9kH2nkjh+PtkWGC9zd7PQunJxejQMo3XlYtlej4hIbnI5zF3tH3txfdf7/SrM3QSHhrl/ysyA/b/Cn9Ng/7K/pweVg7q9raEupPbNBbsLZ+D07r8ee+DsfnD3gkLFoVAI+AVDoWAoXAZCaoIDz7NISc9g7p/HmLx8P6cTUwAoV6wgw9tWolPNEjd9wcSWY+d5e+keVu07k2W6u5uFEgEFCAv0JTktg6hj522vFfPz5r56pXioUWnCgnxz7DOJiDiawlze5pJh7pVXXmHMmDFZpgUHBxMTYz3XyhjDmDFjmDZtGnFxcTRq1IgPP/yQ6tWr2+ZPSUnhmWeeYc6cOSQnJ9OmTRumTJly3ZvQXo3Twtw/nT0AGz6FzbMgJeHv6T5BEN4CyrWCMs2se/bOHYS4Q9b/njsIZ/b9Yw/fTShRB5oPh6pdwM1x57Elp2bw+brDfLTyAHEXrZfUVwnxY0S7SrSrFnzNE3l3RScwcdlelu08BYCHm4XuDcPoXLMEYUG+lAgokGXv2/7YROZtOMa3m05w9q/z+7w83IhoW5HHbi+Hp/bUiYgLUpjL21w2zM2fP59ffvnFNs3d3Z1ixYoB8Oabb/L6668zY8YMKlWqxGuvvcbvv//Onj17bDeafeKJJ/j++++ZMWMGRYoU4emnn+bcuXNERkbe0smDuSLMXZaSBNu+hr0/w+HVkJp488sWLgPFqkCxylC0EpgMSDwFSTGQFAuJMXBqB6T/NRxKkQrQLAJq9QAPL7t8nKtJvJTG9DWH+eT3gySmpANQu1QAQ+6oiI+Xe5YLLo6cu8iWv/a0uVmgW92SRLSpROkiN97LlpqeyW+7TzFj7WHWH7SG3Wol/Jlwfy1qlAy4qVozMg2Ltpzgi3VHqFWqME+2Lk9xP32JiojjuXKYi42NZfTo0SxZsoRTp04RGBhI7dq1eeWVV3jvvfeIj49nyZIltvmXLFnCXXfdxYsvvsjYsWNt08eOHcvUqVM5efKk7fyyq1m3bh2NGze2++fKSS4b5hYuXEhUVNQVrxljCA0NJSIigueeew6w7oULDg7mzTffZNCgQcTHx1OsWDG++OILevToAcDJkycJCwtj8eLFdOjQ4aZryVVh7p8y0uDEJji4wvo4sRF8i1gPwwaGQ1C49eegclC0IngVvPF7XjgLf34Mf3wMl85bp/mFQuMnrId2fYPs+IGyOn8xlWm/H2T6msMkp2Vcd95OtUowvG1FKhT3u+X1GGP4dtMJxv64k/MX03B3szDg9nCGt610zStsMzMNP+2IYeKyveyP/XtYlgKebjzctCyPtyhPYEHHBWAREVcOc7fffjtpaWmMHz+ecuXKcerUKX799Vdq1arF8ePHeeaZZ4iLi8PDwzq4xnPPPcecOXMoXbo0q1evtr1PmzZtCAkJ4csvv7SFuV9++SXLUTuAIkWK4Onp6dDP+F+5bJh76623CAgIwNvbm0aNGjFu3DjKlSvHwYMHKV++PJs2baJu3bq2Zbp27UrhwoWZOXMmv/32G23atOHcuXMEBgba5qlduzbdunW74hDuP6WkpJCSkmJ7npCQQFhYWO4Lc/9mTM6d65aSCJEzYN2HkBhtnebuDdW7QYNHIayRw86rO5OUwtQVB/gu6iT+Ph6EBfpSOsjXNlhx1RL+lC16E0H1Bk4npjDm+x38sNX6ecsW8aVV5eKEBfkS9teAyGFBvvx56CzvLN3LjpPWw90BPp70blyaNfvP2s7HK+TtwaPNwxlwezj+BVzrC0NEXJOrhrnz588TGBjIihUraNmy5RWv7927l8qVK2fZm9aoUSMefvhhhg8fTlxcHL6+vqSmplK4cGHef/99BgwYYAtzmzdvpk6dOg7+VDkvJ8Kcw8eZa9SoEZ9//jmVKlXi1KlTvPbaazRt2pQdO3bYzpsLDg7OskxwcDBHjhwBrAMnenl5ZQlyl+e5vPy1jB8//rphL9fKyXDl7QdNh8BtA2HrPPjzE4jZav156zwoXg3qPwJ1eoF3oZxb71UULeTN6M7VGN25ml3XU8zPm8m96tG1zileXLiNw2cvMmPt4WvOX8jbg/7Nw+n/V2B7pr1h+Z5Y3v55LzujE3j/13188vtByhTx/SsQ/h1AUzMyrxifL/5iGs0rFqVHgzAalyuiO2WISL5QqFAhChUqxMKFC2ncuDHe/xpbtVKlSoSGhrJ8+XIaN25MYmIimzZt4ocffmDy5MmsWbOGdu3asX79epKTk2ndurWTPknu5/Aw17FjR9vPNWvWpEmTJpQvX56ZM2fakvm/T4g3xtxwtOubmWfUqFGMGDHC9vzynrl8ycMb6vWFun3g5CbY+D/Y9g3E7oQlz8LvE6DFSKjfz6Hn1dlTu2rBNCoXxI9bozl85kKWwHX+YhoFPN3o1zScQS3KZTmUarFYuKNKMK0qFefnvw7B7otNYndMIrtjbu7cxu+iTvJd1ElKB/nyQP1S3N+gFCUCfMjINJxKuDxAczIx8cmUKVKQJuWLULRQ7hxUWkTkZnh4eDBjxgwee+wxPvroI+rVq0fLli3p2bMntWrVAqBVq1asWLGCUaNGsWrVKipVqkSxYsVo2bIlK1asoF27dqxYsYKwsDDKly+f5f2bNm2K279GfYiPj8+Re526GqffAaJgwYLUrFmTffv20a1bN8C69+2ft7WIjY217a0LCQkhNTWVuLi4LHvnYmNjadq06XXX5e3tfcVfBvmexQIl61sf7V+37p1bPwXiDltD3brJcMeLUOP+PDEGnn8BTx68rfQV0xMupeHl7nbdu1W4uVnoWLMEHaqHcPBMki0IHjv3dyj08nDLsqcuLMgXDzcL3205yfdRJzl67iLvLNvLu7/spWSgDzHxl0jLuPqZDlVC/GhSvghNyxelUbkgHdYVEZdz33330alTJ1atWsW6dev46aefmDBhAp9++in9+vWjdevWREREkJaWxooVK2jVqhUALVu25IMPPgBgxYoV3HHHHVe897x586hatWqWafkxyEEuGGcuJSWF8uXLM3DgQEaPHk1oaCjDhw9n5MiRAKSmplK8ePErLoCYNWsW3bt3ByA6OppSpUrlnQsgnC09FTbNhJUT4IL1nnEE14C2Y6BiW+fW5sKSUzNYvC2aeRuP8eehv4eU8XCzUDLQGv6K+3uzKzqRXdEJWZb18nBj6B0VGNSyvIZYEclnXPWcuWsZMGAAy5Yt48iRIxw4cIAKFSqwZs0ahg0bxrPPPkv37t2Jjo6mTJkyxMTEULJkST766CMefvhhAJ0zdxUO3zP3zDPP0KVLF0qXLk1sbCyvvfYaCQkJPPzww1gsFiIiIhg3bhwVK1akYsWKjBs3Dl9fX3r16gVAQEAA/fv35+mnn6ZIkSIEBQXxzDPPULNmTdq2VdDIER5ecNtj1vPm1k+FNe/Bqe3w5X3QeDC0exXcnb5T1+X4eLlzX/1S3Fe/FEfOXiA6/hJhQb6E+BfA/V/n0Z27kMr6g2dZe+AMa/af5dCZC7y9dC8/bI1mwv21qFWqsHM+hIjIf1StWjUWLlwIQPny5QkLC2PRokVERUXZLpQoUaIEZcuW5Z133uHSpUs6X+4GHP4v8vHjx3nwwQc5c+YMxYoVo3Hjxqxfv54yZcoAMHLkSJKTk3nyySdtgwYvXbrUNsYcwLvvvouHhwfdu3e3DRo8Y8aMfLt71W68CkKLZ6xXua58E/74CNZ/aL1g4v7pUKiYsyt0WWWKFKRMkWtfqRtU0Iu7apbgrpolMMawMOoEr36/k90xiXT7cA39m4czol1lfLy0zYtI7nT27FkeeOABHn30UWrVqoWfnx8bN25kwoQJdO3a1TZf69atmTJlChUqVMhyAeTlQ63lypWjdOkrT485e/bsFRc+Fi5cOE/svbxVTj/M6kw6zHqLdi6ChU9AahL4l4QeX1jPtROHOJuUwpjvd7Joy0kASgf58tb9tWhUroiTKxMRe3LVw6wpKSm88sorLF26lAMHDpCWlkZYWBgPPPAA//d//4ePjw8AM2bM4JFHHuHxxx9n6tSptuVnzZpFnz596N+/P59++qlt+vUGDZ4zZw49e/a07wfLYS45zlxuojCXDaf3wNyH4Ow+6/h0nd6Ben2cXVW+8tvuU7ywYDvR8ZfwcLPw9gO16Va3pLPLEhE7cdUwJzcnJ8KczqSWW1OsMjz2G1TuBBkpsOgpWPwsZF7/Tg6Sc+6oEszS4S3oXKsE6ZmGiHlRfLrqoLPLEhERJ1GYk1tXwB96zLIOWYIF/pwGX/WFtGRnV5Zv+BXw5P2edXm0mfVQw2s/7mLc4l1kZubbHe0iIvmWwpxkj5sbtHgWus+0Hm7d/QN83hUunrvxspIj3NwsjO5clVEdqwAw7feDPP31FtIyMp1cmYiIOJLCnPw31bpC34VQIACO/QGftYe4I86uKt+wWCwMalmedx6ojbubhQWbT9B/5kYSL6U5uzQREXEQhTn578o0hUd/Bv9S1gsjPmsH0VudXVW+cl/9UnzatwE+nu78vvc07d/9nd92n3J2WSIi4gAKc5IzileFAcugeHVIOgXTO8Lepc6uKl9pXaU4cwY2pnSQL9Hxl3h0xkaGzd3M2aQUZ5cmIiJ2pDAnOcc/FB5dAmVvt45FN7s7rHkf8u/oNw5XJ6wwP0e04LHbw3GzwHdRJ2k7cSULN58gH49CJCKSpynMSc4qEAC9v4V6DwMGlo22DjScdsnZleUbPl7uvNCpGguebEaVED/iLqYRMS+Kxz6PJDlVQ8iIiOQ1CnOS8zy8oMt70PEtsLjDljkwszMkxtx4WckxtcMKs+ip5jzdrhJe7m78susUA7/YyKU0BToRkbxEYU7sw2KBRgOhz7dQoDAc3wDTWsOJSGdXlq94ebgxpE1FZj/WCF8vd1btO8OTX24iNV3Dl4iI2NuKFSuwWCycP3/erutRmBP7KtfKeseIopUh8SR82g6WPAfJcc6uLF9pUDaIzx5uiLeHG7/tjmXonM2kazw6EZE8QWFO7K9IeeuVrtW6gcmAPz6CD+pD5EzdBsyBmpQvwid9G+Dl7sZPO2IY8dUWMnTHCBGR6zLGkJ6e7uwyrkthThyjQID1bhF9Flr30l08C98PhU/ugGN/Oru6fKNFpWJM7V0PDzcLi7ac5LlvtuoWYCJiN61atWLIkCFEREQQGBhIcHAw06ZN48KFCzzyyCP4+flRvnx5lixZYltm586d3HXXXRQqVIjg4GD69OnDmTNnbK//9NNPNG/enMKFC1OkSBE6d+7MgQMHbK+npqby1FNPUaJECQoUKEDZsmUZP348AIcPH8ZisRAVFWWb//z581gsFlasWAH8fWj0559/pkGDBnh7e7Nq1SqMMUyYMIFy5crh4+ND7dq1mT9/fpbPu3jxYipVqoSPjw+tW7fm8OHDOd/Uq1CYE8cq3xqeWAMdxoO3P0RHWQcZfrsS/K8jfDcYVk2EnYvg5GaIOwzJ5yFThwRzSpuqwXzwYF3c3SzMjzzO64t3ObskEblFxhgupqY75XGrwxzNnDmTokWL8ueffzJkyBCeeOIJHnjgAZo2bcqmTZvo0KEDffr04eLFi0RHR9OyZUvq1KnDxo0b+emnnzh16hTdu3e3vd+FCxcYMWIEGzZs4Ndff8XNzY177rmHzL/+nXj//fdZtGgRX331FXv27GHWrFmULVv2lns8cuRIxo8fz65du6hVqxYvvvgi06dPZ+rUqezYsYPhw4fTu3dvVq5cCcCxY8e49957ueuuu4iKimLAgAE8//zzt7ze7LCYfDz4VEJCAgEBAcTHx+Pv7+/scvKfpNPw6xiI+hLMDcKaxc0a/nwCoWxzuO0xKFHbMXXmUd9FnWDY3CgA3n+wLnfXDnVuQSJyVZcuXeLQoUOEh4dToEABAC6mplPtpZ+dUs/OVzvg6+VxU/O2atWKjIwMVq1aBUBGRgYBAQHce++9fP755wDExMRQokQJ1q1bx+LFi/njjz/4+ee/P9vx48cJCwtjz549VKpU6Yp1nD59muLFi7Nt2zZq1KjB0KFD2bFjB7/88gsWiyXLvIcPHyY8PJzNmzdTp04dwLpnLjAwkOXLl9OqVStWrFhB69atWbhwIV27dgWsAbJo0aL89ttvNGnSxPZ+AwYM4OLFi8yePZv/+7//Y+HChezYscO23ueff54333yTuLg4ChcufNUeXe33e9nN5pSb+22I2EOhYtB1MnR4Hc4esD7O/eO/CSete+XSk61h79J56yPuEGz+AkrdZg111bqCh7eTP4zr6VqnJHtPJfLh8gM8/81Wqob4UTHYz9lliUgeU6tWLdvP7u7uFClShJo1a9qmBQcHAxAbG0tkZCTLly+nUKFCV7zPgQMHqFSpEgcOHGD06NGsX7+eM2fO2PbIHT16lBo1atCvXz/atWtH5cqVufPOO+ncuTPt27e/5bobNGhg+3nnzp1cunSJdu3aZZknNTWVunXrArBr1y4aN26cJUD+M/jZk9PD3Pjx4/m///s/hg0bxqRJkwDr7uMxY8Ywbdo04uLiaNSoER9++CHVq1e3LZeSksIzzzzDnDlzSE5Opk2bNkyZMoVSpUo56ZNIthUIgJL1rI+rSU+xhrpL5yH+uHVP3s5FcPxP6+OnUVC/HzQban0vuWkj2lUm6th51uw/y+OzIvnuqeYU8nb614KI3ICPpzs7X+3gtHXfCk9PzyzPLRZLlmmXw09mZiaZmZl06dKFN99884r3KVGiBABdunQhLCyMTz75hNDQUDIzM6lRowapqakA1KtXj0OHDrFkyRJ++eUXunfvTtu2bZk/fz5ubtazy/55UDItLe2qdRcsWND28+XA+OOPP1KyZMks83l7e1/xno7m1G/tDRs2MG3atCypHWDChAlMnDiRGTNmUKlSJV577TXatWvHnj178POz7jmIiIjg+++/Z+7cuRQpUoSnn36azp07ExkZibv7rW1okst5eINfsPVRrDJUaAOJp2DT5xA5HRJOwKq3rYMTd54ElW79L7D8yt3Nwns969L5/dUcOH2B5+ZvZXKvulccmhCR3MVisdz0oU5XUq9ePb755hvKli2Lh8eVn+/s2bPs2rWLjz/+mNtvvx2A1atXXzGfv78/PXr0oEePHtx///3ceeednDt3jmLFigEQHR1t26P2z4shrqVatWp4e3tz9OhRWrZsec15Fi5cmGXa+vXrb/jeOcFpF0AkJSXx0EMP8cknnxAYGGibboxh0qRJvPDCC9x7773UqFGDmTNn2o5JA8THx/PZZ5/xzjvv0LZtW+rWrcusWbPYtm0bv/zyi7M+kjiSXzC0fBaGbYXun0NguDXUzX4Avh0IF885u0KXUbSQNx8+VA9Pdws/bovms9WHnF2SiORTgwcP5ty5czz44IP8+eefHDx4kKVLl/Loo4+SkZFBYGAgRYoUYdq0aezfv5/ffvuNESNGZHmPd999l7lz57J792727t3L119/TUhICIULF8bHx4fGjRvzxhtvsHPnTn7//XdefPHFG9bl5+fHM888w/Dhw5k5cyYHDhxg8+bNfPjhh8ycOROAxx9/nAMHDjBixAj27NnD7NmzmTFjhj3adAWnhbnBgwfTqVMn2rZtm2X6oUOHiImJyXJ829vbm5YtW7J27VoAIiMjSUtLyzJPaGgoNWrUsM0j+YS7h/WcuSfWQpOnrBdKbJ0HH94GOxY6uzqXUb9MIC92qgbA+CW7+fOQwrCIOF5oaChr1qwhIyODDh06UKNGDYYNG0ZAQABubm64ubkxd+5cIiMjqVGjBsOHD+ett97K8h6FChXizTffpEGDBjRs2JDDhw+zePFi2yHW//3vf6SlpdGgQQOGDRvGa6+9dlO1jR07lpdeeonx48dTtWpVOnTowPfff094eDgApUuX5ptvvuH777+ndu3afPTRR4wbNy5nG3QNTrmade7cubz++uts2LCBAgUK0KpVK+rUqcOkSZNYu3YtzZo148SJE4SG/n113cCBAzly5Ag///wzs2fP5pFHHiElJSXL+7Zv357w8HA+/vjjq643JSUlyzIJCQmEhYXpata85PhG6/Amp3dbn1fpDHe+AYXDnFuXCzDGMGxuFIu2nKS4nzffD2lOsH+BGy8oInZ1vasdxfXlxNWsDt8zd+zYMYYNG8asWbOuu1H++5wdY8wNz+O50Tzjx48nICDA9ggL0z/weU6pBjDod2jxLLh5wO4frHvpVk2E9FRnV5erWSwWxt9bk4rFCxGbmEL/mRu4kJK7Rz0XEREnhLnIyEhiY2OpX78+Hh4eeHh4sHLlSt5//308PDxslyjHxMRkWS42Ntb2WkhICKmpqcTFxV1znqsZNWoU8fHxtsexY8dy+NNJruDhDXe8aA11pZtA2kXreHYfNYNDv2fvPTPS4UQkrP0A5jwIH7eARUNg85fWoVTyyHCNBb09+PThBgQV9GL7iQSGzd2sW36JiORyDr8Upk2bNmzbti3LtEceeYQqVarw3HPPUa5cOUJCQli2bJntSpPU1FRWrlxpu1S5fv36eHp6smzZMtuo0NHR0Wzfvp0JEyZcc93e3t62S4glHwiuDo8sgS1zYdloOLMXZnaB6vdCuVbgWyTrw80NLpyFC6fh4hnrfxNj4PgGOPoHpF3I+v7RW6xX1AL4FoWwRtb3rXY3+IU4+tPmmDJFCvJJ3wY8+Ml6ftkVy2s/7uTlLtVvvKCIiDiFw8Ocn58fNWrUyDKtYMGCFClSxDY9IiKCcePGUbFiRSpWrMi4cePw9fWlV69eAAQEBNC/f3+efvppihQpQlBQEM888ww1a9a84oIKyecsFqjzIFTuCL+9Bhs/gx3fWh+3qkBhKNPU+ihc2rqn7ugf1tuOXTwDe360PpaMtO4RrH6Pywa7+mUCmdi9Nk/N3sz0NYcpE+RLv2bhzi5LRESuIlcOUjNy5EiSk5N58sknbYMGL1261DbGHFgvPfbw8KB79+62QYNnzJihMebk6nwKQ6e3oW5v2Pg/SDoFF8/+/bgUb53P29+6l65gMShY1PoIqWUNcMWqWvfeXVbNepsX0lOse+mOrIHdP/61J2+t9XE52NXvZw13Hl6O/uTZ1rlWKEfPXWTCT3t49YedlAr0pW21a5/GICIizqF7s+rerAKQkQaZGeCZA1eKxR+Hnd9Zh0Y5/uff0wsFQ8MBUP8R663MXIAxhlHfbmPuhmP4eLrz9eNNqFFSd9kQcSRdzZq3ueTVrCK5krtnzgQ5gIBS0GQwDFgGw3dYL8bwK2HdG7j8dXi3Oix8Eo5tsF5YkYtZLBbGdqvB7RWLkpyWwaMzNnDyfLKzyxIRkX9QmBOxp4BS1mFShm2F+z6DkvUhI8V6f9nP2sIbYTCjs/V8vn2//H24NxfxdHfjw4fqUSnYOmTJozM2kHjp6vcyFBERx1OYE3EEDy+oeT889hv0/wVq3AfeAdZhUw6vgt/fgi/vgzfKwOfdrOfeZWY4u2ob/wKe/K9fQ4oW8mZ3TCJPzd5Mekams8sSEREU5kQcL6wh3P8/eO4wPLEOOr8LtXpa7y+LgYPLYW4veK82rHoHkk47u2IASgX68tnDDSjg6cbKvad5edEO8vEptyKSy5QtW5ZJkybZnlsslitufO8Ir7zyCnXq1HHoOhXmRJzFzQ2Cq0GDR+Hej2FYFAyNgmbDwCcI4o/Br6/Cu9Xg24Fweo+zK6Z2WGHe61kXiwW+/OMon6465OySRESuKjo6mo4dO97UvM4IYDlJYU4kNwkKh3avwohd0G0qhNaDjFTYOg8+bATfPAZn9ju1xA7VQ3jhrqoAjFuyi5+2Rzu1HhHJO1JTc+62iyEhIfnmRgEKcyK5kWcBqNMLBi63nmdXpTNgYNtX8GFDWPAEnDvotPL6Nw+nb5MyGAMR86L4asMxnUMnIldo1aoVTz31FE899RSFCxemSJEivPjii7ZTNMqWLctrr71Gv379CAgI4LHHHgNg7dq1tGjRAh8fH8LCwhg6dCgXLvx9F57Y2Fi6dOmCj48P4eHhfPnll1es+9+HWY8fP07Pnj0JCgqiYMGCNGjQgD/++IMZM2YwZswYtmzZgsViwWKxMGPGDADi4+MZOHAgxYsXx9/fnzvuuIMtW7ZkWc8bb7xBcHAwfn5+9O/fn0uXLuVwF29MYU4ktytZH3p+CQNXQKU7wWTCltnwQQOY1xsiZ0DcEYeWZLFYeKlzNVpXLsaltExGfrOV9u/+zqItJ8nUvVxF7M8YSL3gnMctnis7c+ZMPDw8+OOPP3j//fd59913+fTTT22vv/XWW9SoUYPIyEhGjx7Ntm3b6NChA/feey9bt25l3rx5rF69mqeeesq2TL9+/Th8+DC//fYb8+fPZ8qUKcTGxl6zhqSkJFq2bMnJkydZtGgRW7ZsYeTIkWRmZtKjRw+efvppqlevTnR0NNHR0fTo0QNjDJ06dSImJobFixcTGRlJvXr1aNOmDefOnQPgq6++4uWXX+b1119n48aNlChRgilTptziL/O/06DBGjRYXM3xSFgxDvb/knV6UHkofwdUaAMV24Ob/e+GkpKewcy1h5m64gBxF63DlVQO9mNE+0q0rxaMxWKxew0ied1VB5VNvQDjQp1T0P+dBK+CNzVrq1atiI2NZceOHbbvg+eff55Fixaxc+dOypYtS926dVmwYIFtmb59++Lj48PHH39sm7Z69WpatmzJhQsXOHr0KJUrV2b9+vU0atQIgN27d1O1alXeffddIiIiAOsfnQsWLKBbt25MmzaNZ555hsOHDxMUFHRFna+88goLFy4kKirKNu23337jnnvuITY2Nsvh2goVKjBy5EgGDhxI06ZNqV27NlOnTrW93rhxYy5dupTlva4nJwYNzpW38xKR6yhVH3p/AyejYO9PcGC59RZi5w5YHxs+gdJN4b5PrOPc2ZG3hzsDW5TnwdtKM2PNYaatOsieU4kM+iKSQF9PPNyz7vz3cndjYItyPNy0rF3rEpHco3Hjxln+sGvSpAnvvPMOGRnW4ZcaNGiQZf7IyEj279+f5dCpMYbMzEwOHTrE3r178fDwyLJclSpVKFy48DVriIqKom7dulcNctcSGRlJUlISRYoUyTI9OTmZAwcOALBr1y4ef/zxLK83adKE5cuX3/R6coLCnIirCq1jfbR63jrY8KFV1mFNtsy13hf2o+bQdQpUucvupRTy9mBIm4r0bVKWT1Yd5H9rDtn21P3by4t2kJ5p6N883O51ieRZnr7WPWTOWncOKlgw616+zMxMBg0axNChQ6+Yt3Tp0uzZY72y/1b2/Pv4+NxyXZmZmZQoUYIVK1Zc8dr1gqMzKMyJ5AUFAqBqZ+uj8ZMw/1GIjoK5D8Jtg6xXyObU7cqu4vKXaoCvJ890qMzAluU4EXflbb9+3BrN5OX7GfvDTrw83OjTuIzdahLJ0yyWmz7U6Wzr16+/4nnFihVxd7/6qSD16tVjx44dVKhQ4aqvV61alfT0dDZu3Mhtt90GwJ49ezh//vw1a6hVqxaffvop586du+reOS8vL9uewn/WERMTg4eHB2XLlr1mLevXr6dv375ZPp+j6QIIkbymSHnovwya/HWy8J8fw6dt4cw+u6/aGIMxBv8CnlQt4X/F4+n2lXiiVXkARi/czlcbjtm9JhFxrmPHjjFixAj27NnDnDlz+OCDDxg2bNg153/uuedYt24dgwcPJioqin379rFo0SKGDBkCQOXKlbnzzjt57LHH+OOPP4iMjGTAgAHX3fv24IMPEhISQrdu3VizZg0HDx7km2++Yd26dYD1qtpDhw4RFRXFmTNnSElJoW3btjRp0oRu3brx888/c/jwYdauXcuLL77Ixo0bARg2bBj/+9//+N///sfevXt5+eWX2bFjRw527+YozInkRR5e0OF1eGg++BaFU9vg4xYQNduuq73RYQ+LxcLIDpV5tJn1EOtz327lu6gTdq1JRJyrb9++JCcnc9tttzF48GCGDBnCwIEDrzl/rVq1WLlyJfv27eP222+nbt26jB49mhIlStjmmT59OmFhYbRs2ZJ7773XNnzItXh5ebF06VKKFy/OXXfdRc2aNXnjjTdsewfvu+8+7rzzTlq3bk2xYsWYM2cOFouFxYsX06JFCx599FEqVapEz549OXz4MMHBwQD06NGDl156ieeee4769etz5MgRnnjiiRzq3M3T1ay6mlXyusQY+PYxOPS79XntB+Gut8G7kNNKMsYw+rvtzFp/FHc3C5MfrEvHmiVuvKBIPnS9qx1zu1atWlGnTp0st9mSrHLialbtmRPJ6/xCoM9CaP0iWNxgyxyY1gpitjmtJIvFwqt31+CB+qXIyDQ88eUm2k1cycvfbeen7TGcv5hzo8CLiOR1ugBCJD9wc4eWz0LZZjC/P5zdB5+0sR6KbTjAejK1o0tys/DGfbXwcHdj7oaj7ItNYl9sEjPXHcFigeqh/vRtXJbuDcMcXpuIiCvRYVYdZpX85uI5WPiEdYw6gDoPwd0fOGSQ4WuJu5DKH4fOsvaA9bE/Nsn22hOtyjOyQ2UNQCz5lisfZpUbc8nDrFOnTqVWrVr4+/vj7+9PkyZNWLJkie11YwyvvPIKoaGh+Pj40KpVqyuuDElJSWHIkCEULVqUggULcvfdd3P8+HFHfxQR1+QbBA/OhQ7jweIOUV/Cd4MhM+PGy9pJYEEv7qxRgle71uCXES358//aMLRNRQCmrjjAs/O3kqZ7v4qIXJXDw1ypUqV444032LhxIxs3buSOO+6ga9eutsA2YcIEJk6cyOTJk9mwYQMhISG0a9eOxMRE23tERESwYMEC5s6dy+rVq0lKSqJz585XjBEjItdgsUCTJ+H+z6yBbsscpwe6fyruX4AR7Sox4f5auLtZmB95nIGfb+RiarqzSxMRyXVyxWHWoKAg3nrrLR599FFCQ0OJiIjgueeeA6x74YKDg3nzzTcZNGgQ8fHxFCtWjC+++IIePXoAcPLkScLCwli8eDEdOnS46fXqMKsIsGOB9Tw6k2G90rXrh0495Ppvv+46xeDZm7iUlkmdsMJM79eQwIJezi5LxGEuH4YrU6YMvr45e/cFcb6LFy9y5MgR1703a0ZGBl9//TUXLlygSZMmHDp0iJiYGNq3b2+bx9vbm5YtW7J27VoGDRpEZGQkaWlpWeYJDQ2lRo0arF279rphLiUlhZSUFNvzhIQE+3wwEVdS/R7AYr1rxJY51mm5KNC1qRrMlwMa03/mBqKOnee+j9Yyq38jQgvf+u15RFyRl5cXbm5unDx5kmLFiuHl5aVzSPMAYwypqamcPn0aNzc3vLyy/0eqU8Lctm3baNKkCZcuXaJQoUIsWLCAatWqsXbtWgDbYHyXBQcHc+TIEQBiYmLw8vIiMDDwinliYmKuu97x48czZsyYHPwkInlE9W7W/14OdMZYA5177rjgvX6ZQOY/3oS+n/3JwdMX6PXJer4a1ITi/joZXPI+Nzc3wsPDiY6O5uRJJ92PVezG19eX0qVL4+aW/TPfnPJNXblyZaKiojh//jzffPMNDz/8MCtXrrS9/u+/OIwxN/wr5GbmGTVqFCNGjLA9T0hIICxMwx6IANZAZ7HA14/A1rmQGA33T4eCRZxdGQAVivsx/4mmdP94HYfPXqTXp38wd2BjihbydnZpInbn5eVF6dKlSU9P1/nheYi7uzseHh7/eU+rU8Kcl5eX7Qa6DRo0YMOGDbz33nu28+RiYmKy3LYjNjbWtrcuJCSE1NRU4uLisuydi42NpWnTptddr7e3N97e+uIXuaZqXaHHLPhmABxaCZ+0gp6zIaSmsysDILSwD3Mea0z3j9exPzaJ3p/+wZzHGuscOskXLBYLnp6eeHp6OrsUyWVyxR0gjDGkpKQQHh5OSEgIy5Yts72WmprKypUrbUGtfv36eHp6ZpknOjqa7du33zDMichNqHIXDPgFAsPh/FH4tB1s/8bZVdmEBfky+7HGFPPzZndMIn3/9yfxyWnOLktExGkcHub+7//+j1WrVnH48GG2bdvGCy+8wIoVK3jooYewWCxEREQwbtw4FixYwPbt2+nXrx++vr706tULgICAAPr378/TTz/Nr7/+yubNm+nduzc1a9akbdu2jv44InlTcDUYuBzK3wHpydZz6Za9nGuGLgkvWpDZAxpRpKAX207E02/6nySlaNgSEcmfHB7mTp06RZ8+fahcuTJt2rThjz/+4KeffqJdu3YAjBw5koiICJ588kkaNGjAiRMnWLp0KX5+frb3ePfdd+nWrRvdu3enWbNm+Pr68v333+PunjuuvhPJE3wC4aH50GyY9fmaSTC7O1yKd2pZl1UM9mPWgEYU9vVk89HzPPTJejYdjbvhcukZmQp+IpKn5Ipx5pxF48yJ3KRt8+G7p6x76YpWhl5zIaics6sCYNvxeHp9up7ES9aA1qZKcUa0r0T10IAs8+2PTeKrjcf4dtNxEpLTGdy6Ak+0Ko+XR64420RE5Ao3m1MU5hTmRG7Oyc0wpxckngSfIOjxBZRt7uyqADged5H3ftnHN5uOk/nXN1qnmiV4vGV5dkUnMG/jMSKPXLnXrnKwH2/eX4s6YYUdW7CIyE1QmLsJCnMityghGuY+aA12bp7QeSLU6+vsqmwOnE7ivV/28f3Wk/z7m83dzULrysXo3iCM5LQMxny/k3MXUnGzwCPNwnm6fSV8vXLHuHoiIqAwd1MU5kSyIfUifPek9TZgAE2egnav5po7RgDsjkng3WV7+XnHKcoW8aV7wzDuq1eK4H8MMnzuQipjf9jJgs0nAAgL8uH5O6vSrlqwDr2KSK6gMHcTFOZEsskYWPkmrBhvfV6iDtz9PpSo7dSy/u1SWgbeHm7XHZBz+Z5YXlywnRPnkwEoUtCLe+uVpEfDMCoU97vmciIi9qYwdxMU5kT+ox0L4Pth1itcLW7Q+EloNQq8Czm7sluSlJLORysOMG/jMU4n/n3/5nqlC9OzYWm61g3F2yP37HkUkfxBYe4mKMyJ5IDEU/DT87DjW+vzgDDo9A5U6uDcurIhPSOT5XtOM2/DMZbviSXjr6spQgMKMKRNRe6vXwpPdx2CFRHHUJi7CQpzIjlo3zL4YQTEH7U+r9rFupcuuLpz68qm2IRLfLPpBDPXHiYm4RIAZYr4EtG2InfXLom723+7l6KIyI0ozN0EhTmRHJZ6AZaPg/VTwfx1t4gqneH2p6FkPefWlk2X0jKY/cdRpqzYz5mkVAAqFC/EC3dVpXWV4k6uTkTyMoW5m6AwJ2Inp3bA72/BjoXAX18xFdrC7c9AmSbOrCzbLqamM3PtET5aeYD45DTcLPBF/0Y0q1DU2aWJSB6lMHcTFOZE7Oz0Xlg9EbZ+9feeugptocN4KFbJubVlU8KlNF5YsJ3vt5ykSEEvfhjanBIBPs4uS0TyoJvNKTqTV0Tsp1gluOcjGBIJ9ftZBxre/wtMbQI/v5Br7vN6K/wLePLW/bWoVsKfsxdSGfzlJlLTM51dlojkYwpzImJ/QeHQ5T0Y/AdUuhMy02HdZHi/HkTOhMwMZ1d4Swp4uvNR7/r4F/Bg09HzjFu8y9kliUg+pjAnIo5TpDz0mgcPfQNFKsLFM/D9UJjWEjZ9DilJzq7wppUu4svE7nUAmLH2MN9FnXBuQSKSbynMiYjjVWwLT66DDuPA2x9itsGiIfB2JfhuMBz7kyturpoLta0WzODW5QF4/ptt7D2V6OSKRCQ/0gUQugBCxLkunIHNX8CmL+Dcgb+nF60MLZ6FWg84r7abkJFp6Pu/P1iz/yzlihbku6ea4VfA09lliUgeoAsgRMQ1FCwKzYdbL5J4ZAnU7gWevnBmD3w7AJa+mKvPqXN3s/B+z7qUCCjAwTMXGDl/K/n4b2QRcQKFORHJHSwWKNMU7pkKT++xjkkHsPYDmPtQrj6frkghbz58qB6e7haWbI/hs9WHnF2SiOQjDg9z48ePp2HDhvj5+VG8eHG6devGnj17ssxjjOGVV14hNDQUHx8fWrVqxY4dO7LMk5KSwpAhQyhatCgFCxbk7rvv5vjx4478KCJiLwX8oc1ouO8zcPeGvUvgf3fC+WPOruya6pUO5MVO1QAYv2Q3fx465+SKRCS/cHiYW7lyJYMHD2b9+vUsW7aM9PR02rdvz4ULF2zzTJgwgYkTJzJ58mQ2bNhASEgI7dq1IzHx75OLIyIiWLBgAXPnzmX16tUkJSXRuXNnMjJy7+EYEblFNe+HRxZDweJwaht8cgcc3+jsqq6pb5MydK0TSkamYfDsTcT+dU9XERF7cvoFEKdPn6Z48eKsXLmSFi1aYIwhNDSUiIgInnvuOcC6Fy44OJg333yTQYMGER8fT7Fixfjiiy/o0aMHACdPniQsLIzFixfToUOHm1q3LoAQcRHnj8GcnnBqu3VP3T0fQY17nV3VVV1MTafbh2vYeyqJ28oG8eVjjfB01xktInLrXOYCiPh46wjwQUFBABw6dIiYmBjat29vm8fb25uWLVuydu1aACIjI0lLS8syT2hoKDVq1LDNIyJ5SOEwePQnqNQRMlJg/iOwelKuHL7E18uDqb3rU8jbgz8Pn+Otn/fceCERkf/AqWHOGMOIESNo3rw5NWrUACAmJgaA4ODgLPMGBwfbXouJicHLy4vAwMBrznM1KSkpJCQkZHmIiIvw9oOeX0KjJ6zPf3kZfoiAjHSnlnU15YsV4q37awEw7feDLNkW7eSKRCQvc2qYe+qpp9i6dStz5sy54jWLxZLluTHmimn/dqN5xo8fT0BAgO0RFhaWvcJFxDnc3KHjG3Dnm4AFImfA7O5wKff9YdaxZgkGtigHwLPzt2pAYRGxG6eFuSFDhrBo0SKWL19OqVKlbNNDQkIArtjDFhsba9tbFxISQmpqKnFxcdec52pGjRpFfHy87XHsWO69Mk5ErqPx49BztnU8ugO/wvSOEJ/7bqc1skNlGoUHkZSSziPTNxCbqAsiRCTnOTzMGWN46qmn+Pbbb/ntt98IDw/P8np4eDghISEsW7bMNi01NZWVK1fStGlTAOrXr4+np2eWeaKjo9m+fbttnqvx9vbG398/y0NEXFSVu6Dfj39d6bodPm0DJzY5u6osPNzd+Kh3fcKLFuTE+WQem7mR5FRdcS8iOcvhYW7w4MHMmjWL2bNn4+fnR0xMDDExMSQnJwPWw6sRERGMGzeOBQsWsH37dvr164evry+9evUCICAggP79+/P000/z66+/snnzZnr37k3NmjVp27atoz+SiDhLyXow4BcoVgUSo2H6XbBjgbOryiKwoBfT+zUk0NeTLcfjiZi3mYzM3Hfhhoi4LocPTXKtc9qmT59Ov379AOveuzFjxvDxxx8TFxdHo0aN+PDDD20XSQBcunSJZ599ltmzZ5OcnEybNm2YMmXKLZ0Hp6FJRPKISwkw/1HY/9fe+tYvWO/reoPzbB1pw+FzPPTJH6RmZPLY7eG88NcAwyIi13KzOcXp48w5k8KcSB6SmQFLR8P6D63Pa9wHXT8ETx/n1vUP30WdYNjcKADGdqtBn8ZlnFuQiORqLjPOnIhIjnBzhzvHQZf3wc0Dtn8DMzrBkXWQnuLs6gDoWqckz7SvBMDL321n+e5YJ1ckInmB9sxpz5xI3nNoFXzVB5L/uuLd3RtK1ofSjaF0EyjdCAoEOKU0YwzPzt/K/MjjFPRy5+vHm1ItVN8/InIlHWa9CQpzInnYuYPw61g4vAounM76mocP3P40NBsKHt4OLy01PZN+0/9k7YGzhPgXYOHgZoQEFHB4HSKSuynM3QSFOZF8wBg4ewCOrrM+jqyBuMPW14LKw11vQYU2Di8rPjmN+6auZX9sEtVK+PP1400o6O3h8DpEJPdSmLsJCnMi+ZAxsG0+LH0Bkk5Zp1W9G+4cDwGlrr9sDjt27iL3TFnDmaRU7qhSnE/6NsDdLfdcgSsizqULIERErsZigVoPwFMbofGTYHGHXYtg8m2w9gPrVbEOEhbkyyd9G+Dt4cZvu2N59fsd5OO/r0UkmxTmRCR/KuBv3Rs36HcIawxpF2Dpi/BZe4jd7bAy6pYOZFKPOgDMXHeE6WsOO2zdIpI3KMyJSP4WUgMeWWId0sTbH05shI9vh1XvQEa6Q0roWLMEozpWAWDsjzuZufaw9tCJyE1TmBMRcXOD+g/Dk+uhYnvISIVfX7Xe7/XUDoeUMLBFOfo0LoMx8PKiHUTMi+JiqmPCpIi4NoU5EZHLAkpCr6+g20fWceiio+DjlrDyLbvvpbNYLLzatTovdqqKu5uF76JO0u3DNRw4nWTX9YqI61OYExH5J4sF6jwIg/+EKp0hMw2WvwaftYPTe+y8agsDbi/HnMcaU9zPm72nkug6eQ1LtkXbdb0i4toU5kRErsYvBHrMgns/se6lO7kJPrrdIVe83hYexA9Dm9MoPIiklHSe+HITb/60W+fRichVKcyJiFyLxQK1ulvPpavQDjJSrFe8Tr/LOhCxHRX3K8CXAxoxsEU5AKauOMBHKw/adZ0i4poU5kREbsQ/FB76Gu7+ALz84Nh667l0+3+162o93N34v7uq8nKXagC8+dNuftyqQ64ikpXCnIjIzbBYoF5feHItlG4CqYnw5QMQOdPuq36kWTj9mpYFYPhXUUQeibP7OkXEdSjMiYjcisKloe93UKsHmAz4fij8MgYyM+262tGdq9G2anFS0zN57PONHD170a7rExHXoTAnInKrPLzhno+h5XPW56snwjf9Ie2S3Vbp7mbhvZ51qVHSn3MXUuk340/iL6bZbX0i4joU5kREssNigdb/B92mgpsH7PgWPu8KF87abZUFvT347OGGhAYU4ODpCwyatZHUdPvuERSR3M8pYe7333+nS5cuhIaGYrFYWLhwYZbXjTG88sorhIaG4uPjQ6tWrdixI+so7CkpKQwZMoSiRYtSsGBB7r77bo4fP+7ATyEiAtTpBb2/Be8A64URn7WFM/vttrpg/wJ81q8hhbw9WH/wHD2nrWPfqUS7rU9Ecj+nhLkLFy5Qu3ZtJk+efNXXJ0yYwMSJE5k8eTIbNmwgJCSEdu3akZj49xdWREQECxYsYO7cuaxevZqkpCQ6d+5MRoZ9x38SEblCuZbQf6n1fLpzB623ATu82m6rq1rCnykP1aOglzubjp6n0/uref/XfdpLJ5JPWYyTR6G0WCwsWLCAbt26Ada9cqGhoURERPDcc9bzUVJSUggODubNN99k0KBBxMfHU6xYMb744gt69OgBwMmTJwkLC2Px4sV06NDhptadkJBAQEAA8fHx+Pv72+XziUg+khQLcx6EExvBzRPuft+6585OTpxP5sUF21i+5zQAVUL8ePO+WtQOK2y3dYqI49xsTsl158wdOnSImJgY2rdvb5vm7e1Ny5YtWbt2LQCRkZGkpaVlmSc0NJQaNWrY5rmalJQUEhISsjxERHJMoeLQ7weo1s16G7CFT8CvY+12pWvJwj78r19D3utZh6CCXuyOSeSeKWt47YedXEy1771kRST3yHVhLiYmBoDg4OAs04ODg22vxcTE4OXlRWBg4DXnuZrx48cTEBBge4SFheVw9SKS73n6wP3T4fanrc9XvW3XK10tFgtd65Rk2fAWdKsTSqaBT1cfosOk31mz/4xd1ikiuUuuC3OXWSyWLM+NMVdM+7cbzTNq1Cji4+Ntj2PHjuVIrSIiWbi5QZuXoOsU6+HWHd/CrPvgUrzdVlmkkDeTetZlej/r1a7HziXz0Kd/8OzXWzSEiUgel+vCXEhICMAVe9hiY2Nte+tCQkJITU0lLi7umvNcjbe3N/7+/lkeIiJ2U/ch6P2N9RZgR1Zb7+maYN/bcbWuUpylI1rSt0kZLBb4OvI4bd9dyZJtug2YSF6V68JceHg4ISEhLFu2zDYtNTWVlStX0rRpUwDq16+Pp6dnlnmio6PZvn27bR4RkVyhXEt4ZDEUCoZT2+Gz9nBmn11XWcjbg1e71uDrQU0oX6wgpxNTeOLLTYycv4XMTKde8yYiduCUMJeUlERUVBRRUVGA9aKHqKgojh49isViISIignHjxrFgwQK2b99Ov3798PX1pVcv61VhAQEB9O/fn6effppff/2VzZs307t3b2rWrEnbtm2d8ZFERK6tRC3r0CVB5SH+qDXQHdtg99U2KBvEj0NvZ8gdFXB3s/DVxuP834JtCnQieYxThiZZsWIFrVu3vmL6ww8/zIwZMzDGMGbMGD7++GPi4uJo1KgRH374ITVq1LDNe+nSJZ599llmz55NcnIybdq0YcqUKbd0UYOGJhERh7pwBr58AE5uAg8fuP8zqNLJIav+YetJhs7ZTKaBvk3KMObu6jc8D1lEnOtmc4rTx5lzJoU5EXG4lCT4qi8c+NX6vOFj0H6s9SpYO/t203Ge/noLxsCA5uG80KmqAp1ILuay48yJiORp3oWg1zxo9IT1+YZPYForiN5q91XfW68U4+6pCViHL3ln6V67r1NE7E9hTkTE0dw9oeMb1itdCwXD6d3WW4Ct/cBuAwxf9uBtpRlzd3UAJi/fzwe/2vdiDBGxP4U5ERFnqdAWnlgLlTtBRiosfRG+6AaJ1x78PCc83LQsL9xVFYB3lu3l5e+2676uIi5MYU5ExJkKFoWeX0KX98DTFw6thI9uh0Or7Lrax1qU47k7qwAwc90Rek5bR3R8sl3XKSL2oTAnIuJsFgvU7weDfofi1eFCLHx+N6x+166HXZ9oVZ5P+zbAr4AHm46ep/P7q3ULMBEXpDAnIpJbFK0IA36B2r3AZMIvr8DcXpAcd8NFs6tttWB+GNKcaiX8OXshlT6f/cGHy/drLDoRF6IwJyKSm3j5Qrcp1sOu7t6wdwl83BJORtltlWWKFOTbJ5vyQP1SZBp46+c9PPnlJp1HJ+IiFOZERHKby4dd+y+FwmXg/BH4XwfY+rXdVlnA050J99fijXtr4uXhxk87YhgyZxNpGQp0IrmdwpyISG4VWgcGrYSKHSD9Enw7AH4ZY7fz6CwWCz1vK82nfRvg5e7GzztOMeKrLWTokKtIrqYwJyKSm/kEwoNzoNkw6/PVE63n0aUk2m2VLSoVY2rveni4Wfh+y0lGzt+qc+hEcjGFORGR3M7NHdq9CvdM+/s8uk/bwblDdltlm6rBfPBgXdzdLHyz6TgvfredfHz3R5FcTWFORMRV1O4BjyyGQiFwehd80hoO/Ga31XWsWYKJ3WtjscDsP47y6g87FehEciGFORERV1KqAQxcDqF1rUOWfHGvdQiTjDS7rK5rnZJMuK8WANPXHGbInM2cTUqxy7pEJHsU5kREXI1/KDyyBOo/Ahjr4MLTO0LcEbus7oEGYYy7pyZuFvhhazRtJ65kwebj2ksnkksozImIuCJPH+gyCR6YAd4BcHyD9TZgOxbYZXW9GpVm4eBmVAnxI+5iGsPnbaHf9A0cj7tol/WJyM2zmHz8p1VCQgIBAQHEx8fj7+/v7HJERLIn7gh8098a6MC6x67jm+DhneOrSsvIZNrvB3nv132kpmfi6+XOM+0r07txGbw8tH9AJCfdbE5RmFOYE5G8ICMNlo+zHnLFQJnm0HOWdWgTO9gfm8Sob7ey4bD1VmNhQT4Ma1OJbnVC8XBXqBPJCTebU1z+/7gpU6YQHh5OgQIFqF+/PqtWrXJ2SSIijufuCW1fhofmg5cfHFkNn7WHuMN2WV2F4oWYN7AJr99Tg6KFvDl2Lplnvt5C+0m/8/2WkxqXTsSBXHrP3Lx58+jTpw9TpkyhWbNmfPzxx3z66afs3LmT0qVL33B57ZkTkTwpZjvM7g4JJ6BgMXhwrvUqWDtJTs3g83WHmbryAOcvWq+qrRLix0ONy9C0fBHKFS2IxWKx2/pF8qp8cZi1UaNG1KtXj6lTp9qmVa1alW7dujF+/PgbLq8wJyJ5VkK0NdDFbAWPAnDvJ1DtbruuMvFSGv9bfZhPVx0kMSXdNj3EvwBNyxehSfkiNAovQslAH9zdFO5yq4xMw5mkFAJ8PCng6e7scvK1PB/mUlNT8fX15euvv+aee+6xTR82bBhRUVGsXLnyhu+hMCcieVpKEsx/BPYtBSxQ/g7wDYICAeDt/9d/C4ElZ8+4uZiawZ+Hz7H/VBJHzl284t6u7m4Q4ONFUEFPAn29CCzohY9Cg3MYSE7LIO5iKnEXUjl3MY3zF1O5/Cvz8/YgsKAXQb6eBBb0wq+AJ9rJ+rfi1VsSXq2h3d7/ZnOKh90qsLMzZ86QkZFBcHBwlunBwcHExMRcdZmUlBRSUv4e7DIhIcGuNYqIOJV3Ieg5B5aMhI2fwYFfHbJaX6DVXw/c/3r8W+pfjziHlCS34p/JIBNI/OshV/gjI82uYe5muWyYu+zf52EYY655bsb48eMZM2aMI8oSEckd3D2g0ztQ8344dxAuxWd9pDj2X2kDXErLJDktnYupGSSnZpCclnHF3jtxHHc3Cz6e7vh6eeDj5YaPpwcFPN1IyzBcTE0nOe3v31Nqeqazy81VChQv7+wSABcOc0WLFsXd3f2KvXCxsbFX7K27bNSoUYwYMcL2PCEhgbCwMLvWKSLidBYLlGlqfTi7FMDnr0eQk2uR6/P661HYyXXIjbns0CReXl7Ur1+fZcuWZZm+bNkymja9+heWt7c3/v7+WR4iIiIirsxl98wBjBgxgj59+tCgQQOaNGnCtGnTOHr0KI8//rizSxMRERFxCJcOcz169ODs2bO8+uqrREdHU6NGDRYvXkyZMmWcXZqIiIiIQ7js0CQ5QUOTiIiISG6V54cmyQmXc6yGKBEREZHc5nI+udF+t3wd5hITrZfk64pWERERya0SExMJCAi45uv5+jBrZmYmJ0+exM/PL8fvG3h52JNjx47pEO4tUN8cQ33OHvXNedT77FPv7M9ePTbGkJiYSGhoKG5u1x6AJF/vmXNzc6NUqVJ2XYeGQMke9c0x1OfsUd+cR73PPvXO/uzR4+vtkbvMZceZExERERGFORERERGXpjBnJ97e3rz88st4e3s7uxSXor45hvqcPeqb86j32afe2Z+ze5yvL4AQERERcXXaMyciIiLiwhTmRERERFyYwpyIiIiIC1OYExEREXFhCnMiIiIiLkxhThxq06ZNtnviiuQ22j7F1WibFVCYuyWnTp3ixx9/RKO53LqTJ0/Svn17WrduTVRUlLPLyZO0fWaftk/n0XabPdpmHcNVtk+FuZs0efJkQkND6dKlCzt27HB2OS5l5MiRlClTBl9fX3bt2sXtt9/u7JLyHG2f2aft03m03WaPtlnHcKnt08h1ZWZmmh9//NG0adPGvP3226ZevXrm/vvvNxkZGc4uLddLTU01Tz31lLFYLGbu3Lm26adOnXJiVXmLts/s0/bpPNpus0fbrGO44vbp4ewwmdtZLBaCg4Pp06cP9913Hw0bNqRVq1b8/PPPdOzY0dnl5VrGGDw9Pbn99tvZtm0bZ86cYffu3YwaNYozZ87g5ubGQw89RL9+/fDy8nJ2uS5L22f2aPt0Lm23t07brOO44vap23n9S0JCAjt37iQsLIySJUtedZ4ePXqwb98+Vq5ciZ+fn4MrzL1SU1PJyMjAx8eHjIwM3N3dSU9PZ/jw4cyfP5+0tDR69+5NuXLl2LFjBzNnzmTcuHE88cQT+Pj4OLt8l6DtM/u0fTqPttvs0TbrGHli+3TqfsFcZty4ccbf39/UqFHD+Pv7m0mTJpnjx48bY4xJT0+37WI9cOCA8fHxMe+//74zy81V3njjDVOpUiXz008/2aalp6cbY4xZt26defjhh82iRYuyLDN06FBTu3Zts23bNofW6qq0fWaftk/n0XabPdpmHSOvbJ8Kc39ZvHixqVq1qlmwYIE5ePCgef3110316tXNo48+apsnMzPT9vOLL75ogoODzbFjx4wxxly4cMEkJSU5vG5nO3v2rHn88cdNrVq1jL+/v7n33nvN6dOnjTFZ+7V161Zz6dIlY4yx/c8RExNjLBaL+eOPPxxfuIvR9pk92j6dS9vtrdM26zh5aftUmPvL0KFDTd26dbNM++CDD0zlypXNtGnTjDF//1VkjDFJSUmmTJkyZujQoebzzz83zZs3N1999ZVDa84NDh48aEaOHGl+/PFHs2rVKmOxWMycOXNsXy7//B/hssvT5syZY4oXL262bNni0JpdkbbP7NH26Vzabm+dtlnHyUvbp8Kcsf5V88QTT5iePXva/tIxxpiTJ0+aQYMGmdq1a5vExETbvJe9/PLLxmKxGC8vLzNq1CiH150bpKenmyNHjtied+/e3dSqVcscOnToqvNf/tLZtWuXad++vXnsscccUaZL0/aZfdo+nUfbbfZom3WMvLZ95vswd/l/hPHjx5uwsLAr/odZtGiRadCggS2lG2NN54MHDzYWi8X079/fxMXFObDi3OlyH8+ePWs8PT3N+PHjs/wPYox1l/SYMWNMv379jK+vr3nooYdMQkKCM8p1Gdo+c4a2T8fSdvvfaZu1n7y4feb7QYMzMzMBiIiIID4+ni+//DLL661atcLNzY2zZ8/app05cwY/Pz9WrVrFp59+SuHChR1Zcq5h/nEhtMViIT09naCgIF544QUmTpzIrl27sszr6+tLYGAgSUlJrFy5klmzZuXOq4Kc4NKlS1edru3zxq7VO22f9hcbG3vVkfG13V7ftfqmbTZnbdy40fb98M/e5snt02kx0kFOnjxp7r//fjNv3jxjTNbj32lpaVnmffvtt42fn5/ZsGFDlul16tQxTz75pP2LzWVutneXp//z9ZIlS5qBAweac+fOmZ9//tnMmDHDGGNy9aCLznDw4EFTq1YtM3r06Cte0/Z5fTfbO22fOe/gwYOmS5cupl+/fmb79u1ZXtN2e2032zdts//NgQMHzN13320sFouZOXNmltfy6vaZ58Pc2LFjjcViMY0bNzYXLlwwxmTd+DMzM83IkSPNF198YYwxpn79+qZNmzbmxx9/NMYYExkZaWrXrm2WLl3q+OKd7GZ7N2vWLNv0y18+3377rXF3dzc1a9Y0FovFfPjhh47/ALlYZmamGTRokPHw8DD333+/7Wq1q82n7TOrW+2dts+ccfnQ1MyZM01gYKDp3r272bBhg4mJicny+uWftd1aZadv2mazJzMz0zzxxBPGzc3N3H333aZw4cLmm2++Mcbk/X/383yY69ixo+nZs6dp0aKFeemll4wxf/9SZ8yYYYKCgkydOnXMpk2bjDHG7Ny509xzzz3Gy8vLtG/f3vj6+pqePXvawkx+crO927p1a5bljh8/bkaOHGksFovp2bOn7TJusdq3b58JDAw01atXt213V6Pt80rZ6Z22z5yTkZFh7rzzTjNx4kTbtH9ve9OnT9d2+y+32jdts7duwYIFxsfHx9x2221m/fr1xhhjmjRpYgYMGGCM+Ts059Xv1TwT5v59ufblXamPPvqoWbBggRk1apSpWrWq2blzpzHGmISEBPPaa6+ZKVOm2P7yufwe8fHxZunSpWby5Mlm9erVDvwUzpETvbssJSXFREREmKCgILN8+XKH1O8K/tnjw4cPm+rVq5tBgwYZY4xZs2aNGTFihHn99dfNkiVLbFdQvfzyy2bq1KnaPnOgd5dp+7w1//5u+PXXX03FihVNWlqaWbNmjbn77rtNhw4dzODBg826deuMMca88sor+f57NSf6dpm22Wv7Z5+nTZtm2wtnjDEXL140vXv3Nvfee68tlCUlJZmxY8fmye/VPHE7r+TkZNzc3PD29gasJzpaLBYAatWqxdy5c7l48SIjR46kZs2avPXWW+zbt48qVarg7u7uzNKdzh69O3XqFMHBwQ77DLndv3ucmZnJwoULuf/++2nfvj27d++mQYMG7N+/n9OnT9OuXTtmzJjh3KJzCXv0Ttvnzfl37wH++OMPOnXqxKeffsrrr79O27ZtKVCgAL/99hvbtm1j586dhISEOLFq57NH37TNXulq3w1ubtZrOi/f+mz48OH89ttvbNmyxfZv2z/ny1OcGiVzwPPPP2/q1atn2rZta9577z0THx9vjLHu1j5+/HiW870mTpxoihYtaiwWi3nvvfdMSkqKM0t3OvXO/q7V43Pnzpm+ffuaZs2amS1bttiGHJg2bZqpXLmymTJlijEmf5/crN45z7V6v27dOtOmTRvTsGFDM2DAAFuPExISTI0aNcwjjzxijDFX7F3KL9Q3x/h3ny8Px3K5r5f/++OPPxp/f/8s4/blVS4b5lJSUsz9999vqlWrZubOnWv69u1rqlWrZjp16mSbJz4+3tx+++3m4sWL5ttvvzVBQUEmICDA1K5d2zbP1UbTzuvUO/u7Vo/vuusu2zy7du0yGzZsMJmZmbYvn7Nnz5rOnTubgQMH5tsvdvXOea7V+44dOxpjrOd53XPPPcZisVwxQv5nn31mSpYsaTvUnZ+ob45xM/92/dN3331nwsPDXfKw6a1y2TC3c+dOU7FixSxXm6xevdr4+PiYCRMmGGOs5ymUKFHC1KhRwxQuXNi8/fbb5uOPPzZ16tSxXQmUH/96V+/s72Z6/G+Xw3GFChXM448/7pA6cyP1znmu1/vx48cbY6z3syxSpIhp165dlmWff/5506pVK3PhwoV894ee+uYYN/puuNy/f/6B5+XlZX744Ycs0/Milw1zkZGRxmKxmLNnzxpjso7oXLhwYXPw4EGTlpZmqlWrZgYOHGgb4fnkyZOme/fupkWLFleMpp1fqHf2d70eBwYGmr179151uSVLlpiGDRuaNWvWOKzW3Ea9c57r9T4gIMAcOHDAGGM9Wb9IkSJm9OjRZu/evWb37t2mZcuW5tVXX3Va7c6kvjnGrX43nD9/3rRo0cI8/fTTDq/V0Vw2zG3evNlUr17dfPDBB8aYv3+pqamppmzZsiYiIsIYY8ypU6eu+Gtnx44d+TqMqHf2d70eh4eH275cMjIyzLZt28xvv/1mBg0aZAICAszzzz+frw8TqnfOc6PeX/5uiImJMdOmTTOFCxc2NWrUMH5+fuaRRx7Jt98N6ptj3Ox3w+URGdLT003FihXN448/blJTU51TtIO4bJg7d+6c6datm+nRo4c5efKkMebvX+A777xjSpQoccUuVe3CtlLv7O9GPQ4NDbX1eObMmaZ169amdevWJioqymk15xbqnfPc6nfDiRMnzPr1682+ffucUm9uob45xq18N1z+o+7zzz83e/bscU7BDpQrr8+NjY3l9OnTpKamAtbLjC9LT08HIDAwkC5durB7926++uorADw8PAAICAggKCiIY8eOZXnfy0Nu5GXqnf3lRI8DAwM5cuQIAPfddx+ffPIJv/32G7Vr13bkR3E49c55cvK7wfw1olVoaCiNGjWiQoUKjvwoDqW+OUZOfTdc/rfr8tBZffr0oVKlSg77HM6Sq8JcWloajz/+OC1atKBLly7cfffdpKSk4O7uTlpaGmD9xV26dIm5c+fy6KOPUqdOHebNm8fy5ctt73P8+HGKFStGmTJlnPVRHE69s7+c7nF4eDgABQsWpHz58k75TI6i3jmPPb4b8sMfd+qbY+jfrhzi7F2Dl3399demfPnypmXLlua3334z06ZNM+XKlbviRrfvvfeeCQoKMl27djXGGLNlyxbz0EMPGS8vL/PEE0+YgQMHGj8/PzN16lRjTP44PKje2Z96nH3qnfOo99mjvjmG+pxzck2YGzx4sBk9erTt+Lcxxjz88MNmxIgRtucffPCBKVu2rPnyyy+vuGnuuHHjzGOPPWbuuuuufHc1m3pnf+px9ql3zqPeZ4/65hjqc85x+u28Lt9aIyYmhrS0NMLCwgA4cuQI9957L7169aJJkyY0bdqU9PR0UlJSKFiwoG1584/bT+U36p39qcfZp945j3qfPeqbY6jPduCMBPnxxx+badOmmRUrVtim/XO36Pvvv28sFotp3ry5admypQkMDDQvvfSSSU5Odka5uYp6Z3/qcfapd86j3meP+uYY6rN9OTTMzZ492xQvXtw0adLE1KlTxxQrVsy89tprxpis96SbMWOG+f33322/6C+//NL4+PiYw4cPO7LcXEW9sz/1OPvUO+dR77NHfXMM9dkxHBbmvvzyS1O7dm3z0UcfGWOs4+x88MEHpmDBgrab5F5rsM9du3YZd3f3LLfwyE/UO/tTj7NPvXMe9T571DfHUJ8dx+5Dk5i/TslLS0ujUaNG9O3bF7COs1O3bl1KlizJrl27gL/Hhfm3hQsX0qZNG5o3b27vcnMV9c7+1OPsU++cR73PHvXNMdRnJ7BXSoyMjDRxcXG25+fPn78igUdFRZmQkBBz7ty5K5Y/cuSI2b9/vxkwYIAJDQ01M2bMMMbkj0uO1Tv7U4+zT71zHvU+e9Q3x1CfnSfHw9z8+fNNqVKlTPny5U3p0qXN6NGjTUxMjO31f15aPHHiRNOsWTNjjMly37S9e/eap59+2pQqVcq0bt06X9yKwxj1zhHU4+xT75xHvc8e9c0x1Gfny9Ewt2HDBlOlShUzadIks2XLFjNlyhRTrFgx88QTT5izZ88aY6y/1Mtjytxzzz1m8ODBV7zPxYsXzYoVK/LVuDHqnf2px9mn3jmPep896ptjqM+5Q46Eucu7QKdOnWpKlSpl4uPjba9NnjzZNG7c2IwdO9Y2LSMjw2RmZpry5cubH374wRhjzJ49e0zPnj3N0aNHc6Ikl6He2Z96nH3qnfOo99mjvjmG+py75MgFEJcH7zt06BCVKlWy3fgWoF+/ftSvX58lS5awY8cOANzc3NiwYQO+vr7Uq1ePiIgIatWqxdmzZylevHhOlOQy1Dv7U4+zT71zHvU+e9Q3x1Cfc5dshblly5YxdOhQ3nvvPf7880/b9GbNmrF27VpiYmIAyMjIoGDBgnTt2hWLxcLSpUtt8y5evJjt27dTuXJlli1bxpo1a1i6dCne3t7/8SPlbuqd/anH2afeOY96nz3qm2Ooz7nbLYW56OhounTpQu/evTl37hyfffYZ7du3t/1i27dvT9myZXnzzTeBv5N7u3btcHNzY//+/bb38vT0pGjRosyYMYMdO3ZQv379nPpMuZJ6Z3/qcfapd86j3meP+uYY6rOLuNnjsRcuXDAPP/yw6dGjhzl48KBtesOGDU2/fv2MMdbB/z7//HPj5uZ2xUmMDz30kGnVqpXteWxs7H89ROwy1Dv7U4+zT71zHvU+e9Q3x1CfXcdN75nz9fXF29ubfv36ER4eTnp6OgCdO3fOMvhf9+7d6dq1KwMGDGDlypUYY4iJiWHfvn307t3b9n7FihXL4Viae6l39qceZ5965zzqffaob46hPrsOizF/DdV8E9LS0vD09ASsIzxbLBb69OmDj48P06ZNs027dOkSHTt2ZOfOndSpU4ft27dTunRpvvrqK8LCwuz2YXIz9c7+1OPsU++cR73PHvXNMdRn13BLYe5qWrRowaOPPkq/fv0wxpCZmYm7uzunTp1i69atbNiwgbJly9KrV6+cqjnPUO/sTz3OPvXOedT77FHfHEN9zn3+U5g7ePAgTZs25ccff7SdyJiamoqXl1eOFZhXqXf2px5nn3rnPOp99qhvjqE+507ZGprkcv5bvXo1hQoVsv1Cx4wZw7Bhw4iNjc25CvMY9c7+1OPsU++cR73PHvXNMdTn3M3jxrNc6fKlx3/++Sf33Xcfy5YtY+DAgVy8eJEvvvhCAwBeh3pnf+px9ql3zqPeZ4/65hjqcy6X3ctgk5OTTYUKFYzFYjHe3t7mjTfeyO5b5Tvqnf2px9mn3jmPep896ptjqM+51386Z65du3ZUrFiRiRMnUqBAgZzMmHmeemd/6nH2qXfOo95nj/rmGOpz7vSfwlxGRgbu7u45WU++od7Zn3qcfeqd86j32aO+OYb6nDv956FJRERERMR5snU1q4iIiIjkDgpzIiIiIi5MYU5ERETEhSnMiYiIiLgwhTkRERERF6YwJyIiIuLCFOZEJN9o1aoVERER+W7dIpK3KcyJiFzFihUrsFgsnD9/PkeW+/bbbxk7dmzOFSgi8hcPZxcgIpIfBAUFObsEEcmjtGdORPKkCxcu0LdvXwoVKkSJEiV45513srw+a9YsGjRogJ+fHyEhIfTq1YvY2FgADh8+TOvWrQEIDAzEYrHQr18/AIwxTJgwgXLlyuHj40Pt2rWZP3/+DZf792HWsmXL8tprr9lqLFOmDN999x2nT5+ma9euFCpUiJo1a7Jx48Ysda9du5YWLVrg4+NDWFgYQ4cO5cKFCzndPhFxIQpzIpInPfvssyxfvpwFCxawdOlSVqxYQWRkpO311NRUxo4dy5YtW1i4cCGHDh2yBa+wsDC++eYbAPbs2UN0dDTvvfceAC+++CLTp09n6tSp7Nixg+HDh9O7d29Wrlx53eWu5t1336VZs2Zs3ryZTp060adPH/r27Uvv3r3ZtGkTFSpUoG/fvly+6+K2bdvo0KED9957L1u3bmXevHmsXr2ap556yh4tFBFXYURE8pjExETj5eVl5s6da5t29uxZ4+PjY4YNG3bVZf78808DmMTERGOMMcuXLzeAiYuLs82TlJRkChQoYNauXZtl2f79+5sHH3zwmssZY0zLli2zrLtMmTKmd+/etufR0dEGMKNHj7ZNW7dunQFMdHS0McaYPn36mIEDB2Z531WrVhk3NzeTnJx8/aaISJ6lc+ZEJM85cOAAqampNGnSxDYtKCiIypUr255v3ryZV155haioKM6dO0dmZiYAR48epVq1ald93507d3Lp0iXatWuXZXpqaip169a95Tpr1apl+zk4OBiAmjVrXjEtNjaWkJAQIiMj2b9/P19++aVtHmMMmZmZHDp0iKpVq95yDSLi+hTmRCTPMX8dlryWCxcu0L59e9q3b8+sWbMoVqwYR48epUOHDqSmpl5zucuB78cff6RkyZJZXvP29r7lOj09PW0/WyyWa067vN7MzEwGDRrE0KFDr3iv0qVL3/L6RSRvUJgTkTynQoUKeHp6sn79elvIiYuLY+/evbRs2ZLdu3dz5swZ3njjDcLCwgCuuNDAy8sLgIyMDNu0atWq4e3tzdGjR2nZsuVV13215XJKvXr12LFjBxUqVMjx9xYR16ULIEQkzylUqBD9+/fn2Wef5ddff2X79u3069cPNzfrV17p0qXx8vLigw8+4ODBgyxatOiKMeDKlCmDxWLhhx9+4PTp0yQlJeHn58czzzzD8OHDmTlzJgcOHGDz5s18+OGHzJw585rL5ZTnnnuOdevWMXjwYKKioti3bx+LFi1iyJAhObYOEXE9CnMikie99dZbtGjRgrvvvpu2bdvSvHlz6tevD0CxYsWYMWMGX3/9NdWqVeONN97g7bffzrJ8yZIlGTNmDM8//zzBwcG2K0bHjh3LSy+9xPjx46latSodOnTg+++/Jzw8/LrL5YRatWqxcuVK9u3bx+23307dunUZPXo0JUqUyLF1iIjrsZgbnVwiIiIiIrmW9syJiIiIuDCFOREREREXpjAnIiIi4sIU5kRERERcmMKciIiIiAtTmBMRERFxYQpzIiIiIi5MYU5ERETEhSnMiYiIiLgwhTkRERERF6YwJyIiIuLCFOZEREREXNj/AwP//5azBPzjAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(2,1)\n", "df_butte_2023['SWE'].plot(label='measured', ax = axes[1])\n", "df_butte_2023['SWE predicted'].plot(label='predicted', ax = axes[1])\n", "df_butte_2023['PRECIPITATION'].plot.bar(ax=axes[0])\n", "axes[0].set_xticks([])\n", "plt.legend(title='SWE')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So it looks like our predicted SWE results in a melt out about 5 days earlier than actual SWE. Not too bad, huh? Also, it looks like a lot of the divergence between measured and predicted ablation coincides with precipitation." ] } ], "metadata": { "kernelspec": { "display_name": "snow-hydrology", "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.10.13" } }, "nbformat": 4, "nbformat_minor": 2 }