Skip to content

Instantly share code, notes, and snippets.

@mauro3
Last active November 12, 2025 09:59
Show Gist options
  • Select an option

  • Save mauro3/04163405f15196a58f38d758bc0009fd to your computer and use it in GitHub Desktop.

Select an option

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",
"**To find out:**\n",
"- how can the frequency of events be increased?\n",
"- how can the max discharge be increased?\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": 34,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: redefinition of constant Main.Q_in. This may fail, cause incorrect answers, or produce other errors.\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG4CAYAAAC+ZBgrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACILElEQVR4nO3deVxU5f4H8M8MMMOOLAIii7gvICrmlpZmqeRSmaa5ZGX9Msv0mi3mvdfydrPbvde0TG27WWlqi5qVqVjmkqmJoCjugoCyyL7PMDPn98cwhx2BWc4M83m/XrxeMnOY+XLwPPM9z/N9nkcmCIIAIiIiIhLJpQ6AiIiIyNowQSIiIiKqgwkSERERUR1MkIiIiIjqYIJEREREVAcTJCIiIqI6mCARERER1cEEiYiIiKgOJkhEREREdTBBIiIiIqqDCRIRERFRHXadIDk6OqJfv37o168fnnrqKanDISIiIishs+fNav38/JCTkyN1GERERGRl7LoHiYiIiKghNpsgHTp0CBMnTkRQUBBkMhl27txZ75h169YhPDwczs7OiI6OxuHDh2s9X1RUhOjoaAwfPhwHDx60UORERERk7RylDqC1SktLERUVhSeeeAIPP/xwvee3bduGRYsWYd26dbjzzjvx4YcfIiYmBklJSQgNDQUApKSkICgoCGfPnsX48eORmJgIT0/PBt9PpVJBpVKJ3+t0OuTl5cHX1xcymcw8vyQRERGZlCAIKC4uRlBQEOTyJvqJhDYAgLBjx45ajw0aNEiYN29ercd69uwpvPrqqw2+xrhx44Q///yz0fdYvny5AIBf/OIXv/jFL361ga+0tLQmcwub7UFqilqtRlxcHF599dVaj48ZMwZHjx4FAOTn58PV1RVKpRLp6elISkpC586dG33NpUuXYvHixeL3hYWFCA0NRVpaWqO9TkRERGRdioqKEBISAg8PjyaPa5MJUk5ODrRaLQICAmo9HhAQgMzMTADA+fPn8cwzz0Aul0Mmk2HNmjXw8fFp9DWVSiWUSmW9xz09PZkgERER2Zjblce0yQTJoO4vLwiC+NiwYcOQmJgoRVhERERk5Wx2FltT/Pz84ODgIPYWGWRnZ9frVSIiIiKqq00mSAqFAtHR0YiNja31eGxsLIYNGyZRVERERGQrbHaIraSkBFeuXBG/T05ORkJCAnx8fBAaGorFixdj9uzZGDhwIIYOHYqPPvoIqampmDdvnoRRE5El3SpWoaJSixAfV6lDIbJKl7OKEejlDA9nJ6lDsTo2myCdPHkSo0aNEr83zDCbM2cONm7ciGnTpiE3NxcrVqxARkYGIiIisHv3boSFhUkVMhFZUJlag5g1h1FUXonfXhqJoHYuUodEZFWOXsnBjE+Oo39oO+yYf6fU4Vgdm02QRo4cCeE228jNnz8f8+fPt1BERGRN9p/PRk6JfnHX02kFTJCI6vjkSDIAID61ACqNFkpHB4kjsi5tsgaJiGhXwk3x3xmFFRJGQmR9ytQa/HE1V/w+q1DVxNH2iQkSNUmrE1Cu1kodBlGLVGp1OHLllvh9QZlawmiIrM+J5DyUV1a37QXlvEbqYoJEjSpVaTD+vcMY8I9YnE4rkDocomY7d7MIFZU68ftilUbCaIisz8mU/FrfF1fwGqmLCRI16n9HknEhsxjllVp8/keK1OEQNVvdhJ6NP1Ftp9MLan3Pa6Q+JkjUoIIyNT46dE38PiG1QLpgiFooOacUAOCq0BedlrDxJ6rFcI24VV0jxRWVUoZjlZggUYM+OZyMYpUGfu4KAEBybinUGt1tforIOlzP1Tf+ER29AADFKjb+RAZqjQ43C8oBAH0M1whvIuphgkT15Jeq8dnv+umfbz4YARcnBwgCcKPqgiKydtdzywAAEUFs/InqSs8vg07Q97B2ae8GgNdIQ5ggUT0fHLiCUrUWfYI8MbZPIEJ89OvHpOWVSRwZ0e1pdQLS8vX/V3sHeQLQTzggIj3DDUSojyvcFPrlEMsqeY3UxQSJarmcVYyNR1MAAC+N7QGZTIYQb/02DalMkMgG3CwoR6VWgMJBLt4dc6kKomopVUPQYb6uYp1emYrXSF1MkEgkCAJe/+EcNDoB9/UOwMge/gAgrkCcVcTF9sj6Ge6Og31c4OGsvzsuZYJEJDJcI5183eBi6EHiNVKPzW41Qqb389lM/H4lFwpHOf42vrf4eICnEgATJLIN1/P0d8edfN3gWtX4sweJqJphEkOoryu0Ov2WXWVqDrHVxR4kAqCf1r981zkAwLy7OiPUt3r3c39PZwBAVhGXoifrZ7g7rjl8oNbqUKnlLEwiALieV6MHyalqiI03EfUwQSIAwIofknCrWIUu7d0wf1TXWs/5e7AHiWyHeHfs4woXRfXmm/wAINJPYkjP089IDvVxhZvSMMTGHqS6mCARfjmfhe3xNyCXAf+eGgVnp9o7OgdU9SBlF7MHiaxfalXjH+brCoWDHI5yGQB+ABAB+htdtVYHR7kMQe1cxJsI3kDUxwTJzmUXV+CV7xIBAE+N6IwBod71jjEkSHmlai4WSVZNEARxOYpQHzfIZDJ+ABDVIE5i8HaBg1wmTvNnnV59TJDsmFYnYNHWBOSUqNAjwAOL7+ve4HHerk5wctDfhd8qYS8SWa+8UjVKqtY8CvbWz77kBwBRNcMNRIiPvs7UUKdXyh7Wepgg2bE1v1zG0au5cFU44IOZA+oNrRnIZDL4exgKtVmHRNbLsFZXoKez+P9Z/ADgYpFE4jUSWpUgsYe1cUyQ7NSu0zfx3i+XAei3E+nq797k8Yap/tlMkMiK1W38AcBVWfUBUMkPAKK614hbjXWQBEGQLC5rxATJDsVdz8OSb04DAOYOD8fkAcG3/RlDHVJmIRMksl5i/VGNZSpcnao+ALhSMJGYIIX51u5B0uoEqLkURi12nSA99NBD8Pb2xpQpU6QOxWLO3ijEkxtPQq3R4b7eAXjt/l7N+jnDVH/OZCNrVnOPKYPqIQQOsRE1VoMEsE6vLrtOkF544QV88cUXUodhMWdvFGLmJ8dRWF6JAaHtsGZ6PzhUTYG+HX9O9Scb0NAQm5uSNRZEAFCi0iC3VA2gOkFycpBD4aBPBbglT212nSCNGjUKHh4eUodhEX9czcWMj4+hsLwS/UPb4fMnB4nbMDRHe/YgkQ1oaIjNxYl7TREBQGpVD6u3qxM8nZ3Exw11euXsZa3FZhOkQ4cOYeLEiQgKCoJMJsPOnTvrHbNu3TqEh4fD2dkZ0dHROHz4sOUDtQI74tPx2P+Oo6hCg4Fh3vj8yUHwqHFxNIc4xMYibbJSKo0WGVX/PxvqQWLjT/ZO7GH1dav1uKuTYaYnbyJqstkEqbS0FFFRUVi7dm2Dz2/btg2LFi3CsmXLEB8fjxEjRiAmJgapqakWjlQ6FZVaLP/+LP6y7TQqtQLGR3bApqcG17pzaC7DNP9b7EEiK3UjvxyCoK+p8HVTiI+7iOu8sPEn+5bWwBA0wKn+jWn+GIuViYmJQUxMTKPPr1q1CnPnzsVTTz0FAFi9ejX27t2L9evXY+XKlS1+P5VKBZWqOjkoKipqedAWdCW7GC9sSUBShj7O+SO7YMmYHpA3s+aoLv+qaf65pWpUanVwcrDZ3JraqOs1Gn+ZrPr/ec1pzET2rLpGz6XW49yPrWFt8lNOrVYjLi4OY8aMqfX4mDFjcPTo0Va95sqVK+Hl5SV+hYSEmCJUk1NptHg39hJi1hxGUkYRfNwU+OyJO/DyuJ6tTo4AwMdVIe5plcPVtMkKpTYwgw2onqXDxp/sXUqNjZxrcnFiD1JD2mSClJOTA61Wi4CAgFqPBwQEIDMzU/x+7NixmDp1Knbv3o3g4GD8+eefjb7m0qVLUVhYKH6lpaWZLf7WEAQB+85lImbNYaz55TIqtQLu6emP3S+MwKge/ka/vlwug5+7oQ6JCRJZn2u3SgAA4e1r11dw+IBI79otfYLUuX3thYENPUic5l+bzQ6xNUfNbnZAn0TUfGzv3r3Nfi2lUgmlUmmy2Ezp+LVcvLP3IuKu5wMA/NyVeH1Sb4yP7FDvHBjD31OJzKIKzmQjq3QtR9/4d/Gr0/grOHxAVFGpxc3CcgBAuF/DNxHcj622Npkg+fn5wcHBoVZvEQBkZ2fX61WyVVqdgNikTHx06BpOpRYAAJyd5Jg7PBzP3N2lVYXYt6Mv1C7kfmxklQx3x+xBIqovJbcUggB4OjvWmsQAAG68RhrUJhMkhUKB6OhoxMbG4qGHHhIfj42NxQMPPCBhZMbLKVFhx6kb2Hz8OlKqai4UDnJMHRiMF0Z3E7cEMQdDoTZ7kMjaVFRqcaNAf3fcuc7dsaEHicMHZM9qDq/VHVlwZS9rg2w2QSopKcGVK1fE75OTk5GQkAAfHx+EhoZi8eLFmD17NgYOHIihQ4fio48+QmpqKubNmydh1K2j0epw+EoOvv4zDbFJWdDo9BsKerk4YfaQMDw2LEychm9OhrWQbhWzB4msS3LV8JqXixN86twdc/iAqPoa6VynhxVgL2tjbDZBOnnyJEaNGiV+v3jxYgDAnDlzsHHjRkybNg25ublYsWIFMjIyEBERgd27dyMsLEyqkFtEo9Xh2LU8/JSYgX3nMsXl4QEgKtgLj9wRggf7dRSL6yzBkISxSJusjTi85udW7+64eqFINv5kv65WTWKo28MK1Bhi40KRtdhsgjRy5EgIgtDkMfPnz8f8+fMtFJHxtDoBR6/mYHdiBvaey0JejaTI29UJD/bviGl3hKBnoKck8XHDWrJWyTlVjX8Dd8euVVuNcJVgsmeNzWADABfDEFslr5GabDZBaosEQcDCrQliYuTjpsDYPgG4P7IDhnT2lXxxxuoaJA6xkXUxNP5dGmz8q3qQKrXQ6QSj1gMjskWCIIjLYDR0E2HoQeJ2PLUxQbIijg5yTLsjBIXllbg/ogOGdPaBoxWtWG0YYsspUUOrE+DADxqyEldzqofY6jIMsQH6JMmSw9JE1iCvVI2iCg1kMqCTb+M1SOxlrY0thZV5ZVxPqUNolJ+7AjKZfigwr1SN9h7WuS4U2Zfb3R07O1YnSGVqJkhkfwxrhAV5ucDZyaHe81wrrGHW0z1BVs/RQS6un8FhNrIWuaVqFDdxdyyXy7jdCNm1pm4ggJrb8bAHqSYmSNQi4kw2FmqTlTDUHzV2dwzUXOeFHwBkf8RV5huo0QMAVyWvj4YwQaIWMRRq3+JUf7ISl7OLAQBd/Btu/AFuWEv27UqWvgepSyM9SG68PhrEBIlaxDDVn9uNkLW4XNX49whoToLEO2SyP5eqbiK6BXg0+Hz1Yqq8PmpigkQtwiE2sjaXsppu/IHqBImzdMjelKk1SMvTb8PTvZFrxFCkrdbooNHqLBabtWvxdA5BEHDw4EEcPnwYKSkpKCsrQ/v27dG/f3/ce++9CAkJMUecZCW4FhJZm0tVPUiNNf4AxJlr5ZUcQiD7ciVbf334uSvqbcNjYOhBAvSLRXpa0fIyUmr2WSgvL8dbb72FkJAQxMTE4KeffkJBQQEcHBxw5coVLF++HOHh4bj//vtx7Ngxc8ZMEuJq2mRN8kvVyCnR/1/s1owaJPYgkb25mKnvYW3qBkLpKBfXteN2I9Wa3YPUvXt3DB48GBs2bMDYsWPh5ORU75jr16/jq6++wrRp0/DXv/4VTz/9tEmDJem1535sZEUMw2sd27k0ub4Rdysne3U5+/Y9rDKZfimM4goNr5Eamp0g/fzzz4iIiGjymLCwMCxduhQvvvgirl+/bnRwZH0MPUi3ilUQBKHexqBElnRJbPwb7z0CWKRN9qu6Ru/214g+QeI1YtDsIbbbJUc1KRQKdOvWrVUBkXUzrJ6t1upQWF4pcTRk7y5n3X74AGCCRPbrcjNq9ICaq2nzGjFoVSXWnj17cOTIEfH7Dz74AP369cOMGTOQn59vsuDI+jg7OcDLRT+8yjokktqlZidI+sa/VMXhA7IfxRWVuFFQNYPNv+lrpHqqP68Rg1YlSC+99BKKiooAAImJiXjxxRdx//3349q1a1i8eLFJAyTrE2CYycY6JJJYs++OlYbdynl3TPbDUH/k76GEl2v9uuGaxB4kFmmLWrVrY3JyMnr37g0A+O677zBhwgS89dZbOHXqFO6//36TBkjWx9/DGZeySjjVnySVU6JCbqkaMhnQtYkZbADgYuhB4t0x2ZHmDkEDgKuSq2nX1aoeJIVCgbKyMgDA/v37MWbMGACAj4+P2LNEbRen+pM1MAyvhXi71lrHpSFurEEiO9ScNcIMWKdXX6t6kIYPH47FixfjzjvvxIkTJ7Bt2zYAwKVLlxAcHGzSAMn6tPfkdiMkvUvi+i5N9x4B3KyW7FN1jR6vkdZoVQ/S2rVr4ejoiG+//Rbr169Hx44dAeiXAhg3bpxJAyTrw+1GyBqcz9A3/r06eN722OqFIjl8QPZBEAQk3dSP6DTnGuGGtfW1qAdp3759GDVqFEJDQ/Hjjz/We/7dd981WWBkvcS1kFikTRJKytA3/r2b0/gbirQreXdM9uFWsb5GTy4DegTefohNrNNjkbaoRT1I8+bNQ/v27TFt2jRs2bIFhYWF5orL7IqLi3HHHXegX79+iIyMxMcffyx1SDajugaJQ2wkjUqtDherhg96B90+QXJxYuNP9uVc1Q1El/bucHZqukYPqO5B4n6F1VqUIF27dg2HDh1CZGQk3n33XQQEBGD06NF47733kJKSYqYQzcPV1RUHDx5EQkICjh8/jpUrVyI3N1fqsGyCvyeH2Eha126VQq3RwV3piBBv19seXz3Nn40/2YeWDK8BNdZB4k2EqMU1SH379sVf//pXnDhxAteuXcPUqVOxZ88e9OrVC1FRUfj73/+OkydPmiNWk3JwcICrq75hraiogFarhSAIEkdlGww9SGVqLUpY00ESSMrQ91736uABufz2292IBaiVWuh0vM6p7ROHoJvRwwpA3MuQNUjVWlWkbRAUFIR58+Zh9+7dyMnJwd///nekpKRg3LhxeOutt4wK7NChQ5g4cSKCgoIgk8mwc+fOesesW7cO4eHhcHZ2RnR0NA4fPtyi9ygoKEBUVBSCg4Px8ssvw8/Pz6iY7YWb0hHuVRdTNmeykQQMBdrNqT8Cqou0BQGo0PAOmdq+8zebX6MHcJp/Q1o1zb8hbm5uePjhh/Hwww9Dp9MZPVxVWlqKqKgoPPHEE3j44YfrPb9t2zYsWrQI69atw5133okPP/wQMTExSEpKQmhoKAAgOjoaKlX9YaB9+/YhKCgI7dq1w+nTp5GVlYXJkydjypQpCAgIMCpue+HvoUSJSoPsYhU6t7/9FFIiU2rx8EGNGowytVbsUSJqi8rUGiTnlgJo/jUibsfDBEnU4lYiNzcXZ86cQVRUFHx8fJCTk4NPP/0UKpUKU6dORa9evSCXy9G+fXujAouJiUFMTEyjz69atQpz587FU089BQBYvXo19u7di/Xr12PlypUAgLi4uGa9V0BAAPr27YtDhw5h6tSpDR6jUqlqJVv2viBmew8lruWUsg6JLE4QhBYPH8jlMrg4OaC8UqvfSoE5PbVhFzKLIQj6G1nDBuO3IxZpc4hN1KIhthMnTqBLly4YPXo0unbtiri4OAwaNAiffvopvvzyS0RHRzc7KTGGWq1GXFycuIK3wZgxY3D06NFmvUZWVpaY5BQVFeHQoUPo0aNHo8evXLkSXl5e4ldISEjrf4E2QCzU5hAbWVhWkQp5pWo4yGXNWiHYwFCoze1GqK0z9LA29wYCYJF2Q1qUIC1btgxTp05FYWEhXnvtNTz44IMYPXo0Ll26hMuXL2PGjBl48803zRWrKCcnB1qttt5wWEBAADIzM5v1Gunp6bjrrrsQFRWF4cOH4/nnn0ffvn0bPX7p0qUoLCwUv9LS0oz6HWyduBYSe5DIwgwF2l3auzVr+rIBVwome2HoYW3u8BrAIu2GtGiILS4uDu+99x48PDywcOFCvPLKK3j66afF55977jlMnDjR5EE2RiarPXtFEIR6jzUmOjoaCQkJzX4vpVIJpbJ5XZX2wJAgcbsRsrSkFhafGrhypWCyE625RlikXV+LepDUajVcXFwAAE5OTnB1da0188vX19ciawn5+fnBwcGhXm9RdnY2i6wtxN+TG9aSNBJv6HuQWjJ8APADgOyDRqvDhcyWD7EZelhVGh20XAoDQAsTpJCQEFy7dk38fuvWrejQoYP4fUZGhkWmyisUCkRHRyM2NrbW47GxsRg2bJjZ35+4HxtJJzFdnyD1DW7Xop/jEALZg0tZJaio1MFD6YhwX7dm/5zhBgLgNWLQoiG26dOnIzs7W/x+/PjxtZ7ftWsXBg0aZJLASkpKcOXKFfH75ORkJCQkwMfHB6GhoVi8eDFmz56NgQMHYujQofjoo4+QmpqKefPmmeT9qWnidiMcYiMLulWsws3CCshkQERHrxb9rGGqP4tQqS07k14AAIgM9mrWIqoGSkc5HOQyaHUCytRaeDg7mSlC29GiBGn58uVNPr9s2TI4ODS/aLIpJ0+exKhRo8TvFy9eDACYM2cONm7ciGnTpiE3NxcrVqxARkYGIiIisHv3boSFhZnk/alphllsRRUaVFRqW1QsS9Rahsa/S3t3cbHS5jL0IJVziI3asNOt7GGVyWRwdXJAsUqDUu6QAMCEC0UCELfuMIWRI0feduuP+fPnY/78+SZ7T2o+T2dHKB3lUGl0uFWsQoiP6f72RI2pbvxb1nsE1JjGzOEDasMMNxFRrbhGXJX6BIl1enqt2mqk5jAb2SeZTFajUJvDbGQZiWLj367FP1u9EB4bf2qbKiq1uJip34anb0i7Fv88l8KorcUJUnJyMoYPH26OWMjGiIXaRSzUJvMTBAFnjOhBqt5KgT1I1DYlZRRBoxPg66ZAkJdzi3/elb2stbQoQTp79ixGjBiBxx9/3EzhkC0RC7U5k40s4EZBOXJL1XCUy1q0AJ6BOM2fRdrURiXWuIFo7pqANbkpWKdXU7MTpKNHj+Kuu+7CnDlz8Nprr5kzJrIRXCySLMnQe9Qj0KNVkwJclRw+oLbtdNUQdEsLtA2qtxthDxLQggRpzJgxmD17Nv75z3+aMx6yIYFe+kVDM5kgkQWcaeXsHANXJw4fUNtmuEaiQlo+BA1U71dYXsmbCKAFCZKbmxsyMjJuO7OM7EeHqjHuzEImSGR+xszOAaobf/YgUVtUotLg6q0SAEb0IDlV1elxGBpACxKkI0eO4OTJk3jiiSfMGQ/ZkEAmSGQhGq0Op9MKAABRrZidA3CGDrVtCakFEASgYzsX+Lm3bt/Q6psI9rICLUiQunXrhiNHjiAuLg7PPfecOWMiG2HoQcoorGDPIpnVxaxilKq18FA6onuAR6teg5vVUlsWdz0fADCwk3erX4M3EbW1aBZbUFAQDh06hPj4eHPFQzYkoGo17fJKLYrK+aFD5mNo/PuHecOhBdsn1MTGn9qyuFT9NRIdZkyCxJuImlq8DpK3tzd++eUXc8RCNsbZyQE+bgoAQEZRucTRUFtmSJCiQ1vf+IvDB5yhQ22MVicg/ropEyTeRACtXEnbxcXF1HGQjQr0rB5mIzIXUwwfGKYwl1VqOSRMbcrl7GIUqzRwUzigRyuHoIHq/QpZpK1n9F5sJSUl0Ol0tR7z9Gz5Im5kmzp4OSMpo4iF2mQ2WUUVSM8vh1zW+gJtoHoRPEEAKip1YsJEZOsMNxD9QtvB0aFV/R4AOMRWV6vOZHJyMsaPHw83Nzd4eXnB29sb3t7eaNeuHby9W3+HR7bHMJMto4BDbGQehsa/Z6An3JWtv6dzqbG4JNdCorYkLsX4IWiAdXp1taq1mTlzJgDgf//7HwICAlq1pDm1DTVnshGZgymG1wBALpfBxckB5ZVa/XYj7qaIjkh6YoF2Jx+jXoc9SLW1KkE6c+YM4uLi0KNHD1PHQzaGq2mTuZ00QfGpgZuyKkGq5AcAtQ23ilW4nlsGmQzoZ8QQNFBjs1rWIAFo5RDbHXfcgbS0NFPHQjaIPUhkTuVqLc7d0G+fMMDI4QOg5l5T/ACgtiHueh4AoLu/B7xcnIx6LUORNrca0WtVD9Inn3yCefPm4caNG4iIiICTU+0/St++fU0SHFk/rqZN5nTyeh40OgEd27kg2Nv42bPcrZzammPX9AnSoHDjhteA6jo9blar16oE6datW7h69WqtbUdkMhkEQYBMJoNWy8bHXhh6kEpUGhRXVMLD2bg7GKKajl3LBQAM7uxjklpHsQeJNRbURhiukaFdfI1+LUMPkkqjg1YntHpR1raiVQnSk08+if79+2PLli0s0rZzrgpHeLk4obC8EpmFFUyQyKT+uFrV+Hc2vvEHqnuQWIRKbUFuiQoXMosBAINN0IPkWmPpizK1xu7b81YlSNevX8euXbvQtWtXU8dDNqiDlzMKyyuRUViBbkYsUkZUU6lKgzPp+vqjISZKkLhSMLUlx5P1w2s9Ajzg28oNamtSOsohlwE6QX+N2HuC1Koi7XvuuQenT582dSwWdfHiRfTr10/8cnFxwc6dO6UOyyaxDonM4eT1fGh0AoK9XRDi42qS1xQTJBZpUxtgyuE1QF8q48a1kESt6kGaOHEi/vKXvyAxMRGRkZH1irQnTZpkkuDMqUePHkhISACgXw28U6dOuO+++6QNykZxJhuZg2F4zVS9RwDgqmTjT21H9TVi/PCagYvCAcUqDQu10coEad68eQCAFStW1HvOFou0d+3ahdGjR8PNzU3qUGxSoKdhLSSupk2mY7g7NmWC5MaF8KiNyClR4XJ2CQBgcLgJrxGlI1Cs4lR/tHKITafTNfplquTo0KFDmDhxIoKCgiCTyRoc/lq3bh3Cw8Ph7OyM6OhoHD58uFXv9fXXX2PatGlGRmy/2INEplai0iDxhqH+yJR3x1WbcTJBIhtnuIHoGegBbzeFyV63erFIXiOt39XOzEpLSxEVFYW1a9c2+Py2bduwaNEiLFu2DPHx8RgxYgRiYmKQmpoqHhMdHY2IiIh6Xzdv3hSPKSoqwu+//47777/f7L9TW8UaJDK1Y1dzodUJCPVxRbC3aeqPgJo9SLw7Jtv2+xXT1h8ZcCJDtWYPsW3duhXTp09v1rFpaWlITU3FnXfe2erAYmJiEBMT0+jzq1atwty5c/HUU08BAFavXo29e/di/fr1WLlyJQAgLi7utu/z/fffY+zYsXB2dm7yOJVKBZVKJX5fVFTUnF/DLhh6kG5yw1oykUOXbwEA7uruZ9LXZZE2tQWCIODQJcM10t6kr80Na6s1uwdp/fr16NmzJ/71r3/h/Pnz9Z4vLCzE7t27MWPGDERHRyMvL8+kgdakVqsRFxeHMWPG1Hp8zJgxOHr0aIteq7nDaytXroSXl5f4FRIS0qL3acsMPUhFFSzsI9M4aGj8u5mn8ecQG9myq7dKcKOgHApHOYaYsP4I4Ia1NTU7QTp48CD+85//4Ndff0VERAQ8PT3RrVs3REZGIjg4GL6+vpg7dy46deqEs2fPYuLEiWYLOicnB1qtFgEBAbUeDwgIQGZmZrNfp7CwECdOnMDYsWNve+zSpUtRWFgofnEvumoezk5wr5odxDokMlZKTimu55bBUS4z+fCBm1Lf+HOrEbJlBy/lANAvDulSY3FHUxBvItjL2rJZbBMmTMCECROQm5uLI0eOICUlBeXl5fDz80P//v3Rv39/yOWWK2uqu4K3YauT5vLy8kJWVlazjlUqlVAqjV+Iq60KaueMS1kluFlQjq7+7lKHQzbMMLwWHeZt8oXqqou02fiT7TL0sN5t4uE1oOZNBHuQWjXN39fXFw888ICpY2k2Pz8/ODg41Ostys7OrterRJYR1M5FTJCIjGGu2gqgukibjT/ZqopKLY5XzWAzxzVSvV8hbyKsdhZbUxQKBaKjoxEbG1vr8djYWAwbNkyiqOxbx3b6tZBuMEEiI6g1OhytWvzOHHfHbPzJ1h1PzoNKo0OgpzO6maG3nitpV2tVD5IllJSU4MqVK+L3ycnJSEhIgI+PD0JDQ7F48WLMnj0bAwcOxNChQ/HRRx8hNTVVXMSSLKujNxMkMt7J63koU2vh565A7w6eJn99sfHnZAKyUQcvVg+vmWOjeBZpV7PaBOnkyZMYNWqU+P3ixYsBAHPmzMHGjRsxbdo05ObmYsWKFcjIyEBERAR2796NsLAwqUK2a2IPUj4TJGq9/UnZAIC7u/tDLjdD419VX1FWqW1xzSKR1ARBwP7z+rrZkT1M38MKsEi7JqtNkEaOHAlBEJo8Zv78+Zg/f76FIqKmcIiNjCUIAmLP6+sK7+ttnlpCQ+MvCEBFpc7kM4CIzOlSVglS88qgcJSbpf4IqFGkXckeJKNqkHJycrhgIgGoHmLLLKyAVtd0YkvUkAuZxUjLK4fSUW7yBSINXJyqEyIOIZCtiU3S30AM7+qn3zPNDAzXCHuQWpEgFRQU4LnnnoOfnx8CAgLg7e2NwMBALF26FGVlZeaIkWyAv4czHOUyaHQCsou5FhK1XGySfuhgRDc/safH1BzkMvEDgEWoZGv2VV0j5uphBSAmXryBaOEQW15eHoYOHYobN25g5syZ6NWrFwRBwPnz5/H+++8jNjYWR44cwenTp3H8+HG88MIL5oqbrIyDXIZAL2ek55fjRn45Oni5SB0S2ZhYCzT+gL4ItbxSy9W0yaZkFlbgTHohZDJgdC9/s70P92Kr1qIEacWKFVAoFLh69Wq99YZWrFiBMWPGYPbs2di3bx/ee+89kwZK1q9jOxd9glRQjoFSB0M25WZBORJvGBp/MydISgfklvIDgGxLbFVxdv+QdvD3aHrvUGNwL7ZqLUqQdu7ciQ8//LDBxRgDAwPxzjvv4P7778fy5csxZ84ckwVJtoGF2tRaht6j6FBv+Lmbd8V6VyfDVH9+AJDt2HdOX380pk+gWd+H0/yrtagGKSMjA3369Gn0+YiICMjlcixfvtzowMj2iGshcao/tdCPZ24CAMZFmLfxB6qn+nOIjWxFbolKXEB1rIUSpIpKnd1PuGlRguTn54eUlJRGn09OToa/v/nGRsm6GXqQuN0ItcTNgnL8mZIPmQyY0DfI7O9nWCySG9aSrdh9NhNanYDIjl4I93Mz63vVnB1n771ILUqQxo0bh2XLlkGtVtd7TqVS4W9/+xvGjRtnsuDItgRxiI1awdB7dEcnHwR6ma+2wqB6uxH7bvzJdvxwWn+NTIzqYPb3UjrKYVij1d5vIlpUg/TGG29g4MCB6NatG5577jn07NkTAJCUlIR169ZBpVLhiy++MEugZP1qDrFxlWJqrh9OZwAAJkaZv/cIqN6wljVIZAsyCsvxZ0oeAMv0sMpkMrgqHFGi0qBEpYE9jwm1KEEKDg7GH3/8gfnz52Pp0qXiStcymQz33Xcf1q5di9DQULMEStbPMMRWqtaiqFwDL1cniSMia5ecU4rEG4VwkMtwvwXqjwDAtWoIgT1IZAt+OpMBQQDu6OQt9tKbm6vCASUqjd3PZGvxamzh4eH4+eefkZ+fj8uXLwMAunbtCh8fH5MHR7bF2ckBvm4K5JaqkV5QBi9XL6lDIitnGDq4s6sffM08e83A3ZAgccNasgG7xOE1y/SwAvprJLtYZffXSKuXq/X29sagQYNMGQu1AR29XZBbqsaN/HL0CWKCRI3T6QR8G5cOAHjAgo2/q1iDZN93x2T9LmYW40x6IRzlMtwfaf76IwM39rICMHIvNqK6OJONmut4ch5S88rgrnS0aOPPHiSyFV+fTAMA3NPT3+zrg9Uk3kTYeZ0eEyQyKc5ko+b6pqrxnxgVJM4sswQ3JkhkA9QaHXbE3wAATLsjxKLvzZsIPSZIZFJcTZuao6iiErvP6mevPTIw2KLvbUiQSuy88Sfr9sv5LOSVquHvocTd3dtb9L15jegxQSKTCq6a6p/O1bSpCT+cvomKSh26B7ijX0g7i763GzfjJBtgGF57ODoYjg6W/aiu7mW172uECRKZVKivKwAgNa9M4kjIWgmCgE3HUgEAjwwMsfh6Wbw7JmuXlleG3y7dAgBMjbZsDytQ8ybCvq8RJkhkUiHe+gSpoKwSRRWVEkdD1uhEch7OZxTB2UmOKRI0/qyvIGv3xR8pEARgRDc/dG7vbvH3502EHhMkMik3pSN83RQA9HdBRHVtPJoCAJg8IBjtXBUWf38OH5A1K1VpsPVP/fDaE3d2kiQG3kTo2UWC9NBDD8Hb2xtTpkxp0XPUOiE++l4kJkhU142Ccuw9lwkAeHxYJ0licKuxF5thNwAia7Ej/gaKKzQI83XFyO7SbPThquRaYYCdJEgvvPBCo3vENfUctU51gsRCbartiz9SoBOAO7v6onuAhyQxGHqQBAEor7TvDwCyLjqdIPawzhnaCXK5NPtZsgdJzy4SpFGjRsHDo+HGuKnnqHVCffQz2VioTTUVllVic1Vx9hPDwiWLw1XhAENduL3XWJB12X8+C1eyS+CudMQUCy9/UZObggkSYAUJ0qFDhzBx4kQEBQVBJpNh586d9Y5Zt24dwsPD4ezsjOjoaBw+fNjygVKzhfpwJhvV9/kfKShRadAjwAP39JRuj3CZTFbjA4A9SGQdBEHABweuAABmDw2Dp7N0m32zSFtP8gSptLQUUVFRWLt2bYPPb9u2DYsWLcKyZcsQHx+PESNGICYmBqmpqeIx0dHRiIiIqPd18+ZNS/0aVINhJltaPhMk0itVafC/35MBAM/d01WyoQOD6q0U7PsDgKzH4cs5OJ1eCGcnOeYOl66HFQDclFwrDDBis1pTiYmJQUxMTKPPr1q1CnPnzsVTTz0FAFi9ejX27t2L9evXY+XKlQCAuLg4s8epUqmgUqnE74uKisz+nrbKUIOUnlcOnU6Q/MOQpLf5+HUUlFUi3M8N4y2471pjuFs5WZu1Vb1Hjw4Ktei+aw1hD5Ke5D1ITVGr1YiLi8OYMWNqPT5mzBgcPXrUorGsXLkSXl5e4ldIiGX3xrElHbyc4SiXQa3VIau4QupwSGLFFZXYcPAaAODZu7vAwQoSZu5WTtbk0KVbOJGcB4WDHP93V2epw6lVpG3PMz2tOkHKycmBVqtFQEBArccDAgKQmZnZ7NcZO3Yspk6dit27dyM4OBh//vlns56raenSpSgsLBS/0tLSWvdL2QFHB7m4aS1nstGHB68hr1SNzu3dMHlAR6nDAVA9hMAaJJKaTifg7Z8vAABmDQlDBy8XiSOqHoLWCYBKo5M4GulIPsTWHHW3IhAEoUXbE+zdu7dVz9WkVCqhVErb7WlLQn1ckZpXhtS8MgwK95E6HJJIVlEFPjmi7z16eWxPi+8p1RjO0iFr8f3pG0jKKIKH0hHP39NV6nAAVF8fgH6YzdnJQcJopGMdrVUj/Pz84ODgUK+3KDs7u16vElmXEM5kIwDvxl5CRaUO0WHeGNvHeq5Z1liQNaio1OI/ey8BAJ4d1QU+bpZfWb4hcrmMExlg5QmSQqFAdHQ0YmNjaz0eGxuLYcOGSRQVNUdI1VpI6UyQ7Nap1Hxsq9qRfGlMT4tvStsUbjdC1uCDA1dwo6AcHbycJV0brCG8RqxgiK2kpARXrlwRv09OTkZCQgJ8fHwQGhqKxYsXY/bs2Rg4cCCGDh2Kjz76CKmpqZg3b56EUdPtcC0k+6bR6rBsx1kIAvDwgGAM7GRdw6zuSu5WTtK6kl2CDQevAgCWT+wNF4V1DWO5KRxwC/Y9kUHyBOnkyZMYNWqU+P3ixYsBAHPmzMHGjRsxbdo05ObmYsWKFcjIyEBERAR2796NsLAwqUKmZjAkSFwLyT59/sd1nM8ogpeLE167v6fU4dTjquAQG0lHEAT8bedZVGoFjOrRHmP7BEodUj0chraCBGnkyJG3nUY4f/58zJ8/30IRkSkYFovMKlKholJrt0V+9ujqrRL8e69+Vs6rMT3hK/GaLg3hXlMkpU3HruOPa7lQOsrxxqQIqxp+NnDjNWLdNUhku9q5OsGj6gJLZy+S3ajU6vCXbQmoqNRheFc/TBtoneuFVd8d2299BUnj6q0S/HP3eQDAK+N6ItTXVeKIGuZWNeRXZsfXCBMkMguZTCbOZEvJYYJkL9775TLOpBfCy8UJ/5kaZbWrqLuxBokkoNbUvoF4fFgnqUNqFIfYmCCRGYX7uQEAUnJLJY6ELOGX81nidglvPhiBQC9niSNqHNdBIim88cM5nEkvhKezI/49ta/V3kAAHIYGmCCRGXXy0/cgJecwQWrrrt0qwaKtCRAEYPaQMEyMCpI6pCbx7pgsbeuJVGw+ngqZDFgzvb9VrJjdlOrteDjERmRynXzZg2QPcktUeOrzkyhWaXBHJ2/8bUJvqUO6LXeu8UIWdPRKDv7+/TkAwIv3dceonv4SR3R7blwokgkSmY84xMYapDarRKXB45/9iWs5pQjycsYHMwdA4Wj9zYq4FxtrkMjMEtML8fQXJ6HW6jA+sgOeG2Ud24ncDmexMUEiM+pUlSDdLCxHRSXv1NuaEpUGczf+icQbhfBxU+DLpwbD38N6645qcuNu5WQBSTeLMOezEyhVazGsiy9WTYuyyin9DakeYmOCRGRyvm4KeCgdIQhcUbutyS9VY+bHx3A8OQ/uSkdsfOIOdGnvLnVYzWZo/HUCUFFpv7uVk/mcTMnDtI/+QF6pGn2DvfDRYwOhdLSd9eDEXlY7HoZmgkRmI5PJxF4kFmq3HddulWDKhqM4nV6Idq5O2PzUYPQNbid1WC3iWmPhUhZqk6n9dCYDsz89geIKDQaGeePLuYPFujdb4cbV5qVfSZvatk5+bki8UYgUJkhtwv6kLPxlWwKKVRoEejrji7mD0D3AQ+qwWkwul8FN4YBStbZqLSTrW+2bbI9Gq8N/Yy9h/W/6Pdbu7t4eG2ZFW90+a83Baf5MkMjMwqtWieVMNttWotLgnz+dx5YTqQCAgWHeWDdrgM3UHDXEVemIUrXWru+QyXSuZBfjxW/O4HRaAQDg/+7qjJfH9oCjg20O1BiGocvseJo/EyQyKw6x2TZBEPBTYgZW7r6AGwXlAIAn7uyEpTG9bGK2WlPclY64Vayy6xoLMl6JSoOPDl7FhkPXoNbo4OHsiH8+FIlJVr4W2O0YapDs+QaCCRKZVSdO9bdJOp2Ag5du4b1fLyM+tQAAEOztgn9PicLQLr7SBmci1UWo9vsBQK1XqtLgm5NpWHvgKnJKVACAkT3a4+3Jfa16FfnmqjvT01Zm35kSEyQyq/CqxSIziypQrtba5Fi8PcktUeGnxAx88cd1XMkuAQC4KhzwzF1d8PRd4XBVtJ0mQ9xuxI6nMVPLXcwsxvZT6dhyIhVFFfr/O2G+rnh1XE+MiwhsM4mEIUHS6ASotTqbmoFnKm2ntSOr5O2mgJeLEwrLK5GSW4peHTylDolqEAQBl7JK8PuVHBy8dAtHruRAq9OvC+SudMSjg0Lw9IjO8Pe0/TviurgQHjWHSqPFqesFOHo1B7FJWbiQWSw+18nXFXNHdMa0gSE2P+Rcl1uNm6FSlZYJEpE5dPJzw+m0AqTkMEGSilqjQ2ZhBW4WliM9vxwXM4twIbMYSTeLkFuqrnVs32AvPNivI6YODIaHs5NEEZtf9X5srEEifa1NRkE5bhZW4HpuKc5nFCEpoxgXMoqg0lSvleXkIMPIHv6YGh2M0b0C4GDFG84aw0Eug7OTHBWVOpSqNPBxU0gdksUxQSKzC/d1xem0Amw8moITKXkAgKYWL665srFQ6/Ea/67xTN3Xauxn0MjP1P/5Ro5rRix1D6z9M7f/vZr7M2ji/dUaHYorNCiu0KBEpUFxRSUKyisbPefOTnLc0ckHd3b1w5jeAehsQws+GsO9qgZpd2IG0vOra+Sas7B2Q6tv132kodep+7dq6m/f+Os08N51X6cV793QKzd1bTT2Qw3/DkIzjrldLLc/5w2HV//3LquavVhSdY0Ullc2WYzs567EnV19cWdXP4ztHQgv17Z741CTu9IRFZVqvLv/ErxdpUmQnh3ZBX7u0izDwQSJzK57oH6dnOPJeTienCdxNPZL4ShHkJczgtq5oJu/O3p28ESvDp7o1cHDLrvP21c1unHX8xF3PV/iaMgaeDo7IqidC4K9XdAj0KPq+vBEZz+3NlNb1BJ+7krklKix/dQNyWKYMTiUCRK1XbOGhAEASiqq79DqtjUyyBp8rtZhNZ6o21TV/plmvFbd5xpp/FoTZ3N/pvb7NPd3a+RnajzhKJfDw9lR/HJXOsHPXQEfN4VdNvKNefzOcCgc5fXWeWnoFMnq/FUaPqbuA/UPaupva8x7NXhMM/7WTf1fbfyYlr9Oc/7bNRRv885XM16nzkPOTg7wrLo2PJwd4e7siABPZ5tb7drc/j0lCj+fzWiwx89S2rlI11snE7hTY6sUFRXBy8sLhYWF8PRkXQ0REZEtaO7nd9squyciIiIyASZIRERERHUwQSIiIiKqgwkSERERUR0s2W8lQ217UVGRxJEQERFRcxk+t283R40JUisVF+uXmw8JCZE4EiIiImqp4uJieHl5Nfo8p/m3kk6nw82bN+Hh4WHStWWKiooQEhKCtLQ0Lh9gRjzPlsNzbRk8z5bB82wZ5jzPgiCguLgYQUFBkMsbrzRiD1IryeVyBAcHm+31PT09efFZAM+z5fBcWwbPs2XwPFuGuc5zUz1HBizSJiIiIqqDCRIRERFRHUyQrIxSqcTy5cuhVEqzOZ+94Hm2HJ5ry+B5tgyeZ8uwhvPMIm0iIiKiOtiDRERERFQHEyQiIiKiOpggEREREdXBBImIiIioDiZIVmbdunUIDw+Hs7MzoqOjcfjwYalDsimHDh3CxIkTERQUBJlMhp07d9Z6XhAEvP766wgKCoKLiwtGjhyJc+fO1TpGpVJhwYIF8PPzg5ubGyZNmoT09HQL/hbWbeXKlbjjjjvg4eEBf39/PPjgg7h48WKtY3iejbd+/Xr07dtXXChv6NCh+Pnnn8XneY7NY+XKlZDJZFi0aJH4GM+1abz++uuQyWS1vgIDA8Xnre48C2Q1tm7dKjg5OQkff/yxkJSUJCxcuFBwc3MTrl+/LnVoNmP37t3CsmXLhO+++04AIOzYsaPW82+//bbg4eEhfPfdd0JiYqIwbdo0oUOHDkJRUZF4zLx584SOHTsKsbGxwqlTp4RRo0YJUVFRgkajsfBvY53Gjh0rfPbZZ8LZs2eFhIQEYfz48UJoaKhQUlIiHsPzbLxdu3YJP/30k3Dx4kXh4sWLwmuvvSY4OTkJZ8+eFQSB59gcTpw4IXTq1Eno27evsHDhQvFxnmvTWL58udCnTx8hIyND/MrOzhaft7bzzATJigwaNEiYN29ercd69uwpvPrqqxJFZNvqJkg6nU4IDAwU3n77bfGxiooKwcvLS9iwYYMgCIJQUFAgODk5CVu3bhWPuXHjhiCXy4U9e/ZYLHZbkp2dLQAQDh48KAgCz7M5eXt7C5988gnPsRkUFxcL3bp1E2JjY4W7775bTJB4rk1n+fLlQlRUVIPPWeN55hCblVCr1YiLi8OYMWNqPT5mzBgcPXpUoqjaluTkZGRmZtY6x0qlEnfffbd4juPi4lBZWVnrmKCgIERERPDv0IjCwkIAgI+PDwCeZ3PQarXYunUrSktLMXToUJ5jM3juuecwfvx43HvvvbUe57k2rcuXLyMoKAjh4eGYPn06rl27BsA6zzM3q7USOTk50Gq1CAgIqPV4QEAAMjMzJYqqbTGcx4bO8fXr18VjFAoFvL296x3Dv0N9giBg8eLFGD58OCIiIgDwPJtSYmIihg4dioqKCri7u2PHjh3o3bu3+GHAc2waW7duxalTp/Dnn3/We47/n01n8ODB+OKLL9C9e3dkZWXhzTffxLBhw3Du3DmrPM9MkKyMTCar9b0gCPUeI+O05hzz79Cw559/HmfOnMGRI0fqPcfzbLwePXogISEBBQUF+O677zBnzhwcPHhQfJ7n2HhpaWlYuHAh9u3bB2dn50aP47k2XkxMjPjvyMhIDB06FF26dMHnn3+OIUOGALCu88whNivh5+cHBweHellwdnZ2vYyaWscwW6KpcxwYGAi1Wo38/PxGjyG9BQsWYNeuXThw4ACCg4PFx3meTUehUKBr164YOHAgVq5ciaioKKxZs4bn2ITi4uKQnZ2N6OhoODo6wtHREQcPHsR7770HR0dH8VzxXJuem5sbIiMjcfnyZav8P80EyUooFApER0cjNja21uOxsbEYNmyYRFG1LeHh4QgMDKx1jtVqNQ4ePCie4+joaDg5OdU6JiMjA2fPnuXfoYogCHj++eexfft2/PrrrwgPD6/1PM+z+QiCAJVKxXNsQqNHj0ZiYiISEhLEr4EDB2LmzJlISEhA586dea7NRKVS4fz58+jQoYN1/p82edk3tZphmv+nn34qJCUlCYsWLRLc3NyElJQUqUOzGcXFxUJ8fLwQHx8vABBWrVolxMfHi0slvP3224KXl5ewfft2ITExUXj00UcbnEYaHBws7N+/Xzh16pRwzz33cLpuDc8++6zg5eUl/Pbbb7Wm65aVlYnH8Dwbb+nSpcKhQ4eE5ORk4cyZM8Jrr70myOVyYd++fYIg8BybU81ZbILAc20qL774ovDbb78J165dE44dOyZMmDBB8PDwED/jrO08M0GyMh988IEQFhYmKBQKYcCAAeLUaWqeAwcOCADqfc2ZM0cQBP1U0uXLlwuBgYGCUqkU7rrrLiExMbHWa5SXlwvPP/+84OPjI7i4uAgTJkwQUlNTJfhtrFND5xeA8Nlnn4nH8Dwb78knnxTbgvbt2wujR48WkyNB4Dk2p7oJEs+1aRjWNXJychKCgoKEyZMnC+fOnROft7bzLBMEQTB9vxQRERGR7WINEhEREVEdTJCIiIiI6mCCRERERFQHEyQiIiKiOpggEREREdXBBImIiIioDiZIRERERHUwQSIiIiKqgwkSERERUR1MkIiIiIjqYIJEREREVIej1AHYKp1Oh5s3b8LDwwMymUzqcIiIiKgZBEFAcXExgoKCIJc33k/EBKmVbt68iZCQEKnDICIiolZIS0tDcHBwo88zQWolDw8PAPoT7OnpKXE0RERE1BxFRUUICQkRP8cbwwSplQzDap6enkyQiIiIbMztymNYpE1ERERUBxMkIiIiojqYIBERERHVwRokImozzt4oxJYTqQj2dsUzd3WGXM4lOIgMBEHArtM3cfhyDsb37YBRPfylDsmqMUEiIpsnCAK2nEjD67vOQa3VAQC8XZ0wfVCoxJERWYcytQavbU/EzoSbAIDvE27g4EujENTOReLIrBeH2IjIplVUavHSt2fw2o5EqLU6GDqNdp/NlDYwIiuRklOKyeuOYmfCTThUXSCVWgH7z2dJHJl1Y4JERDYrNbcMk9cdxbdx6ZDLgFdjeuLbZ4cBAJJuFkkcHZH0YpOyMHHtEVzILIafuxJfPTUYz4/qCgA4n8FrpCkcYiMim3TgQjYWbUtAYXklfN0UeP/R/hjW1Q+F5ZUAgJwSFSoqtXB2cpA4UiLL0+oErIq9iA8OXAUARId5Y93MAQjwdEZ6fjkAIC2vXMoQrR4TJCKyKVqdgDW/XMb7v16GIAD9Qtph/awB6OClr6XwdHaEwlEOtUaHW8UqhPi4ShwxkWXllarxwpZ4HLmSAwB4fFgnvHZ/Lygc9YNG/p5KAPqbCGocEyQishkFZWos3JqAg5duAQBmDwnDXyf0gtKxupdIJpOhvbsSNwrKcauECRLZl4S0AszfFIebhRVwcXLA2w9H4oF+HWsd4+euT5BuFTNBagoTJCKyCWdvFGLepjik55fD2UmOtx6KxOQBDW806euuwI2CcuSVqC0cJZE0BEHAVydS8cauJKi1OoT7uWHDrGj0CKy/35ivuwIAkFemhk4ncDmMRjBBIiKr9/Wfafjr92eh1ugQ5uuK9TOj0Tuo8T0Q3ZX6pq1UrbFUiESSqajU4q87z+LbuHQAwJjeAfjPI1HwdHZq8HjD9SEIQHmlFm5KpgIN4VkhIqtVUanFGz+cw5YTaQCAe3v547+P9IOXS8MNv4GhwS9RMUGiti01twzzNsUhKaMIchnw0tiemHd35yY3YnVxcoBcBugE/TXCBKlhVjXN/9ChQ5g4cSKCgoIgk8mwc+fOWs8//vjjkMlktb6GDBlS6xiVSoUFCxbAz88Pbm5umDRpEtLT02sdk5+fj9mzZ8PLywteXl6YPXs2CgoKzPzbEVFLpOeXYeqGP7DlRBpkMmDJmO74aPbA2yZHQI0eJCZI1IYduJCNiWuPICmjCL5uCnw5dzCeHdnltrvUy2Qy3kQ0g1UlSKWlpYiKisLatWsbPWbcuHHIyMgQv3bv3l3r+UWLFmHHjh3YunUrjhw5gpKSEkyYMAFarVY8ZsaMGUhISMCePXuwZ88eJCQkYPbs2Wb7vYioZQ5duoUJ7x9B4o1CeLs64fMnBuH5e7o1u1bCXWz8tbc5ksj26HQC3o29hCc//xOF5ZXoF9IOPywYjju7+jX7NXgTcXtW1a8WExODmJiYJo9RKpUIDAxs8LnCwkJ8+umn+PLLL3HvvfcCADZt2oSQkBDs378fY8eOxfnz57Fnzx4cO3YMgwcPBgB8/PHHGDp0KC5evIgePXqY9pciombT6QR8cOAKVu2/BEEA+gZ7Yd3MAQj2btlMNPHuuIKNP7UtzZnJ2Rzu7EG6LavqQWqO3377Df7+/ujevTuefvppZGdni8/FxcWhsrISY8aMER8LCgpCREQEjh49CgD4448/4OXlJSZHADBkyBB4eXmJxzREpVKhqKio1hcRmU5heSWe/uIk/hurT44eHRSKr58Z2uLkCADcFPoPizIWaVMbcvZGISa8fwQHL92C0lGO/06Nwj8ejGhxcgQAroZrhL2sjbKqHqTbiYmJwdSpUxEWFobk5GT87W9/wz333IO4uDgolUpkZmZCoVDA29u71s8FBAQgM1O/L1NmZib8/evvYOzv7y8e05CVK1fijTfeMO0vREQA9NuCPLs5Dtdzy6BwlOPNByPwyMCQVr+e0kl/76fW6EwVIpGkvj6Zhr/u1M/kDPVxxfpZA9AnyKvVr2dIqgybO1N9NpUgTZs2Tfx3REQEBg4ciLCwMPz000+YPHlyoz8nCEKtorWGCtjqHlPX0qVLsXjxYvH7oqIihIS0vgEnIr3v4tKxbGciKip1CPZ2wYZZ0Yjo2PqGHwAUDvoEScXGn2xc3Zmco3v6Y9Uj/eDlevvJCk0xrKrNm4jG2VSCVFeHDh0QFhaGy5cvAwACAwOhVquRn59fqxcpOzsbw4YNE4/Jyqq/g/GtW7cQEBDQ6HsplUoolUoT/wZE9kul0eIfPyZh07FUAMDd3dtjzfR+aOeqMPq1lVX7r6kq2fiT7UrPL8P8zadwJr0QMhmw+N7ueG5UV5Ms7KisSpBUGg6xNcbmapBqys3NRVpaGjp06AAAiI6OhpOTE2JjY8VjMjIycPbsWTFBGjp0KAoLC3HixAnxmOPHj6OwsFA8hojMK6OwHNM+PIZNx1IhkwELR3fDZ4/fYZLkCKjuQeLwAdkqw0zOM+mFaOfqhI1PDMKC0c2fyXk7hmFoFXuQGmVVPUglJSW4cuWK+H1ycjISEhLg4+MDHx8fvP7663j44YfRoUMHpKSk4LXXXoOfnx8eeughAICXlxfmzp2LF198Eb6+vvDx8cGSJUsQGRkpzmrr1asXxo0bh6effhoffvghAOD//u//MGHCBM5gI7KAo1dysGBLPHJL1fByccLq6f0wqkf9ukBjiI1/Je+OybbodALW/XZFnKwQ2VE/k9PUewqKNxFMkBplVQnSyZMnMWrUKPF7Q83PnDlzsH79eiQmJuKLL75AQUEBOnTogFGjRmHbtm3w8Kjea+bdd9+Fo6MjHnnkEZSXl2P06NHYuHEjHByqq/w3b96MF154QZztNmnSpCbXXiIi4wmCgA0Hr+Hfey9AJwB9gjyxYVa0WTaTZQ8S2aLC8kq8+HUC9p/Xz86efkcIXp/UB85OLZ+ldjuGIm32IDXOqhKkkSNHQhCERp/fu3fvbV/D2dkZ77//Pt5///1Gj/Hx8cGmTZtaFSMRtVxRRSVe+uY09p7T1/9NiQ7Gmw9GmKXhB1iDRLan7kzOfzzQB9PuCDXb+3GI7fasKkEiorbnYmYx5m2KQ3JOKRQOcrw+qQ8eHRRy2+0QjMECVLIl20+l47Ud+pmcHdvpZ3JGBhs3k/N2xJmevEYaxQSJiMzm+4QbePW7RJRXahHk5Yz1s6IRFdLO7O8rTmHmEBtZMbVGh3/8mIQvj10HANzVvT3WTOsHbzfTTFZoSnWdHq+RxjBBIiKTU2t0eGv3eWw8mgIAGNHND2um94ePBRp+oEYPEht/slIZheWYv/kU4lMLAAAvjO6GhaO7wcFEs9RuhwtF3h4TJCIyqayiCszffApx1/MBAM+P6oq/3NfdYg0/UJ0gsfEna1RzJqensyNWT++He3o2vg6fOSh4E3FbTJCIyGSOXcvF81/FI6dEBQ9nR7z7SD/c29uyDT9QY4YOG3+yIoIg4MND1/DOHv1Mzl4dPPHhrGiE+pp+JuftsE7v9pggEZHRBEHAJ4eT8faeC9DqBPQM9MCGWdHo5OcmSTxOVQWolexBIitRXFGJJTVmcj48QD+T00Vhnpmct+NYdY1otI3PHLd3TJCIyCglKg1e/vY0difqN3t+qH9HvPVQpGQNPwA4OuiH8zQ64bb7LBKZ26WsYsz7Mg7XqmZyLp/UGzMGhUr6/9JJbrhGeBPRGCZIRNRqV7KL8cyXcbh6qxRODjL8fUJvzBoSJnlC4lij3kmrE8SEicjSdp2+iVe+PSPO5Fw3Kxr9LDCT83YcxV5W9iA1hgkSEbXKT2cy8PK3p1Gq1iLQ0xnrZg3AgFDv2/+gBRgaf0Dfi+QoXWcW2alKrX4m52e/pwAA7uzqi/em94evu3Vsem64idDqmCA1hgkSEbWIRqvDv/ZcwMeHkwEAQzv74v0Z/eFnJQ0/ULsHqVKrM9uK3UQNya6ayXmyaibn/JFd8OKYHhadyXk7hl5V1uk1jgkSETVbdnEFnv8qHieS8wAA8+7ugiVjutfqsbEGTjV7kDiEQBZ0/FounjPM5FQ64r+PRGFMn0Cpw6rHUV5VpM0epEYZnSClpKTg8OHDSElJQVlZGdq3b4/+/ftj6NChcHZ2NkWMRGQFTqbkYf7mU8guVsFd6Yj/TO2LcREdpA6rQQ5yGWQyQBCAShahkgUIgoBPjyRj5c/6mZw9AjywYXY0wiWayXk7ToaJDOxBalSrE6SvvvoK7733Hk6cOAF/f3907NgRLi4uyMvLw9WrV+Hs7IyZM2filVdeQVhYmCljJiILEgQBG4+m4J8/nYdGJ6Cbvzs2zI5Gl/buUofWJEe5DJVagTUWZHYlKg1e+e4MfjqTAQB4oF8QVk6OhKvCegdpWKR9e6366w0YMAByuRyPP/44vv76a4SG1t5xWKVS4Y8//sDWrVsxcOBArFu3DlOnTjVJwERkOaUqDZZuT8Su0zcBABOjgvD25Ei4Ka234TdwlMtRqdVyiI3M6kp2CeZtisOV7BI4ymX424TeeGyo9DM5b4dF2rfXqlbuH//4B8aPH9/o80qlEiNHjsTIkSPx5ptvIjk5udUBEpE0rt3SN/yXsvQN/7LxvfD4sE5W3/AbODrIgEoWoZL5/JyYgSXf6GdyBngqsW7mAESH+UgdVrMYEiQOQTeuVQlSU8lRXX5+fvDz82vN2xCRRPaczcSSb06jRKVBew99w39HJ9to+A0MhdosQiVT02h1eGfvRXx06BoAYHC4D9bOGID2HtYzk/N2uJL27Zmsnzw7OxvZ2dnQ1clG+/bta6q3ICIz02h1+M++S9hw8CoAYFAnH6yd0R/+nrY34cIwpZofAGRKt4pVWLDlFI5d08/k/L+7OuPlsT2sbibn7bBI+/aMTpDi4uIwZ84cnD9/HoKgb4hkMpm4vL9Wy43wiGxBTokKC76Kxx/XcgEATw0PxysxPWtNmbcl3EqBTC3uun4mZ1aRCm4KB/x7ahTuj7TOmZy3Y5jmX8ke1kYZnSA98cQT6N69Oz799FMEBATYTH2CPUu6WYSvT6bh1wvZuFFQDl83BR7q3xF/ua87F9SzU6dS8zF/0ylkFlXAVeGAd6b0xYS+QVKHZRTO0iFTEQQBX/xxHf/4MQkanYCu/u7YMCsaXf2teyZnUwwLRbJIu3FGJ0jJycnYvn07unbtaop4yIyuZJfgv/su4uezmbUezy5W4cND13AqNR+bnhoMJfdlsBuCIGDTsetY8WMSKrUCurR3w4ZZ0egW4CF1aEZz5BACmUCZWoPXtidiZ4J+Juf4yA7415S+cLeBmZxNEYu0eX00yui/8OjRo3H69GkmSFZMpdFi3YGrWPfbFVRqBchkwP0RHfBQ/47oFeSJxPQCvPTtGfyZko8PDlzF4vu6Sx0yWUC5WotlOxKxPf4GAOD+yEC8MyXK5ht+A6eqIQTeIVNrJeeUYt6XcbiYVQwHuQxLY3pi7vDwNjFS4sQi7dsyuiX85JNPMGfOHJw9exYRERFwcnKq9fykSZOMfQsyQtz1PLz87RlcvVUKABjVoz1ejemFHoHVPQQd27lAoxPw/Ffx+OTwNTx5Zye0c1VIFTJZwPXcUjzzZRwuZOob/lfH9cRTI9pGw2/gIE5j5gcAtdy+c5l48evTKFZp4OeuxAcz+mNwZ1+pwzIZsYeVNXqNMjpBOnr0KI4cOYKff/653nMs0paORqvDe79ewdpfL0MnAH7uSrw+qTfGR3Zo8ENwfGQHrA28gguZxfjh9E3MHtrJ8kGTRexPysJfvk5AcYUGfu4KrJ0xAEPaUMNvwFk61BpanYD/7ruIdb/pZ3IODPPGBzMHIMAGZ3I2RSzS1gripCqqzejpKS+88AJmz56NjIwM6HS6Wl9MjqSRlleGaR8dw3u/6JOjyf074pfFd2NC36BGLwKZTIZHBoYAgLhqMrUtWp2A/+y9iKe+OIniCg2iw7zx44IRbTI5AlikTS2XW6LCnP+dEJOjJ+7shC3/N6TNJUdAdQ0SALCTtWFG9yDl5ubiL3/5CwICAkwRDxlp1+mbWLY9EcUqDTyUjnjzoQg80K9js372vt4BWPFjEuJTC1Ci0rSZWhQC8krVWLg1Hocv5wAAHh/WCa/d3wsKR9ucwt8cjpzmTy2QkFaA+ZvicLOwAi5ODvjXlL6YFGXbMzmbYhhiA/SF2g5yTs6py+hPwMmTJ+PAgQPo0qWLKeKhVipVafD378/hu1PpAIDoMG+sntYPIT6uzX6NEB9XhPm64npuGU4k5+Kenkx624LTaQWYv/kUbhSUw8XJAW8/HNnspNmWcRozNYcgCNh8PBUrfkiCWqtDZz83bJgdje5tYCZnU2qub8bV5htmdILUvXt3LF26FEeOHEFkZGS9Iu0XXnjB2Leg27iUVYx5m+Jw7VYp5DJgwT3dsOCerq1a2XVgmA+u55bhTHohEyQbJwgCtv6ZhuXfn4Naq0MnX1dsmB2NnoGeUodmETVrLIgaUlGpxbIdZ8Uby7F9AvCfqVHwcHa6zU/avppDbKzTa5hJZrG5u7vj4MGDOHjwYK3nZDIZEyQz234qHct2nEV5pRaBns5479H+GBTe+j2z+gR54rtTwNkbRSaMkiytolKLv39/Fl+f1Df89/UOwH8fiYKnHTT8BizSpqak5pZh3qY4JGUUQS4DXh7XE8/c1dluipUd5DWH2HgT0RCTLBRJlldRqcUbP5zDlhNpAIAR3fywelo/+Lobt1liREcvAMC5m4VGx0jSSMvTN/znbuob/iVje2DeXV0gl9tHw2/ArRSoMQcuZGPh1ngUVWjg66bA+4/2x7Cu9rWpukwmg6NcBo1OYJ1eI6yqQvPQoUOYOHEigoL0s6127txZ63lBEPD6668jKCgILi4uGDlyJM6dO1frGJVKhQULFsDPzw9ubm6YNGkS0tPTax2Tn5+P2bNnw8vLC15eXpg9ezYKCgrM/NuZzvXcUjy8/ii2nEiDTAYsurcbNj4xyOjkCAC6B+iXzs8orECZWmP065FlHbiYjQnvH8G5m0XwcVPgy7mDMX9kV7tLjgDAwVCDxB4kqqLVCVgVewlPbPwTRRUa9Atphx9fGG53yZGBoReJdXoNa1WC9Pbbb6OsrKxZxx4/fhw//fRTs44tLS1FVFQU1q5d2+Dz77zzDlatWoW1a9fizz//RGBgIO677z4UFxeLxyxatAg7duzA1q1bceTIEZSUlGDChAm1lhyYMWMGEhISsGfPHuzZswcJCQmYPXt2s2KU2p6zmbU+AL94chAW3du9VnepMdq5KuDloh+GSc1r3t+YpKfTCVi9/xKe3PgnCssrERXSDj8uGI477bThBwCHqqESjh4QAOSXqvHExj/x3i+XAQCzh4Rh2zND0MHLReLIpGP43GAHUsNaNcSWlJSE0NBQTJ06FZMmTcLAgQPRvn17AIBGo0FSUhKOHDmCTZs2ISMjA1988UWzXjcmJgYxMTENPicIAlavXo1ly5Zh8uTJAIDPP/8cAQEB+Oqrr/DMM8+gsLAQn376Kb788kvce++9AIBNmzYhJCQE+/fvx9ixY3H+/Hns2bMHx44dw+DBgwEAH3/8MYYOHYqLFy+iR48erTklZlep1eFfP1/AJ0f0Q5rRYd5YO6O/WS7uTn5uOJ1WgJScUrsp6LVlBWVq/GVbAg5cvAUAmDUkFH+b0Nvu99SrbvyZIdm7xPRCzNsUhxsF5XB2kuOthyIxeUCw1GFJrvomgtdIQ1rVg/TFF1/g119/hU6nw8yZMxEYGAiFQgEPDw8olUr0798f//vf//D444/jwoULGDFihNGBJicnIzMzE2PGjBEfUyqVuPvuu3H06FEAQFxcHCorK2sdExQUhIiICPGYP/74A15eXmJyBABDhgyBl5eXeExDVCoVioqKan1ZSkZhOR796JiYHD09Ihxb/898dz6dfPVLAyTnsAfJ2p29UYgJ7x/BgYu3oHSU479To/Dmg5F2nxwBgJyNPwH4+s80PLzhKG4UlCPM1xU75t/J5KiKOAzNLqQGtbpIu2/fvvjwww+xYcMGnDlzBikpKSgvL4efnx/69esHPz/Tdu1nZup3oK+7IGVAQACuX78uHqNQKODt7V3vGMPPZ2Zmwt/fv97r+/v7i8c0ZOXKlXjjjTeM+h1a4/DlW1i4NQF5pWp4KB3x76lRGBcRaNb3DPHWJ0gZheVmfR8yztcn0/C3nWeh0ugQ6uOK9bMGoE+Ql9RhWQ1H1lfYtboTWe7t5Y//PtJPLCGgGj1IzI8aZPQsNplMhqioKERFRZkinma9X03N2UOm7jENHX+711m6dCkWL14sfl9UVISQkJDmht1iWp2A93+9jDW/XIYg6Kffr5s5AGG+bmZ7T4MAL/2y+pmFFWZ/L2o5lUaL13clYcuJVADAPT398e4j/eDlyoa/JjkTJLuVnl+GZzedQuKNQshkwIv3dbfbyQpN4TXSNJvZSyIwUN9rkpmZiQ4dOoiPZ2dni71KgYGBUKvVyM/Pr9WLlJ2djWHDhonHZGVl1Xv9W7duNbldilKphFJp/Cyx5sgtUWHRtgRxW4hHB4Vi+cTecHayzLBJgIf+98wqVlnk/aj5bhSUY/6mOJxO1zf8i+/tjudGseFviGGdVDb+9uXgpVtYuDUeBWWV8HZ1wprp/XFX9/ZSh2WVDD1IOg5DN8iqpvk3JTw8HIGBgYiNjRUfU6vVOHjwoJj8REdHw8nJqdYxGRkZOHv2rHjM0KFDUVhYiBMnTojHHD9+HIWFheIxUjqZkofx7x3B4cs5cHFywKpHorBycqTFkiMA4saM2UXsQbImhy/fwoT3DuN0eiHauTph4xODsGB0NyZHjWDjb190OgHv/3IZj392AgVllegb7IUfFgxnctQETvNvmlX1IJWUlODKlSvi98nJyUhISICPjw9CQ0OxaNEivPXWW+jWrRu6deuGt956C66urpgxYwYAwMvLC3PnzsWLL74IX19f+Pj4YMmSJYiMjBRntfXq1Qvjxo3D008/jQ8//BAA8H//93+YMGGC5DPYCsrUmPO/EyhVa9G5vRs2zJJmPyAxQSpWQacT+AEsMZ1OwPqDV/GffRchCEBkRy+smzmgRfvs2SMOH9iPwvJKLN6WgF8uZAOwfK+7rXIQN3TmNdIQq0qQTp48iVGjRonfG2p+5syZg40bN+Lll19GeXk55s+fj/z8fAwePBj79u2Dh0d1EvHuu+/C0dERjzzyCMrLyzF69Ghs3LgRDg7VF8rmzZvxwgsviLPdJk2a1OjaS5bUzlWBpff3wvHkPKycHAl3pTR/Hj93BWQy/QdLTqkK/h7OksRB+ob/xa8TsP+8vuGffkcIXp/Uhw1/M3AKs31IulmEeZvikJpXBoWjHG8+EIFH7jBffWhbIi6FwWukQTJB4JlpjaKiInh5eaGwsBCenqZbK8jw55B6P6Dof8Qit1SNPYtGcC0kiZzP0Df813P1Df8/HuiDaXeESh2WzVj583l8ePAanhoejr9O6C11OGQG38Wl47UdiVBpdAj2dsGGWdHidkl0e6P/+xuu3irF1v8bgiGdfaUOx2Ka+/ndqi4Kw0KNzbF9+/bWvIXdkjoxMvBydUJuqRoFZZVSh2KXdsSnY+n2RFRU6tCxnb7hjwxmw98S7EFqu1QaLf7xYxI2HdPP5Ly7e3usmd4P7VwVEkdmW7iYatNalSB5ebGhbuu8XRUASlFQppY6FLui1ujwjx+T8OUx/dped3VvjzXT+sHbjQ1/S7Hxb5syCsvx7KZTSEgrgEwGvHBPNyzkZIVWMSymyhqkhrUqQfrss89MHQdZmXZVi6mxB8lyMgrLMX/zKcSnFgAAXhitb/hNtc+eveFK2m3P0Ss5WLAlHrmlang6O2LN9P4Y1bP+wr/UPI4OvEaaYpIqYI1Gg99++w1Xr17FjBkz4OHhgZs3b8LT0xPu7u6meAuyMENXdT4TJIs4ejUHC76qbvhXT++He3o2vi4X3V71FGaJAyGjCYKADQev4d97L0AnAL07eGLDrGiE+nImpzHEpTDYg9QgoxOk69evY9y4cUhNTYVKpcJ9990HDw8PvPPOO6ioqMCGDRtMESdZWLuqVZkLyjnEZk6CIODDQ9fwzh59w9+rgyc+ZMNvEtUJEjMkW1ZUUYmXvjmNvef0C/xOiQ7Gmw9GcCanCXApjKYZnSAtXLgQAwcOxOnTp+HrW10F/9BDD+Gpp54y9uVJIt6GBKmUPUjmUlxRiSU1Gv6HBwTjnw+x4TcV9iDZvouZxZi3KQ7JOaVQOMjx+qQ+eHRQiNVMZrF11XuxMUFqiNEJ0pEjR/D7779DoahdRBoWFoYbN24Y+/IkEa+qITb2IJnHpaxizPsyDteqGv7lk3pjxqBQNvwmxJW0bdv3CTfw6neJKK/UIsjLGetnRSMqpJ3UYbUp4k0Er5EGGZ0g6XQ6aLXaeo+np6fXWsCRbAuLtM1n1+mbeOXbM2LDv25WNPqx4Tc5Dh/YJrVGh7d2n8fGoykAgOFd/fDeo/3hw5mcJsetRppm9F5s9913H1avXi1+L5PJUFJSguXLl+P+++839uVJIoZVvEvVGokjaTsqtTq88cM5vLAlHuWVWtzZ1Rc/LBjO5MhMqibo8O7YhmQVVWDGx8fE5Oi5UV3w+ZODmByZCVfSbprRPUjvvvsuRo0ahd69e6OiogIzZszA5cuX4efnhy1btpgiRpKAmyFBUtXvHaSWyy6qwPzNp3Dyej4AYP7ILnhxTA9O4TcjBwf9/Z9Wy8bfFhy/lovnvopHTokKHkpHrJrWD/f15kxOcxLXQeI10iCjE6SgoCAkJCRgy5YtOHXqFHQ6HebOnYuZM2fCxcXFFDGSBNyU+kLhEhV7kIxVt+H/7yNRGNMnUOqw2jyupG0bBEHAp0eSsfLnC9DqBPQM9MD6WdEI93OTOrQ2jz1ITTM6QSorK4OrqyuefPJJPPnkk6aIiayAOMTGBKnV6jb8PQI8sGE2G35LqepA4hovVqxEpcEr357BT4kZAIAH+wXhrcmRcFVY1T7qbRZnejbN6P+F/v7+ePDBBzF79mzcd999kMuNLmsiK2AYYitTa6HTCVzGv4XY8EuPK2lbtyvZxXjmyzhcvVUKJwcZ/jahN2YPCeNMTgtiL2vTjM5mvvjiC6hUKjz00EMICgrCwoUL8eeff5oiNpKQoQcJYKF2S13JLsGDH/yOnxIz4CiXYcUDffDutH5MjiyMM3Ss1+7EDDyw9ndcvVWKAE8ltv7fUDw2tBOTIwsTrxF2ITXI6ARp8uTJ+Oabb5CVlYWVK1fi/PnzGDZsGLp3744VK1aYIkaSgNJRLl48LNRuPn3DfwRXsksQ4KnEtmeGsOGXCBMk66PR6vDPn5Iwf/MplKq1GNLZBz8uGIHoMG+pQ7NL4lIYvEQaZLLxMA8PDzzxxBPYt28fTp8+DTc3N7zxxhumenmyMJlMBjcFC7Wbq/GG30fq0OwWEyTrkl1cgZmfHMfHh5MBAM/c1Rmb5g5Gew+lxJHZL8NSGKzTa5jJ+vwrKiqwa9cufPXVV9izZw/8/f2xZMkSU708ScDD2QlFFRomSLdxq1iF5786hePJeQD0Df9LY3vA0YH1eFLiStrW42RKHuZvPoXsYhXclY7495S+iInsIHVYds+hqmaYNUgNMzpB2rdvHzZv3oydO3fCwcEBU6ZMwd69e3H33XebIj6SkGGqP2eyNS7uur7hzypSwU3hgP9MjWLDbyW4krb0BEHA50dT8OZP56HRCejm744Ns6PRpb271KERqmd68hppmNEJ0oMPPojx48fj888/x/jx4+Hk5GSKuMgKGIqKmSDVV7fh7+rvjg2zotHVnw2/taieoSNxIHaqTK3B0u2J+D7hJgBgQt8O+NfDfcUZsiQ9DkM3zej/qZmZmfD09DRFLGRlXKp2la/QcIZDTWz4bYODg6Hx5/9fS7t2qwTPbjqFi1nFcJDL8Nr9vfDknZysYG3EpTCYIDXI6Bbd09MTV69exWeffYarV69izZo18Pf3x549exASEoI+ffqYIk6SgLOTvv+1opKz2AzY8NsOsQeJ+ZFF7T2XiSVfn0axSoP2Hkp8MGMABoVzsoI1cuRK2k0yuor04MGDiIyMxPHjx7F9+3aUlJQAAM6cOYPly5cbHSBJx9nQg8QECYC+4X9g7e+4mFWM9h5KbHl6COYOD2dyZKXEbRR4d2wRGq0O/9pzAc98GYdilQZ3dPLGTwuGMzmyYqzTa5rRCdKrr76KN998E7GxsVAoqndcHjVqFP744w9jX54k5MIECQAbflvFlbQtJ7dEhTmfncD6364CAOYOD8dXTw+Bv6ezxJFRUxw4xNYko4fYEhMT8dVXX9V7vH379sjNzTX25UlCSjFBst8xipwSFV7YEo+jV/X/l+cOD8erMT3hxCn8Vo89SJaRkFaAZzfFIaOwAq4KB/zr4b6YGBUkdVjUDCzSbprRCVK7du2QkZGB8PDwWo/Hx8ejY8eOxr48SchQg1Rupz1I8an5mL/5FBt+G2Vo/DVs/M1CEARsPp6KN344h0qtgM5+btgwOxrdAzykDo2aqXolbV4jDTE6QZoxYwZeeeUVfPPNN5DJZNDpdPj999+xZMkSPPbYY6aIkSRir0Ns9Rr+9m7YMIsNv63h3bH5lKu1WLYzEdtP3QAAjOsTiH9P7QsPZy7zYksc2cvaJKMTpH/+8594/PHH0bFjRwiCgN69e0Or1WLGjBn461//aooYSSLOdjjExoa/7eBK2uZxPbcU8zadwvmMIshlwCvjeuL/7urMyQo2yFCnx17WhhmdIDk5OWHz5s1YsWIF4uPjodPp0L9/f3Tr1s0U8ZGE7G2a//XcUjzzZRwuZBaz4W8D5Fwl2OR+vZCFRVsTUFShga+bAu/P6I9hXfykDotayYHT/JtkspXtunTpgi5dupjq5cgK2NMQ2y/ns7BoWwKKKzTwc1fgvUfZ8Ns6Nv6mo9UJWLP/Et779QoAoH9oO6ybOQAdvFwkjoyMwWHoprUqQVq8eHGzj121alVr3oKsgNIOEiStTsDq/ZfwflXDPyC0HdbNjEagF6cn2zpOYTaN/FI1Fm5LwKFLtwAAjw0Nw1/H94bCkTM5bZ2ci6k2qVUJUnx8fLOOM/XQxOuvv4433nij1mMBAQHIzMwEoC+ufeONN/DRRx8hPz8fgwcPxgcffFBrNW+VSoUlS5Zgy5YtKC8vx+jRo7Fu3ToEBwebNNa2wFCD1FZnseWXqvHC1ngcvpwDAJgzNAzL2PC3GVwEz3iJ6YWYtykONwrK4ewkx8rJkXioP9vKtsKwWgl7WRvWqgTpwIEDpo6j2fr06YP9+/eL3zs4OIj/fuedd7Bq1Sps3LgR3bt3x5tvvon77rsPFy9ehIeHfgbSokWL8MMPP2Dr1q3w9fXFiy++iAkTJiAuLq7Wa1HNIba2d3txJr0Az246JTb8b0/uiwf7c1mKtqS6SFviQGzUtj9T8bfvz0Gt0SHM1xUbZkWjVwfuu9mWcC+2ptnc7pqOjo4IDAys97ggCFi9ejWWLVuGyZMnAwA+//xzBAQE4KuvvsIzzzyDwsJCfPrpp/jyyy9x7733AgA2bdqEkJAQ7N+/H2PHjrXo72Lt2mqR9tYTqfj79+eg1urQydcV69nwt0nV6yC1vQTfnCoqtVj+/TlsO5kGALi3lz/++0g/eLlwJmdb48B1kJpkc2MJly9fRlBQEMLDwzF9+nRcu3YNAJCcnIzMzEyMGTNGPFapVOLuu+/G0aNHAQBxcXGorKysdUxQUBAiIiLEYxqjUqlQVFRU66uta2t7sVVUavHyt6fx6vZEqLU63NsrAN8/P5zJURtVvZK2xIHYkLS8Mkzd8Ae2nUyDXAa8NLYHPpo9kMlRG8XV5ptmUz1IgwcPxhdffIHu3bsjKysLb775JoYNG4Zz586JdUgBAQG1fiYgIADXr18HAGRmZkKhUMDb27veMYafb8zKlSvr1T+1dW1piC0trwzPbo7D2Rv6tVteHNMDz97dRaxTobaHd8ctc/DSLSzcGo+Cskp4uzrhvUf7Y0S39lKHRWbEIbam2VSCFBMTI/47MjISQ4cORZcuXfD5559jyJAhAOoXhguCcNti8eYcs3Tp0lqz94qKihASEtLSX8GmiENsGtvuQfrtYjYWbUtAQVklfNwUeG96fwzvxin8bR0b/+bR6QSsPXAF7+6/BEEA+gZ7Yf2saHRsxyn8bR2XwmiazQ2x1eTm5obIyEhcvnxZrEuq2xOUnZ0t9ioFBgZCrVYjPz+/0WMao1Qq4enpWeurrVM62vYQm04nYM3+y3hi458oKKtEVLAXflgwnMmRnXCo0TvIIYSGFZZV4qkvTmJVrD45enRQKL5+ZiiTIzvBpTCaZpIE6csvv8Sdd96JoKAgcThr9erV+P77703x8o1SqVQ4f/48OnTogPDwcAQGBiI2NlZ8Xq1W4+DBgxg2bBgAIDo6Gk5OTrWOycjIwNmzZ8VjqJpSLNLWQbCxO4zCskrM/fxP8a54xuBQfD2PDb89cajRK8xhtvrO3SzExLVH8OuFbCgd5XhnSl+snBwp1h5S21e9Wa3EgVgpoxOk9evXY/Hixbj//vtRUFAArVbf29CuXTusXr3a2JevZcmSJTh48CCSk5Nx/PhxTJkyBUVFRZgzZw5kMhkWLVqEt956Czt27MDZs2fx+OOPw9XVFTNmzAAAeHl5Ye7cuXjxxRfxyy+/ID4+HrNmzUJkZKQ4q42q1WwoVRrbqUM6e6MQE9YexoGLt6B0lOPfU/rirYcixR4xsg/yGq0b75Br+y4uHZPXHUVqXhlCfFzw3bPD8MjAtl0yQPWJ6yDx+miQ0TVI77//Pj7++GM8+OCDePvtt8XHBw4ciCVLlhj78rWkp6fj0UcfRU5ODtq3b48hQ4bg2LFjCAsLAwC8/PLLKC8vx/z588WFIvft2yeugQQA7777LhwdHfHII4+IC0Vu3LiRayA1wNmxdoJkC3eW35xMw193noVKo0OIjws2zIpGnyAvqcMiCdQaYmMPEgBApdHiHz8mYdOxVADAyB7tsXpaP7RzVUgcGUmBdXpNMzpBSk5ORv/+/es9rlQqUVpaauzL17J169Ymn5fJZHj99dfx+uuvN3qMs7Mz3n//fbz//vsmja0tcnKQQS7TL7SnqtQCVjzVV6XR4o0fkvDVcX3DP6pHe6ye1h9ertYbM5mXvOYQGz8AcLOgHM9uPoXTaQWQyYCFo7vhhXu6cSanHeNMz6YZnSCFh4cjISFB7MUx+Pnnn9G7d29jX54kJJPJ4OzkgDK11qqn+t8oKMf8TXE4nV4ImQxYNLo7FtzTlQ2/natdpC1hIFbg9ys5WLAlHnmlani5OGH19H4Y1cNf6rBIYuJq87yBaJDRCdJLL72E5557DhUVFRAEASdOnMCWLVuwcuVKfPLJJ6aIkSSkdJTrEyQrnep/5HIOFmw5hfyySni5OGHN9H4YyYafwCJtQL+EyfqDV/GfvRehE4A+QZ7YMCsaIT6uUodGVkDOHqQmGZ0gPfHEE9BoNHj55ZdRVlaGGTNmoGPHjlizZg2mT59uihhJQvq6o0qrm+qv0+kb/v/u0zf8ER09sX4mG36qVrMH0R6H2IoqKrHk69PYl5QFAJgaHYx/PBhhE7WEZBmOXEm7SUYnSAUFBXj66afx9NNPIycnBzqdDv7++jv4K1euoGvXrkYHSdJxtsLVtAvLK/Hi16ex/7y+4Z82MARvPNCHDT/V4yCXQasT7K5I+2JmMeZtikNyTikUDnK88UAfTL8j5LYL4pJ9YQ9S04xOkO6//378+uuvcHZ2hp9f9QJ8Fy9exOjRo5Genm7sW5CElI7WtWHt+YwiPLspDim5ZVA4yrFiUh9MHxQqdVhkpRxkMmgh2FUP0vcJN/Dqd4kor9QiyMsZ62dFIyqkndRhkRWqXihS4kCslNEJkre3Nx588EH8+OOPcHTUv9z58+dxzz334JFHHjE6QJKWNW1YuyM+HUu3J6KiUoeO7VywftYA9A1uJ3VYZMXkcgBa+xhiU2t0eGv3eWw8mgIAGNHND2um94ePG6fwU8O4WW3TjF4o8rvvvkNpaSlmzJgBQRBw9uxZjBw5Eo8++ijWrFljihhJQob92KRcKFKt0WH592fxl22nUVGpw4hufvhhwXAmR3Rb4iydNj6EkFVUgUc/PiYmR8+P6oqNTwxickRNEtdBauPXR2sZ3YPk7OyMH3/8ESNHjsTUqVNx+PBhPPbYY/j3v/9tivhIYlL3IGUUluO5zadwKrUAALDgnq5YdG/3WlO4iRoj1li04TvkY9dy8fxX8cgpUcHD2RHvPtIP9/Zuem9JIoA9SLfTqgSpqKio1vcymQzbtm3Dvffei4cffhh/+9vfxGPsYVPXtkysQZKgB+no1Ry8sCUeOSVqeDg7YvW0fhjdiw0/NZ9DG06QBEHAJ4eT8faeC9DqBPQM9MCGWdHo5OcmdWhkIwxbjbAHqWGtSpDatWvX4GwIQRCwYcMGfPjhhxAEATKZTNybjWyToQdJZcEeJEEQ8NGha/jXngvQCUCvDp7YMGsAwnzZ8FPLOLbRWTolKg1e/vY0didmAgAe6t8Rbz0UCRcFZ3JS83Grkaa1KkE6cOCAqeMgK2XYj81SQ2zFFZV46Zsz2HNO3/BPHtAR/3yQDT+1Tlv8ALiSXYxnvozD1VulcHKQ4e8TemPWkDBO4acWa8s9rKbQqgTp7rvvNnUcZKUMRdqWWAfpUpZ+7ZZrVQ3/8ol9MHNwKBt+arXqGguJAzGRn85k4OVvT6NUrUWgpzPWzRqAAaHeUodFNqot3kCYktFF2gZlZWVITU2FWq2u9Xjfvn1N9RYkAUsVaf9w+iZe+e4MytRadPByxrqZA9CfDT8Zqa3M0tFodfjXngv4+HAyAGBoZ1+8P6M//NyVEkdGtky8gbDx68NcjE6Qbt26hSeeeAI///xzg8+zBsm2KQ0Jkpn2YqvU6rBy9wX873d9wz+siy/ef7Q/fNnwkwm0hSGE7OIKPP9VPE4k5wEAnrm7M14a0wOODkav0kJ2ri1cH+Zk9BW2aNEi5Ofn49ixY3BxccGePXvw+eefo1u3bti1a5cpYiQJiesgmWGILbuoAjM+PiYmR8+O7IIvnhzE5IhMxtbvkE+m5GHCe0dwIjkP7kpHbJg1AEtjejE5IpPgEFvTjO5B+vXXX/H999/jjjvugFwuR1hYGO677z54enpi5cqVGD9+vCniJIkoDUXaJp7mfyI5D899dQq3ilXwUDriP49EYWyfQJO+B5FhuSxb+wAQBAEbj6bgnz+dh0YnoJu/OzbMjkaX9u5Sh0ZtSPUNhMSBWCmjE6TS0lJxc1ofHx/cunUL3bt3R2RkJE6dOmV0gCSt6iJt0wyxCYKA//2egrd2n4dWJ6BHgAfWzxqAzmz4yQxscSG8MrUGr36XiF2nbwIAJvTtgH893BduSpOVjBIBqLkXm+1cH5Zk9BXXo0cPXLx4EZ06dUK/fv3w4YcfolOnTtiwYQM6dOhgihhJQqac5l+q0uCV787gxzMZAIAH+gVh5eRIuCrY8JN52FqR9rVbJZi3KQ6XskrgKJfhtft74Yk7O3EmJ5mFnAtFNsnoT6ZFixYhI0P/gbd8+XKMHTsWmzdvhkKhwMaNG419eZJY9UKRxg2xXcnWN/xXsvUN/1/H98KcYWz4ybxsqQh1z9lMLPnmNEpUGrT3UGLdzAG4o5OP1GFRG2aLPayWZHSCNHPmTPHf/fv3R0pKCi5cuIDQ0FD4+fkZ+/IkMXGIzYhZbD8nZmDJN/q1W/yrGv6BbPjJAmyhSFuj1eE/+y5hw8GrAIBBnXywdkZ/+Hs6SxwZtXUONtbDamkmH9twdXXFgAEDTP2yJBFj1kHSaHX4996L+PDQNQDA4HAfvD+jP/w92PCTZVTP0pE4kEbklKjwwpZ4HL2aCwCYOzwcr8b0hBNnqZEFGDZzFgSI24NRtVYlSIsXL272satWrWrNW5CVEKf5t3AW261iFRZsOYVj1/Rrt/zfXZ3x8liu3UKWZc1DbPGp+Zi/+RQyCivgqnDAO1P6YkLfIKnDIjviUCMh0uoEODpYT4IkCALyStWSLvvSqgQpPj6+WccxG7V9ylYUacddz8f8zXHIKlLBTeGAd6ZEYXxfFuyT5Rk+AKxpiE0QBGw6nooVP5xDpVZA5/Zu+HBWNLoFeEgdGtkZQw8SoB9ms5bpMobNmC9kFGPn83fC09lJkji4WS01qSV7sQmCgC/+uI43f0pCpVZAl/Zu+HB2NLr6s+EnaYizdKykB6lcrcWynYnYfuoGACAmIhDvTOkLD4k+AMi+OdRIkKxlv8K6mzHHXc/HqB7+ksRiLQkjWanm9iCVqTV4bXsidibo1265PzIQ70yJgjvXbiEJWVOR9vXcUjzzZRwuZBZDLgNejemJp0d0Zk87SabWEJsVXCO7EzPw0jfWsxkzP72oSYbF6VQaHSq1ugaLR1NySjFvk77hd5DLsDSmJ+YOD2fDT5IzFGlrtNI2/r+cz8KibQkortDAz12B9x7tj2FdOMuXpCWv0ZxL2ctqrZsxM0GiJnm5OEEu0y9Fn1+mrjcDLTYpC4u/rm74184YgCGdfSWKlqg2sUhbortjrU7A6v2X8P6vVwAAA0LbYd3MaAR6cSYnSa9mD5JUayHdKlbh+a9O4bgVbsbMBIma5CCXwcdNgZwSNXJLqhMkrU7Au7GXsPaAvuGPDvPGBzMGsOEnq+Io4UJ4+aVqvLA1Hocv5wAAHh/WCa/d3wsKR+kbfiKgdg2SFDcRcdfzMH/zKWQVqeCudMS/p/RFTKT1TOhhgkS35eumFBMkAMgrVWMhG36yAVJtNXImvQDPbjqFGwXlcHaS4+3JffFg/44WjYHodmQyGWQy/TpIlryJMEzo+cePSdDoBHT1d8eGWdHo6m9de3IyQaLb8nVXAFlAVlEFTqcVYP5mfcPv4uSAtx+OxAP92PCTdZJiK4WtJ1Lx9+/PQa3VoZOvK9bPikavDp4We3+ilnCUy1CpFSx2E1F3Qs/4vh3wjpVuxmzXt/zr1q1DeHg4nJ2dER0djcOHD0sdklUyZPXLdiZi6oY/cKOgHJ18XbHjuWFMjsiqyS24UGRFpRYvf3sar25PhFqrw329A/D988OZHJFVq15t3vzXSHJOKR764Ch2JtyEQ9WenGsf7W+VyRFgxwnStm3bsGjRIixbtgzx8fEYMWIEYmJikJqaKnVoVsewb1pFpU5s+HctGI6egWz4ybpV7zVl3vdJyyvDlA1H8fXJdMhlwEtje+DDWdHwcuH6RmTdqntZzfs++85lYtL7R3Axqxh+7kp89dRgPGXly1xYZ9pmAatWrcLcuXPx1FNPAQBWr16NvXv3Yv369Vi5cqXE0VmXmIhAPNgvCGduFOKJO8Mxa3CoVf+nJjKwxBDbbxezsWhbAgrKKuHjpsB70/tjeDdO4SfbYO4Na7U6AatiL+KDA/rNmAeGeeODmQMQYAObMdtlgqRWqxEXF4dXX3211uNjxozB0aNHG/wZlUoFlUolfl9UVGTWGK2Jk4Mcq6f3lzoMohYzZ5G2Tidg7YEreHf/JQgCEBXshXWzotGxnYvJ34vIXMw5DJ1XqsYLW+Jx5Ip+Qs8Td+on9NjKZsx2mSDl5ORAq9UiICCg1uMBAQHIzMxs8GdWrlyJN954wxLhEZGJOJhpq5HCskr85esE/HohGwAwY3Aolk/sLa48T2QrzLXafEJaAeZvisPNwgqbndBjlwmSQd1hIkEQGh06Wrp0KRYvXix+X1RUhJCQELPGR0TGMccQ27mbhXh20ymk5pVB6SjHmw9GYOpAtgVkm0xdpC0IAracSMPru/QzOcP93LBhVjR6BNrenpx2mSD5+fnBwcGhXm9RdnZ2vV4lA6VSCaVS2mXPiahlTD3E9m1cOpbtSIRKo0OIjws2zIpGnyAvk7w2kRRM2ctaUanF33aexTdx6QCAMb0D8J9HouBpo5sx28ZAoIkpFApER0cjNja21uOxsbEYNmyYRFERkamZqgdJpdFi2Y5ELPnmNFQaHUb1aI8fnx/B5IhsnqFI29ghtrS8Mjy8/ii+idPP5HxlXE98ODvaZpMjwE57kABg8eLFmD17NgYOHIihQ4fio48+QmpqKubNmyd1aERkIqboQbpZUI5nN5/C6bQCyGTAotHdseCermJxK5EtM0WR9oGL2Vi0NQGF5fqZnO8/2h93drX9mZx2myBNmzYNubm5WLFiBTIyMhAREYHdu3cjLCxM6tCIyETEzWpbucbLkcs5eGFrPPJK1fByccKa6f0wsoe/CSMkkpYxRdo6nYD3f72C1b9UzeQMaYf1MwcgqI3M5LTbBAkA5s+fj/nz50sdBhGZSWsbf51OwPqDV/HffRehE4CIjp5YPzMaIT6u5giTSDLiOkgtvIkoLKvEom3xOHDxFgBg1pBQ/G1C25rJadcJEhG1ba2ZoVNUUYkXvz6N2KQsAMAjA4Ox4oEIODu1nYafyKA1Q2xnbxTi2c1xSMsrh9JRjn8+FIkp0cHmClEyTJCIqM1q6QydC5lFmPdlHFJyy6BwlGPFpD6YPijUjBESSaulRdr2NJOTCRIRtVkOLehB2hl/A0u3J6K8UouO7VywftYA9A1uZ+YIiaTV3B4klUaLFT8kYfNx/X6l9/T0x7uP9IOXq+3OUrsdJkhE1GaJjX8Td8dqjQ7//CkJn/9xHQAwopsf1kzvDx83hUViJJKS2MvaxDVSdybnX+7tjudHtf2ZnEyQiKjNcrzNOkiZhRWYvzkOp1ILAAAL7umKRfd2F4u7ido6cYitkWvk9ys5WLDFPmdyMkEiojarqeGDP67mYsGWU8gpUcPD2RHvPtIP9/ZueCV9orbKcDOgqXONCIJ+Jud/9trvTE4mSETUZrlUzTwrr9SKjwmCgI8PX8O/9lyEViegZ6AHNsyKRic/N6nCJJKMq0KfBpSpNeJjRRWVWPL1aeyz85mcTJCIqM3yqNrmoKhC3/iXqDR46ZvT+Pmsfh/Gyf074p8PRcJFYV8NP5GBh7M+DSiuukYuZhZj3qY4JOeUQuEgx4oH7HcmJxMkImqzPMXGvxKXs/QN/9VbpXBykOHvE/tg1uBQyGSsNyL7ZUiQisor8X3CDbz6XfVMznUzByAqpJ20AUqICRIRtVmeLvoepPjUAjzwwe8oU2sR6OmMdbMGYECot8TREUnPsJnsf/ZdEh/jTE49JkhE1GYFeDqL/y5TazGsiy/ee7Q//NyVEkZFZD1qXiMA8PyorvjLfZzJCQByqQMgIjKXTr6uGNHNDw5yGZ4b1QVfzh3M5IiohvF9O8DHTQEfNwU+nTMQS8b2YHJURSYIrdjCl1BUVAQvLy8UFhbC09NT6nCIqBFanQCdIMDJgfeDRA3RaHWQy2RtfuFHg+Z+fnOIjYjaNAe5DA6wj4afqDUcefPQIJ4VIiIiojqYIBERERHVwQSJiIiIqA4mSERERER1sEi7lQyT/4qKiiSOhIiIiJrL8Ll9u0n8TJBaqbi4GAAQEhIicSRERETUUsXFxfDy8mr0ea6D1Eo6nQ43b96Eh4eHSfdyKioqQkhICNLS0ri+khnxPFsOz7Vl8DxbBs+zZZjzPAuCgOLiYgQFBUEub7zSiD1IrSSXyxEcHGy21/f09OTFZwE8z5bDc20ZPM+WwfNsGeY6z031HBmwSJuIiIioDiZIRERERHUwQbIySqUSy5cvh1LJDTXNiefZcniuLYPn2TJ4ni3DGs4zi7SJiIiI6mAPEhEREVEdTJCIiIiI6mCCRERERFQHEyQiIiKiOpggWZl169YhPDwczs7OiI6OxuHDh6UOyaYcOnQIEydORFBQEGQyGXbu3FnreUEQ8PrrryMoKAguLi4YOXIkzp07V+sYlUqFBQsWwM/PD25ubpg0aRLS09Mt+FtYt5UrV+KOO+6Ah4cH/P398eCDD+LixYu1juF5Nt769evRt29fcaG8oUOH4ueffxaf5zk2j5UrV0Imk2HRokXiYzzXpvH6669DJpPV+goMDBSft7rzLJDV2Lp1q+Dk5CR8/PHHQlJSkrBw4ULBzc1NuH79utSh2Yzdu3cLy5YtE7777jsBgLBjx45az7/99tuCh4eH8N133wmJiYnCtGnThA4dOghFRUXiMfPmzRM6duwoxMbGCqdOnRJGjRolREVFCRqNxsK/jXUaO3as8Nlnnwlnz54VEhIShPHjxwuhoaFCSUmJeAzPs/F27dol/PTTT8LFixeFixcvCq+99prg5OQknD17VhAEnmNzOHHihNCpUyehb9++wsKFC8XHea5NY/ny5UKfPn2EjIwM8Ss7O1t83trOMxMkKzJo0CBh3rx5tR7r2bOn8Oqrr0oUkW2rmyDpdDohMDBQePvtt8XHKioqBC8vL2HDhg2CIAhCQUGB4OTkJGzdulU85saNG4JcLhf27NljsdhtSXZ2tgBAOHjwoCAIPM/m5O3tLXzyySc8x2ZQXFwsdOvWTYiNjRXuvvtuMUHiuTad5cuXC1FRUQ0+Z43nmUNsVkKtViMuLg5jxoyp9fiYMWNw9OhRiaJqW5KTk5GZmVnrHCuVStx9993iOY6Li0NlZWWtY4KCghAREcG/QyMKCwsBAD4+PgB4ns1Bq9Vi69atKC0txdChQ3mOzeC5557D+PHjce+999Z6nOfatC5fvoygoCCEh4dj+vTpuHbtGgDrPM/crNZK5OTkQKvVIiAgoNbjAQEByMzMlCiqtsVwHhs6x9evXxePUSgU8Pb2rncM/w71CYKAxYsXY/jw4YiIiADA82xKiYmJGDp0KCoqKuDu7o4dO3agd+/e4ocBz7FpbN26FadOncKff/5Z7zn+fzadwYMH44svvkD37t2RlZWFN998E8OGDcO5c+es8jwzQbIyMpms1veCINR7jIzTmnPMv0PDnn/+eZw5cwZHjhyp9xzPs/F69OiBhIQEFBQU4LvvvsOcOXNw8OBB8XmeY+OlpaVh4cKF2LdvH5ydnRs9jufaeDExMeK/IyMjMXToUHTp0gWff/45hgwZAsC6zjOH2KyEn58fHBwc6mXB2dnZ9TJqah3DbImmznFgYCDUajXy8/MbPYb0FixYgF27duHAgQMIDg4WH+d5Nh2FQoGuXbti4MCBWLlyJaKiorBmzRqeYxOKi4tDdnY2oqOj4ejoCEdHRxw8eBDvvfceHB0dxXPFc216bm5uiIyMxOXLl63y/zQTJCuhUCgQHR2N2NjYWo/HxsZi2LBhEkXVtoSHhyMwMLDWOVar1Th48KB4jqOjo+Hk5FTrmIyMDJw9e5Z/hyqCIOD555/H9u3b8euvvyI8PLzW8zzP5iMIAlQqFc+xCY0ePRqJiYlISEgQvwYOHIiZM2ciISEBnTt35rk2E5VKhfPnz6NDhw7W+X/a5GXf1GqGaf6ffvqpkJSUJCxatEhwc3MTUlJSpA7NZhQXFwvx8fFCfHy8AEBYtWqVEB8fLy6V8PbbbwteXl7C9u3bhcTEROHRRx9tcBppcHCwsH//fuHUqVPCPffcw+m6NTz77LOCl5eX8Ntvv9WarltWViYew/NsvKVLlwqHDh0SkpOThTNnzgivvfaaIJfLhX379gmCwHNsTjVnsQkCz7WpvPjii8Jvv/0mXLt2TTh27JgwYcIEwcPDQ/yMs7bzzATJynzwwQdCWFiYoFAohAEDBohTp6l5Dhw4IACo9zVnzhxBEPRTSZcvXy4EBgYKSqVSuOuuu4TExMRar1FeXi48//zzgo+Pj+Di4iJMmDBBSE1NleC3sU4NnV8AwmeffSYew/NsvCeffFJsC9q3by+MHj1aTI4EgefYnOomSDzXpmFY18jJyUkICgoSJk+eLJw7d0583trOs0wQBMH0/VJEREREtos1SERERER1MEEiIiIiqoMJEhEREVEdTJCIiIiI6mCCRERERFQHEyQiIiKiOpggEREREdXBBImIiIioDiZIRGTTfvvtN8hkMhQUFEjy/r/++it69uwJnU7X6DGvv/46+vXrJ36/ZMkSvPDCCxaIjohaiwkSEdmMkSNHYtGiRbUeGzZsGDIyMuDl5SVJTC+//DKWLVsGubz5zenLL7+Mzz77DMnJyWaMjIiMwQSJiGyaQqFAYGAgZDKZxd/76NGjuHz5MqZOndqin/P398eYMWOwYcMGM0VGRMZigkRENuHxxx/HwYMHsWbNGshkMshkMqSkpNQbYtu4cSPatWuHH3/8ET169ICrqyumTJmC0tJSfP755+jUqRO8vb2xYMECaLVa8fXVajVefvlldOzYEW5ubhg8eDB+++23JmPaunUrxowZA2dn51qPv/322wgICICHhwfmzp2LioqKej87adIkbNmyxejzQkTmwQSJiGzCmjVrMHToUDz99NPIyMhARkYGQkJCGjy2rKwM7733HrZu3Yo9e/bgt99+w+TJk7F7927s3r0bX375JT766CN8++234s888cQT+P3337F161acOXMGU6dOxbhx43D58uVGYzp06BAGDhxY67Gvv/4ay5cvxz//+U+cPHkSHTp0wLp16+r97KBBg5CWlobr16+38owQkTk5Sh0AEVFzeHl5QaFQwNXVFYGBgU0eW1lZifXr16NLly4AgClTpuDLL79EVlYW3N3d0bt3b4waNQoHDhzAtGnTcPXqVWzZsgXp6ekICgoCoC+k3rNnDz777DO89dZbDb5PSkqKeLzB6tWr8eSTT+Kpp54CALz55pvYv39/vV6kjh07iq8RFhbW8hNCRGbFHiQianNcXV3F5AgAAgIC0KlTJ7i7u9d6LDs7GwBw6tQpCIKA7t27w93dXfw6ePAgrl692uj7lJeX1xteO3/+PIYOHVrrsbrfA4CLiwsAfW8XEVkf9iARUZvj5ORU63uZTNbgY4ap+TqdDg4ODoiLi4ODg0Ot42omVXX5+fkhPz+/VTHm5eUBANq3b9+qnyci82KCREQ2Q6FQ1CqsNpX+/ftDq9UiOzsbI0aMaNHPJSUl1XqsV69eOHbsGB577DHxsWPHjtX72bNnz8LJyQl9+vRpfeBEZDYcYiMim9GpUyccP34cKSkpyMnJaXJxxpbo3r07Zs6cicceewzbt29HcnIy/vzzT/zrX//C7t27G/25sWPH4siRI7UeW7hwIf73v//hf//7Hy5duoTly5fj3Llz9X728OHDGDFihDjURkTWhQkSEdmMJUuWwMHBAb1790b79u2Rmppqstf+7LPP8Nhjj+HFF19Ejx49MGnSJBw/frzRmXIAMGvWLCQlJeHixYviY9OmTcPf//53vPLKK4iOjsb169fx7LPP1vvZLVu24OmnnzZZ/ERkWjJBEASpgyAislUvv/wyCgsL8eGHHzb7Z3766Se89NJLOHPmDBwdWelAZI3Yg0REZIRly5YhLCysRbVRpaWl+Oyzz5gcEVkx9iARERER1cEeJCIiIqI6mCARERER1cEEiYiIiKgOJkhEREREdTBBIiIiIqqDCRIRERFRHUyQiIiIiOpggkRERERUBxMkIiIiojr+H6/gkOi6ATIxAAAAAElFTkSuQmCC",
"text/plain": [
"Figure(PyObject <Figure size 640x480 with 2 Axes>)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"PyObject Text(27.87866577148438, 0.5, 'lake level (m)')"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## This parameter contains the hydraulic friction. For Manning roughness = 0.1 -> 530\n",
"const F1 = 130.0 # kg/m^(-8/3) <---\n",
"## Lake input \n",
"const Q_in = 100.0 # <--- [5, 100]\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 # <--- [100, 500]\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) <--- [100, 10000]\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",
"\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, 500day) # 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",
" sol[1,i] = max(sol[1,i], 1e-14) # trick to not let R-channel close completely as otherwise it will not re-open\n",
"end \n",
"\n",
"# plot solution\n",
"using PyPlot\n",
"subplot(2,1,1)\n",
"semilogy(ts/day, sol[1,:]);ylabel(\"Q (m^3/s)\") # use semilogy or plot\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.11.1",
"language": "julia",
"name": "julia-1.11"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.1"
}
},
"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