2021-02-C3S-climate-extremes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<img src='./img/logoline_12000.png' align='center' width='100%'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tutorial on C3S climate data access and analysis of climate extremes\n",
    "In this tutorial we will access and analyse data from the Climate Data Store (CDS) of the Copernicus Climate Change Service (C3S). We will focus on climate extremes in the area surrounding the city of Lille in Northern France.\n",
    "\n",
    "The tutorial comprises the following steps:\n",
    "\n",
    "0. [Prepare environment](#prepare): 1) Register with the CDS, 2) obtain a CDS API key, 3) install the CDS API and 4) import Python libraries.\n",
    "1. [Search and download](#search_download) data using the CDS API: We will focus on ERA5 reanalysis data of 2 metre (near-surface) temperature.\n",
    "2. [Read data](#read_data): Once downloaded, we will read and understand the data, including its variables and coordinates.\n",
    "3. [View and plot](#view_plot) maximum temperatures in September 2020.\n",
    "4. [Calculate averages](#calculate_averages) of the maximum daily temperatures in September over the period from 1979 to 2019, and compare these with our findings for 2020."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='./img/climate_extremes.png' align='center' width='100%'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='prepare'></a>Prepare your environment\n",
    "\n",
    "Before we begin we will prepare our environment. This includes registering with the CDS and obtaining a unique CDS API key. We will then install the CDS API and import various preinstalled libraries.\n",
    "\n",
    "#### Register with CDS\n",
    "\n",
    "To register with the Climate Data Store (CDS) visit https://cds.climate.copernicus.eu/ and select **\"Login/Register\"**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='./img/CDS_login.png' align='left' width='60%'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Obtain CDS API key\n",
    "\n",
    "To obtain your unique User ID and CDS API key follow the steps at this link:\n",
    "https://cds.climate.copernicus.eu/api-how-to\n",
    "\n",
    "Please make a note of these as they will be needed to download data from the CDS using the API."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Install CDS API\n",
    "\n",
    "To install the CDS API, run the following command. We use an exclamation mark to pass the command to the shell (not to the Python interpreter)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: cdsapi in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (0.3.1)\n",
      "Requirement already satisfied: tqdm in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from cdsapi) (4.47.0)\n",
      "Requirement already satisfied: requests>=2.5.0 in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from cdsapi) (2.24.0)\n",
      "Requirement already satisfied: idna<3,>=2.5 in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from requests>=2.5.0->cdsapi) (2.10)\n",
      "Requirement already satisfied: chardet<4,>=3.0.2 in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from requests>=2.5.0->cdsapi) (3.0.4)\n",
      "Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from requests>=2.5.0->cdsapi) (1.25.9)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in c:\\users\\cxcs\\anaconda3\\lib\\site-packages (from requests>=2.5.0->cdsapi) (2020.6.20)\n"
     ]
    }
   ],
   "source": [
    "!pip install cdsapi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import libraries\n",
    "\n",
    "We will be working with data in NetCDF format. To best handle this data we need a number of libraries for working with multidimensional arrays, in particular Xarray. We will also need libraries for plotting and viewing data, in particular Matplotlib and Cartopy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CDS API\n",
    "import cdsapi\n",
    "\n",
    "# Libraries for working with multidimensional arrays\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "# Libraries for plotting and visualising data\n",
    "import matplotlib.path as mpath\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import cartopy.feature as cfeature\n",
    "\n",
    "# to suppress SSL warnings\n",
    "import urllib3\n",
    "urllib3.disable_warnings()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To enable us to download data through the CDS API, we need to specify our User ID and CDS API key. Please enter them in the fields below by replacing \"UID\" with your User ID, and API_KEY with your API key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your CDS authentication key\n",
    "URL = \"https://cds.climate.copernicus.eu/api/v2\"\n",
    "KEY = \"UID:API_KEY\"\n",
    "c = cdsapi.Client(url=URL, key=KEY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='search_download'></a>Search and download data\n",
    "\n",
    "Now that we have installed the API we can download data from the CDS directly into our environment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Search for data\n",
    "\n",
    "To search for data, visit the CDS website: http://cds.climate.copernicus.eu To facilitate your search you can use keywords, or apply various filters. \n",
    "\n",
    "The data we are going to use in this exercise is the ERA5 hourly data on single levels from 1979 to present.\n",
    "\n",
    "<img src='./img/CDS.jpg' align='left' width='45%'></img> <img src='./img/CDS_ERA5_hourly.png' align='right' width='45%'></img> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having selected the correct dataset, we now need to specify what product type, variables, temporal and geographic coverage we are interested in. These can all be selected in the **\"Download data\"** tab. In this tab a form appears in which we will select the following parameters to download. We will choose a subset area of 1x1 degrees, corresponding to a region of around 111km North/South and 72km East/West in Belgium and Northern France, around the city of Lille:\n",
    "\n",
    "- Product type: `Reanalysis`\n",
    "- Variable: `2m temperature`\n",
    "- Year: `all`\n",
    "- Month: `September`\n",
    "- Day: `all`\n",
    "- Time: `all`\n",
    "- Geographical area: `North: 51`, `East: 4`, `South: 50`, `West: 3`\n",
    "- Format: `NetCDF`\n",
    "\n",
    "<img src='./img/Notebook3_data.png' align='center' width='100%'></img>\n",
    "\n",
    "At the end of the download form, select **\"Show API request\"**. This will reveal a block of code, which you can simply copy and paste into a cell of your Jupyter Notebook (see cell below) ...\n",
    "\n",
    "#### Download data\n",
    "\n",
    "... having copied the API request into the cell below, running this will retrieve and download the data you requested into your local directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-08-20 11:39:36,129 INFO Welcome to the CDS\n",
      "2021-08-20 11:39:36,132 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/reanalysis-era5-single-levels\n",
      "2021-08-20 11:39:36,309 INFO Request is queued\n",
      "2021-08-20 11:39:37,358 INFO Request is running\n"
     ]
    }
   ],
   "source": [
    "c.retrieve(\n",
    "    'reanalysis-era5-single-levels',\n",
    "    {\n",
    "        'product_type': 'reanalysis',\n",
    "        'format': 'netcdf',\n",
    "        'variable': '2m_temperature',\n",
    "        'year': [\n",
    "            '1979', '1980', '1981',\n",
    "            '1982', '1983', '1984',\n",
    "            '1985', '1986', '1987',\n",
    "            '1988', '1989', '1990',\n",
    "            '1991', '1992', '1993',\n",
    "            '1994', '1995', '1996',\n",
    "            '1997', '1998', '1999',\n",
    "            '2000', '2001', '2002',\n",
    "            '2003', '2004', '2005',\n",
    "            '2006', '2007', '2008',\n",
    "            '2009', '2010', '2011',\n",
    "            '2012', '2013', '2014',\n",
    "            '2015', '2016', '2017',\n",
    "            '2018', '2019', '2020',\n",
    "        ],\n",
    "        'month': '09',\n",
    "        'day': [\n",
    "            '01', '02', '03',\n",
    "            '04', '05', '06',\n",
    "            '07', '08', '09',\n",
    "            '10', '11', '12',\n",
    "            '13', '14', '15',\n",
    "            '16', '17', '18',\n",
    "            '19', '20', '21',\n",
    "            '22', '23', '24',\n",
    "            '25', '26', '27',\n",
    "            '28', '29', '30',\n",
    "        ],\n",
    "        'time': [\n",
    "            '00:00', '01:00', '02:00',\n",
    "            '03:00', '04:00', '05:00',\n",
    "            '06:00', '07:00', '08:00',\n",
    "            '09:00', '10:00', '11:00',\n",
    "            '12:00', '13:00', '14:00',\n",
    "            '15:00', '16:00', '17:00',\n",
    "            '18:00', '19:00', '20:00',\n",
    "            '21:00', '22:00', '23:00',\n",
    "        ],\n",
    "        'area': [\n",
    "            51, 3, 50,\n",
    "            4,\n",
    "        ],\n",
    "    },\n",
    "    'NFrance_hourly_Sep.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='read_data'></a>Read Data\n",
    "\n",
    "Now that we have downloaded the data, we can process it, plot graphs and derive information...\n",
    "\n",
    "We have requested the data in NetCDF format. This is a commonly used format for array-oriented scientific data. \n",
    "\n",
    "To read and process this data we will make use of the Xarray library. Xarray is an open source project and Python package that makes working with labelled multi-dimensional arrays simple and efficient. We will read the data from our NetCDF file into an Xarray **\"dataset\"**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = r'./NFrance_hourly_Sep.nc'\n",
    "# Create Xarray Dataset\n",
    "ds = xr.open_dataset(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can query our newly created Xarray dataset ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt, dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:    (latitude: 5, longitude: 5, time: 30240)\n",
       "Coordinates:\n",
       "  * longitude  (longitude) float32 3.0 3.25 3.5 3.75 4.0\n",
       "  * latitude   (latitude) float32 51.0 50.75 50.5 50.25 50.0\n",
       "  * time       (time) datetime64[ns] 1979-09-01 ... 2020-09-30T23:00:00\n",
       "Data variables:\n",
       "    t2m        (time, latitude, longitude) float32 ...\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    history:      2021-02-09 21:41:33 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-bb24dc24-93b4-4ca3-86ec-b8f3b8b6767d' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-bb24dc24-93b4-4ca3-86ec-b8f3b8b6767d' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>latitude</span>: 5</li><li><span class='xr-has-index'>longitude</span>: 5</li><li><span class='xr-has-index'>time</span>: 30240</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-19eba637-7ab6-46b9-adcc-aa811589abe5' class='xr-section-summary-in' type='checkbox'  checked><label for='section-19eba637-7ab6-46b9-adcc-aa811589abe5' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>longitude</span></div><div class='xr-var-dims'>(longitude)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>3.0 3.25 3.5 3.75 4.0</div><input id='attrs-0dc56b8c-2904-4a5d-976c-b780fe6eef72' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-0dc56b8c-2904-4a5d-976c-b780fe6eef72' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-423d42e8-4cbf-4701-8b5e-83402e56246d' class='xr-var-data-in' type='checkbox'><label for='data-423d42e8-4cbf-4701-8b5e-83402e56246d' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_east</dd><dt><span>long_name :</span></dt><dd>longitude</dd></dl></div><div class='xr-var-data'><pre>array([3.  , 3.25, 3.5 , 3.75, 4.  ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>latitude</span></div><div class='xr-var-dims'>(latitude)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>51.0 50.75 50.5 50.25 50.0</div><input id='attrs-321356f2-782b-456e-9640-c61d26577475' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-321356f2-782b-456e-9640-c61d26577475' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-12319f7a-19dd-496b-ba08-42e583d741cc' class='xr-var-data-in' type='checkbox'><label for='data-12319f7a-19dd-496b-ba08-42e583d741cc' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>long_name :</span></dt><dd>latitude</dd></dl></div><div class='xr-var-data'><pre>array([51.  , 50.75, 50.5 , 50.25, 50.  ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>1979-09-01 ... 2020-09-30T23:00:00</div><input id='attrs-4a0b1367-fbe8-472c-9cb7-7f5e05b9c17e' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-4a0b1367-fbe8-472c-9cb7-7f5e05b9c17e' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f780899a-c674-46a1-8d06-9e7b6e04da83' class='xr-var-data-in' type='checkbox'><label for='data-f780899a-c674-46a1-8d06-9e7b6e04da83' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>time</dd></dl></div><div class='xr-var-data'><pre>array([&#x27;1979-09-01T00:00:00.000000000&#x27;, &#x27;1979-09-01T01:00:00.000000000&#x27;,\n",
       "       &#x27;1979-09-01T02:00:00.000000000&#x27;, ..., &#x27;2020-09-30T21:00:00.000000000&#x27;,\n",
       "       &#x27;2020-09-30T22:00:00.000000000&#x27;, &#x27;2020-09-30T23:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-ced5a947-1842-487e-b07a-21ab6e4f47f1' class='xr-section-summary-in' type='checkbox'  checked><label for='section-ced5a947-1842-487e-b07a-21ab6e4f47f1' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>t2m</span></div><div class='xr-var-dims'>(time, latitude, longitude)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-161bbda9-5f3d-4cd9-9f5f-a4c573f3fa78' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-161bbda9-5f3d-4cd9-9f5f-a4c573f3fa78' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-10924021-959d-4882-98ad-dd59942095c2' class='xr-var-data-in' type='checkbox'><label for='data-10924021-959d-4882-98ad-dd59942095c2' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>K</dd><dt><span>long_name :</span></dt><dd>2 metre temperature</dd></dl></div><div class='xr-var-data'><pre>[756000 values with dtype=float32]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-35a265b8-cf0e-4a5a-8926-ede9f0378147' class='xr-section-summary-in' type='checkbox'  checked><label for='section-35a265b8-cf0e-4a5a-8926-ede9f0378147' class='xr-section-summary' >Attributes: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>Conventions :</span></dt><dd>CF-1.6</dd><dt><span>history :</span></dt><dd>2021-02-09 21:41:33 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data4/adaptor.mars.internal-1612904555.1840663-3849-3-fd23b77f-8e8a-49a0-adf3-a29e452afc16.nc /cache/tmp/fd23b77f-8e8a-49a0-adf3-a29e452afc16-adaptor.mars.internal-1612904555.1846411-3849-2-tmp.grib</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (latitude: 5, longitude: 5, time: 30240)\n",
       "Coordinates:\n",
       "  * longitude  (longitude) float32 3.0 3.25 3.5 3.75 4.0\n",
       "  * latitude   (latitude) float32 51.0 50.75 50.5 50.25 50.0\n",
       "  * time       (time) datetime64[ns] 1979-09-01 ... 2020-09-30T23:00:00\n",
       "Data variables:\n",
       "    t2m        (time, latitude, longitude) float32 ...\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    history:      2021-02-09 21:41:33 GMT by grib_to_netcdf-2.16.0: /opt/ecmw..."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the dataset has one variable called **\"t2m\"**, which stands for \"2 metre temperature\", and three coordinates of **longitude**, **latitude** and **time**.\n",
    "\n",
    "Select the icons to the right of the table above to expand the attributes of the coordinates and data variables. What are the units of the temperature data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While an Xarray **dataset** may contain multiple variables, an Xarray **data array** holds a single variable (which may still be multi-dimensional) and its coordinates. To make the processing of the **t2m** data easier, we convert in into an Xarray data array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = ds['t2m']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's convert the units of the 2m temperature data from Kelvin to degrees Celsius. The formula for this is simple: degrees Celsius = Kelvin - 273.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "t2m_C = da - 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='view_plot'></a>View daily maximum 2m temperature for September 2020\n",
    "We will plot the maximum values of 2m temperature over the subset area of Northern France."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we average over the subset area:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_t2m = t2m_C.mean([\"longitude\", \"latitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we select only the data for 2020:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_2020 = Lille_t2m.sel(time='2020')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the max daily 2m temperature for each day in September 2020:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_2020_max = Lille_2020.groupby('time.day').max('time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the results in a chart:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Max daily t2m for Sep 2020')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFNCAYAAADRi2EuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydeVxVZd7Avz92lM0NREBxQUFREJnUMjVNK7PVlqlmpnnbZtpsasapWZppeqdtaqZ1pnqrmcwxK7PVzCzNNcUUwV1J0MAFBUQQ2XneP+7FELlwgXvuuffwfD8fPsq95zznx/c8HH73Oc95fqKUQqPRaDQajUZjPD5mB6DRaDQajUbTVdCJl0aj0Wg0Go2b0ImXRqPRaDQajZvQiZdGo9FoNBqNm9CJl0aj0Wg0Go2b0ImXRqPRaDQajZvQiZdGo+kQIrJfRC50ctuTIjLI/v83ReSvxkbnGkTkThEptMffy+x4NBqN96MTL43GC7EnPTUi0rvZ61kiokQk3pzIWkYpFaKUyu1MGyIyWUQKmr12s4hsFpEyESkQkb+JiF/noj3dtj/wD2C6Pf5iF7Q5QUS+EZETIlIiIutE5Eedj/as47TqRUR6isiHIlIhIgdE5MYm740TkS/t8R0TkYUiEt3kfRGRp0Sk2P71NxERV/8MGo1V0YmXRuO95AE3NH4jIiOBYPPCMYVuwK+A3sBYYCrwGxe1HQUEATvau6M9OfFp9loYsBh4EegJxAB/Aao7H+pZtOXln0ANtp/xJuBlERlhf68H8H9APDAAKAf+02TfO4ArgRRgFDAT+IUBP4NGY0l04qXReC/zgJ81+f5m4K2mG4jIpSKyxT7ykS8ijzR573oRybUnBIjIJSJyRET6tHQwEfmpfXSkWET+0Oy9c0RkvYiUishhEXlJRAKavK9EZEgLbW4XkcuafO8vIkUiktpsu+7A50A/+22/kyLSTyn1slJqjVKqRil1EJgPnNdkv/0iMkdEttpHd94QkSgR+VxEykXkKxHp0UJcQ4E99m9LRWSF/fVzReRb+4jVtyJybpN9VorIYyKyDjgFDGrW7FAApdQCpVS9UqpSKbVMKbW1SRu3iMguETkuIl+IyIBmDmfbz1mRiDzdPLlrpDUvdpezgIeVUieVUmuBT4Cf2vf9XCm1UClVppQ6BbzU1Cm2fvZ3pVSBve2/Az9vKQ6NRnM2OvHSaLyXDUCYiCSJiC9wPfDfZttUYEvOIoBLgTtF5EoApdS7wHrgBfv8pTeA25RSx5ofSESGAy9j++PcD+gFxDbZpB64H9sIy3hsIyx3OfEzvAX8pMn3M4DDSqmsphsppSqAS4BD9tt+IUqpQy20N5GzR6hmAdOwJT6XYUvgfm+P1QeY3bwRpdReoHEEKEIpNUVEegKfAS9g+/n/AXwmZ879+im2EaFQ4ECzZvcC9SIy157knpHw2c/L74GrgT7AGmBBszauAtKBNOAK4JYWHLREUy9DgXr7z9hIdpOft7V9sW+X7eS+Go2mGTrx0mi8m8ZRr2nAbuBg0zeVUiuVUtuUUg32kZUFwKQmm9wNTAFWAp8qpRY7OM41wGKl1GqlVDXwMNDQ5DiblVIblFJ1Sqn9wKvNjuOI/wIzGkfdsCUu85zY7yxE5H+wJSXPNHvrRaVUoX10Zg2QoZTaYv85PgRGO3mIS4EcpdQ8+8+5AJvzy5ps86ZSaof9/dqmOyulyoAJgAJeA46JyCciEmXf5BfAE0qpXUqpOuBxILXpqBfwlFKqRCn1PfAcTW41O6IFLyHAiWabncCWLDbfdxTwJ2BOk5eb738CCNHzvDQa59CJl0bj3cwDbsR2q+et5m+KyFgR+do+SfoE8EtsIz0AKKVKgYVAMrZbRo7oB+Q32a8COD3ZXESGishi+63KMmxJQ++zmzkT+6jVOmCWiERgG9Wa39Z+zbGPFj0JXKKUKmr2dmGT/1e28H2Ik4fpx9mjWAewzdVqJJ9WsCdVP1dKxWJz3g9bAgW2+VTP22/XlgIlgLTS/gH7/g5x4OUkENZs0zBsc7ma7jsE2+jgfUqpNU3ear5/GHBSKaVai0Wj0djQiZdG48UopQ5gm2Q/A/ighU3exjZ/J04pFQ68gu2POQD2uVS3YBsJe6GVQx0G4prs1w3b7bZGXsY2+pOglArDdsvM2RGQudhuN14LrLePTLVEi3/YReRibCNIlymltjl5zI5wCFty1JT+nDnK6HTyoZTaDbyJLQEDW1L1C6VURJOvYKXUN012i2vy//72mFqkFS97AT8RSWjyWgpNbifaR9m+Av5XKdV8BHKHffsW99VoNK2jEy+Nxvu5FZhiH4VqTihQopSqEpFzsI2OASAiQdhu9f0e+B8gRkQczct6H5gptuUQAoBHOfP6EQqUASdFJBG4sx3xf4RtztJ9tDBq14RCoJeIhDf5GaZgGyGbpZTa2I5jdoQlwFARuVFE/ETkemA4ticV20REEkXk1yISa/8+Dtutwg32TV4Bftf4dKGIhIvItc2amSMiPez73ge86+BYDr3Y+8kHwKMi0l1EzsM2X2yefd8YYAXwT6XUKy00/xbwgIjEiEg/4NfYEkiNRuMEOvHSaLwcpdQ+pdQmB2/fhe0PbDm2uTrvNXnvCaDA/gRcNbZRp782GwlpPMYObPPB3sY2+nUcaLqm1m+wJXXl2EZZWkwIHMRfCSwCBtLyqF3jdruxjczl2m/H9cM21ywcWNLkacfPnT12e7Cv4zUTW6JRDPwWmNnCrU1HlGNb2iFDRCqwJVzb7e2hlPoQeAp4x367dju2W69N+RjYDGRhm+j/hoNjteXlLmxLjxzF5vRO+zkGuA3bE5l/brLvySb7vgp8Cmyzx/iZ/TWNRuMEom/LazQasxGRPwFDlVI/aXPjLoqIKGy3cr8zOxaNRtNxXLLCs0aj0XQU+zINt2JfR0qj0WisjL7VqNFoTENEbsc2qfxzpdRqs+PRaDQao9G3GjUajUaj0WjchB7x0mg0Go1Go3ETOvHSaDQajUajcRNeMbm+d+/eKj4+3uH7NTU1BAQEOHxf03G0W2PRfo1DuzUW7dc4tFtjcYffzZs3Fyml+rT0nlckXvHx8Wza5GiZIqisrCQ4ONiNEXUdtFtj0X6NQ7s1Fu3XOLRbY3GHXxFpXl7sNJa41VhQUND2RpoOod0ai/ZrHNqtsWi/xqHdGovZfi2ReB086Ki0m6azaLfGov0ah3ZrLNqvcWi3xmK2X0skXhqNRqPRaDTegCUSr2HDhpkdgmXRbo1F+zWOYcOGUVVVxTnnnENKSgojRozgz3/+MwBz5swhMTGRUaNGcdVVV1FaWmpytN6H7rvGod0ai9l+LZF4+fr6mh2CZdFujUX7NQ5fX18CAwNZsWIF2dnZZGVlsXTpUjZs2MC0adPYvn07W7duZejQoTzxxBNmh+t16L5rHNqtsZjt1xKJ186dO80OwbJot8ai/RrHzp07ERFCQkIAqK2tpba2FhFh+vTp+PnZHuoeN26c6ZNtvRHdd41DuzUWs/1aIvHSaDQaR9TX15OamkpkZCTTpk1j7NixZ7z/73//m0suucSk6DQaTVfDEolX3759zQ7Bsmi3xqL9GkejW19fX7KysigoKGDjxo1s37799DaPPfYYfn5+3HTTTWaF6bXovmsc2q2xmO3XEonXwIEDzQ7Bsmi3xqL9GkdztxEREUyePJmlS5cCMHfuXBYvXsz8+fMRETNC9Gp03zUO7dZYzPZricRr/fr1ZodgWbRbY9F+jeFgaSXzl67l2LFjp59YrKys5KuvviIxMZGlS5fy1FNP8cknn9CtWzeTo/VOdN81Du3WWMz26xUlgzQajcYZGhoU8zd+zxNLdlFTW09Sv++5647bqK+vp6Ghgeuuu46ZM2cyZMgQqqurmTZtGmCbYP/KK6+YHL1Go+kKWCLx0jWtjEO7NRbt13UUHD/Fg4u2su67YhL7hrL7SDmlQf3YsmXLWdt+9913JkRoLXTfNQ7t1ljM9itKKVMDcIb09HTVWpFsjUbTdVFK8e63+fz1s100KMUfLk3iuvQ40v/6FdOGR/HMtSlmh6jRaLoYIrJZKZXe0nuWmOOVkZFhdgiWRbs1Fu23cxw+UcnP//MtD32wjeSYML741URuGjsAf18fknv58PXuozQ0eP6HS29E913j0G6NxWy/lki8KisrzQ7Bsmi3xqL9toyjUj8LFy5kxIgR+Pj48NS8z5j+7Go25pXwl8tH8PZt44jr+cNE+ZE9obiihuwCXQ7ICHTfNQ7t1ljM9muJxEuj0VgLR6V+kpOTee2tBUQNHc1zX+WQ2DeUz+87n5vPjcfH58wlIZJ7++LrI6zYfdSkn0Kj0WjOxhKJ1/jx480OwbJot8ai/bZMS6V+APZWh3Hn4iOcqKzltvMH8s4d44nv3b3FNqZNOpcxA3rw1S6deBmB7rvGod0ai9l+LZF45eXlmR2CZdFujUX7dUzTUj/nTZrCf3L8uO+dLAb36c7I2HCuSI3B18fxwqd5eXlcmBTJrsNlHCrVt25cje67xqHdGovZfi2ReB05csTsECyLdmss2q9jGkv9vL70W9786EuWrvmW312SyMJfnkuwv2+b+x85coQpiVEA+najAei+axzarbGY7dcSiZdGo7EeJRU13P12Jg8tziUmKZ0bokv4xaTBrY5yNWdwn+4M6NVNJ14ajcZjsETiNXz4cLNDsCzarbFovy3zzurtTHn8M5btOMJ9kwfQu2wvE37UvvW4hg8fjogwJTGSdd8VUVlTb1C0XRPdd41DuzUWs/1aIvGqr9cXVKPQbo1F+z0TpRQPvr+VB95cSd7cB/H58Le88ZsfM336NGbOnMmHH35IbGws69ev59JLL+Wiiy5y2Faj26mJUVTXNfDNviJ3/RhdAt13jUO7NRaz/Voi8dqzZ4/ZIVgW7dZYtN8z2Xesgnc35XPr5ZMpOrCbPTu3s337dv70pz8BcNVVV1FQUEB1dTWFhYV88cUXDttqdHvOwJ50D/Blub7d6FJ03zUO7dZYzPZricRLo9FYgw25xQDcMmEg/r6uuTwF+PkwcWgfVuw6ijeUSNNoNNbGEolXTEyM2SFYFu3WWLTfM8nIKyEqLJD4Xt3a3rgNmrqdkhjJkbIqdhwq63S7Ghu67xqHdmssZvu1ROIVGxtrdgiWRbs1Fu33B5RSbMgtZuzAXog4/+SiI5q6vSAxEhG9rIQr0X3XOLRbYzHbryUSL7MLXloZ7dZYtN8fyCuq4Fh5NeMG9XJJe03d9g4JJDUuQs/zciG67xqHdmssZvu1ROKl0Wi8n4y8EgDGDuppSPtTEyPJzi/lWHm1Ie1rNBqNM1gi8Wqs6aZxPdqtsWi/P7Aht5g+oYEMclB7sb00d9u4iv3Xe/SolyvQfdc4tFtjMduvJRKv9PR0s0OwLNqtsWi/NpRSZOSWMHZgT5fM74Kz3SZFhxIdHsQKXTTbJei+axzarbGY7dewxEtEgkRko4hki8gOEfmL/fWeIvKliOTY/+3R2WOtX7++8wFrWkS7NRbt18aB4lMcKatirIvmd8HZbhtXsV+Tc4zqOr1AZWfRfdc4tFtjMduvkSNe1cAUpVQKkApcLCLjgIeA5UqpBGC5/fvOHahaz9kwCu3WWLRfGxl5tvW7xrtwfldLbqcmRVJRU89G+3wyTcfRfdc4tFtjMduvYYmXsnHS/q2//UsBVwBz7a/PBa40KgaNRuMdbMgtoXdIAIP7GDv34tzBvQny92G5vt2o0WhMwtA5XiLiKyJZwFHgS6VUBhCllDoMYP83srPHmTBhQmeb0DhAuzUW7bdxfpfr1u9qpCW3Qf6+nDe4N8t3F+pV7DuJ7rvGod0ai9l+/YxsXClVD6SKSATwoYgkO7uviNwB3AHQr18/Vq5cCcCgQYMIDQ0lOzsbgF69euHj48OxY8cA8PPzY8KECWRmZlJWZlulOj09ncLCQvLz8wFISEggMDCQ7du3AxAZGcnQoUNZu3YtAIGBgYwfP55NmzZx8qRt0G7s2LEUFBRw8OBBAIYNG4avry87d+4EoG/fvgwcOPD0vePg4GDGjh1LRkYGlZWVAIwfP568vDyOHDkC2Cqk19fXn64bFRMTQ2xs7Ok1RkJCQkhPT2f9+vWnh0YnTJjA3r17OXrU9ok9OTmZ6upqcnJyAIiLiyMqKopNmzYBEBYWRlpaGmvXrqWurg6AiRMnsmPHDoqLbbd3UlJSKC8vJzc3F4D4+Hh69uxJZmYmVVVVREdHk5KSwqpVq1BKISJMmjSJ7Oxsjh8/DkBaWholJSXs37/f4XkaMWIEq1ev1uepyXlau3YtQUFBnT5PAD169PDK83TsVAOHTlRxUx//07/nrjhPfn5+jBkz5qzzFB94kuUllbz92ddcf/FEt/4+efN5gjN/n6qqqoiPj/eo3yernKcjR45QXl7ukvMEnnfdM/s8BQUFMWzYMEN/n1pD3PWpT0T+DFQAtwOTlVKHRSQaWKmUGtbavunp6arxRLXEypUrmTx5sivD1djRbo1F+4X3NuXz2/e3suz+iQyNCnVZu47cHjlRxbgnlvPQJYn8ctJglx2vq6H7rnFot8biDr8islkp1eLjk0Y+1djHPtKFiAQDFwK7gU+Am+2b3Qx8bFQMGo3G88nILaFn9wASIt2ztk7f8CBG9AvTy0poNBpTMHKOVzTwtYhsBb7FNsdrMfAkME1EcoBp9u87RXKy03cwNe1EuzUW7de2cOo58a5bv6uR1txOTYxk04ESSk/VuPSYXQndd41DuzUWs/0a+VTjVqXUaKXUKKVUslLqUfvrxUqpqUqpBPu/nX6u2+xHQ62MdmssXd1vfskpDpZWMs6AMkGtuZ2SFEWDglV7j7n8uF2Frt53jUS7NRaz/Vpi5frGiXsa16PdGktX9/tDfUbXLZzaSGtuR8WE0zskQC8r0Qm6et81Eu3WWMz2a4nES6PReCcZucVEdPNnmAsn1TuDj49wwbBIVu45Sl19g1uPrdFoujaWSLzi4uLMDsGyaLfG0tX9bsizze/y8XHt/C5o2+3UpEjKqurYfOC4y4/dFejqfddItFtjMduvJRKvqKgos0OwLNqtsXRlv4dKK8kvqWScAbcZoW23ExL64O8rLN+tbzd2hK7cd41GuzUWs/1aIvFqbY0vTefQbo2lK/ttrM841oCJ9dC225BAP8YN6sXyXYWGHN/qdOW+azTarbGY7dcSiZdGo/E+NuwrISzIj8S+YabFMDUxkn3HKthfVGFaDBqNpmthicQrLMy8C7fV0W6NpSv7zcgr5pyBvfA1YH4XOOd2SqLtlsMKfbux3XTlvms02q2xmO3XEolXWlqa2SFYFu3WWLqq3yMnqthffMqQ9bsaccZt/17dSIgM0YlXB+iqfdcdaLfGYrZfSyRezhSl1HQM7dZYuqrfxvldRk2sB+fdTkmKJCOvmPKqWsNisSJdte+6A+3WWMz2a4nEq7Gqucb1aLfG0lX9bsgtJjTIj6Ro44b8nXU7NTGK2nrF2pwiw2KxIl2177oD7dZYzPZricRLo9F4Fxm5JZwT39Ow+V3tIa1/BOHB/npZCY1G4xYskXhNnDjR7BAsi3ZrLF3R79GyKnKLKgxbRqIRZ936+foweVgfvt59lIYGZWhMVqIr9l13od0ai9l+LZF47dixw+wQLIt2ayxd0e+GxvqMA42b3wXtczslMZLiihqyCkoNjMhadMW+6y60W2Mx268lEq/i4mKzQ7As2q2xdEW/GbnFhAT6MaKfsY90t8ftpKF98PURVuii2U7TFfuuu9BujcVsv5ZIvDQajfewIbeY9Pge+Pl6zuUnolsAYwb00PO8NBqN4XjOla8TpKSkmB2CZdFujaWr+T1WXs2+YxWGLiPRSHvdXpgUya7DZRwqrTQoImvR1fquO9FujcVsv5ZIvMrLy80OwbJot8bS1fyers840NiJ9dB+t3oV+/bR1fquO9FujcVsv5ZIvHJzc80OwbJot8bS1fxm5JbQPcCX5Jhww4/VXreD+3RnQK9uOvFykq7Wd92JdmssZvu1ROKl0Wi8gw25xYyJ74m/B83vakREmJIYybrviqisqTc7HI1GY1E87+rXAeLj480OwbJot8bSlfwWnawm5+hJt9xmhI65nZoYRXVdA9/s06vYt0VX6rvuRrs1FrP9WiLx6tnTPRfyroh2ayxdye9G+/pd7phYDx1ze87AnnQP8NVPNzpBV+q77ka7NRaz/Voi8crMzDQ7BMui3RpLV/KbkVtMsL8vo2KNn98FHXMb4OfDxKF9WLHrKErpVexboyv1XXej3RqL2X4tkXhpNBrPJyOvhPT4Hh45v6spUxIjOVJWxY5DZQDk5+dzwQUXkJSUxIgRI3j++ecBuP7660lNTSU1NZX4+HhSU1PNDFuj0XgJfmYH4Ap69OhhdgiWRbs1lq7it6Siht1Hypk5Ktptx+yo2wsSIxGxLSuRHBOOn58ff//730lLS6O8vJwxY8Ywbdo03n333dP7/PrXvyY83D0jeZ5CV+m7ZqDdGovZfj37o6eTmL0YmpXRbo2lq/h19/wu6Ljb3iGBpMZFnJ7nFR0dTVpaGgChoaEkJSVx8ODB09srpXjvvfe44YYbOh+0F9FV+q4ZaLfGYrZfSyReq1atMjsEy6LdGktX8bsht5ggfx9GxUa47ZidcTs1MZLs/FKOlVef8fr+/fvZsmULY8eOPf3amjVriIqKIiEhocPH80a6St81A+3WWMz2a4nES0+CNQ7t1li6it+MvBLS+vcgwM99l5zOuG1cxf7rPT883Xjy5ElmzZrFc889R1jYDwW+FyxY0OVGu6Dr9F0z0G6NxWy/lki8RMTsECyLdmssXcFv6akadh8pc+ttRuic26ToUKLDg1ixy5Z41dbWMmvWLG666Sauvvrq09vV1dXxwQcfcP3113c6Xm+jK/Rds9BujcVsv5ZIvCZNmmR2CJZFuzWWruB3Y14JSrmnPmNTOuO2cRX7NTnHqKqt49ZbbyUpKYkHHnjgjO2++uorEhMTiY2N7Wy4XkdX6Ltmod0ai9l+LZF4ZWdnmx2CZdFujaUr+M3IKyHQz4eUOPfN74LOu52aFElFTT2vL/ycefPmsWLFitPLRyxZsgSAd955p0veZoSu0XfNQrs1FrP9WmI5iePHj5sdgmXRbo2lK/jdkFvM6P4RBPn7uvW4nXV77uDe+PsKFT0THM4JefPNNzt1DG+mK/Rds9BujcVsv5YY8dJoNJ7Jicpadh52//wuVxDk70tSdBjZ+aVmh6LRaCyEJRKvxjV2NK5HuzUWq/v99vT8LvcnXq5wmxIbwdaCE9Q36KfMmmP1vmsm2q2xmO3XEolXSUmJ2SFYFu3WWKzuNyOvmABfH0b3d+/8LnCN25S4CE5W15F77KQLIrIWVu+7ZqLdGovZfi2ReO3fv9/sECyLdmssVve7IbeEVBPmd4Fr3KbaHwjIstjtRkf1JwFefPFFhg0bxogRI/jtb3/rsA2r910z0W6NxWy/lphcr9FoPI+yqlp2HDrBPRcMMTuUDjOod3dCg/zIyi/l2vQ4s8NxGY7qTxYWFvLxxx+zdetWAgMDOXr0aNuNaTSadmGJxGvQoEFmh2BZtFtjsbLfTftLaFDurc/YFFe49fERUmIjyC6w1ohXdHQ00dG2guVN60++9tprPPTQQwQGBgIQGRnpsA0r912z0W6NxWy/lrjVGBoaanYIlkW7NRYr+83ILcHfVxjdv4cpx3eV25S4cHYfLqeqtt4l7XkaTetP7t27lzVr1jB27FgmTZrEt99+63A/K/dds9FujcVsv5ZIvMxeDM3KaLfGYmW/G/JKSI2LIDjA/fO7wHVuU2IjqGtQ7Dh0wiXteRLN60/W1dVx/PhxNmzYwNNPP811113ncA0zK/dds9FujcVsv5ZIvDQajWdxsrqO7QdPmLKMhKv5YYK9tRKvlupPxsbGcvXVVyMinHPOOfj4+FBUVGRypBqNtbBE4tWrl/df3D0V7dZYrOp30/4S6hsUYwe5tz5jU1zlNjIsiH7hQZZ6slEp1WL9ySuvvJIVK1YAsHfvXmpqaujdu3eLbVi173oC2q2xmO3XEpPrR4wYYXYIlkW7NRar+t2QW4KfjzBmgDnzu8C1blPiIiyxgv0HmQWs2H2UqT2OM2/ePEaOHElqaioAjz/+OLfccgu33HILycnJBAQEMHfuXESkxbas2nc9Ae3WWMz2a4kRr9WrV5sdgmXRbo3Fqn4z8ooZFRtOtwDzPtu50m1qXATfl5yipKLGZW26m60FpTy4aCuLtx7m4Q11vLPxANnZ2WRlZZGVlcWMGTMICAjgv//9L9u3byczM5MpU6Y4bM+qfdcT0G6NxWy/lki8NBqN51BRXcfWghNeWZ/RESn2eV7eOupVVlXLPW9voU9IIJ/eM4HkmDAeXLSNn//nWw6fqDQ7PI2mS2GJxMvPzxJ3TD0S7dZYrOh384Hj9vldJs+jcKHbkTHh+Ih3rmCvlOJ3H2zjYGklL9wwmpGx4bx92zj+cvkINuaVMP3Z1by/ucDh04uOsGLf9RS0W2Mx26+095fNDNLT09WmTZvMDkOj0TjB35bu5tXVuWz983S6B1rnD8jFz60mKiyIubecY3Yo7WJ+xgH+8OF2fnvxMO6afGYVgf1FFcx5P5tv9x/nwqRIHr9qJJFhQSZFqtFYBxHZrJRKb+k9S4x4ZWZmmh2CZdFujcWKfr/evIuy9/9IeurIM+oAPvLII8TExJCamkpqaipLliwxNA5Xu21cwd4bPqw2sutwGX/5dCcTh/bhlxMHn/V+fO/uvHvHeB6eOZw1OUVMe3Y1H2cddOpnXLJkSYv1Ht19nq2IFa8LnoTZfi3xcbSsrMzsECyLdmssVvN7qqaO3UcruP7eP/DSvbPOqAMIcP/99/Ob3/zGLbG42m1KXATvbsrn+5JTDOjV3aVtG0FFdR13v51JRLA//7guBR+flp9O9PERbp0wkMnD+jBnYTb3vZPFkm2HeeyqkfQOCXTYfmVlZYv1HsG959mKWO264GmY7dcSI14ajcYzyDxQigruwazp5wNn1gH0dn5YSNXz53kppczdhBYAACAASURBVPjjR9vZX1TB8z8e3WoC1cjgPiEs/OW5/O6SRL7ec4zpz67ms62HHW7fq1cv0tLSAGudZ43GaCyReKWnt3gbVeMCtFtjsZrfjLxifATS7et3Na0DCPDSSy8xatQobrnlFo4fP25oLK52OzQqhGB/X69IvN7fXMCHWw4ye2oC4wc7/5CDr4/wi0mD+ezeCcT1CObutzO5++3MFpfRaOrXzPNsRax2XfA0zPZrWOIlInEi8rWI7BKRHSJyn/31R0TkoIhk2b9mdPZYhYWFnQ9Y0yLarbFYze+G3GJGxoQTGuR/Vh3AO++8k3379pGVlUV0dDS//vWvDY3F1W79fH0YGRPu8UtK5BSW86ePdzB+UC/unZLQoTYSokJZdOe5zLloGMt2HGH6s6tYuv3IGds0+jX7PFsRq10XPA2z/Ro54lUH/FoplQSMA+4WkeH2955VSqXavzo98zI/P7+zTWgcoN0ai5X8VtbUk51/grGDerVYBzAqKgpfX198fHy4/fbb2bhxo6HxGOE2JS6c7YfKqKlrcHnbrqCypp67386kW4Avz/84FV8H87qcwc/Xh7svGMKn904gKiyIX/53M796Zwulp2yjX/n5+R5xnq2Ila4LnojZfg1LvJRSh5VSmfb/lwO7gBijjqfRaMxly/fHqalvYOzAHi3WATx8+If5Qh9++CHJyclmhNkpUuIiqKlrYM+RcrNDaZG/fLqDvYUnefb6VJctC5HYN4yP7j6P+y8cyuKth5n27GqW7yp0WO/RCudZozEStzzVKCLxwGggAzgPuEdEfgZswjYq1qlJAAkJHRtO17SNdmssVvK7bl8RPgJ1h3a3WAdwwYIFZGVlISLEx8fz6quvGhqPEW5/mGB/nJGx4S5vvzN8nHWQd77N567Jg5k4tI9L2/b39eG+CxOYmhTJbxZmc+vcTaT6fM/HHnCerYiVrgueiNl+DV9AVURCgFXAY0qpD0QkCigCFPC/QLRS6pYW9rsDuAOgX79+Y+bPnw/AoEGDCA0NJTs7G7A9WdO3b1927NgB2FaknTBhApmZmacfGU1PT6ewsPD08GJCQgKBgYFs374dgMjISIYOHcratWsBCAwMZPz48WzatImTJ08CMHbsWAoKCk4/tTNs2DB8fX3ZuXMnAH379mXgwIGsX78egODgYMaOHUtGRgaVlbaSHOPHjycvL48jR2xzJYYPH059fT179uwBICYmhtjYWDIyMgAICQkhPT2d9evXU11dDcCECRPYu3cvR48eBSA5OZnq6mpycnIAiIuLIyoqisYFZ8PCwkhLS2Pt2rXU1dUBMHHiRHbs2EFxcTEAKSkplJeXk5ubC0B8fDw9e/YkMzOTuro6+vTpQ0pKCqtWrUIphYgwadIksrOzT0+cTUtLo6SkhP379zs8TyNGjDhdI0ufJ9t52rBhA35+fp0+TwA9evQw7Tz5BwQw41/fEtVNePrSAR5xniIiIhg2bJhLf59CQ0P5xeclJEUobk0O8JjzdKxKeOSbSuLCfJiT5o+vjxj2+9SA8OTirXyWW0vPYF+euHok/kU5HT5PnnjdM+o8Ofv7VFVVdfp7K173zD5PkZGRREdHG/r36YILLnC4gKqhiZeI+AOLgS+UUv9o4f14YLFSqtWx6LZWrl+5ciWTJ0/uVKyaltFujcUqftfvK+aG1zbw3PWpXDnaM2YUGOX21je/5UDJKb56YJLL2+4IVbX1XP2vbzh0opIls8+nX0SwW477n4+X8/Y+P3KOnuTHP4rjD5cmERrk75ZjWx2rXBc8FXf4NWXlehER4A1gV9OkS0Sim2x2FbDdqBg0Go17WJRZQEigHxeN6Gt2KIaTGhfBvmMnKauqNTsUAB5fsoudh8t45poUtyVdAAPDffn03gncOXkw723K56JnV7M2p8htx9dovBUjn2o8D/gpMKXZ0hF/E5FtIrIVuAC4v7MHioyM7GwTGgdot8ZiBb+naur4fNthZozsS3CAr9nhnMYotylxESgF2wpOGNJ+e1iy7TBvrT/AbRMGcuHwKLceOzIykiB/Xx68OJFFd55LUIAvP3kjgz9+tI2K6jq3xmI1rHBd8GTM9mvkU41rlVKilBrVdOkIpdRPlVIj7a9frpRyvDSykwwdOtQVIWtaQLt1TH5+fou16hp55plnEBGKihyPAljB79LtR6ioqWdWWqzZoZyBUW5TYs1bwb5pnxuaOJxb5zxKSlwEv704EXCuz7mKpn5H9+/Bktnnc/v5A5mf8T0XPbea9fuKnWrHFb9HVsMK1wVPxmy/lli5vnEym8b1aLeO8fPz4+9//zu7du1iw4YN/POf/zw9mTU/P58vv/yS/v37t9qGFfwuyiwgrmcwP4rvaXYoZ2CU2/Bu/gzq3d2UxKuxz2Vv20HCbc9R8u2nzB4dRICfj9N9zlU09xvk78sfLh3Owl+Mx89HuOG1DTzyyQ5O1bQ++uWK3yOrYYXrgidjtl9LJF4ajRlER0c7rFV3//3387e//Q3bVEfrcqi0km/2FXP16FiHRZitSEpcBFn5pRj9VHhzGvvcU0t3s6OolrRRI2ioKAE8p8+lx/fk8/sm8vNz43nzm/3MeH4N3+4vcbi9/j3SdDUskXgFBrZdAFbTMbRb52haq+6TTz4hJiaGlJSUNvfzdr8fbjmIUnjcbUYw1m1qXATHyqs5fKLKsGM44sudhbyxNo/LB/lxcN+udvc5V9Ga3+AAXx65fATv3DGOeqW47tX1/HXxTqpq61tts6O/R1bD268Lno7Zft2ygKrRjB8/3uwQLIt22zZNa9X5+fnx2GOPsWzZMqf29Wa/Sine31zAOQN70r9XN7PDOQsj3abYF1LNzi9165OEB0sr+c3CbIb29GPtK7/vUJ9zFc74HTeoF0vvm8iTn+/m9bV5rNh9lGeuSyGtf4+ztu3M75HV8Obrgjdgtl+nRrxEJF1E7heRp0XkURG5TkQ8ZkJHa2t8aTqHdts6zWvV7du3j7y8PFJSUoiPj6egoIC0tLTTixI2x5v9Zn5fSl5RBdd44GgXGOs2KTqUAF8fsgrcN8+rtr6B2Qu2UFtTQ9XnT/OTn3Ssz7kKZ/12D/Tjf69MZv5tY6mua+Cal7/hyc93nzH61dnfI6vhzdcFb8Bsv62OeInIz4HZQB6wGdgDBAETgAdFZDvwsFLqe4PjbJXG1Xs1rke7PZu6+gaOllcTHR50Vq26kSNHnl61GWyrLG/atInevXu32JY3+12UWUCQvw+XjPTMtbuMdBvo50tSvzCyvndf4vWPL/eyaX8Jg3a8SUpKcof7nKtor9/zhvRm6a/O5/Elu3ll1T6W7yrkmWtTGBUb3unfI6vhzdcFb8Bsv23dauwOnKeUqmzpTRFJBRIAUxMvjcadPPLpDt7ZmM/vxkiLNQlnzJhhcoTGU1Vbz+LsQ1w8om+XXa08NTachZsLqG9Q+Br8YMH3xad4ZdU+zgspYsFniyj20j4XGuTPE1eP5OLkvjz4/laufvkbLu5V2mV/jzRdE8NrNbqCtkoGVVZWEhzsvnkWXQnt9kwOllYy+emvqa1XxEQEs2T2+YR363ji4a1+F289xD1vb2HeredwfoJrCzK7CqPdfpBZwAPvZfPFryYyrG+oYccBePbLvbywIod1D05x65yy1uis3xOVtfx18U4Wbi4gsW8oz1ybQnKMZxUeNwtvvS54C+7w2+mSQSLyNxEJExF/EVkuIkUi8hPXhtlxCgoKzA7Bsmi3Z/Kvr78D4OWb0igsq2LO+9mdWlLAW/0u2lxA37Agzh3subd+jHabGte4kOpxQ4/T0KBYlFnAuYN7eUzSBZ33Gx7sz9PXpvDGzemUVNRw5T/X8dxXe6mtb3BRhN6Lt14XvAWz/Tq7nMR0pVQZMBMoAIYCcwyLqp00rvmicT3a7Q8cKq3kvU35XJsexyUjo3nokkSW7Sxk7jf7O9ymN/o9WlbF6pwirk6LMfwWW2cw2m18r+6EBfmRlW9s6aCN+0soOF7pcUt2uMrv1KQolt0/kctS+vHcVzlc+c917Dpc5pK2vRVvvC54E2b7dTbxaryXMgNYoJRyvBqeRmNRXl65D4C7Jg8G4NYJA5maGMnjS3abXrfPUdmVhx9+mFGjRpGamsr06dM5dOhQp4/1UdZB6hsUs8Z4ViLgbnx85PRCqkayaHMB3QN8uTjZMx9icAUR3QJ49vpUXv3pGArLqrj8pbW8tCKHOidHv9zZ/zWazuJs4vWpiOwG0oHlItIHcP/KgQ4YNmyY2SFYFu3WxpETVbz7bT7XjIkltodtzSoR4ZlrU+gVEsA9CzIpr6ptd7uu8uuo7MqcOXPYunUrWVlZzJw5k0cffbRTx1FKsWjzQVLjIhjcJ8QlsRuFO/pualwEewvL2yyL01FO1dSxZNthZoyMpluAZy27aITfi0b0Zdn9k7g4OZpnlu3l6pe/IaewvM393NX/3YW+7hqL2X6dSryUUg8B44F0pVQtUAFcYWRg7cHX19fsECyLdmvjlVX7aFCKuyYPOeP1Ht0DeOGG0RQcr+R3H2xr93wvV/l1VHYlLCzs9DYVFRWdLr2y41AZewrLvWK0yx19NyU2gvoGxY5Dxtwa+2KHvQC5B/o2ym/P7gG8eMNo/nVTGgXHK7n0hbW8smof9Q2Of7fc1f/dhb7uGovZfp2dXP8z4DLgJvv/rwGmGxlYe2gsqKpxPdotFJZV8fbG75mVFktcz7NXaP9RfE8emDaUxVsPs2BjfrvaNsJv07IrAH/4wx+Ii4tj/vz5nf7E//7mAgJ8fbhsVLQrQjUUd/TdxhXsjVrPa9Hmg8T2COYcDytADsb7nTEymmX3T2RKYiRPfr6ba175hn3H2l5/ycj+7y70dddYzPbr7K3GHzX5Oh94BLjcoJg0Go+i8dP23RcMcbjNnZMGc35Cb/7y6Q5TJwY3LbvS+Gn/scceIz8/n5tuuomXXnqpw23X1DXwSfYhLhweSUS3AFeF7NX0CQ0kJiLYkBXsD5VWsm5fEVenda0C5E3pHRLIyz9J44UbRpNXVMGM59fw+ppch6NfRvZ/jcZVOHur8d4mX7cDowGPufL27WvdSadm09XdHi2r4u2M77lqdEyr9Qh9fIRnr08lLNifu9/OpKLauTk/rvTbvOxKc2688UYWLVrU4fZX7jlKSUUN13jgba+WcFffTY2LINuACfY/FCCPcXnbrsBdfkWEy1P6sez+iZyf0Ie/fraLH//fevYXVZyxndH935109euu0Zjt19kRr+acwrZivUcwcOBAs0OwLF3d7aurc6lrUNzTymhXI71DAnn+x6nkFVXw8MfbnWrfVX6VUmeVXQHIyck5/f9PPvmExMTEDh9jUWYBvUMCmeihC6Y2x119NzUugoLjlRSdrHZZm7aHGAo4J74nA3p1d1m7rsTd14bI0CBe+9kY/nFdCnuOlHPx86t5c10eDQ3KLf3fnXT1667RmO3X2Tlen4rIJ/avxdhqNn5sbGjOs379erNDsCxd2e2x8mrmZxzgitR+xPd27o/fuYN7M3tKAh9kHuT9zW0v0tdZv7uPlPFp9iHWrVvHvHnzWLFiBampqaSmprJkyRIeeughkpOTGTVqFMuWLTv9mH17KamoYcXuo1yZ2g8/345+XnMv7uq7jfO8XDnqtSW/lNyiCmaN8czRLjDn2iAiXJ0Wy7L7JzF+UC8e+XQnN7y2gc++/NrQ/u9uuvJ11x2Y7dfZ55OfafL/OuCAUkovrauxNK+tyaWmroF7p7RvcHf21AQy8op5+KPtpMaFMyTS9eVkausbeHnlPl5YnkNdg+K561NbfKLSVfXuPsk6SG29XrurJZJjwvD1EbLyS5maFOWSNhdtthUgnzHS8x9iMIO+4UH8++c/YuHmAv740Xbie8UY2v81Glfi7ByvVU2+1nla0qVrWhlHV3VbdLKaeesPcEVqDAOdHO1qxNdHeP7Ho+kW4Mvd87dQWVPvcNuO+N1zpJyr/rWOf3y5lxkjo0kf0IM/fLiNXCee+OooizIPMjw6jKTosLY39hDc1Xe7BfgxNCrUZQupVtXW82n2IS7y8ALkZl8bRITr0uO48Zz+LMosIL/klKnxuBKz3Vods/22mniJyFr7v+UiUtbkq1xEPKamQ+NjwxrX01XdvrYml+q6eu6Z0vbcrpaICgviH9ensqewnEcX73C4XXv81tU38M+vv+OyF9dyuLSKV+xPe71442gC/Hy4++0tVNU6TvI6yt7CcrYdPOF1o13u7LupceFk55d2qm5nI8t3HaWsqs7jSgQ1x1OuDb+cNBgfEf618juzQ3EZnuLWqpjtt9XESyk1wf5vqFIqrMlXqFLKYz76ZmRkmB2CZemKbksqapi3/gCXpfTr1Orsk4b24a7Jg1mwMZ+Ps1quDeas3++OljPrlfU8/cUeLhweybL7J3Jxsu02VHR4MH+/LoVdh8t47LNdHY7XEYs2F+DnI1yR2s/lbRuJO/tuSmwEZVV15DV70q4jLMq0FSA/b4jnFiAHz7k29A0P4vofxfH+5gIKjltj1MtT3FoVs/06O7l+sIgE2v8/WURmi0iEsaE5T2VlpdkhuBxHtcfmzJlDYmIio0aN4qqrrqK01Ng6cVZ02xavrcmlsraeezs42tWUB6YNJX1AD37/wbYW/ygfOHCgzfM8+vzpXPS3pXxfXMFLN47mXzeNoVdI4BntTEmM4vbzBzJvwwGWbDvc6bgbqatv4MMtB5k8LJLezY7p6biz76b2t0+w7+R6XkfLq1i19xhXjvbsAuTgWdeGO+31UxvrqXo7nuTWipjt19nHkxYB9SIyBHgDGAi8bVhUGoe1x6ZNm8b27dvZunUrQ4cO5YknnjA7VEtxvKKGt77Zz6Ujo10yKd7P14cXbhiNv58Pd8/PPOtWoK+vr8Pz/MnXGxjyy5fJrQ0jZNdilt0/iZmjHI86zbkokdS4CB58fyvfF7vmk//a74o4Wl7NNR78dJ0nkBAZSrcA306vYP/xlkPUNyjtu530iwjmuvQ43tuUz6FSnbRoPBtnE68GpVQdcBXwnFLqfsBjHrcZP3682SG4HEe1x6ZPn46fn+1h1HHjxlFQYOxzDlZ02xqvr83lVG09s6e6bpm6fhHBPHNNCjsPl/HEkjNvBc6cOfOs85yfX0BBtwQu/9d69h2r4M5rLyY5oo4+oa2POAX4+fDiDaNB4J4FmdTUNXQ69vc3FxDRzZ8LEiM73Za7cWff9fURRsaEk1VwosNtKKVYlFlASlyEIU/CuhpPuzbcZV9rzwqjXp7m1mqY7dfZxKtWRG4AbgYW21/zmMdt8vLyzA7BUJrXHmvk3//+N5dccomhx7a626aUnqph7jcHmJEczdAo1/7hu3B4FLdNGMjc9QdYuv2HW4FN/e7fv59NmzN5ZacPjy7eybmDe7Ps/onsWPmx0+c5rmc3nr4mha0FJ3hq6e5OxXyispZlOwu5PKUfgX7eV7TX3X03NS6CXYfKqK7r2AMOOw6VsftIOdd46Er1zfG0a0NMRDDXjInj3W/zOXzCu0e9PM2t1TDbr7OJ1/8A44HHlFJ5IjIQ+K9xYbWPI0eOmB2CYbRUewxs9cf8/Py46aabDD2+ld02599r8zhZXce9Uzs/t6slfntxIimx4cx5f+vpR98b/ZaVlTP54suQ8T/nu9J6nr5mFG/cnM7rL/693ef54uS+/PzceN5Ym8eXOws7HO9nWw9TU9fg8U/XOcLdfTclLoKa+gZ2HS7v0P6LMu0FyFO84yEGT7w23DV5MA1K8YqXj3p5olsrYbZfZ9fx2qmUmq2UWmD/Pk8p9aSxoWkc1R6bO3cuixcvZv78+Yh49gRcb+HEqVr+s24/lyT3JbGvMQ/sBvj58NKNttuK9yzYcvpWYG7hCYaNn0ZZv7FMueQylj0wkWvT43jrrbc6fJ5/NyOR5JgwfrMwm4MdnPOyKLOAIZEhjIoN79D+XY3UTqxgX1vfwCdZh5iapAuQd4a4nt24ZkwsC77Np7CsyuxwNJoWaWsdr20isrWFr20istVdQbbF8OHDzQ7B5TiqPbZ06VKeeuopPvnkE7p1c1y02VVY0W1L/HtdHuXVde1epb69xPXsxt9mjSI7v5Snlu5md20v0qbNoiY0mn89+Sfm/s+PiA4P7vR5DvTz5aUb0qhvUMxesIXa+vbN98orqmDzgePMSov12uTe3X03OjyIPqGBHUq8Vu45RnFFjVeNLnrqteGuyUOob1C8ssp7R7081a1VMNtvWyWDZrolik5SX+/6RSPNYmtBKfPWH2Bc92PMmzePkSNHkpqaCsDjjz/O7Nmzqa6uZtq0aYBtgv0rr7xiWDxWcuuIE5W1/HtdHtOHRzG8n/HL010yMpqfjR/AG2vzqCrYwYlty0kcPoIn77icJ3HdeY7v3Z3Hrx7J7AVb+MeXe3nwYucLBH+QWYCPwFWjvWO+UUu4u++KCCmxER1awX7R5gJ6hwQwaZh3FCAHz7029O/VjatHx/B2xvfcOWkwkWFBZofUbjzVrVUw22+riZdS6kDj/0VkAJCglPpKRILb2ted7Nmzh+hoj3nIssMopfjTxzvIyi/lswBf5q3fz01j+58x4uDu2mNWcdsab67bT3lVnUufZGyL389IorKmnuDhAfxlXsNZo0quOs+Xp/Rj/b5iXl65j7EDezJ5WNtPJzY0KD7IPMiEhD70Dfe+P1qNmNF3R/eP4KtdhZw4VUt4N+eePzpeUcPy3YX8bHw8/l5SgBw8+9pwz5QhfLDlIK+uzuXhmd43euTJbq2A2X6dXUD1duB94FX7S7HAR0YF1VVZnVNEVn4p901NYMyAHvzxo+389I2NHZ6jo2mbsqpa3liby4VJUSTHuG8uU5C/L09fm8KU/v6G38r782XDSewbygPvZTs172VDbjEHSyuZ5SVP13kSKbHtX0j1k+xDtgLkXnSb0dMZ0Ks7V6bGMD/jAMfKq80OR6M5A2c/Xt0NnAeUASilcgCPWdgnJsb7/0AopXj+q730Cw/i7guG8NYt5/DYVclkfn+ci55dzbvffu+SOnDtxQpuW2Puuv2UVdVxnxtHu5riDr9B/r68dGMalTX1zF6whfqG1vvR+5kFhAb6cdGIvobHZiRm9N2R9gcR2jPPa1FmAUnRYW65ze1KPP3acM+UIdTUNfB/q71vrpenu/V2zPbrbOJVrZSqafxGRPwA92cBDoiN9f5Pimu/KyLz+1LuvGAIAX4+iAg3jR3AF7+ayMiYcB5ctI2f/+dbt69PYwW3jjhZXcfra/OYmhh5+g+mu3GX3yGRIfz1ymQy8kp4fnmOw+0qqutYuv0Il46KJsjf+9buaooZfTc82J/Bfbo7PeKVU1jO1oITXjm66OnXhoG9baNe8zYcoOikd416ebpbb8dsv84mXqtE5PdAsIhMAxYCnxoXVvswu+BlZ7GNduUQHR7Edelndoi4nt2Yf9tYHr1iBBvzSpj+7Gre31zg9OiXo5qPCxcuZMSIEfj4+LBp0yaH+3ua287+PE3bSRs3gZ0v3MrKJ37W4XY6izv9zhoTy6y0WF5ckcM33xWd9X5+fj5jz5tIzj9v573fXWeaE1dhVt9NibNNsHfmd/T9zAJ8fYQrUr0v8fK0a0NL3G0f9Xptda7ZobQLb3DrzZjt19nE6yHgGLAN+AWwBPijUUF1Nb7ZV8ymA8e5c/LgFlcI9/ERfjY+nqW/Op+kvra1mW6bu4mjTszXcVTzMTk5mQ8++ICJEyca8SMZhqt+npp6kLE/5SfPLCJr00av9+Isj14xgkG9u3Pfu1lnzX3x8/Mj5uJfcO5Db3UpJ65mdFwERSdr2pybWd+g+GjLQS4Y1qfNclCajjG4TwiXpfTjrfUHKPayUS+NdWlrHa8+IjJcKdWglHpNKXWtUuoaYD3Q2z0htk1ISIjZIXSYxtGuqLBArkuPa3XbAb26884d43h45nDWflfEtGdX83HWwVY/WTuq+ZiUlMSwYcPajM/T3Hb252lk2f4aqiPimT01oVPtdBZ3++0e6Mc/b0qjrLKWB97LoqHJfK/6oHD21PVmVlosYWFhpjlxFWb13ZTTC6m2Xrdx7XdFFJZVe+2kek+7Njji3ilDqKqr5/W13lOGx1vceitm+21rxOtFoKWFZWKA510fTsdIT083O4QOsyG3hI37S7hz0mCn5tT4+Ai3ThjI5/edz+A+3bnvnSx++d/NTj2546jmY2t4stuO/Dxgm8f02ppcJg7tw+j+PTrcjisww29i3zAeuXwEa3KKeLnJIpMfZB4EbGt3menEVZjVdxP7hhHg50NW/vFWt1u0uYDwYH+mJHnMc0rtwpOvDU0ZEhnKzFH9eOub/RyvqGl7Bw/AW9x6K2b7bSvxGqmUWtX8RaXUF8AoY0JqP+vXrzc7hA7z/PK9RIYG8uNz+rdrv0F9Qlj4y3P5/YxEvt5zjOnPrmLx1kMOt3dU87EtPNVtR38egP9uOEBJRQ33TU3oVDuuwCy/P/5RHJel9OPvy/awMa8EpRQfZBYwblBPegQ0mOrEVZjlNsDPhxH9wlod8SqrquWLHUe8tgA5eO61oSVmTxnCqdp6Xl/rHXO9vMmtN2K237YSr9ZWAHRudUA3UF3tnffuM3KL2ZBbwi+dHO1qjq+PcMfEwSyZPYH+Pbtxz9tbuPvtTEqafapzVPPRGTzRbWd+nlM1dfzf6lzOT+jNqH4hHW7HVZjlV0R4/Kpk+vfsxuwFW/hq11H2F5/iilFRpjtxFWb23ZTYCLYdPEGdg1JNn209THVdA7PGeOdtRvDMa4MjEqJCmTEymrnfHKD0lOePenmTW2/EbL9tJV45InLWEtoicgngHR8dPJjnl+fQOySQG8e2b7SrOUMiQ1l057nMuWgYy3YcYfqzq1i63VZ93VHNR2+lsz/P/A3fU1xRw+wpQyzlpSOEBvnzDbRUbwAAIABJREFU0o1plFTUcPf8TIL8fPj0xT91aSeuIjUugsraevYWnmzx/UWbCxjcpzspugC525g9JYGT1XW84UVzvTTWRFqbmC0iQ4HFwDfAZvvL6cB4YKZSaq/hEQLp6emqtcfY6+rq8PPzmApGTvHt/hKufWU9f7w0idvOH+SydncfKeM3C7PZfrCM+y8cypigQs4//3xGjhyJj48tz3788ceprq7m3nvv5dixY0RERJCamsoXX3xxVnue4vat9ft58vPdlB/YTsFbcwiIjKfxc0PvC25G1dVybNnL1J86gU9QCIFRg4i54a9ntVNb38C5g3txZ2Jdp7y4Ck/wO/eb/fz5kx2MDT7Ke4/cYroTV2Gm2/1FFUx+ZiVPXD2SG5pNI2h877cXD+OuyUNMic8VeELfbS93zd/Mmr1FrH1witMlnczAG916E+7wKyKblVItTiZrNfGy7xwI3Agk21/aAbytlGp7LQMX0VbitXPnTtOrjbeXn7yewe4jZaz57RSCA1w7x6O2voHfvr+Vj7IOMu+WsUxI6PgDqJ7gdvOBEq57dQNjBvQgrX+PTrXlI3DNmFgG9fGMp4Y8wa9tftdBzh/am8hQ763N2Bwz3SqlSH30Sy4e0ZenrjlzOuw/lu3hpa+/45uHpnp1LUxP6LvtZdfhMi55fg2zpybwwLShZofjEG906024w29riVerKZ+IiFKqGvhPG9uYuor90aNHvaqTbj5Qwtrvivj9jESXJ10A/r4+PH7VSLYfPMGv3s1iyX0TOvwH1Wy3padqmL0gi34RQbx+czphQZ77KbUjmO0XbPO9vHmukSPMdCsipMRFnLWCfUODYlHmQc4b0turky7wjL7bXpKiw7hoRBT/WZfHrRMGEh7smdcTb3TrTZjtt605Xl+LyL0icsZYuYgEiMgUEZkL3GxceNbk+eXf0at7AD8ZN8CwYwQH+PLPm9I4WV3L/e9mtVmfzxNRSvGbhVs5Wl7FSzekWS7p0lib1LgI9haWU1Fdd/q1jLwSDpZWco0FE11vYfbUBMqr6nhz3X6zQ9F0UdpKvC4G6oEFInJIRHaKSB6QA9wAPKuUetPgGNskOTm57Y08hC3fH2f13mPcPnEQ3QKMvcc8NCqURy9PZt13xfzr6+861IaZbv+zbj9f7SrkoUuSTi9KaTW8qe96G2a7TY0Lp0HBtoM/LCuxKLOAkEA/pg/37gLkYL7fjjKiXzjThkfxxtpcyqpqzQ6nRbzVrbdgtt9WEy+lVJVS6l9KqfOAAcBUYLRSaoBS6nalVJZbomwDsx8NbQ/PL8+hZ/cAfmrgaFdTrk2P5crUfjz71V4ycovbvb9ZbrcWlPLE57u4MCmKW86LNyUGd+BNfdfbMNttSmzjCva22417vsvj5Tk/4/Abd5I+etTpWpglJSVMmzaNhIQEpk2bxvHjrS+86imY7bcz3Dc1gbKqOuY2G/VyVAvW3efIm916A2b7dbZWI0qpWqXUYaVUadtbu5ecnByzQ3CKrPxSVu45xm3nD6R7oHueWBER/nrVSOJ7dWf2O1vaXa/MDLdlVbXc8/YW+oQE8sy1oxARt8fgLryl73ojZrvtFRJIXM9gsuyJ15p9JYRNvoXP12w6o87ok08+ydSpU8nJyWHq1Kk8+eSTpsbtLGb77QzJMeFcmBTJ62vzKG8y6uWoFqy7z5E3u/UGzPbrdOKl6TwvLM8hops/Pxsf79bjhgT68eKNozl+qpZfL8w+oz6fp6GU4neLtnGwtJIXbxxNRLcAs0PSaDpMSmzE6RGvFd/XkjB8FD+K73FGfdCPP/6Ym2+2TZW9+eab+eijj8wMuctw39ShnKis5a31B06/5qgWrD5HGldiicQrLq714tKewNaCUlbsPsptEwYS4qbRrqaM6BfOwzOHs3LPMV5b4/zat+52Oz/jez7bdpg5Fw1jzICebj22GXhD3/VWPMFtalwEh05UkZVfyvrcYmalxSIiZ9TCLCwsJDo6GrD94T969KjJUTuHJ/jtDCNjw5mSGMlra3I52eQBiEbMPEfe7tbTMdtvuxIvEQkTkZ6NX21sGyciX4vILhHZISL32V/vKSJfikiO/d/OLcwEREVFdbYJw3lh+XeEB/tz87nxpsXwk7H9mTGyL3/7Yg+bDzg3R8GdbnceKuPRxTuZNLQPd7hwUVlPxhv6rrfiCW5T7Q+FPPLJDpSCq9NiTK8P6io8wW9nmT01gdJTtby1fv8Zr5t9jqzg1pMx269TiZeI/EJECoGt2Faw3ww4XtHURh3wa6VUEjAOuFtEhgMPAcuVUgnAcvv3naK1xVU9ge0HT/DVrkJunTCQUBOXRBARnpw1in4RQcxesMWpmmXucnuyuo573s6kRzd//nFdCj4+1p3X1RRP77vejCe4HdEvHF8fISu/lLEDe9I31P+sWphRUVEcPnwYgMOHDxMZGWlmyE7jCX47S2pcBJOG9uH1NXmnl/1oqRasu8+RFdx6Mmb7dXbE6zfACKVUvFJqoP2r1SEJ+0T8TPv/y4FdQAxwBTDXvtlc4MqOhe49vLA8h7AgP37uAU/nhQX589INaRwtr2LO+1sxee1bwDav648fbmN/cQXP/3g0vUICzQ5Jo3EJwQG+JPYNBWyjXS3VB7388suZO9d2SZw7dy5XXHGFKbF2Ve67MIGSihr+u+GAw1qw+hxpXImzidc+4FRHDyIi8cBoIAOIUkodBltyBnT6o4MnD9fvPFTGsp2F3DJhoMcsAJoSF8FDlyTx5c5C/tPGIoLucLtwcwEfZR3iVxcOZdygXoYfz5Pw5L7r7XiK2x/F96R7gC89ynOZN28eK1asIDU1ldTUVJYsWcJDDz3El19+SUJCAl9++SUPPdTpmwBuwVP8dpa0/j04P6E3r6/NY9XqNR5xjqzi1lMx22+btRoBRGQ0trJBGcDp9QiUUrOd2DcEWAU8ppT6QERKlVIRTd4/rpQ6a56XiNwB3AHQr1+/MfPnzwdg0KBBhIaGkp2dDUCvXr0YMWIEq/+/vXsPq+q68z/+XtxVRLl4QECDGgSFiAKJkhA0GtOYTC5eYpNm0jEmk6ZtOpnMpGmmnZnUTPNL21zbNNN2OrFN21ymCTFam5hYDSopIREUxQtqxOhBBBEUUASB9fsDJKjcOevsczbf1/PwKPvss/c6n72Qr/uy1ubNQNvjwJmZmRQWFlJbWwtAeno6FRUVHDlyBID4+HgCAwMpLi4GwOFwMHnyZHJzcwEIDAwkIyODrVu3Ul9fD8DMmTNxOp2UlZUBkJCQgK+vL7t37wYgKiqKCRMmkJeXB8CwYcOYOXMmS3++np2VTTw7ezjXZ11NaWkpx44dA2Dq1Km0tLRQUlICQExMDLGxseTn5wMQHBxMeno6eXl5HeOOZGZmsm/fvo6bO5OTk2lsbOx4PHbcuHFERkZ2nEoNCQkhNTWV3NxcmpvbTqVnZWVRXFzME389yo7jLbx6zzTGB2sOHmy76T4uLo6wsDAKCwsBCA0NJSUlhU2bNqG1RinF7NmzKSoq6hjPJjU1lerqag4dOtSv47T6o3y+u76KSaN9ePPBTKqOV1pynPLz82loaAAgIyPDY47Trl27OHGibfy1lJQU6urqLDlOnvLz5I3HaUz0OHRAMFWH9shx8tDjVFjRzM+3NfKLpUkE1x6S4+Shxwm859+96667buCTZAMopT4FcoGdQOv55VrrV7t9U9v7/IG1wAda6+fbl5UAc7TW5UqpsUCO1jqhp+30Nkl2bm4umZmZvX4Od/P0CVlPnmnipp9twc/Xh7X/lNnlGTmT2TY0tXDby7lUn27ivYevtdUEzX3lqX3XDiRbs+yUb1NzK7Oe3kDGxHBevjvV6ubYKltP5I58e5oku6+XGpu11v+itf6t1vrV81+97FQBrwB7zhdd7dbw5fyO/wCs7mMbum9c86WPAnuCX2w8QHCgn8eOvD56eAAvfW0GZScb+LfsnV3e72Uy2x+u2cX+ynpe/OqMIVl0gef2XTuQbM2yU74Bfj7cmhLN+t0VnDpj/TRCdsrWE1mdb18Lr4+UUg8opcb2dTgJ4BrgHmCuUmp7+9dNwI+B+Uqp/cD89u9tZ19FHe8Vl7Ps6jiPHgQ07bIwHr0hgb/sLOe1/MNu2++728r4v61H+Pacy8mMj3DbfoUQoitL0mJpamllzY6jVjdF2FxfLzWWdrFY9/Zko6v0dqmxtbUVHx/PGQt2+fLlvJn9Lq2BIZSX7iN0RABFRUU8+OCD1NfXExcXx2uvvWb5DX7ntbZq7v3dZ+QdPMG737qGqdEhnV5r5f7772ft2rU4HI6O69kD+TzLly9n7dq1jA6PwP+O50mKHsX3rgrioW9/yyNzcQdP67t2ItmaZbd8tdYs+NkWgvx9effb11jaFrtl62ncke+gLzV2GkKi85fHjHC5a9cuq5twgetvXcroRU8werg/oSPaznbdf//9/PjHP2bnzp0sXLiQZ555xuJWfsnHR/H80hRCh/vz0OuFHePZQFu2y5YtY926dRe8ZyCfZ9myZaxe+xfKT54lwM+Hn901nQe/8YDH5uIOntZ37USyNctu+SqlWJway/YjJ/n8eL2lbbFbtp7G6nx7LbzaR6uf1MXyaWaa1H/nn37wFJ80RDJ85GhCO11iLCkpISsrC4D58+eTnZ1tVfO6FB4cyM/unMGhE6f593eLO+73OnHiBFlZWYSFXXhleSCfJysri9e3V9PY3MJzS1MYO2qYx+dimqf1XTuRbM2yY763zYjG10eRXeC0tB12zNaTWJ1vj4WXUmopsBfIbp/258pOL//OZMO81YHKev684yiLU2Px7TT6enJyMmvWrAHgrbfe6ngk1ZPMmhjOw/Mms2pbGW/18g/PQD7PezvLWbWtjNARAcxNjBzwdoQQwgTHyCCy4iNYta2MllbrB5cW9tTbGa/vA2la6+nAvcAflFKL2l/zmDldUlJSrG5Ch19s3E+Qny93Xhl7wfKVK1fy8ssvk5aWRl1dHQEBnnnD/UNzL+fqSeH85+pi9lfUdZttfz/P4RNn+N7bO5gaHUJEp5HpvSUXUzyp79qNZGuWXfNdnBZL+amz5H1u3VkRu2brKazOt7fCy7fTKPOfAtcBP1BK/RPgMf8dqKurs7oJABw8Xs+aoqPck3EZoSMunPYmMTGRDz/8kIKCAu666y4mTbrk6q1H8PVRvPjV6QQH+vHt1ws5Xn2qy/X683mamlv5zhuFKAU/vCXpgordW3IxxVP6rh1JtmbZNd/rp0QSEuRHdqF1lxvtmq2nsDpfv15er1NKTdJafw5tU/wopeYA7wJJphvXVwcPHmT8+PGWtkFrzc827CfAz4d/vHYip0+UX/B6ZWUlDoeD1tZWfvSjH/Hggw9a1NLeOUKCeOGr0/n6yk/53uq93HKVpqrcSc2Zc6zMbXvAtbamipDQCFpbW1n5/x4nad7ijtcuVni4hiLnKX59TxrRI85e8Jo35WKCJ/Rdu5JszbJrvkH+vtySEk12oZMnb0tipAVTvdk1W09hdb69FV7f5KJLilrrOqXUjcBSY63yMsdOneXxd3aQU3Kcb86ZxD89sIycnByqqqqIjY1lxYoV1NfX8/LLLwOwaNEi7r33Xotb3bNr48fwyPWTeX79Pj58YBmNh3fS0lDLAwuuZFTm3ehzDdQV/gWA4ZOvpqw1iY/W7u52e9+YPZHf/egRr89FCGF/i9NieS3/MO/vPMbSK8dZ3RxhM30dx+snWuvv9bbMlN7G8Tp06BBxcXHuaMoFtNa8U1jGD/+8i+YWzeMLErln1mX4+HjM7W+Dtmvf58TGDu5/BsoHj5kg3NNY1XeHAsnWLDvnq7Vm3nObiBgZyJ++keH2/ds5W0/gjnx7GsertzNe580HLi6yFnSxzBIXD3XgDpW1Z/n+qp38dU8lV8aF8sySFOIiRri9HaaNixpDyHApmkyxou8OFZKtWXbOVynF4rRYnvmghCPVZxgXNtyt+7dztp7A6nx7G07im0qpnUCCUmpHp69SYId7mti787OUu4PWmtXby5j/wma27K/iP/5uKm8+kGHLogvcm+1QJPmaI9maZfd8F86IQSksucne7tlazep8e3uq8XXgFtomtr6l01ea1vrvDbfN7ZYvX47D4SA5Oblj2fbt25k1axbTp09nemoai1f8noff3M6kMSN47+FruS9zwgXjdQkhhPB+0aOHcfWkcN4pLKMvt+QI9+jt93R6ejqffvqp27YzED0WXlrrU1rrQ1rru7TWX3T6qjbSmgEKDQ11yXa6mhrnscce44knnuCp379PXdIi1v32Of5tQSJvPXg1k8YEu2S/nsxV2YquSb7mSLZmDYV8F6fGcrj6DJ8dqnHrfodCtgPV0+/p7du38+STT/LYY4/1uI3Q0FCXbGegbDELp6sGQ+tqapzmVs2L7xXx7dcLGe17jnlpCXxj9qQhc5bL6oHm7E7yNUeyNWso5HtjchQjAnx5u8C9M2oMhWwHqqvf00opamtrATh16hTR0dE9biMlJcUl2xkoWxRemzZtMrLdD3Ydo2LqnWz8/fOcWnk/x9b/hl+88KyRfXkqU9mKNpKvOZKtWUMh3+EBftx0xVje23mMhqYWt+13KGTrSi+++CLf/e53GTduHI8++ihPP/10j+t3l29/tzNQtii8XH39vaVV889vbuMbfyjgdNH7PPfcc9RUlvPCCy9w3333uXRfnk7ubTBL8jVHsjVrqOS7OC2W+sZmPth1zG37HCrZusovf/lLXnjhBY4cOdKn39Pd5dvf7QyULQovpVx32e/jA8c5dOIMa3eU88j1k6nZvp7v3Nf2HMEdd9xh7GY7T+XKbMWlJF9zJFuzhkq+V8WFERs6zK1PNw6VbF3l1VdfZdGitmmk+/J7urt8+7udgbJF4TV79uxBb+NUwzkefauI72XvxM9Hsfqha3j4+niio6M7Tktu3LiR+Pj4Qe/Lm7giW9E9ydccydasoZKvj49iUWosuQeqKD/V4JZ9DpVsXaW/v6e7y9ddv+/7NHK91Xobub6oqGhQNyPm7q/i0beKKP7jk/hU7KGhtobIyEhWrFhBQkICDz/8MM3NzQQFBfHf//3fpKWlDXhf3maw2YqeSb7mSLZmDaV8vzhxmtnP5PDYjQl8a87lxvc3lLLtq7/uruCdbU7W/+L7HN1TwNn6kwwLCefKxQ8wemwcuX94Ft3agq9/AFnLHmfMhCndbisrsoXVr7zQMYWdid/3rhi53qPV1AzuUd/TTc2MDPLjbx+uZlrs6EteLygoGNT2vdlgsxU9k3zNkWzNGkr5XhY+givjQskucPLN2ZOMXwocStn2hdaaJ9fu5uSZJuLu+D5xF71eD0z/zq86vj8JnKyo73Z7iQHneOONN7p8zR2/721ReA3WV5KimJvowN/XFldehRBCuNiStFi+l72T7UdOMmO8jLPlTp8dquFw9RmeuyOFxWmxg95eTk7O4Bs1CLaoNFJTUwe9DSm6uuaKbEX3JF9zJFuzhlq+N10xliB/H7fcZD/Usu1NdoGT4QG+3Jgc5ZLtWZ2vLaqN6mqPGkjfViRbsyRfcyRbs4ZaviOD/PlKUhR/LiqnsdnsmF5DLdueNDS18Jed5SxIHsuIQNdcpLM6X1sUXocOHbK6CbYl2Zol+Zoj2Zo1FPNdnBrLqYZzbNhTaXQ/QzHb7ny4+xj1jc0sTotx2TatztcWhZcQQghh2jWXRxAVEkR2gfvG9Brq3i5wEjN6GLMmhFvdFJexReE1ceJEq5tgW5KtWZKvOZKtWUMxX18fxe0zYsjZd5zjdY3G9jMUs+1K+akGcg9UsSg1Bh8Xzo9sdb62KLxGjhxpdRNsS7I1S/I1R7I1a6jmuyQthpZWzertZcb2MVSzvdiqbWVo3XaJ15WsztcWhVdRUZHVTbAtydYsydccydasoZrv5Y6RpIwbzdsGLzcO1Ww701qTXeAk/bJQ4iJGuHTbVudri8JLCCGEcJclqTHsPVbHrqOnrG6KbRU5T/H58dMuGbfL09ii8AoPt89Nd55GsjVL8jVHsjVrKOd7S0o0Ab4+ZBeYudw4lLM9L7vASaCfDzdPG+vybVudry0Kr6SkJKubYFuSrVmSrzmSrVlDOd/RwwOYN8XB6u1lnGtpdfn2h3K2AI3NLawpOsoNSVGEBPm7fPtW52uLwmvz5s1WN8G2JFuzJF9zJFuzhnq+i1NjOXG6iU0lx12+7aGe7cY9lZxqOMfiVNeN3dWZ1fnaovASQggh3Gl2whjCRwS4ZQqhoebtAieRIYFcGz/G6qYYYYvCy89P5vo2RbI1S/I1R7I1a6jn6+/rw23TY/jrngpqTje5dNtDOdvjdY3k7DvO7TNi8HXh2F2dWZ2vLQqvzMxMq5tgW5KtWZKvOZKtWZIvLE6L4VyL5s87jrp0u0M529Xby2hp1Sxx8dhdnVmdry0Kr8LCQqubYFuSrVmSrzmSrVmSLyRFj2LK2JBLphBavnw5DoeD5OTkC5a/9NJLJCQkkJSUxGOPPdbtdj0t28F+nv7ILixjWuwo4iPNDXJqdb62KLxqa2utboJtSbZmSb7mSLZmSb5tFqfGUOQ8xYHKuo5ly5YtY926dRes99FHH7F69Wp27NjBrl27ePTRR7vdpqdlO9jP01e7j9ayp7zW5SPVX8zqfG1ReAkhhBBWuG16271Ib3ca0ysrK4uwsLAL1vvlL3/J448/TmBgIAAOh8Ot7RwMd32e7EIn/r6KW1OiB70tT2aLwis9Pd3qJtiWZGuW5GuOZGuW5NtmzMhA5kwew6ptTlpadbfr7du3jy1btjBz5kxmz57NZ5991u263pBtfz5PX5xraWX19jLmJUYSOiLARa3smtX52qLwqqiosLoJtiXZmiX5miPZmiX5fmlxWiwVtY18fKCq23Wam5upqanhk08+4ZlnnmHp0qVo3XWh5g3Z9ufz9MWmkuNU1Te5ZYogq/O1ReF15MgRq5tgW5KtWZKvOZKtWZLvl+ZNcTBqmH+PE2fHxsayaNEilFJcddVV+Pj4UFXVdaHmDdn25/P0RXahk/ARAcxJMD92l9X52qLwEkIIIawS6OfLLSlj+WDXMWrPnutyndtvv52NGzcCbZfpmpqaiIiIcGcz+62lVdPU3PWUSK78PCfPNLFhTyW3To/G39f+ZYktRmmLj4+3ugm2JdmaJfmaI9maJfleaHFqLH/85DAfFB/jneceIycnh6qqKmJjY1mxYgXLly9n+fLlJCcnExAQwKuvvopSXQ8Q6inZ/vu7xXx8oIoxn/2Kv+VuHvDn6c2fi47S1NJq/GnG86zO1xaF1/mnKoTrSbZmSb7mSLZmSb4Xmj5uNGNHBbFhTyVvvPFGl+v88Y9/7NO2PCHb2rPneKfQSWNzK/E3P8rRP715SWHV18/Tm7cLy0iMGklSdIhLttcbq/O1xTm94uJiq5tgW5KtWZKvOZKtWZLvhZRSzE10sGX/cRqbWwa1LU/I9i87ymlsbmVJWiwb9lbySm6pkf0cqKyj6MhJFqfGDviMWX9Zna8tCi8hhBDCavOmODjd1MKnpdVWN2XQsgucTBozgmeWTOOGqZH8ZN1eth856fL9vF1Qhq+P4rYZ9h67qzNbFF7eNBCdt5FszZJ8zZFszZJ8L3X1pAiC/H3YsKdyUNuxOtvSqtNs/aKGJWnjUErxzJIUHCODeOj1Qk41dP3wwEC0tGpWbXMye/IYHCODXLbd3lidry0Kr8mTJ1vdBNuSbM2SfM2RbM2SfC8V5O/LNZMi2LC3YlBjWlmd7TuFTnwULJwRA8Co4f689LUZHDt1lsezdwzqs3X28YEqKmob3XZT/XlW52uLwis3N9fqJtiWZGuW5GuOZGuW5Nu1uVMcHKlu4EBl/YC3YWW2ra2adwrLuObyCKJGfXkWKnV8KI/dmMD7xcf4wydfuGRf2YVOQoL8mDfFvWegrO67tii8hBBCCE8wN7GtiNiwd3CXG63ySekJyk42sKSLEeTvz5zIdQlj+NHaPRSXnRrUfurOnuODXce4JSWaIH/fQW3L29ii8LL60VA7k2zNknzNkWzNkny7NnbUMJKiQ9g4iPu8rMw2u6CM4EA/bpgadclrPj6K55ZOJ2xEAA+9Xkh9Y/OA9/PeznLOnmt1yxRBF7O67xorvJRSK5VSlUqp4k7LfqiUKlNKbW//uskV+8rIyHDFZkQXJFuzJF9zJFuzJN/uzUt0sPWLak6eaRrQ+63K9nRjM+8Xl3PzFWMZFtD1WaiwEQH8/K4ZHK4+w/ff2Tng+72yC8qYOGYEM8aNHkyTB8TqvmvyjNfvgBu7WP6C1np6+9d7rtjR1q1bXbEZ0QXJ1izJ1xzJ1izJt3tzp0TSqmHTvuMDer9V2a4rPsaZphaWpPd8FuqqCWH8y/zJrCk6yv991v95D784cZpPD1W7deyuzqzuu8YKL631ZsAtg5nU1w/8JkbRM8nWLMnXHMnWLMm3e9NiRhERHDDgYSWsyvbtAieXhQ8n/bLQXtf91pzLuTY+gifW7KLkWF2/9pNdWIZSsCg1ZqBNHRSr+64VUwY9pJT6OrAV+FetdU1XKymlHgAeAIiOjiYnJweAiRMnMnLkSIqKigAIDw9Ha93xup+fH5mZmRQWFlJbWwtAeno6FRUVHTOSx8fHExgY2DF6rcPhYPLkyR1POgQGBpKRkcHWrVs7DtDMmTNxOp2UlZUBkJCQgK+vL7t37wYgKiqKCRMmkJeXB8CwYcOYOXMm+fn5NDQ0AG2nN0tLSzl27BgAU6dOpaWlhZKSEgBiYmKIjY0lPz8fgODgYNLT08nLy6OxsRGAzMxM9u3bR2Vl2w90cnIyjY2N7N+/H4Bx48YRGRnZUdGHhISQmppKbm4uzc1t1+OzsrLYtWsXJ06cACAlJYW6ujoOHjwIQFxcHGFhYRQWFlJfX09RURFCyRUJAAAWDElEQVQpKSls2rQJrTVKKWbPnk1RURE1NW2HLzU1lerqag4dOtTtcUpKSmLz5s1ynDodp/r6enJycgZ9nABCQ0PlOHU6To2NjTQ0NHjUz5OdjlN9fT179+71qJ8nTzpOU0a1smFPOY3nriDv49x+Hafm5uaO32nu+nkaP2U6eQdPsPByfzZt2tSn4/T80ulc/+wG7v3NFp7IGMa8Odf2epzGOBy8/rcDTA3zofzzPYy14Dg1NTVRU1Nj9OepJ8pV43F0uXGl4oC1Wuvk9u8jgSpAA/8FjNVaL+9tO+np6bqnU4MNDQ0MGzbMFU0WF5FszZJ8zZFszZJ8e7auuJwH/1jImw/MYtbE8H6914psX9qwn+fW72PLY9cxLmx4n9/38YEq/v6VfBanxvLsHSm9rv/JwRPc+T+f8MJXU1g4w/031oN78lVKFWit07t6za1PNWqtK7TWLVrrVuA3wFWu2K7T6XTFZkQXJFuzJF9zJFuzJN+eZcaPwd9XsXEAw0q4O1utNe9sK2PmhLB+FV0A11wewXeuu5y3C5xkF/Te7uwCJyMCfPlK0qVPTbqL1X3XrYWXUmpsp28XAi6ZqfL86VXhepKtWZKvOZKtWZJvz4ID/Zg1MZwNeyr6/V53Z1t4uIbSqtMDHtrh4esnM3NCGP+xurjHgWPPNDXz3s5ybp42luEBVtzp1MbqvmtyOIk3gDwgQSnlVErdB/xUKbVTKbUDuA54xNT+hRBCCCvNTXTw+fHTHKo6bXVTevR2QRnD/H256Yqxva/cBV8fxc/unEGQvy8PvV7I2XMtXa53yx13s/e5O8n+/p2XvPbss8+ilKKqqmpAbfAmJp9qvEtrPVZr7a+1jtVav6K1vkdrfYXWeprW+latdbkr9pWQkOCKzYguSLZmSb7mSLZmSb69m5cYCdDvy43uzPbsuRbW7jjKguQoggMHfhYqalQQzy9NYe+xOp5cu7vLdXwS5jDt/h9fMkbYkSNHWL9+PePHjx/w/vvD6r5ri5HrfX2H1nQD7iTZmiX5miPZmiX59m58+HDiHcH9Lrzcme363RXUnW12yQjycxIcPDh7Eq/nH+bPRUcveO3oyQYO+I7n79IvnaD6kUce4ac//anbxvSyuu/aovA6/8iscD3J1izJ1xzJ1izJt2/mTnGQX3qCurPn+vwed2b7doGT6FFBZPTzycvu/OsNk0kdP5p/e2fnBZdYV20rQ2u4MfnCy5lr1qwhJiaGlJTen4h0Fav7ri0KLyGEEMITzUuM5FyLZst+z7t3qaL2LFv2H2dhagw+Pq452+Tv68NLX0vF10fx0BuFNDa3oLUmu8DJVXFhxIR+OYzDmTNneOqpp3jyySddsm9vYYvCKyrKusdS7U6yNUvyNUeyNUvy7ZvU8aMZNcy/X6PYuyvbd7eV0aphUaprx9OKGT2MZ+9Iobislqff28u2Iyc5WHWaJRddzvz8888pLS0lJSWFuLg4nE4nqampHYPtmmJ137XueU4XmjBhgtVNsC3J1izJ1xzJ1izJt2/8fH2YkzCGnJJKWlo1vn04s+SObLXWZBc6mTF+NJPGBLt8+/OnRrL8mgms/LiUT0urCfL3YcEVUZw49uVQDldccUXHKPfQNir91q1biYiIcHl7OrO679rijNf5aRCE60m2Zkm+5ki2Zkm+fTc30cGJ000UOU/2aX13ZFtcVsu+inoWu/hsV2ePL0hkWuwodpfXcmNSFA/c+3UyMjIoKSkhNjaWV155xdi+e2J137XFGS8hhBDCU82ePAZfH8XGPZWkju99Amp3yC50EuDnwy3Too3tI8DPh1/clcqjbxdx/7UTSb7zjR7XPz+Xot3Z4oyXzBdmjmRrluRrjmRrluTbd6OHB5B2WSgb+jishOlsm5pbWb29jPlTIxk13N/ovsaHD+dP38ggOWaU0f30h9V91xaF18yZM61ugm1JtmZJvuZItmZJvv1z/RQHe8prOXqyodd1TWe7cW8lNWfOscTgZUZPZnXftUXhlZ+fb3UTbEuyNUvyNUeyNUvy7Z+5/RjF3nS22YVOxowM5Np4szexeyqr+64tCq+Ght7/ByEGRrI1S/I1R7I1S/Ltn0ljRnBZ+PA+TZptMtsT9Y18tLeS26dH4+drixKg36zuu0MzdSGEEMKNlFLMTXTw8ecnONPUbFk71hQdpblVu2SKIDEwtii8MjIyrG6CbUm2Zkm+5ki2Zkm+/TcvMZKm5lb+duBEj+uZzDa70ElyTAiJUSHG9uHprO67tii8SktLrW6CbUm2Zkm+5ki2Zkm+/XfVhDBGBPj2+nSjqWxLjtVRXFZrdOwub2B137VF4WV6eoGhTLI1S/I1R7I1S/LtvwA/H7Imj2Hj3gq01t2uZyrb7EInfj6KW1PMjd3lDazuu7YovIQQQghvMDfRQUVtI7uO1rp1v80trazaVsZ1iQ7CgwPdum9xIVsUXlOnTrW6CbYl2Zol+Zoj2Zol+Q7MdYkOlOp5WAkT2W7ZX8XxusYhf5kRrO+7tii8WlparG6CbUm2Zkm+5ki2Zkm+AxMRHMj0caN7vM/LRLZvFzoJHe7P3ESHy7ftbazuu7YovEpKSqxugm1JtmZJvuZItmZJvgM3L9FB0ZGTVNad7fJ1V2d76sw51u+u4NaUaAL8bPFrf1Cs7rtyBIQQQgg3Oj+Kfc7e427Z39qdR2lqbmVJ2ji37E/0zBaFV0xMjNVNsC3J1izJ1xzJ1izJd+CmjB3J2FFBbNjb9Sj2rs42u8DJ5MhgkmOG7thdnVndd21ReMXGys2Cpki2Zkm+5ki2Zkm+A3d+FPst+6tobL70fiNXZnvweD2Fh0+yODUWpZTLtuvNrO67tii8rJ7w0s4kW7MkX3MkW7Mk38GZN8XBmaYW8g9WX/KaK7O98+//gSMv3c0L37y1Y9l3v/tdEhMTmTZtGgsXLuTkyZMu2583sLrv2qLwEkIIIbzJ1ZMiCPL36XFYicFqbdWcuSyTBY/+HD+fL892zZ8/n+LiYnbs2MHkyZN5+umnjbVBXMoWhVdwcLDVTbAtydYsydccydYsyXdwgvx9uWZSBBu6GMXeVdnmHTxBfdhkbpuZcMHyG264AT8/PwBmzZqF0+l0yf68hdV91xaFV3p6utVNsC3J1izJ1xzJ1izJd/DmTnFwpLqBA5X1Fyx3VbbZBU5GBvmRFR/R7TorV65kwYIFLtmft7C679qi8MrLy7O6CbYl2Zol+Zoj2Zol+Q7evPZhJf6658LLja7Itr6xmfeLj/F306IJ9Pftcp2nnnoKPz8/7r777kHvz5tY3XdtUXg1NjZa3QTbkmzNknzNkWzNknwHL2pUEEnRIWy8aFgJV2T7/s5yGs61sCSt66ETXn31VdauXctrr7025J52tLrv2qLwEkIIIbzRvEQHBV/UUHO6yaXbzS50MiFiBKnjQy95bd26dfzkJz9hzZo1DB8+3KX7Fb1TF9/U54nS09P11q1bu329ubm540ZB4VqSrVmSrzmSrVmSr2tsP3KS21/+mBe/Op3bZ7SdnRpMts0trfx680Ge+aCEf50/mb/973+Sk5NDVVUVkZGRrFixgqeffprGxkbCw8OBthvsf/WrX7nsM3k6d/RdpVSB1rrLm8ls8VOzb98+y2cbtyvJ1izJ1xzJ1izJ1zWmxYwiIjiADXsrOwqvgWa7v6KOR98qosh5ipunjWV55gS+M++NS9a77777Bt1ub2Z137XFpcbKSnPjoAx1kq1Zkq85kq1Zkq9r+PgorktwsKmkknMtrUD/s21p1fx60+fc/FIuh6vP8PLXUnn5a6mMCLTFuRWXs7rv2qLwEkIIIbzVvCkOas82U/BFTb/f+/nxeu741d94+v29XJcwhg8fmc3N08YaaKVwFVuUw8nJyVY3wbYkW7MkX3MkW7MkX9fJjB+Dv69i495KZk0M71O2ra2alR+X8swHJQT5+/KzO6dza0r0kHtCcSCs7ru2OONl9aOhdibZmiX5miPZmiX5uk5woB+zJobz1z1tw0o0NjayfPlyHA7HBUXCW2+9RVJSEj4+Ptzwg9/yo7/s4dr4CNY/ksVt02Ok6Oojq/uuLQqv/fv3W90E25JszZJ8zZFszZJ8XWteooODx09TWnWa/fv3s2zZMtatW3fBOlOnJnHPv/+coHHJHDpxmufuSOE3X0/HERJkUau9k9V91xaFlxBCCOHN5raPYn9+0uysrCzCwsI6Xj9SfYb/+vgUvyo6S8gwf37xtVQWp8XKWS4vZIvCa9y4cVY3wbYkW7MkX3MkW7MkX9caHz6ceEcwG/dWXJLta/lfcOOLmykuq+Wni6eRGDWSiOBAi1rq/azuu7YovCIjI61ugm1JtmZJvuZItmZJvq43d4qD/IPVjBjdNrDpsdqzOGsa+MGqYmaMD+WDR7JYeqUUvINldd+1ReHV06j2YnAkW7MkX3MkW7MkX9eblxhJc6vm9+u38n+fHebrr+Rz9lwLTy1M5g/3XUXM6GFWN9EWrO67thhOQgghhPB2qeNHM2qYPyt3NnJ2+06uiAxBhQ/n7pmXWd004UK2OOMVEhJidRNsS7I1S/I1R7I1S/J1PT9fHxYkR6FRRBX8mq0vfYuDB/YTGxvLK6+8wqpVq4iNjSUvL4+bb76Zr3zlK1Y32StZ3XdtMUm2EEIIYQdnz7XQ1NJKSJC/1U0Rg9DTJNm2OOOVm5trdRNsS7I1S/I1R7I1S/I1I8jflx1b861uhq1Z3XdtUXg1Nzdb3QTbkmzNknzNkWzNknzNkWzNsjpfWxReQgghhBDewBb3eLW2tuLjIzWkCZKtWZKvOZKtWZKvOZKtWe7I1/b3eO3atcvqJtiWZGuW5GuOZGuW5GuOZGuW1fnaovA6ceKE1U2wLcnWLMnXHMnWLMnXHMnWLKvzNVZ4KaVWKqUqlVLFnZaFKaXWK6X2t/8Zamr/QgghhBCexuQZr98BN1607HFgg9Y6HtjQ/v2gpaSkuGIzoguSrVmSrzmSrVmSrzmSrVlW52us8NJabwaqL1p8G/Bq+99fBW53xb7q6upcsRnRBcnWLMnXHMnWLMnXHMnWLKvzdfc9XpFa63KA9j8drtjowYMHXbEZ0QXJ1izJ1xzJ1izJ1xzJ1iyr8/XYSbKVUg8ADwBER0eTk5MDwMSJExk5ciRFRUUAhIeHo7XueN3Pz4/MzEwKCwupra0FID09nYqKCo4cOQJAfHw8gYGBFBe33X7mcDiYPHlyx2i2gYGBZGRksHXrVurr6wGYOXMmTqeTsrIyABISEvD19WX37t0AREVFMWHCBPLy8gAYNmwYM2fOJD8/n4aGBgAyMjIoLS3l2LFjAEydOpWWlhZKSkoAiImJITY2lvz8tlGLg4ODSU9PJy8vj8bGRgAyMzPZt28flZWVACQnJ9PY2Mj+/fsBGDduHJGRkR2zr4eEhJCamkpubm7HoHFZWVns2rWr4wbDlJQU6urqOjpjXFwcYWFhFBYWUl9fT1FRESkpKWzatAmtNUopZs+eTVFRETU1NQCkpqZSXV3NoUOHuj1OSUlJbN68WY5Tp+NUX19PTk7OoI8TQGhoqBynTsepsbGRhoYGj/p5stNxqq+vZ+/evR7182SX49Tc3NzxO81Tfp7sdJyampqoqakx+vPUE6PjeCml4oC1Wuvk9u9LgDla63Kl1FggR2ud0IftHAe+6GGVCKBq8C0WXZBszZJ8zZFszZJ8zZFszXJHvpdprcd09YK7z3itAf4B+HH7n6v78qbuGn+eUmprdwOVicGRbM2SfM2RbM2SfM2RbM2yOl+Tw0m8AeQBCUopp1LqPtoKrvlKqf3A/PbvhRBCCCGGBGNnvLTWd3Xz0jxT+xRCCCGE8GS2GLke+B+rG2Bjkq1Zkq85kq1Zkq85kq1ZlubrFZNkCyGEEELYgV3OeAkhhBBCeDyvLryUUjcqpUqUUgeUUi6Zfkh8SSl1SCm1Uym1XSm11er2eDOZu9SsbvL9oVKqrL3/bldK3WRlG72VUmqcUuojpdQepdQupdTD7cul/7pAD/lK/x0kpVSQUupTpVRRe7Yr2pdb2ne99lKjUsoX2Efb05FO4DPgLq31bksbZiNKqUNAutZaxpMZJKVUFlAP/L7TuHY/Baq11j9u/49DqNb6e1a201t1k+8PgXqt9bNWts3btY+5OFZrXaiUGgkU0Dbd2zKk/w5aD/kuRfrvoCilFDBCa12vlPIHcoGHgUVY2He9+YzXVcABrfVBrXUT8CZtc0EK4XHcOXfpUNRNvsIFtNblWuvC9r/XAXuAGKT/ukQP+YpB0m3q27/1b//SWNx3vbnwigGOdPreiXRWV9PAh0qpgvYpnIRrGZm7VFzgIaXUjvZLkXIpbJDaZyOZAeQj/dflLsoXpP8OmlLKVym1HagE1mutLe+73lx4qS6Weed1U891jdY6FVgAfLv9co4Q3uKXwCRgOlAOPGdtc7ybUioYyAb+WWtda3V77KaLfKX/uoDWukVrPR2IBa5SSiVb3SZvLrycwLhO38cCRy1qiy1prY+2/1kJrKLt8q5wnYr2+zvO3+dRaXF7bEVrXdH+j24r8Buk/w5Y+/0x2cBrWut32hdL/3WRrvKV/utaWuuTQA5wIxb3XW8uvD4D4pVSE5RSAcCdtM0FKVxAKTWi/UZPlFIjgBuA4p7fJfrp/Nyl0I+5S0XfnP+Htd1CpP8OSPsNyq8Ae7TWz3d6SfqvC3SXr/TfwVNKjVFKjW7/+zDgemAvFvddr32qEaD98doXAV9gpdb6KYubZBtKqYm0neWCtqmlXpd8B6597tI5QARQATwBvAv8CRgPHAbu0FrLDeID0E2+c2i7TKOBQ8A3zt/XIfpOKZUJbAF2Aq3ti79P231I0n8HqYd870L676AopabRdvO8L20nmv6ktX5SKRWOhX3XqwsvIYQQQghv4s2XGoUQQgghvIoUXkIIIYQQbiKFlxBCCCGEm0jhJYQQQgjhJlJ4CSGEEEK4iRReQoghQSn1Q6XUo1a3QwgxtEnhJYQQQgjhJlJ4CSFsSyn1A6VUiVLqr0BC+7J/VEp9ppQqUkplK6WGK6VGKqVK26duQSkVopQ6dP57IYRwFSm8hBC2pJRKo20qsRnAIuDK9pfe0VpfqbVOAfYA92mt62ibx+3m9nXuBLK11ufc22ohhN1J4SWEsKtrgVVa6zNa61q+nMs1WSm1RSm1E7gbSGpf/r/Ave1/vxf4rVtbK4QYEqTwEkLYWVdzov0OeEhrfQWwAggC0Fp/DMQppWYDvlprmZRYCOFyUngJIexqM7BQKTVMKTUSuKV9+UigvP3+rbsves/vgTeQs11CCENkkmwhhG0ppX4AfB34AnACu4HTwGPty3YCI7XWy9rXjwJKgbFa65NWtFkIYW9SeAkhRDul1BLgNq31PVa3RQhhT35WN0AIITyBUuolYAFwk9VtEULYl5zxEkIIIYRwE7m5XgghhBDCTaTwEkIIIYRwEym8hBBCCCHcRAovIYQQQgg3kcJLCCGEEMJNpPASQgghhHCT/w+ZFfPES96SBwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = Lille_2020_max.day\n",
    "y = (np.around(Lille_2020_max.values, 0)).astype(int)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "ax = plt.subplot()\n",
    "ax.set_ylabel('t2m (Celsius)')\n",
    "ax.set_xlabel('day')\n",
    "ax.plot(x, y)\n",
    "ax.grid(linestyle='--')\n",
    "for i,j in zip(x,y):\n",
    "    ax.annotate(str(j),xy=(i,j))\n",
    "ax.set_title('Max daily t2m for Sep 2020')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum temperature in September 2020 in this area was 32.1 degrees Celsius.\n"
     ]
    }
   ],
   "source": [
    "print('The maximum temperature in September 2020 in this area was', \n",
    "      np.around(Lille_2020_max.max().values, 1), 'degrees Celsius.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which day in September had the highest maximum temperature?\n",
    "\n",
    "Is this typical for Northern France? How does this compare with the long term average? We will seek to answer these questions in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='calculate_averages'></a>Calculate long term average of 2m temperature for September over Northern France\n",
    "We will now seek to discover just how high the temperature for Lille in mid September 2020 was when compared to typical values exptected in this region at this time of year. To do that we will calculate the mean and standard deviation of maximum daily 2m temperature for each day in September for the period of 1979 to 2019, and compare these with our values for 2020."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we select all data prior to 2020:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_past = Lille_t2m.loc['1979':'2019']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate the climatology for this data, i.e. the mean and standard deviation of maximum daily values for each of the days in September for a period of several decades (from 1979 to 2019)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do this, we first have to extract the maximum daily value for each day in the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_max = Lille_past.resample(time='D').max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate the mean and standard deviation of this for the 40 year time series for each day in September:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_m = Lille_max.groupby('time.day').mean('time')\n",
    "Lille_sd = Lille_max.groupby('time.day').std('time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot this data. We will plot the mean plus and minus one standard deviation to have an idea of the expected range of maximum daily temperatures in this part of France in September:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 't2m climatology for Sep from 1979 to 2019')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFNCAYAAADRi2EuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5xU1fnH8c+Z2Qq7dCz0Ioh0EOkiBAUEYoslsaIm9ihoRLCgMXajUUmMwYqKRlF/ARUBEewIKFIEFKV3YWH71DvP7487swzL9p3Zmdl93q/XwM6t584MzHfPOfccIyIopZRSSqnoc8S6AEoppZRSdYUGL6WUUkqpGqLBSymllFKqhmjwUkoppZSqIRq8lFJKKaVqiAYvpZRSSqkaosFLqRgyxkwwxnwZ9jzfGNMhlmUKluOIctXgec81xuwIvg59avr81WGMecAYc8AYszfWZVFKxS8NXqpOMsZsNcacHva8szFmjjFmvzHmoDFmgTHmxJoul4hkiMjmSB83VkGqCv4O3BR8Hb6v7sGMMd2MMQuNMYeMMdnGmO+MMWMjUM7i52kN3AZ0FZHjIn38Cpy/e/Aze8AYc9TgjMaYk4wxi40xOcaYX4wx54atuyQYdEOPQmOMGGNODq5vZIyZaYz5Nfi4r4xytAvum1TF6xhojPk4+G9wvzFmtjHm+LD1xhjzqDEmK/h4zBhjwtb/zRiz1hjjL17O4L53GWO2G2NyjTH/NcY0qEo5laoODV5K2RoBc4ETgWOB5cCcmJaobmoLrKvKjsYYZwmL3wc+xn5PjwFuBnKrXLrStQWyROTXUspWpSBSCT7gbeDqUs49B/gAaAJcA7xujOkMICKzgkE3Q0QygBuAzcDK4CH+AdQD2gH9gcuMMVdG6ToaAzOC52oL5AEvh62/BjgH6AX0BMYD14at/wWYDHxYwrEvBy4DhgAtgHRgekRLr1RFiIg+9FGnHsBrQABwAfnA5BK2aQII0DT4/D5gNvA69pfBWqAzMBX4FdgBjCrjnK2B94D9QBbwz+DyCcCXYdsJcELw51eAZ4GPguX8CjgOeAo4BPwI9AnbdwqwKVi+9cC5weUnAW7ACh4nO7i8IfBqsEzbgLsBRynlGgysAHKCfw8OW9ce+Dx43kXAv4DXg+s+BP5c7LVYA5xTbFlqsGwCFACbwsr+KZCNHcjOCtvnFeDfwLzgPqcXO2az4PEalfG+jAdWBY//NdAzbN3W4Pu7Pvh6vwyklXCM04OfpUDwGl7BDg6CHYS2B18fR/A13hb8zLwKNAweI7T9ldifpUPAdcApwdcrm+BnppzP9gmAFFvWPVguE7ZsIfC3Uo6xBLg37PkB4JSw53cCX5Sy7/bgdeQHH4PKuu4KXE9fIC/s+dfANWHPrwa+KWG/14H7ii17B7i92GfaDdSryf9/9KEPrfFSdY6IXIb9BfFbsX/Lf6yEzYYBe0UkK2zZb7FDW2Pge2AB9pdKS+B+4D8lnS9YE/MB9hdPu+D2/61gcS/E/tJqBniApdg1Ec2wv0ieDNt2E3AqdqD6K3atxvEisgH7S3xp8HobBbefHty2A3Aado3AUTUZxpgm2AHqGaBp8JwfGmOaBjd5A7uGsCl2QL0sbPeZwKVhx+oVvP554ecQEY/YtS0AvUSkozEmGbvGaiF2bdWfgVnFmoAvBh4EMoHiTalZ2DUgrxtjzjHGHFvsuvoCL2HXmDTFfv/mGmNSwza7BBgNdMQO2ncXf31EZBFwJrA7+PpOCFt9GnZ4HI0dZicAI7Bf8wzgn8UONwDoBFyEHbDvwg523YALjTGnFT9/BZhSlnU/aqExbbE/+6+WcYwS9w0aFvy7UfC1WErFrrs0wziyBrQbsDrs+ergsoowHH0dqdivt1I1J9bJTx/6iMUDuzbj9FLWtQJ2AX8IW3Yf8HHY899i/0bvDD7PpJTaFezf+vcDSSWsm0DZNV7Ph637M7Ah7HkPgrVXpVzHKuDsUs7jxA5yXcOWXQt8Wnx77CC1vNixQ1+obQA/YbUG2LUNoRqvVOAg0Cn4/O/As2WUOfz6TwX2EqyFCy57k2BNRvD1ebWc97kV9pf8Juwaqc/DyvJvitX6AD8Bp4V9Rq4LWzeWYE1cCecZDuwMe94ueC0dwpZ9AtwQ9vxE7CbCpLDtW4atzwIuCnv+LjCxnOstqcYrGbvpcHLw51GAF1hQwv73hD4Dxd7P97A/4ycEX0tPKecPXUdSRa67nGvpGfzsnBq2zAK6hD3vFDyfKaHM9xVb9kdgY7CMDbG7FggwqKxy6EMfkX5ojZdSYYwxzbFrWJ4VkTeLrd4X9rMLOCAiVthzsH+bL641sE1E/FUoUvFzFn9edD5jzOXGmFXBTuTZ2LUSzUo5bjMgBbsWLmQbdm1UcS2KbRe+bQvgoIgUhq3bEfpBRDzYfY8uNcY4gD9g1xpWRAtgh4gEyijjDsogIjtF5CYR6YjdZ6iAw7U5bYHbQq9X8DVrHTxvScffVmxdRYTvX/x13IYdusJr4ir8fleUiPiw+0WNww6yt2G/JztL2Pxy7FrKcDcHz/0zdl+xN0vZtzQVue4jGGNOwG5iv0VEvghblQ+Ed4hvAOSLyFE3FJTgJeyyf4pdi7YkuLwy16JUtWnwUnVVSXd+NcYOXXNF5MEInmsH0CaaHayDTUTPAzdh90trBPzA4aaV4td7ALvWoW3YsjbYNX3F7S62Xfi2e4Amxph6YetaF9t2JnaT3UigUOzmp4rYDbQOBrbSyliRL1x7Q5Ed2P3PQs1kO4AHRaRR2KNescAdfi1tgmWqjPDyFX8dQ7WF+4gyEVkjIqeJSFMRGY3d5Lc8fBtjTKjT+TvF9j0oIpeIyHEi0g37e+OIfcM3L2FZpa47+FlehF0bWTykr8PuWB/SiwrejCEiARG5V0TaiUir4H67KPkzr1TUaPBSddU+7C8fAIK3lS8AvhKRKRE+13LsgPKIMaa+MSYt+CUXSfWxv/T2AwTvOgvvh7MPaGWMSQEI1tS9DTxojMkMftndit1EU9w8oLMx5mJjTJIx5iKgK/CBiGwDvgXuM8akGGMGYTfDFgkGrQDwBBWv7QJYhl1DNdkYk2yMGR48doX6xxljGhtj/mqMOcEY4zDGNAOuAr4JbvI8cJ0xZkBwqIH6xphxxpjMsMPcaIxpFezndifwViXKX9ybwCRjTHtjTAbwEPBWFWtCjxAsfxp2LSbBz1hq2PqewWX1jDF/AY7HbqoNdwXwrojkFTt2R2NMU2OM0xhzJvadhQ+UUpT92O91+Fh0Fb5uY0xLYDHwLxF5roTjvwrcaoxpaYxpgV1790rY/snB18EBJAWv2Rlc1yR4LcYY0xW7r+L9xWpUlYo6DV6qrnoYuDvYxPQX4FzsO8iuNEeOadSmuicKhpzfYveP2Y7dtHFRdY9b7BzrsYPNUuyQ1QP7LsiQxdi/4e81xhwILvszdrDZjN0x/Q3s5pjix87CvvvvNux+R5OB8SISOs4l2P3YsrC/kN/C7j8W7tVgmUoKdqVdkxc4C7vj+gHsOzwvF5EfK3gIL3Z/nkXYQ0j8ECzXhODxvwX+hN0H7BB2R/wJxY7xBnYt6Obgo7TAUREvYQfPz4Et2HfU/bkaxwvXFrs5MFT748LurxZyGXb4/xW75vGMYDMwYAc17Bs5ijczApyMfRdvHva/m0tEpMRapmCT84PAV8F/WwOp3HX/ETu03Rv+7zBs/X+wb7hYi/1+fsiRN7U8H7z2P2DfmODi8M0ezTh8B+xHwEsiMqOUcigVNaZiTeNKKVUxxpi3gB9F5N6wZZdjDwMwNHYlqxxjzFbgj2LftaiUUhGhNV5KqWoxxpwSbMJxGGPGAGcD/wtbXw97UE6tXVBK1XkavJRS1XUc9p1i+dhjfV0vwel+jDGjsfv97MNutlNKqTpNmxqVUkoppWqI1ngppZRSStUQDV5KKaWUUjUkagM6RlKzZs2kXbt2sS6GUkoppVS5vvvuuwMi0rykdQkRvNq1a8e3334b62IopZRSSpXLGFN8mrUi2tSolFJKKVVDoha8jDGtjTFLjDEbjDHrjDG3FFv/F2OMBKfxUEoppZSq9aLZ1OgHbhORlcG5z74zxnwsIuuNMa2BM7CnT1FKKaWUqhOiFrxEZA/23GCISJ4xZgPQElgP/AN7vrc50Tq/UkqpxOTz+di5cydutzvWRVGqTGlpabRq1Yrk5OQK71MjneuNMe2APsAyY8xZwC4RWW2MqYnTK6WUSiA7d+4kMzOTdu3aod8TKl6JCFlZWezcuZP27dtXeL+od643xmQA7wITsZsf7wKmVWC/a4wx3xpjvt2/f3+US6mUUipeuN1umjZtqqFLxTVjDE2bNq10zWxUg5cxJhk7dM0SkfeAjkB7YLUxZivQClhpjDmu+L4iMkNE+olIv+bNSxwKQymlVC2loUslgqp8TqN5V6MBXgQ2iMiTACKyVkSOEZF2ItIO2An0FZG90SqHUkopVRnZ2dk8++yzRc9XrVrFoEGD6NatGz179uStt96KYemi65VXXuGmm26q8v6ffvop48ePB2Du3Lk88sgjkSpapV177bV89dVXZW7z0ksv0aNHD3r27En37t2ZM8fuej5hwgTat29Pr1696Ny5M5dffjm7du2KSLmiWeM1BLgM+I0xZlXwMTaK51NKKaWqrXjwqlevHq+++irr1q1j/vz5TJw4kezs7BiWMDGcddZZTJkyJarnmDBhAp9++mmJ65YtW8bAgQNL3Xfnzp08+OCDfPnll6xZs4ZvvvmGnj17Fq1//PHHWb16NT/99BN9+vRhxIgReL3eapc5asFLRL4UESMiPUWkd/Axr9g27UTkQLTKoJRSSlXWlClT2LRpE7179+b222+nc+fOdOrUCYAWLVpwzDHHEOp73K5dO+68804GDRpEv379WLlyJaNHj6Zjx44899xzRx1769atdOnShT/+8Y90796dSy65hEWLFjFkyBA6derE8uXLASgoKOCqq67ilFNOoU+fPkU1MVu3buXUU0+lb9++9O3bl6+//hqwa5qGDx/O+eefT5cuXbjkkksQkTKvc/bs2XTv3p1evXoxbNiwouW7d+9mzJgxdOrUicmTJxctv/766+nXrx/dunXj3nvvLVo+f/58unTpwtChQ3nvvfeKlofXnk2YMIGbb76ZwYMH06FDB9555x0AAoEAN9xwA926dWP8+PGMHTu2aF11bNiwgc6dO+N0Okvd5tdffyUzM5OMjAwAMjIySuwkb4xh0qRJHHfccXz00UfVLhsiEvePk08+WVTNCwQCsS6CUqoOWr9+fUzPv2XLFunWrVuJ65YtWyZdunQRy7JERKRt27by7LPPiojIxIkTpUePHpKbmyu//vqrNG/evMRjO51OWbNmjViWJX379pUrr7xSAoGA/O9//5Ozzz5bRESmTp0qr732moiIHDp0SDp16iT5+flSUFAgLpdLREQ2btwooe/HJUuWSIMGDWTHjh1iWZYMHDhQvvjiizKvs3v37rJz586ic4iIvPzyy9K+fXvJzs4Wl8slbdq0ke3bt4uISFZWloiI+P1+Oe2002T16tXicrmkVatWsnHjRgkEAnLBBRfIuHHjio514403iojIFVdcIeeff75YliXr1q2Tjh07iojI7Nmz5cwzzxTLsmTPnj3SqFEjmT17toiI3HPPPTJnzpwyr+GKK66QJUuWHLX8iSeekBdffLHMff1+v4waNUpat24tEyZMkLlz5x5x3FA5Qm655RZ55JFHjjpOSZ9X4FspJdMkxFyNKjb8ASHJoZ1clVKxM3H+RFbtXRXRY/Y+rjdPjXmq0vvt2bOHyy67jJkzZ+JwHG4wOuusswDo0aMH+fn5ZGZmkpmZSVpaGtnZ2TRq1OiI47Rv354ePXoA0K1bN0aOHIkxhh49erB161YAFi5cyNy5c/n73/8O2Hd6bt++nRYtWnDTTTexatUqnE4nGzduLDpu//79adWqlX2NvXuzdetWhg4dWur1DBkyhAkTJnDhhRdy3nnnFS0fOXIkDRs2BKBr165s27aN1q1b8/bbbzNjxgz8fj979uxh/fr1BAIB2rdvX1QjeOmllzJjxowSz3fOOefgcDjo2rUr+/btA+DLL7/kggsuwOFwcNxxxzFixIii7e+///4Sj7NgwQLuuOMOALZv386XX35JRkYGqampLFu2rGibl19+udRrB3A6ncyfP58VK1bwySefMGnSJL777jvuu+++EreXcmoQK0qDlzpC6INV4LUo8PjJSE2iXopTw5dSqk7Lzc1l3LhxPPDAA0f1G0pNTQXA4XAU/Rx67vf7jzpW8W3C9w9tLyK8++67nHjiiUfse99993HssceyevVqAoEAaWlpJR7X6XSWeO5wzz33HMuWLePDDz+kd+/erFq1qtTjbNmyhb///e+sWLGCxo0bM2HChKJhFCr6/RB+3NB3TVXCzOjRoxk9ejRgN2FOmDCB4cOHF60vLCwkOzubFi1asGPHDn77298CcN1113HdddcdcSxjDP3796d///6cccYZXHnllaUGr++//56RI0dWurzFafCKQ16/RbLT/m2qJgOPiBAQIcflw2fZ/xgKvX7qpZTeRq6UUtFUlZqp6srMzCQvL6/oudfr5dxzz+Xyyy/nggsuqJEyjB49munTpzN9+nSMMXz//ff06dOHnJwcWrVqhcPhYObMmViWVe6xpk6dSv/+/Tn33HOPWL5p0yYGDBjAgAEDeP/999mxY0epx8jNzaV+/fo0bNiQffv28dFHHzF8+HC6dOnCli1b2LRpEx07duTNN9+s1HUOHTqUmTNncsUVV7B//34+/fRTLr744kodo7glS5YU1Zy1bt26KFAWt3v3bvbu3Uvfvn0B++7Vtm3bHrWdiDB9+nT27NnDmDFjqlU2qIEBVFXFiQhef4BDhT4O5Hvx+ANFy6N9XgC3L0BWvrcodAEExF4e7TIopVS8aNq0KUOGDKF79+7cfvvtvP3223z++ee88sor9O7d+4jaoWi555578Pl8RcMc3HPPPQDccMMNzJw5k4EDB7Jx40bq169f7rHWrl3LcccdNVwmt99+Oz169KB79+4MGzaMXr16lXqMXr160adPH7p168ZVV13FkCFDAHvKnBkzZjBu3DiGDh1aYnApy+9+9ztatWpF9+7dufbaaxkwYEBRM+e0adOYO3dupY4H8NFHH1UoIPl8Pv7yl7/QpUsXevfuzVtvvcXTTz9dtP72228vGk5ixYoVLFmyhJSUlEqXpziTCF+o/fr1k2+//TbWxYgqEUGArHwPgbC3JNlpyExLJtnpQEQiXgMWOm+uy1cU9IpLchiaZqSWuE4ppSJtw4YNnHTSSbEuRq0xevRoFixYEOtilCo/P5+MjAyysrLo378/X331VYlBsaL69u3LsmXLKjV/YnWU9Hk1xnwnIv1K2l6bGuOEMYacQu8RoQvAZwkHC7ykJzvJSE0CIhu+vFaAXJfvqPOG8wcEj98ixenQvl5KKZVg4jl0AYwfP57s7Gy8Xi/33HNPtUIXwMqVKyNUsujQ4BUnCr3+UmucAFw+C7fPon6wsztUvf9XqJYz3+On0Ft+/wCAQo9Fan3t6xVp0ajFVEqpRFLaAKi1lQavGBMRrICQ5y777hMAwQ5LLp9FZloSqUnOKn1x+wN2B3qrrGquYrxWAJ8VIMlhNChUU+g9y3f7MMZUO0grpZRKHBq8YihU85Tt8lVqPysgZBf6SEmyA1hSBb6wQ+cq9Frke8oPeSUp8PhpVK/6HQvrsqL3vPDwzRMun0VmahKpyVUL0krVRvpvQSWCqvST17saY8gYQ57bX6map3Bev30XYp7bRyA4Im5J7GEi4FChr8qhC8DjD2AF9A7HqgrdyHAoLHRBMEi7fBwq8BZ9FvQ1VnVZWloaWVlZ+u9AxTURISsr64ix1CpCa7xiRETw+AO4fBXrY1WWQq+Fy2eRkZpEevLhZqvQb4wev92BPhL/hRV4LBqk18ydIrWJiGCJkF3gwyrly8RrBcgK3UiRlgT6G7+qo1q1asXOnTuL5kNUKl6lpaUVzRZQURq8YiBUA5VbySbGso8JeW4/Lq9FZloyKUnGHiai0Iu7jE77leXyWWSkJeHQQFApPkvIdnmpyC/wkbyRQqlElJycXOJkxUrVBhq8YiTH5Y1IDVRx/oBwqNBLitOBPxAoc5iIqir0WsGhLVRFuHxWpUP2ETdSaP8vpZSqNbSPVwwUePxHjA4fDV4rOqEL7KEvtO9F2cKH7KhOzab2/1JKqdpFg1cNCk0JVFDBsbPilYhdi6MBoGSh1yXH5aOgGjczhAv1/wr11dPXXimlEpMGrxoSuqMtx+WNdVEiotCT2OExWkLvc3ahD3cEbpwozuWzOJDnodBrB18NYEoplVg0eNUQY0y5U/MkEit4V6Z+8R8WumniYIEXrxW5GxqOOg92E2ZWgRdvBG+cUEopFX0avGpIoafsKYESUYHHr529g0QEf0A4WOCp8rhslRXq/5Xn9mntl1JKJQgNXlEmIvisAHkR6usTT/wBwetPrL5eoYAS6TJ7/QEOFRw9yXlNKPRaRbMfJNJ7oZRSdZEGrygK72RdWxV4rYSq9fJbQqHXwusP4C9hFP6qBJdCr5/sCA1QW1Vef4CDBd4yZzBQSikVezoYUxQZYyo9GXWi8foD+K0AzjiePDs0/lWBx1/ilElOhyHJYYJ/O0hy2j+HDxIbCjOhawwfLqIwTu5S9QeErAIvjdJTSEmKz/dCKaXqOg1eURKaEigad7bFmwKvn4bp8Tl59uG7DEvviG4FJCwcH36/HAacDgdJoWDmtIOZ02GHmhyXL+767YnYc0E2SEsmPTjqvVJKqfihwSsK7LvbJKJTAsUzty9ARqrgMPE3tY3PEnJcVet7FRAIWAFKys4GYtq0WJ5ctw9/IFA0w0C8vS9KKVVXaR+vKMkujG2fn5pW6I2fOxxDzYAFHj+HCqPT4T0R3ttCr1XUv1D7fSmlVHzQ4BUF+R4//lrcr6skLq9FIA6+3ENNi4cKvSX256prPEWd7jV8KaVUPNDgFUF2vy4rbjpb1yTBDl+x/nL3WUJWvkcHFg0TGl/MH+X5QZVSSpVPg1cEea0AOYV1o19XSQq9salhCo3LlR/FpsVEFxA4WOjFVQd/KVBKqXiiwStCCr3+Otevq7hADCbPtkOX3acuUhNS12a5bnuke9CmR6WUigW9q7Ea4nEsp1gr9FqkJ9fcMAbVuWuxrir0WlgBoWF6ctEYZ0oppWqG1nhVUfio9Bq6DrMC0Z88W5sWq0873SulVGxo8KoCe5wuOFjgjbsBNONBNIeW0KbFyNFO9yqWNPCrukqDVyWJyOEvLK1qKZHPErwRrvUKHctrBcgq8OC1NPBGQqjTvSfBJjtXiU1E8FtS4nypStV2GrwqKbyJRpUuErVe4f8h+yx7JoDsQp++9lGQ6/Ih2uyoakAodB0q9JLnip+Bl5WqKdq5vgLKm2RZHc1TxcmzQ6+1iOCzArh9ATx+S8NWlAXE7q/YuH58zrmpaofw0CXYNdhun0VqkkMDmKozNHiVI1QDkOvy4aoDE15HUoHXomF6crnbhYctrxXA4wvg9lto5UvN8loBCjx+6qfqfwsq8oqHrpA8t4+UjFQ0dqm6Qv+HLUNo+pnsQh8+7VNUaW6fRUZqUomTZ4eHLY/frtXy+AJ1ehy0eJDv8ZOS5CCpkjWVSpWltNAFdm1rgcdPZlr5v6QpVRto8CqFiGCJkF3ow9J2rior9B7+D7V42HL7LLsTfozLqI6UU+ijSUYK6BhfKgLKCl0hhV6LtGSnBn5VJ2jwKoXPErJdXm3uqiaX1ypquvIE+2vpEBzxzRL7RoZG9bS/l6qeioSukDy3jyb1U2ukXErFkgavEri8FrnuujvnYiQJcCDPo7VaCcbjD+Dy+klLdmoNhKqSyoQusH/Z1c+cqgt0OImgoul/3D4NXRGmoSsx5br9WAHRISZUpVU2dIXkefwIOqyJqt00eIXJLvRSoNP/KFUkx6UTasebeH8vqhq67H0h361je6naTZsag/LcPu17pFQx/oCQ5/bToALDgqijhUJSodfCEgEJ1QBLUf9RCf4hoTVyuJZYgitDz1OcDlKTHKQkOXE6TNE54iWoVCd0hbh8FunJTpKc2tFe1U4avIL0xkWlSubyWaQkOXSQy0oKzema4/Lii9B8mPbQKwHAT5LDkJrkIDXJSZKTojuGY/UeRSJ0heS6fTTRwXxVLaXBSylVrlyXj6YZqTiIn9qVeOezhBxX9KYX8wcEv9eiwGthDKQ6HaQkO0lNcuAIhjA4egy9aIhk6AL72gq9FvVStKO9qn2i1sfLGNPaGLPEGLPBGLPOGHNLcPnjxpgfjTFrjDH/Z4xpFK0yKKUiQ7BrblTZDjct+jlUWHNzuoqA2x8g1+Vjf56HgwUeCrwW/rACRKtvWKRDV0iBx09A5w9VtVA0O9f7gdtE5CRgIHCjMaYr8DHQXUR6AhuBqVEsg1IqQnyWUKBzlZYqFBByXD7y3LF9nULv1cECL/vz3OS67D6sxUOMSPXuWhUR/IHIhy6ww36e26c1XqrWiVpTo4jsAfYEf84zxmwAWorIwrDNvgHOj1YZlFKRVeC1+3slO7W/V7hQf67sQu8RtUzxICB2P73QXLNJDoMz9DBhP5cwanxZzZVFoasg8qErJDSdWIp+3lQtUiN9vIwx7YA+wLJiq64C3ipln2uAawDatGkTxdIppSojJ9jfS6cUOsxrBchx+RJipgt/QEoNh8ZwVBgLD2jh73e0Q1dInttPU+1or2qRqAcvY0wG8C4wUURyw5bfhd0cOauk/URkBjADoF+/fgnw35lSdYN9p56PxnV8SqHQHYQFHj/5taQJVgT8UnowcxhwOgwOY2psnlUrIBR4LTJS9V4wVTtE9ZNsjEnGDl2zROS9sOVXAOOBkaI9J5VKOF5/gAKPv2gezrom9N9WdqG3To3/FxAIWEJNz0dR4PGTnmzfram1rCrRRfOuRgO8CGwQkSfDlo8B7gDOEpHCaJ1fKRVd+R4/PuvoDtu1nYhgiZBVULdCV6zl6oj2qpaI5jnEnxQAACAASURBVF2NQ4DLgN8YY1YFH2OBfwKZwMfBZc9FsQxKqSiqi1MKef0BDuZ7seKsE31t5/UHcPusOvVZU7VTNO9q/BIo6deTedE6p1KqZlkBIdfto2F67e7vFfqyL/D4dT7XGMpz+0jNSI11MZSqFp0kWylVLW5fAJe39tZEiNhzJWa7fBq6Yiwg1JobGVTdpcFLKVVteW4fVjUH44w3oWuxAsLBfC9e7c8VFwq9Fv462LdQ1R4avJRS1SZAdqEPK5D44StUfr8l5Lp8ZBV4sRL8mmqbXB3RXiWwunkvuFIq4qyAfadfZloS9VKSisa5SgShslqBAC6f3YlbO8/HL58luLwWack6or1KPBq8lFIRlef24/UHaJCeHNej24fCVkAET3BKHZ+lYStR5Ht8pCbrDAoq8WjwUkpFnMcfICvfQ8P0FFKS4udLMRS2RASvP4DLZ+lYXAkqIJDv9tsBX6kEosFLKRUVAYFDhV7qpziLRriPVc1EKHD5A4Lb58fts9CWxMTn8lmkpzhJKmGCb6XilQYvpVRUFXgtvFaAhukpOKi5ZqHD/bYEd7ApUftt1T65Lh9NdBJtlUA0eCmlos5nCVn5HhqkJ5OW7Ixax/vwpkR3sJO819KmxNrMH7A72qenOLXWSyUEDV5KqRoh2FMMef0BMtMid9dj+HF8ln1Xosdn1fA0ziqW8j1+UpOdNVqjqlRVafBSStUoly/U9JhMsrPqX5JFTYkiuD1+XNpvq84S7EF8G9XTJkcV/zR4KaVqnBUQDhZ4yUhNon5qxWu/woeAcPss3DoEhAry+O2azpQkHdtLxTcNXkqpmMn3+O3ar7RkKKWZSIeAUBWV6/bTLENrvVR80+CllIoprz9AVoHd8T41yVm0vPgQEC6fhc7co8oSECHf4yczTcf2UvFLg5dSKuYCYs/1WC8lQEZqEgGhqCnRrx23VCUUei3Sk504dWwvFac0eCml4kah18Ll1TsSVfXkun00qZ8a62IoVSJHrAuglFLhNHSp6rIn0fYj2jat4pAGL6WUUrVOnsePgIYvFXc0eCmllKp1RCDP7dd+XiruaPBSSilVK7l9Fl5/QGu9VFzR4KWUUqrWynX7Yl0EpY6gwUsppVStZQWEQq8V62IoVUSDl1JKqVot3+PHCmiTo4oPGryUUkrVerna0V7FCQ1eSimlaj2vP4DbZ2mtl4o5DV5KKaXqhDy3T8f2UjGnwUsppVSdEBC7v5c2OapY0uCllFKqznB5LXyWdrRXsaPBSymlVJ2S69KxvVTsaPBSSilVp/gDgks72qsY0eCllFKqzsl3+wmIdrRXNU+Dl1JKqTpHsO9y1I720RUKthpwD9PgpZRSqk7y+AN4/NrkGA2h19QKCPkePz5LA1hIUqwLoJRSSsVKnstP04yUWBejVhARjDEERHD7LFxeC3/ADloFgNNhSE92kpbsxOkwRdvXNRq8lFJK1VmWCLluPxmpSTgMRwWBUA1NogWEmip3KDyJCF4rgMtr4fEHStw2VPuV7/GTkuQgPdlJapKjaP9Ee42rSoOXUkqpOs3ts3D7LAAcBhwOg9MYnA6DI/i3M7jM4Tg6HMRjOPMHBCsgJDsdOINljmQ5w2u3XB4/Lp9FoBKtiF5/AK8/gDGQnuwkPdlJkrNuBDANXkoppVRQQCBgCX5KTxGhAOZ0hP1sDElOR4m1ZjUlFFoKPX7yPP6i5Q4DSU4HyQ67jFUNY+G1W25fAJfvcN+tqpcZCr0WhV6LJOfhpkhHLa4F0+CllFJKVYIlgmUJwUqyIgbITEsiPSWpxkNDKEBlF3qPauoLSLCGCYDDNXuhMJbsdJBURhgLXYvPElw+Px6fVUYsrTq/JeRZfvLcftKSHaQnJ5GSVPsCWIWClzGmH3Aq0AJwAT8Ai0TkYBTLppRSSiUMAXLdfjz+AA3Sk6GGAoOI3ayY7fJhVbC9r6JhDOymWJfPqvCxI8HtC+D2eY/okB/L2sRIKnM4CWPMBGPMSmAqkA78BPwKDAU+NsbMNMa0iX4xlVJKqcTg8QfIyvfgDdY8RXsIBY8/wMECb7WDUSiMFXgtsl0+DuR7OJDvId/jr9HQFS7UIf9A2OuZ6Mqr8aoPDBERV0krjTG9gU7A9kgXTCmllEpUAYFsl490f4DMtMg3PYbCXL7HT6HXKmfr2iHb5aMBdmf8RFZm8BKRf5WzflVki6PiiYhQ6CugfkpGrIuilFIJyeWz8FoBGqYnk+yMTPASEUQg2+Wtduf2RJPr8iEi1EtJ3C7qFRq53hjzmDGmgTEm2RjziTHmgDHm0mgXTsXOhgPrOPvtkXR7rjVvrXs91sVRSqmEZQWEgwVe8j3+YGiqXljyWUJWgafOha6QPLeffLcv1sWosopOGTRKRHKB8cBOoDNwe9RKpWKm0FfIg1/ewxmzBvLLwY10bd6DWxb+iamLJ+K1vLEunlJKJawCj59DhV4CVQhfoe0PHyMaJUwcBV6rqPYr0aYhqmhdXXLw77HAmyJysLy2amNMa+BV4DggAMwQkaeNMU2At4B2wFbgQhE5VPmix5cfD6xn6pKJ7MzdjsM4cBonDuPAhP1sP+zxSY7axuHEgb1Ns3rNubr39fQ9vn+NXsPirQuZ8sktbM/dykVdL2PasIdomNqIB7+8h39/9xQ//Lqa58fP4riMFjVaLqWUqi18lpCV763UsBNlDRVRl9mDtgoN05MTasgJU5GkaIx5BDgHeyiJ/kAj4AMRGVDGPscDx4vISmNMJvBd8BgTgIMi8ogxZgrQWETuKOv8/fr1k2+//baCl1Q1Vf1AByTAf1Y+wyNf3UdGSgNGtD0dQQhIAEssAhJAin4+vFwkgBWwCBAgIAECAXvbAAE2HfyZbM8hhrX5DZMGTGVQq6FRuOLD9uXvYdpnk5mz8R1OaNyZR0dOZ0jrYUdsM+end5j08XVkJGcyY9xrDIxymZRSqrZLTXLQID0ZQ+nDJFRlqIi6JsXpoFE9u34oXsKXMeY7EelX4rqKVtEZYxoDuSJiGWPqAQ1EZG8lCjEH+GfwMVxE9gTD2acicmJZ+8Zr8NqRu42bF/yJpTu/YEzH8Tx++r9oXu+YapelwJvPzDXP8+/vnmZ/4T4GthzCxAFTOK3NyIh+qAIS4NU1L/Dgl/fgtTzc0v8Obux3K6lJqSVu/+OB9Vz1/kVsz93KfcMe5ere18fNh1wppRKRw0CD9GRSk5wl1tq4fcEmtRiVL1EkOw2N6qWUGWJrUrWDlzHm8pKWi8irFSxAO+BzoDuwXUQaha07JCKNy9o/3oKXiPDW+te5+9PbEBEeGPEEv+96WcTfbJffxay1L/Pst0+yO38XvY89mUkDpjCqw7hqn2vd/jXcvugmVu5dwamtR/DoyGfo0PiEcvfL9eTw5/l/ZMHmD/hdl9/z+On/ol5yvWqVpa75tWAvi7YsYOnOLzi9/RjOPvH8WBdJKRVj6clOMtOO7P1Tl4aKiASnw9C4XkpcDLQaieA1PexpGjASWCki5X5jGGMygM+AB0XkPWNMdkWClzHmGuAagDZt2py8bdu2cstZHRUNXgcK9zP5k5uY98tcBrYcytOjn6dtw3ZRLZvH72H2hlk8s/xxtudupWuzHkwccAfjTjgHp6Ny45kU+Ar4+9IHmLFyOo3SmnD/aY9yXpffV+pDGpAATy9/jMe+vp+uzXvw0vj/0rZR+8peVp0hIqzbv4aFm+fx8ZZ5fL/X/iWiXnJ9Cn0F/KHbFTww4gnqJ9ePcUmVUrHkdBgapifjMIacOjhURCQ4jKFx/WScxsQ0fEWkqbHYARsCr4nIWeVslwx8ACwQkSeDy34iQZsaF276kFsX3UCuJ5spg+/j2r43Vzr4VIc/4Of/fnybp5c/yi+HNtKpyYncfMpkzu1yIUmO8u+TWLh5HlMXT2RX3g4u7X4Vd536NxqnNalyeT7ZsoAbP7oSgH+d+TIj24+u8rFqG5ffxVc7PmPh5nks2jyP3fm7MBj6HNePMzqMZVSHsXRuchJPfPMgTy9/jBOadOa5sa/SrXnPWBddKaUSmjHQuF4KSY7Yha9oBK9kYI2InFTGNgaYid2RfmLY8seBrLDO9U1EZHJZ54t18Mr35jHts8m88cMrdGvek3+OeZGTmnWPannKYgUsPvj5/3hq+aNsOPADbRu25+ZTbueCrpeQ4kw5avs9+bu4e8ltfPjLHE5s2pXHRj7DgJZDIlKWbdlbuOqD37N+/1omD57GLf0n4zAVHaWkdtmXv4ePt3zEx5s/4vPti3H5C6mXXJ/hbUdyRodxnN5uNM3rH3vUfl9u/5Qb519JtvsQ04Y9zFW9rot5NblSSiUyAzSql0xKUmxGuY9EU+P7UNS3zwF0Bd4WkSll7DMU+AJYiz2cBMCdwDLgbaAN9lRDF5Q32XYsg9eyXV/x5/l/ZGfedm7sdyt/GXh3qZ3Pa1pAAny8eR5PLnuY1ftW0iKjJTeechsXd59AelI6VsDi5dXP8cjXf8Vv+bh14J1cd/ItJYaz6ij0FXL7oht598f/MqrDOP455kUapDaM6DnikYiw9tdVLNw8j4Wb57Hm15UAtMxszagOYxnVYRyDWw2r0OflQOF+bll4DZ9smc+YjuP5x6j/VKs2UimlFDRMTyYtBlMMRSJ4nRb21A9sE5GdESpfuWIRvDx+D48tvZ9nv/0HbRq2Y/roF+jfcnBUy1BVIsKn2xbxj2UPs3z3UprXO5Yre13Lgs0fsHrfSka0PYNHfvN0VPthiQgvrvo3931+B20atOPF3/6Xk5p1i9r5oi00XVK25xA57myyPdnkuA8V/b3p0M98vHkeewv2YDCcfHx/RnUYxxkdzqRL025VqrESEWZ8P50Hvrib5vWO4dkzX9FhO5RSqpoaBMdMq0kRb2qsaTUdvNbvX8tN869m/YG1XNr9Ku477REyUjKjev5IEBGW7vqSf3zzMF/sWMIx9Y7j/uGPcXbn82us6WrZrq/40weXku/L4x9nPBd3d+wt3/U1vxzaSLb7EDmebLLd2eQUhSv77xxPDjmeQ/gD/lKPUz85g+FtT2dUh3GMbD+aZvWaR6yMq/et5Lp5l7MtZwu3DriTSQOm1GhfQqWUqm0yUpOon1pz4avKwcsY86WIDDXG5MERw4gYQESkQWSLWrKaCl6FXh///u4pHlt6Pw1TG/PEGc8yqsPYqJ43WrZkb6JZenMyU2vkLTrC3vzd/OmDS1ix5xuuP3kidw39W4U6/0dTlusAdy+5jf/76e2iZU7jpGFqIxqmNaJhaiMapTW2n5ewrGhdWmMapTYiM7VBVPuy5XvzmPLJLbzz45sManUq/xrzEi0yW0XtfEopVdvVS3GSmZZc/oYRoDVeFbB6909cN+9qvtn1FeNOOJtHR06PaC1GXeO1vNz72WReXv0fhrQ+jWfHvMyxGcfHpCwf/vw/7lh8CznuQ0waMJULul5Mo9TGZKRkxn0n9rfXz2LK4ltIdaby1Kj/MLrj+FgXSSmlElZN1XyVFbwq9Cu7MaajMSY1+PNwY8zNxphG5e2XKGavm83QV/qxbv9anhn9Ai+Mf1NDVzWlOFN4+DdP8fSo5/lu9zIGvdKDJ755iAJfQY2V4UDhfq798DKu/uAPtMhoyYKLv+LWgVNp3aAtmakN4j50AVzY9RI+vmQprRq04Yq5F3DXkltx+92xLpZSSiUkKw4qmyraVvIuYBljTgBeBNoDb0StVDXs2IxjOfn4U1h82XIu7HpJQnwhJ4qLul3K4stWMKLtGTy+9G8Meqk7r699qcz+U5Hw/sb3OO3Vk5n3yxymDL6PD3//GV2b94jqOaOlY+NOfHDRp/ypz428uOrfjP/vcH45uDHWxVJKKVUFFb2rcaWI9DXG3A64RWS6MeZ7EekT/SLWTFPjoQIPXh0lOKpW7F7KXz+fyrd7lnFi067cPfQBTm8/JqJB90DhfqYunsj7P79Hz2P68PToGTEdcy3SFm6ex8QF1+C23Dw04h9c1PVS/UVBKaUqKD3FSYMa6OdV7aZGwGeM+QNwBfZI9AA100OthuiXV/Sd0mIQ71+0hBfGv4HX8nDZnPM4/90zWbPv+4gcf+7Gdxn2al/mb3qfqUP+yoe//6xWhS6AUR3G8slly+h1bF8mLryGG+dfSZ4nN9bFUkopVUEVDV5XAoOw51vcYoxpD7wevWKp2soYw/hO5/L55d/z4Ign2XBgHaPeGMyNH13Jjtyqzce5v/BX/vjBxVzz4aW0btCGhZcs5Zb+k0l21qrfDYocn9GSd373EZMHTeN/P81mxOun8Pb6WVgBnUxXKaXind7VGFTRSbJVZOV6cvjniieYsXI6gnB17xu4pf9kGqaVf++GiDB347tMXTKJfG8utw+8m+v7TYr50BU1admur7hryW38sH81nZqcyF8G3s1vO59XZ6dtUkqpssRDU2N543it5cjxu4pWYY/jVSMz+mrwqv125e3g0a/vZ/b6WTRKa8ykAVOY0OvaUqc32l+wjymLJ/LhL/+j97En8/To5zmxaalTh9ZqAQkw75c5PPb139h4cANdm/XgjsHTGNVhnDahK6VUmEQIXm3LOrCIVK1tqJI0eNUdP/y6mr99cRefbf+Etg3bc+eQ+zmr8++KAoSIMGfjbO5cfCv5vjxuH3QP1588sU7VcpXGClj876e3eXzpA2zN2UzvY09mypD7OK3NSA1gSilFAgSvYgdpC3QSkUXGmHQgSUTyIljOUmnwqns+3baI+z+/k/UH1tL3uFO4d9jDtG/UkSmLb2HeL3Ppc1w/nho1o87WcpXFZ/mYvWEWT3zzELvydjCw5RDuGHwfg3TeR6VUHZcwwcsY8yfgGqCJiHQ0xnQCnhORkZEtask0eNVNVsBi9oY3ePTr+9iTv5v0pHoExGLy4Glc2/dmreUqh8fvYdYPL/P08kfZV7CX09qM5I7B0+h7fP9YF00ppWIikYLXKqA/sCw0dpcxZq2I1MiIlBq86rZCXyHPf/9Pvt+7gjuH/I3OTbvEukgJpdBXyMw1M5i+4gkOug4wqsM47hg8jW7Na6SLplJKxY1ECl7LRGRAaNBUY0wSsFI71yuVOPK9ebzw/bP8+7unyPFk89tO53H7oHs0yCql6ox4CF4Vvef8M2PMnUC6MeYMYDbwfqQKqJSKvoyUTCYOuIPlV21g0oCpLN66kOGvncxN869ma/bmWBdPKaXqhIoGrynAfmAtcC0wD7g7WoVSSkVPw7RG3DF4Gsuv3sB1fW/hw5//j6Eze/H62pdiXTSllKr1ygxexpjmxpiuIhIQkedF5AIROR9YCjSrmSIqpaKhaXozpg17iG+uXMfQ1sO545Ob+XL7p7EullJK1Wrl1XhNB5qXsLwl8HTki6OUqmnHZhzP8+Nm0aFxJ/704SXa7KiiTkT477rX+MeyR0iE2VOUiqTyglcPEfms+EIRWQDoLVFK1RKZqQ149ex3EBEun3O+TrytomZ/wT4un3M+Exdew6Nf/5WZa56PdZGUqlHlBa+yuv7XzhmIlaqj2jfqyPPjZ7Hp0EZunH+VTrqtIm7eL3M47bV+fLF9MX877XFGtBvFvZ9NZt3+NbEumlI1przg9bMxZmzxhcaYMwFtj1Cqljm1zQjuH/44Czd/yKNf/zXWxVG1RJ4nl1sWXMNV7/+eVpltWHjJUv7U9yaeGf08jdIac+2Hl1HgzY91MZWqEeUN/T0J+MAYcyHwXXBZP2AQMD6aBVNKxcZVva7jxwPreGbF45zUrBvndrko1kVSCezrnV9w84I/sjtvJ5MGTOXWAVNJdtoNJs3rHcO/znyFC945kymLJzJ9zAsxLq1S0VdmjZeIbAR6AJ8B7YKPz4CewXVKqVrGGMODI55kYMuhTFp4Hd/vje7gxap2cvvd/PXzqfxu9miSHcm8f9Fi7hg8rSh0hQxtfRqTBkxl9oZZvL1+VoxKq1TNKXPkemOMkXJuOanINtWlI9crVfMOFO5nzJtD8Vs+Flz8FcdmHB/rIqkEsW7/Gm786Cp+zFrHFT2vYdqwh6ifXL/U7f0BP+e/cyZrfv2ehRd/zQlNOtdgaVVdkggj1y8xxvzZGNOm2AFTjDG/McbMBK6IVEGVUvGjWb3mvHrWO+R6c5nw/oW4/e5YF0lVg8/y8fLq/zB50Z+Zu/Fdcj05ET+HFbCYvuLvjHljKAddWcw65388OvLpMkMXQJIjiX+PfYU0ZxrXfHipftZUrVZejVcacBVwCdAeyAbSsQPbQuBfIrIq2oXUGi+lYmfeL3O46v3fc/5JFzN99AsYY2JdJFUJIsKCzR/wty/uYtOhn0lPqofLX0iSI4n+LQZzevsxnN7+TDo1ObFa7+227C3cvOCPLNv9Nb/tdB6PjnyGJulNK3WMjzd/xGVzzmNCr2t55DdPVbksSpUmHmq8KjRJdvAgydij1btEJDuC5SuXBi+lYuvJbx7msaX3c8+pD3Jjv1tjXRxVQd/v/Za/fj6Vb3Z9SacmXZh26kOMaHcG3+1ZxqIt81m0ZT4bDvwAQJsG7RjZfjSntz+TIa1PIy0prULnEBHe+OEVpn02GafDycMj/sF5XX5f5RB372d38J+Vz/DC+DcY3+ncKh1DqdIkVPCKJQ1eSsWWiHDNh5fywc//x2vnvMfp7cfEukiqDNtztvHI1/fy3o9v0azeMUwedA8Xd59AkuPoG9l35m7nky0L+GTrAr7YvhiX30V6UjpDW4+wa8M6jKFlZusSz7O/YB+3LbqRhZs/ZGjr4Tw9ekap21aU1/Jy1lu/YXP2Lyy6ZBltGrat1vGUCqfBq4I0eCkVewW+As55+3S2Zm/mw99/RuemXWJdJFVMjjubZ1Y8zgvf/wuHcXDdybdwY79byUjJrND+br+br3d+zqItH7Fo83y2524F4KRm3YNNkmM4+fgBJDmS+OiXudy26EYKvfncNfRvXN3nBhymvG7DFbMtewunzxpIpyYnMufCT466E1KpqtLgVUEavJSKDztztzPmzVPJTMlk3h8+p3Fak1gXSWF3nJ+55nme/OYhDrkPckHXS5gy+F5aZLaq8jFFhJ8P/mSHsC3zWb77a/wBP41SG9OlWVe+2fUVPY/pw/QxL3Ji05MieDW2uRvf5ZoPL+WGkycxbdhDET++qpsSLngZYxoQNuiqiBysfvHKp8FLqfixYvdSzps9mkGtTuWNc+eU2HylaoaI8NGmuTzwxd1szv6Foa2Hc++wh+lxTO+InyvXk8On2xaxaMt8vtuzjLM6/45JA6aS4kyJ+LlCJi/6M6+ufYFZ5/yPke1HR+08qu5ImOBljLkWuB9wAaEdREQ6RKyUZdDgpVR8eXPdq0xaeC1/7H0DD4x4ItbFqZNW7l3B/Z9P5ZtdX9G5yUlMG/YQI9uNrlV3nbr8Lsa+OYx9BXtZfOkyjstoEesiqQQXD8Grog3yfwG6iUg7EWkffNRI6FJKxZ8/dLuca/v+mRdWPcustS/Hujh1yvacbVw/7wrGvjmMXw79zOMj/8niy5ZzevsxtSp0AaQnpfOfsa/h8hXqxO2q1qho8NoEFEazIEqpxHLPqQ8xou0ZTFl8C8t2fRXr4tR6Oe5s7v/8TobO7Mn8Te8zacAUvrnyBy7reXWtbu7t3LQLD//mKb7a8RlPLX801sVRqtoq2tTYB3gZWAZ4QstF5OboFe0wbWpUKj7luLMZ+99hZLuzmX/xF7RuoLf+V5Tb7+aQO4uDLvuR5TpQ9PNBd9jPriwOubPYX/grPsvHhV0v5Y7B06rVcT7RiAh/nn817/30Fu+cP5/BrU6NdZFUgoqHpsaKBq/lwJfAWqAonYjIzEgVsiwavJSKX78c3MjY/w6jdYO2zL3wE+qnZMS6SHHB7Xezau+3LN+9lJ+y1gfD1eFQVegrKHXfRqmNaZLeNPhoRpP0pjRNb8a5J15I92N61eBVxI98bx6jZg2m0FfAJ5ctp2l6s1gXSSWgRApeX4vI4IiXrII0eCkV3xZvXcil/zuXXsf25YTGJyKhe3CC/7+Enhf9LXLEz+HrDIbm9Y6hbcP2tGnYzn40aBf3gS7LdYAVu5eyfNdSlu/+mtX7VuIL+ABo1aANzdKbh4WoYKBKOzJcNUlvSuO0JrW66bA6fvh1NWP/O4xhbUby6tnvRGzcMFV3JFLwehDYBrzPkU2NOpyEUgqA19a8yD+/fYKA2P+ODHZH71CHb4M54ufi60LPAxJgb/4eCnz5Rxy/aXpz2jRsaweyBnYgC/3cIrNVjQ6yKSJsyd7E8t12yFqxeyk/H/wJgBRnCr2O6Uv/loMZ0GIw/VoMrPSchap0L676N3ctuZV7hz3M9SdPjHVxVIJJpOC1pYTFOpyEUioqRISD7iy25Wxhe85WduRsZVvOVrbnbmV7zlZ25m3HH/AXbe80TlpktiqqHWvbsD3HZ7QkIyWD+ikZ1E+uT0ZKJvWSM+xlyRmkOlMrfBegz/Kxdv8qlu+yQ9ay3Us5UPgrYDcLntJiEP1bDqJ/i0H0OvbkCs9zqCpPRLj6g9+zcPM85l60mL7HnRLrIqkEkjDBK9Y0eCmlwvkDfvbk72J7zja25WxhR+7WopC2PWcbvxbuLfcYTuMkIyWT+sn1g+HMDmkZyZnUT6lfFM7WH/iB7/euwOV3AdC2YXtOaTGIAS0G07/lYDo1OVGbvGpYtvsQp88aiAMHiy79hgapDWNdJJUgEiJ4BUerby4im4ot7ykiayJXzNJp8FJKVUahr5D9BfvI9+VR4M2nwFdAvjePAl8BBb784LL8w8u8+eT78ikM/h3ap9BXQMfGnRnQcpBdq9VikA7iGSe+3f0NZ799OmNPOJsZ416vdWOYqeiIh+BVZg9OY8yFwFPAr8aYZGCCiKwIrn4F6BvJgiqlVCTUS65H20btY10MFUX9WgxkypD7ePDLe7hzySTuHfaINvGqhFBe/fidwMkiPFcPSwAAHd9JREFU0hu4EnjNGHNecJ3+eqGUUipmbux3K9f0uYmXV/+Hcf89regGB6XiWXnByykiewBEZDkwArjLGHMzh+dsVEoppWqcwzi4f/jjvHr2u+zJ382oWYN544eZJELfZVV3lRe88owxHUNPgiFsOHA20K2sHY0xLxljfjXG/BC2rLcx5htjzCpjzLfGmP7VKLtSSinFqA5jWXzpMk4+vj+3fnwd1827nFxPTqyLpVSJygte11OsSVFE8oAxwFXl7PtKcLtwjwF/DTZdTgs+V0opparluIwWvHXeB9w55H4++Pn/GPn6AL7bsyzWxVLqKGUGLxFZLSK/GGMeLbbcB/QsZ9/PgeIDrArQIPhzQ2B35YqrlFJKlczpcHJz/9uZc+EiBOGst0byzPLHiwb1VSoeVHTwmTNKWHZmFc43EXjcGLMD+DswtQrHUEoppUrVr8VAPrlkGeM6ncNDX03jonfHszdff89X8aHM4GWMud4YsxY40RizJuyxBajKGF7XA5NEpDUwCXixjHNfE+wH9u3+/furcCqllFJ1VcO0Rvxn7Gs8eca/+XbPMn7z+gAWbZkf62IpVfYAqsaYhkBj4GFgStiqvIrM02iMaQd8ICLdg89zgEYiIsYe7S5HRBqUcQhAB1BVSilVdRuzfuS6eZez/sBarulzE3cNfYDUpNRYF0vFQDwMoFpeH68cEdkqIn8QkW1hj6pOjr0bOC3482+An6t4HKWUUqpCOjftwrw/fM7Vva9nxvf/ZPxbw9l0SL9+VGxEbYIxY8ybwFLsZsqdxpirgT8BTxhjVgMPAddE6/xKKaVUSFpSGg+OeJJXz3qHXbk7OGPWIP677jUd80vVOJ0kO0ibGpVSqm7Yk7+LGz+6iq93fs65J17IYyOnk5labq8XVQvEfVOjUkopVdscn9GS2b+bx5TB9zF347ucPmsg3+z8Umu/VI3Q4KWUUqrOcTqcTBxwB/+78GMssThn9hmc9urJPLP8cXbkbot18apFRNh06Gdy3NmxLooqgTY1BmlTo1JK1U15nlz+76e3eXfDmyzb/TUAA1sO5fz/b+9Og+S6zvOO/9+79jYzmAFAEgRJLRRJyFwg0aSSSqyIimOb8iJZEh2TkSiK0RaVLCulUpayP1hOKlUuF+VyWU6lIloMKUcmI5mkZJWiVCSXNrsSiZRCAqAA7qCEhQRIEMQsmF7uPflwb880BjODWXrv51cFoKenp/v06YvuZ8459z2vv4Vfv+ydbClM9riFazNTm+b+/fdxz547+cmLewG4ePxVXH3ebq7cvpurtl/DlduvYefYxWSFBTrjxOmXePLEAZ488ThPnjjAc688y1g0wc6xi9hR2cmOsZ0LlycLUx1ty1L9MNWo4JVT8BIRkedeOciDB/4Hf73/r3jq5SeI/Ihfes3bePfrb+EXX31jX5ah2P/iPu5+9E7uP3AvM7Vprtq+m5uvvJXT9Tn2Hn+Ux47v4ZmXn8KRfd5PFqa4Mg9hV23fzVXnXcPrJq8g9NceSJxzHJk5lIWrl7KQ9UQetl46vVh7sxgUuWTiNczUTvH8zFESl5xxP8WguBDGdlR2cmFlJxeOXcSF+XUXjl3EVGFr28KZgtcaKXiJiEg3OefYc+z/cf/+e3nw8S9zfO4FJuItvP3yd/PuXTfzpp3/CM96t1pnvjHP1598kHv23MkPj/wfYj/m7ZffxPt3f4hrL3jTWUFltjbD/hf3ZUHs2B72HX+UAy8+xnwyD0Dsx1yx9UquOq8ZyK7h57ZdTTEscfDkMwsjWFm4OsBTJ55gtj6zcP9b4kkum7qCy7buyv6d2sVlU7u4ePyShX5qpA2Oz73A0enDHJ4+xNGZwxyZOczR6cMcmT7EkZnDPD9z5KxwFvsxOyo72X3+tfyrn/8Eb7xg2TyzJgpea6TgJSIivdJIG3z/p9/m/v338vWnvsrpxhwXj7+Kd++6mXftupnLt+7qWlsOnnyGL+z5C+77yV9y4vSLvGbLpbzvmg/y2z93K1PFreu6r0ba4OmXn2TfsUfYd3wPjx3fw75jj3Ji/qWF24ReSD2tL3x9QXkHl03t4vKtuxbC1eVTV7CtdF5bRqWSNOHFuWMcnjnE0enDWTibPszh6Z/xnee+xSvVk7z54rfy8es/xZsveeu6H1PBa40UvEREpB/M1mb4xtNf4/799/Ldn/4tqUu55rxruen1N/ObV/wW55UvaPtjNtIG33r2G9zz6J18+7lv4pvPjZf+Ordd82F+4ZIb2jry5pzj+dkj7D32KPuOPcpcfZbXTV3B5VNX8LqpKxiPJ9r2WOs1XT3FX+79PP/1x3/GC7PPs/v8a/n49Z/ibZe+Hd/z13QfCl5rpOAlIiL95tjs83zl8S/z1/vvY8+xH2MY55cvyNYmVS7iwnyNUuvapfPLO9a8luqFmaN8cd/d/Pe9n+fIzGEuKO/gvVd/gPdc/X52VHZ2+Nn1r2qjypf3f5H//PCf8OzJp7l08jI+dt0nuen1/4LIj1b9WQWvNVLwEhGRfvbESwf4n099lYOvPJOtWZo5xOHpQ8zVZ8+4nWce55XOPyOc7WhZUH7h2EUcPPkMd+/5HP/r6a/RSBu85ZJf5LbdH+KXX/trBF7Qo2fYf5I04etPfYXPPnQHe489wo7KhXzk2t/l1qs/QDmqLPszCl5rNCjBqxD61JOUJO3/PhURkc5yzjFdO7W4kHz6EEcW1i1li8mPTB86Y5F602Rhit++8lbed/UHee3k63rQ+sHhnOO7P/1bPvvQHfz9z77LlniSf/mGj/KBN36UrcVtZ9xWwWuN+j14FUKPShzie0YjTTkxU6P/e1VERHptaTg7On2YclTmVy79DYpBsdfNGzg/OvoD/vyhz/CNp79GMSjx3qtv5yPX/i4XjV8CKHitWb8GrzjwqMQBge/hnMPMcM4xX085NV8/9x2IiIhI2z3x0gH+/OHP8MCB+wB4166b+dh1n+QNO67qefDSlkEbEPrGVCliSynC97JTWZuntJoZxcinEKhrRUREeuHyrbv4s1+5k/97+2PcvvsjfO2JB3jLF67lPQ/cxN4X9va0bUoH6xB4xpZiyFQ5JvDPDFytnHOMFUO8Lm6DICIiIme6aPwS/uMNd/DwBx/nk//g9/i7n32f52ee72mbNNWYW22q0TejHAcUwiynrqVgm3OORuI4MVdraztFRERkY1LmuWBsouP7Q6421ajzUlfhGZTjgGKYFWZbzwtlZoRBFthmq41ONVFERETWqByVu7op93IUvJZhQCkOKEfrD1xLlSOfWiOhnvT/yKKIiIh0loLXEqXIpxwHbV2fNVGMeGmmqhITIiIiI07BKxf5HmOFrBZXO9e9mRkejvFiyCunVWJCRERklOmsxlwpDsgrQ7R9/tfMKIQ+hXBtm3iKiIjIcFLwatHJBXfOOcYLwULdr1ESeMboPWsREZGzKXh1STPUTRQ7XzG3X3hmTBRDtlZiJkqj87xFRERWouDVRWZG6GfbDA27cuSzrRIR5xX848CnFGmqVURERtvwJ4A+VI4Dao2UWrKxTbn7WRR4jBUCAu/sTF+JA6qNlCTV+Z0iIjKaNOLVA845Joohw7SjUHNacbIU4a/yxEZpqlVERGQpBa8eMDPMYKILO6R3w9JpxZVOUhilqVYREZHl6BOwR8yMOPQpNlJO15NeN2dDosDLz9RcX34vRT7VRkp9CKdaRUREVqMRrx7L1kMN1pyjb8aWfFpxoxX+J4qhSkyIiMjIUfDqA4O07qkcB2ytRETnmFZcjZnhGYwN0PMWERFpBwWvHjMzAj87E7CfxYHHtkpEJQ7yNWqbG68yM4qhTyHQISgiIqOjvz/tR0gpykpMVBv9te7JN2OsEBCHflv3sITs7M6xYkhtpsqgVJjwPSMOPDwz6klKrZFq83MREVkzBa8+4Vy2kfZLfRJCjGxasbXoaSf2sMQ5JooRL8/V2nrf7RQFHnH+p3kigXMOswDnHI3EUU1Sao2EetIHL56IiPQtBa8+0Qwh48WQk3P1nralFPmU4wCjs/tXkt9/FBilyGeu1h9nd/pmC2ErCjzMLA9ai33RvJxNFUPg+1TiLIg1i+PWGimNfkjRIiLSNxS8+oiZ5VvrpD0JIYXQoxKH+N7ZQaPTKnk1/14FlcjPQlYcegtV91v7YLW+WBrIsvvJRgrTtDkalo2IKYeJiIw2Ba8+NFYIKYQ+p2sJ8/Wk42uIIt+jUggIfW9hHVc3Q1fTRDHkpdnuTDl6ZgvTh+ca1Vqv1p/zPKNgHsXQB0IayeJomNaHiYiMHgWvPhV4xngxZKwQMF9POF1v//qhwDMqhYA48HsauJqPG/hGJQ6YqTY68xhAKQ4oBB6Bv75RrU09bsv9Br6H7xmlKJuWnKslHXu+IiLSfxS8+lRrGCiEPsUooJGkzNUT5mubGwXzLAs4xaj3gWupTm0gXgz9vBTG2VOD3db62pbjAM+MU/O9XdcnIiLdoeA1AJof1L5njBdCxuKA+UbK6Vqyrm13zKAcdfZMxc1aOLtztko7qleEvjFWCBemUfvt+QIUIx/Pg1fm6pp6FBEZcgpeA+SMUbAgWzfUSNOFtWCrLdzu5pmKm2FmeDjGCyGvnN74KJBnUCmEFMP+G9VbTuR7TJYjTs7VtABfRGSIKXgNqIVRMMtGdCpxQDXfcLvWUoS1kE+x+QO0H2RzerXaSJivr3/KsZyHzNb763dmRuDBVDni5bk6idKXiMhQUvAacK2jYHHgUQh9ktQxX0+I80Xk7a443w3OZaNe9UaNZI3tj4Ns66VmkdNBk432ZeHr5FxNxVhFRIbQYH5CybJa14KVIn9hlGsQRnyWarZ5fA0bafuesaUUsqUU4Q3gc21lZhgwWVrciFxERIaH3tmHVDs2su61ZjHS1pMBzvg+MBYHbC1HRHl5iEF/zrD4HLYUw7z+l4iIDAtNNUrfW66qfTH0qRT6/2SBjWoWdB0vhnieMataXyIiQ0EjXjIQJkrZlGPoG1PliPFiOLShq6n53CpxwHhBvyOJiAyDjgUvM7vLzI6Z2b4l13/czB43s8fM7I879fgyPLIz/jy2liOmyjHBAK9d26hiFLBlDevdRESkv3VyxOtu4MbWK8zsrcA7gGucc1cCd3Tw8WXIDPLJAu0QBR5T5YgRffoiIkOhY8HLOfc94MSSqz8K/JFzrprf5linHl+Gz6gGrqZs5C+bavVHvC9ERAZVt9d4XQ682cx+YGbfNbPru/z4IgPNzPAtC1/BABXFFRGRTLeDVwBMAv8Q+DfAl2yFYQwz+7CZPWxmDx8/frybbRTpa1mpkKzQarOMhoiIDIZuv2sfAh5wmR8CKbBtuRs65z7nnLvOOXfd9u3bu9pIkX63UOurlNX60uCXiMhg6PY56l8B/inwHTO7HIiAF7vcBpGh0FrrC0LS1FFPUxqJo56k1BNHOoDbRYmIDLOOBS8zuxe4AdhmZoeAPwDuAu7KS0zUgNvcIG4kKNInWmfqPc+IzCPywSz7r526LIQ1EkcjSamnThtwi4j0UMeCl3PulhW+9d5OPabIqFu6ZNIzI/LPDGPOOer5qFgjzUbGFMZERLpD5bBFhtzSMGZmhD6Evn9GGJueb3C6nvSiiSIiI0PBS2QELRfGxoshDphX+BIR6Ridiy4iQDbqNV4IiAO9LYiIdIreYUUEWBwFmyiGqg8mItIhencVkQWt9cFCX8XBhk0lDthWiSlHPnp1RXpDwUtEzrAYvrQt0TCJAo9yHOAZVAoh28diKvnX3eabUY4DKrGWGcvo0VEvImcxM3COyVLEibla18pNRL7HeDGg1kiZrjZQlb/28AwmCiHOuYVgbWaUIp9S5DPfSJmrNmh0+HUuhB7FMCAKvIW2eAan5hsdfVyRfqLgJSLLyj6g8/A1W6XT2asSB5TjAOcchdAnDnymq3Xm62lnH3gEjBdDzJY/mxWgEHgUw5hqI2GumlBL2tfnoW8UQp9C6OPluy20PnYxCnAOpqsKXzIaFLxEZEVmhodjshzx8mytI+HL94yJYkiYL+hfDAeOiWJEIUyYPt0g0fDXhhSjLMSuptnnke8Rl33qScpcrbHh0OsZFEKfYugT+N5ZI21LleKAFJhV+JIRoOAlIqsyM3yyNV8vz9ZoZ/wphB7jhXDFx4UsDGytRMxUG8zVVGNsPQLPGFvHOqpmnweeMVGMqMSOuVqD07VkTa97FHgUQ5848Bb2Em2939VU8tFOvcYy7BS8ROSczIzAy8PXXG3z9weMFUKKkX/GaMhKj03z9qHPqfk69USjX2sxUVw+1J5Ls889y/q9HAecriXM1RpnjXr6nlHMpxJ9z845urWasUJI6lTEV4abgpeIrElzq6EtpZCTc/UN30/gGROlEH8DH86+Z0yWIk7XE2bmG20dfRs2Y4WAYJP12BYDWMtC/HrCXC0h8LPRrdaF8q0/sxHNIr7OOaoNre2T4aTgJSJrZmbEgc9EEV45vf7wVYr8hRICG/mAXliQ3Vx8P1/XB/QyosCjFLX37X1hIX7oU4wW9/hs/V47HsM5x0QxC/ftXOQv0i9Ux0tE1q0Q+owX1v7B7lk2UjZWCDGzTX9QN8sQbClFTBTDntSi6letpSM6ofW1a1fgWu7+VcRXhpWCl4hsSDFaWwHMyPfYWo7bvg1R8wM6Djy2VWKK0epn7o2KlUpHDBIV8ZVhpuAlIhtWjgPKqwSeShwwWY46GgSa9zteCJkqj/YHdbN0xCCHriYzw4DJUoQ/wq+pDB8FLxHZlEohpLQkfPmeMVWOKG9iPdd6tJZBmCpHI7kVzXpLRwyCbFo6C1/KXjIsFLxEZNPGCiGFMAtfhdBja49Gnprrx8r5ZtClEdoMeqOlI/pdcz3fZFnhS4bDcP16JCI90SwDUAg94uDctbm6obUG1Vxeg2pYi9+PxZsvHdHPOlnEV6Tbhvd/qoh0TWuV+dave6m1BlU58tleiRkrBAv1w4ZFFHiUhmyKcTlZEd+sjttwvYIyahS8RKQt2lEmolOabSuGPlsrWQmKYViE3+nSEf3GzAj8rACvyKAa/l+TRERyrSUoCmFMtZEwV00GtlDneGHwS0es12aL+EJ28kfgWf6vR+BnhVun5xs0OrETvEgLBS8RGTmtU6Nx2aeepMxVG8wPUBX8YugTh6Nbu6wQ+jgHp+ZXDl8L4cr3CFrCVmtQba5HdM4xVR6c7ajMsj1Pm2U3mgF84TJGkjoaaXrW/prSWwpeIjKyWstQTJQiKmnKbDXhdJ9v0ux7xtg6dg4YVsXIx+GYryf43mK4CnzDa5n6Xm1ro6V7TBZDn0LgM12tM1/vbhA3sucU+h5m4GFnBCpawtZ6pM7RSFIaiaOROuppSpK4vg+Xrcyy9Zo4aLbc5X8N0vMABS8RkTMW4o8XQyqFgLlqg7l60pdnQm4Z0tIRG1GKgoV9KVcKWOsJKtltHRPFiGKYMj1f7/j0o1nzefh4+ejbYls2zzMj9D1C/8z7TNKUeh7GGklKI3UkPR4e84zFEO3nQdrz8Na4JrPZd0tDmcu/6IdpeQUvEZFcawCr5KUo5hsp1XrSN5txD3vpiM1o52bdAKGfFeSdqyXMVts//ei1BK5u7YHZyvc8PHPEgNlieM2CWDZNmaSO1GXXOwcpri2/jHgGgefht4Sr5khl00YC6MII5sJfC9/ZfKPbRMFLRGQFZkYh8CiGWW2yaiNlvp5Qa6Q9md4YldIR/aL5IV6KfAqhz8x8vS3rAD0zyrFPMV+j18tRmOVGB0PfCDyH2cprCJ1bPpCl+dcL38+vb57Q0AxbawlY/TA61Qn6HywisorWNUDZ2ZBZCKs1UuYbCdV6d0KYtZSOGNYPpH5lZng4JkoRxUbCqfnGhqbkfM8oR1mh4eb99qtztc3M8FuGlZaWNFnu51c6dvu5HzpBwUtEZI1aQ1gUeMShjys4aklKtZ4FsXauCTPLRwrMKEb+yJWO6CeL04/ZlljrmX4MvGwbqzjo/8C1UWt5TsP4vDdCwUtEZAPOCGF+tlXSmAuoJynz9ZRqI1nTafx+XuLAN1u4vFrZA+mt5mtQjgMKoc/0fH3F9X+h3wxc/bGNlvQHBS8RkU1qDWGh7xEFPhBSz6cj64nD98jDlbcYtpacqbWWsgfSPzzL9o+sNhKmW6Yfo8CjHAVEgdf2MxRl8Cl4iYi0UesHbOAbY8Fi6YdzfQjrw3mwtBbi3VqOOF1LCAOP0FfgkpUpeImIdMionKU16hbOfowDBS45JxWDERERaRMFLjkXBS8RERGRLlHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLrHmhp79zMyOA8+tcpNtwItdas6oUd92lvq3c9S3naX+7Rz1bWd1o39f5Zzbvtw3BiJ4nYuZPeycu67X7RhG6tvOUv92jvq2s9S/naO+7axe96+mGkVERES6RMFLREREpEuGJXh9rtcNGGLq285S/3aO+raz1L+do77trJ7271Cs8RIREREZBMMy4iUiIiLS9wY6eJnZjWb2uJk9ZWb/vtftGTZmdtDM9prZI2b2cK/bM8jM7C4zO2Zm+1qumzKzb5rZk/m/k71s4yBboX8/bWaH8+P3ETP71V62cVCZ2cVm9m0z229mj5nZJ/Lrdfy2wSr9q+N3k8ysYGY/NLNH8779w/z6nh67AzvVaGY+8ATwS8Ah4CHgFufcT3rasCFiZgeB65xzqiezSWb2T4AZ4AvOuavy6/4YOOGc+6P8F4dJ59y/62U7B9UK/ftpYMY5d0cv2zbozGwHsMM592MzGwN+BPwm8H50/G7aKv37z9HxuylmZkDZOTdjZiHwd8AngHfRw2N3kEe83gQ85Zx7xjlXA+4D3tHjNoksyzn3PeDEkqvfAdyTX76H7M1WNmCF/pU2cM4ddc79OL88DewHdqLjty1W6V/ZJJeZyb8M8z+OHh+7gxy8dgI/a/n6EDpY280B/9vMfmRmH+51Y4bQ+c65o5C9+QLn9bg9w+h3zGxPPhWpqbBNMrNXA28EfoCO37Zb0r+g43fTzMw3s0eAY8A3nXM9P3YHOXjZMtcN5rxp//rHzrlrgbcBH8unc0QGxX8BLgXeABwFPtPb5gw2M6sA9wP/2jl3qtftGTbL9K+O3zZwziXOuTcAFwFvMrOret2mQQ5eh4CLW76+CDjSo7YMJefckfzfY8CDZNO70j4v5Os7mus8jvW4PUPFOfdC/qabAnei43fD8vUx9wNfdM49kF+t47dNlutfHb/t5Zw7CXwHuJEeH7uDHLweAi4zs9eYWQTcDPxNj9s0NMysnC/0xMzKwC8D+1b/KVmnvwFuyy/fBny1h20ZOs031tw70fG7IfkC5c8D+51zf9LyLR2/bbBS/+r43Twz225mW/LLReCfAQfo8bE7sGc1AuSn1/4p4AN3Oef+U4+bNDTM7LVko1wAAfBX6t+NM7N7gRuAbcALwB8AXwG+BFwC/BT4LeecFohvwAr9ewPZNI0DDgIfaa7rkLUzs18Avg/sBdL86t8jW4ek43eTVunfW9Dxuylmdg3Z4nmfbKDpS865/2BmW+nhsTvQwUtERERkkAzyVKOIiIjIQFHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLlHwEpGRYGafNrNP9bodIjLaFLxEREREukTBS0SGlpn9vpk9bmbfAq7Ir/uQmT1kZo+a2f1mVjKzMTN7Nt+6BTMbN7ODza9FRNpFwUtEhpKZ/TzZVmJvBN4FXJ9/6wHn3PXOud3AfuADzrlpsn3cfi2/zc3A/c65endbLSLDTsFLRIbVm4EHnXNzzrlTLO7lepWZfd/M9gLvAa7Mr/8L4Pb88u3Af+tqa0VkJCh4icgwW25PtLuB33HOXQ38IVAAcM79PfBqM3sL4DvntCmxiLSdgpeIDKvvAe80s6KZjQG/kV8/BhzN12+9Z8nPfAG4F412iUiHaJNsERlaZvb7wPuA54BDwE+AWeDf5tftBcacc+/Pb38B8Cywwzl3shdtFpHhpuAlIpIzs5uAdzjnbu11W0RkOAW9boCISD8ws88CbwN+tddtEZHhpREvERERkS7R4noRERGRLlHwEhEREekSBS8RERGRLlHwEhEREekSBS8RERGRLlHwEhEREemS/w+QKUBGTN+cBQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y1 = Lille_m\n",
    "y2 = Lille_m + Lille_sd\n",
    "y2 = np.squeeze(y2.values)\n",
    "y3 = Lille_m - Lille_sd\n",
    "y3 = np.squeeze(y3.values)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "ax = plt.subplot()\n",
    "ax.set_ylabel('t2m (Celsius)')\n",
    "ax.set_xlabel('day')\n",
    "ax.plot(Lille_m.day, y1, color='green', label='t2m mean, shading: +/- SD')\n",
    "ax.plot(Lille_m.day, y2, color='white')\n",
    "ax.plot(Lille_m.day, y3, color='white')\n",
    "ax.fill_between(Lille_m.day, y2, y3, alpha=0.1)\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(handles, labels)\n",
    "ax.set_title('t2m climatology for Sep from 1979 to 2019')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the typical range of maximum 2m temperature values for September 15?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now look more closely at the probability distribution of maximum temperatures for 15 September in this time period. To do this, we will first select only the max daily temperature for 15 September, for each year in the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_max = Lille_max.dropna('time', how='all')\n",
    "Lille_15 = Lille_max[14::30]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will then plot the histogram of this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 2.,  3., 11.,  6.,  7.,  5.,  1.,  2.,  3.,  1.]),\n",
       " array([12.53629112, 13.94420452, 15.35211792, 16.76003132, 18.16794472,\n",
       "        19.57585812, 20.98377151, 22.39168491, 23.79959831, 25.20751171,\n",
       "        26.61542511]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEWCAYAAAB/tMx4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQZElEQVR4nO3dfbBcdX3H8fdHYqoYMWguFoF4kYKtZaaC0YJaraAdFBX/cFqottA6zVSnFJk6NlZb2+nUCepgmbHVppWKyuADolIZK5ZKHR2gJhSFGGx8CIJEEnViwCeM+faPPbHr5j7l7ubu/uD9mrmze86eu+dzN/d88ru/3T2bqkKS1J6HjDuAJGlxLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4GpKks1JfnPcOaRJYIFroiTZluS5A+vOS/JZgKr61aq6fp77mE5SSZYdxKjS2Fng0gHyPwZNCgtcTekfoSd5WpKNSXYnuSfJxd1mn+kudyW5L8mpSR6S5A1J7kiyI8l7kjyq735/v7vtO0n+cmA/f53kyiTvS7IbOK/b9w1JdiXZnuTtSZb33V8leVWSrUnuTfK3SY7rvmd3kg/2by8thgWull0CXFJVhwHHAR/s1j+ru1xZVSuq6gbgvO7rOcATgBXA2wGSPAn4R+BlwJHAo4CjBvZ1FnAlsBK4HPgpcCGwCjgVOB141cD3nAE8BTgFeC2wodvHMcCJwDlD/OySBa6J9NFuZLsryS565TqTnwC/lGRVVd1XVTfOcZ8vAy6uqq9V1X3A64Czu+mQlwL/VlWfrar7gb8CBk8SdENVfbSq9lbVD6tqU1XdWFV7qmob8E/Aswe+56Kq2l1Vm4HbgGu7/X8P+ARw0sIfEml/Frgm0UuqauW+L/Yf2e7zCuAE4PYkn0/ywjnu83HAHX3LdwDLgMd2t92574aq+gHwnYHvv7N/IckJST6e5FvdtMqb6I3G+93Td/2HMyyvmCOvNC8LXM2qqq1VdQ5wBHARcGWSR7D/6BngbuDxfcurgT30SnU7cPS+G5I8HHjM4O4Glt8B3A4c303h/AWQxf800oGzwNWsJC9PMlVVe4Fd3eqfAjuBvfTmuve5ArgwybFJVtAbMX+gqvbQm9t+UZKnd08s/g3zl/Ejgd3AfUl+GXjlyH4waYEscLXsDGBzkvvoPaF5dlX9qJsC+Tvgc908+inApcB76b1C5evAj4DzAbo56vOB99Mbjd8L7AB+PMe+XwP8brftPwMfGP2PJ80tfqCD9PO6EfouetMjXx93Hmk2jsAlIMmLkhzazaG/FbgV2DbeVNLcLHCp5yx6T3TeDRxPbzrGP0810ZxCkaRGOQKXpEYt6Ul5Vq1aVdPT00u5S0lq3qZNm75dVVOD65e0wKenp9m4ceNS7lKSmpfkjpnWO4UiSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNWtJ3YurATK+7Ziz73bb+zLHsV9KBcQQuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckho1b4EnuTTJjiS39a17dJJPJdnaXR5+cGNKkgYtZAT+buCMgXXrgOuq6njgum5ZkrSE5i3wqvoM8N2B1WcBl3XXLwNeMuJckqR5LHYO/LFVtR2guzxidJEkSQtx0J/ETLI2ycYkG3fu3HmwdydJDxqLLfB7khwJ0F3umG3DqtpQVWuqas3U1NQidydJGrTYAr8aOLe7fi7wsdHEkSQt1EJeRngFcAPwxCR3JXkFsB54XpKtwPO6ZUnSElo23wZVdc4sN50+4iySpAPgOzElqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckho1VIEnuTDJ5iS3JbkiycNGFUySNLdFF3iSo4A/BdZU1YnAIcDZowomSZrbsFMoy4CHJ1kGHArcPXwkSdJCLLrAq+qbwFuBbwDbge9V1bWD2yVZm2Rjko07d+5cfFJJ0s8ZZgrlcOAs4FjgccAjkrx8cLuq2lBVa6pqzdTU1OKTSpJ+zjBTKM8Fvl5VO6vqJ8BVwNNHE0uSNJ9hCvwbwClJDk0S4HRgy2hiSZLmM8wc+E3AlcDNwK3dfW0YUS5J0jyWDfPNVfVG4I0jyiJJOgC+E1OSGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRQ52NUBq16XXXjGW/29afOZb9SsNwBC5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjhirwJCuTXJnk9iRbkpw6qmCSpLkN+4EOlwD/XlUvTbIcOHQEmSRJC7DoAk9yGPAs4DyAqrofuH80sSRJ8xlmBP4EYCfwr0l+DdgEXFBV3+/fKMlaYC3A6tWrh9idlsq4PtZM0oEZZg58GXAy8I6qOgn4PrBucKOq2lBVa6pqzdTU1BC7kyT1G6bA7wLuqqqbuuUr6RW6JGkJLLrAq+pbwJ1JntitOh340khSSZLmNeyrUM4HLu9egfI14A+GjyRJWoihCryqbgHWjCiLJOkA+E5MSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRg17OlnpAWGcHyO3bf2ZY9u32uYIXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1augCT3JIkv9J8vFRBJIkLcwoRuAXAFtGcD+SpAMwVIEnORo4E/iX0cSRJC3UsCPwvwdeC+ydbYMka5NsTLJx586dQ+5OkrTPogs8yQuBHVW1aa7tqmpDVa2pqjVTU1OL3Z0kacAwI/BnAC9Osg14P3BakveNJJUkaV6LLvCqel1VHV1V08DZwH9W1ctHlkySNCdfBy5JjVo2ijupquuB60dxX5KkhXEELkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1KiRnI3wgW563TXjjiA9oIzrmNq2/syx7PdgcQQuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoxZd4EmOSfLpJFuSbE5ywSiDSZLmNswHOuwB/qyqbk7ySGBTkk9V1ZdGlE2SNIdFj8CrantV3dxdvxfYAhw1qmCSpLmN5CPVkkwDJwE3zXDbWmAtwOrVqxe9Dz/WTBotj6n2Df0kZpIVwIeBV1fV7sHbq2pDVa2pqjVTU1PD7k6S1BmqwJM8lF55X15VV40mkiRpIYZ5FUqAdwFbquri0UWSJC3EMCPwZwC/B5yW5Jbu6wUjyiVJmsein8Ssqs8CGWEWSdIB8J2YktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjRrJR6pJWjw/2mzpjPOx3rb+zJHfpyNwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGjVUgSc5I8mXk3wlybpRhZIkzW/RBZ7kEOAfgOcDTwLOSfKkUQWTJM1tmBH404CvVNXXqup+4P3AWaOJJUmazzAfqXYUcGff8l3Arw9ulGQtsLZbvC/Jl/tuXgV8e4gMS8msB0dLWaGtvGY9OBaVNRcNtc/Hz7RymALPDOtqvxVVG4ANM95BsrGq1gyRYcmY9eBoKSu0ldesB8ckZR1mCuUu4Ji+5aOBu4eLI0laqGEK/PPA8UmOTbIcOBu4ejSxJEnzWfQUSlXtSfInwCeBQ4BLq2rzAd7NjFMrE8qsB0dLWaGtvGY9OCYma6r2m7aWJDXAd2JKUqMscElq1JIUeJJLk+xIclvfurckuT3JF5N8JMnKpciyEDPl7bvtNUkqyapxZBs0W9Yk53enOdic5M3jytdvlt+DJye5McktSTYmedo4M+6T5Jgkn06ypXsML+jWPzrJp5Js7S4Pn+CsE3eMzZa17/ZJO75mzTsRx1hVHfQv4FnAycBtfet+C1jWXb8IuGgpsiw2b7f+GHpP2t4BrBp3zjke2+cA/wH8Qrd8xLhzzpH1WuD53fUXANePO2eX5Ujg5O76I4H/pXfKiDcD67r16ybh93aOrBN3jM2WtVuexONrtsd2Io6xJRmBV9VngO8OrLu2qvZ0izfSex35RJgpb+dtwGuZ4Q1L4zJL1lcC66vqx902O5Y82AxmyVrAYd31RzEh7yWoqu1VdXN3/V5gC713H58FXNZtdhnwkvEk/H+zZZ3EY2yOxxUm8/iaLe9EHGOTMgf+h8Anxh1iLkleDHyzqr4w7iwLcALwG0luSvJfSZ467kBzeDXwliR3Am8FXjfmPPtJMg2cBNwEPLaqtkPv4AaOGF+y/Q1k7Tdxx1h/1haOr4HHdiKOsWHeSj8SSV4P7AEuH3eW2SQ5FHg9vT9JW7AMOBw4BXgq8MEkT6jub70J80rgwqr6cJLfBt4FPHfMmX4myQrgw8Crq2p3MtMZJCbDYNa+9RN3jPVnpZdtoo+vGX4PJuIYG+sIPMm5wAuBl01ouexzHHAs8IUk2+j9KXpzkl8ca6rZ3QVcVT3/DeyldwKeSXQucFV3/UP0znI5EZI8lN5Be3lV7ct4T5Iju9uPBCZiemqWrBN5jM2QdaKPr1ke24k4xsZW4EnOAP4ceHFV/WBcORaiqm6tqiOqarqqpun9451cVd8ac7TZfBQ4DSDJCcByJvdMb3cDz+6unwZsHWOWn0lvqP0uYEtVXdx309X0/tOhu/zYUmcbNFvWSTzGZso6ycfXHL8Hk3GMLdEzuVcA24Gf0PvHeQXwFXqno72l+3rnOJ7FXWjegdu3MTnPks/02C4H3gfcBtwMnDbunHNkfSawCfgCvbnFp4w7Z5f1mfSeTPti3+/oC4DHANfR+4/mOuDRE5x14o6x2bIObDNJx9dsj+1EHGO+lV6SGjUpr0KRJB0gC1ySGmWBS1KjLHBJapQFLkmNssD1gJdkZZJXddefnOSG7gxyX0zyO+POJy2WLyPUA153DouPV9WJ3Zsuqqq2Jnkcvdeg/0pV7RpnRmkxHIHrwWA9cFySW4A/qqqtAFV1N723wk8BJNmW5E3dCH1jkpOTfDLJV5P88fjiSzMb+8mspCWwDjixqp7cv7L78IjlwFf7Vt9ZVacmeRvwbuAZwMOAzcA7lyautDAWuB6UuhNRvRc4t6r29t10dXd5K7CieueAvjfJj5KsdKpFk8QpFD3oJDkMuAZ4Q1XdOHDzj7vLvX3X9y074NFEscD1YHAvvY/DIsly4CPAe6rqQ2NNJQ3JEYUe8KrqO0k+132Y8iPonW/6MUnO6zY5r6puGVtAaZF8GaEkNcopFElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGvV/eaWUTo4aIF8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Lille_15.plot.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the range of maximum temperatures for 15 September in the period from 1979 to 2019. Has the temperature in this period ever exceeded that of 15 September 2020?\n",
    "\n",
    "The histogram shows the distribution of maximum temperature of one day in each year of the time series, which corresponds to 41 samples. In order to increase the number of samples, let's plot the histogram of maximum temperatures on 15 September, plus or minus three days. This would increase our number of samples by a factor of seven.\n",
    "\n",
    "To do this, we first need to produce an index that takes 15 Sep, plus or minus three days, from every year in the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = np.arange(41)\n",
    "days_in_sep = np.arange(11,18)\n",
    "index = np.zeros(287)\n",
    "for i in years:\n",
    "    index[i*7:(i*7)+7] = days_in_sep + (i*30)\n",
    "index = index.astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we apply this index to filter the array of max daily temperature from 1979 to 2019: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lille_7days = Lille_max.values[index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the histogram of maximum daily temperatures in days 12-18 September from 1979-2019:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFNCAYAAADRi2EuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debgcVZn48e8LQYmyBgKyRx1kRFTUuA0qCqIoCuiIIz8XUJRxQXFcIzoC4oy4L+OKyhAVQcQFBBeQIe4iYRMQFMQISIQgIIsCIu/vj3OuFJfuvn1zq+smN9/P8/TTtZ56u7pO99unTlVHZiJJkqTRW226A5AkSVpVmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvFYBEbEkIp7aZ94TI+LXXce0IonifyPi+oj4xXTHo3ZExAsj4pSWynpyRFzZRlkaXkRsExHnRMRNEfG66Y5HaoOJ1youM3+UmdtMtFxEHBIRX+oipmnwBGAXYPPMfMx0BzMmIuZFREbErOmOpQ0RsSgiXt7V9jLz6Mx8Wlfba1tEHBYR50fEHRFxyLh5u0XEjyPihoj4Y0R8NiLWHlDWJhFxYkRcVY+peePmz4mIr0TEtfVxdESsM6C8gyLidxFxc0RcGRFfmeLL7ectwKLMXDszPzaibfQ1wXsQEfH2iLg8Im6MiGOb+ywiLqz7Z+xxR0R8qzH/2RFxQZ3304jYdkAcD4qIEyJiWURcFxHfi4htxi3zH/VY+HNEHBkR927MOyAiFkfEbRFxVI/yXx4Rl9ZYvhsRmy7fHtMwTLw07VaAxGIrYElm3jLNcay0ungPV4DjpGuXUhKPk3vMWxd4N7Ap8GBgc+D9A8q6E/gu8K995r8bWB94APBAYGPgkF4LRsQ+wIuBp2bmWsB84LTBL2W5bQVc2G9mRKw+ou2OGfQevISyH3agvA+zgf8Zm5mZD8nMteo+Whu4HPhqjXtr4GjglcB6wLeAEwcc4+sBJwLbUN6bXwAnjM2MiKcDC4CdgXmU9/HQxvpXUd7jI8cXHBE7Av8N7AHMAX4HHNMnDrUhM33M8AewBHgT8Evgz8BXgDXrvCcDVzaWfSvwB+Am4NeUirwrcDvwN+Bm4Ly67KaUD4PrKB9Qr2iUMxtYCFwPXET58LpyXExvrTHdBsyifHD8tm77V8BzGsvvC/wE+DBwA3AZ8C91+hXANcA+A/ZBz1iB/YBbgb/X13Zoj3UntW1gN+Ac4MY6/5DGvH+r669Tx58B/BGY22O7lwNZ47oZeHyd/rK6T68Hvgds1VgngVcDl9T9eBjli/RnNZ7jgHs133vgIODa+p68sFHWvYEP1DiuBj4NzB637ltr/F+kfHGfBCyrsZ1EaUUE+K+6j2+tr+XjlC+IBGY1trkIeHmP/X4d5Yujb0x93rcfj9s3r6z75nrgE0D0WXc2cFRd7lfAm7n78dvzWK3xXQc8tLHsRsBfgbnAhnW/3FCX+xGw2gT190s0jqE+yzwXOH+Iz4JZdT/MGzf9O8CrG+OvAb7Xp4yPAx8ZsI11gc8DSymfJe8GVh/3nv4P5bPoYmDnPuX837hj5kH1PfkU8G3gFuCplMRzUd2nFwK7N8o4CvhkfX03123fD/hIfW8vBh4xxH67x3sAHA+8uTH+LzXW+/RYf8e6/fvW8QOAkxvzV6vHSM990aO8OfV93KCOfxn478b8nYE/9ljv3cBR46Z9APhEY3zTWvYDh4nFx+QftnitOp5PSaDuDzyM8gF4N7Xp+gDg0Zm5NvB0SkvQdym/iL6S5Rfcw+sqx1C+fDcFngf8d0TsXOcdzF2/vHYBXtQjpr0pScp6mXkH5YvsiZQP7kOBL0XEJo3lH0tJ1DagfNAcCzwa+Kda/scjYq0+r79nrJn5ecqX8c/qazu4z/qT2fYtlF/D69XX96qI2BMgM79CSYI+FhEbUL6gXp6Zy3ps80n1eb0a289qOQdRvmjnUr64x/863RV4FPA4SsJ7BPBCYAtgO8p+H3M/SjKwGbAPcETjFMZ7KV9229fXuRnwznHrzqG0SuxP+fL43zq+JeWL5OP1db+9xnpAfS0H9Hi9vTyWkqhuREneJoppIs+ivG8Pp9SJp/dZ7mBKwvrAusw+4+b3PFYz8zbKsdE83vcGvl/f4zdSjsO5lJaLgyhfclP1JAa0DA3hE8CzImL9iFif0jL2nT7L/hx4SUS8OSLm92h1WgjcQXl/HgE8DWieYh57Tzek7OevR8Sc8RvJzJ24+zHzmzrr/1GOhbWBMyitRadQjpHXAkePOw33fOAddXu3Uerf2XX8eOBDg3bMAFEfzfF7A1v3WHYf4Pi8q1W917pBqZ/DeBIlsfpTHX8IcF5j/nnAxvUzZiK9YmESsWiypjvz8zH6B6Ul40WN8fcBn67DT6b+kqd8UF5D+RW5xrgyDgG+1BjfgvJrdO3GtPdQf01RPlif3pj3cu7Z4vWyCeI+F9ijDu8LXNKY91DKF9bGjWl/ArbvUc5Ese5Lo2Wkx/rLve067yPAhxvj61FabM4HPjNgu/O4Z4vQd4D9GuOrAX+htnrV5XdozD8LeGtj/IPU1or63t9B/RVepx0H/Cflw/cWGr96gccDv2usezu15bRP/NsD1zfGF1Fbswa8vn8sU/f75Y15A2Pq876Nb/F6wrjXuqDPupcBuzbG96dx/E5wrD6W0tK5Wh1fDDy/Dr+LcoronyZRfwe2eFF+2FwPPGiIsvq1eG0KfJ9ySvJO4FRqy2ifcl5Yl7+FcuwvqNM3piQ3sxvL7g2c3nhPrqLR0kg5bfbiPtsZf8wcBXyhMf5ESovrao1px4ztr7r8ZxvzXgtc1Bh/KHDD8rwHlM+039TjeF1Ki3pSW6Yby92H0tr85Ma0f6777snAvSh17k7gbUPEsjmlJXHvxrTfjjte1+jzPvdq8dqZ0uL9MEpL72dqLHtPFIuP5XvY4rXq+GNj+C/APVqGMvNS4PWUJOua2lm0XyfLTYHrMvOmxrTfU1ogxuZf0ZjXHO45LSJeEhHn1g7DN1B+cW3YWOTqxvBfa8zjp/Vq8Zoo1mEMve2IeGxEnF47wv6Z0qL2j9eRmTdQ+npsR0mEJmMr4KONfXQdJSFpvpbxcQ3aR9fn3fu2/Z6yv+ZSvjDOamzru3X6mGWZeevYSETcJyI+ExG/j4gbgR8C602xH07zGBkmpolMWA+q8cfv75szBx2rmXkG5Ut1x4j4Z8oPmhPrqu+nnOo+JSIui4gFk4j9HiLicZQW2OdlbRGKcqXyWIfuYVvBvkpJItYG1qF8kfe9mCbLhQtPpfyIeCXwrtrPaCvKl/7Sxr75DKU1aswfsn7jV2PH3LCa78umwBWZeee48pa3PkzGkZQkbxGltfH0On381a/PpdTTH4xNyMyLKa1gH6eckt2Qcsr6SoBxnfK3HFsvIuZSWvc+mZnNlu6bKe/bmLHh5mdeT5l5GqXl8WuUfbekrudVvCNi4qW7ycwvZ+YTKB+gSTm1A/c8HXIVMCfufiXVlpRfYlA+TDZvzNui1+bGBiJiK+CzlFOdG2TmesAF3L0JfHlNFGvbvkz5ot0iM9el9EP6x+uIiO0p/bSOAQZdqdXrFNQVwL9n5nqNx+zM/Olyxrp+RNy3Mb4lZX9dS/lSekhjO+tm6SjcL743Ujr/PjYz1+GuU6XRZ/mxhO8+jWn3G7dMc51hYmrLUu5+zDa//IY5VhdSTje+mHKK6VaAzLwpM9+YmQ8Ang28oXF6flIi4hGU4+xl9cuTuo0fZe3UnZkPGbK4h1NaX2/JzJspx+wzJ1opM/+WmV+lnIbfjnJ83gZs2HiP1hkXx2YR0dxXY8fcsJrHxFXAFhHR/C4bZd2+K4jMOzPz4Mycl5mbU5KvP/TY9j6UVroct/7xmbldZm5ASXy2As6s89ZqPC4HqKeATwFOzMz/GreNCynv4ZiHA1fnXaciJ3otn8jMrTNzI0oCNotyTGsETLz0D1HumbNTvQz5VsqX3N/r7KuBeWMfcJl5BfBT4D0RsWZEPIzSUf3ouvxxwNtqn5HNKF9Sg9yX8oG6rMbyUlrqYzBErG1bm9LCdmtEPIbSJwWAiFiT0pJwEPBSypfQq/uUs4zS5P+AxrRPU/brQ2p560bEXlOM99CIuFdEPJHSB+qrtQXhs8CHI2Kjuq3NaqtGP2tTjpkbap+d8f3lrm6+lix9nv4AvCgiVo+Il1H6VPW0nDEtr+bxuznlFNWYYY7VLwLPoSRfXxibGBHPioh/qonHjZT69Xd6iIg16vGyGjCrHrur13nbUVr7XpuZ3+q1fo/y1qT0QQK4dx0fcybw8oiYHRGzKadWzxtfRi1n3yi3s1g7IlaLiGdQ+hidkZlLKcnBByNinTr/gfXKuTEbAa+rr28vSuf4bw/zGnoYa118Sy3vyZSE9tjlLO9uJngP5tTXFlFuBfEh4F3N1rd67DyFkoiPL/tR9bifS2kV/FZtCesVxzqUC2l+kpm9Wkm/AOwXEdvWBO0dlNOsY+vPqq9jdWD1+jpm1XlrRsR29XVsSekT+tHMvH6Su0tDMvFS072BwyktC3+kfEAeVOd9tT7/KSLOrsN7U/o3XAV8Azg4M0+t895Faar+HaUvyPGUX8I9ZeavKKfdfkb5gn4o5QqktgyKtW2vppx6uYnS8fu4xrz3UPoKfSpLR+wXAe+Ocnn53WTmXyidiH9ST9s8LjO/QWmFPLaezruAcmXk8vojpX/QVdTL2xsf/m+lnBb7ed3W9yktWv18hNJH5FpKB+zvjpv/UeB5UW5UO9bS9wrKFYN/onx5T9RyN9mYltehlNMuv6MkEl8cmzHMsZqZV1I6cCelg/iYrWvMN9f1P5mZi/rE8FlKIrs38PY6/OI6742UU6yfn8Rpxb/W7UK5mu+vjXkvo9SPKynJ8APocQFOdSPlc+FyypWE7wNelZk/rvNfQum39CvKsXU80LxI5gzKfriWcnw/b9iWmfEy83Zgd0oduJZyBeNL+iUwy2HQe7Ahd11d+R3gyMw8Ytz6L6ZcuPPbHmV/lLL/fl2fXzEgjudQLgp5aa/TkFkugHof5XTn7+uj+cPnHTX2BZTPnL/WaQBrUlrpb6b0t/sZpc+ZRiTGtX5KIxERrwJekJk7TriwOlFbB75UT5OoZRFxJHBVZr5jwoVXERGxL6Wz/BOmOxZpuqxqNyRUR6LcBuIBlF9PW1N+oX98WoOSOhLlzvDPpdxOQZL+wVONGpV7Ufot3ES5EeIJlNMA0owWEYdRTgG/PzN/N93xSFqxeKpRkiSpI7Z4SZIkdcTES5IkqSMrRef6DTfcMOfNmzfdYUiSJE3orLPOujYze/6rxkqReM2bN4/FixdPdxiSJEkTiojf95vnqUZJkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyvFfzVKK4N5C04eaflLDt9tpOVLkkbPFi9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6sjIE6+IWD0izomIk+r4nIg4NSIuqc/rjzoGSZKkFUEXLV4HAhc1xhcAp2Xm1sBpdVySJGnGG2niFRGbA7sBn2tM3gNYWIcXAnuOMgZJkqQVxahbvD4CvAW4szFt48xcClCfNxpxDJIkSSuEkSVeEfEs4JrMPGs5198/IhZHxOJly5a1HJ0kSVL3RtnitQOwe0QsAY4FdoqILwFXR8QmAPX5ml4rZ+YRmTk/M+fPnTt3hGFKkiR1Y2SJV2a+LTM3z8x5wAuA/8vMFwEnAvvUxfYBThhVDJIkSSuS6biP1+HALhFxCbBLHZckSZrxZnWxkcxcBCyqw38Cdu5iu5IkSSsS71wvSZLUkU5avCStHOYtOHmk5S85fLeRli9JKzpbvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpI7OmOwBJw5m34OTpDkGSNEW2eEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdGVniFRFrRsQvIuK8iLgwIg6t0+dExKkRcUl9Xn9UMUiSJK1IRtnidRuwU2Y+HNge2DUiHgcsAE7LzK2B0+q4JEnSjDeyxCuLm+voGvWRwB7Awjp9IbDnqGKQJElakYy0j1dErB4R5wLXAKdm5hnAxpm5FKA+bzTKGCRJklYUI028MvPvmbk9sDnwmIjYbth1I2L/iFgcEYuXLVs2uiAlSZI60slVjZl5A7AI2BW4OiI2AajP1/RZ54jMnJ+Z8+fOndtFmJIkSSM1yqsa50bEenV4NvBU4GLgRGCfutg+wAmjikGSJGlFMmuEZW8CLIyI1SkJ3nGZeVJE/Aw4LiL2Ay4H9hphDJIkSSuMkSVemflL4BE9pv8J2HlU25UkSVpReed6SZKkjph4SZIkdWTCxCsi3hcR60TEGhFxWkRcGxEv6iI4SZKkmWSYFq+nZeaNwLOAK4EHAW8eaVSSJEkz0DCJ1xr1+ZnAMZl53QjjkSRJmrGGuarxWxFxMfBX4NURMRe4dbRhSZIkzTwTtnhl5gLg8cD8zPwb8BfKH11LkiRpEiZs8YqIxcCRwDHA9Zl5C3DLqAOT2jZvwcnTHYIkaRU3TB+vFwCbAWdGxLER8fSIiBHHJUmSNOMMc6rx0sx8O+Vqxi9TWr8uj4hDI2LOqAOUJEmaKYa6gWpEPAz4IPB+4GvA84Abgf8bXWiSJEkzyzB9vM4CbgA+DyzIzNvqrDMiYodRBidJkjSTDHM7ib0y87JeMzLzuS3HI0mSNGNNmHhl5mURsRvwEGDNxvR3jTIwSTPPqK8sXXL4biMtX5Kmapj/avw08G/Aa4EA9gK2GnFckiRJM84wnev/JTNfQrmH16GUm6luMdqwJEmSZp5hEq+/1ue/RMSmwN+A+48uJEmSpJlpmM71J0XEepRbSZwNJPC5kUYlSZI0Aw3Tuf6wOvi1iDgJWDMz/zzasCRJkmaevolXRPS9VUREkJlfH01IkiRJM9OgFq9n1+eNgH/hrrvUPwVYBJh4SZIkTULfxCszXwpQTy9um5lL6/gmwCe6CU+SJGnmGOaqxnljSVd1NeUPsyVJkjQJw1zVuCgivgccQ7mi8QXA6SONSpIkaQYa5qrGAyLiOcCT6qQjMvMbow1LkiRp5hmmxYuaaJlsSZIkTcEwfbwkSZLUAhMvSZKkjvRNvCLitPr83u7CkSRJmrkG9fHaJCJ2BHaPiGOBaM7MzLNHGpkkSdIMMyjxeiewANgc+NC4eQnsNKqgJEmSZqJBd64/Hjg+Iv6z8UfZkiRJWk7D3MfrsIjYnbvu47UoM08abViSJEkzz4RXNUbEe4ADgV/Vx4F1miRJkiZhmBuo7gZsn5l3AkTEQuAc4G2jDEySJGmmGfY+Xus1htcdRSCSJEkz3TAtXu8BzomI0ym3lHgStnZJkiRN2jCd64+JiEXAoymJ11sz84+jDkySJmvegpNHWv6Sw3cbafmSZr5h/yR7KXDiiGORJEma0fyvRkmSpI6YeEmSJHVkYOIVEatFxAVdBSNJkjSTDUy86r27zouILTuKR5IkacYapnP9JsCFEfEL4JaxiZm5+8iikiRJmoGGSbwOHXkUkiRJq4Bh7uP1g4jYCtg6M78fEfcBVh99aJIkSTPLMH+S/QrgeOAzddJmwDdHGZQkSdJMNMztJF4D7ADcCJCZlwAbjTIoSZKkmWiYxOu2zLx9bCQiZgE5upAkSZJmpmESrx9ExEHA7IjYBfgq8K3RhiVJkjTzDJN4LQCWAecD/w58G3jHKIOSJEmaiYa5qvHOiFgInEE5xfjrzPRUoyRJ0iQNc1XjbsBvgY8BHwcujYhnDLHeFhFxekRcFBEXRsSBdfqciDg1Ii6pz+tP9UVIkiStDIY51fhB4CmZ+eTM3BF4CvDhIda7A3hjZj4YeBzwmojYlnLq8rTM3Bo4rY5LkiTNeMMkXtdk5qWN8cuAayZaKTOXZubZdfgm4CLKPcD2ABbWxRYCe04qYkmSpJVU3z5eEfHcOnhhRHwbOI7Sx2sv4MzJbCQi5gGPoPQT2zgzl0JJziKi5z3BImJ/YH+ALbf0P7olSdLKb1Dn+mc3hq8GdqzDy4Ch+2VFxFrA14DXZ+aNETHUepl5BHAEwPz58+3ML0mSVnp9E6/MfOlUC4+INShJ19GZ+fU6+eqI2KS2dm3CEKctJUmSZoIJbycREfcHXgvMay6fmbtPsF4AnwcuyswPNWadCOwDHF6fT5h01JIkSSuhCRMvyh9if55yt/o7J1H2DsCLgfMj4tw67SBKwnVcROwHXE7pMyZJkjTjDZN43ZqZH5tswZn5Y6Bfh66dJ1ueJEnSym6YxOujEXEwcApw29jEsVtFSJIkaTjDJF4PpZwy3Im7TjVmHZckSdKQhkm8ngM8IDNvH3UwkiRJM9kwd64/D1hv1IFIkiTNdMO0eG0MXBwRZ3L3Pl4DbychSZKkuxsm8Tp45FFIkiStAiZMvDLzB10EIs1bcPJ0hyBJ0kgNc+f6myhXMQLcC1gDuCUz1xllYJIkSTPNMC1eazfHI2JP4DEji0iSJGmGGuaqxrvJzG/iPbwkSZImbZhTjc9tjK4GzOeuU4+SJEka0jBXNT67MXwHsATYYyTRSJIkzWDD9PF6aReBSJIkzXR9E6+IeOeA9TIzDxtBPJIkSTPWoBavW3pMuy+wH7ABYOIlSZI0CX0Tr8z84NhwRKwNHAi8FDgW+GC/9SRJktTbwD5eETEHeAPwQmAh8MjMvL6LwCRJkmaaQX283g88FzgCeGhm3txZVJIkSTPQoBuovhHYFHgHcFVE3FgfN0XEjd2EJ0mSNHMM6uM16bvaS5IkqT+TK0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqyKzpDkCSVhbzFpw80vKXHL7bSMtf2eOXZgJbvCRJkjoyssQrIo6MiGsi4oLGtDkRcWpEXFKf1x/V9iVJklY0o2zxOgrYddy0BcBpmbk1cFodlyRJWiWMLPHKzB8C142bvAewsA4vBPYc1fYlSZJWNF338do4M5cC1OeNOt6+JEnStFlhO9dHxP4RsTgiFi9btmy6w5EkSZqyrhOvqyNiE4D6fE2/BTPziMycn5nz586d21mAkiRJo9J14nUisE8d3gc4oePtS5IkTZtR3k7iGOBnwDYRcWVE7AccDuwSEZcAu9RxSZKkVcLI7lyfmXv3mbXzqLYpSZK0IlthO9dLkiTNNJGZ0x3DhObPn5+LFy+e7jBWeaP+nzdJmoj/B6mVQUSclZnze82zxUuSJKkjJl6SJEkdMfGSJEnqiImXJElSR0y8JEmSOmLiJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqyKzpDkDt8b8UJUlasdniJUmS1BETL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkd8apGSdJKY9RXby85fLeRli/Z4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHfGqRkmSKq+a1KjZ4iVJktQREy9JkqSOmHhJkiR1xMRLkiSpIyZekiRJHTHxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOmHhJkiR1xP9q7Mio//9LkiSt+GzxkiRJ6oiJlyRJUkdMvCRJkjpi4iVJktQREy9JkqSOeFWjJEkdWdmvcF9y+G7THcJKzxYvSZKkjph4SZIkdcTES5IkqSMmXpIkSR0x8ZIkSeqIVzVKkqQVQhdXfU73lZm2eEmSJHVkWhKviNg1In4dEZdGxILpiEGSJKlrnSdeEbE68AngGcC2wN4RsW3XcUiSJHVtOlq8HgNcmpmXZebtwLHAHtMQhyRJUqemI/HaDLiiMX5lnSZJkjSjTcdVjdFjWt5joYj9gf3r6M0R8euRRgUbAteuxOV3sQ3Ln9nld7ENy7f8FX0blj9AvHfljh86eQ0AW/WbMR2J15XAFo3xzYGrxi+UmUcAR3QVVEQszsz5K2v5XWzD8md2+V1sw/Itf0XfhuXP7PK72sYg03Gq8Uxg64i4f0TcC3gBcOI0xCFJktSpzlu8MvOOiDgA+B6wOnBkZl7YdRySJEldm5Y712fmt4FvT8e2Bxj1ac0uTpuu7K/B8qe3/C62YfmWv6Jvw/JndvldbaOvyLxHv3ZJkiSNgH8ZJEmS1JFVMvGKiCMj4pqIuKAxbU5EnBoRl9Tn9Vsuf6+IuDAi7oyIKV1N0af890fExRHxy4j4RkSsN4JtHFbLPzciTomITdssvzHvTRGREbFhy/EfEhF/qPGfGxHPbDv+iHht/TusCyPifS3H/5VG7Esi4tyWy98+In5ey18cEY9pufyHR8TPIuL8iPhWRKwzhfK3iIjTI+Kiuq8PrNNbqccDym+zHvfbRit1eUD5rdTjfuU35k+pHg+Iv5V6PCj+NurxgPhbqccDym+zHvfbRit1OSLWjIhfRMR5tfxD6/S26nG/8lurx8slM/m0EqAAAAqVSURBVFe5B/Ak4JHABY1p7wMW1OEFwHtbLv/BwDbAImD+COJ/GjCrDr93KvEP2MY6jeHXAZ9us/w6fQvKhRe/BzZsOf5DgDeN8Bh6CvB94N51fKO2909j/geBd7Yc/ynAM+rwM4FFLZd/JrBjHX4ZcNgUyt8EeGQdXhv4DeUvyFqpxwPKb7Me99tGK3V5QPmt1ON+5dfxKdfjAfG3Uo8HlN9KPR60fxrLLHc9HhB/m/W43zZaqcuU+3quVYfXAM4AHtdiPe5Xfmv1eHkeq2SLV2b+ELhu3OQ9gIV1eCGwZ5vlZ+ZFmdnKTWD7lH9KZt5RR39OuT9a29u4sTF6X3rc+HYq5VcfBt4ylbInKL8Vfcp/FXB4Zt5Wl7mm5fIBiIgAng8c03L5CYz9cl2XHvfXm2L52wA/rMOnAv86hfKXZubZdfgm4CLKP2C0Uo/7ld9yPe63jVbq8oDyW6nHA94DaKEeT1D+lA0ov5V6PFH8U63HA8pvsx7320YrdTmLm+voGvWRtFePe5bfZj1eHqtk4tXHxpm5FMrBBmw0zfFMxcuA74yi4Ij4r4i4Angh8M6Wy94d+ENmntdmueMcUE+zHLm8zdcDPAh4YkScERE/iIhHt1z+mCcCV2fmJS2X+3rg/fX9/QDwtpbLvwDYvQ7vxd1vpLzcImIe8AjKr9nW6/G48kdiwDZaqcvjy2+7HjfLH0U97rF/Wq3H48pvvR73eX9bq8fjyh9JPR63jdbqckSsXk+3XgOcmpmt1uM+5U8rE68ZJiLeDtwBHD2K8jPz7Zm5RS3/gLbKjYj7AG+n5WRunE8BDwS2B5ZSmvnbNAtYn9KU/WbguPqrtm17M4XWrgFeBfxHfX//A/h8y+W/DHhNRJxFOW1x+1QLjIi1gK8Brx/XktOKUZc/aBtt1eVe5bdZj5vl13hbrcc94m+1Hvcov9V6POAYaqUe9yi/9XrcYxut1eXM/Htmbk9p2X1MRGw31Xi7LH95mHjd5eqI2ASgPi/3aaLpEhH7AM8CXpj1pPYIfZkpnCrq4YHA/YHzImIJpZKcHRH3a2sDmXl1rYR3Ap8FlrvTaR9XAl+vzdu/AO6k/O9YayJiFvBc4CttllvtA3y9Dn+VlvdPZl6cmU/LzEdRvnB+O5XyImINypfB0Zk5Fndr9bhP+a3qt4226vIQr2FK9bhH+a3W417xt1mP++yf1urxgPe3lXrcp/xW63Gf96DVulzLvIHS52pXRvB9PK78aWXidZcTKQcs9fmEaYxl0iJiV+CtwO6Z+ZcRbWPrxujuwMVtlZ2Z52fmRpk5LzPnUT78HpmZf2xrG2MVuXoOpbm8Td8EdqrbehBwL9r/I9anAhdn5pUtlwulL8iOdXgnoNVTmRGxUX1eDXgH8OkplBWUX/IXZeaHGrNaqccDym9Nv220VZcHlN9KPe5Vfpv1eED8rdTjAe9xK/V4gmNoyvV4QPmt1eMB70ErdTki5ka9ajciZlP3C+3V437lT6/suDf/ivCgZOhLgb9RPhj2AzYATqMcpKcBc1ou/zl1+DbgauB7LZd/KXAFcG59LPcVhwO28TXKh9wvgW9ROuq2Vv64+UuY2lWNveL/InB+jf9EYJOWy78X8KW6j84Gdmp7/wBHAa8cUR14AnAWcB6lH8ejWi7/QMpVUb8BDqfewHk5y38CpRPuLxvH/DPbqscDym+zHvfbRit1eUD5rdTjfuWPW2a56/GA+FupxwPKb6UeD9o/bdTjAfG3WY/7baOVugw8DDinln8B9QrPFutxv/Jbq8fL8/DO9ZIkSR3xVKMkSVJHTLwkSZI6YuIlSZLUERMvSZKkjph4SZIkdcTES1pFRURGxBcb47MiYllEnDSCbS2KiPkTLPO5iNi2Di+JiKFuWhkRW0TE6RFxUURcGBEHTjHWg6ayflciYt+I2HS645A0OSZe0qrrFmC7emNBgF2AP0xXMJn58sz81XKsegfwxsx8MOVvXl4zlsAtpxUm8YqI1QfM3heYVOJV75guaRqZeEmrtu8Au9Xhu/13XEQ8JiJ+GhHn1Odt6vQ3RMSRdfihEXFB/a9NGuvOjohjo/yR8VeA2Y15n4qIxbV16tDG9Hu0ikXEYc0WrCh/7vy65jKZuTQzz67DNwEXAZs1yvxwRPywtog9OiK+HhGXRMS7x++MiDgcmB0R50bE0XXaiyLiF3XaZ8aSoYi4OSLeGxFnRcT36/5aFBGXRfmj6LFWqRMi4rsR8euIOLixrUHlvisizgAeHxHvjIgz634+IornAfOBo+v6s5uthBExPyIW1eFD6nqnAF+od/P+Wi3zzIjYodeBIWlEurxbqw8fPlacB3Az5c7OxwNrUu5K/WTgpDp/HWBWHX4q8LU6vBrwQ8rdnxcDO/Qo+w3AkXX4YZRWqfl1fE59Xp3y32kPq+OLGsssofw/3jzg7MZ2fwtsMOA1zQMuB9ZplPneOnwg5e9UNgHuTblz9T3KAm5uDD+Ycnf3Ner4J4GX1OEEnlGHvwGcAqwBPBw4t07fl3IH/w0oyecFlIRponKf34hhTmP4i8Czx++v5j6rw/OBRXX4EMqdzGfX8S8DT6jDW1L+Dmbaj0cfPlaVh83O0iosM38ZEfMorV3fHjd7XWBhlP/2S0pSQWbeGRH7Uv6G4zOZ+ZMeRT8J+FhjG79szHt+ROwPzKIkQdvWsnrFtyQi/hQRjwA2Bs7JzD/1WjYi1qL8Hc7rM/PGxqwT6/P5wIWZubQufxmwBdCzvGpn4FHAmeVv65jNXX/Yezvw3UbZt2Xm3yLifEoCOObUsZgj4uuUv2G5Y0C5f6+vY8xTIuItwH2AOcCFlKRtMk7MzL/W4acC29btAqwTEWtnaS2UNGImXpJOBD5Aae3aoDH9MOD0zHxOTc4WNeZtTWkxG9TH6B7/RxYR9wfeBDw6M6+PiKMorW2DfI7ScnQ/4MheC0TEGpRk5ejM/Pq42bfV5zsbw2PjE30GBrAwM9/WY97fMnPsNf6j7JqYNssdvx9ygnJvzcy/A0TEmpTWsPmZeUVEHEL//XUHd3UfGb/MLY3h1YDHNxIxSR2yj5ekI4F3Zeb546avy12d7fcdmxgR6wIfpbRqbVD7G433Q+CFdfntKKcboZy+vAX4c0RsDDxjiPi+AewKPBr43viZUZpuPk85ZfahIcqbyN9qIgflD3qfFxEb1W3NiYitJlneLnW92cCewE8mUe5YAnVtbdFr7uubgLUb40sorWgA/zognlOAA8ZGImL7SbwWSVNk4iWt4jLzysz8aI9Z7wPeExE/ofTHGvNh4JOZ+RtgP+DwsQSi4VPAWvUU41uAX9RtnQecQzlddiQlCZkovtuB04HjxlqCxtkBeDGwU+1ofm5EPHOicgc4AvhlRByd5SrLdwCn1NdyKuX06GT8mNI361xKP7nFw5abmTcAn6WcyvwmcGZj9lHAp8c61wOHAh+NiB9RTlf28zpgfr3w4VfAKyf5eiRNQdzVUi5JK56IWA04G9grMy+Z7ngmo/aFm5+ZB0y0rKRVgy1eklZYUe7HdSlw2sqWdElSL7Z4SZIkdcQWL0mSpI6YeEmSJHXExEuSJKkjJl6SJEkdMfGSJEnqiImXJElSR/4/sHDfmqCIuGkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "plt.hist(Lille_7days, bins = np.arange(10,32,1)) \n",
    "plt.title(\"histogram of max temperature in days 12-18 Sep from 1979-2019\")\n",
    "plt.xlabel(\"Max daily 2m temperature\")\n",
    "plt.ylabel(\"Number of days\")\n",
    "plt.xticks(np.arange(10,32,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even in this increased temporal range, the maximum daily temperature still never reached that of 15 September 2020!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}