02_example_analytical.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2\n",
    "\n",
    "> This example will validate the code of the TEF package by comparing the results with the analytical example of Burchard et al. (2019) and Lorenz et al. (2019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyTEF.tef_core import *\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "import xarray as xr\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining the necessary parameters:\n",
    "\n",
    "Cross section: $ A = 10000 m^2$\n",
    "\n",
    "Tidal period: $ T = 12 h * 3600 \\frac{s}{h}$\n",
    "\n",
    "Tidal frequency: $\\omega = 2 * \\pi / T$ \n",
    "\n",
    "Residual velocity: $ u_r = -0.1 $\n",
    "\n",
    "Tidal velocity amplitude: $ u_a = 1.$\n",
    "\n",
    "Residual salinity: $s_r = 20. g/kg$ \n",
    "\n",
    "amplitude of salinity changes within one tidal cycle: $s_a = 10. g/kg$ \n",
    "\n",
    "phase lag between velocity and salt: $\\phi = \\arccos(-2\\frac{u_r s_r}{u_a s_a}) $\n",
    "\n",
    "N = 10000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 10000. \n",
    "T = 12*3600\n",
    "omega = 2*np.pi/T #\n",
    "u_r = -0.1 \n",
    "u_a = 1.\n",
    "s_r = 20. \n",
    "s_a = 10. \n",
    "phi = -np.arccos(-2*(u_r*s_r)/(u_a*s_a)) \n",
    "N = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u_vel_analytical(t):\n",
    "    return u_r+u_a*np.cos(2*np.pi*t)\n",
    "    \n",
    "def s_analytical(t):\n",
    "    return s_r + s_a*np.cos(2*np.pi*t+phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define analytical solution according to Appendix of Lorenz et al. (2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Q_analytical(S):\n",
    "    t_1 = -1/omega*(np.arccos((S-s_r)/s_a)+phi)\n",
    "    t_2 = 1/omega*(np.arccos((S-s_r)/s_a)-phi)\n",
    "    Q = A/T*(u_r*t_2+u_a/omega*np.sin(omega*t_2)-u_r*t_1-u_a/omega*np.sin(omega*t_1))\n",
    "    return(Q)\n",
    "    \n",
    "def Qs_analytical(S):\n",
    "    t_1 = -1/omega*(np.arccos((S-s_r)/s_a)+phi)\n",
    "    t_2 = 1/omega*(np.arccos((S-s_r)/s_a)-phi)\n",
    "    Qs = A/T*(u_r*s_r*(t_2-t_1)+u_a*s_r/omega*(np.sin(omega*t_2)-np.sin(omega*t_1)) + u_r*s_a/omega*(np.sin(omega*t_2 + phi)-np.sin(omega*t_1+phi))\\\n",
    "    +u_a*s_a/2*np.cos(phi)*(t_2-t_1 + 1/omega *np.sin(omega*t_2)*np.cos(omega*t_2) - 1/omega *np.sin(omega*t_1)*np.cos(omega*t_1))\\\n",
    "    -u_a*s_a*np.sin(phi)/2/omega*(np.sin(omega*t_2)**2 - np.sin(omega*t_1)**2))\n",
    "    return(Qs)    \n",
    "    \n",
    "def q_analytical(S):\n",
    "    q=-A/(T*omega)*1/(np.sqrt(s_a**2-(S-s_r)**2))*(-2*u_r-u_a*np.cos(np.arccos((S-s_r)/s_a)-phi)-u_a*np.cos(np.arccos((S-s_r)/s_a)+phi))\n",
    "    return(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretizing the analytical functions\n",
    "\n",
    "Create the discrete values for time, volume transport, salt, and salt transport for N discrete value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.float64(np.arange(0, 1, 1/N))                  # time array (0 to 1 , so relative to one Tidal cycle T)\n",
    "volume_transport = A*np.float64(u_vel_analytical(t))    # volume transport\n",
    "salt = np.float64(s_analytical(t))                    # salt / salinity\n",
    "salt_transport = A*np.float64(u_vel_analytical(t)*salt) # salt transport\n",
    "\n",
    "data = xr.Dataset({\"salt\": ([\"time\"], salt),                 \n",
    "                   \"volume_transport\": ([\"time\"], volume_transport),\n",
    "                   \"salt_transport\": ([\"time\"], salt_transport)},\n",
    "                    coords={\"time\": ([\"time\"], t)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot analytical functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAEWCAYAAACuSXe8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABE/ElEQVR4nO3dd3zUVfb/8dcJSK8qbUXAElpQwISiYkFFYRUL2EGxgq5f+1rWhqKusrs2XCuKFQsqCvLDtgoIShGkiwgWBJUmvUrI+f1xJ0tkI2Qgmc/M5P18PD6PZD7TTj5Kcubec881d0dEREQkVWREHYCIiIhIPJS8iIiISEpR8iIiIiIpRcmLiIiIpBQlLyIiIpJSlLyIiIhISlHyIiIiIilFyYuIiIiklLiTFzO7qSQCERERESkK21mHXTMbUvAm0MrdM0s0KhEREZE/ULYIj1nj7pfk3zCzJ0owHhEREZEdKsrIy37u/n2B23u6+4oSj0xERESkEDtNXv77QLO93X15CccjIiIiskPxFOwOKrEoRERERIoonuTFSiwKERERkSKKJ3kp2vySiIiISAnSyIuIiIiklHgKdlu4+6wSjkdERERkh4o88pKoxMXMBpnZUjObVeDcnmb2kZnNi32tWeC+v5nZfDOba2YnFDifbWYzY/cNMDOLnS9vZq/Hzk80s0aJ+LlERESkeMS1PYCZnWFmVWPf32ZmQ83skGKO6Xmg83bnbgY+jnX2/Th2GzNrDpwNZMWe87iZlYk95wmgN5AZO/Jf82JgpbsfCDwE9C/m+EVERKQExbu30e3uvtbMOgAnAC8QkoRi4+6fAts3wTsl9l7Evp5a4Pxr7r451khvPtDWzOoB1dx9vId5sRe3e07+a70JHJs/KiMiIiLJryjbAxS0Nfb1ROAJdx9mZncWb0iFquPuvwC4+y9mVjt2fh9gQoHHLYqd2xL7fvvz+c9ZGHutXDNbDewF/E8DPjPrTRi9oXLlytlNmzYtth9IRKQ0mDJlynJ3rxV1HJJe4k1efjKzp4DjgP5mVp5d2Jm6GBU2YuI7OL+j5/zvSfengacBcnJyfPLkybsSo4hIqWVmC6KOQdJPvInHmcAHQGd3XwXsCdxQ3EEVYklsKojY16Wx84uAfQs8rj7wc+x8/ULO/+45ZlYWqM7/TlOJiIhIkooreXH3De4+1N3nxW7/4u4flkxovzMc6BX7vhcwrMD5s2MriPYjFOZOik0xrTWz9rF6lvO3e07+a50OfOJFXS8uIiIikYt32qjEmdmrwNHA3ma2COgL3A8MMbOLgR+BMwDcfbaZDQG+AnKBK9w9vy7ncsLKpYrAe7ED4FngJTObTxhxObsocc2YAQccABUqQPnyUL061K4NdeqEY7/9IDMTGjcO94mIpKuNG2HePPjhh23H4sWwcmU4Vq2CTZsgNzfaOCV9FblJXWlXq1aOd+48mU2bwj/K1athyZJwrF79+8fWqQOtW0ObNuFo2zacExFJNZs2wRdfwOefw/TpMG0azJ0LeXnbHlOxItSrBzVrbjsqVoQ99oBnn7Up7p5T3HFNmTKldtmyZZ8BWhBt7aWUjDxgVm5u7iXZ2dlLt79zpyMvZnYZkE3or9IT+H/uXqzLo1NBw4bw0kuF37dxI3z/PXzzTTjmzIEpU+DDD7f9A8/KgmOPheOOg6OOgmrVEhe7iEhR5eXBpEkwciR8+ilMmACbN4f7GjSAli2he3do0SKMODdqBLVqwR81nHj22ZKJs2zZss/UrVu3Wa1atVZmZGToU3iaycvLs2XLljVfvHjxM8DJ299flGmjY4CzgLHu3sHMnizuIFNdxYrQvHk4Clq/HqZODZ9YPv4Ynn4aBgwIn0aOOy78AjjlFNh772jiFhGBMLry4YcwfDiMGBFGlDMywgjyFVfAkUfC4Ycn3e+qFkpc0ldGRobXqlVr9eLFi1sUdn9Rkpdf3d3NLL8T7ebiCy+9Va4MHTqE48Ybwy+I8ePDL4ehQ+GSS6BPH+jYES64ALp1C4mQiEhJy8uDcePg5ZfhjTdCnUr16tClC3TtGr7WrLnTl4lShhKX9Bb771volGBRkpdHANz93djtocUUV6lToUJIVDp2hH/9K8wdv/UWvPIK9OwZfnGcfTZcfHGolRERKW5LlsDAgfDMM7BgQfiQ1a0b9OgBxxwTRoZFkt1Oi5zc/evtbo8puXBKD7MwJHvPPTB/PnzyCZx8Mrz4YijwPfRQeP112LIl6khFJNW5h1HfHj1g333h9tvhwANDHd+SJeH3zgknKHEpSXPnzi2XmZmZlaj3W758eZn7778/JTob33zzzXXjfU68GzPmmNnbZvalmc2I7do8I943ld/LyAijMS++CL/8Ao88AsuWhVGY/feH/v1hzZqooxSRVJOXF+pYDj0UDjssTFlffjl8/TX85z9hxLdy5aijlJLw66+/lnn22WdrF3ZfbpKsYc/Ly2Pr1q0MGDCgXrzPjXd52WDgOaA70BU4KfZVikn16nDVVWEp4rBhoXfMzTeH1U79+oV5aRGRHcnNDdPRLVuGRQFLlsC//w0//RQ+HDVpEnWEqe/yyy/fp+DIxnXXXfenvn371snLy6NPnz71MzMzsxo3btx84MCB/1M5NGDAgL3OP//8Bvm3O3bseOCIESOqAlSqVKn15Zdfvk9WVlazww47rPGoUaMqtW3btkn9+vUPGjx4cHUIyUefPn3qt2jRolnjxo2b//Of//yfUurrr7++/sKFC8s3bdq0eZ8+feqPGDGiart27Rp37dp1vyZNmmQBHHfccQdkZWU1O/DAA7P+9a9//fc1KlWq1PrKK6/cp0mTJs1btmzZdOHChWUBBg0aVDMzMzOrSZMmzXNycprk/yzHHnvsAUcccURmo0aNWlx//fX/TUTuvPPOOpmZmVmZmZlZ/fr1qw1hBGr//ffP6tmzZ4OsrKzmZ511VqPNmzdnNG3atPnJJ5+8X1Gvf7xN6pa5+/A4nyO7oEyZMI108skweXJIXPr2hQcfhKuvhmuvhRo1oo5SRJJJXh4MGRKmhebPDysgX3opjOKWTbqWpMXnoovYd9YsKhXna7ZowYZBg8ImvoXp2bPnimuuuabBzTffvAxg2LBhNd9///15L774Yo2ZM2dWnDNnzuxffvmlbNu2bZsdf/zx64r6vhs3bszo2LHj2ieeeOKnTp06HXDbbbftM3bs2G++/PLLChdeeOF+PXr0WP3www/vXb169a2zZs2as3HjRmvTpk3Trl27rmnatOlv+a/zwAMPLDrppJMqfv31118BjBgxouqMGTMqT506dXb+4wYPHvxDnTp1tq5bt85at27dvGfPnivr1q27dePGjRmHHnroukcfffSnyy67rP6jjz5a6x//+Mcv999/f70PP/zwm/3222/L8uXLy+S/14wZMyrPnDlzdpUqVfJat27d/JRTTlltZrzyyit7TZkyZY67k52d3ezYY49du/fee2/94YcfKgwcOPCHl19++UeASpUq1cyPs6jiHXnpa2bPmNk5ZtYt/4jzNSROOTlh6PfLL8P0Ur9+odvvQw9t678gIqWXe1jqnJMD55wTVi0OHQozZ4apoXROXKJy+OGHb/z111/L/vDDD3uMHz++YvXq1bdmZmb+Nnbs2KpnnnnmirJly7LvvvvmtmvXbt24ceOKnFjtsccefvrpp68ByMrK2tihQ4e15cuX97Zt22786aefygH85z//qTZkyJC9mjZt2rx169bNVq5cWfarr76qsLPXPvjgg9cXTHD69+9fp0mTJs2zs7ObLV68eI/Zs2dXyI/h7LPPXg2QnZ29fsGCBeUAcnJy1vXo0aPRAw88sHfBqacOHTqsqVu37tYqVar4iSeeuHL06NFVRo8eXeXPf/7zqmrVquVVr14978QTT1w5atSoqgD16tX77dhjj11f1GtSmHj/l74QaArsQeh+B2FHZq1ASoDWreHtt0PvmJtuguuug0cfhfvugzPP/OMmUSKSvr78Em64IRT9N2oURlrOOSeM3pYWOxohKUldu3Zd+fLLL9dcvHjxHt27d18BUJSu9WXLlvW8Ai2KN2/enFHwvoyMcDMjI4Py5cs7QJkyZdi6davF3sMeeOCBH7t37x5XNWSlSpX++6YjRoyoOmbMmKqTJ0/+umrVqnlt27ZtsnHjxoztYyhbtiy5ubkG8Morr/z4ySefVB4+fHj1Vq1aZU2bNm02gG33x8fMdngdCsaxq+IdeWnp7jnu3svdL4wdF+1uEBKf1q3Dp6z334cqVcKQcPv2oYW3iJQOy5eHPlE5OWGEZcCAUIjbs2fpSlyidN55561466239hwxYkTNnj17rgQ46qij1r755pt75ubm8vPPP5edNGlSlSOOOOJ3owwHHHDAb7Nnz660detW5s+fv8eMGTPiKpvu1KnT6ieeeKLW5s2bDWDGjBnl16xZ87u/59WrV9+6fv36P/wbv2rVqjLVq1ffWrVq1bypU6dWmD59+k5jmD17dvljjjlm/cMPP/xzzZo1c7/77rtyAOPGjau2ZMmSMuvWrbORI0fWOOqoo9Ydc8wx60aOHFlj7dq1GWvWrMkYOXJkzY4dO64t7HXLli3r+T9LUcWbvEwws+Y7f5gkwgknhFGY556DhQuhXTu47DL49deoIxORkpKbG4pvMzNh0CC45pqwSeKVV4ZNYyVxcnJyNq1fvz6jTp06vzVs2HALwHnnnbcqKytrY7NmzbKOPvroxnfdddeiBg0a/G55T6dOndbtu+++m5s0aZJ19dVX79u8efMN8bzvtddeu7xp06abDjrooGaZmZlZl156acMtW7b87o9/3bp1t2ZnZ6/LzMzM6tOnT/3tX6N79+6rc3NzrXHjxs1vueWWP7Vs2XKn0zjXXntt/caNGzfPzMzMat++/dr27dtvjF2HdWedddZ+LVq0yOratevKI488ckOHDh02nHvuub8ecsghzbKzs5udd955yw4//PCNhb1ujx49ljVr1iyugt24NmY0sznAAcD3hE67Bri7H1zkF0lROTk5Pnny5KjD+ENr1sCdd4ZPXzVqhOXVF14YlmGLSHr4/PPwAWXmzLBX2oAB/7stSbIxK5mNGadPn/5Dy5Ytlxf360p8BgwYsNfkyZMrv/jiiz+WxOtPnz5975YtWzba/ny8f9o6A5nA8WipdFKpVi2sRJo6Nfwyu+SSsC3BV3HVb4tIMlq7NoysdOgQ2iW8+SZ89FHyJy4iJSWugl13X1BSgUjxOOggGDMmFO1dd12oj+nbNxT0qXumSOp5771Q27JoUUhg7rkHqlaNOiqR4KqrrvoVSHixgiYV0pAZnH9+GHU59VS49daw5cC0aVFHJiJFtXx5KL79859DsvLZZ6HBnBKX/8rLy8vTGss0FvvvW+jKJCUvaax27bA/0ltvhW0H2rQJzat++23nzxWR6Py//wctWoSGc337huXQhx4adVRJZ9ayZcuqK4FJT3l5ebZs2bLqwKzC7i/StJGZHQpM8HiqeyVpdOsGRx8duvLec08Yhn75ZWjaNOrIRKSgdevgr3+Fp54KU8AffggHp/1yiF2Tm5t7yeLFi59ZvHhxC/RBPB3lAbNyc3MvKezOIq02MrMngbbAN8D7wPvuvrg4o0x2yb7aqKiGDoVLL4WNG0OBb58+am4nkgwmTIDzzoNvvw01av36pcfS55JabSSlW5GyVXe/zN0PAe4EagLPm9l4M/u7mR1pZmqJlCK6dQvLLI84Iuwue/LJsHRp1FGJlF5btoTp3MMPD9+PHh1aHaRD4iJSUuIaanP3r939IXfvDBwDjAPOACaWRHBSMv70pzB19MgjYbnlQQeF2yKSWAsXhinde+4Joy4zZsCRR0YdlUjy2+V5Qnff6O4j3f1KDQmmnowMuOoqmDIF6tYNKxr+9rfQvVNESt6770KrViFhefVVeP750K9JRHZORU6lXFZWmGvv3Rvuvz/sWr1oUdRRiaSv336D668PU7YNG4aVRGefHXVUIqlFyYtQsWJY3TB4cOgF07p12PRRRIrX99+HerMHH4T/+7/Q7j8zM+qoRFJPXMmLmfUvyjlJTeeeC5Mnh5qYLl00jSRSnIYNCx8M5s4N7f0ffRQqVIg6KpHUFO/IS6dCznUpjkAkOTRp8vtppOOOgyVLoo5KJHVt3Qq33Ra6XWdmhv3HunePOiqR1Fak5MXMLjezmUBTM5tR4PgemFGyIUqi5U8jvfQSTJoE2dkwUevJROL2669w4olw771hs9SxY2G//aKOSiT1FXXk5RXgFGACYRfp/CPb3XuWUGwSsZ49w5z8HnuE5ZsDB0YdkUjqmDoVcnJg1Ch4+unw70fTRCLFo6hN6la7+3dADXdfUOBYUcLxScRatQrLqY8+Okwl9e4NmzdHHZVIcnvhBTjssFAzNnZs6GotIsUn3pqX8WbWpkQikaS1554wcmQo4B04MIzCaDm1yP/67Te44gq44IKwkeKUKWFHdxEpXvEmLx2BCWb2bazmZaaZqealFChTBv7+97BD9VdfwSGHwKefRh2VSPJYsgSOOQYefzxsrvjhh2FndxEpfkXaVboArSwq5bp1g2bNwsqJY4+Fxx4LU0kipdm0aaHp3PLl8NprcNZZUUckkt7i3dtoAVCDbQW7NWLnpBRp1iysPjr22LAr9VVXqR+MlF5Dh4ZNFd1h3DglLiKJEG+TuquBwUDt2PGymV1ZEoFJcqtRA0aMgGuvDc22unSBlSujjkokcdzh7rtDz5aDD4YvvgjTqSJS8uKdNroYaOfu6+G/3XXHA48Wd2CS/MqWDW3OW7SAyy6Ddu1g+HBo2jTqyERK1oYNcNFF8PrrYTfop5/WMmiRRIq3YNeArQVub42dk1LsootCL4vVq6F9e+2LJOlt0aKw4m7IEOjfPyyLVuIikljxJi/PARPN7E4zu5PQtO7ZYo8qQcyss5nNNbP5ZnZz1PGkssMPD914GzUKHUUfeigMq4ukk4kToU2bsD/R8OFw441g+vgmknDxFuw+CFwIrABWAhe6+8MlEFeJM7MywGOEFVTNgXPMrHm0UaW2hg1DweKpp8J118HFF6uhnaSPl1+Go46CSpXC/l8nnRR1RCKlV7wjL7j7l+4+wN0fcfepJRFUgrQF5rv7d+7+G/AaYQsE2Q1VqsAbb8Add8Bzz4WNHZcujToqkV2XlxcaNJ53XpgWnTgRsrKijkqkdIt3tVEFM7vOzIaa2Vtmdq2Zpeps7z7AwgK3F8XO/ZeZ9TazyWY2edmyZQkNLpVlZMBdd4VixsmTQ4fRmTOjjkokfuvWwWmnhR3W+/QJjef23jvqqEQk3pGXF4EswuqifwPNgJeKO6gEKWym+ndVGu7+tLvnuHtOrVq1EhRW+jjzzLCvy5YtYZ+X4cOjjkik6H78MdRyjRgBAwbAE09AuXJRRyUiEH/y0sTdL3b3UbGjN9C4JAJLgEXAvgVu1wd+jiiWtJWTE/pf5Hfl7d9fhbyS/CZODCOGP/wQ9vW68koV5ookk3iTl6lm1j7/hpm1Az4r3pAS5gsg08z2M7NywNmAxgZKwJ/+BGPGhM6jN98cNq3btCnqqEQK99proTC3cmUYPx5OOCHqiERke/EmL+2Az83sBzP7gdCg7qhU3KDR3XOB/wM+AOYAQ9x9drRRpa+KFeGVV6BfP3jxxbCB3ZIlUUclso079O0L55wTRl0mToTmWn8okpTi7bDbuUSiiIi7jwRGRh1HaWEGt98e/iCcf37ol/Huu9CyZdSRSWm3cWMYERwyJHx98kkoXz7qqETkj8Q78tIWWBHbjPE84CFgL3dfoA0apai6dw/9YNxDQeQ770QdkZRmv/wSponeeAP+8Q8YNEiJi0iyizd5ud3d15pZB+AE4AXgieIPS9Jd69ahI2+LFmEp6n33qZBXEm/atDBF9NVX8PbbcMMNKswVSQXxJi/5+xqdCDzh7sMALR6UXVKvHoweDT16wC23hCZgKuSVRBk2DDp0CMnKuHFwilpUiqSMeJOXn8zsKeAsYKSZld+F1xD5rwoV4KWX4N57YfBgOPpoWLw46qgknbmH6aHTTgudcidNglatoo5KROIRb+JxJmF1zgnuvgqoCdxQ3EFJ6WIWRl6GDg2deNu0gampvPGEJK3Nm8Mu6DfdFJoojh4NdetGHZWIxGtXpo0qAGeY2R1Ab6D9jp8iUjSnnQaffRaSmQ4d4K23oo5I0sny5dCpEzz/fFgS/eqrYQm/iKSeeJOXYcDJQC6wvsAhUixatQodeVu2hNNPh7vvViGv7L5Zs0Jh7hdfhKTlzjtVmCuSyuLt81Lf3dOq14sknzp14JNPoHfvsDv1V1+F5av6lCy74u23QzF4tWphmqhdu6gjEpHdFe/Iy+dmdlCJRCJSQIUK8MILYTff118PfTh+1s5TEoe8vNDRuVu3UJg7ebISF5F0EW/y0gGYYmZzzWxGKm4LIKnDLBRWvvNOGH1p0yb8ARLZmXXrQkFu376hm/OYMWGPLRFJD/FOG3UpkShEduDkk+Hzz8PXI48MBZdnnhl1VJKsvv8+9GyZPRsefBCuuUb1LSLpJq6Rl9gWAGuAOkDDAodIiTr44NCP45BDwu7Ud94ZpgVECho1KozQLVwI770H116rxEUkHcWVvJjZJcCnhF4vd8W+3ln8YYn8r9q14eOPw8Z5d90FZ58N67XWTQgr0v7977AUunbtkOgef3zUUYlISYm35uVqoA2wwN07Aq2BZcUelcgfKF8+rDz617/gzTehfXuYNy/qqCRKGzZAr15w5ZXQpQtMmACZmVFHJSIlKd7kZZO7bwIws/Lu/jXQpPjDEvljZnD99fDBB2FH4Jwc7UxdWs2fD4ceCi+/HKYShw0LS6JFJL3Fm7wsMrMawDvAR2Y2DNACVolEp07w5ZfQpEnoznvzzZCbG3VUkijDh4fEdeFCGDkyrCzK0E5rIqVCkf+pm5kBV7n7Kne/E7gdeBY4tWRCE9m5Bg1g7Fjo0wf69w91DkuWRB2VlKStW+HWW8OKogMOCAlsZ7XOFClVipy8uLsTRlzyb49x9+Hu/ltJBCZSVOXLw5NPhiXU48eHFUnjx0cdlZSEZctCovL3v8PFF4e9sBo1ijoqEUm0eAdZJ5hZmxKJRGQ39eoVkpYKFUI/mAcf1HLqdDJmTNj7auxYeOaZcFSoEHVUIhKFeJOXjsB4M/tWHXYlGbVqFbrwnnRSKOo96SRYujTqqGR3bN0ainGPOQYqVw4J6sUXRx2ViERJHXYl7dSsCUOHwuOPhwSmVauwGuWYY6KOTOK1aBH06AGffho2V3zsMahaNeqoRCRq8Y68/MXdFxQ8gL+URGAiu8MMrrgCJk6E6tXhuONCkadWI6WOESNC4jllStik88UXlbiISBBv8tKpkHMajZGk1bJlmEa66KJQ5HnUUWHvG0leGzbAVVdB165Qv35IXs4/P+qoRCSZFCl5MbPLzWwm0CRW65J/fA+o5kWSWuXKobjz1Vdh1qywT9LAgaGlvCSXSZOgdWt49NGQwEyYEPr4iIgUVNSRl1eArsDw2Nf8I9vde5ZQbCLF6uyzYeZMaNsWeveGE0+En9ViMSls2QJ33AGHHQYbN4Y9rB55RKuJRKRwRUpe3H21u//g7udsV/OyoqQDFClODRrARx+FT/ajR0OLFmFERqMw0fnqq7BH1d13h+LcGTNUXC0iO6Zm2lLqZGTA//0fTJsWpiTOPRfOOCPskySJ89tv0K9fmCb68Ud4661QmFujRtSRiUiyU/IipVbjxqHh2X33hZUtzZqFTr1qbFfy8jsh9+0b9qWaNQu6dYs6KhFJFUpepFQrWzZs6DhzZvhjevnlcMQRMHt21JGlpzVrwqjX4YeH70eMgNdegzp1oo5MRFJJXMmLBT3N7I7Y7QZm1rZkQhNJnMzMUCT6/PMwd26Yyrj1Vli/PurI0oM7DB4cRrcefxyuvDIkiCeeGHVkIpKK4h15eRw4FDgndnst8FixRiQSEbOwP9LXX8M554S+MI0bh+68mkradV9+GUazevaEevXClNEjj6jhnIjsuniTl3bufgWwCcDdVwLlij0qkQjtvXcoHP3sM/jTn0Jb+sMOC916peiWLoU+fSAnB775JvTamTQJ2rWLOjIRSXXxJi9bzKwM4ABmVgvQZ1JJS/kJy3PPwYIFYTlvjx4wf37UkSW3NWvCRooHHACDBsE114Tk5eKLw0ovEZHdFe+vkgHA20AdM7sX+Ay4r9ijEkkSGRlwwQXhj+/f/gZvvw1Nm4YmdwsXRh1dctm0CR5+OCQtd90FnTuHQugHH9TyZxEpXnElL+4+GLgR+DvwM3Cyuw8picBEkknVqqEG5ttvw4qk55+HAw8Mowo//RR1dNHasCE0/WvcGK69NmymOGkSvPFGSPRERIpbvKuNcggjLZcAlwFDzEx7G0mpUa9e+EM9b14oQP33v2G//cLGj3PmRB1dYq1eHXrkNGoU9iFq2DB0L/7oI2jTJuroRCSdxTttNBh4DugGnMS2PY6KhZmdYWazzSwvligVvO9vZjbfzOaa2QkFzmeb2czYfQPMzGLny5vZ67HzE82sUYHn9DKzebGjV3HFL6VHw4bw7LMhienTJ/Qqad4cTj0VxoxJ7+0G5s0LIywNGsAtt0B2Nnz6aWj4d9xxUUcnIqVBvMnLMncf7u7fF9zjqBjjmUVIjD4teNLMmgNnA1lAZ+DxWOEwwBNAbyAzdnSOnb8YWOnuBwIPAf1jr7Un0BdoB7QF+ppZzWL8GaQU2W+/MBKzYEHoFjt2LBx9dEhkHn4YVqTJ7l9bt8K774Y6lsaNw4jTn/8MU6bAe++FpdAiIokSb/LS18yeMbNzzKxb/lFcwbj7HHefW8hdpwCvuftmd/8emA+0NbN6QDV3H+/uDrwInFrgOS/Evn8TODY2KnMC8JG7r4gt9f6IbQmPyC6pVSussFm4MNTD1KgRRif22SesUBoxIuzlk0rcYepUuP56qF8fTj45FOD26xd+zldfDV2JRUQSrWycj78QaArswbYl0g4MLc6gCrEPMKHA7UWxc1ti329/Pv85CwHcPdfMVgN7FTxfyHN+x8x6E0Z1aNCgwW7/EJL+KlUKje569YLp0+Gpp8KU0iuvwJ57QvfucOaZYaSifPmoo/1f+QnLu++GgtvZs2GPPUIn3PPOg65dw20RkSjFm7y0dPeDducNzew/QN1C7rrV3Yf90dMKOec7OL+rz/n9SfengacBcnJy0riKQUpCy5ahFf7DD8OHH4aRildegYEDoXLlUB/SpQscf3woerXC/s9MgOXLw3TXRx+FpGXRohDLYYfBE0+EHbf32iua2EREChNv8jLBzJq7+1e7+obuvislfYuAfQvcrk9Yqr0o9v325ws+Z5GZlQWqAyti54/e7jmjdyEmkSIpVw5OOikcGzbAJ5/AyJHhGBZL1+vVC8nC4YeHjrRZWWGkprht2gRffQXTpsHkyaHQNn8TykqV4IQTwrTQiSdC7drF//4iIsUh3uSlA9DLzL4HNhNGMdzdDy72yH5vOPCKmT0I/IlQmDvJ3bea2Vozaw9MBM4HHi3wnF7AeOB04BN3dzP7APh7gSLd44G/lXD8IkBIEPITGfewj9KoUfD552E7grfe2vbYevVCEnPAAaHmpH79UENTvXroO1O1ang997D3Ul5eSI5WrgzHihVhFOX77+GHH+C778Kmk7m54fWrVAnJ0rnnwlFHheXN5bTZh4ikAPM41nSaWcPCzhfXiiMzO42QfNQCVgHT3P2E2H23AhcBucA17v5e7HwO8DxQEXgPuDKWpFQAXgJaE0Zcznb372LPuQi4Jfa297r7czuLLScnxydPnlwcP6bIH/r551ArM2tWGBGZPTusZFq2bNdfs2rVsCqqUSM46KAwndWqVUiK1K5fSpqZTXH3nJ0/UqTo4kpeSjMlLxKlTZtCYvPTT6E53Nq14diwISQgZcqErxUqQM2a24599gnTT1HV04goeZGSENe0kZndUdh5d+9XPOGISGEqVID99w+HiEhpF2/Ny/oC31cgdNktZU3RRUREJEpxJS/u/kDB22b2L0JhrIiIiEhC7FbNS2zFziR3zyy+kJKTma0FCuv+WxrtDSyPOogkoWuxja7FNroW2zRx96pRByHpJd6al5lsa+hWhrAq6O7iDipJzVXRWWBmk3UtAl2LbXQtttG12MbMtNJBil28NS8nsS15yQWWuPvW4g1JRERE5I8VKXmJTZkUOr9kZrh7tWKNSkREROQPFCl50XwlENvjSABdi4J0LbbRtdhG12IbXQspdmpSJyIiIikl3poXzKwlcETs5lh3n168IYmIiIj8sbh2NjGzq4HBQO3Y8bKZXVkSgYmIiIgUJt6NGWcAh7r7+tjtysD4BOwqLSIiIgLEOfICGFBwafTW2DkRERGRhIi35mUQMNHM3iYkLacAzxZ7VCIiIiJ/IO6CXeD/gNaE5OVCd59avCGJiIiI/LF4p42qAU8CZxA67P5c7BGJiIiI7MAu9Xkxs4OBs4DuwCJ3P664AxMREREpTLwjL/mWAouBXwlLpkVEREQSIt4+L5eb2WjgY8KW75dqmbSIiIgkUrwFuw2Ba9x9WgnEIiIiIrJT2ttIREREUsqu1ryIiIiIRELJi4iIiKSUXWlSVyrtvffe3qhRo6jDEBFJKVOmTFnu7rWijkPSS1onL2ZWAfgUKE/4Wd90975mtifwOtAI+AE4091X7ui1GjVqxOTJk0s2YBGRNGNmC6KOQdJPuk8bbQaOcfeWQCugs5m1B24GPnb3TMKy75ujC1FERETikdYjLx6WUq2L3dwjdjhhQ8mjY+dfAEYDNyU4vOSycSP8+ms41qyB3NxtR9myUKXKtqN2bahcOeqIRUSklErr5AXAzMoAU4ADgcfcfaKZ1XH3XwDc/RczK7RLsJn1BnoDNGjQIFEhl5w1a2DaNJg+HebPh+++g2+/hQULYMOG+F6rWjXYZ59wNG0KzZpB8+bQogXsvXeJhC8iIgKlqM+LmdUA3gauBMa5e40C961095o7en5OTo6nVM2LO8ybB6NGwZgx8MUXIWHJV6UKHHAA7L8/NGoEtWrBXnuFo1o12GOPMOJSpkwYfVm/HtatC8eSJfDTT/Dzz/Djj/D117B27bbX3n9/OPTQcHToAAcfDGYJvwQiEj0zm+LuOcX9ulOmTKldtmzZZ4AWpEcJRB4wKzc395Ls7OylUQeT7NJ+5CWfu6+KbW3QGVhiZvVioy71CHs1pb516+CDD2DYMPj445BcANSrFxKJCy6A1q3DUbdu8SUU7uG9vvoqjOxMmACffAKDB4f769aFE06Azp3h+ONhzz2L531FpNQqW7bsM3Xr1m1Wq1atlRkZGSn/KTwvL8+WLVvWfPHixc8AJ0cdT7JL6+TFzGoBW2KJS0XgOKA/MBzoBdwf+zosuih305o18NZbMHQofPQRbN4ckoNOnaBjx3BkZpbsyIfZtimkTp3COfcwKjNqFLz/PgwfDi+8EEZzOnWCs86CU06BGjVKLi4RSWct0iVxAcjIyPBatWqtXrx4cYuoY0kFaZ28APWAF2J1LxnAEHcfYWbjgSFmdjHwI3BGlEHGbevWMLLywgvw9tuh2LZRI7j8cjj1VDj88JAkRMkMGjYMoz0XXBBi/uKLkGQNGRLOlSsHJ54Il14aRmTKlIk2ZhFJJRnpkrjki/086TAFVuLSOnlx9xlA60LO/wocm/iIdtOyZTBwIDzxBCxaFEYtevUKR7t2yV1XUqYMtG8fjv79QyLz2mvw8sshAWvQAC66CC6+GOrXjzpaERFJYsrwUsHUqXDhhbDvvnDrrWF1z5Ah8MsvIZFp3z65E5ftmUHbtvDggyEJe+MNaNIE7rwT9tsPevYMP7OISBro3r17o+eee64mQL9+/WqvXbtWf3t3ky5gMhs7NkynHHJISFYuughmzw61LWecARUqRB3h7itXDk4/HT78MCzbvvLKUHB8yCFw3HGhXqaUrIgTkfT31FNP1Vm3bp3+9u6mtJ42SknuoZ7l7rvh009DQ7j774c+fdK/uHX//cNozB13wNNPw4AB0KVLmBLr1y8U+qbSCJOIJMZFF+3LrFmVivU1W7TYwKBBC//o7jVr1mScfPLJ+//yyy/l8vLy7MYbb/x57ty5Fd5///0amzdvzsjJyVk3ePDgBRkZ2/KUe+65p/bSpUv3OOqooxrXrFkzd+LEid8Ua8yliLK/ZPL553DkkeGP9LffwiOPwPffw003pX/iUlCNGnDjjaGJ3sCBYXrshBPgiCPCEmwRkYgNHTq0Wt26dbfMnTv3q3nz5s3u1q3bmhtuuGHprFmz5sybN2/2xo0bM1577bXqBZ9z2223La1du/aWMWPGfKPEZfdo5CUZfP013HJLKFytWxcefzxMEZUvH3Vk0SpXDi65BM4/HwYNgnvugWOPDcndAw/AQQdFHaGIJIMdjJCUlEMOOWTjrbfeuu/ll1++zymnnLK6c+fO655//vkaDz74YN1NmzZlrFq1qmzz5s03AqsTHVtpoJGXKC1dCpddFlrq/+c/Yapo/vyw5Lm0Jy4FlSsXrtP8+fDQQzB5MrRqBb17w+LFUUcnIqXQwQcfvPnLL7/86qCDDtp466237vPXv/613vXXX99w6NCh337zzTdf9ezZc/mmTZv0N7aE6MJGITcXHn0UGjeGZ5+Fv/wlTBPddps2PNyRChXgmmtCEnPVVfDcc6EB3333heZ8IiIJ8sMPP+xRtWrVvL/85S8rrrnmmiXTpk2rBFC3bt3c1atXZ7z77ruFbjlTuXLlratXr9bf3t2kC5ho48ZBdnb449uuHcyaFQpTa9WKOrLUseeeYQRm9uwwjXTLLdCyZejmKyKSAFOmTKnYqlWrZk2bNm3ev3//enfccccvPXr0WNa8efOsLl26HNiyZcv1hT2vV69ey7t06ZLZrl27xomOOZ2Umo0Zd9dub8y4bBlcfz289FLo1/Lww3DaaVo9Uxzefx+uuCIU+J53HvzrX2GVlohErqQ2Zpw+ffoPLVu2XF7crxu16dOn792yZctGUceR7DTyUtLc4dVXoXnz0FH2lltgzhzo1k2JS3Hp3DmMYN12W7jGTZqEVUpKzEVE0pKSl5L0009h88Fzzw09TKZOhXvvVV1LSahYMRQ8z5ixrZi3UydYsCDqyEREpJgpeSkJ7vDMM2G05T//Cct6P/8csrKijiz9NW0aesE8+SRMnBhWcj31lEZhRNJPXl5eXloNX8d+nryo40gFSZ28mNlNUccQt19+gT//OeyUfMghYSTguuu0Y3IimYWOxLNmhaLoyy4L2yxoFEYkncxatmxZ9XRJYPLy8mzZsmXVgVlRx5IKkqpJnZkNKXgTaAX0jyaaXfD22yFp2bAB/v3v0K8lI6nzw/TWsGHYB+qpp+CGG0JTu8ceCxs/qt5IJKXl5uZesnjx4mcWL17cgiT/IF5EecCs3NzcS6IOJBUkVfICrHH3//6HM7MnogymyNauDf1HBg0Koy2DB4fpC4meWRh56dw5dOo9/3x4772wG3f16jt/vogkpezs7KXAyVHHIdFItmz13u1u3xpJFPEYPz4UiD7/fFhJNH68Epdk1KhR6ANz991hh+5WrUIdkoiIpJykSl7c/XsAM9s7dntFtBHtQF5e2O35iCPC92PGhJVE5cpFHZn8kTJlwnLqcePCiMwRR8Cdd4aOxyIikjKSKnkpYFDUAezQsmVw4onwt7/B6afDtGnQoUPUUUlRtW8f/pv16AF33QUdO4Zl7SIikhKSNXlJ3mrKceOgdeswBfHkk6EBnWonUk+1avDii/Dyy6H/TuvWobhXRESSXrImL8nZlOP+++Hoo0NDtAkTwnJcrVpJbT16wBdfhL2lTjghjMRs3Rp1VCIisgPJmrwkX0Ywb16YJureHaZMCQWfkh6aNYNJk8IS6jvvhC5dwtSgiOyer76KOgJJU8mavPwt6gD+x9q18PjjYe+catWijkaKW+XK8MIL8PTT8OmnYRpp3LiooxJJXW++CW3bRh2FpKmk3VXazK4r5PRqYIq7T0twOORkZfnk2bMT/bYShalT4YwzQkfehx4KO1ZrelCkaLZuDav67r8f2rfHJkwokV2lpXRL1pEXgBzgMmCf2NEbOBoYaGY3JjyaihUT/pYSkdatw9Rg585w5ZVw8cWwaVPUUYkkv19/Dduj3H9/qAkcPTrqiCRNJXPyshdwiLtf7+7XE5KZWsCRwAVRBialQPXqMGwY3H47PPccHHkkLFoUdVQiyWvaNMjJCQnLwIFhNWb58lFHJWkqmZOXBsBvBW5vARq6+0ZgczQhSamSkQH9+oU9q+bMgexs1cGIFOaVV+Cww2DLllAzdom255GSlczJyyvABDPra2Z9gc+AV82sMqASdkmcU0+FiRPDaEzHjqFwO0lrxUQSKjcXrrsutBxo0yZMt7ZrF3VUUgokbfLi7ncDlwKrCIW6l7l7P3df7+49Ig1OSp/mzcNy6hNOCAW8l1wCmzUAKKXYqlVw0kmhqP3KK+E//4E6daKOSkqJZNtV+r/MLNvdpwBTCpzr6u7vRhiWlGY1asDw4dC3L9xzD3z9dZhSql076shEEmv+/JC4fPstPPNMKGoXSaCkHXkhrCo6KP+GmZ0D3BZhPCKhDubuu+GNN8KS6jZtYMaMqKMSSZxRo0L/luXLw2iLEheJQDInL6cDL5hZMzO7FPgLcHzEMYkEp58OY8eGnhaHHRZWJomku6eeguOPh3r1wjTqUUdFHZGUUkmbvLj7d8DZwFuEROZ4d18dbVQiBWRnh32RmjeH006D++5TIa+kp9xcuPpquOwy6NQJPv8c9t8/6qikFEu6mhczm8nvN2bcEygDTDQz3P3gaCITKUS9ejBmDFx0EdxyS9jLZeBAqFAh6shEiseqVXD22fDBB3DttfDPf0KZMlFHJaVc0iUvwElRByASl4oVQ5+LrKzQ1G7+/FDIW7du1JGJ7J7586Fr1/BVhbmSRJIueXH3BVHHIBI3s7CfS7NmcP75oaBx+HDtPi6pa9SoUNtlFgpzVd8iSSRpa15EUlL37qELrzscfjgMHRp1RCLxyy/MrVtXhbmSlNI6eTGzfc1slJnNMbPZZnZ17PyeZvaRmc2Lfa0ZdaySRlq3Dr/wDzooJDP33qtCXkkNKsyVFJHWyQuQC1zv7s2A9sAVZtYcuBn42N0zgY9jt0WKT716YYO6c88N00nnnaedqSW55XfMHTAgFOa++27YEkMkCSVdzUs+M7uukNOrgSnuPq0or+HuvwC/xL5fa2ZzgH2AU4CjYw97ARgN3LR7EYtsp0IFePnlUMh7662h6PGdd1TIK8mnYGHuwIHaWFGSXjKPvOQAlxGSjX2A3oSEY6CZ3Rjvi5lZI6A1MBGoE0ts8hOcQvu7m1lvM5tsZpOXLVu2Kz+DlHZmYQn1W2/BzJmhkHfatKijEtlm1KiwmeLSpfDRR0pcJCUkc/KyF3CIu1/v7tcTkplawJHABfG8kJlVITS7u8bd1xT1ee7+tLvnuHtOrVq14nlLkd/r1u33hbxvvx11RCJhlOX448OGipMmwdFHRx2RSJEkc/LSAPitwO0tQEN33wgUeTtfM9uDkLgMdvf8pR9LzKxe7P56wNLiCVlkBwoW8nbrBn//uwp5JRq5uXDNNdC7Nxx3HIwfDwccEHVUIkWWzMnLK8AEM+trZn2Bz4BXzawy8FVRXsDMDHgWmOPuDxa4azjQK/Z9L0Ab00hi1KsXhunPOSfUwaiQVxJt9epQmPvIIyGBGTFChbmScsyT+JOfmWUDHQADxrn75Dif3wEYC8wE8mKnbyHUvQwhjO78CJzh7it29Fo5OTk+eXJcby/yx9zDyMttt0H79urIK4nx7behMHfePHj8cbj00hJ/SzOb4u45Jf5GUqokdfKSTJS8SIl4663QkXevvdSRV0rW6NGh7xCE/+8SVN+i5EVKQjJPG4mkv+7dYexYyMsLhbzvvBN1RJKOnnkmNJ2rXRsmTlRhrqS8pEtezOzQWK2KSOlwyCHwxRfQogWcdhrcd58KeaV4bN0aGs5deikceyxMmAAHHhh1VCK7LemSF0IB7RQze83MLjAzFQJI+svvyHvOOaEvzPnnq5BXds/q1aG+5eGHQ8t/FeZKGkm6DrvufhmAmTUFugDPm1l1YBTwPvCZu2+NMESRklGxIgweHDry3nbbto68depEHZmkmoKFuU89FZZEi6SRZBx5AcDdv3b3h9y9M3AMMA44g7BSSCQ9mYUl1G++CTNmQJs2MH161FFJKvnkk9Axd8kS+PBDJS6SlpI2eSnI3Te6+0h3v1JV61IqqJBX4uUepojyO+ZOnAgdO0YdlUiJSInkRaRUyi/kzcoKHXnvvTckMyLb27gRevUKxbldu6owV9KekheRZJZfyHvuuaEOpnt3WFPk7bmkNFi4EI44Al56Cfr1Cz1cqlaNOiqREpW0yYuZ9S/KOZG0V7Fi+MP08MPw7rthZ+o5c6KOSpLBp59CdjZ8801ocnj77ZCRtL/WRYpNMv9f3qmQc10SHoVIMjALy10//hhWrgwJzNChO3+epCf30N7/2GOhZs2w4WfXrlFHJZIwSZe8mNnlZjYTaGJmMwoc3wMzoo5PJFJHHQVTpoQ6mO7d4W9/C43IpPTYsAEuugiuuAJOOCEkLk2bRh2VSEIlXZ8Xwm7S7wH3ATcXOL92Z5snipQK9evDmDFw1VVw//0hmXn11bA/kqS3b76B00+HWbPgjjugb19NE0mplHTJi7uvBlYD50Qdi0jSKl8+NB9r0yZ8As/OhiFDwnSSpKc33oCLL4Zy5eC998Koi0gplXQpu5mtNbM1hRxrzUzLLEQKuuSS0A8GQj+YBx/Uvkjp5rff4Jpr4Mwzw3Th1KlKXKTUS7rkxd2runu1Qo6q7l4t6vhEkk7btuEP2oknwvXXwymnwArNsKaFhQtDndMjj4SC7TFjYN99o45KJHJJl7wUZGY1zaytmR2Zf0Qdk0hSqlkT3n47LKd+/31o3RrGj486Ktkd77wDrVqF+pYhQ8J/23LlIg5KJDkkbfJiZpcAnwIfAHfFvt4ZZUwiSS1/OfVnn0GZMnDkkfDPf6orb6rZsAEuuwxOOw0aNQoF2WecEXVUIkklaZMX4GqgDbDA3TsCrYFl0YYkkgLatIEvvwzTRzfeCJ07w08/RR2VFEX+ZpxPPQU33BBGzxo3jjoqkaSTzMnLJnffBGBm5d39a6BJxDGJpIYaNcLqlKeeCiMxBx0Er78edVTyR9xhwIBQv7RiRdgN+h//0DSRyB9I5uRlkZnVAN4BPjKzYcDPkUYkkkrMoHdvmDYtfHo/+2zo2RNWrYo6Milo4cIwOnb11XDccWH0pVNhDcZFJF/SJi/ufpq7r3L3O4HbgWeAUyMNSiQVZWbCuHFw553w2mtw8MEwalTUUYk7PPcctGgR/vs89ljYu6pWragjE0l6SZu8mNkZZpa/NeqRwIVAVoQhiaSusmVDN9bPP4cKFeCYY+Dyy7VDdVR+/jnsRXTRRWFF0cyZ8Je/hNEyEdmppE1egNvdfa2ZdSBs0vgC8GTEMYmktvyeMNddB08/Dc2bw4gRUUdVerjDCy+EZnOffBL6t4waBfvvH3VkIiklmZOX/N3mTgSedPdhgKrXRHZX5crwwANhJUvNmmEE4NxzYZkW85WouXPDLtAXXBCSl+nTw/5U2ptIJG7J/K/mJzN7CjgTGGlm5UnueEVSS9u2oYfIXXfBm29Cs2bwzDPqC1PcNm0KU3YHHxxGvZ56Cj79NNQiicguSeZk4ExCY7rO7r4K2BO4IdKIRNJNuXJhd+KpU8MU0qWXQvv2MGlS1JGlhw8+CElLv35hN+ivvw4rwDTaIrJbkvZfkLtvcPeh7j4vdvsXd/8w6rhE0lJWVtg3Z/Dg0NCuXbuwg/HSpVFHlprmzAl7TXXuHOpcPvwwXNs6daKOTCQtJG3yIiIJZhZqX77+OnTmffHFMLXx97/D+vVRR5cafv011LEcdFBY/vzPf4a9idS3RaRYKXkRkd+rWhX69w/Ldzt2hFtvhQMPDLUaW7ZEHV1yWrcO7rsvXKfHHgvTb/PmwV//CuXLRx2dSNpR8iIihWvaNOxsPG4cHHBA2CywRYvQ6G7r1p0+vVTYuBEefDAsdb7lFjj88LCK6IknoHbtqKMTSVtKXkRkxw4/HMaOhWHDYI894JxzwsqkQYPgt9+iji4a69aFHi0HHADXXw8tW4al5yNGhARPREqUkhcR2TkzOPnksO/Om29ClSqhoDczEx59FNaujTrCxFi8OEyjNWgA11wTfv7Ro+Gjj8IqLRFJCCUvIlJ0GRnQvXvoDzNyJOy7byhQrV8/bCw4d27UERY/d/jii5CsNWoUals6dgwjLWPGwFFHRR2hSKmj5EVE4mcGXbqE6aQJE0KX3ieeCHUyxx8Pr78e6kFS2dq18OSTkJ0dGvq99lrojjt3Lrz1lkZaRCKk5EVEdp1Z6Anz8suwcGFoxjZnDpx9NtStC5dcEkYnUqVr76ZN8PbbcNZZoSfL5ZeH4uTHHgubKT75pDrjiiQBc/eoY0gJOTk5Pnny5KjDEEl+W7eGOpCXXgojFOvWhZU3XbuGupnjjoNKlaKOcptffw2dcEeOhOHDw4hLrVqhI26vXmHURbs97zIzm+LuOVHHIekl7ZMXMxsEnAQsdfcWsXN7Aq8DjYAfgDPdfeWOXkfJi8guWL8e3n03rFQaORLWrIEKFcIKpqOOCke7donthbJ6NUycCJ99FpKWSZNCXctee8Epp4RRo44doWzZxMWUxpS8SEkoDcnLkcA64MUCycs/gBXufr+Z3QzUdPebdvQ6Sl5EdtNvv4UamXffDSMzM2aEpKF8+dCRtlWrcBx8cOibUq/e7u0BlJcXtjqYPXvbMXly6HjrHl47JyfU7nTpEr4vU6aYfljJp+RFSkLaJy8AZtYIGFEgeZkLHO3uv5hZPWC0uzfZ0WsoeREpZitWhGRm7FiYNi1sDrlixbb7y5eHhg3DiqY999x2VKwYEo8yZcLXDRvC1NTatbBqVUhY8o+CHYHr1IHWreGww8LRtm3oJiwlSsmLlITSmryscvcaBe5f6e41C3leb6A3QIMGDbIXLFiQmIBFSiP3kHDMnAnff7/t+OknWLkyHCtWQG7u/z63UqWQiFSrBvvsE5Zu77NPSH6yssKx116J/5lEyYuUCE3q7oC7Pw08DWHkJeJwRNKbWUg66tf/48e4h+mgrVvDkZcXamg03SNSqpTW5GWJmdUrMG20NOqARKQIzEKiomRFpFQrrX1ehgO9Yt/3AoZFGIuIiIjEIe2TFzN7FRgPNDGzRWZ2MXA/0MnM5gGdYrdFREQkBZSKgt3iYGZrgTTcuGWX7A0sjzqIJKFrsY2uxTa6Fts0cXct65JiVVprXnbFXFXMB2Y2Wdci0LXYRtdiG12LbcxMPSak2KX9tJGIiIikFyUvIiIiklKUvBTd01EHkER0LbbRtdhG12IbXYttdC2k2KlgV0RERFKKRl5EREQkpSh5ERERkZSi5GU7ZtbZzOaa2Xwzu7mQ+83MBsTun2Fmh0QRZyIU4Vr0iF2DGWb2uZm1jCLORNjZtSjwuDZmttXMTk9kfIlUlGthZkeb2TQzm21mYxIdY6IU4d9IdTN718ymx67FhVHEWdLMbJCZLTWzWX9wf6n5vSkJ4u46YgdQBvgW2B8oB0wHmm/3mD8D7wEGtAcmRh13hNfiMKBm7PsupflaFHjcJ8BI4PSo447w/4sawFdAg9jt2lHHHeG1uAXoH/u+FrACKBd17CVwLY4EDgFm/cH9peL3po7EHRp5+b22wHx3/87dfwNeA07Z7jGnAC96MAGoEdvcMd3s9Fq4++fuvjJ2cwKwg+2AU1pR/r8AuBJ4i/Te6LMo1+JcYKi7/wjg7ul6PYpyLRyoamYGVCEkL7mJDbPkufunhJ/tj5SW35uSIEpefm8fYGGB24ti5+J9TDqI9+e8mPDJKh3t9FqY2T7AacCTCYwrCkX5/6IxUNPMRpvZFDM7P2HRJVZRrsW/gWbAz8BM4Gp3z0tMeEmltPzelATR9gC/Z4Wc234teVEekw6K/HOaWUdC8tKhRCOKTlGuxcPATe6+NXzITltFuRZlgWzgWKAiMN7MJrj7NyUdXIIV5VqcAEwDjgEOAD4ys7HuvqaEY0s2peX3piSIkpffWwTsW+B2fcInpngfkw6K9HOa2cHAM0AXd/81QbElWlGuRQ7wWixx2Rv4s5nluvs7CYkwcYr6b2S5u68H1pvZp0BLIN2Sl6JciwuB+93dgflm9j3QFJiUmBCTRmn5vSkJommj3/sCyDSz/cysHHA2MHy7xwwHzo9Vz7cHVrv7L4kONAF2ei3MrAEwFDgvDT9VF7TTa+Hu+7l7I3dvBLwJ/CUNExco2r+RYcARZlbWzCoB7YA5CY4zEYpyLX4kjEBhZnWAJsB3CY0yOZSW35uSIBp5KcDdc83s/4APCCsJBrn7bDO7LHb/k4SVJH8G5gMbCJ+s0k4Rr8UdwF7A47ERh1xPw510i3gtSoWiXAt3n2Nm7wMzgDzgGXcvdAltKivi/xd3A8+b2UzC1MlN7r48sqBLiJm9ChwN7G1mi4C+wB5Qun5vSuJoewARERFJKZo2EhERkZSi5EVERERSipIXERERSSlKXkRERCSlKHkRERGRlKLkRSTBzGzdTu6vYWZ/KXD7T2b2ZjG8byMzO7fA7RwzG7C7rysikmhaKi1SzGKb8Nkf7WFjZuvcvcoOnt8IGOHuLYo5rqOBv7r7ScX5uiIiiaaRF5FiEBvVmGNmjwNfAvua2Q1m9oWZzTCzuwp5ThUz+9jMvjSzmWaWvyPx/cABZjbNzP4Ze+1ZsedMNLOsAq8x2syyzayymQ2Kvd/UAq9V0P2EzrfTzOxaMzvazEbEXudOM3vBzD40sx/MrJuZ/SMW1/tmtkfscdlmNia24eIH2hlYRKKg5EWk+DQBXnT31rHvM4G2QCsg28yO3O7xm4DT3P0QoCPwQGzU5mbgW3dv5e43bPec14AzAWKJw5/cfQpwK/CJu7eJvdY/zazyds+9GRgbe92HCon/AOBE4BTgZWCUux8EbAROjCUwjwKnu3s2MAi4N47rIyJSLLQ9gEjxWeDuE2LfHx87psZuVyEkM58WeLwBf48lNXnAPkCdnbzHEOAjQvv1M4E3CrzfyWb219jtCkAD4ttT6D133xJrZV8GeD92fibQiJCQtSDsjEzsMdqfRkQSTsmLSPFZX+B7A+5z96d28PgeQC0gO5Y0/EBIOv6Qu/9kZr/GdvM+C+hT4P26u/vcXY4eNsfeI8/Mtvi2grg8wu8KA2a7+6G78R4iIrtN00YiJeMD4CIzqwJgZvuYWe3tHlMdWBpLXDoCDWPn1wJVd/DarwE3AtXdfWaB97syNu2EmbUu5Hk7e92dmQvUMrNDY++xR8H6GxGRRFHyIlIC3P1D4BVgfGwa5k3+N3EYDOSY2WTCKMzXsef+CnxmZrPM7J+FvPybwNmEKaR8dxN28Z0RK+69u5DnzQByzWy6mV27Cz/Tb8DpQH8zmw5MAw6L93VERHaXlkqLiIhIStHIi4iIiKQUJS8iIiKSUpS8iIiISEpR8iIiIiIpRcmLiIiIpBQlLyIiIpJSlLyIiIhISvn/1f0ErPgvVisAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2,1,sharex=True)\n",
    "ax[0].plot(t, volume_transport, color='blue', label='volume transport')\n",
    "ax[0].legend(loc='upper left', bbox_to_anchor=(1.05, 1))\n",
    "ax[0].set_ylabel('volume transport / m$^3$ s$^{-1}$')\n",
    "ax[1].plot(t, salt, color='red', label='salt')\n",
    "ax[1].set_ylabel('salt / g kg$^{-1}$')\n",
    "ax[1].set_xlabel('relative time')\n",
    "ax[1].set_xlim(0,1)\n",
    "ax[1].legend(loc='upper left',bbox_to_anchor=(1.05, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup TEF routine\n",
    "\n",
    "Perform TEF analysis using the `pyTEF` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_description = {\n",
    "    \"lon\" : None,\n",
    "    \"lat\" : None,\n",
    "    \"depth\" : None,\n",
    "    \"time\" : \"time\"\n",
    "}    \n",
    "\n",
    "tef = constructorTEF(inputFile=data,\n",
    "                     data_description=data_description)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sort into salinity classes. From the paper and the data we know what s_min and s_max should be:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyTEF.calc import sort_1dim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tef.tracer = tef.ds.salt\n",
    "tef.transport = tef.ds.volume_transport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constructing var_q and var_Q\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1024/1024 [00:00<00:00, 1599.41it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "That took 0.72 seconds\n"
     ]
    }
   ],
   "source": [
    "s_min  = 10.0\n",
    "s_max  = 31.0\n",
    "N_bins = 1024\n",
    "\n",
    "start = time.time()\n",
    "out_volume = sort_1dim(constructorTEF=tef,\n",
    "                       N = N_bins,\n",
    "                       minmaxrange=(s_min, s_max))\n",
    "end = time.time()\n",
    "print(\"That took {:.2f} seconds\".format(end-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "tef.transport = tef.ds.salt_transport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constructing var_q and var_Q\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1024/1024 [00:00<00:00, 1619.11it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "That took 0.70 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "out_salt = sort_1dim(constructorTEF=tef,\n",
    "                     N = N_bins,\n",
    "                     minmaxrange=(s_min, s_max))\n",
    "end = time.time()\n",
    "print(\"That took {:.2f} seconds\".format(end-start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do temporal mean and compute the numerical bulk values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyTEF.calc import calc_bulk_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_mean = out_volume.q.mean('time')\n",
    "Q_mean = out_volume.Q.mean('time')\n",
    "Qs_mean = out_salt.Q.mean('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2/2 [00:00<00:00, 3220.20it/s]\n",
      "100%|██████████| 2/2 [00:00<00:00, 3921.74it/s]\n"
     ]
    }
   ],
   "source": [
    "bulk_values_volume = calc_bulk_values(out_volume.var_Q, Q_mean)\n",
    "bulk_values_salt = calc_bulk_values(out_salt.var_Q, Qs_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_in_num = bulk_values_volume.Qin[0].values\n",
    "Q_out_num = bulk_values_volume.Qout[0].values\n",
    "\n",
    "Qs_in_num = bulk_values_salt.Qin[0].values\n",
    "Qs_out_num = bulk_values_salt.Qout[0].values\n",
    "\n",
    "s_in_num  = Qs_in_num / Q_in_num\n",
    "s_out_num = Qs_out_num / Q_out_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to analytical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q_in: \n",
      "analytical      813.2395113781663,\n",
      "numerical       813.329298933383,\n",
      "relative error  0.01%\n",
      "Q_out: \n",
      "analytical      -1813.2395113781665,\n",
      "numerical       -1813.3292989333866,\n",
      "relative error  0.005%\n",
      "s_in: \n",
      "analytical      28.423565664110257, \n",
      "numerical       28.422911964641628,    \n",
      "relative error  0.002%\n",
      "s_out: \n",
      "analytical      12.747994132742775, \n",
      "numerical       12.748477110828695, \n",
      "relative error  0.004%\n"
     ]
    }
   ],
   "source": [
    "s_div=np.float64(-s_a*u_r/u_a/np.cos(phi)+s_r) #analytical dividing salinity\n",
    "Q_in_ana  = Q_analytical(s_div)\n",
    "Q_out_ana = Q_analytical(10.) - Q_in_ana\n",
    "\n",
    "Qs_in_ana = Qs_analytical(s_div)\n",
    "Qs_out_ana= Qs_analytical(10.) - Qs_in_ana\n",
    "\n",
    "s_in_ana  = Qs_in_ana / Q_in_ana\n",
    "s_out_ana = Qs_out_ana / Q_out_ana\n",
    "\n",
    "print('Q_in: \\nanalytical      {},\\nnumerical       {},\\nrelative error  {:.2f}%'.format(Q_in_ana,Q_in_num,np.abs((Q_in_num-Q_in_ana)/Q_in_ana)*100))\n",
    "print('Q_out: \\nanalytical      {},\\nnumerical       {},\\nrelative error  {:.3f}%'.format(Q_out_ana,Q_out_num,np.abs((Q_out_num-Q_out_ana)/Q_out_ana)*100))\n",
    "print('s_in: \\nanalytical      {}, \\nnumerical       {},    \\nrelative error  {:.3f}%'.format(s_in_ana,s_in_num,np.abs((s_in_num-s_in_ana)/s_in_ana)*100))\n",
    "print('s_out: \\nanalytical      {}, \\nnumerical       {}, \\nrelative error  {:.3f}%'.format(s_out_ana,s_out_num,np.abs((s_out_num-s_out_ana)/s_out_ana)*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: If you are performing more than one TEF analysis at once, it might be convenient to rename the variable `var_Q ` and `var_q`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_volume = out_volume.rename({\"var_Q\": \"salt_Q\",\n",
    "                                \"var_q\": \"salt_q\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot comparison of analytical and TEF numerical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/boergel/miniconda3/envs/xarray/lib/python3.10/site-packages/xarray/core/computation.py:727: RuntimeWarning: invalid value encountered in arccos\n",
      "  result_data = func(*input_data)\n",
      "/Users/boergel/miniconda3/envs/xarray/lib/python3.10/site-packages/xarray/core/computation.py:727: RuntimeWarning: invalid value encountered in sqrt\n",
      "  result_data = func(*input_data)\n",
      "/Users/boergel/miniconda3/envs/xarray/lib/python3.10/site-packages/xarray/core/computation.py:727: RuntimeWarning: invalid value encountered in arccos\n",
      "  result_data = func(*input_data)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEQCAYAAAB4JulQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPzUlEQVR4nO2dd5hTZdbAf2eGoUkfcahKlSpdQBABUUFdFWUVO6gri2LXz45lLWvf1bWiIvaK2MACKmJXkCJIR5SRDgID0oY53x8nYTJDksnMJLkp7+957nOT97Zzkzc59z3vKaKqOBwOh8MRDTK8FsDhcDgcqYNTKg6Hw+GIGk6pOBwOhyNqOKXicDgcjqjhlIrD4XA4ooZTKg6Hw+GIGjFTKiIyVkTWisjcgLY6IjJZRBb71rVDHDtIRBaKyBIRuT5WMjocDocjusRypDIOGFSs7XrgU1VtCXzqe18EEckEHgOOBdoCZ4hI2xjK6XA4HI4oETOloqrTgI3Fmk8Cnve9fh4YHOTQ7sASVV2mqruA13zHORwOhyPBifecSo6qrgLwrQ8Isk9DYEXA+1xfm8PhcDgSnApeCxAECdIWMpeMiIwARgDUzKjQtVGVmhSQUWTZo4VLfkEGuwsy2ZmfST6Z7PEt/stWrAiVKtlSuTJUqWJLVlYsbtWRTMyYMWO9qtaNx7UC+3Vtyeh6YMWqqGRQIBmoZKDY6wIy2EMG+ZrJ7oIK7C7IZNceW6x/V6DA9+yYmWn9u0qVwr5dubL1dQn2q3OkBdHu1/FWKmtEpL6qrhKR+sDaIPvkAo0D3jcCVoY6oaqOAcYAdKtSRadn7wd//QV/5dk6QnZXrMr2ijXJy6jJ11WP4fcttVi5viabqMVmaqLVa1K3ZS2q9jiEjt0r0aMHtGoFGc5/Lm0Qkd/ida0i/bpWLZ3etauvX/uW7dttvW0b7NoV9lwFGZnsqFyLbVm1WCcH8NuO+izbWI9V1Gc19ViXWZ9qLeuT06MJbXvX5tDuQtu27kEqXYh2v5ZYJpQUkSbAB6ra3vf+fmCDqt7j8+qqo6rXFjumArAIGAD8AfwInKmq80q6Xrdu3XT69OmFDaqwY0fRH+O2bbBlC2zeXLhs2mTrJUvgp5/s9ZYtIa+zjKZspiZ/ZdZgTcMuVGjfiiYda9HqmIOo1O+wUn5KjmRBRGaoard4X3effl2cnTsL+7G/L/tfBy5//glr1sDq1eiqVciff+5zqk3UZCnN+S2jGdvrN6NS1/Y0Oq4DHYa2oWqtirG4PYfHRLtfx2ykIiKvAv2A/UUkF7gVuAd4Q0QuAH4HTvXt2wB4RlWPU9V8EbkE+BjIBMZGolBCCFFov8rOLt2xe/ZAXl5RpTNtGrpyFdkrNlPxt400+vlD+H2a3ckk4N8wre4QMg9py4EdatHgmPZkHntMmUR3OCKmUiU44ABbIkTAHrjWrIFVq2DlSgp+XQ4zl1F/7lIa/z6H2ivfJeuP3fAe7BqZxeKqbchr3onqg3rTdFhfKrQ92NnNHPsQ05FKvCnxiS7aqML27WxbuZmf311G8zuHUylvPdX2bCYj2DTQAw/AlVc6m1mSkrAjlViRn8+2mYtYMn42f34xh4oLZtNs00/UYw0Amyrn8GfnI8k+fzA1hh4L1avHX0ZHuYl2v3ZKJQasWVXA968uo8/ovtT+K+R0EIwdC+edFz/BHOUi7ZRKEDZvUr59YTGrXp9GtRlf0HfnxxzAOnZKJVa2O4Zalw+j9rAT3YRMEuGUShgS6ccHQEEB+uZbLJ+zGV58iaYrpu2zy8bJM6jTtCbUqAF14+JY5CgjTqkUpaAAZvywh5+f/Jqs99+m38bxNCaXdRUbsPzoETR96FL2P7iO12I6SsAplTAk6o8PsF/gQw+x5dMfWP37Lg7+5d19dtl9wy1k3X27B8I5IsEplfDMm7OHWXdPotEHT9Bn20dspRoftbqC6v/6P44+pToVEjGAweGUSjiS5ccHwM8/Q4cO+zTfedochvyjNm0OzzYHA0fC4JRKZKjCorfnsuPG2+m46C1W0Ihbaz3CASMGc/4FwsEHey2hI5Bo92s3Y+wVhxwCM2fCkCEUnHra3uab3+hAm2MaQ9WqvPncVnbu9FBGh6MMiECrIe3puPBN8r/8lhpN6jB20ykccv85dGm1lYEDYdIkG7w7Ug+nVLykUyd46y0y3ngdzjxz383nd6ZxY7j5ZlgbLEzU4UhwKhzek5qLpsPtt3OmvMqvdbuzfdZCjj/egocfeSRsSJgjCXFKJVG47759mhodsJteveDuu+Ggg+CSS2D58viL5nCUi6wsuOUWZPJk6rKeL/b05uM7f6RuXbj8cmjUCK67zkJmHMmPUyqJQsOGZoz2L+ecQ5W1v/HOJ1XJr1iFH7MH8cpTebRoAeeeC/PKFg7qcHjHkUfCN98g1atzzD1H8s1/vueHH+D44y2Eq0kTuPRS+P13rwV1lAenVBKV/Hxbb99Oxs4dtP/jY3696n9cdhmMHw/t28PgwTYt43AkDS1awNdfQ04OHH88h1ZfwKuvwoIFZgF+8klo3hwuuAB+/dVrYR1lwSmVRGX48H2aat53Ew89ZE9yt94K06ZB164wbBjk5sZfRIejTDRoAB9/bGmTjzsONm2iZUt49llYuhRGjoRXXrE5l8suc2axZMMplUTlyCNDbsrOhttusye5a6+F11+Hli1tQj8vL34iOhxlpnlzeOcdWLHCskr4QhsOPBD+9z/L7XreefD447brLbe4Cf1kwSmVRCUzEw49NOwuNWvCPfeY6eCUU+Cuu8y68NRThdYzhyNhOewwuPdeUy5PPVVkU8OG1vTLLzbncscd0KwZPPaY69uJjlMqiYoI/PBD4cT9f/9r7ePHw6efWop+H02awMsv2+6tWpn5oGNHiwVIodhWRypy5ZU2Kr/uOli5b568gw+2kfj06danL7kEunSBqVPjL6ojMpxSSRYOPNDWf/87HHWUTaacf36RXQ49FL74AiZMgN277Qnv6KNh1qz4i+twRISIDUl27YIrrgi5W9euMGWKPVNt2QL9+8PQoc5TLBFxSiVZWLVq37bnntunScS8wubNM9v0rFn2ZPfPf1pZGIcj4WjRAq6/Ht5804bbIRAxM+/8+XD77fDee9C6tYV4OZNY4uCUSrJwwQXB23fvDtqclWWmgiVLzMLwzDPQtq2NYhyOhOOqqyxL9/XXl2izrVLFJu4XLICBA81y1r27c69PFJxSSRb++CN4++LFYQ+rVQsefNAeAA84wJ70Tj45+MDH4fCM6tXNffHzz81XPgIOOgjefhveesumYw491HTS9u0xltURFqdUkoVmzcp1eNeu8OOP5mzz0UeWIPm996Ikm8MRDS68EPbfHx56KOJDRGDIEDOJDRtm/btDB/j22xjK6QhLzJSKiIwVkbUiMjeg7X4RWSAic0RkgojUCnHschH5WURmiUjy5Pz2gkceiXjXrCyLa/npJ8u3dNJJcNFF8NdfMZTP4YiUKlXg4ovh/fdh0aJSHVq7tgVPTpliFuHDDzcTWQjrsCOGxHKkMg4YVKxtMtBeVTsAi4AbwhzfX1U7eVG/ImF55ZV92wYOLPVp2rSB776Da66xtBjdujkPMUeCcPHF9vTz6KNlOnzAAJgzB845x2JbevcutX5ylJOYKRVVnQZsLNb2iar6/TS+AxrF6vopyWmnmd25e3erxwKwbl2ZTlWpEtx/P0yebF5h3bvb3IurceHwlJwcm/R7+WXKWkyoRg0YN86cyZYsgc6dzVHFxWzFBy/nVM4HPgyxTYFPRGSGiIyIo0yJTWamPX59/32h6+U//2m2rIMOMqWzbVupTnnUUfZkd9xxNnIZODBoDJrDET+GDYONG2HixHKd5u9/twKrvXrZdM2wYaX+eTjKgCdKRURuAvKBl0Ps0ltVuwDHAqNE5Igw5xohItNFZPq6Mj61JyV33VX4+o8/LArszTdh9uxSn2r//c3V+KmnLIFshw6WOcPhHWnbr8EiduvVg+efL/epGjY0x5Tbb4eXXrIR+fz5UZDREZK4KxURGQb8DThLNfiAVFVX+tZrgQlA91DnU9UxqtpNVbvVrVs3FiInJldfvW/blVfaY1kZEIERI2wS/6CDzAJx+eVuotMr0rZfA1SoYHnwP/wQNm8u9+kyM23S/pNPzFp86KHBpycd0SGuSkVEBgHXASeqalCfIxHZT0Sq+18DxwBzg+2b1tSqZfbnQP7zn3KftnVrc8e8/HJzLDv66DJP2zgcZeeUU+yJ5sNQFvLSc9RRFiDZuTOcdZaZe/fsidrpHT5i6VL8KvAt0EpEckXkAuBRoDow2ecu/KRv3wYiMsl3aA7wlYjMBn4AJqrqR7GSM6np0SMmp61Y0fJXvvCCTd907Vomq5rDUXZ69rRo3SjbYRs2hM8+g1GjzDHlxBOjMhhyBBBL768zVLW+qmapaiNVfVZVW6hqY5+rcCdVHenbd6WqHud7vUxVO/qWdqp6V/grpTHHHrtvWxQrGp1zDnz1lXnNHH64hQ84HHEhM9P+8SdNKrMXWCj8HstPPGEmscMOMy8xR3RwEfXJys6dFrlYnHr1YPXqqF2ma1cbrbRubcGS993nXDMdceLEE63q3FdfxeT0I0eaUlmzxibwv/giJpdJO5xSSVYywnx1Ua6/2qCBpWM67TRL3nfllS6exREH+va1Ecunn8bsEv37W/qinBw45hhzoHSUD6dUkpWsrNDbXnwx6perUsU8Zq64Ah5+GE4/3SXuc8SYGjVs3jCGSgUsrd7XX1tmiaFDrWSEo+w4pZKKZGbG5LQZGZbr7/77LTPsoEFuktMRYwYMsLKPMS4GVKeO5Q076SS47LKIMvA7QuCUSjLzcojY0euui9klRcwV8+WX4ZtvzEIRZWubw1HIUUeZrTUOEx5VqtjD0siRlu34ggucy3FZcEolmTnppODtDz4Y80ufcQZ88IGVc+nXz9VnccSIHj3Mx/3rr+NyucxMePxxuPVWK6x67rmuqmRpcUolmRk9Onj73XdH3Q0zGAMHWmzaihWmWFzOMEfUqVTJ6mHHsUCKCNx2G/z73zaPePrpsGtX3C6f9Dilksy0ahV62/r1cRHhiCPg449NofTrB7m5cbmsI5047DCbV4nzP/v119sc4vjxlpwyDs9pKYFTKslMw4aht8VxoqN370J//8MPh19/jdulHenAYYfBjh2WTjvOXHklPPaYBf6ecoobsUSCUyrJjL+mSjC++SZ+cmC/+88+gy1bzN/fTd47okbPnrb2qEbwxRdbBu9JkyzPpZtjCY9TKslMdralXC3OoYda8Yg407WrlcBYudLmW2LsBepIFxo3tgjc77/3TIQRIyxf6/jx5hXmgn9D45RKMlOtmhXrUrVl4UJr//FHc83ygMMOg7ffhl9+sSwbfwXNRe1wlJIuXTzPanrFFfCvf1mi1UsucXEsoXBKJdn517/MXUWk6MT9HXd4JtLAgRbU/9VXltrF1WRxlJtOnay61o4dnopx881w7bWWjPLmmz0VJWFxSiXZefvt4O3z5xcqm08+ia9MWLqLJ54wc9h55zlzgaOcdOxokYjz5nkqhgjcc4+Zw+6+2+ZaHEVxSiXZmTWrqOkrGFu2xE2cQP75T6t6/PLLZjpw5gJHmenUydYJUNhHxDzCjj/eJvE9sjQnLBW8FsARJfbbr+j7BPkHv+EG2LDB/P0bN4b/+z+vJXIkJc2a2RzirFleSwJYxePXXrPYrKFDYerU4D4z6YgbqaQKDRvC779brVSwgJG+fS0h34QJnoklYgkohw61lGShrHUOR1gyMqBDh4RRKmA67oMPrEDl3/4Gy5d7LVFi4JRKqpCXBwceaEW4wXIlTZtmwSOnnOKpaBkZlkepRw84+2xzTnM4Sk2HDvDzzwkzCgeriffhhxYUedJJsG2b1xJ5j1MqqUL16lC5cujtHueor1IF3n3XiiGdcIJLQOkoA23aWPDTunVeS1KE1q3NFDZ3rjmlJJDO84SYKRURGSsia0VkbkDbbSLyh4jM8i3HhTh2kIgsFJElInJ9rGRMOcK5W3rsiglmJvjgAxtUnXqqS3nhKCWtW9t6/nxv5QjCwIHmFfbmm5aIMp2J5UhlHDAoSPt/VLWTb5lUfKOIZAKPAccCbYEzRKRtDOVMHapXD73t99/jJ0cY2rUzU9jXX5tHmMMRMX6lsmCBt3KE4JprLI3LzTent0dYzJSKqk4DNpbh0O7AElVdpqq7gNeAEIVDHEX4+OPg7Y0bm90pQTjttMIAsmef9VoaR9LQqBFUrZqwSkUEnnnGfGXOOguWLfNaIm/wYk7lEhGZ4zOP1Q6yvSGwIuB9rq8tKCIyQkSmi8j0dQlma407oSo+rlgBN90UX1lK4O67rajfJZfY3KujKK5fByEjw0YrCapUwOYOx483UYcOTc90+fFWKk8AzYFOwCogWIlCCdIWcupLVceoajdV7Va3bt2oCJm0bN8eelsCBI0FkpkJL70EtWpZrQqP4jMTFtevQ9C6dULOqQTSpImZeKdPj2ll74QlrkpFVdeo6h5VLQCexkxdxckFGge8bwS4moKRcMIJobcNHRo/OSIkJ8e8ZpYssej7dPeacURA69bw228Jn6l08GC47DJ4+GF45x2vpYkvcVUqIlI/4O3JwNwgu/0ItBSRpiJSETgdeC8e8iU9o0fbP3RxVBPO/OWnb1/LifnaazZycTjC4k+aGqyfJxj33WflIM47L2H8ZOJCLF2KXwW+BVqJSK6IXADcJyI/i8gcoD9wpW/fBiIyCUBV84FLgI+B+cAbquptFrlkQQSaN7clkAT33b3+eujTB0aNSt/JTUeENGtm6yToKJUqwRtvWFGvdEqqGkvvrzNUtb6qZqlqI1V9VlXPUdVDVLWDqp6oqqt8+65U1eMCjp2kqgeranNVvStWMqYsRxxR9P3q1d7IESH++ZXMTIu437PHa4kcCUsSKRUwcf/zH0ts8eijXksTH1xEfaqxdu2+7sOHH242pgT+tz7wQMv8+u238MADXkvjSFhq14aaNZNGqYBVijz+eJu0T2DHtajhlEoqsWqVKZR77inavmIF3HorfPqpN3JFyBlnWJqyW27xvGyGI1Hxm3iTSKmIwNNPW4jNueemfo17p1RSiTp1wm//7bf4yFFGRCwgslo1GDkyoQdWDi9p1iyplApA/frw5JOWTPW++7yWJrY4pZJKZGVZSbpg3HqrPSYlOAccAA8+aKWI08UG7SglzZrBr78m3cz3qafCkCFmiU4C57Uy45RKKpGRYfVNgwV8DB4cd3HKyrBhZoO+8Ub773A4itCsmXk0rky+8LVHHjGvsJEjUzcuyymVVGPECLMjFadzZytTlwT4zWAZGeZmnKo/PkcZSTIPsEAaNLApz08/hRdf9Fqa2OCUSqrRpk3w9vbtrQRjktC4sZkJPvzQVYt0FMOvVJYu9VaOMvLPf8Jhh8FVV8H69V5LE32cUkk1rrzSHu0zM4u2//yzuRYnEZdeCh07wuWXw9atXkvjSBgOPNCGswnueBKKjAwYM8bq5t18s9fSRB+nVFKVk5K/WkCFCvD44/DHH/t6STvSmKwsq+O7YkXJ+yYo7dubaffppxMu12u5cUolVTnySK8liAq9elltigcecJP2jgAaN076hFq33mqxnFdckVrzhk6ppCq9ehV9n+D5v8Jxzz1mzbvmGq8lcSQMBx6Y1CMVMIVyxx0wdWpqzRs6pZKKfPcdzC2WALpSJRu9JOHMYKNGlnTy7bfh+++9lsaREDRubEolyR/xL7zQTGHXXAM7dngtTXRwSiXV+Owzcy0JFuj4+efQv3/8ZYoCV1wB2dkWu5Lk/yOOaNC4sdVU+fNPryUpFxUqwH//C8uXmxt9KuCUSqrRvn347S+/HB85okz16nD77aYzU8lU4CgjjX11/JJ8XgVgwAArrX333ZCX57U05ccplVTjgAPsUb744/zpp1tbhw7eyBUFRo60wn+jR7u8YGmPX6kk+byKn7vuMsv0f//rtSTlxymVdGHyZK8lKDeZmRYQOX8+vPKK19I4POXAA22dIkqle3fLpPTAA7Bhg9fSlI8SlYqIvC8i7xVbXhSRy0WkcjyEdJSRwOj6ZO+pPoYMsYwzt92W+inEHWHIybF4lRRRKgB33mnmr3vv9VqS8hHJSGUZsBV42rdsAdYAB/veOxKVo4/2WoKok5FhCmXZMjdaSWsyMqBhw5SYU/HTrh2ceaYF/CbzM2AkSqWzqp6pqu/7lrOB7qo6CugS6iARGSsia0VkbkDb6yIyy7csF5FZIY5d7qtlP0tEppf2phw+evTwWoKYcMIJlr7lrrvc3Epa43crTiFuuAG2bYP//c9rScpOJEqlrogc6H/je13X9zZcRN04YFBgg6oOVdVOqtoJGA+E8+Pp79u3WwQyOgIpKDA3qeK16bt0sWx2O3d6I1eUEIGbboJFi2D8eK+lcXhG48aQm+u1FFGlXTvLsPTII8mb7y4SpXI18JWIfC4iU4EvgWtEZD/g+VAHqeo0YGOwbSIiwGnAq6WW2FEyd99tfopXX120feZMy2T32WfeyBVFTjnFktXef7+LW0lb6te3Etop1gFuuMHCb8aM8VqSslGiUlHVSUBL4Arf0gqYoqrbVPW/ZbxuH2CNqi4OdVngExGZISIhShkaIjJCRKaLyPR169aVUZwUY/jw8Nvr1YuLGLEkM9N+fNOnW22KVMP16who0MDC0Ddt8lqSqNKjh8UoP/BAchoVIvH+GquqO1V1tqrOAjKBSeW87hmEH6X0VtUuwLHAKBE5ItSOqjpGVbupare6deuG2i29eP318Nsvuig+csSYc86xsJxU8O0vjuvXEVC/vq1XrfJWjhhwww12WyX9lBORSMxff4jIEwAiUhuYDLxU1guKSAXgFCDkx6WqK33rtcAEoHtZr5eWHHts+O0PPRQfOWJMpUpw8cUwcaLNrzjSjAYNbJ2CSuWooywi4H//Sz7rXiTmr9HAFhF5EvgEeFBVnyvHNY8CFqhq0Bk2EdlPRKr7XwPHAHOD7esIQdu2hVH1tWsXtn/wgbUVz2CcxIwcCRUr2sSmI83wj1SSsFZ9SYjAJZeYeTfZkqiGVCoicop/AX4AegIzAfW1hUVEXgW+BVqJSK6IXODbdDrFTF8i0kBE/Ca1HMwxYLbvuhNV9aPS3pjDR9u2ha8XLvROjhiRkwNnnAHjxqWcad1REils/gLLCVujRvK5F4cbqZwQsPwNUyhZAe/DoqpnqGp9Vc1S1Uaq+qyvfbiqPlls35Wqepzv9TJV7ehb2qnqXWW7NQcAAwcWvp4zxzs5Ysill5pv/6vOlzC9qF4dqlVLyZEK2K0NHw5vvrlvdEAiE06pTAauUdXzgiznx0tARzlp1qzwdarVLfXRpYslZ36pzDN9jqTF71acoowaBbt3w9ixXksSOeGUyoHAmyLypYjcJiI9fPEljmSic+fC1yk6UhGBs8+Gb76x9C2ONCLFlcrBB8MRR8BzzyXPhH1IpaKq96jqkcBxwGzgfOAnEXlFRM4VkZx4CekoAytXWiKhG28sbCsosH/gnJzkGk9HwJln2jpJy8U4ykqDBilr/vJz3nmwZIk9NCUDkXh/5anqBFX9p6p2Bu7E0rS8EHPpHGXn5pttkuHdd/fdtnZtyj3SN24M/frBiy8mzxOdIwqkaFR9IH//O+y3n41WkoFIgh+7BC5AZeAd4PhYC+coB48/btWsiqdqOfpoOOssq3aVYpxzDixeDD/+6LUkjrhRv76VFd6yxWtJYka1anDaaRYIuW2b19KUTCTBj48D3wFjsFT33wKvAYtE5JgYyuYoDy++CHfcAQ8+WLR98mSzEd1/vzdyxZAhQywg0k3YpxEpHAAZyPDhlmDynXe8lqRkIlEqy7H0991UtSvQGQtGPAq4L4ayOcrDkCFw/PFwyCFF26+8Eq67ruhcS4pQsyaceCK89pp5zDjSgBSPVfFz+OFWPuaNN7yWpGQiUSqtVXWe/42q/oIpmdQyyqcaGzZY/pKffy7a/p//WGm5I0KmU0tqzjkH1q2DTz7xWhJHXPCPVFJ8sj4jw+ZWPv448S19kSiVhSLyhIj09S2PY6avSoB7HkxUatUKv71GjbiIEW8GDoTsbGcCSxvSZKQCNq+ycye8/77XkoQnEqUyHFiCpb2/EisvPBxTKP1jJJejvBQUhN8+bVp85IgzFSvC0KFme070JzpHFKhRw770NWu8liTm9OyZHCawSFyKt6vqg6p6sqoOVtUHVPUvVS1Q1SStTZYG5EQQRrR5c+zl8ICzzrIyG++957Ukjpjjj7tKA6XiN4F99BHk5XktTWgiGak4kpUOHUJvO/tsc35PQbp3N6vIhAleS+KICzk5FnuVBpx0EuzaldjFW51SSWXee89q0l966b7bvv/eapamIBUqWEmZTz+F/HyvpXHEnDQZqQD07m15NCeVt0xiDAmX+v4GEekcarsjCTj3XHjqqeC5sxcvTslU+H4GDjTrnguETAPSSKlUrGjxy5MmJW4SgXAjlV+By0VkpoiME5GhvsqPjmRhwgR47DFzIw7k66/hl1/ssSdFGTDAzO3OtTgN8Ju/SnJOSRGOOw5yc2FugpYuDJdQ8jVf7ZPOwMNAM+BtEZkmIreIiCvxm+h88YXlzr7yyqLtvXtb8a4bbvBGrjiQnQ3dujmlkhbk5MCePbBxo9eSxAV/tfBENYFFNKeiqjNV9d+q2h8r0DUP+EdMJXOUn59+Cr/9u+/iI4dHDBxoU0euImSK4/d0TBMTWIMGVj8oUSfrSz1Rr6pbVHW8qo6IhUCOKHLBBTarV5xRo+Af/0j5UonHHGMPsJ9/7rUkjpiSZkoFoG9fs2InYjqimHl/iUhjEflcROaLyDwRudzXXkdEJovIYt866DyNiAwSkYUiskREro+VnClNTo7lg8/MLNr+2GPwzDMwcqQnYsWLnj1NpzoTWIqThkrliCMsY/HMmV5Lsi+xdCnOB65W1TZAT2CUiLQFrgc+VdWWwKe+90UQkUzgMeBYoC1whu9YR2lYvtxyOuzZE3x7ikcHZmVBjx42teRIYdJUqUBiJsaIpJ7KeBE5XkRKpYBUdZWq/uR7nQfMBxoCJwHP+3Z7Hhgc5PDuwBJVXaaqu7BU+yeV5voOCpPthWLgwPjI4SF9+8L8+SmbPMABULu2BSelkVKpV89KDSfiA1MkiuIJ4ExgsYjcIyKlru4kIk2wlPnfAzmqugpM8QAHBDmkIbAi4H2ury3YuUeIyHQRmb5u3brSipba1KxpzuyqFobr58svre2jj7yTLU509/kozpjhrRylxfXrUpCRAQcckFZKBaBPH5tXSbR4lUhyf01R1bOALlhtlcki8o2InCciWSUdLyLVgPHAFaoaaYo/CSZKCPnG+Gq9dKtbt26Ep08jli+32il3313Y1qePBXH06pXyrlHdutn6hx+8laO0uH5dStIoANLPoYdaUozly72WpCgVItlJRLKBs4FzgJnAy8DhwDCgX5jjsjCF8rKqvu1rXiMi9VV1lYjUB4Il7ckFGge8bwSkdsGEWNG0aeht334LH34IZ5wRP3niTJ060KKFi6xPedJQqXTtauvp08P/zONNJHMqbwNfAlWBE1T1RFV9XVUvBaqFOU6AZ4H5qvpQwKb3MGWEb/1ukMN/BFqKSFMRqQic7jvOUVp69Qq/vaR5lxTg0EOTb6TiKCVpqFQOOcScURLNtBvJSOUZVS0SuykilVR1p6p2C3Ncb2xk87OIzPK13QjcA7whIhcAvwOn+s7ZwHet41Q1X0QuAT4GMoGxgdUnHaXgs89sDmXXLjj55MK5lZdegmrVSlY6KUD37haSs3q1TXA6UhB/qhZVM+2mAZUqmWJJRqVyJ1A8IcC32BxLSFT1K4LPjQAMCLL/SuC4gPeTglzXUVomTrR69cU5+2xbX3JJ8ISTKUTHjrb+6SfLm+RIQXJy7IFp8+aSq56mEF27wltvJZYuDZeluJ6IdAWqiEhnEeniW/phpjBHMtA2RHjP6afDqaemdP4vP5062XrOHE/FcMSSNIxVARup/PlnYpWTCTdSGYiVDW4EBM6J5GFmLEcy0Lp1oc/hTz8Vzu6leIqWQGrXtv+cxYu9lsQRMwKVSqtW3soSR1r7AjwWLIis2Gs8CKlUVPV54HkRGaKq4+MokyOabN9uSbC++qpou3+s3K5d4ubQjiItWzqlktKk6UjFr1Tmz7dA30QgpFIRkbNV9SWgiYhcVXx7MY8uR6Kydeu+CiWQ33+Pnywe0rKleU87UpQ0VSoNG1pV8AULvJakkHAuxf4C5tWA6kEWRzJQUuBcXl7KB0CCKZXVq2FLpOG3juQiO9si69NMqWRkmLVv/nyvJSkknPnrKd/69viJ44gJt9wCd91VNLHk9dfbI07PnmnhLXPwwbZesgS6hPVbdCQlmZn2AJVmSgUsuLek0knxpESXYhGpC1wINAncX1XPj51Yjqhy++22AOy/P2zYYN5ffl/bNKBlS1svXuyUSsqShgGQAI0bW8LxRHErjiSh5LtATWAKMDFgcSQj/n/UpUu9lSPOtGhh6zS77fQijZXKjh32rJgIRBL8WFVVr4u5JI7YsWkTXHuteYJNnmxtgQGRc+eaF1gKU7WqJW1evdprSRwxI039xhv7siSuWGGGCK+JZKTygYi4OORk5vHH4emnLTVLMBLJIBtD0jA7enrhH6kkWi74GBOoVBKBSEYqlwM3ishOYDeWekVVtUZMJXNEj+uug4oVYeVK86v1+x8ec4w9wvfp4618cSJNrSPpQ716NhrPy4Ma6fP31KiRrf/4w1s5/JSoVFTVuQ8nO19/Df/3f/u2BxZvnzAhfvJ4RE4O/PKL11I4YoY/W+iaNWmlVLKzbb1xo7dy+AkX/NhaVReISFBfGX+pYEcS0L07nH++5X/Py4PffrP2SZPsyW7QIG/lixM5OfD5515L4YgZfqWyenWhu18aULGiJRxPeKUCXI25Ej8YZJsCR8ZEIkf0qVwZnn3WXqtaxBTAscd6J5MH5OTYD2/3bqtD4Ugx/FH1aeiNUadOEigVVb3Qt+4fP3EcMWP0aLjzzqJtgU7tu3dDhYgKgSYtfs+YDRtcXZWUJHCkkmYkhVIRkVPCHRhQHtiRDBRXKMXZsCFx0pzGCL+ZPS/PKZWUJDvbIuvTUKnUrp0ESgU4Icw2BZxSSSV+/z3llYrf6rd9u7dyOGJERkbauvhVqWIPS4lAOPPXefEUxOExafB0V7myrZ1SSWHq1UuLvlycSpVg506vpTAiMqKLyPFAO6Cyv01V/1XCMY2BF4B6QAEwRlUfFpH7sVHQLmApcJ6qbgpy/HKsINgeIF9Vu0UiqyME774LJ51UtK1+fcud3akTHH20J2LFkzp1bL1tm7dyOGJITo5TKh5TYkS9iDwJDAUuxQIfTwUOiuDc+cDVqtoG6AmMEpG2wGSgvap2ABYB4erZ9lfVTk6hRIETTzTPL1X46CNr+/e/LW3LU08VPsanMH4/hMBkzY4UI41HKjt2eC2FEUmall6qei7wpy8N/mFA45IOUtVV/lgWVc0D5gMNVfUTVc337fYdVq7YEWvmzzdvL5HCuJThw22GLzMTnnjCU/HiQWamrQsKvJXDEUPq1bOC7Wn2JVesCLt2eS2FEYlS8Vug/xKRBliqlqaluYiINAE6A98X23Q+EKoenwKfiMgMERkR5twjRGS6iExft25dacRKL0r6kXXuHB85PMSvVJJhpOL6dRmpV8/c4//802tJ4kpg+JnXRJpQshZwP/ATsBx4NdILiEg1YDxwhapuCWi/CTORvRzi0N6q2gU4FjOdHRFsJ1Udo6rdVLVb3ZKqHKYz/gmFUDzwQHzk8JBEeZKLBNevy0iaBkDu2ZM4YWaR5P66w/dyvIh8AFRW1c2RnFxEsjCF8nJgXIuIDAP+BgxQDZ5SVFVX+tZrRWQC0B2YFsl1HUGoVQsGDIBPPy3aftFF5g51xx1BD0sl/COVRHmic8SAwADIFC/nEEh+fmH/9ppwwY+HAitUdbXv/bnAEOA3EblNVcOG2oiIAM8C81X1oYD2QcB1QF9V/SvEsfsBGaqa53t9DBDW28xRAlWqwJQphe8PPthMXo8/7p1MccbvHZMGPgnpS5pG1efnJ85IJdwz21OY2y8+09M9mIvwZmBMBOfuDZwDHCkis3zLccCjQHVgsq/tSd81GojIJN+xOcBXIjIb+AGYqKoflf72HCHZf/+0S9nr946pVMlbORwxJDBTcRqRSEolnBiZAaORoVicyXjMDDarpBOr6leYC3JxJgVp85u7jvO9XgakTwH1ePHNN9Cvn/2rbt1qbf78X7VrW0GGKlU8Ey/W/OUbFzulksLUqGFfcJqNVLZuhf3281oKI9xIJVNE/EpnAPBZwLYE0YmOUvHKK+YZ41cogfz5Z8pXzNvicxOp7ioEpS4iaRmrsnmzlctOBMIplVeBL0TkXcyt+EsAEWmBmcAcyUbXruG3L1oUHzk8wq9Uatf2Vg5HjElDpbJpU+IolXC5v+4SkU+B+sAnAV5aGVh0vSPZCKdUBg6EDh3iJ4sHrF9vD7JOqaQ49erBr796LUVc2bzZHDwTgbDOlar6napOUNVtAW2LXNXHJCVcbuyPP4brr4+fLB6wfr2F6ySK66UjRqThSCVZzF+OVKOkjHMDB8ZHDo9Yv76wUJcjhcnJgXXrzCUqDfjrL5smTZQYWadU0omBA82vdtUqeP31wnZ/oskBA7yTLQ6sW5c4PzxHDKlXz/pzmqS3+eMPWzds6K0cfpxSSTcqVbIfnT+pZBqxbp0bqaQFaRYA6ZSKIzGoVs1rCeJObi40cjmxU5/69W29cqW3csQJp1Qc3vLWW+YCFThb7U+JP2SId3LFmD//tMnMxiUWbXAkPf4nB/+/bYrjv80GDbyVw49TKunGxx+H3nZQJLXXkpNly2zdvLm3cjjiQL16ljU0N9drSeLC0qWQnW3JBBIBp1TSjaeftknMwDJxCxda20MPhT4uyXFKJY2oUMEUS5qMVBYuhFatvJaiEKdU0pHdu4sqlQED4KabCssMpyCLF9u6RQtv5XDEiUaN0maksmiRUyoOL/nsM6s9Ghh+m5sLd98Nxx6bsk93P/9s/zNp6J+QnqSJUtmyxSIEDj7Ya0kKcUol3ejYMXi6lpEjbRI/UVxIoszixdCmjddSOOJGmigVf7o+N1JxeEflyjBjxr7tTz4Jf/87fPtt/GWKMTt2wLx50L6915I44kbDhvYYn5fntSQxZeZMWydS2j6nVNINkeCZ52rWhL59oXXruIsUa5YuNcXSubPXkjjiRpq4Fc+YYT/dZs28lqQQp1TSjapVC2unfPJJYfumTTB1akqm8P3Jl/60Sxdv5XDEEb9SSXET2IwZ1q8lWDlEj3BKJZ0JjJb67TfYvt07WWLI9OmmS1NwEOYIRRoold27Yc6cksskxZuYKRURaSwin4vIfBGZJyKX+9pvE5E/itWtD3b8IBFZKCJLRCS1c7LHm2eftUebwEmGJk3sn/eWWzwTK1ZMmwY9e7qU92mF/4EphZXK7Nmwaxd06+a1JEWJ5UglH7haVdsAPYFRItLWt+0/qtrJt+xTs15EMoHHgGOBtsAZAcc6ykvFisHb99sPTjopvrLEmI0b7cfXr5/XkjjiSuXKlpI6hZXKF1/Yuk8fb+UoTsyUiqqu8hfzUtU8YD4Qqb9qd2CJqi5T1V3Aa0Bq/dt5yTnnFKa79/POO1aUIdHG0uXkyy/tNp1SSUMaNkx5pdKyZeLk/PITlzkVEWkCdAa+9zVdIiJzRGSsiASbGW4IrAh4n0vkCslRFmbP9lqCmDB1qj20du/utSSOuHPggTZXmILs2WNm3UR8WIq5UhGRasB44ApV3QI8ATQHOgGrgAeDHRakTYO0ISIjRGS6iExflyZFeaJK06a2njfPWzlixNSp0KuXlZFJJly/jgJNm1qteg3615HUzJljWbf79vVakn2JqVIRkSxMobysqm8DqOoaVd2jqgXA05ipqzi5QGCS8kZA0OIIqjpGVbupare6rqxf6Tn8cFt/8423csSAZJ5Pcf06CjRrBtu2WR3pFMOfbLx/f2/lCEYsvb8EeBaYr6oPBbTXD9jtZGBukMN/BFqKSFMRqQicDrwXK1nTmnbtbJ2CtufJk+0h9cgjvZbE4Qn+Ubg/RXUK8f77Nv2ZaPMpENuRSm/gHODIYu7D94nIzyIyB+gPXAkgIg1EZBKAquYDlwAfYxP8b6hqatpnvGL9enjhBZuc99O4seUA273bO7miyDvvmANQz55eS+LwBL9S+fVXb+WIMuvXWzalv/3Na0mCUyFWJ1bVrwg+N7KPC7Fv/5XAcQHvJ4Xa1xEFevUqzAfvJzcXnnoKLr20cASTpOzcCRMnwtChLj4lbWnSxNYpplQmTbIR+AkneC1JcGKmVBwJzjvvwL//bf++n35qExC//GJpWurV81q6cvPZZ5ZL8OSTvZbE4RnVqsEBB6ScUnnnHahfP3Fz2bk0LelK27bw4ovwxhtw772FbT17wu+/eytbFJgwwf5T3HxKmtO0aUrNqWzaZCOV006zismJSIKK5Ygr++1X+Pq33yxZVhKzbZvpypNOshgVRxrjdytOEd5+24wLZ57ptSShcUolndm6FY4+urCH/uc/1mNPOcVbucrJhAnmw3/hhV5L4vCcpk1t5J2f77UkUeGVV6wk9qGHei1JaJxSSWd++w2mTCl8v3Bh6LxgScRzz1mIQqLlRHJ4QMuWplBSILL+jz9srvDMMxMr1X1xnFJJZ9q1M/fhP/6wCYhkCzsPwm+/2Q9v2LDEtTk74oi/3sH8+d7KEQWeeca8vs4912tJwuN+dulOhQoWQbV1Kzz8MNx1l9cSlYsXXrB1ov/wHHHCr1QWLPBWjnKSnw9PPw0DB0Lz5l5LEx6nVNKdTZvgjjsK3998s42tk7CuSkEBjBtnHl/+EAVHmlO7NuTkJL1Sef99MyhcdJHXkpSMUyrpzg03BFcgd9yRdIn4PvzQvEfPP99rSRwJRevWSW/+euIJK2Z5/PFeS1IyTqmkO+Ee6deujZsY0eDBB+2Hd9ppXkviSCj8SiXJHpL8zJ5teewuusis1YmOUyrpjj8/UjB27oyfHOXkp5/g88/h8sshK8traRwJRZs28OefSZut+J57oHp1uPhiryWJjCTQe+Vj9+7d5ObmsmPHDq9FSUwOOQR++MHStPif5GrWtGXbtnKZDSpXrkyjRo3IisO//L332g8vXWJTXL8uBX36mG00NzcqiiWe/XrpUgvkveYaqFUr5peLCimvVHJzc6levTpNmjRBEtm522sCFch++9nTXTlQVTZs2EBubi5Nw42GosDs2fbDu/FG04XpgOvXpWDXLqtq1aCBTdqXg3j2a4D777eR9xVXxPxSUSPlzV87duwgOzvb/fDCkZdXNGPxtm2WqmXLljKfUkTIzs6Oy5P0zTfbU9z//V/ML5UwuH5dCrKybDLir7/Kfap49uslS+DZZ+GCCyyBZLKQ8koFcD+8klizJngaiz17ynXaeHzuM2fCBx/A1Vcnj3kgWrh+HSEiULVqVJSKnS4+n/vo0ZbgYvTouFwuaqSFUnGUQLNmwR+FksDV5PrrLRThkku8lsSR0FStCjt2WDBTEvDTT/Daa3DllclXicIplTTgvffe45577gm9w8aNsGrVvu3LltGkSRPWB5nc3Lx5M+eeey7NmzenefPmnHXWWfz5559RlLpkPvwQPvmk0PzlSE7GjRvHJSU8FYwbN46VK1fuff+Pf/yDX375JfKLVK0Kqkz9+GP+FqJk4ldffUX37t1p3bo1rVq14rHHHov8/FFEFa67DrKzk9Ok65RKipOfn8+JJ57I9ddfH3qnUPnhw5QVvuCCC2jWrBlLly5l6dKltGjRguHDh5dP2FKwc6c9xbVs6UYp6UBxpfLMM8/Qtm3byE9QtaqtQ7jJr169mjPPPJMnn3ySBQsW8PXXXzN27FgmTJhQHrHLxFtvWZ7XW29NTseTmCkVEWksIp+LyHwRmScil/vaXw+oWb9cRGaFOH65r5b9LBFJ2gIfy5cvp02bNlx44YW0a9eOY445hu3btwPQr18/pvtql6xfv54mvkDEcePGMXjwYE444QSaNm3Ko48+ykMPPUTnzp3p2bMnGzduBGDp0qUMGjSIrl270qdPHxb4UlEMHz6cq666iv79+3PdddcVeRJcs2YNJ598Mh07dqRjx4588803kJXF4Guuoes559DutNMY8/bbYe9pyZIlzJgxg9EBxt5bbrmF2bNns3Dhwqh+fqF44AFLqvzwwymRWDkpGTx4MF27dqVdu3aMGTNmb3u1atW46aab6NixIz179mTNmjUAvP/++/To0YPOnTtz1FFH7W33k5eXR9OmTdnte5jZsmULTZo04c0332T69OmcddZZdOrUie3btxf57Xz00Ud06dKFjh07MmDAAAB++OEHevXqRefOnenVqxcLly83c26ICfbHHnuM4cOH06VLFwD2339/7rvvPu6///6ofmYlkZdnD0udOiVHSpZgxNJong9crao/iUh1YIaITFbVof4dRORBYHOYc/RX1ahFLF1xBcyaFa2zGZ06wX//G36fxYsX8+qrr/L0009z2mmnMX78eM4+++ywx8ydO5eZM2eyY8cOWrRowb333svMmTO58soreeGFF7jiiisYMWIETz75JC1btuT777/n4osv5rPPPgNg0aJFTJkyhczMTMaNG7f3vJdddhl9+/ZlwoQJ7Nmzh61bt0LFiowdPZo6NWuyfccODh02jCFHHkl2CJvSL7/8QqdOncgMKP6emZlJ586dmT9/Pq1atYrgkys7S5ZY3suTT4Zjj43ppZICr/r12LFjqVOnDtu3b+fQQw9lyJAhZGdns23bNnr27Mldd93Ftddey9NPP83NN9/M4YcfznfffYeI8Mwzz3Dffffx4IMP7j1f9erV6devHxMnTmTw4MG89tprDBkyhFNPPZXHHnuMBx54gG7duhWRYd26dVx44YVMmzaNpk2b7n3gat26NdOmTaNChQpMmTKFG2+6ifH33AO+B7rizJs3j2HDhhVp69atW+lMbFHgX/+yHF9vvpkUU5pBiZnYqroKWOV7nSci84GGwC8AYi4UpwEpX/C1adOmdOrUCYCuXbuyfPnyEo/p378/1atXp3r16tSsWZMTTjgBgEMOOYQ5c+awdetWvvnmG0499dS9x+wMGNqfeuqpRf70/Xz22We84Evlm5mZSU3f+PqRiROZ8MEHAKxYs4bFK1aQ7ZO5OKoa1ANG45AGo6DAcntVrAj/+1/ML+cIwyOPPLLXPLRixQoWL15MdnY2FStW3Dtv0bVrVyZPngxYbM3QoUNZtWoVu3btChrn8Y9//IP77ruPwYMH89xzz/H000+HleG7777jiCOO2HuuOnXqADbnN2zYMBYvXoyI2OinenWLWQnST0P16XgyfbrVybvgAjjsME9FKRdx0YUi0gToDHwf0NwHWKOqi4MeBAp8IiIKPKWqY0LsFzElPXnFikoBdUoyMzP3mr8qVKhAgc8bpbjfe+AxGRkZe99nZGSQn59PQUEBtWrVYlaIR9T9AksEl8DUqVOZ8sMPfDt9OlWXLKHfhReyY9cu2H//oPu3a9eOmTNnUlBQQIavaElBQQFz5szZaz6IFf/7H3z5pRXiatgwppdKGrzo11OnTmXKlCl8++23VK1alX79+u3tw1lZWXv/oDMzM8n3uatfeumlXHXVVZx44olMnTqV2267bZ/z9u7dm+XLl/PFF1+wZ88e2rdvH1aOUMpg9OjR9O/fnwkTJrB8+XL69etnNYMgqPt8u3btmD59OieeeOLethkzZuwzMooVO3ZYuYZ69SzgMZmJ+US9iFQDxgNXqGpgNN0ZwKthDu2tql2AY4FRInJEiPOPEJHpIjJ93bp1UZM7HjRp0oQZM2YA8NZbb5Xq2Bo1atC0aVPefPNNwH5cs2fPLvG4AQMG8MQTTwCwZ88etmzZwubNm6lduzZVa9RgQdWqfDd3LgQZ5fhp0aIFnTt35s4779zbdueddzJgwAAOPPDAUt1HaVi82JIqH3+8FeFKZRK9X+/tM1WrsmDBAr777ruIjmnoexJ4/vnnQ+537rnncsYZZ3DeeeftbatevTp5eXn77HvYYYfxxRdf8KuvDr3f/BV4rb3m36pVLWYliFIZNWoU48aN2/uQtmHDBm666aYi84axZPRoS2jx7LPmIp/MxFSpiEgWplBeVtW3A9orAKcAr4c6VlVX+tZrgQlA9xD7jVHVbqrarW7dutEUP+Zcc801PPHEE/Tq1Suo225JvPzyyzz77LN07NiRdu3a8e6775Z4zMMPP8znn3/OIYccQteuXZk3bx6DBg0iPz+fDh06MPrmm+nZoYMFPs6caeaCNWv2MRmMHTuWxYsX06JFC+rWrct3333Hk08+Wep7iJQ9ewrNXk89ldjlVKNBovfrIn1m9Gh69uxZ4jG33XYbp556Kn369GH/EKNgYK97+hlnnLG3bfjw4YwcOXLvRL2funXrMmbMGE455RQ6duzI0KE2ZXvttddyww030Lt3b/b4g3gzMqBKlaBKpX79+rz00kuMGDGCVq1a0aBBg73zj7Hm888tw/Y//2lFuJIeVY3JAgjwAvDfINsGAV+EOXY/oHrA62+AQSVds2vXrlqcX375ZZ82Rxh++031xx/3XQoKQh6yYMECbdasmU6cOHGfbdH6/EePVgXVF16IyunKBDBdY/R7CbekW79+88039eyzz47NyVevtv68Y0fY3R599FFt3769bty4Mej2aH3+K1eq5uSotmqlmpcXlVOWmmj361jOqfQGzgF+DnAbvlFVJwGnU8z0JSINgGdU9TggB5jgs5VWAF5R1Y9iKKvDT5Uqwdt37w7pu9uqVSuWLl0aM5EmTrSaYeedB+ecE7PLOBKASy+9lA8//JBJkybF5gI1a8KKFbB5MxxwQMjdRo0axahRo2Ijg4/8fDjjDEuxN2VK4ZRPshNL76+vsNFKsG3Dg7StBI7zvV4GdIyVbI4w1KljtScCk0m2aeNZMMiyZaZIOnUCjwKcHXHkf7F26atcGSpVsjLaYZRKPLjxRvjiC3j+eSjBHyGpcBH1jqJkZsLBBxdtC5ZsMg6sXQuDBtl0zvjxoQdRDkepqF3bogzDZIyINc8+a15eF11kXl+phFMqjn0pnvJ+8WJ7sosjmzfbpGVurmUhbtYsrpd3pDJ16tiTSpxz1fn57DMYORKOOQYeecQTEWKKUyqOfQk2PxIs4WSM2L4dTjgB5s6Ft9+G3r3jdmlHOlC1qg17N2yI+6XnzIEhQ6BVKyssl6xR8+FwSsWxL74cZEUoRTBledi9G047Db76Cl580cxfDkfUyc62YnRRqrESCQsWwFFH2YT8xInJmSwyEpxSSQNKTH1fnIBJ+SYnnsj6TZtsgmP69KjU+A6FPwXLBx/YpPzpp8fsUo4EIi6p731MnTrVUsjsv7/FrRRLahkrli2DAQPskp9+CgcdFJfLekIKDr4cgfhT3wemnyiR/faDrl2D25xjVEZV1bKzvvSSuQ8na4ZWR2wYN24c7du3p0GDBoClvi8XFSqYYlm3zvL9xNC7cfFiOPpoM+t+8cW+fjCphhupxJikSH1PkDTmIsHH52vXAvDcc89x8MEH07dvXy688MISnzTDUVAAV11lk5ZXXAE33VTmUzniSFKlvg9WkiEnx9arVrF9+3ZOP/10OnTowNChQ+nRo8fe85eHWbPg8MPN0jZlChxySLlPmfCk10jFoxzhCZ/6ntBpzPehoIBVq1Zx6623MmPGDGrWrEn//v3p3LlzaT61vezcCSNGwAsvwOWXW7qKVE/BEnU86tdJlfr+xhsZP3580RuoVAnq1oW1a3nitdeoWrUqc+bMiVpi1C+/hL/9zZ7NJk+2yfl0IL2UikckRer7EGnMycoqeoLatfn+++/p168f/pxUQ4cOZdGiRZF9GAGsXm2eMN98Y3Ukbr7ZKZRkIulS3wejfn3YsIFpkydz2Q03ANChQwc6dOhQ+g8kgLFjzW24eXMred24cblOl1Skl1LxKPd9UqS+D5HGHBELEtm2zYYV27bB+vVIOWunTJ8OgwfbtM0bb0CAbnSUFg/6dVKmvg9GVpb94+fnI1GIxcrPh6uvNlPu0UfD668nf9bh0uLmVDwk4VLfh0pjXqNGoTF41y56NGrE1ClT2LBhA7t3794rQ6S89BL06WNzpV9/7RRKMpKUqe9DkZ3NEb168fKLL0JeHnPnzmXOnDkl3k9x/vjDFMkjj5jTyaRJ6adQwCkVT0nI1Pfh0phXrgxA/f3357YLL+Swww7jqKOOitj+vHEjnH225fLq0QN+/NFM947kIylT34dChItuuIGtO3fSoWtX7vv3v+nePWiljZC8+y506GB9+vnn4aGHUjOwMSKimfLY6yXdUoTHnb/+CpoW/7knn9RRo0YFPcT/+b//vmq9eqoVKqjedpvqrl3xFDw64FLfx4WYpr4Px/btqrNmqc6cqX379NEff/wx5K7+z3/9etXzzlMF1S5dVBcujJew0SPa/TpddamjLFSpAt26wcqVtvgJY4suKLCU9ePGmQVt4kSIccVhRxIT89T34ahc2Vy0Fi2ySPsS5lheecUc7zZutIqkt93mWTLvhMIpFUfpqV+/iFIZPmAAw/fbz6K7fKmECwosruyPPyzdyk03WcnUAP8Dh2MfYp76viQqV4Y2bZj60kuwdavlwWvcuIi22LbNAvHPOgu6d7f4k3I6i6UUaaFUNISHiKOMiNiIBcyNC+yXNm8e2q4dm/+qyIqVmezcqWRlWVXidAj6ijeuX8eIrCwLe1+zxh6efAW9dtbO4Y81WWzcqOzcCY8+am7DQTz305qUVyqVK1dmw4YNZGdnux9gLOjYEXJz0Q0bEEDmzaMmsDGrJTXq7iI7uzJBwhEc5cT16xiTkWEj8tq1yf99JZmrV5O1eg01qQVVsqjRpGJq1JOPASmvVBo1akRubi7r1q3zWpSUZM8esxJszatK5p6dVMrMp9aejUjBWipPW0qjiRNh2DCrHulcvaKG69exRdWsuXl5lu6uolSkdlYelXctpvKiRTR94AHo2dOK/vTqZXMxGc6ZFtJAqWRlZQWN3HWUnW3bzAf/1Vfh/fct4OvII2HUKAtozPjx+0Lfyhkz4OOP7cAaNaBRIyuQ0rEjHHqozdqnre9l2XH9Ovrk51vJhbfftr69fr1Np1x0kS21amGaJj/fajK8/z74Y7Rq1bK+3LatPUC1bWsmtHr10k7ZiJYzMjrkiUUqA9OASpjyektVbxWROsDrQBNgOXCaqu6TDldEBgEPA5nAM6paYu72bt26aTSSwDmKogpLllj+osmTLe3EX39Zie9zzrHcXSEzr27aBAsX2oGrVtkvdt06G+KAzc/UrWumhnBLTo45AXho6hGRGarareQ9o4vr17Fj5UrLHDxlCrz3nimSypXh+OPhH/+wYMaQcyaq5in27beWa2jOHPjlFxve+KlQwR6kDjzQ+nHdukWX/fe3h60aNaB6dVvi3M+j3a9jqVQE2E9Vt4pIFvAVcDlwCrBRVe8RkeuB2qp6XbFjM4FFwNFALvAjcIaqhi2g4H585eevv+D33y1d98yZ8NNPNtjIzbXtBx0Exx5rhbSOOKKMk5R79lgCxOnTTdGsWmW/bv/rNWsKlU5xKle2H10k65L2ycrad6lYMXh7hQpI8+ZOqSQp27bB8uWmA+bMgdmzbVm2zLbXqGGK5JRTbBBSrVoZL6RqfXnePDv5778XLqtX2wNVSelgMjNNufgVjb9SZWDf9b+uWLFw8fffYK8zM03BZWYWXSpUQAYOjGq/jpndwRdUs9X3Nsu3KHAS0M/X/jwwFbiu2OHdgSWqugxARF7zHVf6qjyOfVCF9059kV3LctmxA3bshJ3blc2bYdtfINiDhgB994cz6ysHdbbkeNm11R6ivvItgQ8loV6XtO2AA2zp1Mm27dljj4yBCmftWvNT3rEjZjVdHKnDvHlw7rnWjdavL1rgUQRatjRr1ciR0L+/db2oWGFFrD6LL0VMUHbvNqHWrbN1Xp4tW7bs+3rLFhN++3ZTRtu3Fy47dsCuXXa+XbtCP4jFmZgas30jjhlAC+AxVf1eRHJUdRWAqq4SkQOCHNoQWBHwPhfoEeIaI4ARvrc7RWRu1G4geuwPxK5kYtkpWa71vuXneIizl0T9vOKWvNz163IRVi6/1WrRIihlyr3ykqifV1T7dUyViqruATqJSC1ggoiETzlaSDCDYlA7naqOAcYAiMh0L8wTJeHkKh2JLFe8ruX6ddlxcpWOaPfruLglqOomzMw1CFgjIvUBfOu1QQ7JBQIrEDQCVgbZz+FwOBwJRMyUiojU9Y1QEJEqwFHAAuA9YJhvt2FAsNS6PwItRaSpiFQETvcd53A4HI4EJpbmr/rA8755lQzgDVX9QES+Bd4QkQuA34FTAUSkAeY6fJyq5ovIJcDHmEvxWFWdF8E1x5S8iyc4uUqHkysxrlsSTq7SkRZyxcyl2OFwOBzpR3qFejocDocjpjil4nA4HI6okTRKRUROFZF5IlIgIt2KbbtBRJaIyEIRGRjQ3lVEfvZte8QX5Y+IVBKR133t34tIkyjJeJuI/CEis3zLcWWVMZaIyCCfHEt8WQ3iiogs993zLL87o4jUEZHJIrLYt64dsH/Qzy4KcowVkbWBMSBlkaM832Ey9GvfuV3fjuz6rm9Hs4xkLBegDRakMxXoFtDeFpiN5RhrCiwFMn3bfgAOw+JePgSO9bVfDDzpe3068HqUZLwNuCZIe6lljOHnmOm7fjOgok+utnH+LpcD+xdruw+43vf6euDekj67KMhxBNAFmFseOcrzHSZDv3Z92/Xt0nyHSTNSUdX5qrowyKaTgNdUdaeq/gosAbqLxcDUUNVv1T6dF4DBAcc873v9FjAgxk9RZZExVuxNgaOquwB/ChyvCfxOnqfod7XPZxeNC6rqNGBjeeQo73eY5P26rHLGCte3fXjZt5NGqYQhWEqXhr4lN0h7kWNUNR/YDGRHSZ5LRGSOb/jpH16WRcZYEUqWeKLAJyIyQywdCUCR9D2AP31PvOUtrRyx+g4TrV+D69uRkPZ9O6EKWYjIFKBekE03qWqwIEkIndIlXKqXiNPA7HOxMDICTwB3+M51B/AgcH4ZZYwVXlyzOL1VdaVY3rfJIrIgzL6JIC+U4ztMhn4Nrm9HibTq28FIKKWiqkeV4bBQKV1yfa+LtwcekysiFcAq4EZTRhF5GvigHDLGCs9T4KjqSt96rYhMwIb8a0SkvlqS0cD0PfGWt7RylPgdJkO/Lo2crm+HJt36djBSwfz1HnC6z/OlKdAS+ME3vMsTkZ4+u/K5FKaECUwV83fgM5/NsFz4vig/JwN+z4uyyBgrPE2BIyL7iUh1/2vgGOxzCpW+J+hnF0MRSyVHDL/DhOnX4Pp2JLi+7SMangbxWLCOnAvsBNYAHwdsuwnzWFhIgHcC0A37UpcCj1KYQaAy8CY2IfUD0CxKMr6IJYmf4/ui6pdVxhh/lsdhRdCWYiaYeH6PzTBPk9nAPP/1Mdv/p8Bi37pOSZ9dFGR5FVgF7Pb1rQvKIkd5vsNk6Neub7u+XZrv0KVpcTgcDkfUSAXzl8PhcDgSBKdUHA6HwxE1nFJxOBwOR9RwSsXhcDgcUcMpFYfD4XBEDadUHA6HwxE1nFIpAyLSRkSeFJG3ROSicp6rmYg8KyJvRUu+aBLNe/WdL6HvN5BkkjVauL5drvMl9P0GElNZ4xkcFIfgo0ZYxOdiYBkWrFMpzP5PYbl6ynq9DODZKMn+VphtU4Em5Tx/wtxrBPdbBfiCCNKA+++LgBTfERxTEZgGVCivrPFaXN8uuQ8kwr1GcL8p37dTZqTiSyPwNvCOqrbEUg1UwWoIhKIH8F0Zr3ci8BUWmRrJ/oeIyAfFlgNKPjJqxO1efceU537PB95W1T0R7NsDWB2pXABqadE/BYZGQdaY4/p2ibi+7SMh+na0tZRXCzAAmFasrQbwJ1AtyP5tgDeKtTUBFgDPYKkJXgaOAr7GnhC7BznPxCBt+wETsXQNc4Gh5XliIMzTHJZeYSEwBUvNEKyQ0j73Wpb7jda9RnC/3/jvFxjtk3Fy8fvz35fvPub62poBM4FDwx0PdAQmlVdW17dd33Z9u9g54/XDiPUCXAb8J0j7TKBTkPargPODdMR84BBsSDwDGIulgD4Je1IE6Ac8gg1PRwU59xDg6YD3NcPInQ08ieXWuSHEPkF/eEBXLB9TVexPZkmIH94+9xrp/UbzXiO5X2z4vtr3uhswC3sqr+77M7im+H35f3hYBcW933e447EqgevKI6vr265vu76975JQqe/LiRA813+oyncDgfOCtP+qqj8DiMg84FNVVRH5GfuCUdWp2I8hFD8DD4jIvcAHqvplqB1VdQMwMsy5wtEHmKCqf/nkDZWRNdS9Qgn3G817hYjud39gk+/14cC7qrrdJ9/7xfb131dFoC425zBEVeeVdLyq7hGRXSJSXVXzyihrvHB92/XtpOnbKTOngmUF7RbYICI1gBxsCB3YXhWopb7aB8XYGfC6IOB9ARHWn1HVRRQ+af1bRG6J5LgyEjYjaAn3CuW83xjc63Ys2y6E/tMMdl+bsep1vQN3K+FalYAdZZQznri+HQTXt8PiWd9OJaXyKVBVRM4FEJFMrDrdo35tHkB/4PNYCSIiDYC/VPUl4AGgS4wuNQ04WUSq+Oo4nBBkn6S6V1X9E8gUkcrYBOoJIlJZRKoBxwfsWvy+dmH1s88VkTN9bSGPF5FszESwuzzyxgnXt13fHkyS9O2UMX/5hrUnA4+JyGhsyPi6qt4VZPdjgVj6kh8C3C8iBVg9g3L7wAdDVX8Skdcx2+pvQLDheTLe6yfA4ao6xWf2mI3d33TsqQ2C3JeqbhORv2FlXLep6rthju8PTIqCrDHH9W3Xt5Oqb0dzgiaRFqAX9mF3DbLtJyDLaxlLcS9TicCXH7iNYpOZyXavPpk7Ay/6XlfzratiP5wupbmvMMe/DbTy+l7L+Pm4vp2E9+qTOeX7dsqMVIqjqt8AB4XYFqshe8KRjPeqqjNF5HOfmWeMiLTFbNHPq+pPvn0iva99jhcrNfuOqi4s4diExPVtIxnvNR36tqv8mASIyHCso2zyWBSHI6q4vp16OKXicDgcjqiRSt5fDofD4fAYp1QcDofDETWcUnE4HA5H1HBKxeFwOBxRwykVh8PhcEQNp1QcDofDETWcUnE4HA5H1HBKxeFwOBxR4/8BA9SvkHqjWk4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=plt.subplots(1,2,sharex=True,sharey=True)\n",
    "ax[0].plot(Q_mean,out_volume.salt_Q,color='blue',label='numerical Q')\n",
    "ax[0].plot(q_mean,out_volume.salt_q,color='red',label='numerical q')\n",
    "ax[0].legend(loc='lower left')\n",
    "ax[0].set_ylabel('Salinity / g/kg')\n",
    "ax[0].set_xlabel('Q / m$^3$ s$^{-1}$ | q / m$^3$ s$^{-1}$ (g/kg)$^{-1}$')\n",
    "ax[1].plot(Q_analytical(out_volume.salt_Q),out_volume.salt_Q,color='blue',label='analytical Q')\n",
    "ax[1].plot(q_analytical(out_volume.salt_q),out_volume.salt_q,color='red',label='analytical q')\n",
    "ax[1].legend(loc='lower left')\n",
    "ax[1].set_xlim(-1000,1000)\n",
    "ax[1].set_ylim(30,10)\n",
    "ax[1].set_xlabel('Q / m$^3$ s$^{-1}$ | q / m$^3$ s$^{-1}$ (g/kg)$^{-1}$')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "5ddaae1b1aacab30e6f9d299f885d3d1a71d6c802be31530f91bd978b8a5f689"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}