Skip to content

Instantly share code, notes, and snippets.

@mauro3
Last active November 7, 2023 21:51
Show Gist options
  • Save mauro3/04163405f15196a58f38d758bc0009fd to your computer and use it in GitHub Desktop.
Save mauro3/04163405f15196a58f38d758bc0009fd to your computer and use it in GitHub Desktop.
Appl-GL Ng & Björnsson 2003
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Ng & Björnsson 2003\n",
"\n",
"This solves and plots the equations\n",
"\n",
"$\\frac{\\partial Q}{\\partial t} = {\\frac{4 \\overline{ \\nabla \\phi}}{3\\rho_iL}\n",
" \\left(\\frac{\\overline{ \\nabla \\phi}}{F_1}\\right)^{3/8} Q^{5/4}}\n",
" - \\frac{8A}{3n^n} Q N_{l}^n$\n",
" \n",
"$\\frac{\\partial N_{l}}{\\partial t} = \\rho_w g \\frac{Q}{S_{l}}$\n",
"\n",
"**Parameters to play with** are marked with a `<---` in the comment.\n",
"\n",
"As is, the parameters are for a Grimsvötn like flood.\n",
"\n",
"**Technical note:** this will take some time to start-up and for the first execution. To evaluate the cell, after changing parameters press shift-enter."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABpMklEQVR4nO3de1xUdf4/8NfcYbgMd4YRREq8gopQillaKlqZtbVlaWbltllmsWqa2Xe1fhuoPdZubqbWZllGe8myUpM2Rc3wQpIKZl5QUUAEcYbrDDDn98fAkRHQGRw8wLyej8c8YM75zJn3fDLn5ed8zufIBEEQQEREROTG5FIXQERERCQ1BiIiIiJyewxERERE5PYYiIiIiMjtMRARERGR22MgIiIiIrfHQERERERuTyl1AZ2F1WpFQUEBfHx8IJPJpC6HiIiIHCAIAsrLy2EwGCCXtz4OxEDkoIKCAkREREhdBhEREbVBfn4+wsPDW93PQOQgHx8fALYO9fX1lbgaIiIicoTJZEJERIT4Pd4aBiIHNZ4m8/X1ZSAiIiLqZK423YWTqomIiMjtMRARERGR22MgIiIiIrfHQERERJ1Kbb0Ve09egKmmVupSqAuRNBCtWLECAwYMECcqJyYmYtOmTeJ+QRCwaNEiGAwGeHp6YuTIkcjJybE7htlsxsyZMxEUFAQvLy9MmDABZ86csWtTVlaGKVOmQKfTQafTYcqUKbh48eL1+IhERORCVquAqf/cgwff/xlJy7ajpMIsdUnURUgaiMLDw7F48WLs27cP+/btwx133IF7771XDD1Lly7FsmXLsHz5cuzduxd6vR5jxoxBeXm5eIzk5GSsX78eaWlp2LlzJyoqKjB+/HjU19eLbSZNmoTs7Gxs3rwZmzdvRnZ2NqZMmXLdPy8REV2bDb8WYNfxUgBAkakGH/2UJ3FF1FXIBEEQpC6iqYCAALzxxht48sknYTAYkJycjHnz5gGwjQaFhoZiyZIlePrpp2E0GhEcHIy1a9di4sSJAC4toLhx40aMHTsWhw8fRr9+/ZCZmYkhQ4YAADIzM5GYmIjffvsNvXv3brEOs9kMs/nSvzwa1zEwGo287J6ISCIPrfwZe/IuoI/eB78VlUPv64FdL90BuZx3EKCWmUwm6HS6q35/d5g5RPX19UhLS0NlZSUSExORl5eHoqIiJCUliW00Gg1GjBiBXbt2AQCysrJQW1tr18ZgMCAmJkZs8/PPP0On04lhCACGDh0KnU4ntmlJamqqeIpNp9NxlWoiIomdM9Vg78kLAID3H42HVq1AkakGR4srJK6MugLJA9HBgwfh7e0NjUaD6dOnY/369ejXrx+KiooAAKGhoXbtQ0NDxX1FRUVQq9Xw9/e/YpuQkJBm7xsSEiK2acn8+fNhNBrFR35+/jV9TiIiujY/Hy+FIACx3XToEeSFuO5+ACCGJKJrIXkg6t27N7Kzs5GZmYlnnnkGU6dORW5urrj/8pUlBUG46mqTl7dpqf3VjqPRaMTJ3lydmohIervzbHOHht4QAABIiLT9/OVUmWQ1UdcheSBSq9Xo2bMnEhISkJqaioEDB+Ltt9+GXq8HgGajOMXFxeKokV6vh8ViQVlZ2RXbnDt3rtn7nj9/vtnoExERdVy782wjQUOiAgEAMd10AIDDReWtvobIUZIHossJggCz2YyoqCjo9Xqkp6eL+ywWCzIyMjBs2DAAQHx8PFQqlV2bwsJCHDp0SGyTmJgIo9GIPXv2iG12794No9EotiEioo6tpMKME+crIZMBN/WwjQz1CvUGABw/X4F6a4e6Pog6IUlv7vryyy/jzjvvREREBMrLy5GWloZt27Zh8+bNkMlkSE5ORkpKCqKjoxEdHY2UlBRotVpMmjQJAKDT6TBt2jTMnj0bgYGBCAgIwJw5cxAbG4vRo0cDAPr27Ytx48bhqaeewsqVKwEAf/7znzF+/PhWrzAjIqKOJafABACICvKCTqsCAET4a+GhkqOm1opTpZW4IdhbyhKpk5M0EJ07dw5TpkxBYWEhdDodBgwYgM2bN2PMmDEAgLlz56K6uhrPPvssysrKMGTIEGzZsgU+Pj7iMd58800olUo89NBDqK6uxqhRo7BmzRooFAqxzWeffYbnn39evBptwoQJWL58+fX9sERE1Ga5DYGoX9il+ZxyuQw9Q7xx6KwJx4orGIjomnS4dYg6KkfXMSAiIteb+fl+fPNrAeaO641nR/YUtz/7WRY2HizCX8f3w5PDoySskDqqTrcOERERUWtyC4wA7EeIACDcXwsAOFNWfd1roq6FgYiIiDq0KksdTpRUAgD6GS4PRJ4AgDNlVde9LupaGIiIiKhDO1JUDkEAgn00CPHxsNsXwREichEGIiIi6tB+a1hnqI/ep9k+jhCRqzAQERFRh3a84V5l0SHNA1G3hkBkqqmDsbr2utZFXQsDERERdWiN84duCPZqtk+rViLQSw2Ao0R0bRiIiIioQzt+3jZCdGMr6wyF+dnmFZ0z1Vy3mqjrYSAiIqIOy1xXj/wLtpGfG1sYIQIgTrQ+ZzJft7qo62EgIiKiDutUaRWsAuCtUSLYR9Nim5CG7cUMRHQNGIiIiKjDOiGeLvOCTCZrsU2Ir22EqLicp8yo7RiIiIiowzp+vnFCdev3KWscIeIpM7oWDERERNRhHW8yQtSaxkB0niNEdA0YiIiIqMM64cAIUah4yowjRNR2DERERNQhCYJw1UvuASDEt3GEyAyrVbgutVHXw0BEREQdUkmFBeU1dZDJgMhAbavtgrw1kMmAOquAC1WW61ghdSUMRERE1CE1jg5F+GvhoVK02k6lkCNAa1utmpfeU1sxEBERUYd0af5Q6xOqGwV62wLRhUqOEFHbMBAREVGH5Mj8oUYBDfczK63kCBG1DQMRERF1SI2LMjo0QuRlm1jNESJqKwYiIiLqkBoXZXRkhKjxlFlpBQMRtQ0DERERdTg1tfU4U2a7qasjI0SXTpkxEFHbMBAREVGH03hTVx8PJYK9W76pa1OBXo2TqjmHiNqGgYiIiDqcS/OHvFu9qWtTAZxDRNeIgYiIiDqcEyUN84eCrn66DOAcIrp2DERERNThHC9uuOQ+5OoTqoFLp8w4h4jaioGIiIg6nOMNI0Q3ODhC1Dip2lhdi9p6a7vVRV0XAxEREXUogiDYzSFyhJ9WjcapRmW8nxm1AQMRERF1KI7e1LUphVwm3s+M84ioLRiIiIioQ2kcHQr397ziTV0vF+DF+5lR2zEQERFRh9K4QvUNQY6dLmvExRnpWjAQERFRh+LMPcyaunTpPRdnJOcxEBERUYfSuAaRoxOqGzWOEJVV1bq8Jur6GIiIiKhDaRwhutHJEaLGSdW8fQe1hdKZxoIgICMjAzt27MDJkydRVVWF4OBgxMXFYfTo0YiIiGivOomIyA2Y6+qRX1YNwLG73Dfl3zhCVMkRInKeQyNE1dXVSElJQUREBO6880589913uHjxIhQKBY4dO4aFCxciKioKd911FzIzM9u7ZiIi6qJOl1ah3irAS61AiM/Vb+raFK8yo2vh0AhRr169MGTIELz//vsYO3YsVCpVszanTp3CunXrMHHiRLzyyit46qmnXF4sERF1beIVZg7e1LWpS3OIGIjIeQ4Fok2bNiEmJuaKbSIjIzF//nzMnj0bp06dcklxRETkXk6UtO0KMwDw13KEiNrOoVNmVwtDTanVakRHR7e5ICIicl8n2rgGEWA/QiQIgkvroq7P6avMNm/ejJ07d4rP//GPf2DQoEGYNGkSysrKXFocERG5l7auQQRcGiGqrRdQbq5zaV3U9TkdiF588UWYTCYAwMGDBzF79mzcddddOHHiBGbNmuXUsVJTU3HTTTfBx8cHISEhuO+++3DkyBG7NoIgYNGiRTAYDPD09MTIkSORk5Nj18ZsNmPmzJkICgqCl5cXJkyYgDNnzti1KSsrw5QpU6DT6aDT6TBlyhRcvHjR2Y9PRETtRBAEcQ0iZ68wAwBPtQKeDbf6KONpM3KS04EoLy8P/fr1AwD897//xfjx45GSkoL33nsPmzZtcupYGRkZmDFjBjIzM5Geno66ujokJSWhsrJSbLN06VIsW7YMy5cvx969e6HX6zFmzBiUl5eLbZKTk7F+/XqkpaVh586dqKiowPjx41FfXy+2mTRpErKzs7F582Zs3rwZ2dnZmDJlirMfn4iI2klJhQUXq2ohkwFRQc6PEAG80ozazql1iADbHKGqqioAwA8//IDHHnsMABAQECCOHDlq8+bNds8/+ugjhISEICsrC7fddhsEQcBbb72FBQsW4P777wcAfPzxxwgNDcW6devw9NNPw2g04sMPP8TatWsxevRoAMCnn36KiIgI/PDDDxg7diwOHz6MzZs3IzMzE0OGDAEArF69GomJiThy5Ah69+7drDaz2Qyz+dLiXs5+NiIics7Rc7Z/6HYP0MJT7fhNXZsK8FLj7MVqXmlGTnN6hGj48OGYNWsW/t//+3/Ys2cP7r77bgDA77//jvDw8Gsqxmg0ArCFK8A2GlVUVISkpCSxjUajwYgRI7Br1y4AQFZWFmpra+3aGAwGxMTEiG1+/vln6HQ6MQwBwNChQ6HT6cQ2l0tNTRVPr+l0Oi46SUTUzo4W2+YPRYc4f7qsUePijKUVDETkHKcD0fLly6FUKvGf//wHK1asQLdu3QDYLs0fN25cmwsRBAGzZs3C8OHDxavaioqKAAChoaF2bUNDQ8V9RUVFUKvV8Pf3v2KbkJCQZu8ZEhIitrnc/PnzYTQaxUd+fn6bPxsREV3d7w0jRNGhPm0+RoDWtk4eR4jIWQ6fMtuyZQtuv/12dO/eHd9++22z/W+++eY1FfLcc8/hwIEDdlewNbp8cS5BEK66YNflbVpqf6XjaDQaaDTOrZJKRERtd/ScbYSoV+i1jxBd4O07yEkOjxBNnz4dwcHBmDhxIj7//HPx9JYrzJw5Exs2bMDWrVvtTrvp9XoAaDaKU1xcLI4a6fV6WCyWZpf8X97m3Llzzd73/PnzzUafiIjo+hMEAb8XN4wQhbR9hChQvJ8ZR4jIOQ4HohMnTmD79u2IjY3Fm2++idDQUIwaNQrvvPMOTp482aY3FwQBzz33HL788kv8+OOPiIqKstsfFRUFvV6P9PR0cZvFYkFGRgaGDRsGAIiPj4dKpbJrU1hYiEOHDoltEhMTYTQasWfPHrHN7t27YTQaxTZERCSd8xVmXKyqhVwG9HTBHKILPGVGTnLqKrMBAwZgwIABeOWVV1BQUIANGzZgw4YNmDdvHnr16oV7770XEyZMQEJCgkPHmzFjBtatW4evv/4aPj4+4kiQTqeDp6cnZDIZkpOTkZKSgujoaERHRyMlJQVarRaTJk0S206bNg2zZ89GYGAgAgICMGfOHMTGxopXnfXt2xfjxo3DU089hZUrVwIA/vznP2P8+PEtXmFGRETXV+Ppsu4BWnio2naFGQAE8PYd1EZOT6puZDAYMH36dGzcuBElJSX461//ipMnT2LcuHFISUlx6BgrVqyA0WjEyJEjERYWJj6++OILsc3cuXORnJyMZ599FgkJCTh79iy2bNkCH59LQ6pvvvkm7rvvPjz00EO45ZZboNVq8c0330ChuPQ/1WeffYbY2FgkJSUhKSkJAwYMwNq1a9v68YmIyIVcMaEauDRCxFNm5CyZ4OIbvlitVpSWliI4ONiVh5WcyWSCTqeD0WiEr6+v1OUQEXUpL68/iHW7T2PG7TfixbF92nyco+fKMebN7fDTqpD916Srv4C6PEe/v506ZVZaWooDBw5g4MCBCAgIQElJCT788EOYzWY8+OCD6Nu3L+RyeZcLQ0RE1L4aF2Xs5aIRImN1LerqrVAq2nwihNyMw4Foz549SEpKgslkgp+fH9LT0/Hggw9CqVRCEAQsXrwYO3bsQHx8fHvWS0REXYwgCDhSdO1XmAGAn6eq4Zi2UBTozeVTyDEOR+cFCxbgwQcfhNFoxMsvv4z77rsPo0aNwu+//46jR49i0qRJ+Nvf/taetRIRURd0pqwappo6qBSya7rCDACUCjl0nlyckZzncCDKysrCrFmz4OPjgxdeeAEFBQV46qmnxP0zZszA3r1726VIIiLqug4X2u4VGR3iA7Xy2k9xBfD2HdQGDv/Js1gs8PT0BACoVCpotVoEBQWJ+wMDA1FaWur6ComIqEvLbQhE/QyuuWClMRBxhIic4XAgioiIwIkTJ8TnaWlpCAsLE58XFhbaBSQiIiJH5BbYAlHfMNcEIn8tb99BznN4UvXDDz+M4uJi8XnjXe4bbdiwATfffLPrKiMiIrcgjhC5KBAFeHEOETnP4UC0cOHCK+5fsGCB3UKIREREV2OsrsWZsmoArgtEl27wykBEjnNqHaIr0Wq1rjoUERG5icYJ1d38PKHTqlxyTN6+g9rC6en8TU+bERERXYvG+UOumlANXJpUzUBEznAqEOXl5WH48OHtVQsREbmZwy6ePwTwKjNqG4cD0aFDh3Drrbfi8ccfb8dyiIjInRxy8RVmAOcQUds4FIh27dqF2267DVOnTsXLL7/c3jUREZEbqLLU4UiRLRANivBz2XEb5xDxjvfkDIcCUVJSEqZMmYLXX3+9veshIiI3ceisCVYBCPXVQK/zcNlxG0eIKi31qKmtd9lxqWtzKBB5eXmhsLAQgiC0dz1EROQmfs2/CAAYGO7n0uP6eiihlMsAcB4ROc6hQLRz507s27cPTzzxRHvXQ0REbiL7zEUAwKDufi49rkwm4zwicppDgSg6Oho7d+5EVlYWZsyY0d41ERGRG2gcIRrk4hEioOk8It6+gxzj8FVmBoMB27dvx/79+9uzHiIicgMlFWacKauGTAbEhOtcfnz/htt3XOApM3KQU+sQ+fv743//+1971UJERG6icXToxmBv+Hq4ZoXqpsTFGSvMLj82dU1Or1Tt6enZHnUQEZEbaa8J1Y3EQFTFU2bkmGu6l1lFRQWsVqvdNl9f1y2uRUREXVPW6TIAwKAI158uA7gWETnP6RGivLw83H333fDy8oJOp4O/vz/8/f3h5+cHf3//9qiRiIi6EEudFVmnbIHo5qjAdnkP8SozziEiBzk9QjR58mQAwD//+U+EhoZCJpO5vCgiIuq6DhUYUVNrhZ9WhegQ73Z5D/F+ZhwhIgc5HYgOHDiArKws9O7duz3qISKiLm5v3gUAwE09AiCXt88/qv21XIeInOP0KbObbroJ+fn57VELERG5gT0NgejmHgHt9h4BXJiRnOT0CNEHH3yA6dOn4+zZs4iJiYFKZX+55IABA1xWHBERdS1Wq4C9JxsCUVT7B6KyKgsEQeD0DroqpwPR+fPncfz4cbvbeMhkMvEPXH09b6RHREQtO3KuHKaaOmjVCvQ3tN9VyY2nzGrrBVSY6+DTDmsdUdfidCB68sknERcXh88//5yTqomIyCmNp8viI/2hVDg9a8NhnmoFPFUKVNfWo6yyloGIrsrpQHTq1Cls2LABPXv2bI96iIioC9txtAQAMPSG9rncvqkALzXOXqxGaaUZ3QO17f5+1Lk5Hc/vuOMO/Prrr+1RCxERdWG19Vb8fNwWiG6LDm7392s6j4joapweIbrnnnvwl7/8BQcPHkRsbGyzSdUTJkxwWXFERNR1/HKqDJWWegR4qdt1/lAjcXFG3vGeHOB0IJo+fToA4LXXXmu2j5OqiYioNY2ny4b3DGq39YeaCtDa/sHOxRnJEU4HosvvXUZEROSI7UfPAwBu69X+p8sA3r6DnNN+U/yJiIgaXKi04OBZIwDgtuig6/KegY2BqIKBiK7OoUCUlpbm8AHz8/Px008/tbkgIiLqenYcPQ9BAProfRDi63Fd3pMjROQMhwLRihUr0KdPHyxZsgSHDx9utt9oNGLjxo2YNGkS4uPjceHCBZcXSkREndeW3HMAgNv7hFy39wzQ8gav5DiH5hBlZGTg22+/xbvvvouXX34ZXl5eCA0NhYeHB8rKylBUVITg4GA88cQTOHToEEJCrt8feCIi6thqauux7bdiAMDY/vrr9r4cISJnODypevz48Rg/fjxKS0uxc+dOnDx5EtXV1QgKCkJcXBzi4uIgl3NKEhER2dt1vASVlnrofT0woJvuur2vuA4RR4jIAU5fZRYYGIh77723PWohIqIuaEuO7XRZUv/Q63K5faPGQHSxuhb1VgGK6/je1PlIOqSzfft23HPPPTAYDJDJZPjqq6/s9guCgEWLFsFgMMDT0xMjR45ETk6OXRuz2YyZM2ciKCgIXl5emDBhAs6cOWPXpqysDFOmTIFOp4NOp8OUKVNw8eLFdv50RERUbxWQ3jB/KKnf9TtdBgB+nrZ1iASBq1XT1UkaiCorKzFw4EAsX768xf1Lly7FsmXLsHz5cuzduxd6vR5jxoxBeXm52CY5ORnr169HWloadu7ciYqKCowfP95ugchJkyYhOzsbmzdvxubNm5GdnY0pU6a0++cjInJ3+05eQGmlBTpPFYbcEHBd31upkMO/YXHGkgrzdX1v6nycPmXmSnfeeSfuvPPOFvcJgoC33noLCxYswP333w8A+PjjjxEaGop169bh6aefhtFoxIcffoi1a9di9OjRAIBPP/0UERER+OGHHzB27FgcPnwYmzdvRmZmJoYMGQIAWL16NRITE3HkyBH07t27xfc3m80wmy/9D2QymVz50YmI3MJX2WcBAEn9QqFqx7vbtybU1wNlVbU4ZzKjz/UdoKJOpsPOgs7Ly0NRURGSkpLEbRqNBiNGjMCuXbsAAFlZWaitrbVrYzAYEBMTI7b5+eefodPpxDAEAEOHDoVOpxPbtCQ1NVU8xabT6RAREeHqj0hE1KXV1Nbj2wOFAIA/DO4mSQ2hDWsenTPVSPL+1Hm0ORCVlJS066hJUVERACA0NNRue2hoqLivqKgIarUa/v7+V2zT0jIAISEhYpuWzJ8/H0ajUXzk5+df0+chInI3244Uo7ymDmE6DwyNCpSkhlBfDQDgnJGBiK7MqUB08eJFzJgxA0FBQQgNDYW/vz/0ej3mz5+PqqqqdilQJrO/KkAQhGbbLnd5m5baX+04Go0Gvr6+dg8iInLc+v2202X3Dup2Xa8ua0rfOEJUzkBEV+bwHKILFy4gMTERZ8+exeTJk9G3b18IgoDDhw/j3XffRXp6Onbu3Ilff/0Vu3fvxvPPP39Nhen1tpO9RUVFCAsLE7cXFxeLo0Z6vR4WiwVlZWV2o0TFxcUYNmyY2ObcuXPNjn/+/Plmo09EROQaF6ss2Pqb7Wauf4iT5nQZAPE2IedMnFRNV+bwCNFrr70GtVqN48ePY+XKlUhOTsZf/vIXrFq1CseOHYPFYsGUKVOQlJQEne7aF96KioqCXq9Henq6uM1isSAjI0MMO/Hx8VCpVHZtCgsLcejQIbFNYmIijEYj9uzZI7bZvXs3jEaj2IaIiFzrP1lnYKm3IqabL3rrfSSrg3OIyFEOjxB99dVXWLlyZYujKnq9HkuXLsVdd92FhQsXYurUqQ4ds6KiAseOHROf5+XlITs7GwEBAejevTuSk5ORkpKC6OhoREdHIyUlBVqtFpMmTQIA6HQ6TJs2DbNnz0ZgYCACAgIwZ84cxMbGiled9e3bF+PGjcNTTz2FlStXAgD+/Oc/Y/z48a1eYUZERG1ntQr4NPMUAGDykEhJa9EzEJGDHA5EhYWF6N+/f6v7Y2JiIJfLsXDhQofffN++fbj99tvF57NmzQIATJ06FWvWrMHcuXNRXV2NZ599FmVlZRgyZAi2bNkCH59L/9p48803oVQq8dBDD6G6uhqjRo3CmjVroFAoxDafffYZnn/+efFqtAkTJrS69hEREV2bn46X4GRpFXw0Stw7yCBpLY2Tqs+Xm1FXb4VSgkv/qXOQCYIgONKwW7du+OKLLzB8+PAW9+/YsQMTJ05EQUGBSwvsKEwmE3Q6HYxGIydYExFdwZ8/2YctuecwNTESr94bI2kt9VYBvV7ZhHqrgN0vjxJPoZH7cPT72+GoPG7cOCxYsAAWS/Plz81mM/7v//4P48aNa1u1RETUJeRfqMIPh20XskweKu3pMgBQyGUI9m649J6nzegKHD5l9uqrryIhIQHR0dGYMWMG+vTpAwDIzc3Fe++9B7PZjE8++aTdCiUioo5v9Y4TsArArdFB6BUq3WTqpkJ9NSgy1aDIWIMB4VJXQx2Vw4EoPDwcP//8M5599lnMnz8fjWfaZDIZxowZg+XLl6N79+7tVigREXVsJRVmfLHXtojtMyNvlLiaS/Q6D/x6xohCLs5IV+DUvcyioqKwadMmlJWV4ejRowCAnj17IiDg+t6wj4iIOp6Pd52Euc6KgeE6JN4gzcrULQn31wIAzpS1zwLC1DW06eau/v7+uPnmm11dCxERdVLGqlqs2XUSgG106Gp3FLiewv09AQBnyqolroQ6Ml5/SERE12xFxnGU19Shj94HSf061m3lIxpGiPI5QkRXwEBERETX5JypBmt25QEAXhzbW7L7lrUmPIAjRHR1DERERHRN3v7fUdTUWhEf6Y87+oRIXU4zjXOILlbVorymVuJqqKNiICIiojY7dNaItD2nAQBzx/buUHOHGnlrlPDXqgBwlIhax0BERERtYrUKWLghB1YBGD8gDEM60JVll7t0pRkDEbWMgYiIiNrky/1nkXWqDFq1Agvu7it1OVcU0TCPKP8CJ1ZTyxiIiIjIaedMNfh/3+YCAGbeEY0wnafEFV1ZOK80o6tgICIiIqcIgoCX/nsAxupaxHbT4U+3Rkld0lV1D7AFopMllRJXQh0VAxERETnl3/vOYOuR81Ar5fj7QwOhUnT8r5KeId4AgGPnKySuhDqqjv+nmIiIOozfz5Vj4YYcAMCcpF4d5gauV9MYiM6UVaOmtl7iaqgjYiAiIiKHlNfUYvraLFTX1mN4zyBMG36D1CU5LNBLDT+tCoIAHOcoEbWAgYiIiK7KahXw4r8P4ERJJQw6D7z98CAoOtiK1Fcik8nQM7jhtFkxAxE1x0BERERXtWTzb9icUwSVQoZ/TB6MQG+N1CU5rfG02XEGImoBAxEREV3Rmp/ysHL7CQDA4vsHIK67v8QVtQ0nVtOVMBAREVGr/pN1Bq82rDf04tjeeCA+XOKK2u7GhkB09BwDETXHQERERC361958vPifXyEIwGOJkXh25I1Sl3RN+uhtV8SdKKnklWbUDAMRERE182nmKcz97wEIAjBlaCRendC/Q9641Rl6Xw8EeWtQbxWQW2iSuhzqYBiIiIhIZLUKSN10GK98dQgA8PiwHnjt3s4fhgDblWYDwnUAgINnjBJXQx0NAxEREQEAKs11eO7zX7AywzaB+oVR0Vh4T78uEYYaxXRrCERnGYjInlLqAoiISHq/FZkw47NfcPx8JVQKGZb+cQD+ENd5J1C3ZkBDIMrOvyhtIdThMBAREbkxq1XAZ3tO42/f5sJcZ4Xe1wPvTorDTT0CpC6tXcRH2pYMOFZcgZIKM4I64XpK1D54yoyIyE3llVTikdWZ+L+vDsFcZ8WIXsH47vnhXTYMAYC/l1q82mz3iQsSV0MdCUeIiIjcTJWlDqu35+G9bcdgrrPCU6XAi2N74/FhPSDvRLfjaKuhNwTit6Jy7M4rxd0DwqQuhzoIBiIiIjdhtQr4cv9ZvPH9bzhnMgMAbukZiMX3D0BEgFbi6q6foTcEYs2uk9j++3kIgtClJo1T2zEQERF1cbX1VnzzawHe23ZcvLFpuL8n5o3rg/EDwtwuEAyPDoJaIcfJ0iocK65AdKiP1CVRB8BARETURRmra/HlL2fw4c48nCmrBgD4eCjx3O09MXVYD3ioFBJXKA1vjRLDegZi25Hz2JJ7joGIADAQERF1KYIg4NczRqzbfQobfi1ATa0VABDopca0W6Pw6NBI+HqoJK5Sekn99Nh25Dy++bUAz4680e1Gyag5BiIioi7g93Pl+ObXAnzzawFOllaJ23uH+uDRod3xx/gIeKrdc0SoJXfHhmHRNzn4ragcB84YMTDCT+qSSGIMREREnZC5rh5788qw7Ugxth4pxvHzleI+D5Uc4/rr8ejQSMRH+nP0owU6rQp3xejxVXYBPs08xUBEDERERJ2Bpc6KnAIj9p0sw+68Uuw6Xooqy6U7tqsUMozoFYJ7BoZhdN9QeGn41/vVTEnsga+yC7B+/1m8MDoa4f7uc6UdNcf/Y4iIOhirVcCpC1XILTAhp8CIX06XITv/ojgfqFGwjwYjegXj9t4hGB4dBJ0n5wY5Iz7SH7f0DMRPx0qxbMvvWDZxkNQlkYQYiIiIJFJvFVBwsRp5JZXIK6nEseIKHC404XChCZVNRn8a+WlVSIgMQEIPfwzvGYR+Yb5usZBie3pxbB/sOv4Tvtx/FvcPDsfw6CCpSyKJMBAREbUTq1VASYUZZy9Wo9BYg4KGn2fKqpBXUomTpVWw1FlbfK1GKUefMF/0C/PFgHAdburhjxuCvBmAXGxQhB8eHRKJtZmn8ELafnw14xa3WqSSLmEgIiJygqXOCmN1LS5UWlBaYUZJw88LlRaUVNh+L6204JypBudMNaitF654PLVCjh5BWkQFeSEqyBt9w3zQL8wXUUFeUCp4u8nr4eW7+iLrVBlyC02YuPJnfDLtZvQM4dpE7kYmCMKV/2/tQt577z288cYbKCwsRP/+/fHWW2/h1ltvdei1JpMJOp0ORqMRvr6+7VwpEblKXb0V1bX1qLbUo7q2HlUNP2ssl36vttSjylKHCnMdymvqYKqphamm4ffqWpSLz2ubzeO5GrkMCPX1QJjOA2F+nujm5wmDzgNRwd64IcgLBj9PKDjqI7lCYzUmf7AbJ85XwkMlx19G98KUxEho1Rw36Owc/f52m0D0xRdfYMqUKXjvvfdwyy23YOXKlfjggw+Qm5uL7t27X/X1DETUFQmCAEEArIIAAQ0/BdgeEGAVbPNc7B6CAKtVQN3l2xv21VutqLcCdVYrrFbYbbPbJwioqxdsP622Y1rqBVjqrKittz0sdVZYmvxeWy/AIv7etI2A2iZta5oEoKuN0LSFTAb4a9UI8FIj0EuNIG8NAr3VCPRq/KlGsI8GYX6eCPXRcKSnkyipMOOFtP346VgpAMDXQ4mk/noM7xmEmG46dA/QQq3kf8vOhoHoMkOGDMHgwYOxYsUKcVvfvn1x3333ITU1tVl7s9kMs9ksPjeZTIiIiHB5IPpwZx5+LyoHYPsCAmxfRo2a/sex3y40a2DfVmi2vfXjNm+Llt7r8mO00MaR2nHVtm2r/Up1XUvtLfWFs7ULDb9YG4OG1bbt8kDS+LxpOLkUUoTLttvexdqwr/Fni69vKOjy4ONOZDJAq1LAU93wUDU81Jd++mhU8PFQwsdDBV/Php9Nnvt6qODroYK3h5KjOl2U1SrgP7+cwfIfj+H0hSq7fTIZEOStQZC3Bt4aBbRqJbRqBdRKORQyGWQyGRRyQCGXQS6TiT/bqi0vlaFt79dRlqqaNKQ7bgz2dukxHQ1EbjEWaLFYkJWVhZdeeslue1JSEnbt2tXia1JTU/Hqq6+2e20Zv5/H9t/Pt/v7ELmCUi6DXC6DUi6DQiaDQtHwU37pL3/lZdvsHi1tb7JNrZRDpbA9NEo5VAoZVAq5uF3d5HeVwtZe3dBebKOUQaO0BRxtQ9jxUCmgUcq5QCFdlVwuw0MJEXhgcDj2nbyALbnnkJ1/ETkFRtTUWnG+3Izz5earH4jaZESvYJcHIke5RSAqKSlBfX09QkND7baHhoaiqKioxdfMnz8fs2bNEp83jhC52gODu2FIVECz7U3/3m6a+Bu3y67StrW/95t+IcjEbWi2rWnb1vY37rhaLc7UjavV0trr29AHzvRha21x1baX3ksutx1bJrNtl8tsz+UyW4PG38V9Db/LAMhlttc1/muz6XPbT9u7yC87tqzJcS5/P/llx5Y1PWaT7Y2hhVc3kTtRyGUYckMghtwQCMA2+nqh0oJCYw1KKy2ottSh0lyPSksdLHVWCELj6WHb6d/6hlFbq7VtQ7ECnH9dW0d9O9JgsZRX+LlFIGp0+b8OBUFo9V+MGo0GGo2m3Wu6d1C3dn8PIiK6NjKZDIHeGgR6t//3AknDLWaHBQUFQaFQNBsNKi4ubjZqRERERO7HLQKRWq1GfHw80tPT7banp6dj2LBhElVFREREHYXbnDKbNWsWpkyZgoSEBCQmJmLVqlU4ffo0pk+fLnVpREREJDG3CUQTJ05EaWkpXnvtNRQWFiImJgYbN25EZGSk1KURERGRxNxmHaJrZTQa4efnh/z8fC7MSERE1Ek0XiV+8eJF6HS6Vtu5zQjRtSovty2e2B6X3hMREVH7Ki8vv2Ig4giRg6xWKwoKCuDj4+PSxd0akytHnq6OfeUc9pfj2FeOY185jn3luPbsK0EQUF5eDoPBALm89WvJOELkILlcjvDw8HY7vq+vL/+HcRD7yjnsL8exrxzHvnIc+8px7dVXVxoZauQWl90TERERXQkDEREREbk9BiKJaTQaLFy48LrcJqSzY185h/3lOPaV49hXjmNfOa4j9BUnVRMREZHb4wgRERERuT0GIiIiInJ7DERERETk9hiIiIiIyO0xEBEREZHbYyCS2HvvvYeoqCh4eHggPj4eO3bskLokyW3fvh333HMPDAYDZDIZvvrqK7v9giBg0aJFMBgM8PT0xMiRI5GTkyNNsRJLTU3FTTfdBB8fH4SEhOC+++7DkSNH7Nqwv2xWrFiBAQMGiCvhJiYmYtOmTeJ+9lPrUlNTIZPJkJycLG5jf9ksWrQIMpnM7qHX68X97Cd7Z8+exaOPPorAwEBotVoMGjQIWVlZ4n4p+4uBSEJffPEFkpOTsWDBAuzfvx+33nor7rzzTpw+fVrq0iRVWVmJgQMHYvny5S3uX7p0KZYtW4bly5dj79690Ov1GDNmjHgDXneSkZGBGTNmIDMzE+np6airq0NSUhIqKyvFNuwvm/DwcCxevBj79u3Dvn37cMcdd+Dee+8V/7JlP7Vs7969WLVqFQYMGGC3nf11Sf/+/VFYWCg+Dh48KO5jP11SVlaGW265BSqVCps2bUJubi7+/ve/w8/PT2wjaX8JJJmbb75ZmD59ut22Pn36CC+99JJEFXU8AIT169eLz61Wq6DX64XFixeL22pqagSdTie8//77ElTYsRQXFwsAhIyMDEEQ2F9X4+/vL3zwwQfsp1aUl5cL0dHRQnp6ujBixAjhhRdeEASBf66aWrhwoTBw4MAW97Gf7M2bN08YPnx4q/ul7i+OEEnEYrEgKysLSUlJdtuTkpKwa9cuiarq+PLy8lBUVGTXbxqNBiNGjGC/ATAajQCAgIAAAOyv1tTX1yMtLQ2VlZVITExkP7VixowZuPvuuzF69Gi77ewve0ePHoXBYEBUVBQefvhhnDhxAgD76XIbNmxAQkICHnzwQYSEhCAuLg6rV68W90vdXwxEEikpKUF9fT1CQ0PttoeGhqKoqEiiqjq+xr5hvzUnCAJmzZqF4cOHIyYmBgD763IHDx6Et7c3NBoNpk+fjvXr16Nfv37spxakpaXhl19+QWpqarN97K9LhgwZgk8++QTff/89Vq9ejaKiIgwbNgylpaXsp8ucOHECK1asQHR0NL7//ntMnz4dzz//PD755BMA0v+5Urb7O9AVyWQyu+eCIDTbRs2x35p77rnncODAAezcubPZPvaXTe/evZGdnY2LFy/iv//9L6ZOnYqMjAxxP/vJJj8/Hy+88AK2bNkCDw+PVtuxv4A777xT/D02NhaJiYm48cYb8fHHH2Po0KEA2E+NrFYrEhISkJKSAgCIi4tDTk4OVqxYgccee0xsJ1V/cYRIIkFBQVAoFM1Sb3FxcbN0TJc0Xr3BfrM3c+ZMbNiwAVu3bkV4eLi4nf1lT61Wo2fPnkhISEBqaioGDhyIt99+m/10maysLBQXFyM+Ph5KpRJKpRIZGRl45513oFQqxT5hfzXn5eWF2NhYHD16lH+uLhMWFoZ+/frZbevbt694IZHU/cVAJBG1Wo34+Hikp6fbbU9PT8ewYcMkqqrji4qKgl6vt+s3i8WCjIwMt+w3QRDw3HPP4csvv8SPP/6IqKgou/3srysTBAFms5n9dJlRo0bh4MGDyM7OFh8JCQmYPHkysrOzccMNN7C/WmE2m3H48GGEhYXxz9VlbrnllmbLgvz++++IjIwE0AH+vmr3advUqrS0NEGlUgkffvihkJubKyQnJwteXl7CyZMnpS5NUuXl5cL+/fuF/fv3CwCEZcuWCfv37xdOnTolCIIgLF68WNDpdMKXX34pHDx4UHjkkUeEsLAwwWQySVz59ffMM88IOp1O2LZtm1BYWCg+qqqqxDbsL5v58+cL27dvF/Ly8oQDBw4IL7/8siCXy4UtW7YIgsB+upqmV5kJAvur0ezZs4Vt27YJJ06cEDIzM4Xx48cLPj4+4t/j7KdL9uzZIyiVSuH1118Xjh49Knz22WeCVqsVPv30U7GNlP3FQCSxf/zjH0JkZKSgVquFwYMHi5dLu7OtW7cKAJo9pk6dKgiC7dLMhQsXCnq9XtBoNMJtt90mHDx4UNqiJdJSPwEQPvroI7EN+8vmySefFP9fCw4OFkaNGiWGIUFgP13N5YGI/WUzceJEISwsTFCpVILBYBDuv/9+IScnR9zPfrL3zTffCDExMYJGoxH69OkjrFq1ym6/lP0lEwRBaP9xKCIiIqKOi3OIiIiIyO0xEBEREZHbYyAiIiIit8dARERERG6PgYiIiIjcHgMRERERuT0GIiIiInJ7DERERETk9hiIiIiIyO0xEBEREZHbU0pdQGdhtVpRUFAAHx8fyGQyqcshIiIiBwiCgPLychgMBsjlrY8DMRA5qKCgABEREVKXQURERG2Qn5+P8PDwVvczEDnIx8cHgK1DfX19Ja6GiIiIHGEymRARESF+j7eGgchBjafJfH19GYiIiIg6matNd+GkaiIiInJ7DERERETk9hiIiIiIyO1xDhERdQmCIKC2XkB1bT1qGh7VtfWottSL26otVrv95jor6uoF1FmtqK0XUFdvRZ1VQG29FfVW2/HqrLY2AgTI0DAHQQbIYJuTYPtpey6Xy6BRKqBRyqFRyaFRKuDR8FOjlEOrVsDXQwVfTxV8PZXi715qBZfzIJIYAxERScpcV4+KmjpUmOtQ3vCzoqYOlRb75xVm+9/LzXWoqKlFhbkOlWZb6Km3ClJ/nDZRyGXw16oQ4uOBEF8NQht+hvh6oJufByIDvRDu7wmNUiF1qURdFgMREV2zmtp6GKtrLz2qLv1+sboWpqb7LntY6qwur0chl8FTpYCHSgFPtVz83UOlgGfDo3HkRqGQQSWXQamQQ6mQQSW3/VQ2bpPLoFLIIZMBgmAbiRLQ8Dtsz9HwvF4QYKmzwlxXD3OtFTUNP80N26os9TDV1MHUpE/qrALqrQJKKiwoqbAgt7DlzySXAQY/T/QI9EKPIC366H3Rz+CLPnofaNX8q5zoWvH/IiICYPtir7LU40KlBReranGhyoKLVRaUVVpQVtV6oHFVqPFSK+ClUcLbQwmfhp/eGiW8NSr4NP7uoYSXpmG/XRsltGoFPNQKeCgVUClkneIUlCAIqKm1wlhdi9JKM4rLzSg21aDYZPv9nKkGZ8qqcaq0EpWWepwpq8aZsmrsPHbpGDIZEBXohdhwHRJ6BODmHgGIDvGGXN7xPz9RRyITGv95Q1dkMpmg0+lgNBq5DhF1eIIgwFRTh4tVFjHglDX5vTHsNN1XVlkLS33bg41cBvh6quDnqYLO0zY3RneFh2+Tn94aJRT8Am+VINhGkE6VVuJkaRWOFVfgcKEJuYUmnC83N2uv81Thph7+GNE7BHf0CUE3P08JqibqGBz9/mYgchADEUnFahVgrG4ILQ3B5VKgqbWN4jRsb2xzscp2KqYt1Eo5ArRq+GlVCPBSw7/h9yuFGp1WBW+1kqMSEjhfbkZuoQn7T5dh78kL+OXURVTX1tu16aP3wR19QnDPQAP6hvHvL3IvDEQuxkBErlJvFcTRmtIK288LlWaUiL9bUFppFn8vq6pt82RhrVoBf60a/l4q20+tGv5aFfy0agR42Ycefy/bPk8Vr3jqzGrrrcgtMOGn4yXY+lsxsk6Voekfnz56H/whrhvui+uGUF8P6Qoluk46RSDavn073njjDWRlZaGwsBDr16/HfffdJ+4XBAGvvvoqVq1ahbKyMgwZMgT/+Mc/0L9/f7GN2WzGnDlz8Pnnn6O6uhqjRo3Ce++9Z3cDt7KyMjz//PPYsGEDAGDChAl499134efn53CtDETUmtp6K8oqLSgVw4wFFyrM4u+Noacx5FysrkVb/q/z8VCKgca/McQ0BhwvNQIu2+enVcFDxauS3F1ZpQUZv5/HpkOF2PrbefG0qEIuw7j+ejxxSw/ER/ozBFOX1SkC0aZNm/DTTz9h8ODBeOCBB5oFoiVLluD111/HmjVr0KtXL/ztb3/D9u3bceTIEfEmbc888wy++eYbrFmzBoGBgZg9ezYuXLiArKwsKBS2L4M777wTZ86cwapVqwAAf/7zn9GjRw988803DtfKQOQ+zHX1KKusRUmFucmIjW0Ux35Ux4KSCjNMNXVteh8/rQqBXmoEemkQ4KVGgLcagV62kZuAhu2BDdv8vdRQKbiOKl2bi1UWbDxYhPX7z2DvyTJxe2w3HZ4ZeSPG9dfztCd1OZ0iEDUlk8nsApEgCDAYDEhOTsa8efMA2EaDQkNDsWTJEjz99NMwGo0IDg7G2rVrMXHiRABAQUEBIiIisHHjRowdOxaHDx9Gv379kJmZiSFDhgAAMjMzkZiYiN9++w29e/d2qD4Gos6rpra+YdTGNkpzacSmSchpCDgXKiwoNzsfcOQywF/bJMx424eapiEnoOHUlJIBhyT0W5EJa346ifX7z8LccJVgvzBfzE7qhTv6hHDEiLoMR7+/O+xl93l5eSgqKkJSUpK4TaPRYMSIEdi1axeefvppZGVloba21q6NwWBATEwMdu3ahbFjx+Lnn3+GTqcTwxAADB06FDqdDrt27Wo1EJnNZpjNl67eMJlM7fApqS2sDXNwShtGaEorLCitsM3BKW2Yi1NaYbbtLzej0lJ/9YNexrZQnm105lK4USPAS2M3khPkbdum81TxKinqVProfbH4gQGYO64P1uw6iX/uzENuoQnTPt6Hm3r447V7YzgBm9xKhw1ERUVFAIDQ0FC77aGhoTh16pTYRq1Ww9/fv1mbxtcXFRUhJCSk2fFDQkLENi1JTU3Fq6++ek2fgRxXbalHSYX5UsARg01D6GkY2bFNPDbD2TnGKoWsYZRGc9lpKTUCvTWXjeqo4euh4qkDcgsBXmrMGtMLTwzrgZXbT+DjXSex92QZxr+7E1OGRmJWUi/4eqikLpOo3XXYQNTo8mFbQRCuOpR7eZuW2l/tOPPnz8esWbPE5yaTCREREY6W7fYar6RqDDRNg86lbZeeV7VhFEecg+OtQZC37XRUkLftFFWQt217435fDyVPARBdgb+XGi/d2QePJUbi9e8O47uDhViz6yTSc8/h7w8NxNAbAqUukahdddhApNfrAdhGeMLCwsTtxcXF4qiRXq+HxWJBWVmZ3ShRcXExhg0bJrY5d+5cs+OfP3++2ehTUxqNBhqNxiWfpStoXMW4tMKCkiah5tKpqsbTVg1zdKosTl9JpVbKEex9aSJxYMPvl7Y1hh3biA4nGRO5nsHPE/+YPBgPHz2PBesP4fSFKjyyOhN/Gh6FOWN7835q1GV12EAUFRUFvV6P9PR0xMXFAQAsFgsyMjKwZMkSAEB8fDxUKhXS09Px0EMPAQAKCwtx6NAhLF26FACQmJgIo9GIPXv24OabbwYA7N69G0ajUQxN7qqu3ooLDaM4TU9TNQYdW/i5FHRqap1fxdhfq7o0guOtQVCToBPUZGQn0FsNbw1HcYg6ilujg7HxhVvxt29zkbY3H6t35GHfqTK8/2g81y+iLknSQFRRUYFjxy7dlCcvLw/Z2dkICAhA9+7dkZycjJSUFERHRyM6OhopKSnQarWYNGkSAECn02HatGmYPXs2AgMDERAQgDlz5iA2NhajR48GAPTt2xfjxo3DU089hZUrVwKwXXY/fvx4h68w6ywEQUClpR4l5Wa7OTiNE4zPNwk6pZW2FY2dHcXRKOWXgkyTn4Fe6ianqxqurNKqeSUVUSfmrVFi8QMDMKpvKGb/Kxv7T1/E+Hd34v1HByM+MkDq8ohcStLL7rdt24bbb7+92fapU6dizZo14sKMK1eutFuYMSYmRmxbU1ODF198EevWrbNbmLHpfJ8LFy40W5hx+fLlnWJhxsZF/0pamFwsnqKqvDRXx+zkTTZlMiBAq252SuryOThBDdu1aq5iTOSOTpZU4um1WThyrhxqpRzvPhKHsf31UpdFdFWdbh2ijq69AtGuYyU4UVLZwoRjW9C5WFXr9DE9VQoE+TROMrYPOk1/Nq6Jw8vFicgRleY6JH+RjfTcc5DLgNT7YzHxpu5Sl0V0RZ1+HSJ3sXzrMew6XnrFNnIZxEX9GoOOGGxaOFWlVfM/KxG5npdGiRWTB2PB+kP4Yl8+5v33IKos9XjiliipSyO6ZvzmlFhCjwBo1UoENwk6jZOPg3xsgcdPy1EcIuoYlAo5Fj8QC38vNd7POI5Xv8mFVq3gSBF1egxEEps1ppfUJRAROUUmk2HeuN6wCgJWbT+Bl748CA+VAvcO6iZ1aURtxkuAiIjIaTKZDPPv7INHh3aHIABz/v0rMk9c+fQ/UUfGQERERG0ik8nw2oQY3D0gDLX1Ap75NAunSiulLouoTRiIiIiozeRyGf7+4EAMDNehrKoW0z7ehwpzndRlETmNgYiIiK6Jh0qB1Y8lQO/rgWPFFXhl/UFwRRfqbBiIiIjomoX4emD5pDgo5DJ8lV2Af2edkbokIqcwEBERkUsk9AgQr5xd+HUOjhVXSFwRkeMYiIiIyGWeGXEjhvcMQnVtPeb99wDqrTx1Rp0DAxEREbmMXC7Dkj8OgJdagaxTZVj780mpSyJyCAMRERG5VDc/T7x0V18AwNLvjyD/QpXEFRFdHQMRERG53OSbu+PmqABUWerx6je5UpdDdFUMRERE5HJyuQwpf4iBQi7DD4fP4adjJVKXRHRFDERERNQueob4YMrQSADA//s2lxOsqUNz+uauJ0+exI4dO3Dy5ElUVVUhODgYcXFxSExMhIeHR3vUSEREndQLo6Kxfv9Z/FZUjrS9pzF5SKTUJRG1yOFAtG7dOrzzzjvYs2cPQkJC0K1bN3h6euLChQs4fvw4PDw8MHnyZMybNw+RkfwDT0REgL+XGsmjo/HqN7l464ejeGBwODxUCqnLImrGoVNmgwcPxrJly/Doo4/i5MmTKCoqQlZWFnbu3Inc3FyYTCZ8/fXXsFqtSEhIwL///e/2rpuIiDqJyUMi0c3PE+fLzfhs92mpyyFqkUxw4IYz3333He6++26HDlhSUoK8vDzcdNNN11xcR2IymaDT6WA0GuHr6yt1OUREncrne05j/pcHEeStwY65t8NTzVEiuj4c/f52aITI0TAEAEFBQV0uDBER0bX5Y3w4wv09UVJhxme7T0ldDlEzbb7KrLi4GIcOHcKBAwfsHq5UV1eHV155BVFRUfD09MQNN9yA1157DVarVWwjCAIWLVoEg8EAT09PjBw5Ejk5OXbHMZvNmDlzJoKCguDl5YUJEybgzBneeJCI6HpRKeR47vaeAIAPd+ahtt56lVcQXV9OB6KsrCzExMQgLCwMAwYMwKBBgxAXFyf+dKUlS5bg/fffx/Lly3H48GEsXboUb7zxBt59912xzdKlS7Fs2TIsX74ce/fuhV6vx5gxY1BeXi62SU5Oxvr165GWloadO3eioqIC48ePR319vUvrJSKi1v1hcDcEeWtQaKzBdwcKpS6HyI5Dc4iaGjBgAHr27Il58+YhNDQUMpnMbr8rrzAbP348QkND8eGHH4rbHnjgAWi1WqxduxaCIMBgMCA5ORnz5s0DYBsNCg0NxZIlS/D000/DaDQiODgYa9euxcSJEwEABQUFiIiIwMaNGzF27FiHauEcIiKia/fu/47i7+m/o7/BF9/OHN7sO4TI1Vw6h6ipvLw8LF26FEOGDEGPHj0QGRlp93Cl4cOH43//+x9+//13AMCvv/6KnTt34q677hJrKSoqQlJSkvgajUaDESNGYNeuXQBsI1q1tbV2bQwGA2JiYsQ2LTGbzTCZTHYPIiK6No8OjYSHSo6cAhN+PlEqdTlEIqcD0ahRo/Drr7+2Ry3NzJs3D4888gj69OkDlUqFuLg4JCcn45FHHgEAFBUVAQBCQ0PtXhcaGiruKyoqglqthr+/f6ttWpKamgqdTic+IiIiXPnRiIjckr+XGn+MDwcArP2Zk6up43B6peoPPvgAU6dOxaFDhxATEwOVSmW3f8KECS4r7osvvsCnn36KdevWoX///sjOzkZycjIMBgOmTp0qtrt8yFUQhKsOw16tzfz58zFr1izxuclkYigiInKBR4dG4tPM00jPPYfi8hqE+PAuByQ9pwPRrl27sHPnTmzatKnZPplM5tKJyi+++CJeeuklPPzwwwCA2NhYnDp1CqmpqZg6dSr0ej0A2yhQWFiY+Lri4mJx1Eiv18NisaCsrMxulKi4uBjDhg1r9b01Gg00Go3LPgsREdn00fsirrsf9p++iH/vO4MZDVefEUnJ6VNmzz//PKZMmYLCwkJYrVa7h6uv2qqqqoJcbl+iQqEQL7uPioqCXq9Henq6uN9isSAjI0MMO/Hx8VCpVHZtCgsLcejQoSsGIiIiaj+Tbu4OAEjbexpW3vSVOgCnR4hKS0vxl7/8pdm8nfZwzz334PXXX0f37t3Rv39/7N+/H8uWLcOTTz4JwDYilZycjJSUFERHRyM6OhopKSnQarWYNGkSAECn02HatGmYPXs2AgMDERAQgDlz5iA2NhajR49u989ARETNjR9gwGvf5iL/QjV+Ol6CW6ODpS6J3JzTgej+++/H1q1bceONN7ZHPXbeffdd/N///R+effZZFBcXw2Aw4Omnn8Zf//pXsc3cuXNRXV2NZ599FmVlZRgyZAi2bNkCHx8fsc2bb74JpVKJhx56CNXV1Rg1ahTWrFkDhYJLxxMRScFTrcC9gwz4NPM01u8/y0BEknN6HaLXX38db731Fu6++27ExsY2m1T9/PPPu7TAjoLrEBERuda+kxfwx/d/hpdagX2vjOH9zahdOPr97XQgioqKav1gMhlOnDjhzOE6DQYiIiLXEgQBty7dijNl1Xj3kTjcM9AgdUnUBTn6/e30KbO8vLxrKoyIiAiw/SN6wkAD3tt2HF9nFzAQkaTafHNXIiKia3VfXDcAQMbvxbhYZZG4GnJnDgWixYsXo6qqyqED7t69G9999901FUVERO6hV6gP+uh9UFsvYEvuOanLITfmUCDKzc1F9+7d8cwzz2DTpk04f/68uK+urg4HDhzAe++9h2HDhuHhhx/mHBsiInLYnTG2hXW35DAQkXQcCkSffPIJfvzxR1itVkyePBl6vR5qtRo+Pj7QaDSIi4vDP//5Tzz++OP47bffcOutt7Z33URE1EUk9beta7fj6HlUWeokrobcldNXmQmCgAMHDuDkyZOorq5GUFAQBg0ahKCgoPaqsUPgVWZERO1DEATc9sZW5F+oxvuPDsa4mLCrv4jIQe12lZlMJsPAgQMxcODAayqQiIgIsH2vjO2nxwc78/B9zjkGIpIErzIjIiLJJfW33az7f4fPobbeKnE15I4YiIiISHLxkf4I8FLDVFOHrFNlUpdDboiBiIiIJKeQy3BbtG0u6vbfz1+lNZHrMRAREVGHcFsv2w1eMxiISAIMRERE1CE03vE+p8CE8+Vmiashd+PQVWb333+/wwf88ssv21wMERG5r2AfDfobfJFTYMLOY+fxh7hwqUsiN+JQINLpdO1dBxEREW7rFYycAhMyjjAQ0fXlUCD66KOP2rsOIiIi3BYdjBXbjmPH0RJYrQLkcpnUJZGbaNMcorq6Ovzwww9YuXIlysvLAQAFBQWoqKhwaXFERORe4iP94alSoLTSgqPF/E6h68fplapPnTqFcePG4fTp0zCbzRgzZgx8fHywdOlS1NTU4P3332+POomIyA2olXIk9PDHjqMlyDxRit56H6lLIjfh9AjRCy+8gISEBJSVlcHT01Pc/oc//AH/+9//XFocERG5nyFRAQCAzBOlEldC7sTpEaKdO3fip59+glqtttseGRmJs2fPuqwwIiJyT0NvCAQA7M67AEEQIJNxHhG1P6dHiKxWK+rr65ttP3PmDHx8OLRJRETXZkC4HzxUclzgPCK6jpwORGPGjMFbb70lPpfJZKioqMDChQtx1113ubI2AMDZs2fx6KOPIjAwEFqtFoMGDUJWVpa4XxAELFq0CAaDAZ6enhg5ciRycnLsjmE2mzFz5kwEBQXBy8sLEyZMwJkzZ1xeKxERXTu1Uo74SH8AwG6eNqPrxOlA9OabbyIjIwP9+vVDTU0NJk2ahB49euDs2bNYsmSJS4srKyvDLbfcApVKhU2bNiE3Nxd///vf4efnJ7ZZunQpli1bhuXLl2Pv3r3Q6/UYM2aMePUbACQnJ2P9+vVIS0vDzp07UVFRgfHjx7c40kVERNIbGmU7bZZ54oLElZC7kAmCIDj7ourqanz++ef45ZdfYLVaMXjwYEyePNlukrUrvPTSS/jpp5+wY8eOFvcLggCDwYDk5GTMmzcPgG00KDQ0FEuWLMHTTz8No9GI4OBgrF27FhMnTgRgWyIgIiICGzduxNixY1s8ttlshtl8ael4k8mEiIgIGI1G+Pr6uvRzEhGRvT15F/DQyp8R5K3B3gWjOI+I2sxkMkGn0131+9vpEaKqqip4enriySefxPLly/Hee+/hT3/6k8vDEABs2LABCQkJePDBBxESEoK4uDisXr1a3J+Xl4eioiIkJSWJ2zQaDUaMGIFdu3YBALKyslBbW2vXxmAwICYmRmzTktTUVOh0OvERERHh8s9HREQtGxCug0ohQ0mFGfkXqqUuh9yA04EoJCQEjz76KL7//ntYrdb2qEl04sQJrFixAtHR0fj+++8xffp0PP/88/jkk08AAEVFRQCA0NBQu9eFhoaK+4qKiqBWq+Hv799qm5bMnz8fRqNRfOTn57vyoxER0RV4qBTob7DdNuqX02USV0PuwOlA9Mknn8BsNuMPf/gDDAYDXnjhBezdu7c9ahNPx6WkpCAuLg5PP/00nnrqKaxYscKu3eVDqY5cpnm1NhqNBr6+vnYPIiK6fgZ3t/1DloGIrgenA9H999+Pf//73zh37hxSU1Nx+PBhDBs2DL169cJrr73m0uLCwsLQr18/u219+/bF6dOnAQB6vR4Amo30FBcXi6NGer0eFosFZWVlrbYhIqKOZ3CkHwAg6xQDEbW/Nt3LDAB8fHzwxBNPYMuWLfj111/h5eWFV1991ZW14ZZbbsGRI0fstv3++++IjIwEAERFRUGv1yM9PV3cb7FYkJGRgWHDhgEA4uPjoVKp7NoUFhbi0KFDYhsiIup4Gi+9/62oHFWWOomroa6uzYGopqYG//rXv3Dfffdh8ODBKC0txZw5c1xZG/7yl78gMzMTKSkpOHbsGNatW4dVq1ZhxowZAGynypKTk5GSkoL169fj0KFDePzxx6HVajFp0iQAgE6nw7Rp0zB79mz873//w/79+/Hoo48iNjYWo0ePdmm9RETkOmE6T4TpPFBvFfBrvlHqcqiLc/rWHVu2bMFnn32Gr776CgqFAn/84x/x/fffY8SIES4v7qabbsL69esxf/58vPbaa4iKisJbb72FyZMni23mzp2L6upqPPvssygrK8OQIUOwZcsWu1Wz33zzTSiVSjz00EOorq7GqFGjsGbNGigUCpfXTERErjO4uz++O1iIX06XIfHGQKnLoS7M6XWItFot7r77bkyePBl33303VCpVe9XWoTi6jgEREbnOBztO4G/fHcaoPiH48PGbpC6HOiFHv7+dHiEqKipiICAiousiruFKs1/PGHmjV2pXTs8h8vX1xfHjx/HKK6/gkUceQXFxMQBg8+bNze4hRkREdC36hflCIbct0FhkqpG6HOrCnA5EGRkZiI2Nxe7du/Hll1+iosJ2J+IDBw5g4cKFLi+QiIjcl6dagegQbwDAgTOcWE3tx+lA9NJLL+Fvf/sb0tPToVarxe233347fv75Z5cWR0RENCDctmL1QQYiakdOB6KDBw/iD3/4Q7PtwcHBKC0tdUlRREREjWLD/QAAB84yEFH7cToQ+fn5obCwsNn2/fv3o1u3bi4pioiIqNGAbo0jRBfh5IXRRA5zOhBNmjQJ8+bNQ1FREWQyGaxWK3766SfMmTMHjz32WHvUSEREbqxPmA9UChnKqmpxpqxa6nKoi3I6EL3++uvo3r07unXrhoqKCvTr1w+33XYbhg0bhldeeaU9aiQiIjemUSrQW29bbPcgT5tRO3F6HSKVSoXPPvsMr732Gvbv3w+r1Yq4uDhER0e3R31ERESI7eaHQ2dNOHjWiLtiw6Quh7ogpwNRoxtvvBE33nijK2shIiJq0YBwHT7fwyvNqP04FIhmzZrl8AGXLVvW5mKIiIhaEtswsfpAw8RqrlhNruZQINq/f79DB+MfUCIiag+9Qm0Tq001dThTVo2IAK3UJVEX41Ag2rp1a3vXQURE1Cq1Uo7oEB/kFpqQW2hiICKXc/oqMyIiIin0M9huLJ5bYJK4EuqKGIiIiKhT6BfWEIgKGYjI9RiIiIioU+AIEbUnBiIiIuoU+jaMEJ29WA1jVa3E1VBXw0BERESdgs5ThXB/TwA8bUau16ZAtHbtWtxyyy0wGAw4deoUAOCtt97C119/7dLiiIiImuI8ImovTgeiFStWYNasWbjrrrtw8eJF1NfXAwD8/Pzw1ltvubo+IiIiUeM8opwCrlhNruV0IHr33XexevVqLFiwAAqFQtyekJCAgwcPurS4y6WmpkImkyE5OVncJggCFi1aBIPBAE9PT4wcORI5OTl2rzObzZg5cyaCgoLg5eWFCRMm4MyZM+1aKxERuV5/g23Fak6sJldzOhDl5eUhLi6u2XaNRoPKykqXFNWSvXv3YtWqVRgwYIDd9qVLl2LZsmVYvnw59u7dC71ejzFjxqC8vFxsk5ycjPXr1yMtLQ07d+5ERUUFxo8fL45uERFR59A4QnSsuALmOv4dTq7jdCCKiopCdnZ2s+2bNm1Cv379XFFTMxUVFZg8eTJWr14Nf39/cbsgCHjrrbewYMEC3H///YiJicHHH3+MqqoqrFu3DgBgNBrx4Ycf4u9//ztGjx6NuLg4fPrppzh48CB++OGHdqmXiIjah0HnAZ2nCnVWAUfPVUhdDnUhTgeiF198ETNmzMAXX3wBQRCwZ88evP7663j55Zfx4osvtkeNmDFjBu6++26MHj3abnteXh6KioqQlJQkbtNoNBgxYgR27doFAMjKykJtba1dG4PBgJiYGLFNS8xmM0wmk92DiIikJZPJOLGa2oVD9zJr6oknnkBdXR3mzp2LqqoqTJo0Cd26dcPbb7+Nhx9+2OUFpqWl4ZdffsHevXub7SsqKgIAhIaG2m0PDQ0Vr34rKiqCWq22G1lqbNP4+pakpqbi1VdfvdbyiYjIxfoZfPHziVLOIyKXcnqE6OLFi3jqqadw6tQpFBcXo6ioCPn5+Zg2bRqOHTvm0uLy8/Pxwgsv4NNPP4WHh0er7WQymd1zQRCabbvc1drMnz8fRqNRfOTn5ztXPBERtQuOEFF7cDoQ3XXXXaipqQEABAUFISQkBABw5MgRjBw50qXFZWVlobi4GPHx8VAqlVAqlcjIyMA777wDpVIpjgxdPtJTXFws7tPr9bBYLCgrK2u1TUs0Gg18fX3tHkREJL3GidWHC0wQBEHiaqircDoQ+fv747777kNdXZ247fDhwxg5ciQeeOABlxY3atQoHDx4ENnZ2eIjISEBkydPRnZ2Nm644Qbo9Xqkp6eLr7FYLMjIyMCwYcMAAPHx8VCpVHZtCgsLcejQIbENERF1Hj1DvKFWyFFursOZsmqpy6Euwuk5RP/9738xZswYTJo0CV988QVycnIwatQoTJ48GcuWLXNpcT4+PoiJibHb5uXlhcDAQHF7cnIyUlJSEB0djejoaKSkpECr1WLSpEkAAJ1Oh2nTpmH27NkIDAxEQEAA5syZg9jY2GaTtImIqONTKeTopffGobMm5BQYERGglbok6gKcDkQeHh749ttvMXLkSDz44IPYsWMHHnvsMbzxxhvtUd9VzZ07F9XV1Xj22WdRVlaGIUOGYMuWLfDx8RHbvPnmm1AqlXjooYdQXV2NUaNGYc2aNXYLSxIRUefRL8y3IRCZMC4mTOpyqAuQCQ6cgG3pkvOioiKMHj0a48ePx+LFi8XtXXWujclkgk6ng9Fo7LKfkYios/h410ks3JCDO/qE4J+P3yR1OdSBOfr97dAIkZ+fX4tXZAmCgPfffx8rV64Ur9ri6s9ERNTeYrrZvtgOneU9zcg1HApEW7dube86iIiIHNY3zBcyGVBcbkZxeQ1CfFpfmoXIEQ4FohEjRrR3HURERA7TqpW4IcgLx89XIqfAhJDeDER0bZyeVN2oqqoKp0+fhsVisdt++c1XiYiI2kNMN50tEJ014vbeIVKXQ52c04Ho/PnzeOKJJ7Bp06YW93MOERERXQ8xBh2+zi5ADm/hQS7g9MKMycnJKCsrQ2ZmJjw9PbF582Z8/PHHiI6OxoYNG9qjRiIiomb6N6xYfaiAE6vp2jk9QvTjjz/i66+/xk033QS5XI7IyEiMGTMGvr6+SE1Nxd13390edRIREdnpb9ABAPIvVMNYVQudViVxRdSZOT1CVFlZKd6/LCAgAOfPnwcAxMbG4pdffnFtdURERK3QaVWICPAEAOQUcpSIro3Tgah37944cuQIAGDQoEFYuXIlzp49i/fffx9hYVwtlIiIrp/+YbZRopyznEdE18bpU2bJyckoLCwEACxcuBBjx47FZ599BrVajTVr1ri6PiIiolbFdPPF5pwi5HAeEV0jpwPR5MmTxd/j4uJw8uRJ/Pbbb+jevTuCgoJcWhwREdGV9O9mGyE6xCvN6Bq1eR2iRlqtFoMHD3ZFLURERE5pvNLs+PkKVFnqoFVf89cauSmH/uTMmjXL4QMuW7aszcUQERE5I8THAyE+GhSXm5FbYEJCjwCpS6JOyqFAtH//focO1tINYImIiNrTwAg/pOeew/7TFxmIqM14c1ciIurUBnf3twWi/DKpS6FOzOnL7omIiDqSuO5+AIBfTl2UtA7q3BiIiIioUxsQroNCLkORqQaFxmqpy6FOioGIiIg6Na1aib5hPgA4SkRtx0BERESdXlyEPwBg/2nOI6K2YSAiIqJOb3CkHwBg3ykGImqbDh2IUlNTcdNNN8HHxwchISG47777xPuoNRIEAYsWLYLBYICnpydGjhyJnJwcuzZmsxkzZ85EUFAQvLy8MGHCBJw5c+Z6fhQiImpHNzVcbn/wrBEV5jqJq6HOqEMHooyMDMyYMQOZmZlIT09HXV0dkpKSUFlZKbZZunQpli1bhuXLl2Pv3r3Q6/UYM2YMysvLxTbJyclYv3490tLSsHPnTlRUVGD8+PGor6+X4mMREZGLhftrERmoRb1VwJ68UqnLoU5IJgiCIHURjjp//jxCQkKQkZGB2267DYIgwGAwIDk5GfPmzQNgGw0KDQ3FkiVL8PTTT8NoNCI4OBhr167FxIkTAQAFBQWIiIjAxo0bMXbsWIfe22QyQafTwWg0wtfXt90+IxERtc38Lw/g8z35+NPwKLwyvp/U5VAH4ej3d4ceIbqc0Wi7m3FAgG1oNC8vD0VFRUhKShLbaDQajBgxArt27QIAZGVloba21q6NwWBATEyM2KYlZrMZJpPJ7kFERB3XsBttNxj/6ThHiMh5nSYQCYKAWbNmYfjw4YiJiQEAFBUVAQBCQ0Pt2oaGhor7ioqKoFar4e/v32qblqSmpkKn04mPiIgIV34cIiJysaE3BAIADheacKHSInE11Nl0mkD03HPP4cCBA/j888+b7bv8HmqCIFz1vmpXazN//nwYjUbxkZ+f37bCiYjougj20aB3qG09op+OlUhcDXU2nSIQzZw5Exs2bMDWrVsRHh4ubtfr9QDQbKSnuLhYHDXS6/WwWCwoKytrtU1LNBoNfH197R5ERNSxjewdDABIzz0ncSXU2XToQCQIAp577jl8+eWX+PHHHxEVFWW3PyoqCnq9Hunp6eI2i8WCjIwMDBs2DAAQHx8PlUpl16awsBCHDh0S2xARUdeQ1N/2D+WtvxXDUmeVuBrqTBy6271UZsyYgXXr1uHrr7+Gj4+POBKk0+ng6ekJmUyG5ORkpKSkIDo6GtHR0UhJSYFWq8WkSZPEttOmTcPs2bMRGBiIgIAAzJkzB7GxsRg9erSUH4+IiFwsLsIPwT4anC83I/NEKW7rFSx1SdRJdOhAtGLFCgDAyJEj7bZ/9NFHePzxxwEAc+fORXV1NZ599lmUlZVhyJAh2LJlC3x8fMT2b775JpRKJR566CFUV1dj1KhRWLNmDRQKxfX6KEREdB3I5TKM6ReKdbtPY+PBQgYiclinWodISlyHiIioc/j5eCkeWZ0Jb40SexaMglbdof/tT+2sS65DREREdDVDbwhAZKAWFeY6bDzY+vIqRE0xEBERUZcik8nwYLztiuTP95yWuBpyhKXOiiJjDWpqpbulFscRiYioy3kwIQLv/O8Ysk6VYfeJUgxpWLSRrp+a2nqUVJhRUmFBSbm54Xfb8/MV5ibbLDBW1wIA/j09UbxR7/XGQERERF1OqK8H/pgQjnW7T+PdH48xELlIlaUOJeUNgabxUW5pEnYuBaByc51Tx1bIZTBW1bZT5VfHQERERF3SMyNuxL/25mPnsRL8+Ns53NGn9cV43ZUgCCg31zWM1jQJNuVmnG94XlpxaV+VxblTWiqFDEHemoaH2vbT59Lz4CbP/TxVkMuvfJeJ9sRAREREXVJEgBZP3NIDq3fk4ZX1h7BlViC8NV3/a08QBBira1FSYcb5y0dv7J7bRnqcXcDSQyVvEnI0CPaxBZ1AL3WTsKNBsLcGvp7Kq95Kq6Po+n8yiIjIbf1lTC9szilC/oVqJKdlY+WUeCgkHIVoq3qrgLIqixhqSivNOH/5qE6TfbX1zq2o46VW2I3ciIHHR4Pgy557qRWdJuQ4g+sQOYjrEBERdU6/nC7Dw6syYamz4oHB4Ui9PxZqpfQXWdfWW3Gh0tIQbJqEm4bnpeI+Cy5UmmF18tvax0NpOyXlrUGQj9puVCfIW90QdmzPPdVdd6FiR7+/OUJERERd2uDu/lj20EA8//l+/PeXMzh+vgJ/uy8GMd10Ln+vlq6sKrULPZeCz8U2TCD216oQeNkoTrBP81GdQC81PFRdN+S0B44QOYgjREREndvWI8WYuW4/Khqufhp6QwBG9w3FwAg/RAZo4adVQ6WQQSaTwWoVYK6zoqa2HuU1dSirsuBClQUXqyy4UFnb8NOC0ianrEorLE5fWSWXAQFeDROMLztl1TT4BPtoEOClhkoh/chWZ+Po9zcDkYMYiIiIOr9CYzVSN/6Gbw8UtHgKSi4DlAq50xONm2rpyqrGcBPsY3/ayl+rlvTKKnfAQORiDERERF3H2YvV2HSwED8dK8HR4goUXKxudY6ORilHgJctvPh7qeCvVSPASw0/rdp2ZVWTOTlBXp3ryip3wEDkYgxERERdV229FVXmetTU1aO23goPlcL2UMqh5GmqTo2TqomIiBykUsih08qhg0rqUkgijL1ERETk9hiIiIiIyO0xEBEREZHbYyAiIiIit8dJ1Q5qvBjPZDJJXAkRERE5qvF7+2oX1TMQOai8vBwAEBERIXElRERE5Kzy8nLodK3froXrEDnIarWioKAAPj4+Ll1wy2QyISIiAvn5+Vzf6CrYV85hfzmOfeU49pXj2FeOa8++EgQB5eXlMBgMkMtbnynEESIHyeVyhIeHt9vxfX19+T+Mg9hXzmF/OY595Tj2lePYV45rr7660shQI06qJiIiIrfHQERERERuj4FIYhqNBgsXLoRGo5G6lA6PfeUc9pfj2FeOY185jn3luI7QV5xUTURERG6PI0RERETk9hiIiIiIyO0xEBEREZHbYyAiIiIit8dAJLH33nsPUVFR8PDwQHx8PHbs2CF1SZLbvn077rnnHhgMBshkMnz11Vd2+wVBwKJFi2AwGODp6YmRI0ciJydHmmIllpqaiptuugk+Pj4ICQnBfffdhyNHjti1YX/ZrFixAgMGDBAXfktMTMSmTZvE/eyn1qWmpkImkyE5OVncxv6yWbRoEWQymd1Dr9eL+9lP9s6ePYtHH30UgYGB0Gq1GDRoELKyssT9UvYXA5GEvvjiCyQnJ2PBggXYv38/br31Vtx55504ffq01KVJqrKyEgMHDsTy5ctb3L906VIsW7YMy5cvx969e6HX6zFmzBjxfnPuJCMjAzNmzEBmZibS09NRV1eHpKQkVFZWim3YXzbh4eFYvHgx9u3bh3379uGOO+7AvffeK/5ly35q2d69e7Fq1SoMGDDAbjv765L+/fujsLBQfBw8eFDcx366pKysDLfccgtUKhU2bdqE3Nxc/P3vf4efn5/YRtL+EkgyN998szB9+nS7bX369BFeeukliSrqeAAI69evF59brVZBr9cLixcvFrfV1NQIOp1OeP/99yWosGMpLi4WAAgZGRmCILC/rsbf31/44IMP2E+tKC8vF6Kjo4X09HRhxIgRwgsvvCAIAv9cNbVw4UJh4MCBLe5jP9mbN2+eMHz48Fb3S91fHCGSiMViQVZWFpKSkuy2JyUlYdeuXRJV1fHl5eWhqKjIrt80Gg1GjBjBfgNgNBoBAAEBAQDYX62pr69HWloaKisrkZiYyH5qxYwZM3D33Xdj9OjRdtvZX/aOHj0Kg8GAqKgoPPzwwzhx4gQA9tPlNmzYgISEBDz44IMICQlBXFwcVq9eLe6Xur8YiCRSUlKC+vp6hIaG2m0PDQ1FUVGRRFV1fI19w35rThAEzJo1C8OHD0dMTAwA9tflDh48CG9vb2g0GkyfPh3r169Hv3792E8tSEtLwy+//ILU1NRm+9hflwwZMgSffPIJvv/+e6xevRpFRUUYNmwYSktL2U+XOXHiBFasWIHo6Gh8//33mD59Op5//nl88sknAKT/c8W73UtMJpPZPRcEodk2ao791txzzz2HAwcOYOfOnc32sb9sevfujezsbFy8eBH//e9/MXXqVGRkZIj72U82+fn5eOGFF7BlyxZ4eHi02o79Bdx5553i77GxsUhMTMSNN96Ijz/+GEOHDgXAfmpktVqRkJCAlJQUAEBcXBxycnKwYsUKPPbYY2I7qfqLI0QSCQoKgkKhaJZ6i4uLm6VjuqTx6g32m72ZM2diw4YN2Lp1K8LDw8Xt7C97arUaPXv2REJCAlJTUzFw4EC8/fbb7KfLZGVlobi4GPHx8VAqlVAqlcjIyMA777wDpVIp9gn7qzkvLy/Exsbi6NGj/HN1mbCwMPTr189uW9++fcULiaTuLwYiiajVasTHxyM9Pd1ue3p6OoYNGyZRVR1fVFQU9Hq9Xb9ZLBZkZGS4Zb8JgoDnnnsOX375JX788UdERUXZ7Wd/XZkgCDCbzeyny4waNQoHDx5Edna2+EhISMDkyZORnZ2NG264gf3VCrPZjMOHDyMsLIx/ri5zyy23NFsW5Pfff0dkZCSADvD3VbtP26ZWpaWlCSqVSvjwww+F3NxcITk5WfDy8hJOnjwpdWmSKi8vF/bv3y/s379fACAsW7ZM2L9/v3Dq1ClBEARh8eLFgk6nE7788kvh4MGDwiOPPCKEhYUJJpNJ4sqvv2eeeUbQ6XTCtm3bhMLCQvFRVVUltmF/2cyfP1/Yvn27kJeXJxw4cEB4+eWXBblcLmzZskUQBPbT1TS9ykwQ2F+NZs+eLWzbtk04ceKEkJmZKYwfP17w8fER/x5nP12yZ88eQalUCq+//rpw9OhR4bPPPhO0Wq3w6aefim2k7C8GIon94x//ECIjIwW1Wi0MHjxYvFzanW3dulUA0OwxdepUQRBsl2YuXLhQ0Ov1gkajEW677Tbh4MGD0hYtkZb6CYDw0UcfiW3YXzZPPvmk+P9acHCwMGrUKDEMCQL76WouD0TsL5uJEycKYWFhgkqlEgwGg3D//fcLOTk54n72k71vvvlGiImJETQajdCnTx9h1apVdvul7C+ZIAhC+49DEREREXVcnENEREREbo+BiIiIiNweAxERERG5PQYiIiIicnsMREREROT2GIiIiIjI7TEQERERkdtjICIiIiK3x0BERJ3Ktm3bIJPJcPHiRUne/8cff0SfPn1gtVpbbbNo0SIMGjRIfD5nzhw8//zz16E6ImorBiIi6rBGjhyJ5ORku23Dhg1DYWEhdDqdJDXNnTsXCxYsgFzu+F+fc+fOxUcffYS8vLx2rIyIrgUDERF1Kmq1Gnq9HjKZ7Lq/965du3D06FE8+OCDTr0uJCQESUlJeP/999upMiK6VgxERNQhPf7448jIyMDbb78NmUwGmUyGkydPNjtltmbNGvj5+eHbb79F7969odVq8cc//hGVlZX4+OOP0aNHD/j7+2PmzJmor68Xj2+xWDB37lx069YNXl5eGDJkCLZt23bFmtLS0pCUlAQPDw+77YsXL0ZoaCh8fHwwbdo01NTUNHvthAkT8Pnnn19zvxBR+2AgIqIO6e2330ZiYiKeeuopFBYWorCwEBERES22raqqwjvvvIO0tDRs3rwZ27Ztw/3334+NGzdi48aNWLt2LVatWoX//Oc/4mueeOIJ/PTTT0hLS8OBAwfw4IMPYty4cTh69GirNW3fvh0JCQl22/71r39h4cKFeP3117Fv3z6EhYXhvffea/bam2++Gfn5+Th16lQbe4SI2pNS6gKIiFqi0+mgVquh1Wqh1+uv2La2thYrVqzAjTfeCAD44x//iLVr1+LcuXPw9vZGv379cPvtt2Pr1q2YOHEijh8/js8//xxnzpyBwWAAYJv4vHnzZnz00UdISUlp8X1Onjwptm/01ltv4cknn8Sf/vQnAMDf/vY3/PDDD81Gibp16yYeIzIy0vkOIaJ2xREiIur0tFqtGIYAIDQ0FD169IC3t7fdtuLiYgDAL7/8AkEQ0KtXL3h7e4uPjIwMHD9+vNX3qa6ubna67PDhw0hMTLTbdvlzAPD09ARgG80ioo6HI0RE1OmpVCq75zKZrMVtjZfKW61WKBQKZGVlQaFQ2LVrGqIuFxQUhLKysjbVeOHCBQBAcHBwm15PRO2LgYiIOiy1Wm03EdpV4uLiUF9fj+LiYtx6661OvS43N9duW9++fZGZmYnHHntM3JaZmdnstYcOHYJKpUL//v3bXjgRtRueMiOiDqtHjx7YvXs3Tp48iZKSkisuhuiMXr16YfLkyXjsscfw5ZdfIi8vD3v37sWSJUuwcePGVl83duxY7Ny5027bCy+8gH/+85/45z//id9//x0LFy5ETk5Os9fu2LEDt956q3jqjIg6FgYiIuqw5syZA4VCgX79+iE4OBinT5922bE/+ugjPPbYY5g9ezZ69+6NCRMmYPfu3a1eyQYAjz76KHJzc3HkyBFx28SJE/HXv/4V8+bNQ3x8PE6dOoVnnnmm2Ws///xzPPXUUy6rn4hcSyYIgiB1EUREncXcuXNhNBqxcuVKh1/z3Xff4cUXX8SBAwegVHKmAlFHxBEiIiInLFiwAJGRkU7NbaqsrMRHH33EMETUgXGEiIiIiNweR4iIiIjI7TEQERERkdtjICIiIiK3x0BEREREbo+BiIiIiNweAxERERG5PQYiIiIicnsMREREROT2GIiIiIjI7f1/DMey5k7A93EAAAAASUVORK5CYII=",
"text/plain": [
"Figure(PyObject <Figure size 640x480 with 2 Axes>)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"PyObject Text(24.0, 0.5, 'lake level (m)')"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## This parameter contains tuphe hydraulic friction. For Manning roughness = 0.1 -> 530\n",
"const F1 = 130.0 # kg/m^(-8/3) <---\n",
"## Lake input \n",
"const Q_in = 50.0 # <---\n",
"\n",
"# these are fairly \"fixed\" constants\n",
"const L = 3.35e5\n",
"const ρ_w = 1000.0\n",
"const ρ_i = 910.0\n",
"const g = 9.81\n",
"const n = 3\n",
"const A = 5e-25\n",
"\n",
"## Geometry\n",
"ice_thick = 1000.0\n",
"Δx = 50e3\n",
"Δz_b = 0.0\n",
"Δice_thick = ice_thick - 0 # zero ice thickness at outlet <---\n",
"const Sₗ = 1000.0^2# lake area (constant with elevation) <---\n",
"\n",
"## Model\n",
"# gradient of hyd pot\n",
"const Ψ = (ρ_w*g*Δz_b + ρ_i * g * Δice_thick) / Δx\n",
"const day = 3600*24.0 # day in seconds\n",
"\n",
"\"Lake level as function of N (assuming z_b=0)\"\n",
"z_l(N) = (ρ_i*g*ice_thick - N)/(ρ_w*g)\n",
"\"Effective pressure at lake in terms of lake level\"\n",
"N(z_l) = ρ_i*g*ice_thick - z_l*ρ_w*g \n",
"\n",
"# define ODE\n",
"\"This function returns [dQ/dt, dN/dt] for given Q, N\"\n",
"function ode(u)\n",
" Q,N = u\n",
" dQdt = 4*Ψ/ (3ρ_i*L) * (Ψ/F1)^(3/8) * abs(Q)^(5/4) - 8A/(3n^n) * Q * abs(N)^(n-1)*N\n",
" dNdt = ρ_w * g * (Q-Q_in)/Sₗ\n",
" return [dQdt, dNdt]\n",
"end\n",
"Q0,N0 = \n",
"\n",
"## Integration time (adjust such that all the flood is captured) <---\n",
"#tspan = (0.0, 500day) # going to 500day will show cycles\n",
"tspan = (0.0, 60day) # short time span will show one outburst clearly\n",
"dt = 3600\n",
"ts = tspan[1]:dt:tspan[end]\n",
"\n",
"# solve it\n",
"sol = zeros(2, length(ts))\n",
"sol[:,1] = [1, N(ice_thick*0.95)] # inital discharge and lake level <--\n",
"for i = 2:length(ts)\n",
" # using the explicit midpoint method (2nd order)\n",
" # https://en.wikipedia.org/wiki/Midpoint_method\n",
" sold = sol[:,i-1]\n",
" sol[:,i] = sold .+ dt*ode(sold + dt/2*ode(sold))\n",
"end \n",
"\n",
"# plot solution\n",
"using PyPlot\n",
"subplot(2,1,1)\n",
"plot(ts/day, sol[1,:]);ylabel(\"Q (m^3/s)\")\n",
"subplot(2,1,2)\n",
"plot(ts/day, z_l.(sol[2,:]) ); xlabel(\"time (d)\"); ylabel(\"lake level (m)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: the water level can go negative because the hyd. pot. gradient is constant, i.e. not dependent on $N$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.9.3",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
[deps]
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment