{
"cells": [
{
"cell_type": "markdown",
"id": "27b5cc66-b11d-42dd-b46d-ba8169bd0ada",
"metadata": {},
"source": [
"## Tutorial: Combining adjoint sensitivities with emissions data to calculate air quality impacts\n",
"### Example: LA County\n",
"\n",
"In this tutorial we will demonstrate how to combine results from the GEOS-Chem adjoint, known as sensitivities, with emisisons data from the US EPA's National Emission Inventory (NEI) to determine how different emission scenarios would impact air quality in Los Angles, CA. \n",
"\n",
"First, we will confirm that you have the correct files needed to perform this tutorial. You should have four netcdf files in this folder. Three sensitivity files of the form:\n",
"\n",
">la_sens_no.nc \n",
"la_sens_o3.nc \n",
"la_sens_pm.nc \n",
">\n",
"\n",
"And eleven emissions files each of the form:\n",
"\n",
"> nei_emis_onroad.nc\n",
"\n",
"For each of the emissions files the component of the filename following \"nei_emis_\" will indicate the sector. For both the sensitivities and files there will be date and time information in the filenames to indicate when the file was created. These files should be stored in the same directory as this tutorial.\n",
"\n",
"We will present a quick example of how these files can be used. We will check to make sure these files are in the current directory but first we need to load in the python libraries we'll be using."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "3f91f282-a5b0-4ed6-a405-53740b4947a3",
"metadata": {},
"outputs": [],
"source": [
"\n",
"# for this tutorial we will use netCDF4, xarray, numpy, os, and matplotlib.pyplot\n",
"import netCDF4 as nc\n",
"import xarray as xr\n",
"import numpy as np\n",
"import os\n",
"import matplotlib.pyplot as plt\n",
"import datetime\n"
]
},
{
"cell_type": "markdown",
"id": "6b2e43d3-d1d8-4f70-b29b-c5177901d592",
"metadata": {},
"source": [
"Let's take a look at the files in the current directory. We'll also store specific filenames used in this tutorial in variables."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "fdef4c72-3089-4734-bdbb-c04c7bf99063",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"la_sens_no_19112021_080520.nc\n",
"la_sens_o3_19112021_080452.nc\n",
"la_sens_pm_19112021_080507.nc\n",
"nei_emis_ag_19112021_082921.nc\n",
"nei_emis_c3marine_19112021_082923.nc\n",
"nei_emis_egu_19112021_082925.nc\n",
"nei_emis_egupk_19112021_082933.nc\n",
"nei_emis_nonroad_19112021_082936.nc\n",
"nei_emis_oilgas_19112021_082938.nc\n",
"nei_emis_onroad_19112021_082940.nc\n",
"nei_emis_onroad_catx_19112021_082946.nc\n",
"nei_emis_othpt_19112021_082948.nc\n",
"nei_emis_ptnonipm_19112021_082951.nc\n",
"nei_emis_sf_19112021_082953.nc\n"
]
}
],
"source": [
"# get a list of the files in the current directory\n",
"files = os.listdir('.')\n",
"\n",
"# loop through these files\n",
"for name in files: \n",
" if name.endswith('.nc'): \n",
" # print all files ending with the netcdf extension (.nc)\n",
" print(name)\n",
" # save the filenames into variables for the pm2.5 sensitivity file and the onroad_catx emission file\n",
" if \"_pm\" in name: pm_sens_filename = name\n",
" if \"_onroad_catx\" in name: on_road_emis_filename = name \n",
" "
]
},
{
"cell_type": "markdown",
"id": "885481bf",
"metadata": {},
"source": [
"You should see 14 files listed above: 11 emissions files and 3 sensitivity files. For the purpose fo this tutorial we will focus on just two of these files, the pm sensitivity file and the onroad catx emissions file. Let's read in data from the latter first and print out info from the file:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6cc914bf",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xarray.Dataset {\n",
"dimensions:\n",
"\tlat = 89 ;\n",
"\tlon = 91 ;\n",
"\ttime = 396 ;\n",
"\n",
"variables:\n",
"\tfloat64 lon(lon) ;\n",
"\t\tlon:units = degrees_east ;\n",
"\t\tlon:long_name = longitude ;\n",
"\t\tlon:comment = centre of grid cell ;\n",
"\tfloat64 lat(lat) ;\n",
"\t\tlat:units = degrees_north ;\n",
"\t\tlat:long_name = latitude ;\n",
"\t\tlat:comment = centre of grid cell ;\n",
"\tdatetime64[ns] time(time) ;\n",
"\t\ttime:long_name = reference time of sst field ;\n",
"\t\ttime:comment = ;\n",
"\t\ttime:axis = T ;\n",
"\tfloat64 BC(time, lat, lon) ;\n",
"\t\tBC:long_name = emissions of BC ;\n",
"\t\tBC:units = (kg / m2 / s) ;\n",
"\tfloat64 OC(time, lat, lon) ;\n",
"\t\tOC:long_name = emissions of OC ;\n",
"\t\tOC:units = (kg / m2 / s) ;\n",
"\tfloat64 NH3(time, lat, lon) ;\n",
"\t\tNH3:long_name = emissions of NH3 ;\n",
"\t\tNH3:units = (kg / m2 / s) ;\n",
"\tfloat64 NOx(time, lat, lon) ;\n",
"\t\tNOx:long_name = emissions of NOx ;\n",
"\t\tNOx:units = (kg / m2 / s) ;\n",
"\tfloat64 SO2(time, lat, lon) ;\n",
"\t\tSO2:long_name = emissions of SO2 ;\n",
"\t\tSO2:units = (kg / m2 / s) ;\n",
"\tfloat64 SOAP(time, lat, lon) ;\n",
"\t\tSOAP:long_name = emissions of SOAP ;\n",
"\t\tSOAP:units = (kg / m2 / s) ;\n",
"\tfloat64 SO22(time, lat, lon) ;\n",
"\t\tSO22:long_name = emissions of SO22 ;\n",
"\t\tSO22:units = (kg / m2 / s) ;\n",
"\tfloat64 ALK4(time, lat, lon) ;\n",
"\t\tALK4:long_name = emissions of ALK4 ;\n",
"\t\tALK4:units = (kg / m2 / s) ;\n",
"\tfloat64 ACET(time, lat, lon) ;\n",
"\t\tACET:long_name = emissions of ACET ;\n",
"\t\tACET:units = (kg / m2 / s) ;\n",
"\tfloat64 MEK(time, lat, lon) ;\n",
"\t\tMEK:long_name = emissions of MEK ;\n",
"\t\tMEK:units = (kg / m2 / s) ;\n",
"\tfloat64 ALD2(time, lat, lon) ;\n",
"\t\tALD2:long_name = emissions of ALD2 ;\n",
"\t\tALD2:units = (kg / m2 / s) ;\n",
"\tfloat64 PRPE(time, lat, lon) ;\n",
"\t\tPRPE:long_name = emissions of PRPE ;\n",
"\t\tPRPE:units = (kg / m2 / s) ;\n",
"\tfloat64 C3H8(time, lat, lon) ;\n",
"\t\tC3H8:long_name = emissions of C3H8 ;\n",
"\t\tC3H8:units = (kg / m2 / s) ;\n",
"\tfloat64 CH2O(time, lat, lon) ;\n",
"\t\tCH2O:long_name = emissions of CH2O ;\n",
"\t\tCH2O:units = (kg / m2 / s) ;\n",
"\tfloat64 C2H6(time, lat, lon) ;\n",
"\t\tC2H6:long_name = emissions of C2H6 ;\n",
"\t\tC2H6:units = (kg / m2 / s) ;\n",
"\tfloat64 CO(time, lat, lon) ;\n",
"\t\tCO:long_name = emissions of CO ;\n",
"\t\tCO:units = (kg / m2 / s) ;\n",
"\n",
"// global attributes:\n",
"\t:file creation time: = 19-Nov-2021 08:29:48 ;\n",
"\t:description: = This file contains emissions from the US EPAs National Emissions Inventory (NEI) for 2011 v2.1. This data has been regridded to the 0.5° x 0.667° resolution. This data is for the onroad_catx sector. ;\n",
"}None\n"
]
}
],
"source": [
"# read in emissions data from the onroad_catx file\n",
"emissions = xr.open_dataset(on_road_emis_filename)\n",
"\n",
"# print out all of the metadata of this xarray dataset\n",
"print(emissions.info())"
]
},
{
"cell_type": "markdown",
"id": "0ee6a372-4c85-4fb5-a3f9-0822925df286",
"metadata": {},
"source": [
"\n",
" \n",
"Each of the variables stored in this file represent a chemical precursor species of one of the three pollutants ($PM_{2.5}$, $O_3$, $NO_2$). Each of these species have three dimensions, one of latitude (89 grid cells), one of longitude (91 grid cells), and one of time (396 days). In this tutorial we will only consider NOx emissions.\n",
"\n",
"In the previous step, you'll notice we read in emissions from the \"onroad_catx\" file. So all of the emissions in this file are specifically emissions from the \"on-road\" sector in California and Texas.\n",
"\n",
"Next we will store the coordinate data, time data, and NOx emissions from this dataset into seperate variables. For this first step of the tutorial let's just look at emissions from a specific day. I've arbitrarily chosen July 1st, 2011 but feel free to try out any day in the date range."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "b27b7143-94ac-4d68-a3ad-7338bb713505",
"metadata": {},
"outputs": [],
"source": [
"# choose a year, month, and day from December 1st, 2010 to December 31st, 2011\n",
"year = 2011\n",
"month = 7\n",
"day = 1\n",
"# create a datetime variable based on this information\n",
"date = datetime.datetime(year,month,day)\n",
"\n",
"# store time data into variable\n",
"time = emissions.time\n",
"\n",
"# get index for given day\n",
"day_index = np.where(time == np.datetime64(date))\n",
"\n",
"# if day index is out of range replace with July 1st, 2011\n",
"if (day_index[0].size) == 0: day_index = 212\n",
"\n",
"# store coordinate data into variables\n",
"lat = emissions.lat\n",
"lon = emissions.lon\n",
"\n",
"# store emissions data into a variable\n",
"nox_emis = emissions.NOx[day_index[0],:,:].squeeze()\n"
]
},
{
"cell_type": "markdown",
"id": "fde390ea",
"metadata": {},
"source": [
"Now that we've stored the coordinate data and emissions data into seperate variables let's plot our emissions data."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "dfb04e07",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(30.0, 45.0)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAAGdCAYAAACo8fERAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAxOAAAMTgF/d4wjAABn3ElEQVR4nO3deZxkdXnv8c9T1VU93dPDDCAEccBxR9yIiluMohKXe9HrFsQdl6hxS0LMZkz0uiTGhSTuIhrEJUQxQVRUoldijLsCgor7CMMiAg7MTHfX+tw/fudUnzp9qurU1l3V/X2/XvWaqrPVqZ7uPt9+fssxd0dERERkUIX1PgERERGZbgoTIiIiMhSFCRERERmKwoSIiIgMRWFCREREhqIwISIiIkNRmBAREZGhKEyMmJmdamZ7Eq/PMrMPJ17f0cy+ZmYVM7tojOdxgpm5mc2M6z0mmZk938x2r/d5rDUze42ZfWUEx3EzO3EU5yQiG5/CRIqZ7YoCwDVmtmxmPzazt5nZzgEP+UfASxKvXwksAncGnjjs+XbxVeDW7l4f43tMrSj0uZl9NLX8RDPz1LLtZna6mf3SzKrRv281s4PW9qxzeQvwuBEc59bAl0dwnIlhZs8xs5+b2ZKZ/ZeZ3TmxLg7f6cdixnHONrO/ip4fYWbnmdmimV1rZq9MbfsQM7vAzH4dHe+OHc7t9ma2z8zKZvYHZvZVM7s52u8TZnb71PZDv6+Zvd3MLjWzevIPHpFBKEwkmNldgG8DhwJPIVzwnw3MAH8yyDHd/WZ3vzmx6PbAV9z9l+5+0wDnWMhTbXD3qrtf1+/x15qZza7j21eAk83sHp02MLN54L+ARwB/ANwReF70+iIzm1uLE83L3fcP8n2VcZzr3L06inOaBGb2cOAM4O+B44HrgM+YWTna5KuEAJV8fBX499RxCsBjgM9Ei/4NOAR4EPBi4K/M7LmJXbYSfqe0Xewz/G/gP6Ov+UOBDwK/S/g+2wJ81sxKie1H8b5N4F3AF3qcm0hv7q5H9CD8UH0DsIx1O6J/HwR8CdgL/Br4V+BWie1OBfYkXp8FfDh6vhvwxOM10fLHAJcRLm4/BZ6V2H9XtO2TgW8CVeC+wEXAm4D3AvuiY5+S2O+EaL+ZPOfd4etxP+Br0XldBfx5ar1Hn/cLhGrLd4B79jjmRYS/nt8H3AK8I1r+rOizV6KvxWMS+/wWcC7hArCP8BfzcanjngD8EFgCPg38BbC7y3mcCuwBPgacl1h+YvixaL3+2+g9D0/tf3i0/FXR6/8VfQ3uktjmP4DPdTmHeeCd0f/H3ui8d6W+dz5CuADeBPwKeA6wA/g4sB/4XvJrAbyGEFbj108FrgCWo6/fGYl1fwz8Ivqa7yH6fkz8356YeJ3ne/TxhO/RA9H/89GJbX4PuDj6/7kB+EyP75M/I3zPVYCvA/fL+L97cnT+e4EPALNdjvfvwEcSr7dG/1+P77D9UUAD+L3U8gcCV0bP7xl97jsn1r8WuCTjePHX6I4d3u/zwPM6rLt1tO89R/2+6d9Reugx6GPdT2BSHsCtCEn9KT22eyRwMuEv1PsCXwE+llh/Kp3DxGGEsPIW4AhgIfphrwCvA+4CvBSoA78T7RP/Mvhh9N53BLZHv6xvJlRM7ki4iCwRXfRYHSa6nnfG59xG+KV/JnBXwkXpAPC0xDYO/Bz4P4QqzqeA7/T4+l1EuAj/ZXQutyMEnTrw8uhr8Nroa7Ir8TV4OXCP6H3eA1wJbInWbwd+Q/gr6xjgRdHr3V3O41TCBemY6L3vFy1Ph4nLgPd0OMZ7ge8lXr8v+v8tEsLRXmBnl3M4G7gw+v+4C+GCeBlQTHzv3EIIE3cG/poQJi8Ang7ciXCR/E7imK8hChOEi9Ayocp22+h9XhitOz76/nkUcHT0f/CM1P/tiYmvf57v0UuAhwF3I4SKT0TrZ6L3+qPoPO4J/HGXr8vTCN9rTyd8751B+F48KPF/t0T4frtH9J43Ai/rcsw9pC7WhHD99x22/+voe6yQWv564N3R8+cDV6XWP5QQQuZSy+Ov0aqLOiHYLAFHdDiXe0T77hzl+2b9jtJDj0Ef634Ck/IA7h/90B3X534PAGqsXABOpUOYiF5/hfa/AN8IfDN1zHOAj0fP418Gz05tcxFwQeL1TPQL+KTo9QkkwkSv885Y/yLg6uT+0bl+K/HaSVQrCH+1ObDQ5et1EfCljM/7sdSyrwNv7nCMIuGv8odEr/8wuljMpI65u8t5tP6fov+jC6Pn6TCxRIcLHyHILSZebyP8pfx2Qph5dpf330W4QB+cWFaK/g8fnDiv72d87nek/h8d2Ba9fg0rYeI+hIv4qv8P4EnAj7p8fyTDRN7v0ZMT658K3BA9PzRaf1TOn6mvA29KfW9fBbwk8X/XBH4rsc17gXO7HLMKPDa17GPA+zts/yPgDRnLL2HlZ+yVpMIzcGzWZ6V7mHg88O0O52HAJ0lUuEb1voltzkJhQo8hH+oz0Scz22lmH4o6cu0Dvkj4ZXfEgIe8C+GXZ9LXouVJF2fse1n8xENHyxsI5fdRnPddCL+wkh04s87rssTzuI/G4WZ2tJntTzye3uWzdP0amFnJzP7OzH5oZnsJF8h5Qik63v+7qXP9ZofPleX/AieY2UP72GcVd99H6FfxUuCr7v7BLpvfjRAeroq/RoQAMkfoVxO7PHH8BuEv8O8n1v8q+vewjPe4lNAM8vOoU/HJiT4CXyBcZH5mZu8xs/9tZtbhXPN+j6a/Fw41s6K730gIH5eb2TlRR8iFDu+16v2i/9dvp97v1+7+q8Tr6+jwvd8vM3sQoRL0wdTy20TLvxgvGsX7EfpLfLrDurcSKhPPSZ7KiN53Q4g6yO+OOprefYTHfa6ZXRZ1UH1pal0h6sD6MzP7qZm9eFTvO60UJlb8jPDLNf0LMu0sQqn2BYRS8ZOj5aVOO/SQ9xfDql7lhMpCktP5//Qs+jvvvOeVPAeP/i0A1wDHJR7nJ7ZLf5Ze7/UXhI6wrwIeHB1vLyvnbon37pu7/4LQnPOGjNU/I5TasxxD6D+Q9DuEcvNRiQt3lgVC1eO41OPOwHmJ7bL+jzt9zds3DBfhEwjNHL8i9LH5qpmVPXQKviehqlMlNLF8ssO5DvO9YNG5PJXQ1PYj4BWEYHFozuP2eq/4/br9Prue1WHjsGh52qnA19z9x6nlJxGqakvR6191OGaTEOx7igLc/2KlQ2dy3d8RmiYf4e7XJlYN/b4bzLmE3wu/HPFxv0P4+n80Y90zCNWgOxP6lv25mR0z4vefKgoTEXe/gdCG+sdZf6GZ2fbo6QOA0939C+5+BaGvxTCuiI6Z9MBo+Sj1e95XAPdJjRzJfV7uXnf3nyYe+3q8V7evwQMIJfVPuPvlRM0DiW1/BNzbzIqJZcfnOc+E1wP3JnQ0TDoXeJqZtf3yjl4/nVAqj5f9NvBXwGMJFYZXd3m/SwnVlbnU1+mn7n5Ln+fekbs33P1L7v4XhF969yGEFjyM+LnA3V8enfNj058zMpLvUXf/hru/GvhtQifSR3TY9EfJ94u+B+/b7/ulfJPQtyI+5jyhafMbyY3MbAvhAnJWxjHSFYRvAjvN7E6JZQ8HLksEjl5+m9B89e3Uebya0Dfi96Kwm/4sw77vhuHuX3b3PenlZnYnM/uMmX0rGgLbV/XA3S919x8SQlraUwh9qRoeRk99DDhloA+wQShMtHspoTLxBTN7pIU5J+5vZm8n9OqH8JfqM6Nv1EfTe8hXL+8G7mVmrzWzO0fltCcD/zTkcdP6Pe+PALPAu83sGDN7KvCyMZwXwNuAJ5rZS6OvwWsJv2TflTj3R5vZvc3s3oTy83Ji/48CBwH/bGZ3MbMXEDoW5ubu10Tvl/6F82bCSJn/jL4njrIwmdPngR8TytBEVYgPEjrnfZZQSXmFmd23w/tdQeg8eY6ZPcrMbhfNDfD2If9ib4m+d/8i+rrdltAptAL80sxOMrOXmNk9LMxh8BTCX7U3ZhxqqO/R6LO9ITqf2wK/T6jM/KTDLv8MvNjMnhb9tfcuQjgbZi6EdwJPMbPnmdndCJWYawidWZOeQPi+/7fUZ9hCuGC3Kgju/j3CyKL3mdm9zOz/EPrRvC2x34KZHUf4KxbgrmZ2nJkdEr0+idD3yRP7/CWhGvcs4DcW5pQ4Iq50jeh94wn0jiMMMT04Wh9vP9WiPyw+Cvypux9PCL8vin5/jMLRtFdCdkfLNq/17rQxaQ9Ce/UHgWsJF6yfEH5IbxOtP57wV+Uy8C1CRzZnZeTBqfTRATNaFg+7qxLK5s9OrNtFRgcqQkfG16eW7QaeHz0/gfbRHF3Pu8PX4n6Etut46GDW0NATe51rr/OOlsdDQ6usHhp6GPBZQvPILwh/AewBTk1s83DCX7TL0bZ/Sc4OmIlltyKMnvDU8h2EC+dVhPL6lcDpRKMLom3+Lnr/ucSyNxP6N2QOWSTMH/BWQkfXCmFkzHviY6S/d9L/x1lfc9o7YN6VMFrkBlaG7v6vaN2DCRekvYROnf8N3L/L/21f36Mkvv8IQ3s/SejXsAz8gMTQ0g5fmz+L/o87Dg1Nbd/63F2O+dzo67dMmDvkLhnbfA7414zljyExciex/Ijosy1Gn++vU+vjr0P6cWq0/hvAEzP+j7P2OWFU75v4WUyv7/gzM+mP6Ot29+j5sdHX5pLE4xdEI5ai7/cbOjzSnVjPAl6aWnYZcHzi9UuAD6z312A9HxZ9IUREpAMzeydwi7v/1QiPeTjhr9vDvXszoORgYfr8k9z98qj69Fl3H7paYGZnEUbbvCOx7DPAWe7+8ej1mwgju14z7PtNKzVziIj0dinZ/SiGsYMw7FhBYvR+BCya2bPiBVGzziFd9unHx4EXmlkxOuZTSDWNbTaqTIiIyNSKqkb/h9D0cwOw393vGHVQ/UdCX4YiYabZp7v71TmP+wzCHCsHE5r3DhDmKrk46pPxNuDR0eb/mKxcbEYKEyIiIjIUNXOIiIjIUBQmREREZCg9b2U9LrOzs37YYVkzAIuIiLS7+uqrq+4+O873uO9xs37d9Y2hj3P1tY1ve5jfYtNYtzBx2GGHsWfPqknLREREVjGzX4/7Pa67vsHPvjP83FNbjvzFrUdwOlNFzRwiIiIyFIUJERERGYrChIiIiAxFYUJERESGojAhIiIiQ1GYEBERkaGs29BQERGRSdJ0Z79X1/s0ppIqEyIiIjIUhQkREREZisKEiIiIDEVhQkRERIaiMCEiIiJDUZgQERGRoWhoqIiICODAzc3hb0G+GakyISIiIkNRmBAREZGhKEyIiIjIUBQmREREZCgKEyIiIjKUvsKEmb3azNzM7p5a/uxo+UmjPT0RERGZdLnDhJndG3gAcGVq+U7ghcDXR3tqIiIiMg1yzTNhZrPAO4GnAV9KrT4D+BPgH0Z7aiIiImunibGvqemXBpG3MvFa4MPu/ovkQjP7Q+D77v6NXgcws9PMbE/82L9//wCnKyIiIpOmZ5gwswcCxwPvSi2/HfAHwN/meSN3P93dd8aPhYWFQc5XREREJkyeysRDgWOAX5jZbmAn8HngQcCRwA+j5Q8A3m9mfzCeUxUREZl+ZrbFzM4zsx+b2SVm9jkz25Wx3bOi9fHjBjP792jdLjOrp9bfIbGvm9n3Eut+d5yfqWfjkLu/EXhj4gR3Aye5++XARxLLLwLe4u6fHv1pioiIbChnAJ91dzezl0avH5ncwN3PBs6OX5vZZSSuu8Bedz+uy3s8yN3XpE+B5pkQERFZQ+6+7O4XuLtHi74O3L7bPmZ2P+C3gPPHfX6D6Lvbqrvv6rD8hGFPRkREZANYMLM9idenu/vpXbZ/OfCpHsd8HvAhd68llh1kZt8CisB5wBvcPXnb04vMrAR8Efgbdz+Q+xP0SWNgRERECLcgv6U5O4pD7Xf3nXk2NLNXAncCXtRlm3ngKYS+irFrgZ3ufr2ZHQL8G/CnwJui9bd19yvNbCvwHuDNwIv7/iQ5qZlDRERkHZjZK4AnAo9x98Uumz4Z+KG7/yBe4O4Vd78+en4T8AHgdxPrr4z+PUAYjTnWDpgKEyIiImvMzE4Dngr8nrvv7bH5c4H3p/Y/PGrCiCeWfCJwcfT64KiagZkVCFWNi0f6AVIUJkRERNZQdBuKtwI7gC9FQze/Ea0708wel9j2DsB9CM0YSQ8GLjazS4HvAtcBb4jWHQN8PVp3GXAo8Mdj+0Coz4SIiMiacvc9gHVY9/zU658B2zK2+3fg3zsc42vAPYc/0/xUmRAREZGhKEyIiIjIUNTMISIiQnTXUN+y3qcxlVSZEBERkaEoTIiIiMhQFCZERERkKAoTIiIiMhSFCRERERmKwoSIiIgMRWFCREREhqJ5JkRERAB3Y19jbr1PYyqpMiEiIiJDUZgQERGRoShMiIiIyFAUJkRERGQoChMiIiIyFIUJERERGYqGhoqIiBBuQb63Mb/epzGVVJkQERGRoShMiIiIyFAUJkRERGQoChMiIiIyFIUJERERGYrChIiIiAxFQ0NFREQAx9jX3LLepzGVFCY2gEeWTlnvU5hIF9bOWe9TEBHZFBQmpoQCQ//G+TVTUBERWaEwMaEUHiZb/P+jUCEiog6YE0lBYnro/0pE+mVmW8zsPDP7sZldYmafM7NdPbb/gZl9O2OdmdkXzeyG1PJnmtmlZnZ5tP7oMXyUFoWJCaOL0/TR/5mIDOAM4C7ufhzw6eh1J28AvtZh3UuB3ckFZnYM8A/AI9397sDZwLuHPN+uFCYmiC5K0+uRpVP0/yciubj7srtf4O4eLfo6cPusbc3sd4E7AR/KWHcn4BTgjalVdwcucfdfRa8/DTzGzA4dxflnUZiYALoQbRz6fxQRYMHM9iQep/XY/uXAp9ILzWwr8E/AH2asKwDvA14C1FKrLwHuY2Z3jF4/CzDgtv18iH6oA+Y608Vn43lk6RR1zBSZQk039jdGMs/EfnffmWdDM3slofLwoozVbwbe6e5XR1WIpFcAX3b3S9L9Ldz9p2b2h8CHzKxIqEzczOrQMTKqTKwjBYmNS/+3ItKLmb0CeCLwGHdfzNjkwcDfmtlu4BzgHmb2/WjdQ4BTo3VfAQ42s91mdjCAu/+7uz/Q3e9H6I+xBfjZuD5LX2HCzF5tZm5md49ef8DMfhT1Rv2ymR03lrPcgHSx2fjUfCUinURNH08Ffs/d92Zt4+73dPdd7r6L0DfiMne/W7TuJHc/Olr3YOA30ba/iY5/6+jfIqEz5js7BJaRyB0mzOzewAOAKxOLzwPuFvVGfRPwsVGe3EalC8zmov9vEUkys53AW4EdwJeiP8i/Ea0708weN4K3+Rcz+wHwI0ITxytHcMyOcvWZMLNZ4J3A04Avxcvd/fzEZl8HbmtmBXdvjvQsNwhdVDYv9aMQkZi77yF0iMxa9/wOyy8C7tth3W7gVqlljx7qJPuUtzLxWuDD7v6LLtv8EXBBpyBhZqcle7fu37+/33MVERGRCdQzTJjZA4HjgXd12eYZwMnACztt4+6nu/vO+LGwsDDI+YpMLVWmRGSjytPM8VDgGOAXZgawE/i8mT3f3T9rZk8BXg08wt2vH9+pTr8La+ds6AuKzc7ilcp6n4aIyECckQ0N3XR6Vibc/Y3ufmSiR+ke4FFRkDgZeD1wortf2fVAAmzcG0PZ7GzbvyIisnkMO2nVR4DrgE9GVQsIFYobhzyuTJF0gFCFQkRkc+k7TETVifh5aaRns0lspOaOTpUIBQoRkc1DM2Cuk43a3JGkJg8Rkc1BYUIGlicsKFCIiGx8ChPraJqrE/2EBAWKFRuleUtEJElhYp1NY6AYJBwoUIiIbFy6Bbn0ZZhQoE6ZIjLJHOOWmuaZGIQqExNgGqsTg1KFQkRk41GYkNx6BYHC3NxIjiMiItNFYULWhQKFiMjGoTAxATZKD//m0tJ6n4KIiKwDhQlZN5uxOrGZ+seIyOahMLHOpqkqoZEYIiKSRUND19E0BYlx0XBREZkUTYcDjfJ6n8ZUUmVCZI2oiUNENiqFiXWiqsSKzdh3QkRkI1GYWAcKEiIispEoTGxANjs7dX/tT9v59ktNHCKykSlMrLG1rEpM2wV62s5XREQChYkNJn1BnrYL9LSdr4iIaGjomhp3VaLThXjahl9O2/n2oiYOkengGIt1DQ0dhCoTkttaVg02UoVCHW5FZKNTZWKD6HXxHfdf+53uGKr7dYiIbHwKE5vIejQfpENGP+FiozV3iIhsVGrm2ADWokmg23t0qkp02jb5GOZ9p4maOkQkZmZbzOw8M/uxmV1iZp8zs10Z2+0ys4vM7GYz+3aHY5mZfdHMbkgs22pm3zCzS6NH5vFHSWFikxn1xbmfINFp/17BYqMEChGRhDOAu7j7ccCno9dptwCvAp7W5TgvBXanli0BJ7r7vdz9XsDngNOHPeFuFCam3HpWJYYNEv0cT4FCRDYKd1929wvc3aNFXwdun7HdTe7+FeBA1nHM7E7AKcAbU/s13X1ftI0BBwHNEX6EVRQmNqGxXZhny+2PAYw6oEwSNXWIbBoLZrYn8Titx/YvBz7VzxuYWQF4H/ASoNZhmy8A1wEnR+8xNuqAuUbGcSGZqKpEVnhILqtUc79nYW4us6OmOmSKyDi5G4v10igOtd/dd+bZ0MxeCdwJeFGf7/EK4Mvufkmn/hDufmIUOv6a0Fzy4j7fIzdVJjapYYJIriCR1mfFolOFQs0dIrJRmNkrgCcCj3H3xT53fwhwqpntBr4CHGxmu83s4ORG7t4kVDCeOYJT7khhYkr1Gl0xrpESAwWJtJzNIRsxUKipQ0QAoqaPpwK/5+57+93f3U9y96PdfRfwYOA37r7L3X9jZr9lZockNj8F+N4ozrsTNXNsYJ2aC/JKX7RHEiSyxMfJaApRk4eIbDRmthN4K/Bz4EuhjyQVd7+/mZ0JnO/u55vZLPAzYBbYbmZ7gA+5+1/1eIudwPvMbAaw6BjP6HFOdwAeEe27BFwK/D93X87zmRQmNph17cBYLkM1f9+IVWbLffWtEBGZRu6+h3CRz1r3/MTzCuHi3ut4u4FbJV5/B7h3nnMxswcAfw/8FvA1QofNWxGGnL7LzD4IvNHdu/5lqjCxRi6snbMmJe7m0tLIAoVXKm3ViVXHrlTbqxPDBIkuhq2wiIhIR38KnObuF6dXmNk8oa/F04Ezux1EfSbW0CjvHtmtxB9feHtdgPM0E6S3WXXMUVcSRtV0IiIiPbn772cFiWjdoru/1927BglQmFhzax0oBtm317ZjDxQiIuvAgaV6aejHNDKz46NKBGZ2spm9xcyOzLu/wsQ6WKtAMc59FChERDaUM4FKNKvmGwgTYf1L3p0VJtbJhbVzRhYqhqkyDPM+YwsUGU0dG3lmTBGRCdBw9wbwGODd0YiRw/PurDCxzkZZpehlFEMpvVJpO05zaak9VKhCISIyjWbN7AjgJOCiaFkx784KExNgFIGiV1AY9ZwMXasUgwSKcqoaoY6YIiJr6R+BK4B97v7daN6JvXl3VpiYEOMMFOOa3GlkgSIdJHKa5pkwRUQmgZkdD+DuZ7r7Dnd/UrRqN3Bi3uP0FSbM7NVm5mZ29+j14Wb2OTP7iZldbmYP7ud40m4cgWLcs0SOtELRozqhfhMiIiP3ajP7qZm9w8xONLMigLs33D33L/HcYcLM7g08ALgysfiNwNfd/U7Ac4CPRNN3yoBG0TEz7tewVtNNDxUoBqxKiIjI8Nz9JOBehH4SzwF+YmYfMrMnxkNF88h14Y/mB38n8DTgS4lVJwO3i07oW2b2K8INRy7KewKSba1mzByVOFDETQ9ts2XGgSJPP4hhp+QWERmQu7FUnc55Iobh7geAc4Fzo4LAw4HHA282s++7++N6HSNvFeG1wIfd/RfRDUkws0OBgrv/OrHdbuDo3J9ANpzkFNw9p9/OQ/frEBFZM+5eBy6MHpjZ/fPs17OZw8weCBwPvCvrfdObdznOaWa2J37s378/z/nJlBv1PTXS/SbUCVNEZDhm9kAze1o0NDS5/Nnu/o08x8jTZ+KhwDHAL8xsN+EOZp8H7he92WGJbW9Le5+KFnc/3d13xo+FhYU85yeijpgiImNiZi8FzgaeClxmZk9MrP6jvMfpGSbc/Y3ufqS773L3XcAe4FHu/lng48BLohM6HjgC+EruTyGSJatTZpfmEZudVYVCRGQwfwDcx90fSygevNHMnhGt69jakDbsyIu/AD5kZj8BqsAzo/YWkbHKui15MlCs1UgWEZEpZ+5+C4C7/8DMHg78ZzRENN2VoaO+w0RUnYif/wp4ZL/HkI0t2QlzpFKdMbMCRSx+f4UKEZGu6mZ2uLtfD+Due8zsEcAXgKPyHkRzQsiGplAhInk5UKlvusviG4E7ANfHC9z9mihQ/N+8B9l0XzVZZ4MMD03qozqRpCYQEZHV3P1jHZZfC7wg73EUJmT6DBgoYqpWiIi0i2a7fBpwexLZwN3/PM/+utGXbAiFubm+h4zGo0C6PcZhmmY2FZFN4z+AJwB14EDikYsqEzKZek2r3WFmzDhQjGqyrG6BQpUNEdlAdrr73QbdWZUJ2ZAGqVT0S3NbiMgGcpmZ3XrQnVWZkMkVT17VqUIRd+Tscu+OfvtT9MtmZ1WhEJGN4HXAN8zsEmA5XujuJ+fZWWFC2oxi1MPI/2LvdJvyaupupB1CxaibPkRkg3KjWtu0l8UPAucD3wUa/e68ab9qslo6BKzJcMpqtXNY6CW5X7Xa8w6j465SiIhMsbK7v3TQndVnQoDe1YS8IxzS6zP7LYzyluKzpfCIg0WPOSx0kzARkUz/Y2b3GHRnVSak72aJkc/TMGh1YraU/RzWrEKh/hIiskE8AHiumf2I9j4T98uzs8LEJjdM/4Z0M0hfxxp6JswQHrwcvoWtWg/LKrX2YNKlH4WaPERkvZjZ24DHAbcF7uHul2dsY8CbgP9F6MdwI/AH7v5TM7sdcC5QjB5XAC9w99+Y2bHARxOH2gEc5O6HdDmlPx7m8yhMbGKj7CiZdayxNSnMllohohmFiri9zmAlUPToR6GOmSKyjs4lBIWvdNnmccBDgOPcvWZmrwL+DjgZuAZ4sLsvAZjZPwF/A5zm7j8AjosPYmbvoMMdQM3sncB5wJeGueu3+kxIm3h+hjXpW9Bv34mof0QySHi5iJeLrVDh5Rn1oxCRiefuX3b3PTk2nQW2RFWKg4A90f6VRJAoAgtAM72zmc0Spsl+f4fjfwl4NvBjM/uwmT0pmlq7L6pMbFJ5KgnJ1/3+9d73BbpXv4lUs0YcJGJeLtIECpUaXp5ZafaIj62RHiKydhbMLBkUTnf30wc4zqeAE4DrgH3A1cBD45VmVga+SWgquZRQyUh7IvALd78k6w3c/VzgXDObAR4GPB74BzP7AaFicb6739DrRBUmNqFBmiTG1iSQp+9EIkisVCBCkGiWQ3GtUG22AgWEkltboEi+Xwf9Bgp1vhTZWByo10dSsN/v7jtHcJx7A8cAtwFuIdwu/B3AqQDuXgWOi0LF24EXEZpOkp5L56pES9TE8Z/RAzO7H+FeHV8Ceo7yUDPHJjNs34Y1awKJRc0a6SDRLBdolgs0okccKlohI7FfK1CUyyG4dAkvavIQkQlyKqEvw153bxImlnpYeqMoVPwL8MzkcjO7LfAg2jtjZjKz+eQDuNzd/8rdcw0XVZjYREbZSbJb34q+j5lVLcjoHwErQQKgUV759o0DRbNcaAsUYZ9UoAAFChGZBj8HHmFmcYn1scDlAGZ2tJltjZ4XCJ0yv5fa/znAf7j73hzvtZ/QlNJ6mNmymX3ZzO7Sa2c1c2wS4xxtMbKLb9xvokv/iGSQaM6uhIlCpUmjXKBYbdIsF6JmjgbN2dJKP4r0e/UY6aE+FCIyLtEoiv8DHAF8wcz2u/sdzexMQj+F84F3Ancl3ISrClwLvDA6xN2BN4Z+mRQI02C/PHF8I1Q2npPzlP6GECj+hTAw7tnAPKG/xnsJfTc6fx73zNEiY7dz507fsydPR9bN7ZGlU0Z2rDUdvpmlV9+IVtUgX0WiOVugUQoRoVhzCpXQY6JYbVKohudWbVCo1KLn0ainSm3lvh5d+lD0ChPD9pm4sHbOUPuLbCZmdvWI+iF0NHPodt/17j8b+jg/e8rfjP1cR83MvuPu90kt+4q7P9jMLuvV3KFmjgk3ygtO1sVv5H99z5Y7P7oZMEg0y+HRKFlbpaKZaAJppjth5jDuIAGjDYoiIkOaN7Pbxy+i54dGL3vOP6EwscmMLVDkCQwQQkPWA9pGXqRHbUB2kGiUrPU8vV16/zikTBIFChGZEK8CvmlmnzezzwHfAF5lZgvAx3vtPHm/XWWVC2vnjPSikzX1dXNpafAmj2SIGPQOoD2Gf3YKEs1W/jAKJQMKreaOZrnQau4YpVEPCY3/b9XsIbLOHBrVYu/tNiB3/4SZfZlwjw4Dvu7u10er/67X/qpMTIlRX2hGVqGIg0SywtB1+1L2g/bKQd4g0SgbjXJ4nmzuSFYnwr6ppo5BQ88YqUohIuvFzB7l7r9290+5+/nufr2ZvSjv/goTU2TiAkUySLSWdQgLidCQeS4dRm70ChLNEiuhIqO5IzlUNK/1HMWhQCEi6+TNyVuQm9kzCRNe5aIwMWXWIlDkkg4SPcJC6/2iiaTSD2ivHuQOEjPRIw4VpdWdMdPvn+c8O57/Gsx6qUAhIuvgFOCjZnakmT0ReAXwmLw7K0xMoXUPFFlBIj5Wh7CQDA1ZskZuhOU9gkQUIpozdGzuiKsTg4zqWC8KFCKylqI7jb4cuBB4HfAod78x7/7qgDmlxtkps2tnzA5BIu9IiW4X9PTIjVxBItHC0qwDGIUwrQTFWv9zqEzSRFXqmCki42Zm6Xt51IGfAKeZGe7+53mOo8rEFFvzCkWOINGcLXV9dHzvVD+JvEGiVZkod27uGLTvxKpzXKcbe6lKISJjdCD1+A/ClN3x61xUmZhy46pQrKpO9AgS/TQhdLqoDxQkEt/B8TDRRtkAT3TGDNNsx+dZrPacf2XiPLJ0iioUIjJy7v5/R3EcVSY2gLFfZPoMEl4udn1kGTRINEtOs+TRMeg4uqP9zqLTmaFVoRAZP68Xhn5MEzM7ucf6I83sQb2OM12fWjqahL9au4WFtPgOn81ygdrCDLVtM9QWitS2FsJj3qhtDY/6HOGxtb2PRJa4v8QgJqm/RCcKFCIyYvc3s++b2f81s0eZ2b3M7AFmdqqZfQz4JLDc6yDT+SeaZIoDxXpccLqFiGZ5dWbNPeyzlB0gVkLDyrwSherKumLVw82/qqFqkbz517RTx0wRGRV3/1MzeythTok/A3YCi4TbmZ8NfMZz3BFUYWIDGnU/im6y+kpkhYdYMkQA1LYW+goRSckqRKEWPeoZ6yrtIcI69JmYttuOqx+FiIyCu18DvD56DETNHLJKnlELWf0Omol+CUmNqD9EuhoxTHNGWqeqRHp4aHw78o1CzR4iMgkUJqQ/GZWITk0cyQABIURkBYn6HFS3Dxgiau1ViWLVV/WbiJs4rNro7+BTQoFCRNabwsQGNfLyd8aNsdLTX8c6hQhgpYNlshqxdWVY5yCSVQmgrSqRbuLYqBQoRGQ9qc+E9C3r7p4A6Tt1wkrfiLiTZW0+3OWzPhcCRH0r1OdCU0RhxijU+x+R0a0qMUo2O7tuE1floY6ZIkNyg5r+xo6ZWdndq3m21VdNhpbVTyKrGrG8I1QkqgeFEFHdHoJEY74ZPZz6nPdVpehVlShukBEc/VCVQkT6ZWafNrODE6/vAHwt7/4KExvYSP5CnU3dXjySHsWR7FyZrEak+0bEzRr1OWjMhyDBfAPmGzTnGzTmm9QO8lbloptCNV9VIu4vUajUwkiODdYJM4sChYj06b+Ab5vZA6OJrC4C0vft6Ch3M4eZXQgcATSBfcDL3P0SM7sv8HZgS/T4F3fPfQIyvdIdL5O3/Y6DRLMEta1Gc2alb0R9zvESNOcbUGoyMx8u7o1qES8VaNYcKNKshaaPmS6jNZPNIn3d2Kuaq3I31TR0VETycvc3m9k3gS8BNwIPcfef5d2/nz4TJ7v7XgAzezzwAeDewPuAV7v7+WZ2CHCFmX06up2pbDDpIaHNVGdLWJk7orbVovAQRmmEakNUjSg5NlenWG4wPxcu7NWZGer1Ao1qkSZAzfBaAVjdlyKuSkB7VWKzN3GkqR+FiORhZruANwMfBO4OvNLMXuLuPWe/hD6aOeIgEdlOqFDEdkT/bgWqwE15jytTIDWSozlbyqxK9NOsYXN1ZudrzM9V2TG/xI75JbbNLTM/V2V2vobN1RNNH6EvRT1x37G4eaNXVWLUs17Gt2mfNmr2EJEe/hs43d1fCDwE+A3wzbw79zWaw8zOBh4WvXx09O9zgE+a2euBw4AXuPt1GfueBpwWv96+fXs/by0DWsvZMLuJb8ZF9G+xHOZ8mJ2pMzfT3oehWpuhWG5Qj26Y0yw1KSZ6WKcnqEpXJbrx8szKBNwZTR3TNgtmP1SlEJEuTnT3HwG4ewN4hZn977w799UB092f5e5HAa8ilEMgzOX9Z+5+NHA34A1mdpeMfU93953xY2FhoZ+3liEMffFIXXRDR8ZGall2BSBrEiqbCduWS2Fa6/mZGvOJQBEvj7frJXkPjvhc0vfi8HKx1WnUyzOhM2m53N7BNN4/eev1DWgSwqXIxKrZ8I/pVDOzJ0WP2wO4+2fy7jzQaA53/yDwMDP7LeAJ7v6xaPnPgW8APW9XKmtr4EBRaQ8S6ftaFKorfRN6VQWapZWL+0wUFObKNeZnqszPVFdVKLopZNxeI9lXIqt5o2ugyAgVG9kjS6coVIgIZlYys48AlwF/A/wtcJmZfcTMcg/UzxUmzOwgMzsy8foJhN6eNwDLZvbQaPmtgAcAl+f+JLJmhq5QpIZU5p2euplsTCs1VzVxbC1W2VqsMj9TY64c3iMOG0QBxDt8Sxdq9DXbZcdAAW2BYqNXJ2IKFSKb3muBWeAodz/O3e8FHB0te13eg+TtM7Ed+ISZzRE6Xv4aOMndG9F41NPNbAYoAW9x92/18UFkDY3qNuWFSq1trolitdk2NDQtHQbKpXpUlahxUCl0Fj7QKLM4U2JpprSq30Q3yYpInluNezmMFilUaqv7UMyWW9WYjdx/Ik39KUQ2rScB90yO2nD3G83sWcClwF/mOUiuyoS7X+Xu93P3e7j7vdz9RHe/JFr3BXe/T7T8WHf/5/4/i6y1vi8aiX4TyaYOqzZWXbw7TmldcmymuVJ1AOZnqiwUl1koLrdVJ+J+E90UEx0uuzVvZBmmQjGtIzryUJVCZG2Y2dvMbLeZuZndvcM2J5jZopldknjMRetuZ2bfiZZdZmYfT81gea6ZXRMdv1snxVrW8E93XySMzsxFM2BuYr0CRa+/ytO38y5Umm0dIbtJNnHEYeKg0jLzM6k+GjPN1giQ1vukwkq/QSKWN1BsNmr6EFkT5wIPBn7ZY7sfRM0P8SP+xXwN8OBo2T2Aqwl9HmLvAY7LcR4NMzs6vTCadyL3L1WFiU2u7wpFh6moc13Io/4S6SaObYXwiKsTczM1ZmfqbRWMVe9XS/SXGGIuiTyBYrP0n0hTqBAZH3f/srvvGWL/ShwszKwILJC4+EetBtfnONTpwGfN7FFmtj16PBr4DPCPec9HYUJ6B4rK6kpXuqmjk2YpzDGRHMkBMDdTazVx7CgusqO4mKhO9NHUETdz9BkomuVC66FA0Z0ChUjfFsxsT+JxWu9dOrqLmX3XzL5lZi9OrjCzspldQhgMcUdCZ8q+uPtZwD8TZrO+KXqcAbzd3T+Q9zi6BbkAg3fMTHbELFab1Ch23DbuLzE7U2c+0cSxrRiqdtsay+wvhqaOG5kPx2x1wlx93EGCRPIOp41yIXQcLRcoQM9OmZupQ2aa7vMhm4JDYTS3IN/v7jtHcJzvAjvd/WYz2wlcYGY3JKZjqALHmVmZcI+sF9HHzbli7n4GcIaZHRa9/nW/x1BlQvLrcnOsuDrRca6JqN/DShNHtdXEAbDNllvViWRTRyxd2ShWVzpe9pKsQkAIEY3E83gb9aHoTs0eImvL3W9x95uj53uAfwV+N2O7KvAvwDP7fQ8zOzp+AHPAXGpZLqpMSP8qNZgtYdV6242/0hWC+Bbi8bDQeH6JZMfLHcVFttkyBxUq7Csusa2x0tSxVK6xb6ZJY6aJJyoT8YRVcWdPqzZW3SsEVlchWsszhrDmrVAUCB1TbXYWr1R6f602IFUpRNaGmd0a+JW7N81sG3AS8P5o3dHAje5+wMwKwMnA9wZ4m+8ADiSn7nTCPBMLZJWFM6gyIUNLj+poNT/kn9CyK69n35cDaHW+tGoj9yRaPd8vqlCEG5pFFYptW1uzZRbm5ijMzWGzsxt6mGg3qlCIDMfM3mlme4CdwBfM7KfR8jPN7HHRZk8izEZ5KfB14D8JFQgId/b8mpl9jxAibgW8PHH886PjA/zIzC7KOg93P8zdD4/+PQw4Eng9sJh4r55UmZA28Y3B4r+641pDAUKpv1oNF9VUdaJQqdEEZm+sUqjMADPUqkYIu0a9VqTOLPuqKyH3QCM0HewtzbOjuMhVtUO4rrqDn+2/FXtu2c6+pS0s3zJL4eYSxUWjfDOUb4HSAWfL3galffVV1ZBklaIQVRsghI64OlGoZE+wlb6XR/z54y3TM+4nqxSt/TZRtUIVCpHBuftLgJdkLH9+4vk7gHd02P8C4IIux39cp3WdmNlTCbNeXg483N1/mHdfhQnpyisVbHaW5tJSdhkrI1CU9odVhWqRcMk1CjWjSpFGzdgX7bpUD+0f+7dsYaG4zDXLO/h1ZYE9t2znN/vmqS+WKNxconSLMXMgBIkte5uUDjRbQSKrGtFPoOjV5yKeLZPWJ2kXB4rWe2/SYCEi08vMHgm8ETgAPNvd/6ffYyhMSE9tgWJuLgwVjasUEYNVgSI0QZSAAoX6SoWiBuwj3GocYLFeZn6myg3LW7lpcT4EiZtnKSwWKd1ilPfCzFIIErM31XNNUJU3UMTLuh1vVZWiPNM2NDYOWemRHpshWKg6ITLdzOxC4A6ECa/Oi5bNx+ujmTB7UpiQXLoGiqjZIx0oCsCWm2pAKRpuZVGfhxAolusFrgeWqiXmyjX2Ls6xuFSmfvMsxZuLzCyFIDF7s1Na9NxBIpYnUPT1NeizStF2Lhs4WChQyEZiI+rrNUVOjP79MKHjZcyi17k6YCpMSG7DBIpibQYoUqjFl+AijVqBZaBeL7BvpkllsYQvzVC8uUg5atqYvdmZ3duktL+xaqbLPB0uewWKfmfPTFcpDEIHzagTardA0TqnDTgSRIFCZDq5+0gGYihMSF/ydMzMChSlfaFZoFgrxHtRrxs1StRrBeqlJiwWKS4WKN8SOluWDoQgEaob0fsNMHV2t0Ax6DFbdx6N3yOxrlOzR9s5KVCIyAaiMCED6dgxs0egKFQKhG+7QjS8M/Sj8FKR4qIxswTlm2F2r1M60DlI9DsMtFOgGES8bzwnRfy8n2YP2JiBQkQ2J4UJGVhmswd0DRQAszfVSQ4dLdSMZglmDoSOlrN7V4Z+xgapHqRlBYq8x+00ARbRMQZt9thogULVCZHNSZNWySr9XAxazR7xBbNSXRnlUam1Rj0UKjWs2qBQDbcLn72pzuzeJrM3exj2uTcM/RxXkIglKxp5jhtPw91IPJqz4dF6nZiK28sz7VNxJya56nhOG2zyK01oJTI9zOyvon8fPsxxFCZkaF6p4JVKX4GitL/OlptqrUAxe7OzZW8zV5AYdqbLXvsn7+WRDhDN2QKNktEoWVuoSAaKtpkzo0ARh4qu56VAISJr7/ejf98yzEHUzCEj45VKKPfnaPLwcpFCtdkaOgpkjtgYl6z7eSRvBNZaFs2U2ShZtI0ltrfo/iBq9siiJg+RqXDAzD4F7DKzj6VXuvvJeQ6iMCEjNXigWG2UzRtZkoEiq0Nm1pTbeSVHe7QmuIruQJpn+KiIrAMnMXx903gc8EjgnsBnBj2IwoRkiu/RMahWp0xYGTYKq6bejgNFa4REhwDRqWkifZOx+BbiecWBImsOiuQ9PIo1p1EKlYi4OlGo+qpbrqfv79FWoajUcgWKjVadgP764YjI2nH33wD/ZmY3uPsXBz2O+kxIRyO5AMTViWq1VepPTkUdh4RCambL+C6gWXcDLVRqrUdacl2nR1ryHGKtOSgS9+5o3Q216q3bn8fbdGqeSXbMJA465Xx9KDYS9aEQmXj/ZWZ/amafNbMLzOxPzCx3wUFhQrq6sHbOwKGirUMmrAoUhdbzRuvfXuEhKwz0K+tYWaEmGSha9/FIVCK6VSXS+g0UG6kzZkyBQmSivRV4OPBe4Izo+el5d1aYkFwGGS7aplJNPO8cKGKjDA/dZAUKYFWgCNuuBIpWlSKqSuS6UZgqFAoUIpPrYcBj3f08dz8PeAJwQt6dFSYkt0EqFG39AlJDRmN5miKSrFof6pE2SKBIP+8lOXIkb6DYiNUJCIFCoUJk4hjtmcBYPbFvRwoT0peh+1F0mIMiS94w0K+sY6UDRbofRbHazG726OMupnH/ifB8c1coQFUKkQnzeeDzZvY0M3sq8Gngs3l3VpiQvuUJFMmmjsxRC6lAMY7QkEc6UHTqRxErdnieV7+BYqNWJ2IKFDJJDCjUh39MqT8HPg48EXgy8B/AX+bdWUNDZSBxoMh7MWgbKpq+dfmghu1PEV3I4+GqsUKl1rrgdxs6CkPecbRSw8szfQ0b3Yg0uZXI+nP3JvCe6NE3VSZkKANfBNIjPAZ5DKutCtG92QOyR3r0q5maeht6Vyg2enUC1I9CZNopTMjQ8gaKVX9tJwPFekkFk16BAlaqEf3cKCx+JJerD8VqChQi00lhQkYiK1DkmsUxGSi6PUYp65gDBoq0TuEhedfRuImkn0CxGaoTMQUKkemjMCFrKrMvQCVHWOgVNvp5pI/ZOo/+AsWqz5Zxfw9ov3FYx327TAOuCoWIjJuZHW9m89Hzk83sLWZ2ZN79FSZkrLKqEx0DRfyYEHmGjialpwSPJYeVprdPHr9tBEsi5GymzphJ6kchsqbOBCpmdifgDUAN+Je8OytMyNjlDhSxZLBYi3DRoTrRS8fpv6NQkQ4WcahIBgurNton6RrzjJ/TSIFC1oxDoTb8Y0o13L0BPAZ4t7v/FXB43p0VJmRN9B0oktLhYg2rF52qE6u36y9YFKrNtm0zqxITVKVZbwoUImM3a2ZHACcBF0XLip03b6cwIWtmqECRlhUw+n0kDVidyNJPsGhr3lBVoisFCpGx+kfgCmCfu3/XzO4A7M27s8KErCmvVFaFiubSUuuxpgb4y7/fG491ChYdmzdUlehKgUJkPNz9THff4e5Pihb9Ajgx7/6aAVPWhVcqmcMd04FiTUcyJGfkrNQ6zpA5qKxOm72mDd+snS+7iQOFZs2UaWZmbwMeB9wWuIe7X56xzcOBvwe2AU3gk8Cr3N3N7B7AOwn9GmrA14CXuXsl2teBy6L9iNb9d49zuj9wB9qzwdl5Po8qE7Ju8sxDkaxajOXCugbViU7HWNW8sZ6Td00hVSlkyp0LPBj4ZZdtfgM81d2PBe4LPBR4arRuGXipux8DHAdsB/40tf+D3P246NErSLwb+FfCfTkeGz1OyvthVJmQddWpQtHJ2CsXOasTowgUrfdYtUyhIi/d10Omlbt/GcCs812+3f3ixPNlM7sEuH30+ieJdQ0z+xZwzBCndCJwrLsvD7Jz7sqEmV1oZt8zs0vM7L/N7LhouZnZa8zsx2Z2uZldNMiJyOaV1Y8ir3TlYqBKxjpcvK1abw8SqkoMTBUKmUALZrYn8Tht2ANGIy2eDFyQsW4r8HzgU6lVF5nZpWZ2erRNN9cOGiSgv8rEye6+F8DMHg98ALg38HLgHsDd3b1qZrce9GRkc0sHilFMId0tUHSsaoy978T03qN4UqlCISPhUBhNrt/v7jtHciTAzA4iBIU3uft3U+tKwL8BF7r7JxOrbuvuV0Yh4j3Am4EXd3mbr5rZx4BzCE0oALj7qvCSJfdvxjhIRLaz0qnjz4AT3L0abXdt3mOKdDOOcJGUeVv0EesYHDpVJRJVEnW+7I8ChWxEZrYN+BxwvrufnlpXAj4GXAv8UXKdu18Z/XvAzN4FnNHjre4f/fuy5GHIqIRk6evPLDM7G3hY9PLRUVo6DHiCmcXDSf7R3f8tY9/TgFapZ/v27f28tcjYw0WbEVUn0tv2U5UozM0pUPQp2eShYCHTzswWCEHi8+7+utS6GUIV4SbgBe7uiXUHAxV3XzSzAvAU4GK6cPeHdVvfS1+jOdz9We5+FPAqQsmkBJSBOXd/AHAycLqZ3T1j39PdfWf8WFhYGOa8ZQKt9S/vYfpa9K3DTcD61dg2R2PbHL4t0cRSTlREUtWRZFPMZrpz6CioL4VMMjN7p5ntAXYCXzCzn0bLzzSzx0Wb/RFwP8If7JdEj7+O1j0FeCJhlMfF0bp3RuuOAb5uZpcShoceCvxxjnN6kpm928zeZWZP6OvzJMJMX8xsifBF2A3cy91/Hi3/GHCBu5/Vbf+dO3f6nj17BnpvmWxr/Ut80ItsW5+JTk0cbRf6lTt7DtJ3ojlbwsthdtp40irbl6g85GzuWLMAtUGoQrExmNnVo+yHkKW0bYff+SWvHvo4P/iH08Z+rqNmZn8LPJ6VeSWeAZzn7q/Ps3+uyoSZHZS8FWmUWG4klFf+FXh0tPxgQor6Xs7zlw3owto5G+cX+Iim2Y6DRLNcoLYwE57PlkKFIg4pHSoU6Y6iqlD0RxUKkVyeDDzY3f/J3f+JMKfFyXl3ztvMsR04z8wui8omLwFOitpoXgk8xswuB/4b+Pt0b1PZnNYqUAzyl3pf81NkBIp+mjriIFFbmKG2bYbmbAgUjYVStG4mO1B0OV8Fiv4oUIj0ZO6+GL9w9wNA50kwUnLVat39KkLFIWvdDYSZskRWubB2zsT9Il8VJPKM4ujQIbOXZJBozhZolOKfzSjHL5Raryw+drkc3m+23Nbcke6QabOzavLog6bhljym+Bbiw/pmNMjiPYRRHH8AfCvvzppOW8Zuw/zy7nNiqawgUdtaaD2aswUa5QKNhRKNbVuyKxQ9go7NzqpK0adJC7ciE+LlhCGmbwPeAVxP+zDRrjSdtqyJOFCs9y/ygaoSSckKRRfN2VIICeUQGmpbQ5hotgoaBlsLFEoG+6NF27ZQ3LfctUIRn3+6U6aqFP3RnBQi7aJmjb8YdH+FCVlT42r2yHOPj5HdxyOuUHQY0dEtSNS2hmaOQg04ELavLRTbAoVVGxSha4fPrDkoFCj6o0AhAmb2++7+cTPLnB3T3d+V5zhq5pA1NzG/wIec8TKrE2avINGcITyi142S0SwbtYVi6JxZLuDlIo14lEeX5o6scKQmj/6sd6VMZALE80Idn/G4b96DqDIh62KtO2aO/O6iGRrbttAsh34QtYUizbJF/SRCmGjOQD261c5MXJXYajSrYfY3ALbNUNpXp0CY4KpQnsH2kdkhM/5cqlAMRx0zZTNz91dH/z5nmOOoMiHrZtS/vDtdQDODxLD34Uh1xswbJJql8KhvhfpceN4oG7V5W9ln28zquShUoRg7VSlkMzOzF5rZ9uj5O8zs22b2kLz7K0zIupqKCa4q1dUPaPVpaEYjMOKmjVgzMXq07flM+/N4Xdzk0VpXTvx4qsljTShQyCb2Ene/2cx+h3An8L8G3pJ3Z4UJmQgTFyjSwSFrfbWKVethWuxqg2K1SaHSpFhzClWnUINi9G+hFm5tXKhBoZ56Hq1v7VdphmNVmxQ6dcIcwx1OJVCg2LzMo5/JIR9TKj7zhwNnu/vn6aMrhMKETIyJCxR57DuA7VuiUKlR3F+jtL/eChSlA81WoJhZWgkUMwdWgkTpgFM64K3tS/sbrSBR3Lfceps89wJRdWJ0Hlk6RaFCNpummZ1CuIHYF6Nluf9qUZiQiTJVgSKuWlSrrUBh1UbvQFGDmaUQJOKKROlAs1WRKO6vtQWJlmRTB6g6sQYUKGQTeSlwCvA+d99tZncGvpR3Z43mkIkz9pEevS7CnZo2Om07W44CBRSrdRrb5ijth2a1ANtm4EAz0Rci/Bs3f5QORM0iiSDRsWkjFk9mFX+WLlNuD2tUlY1pHl2i+ShkM3D3rxPuGoqZGXCtu+eeAVOVCZlIU9ExM5aoUFCpUdy3FG4xXm1S2ldv9YUoLTrFaniUFr0tSJT21/MFiTUy6mm6p725RRUK2ejM7P1mtsPMysAlwK86TWSVRWFCJtokBorm0lLnv/4TgaK4b7kVKEr7G61AUVoM4aJY82hdPYSPHkGi491FU5WWYebUGOe9Pqb9PiIKFLLB3cfd9wKPAi4GjgBemHdnhQmZeJMYKCB1f4xk00gUKKxabwWKYrXZChTJIFGoNnMFiUw5+0/kuYCv5YV+mkOFOmbKBha3xT4E+LS73wI08+6sMCFTYc0CRY/+El37I+QMFKX9jeGCRPL25x1uOtZPdWIUF/bC3FzHR7f3neZQIRtTayj3EI8pdZ2ZvQf4feALZlaCcJugPBQmZGqsdz+KrCCxalmXQFHcXwudLBNDPwfuI5EMFK1l/Y3uyHsx7xYUegWG5P7Dnse49XseChSywTwduAI4JWruuA1wet6dFSZk6qxHoOhWkeharUgEinguio5DP3Nom2+iR/+JThfwfkPEqExaqIjfL/2+ChSyGbn7De7+T9GoDtx9t7uflXd/hQmZSmsZKPoeapluKokDRTQXxbAjNjInsOrQ3BGLL5DrESL6Pf5adALtdfx+A4VChUwrM/tQ9O+3zOyb6Ufe42ieCZlaa33n0SSvVNouOM2lpfYLZDz/RKxahXIZ27eU3URB71ku43uAFCo1vDwTboE+W2rdI6Q1/0Ri7onkvBPdLpBrcVfVfo367qf9BhTdfVU2iX+K/n3FMAdRmJANJX2R70uHzpedbvOd3ibzeFn9GJKViUSwsOrKpP7pYNGcLeHl0BeqyeCBouO5pqXPu5+JvHLIU+0Z9YW83+8NBQnZDNz9O9G//zXMcRQmZMPpu2oQL+ui28ySuf+qj6oTvXQKEvFdRAt0CBQQQkUyUABUqr3PsVPnzeT5jihQ9AoS47yIJ4/dLVgoSMhmY2Z3Jdwp9PYksoG73y/P/goTMtXyNnV0DRQ5L5Lx/n39ld+huSOvZJBoJG5JnhkoYKVK0WXK7dayTrLOr9zf1yrzs6xjiOj2fulQoSCxiXmY6n6T+hhwNvABoNHvzgoTsiHlKmkPeGGMqxS5KxKdmjsypKsSySDRnG3vL50rUEB7lSJLp3CT7NuRPN4AX7dJCxKT8t4iE6Tm7m8edGeFCdmwejZ3DGEtOiw2Z0ttQaJRiieoW12hSC7t2I8iKU+AYCXcGAwUKCY5RIhIm8+Z2aPd/XOD7KwwIZvKKANFX3JUJ5JVibh5IxkkVu48CulAATk6ZnaSMbpkdYVkZuBA0YmChMhE+SLwSTNrABXC3xDu7ofn2VnzTMiGlnXBGuUtuvvS5eKbFSSaqSDRKFnb8+ZsVLWIHl4utoaPto43W+owW2Zp1Tovz7Qe8XnEj9YxkxNl9TnjZut9KhUFCdn0zOxOZvZVM/txNKfDsRnbFMzsLWZ2uZldEd3ZsxytWzCzz5vZDWZ2Q8a+zzCz75nZJWZ2sZk9pscpvRc4Ffht4HjgvtG/uagyIRteVv+JdatQJGVc5NMdLtvCQ2vz7ApF/KpjP4oMWXNbNFc1dRTbjwm5KhRZoU0hQqTlvcAZ7n6WmT0ZeD/wwNQ2zwPuCdwbqAFnAn8EvDl6/SbgRuALyZ3M7BDgXcBd3P1aM3sw8O9AtyrDje5+7qAfRpUJ2bSGqlAM+Fc5leqq/gvJSkC6w2UySDTKRqMcPU9VKJJViswKRUqyAhG/d3sVoth6xK+HqVCoGiGywswOJwSED0eLPgHczsx2pTa9F/AFd6+6uwMXAM8EcPeKu38R2JvxFgXCXx0L0esdwJ4ep/UfZvYiMzvEzObjR97PpMqEbAqdRnf0XaFIXjiTz4fsP5AVJFrNGlGQaK9MOO0VCujVMbPT+8bi4NBalxyKWm0OXKFQiJBNaMHMkhfv0909edOso4Br3L0OoWOCmV0JHA3sTmz3LeAFZvZuQj+GU4Bdvd7c3W8wsxcB3zWzm4A54MQeu/1d9O+7WPkF4+S8c6jChGwaQweKbn+B9zlnReucyjMdg0Sz3B4kmm0/rfkCBaw0e8SymjGSkiEiuWyQQKEgIdNmRLcQ3+/uO3tsk57QIv3DDGHeh9sCXwYOEJozHt7rzc3sIODFwH3d/Udm9ljgXDM7Ng4wq07GfaiWCoUJEXIEirzNGgOGiqwgEQ8FjYNEXJlolgh/Z2AUalCoh3/jyXYKtSLFmlOoevi30gSy59xIz13RSTjG6gpF+MyJUSPRBFnr1slVZDpcBew0sxl3r5uZEaoVVyY3ipo2Xhs9MLNTgB/kOP4jgZvd/UfRcT5lZh+I3uMXyQ3N7Bh3vyJ6PpMMG2b2O+7+P3k+kMKEbCrdJrOKL4B93Wcja7t+z6m8uoqYbN6AECDqW7P3X/lLKmxbqCae16BQb/9rq9MMf93+IisdMNgfbVcNwaI5WwoNs3Enz8R8FvHEXrpZlshq7n69mV0MPAM4C3gSsNvddye3M7MtwBZ332tmtwL+EvibHG/xc+DeZnZ49F4PJBQrr87Y9qOE/hsA30w8B3h76nVHChMy9fq9e2iv2TE7hopKYnrqrOUDSDc5JOeTSDZv1LdCfa7DNL85u3wUakahDnFFA+LgsTpIFFYVQqMKxv6VV1YNM+62NXfAqum7B73x2rSFkEeWTuHC2jnrfRoyPV4InGVmrwRuAZ4NYGZnAue7+/nAduC/orkfisA/ufun4gOY2XeBWwMHR300vuTuz3T375rZ3wMXmVmNMPLjZHfP+mVlHZ5nve5IYUI2hPiXeN5Q0em+DEkdmz5GdMOr5GiK5HTZyWGgySDRmG9mHQZKqZBRCtvZTPi3WG4wM9OkXi9QWSxBrQA1o1AL72e1EDSgPUS0BwwjK1Ck+2NkVScGNY1VDQUKyStqgkgPBcXdn594/ivgmC7H6Fg1cPd/Bv45z6l0eJ71uiOFCdlQBqlSQOdQ0bXpI4esi2khVdlIdnhcVZUoQX0OvATMt997Jw4LEAIDwEy0rFyqMztTZ65c49Ati8zPVFmsl7lxeZ6laol9S1uo1ws0qkWa9QKNKGAAbSEjvI7/ODGapQ6BArpWJwalQCEydluiO4Za6jnAlrwHUZiQDaffQAGjCxXD/CWeHAqabN5ozDvN+QYz8ysVgJlEkIiDA8BcucbcTI35mRqHze7nyC17uXVpL0eVb+Sq6qFcW9vBNcs7+HVlgcV6iZsW56nUZ6jWZlrhwgGvF0IFA2jU4hBhdK1QkKiJJqoTw35dFChExmqeMH9FLPlclQnZ3AYJFNDeTt9pGOkoxZ0v4wmqIKpKJJo3mqUmNldnZqZJuRRCQzo8AMzP1JifqXLrLbewUFzmrluu4ZjydRxTitpMtixyRe3nXLHlCK6qHcJ11R3csnVLK1gs1UssVUtt4QKgUS3SoESVzoGiWM24Y3GiOjFsqFCgkLVgDsXa5roFubvvGsVxFCZkw+q3H0Vann4VA0t0vIynzQZWqhKp5o3Z+Rrb5paZK4fgEFcfAOZnqmwtVjmotMwR5b0cVbopESLC+8ze+udUrr09x5RKHFO6kasbe7iidmioVmzZwf7GFm6pZQcL5mAf0CjNUKVIVqAoxMNFWV2dSBpFpUJEJo/ChGx4g1YpYuMKFfEtxpMdL5NViWTzxra5ZQ6ZX8wMEAvF8Dh2y9UcU7qR2xTnSIaIWPy8cu3tuU1xjtsUF9k/u5dvVHawrzHHVbVDOHLLFvY3tnDt8kEs1sutYAGwjzka0AoUzZlEoABK+0NzR7Ga6MXZYd6NQUKFqhMikyt3mDCzC4EjCL8v9gEvc/dLEuufTRgv+1h3//RoT1NkOMMGChhvpSKe7bI+t7p5Y36uyiHzi9xqy4FWgAA6hIhwkU6GiLRkqFiwMo/Yssh+38sPqr/hmvrB7G3Mc0R5b2gGqW3hQKPcakpJB4q2CgUhULSGisJKdWJEoUKBQmQy9VOZONnd9wKY2eOBDxBNZmFmOwljZr8+4vMTGZlRBApYPf/BMOEibuJolIza1pVA0Zhvtpo3dswvcastB1p9IRaKy2wrLHNU+ca+QkRaOlTcb7YJszdyRe06rm5sZ1thmX3lUKnYWtzR2m+xXKbObMdAAVso7lvGKhmzYHUJFQoUItMrd5iIg0RkOyv3EQI4A/gT4B9Gc1oi4zGqQJGUdWHrFTDCXTuLrY6Xta3tzRuUvNW8ceiWRbYWqxxR3su2wjLHzl7NseVlFqzMICEiLbnvSr+KRa4u3cgVtUPZ12gfwXItB7EIrUDRLK0OFFYttTd3pPtPZISKfqoUChQik6WvPhNmdjbwsOjlo6Nlfwh8392/EaYXF5lsw3bMzKOf6kV6OGhy9MZcOYzQOHLLSpC432wTCBfjYUJElnS/iu2Fvfygusze0jz7G1u41ZYDLNVLVGszNObqNGpGoVakvjW+R0gIFMVqiUKltlKdSExm1f6Gg4eKaQwUIhtVX2HC3Z8Frf4RbzazlwB/APxOr33N7DTgtPj19u3b+ztTkRHr9lfiOKsXzb03UyiXsdkShUqT0oEwj0NzJkx3XVws0CjNUC83WKqWuGF5K1uLO1goLnNN/WCuKKwM96xce/uRBorKtbdvPd/vVX5Q3cI+38K+ZhjtsVgvs1QtUa8X8HqBQq0Q7v1RXZkxM9xkrMNsnZ10CBUa9SEyHQYazeHuHzSz9wAPAo4EfhhVJY4A3m9mr3L396X2OR1o3c99586dm2swr0yVcTSHxLxSoXn9rylWq2ypbKd46DxQii7MRnW7UaVIvbaFG6rFMDwTONAoc92WHVxVO4SrE50u4wAwTKjIDhFhlMfexjzXVcNEVzcuz4fZMxdLsFikuGjMHICZpXADsdKBZusOo9atmaOTVKjoFSimqTqhJo7pUOhwIzzpLleYiO6NvuDu10SvnwDcCHzU3T+S2O4i4C0azSEbwbgDReNX11MkDOIsVLdQObQMUafGQtWobi9SIwydqtZmWNpeYrFe5tqZg7iuuoMflPe2jeQYNFQkg8QVtRpXN3ZwVfVQ9jVD58t4qOiNy/PsXZwL9/dYLFJcLDCzFCoSpQNOaTHc8rxYbbZuAgbkDxJJqQmvprlCoRAhm0HeysR24BNmNkfoePlr4KToXusiG9Y4AwUQAkWlysytDgYWKFRmKNZmWN4R+h4U6kVqtQLL9QLX1AvsnZtjx/wSNyxv5dotIVRcteW6tk6ZeUPF6hCxfVWIiIeG3rC8lb2LcywulfGlmShIrFQlCrXwF11clWgZJEjEcgaKSa5OKEjIZpErTLj7VcD9cmx3wrAnJDJpxh4o9u6lCMxUatitDoqWzlCsGpVaGClRo0S9VmBftcjiUpn5uSpL9VIrVFy7dQc/KO1dFSqyAkU/ISKeuOqmxXkWl8rUF0sUFoutIFGoh+aNYs1b0xAXqs1wJ9GsoaH9muIKhYKEbCaaAVMkh3GPAGns3YstLVEEtlTmKFbnWT6kBBSiURJGtVbEF4vU5xvcslhicanMzEyTm+bmuWF5K7fd+ptwI6/EPTmSVYpOIeLa2o7MAAG0ptWO+0kUbi619ZMo1FaqEsDqJo5BlKN+E8kJr3oEikmrTihIyGajMCHSh7H3o7hqD8XfOjyaDHs+NHnUVvpR1LdCc3GGZsmpz4dgUVkMF/ubFuf55fzBXLN1Bz/ccmTbjb7iIHF1Y2nlnhxRiMiaOrt1Xw5Cf41KssNlIkTEVQmgrYmjrfNlP8rl9ud9BIpJoSAhm5HChEif1qofRbl6MIXqAlCmWCuwvCN0eAyTWxnNpSL1xQJeguWlGSpzq0PFVVvDyI9ttsw+38IPlu/UmiY786ZehPAAtN01dKWfBG3NG3FVIu542feQ0FgyRKSX5wwUk1CdUJCQzUphQmQAyYvGOIJFqx8FYNU5CofMAjOtCa5mluKpt41mGeqLMzRLRZbnSyzPz/Kb+Xlu2jbPLw8czI+3HsGRW/ZyzfKOXLcbB/DoNdEkVOl+Eq3KREZVotBvX4l0kEjcUZVKbWoqFAoSG8AmvAX5qChMiAwpfREZVbho9aM4/DDKHESx2qRRLlBbKNIsh/t5lA4Q7jQaBYv61iLNm4s05kv8erHEb+bnuXbuIHbMH74qPDSqxVWhgVr7LLatGS0X24NEsRp3tvS2X7599ZfIqkYkg0T8uo9AMQnVCZHNSGFCZMRGWbVITnBVOHQ7xdkSpf1FmuUCtW1RpSIZLG4x6nPQLBn1W8rU50rcsn2WW0pbwwFrhVZgKNQKrbt7WqqYUEiFirifBKw0byRDRGYTR7dhod2qEWl9BgoRWXsKEyJjNIpg0ZrgqlKluG0BZkt4eYbi/hJeLlJbmKE5W1gJFgdCpaK2ZDRnjPotMzQ7XKsLOVskklWJ1r5RVWLV3BLQfVhojyDh5fBrKbMTZ45AoeqEyNpTmBBZI8MOL42bPQpzc9hsmWK5jG+bw6ohVDTLhbZmkEKNVv+K+lzv46d1CiBDq1bbA0Wl1hYorFpvBYquKkNMiCUiI1XovYmIjNIwHfW8Ugl/gVeqUK1i+5Yo7lvCqg0K1SbFapPS/galA81wn4yof8PMEK0AzQ4DLeLmlUzdmi1gdRNIqpLRVpWI11WrK/tlBIn4bqMQqhO9bgMvIqOjyoTIOhimSuGVCo1KBZudpbAj3H23SGgesNkSLMQX8gIcaNIst3pG9F2hiINEXKXIahZpzhYoZg0J7XTb8Vi8Lq5SxKEhPZojuS30VZFQk4fI2lCYEFlHw8xZ4ZVKuJ353BywgLFSaiyUizSrBdgWfsTjy3OzZDSH+KkPocIopIbPNcqFweeY6BQq0uuhZ5CY1PknZEp4hz5A0pOaOUTW2UiaPfbth30HsH1LFCq1VrNHaV+d0v4GhWq4q2fpgFPoc3LKZml1U0ezxKomjma5gJeLA3+WzCrGABWJZHNHTE0eIuOlMCEyAdaiH0UhGtLZT6BIdsLM6pC50oSSWNarv0Q3cb+IZP8IGElnSwUKkfFRmBCZEMMGisbevTT33hwuwpUaxX1LFPctU9xfa+uYmSdQJKsRyWaR5szqUNGxE+aoDBAkCnNzqlCIrCGFCZEJcmHtnOGrFHtvDs0elRpWrVOo1FqBoljz3IECVoJEuqmjkapINGfbf5XkGtqZx5AVCQUKmVRmdicz+6qZ/djMvmlmx3bY7h5mdpGZ/dDMfmRmT4yWL5jZ583sBjO7IbXPkdG6H5nZ98zsY2Z2yDg/j8KEyAQaRz+K4v4apX31tkBRqNFXH4puzR6NcurXyWyp8w28eqlURzaPhAKFTKj3Ame4+52BNwHvT29gZvPAecCr3P2uwN2A/45W16L9Tsw4dgN4nbvfxd3vCfwSeOPIP0GCwoTIhBp1P4pCpdbqlFmseRQqvGOgaJagWfLwSDV1xOJZN1vr0oGiXyMMEUkKFDJJzOxw4N7Ah6NFnwBuZ2a7Ups+Dfiau38FwN3r7v7r6HnF3b8I7E0f391/Fe8T+QZw+5F+iBSFCZEJNsp+FBb1oWiN8jjQbAsUSVkTVXWtSiT6Tawa0ZGnOjGqEDFbXnmkKFDIGlowsz2Jx2mp9UcB17h7HcDdHbgSODq13bHAspl92swuMbOzzeywfk7EzIrAS4BPDfZR8tE8EyITbpi5KCAxH8WO7Vi1RHHfMmzbQqHSDPNPbC3AAUhPahVXILxEKKgSQkYcPBrl1fNN9C0jQKTnicgKAZkyAkSa5qGQbgyyJ2Dr335339ljm/QPT1ZP5hLwKOABwDXA64F3AifnOQkzM+BdhOrF2/PsMyiFCZEpkLdC0Sl0eKXSunDbtq1YtUFpP9SYaQWKZhWSgaJQDxUIa7u518rz+DbkwGC3Ic8RJNLLcgcLkcl2FbDTzGbcvR5d9I8iVCeSfgl8yd2vBjCzjwAX9PE+b4uO+3h3H+tsXAoTIhtIVuhYFTAqNQrlGZpAsRpaOgsla5slc+YA1LdCYSa6XXl0x9C4f0X8vFjzMH9Fpbn6NuTd7hyaM0j0pVLtWZ3QrcplErj79WZ2MfAM4CzgScBud9+d2vRjwPPM7CB3vwV4NHBpnvcws7cBdyQEibHfFU99JkQ2uDhgJC+k8Y204kmtWqEgGRiqMLO0Eh6StyEvHfDWPqX9jVaQsGojmoEz0aNz1U29Bv+91jMMxMceR1gRGa0XAi80sx8Dfwk8D8DMzjSzxwG4+5XA3wNfM7NLCSM3XhIfwMy+C3wNODjqm/GhaPnvAC8DdgHfiPpb/Mc4P4wqEyKbSWXl9t+FSo0m4S+K0r46bJuhWXM4ALWtRiHx26GfINGm242+Evq50Dej27B3/YwiE87dfwQ8MGP581OvzwbO7nCMe3dY/j9k98EYG4UJkc0musDHNwZrRKMvCpUmhWhURqv/BO1NG8XqSpAo1nx10wZR1aNTE8eILvQ9A0XG9iIyPgoTIpucVRsr7Z37obZQjEKF04zCRRwkSosrQaK0r95+jG59JKBjkBj0Qp83UPQ6vkZyiAxPYUJkk0s2dwA0awU40KS2tdAaBlqM+1PElYnoNs1x80asa1+JMei3QiHSlbOq0ib5KEyIbEaJpo70fTRCUCjQrHnU3BH3lWi2gkS6eaOtKpFVoRhxVSJ9jE6BQs0bImtDYUJks8kYQpmuTsDKcNFGyVrTb6eDRLJ5o60qkX6/McsKFAoSImtHQ0NFNpFVF9j0ME5WDxdt3RgsI0isElclqp2HaHY8lyEpPIisH4UJkc0q1achzA+xEhDa7jAaBYlYvF3PqsQ6UbAQWVtq5hARrFpv9Z2IR3c0y4VWR8s4SIyqc9q4LvYKESLrQ5UJkU2ibQhksvkh0WEy2ZGy1dyRChIdqxJ9NHFMCg0LFRkNVSZENoF4Su1Hlk4hri2k/5KwxPKsbYYNEZNUNRh1iBjmVvEiG4HChMgmEt/O3GZnwwgICCM74iAwW8Kq9RAoZkurOlkOEiT6DRHpC73Nzva1fz/HHpZCxAbjfdz1VtqomUNkk7mwdk7rotpcWloJANVqKxxYtd7e5FGptW7g1Zouu1IL+1Sr4RipINFcWsodJLxSaT26rRs0DAyzbycKEiIrVJkQ2YTiCgWwan6JmCWWj6sSMcgFPrlPr6rFuPpEKEiItFNlQmSTii+IXqmsBIBKNREUEpWI6DXQMUjkrUQMW2XodKzk8cZRiYgpSIispsqEyCaWp2NmS5cQkcdajJwY53soRIh0psqEiHTvR1GtZgaJfisR00xBQqQ7hQkRAVYCRavZI9mpMvG8V4gYRYfJSaIgIdJb7mYOM7sQOILQX2sf8DJ3v8TMPgD8DrAE3AK83N0vGcO5isiYxR0zvVIJHTPn5qBS7RkeNioFic3G2++AK7n1U5k42d3v6e7HAW8FPhAtPw+4W7T8TcDHRnmCIrK20h0zk0EiXXVQkBAR6KMy4e57Ey+3E02S5+7nJ5Z/HbitmRXcfTST+IvImkt2zNxsFCJE+tfXaA4zOxt4WPTy0Rmb/BFwQVaQMLPTgNPi19u3b+/nrUVkHejCKiJ59NUB092f5e5HAa8C3pxcZ2bPAE4GXthh39PdfWf8WFhYGPScRUREZIIMNJrD3T8IPMzMDgUws6cArwZ+z92vH+H5iYiIyITLFSbM7CAzOzLx+gnAjcBNZnYy8HrgRHe/cjynKSIiIpMqb5+J7cAnzGyO0PHy18BJ7u5m9hHgOuCTZvFNjHmEu9848rMVEREZF0/ch0b6kitMuPtVwP06rCuN9IxERERkqmgGTBERERmKwoSIiIgMRWFCRERkjZnZnczsq2b2YzP7ppkdm7HNA83skujxfTN7r5nNJtafZGZXmNlPzewTZraQWHewmX3EzH5iZj80szeO8/MoTIiIiKy99wJnuPudCbeieH/GNpcCx0e3q7gHcBjRXE5RcHg/8Hh3vyNwLfDXiX0/AFzs7ndy97sC/zyuDwIKEyIiImvKzA4H7g18OFr0CeB2ZrYruZ27L7p7fOexMhCPqAR4DPBtd78iev0u4KnR8e8YHf/0xLGuHf0nWaEwISIiMloLZrYn8Tgttf4o4Bp3rwO4uwNXAkenD2Rmu8zsEuAGwp25z4hWHQ38MrHpbuA2ZlYAjgWuAt5jZt81swvN7LdH9/FW6+veHCIiIhuWO4zmFuT73X1nr3dLvbbMjdx3A8dFzRofBp4IxDfNSR8jVgIeCPyNu7/AzB4FfMrMdsUBZtRUmRAREVlbVwE7zWwGwMKMj0cRqhOZ3H0/IUQ8PVp0JbArscku4OroRpu/jJ5/Kdr384Rmkl4BZ2AKEyIiImsouofVxcAzokVPAnZHVYgWM7uDmZWi52VCVeJ70erPAceb2THR6xezUrH4DnCLmd0z2ve+0fKrR/9pAjVziIiIrL0XAmeZ2SsJfSGeDWBmZwLnu/v5wAnAn5hZg3C9/n/A6wDcfZ+ZPR84L6pwXBYfI7rVxanAmWa2BVgGnpTozDlyChMiIiJrzN1/ROjXkF7+/MTz95M9ZDRefz5wfod136bDbTDGQc0cIiIiMhSFCRERERmKmjlEREQgDA2tVtf7LKaSKhMiIiIyFIUJERERGYrChIiIiAxFYUJERESGojAhIiIiQ1GYEBERkaEoTIiIiMhQNM+EiIhIrKJ5JgahyoSIiIgMRWFCREREhqIwISIiIkNRmBAREZGhKEyIiIjIUBQmREREZCgaGioiIgLgTnNpab3PYiqpMiEiIiJDUZgQERGRoShMiIiIyFAUJkRERGQoChMiIiIyFIUJERERGYqGhoqIiAC445XKep/FVFJlQkRERIaiMCEiIiJDUZgQERGRoeQOE2Z2oZl9z8wuMbP/NrPjouWHm9nnzOwnZna5mT14bGcrIiKyAZjZnczsq2b2YzP7ppkd22G750XX15+Z2RlmNpNYd5KZXWFmPzWzT5jZQmLd/aPr9Y/N7Itmdutxfp5+KhMnu/s93f044K3AB6LlbwS+7u53Ap4DfCT5YUVERGSV9wJnuPudgTcB709vYGa3A14HPBi4I3AE8Lxo3UK0z+Pd/Y7AtcBfR+sM+Ajwx9HxPwucPs4PkztMuPvexMvtQDN6fjLwzmibbwG/InxwERERSTGzw4F7Ax+OFn0CuJ2Z7Upt+mTgP9z9V+7uwHuAp0brHgN8292viF6/K7HuvkDF3S+KXr8XeLyZlUb9WWJ9VRDM7GzgYdHLR5vZoUDB3X+d2Gw3cHTGvqcBpyUWNc3s2v5Od10tAPvX+yT6oPMdL53v+E3bOet8x+uIcb9BhaVv/2f930bRHNA0sz2J16e7e7IycBRwjbvXAdzdzexKwrVzd2K7o4FfJl7vZuX6mrXuNmZWSK9z931mtg+4NXDl4B+rs77ChLs/C8DMng28GXgm4KnNrMO+p5Mos5jZHnff2dfZriOd73jpfMdr2s4Xpu+cdb7jlbo4j4W7Hz/u90i+Xep15rUztV16m/QxBjn+SAw0msPdP8hKhQIzOyyx+raMKfmIiIhsAFcBO+P+hVEfh6NYfe28EtiVeJ28vqbX7QKudvdmep2ZbQO2EfpVjEWuMGFmB5nZkYnXTwBuBG4CPg68JFp+PKEU9ZXRn6qIiMj0c/frgYuBZ0SLngTsdvfdqU0/ATzBzH4rChwvAs6J1n0OON7Mjolevzix7jvAFjM7IXr9QuA8d6+N+KO05G3m2A58wszmCB0vfw2cFLXz/AXwITP7CVAFnhm3A/Uw1p6lY6DzHS+d73hN2/nC9J2zzne8pu18e3khcJaZvRK4BXg2gJmdCZzv7ue7+8/N7NXA/xD++P9/RKM+on4QzwfOiyocl8XHcPemmT0DeE903b6aleAyFhY6iIqIiIgMRjNgioiIyFAUJkRERGQoChMiIiIylLGHCTN7rpldZmZ1M3tpat3Lo/t5xPf8eEpi3RYzOyva93IzO9/MbjWp5xutf6iZfcvMvh/Nl/7AST7faJvDzOxXZnbuuM91mPM1s6eY2cXR+svM7GWTfL7R+ldZmE//Z2b2ugk43/9tZt82s4qZvSW1bhJ/3jqeb7R+zX/ehj3naJtJ+pnr9j0xiT9zvb4n1vxnToK1uIfGdwhTbv9VxrrvA7/j7jeb2VHAd83s6+7+S0JP1wXgntGokfcBfx49Ju58LQyd/SDwGHf/oZltAbaM+VwHPt/ENu8CLiCMQV4Lg57vHsLX9joz2w58x8y+6+7/M4nna2YPIUxte0+gDvyPmX3F3T+/juf7E8K8/r/P6u/NSfx563i+6/jzBoN/jWOT9DPX7Xwn8Weu2/fEev3MCWtQmXD3S939h6zcyyO57ovufnP0/CrCfT2OSmwyD5QsDHtZIHxzT+r5vhj4cLQv7r6cup/JpJ0vZvb0aNl/jfs8E+c00Pm6+/+4+3XR85uBK4DbTer5Ak8BznL3A+5eIdwY76npY6zx+f7Y3S8l/KLNMmk/b93Od11+3qL3GvhrPIE/cx3Pd0J/5rp9fdflZ06CiekzYWYnAgcTUimEG5PcAlxP+OHbDrxjfc5utYzzPRaYM7MvRCXvt5vZ/PqdYbv0+UZ/2Z0G/OV6nlcnGV/f5LpjgQcSxlxPhIzz7Tan/iSa6J+3DBP985Zl0n/mupnEn7kM0/Yzt6EMHSbM7L/N7IYOj6N6HwHM7B7AvwBPcfelaPGJhLnFjyDcnGQv8LcTfL4l4ARC+e2+hF/Gr5ng830f8OfuPtIb/YzxfON1O4FPAi9y92sm/Hy7zam/bufbwcT+vHUwlp83GOs5T+zPXI/jT9zPXBcj/5mTfIbuM+HuvzvM/lHi/TTwXHdPTsP9IuBsd1+OtvsIof32NcO83xjP95fAxe7+m2i7cxhBe/MYz/eBwPvNDEJJe87MPu/ujxrm/cZ4vvFfdl8AXu/uHx/mfWJjPN9uc+oPbNjz7WIif966GMvPG4z1nCfyZ66bSfyZ62IsP3OSz7o2c5jZXQkdkV7g7v+ZWv1z4FEWAU4CLl/rc0zqcb4fBR5mZrPR60cDl67l+aV1O193P8Tdd7n7LuAVwGeH/aU2rG7na2a3Br4I/IOHG82tux7fDx8Hnm1mW6PvieeyMm/+JJq4n7ceJu7nrZdJ/JnrZhJ/5nqYtp+5jcXdx/ogzAe+BzgA/CZ6/tvRuv+Mll2SeDwqWncIcC7wA0Kv+Y8Dh0zq+Ubr/xz4IWGO9H8Ftk/y+SaOcSpw7rjPdcjvh/dF+yTXPWdSzzda/7eEi/TPgb+bgK/vCdHrW4B90fPHResm8eet4/lG69f8523Yc04cY1J+5rp9T0ziz1yv74k1/5nTIzx0bw4REREZysSM5hAREZHppDAhIiIiQ1GYEBERkaEoTIiIiMhQFCZERERkKAoTIiIiMhSFCRERERmKwoSIiIgM5f8DAtIChPkR1KEAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# create a meshgrid based on the coordinate data for the purpose of plotting\n",
"X, Y = np.meshgrid(lon, lat)\n",
"\n",
"# define contour levels\n",
"levels = np.linspace(1E-14,3E-10,50)\n",
"\n",
"# create figure\n",
"plt.figure(figsize=(8, 6), dpi=80)\n",
"\n",
"# plot contour\n",
"cf = plt.contourf(X,Y,nox_emis,levels=levels)\n",
"# plot colorbar\n",
"cb = plt.colorbar()\n",
"\n",
"plt.title('California on-road NOx emissions on '+date.strftime('%m/%d/%Y'))\n",
"cb.set_label('Emissions of NOx (kg/m2/s)')\n",
"\n",
"# limit figure extent to california\n",
"plt.xlim(-128,-110)\n",
"plt.ylim(30,45)"
]
},
{
"cell_type": "markdown",
"id": "225ace44",
"metadata": {},
"source": [
"We can see a map of California with large on-road emissions in high population areas like Los Angeles and the bay area.\n",
"\n",
"Let's now load in the adjoint sensitivities just for $PM_{2.5}$. This filename was stored earlier in the \"pm_sens_filename\" variable."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "7d101fbf",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xarray.Dataset {\n",
"dimensions:\n",
"\tlat = 89 ;\n",
"\tlon = 91 ;\n",
"\ttime = 396 ;\n",
"\n",
"variables:\n",
"\tfloat64 lon(lon) ;\n",
"\t\tlon:units = degrees_east ;\n",
"\t\tlon:long_name = longitude ;\n",
"\t\tlon:comment = centre of grid cell ;\n",
"\tfloat64 lat(lat) ;\n",
"\t\tlat:units = degrees_north ;\n",
"\t\tlat:long_name = latitude ;\n",
"\t\tlat:comment = centre of grid cell ;\n",
"\tdatetime64[ns] time(time) ;\n",
"\t\ttime:long_name = reference time of sst field ;\n",
"\t\ttime:comment = ;\n",
"\t\ttime:axis = T ;\n",
"\tfloat64 BC_PM(time, lat, lon) ;\n",
"\t\tBC_PM:long_name = daily sensitivity of BC_PM ;\n",
"\t\tBC_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 OC_PM(time, lat, lon) ;\n",
"\t\tOC_PM:long_name = daily sensitivity of OC_PM ;\n",
"\t\tOC_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 NH3_PM(time, lat, lon) ;\n",
"\t\tNH3_PM:long_name = daily sensitivity of NH3_PM ;\n",
"\t\tNH3_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 NOx_PM(time, lat, lon) ;\n",
"\t\tNOx_PM:long_name = daily sensitivity of NOx_PM ;\n",
"\t\tNOx_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 SO2_PM(time, lat, lon) ;\n",
"\t\tSO2_PM:long_name = daily sensitivity of SO2_PM ;\n",
"\t\tSO2_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 SOAP_PM(time, lat, lon) ;\n",
"\t\tSOAP_PM:long_name = daily sensitivity of SOAP_PM ;\n",
"\t\tSOAP_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\tfloat64 SO22_PM(time, lat, lon) ;\n",
"\t\tSO22_PM:long_name = daily sensitivity of SO22_PM ;\n",
"\t\tSO22_PM:units = (ppbv) / (kg / m2 / s) ;\n",
"\n",
"// global attributes:\n",
"\t:file creation time: = 19-Nov-2021 08:05:18 ;\n",
"\t:description: = This file contains sensitivity data calculated from the GEOS-Chem adjoint. We calculate the sensitivity of pm with respect to 7 precursor species at the 0.5° x 0.667° resolution from December 1st 2010 to December 31st 2011. ;\n",
"}None\n"
]
}
],
"source": [
"# read in adjoint sensitiviyt fo PM2.5 to precursor species\n",
"pm_sens = xr.open_dataset(pm_sens_filename)\n",
"\n",
"# print out datset information\n",
"print(pm_sens.info())"
]
},
{
"cell_type": "markdown",
"id": "fc5a9eb7",
"metadata": {},
"source": [
"You'll see a list of variables here indicating the sensitivity of $PM_{2.5}$ to various chemical species. Note these variables are in the same units as the NEI emissions.\n",
"\n",
"Though for this example we will focus just on the sensitivity of $PM_{2.5}$ to NOx, you can consider the sensitivity of $PM_{2.5}$ to other chemical speciies by using the corresponding variables in this file. Let's isolate the NOx sensitivity for the same day that we chose earlier and plot the sensitivity."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "f68b7d53",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(30.0, 45.0)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAAGdCAYAAACo8fERAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAxOAAAMTgF/d4wjAABuqklEQVR4nO29eZxkVXn///7cqu6eGXoYdkEGHP3ihkKQgIoxivtGjEYFTVBwiRiMMZKfGqMRv+7RhG9EIYqKiGBcIBLcQCSggluU1QUQZYBhkc2BWbu7qp7fH+fe6lu3b1Xd2rqrup7363VfXffcc849VTO3zqee5znPkZnhOI7jOI7TLdFSD8BxHMdxnNHGxYTjOI7jOD3hYsJxHMdxnJ5wMeE4juM4Tk+4mHAcx3EcpydcTDiO4ziO0xMuJhzHcRzH6QkXE85IIelwSSap3EudTP0zJJ3Vv1E6juOMFy4mhgBJl8aTn0naLOmnkp6Tc/01mXYrJd0fX9svLjtI0lck3S5pi6QrJb20zf0PT90/OTYO5M32zg+BvcysAiDp/ZIubVWnAG8G3picSNog6dheBtlkXN30s15SRdKjM+WXSXpPpuw5kr4X/x/aHL9+Vq9jGASSXiDpV5K2S/q5pCdmrr8rfg5mJF3Wop/3SvpU/Hpa0uckPSDpXkn/Ly0oJT1G0tck3Rb/H39mkz5XSdoq6cGS/kLSxZLuk/QHSd+RdFCmfs/3Lfp+HWdYcTExPPw7sBfwOOAK4L8TgRCzAXhlps2LgQcyZY+L6x4FHAB8DviSpMMLjGFtPIa9gEd0NPpFwsxmzezOXutk6t9vZvf3PrqBMQe8t1UFSa8EvgF8DzgEOBj4H+Bbkl4x8BF2gKRHAf8FfJEwzssJ49w1Va0MnAV8uU13LwC+Gb8+BXgi8CzgZYRn4N2puquAG4C/a9PnM4DrzOx24CnAt+I+n0h4tr6TGWs/7lv0/TrOcGJmfizxAVwKvD91Xga2AG9IXf8YsBnYN1XvQuBDgAH7tej/QuCkFtcPj/sodzDmVwDXAduBO4HTUtdWEb5g7wY2Eia5danrZxC+ON8P3AfcDpyQur4C+DRwF7Atvs+LsmMFjo1fp491mTqPA6oES0V6/N8APp4eT+qzTvd3KfDyeIylVHsBNwPH5nw2ueOKrz0PuBaYAW4EXtXmc15PEJoV4HGp8suA98SvV8ef87/ktP8I8AdgGigRrDanp67/cTyWP2kxhlfFY52Jx/68nP87zwB+BWwCzgN2btHfScDlOZ/l3+fUfQ9wWZN+9gK2AjsAO8ef0bNS118D3JP+d0tdM+CZTfr9FPC+JteSZ/OF8Xnf7tvu/frhxzAfbpkYQiyY5+eAiVTxJuB84GgASQ8GnkyxXzK7ESbtdvwmNvGfF/96zEXSXgSLx4nAI4EjgJ+nqnwSeDhh4nwCQRR8XVIpVeeFhPf3RMIX6L9JOjC+9neESe55wP7AW1hogYHw3v8d+BHzFpVb0xXM7ErCRFh39UjaifArMu+z+wvgDuDv4/7+gjA5rgTS5umnEj7Xc4uOS9K6uK/zgAPjOqdL+pOcPtL8Iu7z/U2uPxtYA/y/nGv/D9gJeLaZVYFjgKMk/ZmkKeBMgqi6PK9jSU8CTgdOjsf8NeC8+L2keRdBRD2NYBF7V4v383iC1QQAM7P4/Akt2uTxAuB7ZraF8P9FBPGXcDGwK7DfwqYteT5BbOYxTRC7yfPUz/s6zsjiYmLIkDQh6R8JvzZ/kLl8JvOujqMJ5teNbfp7CfBo4OwW1e4AXkdwmyQm8csl7dGk/oMJv1K/aWY3m9nPzCzxW68jmHmPisuvB44DHgYclurjVjN7u5ndYGanEczAT4mv7QNcaWY/N7ObzOzbZvY/ZDCzbYRfibNmdmd8VHPG++V4TAkvIgicBROomd0H1ID74/7uM7PtcR+vSlV9FfA1M9vUwbjeAFxtZv9sZteb2SeAcwjCpR0nAs+OJ/cs+yXjzRnLHcD9cR3M7DfAPwKnAZ8gfAe0mvj/DvgvMzs5HvO7gStJxZjEvN3MfmpmPwM+QxBbzdiD8PmnuTsu74QjmHdx7AFsNLO5TJ/JtULE8RBTwP82qfI+4NcEodi3+zrOqONiYnh4m6TNBLPtW4G/MbOrMnUuAnaSdChBVJzZqsN44vkc8Dozu6lZvXiS+KyZXWVmPyD4fTfSOHmmuRq4BvhdvBLiSEmT8bXHECwOtyaBgAQz+0qCoEj4RabPO5n/8v0C8NI4MO+Dkv641fsswJeAJ0naJz4/EvhK/Iu4KGcAL5K0WtJK4CW0+fxzeCTw40zZj+LylpjZjfEYPpBzWW2aZ69/gmCteS1wTCyWmlF0zNemXqf/LYuMp2Niq8ozmLcg5PXZzZbILwC+ZWa1nHu+mSC2j0qJ1n7d11lCJJ0cBzubpMcWqL+TpKtSxw1xoPQuizHeYcTFxPDwaeAgYB8z2zX+td5A/AX2ReDfgAcBFzTrLBYc3wLeamZf7GQg8a+sa4CHNrleIfjKjwJ+T/DL/zAWFNOEOIeDMscjCOb9hPQvOQhfwFHc/0/je/878BCCleT/6+Q9ZMb7K+CXwMsk7UxwV3ypwz5+DNxCEBEvIridvtvhUHqdRN8LHJazGuBGYE3sfmq8YSjbMa6TsBvh36NGo8DLo+iY0/+e9X/LJvyehWJjdxZaK1pxOHCLma1P9bmTpLRrMLlHJ/2mrR11JB0H/F/guWb2y9Slft3XWVrOIbiNby5S2cw2mtlByUGw9H07tmyOJS4mhoc/mNmNeabqDJ8H/hT4Usa0WkfS4whBl+9P3A+dEMc2PIYQ/JeLmVXN7BIzezvBB/7HBNFwNSEAc2X8ftJHXtxDs/7vM7MvmNlfESLjX9Ok6hwhsLAdXyKIn78guFiambFb9XkGwSL0SuDsvF+vbfq4jhAjkuawuLwtZnYrITgwGzvxHUJMyVtymr2FYGX6TqrsVIIl4XjgE5Ie1OK2PY25CT8lxFakeRrwkw76eAGNcQ1XEERM2r3ydOBeGoVUUyTtRgjYvTBT/mqCgH9h7MZJ0/N9naXHzL5vZhuy5ZIeLumbkv5X0tWSjm/SxauBzw52lMNNoaQ+zvBgZtfEX3qb867HJrqLgP8EzpK0Z3xpm8XLHyV9CNjbzF4Vn78Z+C3BF7ya4GbZnSZxFpKeQPhleBHhS/NlhBiKm83s95L+i7Ac9QRCLMQ+cZ33mNm97d6jpLcQluBdRQh2ezZwfZPqNwOPjANG76F5oOmXCf7uVbQPWr0ZeIqkb5L63Ajul/cRRPhbC/SRHdd/AG+R9F7CapZnEwJDn9K0l4V8kPBvVSK2jJjZA/G/4WckbSNYr4xgkn8LwZWxGSBeJvps4AAzu0XSiwm/qv68yf1OBr4v6W8JguRowoT78g7GnOU04I2S3kGwVh1HsJ58IakgaV9gF2BPYIckt0PK9fcCwhc4cfl9kr4IfCye/HcgiK5TE5dEbDnbPzWO/yPpHuDOWMQ/H/hhWvRK+iuCgHsdcEPqebrfzLb16b5F3q+zyMQ/qr4IvNLMrpO0CvixpB+b2RWpeocRAm6bBe2OB82WefixeAeZpaGdXCcshawvDSWsjMguSzTgjFSbM4BLU+dvI0xQ2wlm228AB7YYz6MJE8s9hBiPnwPPT11fQfgldxtBZPyOsMJjZer+ZzV7j8DrCW6WrQSx8hVgz/ja4aSWsRLcKt8guB0WLA3N3ONncfkfZcobxhO3v45gXbg0U/cC4OcF/k0XjCsuT5aGzhJ+uR7Tpp/1hJiXdFmyHPg9mfLnAd8nBH9uiV8/O3V9z/jzfE2qbG9CTMsrW4whWRo6S/OloeVU2bHAhjbv6wiCeJ0h/Lp/Ys6/yYL/x6n/f/eRWXoZf+ZnEKw09xHcZOlxrcvrk/kltl8mtUQ59f8yr82x/bpvu/frx+Id8fP22Pj1/oTvoKtSx03A0Zk2nwY+stRjX+pD8YfhOE4BJF0NfNbMTl7qsYwrkt4KHGxmfUvGpZCx8h7gCRZWIDljiKT1wBFm9gtJjyHEQezbov4OhNVwjzezXlx/I4/HTDhOASTtIun1hCWWna7icPrLzYRf//1kF+BDLiScFNcDWyXVV7VJ2i+zYuNlwDXjLiQAt0w4ThHiXyw7Am8xs88v8XAcx+kjkk4hxA3tSbBQbTaz/SQ9nJD4bV9CnNLdwF+Z2W1xux8QMsp+bmlGPjy4mHAcx3EcpyfczeE4juM4Tk+4mHAcx3EcpyeWLM/E1NSU7b777kt1+5HkntvGNrna0LDb3mObLddxlpTbbrtt1symBnmPQw5aYXfeVem5n9vuqP7MzA7tw5BGhiUTE7vvvjsbNixIOOY04Vmlo9yONAzcARdVi2zU6jhOP5F0d/tavXHnXRXWX7Gu534m9vrtgtT2yx2fnhzHcRzH6QkXE47TIc8qHdW+kuM4zhjhe3M4juM4DiGH+Vx9d3mnE9wy4TiO4zhOT7iYGAHcrD58+L+J4zjOPC4mHKdLXFA4juMEPGbCcXrgWaWjfKmo4ywXDGas9zwT44hbJhzHcRzH6QkXE47TI+7ucBxn3HE3x5DjE9Vo4O4Oxxl9DGO7Lw3tCrdMOE6fcOHnOM644mLCcRzHcZyecDHhOH3ErROO44wjLiYcp8+4oHCc0cSAOaznYxxxMeE4juM4Tk+4mHCcAeDWCcdxxgkXE47jOI7j9ITnmXCcAeG5JxxntDBgu2mphzGSuGXCcQaIuzscxxkHOhITkk6UZJIemyk/Ji4/or/DcxzHcRxn2CksJiQdDDwRuCVTvhY4Dvhxf4fmOMsDt044zugwZ+r5GEcKiQlJU8ApwPGwYBHtacBbgJn+Ds3xSchxHMcZBYpaJt4LnGVmN6ULJf0N8Esz+0m7DiSdIGlDcmzevLmL4TrOaOLC0HGc5UxbMSHpMOBQ4NRM+UOBvwbeXeRGZnaSma1Njunp6W7G6ziO4zjOkFFkaehTgUcBN0kCWAtcCLwNeDDw67h8T+Czkt5lZp8ezHAdx3EcZzAYYruVlnoYI0lbMWFmHwY+nJxLWg8cYWa/AM5OlV8K/KuZfaP/w3Qcx3EcZ1jxPBOO4ziO4/REx2LCzNbFVols+eFulXCc5ngQpuM4yxW3TDiO4zgO8RbkVur5KIKkkyWtz0sEmaojSR+V9EtJ10i6RNJ+8bVpSRdKukfSPTltj47bXCXpSknP6+WzaYeLCcdxHMdZfM4Bngzc3KLOC4GnAAeZ2YHAxcAH42tzwEeAZ2YbSdqFsALzOWZ2EPAm4PN9G3kOvtGX4ziO4ywyZvZ9gHg1ZCumgBWSKsCOwIa4/QxwsaR1OW0iQECSg2GnpN2gcDHhOIuI7yTqOMNLWBo60Y+upiWlJ++TzOykLvr5OnA4cCewCbiNkK6hJWZ2j6Q3AFdIug9YSY4Fo5+4mHAcx3Gc/rLZzNb2oZ+DCXme9gYeIKRp+ARwbKtGknYkbH9xiJldL+nPgHMk7W9mlT6MawEeM+E4juM4w8mxwCVmttHMaoS4h6cVaPds4H4zux7AzL4O7AzsM6iBuphwHMdxnOHkd8AzJCW+lz8DFqRmaNLuYEl7QH1bjIjgJhkI7uZwnEXG4yYcZ3jZboszLUo6BfhzwlYU35W02cz2k/QZ4HwzO5+wW/ejgWslzQJ3AMel+rgC2AvYOY7RuMTMXmlmV0j6EHCppDnCyo8jzWx2UO/HxYTjOI7jLDJm9kbgjTnlr0u9niFsqNmsj4NbXPsY8LEeh1kYd3M4zhLg2TAdx1lOuJhwHMdxHKcn3M3hOEuEx044znBhwOwixUwsN9wy4TiO4zhOT7iYcJwlxGMnHMdZDrg9x3Ecx3EAs76l0x473DLhOEuMWyccxxl1XEw4juM4jtMTLiYcZwhw64TjOKOMx0w4juM4DmFp6JwvDe0Kt0w4juM4jtMTLiYcx3Ecx+kJFxOO4ziO4/SEO4ccZwjwtNqOMwyI7TXPM9ENbplwHMdxHKcn3DLhOI7jOA6SVgJ7AdvM7I5O2rqYcBzHcRzC0tBxS6ctKQJeCfw18DhgI7BC0hzwNeD/mdkN7fpxN4fjOI7jjC8/JIiIfwDWmNneZrYrcCDwI+Azkl7erhO3TDiO4zjO+PLnZvb7bKGZ3QWcCZwpaY92nbhlwnEcx3HGlDwhIWmNpMem6tzVrh+3TDjOEuPLQh1nODDEnJWWehhLgqQLgJcDFeDquOxMM3t3kfZumXAcx3Ec50FmthF4PvDfwMOBFxVt7GLCcRzHcZxkGctTgIvMbA6oFW3sYmJI8S2pHcdxnEXkF7Gr4wjgfySt6qSxx0wMIS4kxgePl3Cc4cGAmfFNp30s8FzgajPbKmlv4B+LNnYxMWS4kHAcx3EWC0nXAucDXzOz85JyM7sNuK1oP+7mGCJcSDiO4ziLzNOBm4D3SPqNpI9Leoakjpa1uJgYElxIOI7jLDVhaWivR6E7SSdLWi/J0jkdcuodIOlSSb+WdL2kv4jLpyVdKOkeSffktNtZ0tmxQPi1pA/n9W9md5vZZ8zsCEImzO8DrwV+I+nM5H7tcDfHEOBCYjzxeAnHGWvOAT4CXNasQhwEeR5wjJldJqkM7Bxfnovb3wt8N6f56cDlZvZXcV97tRuQmW0Gvgp8VdIE8Azgz4H/atfWxcQS40LCcRxn/DCz7wNIalXtL4EfmdllcZsKcHf8ega4WNK6bCNJ+wEHAy9J3a+jXUDNbE7ST8zsgiL1O3JzSDoxbZKRdHpsdrlK0vclHdRJf+OOCwnHcZxlybSkDanjhC772R/YLukb8Tx7pqTdC7a7FfikpCskfUfS4/IqSvqjeB7fJulcSbulLl9cdKCFLROSDgaeCNySKj4PeL2ZVSQdAXwFeETRPscZFxKO4zjDhRls78/S0M1mtrYP/UwAzyHMvbcD7wdOAY4s0O4w4J/N7PWSngN8XdK62LqR5mPACcCPgb8HfiDpmfFqjpZmkzSFxISkqfgN/CVwSVJuZuenqv0YeIikyMwKZ80aN1xEOODxEo7jFOJm4JJ4YkfS2cC3Cra7zcwuATCzCyVNAmuB9Zm6O5rZN+PX/yzpekLSqmcSUm8Uoqib473AWWZ2U4s6bwa+1UxISDohbfbZvHlz0TE6juM4zjjyFeBQSTvG588l3oSrDT8HHpB0IICkQ+LyvLwRqyTVtYCZnQW8m+Di2LXoQNuKCUmHAYcCp7aoczTB7HJcszpmdpKZrU2O6enpomN0HMdxnGWFpFMkbSBYC74r6ca4/DOSXghgZrcAHwJ+JOlq4JnAG1N9XAH8CNg5/qH+hbidETJafkbSNYT5+yXxfhtZLids7lXHzL4MvAtouwIkoYib46nAo4Cb4qjTtcCFkl5nZt+WdBRwIvCMInuejzPu4nAcxxleDDFTW5xFjmb2RlLCIFX+usz5mcCZTfo4uEX/PwMeX2Acr21S/hWCZaQQbT81M/swUE92IWk9cISZ/ULSkYSAkGfGCspxHMdxnBEkzmuxlpQ2MLNfFWnbqwQ7G7gT+O/UWtlnmNm9PfbrOMsWD750HGfYkPQWQnzkfcxvPW7Aw4q071hMmNm61Oux3V7NcRzHWX4UTYe9DHkT8Egzu72bxr43xyLh8RKO4zjOEHNrt0ICPJ224ziO4zhwoqTPEPJYbE8KzaxIXgsXE47jOI7j8CLgzwhZrKtxmVEsSZaLCcdxHMeBMHMu1tLQIeTPgXVmtq2bxmP7qS0mix0vodLCACKrVnNqOo7jOA4AvyVsa94VLibGhKzAcHHhOI7jpPgNYU+O82iMmWia/TqNi4llRp5VolU9FxWO4zgOsIJgnTggVVZ4oy8XE2OOSiUXFI7jOACISm288kxIOtTM/tfMXt1LP55nwilszXAcx3GWHSdKulHSJyQ9U1JXE4KLCcdxHMcZU8zsCOCPgEuBVwO/kfQFSX8R79VRCHdzLCN6sTC4u8NxnHHHbDyXhprZFuAc4BxJZeDphLwTH5X0SzN7Ybs+xu9Tc5rigsJxHGe8MbMK8J34QNITirRzMeE4juM4Y46kr7Jw9cb9wI8k/a+Z1XKa1fGYiWVCv4IoPRjTcRxnLLkL2Ae4LD72JuSbOBL493aN3TLhOI7jOISf5RUb29/YBwKHm9kMgKTTgK8DLwCuatd4bD81pzlunXAcxxk79gBmU+dzwFozmwVm2jV2y8QywCd/x3Ecp0e+B3xT0hcIRpqjgcskTeNiwnEcx3GKImarYzstvhE4DngpIOBC4JNmNgc8sV3jsf3UHMdxHMeps6OZfQL4RFIgaT/gxiKNPWbCcRzHcZyvSZpMTiTtSwjALISLCcdxHMdxzgHOBpD0IODbwD8UbexuDsdxHMchRB3OjtmuoQlmdrKkdZJOAZ4EvMfMvlW0vYuJEcdXcjiO4zjdImn/1OnngNOBi4BfStrfzH5VpB8XE47jOI4zvnwzp+xl8WHAw4p04mLCcRzHccYUM3toP/rxAEzHcRzHianUop6PIkg6WdJ6SSbpsW3qrpD0K0k/S5VNS7pQ0j2S7mnR9vT4HtNNru9QYKxt67iYcBzHcZzF5xzgycDNBep+APhRpmwO+AjwzGaNJP0ZC3cCzfIDSe+S1GChkDQp6XmSzgeOajdAFxOO4ziOs8iY2ffNbEO7epL+FHg48IVM+xkzuxjY2KTdrsCJwAltbvEnhHTZ35F0p6QrJV0H/B54DfB+Mzu93Tg9ZmKE8ZUcjuM4/cNM/VoaOi0pLRROMrOTOu0kdi/8O/BCgqDohFMIyzvvl9S0kpltAz4KfFTSWmAtsBW4PtlBtAguJhzHcRynv2w2s7V96OejwClmdpukwmJC0suAWTP7RoG6pwDnAZfElpK21pI83M3hOI7jOMPJk4F3S1oPfAk4QNIvC7R7GvD0OMBzfVz2S0kH5NS9BDgG+I2ksyS9RNKqTgfqYsJxHMdxhhAzO9DM1pnZOuDlwLVm9pgC7Y43s7WptgCPMbNrc+qeY2ZHE9wonweeDlwj6XxJr5G0W5GxupvDcRzHcQjLHoou7eyV2L3w58CewHclbTaz/SR9BjjfzM4v0McVwF7AznGMxiVm9spuxmNmFULmy4vivh8PvJhgucizaDTgYsLJRaUSVq0u9TAcx3GWJWb2RuCNOeWva1L/UuCQTNnBBe/VPAIzg6SSmVXN7KfAT4F3FGnnbg7HcRzHGWMkPSWOl7gDmJO0XdJPJb1d0s5F+nAxMaL4slDHcRynVyR9k5CL4uvA44EpYHfgb+PX/yPpue36cTfHMkUTZWyu0lsf7upwHGfMmBu/LcjfYWbXZMrmCC6On0r6ENB2/w4XE8sQTfg/q+M4jtOeHCGRvT4H3NCuH3dzLGNcVDiO4zjtkHSEpOfHr58s6WOSXttJHx3NNpJOBN4DHGBmv5C0B3Am8H8Iub3fYGaXddKnM9y4q8NxnHHBEJXqeLk5JL0PeDYwKenphLiJbwLHSnqwmb2vSD+FxYSkg4EnArekij8M/NjMnivpUOAcSf8nXq/qDAH9iJ1wHMdxli0vAh4HrALuBPYxs3slnQpcDhQSE4XcHJKmCJuGHE/jdqZHxuWY2f8Sdhl7crHxO47jOI6zxMyZWcXMHgB+Y2b3ApjZJqCwWbpozMR7gbPM7KakIN7eNDKzu1P11gP7Fr2503/y4iR6jZ3wZaiO4zjLlvQX/BuSFwpbjU4W7aStmJB0GHAocGrOZctWb9HPCZI2JMfmzZuLjtHJ4JP7aPOs0lFLPQTHcfIwqNbU8zFinBhvdY6Z/ShV/nDgi0U7KWKZeCrwKOCmePextcCFhCANJO2eqvsQGmMq6pjZSfHGI2vNbO309HTRMTrOssMFheM4w4CZnWdmW3LKbzCzDxTtp63928w+TAi0BCAWFEfEqzm+Ssgt/p44AHNPwFdzpBiWSaPXQExf1dF/kv8bF1W/vMQjcRxn3Im3Hf9L4GGktIGZva1I+17zTLwdeJKk3wBnAK/0lRyO0xnDIjgdxxlrvkbYJbQCbEkdheg4Mi+1Nzpm9nvC+lRnCBh0kiq3TgwOt1I4znAwhum0E9aa2WO6bewZMMcIz4g5/LiVwnGcJeJaSXt129hnF8cZMp5VOsotFI7jLDbvA34i6Spge1JoZkcWaexiwukId3UsDu72cJzFx4Da6C3t7BefB84HrqCDZFUJLibGDE+vPVq4lcJxnEVi0sz+ttvGHjPhdIwnzVpcPI7CcZxF4HJJB3Tb2C0TI0azibyT4Eq3Towe7vZwHGfAPBF4jaTraYyZeHyRxi4mnK7w2Imlwd0ejjNYKtWxNdj/fS+NXUyMKW6dGF3cSuE4Tr+QdApwHnBJL0knx1aCOb6b6KjjsRSO4/SBS4BjgBsknSXpJXFq7Y5wMTHmuKAYbZ5VOspFheM4XWNm55jZ0cAjCMtDnw5cI+l8Sa+RtFuRftzNMWJYtZo7gdtcpWth4JuAjT7u+nCcPmCiVhvP39ixi+Oi+EDS4wl7dVwCtF3lMZ6fmrMAt1AsD9xK4ThON0halT6AX5jZO8ys0HJRFxNOHRcUywN3fTiO0wWbgU3pQ9J2Sd+X9Mh2jd3N4TTQD5dHGnd/LB3u+nCczlmsdNqSTgZeCDwEOMDMfpFT5+nAh4DVQA34b+BdZmaSpoFzgT8GMLPdUu0eDHwOWAfMANcBbzCz+1oM6Z8JguJzgAhBmauAO4FPAYe3ej9umXAW0M/dRVUqLTj6SV7/g77nqOFWCscZSs4Bngzc3KLOH4BXmNn+wCHAU4FXxNfmgI8Az8xpVwXeZ2aPNLMD43t8uM14/sLMPmZmD5jZ/WZ2MnCEmZ0B7NruzbhlwsllkHko2lkvBjH5t+pzHKwnbqVwnOHCzL4PIDW3hJjZlanX2+MdPR8Wn88AF0tal9Pu98DvU0U/Ad7QZkirJD3MzH4Xj+thzIuItpOBWyacpvTTQtHyPktsRRgny4VbKRxnNJG0J/BS4FsdtisBbwS+3qbqu4CfSrpQ0gUEAfKu2J3y1Xb3ccuE05JxyZSZCIpxsVK4hcJxFmJAtT/ptKclbUidn2RmJ3XbmaQdCWLgI2Z2RQftBJwKbAQ+3qqumZ0r6QfAEwgxEz82s7viyx9sdy8XE05bxkVQwPjkzHC3h+MMlM1mtrYfHUlaDVwAnN+FIDkZ2Ad4kZnV2lWOxUM7C0Yu7uZwCrFYLo9hwN0ejuMMA7GL4QLgQjN7X4dtTwb2A15sZrMF6tckVbNH0fu5mHAKM26CYlxEhQsKx1l8JJ0Su0LWAt+VdGNc/hlJL4yrvRl4PPBiSVfFxztTfVwB/AjYWdIGSV+Iy/8EeBNhaehP4nZfazOk1cCO8fEg4K3AO1u2SL8fMytat6+sXbvWNmzY0L7iiDOIL+pmk9xiTfbj4vJIGAe3B7jLwxluJN3WL9dBM8q7rrF9P/H2nvu56S/fOfCxLgaSLjWzw4vUdcvEMmKxJvlxslDA+Lg93ELhOE6CpIcT4i0K4WLC6YpxFBTjICpcUDjOeCLpbkl3xce9wM+A9xZtP14zgtMWlcN/Cau0t3KM0yqPhHFY7eFLR52xxcCqi5NOewg5JPW6AtxpZh6A6XROIiSyr1u2GTMLBYyHlcItFI4zXpjZzanjtk6EBLhlwqG5cFC5PDQWikGKlm7HvtytFJ6LwnGWP5K+AXzIzC7PubYT8Fpgi5l9slU/LibGnHYWiH4JimG2YCRj60ZULHdBAe72cJxlzjuBD0raH/gxYZfQlcCjCQGYpwKnt+tkeL/hnYFT2JVRMI5imAVDEbq1sLigcJzlgy3SFuTDgpldDbxA0j6EXUnXAluBLwGXmVmhL8XR/vZ3uqaokMi2KWKlGGW6tVK4oHAcZ5Qxs1uBs7pt7wGYY0g3QqIfbUeJbqwsyz0oEzww03GcfFxMOB2x3C0TaTRR7lhUjIugcFHhOE4aFxMjSCtzehHz/DgJgn7ggiIfFxTOsqSm3o8xxMXEmNKNoBhGEaJyOffo+306tFKMk6BwUeE4o4+k/5NTdljR9i4mliFFgwc7EQfDIiSKioZmIqNX4THqK1YGhYsKxxl5LpH0suRE0tvoICDTvxlHFKtW+/Lr1yqVthPqUgqJxQj47HSVStEVH+OwwiOLr/hwRp0xTqf9NODLkp5OWB4qwvbnhXAxsUyxuUpffkUvhZBYihUj3Sx7LZKXYlwFBXjmTMcZJczst5L+AfgOcA9wsJndW7S9uzmWMYNwdwyCQcc7dDoOpz+468NxRgdJJwCfJVgoPgT8SNJTi7b3b84Rpoiro6iFIs/dMWiR0ZeJO3n/ffz134mVwq0T7UkLCrdWOM7Q8izgCbE14oeSfgh8EXhUkcYuJpw6aUExKCHRt1/+WRGVPe9x8u73VuzjLigSXFg4Q40JjenSTuD5ZmbJiZldIenQoo0Lf7NL+g6wJ1ADNgFvMrOrJB0CfBxYER+fM7OPFO3XGTydxE+MnIgoWq/LiXyYdk5dbnhsheMMB5IeBBwv6eC46ArgP8zszqJ9dPINf6SZbYxv/CLCLmIHA58GTjSz8yXtAlwn6Rtm9qsO+na6pOiqjn4FZHbKoouITtoXFBidbHTWcudUt07k4tYKx1k6JD0K+B5wKXAxYRXHE4GrJT3VzK4r0k/hb/pESMSsIVgoEnaK/+4AzAL3Fe13uXNR9ctjF4TW1yDGQSZ/6jDeooiVwgVFb7iwcJxF56PAm83sS+lCSa8A/g14QZFOOvrWl3QmIdIT4Lnx31cD/y3p/cDuwOvzTCNxpOgJyfmaNWs6ubXTB5bKOtE1i5VFMn2fNhN9Wig1ExYuKPqDCwtnSRi/PBOPygoJADP7T0nvLdpJRzOLmb0KQNIxBDXzfOCtwFvN7CuSHgZcKumnZnZ9pu1JwEnJ+dq1aw2nL3SSwGrkBEUbOnkvhWIaOrBWtLJUFBEU0HqfFWceFxaOMzBaqafCyqqrPBNm9nngaXHQxovN7Ctx+e+AnwBP6qZfp3s6mZSWS6BgxxtwxftrZI9cSqX5o1WfrVJ6FxjfuOzh0U+S/BXj5j50nAFxnaQFD5OklwM3FO2k0LexpB2BaTO7PT5/MXAvIUvW9jhI43uSdiMEbvhqjiWgXym2h4JFfB9t02O3sVb0YqEAt1L0glssnL5Ta19lmfFWgkfhxcCPAAP+BDg8PgpR9KfdGuBcSSsJH/XdwBFmVpV0JHCSpDIwAfyrmf1v0QE4S8OouzsGMfbCogIWCItu0nEvuL+Lip5wYeE4nWNmv5b0OOBvgGfHxVcAb0kMCEUo9I1sZrfSZMMPM/su8MdFb+gMlnGIn+h6zFELr15t/udIoY28SqXCgqLTHBQeoNk7vuGY4xQnFg3/3EsfozeTjCCLvTx0HARFIVqJh2Z1OxEVAxYU4FaKXnBB4TjtkfTuFpfNzN5XpB/f6GuZMtKTUAsh1FT4RNHCoxXl8vyR7Sfnni0DNbP1mwRldiPaVCotnziYJcCDNJ2OMFBVPR9FkHSypPWSTNJjW9R7raTfSPqtpNPikAIkTUu6UNI9ku7JafcESVdJukHSxZL2anKL1TnHjsCxwHsKvRlcTCxrigqKUVndkTsZFxEOafIERJ6g6ERU5Kz66KegAF/10Qu+8sMZUs4Bngzc3KyCpIcC74vr7UfY0uK18eU5wmKHZ+a0E3A28Pdm9gjg26RSM6Qxs7emD+BCQuDlnUDhXUNdTDjA8AuKpkKiCM2sEHl1sv13KirS9QYgKFxUdI8LCmeYMLPvm9mGNtVeCnzNzH4fb8L1SeAVcfsZM7sY2JjT7hBgxswujc8/BbxI0kSzG0k6KN6D62TgA2b2JDO7rOj7cTGxzBm5/BP9nCybCYhS1Hi0a9NCVCzsu7igcCvF4uOCwlkkpiVtSB0ntG+Sy740Wi7Wx2UdtTOzTYQNOhe4OiStk3Q2cD7wFeCxZva1Tge6TCPvho+l3KNj1PNPdGSVaGV9yAqHbHk1tcA86ScdTBlFDQGaydgWiLBMYGaveShy2/mKj67xwEynFepPbubNZra2Lz2FvA8JneT6zr6TZm2vA24jWCQmgTcEL0ncidmpRW7mYmJMWMzdRa1S6e9mX/2gmZBoVicRFun3Uak0iphYWOSu+sgkumq182j28y4qLlxQdI8LCmdEuAVYlzp/SFzWUTtJSWDlHTl1v0QQHo/JuVZYWg3ZN/7yZlR2EB225aK546nVOgu8zCMreNITfVZ8VGvNhUWfREVC+r0WyZ7pgqI7XFA4I8C5wGXxhlt3AW8gTP7t+DmwQtLhcdzEccB5ZjaXrWhmx/ZjoB4zscgs5ZfXSEw6izXGPMtJOlAzez0bY5Guk4mpyI2HyKz4ULlcP1rRdg8RPDDTcfpKVb0fBZB0iqQNwFrgu5JujMs/I+mFUN/v6kTgcuC3BEHx2VQfVxBSYO8cx2Z8IW5XA44GPibpBsI24v/Qnw+oyfsJAaKLz9q1a23DhnaBrMuXpbJQdDLp9GKd6MnNkZe7oVV+iTTNAi6zdDO+PGtC4g5JX8vEVeRaF5qIpiIpuVtZK0ZCMA4Zbp0YDSTd1sc4hFzKO+1k+57YKodTMW464R8GPtZhwy0TS4R/gXVGX1eadCt08iwXfbJU1OsWtFY0veYWCsdxlgAXE0vIUgiKxVoq2uumV4Wp9bDFX3qb8aJHmkRAZN0fCX0QFd3kqnBB0RlLGcfkCbWcYUPSgyTt0Wk7FxNLjFsocmjmAuhG3OTlkYDu81nkiYtWVgpY4IrpRFRAc2tFq1gKFxSdsRQTevqeLiqGB1V7P0YRSY+WdC1hqej1kq6R9Kii7YcnZH+MGeaNwJYVPcVx5Ojuai1HAFTmr6VzVRTdSCyz+iNNs5UgzXJV+EqP4aTVs+7bqDtLyKnAh8zsiwCSXg78B/C0Io3dMjEkDOsXx1BkxUzRl/GkBUA2G2azIzuOcuZaMytFugxys2l2Y6lYUNbCQjGWwrELBi3oO7U+uLXCWWR2ToQEgJl9CdipaGMXE0PEsAqKblm0uIlm5E3oWSFRACtHC46kvN5PKZoXAM1iKfooKjoRFOBuj6WkV1HgosJZJKqS9k9OJD0SKByU5m6OIWOxXB6duDqGPolVNwmsWgiJukjIknMPi4ehSi30mbg+qtUgHiqVxnTd2TTdPbg/8twerdJzu9ujPdlnrxeB3+/nOOlvuf3oGCoMVOskY/Wy4p3A9yRdSch8eRDwyqKNh2eGcOqMSqbMgVKt9m/TrzZBl52Ih3R9VWr1/TrqZc3GkIiKdJru7L4fsEBUtNv3Axbu/ZErRpJrLig6opsYhsVwlyS4sHD6hZldEFsmnkD4KvuRmd1TtL2LiSFlMQTFUFsn2oyrp/GkrBKdWSGiJudREBakLBVNb54ToAlNLRW5qbULCIqkrQuK/jGMAt9Tgjv9QtKpZnY88I2csra4mHCWF5VK+1Ub1VrreIk2LpOsqFAlZZmoJHXia80HOj8WKOz+cEHhZHH3R38Z1aWdfeCJOWWHFW3sYmLMGcplot1YS/qx8Vere6WCLi2diCrl4kiTFRWUo8a4ijqdiYoGF4YLCieFWymcbpD0MuBIYJ2kr6QurQG2FO3HxYSzqNuTt2WQwqbLOIwFqzcAKwtVrICwKOoC6VxUuKBwsrigcLrgBuCbwOPjvwkPABcX7cTFhAMMiYWim4m+iMAp4vpIiPLjIup/owgrKy6blwVZYZH3ZDVYK8rzIqPpyg/IFxUuKJwWuKBwOsHMrgaulvRNM7u72348z4QzMDraOXSphUwBEiFhUeMBQVgkR71ulHKNpHNUxDkm5hNfZfJTwMJkWdmNxGLqQionH0WnuSic5cMwBouOAgJU6/0YRXoREuBiwklR5NfpsGXEbEm1t6c6GyeRCAkAK0VYKbFWZARGSlgkosKinMRXSdlUeV5UQHtRAQ1JrxoEQpeCYsmtUk7f8WRXzmLiYsJpYEnM3f2eyLrNvBmlYyJaC4l6vVhUNBUXOcKiNllub6WARtdMq51JcUHhNMcFhbMYuJhwlpalmsBaWC3yhET9WilqFAhZ90YnwqJTKwW4oHC6wgWFM2hcTDgLaGed6Juro98TV60Ht0ZO4GU6oNIizQuJnJiJ0CYjGDJWiwXColMrRautzl1QOG1wQVGMcYuZiJeGJq93k/RNSfdLulTSvkX7cTHhjDwdiZt2QimxFtSDJ9UgJObrpYMt2wdlLhAXUXsrRcMGYgmtdiVNCYqGwMwULijGGxcUTg7vSL3+IHAt8EjgfOBjRTvx0G4nl16XirZdybHYk1WB5aHpHUHD+UIhkYiAdJ00qli9XkN5zTJ1I+ob8pWBSoRNRihZ9lkBmwzLQBt6yi4jze71UWDpaNFlo75kdHniGTOdDOmvmCcAB5tZFThJ0jFFO3HLxBAzzA/7yKzq6GJFh0X5j0WDSyNaeECONaKZ9SJtqYjS9Rq3OW9qpShgoWjALRSO0x4L6bR7PUaMKUmPjjf5qsVCIsGKduKWCWcgWKXSWZ6JHsjNnVAwUVV6X41esWbdlHNyXlaAcpzsqhRbKZLhVqKGJ1OVkOiqno4bmlsoEmILRav9PPIsFI7jjB2rCJkvBSBprZltkLSGuvm0PS4mnKa0c3Us+k6i3ZJs7JUIjCStdrsNv3LIujUKtcm7RVl1QREIgqLuDokFRSJ2CguK9NbmXQgKd3c4znhhZuuaXJoDXlK0H3dzOMuHIht9JRNni8lRvawKKYhF8zEZwLzbI7PaIzcwE1q7PAaQ+dLdHcsPD8Z0WmFmW83spqL1XUwMOUsdN7Foy0QHTcHYifp+GS1o6s5oQa2kBUdaUKTjKMI95gUFUFxQJLRLvZ0i645qJjxcUCw/XFAsZNyWhrZC0mlF67qYcEaWhiWQeWTjARJBkWedqNUKCYlOSYQDLAzYTAuKcL25oGgIzGwmKJLzLF0EZLqgGB9cUDgt+HrRii4mnLZ06y8fVHBfWxHRCdVFFBEl1Y+kvMFCUe5gpUeeoGjl7iiQ1KooLiiWHy4onDzMzMWEs3gslqujpYgoEi8BC60TfSTtwoCFIgKgVg5HXVikBAXQICiAlgmuOhYUMc0ERVHrBLigWI64oBhfJO0j6UmSpjLlzyrah4uJEWCp4yaGga4tEYuw9DEtIBLq4iArIuL4iFp5/nqDhSJqFBTpwEygN0GRFz8BLQVFK1xQLD9cUCxenglJD5f0Q0k3SPppnOchW0eSPirpl5KukXSJpP1S14+QdJ2kGyWdK2k6de3ouM1Vkq6U9Lwm4zgS+DnwKeBGSYelLv9LsXfTgZiQ9J3UwH4g6aDUm31P/IH8QtKlRft0RoelCsTsq0sjIW2d6NKFk7VCJGStEVkRYSXCEanRShELinofTVZ61M+LCooszQRFE0Zi6a/TV1xQLBqfAk4zs0cAHwE+m1PnhcBTgIPM7EDgYkLKa2Lh8FngRWa2H3AH8M742i7AqcBzzOwg4E3A55uM4x8JWS8PAF4NfFnS0+NrhdfCd2KZONLMDowH9m/A6XH53wEHAI81s8cCr+igT8fJZSAiIk16dUcSN9HDktDiIiKpN59RMysoWq306FhQtIifgPwVHm6dcFxQDBZJewAHA2fFRecCD5W0Lqf6FLBCkoAdgQ1x+fOAn5nZdfH5qczPvxFBCCSWip1S7bJEZrYBwMy+C7wA+Gzs4iicAbOwmDCzjanTdGastwJvN7PZuN4dRft0ijMMro6+WSdaTEBdiYh28RIduDpUCcKiaK6JTkUEpFwapUa3h6ViJ1qt9Oi3oKjTRFC0+/dwQbE8GVdB0aelodOSNqSOEzK32Qe43cwqAGZmwC1AdpfOrwOXAHcSLA/PAN4dX9sXuDlVdz2wt6TIzO4B3gBcIelmwo//Y5u85ZKk1cmJmV1LEBSnAQ8p+LF1FjMh6UxJtwLvB46RtCOwO/BiST+Oj9z/gZJOSH+4mzdv7uTWzojSyYqOriwR7YREepvuLG1yTySbdvWTBW6ReB7OCoq2Kz0yG4nlCoo0XQRkNq2Td90FxbJkXAVFH9hsZmtTx0k5dbJfMHkuhYOBRwF7Aw8muDk+0aKP0FGYm48HDjGzhwCvBc6RlPcgnwYc0jAws18BzyfEUhSiIzFhZq8ys32AdwEfBSaASWClmT0ROJKw09hjc9qelP5wp6ens1WcEaAvaZX7lZq5mZBIBETbnUvb//dXrf+CIkueoAjn+Ss9QpvUKo9Wgqpg/EQei7W3iuOMIbcCa5PJPXZh7EOwTqQ5FrjEzDaaWY0Q9/C0+NotwLpU3XXAbXG9ZwP3m9n1UF/iuXN8jwbM7ONmdklO+a/N7NlF31BXqznMLP2GNhP7fczsFuByMirH6Q/D4OqA1oJiyTJiFhUQXaBqrWcrRZIVL6oaUXVhX0UFBbTeH6TBOtFMLHWZbtuDMR2nP5jZXcCVwNFx0UuA9Wa2PlP1d8AzJE3E538G/CJ+fQFwqKRHxefHA19KtTs4js0gXqERAbdlxyLp+a2Oou+p0LdDbDKZNrPb4/MXA/cC9wH/CTwXOFXSzsDjgQ8XHYAzmrTbBKyh7iB2EI2am+0HhWoGFbBJEVWNWkmoFiZ8VQ0riagSYiZUMywSqs4LhXZYKSwrC4LCiFDoN4IaIsIyG4RF1HcbTW0M1pQWu4s21it1ZT3yDcGckceKL+3sA8cBZ0j6J+AB4BgASZ8Bzjez84FTgEcD10qaJcRNHAdgZpskvQ44L7ZwXJv0YWZXSPoQcKmkOcKmXUcmsY0ZvgFcQ5jPs79UDPhWkTdT9Jt4DXCupJWEwMu7gSPMzOIP4nOSjo/rfsjMrijYrzPCNBMUhXYTTXbu7LQdLImQSJMIiF77yKOooNBsB5aS5HNutrtoMqa8HUML7CrqOE7nxC6Iw3LKX5d6PQP8dYs+zgfOb3LtY8DHCgzl/cBRwO+B083sogJtFlDo29jMbiVYHPKu3UMwvThjSFELRd+sE0ssJNJ0ap0oKkLaCYpS3G8z6wQwv2U5tN/krEPrhAsKx1k+mNm7JZ0IPBN4jaRPEDwOn+pkdaZnwBwxhiVuIk2eaXsgk03ucsZo4bEMaBZDEcoagzEL0+VW5cOed0KlUu7hOE4xLHCRmb0COIEQf3F0m2YNLP3PO2dZ0EkMRZ0cV0dT0kIiOylmaSUoCm5Fnias6KgB4Ze/MW+NaGad6Ad5FgqwurtDJMGYbawT1VpzS0Pi7sizTrRgKawTnf7/Sup7HIfTCYsYMzE0SNqdEG9xDCFI803A1zrpw8WE0zeygiIbA9G1q6MTIdGOgu1UmZ+cs6snVLGWKyqAlq6OTmgmKKjGQiIVjKlak2/BdOBlIh7ygjGT95cIhZQAycZODJp+WhZcVDhOcyT9F7A/YVXm85JsmJ2yPGzCY8YwujoWhTwhMoAloUW3JE8v9xwkWZcHLMyQmVxfkBUzTTJBF9gIrAiDWCo6SBeFuz8cJ5cXAbsBbyZkzLwrPu6WdFfRTtwy4SwqC6wT7VwdeQGXzXz+6fNuf0XXah1PqgntAjF7IbFQhNfz7o68YEylk+JF0byroxvaWCdyV4B0uUR0sSZ6t1R0xrNKR43VDxgNPk/dsPHQfnTilgmnr2S/oHvyq3ciJLIMKIFVQpLEqhPrRK++2PRuo9Zkh1Fo3LOjTvK5FbRONLM6DCIr5lJZDNxK4ThgZjcDq4E/BspmdnP6KNqPi4kRZZx+KbSlVGo80qTdIEUERrO9LdqQ3XNjUDQIirz9O1JptoHW+3V0SVZQ5K4AabWZ2xCtuljq+zvOUhPniPoB8Hbg53FSyo5xN4fTdzoOxOxkVUcR8vpKW0waLB0lKC2cbK0cxb/uk90/GwMxQ3l4XUtZCWDeapDsHhrKFg6pt8RXivucX92RdneoVg17diQrO2Zr80GXiesiG4yZWdnREIgJDZ9hXjIraLRE+US9PPAfLsue44EDzGyDpAOA/6DDlRzglomRZpgf8nbujmYrA3LdInl1O42JyFov8qwYqcDFtJBoGF89yDGc5wmJ9DbkoSzdfv5oRa2kpkf93qldRmHe3dH4N7ObaC/ujsznpXK5kJXCcXplsXYvVZxOu9djxJhLVnDE24/v0E0n/uQ7i0ZLC0WedSIvGLJa629iqlLjZJvnDkisEnlCIs8aAQuFRJ54KOIaybNo1FCc+2J+uWitRsPeHVZqzD1h5WhhVsyshaIIOfkqek25XUSADDKnhe8p0pxh+MHi26APnClJj2Z+X46G83g78ra4mBhxLqp+eWgftrxEVkVyT+Tu0VGpNI936NZNkhISDffPcW+k4xCgvZDIWiPS7VpR1B1i9bkvCApVG/fusFpjMiurxZaGZoIi/jxaujsSunB79GqxyLb3dN7jwbB+ty0zVrFwM6/k3ICHFenExYSz6DQVFEWtEwmtBEY7MtaNPPeG5VhAsluCt4uPaCUkWu0mmhUQ6bqqhb4i5leUtFouuiB+Ik9QQMtkVnkxEVkrRSIK80RFv3FxMXiW2irhQmJxMLN1/ejHYyaWAUv90LeiW/Nxy8mh25wJCbFIyLo3EiHRMI6MeyMrJJrFR6TjIpJYh+R6cmTJxlOk66avWRTK6is5UstF69dbxE/kxlBAozBrkshqgTjIiT0ZxPLRdmii3HA4o81SCgnVej/GEX/qnIHTkbujlcsizxLRiXUiLSKg6ZLJPPdGMyGR1A9lyfl8X2lrRFZANAvCzLNsNJSlvqy6iZ/ItVDAfPwElVx3R0JuTEQBK0UzmomPxUzf7TQyzD9QnOHELRPLhFF8+NuZpuvXi2xA1c4CkiMk6vdpsgw0GydRr58REmlLQzsh0Ww1R9ZikbZAWAlMqSMpj1J1SylRkxU+mfwTTS0UaRGXdfHkWChyrRQZkhUfrY5mtLveCrdOdM9Sf5e4e2M0cTHhLArN3B1pQVH/JdpKGNTrdGBLbCIkuomTSISERaJIoGUzoZG+lic0siIiqV9vl6ozP7bG5aK1UsbNUVRQ5Lk7ks+tC7dHr/QiKpzRYhiExLgtDZX08fhvR1uOZ3ExsYxY6l8U7SgiKJrSwfbYDbQRErnjaeLeSAsJKBYfUe8zdatszMQCAZG2RJRyjqS8m/iJTgRFs9TlXVopeqVTUeHWic5Zyu+QYRASY8qfxn9P6KUTf9qcoSIbO5G7TDQhyTmRxE1k4y3SmS2jRlFR/9vCvZF1FzSMs4P4iKyQyJY3lCmnLM+yUQVKjfETVJO6IX6iVE2t+shsV14ohgJYED+RkImjgJxYim4ERYGA3U7iMZzRwIXEknKrpGuBh0r6afaimT2+SCcuJpYZw5x3AsgNxgRai4aEZJlokaDLFomtigqJhjZp90YqXqGdgOiXeMhb/aHKQjET1XcMXZhuu+fdQ5LPPJnEEwtFJjgTClqb8sj+32ghLvJ2MXW6Z6msEsP8fTUm/AVwMHAW8NZuO3ExsQwZVUGxgCLWiVb3aWKVgGJCIjdOogTVieaxEq3EQ17sQ9E+0ufRXPw6tk7Ur9XmE1epGuInIkL+CarB3dGTdSIhz0oB/RUVzrJnWL+jxm1pp5nNAT+R9EIz+3W3/biYcJaEVoIiLytmU/LSaxexSrRxbQAthUQSq1DvN0cMNBMPrawW2T5y30MtVScWFKqlrCWZ5aJ57g7Rfsv0BpLPNBEVWSsFDM5SMYhcJU4Di22VGFYhMeZsiIMxnwkY8F3gnWa2qUhjFxPLlGG3TsBCQdF3Vwc0t0p0KCRqE7GlogS1iWIuizzx0E405MZnpMpUC/ePSMVLVDN1Gs4b3R2lmsWfR2O67bbWiUQ0LJWoGCC+L8fiMezfSWPOqcBW4BXx+evjslcWaexiwllSiro8unV1NPSRXgbahZCoTVAXE1Dc8tBSTDQxojSzTFgpuDksmhcVAJYyNLRyd1jEgmBM1QpMplnXxpCLimEULeOOC4mh50Az+6PU+fGSri7a2MWEs+S0FBRdbuLVYIlILQNNbyverZBIlmZCZwIi6xqpjy/qzOWQjCsRFdAYP1HE3aFMfwusE7M5juM8wZDdy6OZqMhZ+QEFJ/0eXB3tcKvE4rg4RkZI2OjliegjJUmrE7eGpB3oIH2Ei4llzCi4OprRMm6iqKsjbzfQtHuj1X4bE/PXskKiNtlaPDScR+kyaxkLUafUQlxURQ0jmotdLgRRUZuEaJaG5aKt3B2Ks2angzEb4iiiaKGrI5l4s4Iha6XIq5NjpYDhdn84/WFUv4PGkDOBH0s6mxAz8XLg80Ubu5hwhopW7oy2ro68JafZLJdp90YbIVGdTGWWzAiJBotDUQFRsrprpGvKwdIQXBmNsRNJQGZ9HE3cHao1ujsSERFWe7SwTmS3Hi8qKtpYKaCAqOjAOlFUmLhVYvBWCRcSo4OZfUTSNYQATAFvN7MLirZ3MbHMGRXrRN9cHelEVen+c9wbWSFRnZxfqZEIibpVojwvIJKy0Nf837SACH3SICCsZL3lnE0sDpYIiVhUZAIyk/iJulUik8xKVRpyT6SDMZtaJ9ICITuxtxMVBa0U4JaK5cQofO/kEVU7XOnUJZIeTvjlvxuwETjWzH6VU+8A4OPAgwgP6jvM7L/ia0cA/0qYy68GjjGzzfG1nYFPAI8n2CH/28z+sdWYYvFQWECkcTHhDC2FXR2lyYXXs1uL57g3OhUSySqObCprWCggIImtiAVEUtYiPkK1NmaLEsHNUbYgHqrBipINyEziJxJxkU5mpVp4v+TknlD9PFgnVEkLiIxAyFopoPMgzSZWCmiyM6kzEoyqiFgCPgWcZmZnSHop8FngsHQFSauA8wgi4TJJZWDn+Np03OapZnadpE8A7wTeETc/HbjczP4qrr/XIN9ML7+THGcgtJpEOp1gWrk3oHMhUZsgdoVAbcKwyMLfEtikYWWjVjZqk4ZNGEwYNlGD+FCzo2T1OtnDIps/JmoQxUJlMilPCZ0odsWkdxdtsndHvaycfBbz+3YAjXt2JJSixjwe2Y29yuWFcSzZvB9tNg9LGOR+H+7iGIyLw4VEMSTtwXzWSYBzCems12Wq/iXwIzO7DMDMKmZ2d3ztecDPzOy6+PxU4mWdkvaL+z8p6cjM7hjAW6njYmIMGPYNwArR6ss/u4No2irRwr2RFhLVifA6iITmQiL8zREQUzVq5UYBYWkBUa4RTVWJJmq5h8qNB5HVjwWioxTuk6wqCaJmofUkrEjJWFBKpD6P9CqW+Y3AwnmjiFjgOupUVGTrZ6+3EBXOaOBCooFpSRtSxwmZ6/sAt5tZBcDMDLgF2DdTb39gu6RvSLpK0pmSdo+v7QvcnKq7HthbUhS3uxX4pKQrJH1H0uP6+xYbcTeHMzQUTrOdkLg6mvWXJKyKrRKhLL6WLP9MZbNs2IWzhZCwEljZcl0YKgUBEMXuDJUMRUapVCOKjMmJYFnZPhOSVVSrYUBK3CSxq0NJn9VG14cl9yG4MWRgceyE1WJXRSnndS21LLRJMKZmk/ehVFBmVA/EtHLU6PqAhUtD8+Ip2i0lLRCgWZQiliu3SgR88s+nT+m0N5vZ2jZ1sj7PPD/nBPAc4InA7cD7gVOAI5v0kW53GPDPZvZ6Sc8Bvi5pXSJgskj6T+DjZvbDNuPOxX8KjAmjZp1IJoWGjZziSSB3wmix4ZOylot2904CJjO5Iuoug5U1aquqsKoKqyowVUUrK5RWVSitrDC5ao6pVXOsnJ5heno7q6e3s8vqrey10/3ss2Yj+6zZyJodtrFiao6JiSqlUo1SKYxRkdUPiMVI+kgES2as9fM2Wiz9RalBBZoNYPtxx1lm3AqsjWMgkCSCteKWTL2bgUvM7LbYenE2IaCSuO66VN11wG1mVovb3WZmlwCY2YXAJNBK4FwK/IekKyW9VtKKTt6Qi4kxYhQERfpXYyFBkfcrNi5Ti1+4eRNpdufOrLXCouDOYEW1LhxKU9W6eFixcpYdVs2wauUsO+2wjV2mt/Kg1ZvYe/X97LN6I/tOb+Qh0/fxiB3vYt2aP7Bqco4VU+Eol6sNwqJUqjUIC+UFbnb49Kaj1KOM9lJl/ppqfRYZRfdZWQTcKuEMA2Z2F3AlcHRc9BJgvZmtz1T9CnCopB3j8+cSVm1AWHVxqKRHxefHA1+KX/8ceEDSgQCSDonLb2sxpk/FGTD/nmANWS/pI5IeUuQ9Dc9T7iwKo7BUNPnCV2rH0IaVHdndRJMVHdXqgmA/VQyb1PzfGg07bcLCbcSzOSSCoDBqU4ZN1SivrDA1FbbuLMfui4moytREhYmoSjmqsao8x2RUZVV5lh3KM0yXZthtYjPTpe2s0BzTpRkmSxVu37Ijm2emiKIStVqNSjWiVkvew7wYqlajlCskBF4q5QIJgZgqZKLNioWezbpZ1wUMNGtlt7iQcIaM44AzJP0T8ABwDICkzwDnm9n5ZnaLpA8BP5JUIYiB1wOY2SZJrwPOiy0c1yZ9mJlJOhb4TGxh2A68JN4htB2/JAidJwKPAn4g6RQz+5dWjVxMjCGjIChgPoairaCYmgyTWZn4b/Dt5+WbaHqvxAKhlDsjvRR00rAVNUqrKuywaoapcoVSZJRLVVbGwqEc1ZiemKGsWl1ArClvY7q0ndXRdnYvP8AepU1MqMqKaI4V0Rw7lGa5besa7t22A9srZcrUoBRERaPpIUzWVYKoMAiBmHHcBNX8ZaXZ1MDZ8+ya+qxLaEGMRD/IEx+OMyQsVjptM7uezFLQuPx1mfMzCdkp8/o4Hzi/ybWfMe8SaYukJwB/CzyNkP/iMDO7LV6C+ivAxYSzkOUkKOqujkoFSqUgJCaj+b+1EDSpmmGIqGpUo+Y5HerCIppfAlqbNLSywoqVs6yanGP11PZcATEVVZguzbC6tJ3p0nZ2Km1l19Imdi9tZvfSHLtFK4ASK7SB1VEQGjuUZ/hdtBu/37Kaak3M1UoZUQHzwqJGbS4KcRPt8lLEqBofmfm7T4FmxcgGWg6AlkuK3SrhOO34NPAx4K/NbHtSaGabJX2gXWMXE2PMKAsKIIiKahV63E0UUu6N9JEsA11Zg1UVplbNsePKGXZduYWdprZRVo2pqMIO5RkmVG2wQqwubWPHaFuDiJjQfHKtR0ysZHV0NztGof6EqkxPzNTdHnOxL6YcWyTmp8kgJOqujhaiqBV5MSPpuImOWQJXhxUUJy4kHKcQbzezb6cLJD3XzC4ws0+1a+xiYswZVUEB1K0UNldBUdQYN1GLE07XahglVK1h5VI9dqKh7yQHg2gIuqxNzAdcTq6aY/XK7aye2s5uK7awy+SWuoCYULVuhVgdbasLhDwRUd7zRgAqd+7HXqUdWK1trNBtrIjmWFPexg6lWW7avAubZlYEK0WOqMgGY5qoZ7y0AtaGbPAlDCDochC4KHAGjY3IszAYPgB8O1P2QQqm13Yx4Yy+oKindi4eN5E18TckfYr/pgMuV62cZXpqhj1WbuZBUw80WCFWRHOFRUT6vHLnfkxHUxw4WWV1NO/2KEdVbtu6E7/fshoIE+gcJWqxW6NUqlGLopDfomYw1946kcRG5H1RNrNI9CVeYokDMd0q4TitibNlPgLYUdLzU5fWAKuK9lNYTEj6DrAnIRJsE/AmM7sqdf0Y4Azgz8zsG0X7dYaDURYUwIK4iTTZFR0N/aWXfiaujTLUVtQaAi6np2bYbcVW9pjaxG4Tm9mtvKnuytgp2sbqaJbVUa2tiEiTCIoJlRrcHsHSMcP0xAw3b9qZmbkyE1GVOUIsRa0WuzpiC4OVLARhzjUm5sqSF1jW17iJdoGV6biJdN0BxVO4kHCcQvwJcCxhI7G3psofAP6haCedWCaONLONAJJeRNhE5OD4fC1hmcuPO+jPGTJGVlAkcROpTb+yQZhpkuWhUQWqU3Gfop5Ku7piYcDlriu2svuKICTWTt7HrqVNDSJir9IODfdoJSLy6jW6PdazU2kra8ohLuP2rTuyaWYFUK27PRSn2g7ZMBMREbYnz5InGFRrvTtiqxwdi0Un+7D4pmCO0x1m9nng85Jea2af7bafwmIiERIxa0gvgofTgLfQZumIM/yMpKBI4iamYotAJm4CWLCio6GvUtgcqzYB1ak4TiIOuNxhapZdV27hQSs2sevEFvaa2Mje5T+we2lLTyIiSzO3x4Sq7FCeaXB7zEYlSiVRi+K9QFusgE2ERPI3qrTPfNlTIGYeS+DqcKuE0w0iP6ZoOSPpoWZ2EyGXxf7Z63nboufRUcyEpDMJa1AhZOJC0t8AvzSzn4SMoI6zONQFRYu4CaAehJkmqhq1idSmX8kxGeIkShM1JicqTE/NsKo8xw7lGXabCK6N1dF2domM3VJColsRkSbt9nhwKWK7bWT71ARzVqJSK1GpRdy7bQfKpRqVSqmeYpu5YFkRzV0c4XNoLxJ6Dj7rxNXRA0VXcjiO05aPA0cA38y5ZsDDinTSkZgws1dBPT7io5LeCPw1wefSknjXtPrOaWvWrOnk1s4iMirWCUhlyyyXsZnZYHOIs2CqEkGU7HkREVGlRim4CMoimjPK262+6VdtAqLtUdhCfCJi+8wEm8tTrCrPcd/sDiF/RLSdu6NpVuh+VmiG6Sj4SSp37tezoKjcuV/99d21KndXd+Luyo7cMzcNQDman6SjvPTaNHFpVNOvre7iUC1YIVSz8DcvYVW/3R3tREAf7udWCccpjpkdEf99aC/9dLWaw8w+L+mTwJOABwO/jq0SewKflfQuM/t0ps1JpPZWX7t27diuvxkFRk1Q1LZtrwdkCuq/jgVQi4gAq83/bA+7c0awpUZpVqgaIYNoDlQtUamKmbmI+6oRc7USWysTVCxiplZmU20FD9RWcndpE7uXtvDgUsR0NFUXA92IiqTtnFW5qTLLb+d249fb92b99l2p1IJVZTKqUi5ViaLWj23dtWFBSETVsFNoaaa1kIhma6hWQ5VaUyGRu3NomkFltnSB4DgDRdK/AqcXdWtkKSQm4k1Gps3s9vj8xcC9wBfN7OxUvUuBf/XVHMuDZGOwURAVVq1i1Wo9fCAdQKhqHLRYSYmJOP20KqJWFaoaUTUimgv7W6gWUamIak08UBOzcyXmaiVmq2U2V6fYPLWCu8s7snv5ATaV72OX0lZ2j0pdiYqk/ubaDL+riOtm13LTzO7cvG3XhnqTpbD3RxGSrJcJpdmFQiKaraFqbV5QtBESXdNOCLjLwhkiBrab7vCzBfiGpHsICyz+08zuL9q4qGViDXCupJWEwMu7gSPiLVGdZc4oWSlq22eCoKhW58MsS1Hd/ZHICYvCKyuFPTBUA22sEc1FqCqiWYjmxFy1RHVObJ0LG3DNzJXZWplgS2WKXSa38IeJHbi7vCP7TNzLfaVN7FLaytrSFBMqFRIVSZ17qlu5uTLFL2f25qaZ3fn9zI5N20xEIQgzjZo8iUEYhS/I0pwtsEYMXEg0o40Fo2Vq7Iz48JUcjtM7ZnYicKKkpxOWir5f0oVm9ldF2hcSE2Z2KwU2DDGzw4v054weoyYoVKoQEbs5YuFAuUzs3IiTWZWxcrITZ7gSPVBFtSiICgPVRGUuuD2218RcbKHYWplg4+RK7pvcIYiKyg7sNbGR3csPFBYVSfkd1S38anZnfjv7INZv342Nc/l5YsqqNcRNJMtD6+e12KqS2ouj4XULIRFtryyOiGhFs/u7i8NxFg0z+x9JiUXiKKB/YsJxYLQERRJHEUHY/Ct1LS+OwqohqVVtMmJyc43STHB7qKoQR1GLqFiZ6lxUd3tsnprigampuqi4Z3KavSbv5+7yprqoeHB5S5zIqtQQYJlww9w2fju3O7+b3YPrt+5Zj4/IYyqqMBlVKTUJvkzTEC9Ri60SHcZH9B13ZzjDji3yBnhDhKQ9gFcSrBIR8DngbUXbu5hwOmLU4iiaBmZm4igsiupWivALP0Iba7GVQvGv+4jKrKhWg9tjZmqCLVOTC0TFg6Ye4J6J1dxd3sRd1QfYu7yR3UvbG3JSJPERv53bk5tm9uDGrXsUek+JZaJcqjGTcz29O2ji4ijNBBERzVQH69boJviyncBwq4TjLBa/As4h7BracQJKFxNOV4yKlaJpYGYmjsLK0QIrhSYjoqoFl0ctCr/058IW4dU5UV1RYutUaaGomF3JxhWruGdiOhYVO7L3xH3cXdrMg0th8ry5MsVv5/bg19se3DI+Io9yqQpMUCrVqKbyZzQuAZ3/qxqxJWIJ4iO6EAOdxEs4jtM39jGzbd02djHhdM2oCArICczMxFFQi5paKaZqSebIKP4LlbmI6oyoTUX5omJ2ih0nZ+ZFRWU1ayfv4+7SJoC28RHNmIyqTETVpnkmIBMvkbg4qrV4BUdl8dwazcizYCx1vIbjjCmSXmZmXwVenZd40sxOLdKPiwmnJ0ZNUCSBmfU4irTbo5WVomJEc8bsXIRqIY6iukJUp9RWVOy2YisbJ1dxz9w0e0/9AYD123djS2Wq4/cwWao0BGFmSQddRlWruzhU65+QaJtrol8UsGr4Sg6n37Tas2aZ8ljgq8ChOdcKfxguJpyeGek4ipTbo52VYqJmcQKoEtGcqM6EjcLaioqZKe6ZWhVERYeWiCzlOCAisU4osvC0V1MrOeIgshA4SiwkqktrkchzT7jLwnGWnHhJKMBbzeye9DVJuxXtx8WE0zdGxUrREEeRdnsUtVLUoLxNzE1HVFZ0Jio2rlzBZJx4qhzVmIyqTMZxFGXVmCqwy9BkgcRViYsjqsaBl0shJDxewnFGie8Q7wTepiwXFxNOXxkVQQE5bo+iVoqtEM2K0pxRWRFRWanComKmUq7HO0yklnmWS9V6dstEZKRfJ4JjthqSZs1UyszNlajNBCtJNCdKM1DaHo7ydqO8zShtr1HaXiHaOrvg/Vu5xa5gLch1cQwqjbbjOANFUhmYBKI4MWUSOLEGKGxKdTHh9J2RdnvAfHAmwUphU2Wi7RWsHFGjTES1bqUoba8RzcVuj1lRmRNRJbgYojmFVR8VUZsTlbmITXNR2O0zFhGlUq0uLqKoRrlUi19bXVyUIotXcEClWuKBbVNs2zyFbStT2hZR3hRR3goTmwl/t9aYeqBGedMc5c2zuUICFoqCIuKikJDw5ZzOCDOGeSbeCZxIiI/Ykip/APi3op24mHAGxqhYKRp2maxW0YopmJkNk+TUJJoh3om0TEQFq0V1K4VVxWTFKG+LqKyMKK0Uc3NREBOzsaioiOpciVpF1OYirBRn3YyMavw6yWSpKN5WPH4NjYKjVhOzWydgc5nS1ojyNjGxCcrbYGKLMflAjYktNSbun6W0efv85F9NWV2akBYKecJigZDIs0a4kHCckcLM/i/wfyX9h5n9Tbf9uJhwBsooCQqrVolWANtpzJpZimBqCs1WoBZEBWXQ7LzrI6SnLlGaiyjNGJWVcTzFLFS3Q2UHqM5FVKcMS1JDJKIivpGVQjClRfPXLTIqUBcYVhXR1hLlzRHlbVDeAhNbgpCY2FJj8v4K5S2zlB7YHvpYYDVInRcUFrm4NcJxlhW9CAlwMeEsAqMiKCAVR7FyBTYzG+IoJidhZiaOo5gAKqgShSWS1OqWipAUqlSPpyjNRlTmQjxFVJmPp7CIWFCoLhzqf0vJuTWcJ8KDqihvVd2lUd4KUw/UmNga3BoTG7cH0QPt4xgKCoumbeplLiScZYKFpdTjhKSLzewZku6mcSmoADOzQul5XUw4i8Jyi6OgFAUrRRRyLiSrPlQxrKxYVITsmdl4inlLROZvg6BQ7nk0F4uIbTC5yZjYXGNyU5Xyprl5t0Y7q0EpZ/+PbJs8ceFCwnGWI0fHfw/ppRMXE86iMipWinZxFACq1uYn3XjVR1pURHMimi1RmitRmRLRXLBUWGouX2iRSF1LylLiI6rEQZZbjKkHqpQ3V/PdGq0m+fS1PGGR7qdIH47jjCxmdkf88g4zmwWQ9DBgf+DbRftxMeEsOqMkKFrGUQCUY+tFbK3QbBy8GMVukEkjmqsSrSxTnjEqU8JKopZ68iwKveaJjFppPr2tlaA0ayE+YmsIsiw/MBPcGq1ERJKfoZzzuBcRFq3aFKVc7ixJlYsVx1lsLpf0dMIy0R8A64EjgDcUaexiwlkSRsntkRtHkbg9sqskyuX5HUlLEZqtYJNlotkq1RVlotnSvHgop4RCyiLRICDqlotQpqrV3RrljVubi4h2GSdbCYum1gqf4J3lj2qLEzMh6eHA54HdgI3AsWb2qyZ1VwBXAFvN7JBU+RHAvxLm8quBY8xsc6bt6cCrgdXZaxnKZrZJ0quAz5vZP0m6puj7cTHhLCmjZKVIx1EA4fVsnMOhibhQuYxmKiGj5vYytRVlLBV/kRYUMG+lALDSwnqqGBP3bw+5I6q1+Qm+0+yQnVorXEg4Tr/5FHCamZ0h6aXAZ4HDmtT9APAj4I+SAknTcZunmtl1kj5ByBnxjlSdP6P4/hrJZkGHA1+MXxfOuuFiwllyRklQJKhcboyrSCbdmSAu6qIjiuZdIaWI0lT8vGYCHBfkdYjyr6tSQ9tmw+SeJyCKZKJM37uotcJxnL4haQ9Cmupnx0XnAp+QtM7M1mfq/inwcOAkUmICeB7wMzO7Lj4/FfgWsZiQtCshGdUzgNcUGNb/SPoVQRccJ2lnoPCvFBcTzlCQuD1guF0fSRyFSgufsbTVwmZmQll6gi6V0MT2xkaxaGiwT2Qn9VLUeL3X1NXNEli1slZ0Q6nkYsQZV6YlbUidn2RmJ6XO9wFuN7MKhPWXkm4B9iXEKgAgaQfg34EXEgRFmn2Bm1Pn64G9JUVmVgNOAd5jZvfnbS2ew5sIYuV3ZjYnqQT8dZGG4GLCGUJGIZ7CcibJbJlKpbZbZKcFSL0sbzJvEseQ1z5r1aiT7XcxREUrQdEqKLMHIaJSKfffx3HaIehXnonNZra2TZ3sjfJm/I8Cp5jZbXGMRbs+QkfSy4BZM/tG+6HGHQVBcw2wp6Sd4uJ7i7bvbqcfx1kE0taKUaTIhJYnNnJ3x2zSV65YqTXZIbTZxN3M0jFku3TmiaxcMUUQFI4zxNwKrI032ULBdLAPcEum3pOBd0taD3wJOEDSL+NrtwDrUnXXAbfFVomnAU+XtD5uC/BLSQc0G5CkYwmBoNcCP4+PnxV9Qy4mnKFmOQiKdqKiI0GRZxGZqzQXFVkWS1C0msxbpfFuIg66wQWFM6yY2V3AlcwnjHoJsD4bL2FmB5rZOjNbB7wcuNbMHhNfvgA4VNKj4vPjCYIDMzvezNam2gI8xsyubTGsfwYeb2a7mtnu8VEo+yW4m8NxFoUQZ9F8crO5yoKJNBEUC36RV6u5k3VeH3VBEWWCLnODLWv5E32z+p3Syq0RRfnip0fc5eF0hMXJ6BaH44AzJP0TYYfOYwAkfQY438zOb9U4Xsb5OuC82MJxbdJHl9ydCubsGBcTztAzKqs92lFEUMDCX+dWqeQLClggKpr1QTo3BnQuKDoh+x67jX/IaadyeYHVRhPllrEpLiicYcTMridnKaiZva5J/UvJpLyOBUdL0RHXKxKB+V+S/pawLLQeKW5mWwu0dTeHMxqMursjoa9xFNB5LEWaTlweg4qf6NfKkTaoVHK3h+O05sPAycA9wCZgc/y3EC4mnJFhOQmKbuMoOo2lWEAvgqKf9Gr96BIXFI6Tj5lFqaOU/C3a3sWEM1IsF0EB7a0UzQIrO7FS5PbRraAoYp3ow2Td4KIZwOTvgsJphWrW8zGqSDpI0l/Gr3eStFfRti4mnJFjuQmKJbFSDFJQtKKVW6NZfowMnSwRbdqHCwrHaUDSGwh7hbwvLtoVOLtoexcTjjMEdCMooHMrRQNFBUVRFgReRs2vDQEeR+E4DRwHPJGwsgQz+y1QeGmoiwlnJFlO1omEXtweRa0UXQmKXqwTRWIj2gVhDnjCd0HhpFHFej5GlFkz25YpK/ywu5hwRpblKigGbaUYuIUCiomIHoIwc1OO94ALCsfhbkmPIE7RLemVhEydhXAx4Yw0y1FQwICsFJn2DbQTFJ1aJ3pcDbIgBqLAZN/P7JmOM4b8PXAW8Mg4Bfc74rJCuJhwRp7lLCj6aqVYbEGR16bdtYJBmNCfQMyGtm6dcMYYM7uREDNxCPB84IA4bqIQLiacZcFyFRTQZytFJo5iQdvsJmHtBEW272Z1+5GB0id7Z9AYqFbr+RhV4k3CZoBnA8/rpK2LCWfZcFH1y8tWVHSbORP6YKVoZYHIEyut6CAeo4iVwa0TjtMbki6SdFD8+sGEnUKfA/yrpLcX7cfFhLPsWM6Coq+JrroVFP3OPdEJOZN9vwWF44wZe5vZVfHrvwS+Z2bPI+wb8ldFO3Ex4SxLlquggAFsa96qXVFB0aLPjl0cqbiJXFHg1gPH6SfbU6+fBHwLwMz+gC8NdZzlLSigvajoKI4i025BHEVCM0FRxN3RyoLRYqnnUrg7nHHFUKXW8zFi1CStlbQD8FTge6lrq4p24mLCWdYs5ziKhFaiorDbo12CqyJBZQUsH32joHWiG0HhcRPOmPFB4OfA9cAlZnYDgKQnAeuLduJiwhkLlruggNZBmn1xeySColt3Rydkloh2a50o2tZxxhUz+y/gQOAI4GWpS+uB1xftp7CYkPQdSddIukrSD1LRn6dLuj4u/35S7jjDxrgIilZWigVl3cRRQHF3R7vyXujAguCCwinKGLo5MLPfm9lVZmapstvN7JaifXRimTjSzA40s4OAfwNOj8vPAx4Tl38E+EoHfTrOojIOggKai4o8t0fROAqgu/iJPLrIkNmLdaJo+3pdd3U4TkcUFhNmtjF1ugaoxeXnm1ny7fFj4CGS3H3iDC3jEEeR0EpULCgrEEcBFM9B0StFsmF2OOm7hcJxBkNHk76kMyXdCrwfOCanypuBb8VZtLJtT5C0ITk2b97c3Ygdp0+Mi6CA/HiKTt0ebVd4NLNOFBUcBTbv6tU64TjOYOhITJjZq8xsH+BdwEfT1yQdDRxJ2BM9r+1JZrY2Oaanp7sds+P0DRcUAxYUeem4W7k4yuVGUbFE1okiGUedZYgxn1K+l2MM6codYWafB54maVcASUcBJwLPMrO7+jg+xxk44yYoshNlp3EUHQmKets2IiJLC0HhrgrHGT4KiQlJO8Y5u5PzFwP3AvdJOpLg9nhmJ5GfjjNMjJOggB6sFJ0IivTRDZ1aKFL04upwq4TjdE7RJ24NcK6klYTAy7uBI8zMJJ0N3An8t6Sk/jPM7N6+j9ZxBkgiKJ5VOmqJR7I4WLW6YNWCzVUW/PK3SqVxcq5WoVRqrFurzU/4lUqh+IdClMu5MReaKDcKmlJpcAmyHMdpS6En3sxuBR7f5NpEX0fkOEvMRdUvu6DoQFBA7HoYlKBIiKIGK8gCQZFB5XJ+/EcL3CrhjGKeiGHAl3A6Tg7j5Pbo1eXRUL9fy0ZLqa+mHtwdjuMsDv5kOk4Txk1QFA3MbGAQgiIREs0ERYoGC4onmnKcJcPFhOO0YJwSXEExK8VABUUpan0OHVknPOeE0xm2MHi4m6MAkh4u6YeSbpD0U0n759R5uqSfSPqVpF9I+oBSwYmSjpB0naQbJZ0raTouf7CkC+OtLq6R9BVJu/TtY8rBxYTjFMAFRZeCIm9zsGbkCYd0eRN3h1snnBHlU8BpZvYIwlYUn82p8wfgFWa2P3AIYYvwVwDEwuGzwIvMbD/gDuCdcbsq8D4ze6SZHQjcDHx4kG/GxYTjFMQFRReCAhoFRTNRkXVpZJNXpa/1Ac9V4SwlkvYADgbOiovOBR4qaV26npldaWa/i19vB64CHhZffh7wMzO7Lj4/lVhoxBt3XZbq6iepdgPBxYTjdMA4uT2axVE0nHcqKKC3OAqYFxTNrBMZ3NXhDCH7ALcn+1rFu3XeAuzbrIGkPYGXAt+Ki/YlWBwS1gN7Z/fGklQC3gh8vV+Dz8PFhON0wbgIClhopVgyQTFgfFmogzG/wV0vB0yn96KSdEKTu6VRTp1wQdqRIAY+YmZXtOgj204Ei8VG4OPtP4DucTHhOF0yToJiqGgnRDxuwll6Nqf3ojKzkzLXbwXWSipDfdLfh2CdaEDSauAC4PxMP7cA61Ln64DbMhttnhz3e1TeBpz9xMWE4/TAOLk90vTFOtEJ3abkdpwhJN7D6krg6LjoJcB6M1ufrhcHWV4AXGhm78t0cwFwqKRHxefHA19KtT0Z2A94sZnN9v1NZHAx4Th9YFxFRU+kBUiz7csdZ/lyHHCcpBuAfwReCyDpM5JeGNd5MyH79IslXRUf7wQws03A64DzJN0I7A18MO7jT4A3EawVP4nbfW2Qb8YjkxynjyzX/T2KpN1ulnI7r24hXFQ4S8Ei/b8zs+uBw3LKX5d6/QHgAy36OB84P6f8clrEYAwCt0w4zgAYVytFob0wigRi9vkL3Vd0OM5gcTHhOANiHFwfrTbaAjrbybNZXEQ38RIehOk4i4rLdccZMMvF9ZHn6sitl3V3JOXduDocZ1ExD/btErdMOM4isVytFIV2GM3Sbc4Jj6NwnKHEfyY4ziKSJyhG3WLRlk4CMau15nt0OI4ztLiYcJwlJiswhllcNHN15ImEZu6OOrVa+x1AOzA5a6LcPobDcZyB4GLCcYaMZWm9SFkncqlU+raJV51SqbMAUMcx3JXWJW5PdJwRYBRWhhS1CrRfAeIBcI4zarhlwnGcjkg2xCri7kgCMVUuz1sJSqX5ekkgZhQ1WieygqLSPCV3U3HSoVVCpZJv9uU4XeKWCcdxuqLZxNt2dUfczuYq83UTgVCpzAuH5HUlVSclJBraO46zpLhlwnGcruk0IBM6sFIkFLVEOE4/6HYzujHHLROO4/REKwtF11aK9NGmv1zcXeE4i4qLCcdxeqZVrEGeCLBKZV5UVKsNoiKvba/WiLwkWp6N03H6h4sJxxkhhnlFh1WrXYmKOhkrRdcCogerRJF04c4yxqzh/1+3xzji0txxnL7Sbg+P5Ms2sQw0i6VYwCK5LnxVh+N0josJx3H6TpFNwfJERT1jpk/mjjNSuJvDcZyBkLg92v3KT5uGG2IplhB3dzhOZ7hlwnGcgdMq0VW9Tmo5aVpQtNzfo8i9uxQn7u4YU/zfvCtcTDiOs2ikJ+dm+SmA3Cya9Xb93sOjBS4oHKcYLiYcx1kSWgmLbET8sIgLx3Hy8ZgJxxkxhnl56KBotdyulRuj5bWCS/g8fsJx2uOS3nGcJafo6o9miaYaVoKkylr15Th5DEMA8CjilgnHcYaCIrEJ3VooivbhOE53uJhwnBFkHF0dCe0ERavlpd0KCXd1OE5r3M3hOM7QUMTdAa1dHq3aOE5L4nTaTue4ZcJxRpTlap0ouhRzsb/03TrhOM1xMeE4zshSVFD4r03HGSwuJhzHGTo6SRTVTii4kHCcwVNYTEj6jqRrJF0l6QeSDorL95B0gaTfSPqFpCcPbLSO4zSwXF0d0JmgaNpHn4WEuzqWP+k9Zbo9iiDp4ZJ+KOkGST+VtH+Teq+N59ffSjpNUjl17QhJ10m6UdK5kqZT154Qz9c3SLpY0l49fzgt6MQycaSZHWhmBwH/Bpwel38Y+LGZPRx4NXB2+s06juMMmjzR4BYJZ8j5FHCamT0C+Ajw2WwFSQ8F3gc8GdgP2BN4bXxtOm7zIjPbD7gDeGd8TcDZwN/H/X8bOGmQb6awmDCzjanTNUAtfn0kcEpc53+B3xPeuOM4Tk906+4YpJBw64TTK5L2AA4GzoqLzgUeKmldpupLga+Z2e/NzIBPAq+Irz0P+JmZXRefn5q6dggwY2aXxuefAl4kaaLf7yWhIwuCpDOBp8Wnz5W0KxCZ2d2pauuBfXPangCckCqqSbqjs+EuKdPA5qUeRAf4eAeLj3fwhDHX2tabZ25QQynEqH3GozbePQd9gxm2/ey7ta/0wx1Qk7QhdX6SmaUtA/sAt5tZBcDMTNIthLlzfarevsDNqfP1zM+vedf2lhRlr5nZJkmbgL2AW7p/W83pSEyY2asAJB0DfBR4JWCZamrS9iRSZhZJG8xsbUejXUJ8vIPFxztYRm28MHpj9vEOlszkPBDM7NBB3yN9u8x57tyZqZetk+2jm/77QlerOczs88xbKJC0e+ryQxiQ8nEcx3GcZcCtwNokvjCOcdiHhXPnLcC61Hl6fs1eWwfcZma17DVJq4HVhLiKgVBITEjaUdKDU+cvBu4F7gO+CrwxLj+UYIq6rP9DdRzHcZzRx8zuAq4Ejo6LXgKsN7P1marnAi+W9KBYcLwB+FJ87QLgUEmPis+PT137ObBC0uHx+XHAeWY2MEdgUTfHGuBcSSsJgZd3A0fEfp63A1+Q9BtgFnhl4gdqw0AjSweAj3ew+HgHy6iNF0ZvzD7ewTJq423HccAZkv4JeAA4BkDSZ4Dzzex8M/udpBOBywk//v+HeNVHHAfxOuC82MJxbdKHmdUkHQ18Mp63b2NeuAwEhQBRx3Ecx3Gc7vAMmI7jOI7j9ISLCcdxHMdxesLFhOM4juM4PTFwMSHpNZKulVSR9LeZa38X7+eR7PlxVOraCklnxG1/Iel8SbsN63jj60+V9L+SfhnnSz9smMcb19ld0u8lnTPosfYyXklHSboyvn6tpDcN83jj6+9SyKf/W0nvG4LxvkDSzyTNSPrXzLVhfN6ajje+vujPW69jjusM0zPX6v/EMD5z7f5PLPoz5wQWYw+NnxNSbr8j59ovgT8xs/sl7QNcIenHZnYzIdJ1GjgwXjXyaeBt8TF041VYOvt54Hlm9mtJK4AVAx5r1+NN1TkV+BZhDfJi0O14NxA+2zslrQF+LukKM7t8GMcr6SmE1LYHAhXgckmXmdmFSzje3xDy+r+Mhf83h/F5azreJXzeoPvPOGGYnrlW4x3GZ67V/4mleuYcFsEyYWZXm9mvYWFSXDO72Mzuj1/fStjXY59UlVXAhMKyl2nCf+5hHe/xwFlxW8xse2Y/k2EbL5L+Ki773qDHmRpTV+M1s8vN7M749f3AdcBDh3W8wFHAGWa2xcxmCBvjvSLbxyKP9wYzu5rwRZvHsD1vrca7JM9bfK+uP+MhfOaajndIn7lWn++SPHNOYGhiJiQ9E9iZoEohbEzyAHAX4eFbA3xiaUa3kJzx7g+slPTd2OT9cUmrlm6EjWTHG/+yOwH4x6UcVzNyPt/0tf2BwwhrroeCnPG2yqk/jAz185bDUD9veQz7M9eKYXzmchi1Z25Z0bOYkPQDSfc0OfZp3wNIOgD4HHCUmW2Li59JyC2+J2Fzko3Au4d4vBPA4QTz2yGEL+P3DPF4Pw28zcz6utHPAMebXFsL/DfwBjO7fcjH2yqn/pKNtwlD+7w1YSDPGwx0zEP7zLXpf+ieuRb0/ZlzitFzzISZ/Wkv7WPF+w3gNWaWTsP9BuBMM9se1zub4L99Ty/3G+B4bwauNLM/xPW+RB/8zQMc72HAZyVBMGmvlHShmT2nl/sNcLzJL7vvAu83s6/2cp+EAY63VU79rul1vC0YyuetBQN53mCgYx7KZ64Vw/jMtWAgz5xTjCV1c0h6NCEQ6fVmdlHm8u+A5ygGOAL4xWKPMU2b8X4ReJqkqfj8ucDVizm+LK3Ga2a7mNk6M1sH/H/At3v9UuuVVuOVtBdwMfAvFjaaW3La/H/4KnCMpB3i/xOvYT5v/jAydM9bG4bueWvHMD5zrRjGZ64No/bMLS/MbKAHIR/4BmAL8If49ePiaxfFZVeljufE13YBzgF+RYia/yqwy7CON77+NuDXhBzp/wmsGebxpvo4Fjhn0GPt8f/Dp+M26WuvHtbxxtffTZikfwd8cAg+38Pj8weATfHrF8bXhvF5azre+PqiP2+9jjnVx7A8c63+TwzjM9fu/8SiP3N+hMP35nAcx3EcpyeGZjWH4ziO4zijiYsJx3Ecx3F6wsWE4ziO4zg94WLCcRzHcZyecDHhOI7jOE5PuJhwHMdxHKcnXEw4juM4jtMTLiYcx3Ecx+mJ/x+pJ/Es4rHWBQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# store NOx sensitivity and remove empty index\n",
"nox_sens = pm_sens.NOx_PM[day_index[0],:,:].squeeze()\n",
"\n",
"# create figure\n",
"plt.figure(figsize=(8, 6), dpi=80)\n",
"\n",
"# specify contour levels for plot\n",
"levels = np.linspace(1E2,2E7,100)\n",
"# plot figure and colorbar\n",
"cf = plt.contourf(X,Y,nox_sens,levels=levels)\n",
"cb = plt.colorbar()\n",
"\n",
"plt.title('PM2.5 sensitivity to NOx on 01/01/2011')\n",
"cb.set_label('Sensitivity of PM2.5 to NOx (ug/m3)/(kg/m2/s)')\n",
"\n",
"# limit domain to California\n",
"plt.xlim(-128,-110)\n",
"plt.ylim(30,45)"
]
},
{
"cell_type": "markdown",
"id": "7d1e26ef",
"metadata": {},
"source": [
"You can see there is a high sensitivity area in the grid cell containing LA. These hotspots indicate areas where emissions contribute more to LA's annual average $PM_{2.5}$ exposure (on a given day).\n",
"\n",
"If we combine the sensitivities here with the emissions discussed previously we can calculate how much emissions on our chosen day contributed to $PM_{2.5}$ in LA.\n",
"\n",
"Let's do that now:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "728dd4dd",
"metadata": {},
"outputs": [],
"source": [
"# combine the sensitivities and emissions\n",
"cont = np.multiply(nox_emis,nox_sens)"
]
},
{
"cell_type": "markdown",
"id": "204b31bf",
"metadata": {},
"source": [
"So how much did emissions from this day contribute to annual average PM2.5 exposure in LA? Let's find out:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "b6132c45",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"On-road emissions of NOx on 2011-07-01 00:00:00 contributed 0.029 ug/m3 to annual average PM2.5 exposure in LA.\n"
]
}
],
"source": [
"# get the total number of contributions by summing over space\n",
"total = cont.sum().data\n",
"\n",
"print(\"On-road emissions of NOx on %s contributed %0.3f ug/m3 to annual average PM2.5 exposure in LA.\" % (date,total))\n"
]
},
{
"cell_type": "markdown",
"id": "6989229f",
"metadata": {},
"source": [
"Naturally, let's extend this calculation to all days of the year in order to estimate how much on-road emissions from California contributed to LA's $PM_{2.5}$ exposure."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "d857a6ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In 2011, on-road emissions of NOx contributed 1.311 ug/m3 to annual average PM2.5 exposure in LA.\n"
]
}
],
"source": [
"# get NOx emissions from all days\n",
"nox_emis = emissions.NOx[:,:,:]\n",
"\n",
"# get PM2.5 sensitivities to NOx for all days\n",
"nox_sens = pm_sens.NOx_PM[:,:,:]\n",
"\n",
"# combine the sensitivities and emissions and sum over space\n",
"cont = np.multiply(nox_emis,nox_sens)\n",
"total = cont.sum().data\n",
"\n",
"print(\"In 2011, on-road emissions of NOx contributed %0.3f ug/m3 to annual average PM2.5 exposure in LA.\" % (total))"
]
},
{
"cell_type": "markdown",
"id": "bee4a885",
"metadata": {},
"source": [
"## Running emission scenarios\n",
"\n",
"In the previous section we looked at how we can estimate the amount different sectors contribute to a pollutant exposure. In this section we introduce a simple extension and demonstrate how to consider the impact of an emission scenario on a specific air pollutant.\n",
"\n",
"Let's generate a test emission scenario where on-road NOx emissions from California are reduced by some amount for every day in the year. Specify some percentage to scale emissions between 0% and 200%."
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "d3db1286",
"metadata": {},
"outputs": [],
"source": [
"# specify emission scenario scaling\n",
"user_scaling = 150\n",
"\n",
"# apply scaling to emissions by converting percentage to relative scale\n",
"nox_emis = emissions.NOx[:,:,:] * user_scaling / 100"
]
},