{ "cells": [ { "cell_type": "markdown", "id": "360075b1-2ac9-4ca1-9615-403e5481d250", "metadata": {}, "source": [ "# Calculating zonal means" ] }, { "cell_type": "code", "execution_count": 1, "id": "ec783093-5325-45d2-8712-a741e8faa220", "metadata": {}, "outputs": [], "source": [ "import intake\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from gridlocator import merge_grid" ] }, { "cell_type": "markdown", "id": "de341c08-79f1-4a4c-a93e-997fbf8e87cc", "metadata": {}, "source": [ "In this example we will compute the zonal-mean air temperature.\n", "We retrieve the 2m-air-temperature from the `dpp0067` NextGEMS simulation using `intake-esm`." ] }, { "cell_type": "code", "execution_count": 2, "id": "06aee1d5-5f60-4586-9b10-87a5315164b3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "--> The keys in the returned dictionary of datasets are constructed as follows:\n", "\t'project.institution_id.source_id.experiment_id.simulation_id.realm.frequency.time_reduction.grid_label.level_type'\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [1/1 00:00<00:00]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "catalog_file = \"/work/ka1081/Catalogs/dyamond-nextgems.json\"\n", "col = intake.open_esm_datastore(catalog_file)\n", "cat = col.search(\n", " variable_id=\"tas\",\n", " project=\"NextGEMS\",\n", " simulation_id=\"dpp0067\",\n", ")\n", "cat_dict = cat.to_dataset_dict(cdf_kwargs={\"chunks\": {\"time\": 1}})\n", "\n", "ds = merge_grid(list(cat_dict.values())[0]) # Include the grid information!" ] }, { "cell_type": "markdown", "id": "d844c228-de01-4188-8a39-e082385a8843", "metadata": {}, "source": [ "## The general idea\n", "behind the zonal average is to calculate a weighted average of values in a certain latitude bin.\n", "Therefore, in a first step, we count how many cells are in given equidistant latitude bins." ] }, { "cell_type": "code", "execution_count": 3, "id": "a8430e1a-0f1b-4148-b15f-4dac307fbe86", "metadata": {}, "outputs": [], "source": [ "hist_opts = dict(bins=128, range=(-np.pi / 2, np.pi / 2))\n", "cells_per_bin, lat_bins = np.histogram(ds.clat, **hist_opts)" ] }, { "cell_type": "markdown", "id": "a10c6bda-54c4-491c-9925-5b308ddfb3ff", "metadata": {}, "source": [ "Now comes the trick! In a next step, we will repeat the histogram but account a weight to each cell.\n", "Usually, `histogram` will weight each data point with one, i.e. it will count the values in a certain bin.\n", "Here, we will weight each data point with the 2m-temperature. Thereby, we will compute the cumulative sum of all temperatures in a given latitude bin.\n", "\n", ".. tip::\n", " The `np.histogram` function is more efficient when passing a range and a number of bins.\n", " This is because, when constructing the bins internally, the function can assume equidistant bin sizes.\n", " This is not the case when passing a sequence of bins directly." ] }, { "cell_type": "code", "execution_count": 4, "id": "5bbb275b-3389-4008-b37b-22fb37b942c1", "metadata": {}, "outputs": [], "source": [ "varsum_per_bin, _ = np.histogram(\n", " ds.clat, weights=ds.tas.isel(time=1, height_2=0), **hist_opts\n", ")" ] }, { "cell_type": "markdown", "id": "5d43cca2-88b3-4d44-8d7c-08d65ffd0a07", "metadata": {}, "source": [ "The zonal mean can now be computed by dividing the cumulative values of the temperature with the number of cells in each bin." ] }, { "cell_type": "code", "execution_count": 5, "id": "70c53f39-1422-48eb-9a78-aad2b8f83d9f", "metadata": {}, "outputs": [], "source": [ "zonal_mean = varsum_per_bin / cells_per_bin" ] }, { "cell_type": "markdown", "id": "9ad77d8e-2bea-4bb2-a813-c21e2efadae3", "metadata": {}, "source": [ "We can check our result by plotting the zonal mean as a function of the latitdue bins.\n", "While doing so, we will scale the bins by their area so that the visual appearance of each latitude bin represents their actual proportion." ] }, { "cell_type": "code", "execution_count": 6, "id": "718e97e5-2ada-480b-b4d7-3df5e41e705e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 0, 'latitude')" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvtElEQVR4nO3deXxV9Z3/8dcnO9kIhARCFhL2fY2gUFCsK+4LiFrr0inVsa3btGrtdGw7dmw7dX5Oq2Np1Vo3REHFui8oymoSdsISIISQQEIgGyH75/fHPaQRAgTIzcm99/N8PPLg5pxzT943Xu8n3+/5nu9XVBVjjDEGIMjtAMYYY7oOKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoXXioKIRIjIKhFZKyIbReSXzvaeIvKxiGxz/u3hbE8XkcMissb5esZb2YwxxrRNvHWfgogIEKWq1SISCnwF3ANcCxxQ1cdF5CGgh6o+KCLpwD9UdaRXAhljjDkpr7UU1KPa+TbU+VLgKuAFZ/sLwNXeymCMMebUhHjz5CISDGQDA4GnVHWliPRW1WIAVS0WkcRWT8kQkdVAJfBzVf2yjXPOAeYAREVFTRg6dKg3X4Ixxvid7Ozs/aqa0NY+r3UffeOHiMQBbwI/Ar5S1bhW+w6qag8RCQeiVbVMRCYAbwEjVLXyeOfNzMzUrKwsr2Y3xhh/IyLZqprZ1r5OGX2kquXA58AlwD4RSXKCJQElzjF1qlrmPM4GtgODOyOfMcYYD2+OPkpwWgiISDfgAmAzsAi41TnsVuDtVscHO4/7A4OAHd7KZ4wx5ljevKaQBLzgfNAHAfNV9R8ishyYLyLfAwqAmc7x04BfiUgj0ATcqaoHvJjPGGPMUbxWFFR1HTCuje1lwLfb2L4AWOCtPMYYY07O7mg2xhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGlhRcEYY0wLKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGlhRcEYY0wLKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGkR4q0Ti0gEsAQId37OG6r6HyLSE3gNSAfygVmqetB5zsPA94Am4Meq+qG38pnAVV3XSEFZDQUHaqg83EBdYxN1jc2er4Z/Pq5taCIiNJgRfWMZmdydgYnRhAbb31HGv3mtKAB1wPmqWi0iocBXIvI+cC3wqao+LiIPAQ8BD4rIcGA2MALoC3wiIoNVtcmLGY2fq65rJHvXQVbuKCMr/yDbS6spO1R/3ONFIDwkiPCQYCJCg6iqbaSm3vMWDAsJYkBCNL2iw+gZFUaPyDASYsIZ2ieG0SlxJMSEH/e8VbUN5JVUs62kmrLqehqbmmloVhqbmokKD2F431hG9u1+wnMY0xm8VhRUVYFq59tQ50uBq4DznO0vAJ8DDzrb56lqHbBTRPKAicByb2U0/qXicANb9laRW1zJ5r2VbCzyfDU1K8FBwqjk7lw4vDdp8ZH06xlFWs9IekSFthSA8JBgQoMFEWk5Z1Ozkl92iA17KthYVEleSTUHDtVTcKCGA4fqqaptbDk2qXsEI/rGEhYSRP2RlkdjM4UHaiiqqG0zc0iQ0NisLd/3jg1nVHIcZ/fvyZSBvRjSO4agIGnzucZ4gzdbCohIMJANDASeUtWVItJbVYsBVLVYRBKdw5OBFa2eXuhsO/qcc4A5AGlpad6Mb3xAxeEGPtywl7fX7mHZ9jLU+XyNiwxlWJ9Y7jp3AJP692R8Wg+iwk/97R4cJAxIiGZAQjRXjT3m7cihukY2FlWyrrCc9XsqyC32FKGwkGDCQ4IICwliYkZPBvWOYVBiNIN7x5AYG05ocBAhQZ4CVFnbwKaiypbCs7rgIJ/k7gMgPiqMyQN7MfusVCYPiP9GwTLGG7xaFJyun7EiEge8KSIjT3B4W+92PWaD6lxgLkBmZuYx+43/O3ionsVbSvhw414WbymlvrGZ9PhIfjh9IOP79WBYn1h6x4Z3ygdoVHgIEzN6MjGj52mfIzYilLP7x3N2//iWbXvKD7Msbz/Lt5fx+dZS3llbxODe0dw2OYNrxiXTLSy4I+IbcwyvFoUjVLVcRD4HLgH2iUiS00pIAkqcwwqB1FZPSwGKOiOf6ZqampWy6jr2Vtayt6KW7aWH+GzzPrJ3HaRZITEmnJsnpXH12GRGp3T3q7+ik+O6MTMzlZmZqdQ2NPHO2iKeX5rPz95cz28/2Mz3p2Zw+5SM02r9GHMiouqdP7ZFJAFocApCN+Aj4LfAuUBZqwvNPVX1pyIyAngFz3WEvsCnwKATXWjOzMzUrKwsr+Q3nauhqZn1eyrYsKeC9YUVrN9TQV5J9Tf62wGGJ8VywbBEvj2sN6OSuwdUf7uq8nX+QeYu2c4nuSX0ig7j7ukDuWlSGuEh1nIw7Sci2aqa2eY+LxaF0XguJAfjuR9ivqr+SkTigflAGlAAzFTVA85zHgHuABqBe1X1/RP9DCsKvq3icAOfbynh09wSPt9SQqVz0TY+KoyRyd0ZlhRLclwEvWMj6NM9guS4bsRH2+gcgJyCg/z+gy0s31FGclw3vj81g+smpBATEep2NOMDXCkKncGKgu9QVbbuq2bt7nLWFJazdnc5m/dW0dSsxEeFMX1oIucPTWRsahxJ3SP8qivIW1SVpXllPPHxFnIKyokKC+b6CSncck46AxOj3Y5nujArCsYVBw/V82Xefr7YUsqSbaWUVtUBEBMRwtjUOMamxnHeEE8hCA6gbiBvWLu7nBeW5fOPdcXUNzXz7aGJ3HPBIEanxLkdzXRBVhRMp9pVdoinF29nQU4hjc1KXGQoUwclMHVQLzL79SA9PiqgrgV0pv3Vdby8ooDnlu6k4nCDFQfTJisKplNsL63mqc/yeHttEcFBwuyzUrlmXDKjU6wl0Nmqahv4+/Jd/OXLHZTXNHDZqCQevXKE3TFtACsKxsuampU/L9nO/3y8lZCgIG6elMacaf1JjI1wO1rAq6pt4NmvdvL04u1EhQfz6JUjuHJMX7tmE+CsKBiv2X2ghgfmr2VV/gEuG5XEL68aQS8bIdTl5JVU8W+vr2PN7nIuHN6bx64eaUU7gJ2oKNiUj+a0LVpbxKVPfsmm4kqemDWGP900zgpCFzUwMYYFd03mkRnDWLK1lOufWc7BE0wMaAKXFQVzWl77uoAfv7qaoX1ieP+eqVw7PsW6JLq44CDh+9P688r3z2ZvRS0/fDWHxqZmt2OZLsaKgjll87/ezUML1zNtcAIv/cskUntGuh3JnIIJ/Xrw2DUjWZpXxn++m+t2HNPF2MQp5pS8nrWbBxeuY+qgBObeMoGIUJtewRfNzEwlt7iK55buZHhSLLPOSj35k0xAsJaCabfXvi7gpwvW8a2Bvawg+IGfzRjK1EG9+PlbG8jedcDtOKaLsKJgTkpVeWpxHg8uWM+0QQn85buZVhD8QEhwEH+8cRxJcRF874Us1hdWuB3JdAFWFMwJNTcrv3xnE7//cAvXjEu2guBn4iLDePGOSUSHh3DTX1eQveug25GMy6womOOqa2ziR/NW87dl+Xx/agZ/mDmGsBB7y/ibtPhI5v/gHOKjwvjusytZuaPM7UjGRfZ/uGlTXWMTd72Uw7vrivnZjKE8ctlwm6/Ij/WN68b8H5xDUlw3bn1+FV9t2+92JOMSKwrmGLUNTfzgxWw+21zCY9eMZM60AW5HMp0gMTaCeXPOJj0+ijkvZrGpqNLtSMYFVhTMN9Q2NDHnxWw+31LKf107ipsn9XM7kulEvaLD+fsdE+neLZR/eeHrlunOTeCwomBa1Dc2M+fFbL7cVsrvrhvNjRPT3I5kXJAYG8FfvpvJgZp6fvBiFrUNx10R1/ghKwoG8Aw7/cXbG1iytZTHrx1lNzMFuJHJ3fmfWWPJKSjnZwvX48sTZ5pTY0XBAPDsVzuZ9/Vu7p4+gBvOshaCgUtHJXH/hYNZuHoPf16yw+04ppNYUTB8mruPx97L5dKRfXjgwiFuxzFdyI/OH8hlo5L4/YdbWF1g9zAEAisKAW7z3kp+/OpqRvSN5Q+zxtiwU/MNIsJ/XTeKPrER3PfaGg7VNbodyXiZFYUAVtvQxN0v5xAVHsJfv3sWkWE2P6I5VmxEKE/MGsOuAzX8+h+b3I5jvMxrRUFEUkVksYjkishGEbnH2T5GRJaLyHoReUdEYp3t6SJyWETWOF/PeCub8fjtB5vZXnqIJ2aNpU93W4XLHN+k/vHcee4A5n29mw837nU7jvEib7YUGoEHVHUYcDZwt4gMB/4KPKSqo4A3gZ+0es52VR3rfN3pxWwBb2nefp5fms9tk9P51qBebscxPuC+CwYzom8sDy1YR0llrdtxjJd4rSioarGq5jiPq4BcIBkYAixxDvsYuM5bGUzbKg438G+vr6V/QhQPXjLU7TjGR4SFBPHk7LHU1DfxsA1T9Vudck1BRNKBccBKYANwpbNrJtB6QHyGiKwWkS9EZOpxzjVHRLJEJKu0tNSbsf3WLxdtpKSqjv+ZNZZuYTbjqWm/gYkx/OTiIXy6uYSPNu1zO47xAq8XBRGJBhYA96pqJXAHnq6kbCAGOLJ6eDGQpqrjgPuBV45cb2hNVeeqaqaqZiYkJHg7vt95e80eFq7ew93TBzImNc7tOMYH3TY5naF9Yvjloo3U1NtoJH/j1aIgIqF4CsLLqroQQFU3q+pFqjoBeBXY7myvU9Uy53G2s32wN/MFmrySKh5euJ6z0nvwo/MHuh3H+KiQ4CB+ffVIiipq+eNneW7HMR3Mm6OPBHgWyFXVJ1ptT3T+DQJ+DjzjfJ8gIsHO4/7AIMBuo+wgh+oaufOlHCLDgvnTTeMJDbbRyOb0nZXek+snpPCXJTvIK6lyO47pQN78ZJgC3AKc32qY6QzgRhHZCmwGioDnneOnAetEZC3wBnCnqtrCsR1AVXnkzfVsL63mydnj6B1rw0/NmXv40qFEhYfw729ttIvOfsRrdyup6lfA8W6PfbKN4xfg6WoyHeyVVQW8taaI+y8czJSBNvzUdIz46HB+cvEQfv7WBhatLeKqscluRzIdwPoQ/Fz2roM8umgj5w5O4IfT7TqC6Vg3TkxjTEp3/vPdXCprG9yOYzqAFQU/tq+ylrteyiapezeenD3W5jUyHS44SPj11SPZX13HEx9tdTuO6QBWFPyUZ43lbKrrGpn73QnERYa5Hcn4qdEpcXxnUj/+vjyfjUUVbscxZ8iKgp96dNEmcgrK+f31Yxja55jbPYzpUP920RB6RIbx87c20NxsF519mRUFP/T2mj28uqqAu84bwGWjk9yOYwJA98hQfjZjGKsLynk9e7fbccwZsKLgZ4orDvPvb21gfFocD1xo9/6ZznPt+GQmpvfk8fc3c/BQ/cmfYLokKwp+RFX56RvraGhSnpg1lhC7Qc10IhHhV1ePoLK2kd99uNntOOY02aeGH3lpxS6+3LafRy4bRnqvKLfjmAA0tE8sd0xJ59VVu8mx5Tt9khUFP7GjtJrH3svl3MEJ3Dwpze04JoDdc8FgeseG8/M3N9DY1Ox2HHOKrCj4geZm5SdvrCM8JJjfXT8az7RTxrgjOjyEX1w+gk3Flby0YpfbccwpsqLgB15csYvsXQf5xeXDbV4j0yXMGNWHqYN68YePttoqbT7GioKPKzxYw+8+2My0wQlcO97mnjFdg4jwyytHUNfYzG/ey3U7jjkFVhR8mGf20w0o8JtrRlq3kelS+idE84Nz+/PWmiKW5u13O45pJysKPuytNXv4YmspP714CCk9It2OY8wx7p4+kP69orjvtTXsr65zO45pBysKPqqipoFfvbOJ8Wlx3HJOuttxjGlTRKhnUafyww3c99oamwLDB1hR8FF/WryN8sMN/PrqkQTb7KemCxveN5ZHrxjBl9v28/TntnxnV2dFwQftKjvE35blM3NCCiP6dnc7jjEndePEVK4c05cnPt7Kyh1lbscxJ2BFwQc9/v5mQoODeOCiIW5HMaZdRITfXDuKfvFR/OjV1ZTZ9YUuy4qCj1m18wDvb9jLnecOsHsSjE+JDg/hqZvGU17TwIML1tu6zl2UFQUf0tys/Oe7m+gTG8H3p/Z3O44xp2x431h+eskQPsndx7yvbYrtrsiKgg/5OHcf6wor+MnFQ+gWFux2HGNOyx1TMvjWwF786p1N7Nx/yO045ihWFHzI80t3khzXjavG9nU7ijGnLShI+O+ZYwgLCeLeeatpsEnzuhSvFQURSRWRxSKSKyIbReQeZ/sYEVkuIutF5B0RiW31nIdFJE9EtojIxd7K5otyiytZseMAt5zTz9ZJMD6vT/cI/uvaUawtrOCPn25zO45pxZufLo3AA6o6DDgbuFtEhgN/BR5S1VHAm8BPAJx9s4ERwCXA0yJifSSOvy3NJyI0iNlnpbodxZgOMWNUEteNT+Gpz7dTeLDG7TjG4bWioKrFqprjPK4CcoFkYAiwxDnsY+A65/FVwDxVrVPVnUAeMNFb+XzJgUP1vLVmD9eMSyEuMsztOMZ0mPsv8iwZ+6JNsd1ldEo/hIikA+OAlcAG4Epn10zgyJ++yUDr4QiFzraAN+/rAuoam7ltcrrbUYzpUMlx3bhoeG/mrdrN4fomt+MYOqEoiEg0sAC4V1UrgTvwdCVlAzHAkRW+25qr4ZiBzCIyR0SyRCSrtLTUW7G7jIamZl5cvovJA+IZ0ifG7TjGdLjbJqdTcbiBt9fscTuK4QRFQUSOO3+CiJzVnpOLSCiegvCyqi4EUNXNqnqRqk4AXgW2O4cX8s9WA0AKUHT0OVV1rqpmqmpmQkJCe2L4tI827qO4opbbp2S4HcUYr5iY0ZOhfWL427J8u6GtCzhRS+FTEelx9EYRuQhYeLITi2dy/2eBXFV9otX2ROffIODnwDPOrkXAbBEJF5EMYBCwqr0vxF/9bdlOUnt24/yhiW5HMcYrRITbp6SzeW8VK3cecDtOwDtRUfgzsFhEWv4cF5GbnO2XtePcU4BbgPNFZI3zNQO4UUS2ApvxtASeB1DVjcB8YBPwAXC3qgZ0J+OGPRV8nX+QW89Jt5lQjV+7amwycZGh/G1pvttRAl7I8Xao6l9EpBb4zGkd3ADcCUxX1fyTnVhVv6Lt6wQATx7nOY8Bj53s3IHi+aX5RIYFMzPThqEa/xYRGszss9KYu2Q7e8oPkxzXze1IAeuEF5pV9UXgV8Bq4CZgSnsKgjlz+6vreGdtEdeNT6F7t1C34xjjdd85Ow2AF5fb8FQ3nehC83oRWQf8AogE4vF0Jx3Zbrzo1ZUF1Dc1c+vkfm5HMaZTpPSI5PyhvVmYU2grtLnouN1HwOWdlsJ8Q31jMy+u2MXUQb0YmGjDUE3guGJMEp/k7iNr10EmZvR0O05AOtE1BWvDueT9DcWUVNXx2+tGux3FmE717WG9CQsJ4r31xVYUXGIzq3VBL68ooF98JOcO9v/7MIxpLTo8hPMGJ/D+hmLrQnKJFYUuZkdpNavyD3DDWakE2TBUE4AuG53Evso6sgsOuh0lIJ3oQvNcEblGRKxTuxPNzyokOEi4fnyK21GMccWRLqR31xW7HSUgnail8BwwBnhPRD4VkQdFZEwn5QpIDU3NLMgpZPqQBBJt/WUToKLDQzjXupBcc9yioKorVPVRVZ0KzAIKgAdEZLWIPCciszotZYBYvLmE0qo6bjgrze0oxrjqslGeLqQc60LqdO26pqCqZar6qqp+V1XHAU/hmZvIdKD5WbtJiAln+hC7wGwC27eHJXq6kNZbF1JnO60Lzaqa7UxJYTpISWUti7eUct34FFtu0wS8mIhQpg1K4P31e60LqZPZp08X8UZOIU3NyqxMu8BsDMDlo5PYW1nL6t3WhdSZrCh0AarK61mFTEzvSf+EaLfjGNMlfHtYIhGhQSzIscV3OtNJi4KIzDwyLFVEfi4iC0VkvPejBY5VOw+wc/8hZp1ls6Eac0RMRCgzRiWxaE0Rh+oa3Y4TMNrTUvh3Va0SkW8BFwMvAP/n3ViB5bWs3cSEhzBjVB+3oxjTpdw0MY3quka7Z6ETtacoHFno5jLg/1T1bSDMe5ECS2VtA++tL+aKsX2JDDvR/ITGBJ4J/XowKDGaV1YVuB0lYLSnKOwRkT/juVfhPREJb+fzTDu8s7aI2oZmbrCFdIw5hogwe2Iaa3aXk1tc6XacgNCeD/dZwIfAJapaDvQEfuLNUIFk/te7GdonhtEp3d2OYkyXdO24ZMJCgphnrYVOcdKioKo1qroQqBCRNCAUz/rK5gxt3lvJ2sIKZmWmImKT3xnTlh5RYcwY2YeFq/dwuD6gl23vFO0ZfXSliGwDdgJfOP++7+1ggeC1r3cTFhzENeOS3Y5iTJd248Q0qmob7Q7nTtCe7qNfA2cDW1U1A7gAWOrVVAGgoamZRWuKuHB4b3pE2XV7Y05kYkZP+idEWRdSJ2hPUWhQ1TIgSESCVHUxMNa7sfzfl9tKKTtUb60EY9pBRLhpYhpZuw6SlX/A7Th+rT1FoVxEooElwMsi8iRw0jtJRCRVRBaLSK6IbBSRe5ztY0VkhYisEZEsEZnobE8XkcPO9jUi8syZvLCu7s3VRfSIDGWara5mTLvcPKkfvWPD+c17uajafEje0p6icBVQA9wHfABsBy5vx/MagQdUdRie7qe7RWQ48Dvgl6o6FviF8/0R21V1rPN1Z/tfhm+pqm3go417uXx0X8JCbHSvMe3RLSyY+y8cTE5BOR9s2Ot2HL/Vnk+kX6hqs6o2quoLqvq/wIMne5KqFqtqjvO4CsgFkgEFYp3DugNFpxfdd324cR91jc1cM966jow5FdeNT2Fw72h++8FmGpqa3Y7jl9pTFC5sY9ulp/JDRCQdGAesBO4Ffi8iu4H/Bh5udWiGs4jPFyIy9VR+hi95c3Uh/eIjGZca53YUY3xKSHAQD106lPyyGl61i85ecaI1mu8SkfXAEBFZ1+prJ7CuvT/AuR6xALhXVSuBu4D7VDUVT5fUs86hxUCas4jP/cArIhLbxvnmONciskpLS9sbo8vYW1HLsu1lXD022e5NMOY0TB+SyNn9e/LkJ9uoqm1wO47fOVFL4RXgCmCR8++Rrwmq+p32nFxEQvEUhJedG+AAbgWOPH4dmAigqnXOKCdUNRvPtYvBR59TVeeqaqaqZiYk+N5F2kVr96AKV9uoI2NOi4jwsxnDKDtUz5+/2OF2HL9zojWaK1Q1X1VvVNVdrb7aNR5MPH8GPwvkquoTrXYVAec6j88HtjnHJ4hIsPO4P57lPv3qv3hzszJv1W7GpcWR0SvK7TjG+KzRKXFcPjqJ55buZH91ndtx/Io3h75MAW4Bzm81zHQG8H3gDyKyFvgNMMc5fhqwztn+BnBnewuQr/hiayk79h/itsnpbkcxxufdd+FgahuaeHrxdrej+BWvzdWsql8Bx+s0n9DG8QvwdDX5reeW7qR3bDiXjkxyO4oxPm9AQjTXT0jhpRW7+JepGfSN6+Z2JL9gg+Q7ybZ9VXy5bT+3nN3P7k0wpoP8+NuDUJQ/frbN7Sh+wz6dOsnzy/IJDwnixolpbkcxxm+k9Ijk5kn9mJ9VyM79h9yO4xesKHSC8pp6FuYUcs24ZOKjw92OY4xf+dfpAwgNFv7fJ1vdjuIXrCh0gldWFVDb0MxtU9LdjmKM30mMieD2KRksWltEXkmV23F8nhUFL6ttaOL5pflMHdSLoX2OuRfPGNMB/uVbGYQGB/Hc0ny3o/g8KwpetjBnD6VVddx57gC3oxjjt+Kjw7l6bF8W5hRSXlPvdhyfZkXBi5qalblLtjM6pTuTB8S7HccYv3b7lAxqG5p5ddVut6P4NCsKXvTBhr3kl9Vw57kDbJ4jY7xsWFIskwfE8/fl+TaD6hmwouAlqsozX2wno1cUF4/o43YcYwLC7VMyKK6o5cONtt7C6bKi4CXLtpexfk8Fc6b1JzjIWgnGdIbzhybSLz6S5+2C82mzouAlzy/dSXxUmK3BbEwnCg4Sbj0nnexdB1m7u9ztOD7JioIX7Co7xKebS7h5UhoRocFuxzEmoMzMTCE6PITnl+50O4pPsqLgBS8s20WwCN85u5/bUYwJODERoVw3Ppn31u/lwCEbnnqqrCh0sOq6Rl7P2s1lo5NIjI1wO44xAenms/tR39TM61k2PPVUWVHoYAuyC6mqa+T2KRluRzEmYA3uHcPE9J68sqqA5mZ1O45PsaLQgZqblReW5TM2NY6xqXFuxzEmoN18dhq7ymr4Km+/21F8ihWFDrRkm2dltdtt4jtjXHfJyD7ER4Xx0opdbkfxKVYUOtDfluWTGGMrqxnTFYSHBDMzM5VPN5dQXHHY7Tg+w4pCB9lRWs3nW0r5jq2sZkyXcdPENJpVmWfzIbWbfXp1kL8v30VYsK2sZkxXkhYfybRBCcz7uoD6RpsPqT2sKHSAqtoGXs/azeVjkkiIsZXVjOlKbpuSzr7KOt7ILnQ7ik+wotAB3sgu5FB9E7dPtmGoxnQ15w1OYGxqHE8tzqOuscntOF2eFYUzdGQY6oR+PRiV0t3tOMaYo4gI9184mD3lh5mfZa2Fk/FaURCRVBFZLCK5IrJRRO5xto8VkRUiskZEskRkYqvnPCwieSKyRUQu9la2jvTF1lLyy2q4bXK621GMMccxdVAvMvv14KnP8qhtsNbCiXizpdAIPKCqw4CzgbtFZDjwO+CXqjoW+IXzPc6+2cAI4BLgaRHp8rPJPb8sn96x4Vwy0tZMMKarOtJa2FtZy7xVBW7H6dK8VhRUtVhVc5zHVUAukAwocGQF++5AkfP4KmCeqtap6k4gD5hIF5ZXUs2SraV8Z1I/QoOtJ86YruycAfFMyujJU59vt9bCCXTKJ5mIpAPjgJXAvcDvRWQ38N/Aw85hyUDrwcSFzrajzzXH6XbKKi0t9Wbsk/r78nzPMNRJNgzVmK7uSGuhtKqO+TZR3nF5vSiISDSwALhXVSuBu4D7VDUVuA949sihbTz9mJmsVHWuqmaqamZCQoK3Yp9UZW0Db2QXcsWYvvSKtmGoxviCSf3jGZPSnb8v34WqTZTXFq8WBREJxVMQXlbVhc7mW4Ejj1/nn11EhUBqq6en8M+upS7npRW7qKlvsnmOjPExt5yTTl5JNct3lLkdpUvy5ugjwdMKyFXVJ1rtKgLOdR6fD2xzHi8CZotIuIhkAIOAVd7KdyYO1TXy1y93ct6QBEYm2zBUY3zJ5aOTiIsM5cXlNlFeW0K8eO4pwC3AehFZ42z7GfB94EkRCQFqgTkAqrpRROYDm/CMXLpbVbvk1aCXV+7iwKF6fnT+ILejGGNOUURoMDdkpvLXr3ZSXHGYpO7d3I7UpXitKKjqV7R9nQBgwnGe8xjwmLcydYTD9U3MXbKDqYN6MaFfD7fjGGNOw3fO7sfcL3fw6soC7r9oiNtxuhQbR3mKXllVwP7qen78bWslGOOrUntGMn1IIq+s2m0T5R3FisIpOHConj99to1z+sdzVnpPt+MYY87ALef0Y391HR9u3Ot2lC7FisIpeOzdXKpqG3n0yhFuRzHGnKFzByXQLz7SVmY7ihWFdlqWt58FOYXMmdafIX1i3I5jjDlDQUHCrMxUVu48QP7+Q27H6TKsKLRDbUMTj7y1gX7xkXYtwRg/cv2EFIIEu8O5FSsK7fDCsnx27j/Er68aSURol5+jzxjTTr1jI5g+JJE3sgtpbLILzmBF4aTKa+p5anEe5w1JYNpg96bVMMZ4xw1npVJSVcfnW9ydS62rsKJwEk8tzqOqrpGHLh3qdhRjjBdMH5pIr+hwXrMuJMCKwgntPlDDC8t2cd34FIb2iT35E4wxPic0OIjrJiTz2eYSSqpq3Y7jOisKJ/DEx1sRgfsvHOx2FGOMF83KTKWpWVmYs8ftKK6zonAcG/ZU8NaaPdw+JYO+cTY3ijH+bEBCNGel92B+1u6An1LbisJx/PaDzXTvFspd5w1wO4oxphNcOz6FHaWH2FRc6XYUV1lRaMOSraV8uW0/P5w+kO7dQt2OY4zpBBeP6ENwkPDe+mK3o7jKisJRmpuVx9/fTEqPbtxyTj+34xhjOknPqDAmD4jn3XXFAd2FZEXhKG+v3cOm4kp+cvEQwkPsRjVjAsmMUUnkl9UEdBeSFYVWahua+O8PtzIyOZYrRvd1O44xppMd6UJ6d13gdiFZUWjlxeW72FN+mIcvHUZQ0PHWBzLG+KuWLqT1gduFZEXBUVHTwJ8W53Hu4ASmDOzldhxjjEsuG5XErrIaNhYFZheSFQXH05/nUVnbYNNZGBPgWrqQAnQUkhUFYE/5YZ5fls8145IZlmTTWRgTyHo4XUjvBWgXkhUF4A8fbQHgAVvA2xgDXDG6L7vKalizu9ztKJ0u4IvCpqJK3ly9h9snp5Ns01kYY4BLRvUhPCQoIOdC8lpREJFUEVksIrkislFE7nG2vyYia5yvfBFZ42xPF5HDrfY9461srT3+wWZiI0L51/MGdsaPM8b4gNiIUC4e0YdFa4uoa2xyO06nCvHiuRuBB1Q1R0RigGwR+VhVbzhygIj8Aaho9ZztqjrWi5m+4att+1mytZRHZgyje6RNZ2GM+afrJqSwaG0Rn+WWcOmoJLfjdBqvtRRUtVhVc5zHVUAukHxkv4gIMAt41VsZTpKPxz/IJTnOprMwxhzrWwN7kRgTzoKcQrejdKpOuaYgIunAOGBlq81TgX2quq3VtgwRWS0iX4jI1OOca46IZIlIVmnp6S+f92luCRv2VHLvBYNs3WVjzDGCg4RrxiXz+ZZS9lfXuR2n03i9KIhINLAAuFdVW98NciPfbCUUA2mqOg64H3hFRI4ZH6qqc1U1U1UzExJOb81kVeV/P9tGWs9Irh6XfPInGGMC0nUTUmhsVt5eU+R2lE7j1aIgIqF4CsLLqrqw1fYQ4FrgtSPbVLVOVcucx9nAdsArS559vrWUdYUV/Ot5AwgNDvgBWMaY4xjcO4ZRyd1ZkB04XUjeHH0kwLNArqo+cdTuC4DNqlrY6vgEEQl2HvcHBgE7vJHtqc/ySI7rxrXjU7xxemOMH7lufDKbiivZsrfK7Sidwpt/Jk8BbgHObzXMdIazbzbHXmCeBqwTkbXAG8Cdqnqgo0Ot3V1O1q6D3PGtDMJCrJVgjDmxy0b3JUjgnbWB0YXktSGpqvoV0OZUo6p6WxvbFuDpavKq55fuJDo8hFmZ1kowxpxcQkw4kwf04p11RTxw0WA8nSD+K6D+VN5XWcu764uZmZlCTITdl2CMaZ8rxnhmTt2wx/9nTg2oojD/6900Niu3TU53O4oxxodcPKIPocHCO+v8vwspoIrCuj0VDEyIpl98lNtRjDE+JC4yjGmDEvjH2iKam/175tSAKgrbS6sZkBDtdgxjjA+6YkxfiipqySk46HYUrwqYotDQ1ExBWQ39E6yVYIw5dRcM7014SJDfj0IKmKJQcKCGxma1loIx5rREh4dw/tBE3l2/lyY/7kIKmKKwo/QQgLUUjDGn7YoxfdlfXcfKHWVuR/GagCkK20urAehvLQVjzGmaPiSRqLBgvx6FFDBFYUdpNb2iw+neze5PMMacnm5hwVw4vDfvb9hLfWOz23G8ImCKwvbSQwywriNjzBm6YkxfymsaWJq33+0oXhEwRWFHabV1HRljztjUQQnERoT47SikgCgKBw7Vc7CmwVoKxpgzFhYSxKUjk/ho0z5qG/xv/eaAKAo793suMttwVGNMR7h8TBLVdY18vuX0V3/sqgKiKJRWeZbSS4wNdzmJMcYfnNM/nvioML8chRQQReFgTQMAPSLDXE5ijPEHIcFBzBiVxKe5+6iua3Q7TocKkKJQD0BcpA1HNcZ0jGvHJ1Pb0Ox3S3UGRFGoqGkgLCSIbqHBbkcxxviJcWk9GJcWx/NLd/rVzKkBURQO1tTTIzLU71dMMsZ0rjumZJBfVsPiLSVuR+kwAVIUGux6gjGmw10ysg9J3SN49qudbkfpMAFRFMpr6u16gjGmw4UGB3Hr5HSWbS8jt9g/luoMiKJgLQVjjLfMPiuVbqHBPL/UP1oLAVEUPC0FKwrGmI4XFxnGdROSeWtNEfur69yOc8a8VhREJFVEFotIrohsFJF7nO2vicga5ytfRNa0es7DIpInIltE5OKOyKGqlNc00MO6j4wxXnL7lAzqG5t5eUWB21HOWIgXz90IPKCqOSISA2SLyMeqesORA0TkD0CF83g4MBsYAfQFPhGRwap6RpOLVNU10tis1n1kjPGaAQnRTB+SwIsrdnHnef0JD/Hd4e9eaymoarGq5jiPq4BcIPnIfvGMD50FvOpsugqYp6p1qroTyAMmnmmOCuduZrvQbIzxpju+lcH+6jreWVvsdpQz4s2WQgsRSQfGAStbbZ4K7FPVbc73ycCKVvsLaVVEWp1rDjDH+bZaRLa0J8Os355i6H/qBfjnxOnHF4iv2Rvs93jmfO53OPP0P2ta8/br7ne8HV4vCiISDSwA7lXV1mO2buSfrQSAtu4sO+Y2QVWdC8zt0JAnICJZqprZWT+vKwjE1+wN9ns8c4H6O3TzdXu1KIhIKJ6C8LKqLmy1PQS4FpjQ6vBCILXV9ymA/01BaIwxXZg3Rx8J8CyQq6pPHLX7AmCzqraeSWoRMFtEwkUkAxgErPJWPmOMMcfyZkthCnALsL7VsNOfqep7eEYZte46QlU3ish8YBOekUt3n+nIow7SaV1VXUggvmZvsN/jmQvU36Frr1tU/Wd2P2OMMWcmIO5oNsYY0z5WFIwxxrSwonAUEZnpTMvRLCLHHRImIpc403HkichDnZmxo4lITxH5WES2Of/2OM5x+SKy3pmiJKuzc3ZFJ3sfiMf/OvvXich4N3J2de34PZ4nIhWtpsj5hRs5O5KIPCciJSKy4Tj7XXnvWFE41gY8w2WXHO8AEQkGngIuBYYDNzrTdPiqh4BPVXUQ8Knz/fFMV9WxgTh2/GjtfB9cimck3SA8N13+X6eG9AGn8P/Tl857b6yq/qpTQ3rH34BLTrDflfeOFYWjqGquqp7sLumJQJ6q7lDVemAenmk6fNVVwAvO4xeAq92L4lPa8z64Cvi7eqwA4kQkqbODdnH+9v9Tu6jqEuDACQ5x5b1jReH0JAO7W33f5pQcPqS3qhaDZ84qIPE4xynwkYhkO9ONBLr2vA/87b3iDe39HZ0jImtF5H0RGdE50VzlynunU+Y+6mpE5BOgTxu7HlHVt9tzija2demxvSd6zadwmimqWiQiicDHIrLZ+WsnULXnfeBz7xUXtOd3lAP0U9VqEZkBvIWnW8WfufLeCciioKoXnOEpfG5KjhO9ZhHZJyJJqlrsNE/bXIVcVYucf0tE5E08zf5ALgrteR/43HvFBSf9HbWeN01V3xORp0Wkl6r61GR5p8iV9451H52er4FBIpIhImF47tBe5HKmM7EIuNV5fCtwTGtJRKKcdTEQkSjgIjwX5QNZe94Hi4DvOiNJzgYqjnTVmRYn/T2KSB9n6hxEZCKez66yTk/auVx57wRkS+FEROQa4I9AAvCuiKxR1YtFpC/wV1WdoaqNIvJD4EMgGHhOVTe6GPtMPQ7MF5HvAQXATIDWrxnoDbzp/H8ZAryiqh+4lLdLON77QETudPY/A7wHzMCzPkgNcLtbebuqdv4erwfuEpFG4DAwW318OgYReRU4D+glIoXAfwCh4O57x6a5MMYY08K6j4wxxrSwomCMMaaFFQVjjDEtrCgYY4xpYUXBGGNMCysKxrQiItUn2R8nIv/a6vu+IvKG83isc7ftqf7MR0Xk3049rTEdz4qCMacmDmgpCqpapKrXO9+OxTOu3BifZUXBmDaISLSIfCoiOc4aEkdm7XwcGODM6f97EUkXkQ3Onbi/Am5w9t1wdAvAOS7defyIs37AJ8CQVscMEJEPnEkHvxSRoZ33qo2xO5qNOZ5a4BpVrRSRXsAKEVmEZ62Jkao6FuDIh7yq1jsLv2Sq6g+dfY+2dWIRmYBnKodxeP4fzAGynd1zgTtVdZuITAKeBs73yis0pg1WFIxpmwC/EZFpQDOeKYt7d9C5pwJvqmoNgFNsEJFoYDLwujOdCEB4B/1MY9rFioIxbbsZz/xXE1S1QUTygYhTPEcj3+yibf38tuaXCQLKj7RCjHGDXVMwpm3dgRKnIEwH+jnbq4CY4zzn6H35wHgAZ33dDGf7EuAaEenmzDx7BbRMD71TRI5MSCgiMqbjXpIxJ2dFwZi2vQxkikgWnlbDZgBVLQOWOheNf3/UcxYDw49caAYWAD1FZA1wF7DVOUcO8Bqwxjnmy1bnuBn4noisBTYSAMtSmq7FZkk1xhjTwloKxhhjWlhRMMYY08KKgjHGmBZWFIwxxrSwomCMMaaFFQVjjDEtrCgYY4xp8f8Bbz9OB0ckyJkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(0.5 * (lat_bins[1:] + lat_bins[:-1]), zonal_mean)\n", "ax.set_ylabel(\"tas / K\")\n", "ax.set_ylim(270, 305)\n", "\n", "# Scale the x-axis to account for differences in area with latitude.\n", "ax.set_xscale(\"function\", functions=(lambda d: np.sin(d), lambda d: np.arcsin(d)))\n", "ax.set_xlim(np.deg2rad(-80), np.deg2rad(80))\n", "ax.set_xlabel(\"latitude\")" ] }, { "cell_type": "markdown", "id": "2e87d916-0df4-4086-b8a8-b861adbb6edb", "metadata": {}, "source": [ "## Multi-dimensional input\n", "\n", "Until now, we used data at a single time step to illustrate the general idea of calculating zonal means by using histograms.\n", "However, most real-world data has several other dimensions like time or height.\n", "A straight-forward way to calculate zonal means for this kind of data is to unravel it, i.e., to get rid off every dimensions.\n", "This approach, however, will loose all information of the thrown away axes; data at different heights or times will be mixed.\n", "Fortunately, there are alternatives that allow us to calculate zonal means while maintaining the dimensional structure of our dataset.\n", "We achieve this by using `xr.apply_ufunc` which lifts a function (in this case `_compute_varsum`) from (numpy) arrays to (xarray) DataArrays.\n", "This lifting into the world of DataArrays involves describing the dimensions, shapes and data types which the function cares about. Afterwards, xarray applies the usual looping and broadcasting rules over the dimensions the functions does **not** care about. Any necessary looping may then be carried out in parallel (e.g. using dask)." ] }, { "cell_type": "code", "execution_count": 10, "id": "1ce575fa-e642-46a1-8a95-01b86a45a90d", "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "\n", "\n", "def calc_zonal_mean(variable, **kwargs):\n", " \"\"\"Compute a zonal-mean (along `clat`) for multi-dimensional input.\"\"\"\n", " counts_per_bin, bin_edges = np.histogram(variable.clat, **hist_opts)\n", "\n", " def _compute_varsum(var, **kwargs):\n", " \"\"\"Helper function to compute histogram for a single timestep.\"\"\"\n", " varsum_per_bin, _ = np.histogram(variable.clat, weights=var, **kwargs)\n", " return varsum_per_bin\n", "\n", " # For more information see:\n", " # https://docs.xarray.dev/en/stable/generated/xarray.apply_ufunc.html\n", " varsum = xr.apply_ufunc(\n", " _compute_varsum, # function to map\n", " variable, # variables to loop over\n", " kwargs=hist_opts, # keyword arguments passed to the function\n", " input_core_dims=[[\"cell\"]], # dimensions that should not be kept\n", " # Description of the output dataset\n", " dask=\"parallelized\",\n", " vectorize=True,\n", " output_core_dims=[(\"lat\",)],\n", " dask_gufunc_kwargs={\n", " \"output_sizes\": {\"lat\": hist_opts[\"bins\"]},\n", " },\n", " output_dtypes=[\"f8\"],\n", " )\n", "\n", " return varsum / counts_per_bin, bin_edges" ] }, { "cell_type": "markdown", "id": "35e6f478-b217-48c2-89b6-c6d5a74abae4", "metadata": {}, "source": [ "Using this function we can calculate the zonal means along the time dimension." ] }, { "cell_type": "code", "execution_count": 11, "id": "0e6d8689-6154-4498-8fd4-2eb37d343884", "metadata": {}, "outputs": [], "source": [ "zonal_means, lat_bins = calc_zonal_mean(\n", " ds.tas.isel(time=slice(24, None, 48), height_2=0), **hist_opts\n", ")" ] }, { "cell_type": "markdown", "id": "a79344f8-00ca-4fb8-8166-63fb0428b1ba", "metadata": {}, "source": [ "We can now either plot the zonal means for individual timesteps or the whole dataset." ] }, { "cell_type": "code", "execution_count": 12, "id": "3d7a4f19-b6f4-4fdc-bc18-8834b88331e8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 0, 'latitude')" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABkvUlEQVR4nO3dd3zV1f348de5Izd7DzJJAglJWAHC3sgeAop7W7Va22q33d+2v9rWtrZ11TrqFgeishFkb8gmi0zI3nvceX5/5HpDFRCVcCE5z8cjD+7nfM793Pe53Nx3Pp/PGUJKiaIoiqIAaJwdgKIoinLlUElBURRFcVBJQVEURXFQSUFRFEVxUElBURRFcVBJQVEURXHot6QghHAVQhwTQmQKIXKEEL+zl/sLIXYIIQrt//rZy6OFEN1CiAz7z/P9FZuiKIpybqK/xikIIQTgIaXsEELogQPAI8B1QJOU8s9CiMcAPynlz4QQ0cAmKeWofglIURRF+VL9dqYge3XYN/X2HwmsBF6zl78GrOqvGBRFUZSvRtefBxdCaIFUYDjwrJTyqBAiREpZDSClrBZCBJ/1lBghRDrQBvxKSrn/HMd8AHgAwMPDY0JCQkJ/NkFRFGXASU1NbZBSBp1rX79dPvqfFxHCF/gQ+B5wQErpe9a+ZimlnxDCAHhKKRuFEBOAj4CRUsq28x03JSVFnjhxol9jVxRFGWiEEKlSypRz7bssvY+klC3AHmAxUCuECLUHFgrU2esYpZSN9sepQDEQfzniUxRFUXr1Z++jIPsZAkIIN2A+kA9sAO6yV7sL+Pis+lr741ggDijpr/gURVGUL+rPewqhwGv2L3oN8J6UcpMQ4jDwnhDiW8AZ4AZ7/VnA74UQFsAKPCilbOrH+BRFUZTP6bekIKXMAsado7wRuOYc5R8AH/RXPIqiKMqXUyOaFUVRFAeVFBRFURQHlRQURVEUB5UUFEVRFAeVFBRFURQHlRQURVEUB5UUFEVRFAeVFBRFURQHlRQURVEUB5UUFEVRFAeVFBRFURQHlRQURVEUB5UUFEVRFAeVFBRFURSHfl2jWVGUS0tKSVljFyaLDa1GoNMItBpBsLcBg07r7PCUAUAlBUW5CpTUtXPgwG568t7F6pJJp9ZCl1bSpQEbEm+bDVerDq3ZH2ELwkNq8RRaPLRaXJBo3LwJjRhO7LAEPIJjITAedAZnN0u5AqmkoChXsO2p+zi99WWGc5Qynx6ypRsziocSpTcT4drKUEMD4brms55R7XhkkRq6hIEeXPDu7MK10QyZvfuaDGH43fk2IvwL62Apg5xKCopyBWpsO8N7rz7O+LZNHMKb+mI3NM1GYpvKOdFajKcevDxc6fYIpclzCraoBAzRIXj6tBPjW0SAWz1CI9FiBZugRgZR2JhAd4M3/s3t/KRnLZ4vzodFj+My5QEQwtlNVq4QKikoyhWks7OUDz99iehD68nP6+C19HqOlleep7YZaAdOAZ8Q6OvFyJhYjOEJjBgyB1etC9aeNqymJlw1TSwMasI7sAdzcA/3yJ/yWOVa5m3/KR0lB/C8/llw9b58DVWuWEJK6ewYvraUlBR54sQJZ4ehKF+bxdJJbfU+TqRtoulkHjVZ5RzIrWdXfis9lq9/XCEEI8JjSI6fyIj4ZIZofTF2dmM15qOx1jEqMYHn4/2ZkFfMTzTvYvIMw/26Z2DY3EvXOOWKJYRIlVKmnHOfSgqKcvk0NDeQkbuJvPLjdNVUUXagjLKiSgrKWyhtMp/zOUKjIcTHn0B3A2HeIYS7+SJkO922ejpsLbR295BW3kOH8fyvmxAezcyZ8xkXPRwagunu2ENsuCf7p+kpy47mCevzDNNUYxlzK7rFfwR3/356B5QrgUoKyqBn6+ykdNtWis800tgmMfcIdCaBi82CGRPtsot6Yxsdli4sNhsWKbHYLOgMegJiXVmybBnjR85BfI1r72armc3p71CYtZVTuZJTB0s4U3aKspqWCz4vNCyMWaGxjB8SRqDBC43Wh3bvYMxaFzS2TjB1oNW1YonSoHczUtBZTdnhM5zJKSa7qptz/WYPGxLJ3LmzmeQ5iY7OPAK8W2lbauHfp2/hOzUf8G3dJmyufris/BckrvjKbVWuDk5JCkIIV2AfYKD33sU6KeVvhRD+wLtANFAG3CilbLY/5+fAtwAr8H0p5fYLvYZKCsq5NHQY2Z5dSeuu9/AvLsLV5kaxqOVodzUNbV20t3XT3tZFV1snHe1dtLcbsZit5z2ewd1A5NAhjBgSR1hEBD5D3PH0BT83DZG+XkybdwdBkSPQaPrGgpqtZj7O2UR+2uuE9XSweZOWfanZNNRUnfd1tHotkbHBjAmLZJzLcJIMroQtXUzAigXUtRipru2kobGLpvoejOXNeNdX4mquBVs7wtCDHFKF8PWicEQA+1uHEvjWWzRnniSzuhnr537NoyPCeWjmvbjorGg5SeA8HS94LsPlpIm/2p5njKYU6/Qfob3mV6BRY1wHGmclBQF4SCk7hBB64ADwCHAd0CSl/LMQ4jHAT0r5MyFEErAWmASEATuBeCnleX9bVVJQPtPeY+ajjCr2p2WSWLyHrpJqTtSXk1p/msqSGmyWb3CB/hw8DHrCfLwI9fEhMTiOmNBI/KK0dAcacaeZnjojR6q9SK/upCTtIMaOti8cQ6MRJMT7kRIdQYx7HJEmP0IRxET5E/3tW3EbM+aCMdhsNooLm0k7XEX10QK0nel4u1rpDstEugfSFSPIDogn+6Qv3s+8wMGSXMzWvl+nQH9ffrrw27i4hmDt2oN7WDv5CyawsWQ+/6/m39yi2405cgH6218Bg9clff8U53L65SMhhDu9SeEh4HVgjpSyWggRCuyRUo6wnyUgpfyT/Tnbgf+TUh4+33FVUlCqWrp55WApJ0/uYmRGJkdKCjlUlE5z0xe/hC9Ep9Pg4e6Gi8EFjUaDRiMQQtDW0k5HV89FHUOv0zAkNIT29k5aWs79+i56DdNHDWNsUBwjDCFEB0YyavEkglKScImMQOj1Xynuz1hMVtL3VXJ8czbmpjSiXSrRhp3ktEcIXqFmMsKGcfLoEGyvf8yewqNYbTYAQoL8+dmC29F7LMLWk4bZdhzLLAOv6u/khqJP+I3+DcwuQzHc/T6asLivFZty5XFaUhBCaIFUYDjwrP2MoEVK6XtWnWYppZ8Q4hngiJTyTXv5y8BWKeW6zx3zAeABgKioqAmnT5/ut/iVK9fJylZe3F9CdfYevI/uY3vOESpqas5bPzrAQEhEFMbQOFoCo2kLiSIgWE+0ZxeRXkaGepgI1mtw0xvQ6AxIjQ5p1VJSWMWpmk7yT7dQmV9K55lyjE0NWE2mrxRvoI87y0fMZk6QL/Ez5zP2vhtxD/D8hu/CF9lskpL0OvZvyKXtdCox+hOEhuVwwjMZn7gz7PaeQ9MzRXz88TuOxBAWGshPb04koPJ6WnwTsBmPYfY/wSfjlxCQ38azun/hhkQz7xlcZl13yWNWLr8r4UzBF/gQ+B5w4DxJ4Vng8OeSwhYp5QfnO646UxhcbDbJnlN1vLivmLCj75KbmcnenCN093yx242vK0xPCMZ12ETKRi1DxsUwOkrPRD8fJnh7MdrLEw+DG7j6XtQ1c5vNxt69G0nfthZdbTOyWYuxzZXmLh05zdVk1OVR3troqK/TaokPjmKE/3DGeYcwyaeK0Y/9kbDZUy7lW3JBLXVd7HjrOLUZxxnlvZ2mIUZaAyKxxday/WkP1m3eyGe//1Fhofz0Z/4E7YjB6LWKVs9Qeqx7KBynIbduBM+LfzBSnKYn/AEMdz/+tc9olCuD05OCPYjfAp3A/ajLR8oFGI1GcvPyycguwGo1YrZ10mNtpaUrj7bddezKLCWrLP8Lz3PVwdwRvviPmkpN8g3MnjKFOyfEEOVp+Fq9hi5Ge30de576E9bcg5hrXKkwBxEgTUTqa9H4WNBHxpH04E/xHT+hX17/YphNFj760wcYSw8TH7WHDI8ZWMaXse8JI+/t3OnopTTEP5gf/CIaj0Ij8QXjKR62GrPoojU8jXe1E/m18TVu1O2jU5+C4c7X0UWGO61NyjfjrBvNQYBZStkihHADPgH+AswGGs+60ewvpfypEGIk8DZ9N5o/BeLUjebBwWw2k56ewY7cQkraCmg4UgpF1Zh7ejCaTXSbjRTWl9PQ1vSF5w710zF+YjLmmTcwfeICvn/NGNydNGNoV2UNhkA/tIYrb7K5+so6tv/sGfxDjlAZGUPH6HqO/L6B9/btd9TxcvXgoe9NICi8g3HvaWgOvp6a4IlI1wxeDwthflsav9O/hg1PGP8bDMvvQajeSVcdZyWFMcBrgJbedRvek1L+XggRALwHRAFngBuklE325/wSuBewAI9KKbde6DVUUri6tbS0cPzQMXacLueglw+WgpOYdx2hNOMQLV3tF3yuAMbHh+A67zpCp8/jT0sWMDzA5/IEfpXb+9rHNBY9TUtkENWxgqI/VfHW3oOYrb09tHQaLbfeNpGZ13XB64LE6mhykm7HqNXxaUI55noNf3J9iTGyBKP7RPS3PYMmPMHJrVK+iivi8lF/UEnhymaz2ejo6KC1tZXW1lYaG5soKT9NpbGGEi8DJz2iaGow0711H7rdH3OmrvxLjxngrmHsxNF0rbyT5eOT+e60Wfjo1RReX9W+19ZxuvxVRLSV/JAw2l4u4bUNqbR19yXjVStTuOPhHvZ+Es6d22spHLGGxoCpFA7NYHNXLDfbdvNTl3dxl2bkmAfRLXtMdV29SqikoFw2UkrKy8tJO/ARB9s7KfcOpkbvRaPRg3aTHh/PJjpN0LH7JEP2bCanNBuL7YtXCMO8NIwfHYtLUAi4uYO7B9LHC7eJCTw8dRHTop13jX6gOPDiK6TV7SUsoYg070S6N1aw7s1MKpr7pt9evWoK93+nhbeyx/HdN9Pp8B5DScIttOklW6OaaGt05Reub7OaA1h0vmgW/hrNhLtBqxL1lUwlBaXf2Ww2jhw9xoas42QEBVLQE8qsg9sJ7K6kOd6fyjMW2o7mU12Wz+mmc8/6adDC0gQ3AifMpmbRYqb5pRGCG8OC40iIHI+Pz0gMhuDL3LKBbf8zf+FQazXDEk6Q6ZdIw4Eedj+3n4Lavq7e1103k/seaOS/NTO4/8X9BLW5U54yn3KX+ZQOKWWjCCSh5wy/9FjLJEseRo9IXG56GRE12YktUy5EJQWlX7W0tPDkujf5MDKO5kbJ9Ydexb9dT2pRC0WlxRTVFWG2nnuyN4BJ4RpWjvWnadwaqkZHco3LfkYGr2LytPsRQt3E7G8Hnvklu1vbGRGfycmAOEqPuHLknxspqj/jqHP99XO57756/tsxj5Wv7GJ8sYW60CBqR6ykUoxkT2wZmU3hzNdn8Dv9awTYOjF87xDCP8aJLVPORyUFpV/YbDa2bdnEe13F7NbHsvrYM/ie9mdXeimHT6VjsZ17agmtBsaHG1gWC2OTotkXuxo5vJMU80n0pmtYtupRPDw8LnNrBrfD6/7OlsJCRsacIisonpr9ZvY8u5PSht7EIBDccONC7r2rhjc11yCOlPKtrXkEt0L5iBiqI1ZT5A2bPDzx7uxig+uv0Rr88PzhMXBR/5dXGpUUlEuurq6OV975B1sTEtAUHmTh/ho2HC/iaGkBNmn7Qv24QC3R8TFcP7SbW4a2UeYSy6tBixgam0dQazetNSksvPY+hg0b5oTWKABHtz3DBxlpJA89w7bAORh25/DJC4c409ibGDRCwy23XctdN1Zyyi2I19qWMf/T91izvxa9RdI5MoysYffwgX8XIS2NvOLyBN2hM/F8YINa2e0Ko5KCcslIKTmx7whbsz7k/egoFh14l+5Del45cohuy/+OLJ4crmX+1GisY2dzi0xjjLWEQhnO2oDZREWeglYPGqtGMnPOalKmpPzPLKOKc+zf+iQfZaUzJb6EfxseYMTeDWx88RiVzRUAaLVa7rzzBm5dXo3Jp4bNzWs4YArigY9fZk6mGelqo2TSVH4etZzbTJ/yC83bdE58FI9lv3Nyy5SzqaSgXBLdnV2s/++7FMkcPgly4dsfV/Hvfds5Uv2/cw5NifFgzsoUqsct5cEzHzLFkku5DGKd93R0Xj3QHIVbRyTDJiQybcksDFfgQK/B7KN3f0JmSTNJY7P5vek3zNz7AutePkJ9ex0ALnoX7rr7Vq6fV4cuOJ+K+jj+n/tDDC/8D9/aVsrIM7BtxWL+pZ/H027PsdxyCPN1r+IydrWTW6Z8RiUF5RsryShg/ccf0e3TQIFHG1M2a/j1lpdpMffdQPb29+OW++dSM3oR3y9byzxzGrXSlw366bhYk4kxxuJrcydq9Sh8J0c5sTXKBUnJ88/MwdaRiGdKLr9t+B3X7vs7r71+kNauFgDcXN24/Y7bWDG9BffINCyNfjxu+yEV+izu2rKBNQckT3z7Ho7UDucjj/9juLUK3bd3QeiFpwNXLg+VFJSvzWqxsvOtzRwpScPTqxtogCNu/OrD5+k5a3xB0vSJhNz9AD+oWccKy2GapQc75UKijTfiYdaj9TQTumQkXpMiEVp1mehKZ+ts5vcvTiXONoPm5Fr+UfEI1376FC++u4tOYyfQe8aw5oY1rJihJyjuIMKo553GW/jUU/KX59cytEHHd2/9GbZG2Gb4Bd4GLfrvHAAfNWeSs6mkoHwtbTVNrH35LarNjSR5nMG9cwhHc2t44oPXMNmnXfbyciH5rnuYFeHJr7tewYKWg8zDo+NmulvaiZ5gIPaeFWjc3JzcGuWrasrbzU8+/Q5LXcdQEguvlN7Bom2v8Z8Nm+k2dQOgEYKly5axdGYsI0YdQLq0sat+IYeai3nihdOcHhnNz+IfJs5Szsfi1+iCh6O9f4ca+exkF0oK6k825Zy669t54z+vUW9qZLn2GKEtU9mef4q/rHvVkRB8/QwkfffX/NivjP/X/SIF2mFktv+LjvLrcTOdZM5flzH8OzeqhHCV8k+cyyMR17Kh+TTx9U3cGP0BB69ZxaPX3oS/ZwAANinZtGkT73ycysHU+RibhzIvZBshoeN4cbGOkallrAnaTx5D+S6PIhryke/fA9ZLuxKecumopKB8gbmhi3eefQOjrZx7ZT41PTfyj4x3+Pd76zHbzywDAlwZf/fPeEv3Jss0R0jTXENJ9c/oLDzBvOs9Gf/0r9D5+zu5Jco3Nebax7nVVUtmjg8TO08yI/Y4h2dN49GVc4gM6BuYtu/QQTa9s4WjWVOprxnOt7zf58CU69mZLLj1PxsZMbqST8zj+LvpVkTRDtj/Nye2SrkQlRSU/9FS3sb6f71DEyWstJr40DqR//fJP/lwyy6s9oQQHOjBnPse5kOvF/DTtZNte4CSkiUMb9vB3P8+hu+yJU5uhXLJaLQsvHs9sT4ZdB+ZxSL5CSPiytk0+3oeXjmU0eHJjqrHck+y5ZU3OJMzHrPRjd9rtrF2UTSNHpKHDq5FM8qNZ23L2GUch9z7Fzhz1HntUs5LJQXFoam2g0PPb+KMNo815hbe7Qnk72//i6MnMh11hkb4sfJ7t/KS4W3ate7k93yH01kapgSXkPTKU+iCgpzYAqU/CK8QVl37GqXB6zEcX8P1ureZOiSDt+Y+yp0rdMwY2jfH0YnK07z58psU5EzEz7Oekf5xPH2tjuE7KlntsQ3zcC++Lx+m0eaDXH8/9Hy1tbSV/qeSggKAzWoj45k95OiyuNFSzVvdEfzr1acoLeubznrC6Gge+OV8/mzdTKfGlTSP1bQfz2PG8uGE/+VxNC4uTmyB0p8846azJv47bPXYhPvJG1ju/gGrfbbw0jV/5Npru1kSM8lRN7u2gmf/sYPy8qHc5LGFssTVvD1Hx63PbiI8poHOQC8eMn0PWs7A1p86sVXKuaikoABw4KV0SsjjBmsJ7/ck8MxbT1PX2OzYv2JRMg/9Ip4Hqg9jEnrWDl2N7bCFqT+6g4D77uu35S6VK0fMkh9zoz6Sj00ZuGfewQy3Xdzv+l9envMkM2/o4cZhfdOZF1dX8PZ/ejBbXPi12MG+qUMp0Uruy3mNnlEBpGtH8JJpKWSuhex1TmyV8nkqKSgUHK3CerqYUZxgg2kCz659iuq63mUvhRDcd9d0ljwYwbK8AmwIfhn/beL2uTPrjmvwXrzIydErl40QjLvvLW6xVLDReBzPtPsYbUjnh4a/83bKXxh/h4Vb40Y6qn965AgHdsXj41nHaJ94XlyuJ/ndfJZqttI1NpA/227mlG0ocvMPwXjhlfaUy0clhUGuo6Ub64c5VGv2UdEzk+fefZKyqnrH/nvunMTwNZEsTj+FOz3cN+xXrDhpZeLcEfitXOHEyBVn0Lh6Ev/AZu6yFrHdfACvEw8xVFfKzzz/zPuJ/4+xd+pI9A8EerurvvHKYWqqI7nefTOd4dfy7Fwd9731LhEB1dii3Pm5+W5ETyvkbXJyy5TPqKQwyOU+k0aeJpd4awx/2/pPCs7UOfbddnMKITcns/xEDmE08ED4L7ixuZFYfy/C7rjdiVErzuQRGEXY/du53XqKvabdeJ74DoGaWn7q8ThbEv/Aslv90Gu0AJQ31LL+vzqsFj2P6A9QGRfJbq2N75U+iznem2xDPGdsQZhT33Byq5TPqKQwiJWm16DvaCZMVPJq/hGO5hQ79l137STC75nMdUdSSZRn+H7AD5klO4htk4z58SNOjFq5EniHDCXwge3cSD5pPbvQpt6Lj76JR11+z+mJv2ReSt8YlY1793LscCKB7uVEByezdYYejx0V3CjX0jEqkPW2mejKD0HruVfkUy4vlRQGsY73sinTHqa8NZQ3tq13lE+ekETyd0cw80gRKZziN173E+WtIUXjz7Rf/8KJEStXEt8h0fjct5UV2hzONGfQnX4bboZ2Fvq+w4jVi4jycgfAbLXw1itp1NeHc5PuQ2yBt/O3OVqWb93CCN9C1ouZCCRkv+fkFimgksKglffxKWqoIMo6mn9u/T/aenqnrvD39WHVz4fhnuXBEusR3tfNoTsihHm+kUy6/0EnR61caQLCY3G5ZxPXuB+lq+IMjVnXEeVWQsPwMSy7NQxBb6+0U+VlfPiuJ9Km5SH9bgK8JvLXcB2rmz6kZEgsqbY4zMffgKt4LraBot+SghAiUgixWwiRJ4TIEUI8Yi8fK4Q4LITIFkJsFEJ428ujhRDdQogM+8/z/RXbYGfpMmE6kI9eU807mX/ncEnvACIB3P1YEmeaJ3FD2w6qRCAfJi5kSeAQpq5S9xCUcwsbGk/3zR8xzX83hqIuaoqmc6v/KxSN/xlTx/k66m3cspu0jCRCDacICwnH28+L+pICdCGw3joTfWsxVGee/4WUy6I/zxQswI+klInAFOBhIUQS8BLwmJRyNPAh8JOznlMspUy2/6g/S/tJwYvpnHEpoL1By7M7+2aZXbhsGBUxt7HkzEHCaeDXwx5iqbab+fNVQlAuLDFxFBXL32N80Ca8cuLoafFlztC9jFu9iED33kGN3cYe3n41l6baYdwg3qVnyA2c8JeM9spjk3YqZqmDrHed3BKl35KClLJaSplmf9wO5AHhwAhgn73aDuD6/opB+aKuomYaK4sYagvliU//TLd9ssrwUD8C77uDCRn5rNAc4jm/NSTo6rh91Y+dG7By1ZgxcTypM18hwvN1jDmLSdEe4WTiChbeFOaok52fx4cbfbB1BvItwxsM8xxLVF0mjUOC+dQ2DkvaWjWDqpNdlnsKQohoYBxwFDgJXGvfdQMQeVbVGCFEuhBirxBi5nmO9YAQ4oQQ4kR9ff25qijnIc1Wil4+iFlfwfaCrRwobgF6B6it/OVkEo4auU+/nkzdMIrCg/nZzX92bsDKVWfNwpmkh92Ha3kx1WXjecD/GSom/5pRSb6OOps3budIWgruVhOJQa009KRhCzHwgXUmOlMzFO9yXgOU/k8KQghP4APgUSllG3AvvZeSUgEvwGSvWg1ESSnHAT8E3v7sfsPZpJQvSClTpJQpQWryta+kdsMpqmQ5AV3x/GNnX0+PeYsj0HTOJll/Ak/RzXMjVvPnlb9UU1coX8ut334EQ1AlhpwJeHSYGBVdSMoNs/GyX0Zq62hn86ajVOdfQ4K2gECbK9Ee5ezTj6VNeiIz1jq5BYNbvyYFIYSe3oTwlpRyPYCUMl9KuVBKOQFYCxTby41Sykb741R7eXx/xjeYmCo7yD94jAAt/OPI/1Hd3nuK7uNtwO26hxjSXsNyzVGeD1jD92avxs3Fw8kRK1crT4MO/6V/JNK2jsa8eazQrSd73LeYsirEUedY6gl2Z3RiPT2N2AAPfDrS6QzxZZt1Atb87eoSkhP1Z+8jAbwM5EkpnzyrPNj+rwb4FfC8fTtICKG1P44F4oCS/opvMJFWGznP76DRvZKGGgtvHc1y7Jv1wBRSSszcpdlIoSacyvHTGBU8yonRKgPB1OmTKfYdh2uxluqyUTzo9RRds3/JsHhfR51NH28g9XQY8T7luLWnYg11Y7ctGZ2tEyrVMrvO0p9nCtOBO4B5Z3UzXQrcIoQ4BeQDVcAr9vqzgCwhRCawDnhQStnUj/ENGnXbC6g0NRBnGcEvd/wOq70reHxSAEkuS0h0OUGEaOCPcd/i8Wl3OTdYZcBY/aO/E+O/G5k/BZ+mHuJjaxh70yQMBh0A9U2NbN99BE1TBMON4O3RxhGXBKxSgzy1w8nRD1792fvogJRSSCnHnNXNdIuU8l9Synj7z2NS9o5WkVJ+IKUcKaUcK6UcL6Xc2F+xDSY2k5WMrXvoNtTyxokNpJc3AqDVCAJWPoqnvorr5T5e91rK7fNuQqvuIyiXiJe3J41xdzCkLY2KvDlca95I0aTvMmVu32WkAwcOsD/dnYgwNzyN6TQFBpMhh9GTpSbIcxY1onmAK/34GMU+jXhXh/LU3jcd5fEr57JUY2Gl3EOD8GXHxDuYHxLhxEiVgWj5tx6mxaMbjxpXSvImca/ry9hWX09chC8AVpuNN988RpSmkmGNGZgD3NltTcatLR86VO9CZ1BJYQCTNknG/sPEmoP54dbHMVl6p7IIiQxkcfwSgrxyGS3KeCLsTv42eYGTo1UGqkWP/YehLmno6iPwLDXgFRvPqGWBaHW9Xz+VlZWsf8dMvNmE1hf22sb2PrFkt/OCHsRUUhjATm05RJVfPZ8eO0ZuZe89e61GQ+SKHxEWUsjSrjRyNUMJn3cXIQa9k6NVBqqA4BBqx87Ds7GDyopErmvbTcG8m1kzpu/MdOPmMoJCLQTa8sh1j6ZFemAt+MSJUQ9eKikMYIe370RTBf/a/5qjLGjhDdwU1kqYpZxIUc8/4+7ne7GxToxSGQzuvPth2jxr8WnVUFcQw3g/PZolvrh79I5daGlpozOzgqjmkxgDvNhrG4ut4BOw2Zwc+eCjksIAlblpFx2+bbz96XbM9j7f7hEx3Js8jsCoVJZ2ZLBHn8xtC+9Cr1E3l5X+t+Lnj+NtNKHp8mZqcRWpE1YwLTzAsf+TTzqI7W7A5m9gtzUZvaUVatQEeZebSgoDkJSSg9s+oTy3mkOl6b2FQrB4we3EjTiITwV408WGlO8zO+ALg8YVpV8MCwrhzJRIPBsFdTUhXKO14TLezbH/RGoFER7deHo2sd9mHytT9KmToh28VFIYgPa/v55O707e27PPURYxdi5zkivQ6VtY2pPK++7z+MFstcaycnn98vb7adNJAo2uRJyEhmUL8PXuTQydnZ1UZNUSYkynwSuQAiKw5Kv7CpebSgoDjM1mI3X/ETIP5VDW1Lu8ocbFjfuXJhA9/DjDT5kwCj21C39BpKuLk6NVBhudVsO4B9agadRhNboyVvoyMzzQsf/wnhqGtxRh9Xdlp2U82qrj0NPqxIgHH5UUBpgNz79EtUs5G44cc5QtnzGH5MlHKM0ZwxSRx0sxD/Lg2PFOjFIZzBaPSqQ0zJfwTgMhZ3QEJ/VdQsrIKCfK2ILw07DbmozABiV7nBfsIKSSwgDS09PNyeLTpH9yinZjJwCePsGsut1MYVc0d3bs5oQ+iXtv+QUuGvVfrzjPjx+5jc5WT3yNnrTPSmRIgCcARqORmqIqvF2LSZdxdOIKBducHO3gor4ZBpC3nnia8q6D7M7tm/Du3iVzcA1vIy6/EzdhxP+Wf+Proi4bKc4V4ueOaWISIW06PD2GsCCg7xLSiX1VxHamYvZx5YAYiS1vi5o19TJSSWGAaGlqoshcQurmGqw2KwBx4SOYfG09B04u4loOk5f0ILGxyc4NVFHsvn/PLFq7PImpCCZymM5RnpVZQ2xdMTZ/AxtM09CYWqD8qPMCHWRUUhggnnn8X1SezCe1vMBeInhg1TXs1k7mZy1vU+ESRfJ1v3VqjIpyNr1ei37aKMJ7/OmYGEHsEF8ALBYrLTnVGLy62GNNxoIWCrY4N9hBRCWFAaCmqopG36Ps3VvgKFucNAPdzAZic5uJ0tQTcN2ToFOXjZQryy2rJ+LaqcUYHMpSj75LSPkna/HTpdOhcSNHOxRrzkaQ0omRDh4qKQwATz/1BNWbLZxprgHARe/Kratm8lHdYr4rP6I2fC5uCWrCO+XK4+VjoN4E0Q0jCB9idZTnnGxidGsmtgBXPrZORdt2GurznRjp4KGSwlXuVF4RZksWWzKOOMoemLSaLUk+fLv2Y9yEmZDr/u7ECBXlwgJmTSOqOwDjuAB8vVwB6Ow04lVcCkFaNvVM662Yv9mJUQ4eKilc5d59/zGytplo7+ntghrsFUDCqmRaCr24SbcHOel+CBjm5CgV5fxWLR+FoUNLW1goUwL65kKqLm7BxSOLOuFHqTYES9bHToxy8FBJ4Sq2Y8sWmops7MjrO0v48dx7eM4Qx29tr2Fx8UE/9zEnRqgoX87D20CdTUOwKYUx+r7JGU9mNzOzbTdWPxc22yaja8iCtmonRjo4qKRwlZI2G9npT7BvfyE2exfU0UOGcXzZCOZWpzJFk4fL/F+Bm69zA1WUixAyexrDu0LwG+XlKCssamJMfTUiEDb0TO8tPLXVSREOHiopXKXeevlxivZHkFZ2EgCB4OZVD5BX6cNvXV7HNiQZUu51bpCKcpEWLY5H16HFGB/E0BAfAKxWG/V13bi7p3NKRlCv9cKYut7JkQ58Kilchdo6WmkrPcSHaXscZUsTZ/D68ER+bnubANGOZuVToNE6L0hF+Qo8fAxUS4HeayQzPPrOFtKy21jQvgObjwvb5UT01YfUBHn9TCWFq9CmF37C/n0e1NT3zoJq0LmgvfkuRjWc4lbtLjRTvgOhY50cpaJ8NZHzpjC8bTTD+3ICWZmtjOoxQoCF97vnoMEKOR85LcbBoN+SghAiUgixWwiRJ4TIEUI8Yi8fK4Q4LITIFkJsFEJ4n/WcnwshioQQBUKIRf0V29Usq+QkbSVmPsrsm2d+zbQ15JqD+JPrS9i8w2HuL5wYoaJ8PXMWDMfQrcNrjC86Xe9XU01tMxZhxtfjBJlyGNU6H3r2v+jkSAe2/jxTsAA/klImAlOAh4UQScBLwGNSytHAh8BPAOz7bgZGAouB54QQ6vrHWaxSUvHWb3jncCU9HW0ABHn6kb1wFQ/Kj4mTVWiW/xNcPJwbqKJ8DR4+Blp6THRFBzImwt9RnpHZydK2Hdg89bzLXFxbsqGhyImRDmz9lhSklNVSyjT743YgDwgHRgCfLQm2A7je/ngl8I6U0iilLAWKgEn9Fd/VaMuhTWQVDGVf5m5H2YzF3yKqo4ZHdOuRo9dA/EInRqgo34wcHo2XYQLT9X1TsqSn9jDcpgf/bt7qmI8NsBx5zXlBDnCX5Z6CECIaGAccBU4C19p33QBE2h+HA+VnPa3CXqYANikxbd/AK6kHkPZphEeExJA3Ygr/NDyHzd0PsfRvTo5SUb6Z+WumEdkYS2hwX1I4lV+FWxCEuB2jHn9y3SKwpL4KNuv5D6R8bf2eFIQQnsAHwKNSyjbgXnovJaUCXoDps6rnePoXZsASQjwghDghhDhRX1/fX2FfcT4+mc6BLFdO5Z9wlPmu/A6PibUMowqXNa+Am58TI1SUby482hdbpw2XUV54eRkAaO/upCC3m5XtnyBdNbwml+Aq27AV7HJytAPTeZOCEMLnAvsmXszBhRB6ehPCW1LK9QBSynwp5UIp5QRgLVBsr15B31kDQARQ9fljSilfkFKmSClTgoKCLiaMq56Uksp1H/Lm8U2OsvGJ0xkRYOIu7Q56xt0BsbOdGKGiXBpCI2i29WCOCGRmZN8fOUf2GIlwCUT4N/Fx6wyMWg0tm//hxEgHrgudKXwqhPjCn55CiIXAl44gEUII4GUgT0r55FnlwfZ/NcCvgOftuzYANwshDEKIGCAOOIbCpuJiso7U0VBVBoBeq8O28Eb+6vIfmjyCcVWXjZQBJGLqRPzkRObq9Y6ytPRS3EI6GaXdhUnq2ec5Gp+Ow8iOZidGOjBdKCn8B9gthHD8OS6EuNVevuwijj0duAOYJ4TIsP8sBW4RQpwC8uk9E3gFQEqZA7wH5ALbgIellOqiIZD25ibWZ/UtMhI1eTn/8noPX9GB923vg97VidEpyqU1c/UkAutC8Yv2wM2z995CY1sLuZlmrjemIl3gect1aLHRuP4ZJ0c78OjOt0NK+aIQogfYZT87uAl4EJgrpSz7sgNLKQ9w7vsEAP86z3P+CPzxy449mGwvzOfMwWIaayoA0OtcuH2SJ7O0u6hKeZiwsGTnBqgol5iXnxvdnSZao12YkxDI1hO9V5EPfWpm2v0uGGwVnGiOo8nHDVvpO0j5K3ovTCiXwgVvNEsp3wB+D6QDtwLTLyYhKJfOjrUH+fjkB47tuMkz+LXPRk75xxG2VOVPZWDq0VsQoYFce9ZQpbSMYrrdrMyXnyCssNVtKkGygtasHCdGOvBc6EZzthAiC/gN4A4E0Hs56bNypZ+VNVbTuD+dxpre6YJ1Lq68P/UUtXo/Qm75ANRfR8oANXbpfPy7JhDh6oabd+/l0aaOVvYfHMJyQwlSa+N18zIEUP6JuoR0KZ338hGw/LJFoZzTM29vY2POe47t1ZPCGO7RTPa4HzEhaKgTI1OU/pW8ZDylW/ZQNFzLDFswOw6eAeDkgQ6un9xFQEcR+c3xNHh64m/aibnbhN5NrUF+KZz3TEFKefpCP5czyMHIZOmhcdsJGqt7x2K4GAw8P62e/QHJJC9SC+coA5veRUd7Vwcdw924Ttc3XOlYbipZLWNYwT6EGTa6zyLUXE/2lh1OjHZgUbOkXqGe27iOT3P7lh98MMVAo2cEkQv+hlarpoRSBj43H3f89KOIM7vh5usOQEt3J6VHAlkakAMaK+9Zr0EDdBa+4NxgBxCVFK5AUkryPjhJeVnv1Ng6rYafTIassGuITxjv5OgU5fKYfusy/AuHkh2jZcqEUEd5RloGZ7TRDHXNJbcxlHqDN7HadMpzSp0Y7cBxoRvNLwghVgshvM5XR+kf+zO2cDitr8fR3WO0HAmeweyVv3ViVIpyeUVNHUVrNzQN17PGYHOUHziVTlX5RJbqjyHM8LH7HEK7Wjnx6QcXOJpysS50pvBfYCywRQjxqRDiZ0IItXLLZfDOq8fJzuubGvimKWG4JtxKYGCgE6NSlMtLCIG3iwtuvuGMbdXjGdo7nXaXyUjR0TKm+pQitEY2WaaikeAvPqCzrcPJUV/9LnSj+YiU8v+klDOBG4EzwI+EEOlCiP8KIW68bFEOIpVVxzi5r28WkRXxOjJjF3LNwuucGJWiOMfkmxcTUJRAZqiWa+ZFOcq3pH5KQfdkRrqcIKMxkga9DyNsJRzdvsd5wQ4QF3VPQUrZKKVcK6W8U0o5DniW3rmJlEvs6f9u4VBOrmN75dQYxibfgJubmxOjUhTniL1mAuYWH07HabjP3YjQ9/aiL2mopf6kN0sNGWDTsMF9FkGtRipq12Gz2S58UOWCvtaNZillqn1KCuUS6uwoI3PTNqzm3imfJoVrODVyNbNmqRlQlcFJaAQBri5ohvgSmmMjaWbf+JxdRz8lxtCFXtfGVtMktDZJoscRCk/kOTHiq5/qfXQF+e8b/2HPyQzH9uKpCSwcvxwXFzUoRxm8Jt14Db4Vw8gO1HBHct8Z84HCbHLrpjLB5SgnWqKp1gcQ3d5I2ol1Toz26qeSwhXCZGph73vH6ek0AjDMT0PF5DuYNWuKkyNTFOeKnj8R9+o4MuMF15RY8IvunbjZaDGTc7yQJS65SKnlQ/c5+DebcfHeTUttk5Ojvnp9aVIQQtzwWbdUIcSvhBDrhRCqs/wltumjv7I7u2/5iNlTRnJbyhz0Z80pryiDkUYjCHJzoSfaFUok1y7sW+ZlZ9o+9FZv3PU1bOyahFZKEmURx7ZvucARlQu5mDOFX0sp24UQM4BFwGvAv/s3rMHFZjOz/vUsmhrbAQhwEzTPuofZMy9qgTtFGfBSVs3CtzqKw4kabnD1ROvae0m1vLmejMJQpuqPkdcewSnXKIIbTNRaPnbcm1O+motJCp+9s8uAf0spPwbURe5L6MjB/7I/87Bje9LEUTw6dYaazkJR7GIWTiKkYhJHR2mIOtjJpLlDHPtSU4+xTFcMSD7SzSSw1UiEfzY5u9XCjV/HxSSFSiHEf+gdq7BFCGG4yOcpF+ntf2/hTEUjAK46kPPuZMbMCU6OSlGuHFqdhnBXd9x9PKh3Eaweedb6zacyqe2KJtotg43NKQBEtLaSW7rWWeFe1S7my/1GYDuwWErZAvgDP+nPoAaTirLDHDiR4dgeNSaJX8+Zi0aj8q6inG3CjXMILkph92gNs4o1eA/xBnpHOJ/Mt3Gz9gDl1hBSDQn41QgMIUepPVXh5KivPl/6zSOl7JJSrgdahRBRgJ7e9ZWVS+DVP/+dzKIzjm3Pa25jysxxToxIUa5MMXPHEN2aTEGSBo+8Tq6Z5evYdyz9CFEC3Fyq2WCeQmB3J6HaGo7ue8d5AV+lLqb30bVCiEKgFNhr/3drfwc2GHR21LHrYB7Yp4uPjYvlT4sXqbMERTkHIQTjxrkT2BlBdrSWpSFhjn0ZZadIa01iuctONndNxoKGwBobXe6f0NPa5cSorz4X8+3zB2AKcEpKGQPMBw72a1SDxJuP/4IDp045tkPmrWLyXNXbV1HOZ+yDNzGsbCF7RwvG5fQQMaJ3kjyLzUphdgVLOE2TxpWduomEVRsJ9i/m+Fb1N+xXcTFJwSylbAQ0QgiNlHI3kNy/YQ18NquV7duzMJssAASFBPHEyusQat1lRTkvnUHHaEM3rUMN0NDDomnejn1Hs05QLocz2m0/z3auwN1mJLKmi/Kuddisaj6ki3UxSaFFCOEJ7APeEkL8C7B82ZOEEJFCiN1CiDwhRI4Q4hF7ebIQ4ogQIkMIcUIIMcleHi2E6LaXZwghnv8mDbvSbfzP39lV1DfxXeS8BUxdqEYvK8qXmfaDWwitnMbREYIF3hEI++XWgqoychuCuJ/9ZMsYdhvGE15uISA4k4L96U6O+upxMUlhJdAF/ADYBhQDyy/ieRbgR1LKRHovPz0shEgCngB+J6VMBn5j3/5MsZQy2f7z4MU34+qz8b1ttLZ1AuDh6cFvl6xQ4xIU5SJ4xsYwpdyXtBGC+NR2Rk7wdewrzMhDq/EnxJDF0x3X4m41EdPUQlbB69QUV9La0orJZHJe8FeBi0kKv5FS2qSUFinla1LKp4CffdmTpJTVUso0++N2IA8Ip/e26mfnfD5A1dcL/eqVtesAO0/1nSXEzJrGojWrnRiRolxdZlybhDUgAEu7kYWTAxzlR7LSyDElcqf7+6Rbh3NAP4bwcjMeQ/ay7cBDbNpyAx+sX8Lura86L/gr3MUkhQXnKFvyVV5ECBENjAOOAo8CfxVClAN/A35+VtUY+yI+e4UQM7/Ka1xN3vrnk5yurgVAp9PxyILpGNwMTo5KUa4eoWtuIb54NMfjBbN1EejdeucIq26po+hMN6PNggh5nGc6r8XDbCKurY7IiHyGBNQT5F1Nl/wrR45sc3IrrkwXWqP5ISFENjBCCJF11k8pkHWxL2C/H/EB8KiUsg14CPiBlDKS3ktSL9urVgNR9kV8fgi8LYTwPsfxHrDfizhRX19/sWFcMWqKqtiVk+PYjp48ljV3fN+JESnK1UcIwTWYyBihYejhBmbP6VtKPutEGqcZQVLIRk5Y4zmmTSSmyMjsVBtzj5Qx/1glJXVTyKx/huzsVCe24sp0oTOFt4EVwAb7v5/9TJBS3n4xBxdC6OlNCG/ZB8AB3AV89vh9YBKAlNJo7+WElDKV3nsX8Z8/ppTyBSllipQyJSgo6GLCuKK88bs/k1ra1w31vnkj8Q3wu8AzFEU5l/E//jHCNwhTj5nFU/uW6kzNzyK3I4I51lZSWvbzh67b2OWZwjv6CTzFao7b4rmnbBMvar7DE6d3sef4CXp6epzYkiuL7nw7pJStQCtwy9c5sOjtW/kykCelfPKsXVXAbGAPMA8otNcPApqklFYhRCy9y32WfJ3XvlIZu0zsSj+CtA9WixgVy513qwXsFOXr0A0dweQcN47FC8Y2+BMxzIuK4nZMFjN5mflMnhGHIfBT0i2/4dstjyCsEm9LOwXd/vwz8FX+lP5P7pj2F3Z0CNwPpBFkM7EsMIDfjB/t7KY5VX8OnZ0O3AHMO6ub6VLgfuDvQohM4HHgAXv9WUCWvXwd8KCUckCtlPHRX15lT2GaY/v2uQmExkY4MSJFubotGpNI9ggtIXvruHaRu6M8NTWNQlMiwSFt3HNyPf6WHn65KIG/JI0koSOA9R0TWWA5yi/T3uBW1jFFcwB0TTzXaiW3cUB97Xxl5z1T+KaklAeA843E+sIUoFLKD+i91DQgSSnZsm0dPabemcj9IwK5+5a/OjkqRbm6hd7/WwxPH6WHRqZETOAFw3YsRisV9VWcqupkgqcv3T6pvLE5ndCEH+N/+61opaBgfzOV7gUsb9tM6pH/MH9FMOmW7dzaNJRnTx7m2dnLnN00p1GT7FwmuZ9msbnggGP7hrnxjJia5MSIFOXqp/HwZm6BkWPxgqgDVmbN8nTsyzmRhaEnke1ToSDKQN1f/kLJylWMi2jE4LOEY62TiNHU0NX5DH9/8zTD909jjMzmU4snFtvgHQGtksJl8sqf/0pjazcAbt7u3Dj7l06OSFEGhrnXPcjJBB1BxWeYMzHOUZ6ek0lhVwSjhZU/3GjiiRtCqGyrpuGRB4n1LafSeDcnexK4U7eDW11+yhPVx0murKNF48N7edlObJFzqaRwGXTUd/JR3ieO7WVzhjHjtoVOjEhRBg7vpXfjo9VRGSBJakomLMoNgB6zkZzsXEa2RbOqvQN9WCV/usPM4QRByHt/wd/fhePyb5wZ80vCRRP/1P+K6FMFeFpbeKu61Mmtch6VFC6DF37wJ4qresdUaHVaVifciotrv93OUZRBRWi1LGzyZfNEDVHHSlgyr6+retqxVCp7pjHNMp1fB0/m791uRE5soSrARtjBf9LVaiarahHdtxzmjFcK39Z+wMzKI2RoIqhs73Biq5xHJYV+ZrXa+CD1fcf21BlDWfDQt50YkaIMPLO+/yQ1cXqwVTE+ZC56Q+9XW0V9FcUVpykpjsAU/T0SHjlCdcLdnF7QhVdLIREtu6kqbGbdP4tI73gIjZAsOJ2NVeh4NvPAl7zqwKSSQj/76E9vc7igb7DabYkLCIpWg9UU5VLSDUtmboGN7eMEMenuzJ3dNxlC+vE0zN6+fPjEH3jx+/cRrp9ItE882xdaic9Yx5Qz/yIlqo5uUxBZ5uHMsBwnxlzIth6B/GxQ0SCikkI/e/WjfzkGqyWMGsLUWbc6NyBFGaDWrHiMw+M0BNemMX1U32JV2bknOW1qQMbPoFsTQto7r3HYdxVrfE28uUhQ2ViA9+u/Y8S235LTNZMoTS1LC49RpQ3h09IiJ7bIOVRS6EfHPj7E9qy+wWr3TJ3OmBsH7Dx/iuJUAQtuZXSXjhPD2oiun0VCoisAVpuVmgP5uKOnxd+VrrhxeB7J4tXER3k0sJoP17jz0Hf1fDqhGWvNcDqkKzNqc9FJM6+XZjq5VZefSgr96J9//Slmc29/5+ihvoz3n4vQqJXVFKW/rNZPYvNEwdDSAuZNG+Yo33pyB8nb3mD5yVaCbV50hoRhO1TG8xH38Y+mM/y3tZK8STq8mtM4bBzHJHGciZ0nOCyCMVqsTmzR5aeSQj8pzy1jQ/oRx/at10xmwnfUmgmK0p8mfe9ZvP00NLtmkui2HB+f3q+4pk4Tm6dXUha6jXE5eXhKVzzdfag1JzJ54jtstF3LY9W1pA1No6xjDh7CyIr8dNqFN29m73Nyqy4vlRT6yR8evZfOrt6/MMKD3JjqOgu/qDAnR6UoA5vQ65le7sbHk2BokZV5c0Id+/70rwZ+taOEl00bSalyRWh0hJQW8n16eHbKd/iR+Dmm+HZ8K0yU2UKY3J6Bl2zjo4ZaJ7bo8lNJoR+0N7aw7sRBx/byheNIWjXLiREpyuBx9/fWczpG4tlxkJShM9Hr+y7ZFhSYeGlzLbd8+BgTW4diczFQu2cnjxUdoyopkQSjP9UeRzjUPZWRmgJWVe4kQxtLdWujE1t0eamk0A9+fe+dNDf3rgMb5KVlsfc8YhfMcHJUijI4uIVEcE2Dni3JrURUjeXPfxrOzOl+6LV9A0armnr46yd/YpY5CZu7NxWnz7C4YAd7u5eTM6qY5vpZNEtP7izZgVm48GLWDie26PJSSeESs5otvH98l2N76TUJRI4f6cSIFGXwWT3mB6SNlITXpeNhnMBvfufPO++Hc/fdfWOEdublk33yIIvMyWgNnnh1u3IgMIVkvcC9K413e5YwWuayouETdhj1TmzN5aWSwiX2+3vvpKq6EwBvAyyOXMS4e290clSKMrgkLb2b+d2CE9E52ArmkP7eTCqzJzFj2ixmzuxLDL/f/QzZZa8zrSUenUZHnGcNZ1qu4djwo1S0LqNK+vPjgjcp1MTyfvpmJ7bo8lFJ4RKy2WysPbTdsX3tjEiCg6LQaNTbrCiX2+TW0ewfZyOqJh93cQey8haaj09k9sgV+Pr2/uXf2GFjXc4OjjU9SUi1P4mNVXyomU9YSAsjGk/yUs9KRlhOc0PdZv6v2YWjuUe+5FWvfurb6hL61/fupbCkd9UmFy0sGT2fGT+8z8lRKcrgNPHHLzJDY6ba7RMSc14ipCUfLzmaQNMcls2e76j3SVYHlRWFNLd/SnC7Dn2QQN8ynryQfVR3LaDIFsb/FbxBl82Dn1cUk5t/0omt6n8qKVwiNquRV/f2nSWsmuCHv0ciLh4eToxKUQYvL08vpp/xZf1sycnQUuo7c9E3rCW0qZLx3kmMHZnoqPv0nmY+yfkE12YvEg1nWNe9CjG8lAlNVfzDfBMBsobtaT+hSBvHnwu2U19f78SW9S+VFC6Rv9x3HVk5NQBoBCybNpMZP1ZnCYriTF6Lf80cbRcvLuvg36uL+NtcPdsT6vGRSaxJSiBueN8a6Vuz23hp3ZPEnK6m3tUPf9NoWgz7aOiawq/NdxPfdZKNmY+w330a7xxa58RW9S+VFC6BzsYCXtl/1LF9U7IbQYbJePqq2VAVxZnGzF3OmjIDv6trJrFbj9Yzn2MJu+gx78HLfRr3pKQweUqMo352bSmHtx8mKLCVt5pvpH5EKtPau9lnns9Pzfczui2PN04+xgFhoqmpyYkt6z8qKVwC/+++eyks7h3c4qKFGxZMJuU7ajZURXE2jVZDZexy5jdoWdeQyQeFJjxa43h9dhahLS0YhD9rokeydPVwx3PSTp5gtO0UTdKbaO8kDKaDXNfmTmvQKn5qfoAZrRmMbSxj4+Z3nNiy/qOSwjdUfvQF3jqe49i+d6IbrpaZBEVFOy8oRVEcxtz1B4rcQmm0+BCrP8P7zQfpdhV8Ep+JPxPR4Mtcj+H4BRoA6DR2YjtYgMYbXm24ntL4nQQ3ZjMm3wSG5XxqHcd3z7xDruYMjfV1Tm7dpaeSwjdgs5r442P/obyyFQAPPSyZt5CE2691cmSKonzGx8ud9mtfxqjV0uXqj59JcFtLD0dHltGqPYKvZiYaSwCJo2Idz8nOLCE2oJraTm8mx4/nYOSLJOW+Snx1N2/03IG7rYdJrXl8+OZrTmxZ/+i3pCCEiBRC7BZC5AkhcoQQj9jLk4UQR4QQGUKIE0KISWc95+dCiCIhRIEQYlF/xXapHHn5et7PzHNsf2uGHy7GaUQnT3BiVIqifN7s8SNZG/04BmMLXX6RPNJWjbvZhdemV9BNKn5iOhO9hzrq5xalMa41F6mB1xpXkDjahyOxJ5h08LcMM0bwpnU+q+r20aUvoKqs1Iktu/T680zBAvxISpkITAEeFkIkAU8Av5NSJgO/sW9j33czMBJYDDwnhND2Y3zfSHdLHn9+oZym5m4Agj0Es2fdQPzNcxBCrZmgKFeae2+8nj9r7ye44xSdJj9+2VSH2a2W9HGBdIl0El0nExzhBYDRaqJidx2eISYyq124e80WmuZ4kzm0lck569nQc0PvYjzGDPZsHVgjnfstKUgpq6WUafbH7UAeEA5I4LMFVH2AKvvjlcA7UkqjlLIUKAImcYXa9Meb2XGy717C/Ysi8e0aRcyEiU6MSlGU8/HzcGHydd9np3UcHvoelra3E9PtwtGArXSPmgjaVsZM6DtbyMnJY6xXIVjgqb1l/OHm7exZ4YlfwwGmd2h42rKK8e2FWBpSsVkHzkI8l+WeghAiGhgHHAUeBf4qhCgH/gb83F4tHCg/62kV9rLPH+sB+2WnE84aQFJeuJYnt1joMVoASArSkDT+HoaunqjOEhTlCrZ4VCgfB30bPT10imheqi9FZ9PxZuAbdEZFkuwT76ibW5FBQkEx0k3Du9kVuLr68Jcb1/HJRMmUjDc5qp1DlzQw1K2YjIP7ndiqS6vfk4IQwhP4AHhUStkGPAT8QEoZCfwAePmzqud4uvxCgZQvSClTpJQpQUFB/RX2edlsJt7/45Mczeu7l3DDdeOJbAsjdsrkyx6PoihfzaI5s3nHMhdPcQbfTi3/rKrBqm/l/YhdxHUNIXKEPwBWaSP3kImw0AZaWgWHMqoI9h1K1MLxuPQUsKC6g83WyaR05pOx74CTW3Xp9GtSEELo6U0Ib0kp19uL7wI+e/w+fZeIKoDIs54eQd+lpSvG4Z0P8+zeFqTszVfzYvWEDruX6Bsnq7MERbkKLB45hHc8bsOkcYGQccwwtvFAvZEOr9P0uHQzflLfKOec3F3M1x5Fumq576NM2nvM3Lz0SQ6mwPQTr7HZMBMPacSXTDpaW53YqkunP3sfCXrPAvKklE+etasKmG1/PA8otD/eANwshDAIIWKAOOBYf8X3dTS1F/L+UyWUlBUBvac2k29ZzvguXyInjHVucIqiXBSdVsOK6ck8a1qOS+tRjLH38p32aia2azkek87okOjeuWqAnKYzuBxwIXlUPp09Nr7132PoPIKYsmAKVlmPR487xbZQEjQF7P9oYNxw7s8zhenAHcA8e/fTDCHEUuB+4O9CiEzgceABACllDvAekAtsAx6WUl4xd2+klHz6zr28cazvstEt4zzxCbyB+LvnOC8wRVG+spsmRrFWu4IWXRDu3Qfp8r+R77dXkh1rJqQigGFjQ4De69enjuezrGwfbsPg2JkW3jx6mimL/0D+OCu3b3mXd+VckkynaShOc26jLpH+7H10QEoppJRjpJTJ9p8t9vIJUsqxUsrJUsrUs57zRynlMCnlCCnl1v6K7evYcfS3vP5GAE311UDvQLXIm29mqc0Hn6GhX/JsRVGuJD5uelakDOfRnvuQ9QV4RgqCW+MJlFaErYspk/vuVx4t248x/ya+HfAGNj8Xfrshh2oCmDl3Eu6mBk54jcQstQxxOUlVYYkTW3VpqBHNF6GmrZhjH+Ww4/gnjrLvTPfFzX0J8XfMdGJkiqJ8XfdMj2avdTQHI+5DZL+L/8jJLOjs4mRUAXGxw9Doe9d0LuzsQFO4C8OeOGYlHcNqg79szyVm/o8ojLexNHc3O23jSbFmk71l15e86pVPJYWL8MH7P+WVjZUYe3oHqsX6CVyuvZPbAyMxBHk5OTpFUb6OoQEeLB8TxgOn59IdNRtDyRtMatFyPN6MzylPRk4OdtRNLfkYbfs45hdnQYCezSdrsASPIWCUgWlHc1hrmIef7MCrQSWFAS89dyPpH+soye2bGvv/LQ3G6DaPmOvGOzEyRVG+qZ8vSUCi5Zfi+wjPQIJaPPDVWNG3GJkxMcBRb297PcNKP6QjfzXx4WWYTbA5r5bpc9ZQGSgQbmaO2+IZoTtIU/7VfQlJJYULsNlsvLllE+8f6sv+t4zSkTXlHh4ZPw6NyxU7C4eiKBchzNeN718Tx/oCI7lJPyBaVLOgs4viIUV4x41D5+4KQEmrjWbXw/g1t7DImo40aHhyZy4uY26mLcHIfYc+5Gl5HQGyhdYNTzm5Vd+MSgoX8PzLT3DqoyLaGhsA8DJo+N6iaMxe8wmfGPklz1YU5WrwrRkxxAZ58MPsCDz8dExs03BglAn3Sg3jZ/TdcD7YUk6seTfemYF4hJs5XWukXBvBpNEheBT2UBQaRqYtloDOjUir2Ykt+mZUUjiP5uo6Thak88mJvpGKf5ij5+2Ye/jdsqlqoJqiDBAuOg2/u3Yk+U2SAr85GFp98NbZ8C1u4ZoJ3o56b2RZiY/ag6EyjFmBJ5DAv/YXETn+Joribaxq38VT1tV4iwY6Nz/vtPZ8UyopnMdP33+d47vzMBlNAMQF6pk8cSRLJt2IR6C7k6NTFOVSmhkXxOKRQ3iyNpkwXTMLOjvJDzuNLnwS7sE+ADR1W3mtoJNAayqTyqshUM+G1ArMSasJie9g/rZj7AlMIU9Gocn4D8gvzNJzVVBJ4Ry2bdyBS/p6UtOzHWVPL9LzXuhdLJ04/ALPVBTlanXL5Cg+NSUR5q1lUrvgcKKJGlMIq28KcdR5fL+FsX7b6D45kuERpzGZJJvKdUwfPppcH8GYIYW8aZmPu60cWZXhvMZ8AyopfI7JaOb54hy2HD7jSPTz4j0IHjaCH995v3ODUxSl30yNDcDd1cARj2uwtIfgL2wMKz3Fsom+uPj2dj1v7DSyta6TwNYW5pKBdNXy7/1F6MbchHV0J9fl7GCD5yxMUodp78tf8opXJpUUPufRF98i5PjLlBVUAqDVCp5bADv9VxHm6ebk6BRF6S8uOg0LEkN4rikFV30PCzs7KfAopbpxOKvX+Dvq/XG/kRjPnYTmGdCFazlV0cHp8BUs8gR5qhhLuAs7bOOh8GO4Cm84q6Rwlrz8UrqaCli3p9JRdl1KEB4BIdx+9w+dGJmiKJfD4lFDONoTzhhfPVM74ES8icPWZG5cqkfv0Xsv8XRjN8UtaViLhjAmLBcp4PW0eoYk305ZkoUJ7jmss83CINugcIeTW/TVqaRgJ6XkwYOZmLLW01TVDICruwvPzepiq8dyQn08nRyhoij9bVZ8EO4uOg57zKOnM4xAm40gWUVH3XBWrPZ11HviqJUh7qcZ31GMLciVd0+UYxp/H3PdO0iozGaP3yQapTfWtLec15ivSSUFu6de38qkk2/z4Sd9i7/dMicaFzcPpt/wYydGpijK5eKq1zI3IZjn6seg19hY1NlJj2kfh4snccsqAxq9HoCMKhP+tg1ElUj0ERo6ui1srzIwLXIu5uaTyDADH1mnQ+F26Gpycqu+GpUUgJbTTazt7iQrP42e9t75jfyDPHlqXC0btdeQFKMGqinKYLFk1BCyOv1ICPBgRoeN42Emjg1NwNIRTcqUvnEL2wrrMBX7Mt4vE1w1vHn0NNopDxLrXk+CfxHr5Uy00gwnP3Bia766QZ8UbCYrj2w/yvy0p9m5s8xRfuvSJAxaDYHzf+C84BRFuezmjgjGoNOQ5TGVrq4ogqxW4nxKKE+PZel8naPe2lyJn76eUW0lmCM9OVrSRIZ2NNd5RjCivYCsoASybbHIQ8+A1eLEFn01gz4pHPswjfbWXDZmncFm6V3TJ3ZYMH8cWsY2OZl5k9Wkd4oymHgYdMyKD+LZujHYhJ5FnZ24y8N8IiYxPGY4encDACWNZny6NhJTakEfKXAxCP76SQFeE+4ltjYba5gHT1tWIlrKIGf9hV/0CjKok4KpqoPfdNQRdvR9so6edpRfv2Ys3qKH5uQH0GsH9VukKIPSklFDyGr3JDrQi/ntZopEGUVJMVhOBjNxVqCj3paCBmynAxmvTcMS7cbBokaOes5lQccphvjXs0s7njIZidz/d7DZnNiiizeov/E2fXSUEXXv815GlaNs0oRhPOJezBFbIksWLHVidIqiOMs1iSHotYLTXhNo64kltNNCTHA5eZXRLJ3VNzvyu7lW9JYGkltz6Yr0xcfThcf31JEYO4OR1hyMw3150rwSUZ8P+Zuc2KKLN2iTgqmmgxcNddTtPEV9cW+PI71Ww5KbRhFua+BE5J0EeRmcHKWiKM7g46Zn2rBAnqsbRZv05K62VkK1WRwKG8PQ8OG4+PSOcK5ut+LSuZOYM5IwTQX64a5klreQ5r+UxKZTmCO8OCamUU4otn1/vSrmQxq0SWHbrnQm1jzP5pN9C2Isnj2Om43F5NsiWLLyDidGpyiKsy0ZNYSMFjcignwZ2eGKrieH4oRhmCqCSJnbNx/SzpwybKfDWWDdTkWQJ2H+bvw+N5QZLaUYhJHgofAv87VoarJg/QNg7nZiq77coEwK0ibZ0XSUwk80dNbVAuBp0DHtxmEkWM6wJ/gWhgWrZTYVZTBbkBSCRkCz72j2WccxozYf9yALbcVerJxtddRbn2dBa2xmRE0lHqIT96EeZFZ1EBE4gdFkUBDtRallNk9zM2S/B/9dDK0VTmzZhQ3KpFCSXkZEw1o25uU7ypbOnMycxnKqpD+zVz3oxOgURbkSBHgamBTjzzP1YxiqbWJ+o5U41zzKZSBDw8Jx8fcDoLFb4mI5hEt5EnPZQa5X79fqMfeFjO7KpVXvzVIvL542XssTfr9FNhbDu7dfsZeS+i0pCCEihRC7hRB5QogcIcQj9vJ3hRAZ9p8yIUSGvTxaCNF91r5+W6Viw75/k/epgc76GgA8XXQk3hjNFFMOn3hdS2JE4JccQVGUwWDJqFAyGjQkx0VxxDyZGGMWaRHxNDdGkTSl7xJSWm4GLjWeTOs4Aq4aAoPcWVvmQUprBVppISOsm4eljeeqR7B5yINQlQ7lRy/wys7Tn2cKFuBHUspEYArwsBAiSUp5k5QyWUqZDHwAnN2Bt/izfVLKfvlz3WKxoOnYzqbTfT2Olk6dTlxbD0apI3n5d/vjZRVFuQotGjkEgJ1uS3CX3cxozKJw+HDMpX7Mm9LXxXTTKQshYRJdVQgTOE6rv5bM8hYio1aSxEn2hBq4Dl/uHBXKT04lYtJ7w9Erc3W2fksKUspqKWWa/XE7kAeEf7Zf9K5neSOwtr9iOJe1a9fSsF9DY3nvuAQXrWD4bVEsbz/IXt0UkhPiLmc4iqJcwYb4uDI+ypeXysNZFNSIocoDXSCY6yUpcd4IvQsAOfU2gvVHcK9IZqrtAB1BvTOqZlqnMc6SRr2rLzmu3dx2qonR0aG8ZZqNzN0ArZUXenmnuCz3FIQQ0cA44OzzpZlArZSy8KyyGCFEuhBirxBi5nmO9YAQ4oQQ4kR9ff1XjqW2bj2f1Lo6tqcmjsbd6o0PXbhP//ZXPp6iKAPbklGh5Na005WwBnOnK7GGUmo1nnR2DCVobN9KjDv37MPb6kF4fQN6dwu+/q5sy6tnvLULIW287VOLh8nAL4P9ed26ACklHH/JiS07t35PCkIIT3ovEz0qpWw7a9ct/O9ZQjUQJaUcB/wQeFsI4c3nSClfkFKmSClTgoKCvlIsJVXN2IoKOZGf7iibcG04y2oPU0QEU2cv/0rHUxRl4FsxNgyNgHfNM5irO8nwnlOkh8XT1BTOhOm+jnqb8joYkdSJtno040ilw1+QerqZyICpxFFAdmwIpe3Z+B+vZ8X4Ceywjsd8/JUrrotqvyYFIYSe3oTwlpRy/VnlOuA64N3PyqSURillo/1xKlAMxF/KeF7/71PkpntiMxkBGBoYgGdCBMmyiLLhN6NVU1ooivI5Q3xcmREXxBvZPQQlTCOupZLchHg0NQYWTujrQbTntA2/ipfxqx3BeFMqnUN6u7UX9kxkgjxGtXcwO225GG1d3FhjYZvHKvTGZswZ7zmraefUn72PBPAykCelfPJzu+cD+VLKirPqBwkhtPbHsUAcUMIl5O6azt76vrnNZ88ZxsiKGrqlCxOXPnQpX0pRlAFkzYQIqlp7yA1dTVJdCZYAdzStnYToPDAMjQXAZJHszCgl3E9LaHUHLu5WvL00fJjdzmTLGXTSzNGx15DfdATOtLNq3irybJG07nnmiuqe2p9/Gk8H7gDmndXN9LPJhG7mizeYZwFZQohMYB3woJTykq1OYbXYMJWc4czpvjwzaoGBZT1H2O85Fx9/1Q1VUZRzW5gUgperjperhjLd3ECwSx2Nwp3mliEMnd53X2FDmSsJHh/gXjmOZFLpDNGTWd5CkPsY5rKT7LgoMmUlEitj68wcD7mBwM5TNOfudmLr/ld/9j46IKUUUsoxZ3Uz3WLfd7eU8vnP1f9ASjlSSjlWSjleSrnxUsZTmV9N9akgbNbekYgxoUF4GcPxEEa8ZqkbzIqinJ+rXsu1Y8PYktuA26iVxJmKyAqNpbUplOnT3Rz1PsjuwFazjyHGcEZ15NIZ4YdWQHr7LJbJj0HAwXHzqegsojO1mtkrH6JZelKx7R9ObN3/GjQX0bfv2M+ehgbH9sQpfsxqPMlJEcuUSfOcGJmiKFeDG1Ii6THb2GW4hpEdhaQnjsa1TjIjsgdtSO/qjB3dZtYXu5Lsf4rQChOuLiZ8PI28n2El0m0Y0217yRwxjAxLIbLbRlCzIC/sOpLa9lOQn+PkFvYaNEnhdFMhxQV9b3rCRF+SOEN2xHJ6b38oiqKc39gIH4YHe/LyKTdm9ZRjCXBDNLWjbdHit3i+o95L+Qb8W17Cv2Yss2y7qYsNprbNiNH9ZlZp1mHTCDYnJdNlaaP9cCWjV/0QKQRl255yYuv6DJ6kUJ2Osae311GAjwdhrkMBSJi+xplhKYpylRBCsGZCBCdON5MQMgm93kq11pvW5lAmLhkKovfrdG9mBWdau4h382Bs5SlsQa7ohYV1x7yJ0Lsw1XyItKQR7NWUYippw80QRp7PLKY0b6Sns+1Louh/gyIpdLYYaak949geNdqTkW2VFItwxiWMdWJkiqJcTVaPC0cI2GieR4ytlKyoaLoahjAp8DSGMZMc9V45E8ZIuRFRFssUDmOK8GRHcSshQ27lTv2LGKxGnk9JxIagPasac8q38RGdnN7zqvMaZzcokkJlQRVd9SbHtm+oC5NseeQEzHBiVIqiXG1CvF2ZPiyQtSe7GG0uJTcuCZ+6Tka0leGxvG/w64sHq3Dp2cUIy1AmVOdiivTGIjS8+oYZb2Hhuo4PqQwK4I1IG02pJSRNXkiujMY3+79O7546KJJCTWYGjR19k1e56MNxEVYCxl/rxKgURbkarRoXzpmmLhL1Lpj8PTE2t9HV7MvkGd4IN08Aquva2VfvxXivBrSlgYzyyEFEGHi3JxJt4wgWeW9gaGMNL8Z7Ut8qMei1HA2+gZCeUmTpXqe2b1Akhaqq0zR1tDu2k1ystOLBlIkLnBiVoihXo0UjQ3DVa6jsGIN001KmD6KhbijTXY7hNm+xo96/81wJ7HiJ4aZoJtdk0hUXgFYnWJsah1Zj5trmTZg1gneHetFZ1oBXys00SG/a9z7jxNYNkqSQRxetLc2O7TkeNWS4j0OrV2swK4ry1Xi56lmQNIStWTb8aCY/JhqXJj8S28vwWrbUUW/DkVJ6zFUk6jxwLXYjRXcc4whvtssZmFsNJPhnMbyuko8i9Jw6kMuckZG8a52D5+lPoaPOae0bFEmh22amva3DsT3FtwVj2HgnRqQoytVs9bgwWrrMxMkaymIj8as+Q2PdEKbEd6CLiAHAaLTyTsNwosRmIszhzKs6gIww4OKqY1/NNHx8q5hQm023TrDOaCXQ00Bh8CI02CBvg9PaNiiSgrm9rneaWsDDyw1XnSB+1CwnR6UoytVqZlwQ/h4uGHp0tPoH09JqoqEhmqniIK5LVzrqPXuoHA/bdkZbIug5E8pSNtGe6Mun9dMRQjLSPZPYhgY+DvehvaOTuFGTKbKFYcpaf4FX71+DIik0d5Y7Hgf6uGNFEJ041YkRKYpyNdNrNawYE0phsQ7ppSPTZziGVi/iWmrwW3ANaHq/WjNz6ynSD2OI7hSTO8YxuTYDn8AOGmyhlDdFEhJczMiKUzS6ann50EnmJASz2TYFXflhaK91StsGRVJoqe+b3mKol+C0NgoMnk6MSFGUq92qceE017oikOQmxBBxpoz6+lAm+xfiOn6ao94zJ9oIEu8yQkagzY1nHjvpjvNmb+0U3L2bibfmEdlhYlO3haRQbw67znTqJaQBnxSkTdLeoHNsj/LpoTlgjBMjUhRlIEiO9CXG140htnrMIQaq2o00NAxlmu0AhhWrHfXeOlCE1teI3rWYGebFjGwuwjbEjbSaiRjNeiIjchleX0uepzsdFitRCRMoJhxbzkdOadeATwo9nWYaO/tWNkrwtRASP8WJESmKMhAIIVg1LhxNawfNYaHsCU7Bq92TqOZmgiaPQePpC0BTo4mPG6Pxtb2AQbgRludNpDiDJdKdnRUzCQwqY0hXPlaNYHtpOXNGBLPRMhlx+qBTLiEN+KTQXt9Oe0ffsgxDfQVhidOdGJGiKAPFquRwuuv1dLt50RTlTXhlOeVVw5mmP4rbNUsc9f6x4xCGYYG4aE4yqmcNE61HaIwdwsGyeUiblmS/3bhZbGwvrWH68EB2yEkIJBTtuOxtGvBJobWskvaWFsd2qI8ezZBRzgtIUZQBIzrQgyhdAAChkY0UdJppbQplYns2rqtvcNQ7lNpEyehH8Q48gqvVm/FVp0ErMHi6Ul0UTVhIISNamjgqtXi76vCIGEuj8IeinZe9TQM+KdSVnqGltW/mQVefANDqnRiRoigDyc1JYRhs3cgILQdCJhFkdEN32ovwSFfcR/T+AWqzwd///iiGbz2OXlfC5KoJRMkyuhJd2FecjEAyQhyhzuBCYWcPcxKD+dQ8GlvRbrBaLmt7BnxSKCk+hdnU+6a66ARmdz8nR6QoykBy7dhwPLvaKNPF4hZiJqq8kpqmMCYZj6NffZOj3gef5mNx98N78Uh8W1NIMR+nwiucLkbTVhvNeI9dAGw9U8nM4UHstY1FY2yBytTL2p4BnxQK2/rGKPi5aWnXeDsxGkVRBhp/DxeCbZ6UE0VyaA6HLSYirAFElTdhmLMAnbsHALW1Zl7/+3dxnTYBraeFuVVdAFgSNPSUhBJhOE14dxc7q5pICvMmQ5+MDc1lv4Q04JNCg7HvJnOAG3TqVFJQFOXSuikqCim0tEW4cSokhZhmV6yV3kRrqwhcMMdR76X1HyM0Ap9FIxlnv4RUGuZDR5E7UgpGGnNJFzrqzRbioyPJ1cSrpHCp9Zj6uqMGuUmkm7p8pCjKpfWt0WF4WtrJchuFt0835voi4mxhJDYVY151r6Pe0RNNbHr1H7iPC8bFEsUUYwblbhF02SRtzWHM063DhuDZ07VMjg1gu3EUsiodOhsu8OqX1oBPCl3GTsdjfzeBd0iME6NRFGUg0uu0JOvgJGOIHZLHLl0zo4zhDCnpRh8VTfi4RKB3/ZyfP/43bEh8Z0ezqLp3YG35iFA6KoYR6VHIvNp6XqusJy7cm322Mb1dU8v2X7a29FtSEEJECiF2CyHyhBA5QohH7OXvCiEy7D9lQoiMs57zcyFEkRCiQAix6FLE0dnd6njs7yYYEpV4KQ6rKIryPx5KHIVV6GkO96Q8aDIdHc2M7ggi1liK5/cfcazhfLKwiue+/XM8J4WSWJ5ElCzl2KjR6HJ6r2Is79mHWcI+m5Ez+pje+wp1+ZetHf15pmABfiSlTASmAA8LIZKklDdJKZOllMnAB8B6ACFEEnAzMBJYDDwnhNB+0yC62v73TGFIqDpTUBTl0psb4ou37CDXK55gD6htL2OsNZroyjpahk5m5rxoR90nNrzFmQPHcY+IYVp3PoU+w/GoaKCj3R8/vz1MO13L61UNxEcGU60JgfoBkBSklNVSyjT743YgDwj/bL8QQgA3AmvtRSuBd6SURillKVAETOIbsJptdPf0NdHfTeAZEPZNDqkoinJOGiGY7WEhW4whyD+H47oGDGZIKe9EKy2EPbAGrYsrABUNVbz95HsYRviwuNIFgHb3VuorR6LzqeYW/Q66JegCDOSaw7BexjMFIS/DItFCiGhgHzBKStlmL5sFPCmlTLFvPwMckVK+ad9+GdgqpVz3uWM9ADxg3xwBFPRz+IHA5bvLc2UYjG3uD+p9/OYG63vY3+0eKqUMOtcO3bkKLyUhhCe9l4ke/Swh2N1C31kCgDjH07+QsaSULwAvXNIgL0AIceKzxDVYDMY29wf1Pn5zg/U9dGa7+zUpCCH09CaEt6SU688q1wHXARPOql4BRJ61HQFU9Wd8iqIoyv/qz95HAngZyJNSPvm53fOBfCllxVllG4CbhRAGIUQMEAcc66/4FEVRlC/qzzOF6cAdQPZZ3U5/IaXcQm8vo7MvHSGlzBFCvAfk0ttz6WEppbUf47tYl+1S1RVkMLa5P6j38ZsbrO+h09p9WW40K4qiKFeHAT+iWVEURbl4KikoiqIoDiopfI4Q4gb7tBw2IcR5u4QJIRbbp+MoEkI8djljvNSEEP5CiB1CiEL7v+ecNdA+LUm2fYqSE5c7zivRl30ORK+n7PuzhBDjnRHnle4i3sc5QojWs6bI+Y0z4ryUhBD/FULUCSFOnme/Uz47Kil80Ul6u8vuO18F+/QbzwJLgCTgFvs0HVerx4BPpZRxwKf27fOZa5+mZND1Hf+8i/wcLKG3J10cvYMu/31Zg7wKfIXfp/2fTZEjpfz9ZQ2yf7xK75Q+5+OUz45KCp8jpcyTUn7ZKOlJQJGUskRKaQLeoXeajqvVSuA1++PXgFXOC+WqcjGfg5XA67LXEcBXCBF6uQO9wg2036eLIqXcBzRdoIpTPjsqKXw94UD5WdsVnDWv01UoREpZDb1zVgHB56kngU+EEKn26UYGu4v5HAy0z0p/uNj3aKoQIlMIsVUIMfLyhOZUTvns9Ps0F1ciIcROYMg5dv1SSvnxxRziHGVXdN/eC7X5KxxmupSySggRDOwQQuTb/9oZrC7mc3DVfVac4GLeozR65+vpEEIsBT6i97LKQOaUz86gTApSyvnf8BBX3ZQcF2qzEKJWCBEqpay2n57WnecYVfZ/64QQH9J72j+Yk8LFfA6uus+KE3zpe3T2vGlSyi1CiOeEEIFSyoE8WZ5TPjvq8tHXcxyIE0LECCFc6B2hvcHJMX0TG4C77I/vAr5wtiSE8BBCeH32GFhI7035wexiPgcbgDvtPUmmAK2fXapTHL70fRRCDLFPnYMQYhK9312Nlz3Sy8spn51BeaZwIUKI1cDTQBCwWQiRIaVcJIQIA16SUi6VUlqEEN8FtgNa4L9Syhwnhv1N/Rl4TwjxLeAMcAPA2W0GQoAP7b+XOuBtKeU2J8V7RTjf50AI8aB9//PAFmApveuDdAH3OCveK9VFvo9rgIeEEBagG7hZXuXTMQgh1gJzgEAhRAXwW0APzv3sqGkuFEVRFAd1+UhRFEVxUElBURRFcVBJQVEURXFQSUFRFEVxUElBURRFcVBJQVHOIoTo+JL9vkKI75y1HSaEWGd/nGwfbftVX/P/hBA//urRKsqlp5KConw1voAjKUgpq6SUa+ybyfT2K1eUq5ZKCopyDkIITyHEp0KINPsaEp/N2vlnYJh9Tv+/CiGihRAn7SNxfw/cZN930+fPAOz1ou2Pf2lfP2AnMOKsOsOEENvskw7uF0IkXL5WK4oa0awo59MDrJZStgkhAoEjQogN9K41MUpKmQzw2Ze8lNJkX/glRUr5Xfu+/zvXgYUQE+idymEcvb+DaUCqffcLwINSykIhxGTgOWBev7RQUc5BJQVFOTcBPC6EmAXY6J2yOOQSHXsm8KGUsgvAnmwQQngC04D37dOJABgu0WsqykVRSUFRzu02eue/miClNAshygDXr3gMC/97ifbs559rfhkN0PLZWYiiOIO6p6Ao5+YD1NkTwlxgqL28HfA6z3M+v68MGA9gX183xl6+D1gthHCzzzy7AhzTQ5cKIT6bkFAIIcZeuiYpypdTSUFRzu0tIEUIcYLes4Z8ACllI3DQftP4r597zm4g6bMbzcAHgL8QIgN4CDhlP0Ya8C6QYa+z/6xj3AZ8SwiRCeQwCJalVK4sapZURVEUxUGdKSiKoigOKikoiqIoDiopKIqiKA4qKSiKoigOKikoiqIoDiopKIqiKA4qKSiKoigO/x/enQyUIPym1QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "for zonal_mean in zonal_means:\n", " ax.plot(0.5 * (lat_bins[1:] + lat_bins[:-1]), zonal_mean)\n", "ax.plot(0.5 * (lat_bins[1:] + lat_bins[:-1]), zonal_means.mean(\"time\"), lw=3, c=\"k\")\n", "ax.set_ylabel(\"tas / K\")\n", "ax.set_ylim(270, 305)\n", "\n", "# Scale the x-axis to account for differences in area with latitude.\n", "ax.set_xscale(\"function\", functions=(lambda d: np.sin(d), lambda d: np.arcsin(d)))\n", "ax.set_xlim(np.deg2rad(-80), np.deg2rad(80))\n", "ax.set_xlabel(\"latitude\")" ] }, { "cell_type": "markdown", "id": "6b3df69e-d608-4366-a278-e1fbb398ecf2", "metadata": {}, "source": [ ".. tip::\n", " The concept of using histograms to compute zonal means can also be generalized across other dimensions,\n", " i.e., a meridional mean or to compute distributions in temperature or humidity space." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (based on the module python3/2022.01)", "language": "python", "name": "python3_2022_01" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.9" } }, "nbformat": 4, "nbformat_minor": 5 }