
{
"cell_type": "markdown",
"id": "b7f75dd4",
"metadata": {},
"source": [
"Now, as we did previously, let's combine these emissions with sensitivities to estimate the impact of the proposed emission scenario:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "99fcb0c3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"An increase of on-road emissions of NOx by 150% increases PM2.5 by 1.967 ug/m3.\n"
]
}
],
"source": [
"# combine the sensitivities and emissions and sum over space\n",
"cont = np.multiply(nox_emis,nox_sens)\n",
"total = cont.sum().data\n",
"\n",
"print(\"An increase of on-road emissions of NOx by %d%% increases PM2.5 by %0.3f ug/m3.\" % (user_scaling,total))"
]
},
{
"cell_type": "markdown",
"id": "819fdfc1",
"metadata": {},
"source": [
"We can also consider scenarios where emissions are reduced. Specify some percentage between -100% and 0%:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "bfa7c736",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A decrease of on-road NOx emissions by -50% reduces PM2.5 by -0.656 ug/m3.\n"
]
}
],
"source": [
"# specify emission scenario scaling\n",
"user_scaling = -50\n",
"\n",
"# apply scaling to emissions by converting percentage to relative scale\n",
"nox_emis = emissions.NOx[:,:,:] * user_scaling / 100\n",
"\n",
"# combine the sensitivities and emissions and sum over space\n",
"cont = np.multiply(nox_emis,nox_sens)\n",
"total = cont.sum().data\n",
"\n",
"print(\"A decrease of on-road NOx emissions by %d%% reduces PM2.5 by %0.3f ug/m3.\" % (user_scaling,total))"
]
},
{
"cell_type": "markdown",
"id": "32765f11",
"metadata": {},
"source": [
"In the prior examples we considered scenarios using emissions data from all days of the year. But we can also assess air pollution impacts based on emission scenarios across different time periods. For example, what would happen if on-road emissions in California on Saturdays increased by 15%?"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "b13b5132",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"An increase of on-road NOx emissions on Saturdays by 15% increases PM2.5 by 0.029 ug/m3.\n"
]
}
],
"source": [
"# specify emission scenario scaling\n",
"user_scaling = +15\n",
"\n",
"# get NOx emissions from all saturdays (December 4th 2010 was a saturday, indexing every seven days)\n",
"nox_emis = emissions.NOx[4:-1:7,:,:] * user_scaling / 100\n",
"\n",
"# get PM2.5 sensitivities to NOx for all days\n",
"nox_sens = pm_sens.NOx_PM[4:-1:7,:,:]\n",
"\n",
"# combine the sensitivities and emissions and sum over space\n",
"cont = np.multiply(nox_emis,nox_sens)\n",
"total = cont.sum().data\n",
"\n",
"print(\"An increase of on-road NOx emissions on Saturdays by %d%% increases PM2.5 by %0.3f ug/m3.\" % (user_scaling,total))"
]
},
{
"cell_type": "markdown",
"id": "e4ab70c5",
"metadata": {},
"source": [
"Lastly, we can also consider scenarios from specific locations. Emission scenarios can be run over any subset of grid cells, for the below example I will consider Los Angeles but feel free to choose any grid cell in the domain by changing the lat_location and lon_location variables."
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "eabf52b9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"An increase of on-road NOx emissions from (34.05°,-118.24°) by 25% increases PM2.5 by 0.014 ug/m3.\n"
]
}
],
"source": [
"# choose latitude and longitude coordinate\n",
"lat_location = 34.05\n",
"lon_location = -118.24\n",
"\n",
"# find the grid cell that contains the above point.\n",
"lat_index = np.abs(lat - lat_location).argmin()\n",
"lon_index = np.abs(lon - lon_location).argmin()\n",
"\n",
"# specify emission scenario scaling\n",
"user_scaling = 25\n",
"\n",
"# apply scaling to emissions and isolate specific grid cell\n",
"nox_emis = emissions.NOx[:,lat_index,lon_index] * user_scaling / 100\n",
"\n",
"# get PM2.5 sensitivities to NOx for all days\n",
"nox_sens = pm_sens.NOx_PM[:,lat_index,lon_index]\n",
"\n",
"# combine the sensitivities and emissions and sum over space\n",
"cont = np.multiply(nox_emis,nox_sens)\n",
"total = cont.sum().data\n",
"\n",
"print(\"An increase of on-road NOx emissions from (%0.2f°,%0.2f°) by %d%% increases PM2.5 by %0.3f ug/m3.\" % (lat_location,lon_location,user_scaling,total))"
]
},
{
"cell_type": "markdown",
"id": "b9f72707",
"metadata": {},
"source": [
"This concludes the tutorial on combining adjoint sensitivities with emissions. All of the information given here can be extended to any of the eleven sectors or the six other precursor species of $PM_{2.5}$. Additionally these examples apply to the sensitivities calculated for $O_3$ and $NO_2$ and their respective precursor species. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9ccba8a3",
"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.8.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}