{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a6fc91cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Stommel/ Marotzke model\n",
    "## Marotzke, J. (1990). Instabilities and multiple equilibria of the\n",
    "## thermohaline circulation. PhD thesis, Berlin Instit Meereskunde, Kiel.\n",
    "## Marotzke, J. (2000). Abrupt climate change and thermohaline circulation:\n",
    "## Mechanisms and predictability. Proceedings of the National Academy\n",
    "## of Sciences, 97:1347–1350.\n",
    "\n",
    "## Eli Tziperman, 201910\n",
    "## Edited by Kirstin Koepnick, 202603\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "88cdc18a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# font required by PUP:\n",
    "plt.rcParams['font.family'] = 'Myriad Pro'\n",
    "\n",
    "########################################################################\n",
    "def q(DeltaT,DeltaS):\n",
    "    ########################################################################\n",
    "    # MOC transport in m^3/sec as function of temperature and salinity\n",
    "    # difference between the boxes\n",
    "    flow=k*(-alpha*DeltaT+beta*DeltaS);\n",
    "    return flow\n",
    "\n",
    "########################################################################\n",
    "def steady_states(Fs,X):\n",
    "    ########################################################################\n",
    "    # 3 steady solutions for Y:\n",
    "    y1=X/2-0.5*np.sqrt(X**2+4*beta*Fs/k)\n",
    "    if X**2-4*beta*Fs/k>0:\n",
    "        y2=X/2+0.5*np.sqrt(X**2-4*beta*Fs/k)\n",
    "        y3=X/2-0.5*np.sqrt(X**2-4*beta*Fs/k)\n",
    "    else:\n",
    "        y2=np.nan\n",
    "        y3=np.nan\n",
    "    Y=np.array([y1,y2,y3])\n",
    "    return Y\n",
    "\n",
    "########################################################################\n",
    "def rhs_S(time,DeltaS,time_max,DeltaT,hysteresis):\n",
    "    ########################################################################\n",
    "    # right hand side of the equation for d/dt(DeltaS):\n",
    "    # Input: q in m^3/sec; FW is total salt flux into box\n",
    "    rhs=-(2*np.abs(q(DeltaT,DeltaS))*DeltaS+2*Fs_func(time,time_max,hysteresis));\n",
    "    return rhs/V\n",
    "\n",
    "########################################################################\n",
    "def Fs_func(time,time_max,hysteresis):\n",
    "    ########################################################################\n",
    "    # total surface salt flux into northern box\n",
    "\n",
    "    ########################################################################\n",
    "    # Specify maximum and minimum of freshwater forcing range during the\n",
    "    # run in m/year:\n",
    "    FW_min=-0.1;\n",
    "    FW_max=5;\n",
    "    ########################################################################\n",
    "    if hysteresis:\n",
    "        if time < time_max/2:\n",
    "            flux=FW_min+(FW_max-FW_min) \\\n",
    "                *time/(time_max/2);\n",
    "        else:\n",
    "            flux=FW_max-(FW_max-FW_min) \\\n",
    "                *(time-(time_max/2))/(time_max/2);\n",
    "    else:\n",
    "        flux=2;\n",
    "\n",
    "    # convert to total salt flux:\n",
    "    return flux*area*S0/year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7e359091",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finding steady states, including unstable ones, as function of F_s...\n",
      "calculating dDeltaS/dt as function of DeltaS for stability analysis...\n",
      "doing a hysteresis run...\n"
     ]
    }
   ],
   "source": [
    "########################################################################\n",
    "## Main program:\n",
    "########################################################################\n",
    "\n",
    "# set parameters:\n",
    "meter=1;\n",
    "year=365*24*3600;\n",
    "S0=35.0;\n",
    "DeltaT=-10;\n",
    "time_max=10000*year;\n",
    "alpha=0.2; # (kg/m^3) K^-1\n",
    "beta=0.8; # (kg/m^3) ppt^-1\n",
    "k=10e9 # (m^3/s)/(kg/m^3) #  3e-8; # sec^-1\n",
    "# area of each of the ocean boxes, rediculous value, needed to get\n",
    "# all other results at the right order of magnitude:\n",
    "area=(50000.0e3)**2\n",
    "depth=4000\n",
    "V=area*depth # volume of each of the ocean boxes\n",
    "Sv=1.e9 # m^3/sec\n",
    "\n",
    "print(\"finding steady states, including unstable ones, as function of F_s...\")\n",
    "# -----------------------------------------------------------------------------\n",
    "Fs_range=np.arange(0,5,0.01)*area*S0/year\n",
    "DeltaS_steady=np.zeros((3,len(Fs_range)))\n",
    "q_steady=np.zeros((3,len(Fs_range)))\n",
    "\n",
    "# test:\n",
    "Fs=Fs_func(0.0,0.0,False)\n",
    "y=steady_states(Fs,alpha*DeltaT)\n",
    "\n",
    "i=0\n",
    "for Fs in Fs_range:\n",
    "    y=steady_states(Fs,alpha*DeltaT)\n",
    "    DeltaS_steady[:,i]=y/beta\n",
    "    for j in range(0,3):\n",
    "        q_steady[j,i]=q(DeltaT,DeltaS_steady[j,i])\n",
    "    i=i+1\n",
    "\n",
    "print(\"calculating dDeltaS/dt as function of DeltaS for stability analysis...\")\n",
    "# ------------------------------------------------------------------------------\n",
    "hysteresis=False\n",
    "time=0\n",
    "DeltaS_range=np.arange(-3,0,0.01)\n",
    "rhs=np.zeros(len(DeltaS_range))\n",
    "for i in range(0,len(DeltaS_range)):\n",
    "    DeltaS = DeltaS_range[i]\n",
    "    rhs[i]=rhs_S(time,DeltaS,time_max,DeltaT,hysteresis)\n",
    "    i=i+1\n",
    "rhs=rhs/np.std(rhs)\n",
    "    \n",
    "print(\"doing a hysteresis run...\")\n",
    "# --------------------------------\n",
    "hysteresis=True;\n",
    "y0=[0]\n",
    "teval=np.arange(0,time_max,time_max/1000)\n",
    "tspan=(teval[0],teval[-1])\n",
    "sol = solve_ivp(fun=lambda time,DeltaS: rhs_S(time,DeltaS,time_max,DeltaT,hysteresis) \\\n",
    "                ,vectorized=False,y0=y0,t_span=tspan,t_eval=teval)\n",
    "T=sol.t\n",
    "DeltaS=sol.y\n",
    "\n",
    "hysteresis=True;\n",
    "FWplot=np.zeros(len(T))\n",
    "qplot=np.zeros(len(T))\n",
    "i=0;\n",
    "for t in T:\n",
    "    FWplot[i]=Fs_func(t,time_max,hysteresis);\n",
    "    qplot[i]=q(DeltaT,DeltaS[0,i]);\n",
    "    i=i+1;\n",
    "    \n",
    "N=len(qplot); N2=int(np.floor(N/2));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a8bd8853",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plotting...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGFCAYAAABg02VjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAA0RRJREFUeJzs3XdcE/cbB/BPmLIEZCgoiFtxACou1ELVat174N6jtc7WukWt1llH694VcSF17z1x4xaVIQgOUJZs8v39cb+cBBJIICHrefeVl83l7vKc5r53z32XgDHGQAghhBBCCCHFoKfqAAghhBBCCCGajxILQgghhBBCSLFRYkEIIYQQQggpNkosCCGEEEIIIcVGiQUhhBBCCCGk2CixIIQQQgghhBQbJRaEEEIIIYSQYqPEghBCCCGEEFJslFgQQgghhBBCik3nEov379/D3t4e9+7dK3Td8ePHY/DgwQWuIxQKsXbtWtSrVw+mpqYoV64c2rZti7S0NH6dlJQUfP78udixy0IgEOCPP/5Q+H6FQiGioqLk3u7z589ITk5WeDyEKEruMiEyMhICgQD+/v5S1w8MDEStWrWQmppa4H5PnTqFli1bwsrKCtbW1mjUqBFCQkL4z7OzsxETE6Ow4yiIj48P2rRpo5R9R0ZGyr1Neno6Pnz4oIRoCCk+We8T7t69i7Jlyxb6Ww4ODka7du1gb28PCwsL1KtXD2fPnuU/L8r1dejQoahWrVqB6zx8+BCOjo44fvw4AODy5csQCAS4fv06AODkyZNwdHTEgwcP+G1K8n6FaCedSyx+/fVXtGrVCg0aNCh03Tlz5iAoKAjnz58vcJ158+Zh/vz5eP36NU6fPo0+ffrAxMQEABAdHQ1ra2ucPHlSYcegCu3bt8ewYcPk2ubmzZuwtbXFo0ePlBQVIcUnT5kAAD169IClpSUWL14sdZ0zZ86gY8eO6NKlCx49eoTg4GCMGzcOzs7O/Dp16tTB3Llzix2/Kv3555+oWrWq3NtZWVlhy5YtSoiIkOKTtUxo2LAhfHx8MHXqVKnrPH78GD4+Pqhbty5u3bqFkJAQzJo1C1WqVOHXKcr1VRbm5uaoUaMGbGxsJH5uZWWFGjVqoHTp0gC0536FqJaBqgMoSREREThw4ACePXsm0/p2dnaYPHkyFi1ahFatWklcZ+vWrRg2bBi6du0KAHB0dISbmxv/eU5ODrKzs4sdu6plZGTAwEC+n0tmZiYYY0qKiJDik7dMEPnzzz/RpUsXTJs2Debm5vk+3759Ozw8PDBlyhR+WfXq1cXWycjIKFrQaqSox6ANx060k7xlwh9//AFXV1fMnz8flSpVyvf5nj17YG5ujmXLlvHLKleuLLZOUa6vsqhatSouXrwo9fOmTZuKfa4t9ytEtXSqxuLw4cPw8PAQO6lnzZqFGjVqwNTUFLVq1cLy5cvFbob79OmDS5cuIT4+XuI+DQ0NERwcLPFk3LlzJ1xcXAAAAwYMgEAgwM6dOwEAqamp+OWXX1C1alXY29uja9euYk0Kzp49i5YtW6JMmTIoV64cevfujY8fP4rt/9SpU2jcuDEsLS3RtGlTsacMQqEQLi4u6Natm9g2V69ehUAgwKlTp/LFGxoairZt26J06dIoW7Ys/Pz8AHDNqy5duoRz585BIBDAx8en0Bj9/Pzg7e0NAGjevDkEAgEuX74MAIiLi8PgwYPh5OSE8uXLY+DAgWJVr7dv34aXlxfMzMxQoUIFbN68WeLfPSHFJalMALgL/fjx41GlShWULVsWgwcPxpcvX/jPvb29YWJiIrU209DQEG/evMH79+/zfSZqbhUREYEtW7ZAIBBg6NCh/Odr166Fh4cHLC0t8d1334l9x5s3b9CjRw84OjrC0tISLVu2zNdc49WrV+jYsSPKlCmD2rVrY9myZWJl2sCBA2Fvb4+srCx+WVpaGiwsLPD777/ni/fr168YMWIEypYtC0tLS3To0AGMMfj4+GDevHnIzs6GQCCAQCAoNEZRUwyAK3sFAgFfzmRnZ2PevHlwdXWFtbU1fvjhB7GmY+/fv0evXr1gbW0NW1tbjBo1SuLfPSHFIalMEDV5bty4MSwsLGBnZ4dFixYBAKpUqQJ3d3ccOXJE4v4MDQ3x+fNnqYlKUa6vuV28eBHe3t6wtLSEh4cHDh06xH92/fp1sWtvXv7+/hAIBIiMjJR6vyLvfQQhYDqka9eubPLkyWLL1qxZw+7cucOioqLY6tWrGQAWEBAgtk6ZMmXY4cOHJe5zy5YtTCAQMHd3dxYQEMCysrL4zzIzM9mjR48YALZx40aWkJDAMjMzmVAoZM2bN2dt2rRhjx49Ys+ePWNdunRhTk5OLD09nTHG2NmzZ9l///3HoqKi2NWrV1nFihVZ27Zt+X2fOHGCCQQC9vvvv7NXr16xI0eOMDc3NwaALVy4kDHG2KJFi5iBgQF7//49v93YsWOZs7Mzy8nJyXcsDRo0YG3btmUvX75kISEh7ObNm4wxxhISEljz5s3Zd999xxISElhKSkqhMaanp7Njx44xAOzkyZMsISGBZWdns7S0NFa1alU2YMAA9vLlS/bgwQPWtGlT1qhRI8YYY9nZ2cze3p4NGjSIvXnzhgUHB7NHjx7J8K9LiPzylgkREREMALOxsWFz585lL1++ZMeOHWMODg6sadOmYtt269aNTZo0SeJ+Q0JCmJmZGStbtixbtmwZ+/LlC/+ZUChkCQkJzNnZmQ0cOJAlJCSw1NRUxhhjM2fOZC4uLuz8+fMsPDyczZs3jxkaGrKHDx/y8W3YsIG9fv2aPXv2jLVt25bZ2Njw52RsbCyzsrJiP/74I3vw4AG7desW69KlCzMwMGCtW7dmjDF248YNBoAdOnSIj2nfvn1MIBCw0NDQfMcyZcoUVr58eXbr1i325s0bduLECcYYYykpKez3339n+vr6LCEhgSUkJBQaY3Z2NktISGAA2MyZM1lCQgJf5vXv3595eHiw4OBg9vr1azZ69GhmYWHBl189evRg7u7uLCQkhD1//pxduHBBxn9lQmQn6T6hV69ezN7envn7+7O3b9+yp0+fil2XJk2axLp16yZxf+/evWNly5ZlFhYWbM6cOSwmJkbs86JcXxljbMiQIczMzIw1btyYnTx5kr18+ZKNGTOGv+Yyxti1a9cYAHbp0iXGGGOXLl1iANi1a9cYY4zt3r2bAWARERFS71fkvY8gRKcSi0aNGrHly5cXuE6FChVYv379xJbVrl2brV+/Xuo2Z8+eZc2bN2cAmJOTEzt69Cj/mehGZffu3fyy//77j+nr67OIiAh+2YsXLxgAtm/fPonfMWvWLCYQCPiLcN26ddmPP/4ots6dO3fEEosPHz4wIyMj/pizsrKYra0t8/Pzk/gdVlZWbNSoURI/8/b25m9MpMkbY95CjDHGVq1axaysrFhycjK/7NSpUwwACw4OZl++fGECgYAtWrSowO8iRBHylgmi87V///5i6x06dIgBYBcvXuSX/fTTT6xPnz5S9/3q1SvWp08fZmhoyMzNzdkff/zBhEIh/7mLiwsbMWIE//7Dhw/M2NhYrKwQCoXMwcGBjR07VuJ3nDt3jgFgp0+fZowxNn78eGZra8vfnDDGWHJyMitXrpzY+evu7s46duzIv+/SpQvz8fGR+B1du3ZlderUYRkZGfk+mzdvHjMwMJD6dyApRsaYWDnFGGMPHz7MV1YkJiYyIyMjtmTJEj7mdu3a0c0MUaq8ZcLJkycZAHb+/Hmp2yxbtow1btxY6ufv379nI0eOZKampszIyIhNnDhR7HwqyvV1yJAhzMDAgD19+lRsPTc3N/bdd98xxuRLLBiTfL8i730EITrVFOrLly+wtLQUWxYREYFZs2ahadOmqFy5Mt6/f49Pnz6JrWNpaSnWDCKv1q1b4+rVq3jw4AGcnJzQpUsX3LlzR+r6wcHBYIzB3d0d1tbWsLa2RuPGjQEAYWFhAICsrCxs27YN7du3h6urK1avXg3GGOLi4pCeno6nT5+ibdu2YvvN29HM3t4ePXv2xLZt2wBwHUo/f/4stZPYsmXLsHXrVjRs2BD79u0rtH9EQTEWdOwpKSlwcnLij7137978sVtZWWH27NmYMWMGWrVqRVWtRKkklQkA8MMPP4i9b926NQCIDURQWLlQtWpV7N27FxEREejTpw9mzpyJv/76S+r6Dx8+REZGBsaNG8efG2XKlMGnT5/4cgEATp8+jf79+8Pd3R39+vUDAL7Munv3Lpo3bw4zMzN+fXNz83z9O8aNG4eTJ08iNjYWCQkJOHnyJEaOHCkxrpkzZyI2NhaVK1fG8uXLkZKSIvUYZIlRkuDgYABAhw4d+GOvWLEisrOz+WNfuHAhrly5AldXV2zatEmsKRchipK3TLh27RosLS3x/fffS92msLKgbNmy2LRpE6KjozFlyhSsXr26wA7fgGzXVycnJ7i6uopt17p1a4UOmCLvfQQhOtV529raGomJifz758+fo0WLFmjSpAl++eUXVKpUCRMmTMi3XUJCAqytrQvdv7u7O06cOIGyZcvi4MGD8PT0lLieqD3zgwcP+PbGuWMEgK5du+LRo0eYNm0aatSogTt37mDmzJkAuPbfQqEQhoaGhcY0btw4NG/eHLdv38bevXvx448/okKFChLXHTFiBHx8fLBs2TIMHDgQR48exe7du6Xuu6AYpUlLS0ONGjX44e9ys7OzA8D1z+jWrRsWL16M9u3bY8aMGVi4cGGhx0qIvPKWCSJ6euLPXIRCIQDA1NSUXyZrueDo6IgtW7bg6dOn2LdvHyZPnixxPdEQ1fv27UOtWrXEPhONMrds2TLMnTsXU6ZMQd++fWFoaIgff/xRbB+ylAu+vr749ddf8e+//8LW1hbm5ubo3r27xHUbNmyI169fY/369Vi8eDG2bt2KO3fuSOy0LkuMBR37tWvXYGFhIfaZ6H2HDh3w+vVrrFq1CpMnT4a/vz8uXLgAfX39Qo+XEFnlLRPS09ML3UbWssDa2hqLFi1CWFgY9u3bhzVr1khdV5bra95yCuA6YOcupxRBnvsIQnQqsXBwcMC7d+/49//++y8yMjIQFBTEX4wzMzPFtmGMISYmBo6OjjJ9h5mZGQwNDfnh20T7zT3mvZubGxITE5GUlIR69erl20doaChOnDiBgIAA9O3bFwDEOn5ZWlqiYsWKuHDhAsaNG8cvv3LlSr59eXl5oV69etizZw+OHj2KHTt2FBh/lSpVsGHDBpQpUwaLFy/GunXrULp0aRgaGuLr168yx1jQsZ8+fRqlSpVC2bJlpcbh7u6Offv2QV9fH+vXr6fEgihF3jJB5MaNGxgwYAD/XtSB2sPDg1/27t27fJ2+C2JlZSU2yIOhoaHYuVGvXj0IBAKEhoaiXbt2Evfx119/YdCgQViwYAEAiHVuBrjz6+zZs8jMzISRkREAbi6Zp0+fisVuZmaGwYMH4+DBg7Czs8OgQYNgbGxcYOzTp09HgwYN0LZtW1y8eBGdOnWCoaEhcnJykJGRwW9fWIySjl00kl5YWBg6d+4sNQ4HBwcsWbIEVapUwejRo/Hs2TPUrVtX6vqEyCtvmVC/fn0kJibi0qVL/IAkeb17907mewSASzBE9wgAinR9FX3v27dv+WGsGWO4ePGi2LkuD0nXbED++wii41TYDKvErVy5ku8kzBjXuRn/b0/44sULNnXqVGZkZMRatmzJr/P48WOmr6/PPn/+nG9/mZmZ7KeffmLnz59n0dHR7PHjx8zX15dZW1vzbRazsrKYubk5a9euHQsPD2epqaksIyOD1a5dm9WoUYOdPHmSxcbGspcvX/J9M6Kiopi+vj7r06cPe/36NQsMDGTlypVjAFhYWBhjjLFt27YxgUDA5s+fz968ecOOHDnCPDw8mLW1NZs3b55YnOvXr2dGRkasQoUKLDs7W+LfTU5ODvvnn3/4jtstW7ZkTk5OfJvw0aNHMzMzM3bz5k0WFxcnU4zR0dFMT0+PDR06lEVFRbHMzEwWHx/PypYtyxo3bsyuXLnC3r9/z548ecLOnj3LGGPsy5cvbPPmzSwsLIzdunWL1apVizVv3rxI/96EFCZvmSBqY1yxYkW2dOlS9vr1a3b06FHm6OjIunfvzq8nFAqZra0tO3LkiMT9LliwgB08eJCFhYWx8PBw9ueffzI9PT2xQSDatm3LnJ2dWUhICN+5e+TIkczCwoJt3ryZRUVFsYiICHbo0CG+XXX16tVZrVq12KNHj9jly5eZl5cXA8C2bdvGGGPs+fPnzMjIiPXs2ZM9efKE3bx5k7Vr1445Ozszb29vsRifP3/OBAIBMzIyYi9evJD6d7R7925279499ubNG/bzzz8zfX19fv2AgAAGgK1fv57FxsbKFCNjjNWoUYPVr1+fhYaGsqSkJMYYY23atGEODg5s//79LCYmhr1+/ZodPHiQ32bTpk3s8ePH7MWLF6xnz57MzMyM7zBOiKLkLROys7NZw4YNmYODA9u7dy9/rb9//z6/ToMGDdiaNWsk7u/vv/9mu3btYi9fvmRRUVFsy5YtzNjYWGz9olxfhwwZwsqXL8+8vLzY+fPn2cuXL9no0aOZgYEBH5u8fSwk3a+IyHIfQQhjOtZ5OywsjBkaGvIjn3z9+pX179+fWVpasho1arA//viDjRs3jrm5ufHbzJw5U2qnqvj4eDZ48GDm4uLCjI2NmbOzM+vduzd7/vy52Hrbtm1jtra2zNramt24cYPfdsyYMaxKlSrM2NiYOTo6st69e/PbbN++nVWsWJGVKVOGde/enR05coQBYA8ePBDbb+3atZmVlRVr27Ytu337Nqtfvz77/fffxb7//fv3DECBHaI/fvzIvL29maWlJStdujRr1aoVPxKN6O+ufv36zMTEhE2cOFHmGBcsWMAsLS2ZnZ0di4yMZIwx9vbtW+br68ucnJz4vzfR6DrPnz9njRs3Zubm5qxMmTKsa9euLDw8XGrchBRH3jJBlFg8evSI9enTh9nb2zMnJyc2ZcoU9vXrV367s2fPMktLS7FO0rnNnj2bubq6MjMzM2ZnZ8e8vb3ZqVOnxNa5c+cOq169OjMzM2MrV65kjHEX9qVLl7K6deuyUqVKMVtbW9aiRQv29u1bxhhjwcHBzMPDg5mbm7MWLVqw69evMzMzM/bXX3/x+7116xbz8vJiFhYWrEmTJmzXrl1s8uTJrEmTJvnirFmzJvvhhx8K/DsaPnw4s7e3ZyYmJszNzU3sZj8zM5N169aNmZqaMk9PT5ljPHHiBKtQoQKzsLBg+/fvZ4xx5fHvv//OatasyYyNjVnZsmVZmzZt2NevX1lmZibr2LEjK1OmDDMzM2NNmjQR60hPiKLkLRMYYywpKYlNnDiRVatWjRkbG7MKFSrw5+zLly+ZkZGR2GAsua1du5Z5eHiw0qVLMysrK9akSRO2Z8+efN8p7/V1yJAhbNiwYWzv3r3Mw8ODWVhYsGbNmrGrV6/y+5U3sWBM8v0KY7LdRxDCGGMCxnRrBrP+/fsjOzsb+/btK3TdDx8+oGrVqjh69KjUKlBNsGDBAqxatQrh4eFi1a+EEPnKBJFGjRqhffv2mDdvnvICU7IrV67A29sbV69ehZeXl6rDIURtyFMm9OrVC2ZmZlrfPIjuI4jMVJ3ZlLSYmBhma2vL7ty5U+i6P//8Mxs8eLDyg1KSzMxM9tdffzF9fX22a9cuVYdDiFqSp0xgjLGDBw+ymjVrijUT0DSnTp1itra2bNiwYaoOhRC1I2uZcOfOHWZvb88+fPhQQpGVPLqPIPLSuRoLXTJr1iwsW7YMixcvljoSDSFEt5w7dw5t2rSBr68vtm/fznfyJoSQvOg+gsiLEgstFhERAcYYKlWqpOpQCCFqIi0tDSEhIWjSpImqQyGEqDm6jyDyosSCEEIIIYQQUmw6NfM2IYQQQgghRDl0YoI8oVCImJgYWFhY5JvpmhBdxhhDcnIyHB0dJc7iqq2oTCBEOl0sF6hMIEQ6ecoEnUgsYmJi4OTkpOowCFFbUVFRqFChgqrDKDFUJhBSOF0qF6hMIKRwspQJGpNYpKSkYNSoUTh37hyqVauGLVu2oFatWjJta2FhAYD7C6Hxlwn5JikpCU5OTvw5om4yMzOxcOFC7N+/Hx8/foSnpyc2btwIFxcXZGdnY9KkSTh48CAcHBywevVqtGjRQqb9UplAiHTqXi4oA5UJhEgnT5mgMYnFkCFDYGZmhoiICKxduxbt2rXDmzdvYGBQ+CGIqjVLly5NBQYhEqhr1b+BgQGysrJw4cIFlClTBiNGjICvry9u3LiB33//Hc+fP0doaCiOHj2KDh06IDw8HDY2NoXul8oEQgqnruWCMlCZQEjhZCkTNGJUqMjISFSpUgWRkZEoX748GGOws7PDjh070LFjx0K3T0pKgqWlJRITE6nAICQXTTs3rl+/jubNmyMhIQHOzs44efIkmjVrBgBo0KABhg4dip9//rnQ/WjacRNSknTx/NDFYyZEVvKcHxpRYxEcHAwXFxeUL18eAJcxeXp64uHDhxITi4yMDGRkZPDvk5KSAAAP3z9Ey9ItSyZoQojCGRkZQU9PDw8fPkR6errYXAyNGjXCw4cPJW5HZQIhhGiXqKgo3L59Gw8ePMCzZ88QHR2NmJgYpKSkICsrCwBgaWkJa2trVKxYETVq1EDt2rXh5eWFmjVr6lSNXEnSiMTi3bt3sLe3F1tmY2OD2NhYiesvXrwYfn5++ZZ7d/CGgaEBDPQMkC3M5v+0amoFQU0BfqjyA2Y2monfxv0mcb8ODg7YuHEj/37JkiW4fv26xHXnzZuH+vXrAwAeP36MmTNnSlyvUaNGmDVrFv/+p59+QlRUlMR1t23bBltbWwDA8ePHxWLJrWfPnhg0aBD/vnPnzhLXMzExwb59+/j3GzZswIkTJySuO3nyZHh7ewMA3r59K/WpcI0aNbBs2TL+/YwZM/DkyROJ6/7111+oUqUKAODatWtYunSpxPXatGmD8ePH8+8HDBjA3xjmFRgYCENDQwBAQEAAAgICJK43bNgwdO3aFQCQmJiIgQMHSlxP2//Nt27dKnE9dRUYGIgmTZrg8+fPsLGxERudwsbGBg8ePJC4nbQyod2/7WBqbgoLYwskZyTn+7N86fLY3nU76jvUV9oxEUIIKRxjDDdu3MDBgwdx6tQpvHjxotBtUlNTERsbi2fPnuHkyZP8chsbG7Rr1w7du3dHu3btYGpqqszQdYpGNIVasWIFgoKCcO3aNX5Z//79YWdnh1WrVuVbX9LTSScnJ+B3AKUK/i5TQ1PYm9rTTQbRCZpU/f/ixQs0bNgQQUFBSEpKwoQJExAdHc1/PnPmTDx9+hT//fdfvm2LUyYYCYzgWNpRYplgb2aPli4tMee7OXC0cFTQkRKiWppULiiKLh6zpvj48SM2btyIHTt2ICwsTOp6AoEAZcuWhaWlJQwNDcEYQ0JCAuLj45Geni51OwsLCwwYMABjxoxBvXr1lHEIGk/rmkI5Ojri06dPYsvi4uKk/gCMjY1hbGycb3m1MtXwWfhZ4tPJT6mfkJadhtSsVEQkRgAA4tPj8/3ZdHNTqTcZlHgQohzJycno2bMnRowYgTZt2uDmzZuIj48HY4yvzo6Li4ODg4PE7aWVCa62rviQ/UFqjUVSZhIyhZkFlgkv4l9g/5P90BPoUeJBiJZLS0vD169f+Zpkojzh4eH4888/sXPnTrEHQwCgp6eHJk2a4Pvvv0f9+vXh5uYGJycnvsVCbowxxMbG4uXLl7h79y6uXbuGy5cvIzExEQB3fVm/fj3Wr1+PNm3aYO7cufDy8iqRY9RGGlFjERoailq1auH9+/ews7ODUCiEra0tDh06xDfPKYgsmVZMcgx+O/MbzoadhamhqdSbjCxhVoHfVdDTTUo8iLrRhKd0WVlZaN++PQwMDHD06FEYGBggMTERDg4OuHr1Kho0aAAAcHNzw6RJkzBkyJBC9ynrcZ9+fRpD/xsKYwNjied0WnYaUrNTC/wuAQSwMraixINoDE0oFxStsGN+8+YNxo0bhytXrmDkyJFYs2aNCqLUDZ8+fcIff/yBdevW8X0lAK5GolWrVhg4cCA6d+4MKyurIn9HZmYmLl68iP3792Pv3r1ITRUvx9u0aYOVK1eiTp06Rf4ObSJPmaARiQUAtG7dGpUrV8bq1auxcuVK7NmzB0+ePJGp842iCsnCbjKKk3ikZqXihyo/YEmbJXSDQUqMut9AMMYwcOBAvHnzBmfPnoW5uTn/2YgRIxAbG4t9+/YhKCgIkydPRkREBMzMzArdr6KO+37sfYw6MgoJGQlISEtQeOJBDyOIKqh7uaAMhR1zYmIibGxskJOTgxo1asjUvp/IRygUYtu2bfj111+RkJDAL7ewsMCIESMwYcIEVKxYUeHfm5iYiN27d2PFihUIDw/nl+vp6WHMmDFYuHAhrK2tFf69mkQrE4u4uDgMGjQIt2/fhoeHB7Zt2ybzLJklVUgqIvGQ1scjIycD27tsR9uqbZUWP9E96n4DcfDgQfTq1QumpqYwMjLil69btw5du3bF8OHDcfbsWVSrVg0bNmyQuX1sSR23IhIPehhBSpq6lwvKIMsxN2/enB+8IzIyEs7OziUZolZ78eIFRo8ejStXrvDLTExMMGnSJPz666/Fqp2QVVZWFvbs2QM/Pz+xBMPR0RFbt25Fu3btlB6DutLKxKI41KWQLCjxEPXxEDKh1O0NYYjSpUpTcwqiMOpybpQ0dTnuwhKPoj6MsDaxhmUpS2zqtIlqO4jc1OX8KEmyHPP8+fMxd+5cAMCWLVswfPjwkgxRKzHGsHnzZkyYMEGsg/XAgQOxZMkSqf3mlCk9PR2rVq3CwoUL8fXrV375iBEjsHLlSp2akV6EEos8NKGQLKyPR0xKDDJzMqVuL605BT3VJAXRhHNDGTTluIv7MKKUfimYGZpREysiF005PxRJlmO+desWmjZtCgDo3bu32NDdRH6JiYkYNWoU9u/fzy+rXLkyNm7ciNatW6swMk5UVBRGjBiBM2fO8Mtq1qyJQ4cOoVatWiqMrORRYpGHNhSS92PvY0jQEMQkxxSpOQU91SSSaMO5URTacNwFPYxgYEjISCgw6ZDWxIqaXRJtOD/kJcsxZ2dnw87ODgkJCShTpgw+fvwIfX39Eo5UOzx//hydOnXCmzdv+GXjxo3DsmXL1GpOCVGNypQpU5CSkgIAMDc3x/bt29GzZ08VR1dyKLHIQ9sLyYKaUxT1qSbdXOgGbT83pNH2445JjoHfJT9cibyCT18/yd3ESlKzS6r91B3afn5IIusx9+zZE4GBgQCA27dvw9PTs6RC1Brnzp1Dz549+eFeLS0tsW3bNnTv3l3FkUn36tUr9OjRA48fP+aXzZ49G35+fjoxgzclFnnoYiEpUpynmnRzof109dzQ1eMWKaiJVWHNLiXVflLTKu2ii+eHrMe8adMmjB49GgCwcOFCzJw5s6RC1AqbNm3CuHHjkJOTAwBwd3dHUFAQXFxcVBuYDL5+/YrRo0fD39+fXzZ06FBs3LhR4vwZ2oQSizx0sZCURUFPNenmQjfo6rmhq8ctC2nNLgur/ZTUtIoGltBMunh+yHrM4eHhqFy5MgCgZcuWuHz5ckmFqPH++OMPzJo1i3/fuXNn+Pv7iw0jru4YY1i9ejUmT54M0e1zu3btcODAAY06DnlRYpGHLhaSxaXImwtKONSXrp4bunrcxVFQ7WdBTaskDSxBCYd608XzQ55jrlatGl6/fg0DAwN8/vxZJ0cJkgdjDDNmzMCff/7JL5syZQqWLFmisX1UDh48iAEDBvAzgjdt2hSnTp3S2vOFEos8dLGQVJai3lxQwqGedPXc0NXjVhZpTasKGliCEg71pYvnhzzH/NNPP2HdunUAgCNHjqBTp04lEaJGYoxhwoQJWLt2Lb9s2bJlmDp1qgqjUowrV66gS5cu/GR+2pxcUGKRhy4Wkqog7eZCnoSD+nCULF09N3T1uEuatIEl5E046CFEydLF80OeYz58+DC6du0KABg/fjzWrFlTAhFqHsYYpkyZgr/++otftm7dOowdO1aFUSlWSEgIWrVqhfj4eABAkyZNcOrUKVhaWqo4MsWixCIPXSwk1UlREo7cfThoWFzl0dVzQ1ePW10UJeHI+xCCRq5THl08P+Q55qSkJJQpUwY5OTmoUaMGXrx4UUJRahY/Pz/MmzcPACAQCLBjxw4MGjRItUEpQd7konnz5jhz5gxMTExUHJniUGKRhy4WkppAUsJRUB+OvMPiUrOJ4tPVc0NXj1vdSUs4pD2EyDtyHT2EUAxdPD/kPeYWLVrg2rVrAICIiAhUrFhR2SFqlFWrVmHSpEn8e22fqfzRo0do1aoV4uLiAHAd0wMDA2FgYKDiyBSDEos8dLGQ1FSS+nAUNCyuPvRhbGjM125QsiEfXT03dPW4NZWkhxAFjVyX+yEENa+Uny6eH/Ie84IFCzBnzhwAwObNmzFixAhlh6gxAgIC4Ovry79fuXKlWJKhre7fv4/vvvuOn0hvxIgR2LRpk1bMc0GJRR66WEhqE2nD4goEAnxO/5wv4cjbRpuaTEinq+eGrh63NpE0cl1BDyFyN6+kBxAF08XzQ95jDg4ORpMmTQAAvXr1wv79+5Udoka4du0aWrVqhcxMLumfO3cu3xxKF5w/fx4//vgjsrK4GtY5c+bAz89PxVEVHyUWeehiIakLJNVuSGujnbfJBHUG5ejquaGrx63tJD2EkNa8kh5ASKeL54e8x5yTkwNbW1skJCTA2toanz590tihUxXl9evXaNKkCd/XQJue2Mtj//796Nu3Lz/Pxd69e9GnTx8VR1U8lFjkoYuFpK6S1EZbWpOJ3J1BdTXR0NVzQ1ePWxcV9QGELicaunh+FOWYe/bsicDAQABcDUajRo2UGaJa+/LlC5o0aYLQ0FAAQOvWrXHixAmtn5FampUrV2LKlCkAgFKlSuHq1ato2LChiqMqOkos8tDFQpJ8I6nJhKTOoLkTDV1pKqGr54auHjfhyPoAIneioUt9NXTx/CjKMW/atAmjR48GAMyfPx+zZ89WZohqSygUolOnTjhx4gQAwNXVFTdu3NC6IVflwRjDiBEjsG3bNgCAo6Mj7ty5A0dHzSw7KLHIQxcLSVKwvJ1BJSUauZtKaGuioavnhq4eN5Eu7wMIaTWdor4a2pxoaPL5kZWVhVWrVmH+/Pn4/fffMXPmTJm2K8oxR0REoFKlSgAALy8vfpQoXZN7WFkbGxvcvXsXLi4uKo1JHWRkZKBVq1a4fv06AG6Oi8uXL8PIyEjFkcmPEos8NLmQJCUjb6IhqalE7kRDW5pO6eq5oavHTWSXN9GQ1ldDlGhoU9MpTT4/li1bhosXLyIyMhK+vr5KTSwAoFatWnjx4gX09PQQFxcHa2vrooaukU6cOIGOHTuCMQY9PT2cPn0arVu3VnVYauPjx4/w9PTE27dvAQCTJk3CypUrVRyV/CixyEOTC0miGnmbSkhKNERNpzT5yaWunhu6etyk6PL21ZCUaIiaTmn6fBqafH6kpqbC1NQUPj4+aN26tdITi8mTJ/MzS2tDJ115REREwMPDAwkJCQCAxYsX4/fff1dtUGro7t278PLy4kfKOnToELp166biqORDiUUemlxIEvWQN9GQ1HTK1NAUjuaOGnVDoavnhq4eN1GcvImGpKZTovk0NK2GUxvOj8ISi4yMDGRkZPDvk5KS4OTkJPcxnzt3Dm3atAEADB48GDt27ChW3JoiOzsbLVu2xM2bNwEA3bp1Q2BgoM6NACWrdevW4aeffgIAlC5dGvfv30eVKlVUHJXs5CkTtGNKQEKUrL5DfdwdfZd/n7vplOjJZWpWKl5/eQ0A8NrqpZE3FIQQ2ThaOGJ3j938+9xNp0TzaaTnpCM9Jx3x6fFourkpHEs7alWTKU22ePFihcwv0KJFC5iZmeHr1684deoUhEIh9PT0FBCheluwYAGfVFSuXBk7duygpKIAY8eOxZUrV7Bv3z4kJSWhT58+uHHjhkb2tygM1VgQUky5n1xmC7PzTdBlJDCCs5WzWtZk6Oq5oavHTUpG3vk08tZwippMqeuDB204P0qqxgIAOnfujKNHjwLgmr00aNCg6IFrgKtXr8Lb2xtCoRD6+vq4fv06GjdurOqw1F5ycjIaNmzID8k7Y8YM/PHHHyqOSjbylAnan1YTomSiJ5cffv2Ax+MeY4THCNS0qQmbUjYw1DNEJsvE6y+vcS/2Hry2eqHssrIYeGggYpJjVB06IUQJHC0csbHTRjz/+TnipsXhaL+jcDBzgIulC4z0jZCFLMSnx+PRx0dourkpqq2phoabGuJ+7H1Vh64zjI2NUbp0abFXUbVv357//5MnTyoiPLX15csX9O/fH0Ih9/Bs/vz5lFTIyMLCAgEBATAw4BoL/fnnn1o5kpjcicX79++xaNEi+Pj4wMnJCSYmJrC1tUWdOnUwatQorT+pCCmItBuKMqXKQE+gh/ScdHxM/Yjdj3ej2tpqcFvvRjcThGi5tlXbImZqDMInhuPm8Juoa1dX6oMHxxWOOP36tKpDJnL48ccf+f8XzeWgrcaMGYOoqCgAwHfffYdp06apOCLNUr9+fSxYsAAAN//HwIEDkZSUpOKoFEuuxGLmzJlo0KABMjMzMX/+fNy7dw/JyckICwtDYGAgvLy8sHz5ctSvXx9PnjxRaKBHjhxBpUqV+E5ShGgC0Q1F7poMM0Mz6An0kJqVyj+xdP3bFWOOjaFaDEK0XH2H+ng07pHUBw+xKbHo5N9Jq2o2mzZtqnU3T7lVrFgRtWrVAsDNwP3582cVR6QcBw4cwP79+wEA1tbW+Pfff6Gvr6/iqDTPr7/+ihYtWgDgRtaaMGGCiiNSLLkSi1KlSuHFixeYN28ejI2NYW9vDwMDA5QuXRo1atTgR0RYvXo1wsPDFRbknTt3MG3aNLi6uipsn4SUpNw1GaHjQ9Gvdj/Ym9rzTyyfxz/Hpnub4LHBgxIMQnSEpAcPoqZSoprNGn/X0PgEIzg4GOnp6fmWJyUl4ddff1VBRIonag4lFApx5swZFUejeHFxcfyoRgCwfv16ODk5qTAizaWvr49du3bxze927NiBU6dOqTgqxSly520bGxvEx8eLLcvOzkb58uXx4cMHhQSXe7+MMSxatAjXrl3D2bNn5dpeGzqiEe0kGl0qMSORnydDAAHsTO3QrVY3pc/0ravnhq4eN1F/otGlwhLCxObJMDcyR9caXUtkvhxFnR/du3dHo0aNMHPmTISEhKBOnTpin8fGxqJChQrIyckpbsjFVtxjPn/+PD8x3MCBA7Fr1y5Fh6hSvr6+CAgIAEBDyyrKtm3bMHz4cACAs7Mznjx5AgsLCxVHJZlSh5tds2YNvn79irS0NCxevJhfnpOTg9u3b8PExET+iAsh6ugiK0mjPRCijkRPLEXzZLyIf4GvWV/xMfUjNt3bhCMvj+CY7zG1GzWGEKIcoqZSotHmDoceRmpWKlIyU7D78W4cCT2C/T33a8RwtRUrVsSxY8fAGIObmxtsbGzg5uYGNzc31KtXD48ePYKDg4Oqw1SI5s2bw9zcHCkpKVo37Ozhw4f5pMLa2hr//PMPJRUKMHToUOzZswfnz5/H27dvMX36dPz999+qDqvY5P7VR0ZG4uzZs8jJycHZs2f518WLF2FmZoYDBw4oI065LF68GJaWlvyLquuIuhPNk3Fl6BU0KNcAZoZmYGCITYmF11Yv6sxJiI4RjTb38ueX6Fe7H8yNzKEn0ENSRhK6BHTRiCaTf/31F65duwZDQ0PcunULmzdvhpeXF16/fg0/Pz8EBgZi6dKlqg5TIYyNjdGqVSsAwKdPn3D/vnYMyvHlyxeMHTuWf79q1SqtSQZVTSAQYNOmTTA1NQUA/PPPP7h69aqKoyq+IjeF+v7773HhwgVFx1MgPz8/mZpCKXJ8akJUQVSD8fjTY2TmZMJIYIS65eoqfB4MXW0SpKvHTTSXqAbj4POD/Azf5czLKaVGU9HnR3Z2NvT19dX6Kbcijnnjxo0YM2YMAO5+Zc6cOYoMUSVGjhyJLVu2AOD6kRw7dkyt/x010apVqzBp0iQAQPXq1RESEoJSpUqpOCpxJTKPxYULFyAUCnHz5k0EBgYiJCSkqLtCQEAA6tatK/E1bNgwufenyPGpCVEFUQ3Gkb5HUEq/FDJZJu7F3kOXgC5q/5SSEKJ4ohqMG8NvwMLIgq/RbLWzlVoOWd2sWTOMGzcOmzZtwv379yV23tY2uYed1Yah92/evMknFRYWFtiwYQMlFUowfvx4fi6Q0NBQjZk0T5oi11i8evUKnTp1Qnh4OOzt7fHhwwd4eHggMDAQFSpUUHScAGSvsciLnk4STXY/9j467umIT18/wUDfAB7lPHCw90GFdOBU1Lnx/v17bNu2DefOncOrV68QFxcHMzMzlCtXDs2aNUO3bt3ELrqqRmUC0WSiGs3n8c+RkZUBBwsHBI8MVlinbkWcH0uWLEFISAhCQkL4mYarVq0KNzc3uLu7w93dHW5ubmrTrEZRZUKdOnXw9OlTCAQCfPz4Eba2tgqMsuRkZ2fD09MTDx8+BMA9Vde2YVHVydOnT+Hh4YGsrCwYGRnh6dOnqFq1qqrD4pVIjcXo0aPRsmVLxMfHIyoqCnFxcahTp47YcGSEkOKr71Afd0fdhWd5T4ABd2Puwu+Sn6rD4uWe38bPz0+h89tkZWVh2bJlsLCwEHuKk52djfHjx8PBwQH169fXinaphMhKVKM5sN5A6OnpIS4tDj3391Sr2sxp06Zhz549ePr0KW7duoWyZcvCw8MDxsbG8Pf3R/v27VGhQgWULVtW1aEqlOgBCmNMo4edXbduHZ9UuLm50b2dktWuXRuTJ08GAGRmZmLChAko4nN/lStyYvHo0SOsWbMG5ubmAIDSpUtj7dq1uH79usKCI4RwHC0ccbD3Qdia2UIoFOLE6xNqcxORe36bFi1aSJzf5vz580Wa32bVqlW4ePEinJ2dxZb//vvveP78OUJDQzF16lR06NAh3/DXhGi7Od/NQUPHhhBAgEcfHmHFjRWqDkmiUaNG4Z9//sGePXuwc+dOPH78GMeOHYODgwNGjhyp6vAUSjSfBQAcO3ZMhZEUXWxsLGbPns2/X79+vdyjcxL5zZo1i2/xc+LECRw9elTFERVNkROLevXq5btJePfuXb4bAEWaO3eu3M2gCNEWjhaOONz3MJwsnWBhZIHN9zarOiQAwOzZs/mxt/fu3Ss2cEJuLVq0QKdOneTa908//YQTJ07A3t6eX/b161ds3rwZ8+fPh4WFBXx9fVGtWjV+OERCdIXogYN7OXeYGZmpbfv358+fo169emLL2rdvj3Xr1iE4OFhFUSlH8+bNYWlpCYDrZ5GVlaXiiOT366+/8sP0Dx8+HE2bNlVxRLrB3NwcK1Z8ezgwYcIEpKWlqTCioilyYjFv3jz06dMHe/fuxc2bN/Hff/+he/fuGDBgAG7evMm/CCGKU9+hPq4Pv44+tftgZAP1e9I3depUlCtXDmPHjlXIDYNoGL7cnj59ivT0dDRp0oRf1qhRI77aXpKMjAwkJSWJvQjRBqLkYlzDcZjcdLKqw5GocePG2LBhQ77ldevWxYMHD1QQkfIYGhryzaESEhI0rhXH1atX4e/vDwAoU6YM/vzzTxVHpFt69eqF77//HgAQERGhkX//RU4svL298eTJE/j6+sLLywvdu3fH06dPMXXqVHh5ecHLywvNmzdXZKyEEABIdgQuz+X+VDPR0dE4evQojI2N0aNHD7i6umLp0qWIjY1V2He8e/cONjY2YpNP2djYFPgdNLcN0WpqXCYAXHv9DRs2YMiQIXj06BGEQiHS09OxfPlymJmZqTo8hctdM6tJzVmEQiE/7CkALFq0SGM7n2sqgUCAv//+m296tmTJEkRERKg2KDkVObEQCoWFvnJychQZKyE6LyYG6N8fCAwENqtHS6h8mjdvjlWrViE6Ohrr16/H/v37FdpEUtKMtowxGBoaSt1m+vTpSExM5F9RUVEKi4cQVVu5Eti4kftTHdWqVQvBwcGIioqCu7s7TExMYGFhgW3btmHx4sWqDk/hfvzxR+jr6wPQrMTC398f9+7dA8A1dx8xYoSKI9JNtWrV4hO8jIwMzJgxQ8URyadY880rah4LQohs5s8Hbt8GjIwAde3zmJGRgaCgIPTt2xedO3dGTk6OQmfXdXR0RHx8vNiIGXFxcQUOW0lz2xBtdf8+sG8fkJUFqNsgMjNmzMDt27cBADVr1sT58+cRHh6OAwcOICgoCGFhYfD19VVxlIpnbW3Nt9h49eoVXr58qeKICpeamorp06fz71esWMEnR6TkzZw5k68tCggI4M8jTVDkxCI0NBSurq7w9vbGxIkT4enpicaNGyM6OlqR8RFC/u/+fSAgAMjM5BILRzVs9dC/f3/Y29vjl19+QcWKFXH9+nU8ePBArHq9uFxdXSEQCHD//rdJwW7dusVPMESIrrh/H2jVCvj0CTA3B6ZMUXVE4mJjY9GxY0c4ODhg1KhROHHiBMqVK4fOnTvzy7WVpjWHWrFiBd69ewcA6NixI1q3bq3iiHSbpaUl/Py+DSs/ZcoUjRl+tsiJxZgxY2geC0JKiOgGIj2du4H4+29VRySZoaEhDh06hLdv32LJkiWoU6eOwr/D0tISvr6+mDNnDlJSUvDvv/8iJiYGvXr1Uvh3EaKuYmKALl2ApCTA2JhrHqluDxu2b9+ODx8+YP/+/bCyssLkyZNha2uL7t27Y8eOHYiLi1N1iEqTO7E4cuSICiMpXExMDN9JWF9fH8uWLVNxRAQARo4ciRo1agAArl27hv/++0+1AcmI5rEgRI3FxACjRwPffQckJgKlSgHnzwP166s6Msl27NiBVq1aQSAQIDo6Gi9evFDKU5a1a9fC0tISlSpVwvr163H+/Hmt7ARKiCSnTwO1anE1Ffr6QN++6lsmCAQCtGjRAkuXLsWLFy9w+/ZtNGnSBJs3b0b58uXRsmVLLF++nH9ari2qV6+O6tWrAwCuX7+u1vPszJ49G6mpqQCAsWPHombNmiqOiADcg7rczYh/++03ZGZmqjAi2WjUPBaE6JKYGKBzZ2DLFiAlBdDTU88biMmTJ+ebUXvUqFGoWLEiXF1dUbduXURGRhbrOy5evIiZM2fy701MTLBnzx58+vQJN27cyDdGPiHaSPSgQVRTkZUFNGwIzJ2r6shkV6tWLfz222+4fv06oqOjMXjwYFy9elUr56Hp3LkzAK4/6smTJ1UcjWRPnz7F9u3bAXC1wXM16cekAzp16gRvb28AwOvXr7F+/XrVBiQDmseCEDUTEwMMGADUqAE8fAgIhVzzpxEj1PMGwt/fHxUrVuTf79q1C7t27cLu3bvx5csX/Pjjj5g6daoKIyREs4kSCg8P7kFDRgbX/MnXFzh4UP2aQIk8ffoUr1+/lvq5nZ0dhg8fjsOHD2tlGZG7OdSlgweBiAggLk6lvewZ40IQhTJ79hy+Vnn69Ok0vKyaEQgEWL58Of9+4cKFSEpKQlxqHCISIhCXGqd2fS8ErIgR5R3uUeLOBQK1GHI2KSkJlpaWSExMpNFgiNqKiQF++w04fBhITeUSCoEAKF+eW6aMmgpFnBt2dnb4+PEjBAIBGGOoVasW2rdvj5X/H/vy/fv3qF+/PmJiYhQZerFQmUA0QUwM4OcH/PcfdxMoFHLLS5cGLl5UXu2los6Ptm3bonbt2nxZAACHDx/Gnj17YG9vj0mTJqFy5cqKCLnYlFEmZMfFYbazM4anpaFq7g+qVAHGjwcGDwasrBTyXYVJSAB27gTWrwfEB6l6AWA97O1PITz8gcRJSYnq9e/fH3sO7QHcgDI/lsFnfOY/q2JdBeMbjcdg98GwKmWllO+X5/ygeSwIUbHcNRQBAVyzJ1EtRf/+QHCw+jV/yq1+/fr466+/EBcXh2nTpiEyMhKTJ3+bAdja2loj2oUSoi7y1lB8/MiVCaamQIMGyk0qFCkkJAQ9evTg3z9//hy9evXC9evXsXfvXjRu3FitHjgo1OnTMHBxwR9paciXOoWFAZMmARUqcB1mlB8KKlYEpk4F3N2B/fuBs2eBevUWAHgIYDmSkh7h6lVKKtRV69GtgckA2gGfhZ/FPgv7EoZJpyehwsoKOP1a+b+nwsiVWPj7+8tU5RIREYGrV68WOShCdMH9+0C9ekD16vkTigEDuKdK//6rvs0cRFauXInVq1ejbNmyWLNmDdatW4cKFSrwn1+4cAGurq4qjJAQzXD/Ptdnws1NckJx9Spw965mJBUAkJiYKDbL/a5du1C5cmVERkYiOjoa7u7u/GhEWuX0aaBDByAtDXqQcKPFGPdKS+PWU2Jycfo00LEj0KIFEBUF7N0L9OoF6OtfxKNHcwD0g5OTF3x8DNGxY4nkOUROp1+fxshLIwFDAALk+0Gx//+XlpWGDns6qDy5kCuxOH78OOrUqYOdO3fmG8EhKysL165dw8iRI+Ht7Q2BQKDQQAnRBqInkbVqAU2bAo8fA1+/amZCIVK7dm28evUK9+7dQ1RUFIYOHSr2ea1atbBmzRoVRUeI+jt9mjvfW7QA7t371uxJUxMKkQoVKiA2NpZ/f+7cOfTu3Rv6+vowNjbG9OnTcebMGYV+Z0pKCnx9fWFvbw8vLy88f/5cofsvVEIC0KMHlziI2q5JIxRy6/XowW2nhFB69wbatuWa05Urxy1njIkNhvHHH+Nx5Ige2rbl1ldCKKSIEtIT0GN/D+6hfiG31UIIwRhDj/09kJCeUCLxSSJXYrFnzx7s2rULN27cQIsWLWBubg4nJyfY2trC2toa8+bNg7u7O549e8bPOkkI+XbjULcu9yTyxQtuojtDQ8DeXjMTChF/f38YGRnB3d0ddnZ2+T53cXGBu7s71WQSkouoxtLWFujUCYiN5fpWCQSAjY1mJxQibdq04ftXREZG4sGDB2jTpg3/eZUqVRAVFaXQ7xwyZAgMDQ0RERGBzp07o127dsjOzlbodxRo585vneRkIRRy6+/apbRQtmwBDAy+LT927Bg/uE7t2rXh6+sLAwNg82alhUKKaOfDnUjNSoUQsv2ehBAiNSsVu0JU949Y5M7bAFdLER8fDxMTE1haWioyLoWijppEFe7fB4YM4WopRENDipQqBVhbA9u3c0+TVEUR54avry9CQkLw22+/oXXr1ihfvjz/WVZWFoKDg7Fz506cPXsWu3fvVouHDlQmEFW4fx8YNYp7IhwVxT1cEDEy4vr0tmwJzJmj2gcMijo/3r17Bw8PD5ibmyM9PR3GxsYICwvjWzTcvn0bnTp1wocPHxQSd2RkJKpUqYLIyEiUL18ejDHY2dlhx44d6NixY4HbKuSYGQOqVeP6UMhzayUQAJUrA69ecf+vAIxxNePu7lzzp2/LGRo2bIj79+8DAIKCgtC1a1f+8z59gJAQ4PlzhYVCiogxhmprqyHsSxgYZP89CSBAZevKeDX+lcJaD8lzfhgU+GkhDA0NUU5Ut0YIwenTwNCh3FCQMTH5bxycnQFLS2DTJs19CpnXnj17cO/ePWzatAl+fn74+PEjrK2tkZaWhvT0dDRp0gTdunXD6tWracQRonNyP2BIS+OeCIsYGnIjPJUvzz1k0JYyQaR8+fK4c+cO1qxZg4SEBPz8889iNzoXLlzgJ5FThODgYLi4uPAPNwQCATw9PfHw4cN8iUVGRgYyMjL490lJScUPID4eePNG/u0Y47b7/JmrrlKA+HiuFnzBAvHlJ06c4JMKDw8PdOnSRezzHj24zt0KDIUUUXxaPN58kf/3xMDw5ssbfE77DBvTkv9HLFZiQYiuy/0EMiEhf82Ett84iDRo0AAbN24EoDk1mYQog2jY6LNnuT4SeR8wCATcCKOlSqm+xrIkVKxYEStWrJD42bNnz9CzZ0+Ffde7d+9gb28vtszGxkasn4fI4sWL4efnJ/d3BAQESJ3M7+eOHfGD3Hv8ZkSfPviY5+HLvHnzUP//F47Hjx+L9Y3IrVGjRpg1axb/ftKk2QAWYOXK2fj33xB+ee7mqHPmzMGJEyf4shsAPn1yBzAfffqMwKBBLTFo0CD+M9GEf3mZmJhg3759/PsNGzbgxIkTEtedPHkyP+Hb27dv8fPPP0tcr0aNGli2bBn/fsaMGfkmYhX566+/UKVKFQDAtWvXxGarzq1NmzYYP348/37AgAFSE8rAwEAYGhoCKPjffNiwYXyNT2JiIgYOHChxPQcHB7G/5yVLluD69esS1xX9m6dkpkj8XFZ9BvWBaSb3e9q2bRs/R8nx48fFYgGAfv36oV+/fsX6PpFiNYXSFNTsgSiKaFz5K1eAT5/yP4EEuJoJR0duEit1v3HQ1XNDV4+bKEfumkpRuZC7ib3oAYO9vXo0dSpMSZ4fDx8+hLu7u0L2tWLFCgQFBeHatWv8sv79+8POzg6rVq0SW1dSjYWTk1PxjjkuDpDQz0yu7RVUTSAKZf9+bhQoADhz5gza/v+CVK9ePTx48CDfnGT793PNoRQYCimiuNQ42C0r+u8p7tc4hdVYKK0p1NOnT2FsbIyqVasWvjIhWiJ3UwbGuJqJ3DcNoieQenraXzNBCCm8phLgaivs7TXjAUNJS0xMhL+/P7Zu3YqQkBCFda52dHTEp0+fxJbFxcWhXr16+dY1NjaGsbGxQr6XZ2PDdZQpah+LMmUUGkqNGkBgIJdYMMbEamhmz54tcaLjwEBuOwWGQorIxsQGVayrFLmPRRkT1fwjyjUq1OTJk7Fu3TqxZYcPH0afPn0wfvx4hIWFKTQ4Qkpa7uFgbW2BSpW+DQsbH8+1OxUKuWYMNjbceqNGAU+ecE94QkJ0N6l48OABTYRHtJJoVLdKlbhyQTQs7Js3XLmQlcXVVLq4fBvl7dUrIDycK1MoqeBcuHABAwYMgIODA/z8/ODi4iLT3FiyatCgAV6/fs0nF0KhEHfu3EHjxo0V9h0FEgi4GbWL4pdfFNpbWiAAxo7lEoX374GLFy/ixo0bAABXV1d079493zaxscChQ8C4cdRxWx0IBAKMb1S039MvjX9R2bQPcjWFKleuHAIDA+Hl5QWAm0XTzc0N9vb2fJViSEgIHNWsjpeaPRBp8jZtklQjAXxrymBtrV2drxV5bujr60NfXx81atSAm5sb3N3d4ebmxpcR6oTKBFKQ3E2bkpMl10hoc02lIs+P6Oho7NixA9u3b8eHDx/QpUsX9O/fH23btuXvIXJychQUOdC6dWtUrlwZq1evxsqVK7Fnzx48efKk0JsshR1zQgI3o3be9nDS6OkBJiZAdDT3g1KghARuxu0WLYDk5Fa4cuUCAG7Ajbzt6bOzga5duSGOIyMVHgopooT0BFRYWQFpWWkyDTmrJ9CDiYEJoidHw6qUlcLiUFpTKGmzaD59+hTZ2dno2LEj/vzzT5oMi6itvE0YpCUSpUoBZmaAhQU1ZZDViRMnMGjQINSuXRt6enrw9/fHb7/9BoFAAHt7e7FEQ1GdxAgprrydraUlEqK+U8nJmtNXQtXat2+Pixcv4vvvv8f8+fPRtWtXmJmZ8Z8r44nq3r17MWjQIDg5OcHDwwOnTp0q2Se3VlZcNUGHDlzSUFByoafHZaiHDinlTt7Kiusz0aGDEDk5kwA8Q/XqpdG7d2+x9WJjgZEjuYT6+HFKKtSJVSkrBPYORIc9HaDH9ApOLhh3Th3qc0ihSYW85EosRLNoOjs7AxCfRVNfXx/Tp0/HuHHjFB6kv78//vrrL7x584YfYaJVq1YK/x6iXSTdMEjqbA18SyS0rUaiJE2YMAFbtmxBp06d+GWnT5/GTz/9hDFjxuDTp08ICQmBv78/JRZEJfLWUFpYSO5sDYgnEtpWI1FSTp06BV9fX0ycOBENGzYske+0tbWVOiJRiWnblrtD79Hj2wUnV+MQ0U9NYGICwaFDwA/FGUuq8FAaNvRDcPBkAG9hZRWNwEB9WFsDX75wOdChQ9w18vhxpYZCiqht1bY47nscPfb3QGoW93sS63MhBDcrdxZwoO8B/FBFtf+IciUWolk09+3bx8+iKZpVE1DOLJoAEBsbi61bt6JOnTrYtGkTunXrhvDwcNjQkAXk/+S5YcjdhIESCcWJjIxEnTp1xJa1bdsWS5cuxb59+8SGIySkJORt0pS3hjI+/tu6os7WlEgozvXr17Ft2zZ8//33cHBwQP/+/eHr66sbA8C0bcs1b9q1C1izRmx+izAAawD4/P03uin5Tv7Zs2cIDp4PYDWsrCYgIWEe+vT59nmNGsCKFcDgwdy1kKintlXbInpyNHaF7MKa4DVi81tY5Fgg+Vwy8BAItwkH6qouTgAAk0N0dDSzs7NjlSpVYg4ODszFxYUJhUL+8+DgYGZvby/PLuWWkZHBzMzM2H///SfzNomJiQwAS0xMVGJkpKTcu8dYgwaM1azJWJUqjNnYMFamDGN6eoxxtw7iL1NTxlxcuPVq1WJs9GjG3r1T9VGoB0WeGz4+PmzatGn5loeFhTELC4ti71+RqEzQLu/eMTZqFFcm2Nh8O98NDSWXCaVKfVvP3p6xAQOoTMhN0efH169f2datW5mXlxfT09NjjRs3ZmvWrGHnz59nenp6CvmO4lJqmSAUMhYXx1h4OLsUGMgAMACsT58+iv+uPIYMGcJ/38qVK3OHwuLiuNCIZhEKhSzuaxwL/xLO4r7GsadPn/L/xhUqVGCZmZkK/055zg+5aixKehZNSQwNDSEQCGBiYiJ1HaXMqElUIvdQrxYWBTdnAsT7RqSmctW6S5ZQO+iSsG7dOjRu3BhxcXGYNGkSatWqhezsbKxduxbW1taqDo9oCUlNHAuqicjdpIlqKFXD1NQUw4YNw7Bhw/Dy5Uts3boVixYtwocPH1Q2ck2JEgi4YQRtbNCsfHlYW1vjy5cvOHHiBDIyMhQ/7O3/RUdHw9/fHwBgZWWFESNG5A6FaCiBQAAbUxt+jgobVxt06tQJR48eRXR0NPbu3St1or4SociMZuDAgWzlypWK3GU+J0+eZGZmZuzz589S15k7dy6fveV+0dNJ9fXuHWP9+3NPD0VPG11cGDMykvzEUSDg1hHVWFStytVi3Lun6iPRLIp+SvfixQvWtm1bJhAImLGxMTMwMGAmJibM399fIftXFKqx0AynTjHm4CBeJpiZSa+dzF0TYWPDWL16VCYUhSLOj6ZNm7KxY8eyjRs3slu3brHU1FSxz7Ozs1lQUBDr1KlTccNViJIsEwYPHszflxw/flxp3zNlyhT+e2bOnKm07yGqd+XKFf7fum7dumKtiRRBnvNDITNviya72bJlC0JCQhQ6dFxuaWlp8PT0RNeuXbFw4UKp6yllRk2iMHnbPRfUHwL4NtSrqMaCRmRRHGUNuxodHY379+9DT08PDRo0gIODg8L2rQg03Kx6yTtam+hclzQ6k0juPhFUE6FYijg/lixZgpCQEISEhCA0NBQAULVqVX4oag8PD9SrV09tyoaSLBMOHz6Mrl27AgCGDx+OLVu2KPw7EhIS4OTkhJSUFBgbGyMyMhJly5ZV+PcQ9cAYQ9OmTREcHAyAGzihrQKHslTacLN5XbhwAdu2bcOhQ4dgYWGB5s2bIyQkRO79BAQEYNGiRRI/8/T0xLZt2wAAY8aMgampKebOnVvg/pQyoyaRm6QEQtLNgrQOlDTUq+aqUKECKlSooOowiJqRlkDkbd4orSkTNXHUHNOmTeP//969e+jSpQs8PDxgaGgIf39/zJgxAwKBAHZ2dnj//r0KIy15P/zwA0xNTZGamoojR44gJycH+vr6Cv2ODRs2ICUlBQAwZMgQSiq0nEAgwK+//oqePXsCAJYuXarQxEIu8laHREVFsQULFrDKlSszMzMz5uvry44fP86ys7PZ48ePldYRa+7cuczZ2Zm9K0IPO2r2oFySmisU1HES4Jo45V6fOlCqhq6eG7p63CVFNMCCqKmi6Fw3NZVeJggEjFlbU1MmdaDo86N+/fr5Blw5fvw4K1++vNo00SnpMqFHjx5805WLFy8qdN8ZGRmsXLlyDAATCATs1atXCt0/UU/Z2dmsatWq/O/qngILT6V13lbFZDcAsGnTJqxfvx6XLl1Su1m9dYWkDpOy1EDkfdpIQzkSoj2k1UoWVAORe7hnat6oG54/f4569eqJLWvfvj3WrVuHtWvXqigq1erRowcCAwMBAAcOHIC3t7fC9r1//36+Fqhbt266Mbwvgb6+PqZMmYKxY8cCAJYtW4aAgIASj0OuPhZ6enoFTnbz9OlT1KtXT6F9LN6/f48KFSpAX18fpqam/HJfX1/8888/Mu2D2lPLTtIoTIX1gQAogdBUunpu6OpxF4U8DxVyowRCcyn6/PDx8UGjRo2wZMkSseXh4eHw9PREXFxcsb+juEq6TEhOToa9vT3S09Nhb2+PmJgYhTSHYoyhUaNGuHv3LgDgypUraNGiRbH3SzRDWloaKlasiE+fPsHAwAAREREoX758sfertD4Wqpjsply5csjOzlba/nWRtHbOFhbcTURmJree6CljQX0gKIEgRPNJmmCyoIcKhdVKUgJBclu3bh2aNGmCDx8+YPLkyahTpw4yMzOxfPlysVYPusTCwgLt27fHoUOH8PHjR1y5cgU+Pj7F3u/Nmzf5pMLDwwPNmzcv9j6J5jAxMcHo0aOxcOFCZGdnY8OGDViwYEGJxiBXYtG0aVM0bdoUq1evxt69e7Ft2zb4+fnB09MT/fv3R+3atZUVJ5GTtJqHgpopiP7MOwoTdZgkRPMVVPMgbR4IeqhAFKFWrVoIDg7Gzz//DHd3dxgaGkIoFMLAwABbt25VdXgq06tXLxw6dAgA1xxKEYnF6tWr+f+fMGGCbswTQsSMGTMGf/75J7Kzs7Fx40bMnDkTpUqVKrHvL/Zws6LJbv79919+shtlDTdbVNra7EHWmgdJJDVToFGYdI+2nhuF0dbjlrfmIa/cE0zSQwXdpczz4+3bt3j48KHaDUWtijIhJSUFdnZ2fHOod+/ewcCg6IN1RkVFoVKlSsjJyYG9vT3evn1LI2TqqH79+mHv3r0AgO3bt2PIkCHF2l+JDTcLADVq1MDSpUuxePFiHD16lB8alhRfQTcJRa15oGYKhGi2gmoji1rzQPNAkJLi7OwMZ2dnVYehFszNzdGhQwcEBgbyzaG+//77Iu/vn3/+4R/sjh07lpIKHfbLL7/wicWaNWswePDgEqu9UsgEeepOnZ9OynOTIAnVPJDiUOdzQ5nU+bgLarIkS20kQDUPpHjU+fxQFlUd8/79+9GnTx8AXBOW9evXF2k/qampqFChAr58+QJDQ0O8ffsW5cqVU2SoRIPk7cR/9erVYvW3KdEaC1IweW4SJD1dBCTfJFDNAyGaS9owrYV1li6sNpJqHgjRLB06dICJiQnS0tIQGBiItWvXFqk51O7du/HlyxcAQN++fSmp0HECgQATJkzAwIEDAXB9b0qqIz8lFsVUWOJANwmEKE9KSgpGjRqFc+fOoVq1atiyZQtq1aql6rAKTBwKmvulsCZLVBtJiHYxMzNDhw4dcPDgQXz69AmXL19Gq1at5NoHYwxr1qzh30+YMEHRYRIN1KtXL0ydOhUfPnxAUFAQ3r59WyLNECmxKERBTZVkTRwAukkgRBmGDBkCMzMzREREYO3atWjXrh3evHlTrA6QhSnsYYKsiYOkYVqpyRIhuqd37944ePAgAG50KHkTi8uXL+Pp06cAgObNm6NBgwYKj5FoHmNjY4wZMwZ+fn7IycnB+vXrsXjxYqV/r073sSisc7Ss7ZkB6YkD3SQQdabJbakjIyNRpUoVREZGonz58mCMwc7ODjt27EDHjh0L3Lag4y7KwwRJpCUONEwrUXeaXC4UlSqP+evXr7Czs0NaWhpsbW0RGxsr18ORPn36YP/+/QCAgIAA9O3bV1mhEg3z/v17ODs7IysrC2XKlEF0dDRMTEzk3g/1sZCibpkofDWoAAs7EySnGRU6gkphTZUocSBEdYKDg+Hi4sLPKioQCODp6YmHDx/mSywyMjKQkZHBv09KSgIAuJcJRxKcYFHWBMkZJoX2e5KlFpISB0KIPMzMzNCxY0ccOHAAcXFxuHjxItq0aSPTth8+fODnwrC3t0f37t2VGSrRMOXKlUPv3r3h7++Pz58/4+DBg3y/C2XRqcTibU55IKc04qMZANELKIU0mCEVFnpfkSw0g0VZUyRnm1FTJULU2Lt372Bvby+2zMbGBrGxsfnWXbx4Mfz8/PItD8+ZBsAI8TEMgADx8dyf+hgAK/jAQu8rEoU5SNcfh2w9IxgYCJCTA9jZAa6uQKVKDti4cSO/vyVLluD69esAgHnzxL9r3rx5qP//TOPx48eYOXOmxONq1KgRZs2axb//6aefEBUVJXHdbdu2wdbWFgBw/PhxsVhy69mzJwYNGsS/79y5s8T1TExMsG/fPv79hg0bcOLECYnrTp48Gd7e3gC4uQl+/vlnievVqFEDy5Yt49/PmDEDT548kbjuX3/9hSpVqgAArl27hqVLl0pcr02bNhg/fjz/fsCAAXyymFdgYCAMDQ0BcE9zAwICJK43bNgwdO3aFQCQmJgo9eLr4CD93zwvTf03P3z4sMR1ifL07t0bBw4cAADs27dP5sRi69atyM7OBgAMHz4cRkZGSouRaKYxY8bA398fALBx40alJxY61RTKGQ/wFQ6wQDKSYQFrfIElkrEJo1AfD79tUFAVBQ3HRLSIJjd5WLFiBYKCgnDt2jV+Wf/+/WFnZ4dVq1aJrSupxsLJyQmVcB9JcOTLBAskIwMm2I6haIuz4l9YUBUFjbRAtIgmlwtHjhzBhAkTULVqVZw9e7bwDf5P1ceclpYGe3t7pKSkwNLSEu/fvy90tuScnBy+OahAIEBYWBhcXFxKJmCiMRhjqFOnDp49ewaAe8hRp04dufZBTaGkeHw5CaUnj5A8VXWyjXinivh46e0g4uOBFy+A/fvzTyBBbaQIKRGOjo749OmT2LK4uDjUq1cv37rGxsYSJ4t6+MIUpVfOkdLRyka8U0VqKhARwW0orWzw8pI8NjS1kSJE6e7cuYNp06bB1dUVmbJ0jlQjJiYm6N69O3bt2oXExEScPHkS3bp1K3Cb06dPIzIyEgDw448/UlJBJBIIBBgzZgx++eUXAFytxdq1a5X3fbpUYyHTk4jCem5KmvK6IPSUk6gxVT+lK47Q0FDUqlUL79+/h52dHYRCIWxtbXHo0CG+iY40ch23LMNAyTKbpUhBvbqp/SVRA5paLmRnZ4MxhkWLFuHatWsaVWMBcIlCu3btAHBN2URNo6Tp1KkTjh07BoCrqenUqZPSYySaKSEhAY6OjkhLS4OlpSXevXsHMzMzmbeX5/ygxKIo7t8HRo2SUvMh59AxgPQZ8CgBIUqmDhfT4mjdujUqV66M1atXY+XKldizZw+ePHkCgUBQ4HYKP25ZhpiTNA6tNIU1x6QEhCiRppcLfn5+hSYW0ppHqvKYs7OzUb58eXz8+BHGxsb4+PGj1FgiIyNRqVIlMMbg5OSE8PBw6Ovrl3DERJMMGzYM27dvB8D1zRk2bJjM21JTKGWrXx/4/zTpUsnzlDM9nXsV1PQKKLyZBfX9IDpm7969GDRoEJycnODh4YFTp04VmlQohaMjIKUTrZjCZs6TpzkmAHTqVHACQg8lCJFI2oAOqmRgYIA+ffpg7dq1yMjIQFBQEAYPHixx3c2bN0P0XHjUqFGUVJBCjRkzhk8sNm7cKFdiIQ+qsVAlWZ9yytPMQiAArKyk9/2gmw2Si9qeG0qm1sctS3NMeSbZAQqvFaUaEJKLWp8fMtDUGgsAuHXrFpo2bQqAG/3szJkz+dbJysqCk5MTPnz4AAMDA7x9+xYODg4lHSrRMIwx1K9fHw8fPgQA3L9/Hx4eHjJtSzUWmkLWp5yyJiCivh9fvnDbFbcGhG42CCl59esDjx4Vvp4sCYi8taKF1YDQwBREhQICArBo0SKJn3l6emLbtm0y70vagA6q1rhxY1SqVAnh4eE4f/483r9/j3Llyomtc+TIEXz48AEA0LVrV0oqiExEnbjHjBkDgKu12LBhg+K/h2ostIgsfT/krQEBZGvvTTcbGklnzo08dOa4ZX0oIW8NCFDwwBQ0CpZG0/TzQ5Yai7zU6ZhnzZqFP/74AwCwevVqfjQfkfbt2+PkyZMAgDNnzsg85wUhycnJcHR0REpKCszNzRETEwMLC4tCt6PO23moU4GhFlR9s0F9QdSGrp4bunrcBZK1CZY8A1MAhY+CRc0z1Y6mnx+anlg8e/YMtWvXBsDVYNy6dYv/LDo6GhUrVoRQKETFihURFhYGPT09VYVKNNCYMWP4iTU3b96MESNGFLoNNYUiBZO1CRYg/82GLGP9yzIPCDXFIqRkydoES5aBKXKPgpWZKVuZAMjWPJNqR4mWc3V1hZubG0JCQhAcHIw3b97wM9Lv3LkTwv8n9UOHDqWkgsht5MiRfGKxdetWmRILeVCNBVEMeW425JkHBJDeFMva+luzLrrhKBJdPTd09bhLnCyjYBWleSYgW+0o1YQUiS6eH+p2zMuWLcNvv/0GAJg/fz5mz54NoVCIatWqISwsDAKBAOHh4ahYsaKKIyWahjEGd3d3PPr/g6SnT5/C1dW1wG2oKVQe6lZg6DxZ+4IUpSkWINsNB9WGANDdc0NXj1ttydM8U96mWIBso2JRM02eLp4f6nbMUVFRcHZ2BgBUr14dL168wOXLl+Hj4wNA+ohRhMhizZo1mDBhAgBg6tSpWLZsWYHra11i8c8//+Cff/5BdHQ0KlWqhL/++gvff/+9zNurW4FB5FBQU6y8NRZFueGQpWO6Fj/51NVzQ1ePWyvIUztalJoQWYbs1vLO6bp4fqjjMX///fe4ePEiAODmzZtYt24d/v33XwDcHD59+vRRZXhEg8XHx8PR0RGZmZmwt7dHdHQ0DA0Npa6vdYnFqVOnUKdOHVhaWmLWrFnw9/fHhw8fZJ4QRh0LDKIE8txwFLU2RMuefOrquaGrx61z5KkJKUozTUC2zuka1lRTF88PdTzmHTt2YOjQoQAAX19fBAUFIS0tDdbW1oiJiUGpUqVUHCHRZL1798aBAwcAAEFBQejatavUdbUuscjt3r178PLywufPn2FqairTNupYYBA1IGvH9JJ48qmiRERXzw1dPW5SCHmaaYo6p2dlyfcdsjbVVGEioovnhzoec3JyMsqWLYu0tDSx5ePHj8eaNWtUFBXRFqdPn0a7du0AAB07dsTRo0elrqu1icXLly8xY8YM2NnZFTiph7rOqEk0WEk8+ZQnEVFQ0yx1vJiWBF09bqJgsnZOL2pTTUD2RESBTbN08fxQ12MeMGAA/P39xZY9ePAA7u7uqgmIaI2cnBy4uLggOjoaenp6iIqKgqOUhxhamVg0b94c169fR6NGjXDp0iWYmJhIXXfevHnw8/PLt1zdCgyixeR98lmURASQvWmWlM7q6noxVTZdPW6iQvI21SxqIiJr06wCakp18fxQ12M+c+YM2uYqsz08PHD//n0VRkS0yezZs7Fw4UIAwOLFi/H7779LXE8rEwsA+PLlCxYtWoTAwEDcu3cP1tbWEtejGguiceRNRIrSNAvI11k9ycwMlm/f6ty5oa43EYTw5E1Eito0S0JNaZKNDSxDQ3Xq/FDXMiEnJwfOzs6IiYkBAPz999/46aefVBwV0RZhYWFwc3NDr169MHbsWHh6ekpcT+MSi4CAACxatEjiZ56enti2bZvYskqVKmHGjBkYOXKkTPtX1wKDkCKTt2mWlM7qSQAsoXu1eVQmEK0kT9OsAmpKdbFcUOcyYcWKFZg6dSoqVKiAR48eSX2oSkhRpKWlFdgKCNDAmbf79euHfv36yby+mZkZBAKBEiMiRM3JM3u6iKTO6mZmwNu3SgmREFLC2rblzm9ZSasptbEBQkOVFSWR0+TJk9GkSRO4uLhQUkEUrrCkQl5qkVgUZufOnejcuTPMzMywdetWREdHo3379qoOixDNUr8+8P+ZNnlJSVwncEKI7qlfH7h7N/9yKhfUikAggJeXl6rDIEQmeqoOoDApKSk4duwYXF1d4ejoiAMHDuDcuXNSe64TQgghhBBCSp7a11iYm5vzE3gUlagbSVJSkiJCIkRriM4JNehqVaKoTCBEOl0sF6hMIEQ6ecoEtU8sFCE+Ph4A4OTkpOJICFFPycnJsNShpg9UJhBSOF0qF5KTkwFQmUBIQWQpE3QisShTpgwA4O3btzpTSCqSaLjeqKgotRstQxOo898fYwzJyck617SQyoTiUefftCZQ978/XSwXHB0dERUVBQsLC6mDw6j7v5sy0DHTMQPylQk6kVjo6XFdSSwtLXXmR6IMpUuXpr+/YlDXvz9dvLGmMkEx1PU3rSnU+e9P18oFPT09VKhQQaZ11fnfTVnomHVDQccsa5mg9p23CSGEEEIIIeqPEgtCCCGEEEJIselEYmFsbIy5c+fC2NhY1aFoJPr7Kx76+1M/9G9SPPT3Vzz096eZdPHfjY5ZNyjymAVMl8aTI4QQQgghhCiFTtRYEEIIIYQQQpSLEgtCCCGEEEJIsVFiQQghhBBCCCk2SiwIIYQQQgghxab1iUVKSgp8fX1hb28PLy8vPH/+XNUhaYzMzEzMmTMHNWvWRJkyZdC2bVtERESoOiyNExUVBVNTU/j5+ak6FAIqE4qDygTFoXJBc/3zzz9wdXVF6dKl4ebmhgsXLqg6JKXz9/dHw4YNYW1tDXd3d5w/f17VIZWII0eOoFKlSmjTpo2qQ1EKZVwPtT6xGDJkCAwNDREREYHOnTujXbt2yM7OVnVYGsHAwABZWVm4cOECYmJiYGdnB19fX1WHpXFmzJihU8PWqTsqE4qOygTFoXJBc1WpUgVnzpzBu3fv4O3tjd69eyMnJ0fVYSlVbGwstm7diri4OIwePRrdunVDfHy8qsNSqjt37mDatGlwdXVVdShKo4zroVYPNxsZGYkqVaogMjIS5cuXB2MMdnZ22LFjBzp27Kjq8DTO9evX0bx5cyQkJMg8tbuuu3v3Lnr06IHmzZujevXqmDt3rqpD0mlUJigWlQlFQ+WC9rh37x68vLzw+fNnmJqaqjqcEpGZmYkyZcrA398fXbp0UXU4SpOdnQ3GGBYtWoRr167h7Nmzqg5JoZR1PdTqGovg4GC4uLigfPnyAACBQABPT088fPhQtYFpKCMjI+jp6cHIyEjVoWiMKVOmYPr06fR3piaoTFAsKhOKhsoF7fDy5UssWrQIQ4YM0ZmkAgAMDQ0hEAhgYmKi6lCUysDAAIaGhqoOQ2mUdT3U6sTi3bt3sLe3F1tmY2OD2NhYFUWk2QIDA9GkSROtL0wUJSgoCPHx8Rg5cqSqQyH/R2WCYlGZID8qF7RD8+bNUbNmTURHR+Ovv/5SdTgl6vTp02CMwdPTU9WhkGJQ1vVQqxMLoVAIPT3xQ2SMaXUGqiwvXrzA33//jXnz5qk6FI2QlZWFadOmYeXKldDX11d1OOT/qExQHCoT5Eflgva4du0aPn/+jJYtW6J27dr48uWLqkMqEWlpaZg6dSomTpwIa2trVYdDikFZ10ODYm2t5hwdHfHp0yexZXFxcahXr56KItJMycnJ6NmzJ0aMGKG1IyMo2j///IPq1avjhx9+UHUoJBcqExSDyoSioXJBMwQEBGDRokUSP/P09MS2bdsAANbW1li2bBkOHjyIgwcPanQtlKzHPGbMGJiammpFvyBZj1lbKet6qNWdt0NDQ1GrVi28f/8ednZ2EAqFsLW1xaFDh+Dt7a3q8DRCVlYW2rdvDwMDAxw9ehQGBlqdiypMo0aNEBISwj8NyMrKgkAgQMuWLXVmmD51RGVC8VGZUHRULminOnXqYOLEiRgxYoSqQ1GqefPmYfv27bh58yYcHR1VHU6J8fPz08rO28q6Hmp1U6jq1avDx8cHM2fORFpaGhYvXgwHBwd89913qg5NIzDGMHToUKSkpODAgQN0AyGH27dvIyMjA2lpaUhLS8PAgQMxa9YsunlQMSoTiofKhOKhckE77Ny5E1++fEFmZibWr1+P6OhotG/fXtVhKdWmTZuwfv16nDp1SqeSCm2mrOuhVicWALB3715ER0fDyckJly5dwqlTpyAQCFQdlkYIDAyEv78/Hj16BCcnJ1hbW8Pa2hoBAQGqDo2QIqMyoeioTCC6LiUlBceOHYOrqyscHR1x4MABnDt3Tqtvtt+/f49x48YhISEBzZo148/7n376SdWhkWJSxvVQq5tCEUIIIYQQQkqG1tdYEEIIIYQQQpSPEgtCCCGEEEJIsVFiQQghhBBCCCk2SiwIIYQQQgghxUaJBSGEEEIIIaTYKLEghBBCCCGEFBslFoQQQgghhJBio8SCEEIIIYQQUmyUWBBCCCGEEEKKjRILHdayZUsIBAKxl76+PpKSkmTavl+/fvj333+L/P379+/H8OHDi7y9PObPn48ZM2aUyHcRosmoXCCEyOvGjRsQCARo166dxM+jo6MxduxYVK1aFaVKlULZsmXxww8/4PHjxxLXL245Ig2d88onYIwxVQdBSh5jDJaWlpg7dy769+/PL9fT04O9vX2h29+7dw8DBgzAkydPoK+vX6QY2rZti3nz5qFp06ZF2l4eSUlJqFGjBu7duwdHR0elfx8hmojKBUJIUYwYMQKpqakIDAzEq1ev4OzszH8WEREBT09PeHt7Y/LkyXBwcEBUVBQCAwMxbtw4VK9eXWxfiihHcps8eTKcnZ0xceJEOudLANVY6KhXr14hOTkZLVu2RLly5fiXLDcPALBx40b89ttvRT7p3759i3fv3pXIzQMAlC5dGiNHjsS2bdtK5PsI0URULhBC5PX161fs27cPEydOxPfff48dO3aIfb527VqYmZlh3759aNq0KVxcXNCiRQusWrUqX1IBFL8cyevBgwf4+vUrADrnSwIlFjrq3r17MDAwQL169Yq0/dmzZ9G6dWsAQFpaGgQCAY4ePYpmzZqhXLlymD9/Pl69eoVOnTqhTJkyGDp0KIRCIb/99u3bMWzYMABA9+7dMXjwYP6zr1+/wtTUFBcuXAAAhIeHo3379rCzs0OjRo1w6dIlft3x48fDyckJdnZ26NevH99cQyAQ4MiRI3B1dcXBgwcBAK1bt8a5c+eKdLyE6AIqFwgh8tq3bx/KlSuHRo0aoX///ti+fTtyN4b58uUL0tPT8fbtW5n2l7scAbjz9vLly/Dx8YGNjQ28vb3x8uVLmT738fHBtWvXsHDhQlhbW+PWrVt0zisbIzpp6tSpTCAQMDMzM/7VpEkTmbbNyMhgJiYm/PucnBwGgHl7e7OYmBh2+/ZtBoDVq1ePvX79moWHhzNzc3MWEBDAr1+tWjX28eNHxhhjJ06cYGZmZiwlJYUxxtjevXtZxYoVmVAoZOnp6axSpUps2bJlLDMzkwUGBjJzc3P2+fNnxhhjz58/Z1+/fmXp6enshx9+YHPmzGGMMQaANW3alL1+/ZqP8/3796x8+fLF/8sjREtRuUAIkVezZs3Y3LlzGWOMJScnM1NTU3b27Fn+83v37jFnZ2cmEAhYgwYN2LRp09jTp08l7itvOcIYd97Wrl2bPX78mH39+pWNGDGCVa9enWVlZcn0ube3N1u4cCG/PzrnlYsSCx3l4+PD+vTpw169esW/YmJiZNo2JiaGOTo6ii0DwI4dO8a/r1GjBps+fTr//vvvv2eTJ09mjDF2+vRp1rNnT/6znJwc5uzszHbt2sUYY6xHjx78jUBgYCArW7YsX0AwxsTWffjwIevXrx+rWrUqK1WqFPvxxx/5eHbu3CkWY2ZmJjMyMpLpGAnRRVQuEELk8eLFCwaAvXz5kl/Wr18/1q9fP7H1srOz2aVLl9jMmTNZtWrVmIGBATt06FC+/UkrR3bv3s2/T0pKYoaGhuzatWsyfZ43saBzXrmoKZSOevDgAVq0aIGqVavyLwcHBwBc+8b69eujTp068PX1zbetiYkJ0tPT8y23srLi/9/BwQFmZmb8+0qVKuHz588AgK1bt2LEiBH8Z3p6ehg2bBh27dqF1NRUnDp1im8C8ebNG8THx8POzg7W1tawtrbGp0+f8OnTJ4SHh8Pb2xuenp44efIkpk+fjrS0NH6/lStXFosvNTVVLCZCiDgqFwgh8ti6dSs8PT3F+kr0798fhw4dwpcvX/hl+vr6+O6777Bw4UI8ffoU9vb22LNnT779SStHXFxc+P+3sLCAvb09YmJiZP48NzrnlctA1QGQkhcWFoaEhATUr18/32dfvnzBpk2bcOfOHejr6yMhISHfOlZWVkhLS0NGRgaMjY1l+k6BQAAAiI+Px/379xEQECD2+bBhw/DHH39g165d8PT05C/+5cuXR/Xq1fH06dN8+/z777/h7OyMSZMmAYDYzUPu7xSJjY1FuXLlZIqXEF1D5QIhRB7Z2dnYtWsXfv/9d7Hlbdu2hYWFBfz9/fHzzz/n204oFCIjIwN2dnb5PpNWjnz48IH//6SkJHz48AFVqlSR6XOBQICcnBz+czrnlYtqLHTQvXv3oK+vDzc3t3yfGRgYID4+HtOmTcPTp0/Fnjbm1qRJEwQHB8v93bt370b//v2hpyf+03NyckLr1q0xc+ZMsTHsO3TogM+fP2PlypXIysoCAERFRQHgnkqEhYXh3bt3ePjwIbZs2SJWeOR169YtNGvWTO6YCdEFVC4QQuRx7NgxfPjwAXXq1MGTJ0/414sXL9CiRQts3boVAwcOxOLFixEcHIyIiAhcuHABHTp0AGMMkydPlrhfSeXI7NmzERoaitTUVEyaNAnu7u5wd3eX6XN7e3vcvn2bLyvonFcuSix00P3791GjRg2Ymprm+8zCwgKPHz+Gm5sbevbsiWPHjkncR/fu3bFr1y65v3v79u0YOnSoxM969OgBgUCAnj178sssLS1x7tw5nD59Gg4ODrC3t4evry9ycnLQt29ftG3bFnXq1MGvv/6KqVOnFjiJ165du9CtWze5YyZEF1C5QAiRx9atWwEAbdq0Qd26dcVeQUFBePjwIdzd3XHs2DF07NgRtWrVwpgxY1CzZk2EhISgatWqEvcrqRwZOnQo+vTpgwoVKuDjx48IDAwUexBR0OeTJ0/G48eP4ejoiLdv39I5r2yq7uRB1EtoaCj//yNHjmQHDx6UuF5KSgpzcnJib968Ucj3CoVC5uPjwxYsWKCQ/eV1+fJl5ubmxoRCoVL2T4g2o3KBEFJS8pYjAPiO2JIU9nludM4rH9VYEDELFixAjRo14OHhgVKlSknN6s3MzLB7926++UFxvHv3Dj179sTHjx8xZcqUYu9Pki9fvmDXrl352lcTQgpH5QIhpKQoshzJi8555RMwlmsWE0JUICgoCLt378a6detQtmxZVYdDCFEDVC4QQgCu8/W1a9fg5eVVpM9JyaLEghBCCCGEEFJs1BSKEEIIIYQQUmyUWBBCCCGEEEKKjRILQgghhBBCSLFRYkEIIYQQQggpNkosCCGEEEIIIcVGiQUhhBBCCCGk2CixIIQQQgghhBQbJRaEEEIIIYSQYqPEghBCCCGEEFJslFgQQgghhBBCio0SC0IIIYQQQkixUWJBCCGEEEIIKTZKLAghhBBCCCHFRokFIYQQQgghpNgosSCEEEIIIYQUGyUWhBBCCCGEkGKjxIIQQgghhBBSbJRYEEIIIYQQQoqNEgsFOnHiBAQCgcRXr1695N7f8ePH0bBhQzg5OeHAgQNKiJgQQgghhaHrMSGy0cnE4v3797C3t8e9e/ckfn758mUIBAJcv34dALBixQr4+PhAKBQWuN/79+8DAA4fPoybN2+KvVasWCFXjG/evEGPHj3g4+ODY8eOoWXLlnJtXxAfHx+0adNGYfvLKzIyUmn7loePjw8GDRqk6jCIBiisTMhNlvLg5MmTEAgE8PLykvj5sWPHIBAI0KJFC7HljDFs2LABrVq1go2NDSpXrowePXrg7t27Ur/r2rVr6NOnD2rXrg0TExM4OjqiTZs2yMrKkri+upyfJUGXjpUoVu4yIe/1OCwsTKZ7gqJQ9vUZUI/zgq7P2ksnE4tff/0VrVq1QoMGDWRa/+eff0ZYWBi2b99e4Hr3799H6dKl0alTJzRp0kTs5ezsLFeM+/btQ2ZmJmbPng03NzeULVtWru1V5c8//0TVqlVVHQYAoEqVKqhSpYqqwyAaQJ4yQZby4PPnz7CyskJwcDDCwsLyfb5jxw7Y29sjISGBX5aZmYl27dphwYIFGDJkCG7fvo2dO3fCysoKTZs2xaZNm/Lt56effsIPP/yAmjVrYvPmzQgPD8eZM2cwYsQIGBoa5ltfnc5PZdOlYyWKl7tMyHs9njhxokz3BOpIXc4Luj5rLwNVB1DSIiIicODAATx79kzmbYyNjeHn54cFCxZg+PDhUte7d+8e3NzcIBAIih1nZGQkKlSogNKlSxd7XyUpIyND1SHwtmzZouoQiAaQt0yQpTz4/PkzypcvD0tLS+zevRtz5szhP4uPj8fRo0cxduxYHDx4kF++ZMkS3LlzB48ePUKFChUAcBffFi1awNLSEhMnTkTbtm1RsWJFAMCmTZuwZcsW3Lp1Cx4eHvx+ypUrhzp16kiMS53OT2Ur6rFu3rwZa9asQWhoKCpUqIB58+bh7NmzCA4OxsuXLxUcJVFHecuEvNdjWe8JAPX7PalLGUDXZy3GdMyqVatYkyZN+PdCoZAtX76c1apVi1lbW7Mff/yRbdmyhQFg165d49dLSkpiBgYG7OHDhxL3GxcXxwCwcePGsaysLLGXUCiUK0Zvb28GgH95e3szxhjbsWMHa9CgATM3N2e1atVic+bMYRkZGWLburi4sI0bN7KffvqJ2drastevX0vcf+vWrdmuXbtYgwYNmIWFBevevTv7+PEjY4yxAQMGMDs7O5aZmclvk5qayszNzdm0adMYY4wFBwezZs2aMVNTU1a+fHm2adOmfHHn/XmtWbOGubu7s9KlS7OWLVuyc+fOyRR3QdtJikOkatWqbMiQIYWul5ekWFxcXNiIESPE1ps3bx4zMDDg3zdq1IhNnz6drVmzRuLfK1FPecsExgr+zTFWeHkwb9485unpyf78809WrVo1sc/Wrl3L3N3d2ebNm5mZmRljjLG0tDRWqlQpNmfOHIn7+/LlC7OwsGCTJ09mjDGWlZXFSpcuzX799VeZj1Pa+Snt3AsICGANGzZkFhYWzMnJiY0cOZIlJyfz+5Pl917Qeefi4sLWrl3L/v77b9agQQNmaWnJWrVqxZ48eSIWd2HlnqT4CyuLpJk0aRL/73DhwgW2YcMGZmNjwxwcHJivr6/Mf9dEs+UuE6RdjwsrAxgr2u+psOszY+p1jS6srJRWBuS+Phe0HtE8OpdYdO3alb84M8bYtGnTmJmZGdu2bRsLCwtjW7ZsYVZWVvkSC8YY8/DwYKtXr5a43zNnzuQ7YUWvggoeSVJSUpivry+rVKkSS0hIYCkpKWzevHnM1NSUbd26lUVERLDAwEBmY2PDunbtKrati4sLc3Z2ZlOmTGFRUVEsOzs73/69vb2Zubk58/X1ZY8ePWLnzp1jtra2rFu3bowxxm7cuMEAsEOHDvHb7Nu3jwkEAhYaGsqys7OZvb09GzRoEHvz5g0LDg5mjx49YikpKez3339n+vr6LCEhgSUkJPDbz5w5k7m4uLDz58+z8PBwNm/ePGZoaMj/3UiLu6DtpMUhIiq4ClsvL0mxyJJYNG/enNnZ2bF+/frxf692dnb5/o2IeslbJhT2WxUpqDz45Zdf2Pfff8+ePn3KALDg4GD+swYNGrDFixezPXv2MAAsMzOTPXjwgAFgly9flhqnt7c3++GHHxhjjD158oQBYKdPn5b5OKWdn9LOvYMHD7IzZ86w6OhoduLECWZlZcVGjx7N76+w33th552LiwurUqUK69atG7tz5w67c+cOa9asGbO0tGQfPnxgjDGZyj1J8RdUFklz+PBhBoAdPHhQbPkff/zBALCVK1fK/HdNNFvuMkHS9VikoDKgqL+nwq7PjKnPNbqwbQoqA3InFvJeo4l607nEolGjRmz58uWMMcY+fPjADA0N2ZIlS8TWWbZsmcTEokOHDvzTgLwWL17MALCgoCD+Iil65eTkyB3nkCFDWNWqVcXiXLp0qdg6AQEBDAC7efMmv0x08ud+kpGXt7c3q1ixothTvzFjxrBSpUrxsbq7u7OOHTvyn3fp0oX5+PgwxrinpwKBgC1atCjfvvPebIviNzY2Zrt37+aXCYVC5uDgwMaOHSs17sK2KygOxr4VXIWtl5ekWGRJLL777jvm7OzM0tPT+WVjx45lpUqVkpjgEfWQt0wo7LcqUlB50L9/f9a5c2fGGGNVqlRhEyZMYIwx9uzZMwaAvXnzhh09epQBYB8/fmQHDhxgANirV6+kxunr68tq167NGONu+gGwly9fynWsks5PWcoMxrinpA4ODvz7wn7vhZ13Li4uzMXFhX39+pVfFhsby4yNjdmcOXNkLvekxS/pWAvSoEED1rRp03zL9+/fzwCwK1euyLwvotlylwmMiV+PcyuoDCjq70mW6zNjqr9Gy7JNQXHkTizkvUYT9aZznbe/fPkCS0tLAMDDhw+RlZWFtm3biq0jrQOnpaUlvnz5IvGz+/fvw8zMDJ07d0bDhg3FXnp63F/zixcv0Lp1a5QpUwbW1tYYP368TDGL4vzxxx/Flnfo0AEA8o1k4+PjI7HjZm5Vq1aFkZER/75Ro0ZIT0/Hu3fvAADjxo3DyZMnERsbi4SEBJw8eRIjR44EAFhZWWH27NmYMWMGWrVqhVOnThUaf0ZGBsaNGwdra2tYW1ujTJky+PTpk1jH1rxxF7adrHHIG6+kWGQhEAhQrVo1GBsb88s8PT2Rnp6OmJgYufZFSk7eMkGW3ypQcHnw+fNnvj12r169cPjwYQBAYGAgmjVrhsqVK8PCwgIAkJCQAEdHRwAo8HcSHR2N6tWrAwCcnJz4ZYog6feenJyMtWvXonXr1qhZsyYOHTqET58+8Z8X9nuX5bzz8fGBqakp/17UP+TRo0dylXtFOV9z+/DhA+7du4e+ffvm+ywyMhJ6enpi/ViIdstdJhREWhlQ3N9TYddnQPXXaFm2UeY1mqgvnUssrK2tkZiYCABIS0sDAJkvSAkJCbC2tpb42f379+Hu7s4nEZL0798fw4YNQ3x8PCIjIzF48GCZvlc0bGROTo7YctH7UqVKiS2X5XjydjAXFWKi4fN8fX1hbm6Of//9F4cOHYK5uTm6d+/Or+/n54cHDx7A1tYW7du3x6xZs6R+l+jved++fXj48CH/ev36NXbt2iU1blm2kzUOeeKVFIu+vj5SU1PFljHG8m1X2N8rUT+SyoTCfqtAweXB58+f+c/69u2LiIgIPHr0CIcOHeKHWBTduCQkJKB27drQ09PDhQsXJO4vISEBDx48QN26dQGAX//SpUvFO/j/k3Tu+fj4YP369ejbty/Wr1+PoUOH5tuusN97YeedpPIyJycHpqamcpV7xUkqACAqKgoAl9jkdfz4cdSsWRPm5ubF+g6iOXKXCQWRVgYU9/cky3VE1ddoWbdR1jWaqC+dSywcHBz4rL9evXoAkO9ifuXKFYnbvnv3jn+ymFtiYiLCwsIKHaoyLCwMOTk5EAqFKF26NBo2bChTzB4eHtDT08OJEyfElh8/fhwA95RQ0czMzDB48GAcPHgQgYGBGDRokNiTSQBwd3fHvn37+BsPgCt4cnJyxEaeqFevHgQCAUJDQ1GxYkWxl729vdQYZN1OUhySyLqeJC4uLvlqhkJDQyUmF0Sz5C0TZP2tSisPAPHEws3NDTVr1sSqVavw7Nkz9O7dG8C3xEL0dHTMmDFYvXo1f1OS24IFC2BgYIBx48YB4M7P8ePHY/ny5XKNLCPp/JTk6tWruHfvHrZt24YRI0YUq0agoPPu9u3byM7O5t/Hxsbi+fPn8PDwKHa5J+1YJZ2zdnZ2ALha5dwOHjyIS5cuyVxWE+2Qu0woiLQyoCR+T6q+RsuzTUlco4n60LnhZr/77jvs3bsXAFCpUiUMGjQIs2fPhq2tLZo0aYKzZ89iz549ACB2wUtMTMSTJ0/g4+OTb5/3798HYwz169cv8LsDAgLwxx9/4Ndff8XAgQOxaNEimS7Wjo6OmDhxIhYuXAg7Ozv88MMPCA4OxoQJEzB8+HA+QVK0sWPHYu3atTA0NMTKlSv55QkJCTh48CBatWqFjx8/4uHDh3B1dQUAVK5cGYwxbN++HV27dkW5cuVQqVIljBgxArNmzYKpqSnatWuHnJwc3L9/H+3bt89XGIoUtl1aWprUOHIrKF5ZtWvXDr/99ht27tyJ5s2bY+/evTh79ixycnLAGFPIEMNENfKWCbL8VgsqDwDxxALgai3mzZuHnj178stz11gAwB9//IEbN26gSZMmWLJkCby8vBATE4OtW7ciICAAO3bsEJvPZvHixbh79y4aNGiAWbNmwcfHBy4uLvj06RMeP36Mfv365YtL0vkpiaiZ1r59+2Bvb4/Lly9j69at/IORgmpmRWQ57168eIFRo0Zh4sSJyMzMxMSJE2FnZ4cxY8bA3Ny8WOWepGP18fGBubk5jh49Kraus7MzGjRogJUrV6JcuXKoVq0azp8/j/379wMAJRY6JneZIE1BZUBJ/Z5UeY2WZRtZr72KuEYTNaKy3h0qEhYWxgwNDVloaChjjLH09HT266+/MmdnZ34IuDdv3jAA7NSpU/x2mzdvzjdspMjy5csZAJlHMYiIiGDOzs7s+PHjUtfJ21ksJyeHrVmzhtWtW5eZmpoyV1dXtmLFinwdFiV1Ms5LNJxdbrt372YAWEREhNjymjVr8iPRiDx//pw1btyYmZubszJlyrCuXbuy8PBwxhhjmZmZrFu3bszU1JR5enry22RlZbGlS5eyunXrslKlSjFbW1vWokUL9vbt2wLjLmi7guJg7FvnsMLWy0tSLFlZWWzcuHHMysqKVatWjU2cOJHvQCvqvCrP3ytRH3nLhMJ+q4wVXB4IhUKmp6fHduzYwS97/vw5A8AOHz7ML8vKymIA2IYNG/hlGRkZbP78+axJkybMwsKCVatWjfXt21dqJ22hUMi2bt3Kvv/+e+bg4MCMjY1Z+fLlWZcuXVhWVla+9SWdn9LOvUWLFrFy5cqxcuXKsaFDh7KtW7cyAOzLly+MscJ/74Wddy4uLmzo0KFsxowZzMXFhZUpU4b16NGDRUZG8uvIUu5Jiz/vsX79+pXp6+uzXbt2Sfy7DAsLY23btmXm5uasbNmybOzYsXxH2xs3bkjchminvGWCpM7bBZUBon0U5fck73VEldfowrYpKI7cnbflvUYT9aZziQVj3OgqvXv3lnn9tLQ05uTkxLZv317k7wwMDGRhYWGMMcYePHjAypYty79XV5cvX2YCgSDf6FiEaBt5ygRFlAdEtocginT69GlmaWkpNhdHYRYvXswMDAxYamqqEiMj6qigMqGoZYCif090jSbqSOf6WADA8uXLceHCBdy9e1em9f/55x9UqVJF5s7Wkly+fBlNmzaFubk5Bg0ahE2bNqFSpUpF3p+ynT59Gj169MDQoUPh5eWl6nAIUSp5ygRFlAek5J05cwZTpkyRqxP2nTt34OrqChMTEyVGRtRRQWVCUcsARf6e6BpN1JWAMep9SsSdO3cObdq0ga+vL7Zv3y427B0hhChCpUqV0Lp1a2zevFnVoUjl7OyMNm3aYOvWraoOhWgBRf2e6BpN1BklFiSftLQ0hISEoEmTJqoOhRBCCCG50DWaqDNKLAghhBBCCCHFppN9LAghhBBCCCGKpRPzWAiFQsTExMDCwoLmGyAkF8YYkpOT4ejoKNPcBNqCygRCpNPFcoHKBEKkk6dM0InEIiYmBk5OTqoOgxC1FRUVhQoVKqg6jBJDZQIhhdOlcoHKBEIKJ0uZoBOJhWgW2aioKJQuXVrF0RCiPpKSkuDk5MSfI7qCygRCpNPFcoHKBEKkk6dM0InEQlStWbp0aSowCJFAXav+MzMzsXDhQuzfvx8fP36Ep6cnNm7cCBcXF2RnZ2PSpEk4ePAgHBwcsHr1arRo0UKm/VKZQEjh1LFcoDKBENWRpUzQjcaTRK3E3IsFcnJUHQbRAAYGBsjKysKFCxcQExMDOzs7+Pr6AgB+//13PH/+HKGhoZg6dSo6dOiA+Ph4FUdMCFEmKhMIUW86MdxsUlISLC0tkZiYSE8iVGzv7w/x3dL2yGzfDRWP/g2o4RMxXaJp58b169fRvHlzJCQkwNnZGSdPnkSzZs0AAA0aNMDQoUPx888/F7ofWY571Spg/XqgVi3Axwfw9gbq1gV0pC8r0WGaVC6UZJnw7Bnw00/A27eAoyPg4gJUqvTtz0qVgAoVAAOdaAtCdIk8ZQL9/EmJCQwEDi55g954D73j64DF5YEZM1QdFtEgRkZG0NPTw8OHD5Geni42QVSjRo3w8OFDidtlZGQgIyODf5+UlFTg94SGAlOnchVroaHA4cPccnt7oFUroE0b7qUj/VoJUVslVSYIhUDfvsDjx9z7sDDg2rX86xkYcIlG5crcq0qVb/9fuTKg5nkaIcVGiQUpEVevAv37AxnogYPNV6P3tV+AmTO5xz5Dhqg6PKIhAgMD0aRJE3z+/Bk2NjZiw97Z2NjgwYMHErdbvHgx/Pz8ZP6eqlWBjRsBf3+gbVvg4kXuN/zxIxAQwL0AoGbNb0mGtzegQ31diSYKD+cy5QoVgNq1VR2NQpRUmaCnBwT8HY/dUx+iax9jhJdvjogI8K/wcO7PzEzg9WvuJYmtbf6EQ/T/5ctTjSjRfJRYEKV7+hTo3BnIyAC6dgV6HBwPzIgGli4FRowAypUD2rVTdZhEzb148QJ///03goKCkJSUlG8sbcYYDA0NJW47ffp0TJ48mX8vGuFCGj09YPhw7gUA06Zxv99bt4CzZ7nX3bvAixfca+1a7kll48bfEo1GjahJBFEzQUHAlCncU57du1UdTbGVZJkAALUTrmPxnS4APNH49u18nwuFwLt3XG1GWBjw5o34/8fFfXtJ2BxGRlxzqrwJR7Vq3P8bGxf+d0KIqtFljyhVdDSXMyQkAM2aAXv2APr6ABYvBmJiuItbz57ApUtAw4Yqjpaoq+TkZPTs2RMjRoxAmzZtcPPmTcTHx4Mxxo9SERcXBwcHB4nbGxsbw7iYV2VjY+C777jXwoXAly/AhQvAuXNcovHmDXD9OveaNw+wsgJ++AFo3547B8qWLdbXE1J86encn6VKqTYOBVBJmWBpyf2ZmCjxYz09wMmJe333Xf7Pk5K4mo28CUdY2LfajpcvuZekfbu4ANWrAzVqiP9JNR1EnVBiQZQmIQH48UcuuahZEzh6FDAx+f+HenrA1q3Ahw/cXVmHDsCNG9xjGUJyycrKQvfu3eHk5ITly5cDAFxdXSEQCHD//n00aNAAAHDr1i1MmjSpxOKytgZ69OBeAHfDcPYsl2icO8clHvv3cy8AaNCAOx86dOBqM+hGgJQ4UZ8CDX/0rbIywcqK+1NKYlGY0qUBNzfulVd2NnetzJtwvHnDtV5LTv722alT4tuamnK1GnkTjho1vuVChJQUSiyIUqSnc82enjzhulGcOgWUKZNnJSMjrkf3d98BDx5wjdlv3OB6yBICrinD0KFDkZKSgrNnz8Lg/22LLC0t4evrizlz5mDfvn0ICgpCTEwMevXqpbJYK1UCRo3iXjk5XFOHEyeAkyeBe/e+vRYu5M6Jbt2A7t2Bli2pyRQpIVpQY6HSMqGQGoviEHX6dnEBvv9e/DPGuGdwoaFcbUbuP9+8AVJTgZAQ7pWXvb3khKNyZe4STIii0eWMKJxQCAwaBFy+zD2hOXkSqFhRysoWFtzdV9OmXAnZoQPXU9bcvERjJuopMDAQ/v7+MDU1FWv/vG7dOqxduxbDhw9HpUqVUK1aNZw/fx5mZmYqjPYbfX3uJ920KbBgAfD+PXD69LdEIyYG+Ocf7mVjA3TpwtV8tGql8Q+TiTrTgsRCpWWCKLFIT+dqf0roZBUIuK6I5cpxDyJyy8r61ic/b9IRG8sNOPHxIzf4RG76+tzDEElNqxwdaSR4UnQ0jwVRKMaASZOA1asBQ0OupiLv0xeJQkO5Thjx8Vx7kcOHuR0QpdLVc0OVx52RAZw/z1XWHT7M/eRFLCyAjh25JKNdO0BN8iSiLcaM4YY7mz8fmD1b6mq6WC7IdMw5Od+qFz98UPva9aQk4NWr/AlHaCiQkiJ9OzMzyQlH9eo0XK6uonksiMqsWMElFQCwa5eMSQXAlVjHjnEbnDwJjB7N9cGgxyZEyxgbcx2627fn7vGuXuWSjKAgriZDNJxtqVJcctGjB5dsiJp3E1JkohoLqhYrGn19LvtPTuaaQ6l5YlG6NNe36/9dTniMcbUZeROOly+52o+vX7nWyZJG6i1Xjks0atbkRix2deX+LFuWLteEQ4kFURh/f+DXX7n/X7GCm0xILk2aAPv2cZ0ztm/nhrpYsEDRYRKiNgwMuFm9fXyANWuA4GDg0CEu0QgPB/77j3sZGnLNpLp355pNqfn9DFFX2dncn1QbXHQmJlxiIUrSNJBAwDV3cnTkyp7cMjO5DuKSmlZ9+MA163z/nmvqnFuZMlyCIUo26tThXnZ2JXdcRD1QYkEU4tw5YOhQ7v8nT+ZeRdKpE7BhA9cDduFCLrkYM0ZhcRKirvT0vvXLWLqU64gpSjKePeOaFZ46xZ0OLVpwSUb37jT7N5GDUMj9SUOSFZ2oKZQoSdMyRkZcbUTNmvk/S0z8lmg8e8a9nj7lukd+/szVvubty2FvzyUYdet+e9WuTc08tRklFqTYHj7kbnCysrhaimXLirnDkSO5WYb8/ICffuLqXrt2VUCkhGgGgQBwd+de8+dzk/AFBXFJxr173NPCy5eBCRO4zpwDBgC9elFzKVKInBzuT0osik5U25OVpdo4VMDSEvD05F65paVxZZQo0RC9wsK4juMXLnAvEYGAG1m+Xr1vw+/Wq8eNiEXNqTQfJRakWMLDub7WyclcleqOHQq6Zs2dyyUXW7YA/fpxvV2bNVPAjgnRPDVrAtOnc6/ISK4m49AhbjK+K1e41/jx3Az3Q4dyE/Pp66s6aqJ2RDUW9OMoOh1OLKQxMQE8PLhXbl+/As+fA48fc0PPP3rE/f+HD8Dr19zr0KFv65cuLZ5suLlxNRz8/FdEI1BiQYosLo7rXPr+PVcYBAUpsE+gQACsX8/t/NgxronU9euS62cJ0SEVK3Ijr02aBERFcbPZ//sv94TwwAHuVbEi15pw2DCuwo8QANQUShG0vCmUIpmZAQ0bcq/cPn7kkoxHj77Nv/HsGTeK1bVr3EtEX5/rs9GgAVC/Pvdyd6emVOqMEgtSJKmp3L1+aCjg7MwN5KTwGT4NDIC9e7leq8HB3AR6N29yPc4IIXByAqZNA377jWuSuGMHl2RERgIzZ3IVf127cv0yfHzoflLnUVOo4qMai2Kztwdat+ZeIpmZXN8NUaIREsKVaZ8+cbUcjx9z5RvAPXesWVM82fDwoKFw1QUlFkRu2dlc66RbtwBra65DqdLu9c3MuBqLZs24Abl//JFr96HwLIYQzSUQfGuK8OefXK3Fhg1cHn7wIPeqVo0bxXnwYMDWVtURE5WgGovioxoLpTAy+ta5e8AAbhljXIvo+/e5vmWiP2NjuSZWz58Du3dz6woEXCdxLy/udqFZM252ceqzUfI0rnSJioqCqakp/Pz8VB2KTmKM60995Ag3zv7Ro0CtWkr+UltbbtrismW5utPu3blZxggh+ZiYcDPf37jBPfUbN44bev/VK2DqVG6gtQEDuOYG2j89KhFDfSyKj2osSoxAwI1617kzN5bL0aPcXD8xMcDx49zAFl27cjW3jHG1Ghs2cOVf1aqAgwN3u7B8OVce0m1DydC4xGLGjBkwpsl9VGbhQmDTJu6BV0AA93SgRFSqBJw4AZibc8NLDB367SJJCJGoXj3gn3+4C/GmTVyNRmYmN+dMixbc6C7799PDV51BNRbFRzUWKufgwE0wOns217fz7VuujAsMBKZM4YbsNjLiOokHBXHza3l5cU2lvLy490FB3OdE8TSqdLl79y6uXLmC9u3bqzoUnbR1KzBnDvf/f/+tghFg69fnhpAwMOCymt9+K+EACNFM5ubcKM737gG3bwPDh3M1jvfuAX36cDPpbtjADRtJtBj1sSg+qrFQS3lrJxITuVrZpUu5exU7O+6hyo0b3Drdu3MDW1StytVwbNjA1XiIThFSdBpVukyZMgXTp0+HkZFRgetlZGQgKSlJ7EWK5/hxrn02wHUKHTtWRYG0aQNs28b9/4oVwF9/qSgQQjSPQMDVUmzZwj3lmzsXsLHhxpsfO5YbW37NGo2eVJgUhGosik9UY0GJhVorVSp/7cSrV8DOndy9TN26XHn45g034MXYsVwNb5ky3Dgxfn7cxL/Jyao+Es2jMaVLUFAQ4uPjMXLkyELXXbx4MSwtLfmXk5NTCUSovW7fBnr35jL5IUOABQtUHNDAgVwPVYCb4nvfPtXGQ/D+/XssWrQIPj4+cHJygomJCWxtbVGnTh2MGjUKJ0+eVHWIJA87O2DePG4EqTVruHbKsbHcpHtVqnCjPWdmqjpKolDUx6L4RDUW1BRKowgE4rUTjx4BX75wg8/MmcONUGVuzg15e+YMVza2acNNOuruziUoFy9SmSgLjUgssrKyMG3aNKxcuRL6MhSI06dPR2JiIv+KiooqgSi1U2go0KEDN7xsu3ZcO221GGXht9+4GcEArqS4eFG18eiwmTNnokGDBsjMzMT8+fNx7949JCcnIywsDIGBgfDy8sLy5ctRv359PHnyRNXhkjzMzLhT6fVrYONGbvjomBiu03fNmtw8GdSdSUtQU6jio6ZQWsPS8lvtxNmzQEIC8OAB1y+tf3+ua6dQyA2CsXw58P33XA1v9+7A5s1AdLSqj0A9aUTp8s8//6B69er44YcfZFrf2NgYpUuXFnsR+X34wCUTcXHcBDcHDnwrU1VOIOCaQfXsyT1C6NqVewRBSlypUqXw4sULzJs3D+/evYOVlRUMDAxQunRp1KhRA4MHD8b58+exevVqhIeHqzpcIoWRETepXmgo14eqXDkgPJy7wDZpwg0vTTQcNYUqPuq8rbX09bnaiXHjuGFsw8K4hyx793LDdJctC6SkcE2rRo3iannr1QN+/x24fJlyTREBY+o/4GCjRo0QEhICvf8XhllZWRAIBGjZsiXOnz9f6PZJSUmwtLREYmIiJRkySk4GvL25caOrVOE6PNnbqzoqCdLTuUcOV65wk2ncvMk9ciUyUfS5UaFCBXz9+hV9+/bFkCFD0LhxYwVEqXhUJhTu61dg1SpgyZJv7YwHDQKWLVPTsoAUrnlz4Pp1bhCMbt2krqaL54fMx9ynDzeU2po132rNiU4QCrkajRMnuEmBb90SH7K7dGmu+VT79txDWW2ay1eeMkEjHlvcvn0bGRkZSEtLQ1paGgYOHIhZs2bJlFQQ+WVlcRUB9+9z7bBPnVLjG4lSpYD//gNq1+YeLbRrB3z+rOqodFZ0dDSOHj0KY2Nj9OjRA66urli6dCliY2NVHRqRk5kZN1BDaCg3ujMA7NrFjSC1ZQvNgaGRqMai+KjGQmfp6XGzfc+ezT1s/fSJayo6YAA33VZSEjfk7fDh3HxBHh7AjBlAcLBulZdUuhAxjAEjRnCdl0xNudGgqlZVdVSFsLbmHh+UL89Nxdm5M42bqULNmzfHqlWrEB0djfXr12P//v1wplokjVWuHDcQW3Awd6FMSOCGru3ShbuwEg1CfSyKj/pYkP+zsQH69eNGlXr/nisj584FGjXiWms/fAgsXsw1Ja1UiesAfueO9icZGlm6bN++HXPnzlV1GFpp5kzuqaS+PnDwIDc0pUZwcuKqViwtuar+/v1pQGoVycjIQFBQEPr27YvOnTsjJycHS5cuVXVYpJgaNeJGiFu+nOuPcfQo1774zBlVR0ZkRjUWxUfDzRIJ9PW5MnLePC7BeP+eu5fq04er/Y2M5MrORo24h7WLFnGj8GkjKl0I7++/uewa4Jo6/PijauORW506wOHD3F1PUBDwyy/a/2hAzfTv3x/29vb45ZdfULFiRVy/fh0PHjzApEmTVB0aUQADA25m2zt3AFdX7uLZti0waRLNfaERaLjZ4qPhZokM7O25kfH37uVqdgMDuSTD1JTrFD5zJvc8tHt3rsGFNj0HVXhiQePZa6bAQO4+HAAWLuTmq9BI333HDecgEADr1n2b74KUCENDQxw6dAhv377FkiVLUKdOHVWHRJSgXj3g7l3g55+596tWAS1bciPJETVGTaGKj2osiJxMTLgEYu9e4ONHbpI+Ly/udAwK4jp7V67MzRH27p2qoy0+hZYuNJ69Zrp6lWs5xBgwZgzX2Uij9erF3ekA3MHs3KnScHTJjh070KpVKwgEAkRHR+PFixfQgIHnSBGYmABr1wLHjnGz1d65AzRtCrx8qerIiFTUFKr4qMaCFIOZGTe63rVrwJMn3ISk1tbA27fcRH3Ozlw30cuXVR1p0Sm0dMk9nv2nT59gb29P49mruadPuR9xRgY3FcTff6vJBHjF9csv3CR6ADdEw6lTqo1Hi02ePDnfg4JRo0ahYsWKcHV1Rd26dREZGami6IiydejAjfJcuTI370WzZsC9e6qOikhEiUXxUedtoiC1a3PPQN+94zqAt2zJnaJHj3LD/fv4AJcuqTjIIlBo6fK/9u47LKpr6wPwb+gdAUEGHAUbggXFrhjFxBZjRxNr7DFGzbXEmmuLkcRYo7EkxiiK2IgxBvV+ltjQqBGx1ygoiApIFUTK+f5YmYGRImVmzpT1Ps95ph3mLIaZzayz9177v//9L2xtbRXXo6Oji92vffv26NmzpyoPzSogNpaqs6ak0JeBHTv0bOhtUBDVgcvLo/q5f/8tdkR6KSQkBDVr1lTcDg4ORnBwMLZv347k5GR0794d06dPFzFCpm716lFy0bIlVXvu3JnqvTMtw3MsKo/LzTIVs7SkryonTwI3b9LIEVNTSioCAijJ0KUEQ22nLb755hsMHToUISEhOHv2LM6dO6fYmPhSUmhydmwsUL8+ZciWlmJHpWJGRsDPP9O3nJcv6dTqP/+IHZXeyc/Ph42NDQBAEAQsWbIEEyZMwKBBg2Bvb49p06YhIiJC5CiZurm4AEeP0nCo5GRqXx49EjsqpoTnWFQe91gwNfL2BtavB+7fBz79lGrRnDxJCUbfvtQrrO3U1rr07t0bZ8+exbBhw+Dv74927dqhXbt28Pf3V9chWRnJhz1dvw5IpTRKyNFR7KjUxMyMZqY3bUqzprp2pUumMn5+fli5ciUSExMxc+ZMxMTEYOrUqYrHHRwc8Pr1axEjZJpia0sVTho3poncH3wAZGSIHRVT4KFQlceTt5kG1KhB9WfkCYaxMa0F7ONDpWq1uYqU2lqX/Pz8Yrc8bX41DEB+Pk0cOnmSlp8/dAgoNIpFP9naAgcPAh4e1GPB33ZUasWKFVi9ejWqVauG77//HuvWrUP16tUVjx8/fhw+Pj4iRsg0yd6eekBdXYFr1wAeBadFOLGoPJ68zTRIJqME48oV6rV49YpK1XbooL29F2prXU7o0oAwAyEIwNSpwO7d1Dbu2wf4+oodlYa4ulLXjJMTla8ZOJDPOKlIgwYNcO/ePVy6dAmPHz/GyJEjlR739vbG999/L1J0TAw1agAhIXR940aqHMW0gPzEHs+xqDjusWAiaNAAOHaMilza2tI6wL6+tHSXtlFbYtGjRw/UqlULCxcuLHESN9Os5cuB1avpenAw0KmTuPFonJcXfcOxtKSumvHjeQE9FQgJCYGZmRmaNGkCZ2fnIo97eHigSZMmiI6OxunTp0WIkImhUyc6kQEAkybxAnpagXssKo97LJhIJBIacXL1KuDvD6Sn07yLlSvFjkyZ2lqXhIQEBAUF4fr162jUqBHeffddhIaG8lhrkYSEAF98QdeXLwc++kjceETTujWwaxf9Y928GZg/X+yIdF54eDgaNGiArVu3Iu6N1X1ycnJw5swZjB07Fh07doREL2oZs7L66ivA3R2IjgZ++EHsaJimEwv5grmdOnXSnwVzuceCiczDAzh+HPjkk4KRKN9+K3ZUBdTWulhZWeHDDz/Enj178Pz5cwwcOBCTJk1C9erVERQUhOzsbHUdmr3h6FFAPjplypSCs4gGq2dPYMMGuv7VVzRWg1XYjh07EBwcjLNnz6J9+/awsbGBTCZD1apV4eDggAULFqBJkya4efNmuYs35OTk4LvvvoOtrS2+/vprxf25ubmYNGkSpFIp/Pz8uCdES1lZAQsX0vVVq/gkr+g0WG628IK5CxcuVOmCuaK2C9xjwbSAqSlVj1q8mG7PmkVFMLWBiTqfPCsrC7/++iuCg4Nx4sQJvPfee+jfvz9CQkIQGhqKY8eOFTt0gqlOVBQtJZ+TA3z4IbBsmdgRaYmxY2lVmoULgQkTgGrVqFQWq5BmzZph478JWk5ODpKSkmBpaQl7e/tKPe+qVavw559/okaNGkr3z5o1C7du3cLdu3dx4MAB9OjRAw8fPoSTk1OljsdUb+hQYM4cKm3922+0pAwTiQbLzcoXzJWvbSVnZ2entGju6dOn8fDhQzRs2LDMzy1qu8DlZpmWkEhoIvfLl7Rs16efAg0bAq1aiRyYoCYff/yxYGtrK9SsWVNYuHCh8PjxY6XH+/XrJwQGBqrr8EpSU1MFAEJqaqpGjqctHj4UBFdXQQAEISBAEF69EjsiLZOfLwhjxtALZGEhCBERYkekcdr+2Xj58qUgCILQsWNHYfHixYIgCEJGRoZgZ2cnRBT6e/n5+Qlr1qwp8/Nq+++tb+bMoY9Z375iR2LgqlalP8T166XupurPR2hoqPBKhf+A1NEulPl33raNXsP33qtY8IypWH6+IPTrR2/L2rUFITNT9ccoT5ugttMW6enp2L17Nx4+fIh58+YplZ8EgDVr1mDIkCHqOrzBS0qiVbWfPqWa8vv2AebmYkelZSQS6kv84AOaWdqzJ3D7tthRsUKsrKyK3Hfjxg28evUKrVu3VtzXsmVLREVFlfg82dnZSEtLU9qY5gwcSJeHD9PZNSYSkSZvT58+Ha6urvj0009x/vz5Sj+fKtqFCrcJPBSKaRmJhKaMurlRRf0lS8SNR6Wty9ixYxWNRlhYGLp161biZE03Nzf04aEnapGZSd+R79yhso+HDlFteVYMExNg507qO3zxghbQe/JE7KhYKeLi4uDk5ASjQl+OnJycEB8fX+LPBAUFwd7eXrHJZDJNhMr+1bgx4OkJZGXRGjpMJCKVm42NjcWBAwdgbm6O/v37w8fHB0uXLi31M1te5W0XKtwm8ORtpoXs7QF5Vffly4GEBPFiUWli0aVLF8yfPx+NGzfGqlWrkJSUpMqnZ2WQmwsMGgScOwc4ONAZQjc3saPSctbWtKJX3brAo0fA++8DqaliR6VzLl++rJGqb/n5+UpfHgBAEASYys8kFmP27NlITU1VbI8fP1Z3mKwQiQRo356u//WXuLEYNBHLzfr7+2PVqlWIjY3F+vXrsXv37iJzJCqjvO1ChdsE7rFgWqpfP6B5czqBs3ateHGotHUZMGAADh8+jPDwcKSlpaFt27b48MMPceTIEVUehpVAEIDPPgN+/x2wsKDvyt7eYkelI5ydKQurVo2WuOzXD+DSyOXSvHlz2NjYoFGjRhg6dCiWLVuGI0eO4Pnz5yo9jpubG5KSkiAUWoMkMTERUqm0xJ8xNzdXTBqVb0yz2rShS04sRCRiYpGdnY19+/bho48+Qq9evZCXl4elS5eq7PnL2y5UuE3gHgumpSSSgqqfwcEFH3dNU0vrIpPJMG/ePNy+fRtjxozBzz//DC8vLyxatIjPFKrR4sXAjz/S/4wdO4B27cSOSMfUqgUcPAjY2FCR6BEjxPtk6qCDBw/CwcEBDRo0gJGREUJCQtC1a1dIpVJIpVJ0794ds2bNQmhoaKWO4+PjA4lEgsjISMV9f/31F1qJXgqDlaZpU7q8eVPcOAyaSInFkCFD4OLigsmTJ6NmzZqIiIjA5cuXMWXKFJUdQ2PtAvdYMC3WuzetzB0dTSNXxKDW1kUikaBz587YuXMnzp07hypVqqB3797qPKTB+vlnYN48ur52La3GyCrAzw8IC6OzUqGhwMyZYkekMz7//HNs2rQJO3fuRHBwMC5fvoxDhw7B09MT06ZNQ+PGjXHlyhXMmDGjUsext7fH4MGDMW/ePGRkZGDbtm148uQJBgwYoKLfhKmDlxddxsXRirFMBCLNsTA1NcWvv/6KR48e4dtvvy1Xadmy0li7wD0WTItZWQE9etD1//s/cWJQ6zoWu3fvxsWLF/HixQvFfU3lp62YyoSH0wqMANWL//RTcePReV26UImF4cNp4Q93d+A//xE7Kq0XExNT5AtD165dsXTpUuzatQu7du1S2bHWrFmD0aNHw9PTE3Xr1sWxY8dgbW2tsudnqlelCuDkRBXrHj6kCd1Mw0TqsdiyZYviemxsLDIyMuDl5VVicZeK0ki7wOtYMC0XEEA1aU6cEOf4akssRo0ahb1796Jly5bw9PRUNCCqbkgM3YULVMoxLw/4+OOCVRhZJQ0bRtWhZs2i5cqlUlphkJWoTZs22LhxI7755hul+5s2bYoRI0ZU6rn//PNPpduWlpbYsWNHpZ6TaZ6bGyUW8fGcWIhCg4nF1KlTMWrUKKWTDePGjcPPP/8MQRDg4+OD8PBw1KxZs8LHEKVd4KFQTMsFBNDlX39RJX0LC80eX22Jxf79+3Ht2rVKNRqsdPfuUZdXZiatWfHTTzR5h6nIjBk0bmPNGuq9cHEp+MSyItatW4dWrVohMTERU6ZMgbe3N3Jzc7FmzRo4ODiIHR7TAlIpcO0aJRZMwwShILHQwFCokJAQLFy4UHE7ODgYwcHB2L59O95//30sXrwY06dPx549e9Qei0rxUCim5erUARwdqYL+7dtAkyaaPb7aTlu0a9dOqToDU61nz2jJhcREoFkzYM+eghMpTEUkEmDlSiAwkCpE9ekDXL0qdlRaq379+rhw4QJiY2PRqFEjWFlZwdraGhs2bEBQUJDY4TEtIC/Qw4mFCAr/P9ZAj0V+fj5sbGz+PbSAJUuWYMKECRg0aBDs7e0xbdo0REREqD0OleMeC6blJBJA3lF47Zrmj6+2HosVK1agb9++6NGjR5ExjrNnz1bXYQ1CejottfDwIVC7Ns2x+Lf9ZqpmbAxs20aZ3OnTQPfuVGpBhfXX9YmXlxcOHz6M2NhYREZGwsjICM2aNSu1FCwzHFWr0mVysrhxGKTCFe40kFj4+flh5cqVGD58OJYuXYqYmBhMldfCBODg4KCRdW9UjnssmA6oXx84dQq4f1/zx1ZbYhEaGoorV67g1atXSl8qJBIJJxaVkJNDJ9AjI5WXXmBqZGEB7N9PK3zduEHjzs6cob5GVqzq1aujevXqYofBtIytLV2mpYkbh0HScGKxYsUKvP/++/jiiy9gamqK9evXK7UJx48fh4+Pj9rjUDnusWA6wNOTLh8+1Pyx1ZZYrF+/HhcvXkSzZs3UdQiDIwjAmDFUQszKinoq6tQROyoD4eAAHDpEq3zdugX06gUcOQJYWoodGWM6Q74GGScWIpCXmgU0MseiQYMGuHfvHm7evAl3d3c4OzsrPe7t7Y3vv/9e7XGoHPdYMB1QqxZdPnig+WOr7bSFp6cnZDKZup7eIM2dS6spGhsDe/cCLVqIHZGBkcmoi8jeHoiIAIYMUf5nzRgrlTyx4HUsRKDhHouQkBCYmZmhSZMmRZIKAPDw8ECTJk0QHR2N06dPqz0eleFys0wHuLnR5dOnmj+22lqX2bNnY9CgQTh8+DDOnTuntLHyW7sWkM9//eknGurPRNCwIQ2LMjMD9u0DJk9WnhTJirh48SKys7PFDoNpAe6xEJGGE4vw8HA0aNAAW7duRVxcnNJjOTk5OHPmDMaOHYuOHTvqVhl6eY9Fbi63/UxryYfIP3um+WOrrXXp1asX/vzzT7z//vto166dYvP39y/3c71+/Rrz5s1D/fr14ejoiK5duyI6Olr1QWupX3+l768A8NVXwMiR4sZj8Dp0ALZvp9IL69YBb6zbwJQtXboUzs7OGDVqFG7duiV2OExEPMdCRIV7VzWQWOzYsQPBwcE4e/Ys2rdvDxsbG8hkMlStWhUODg5YsGABmjRpgps3b1boe4FoCpdf5B5rpqXkiUVGBi1JoElqm2ORmZkJc3NzlZyJMDExQU5ODo4fPw5HR0eMGTMGgwcPxtmzZ1UQqXY7fRoYPJhOjIwfT8OhmBYYMIBqZn7+OS137uZGKxSyInbt2oVbt27h9OnT+OOPP+Dt7S12SEwk3GMhosI9FhqYYwEAzZo1w8aNGwFQL0VSUhIsLS1hb2+vkeOrReHEIje3oAeDMS1ia0t1Z169ol4L+WRuTVDpaYuQkBDF2hUWFhYlJhXlHVNpZGSEoKAguLm5wcLCAp9++inOnTuH1NRUlcStrW7epDnC2dlA7940HEqXeoz13uTJtIgeAIweTfMvDNDu3bsRX8rCBM+fP8edO3cwfvx4fPHFFxqMjGkb+QqwPDJOBBoeCvUmU1NTuLq66nZSASgnEjzPgmkpiQSoUoWua/qrskpbl/DwcDRs2FDtYyrNzMxgZGQEMzOzYh/Pzs5GWlqa0qZrYmOpqmlKCtC2LRAaqrGTTKw8goIKJnEHBgJ//y12RBo3Y8aMUj/PxsbGmDlzpgYjYtrK3JwuObEQQeHEgs9QVdybPRaMaSmximWoNLHQ1JjKsLAwtG7dGpYllPoMCgqCvb29YtO16lQpKTQ5+/FjWuTk99+5qqnWMjICNm8G3nsPePkS6NED+OcfsaPSqNTUVLi6upb4uLOzM5KSkjQYEdNWnFiISD4fwMhIo4nF5cuXdXMhvJJwjwXTEWLNaVP54EB1j6m8ffs21q5di3379pW4z+zZs5VW+ExLS9OZ5CI7G+jTB7h+HZBKaXSNk5PYUbFSmZkBYWE0qTsqirqaIiIAFxexI9MId3d3PH36tMTkIiEhAdbW1hqOimkjeSezPn3P1BnyHgsND4Nq3rw5jI2N4eXlBV9fXzRp0gS+vr7w9fWFiy62kRIJDR/Iy+PEgmk1veixeJOqx1Smp6cjMDAQY8aMQefOnUvcz9zcHHZ2dkqbLsjPB4YPB06epEzz0CGgZk2xo2JlYmdHfzAPD+D+feCDD6gHwwC0bt0amzdvLvHxkJAQNGnSRHMBMa1VuMeCK3VqmEiJxcGDB+Hg4IAGDRrAyMgIISEh6Nq1K6RSKaRSKbp3745Zs2YhNDRUo3FVSuGSs4xpKb3psVCXnJwc9OvXDzKZDMuWLRM7HJUTBGDqVGD3bhrC+dtvgK+v2FGxcnF1pS6mdu2AixeBgQPpD1l4TK4e+uqrr9C0aVNYWFhg8uTJMPn3n25+fj5+/PFHzJs3DydPnhQ5SqYN5IkFQCd7S5gmx9Sh8FAoDfr888+xadMm9OzZU3Hf//73P3z22WcYP348EhIScOXKFYSEhGDQoEEaja3CTE0pO+YeC6bFxOqx0InEQhAEjBw5EhkZGThy5Ijii4s+Wb4cWL2arm/dCnTqJG48rIK8vIA//qA/4MGDVCN40ya9niwplUpx4MABjBs3DvPnz0fdunVhYmKCu3fvwtraGlu3bkXTpk3FDpNpgcKJRHY2JxYaJe+x0HAVkJiYGDRs2FDpvq5du2Lp0qXYtWsXdu3apdF4VIJ7LJgO4B6LUoSFhSEkJARWVlZKcyXWrVunO2c4SrFjByCvwrlsGaAHv5Jha90a2LWLJsts3gy4uwOLFokdlVq1aNECf//9N86fP48bN24gOzsb3t7eaNWqFWxsbMQOj2kJXltMRCINhWrTpg02btyIb95YSLRp06YYMWKERmNRGfkbmXssmBaT9xBrek6bRhKLoUOHYt26dbCzs8PBgweRnZ2Nvn37lvnnAwMDFetj6JtjxwB52zplCjBtmqjhMFXp2RPYsAEYN46WS3d3Bz75ROyo1CY7OxurV6/Gvn37cPfuXVhYWMDLywuBgYEYP348jESom8+0T+GT5XyyV8NESizWrVuHVq1aITExEVOmTIG3tzdyc3OxZs0aODg4aDQWlZEnFvwmZlpMrGIZGkksrl69Cjs7O9y8eRPTp09H+/btcfLkSaxatUoTh9daUVFA37500uPDD6m3gumRsWOBuDhg4UJgwgSag9G7t9hRqVxGRgZatWqFvLw8jB49Gg0aNEBeXh5u3bqFVatWYevWrTh16hTMCw+wZwbJyIi2/Hz+TqZxIs2xqF+/Pi5cuIDPP/8cjRo1gpmZGfLy8mBqaopNmzZpNBaVkQ+F4h4LpsX0OrEwNTWFIAjYsmUL5s6diyFDhqBZs2aaOLTWio6mtSrS04GOHWleBZ/U1UPz51NysWkT8NFH1EXVtq3YUanUzJkz4eXlhd27dyvNf+rZsyemTp2KgQMHYv78+UWGQjDDZGJC/+g4sdAwkeZYAICXlxcOHz6M2NhYREZGwsjICM2aNYNUKtV4LCrBQ6GYDtCLxCIjI6PY8dSffPIJWrRogRcvXmD+/PkAgJcGUoqzOElJtNTB06dAo0ZUOIhP5uopiQRYvx6IjwfCw2mIVEQErXyoJ44dO4ZDhw4VW1TBxMQEK1euxPvvv8+JBQPAiYVoRBoKVVj16tVRvXp10Y6vMjx5m+kAsRILlbYwDg4OSExMLHL/uHHjcPToUVy9ehXW1ta4f/8+WrVqpcpD64zMTPpueecOUKMGLX2gomU+mLYyMaHJ3C1bAi9eUFb55InYUalMamoq3NzcSnzc1dUVz58/12BETJvxdzKRyF9wEXos9I78TOCrV+LGwVgp5ImFpjvWVJpY5OXlIV9+VgRAu3bt8OzZMwBAlSpVFL0ZderUwdatW1V5aJ2Qm0sVn86dAxwcaMkDd3exo2IaYW1NZWjr1gViYoD339d8DTg1adq0KcLDw0t8PDw8vEi5SWa4OLEQifzbhZ6vq6MRVlZ0mZUlbhyMlUIveizedPXqVYMe8lSYIACffQb8/jtgYUGX3t5iR8U0ytmZsslq1YArV4B+/TT/iVeDhQsXYty4cdixYwfyCtUQzc/PR0hICEaPHo3Zs2eLGCHTJvIT5pxYaJg8seDFQyrP2pou+fsN02Lycwh6lViwAosXAz/+SEPud+wA/P3FjoiJolYtWjjPxoYmco8cWTD2WUe1aNECW7ZswX//+184OzujZcuWaN26NZydnTF9+nR8//336NKli9hhMi3BPRYi4R4L1ZH3WGRmihsHY6XQi8nbALBjxw688847aNSoEQBAoscrDpfV5s3AvHl0fe1aKjHLDJifHxAWBvToQVmmmxvw3XdiR1UpH3zwAbp06YKIiAjcuXMHr1+/Rr169eDv788L5DElnFiIRP7tghOLypP3WHBiwbSYXiQW/v7+mD9/PtLT02Fqaorc3FzMnTsX7du3h5+fH3x9fWFhYaHKQ2q98HBaIw0A5syh5QwYQ5culHEOH04LmLi7A//5j9hRVYqZmRkCAgLQtm1bpKamIjU1Fbdv31Zc79evn9ghMi3AiYVIuMdCdeQ9FjwUimkxvUgsTp06BQC4d+8eLl26hMjISFy6dAlz585FSkoKTExMUL9+fVy9elWVh9VaFy4AAwfSukQff0zDoRhTGDaMqkPNmkXLrkultFKijmnevLkieUhNTUXuv98YBUGAtbU17OzsUKVKFU4sGABOLETDiYXqcI8F0wF6kVjI1a1bF3Xr1sVHH32kuO/hw4f4+++/cfnyZXUcUuvcu0cjXTIzga5dgZ9+ovkVjCmZMQOIjaUxcsOHAy4uQECA2FGVS0xMDKysrDBs2DD069cPHh4eioTCiFd9ZG/gxEIkPHlbdXiOBdMBelFutjSenp4YMGAAlixZoqlDiubZM0omEhOBZs2AvXv5JBErgUQCrFoF9O9PpxX69AGuXRM7qnJ5+vQpfvzxRyQlJaFnz54YPHgwtmzZgqdPn4odGtNCnFiIhHssVIeHQjEdoJflZg1RejotUfDwIVC7Ns2x4LmrrFTGxsD27UD79rS2RbduwKNHYkdVZsbGxujatSu2bNmC+Ph4TJkyBVFRUWjWrBnat2+PtWvX4omaFgTMyMjA4MGD4eLignbt2uHWrVtqOQ5THXliUagyMdME+dl1S0tx41AzjbQJTk50mZCg+udmTEW43KweyMkBBgwAIiOVlyxg7K0sLID9+wEfH5p30a0brdKtY8zNzdGnTx9s2bIF169fR4sWLTB16lTUqFFDLccbMWIETE1NER0djV69eqFbt26KOR5MO3GPhUhSUujSwUHUMNRNI22ChwddRker9nkZUyHusdBxggCMGQP873/US/rHH0CdOmJHxXRK4eXYb90CevfWuZVd79+/j2+++Qb+/v5wc3PDhQsXsGTJEty5c0flx4qJicFvv/2GJUuWwMrKCjNmzMDLly9x+PBhlR+LqQ4nFiJJTqbLKlVEDUOdNNYmyBOLe/d0fh0ipr/0avK2IZo7FwgOplEte/YALVuKHRHTSTIZJRf+/sCZM8CQIfSGki9XrIXu3LmDvXv3Yu/evXjw4AECAgIwbNgw7N27F66urmo77vnz5+Hh4QF3d3cAtGZOixYtEBUVhQ8++EBp3+zsbGRnZytup6WlqS0uVjpOLDREEIC7d4ETJ4A//wR27aL71dR7qA001ibUqUMJWkoKDS3z9qZkw8MD8PQsuO7hAdjbV+6XYqyCOLHQYWvXAkFBdP2nn2iOBWMV1rAhDYvq0gXYtw+YPJneZFpaVszb2xtSqRS9e/fG3Llz4eTkBFtbW6Snp8PIyAj29vYwNzdX+XHj4uLg4uKidJ+TkxPi4+OL7BsUFISFCxeW+xihoaEIDQ0t9rFRo0ahT58+AIDU1FQMGzas2P2kUik2btyouP3tt98iIiKi2H0XLFgAPz8/AMC1a9cwd+7cYvdr2bIlvvzyS8Xtzz77DI8fPy52382bN6Nq1aoAgPDwcKVYCgsMDMTw4cMVt3v16lXsfpaWltgl/5IKYMOGDTh48GCx+06dOhUdO3YEADx69AgTJ07EzZv0WFAQsHUrXffy8sJ3hRaJnDNnDq5fv17sc65cuRK1a9cGAJw5cwZLly4tdr/OnTtj0qRJittDhw4t8ctjWFgYTP8dkKzTf/OMDCApCZvr1EHVc+eA+HiEA1BEYm4O/O9/GOTujkGDBhV7HF2miTYBAL2Oy5cDn35K39quXKGtOFWqFJ9wyDc7u4rFwNhbiJVYSARBEDR7SM1LS0uDvb09UlNTYafiD/GvvwKBgXRy6KuvgELtPmOVs2cPrWshCMCSJcDs2So/hCo+GxMnTkRqaiqSk5ORkpKiuExJSUHWv0O5TE1Nlc4OqsLy5cuxb98+nDlzRnHfkCFD4OzsjFWrVintW9zZSZlMppY2gZWuSxfgyBGqVzBkiNjR6LgbN4CTJ4GICOqZeLNIgrk50KYN0LEjlbFu1Yruewt1/s9UJ423CZmZNM+iuO3hQyoN+TaOjkWTjerVqVekXr0y/b0YK05MDL2dLCwqP6q6PG0C91hUwunTwODB9L3vk09oOBRjKjNgABAfD3z+OS3b7uZGKy1qmbVr15b42OvXrxVJhqq5ubkh4Y2qLImJiWjcuHGRfc3NzdXSa8LKj4dCVVBuLpWijogAzp6lLSZGeR8zM6B1a0oiOnak6xYWooQrBo23CVZWVHDDx6f4xzMy6G/0ZsIhv56UREU6Xrygqi/FcXGhIbIyGQ1jk1+Xb1JpwYeKsULEWseC340VdPMm0KsXkJ1Nc2x/+EFrR6owXTZ5Mi2g9913VB2gWjWqGKUjzMzM4OLiUmR4gio0a9YM9+/fR0JCApydnZGfn4+LFy9ithp6dpjqcGJRRqmpwOXL1CNx8CBw/XrxC7LVqEGLawYEUO+EnpeTLY3WtQk2NkCDBrQVJz29IPGQJxwPHwJ37tBlVhbw/Dltly4V/xzGxtTr4eoKNG5M7wdHR6B+fSoE4upKyYkWz9Nj6iEvN5uXR5um3gKcWFRAbCx9t0tJoXZ8xw7+zDI1+uYbGuIQEkLj7k6epJUXDVy9evUQEBCAuXPnYvXq1VixYgWkUik6dOggdmisFJxYFEMQgKgoGs508iStY3PlStGKQ3Z29E+nXTugbVuqEmJrK0bEWknn2gRbW5pT17Bh0ccEgXo0Hj8u2B49Ur4dF0cfpIQE2kpaXNXWluZ4uLhQLXwvL+rpcHICGjWioVeWlnx2VM/IeywA6rXgxEJLpaQA3bvTZ7p+feDAgYJFOBlTCyMjYPNmWtL96FGqDnD2LK3AaOB27tyJ4cOHQyaToWnTpjh8+DAk/M9Rq8n/uRl0YpGWRonDmTPA8ePUG1HcSvVOTrRw5jvv0Na0KbUHrER60yZIJEDVqrQ1bVr8Pnl59L55/pySimfPqAfk/n26Ln8sPR24erX045mZUfIhk9Gw25o1qfejRg3q+ahWjSpc8crtOqNwYvH6teZGRXJiUQ7Z2UDfvvQ/QCqlqqDyBTgZUyszMyAsDOjQgc5sdutGyYWzs9iRiapq1aolViRi2sngeizy84Hbt4ELF2h+xMGDRSdZA/RFsmbNgrkRAQE0eZeVi0G1CcbG9KXf3b3k5CM3l5LYpCR63926Rb0biYl0PTYWePWKvnneuUNbSSQSek/KZJTwNGlCCUfVqpSAVK9OyQcP4dAKhXNATVaG4sSijPLzaRjriRPUq3jwIP0PYExj7Ozojde2LZ2R6tGD6tNbW4sdGWNlpveJxZUrNBE3IoLmSDx4ULDqdWFSacGQprp1KaGwsdF0tEzfmZiUPnRWEGiSeWwsDcWIjaVejkePCrbYWFpgURCUk4+dO4s+n7k59Xw4OtKZ19q1aXN1pS9PjRvTe5+pnbExbXl5nFhoHUEApk4Fdu+mDHDfPkrUGdM4eVdZu3bAxYvAwIHAb79x9zTTGXqVWLx+TZNqDx2iHsTbt2nce3H8/Oissr8/zZPw8tJsrIwVRyKhL/ze3rSVJCenoLcjLg745x9aeTw5mW7fv0+TzbOz6XNQGktLSjwcHam3o25dGn4lk9EYc6mUNgOqaKYuZmb0Z+HEQsssXw6sXk3Xt24F3n1X3HiYgfPyAv74A+jUiXowxo8HNm3iiXdMJ+hsYpGdTePYL12iHonISBq3XtJ/7L596bPauTPNk+Dkn+kyU1PqbSjN69fU63H/PvXSycvqJiTQnI87d+h6VhYlI/Ik/M8/i38+B4eCJEO+ubkVvY97+kokTyw0WXKWE4u32LED+OILur5sGaCHi5UyXdS6NXVD9+1LE7vd3YFFi8SOirG30onEIiurYEiTPJG4fr34oKtUoTO+vXtTYYWWLelMLCf6zNCYmRUMfSpJejrN73jxgpKNe/do/sedOzTs6ulTWr8pO5t6Q5KTqb5/aWxt3558uLnRcGID+1zKz2dwj4WWOHYMGDGCrv/nPzQcijGt0asXsH49rc741VeUXHzyidhRMVYqrUss0tNptdO//6YvNxcu0BoCeXlF93VyovHqfn4Fl56eBvdlhbEKs7UtKH9bEkGgHo/4eJpwHh+vvBW+7+VL+gynpwN375Z+bEtLmmDu7k6JhpMTTT6vVatgLZBatagHRE8+0/LKUJxYaIGoKDoZnJNDw9iXL9eb9xnTJ+PGUXfyokXAhAnUMPbuLXZUjJVI1MQiLY0mVF+6VLCVVAXHxYWSh8KJhEzG/wgYUzeJhIZBOTiUvKq5XHp62RKQ1FTqibx3j7bSmJgUlNgt3APy5qWzs9aXf+bEohQZGRkYN24cjh49irp162LTpk3wLm2iUSVER9NaFenpVKgjOFjr3zvMkC1YQMnFzz8DH31EdfHbtBE7KsaKpZHEQhBoSMWBA7T2y5MnBQuMlaRdO0oeWrWiNSPc3TmJYEzbyXtA3laaOTOT2oQHD6g9ePasIPl4+rTgMj2dGqeYGNpKY2JCyYerK/V41KxJsbi4UI+MVFqw+rlI1RstLeny5UvNHVNnEosRI0bA2toa0dHRWLNmDbp164Z//vkHJiaq/RWSkmiJgKdPaUHK336j6mmMaS2JBNiwgd604eHABx9Qqcv69cWOjLEiVJ5YvHpFZ4OiomguxA8/0JeIktSoUdAT4eMDNGjA60Uwpu+srGiYU61ape+XkUFJRmwsJRnyXo83L58/p0as8CT0t5H3gjo60nofTk7UHjk6Ulx2dpSoqDAJqVaNink9e6ayp3wrnUgsYmJi8NtvvyEmJgZWVlaYMWMGvvvuOxw+fBgffPCByo6TmQn07Ek94zIZVRC0t1fZ0zOmPiYmwK5dVCnqwgXKjs+d43rhTOtUKrHIyqIyy8HB9M89KQk4f77k/WvXpn/W/v5U0aZ3b4NfVJIxVgobGzrR8LaTDbm5Bb0eMTH0xTEvr2BCeno6JR9PnlC7BdDtQ4feHoODA+DhQaV4nZyoDXNyorK8Tk500tDFhb6gvqVX1dWVLuPj335YVdGJxOL8+fPw8PCAu7s7AEAikaBFixaIiooqNrHIzs5Gdna24nZaWtpbj5GbSxWfzp2jv+nhw9SDxZjOsLamMrTt2lHD1r07cOoUNUqMaQl5YlHqmN+cHBqyEBxM1Zmio6lKzIMHtFppcVq3BurUAXx9qRfj889pWAJjjKmafB6GuzvQvHnJ+wkCze2Kjqa27PVrmu+RlEQJSUoKJSUpKTRe6dWrgmpYb2NqSidKnJ2pa8LUlIZjSaXU9tWpg+ZGUpyAFA+u2wOwUs3v/hY6kVjExcXBxcVF6T4nJyfEl5CCBQUFYeHChWV+fkEAPvsM+P33UBgZhcLHB5g1S3mfUaNGoU+fPgCA1NRUDBs2rNjnkkql2Lhxo+L2t99+i4iIiGL3XbBgAfz8/AAA165dw9y5c4vdr2XLlvjyyy8Vtz/77DM8fvy42H03b96MqlWrAgDCw8OVYiksMDAQw4cPV9zu1atXsftZWlpi165ditsbNmzAwYMHi9136tSp6NixIwDg0aNHmDhxYrH7eXl54bvvvlPcnjNnDq5fv17svitXrkTtf0vXnTlzBkuXLi12v86dO2PSpEmK20OHDi0xoQwLC4PpvzXYQkNDERoaWux+Ovk3v3ABG2vUoAbryhUMatsWgyIjC2ZwMSayf9+q+PFH4MgR6lSYYbIC3pmXUP3UDjoDJwglP4GjI60P4elJw/6cnKinTsXDYhljrNIkEupZ8PWlrTT5+dTj8egRleNNSKAE5MEDup2UVLBCeloanYB58oS2Ekz7d8MvQM4vJog1rolUS1cYWZghw0GGLGltmNuYwMTGEq12TlHJtDKJIJTWgmuH5cuXY9++fThz5oziviFDhsDZ2RmrVq0qsn9xPRYymQypqamwK+bs7cqVVEpWIgH27gX69VPLr8GY5kRGAh060HjRxYuBEhKYtLQ02Nvbl/jZ0FeG+ntrg6QkemveuEG3HZGEh/CEHdKVdzQ1LaiH36YNFSbw9aUhADypWq0M8fNhiL8z02GvXlHikZBAQ6yePqUy2fLKV+npQHIyhCfxyH/6DMYooaf3Xxmwhs3r5BIX8izP50MnTvG4ubkhISFB6b7ExEQ0LmEVSHNzc5iXY8Z1//60cPGECZxUMD3h5weEhQG//AJMmyZ2NIwpODnRgtWXLtH/v6P7LXDw9mw0fXkGXu4ZgLc3MGoUTXI0NhY7XMYY0z4WFjQZWCYrdTcJAOO8PGQ+TUP8ybt4mZCJvLv/IC8mFq/SXsMoKQFmqc8Rb+uFnpmZKplYrBM9Fnfv3oW3tzeePn0KZ2dn5Ofno2rVqvj1118VQ29KU5ZM69Ur+jsxZkgM9Sydof7ejJWFIX4+DPF3ZqysyvP50InVGerVq4eAgADMnTsXWVlZCAoKglQqRYcOHVR2DE4qGGOMMcYYqzidSCwAYOfOnYiNjYVMJsOJEydw+PBhSHicLWOMMcYYY1pBJ+ZYAEDVqlVLrEb0NvLRXmUpO8uYIZF/JnRgRKRKcZvAWMkMsV3gNoGxkpWnTdCZxKIy0tOp2ojsLZNcGDNU6enpsDeg1SC5TWDs7QypXeA2gbG3K0uboBOTtysrPz8fT548ga2tbYnDp+QlaR8/fswTt9SMX2vNKMvrLAgC0tPT4ebmBiMjnRkZWWncJmgXfq01h9uF4pWlTXgTv29Lxq9N6XTt9SlPm2AQPRZGRkaoXr16mfa1s7PTiT+yPuDXWjPe9jobyhnJwrhN0E78WmsOtwvKytMmvInftyXj16Z0uvT6lLVNMIxTEYwxxhhjjDG14sSCMcYYY4wxVmmcWPzL3Nwc8+fPL9eK3axi+LXWDH6dK4dfP83h11pz+LVWHX4tS8avTen0+fUxiMnbjDHGGGOMMfXiHgvGGGOMMcZYpXFiwRhjjDHGGKs0TiwYY4wxxhhjlcaJBWOMMcYYY6zSOLEAkJGRgcGDB8PFxQXt2rXDrVu3xA5Jb+3fvx/m5uZwcHBQbGlpaWKHpRdycnLw3XffwdbWFl9//bXi/tzcXEyaNAlSqRR+fn44ffq0iFHqBm4TNIfbBPXhNkF9uI0o3uvXrzFv3jzUr18fjo6O6Nq1K6Kjo8UOS+s8fvwYVlZWWLhwodihqBwnFgBGjBgBU1NTREdHo1evXujWrRtyc3PFDksvPX78GBMmTEBycrJi05VVJ7XdqlWr8Oeff6JGjRpK98+aNQu3bt3C3bt3MX36dPTo0QNJSUkiRakbuE3QHG4T1IfbBPXhNqJ4JiYmyMnJwfHjx/HkyRM4Oztj8ODBYoeldebMmaOXpWYBAIKBi46OFoyNjYXY2FhBEAQhPz9fcHJyEg4cOCByZPppxowZwrJly8QOQy+9fPlSEARB6Nixo7B48WJBEAQhIyNDsLOzEyIiIhT7+fn5CWvWrBElRl3AbYJmcZugPtwmqAe3EWV35swZAYCQkpIidiha4+LFi0KNGjWEwYMHCwsWLBA7HJUz+B6L8+fPw8PDA+7u7gAAiUSCFi1aICoqStzA9NTjx49RvXp1scPQS1ZWVkXuu3HjBl69eoXWrVsr7mvZsiW/v0vBbYJmcZugPtwmqAe3EWVnZmYGIyMjmJmZiR2K1pg2bRpmz56tt6+JwScWcXFxcHFxUbrPyckJ8fHxIkWk3x4/foxx48bB1dUVgwYNwvPnz8UOSa/FxcXByckJRkYFH3V+f5eO2wTN4jZBs7hNqDxuI8ouLCwMrVu3hqWlpdihaIV9+/YhKSkJY8eOFTsUtTH4xCI/P1+pgQUAQRBgamoqUkT67ejRo0hOTsapU6eQkpKCbt26iR2SXuP3d/nxa6ZZ3CZoFr+/K49fw7K5ffs21q5diwULFogdilbIycnBzJkzsWLFChgbG4sdjtoYfGLh5uaGhIQEpfsSExMhlUpFiki/mZubw8jICPXq1cNPP/2Ey5cv4+HDh2KHpbfc3NyQlJQEQRAU9/H7u3TcJmgWtwmaxW1C5XEb8Xbp6ekIDAzEmDFj0LlzZ7HD0Qo//PAD6tWrhy5duogdiloZfGLRrFkz3L9/X9FI5Ofn4+LFi2jVqpXIkek/ecaur+MMtYGPjw8kEgkiIyMV9/3111/8/i4Ftwni4TZB/bhNqDxuI0qXk5ODfv36QSaTYdmyZWKHozV27NiBI0eOwNLSEpaWlti2bRsWL16Md999V+zQVMrgE4t69eohICAAc+fORVZWFoKCgiCVStGhQwexQ9M7qampOH/+PPLy8vDixQtMnDgR3bt3V0yAY6pnb2+PwYMHY968ecjIyMC2bdvw5MkTDBgwQOzQtBa3CZrDbYLmcZtQedxGlEwQBIwcORIZGRnYs2cPTExMxA5Ja1y4cAHZ2dnIyspCVlYWhg0bhi+//BLHjh0TOzSVMvjEAgB27tyJ2NhYyGQynDhxAocPH4ZEIhE7LL0THx+PyZMnw9nZGY0bN0a1atUQGhoqdlh6b82aNbC3t4enpyfWr1+PY8eOwdraWuywtBq3CZrBbYI4uE2oPG4jihcWFoaQkBBcvXoVMplMseglf64Nh0QoPNCSMcYYY4wxxiqAeywYY4wxxhhjlcaJBWOMMcYYY6zSOLFgjDHGGGOMVRonFowxxhhjjLFK48SCMcYYY4wxVmmcWDDGGGOMMcYqjRMLxhhjjDHGWKVxYsEYY4wxxhirNE4smGJlTHt7e0gkEtjb28PBwQG+vr6Ij49HkyZNEB8fr9YYsrKyUL9+fTx48ABbt25F3bp1K/V8wcHBGDx4sIqiY8ywcJvAGGOsIjixYEhOTkZycjKuXr0KALh58yaSk5Nx5coVSKVSREVFQSqVqjWGdevWoUOHDqhVq5ZKnm/IkCGIjIzE33//rZLnY8yQcJvAGCuLd955BxKJRGkzNjZGWlpamX5+0KBB2LZtW4WPv3v3bowePbpcP7No0SLMmTOnwsdkpePEgmmFjRs3YtasWSp7PmNjY8yYMQMbNmxQ2XMyxjSH2wTGtJsgCIiKisKyZcsQHx+vtNnZ2b315y9duoSoqKhK9ST+/PPPGDNmTLl+5j//+Q9++eUXPHnypMLHZSXjxIKVKi4uDhKJBDExMQCAatWqITQ0FJ06dYKLiwsmTJiA58+fY+DAgXB2dkbPnj2RmZmp+PmjR4+iefPmqFq1Krp3717sBzkmJgZ5eXnw9PQs8lhCQgI8PDywZMkSTJ06FR06dFB6vGbNmti8eTO2bt2Kfv36Yfbs2ahRowYA4L333sPRo0dV+XIwZvC4TWCMAcC9e/eQnp6Od955B66urorNxcWlTD+/ceNGzJgxA8bGxhU6/qNHjxAXF4c2bdqU6+fs7OwwduxYbN68uULHZaXjxIKVi7m5ORYvXoyffvoJFy5cQHBwMFq1aoUZM2bg9u3buHbtGr7//nsAwK1bt9CnTx8sX74cT548gY+PDwYMGFDkOe/evQsvL68i9+fk5KB///7o2bMn5syZg7Fjx+L06dN4+PAhAOD8+fNITExUPOeRI0eQk5ODO3fuAABq1KiB58+fIycnR10vB2MGj9sExgzTpUuXYGJigsaNG1fo548cOYL33nsPAM2pkkgkOHDgANq2bQtXV1csWrQI9+7dQ8+ePeHo6IiRI0ciPz9f8fO//PILRo0aBQDo168fPv74Y8VjL1++hJWVFY4fPw4AkEgk+P333+Hj44O9e/fySQY14sSClYuxsTEmTJiA2rVrw8PDA/7+/mjTpg2aN28OJycnvPvuu7h48SIAYP369ejWrRs6dOgAMzMzTJw4EWfPnsWjR4+UnjMxMREODg5FjjVx4kS4uLhg9erVAABvb2+0bdsW27dvBwDs2bMHgYGBsLW1VfzMokWLYGlpqbjt4OCAxMRElb8OjDHCbQJjhikyMhJ5eXlwcnKCjY0NbGxsytx78Pr1azx79gwymQwAnaAAgBUrViAsLAwHDhzA/PnzERgYiFWrViEyMhJ79+7F7t27AQD5+fkICQnBsGHDAABjx45FWFgYXr58CQD4448/4OLigoCAAMUxv/nmGxw4cACBgYHw8vLC/fv3VfZasAKcWLByq1KliuK6VCqFtbW14ranpydevHgBAPjnn38QHh6uqDDTtGlTWFhY4Pnz50rPZ2lpiVevXindFx0djQMHDuDZs2cQBEFx/9ixYxEcHAxBELB3716MHDlS8Vi1atVgZWWl9DyZmZlK8THGVI/bBMYMz6VLlzBw4EBERUUptl9//bVMP5uUlKR08sDIiL6OTp8+HVKpFC1atICXlxd69OihOGnRsmVLxUmKo0ePwtfXF87OzgCArl27wsnJSXH8PXv24OOPP4ZEIlEcY/z48ahduzYAwNHREQkJCZV/EVgRnFgwlSr8IXZ3d8eoUaMUFWZSUlKQlZWF5s2bK/2Mq6trkdKVzs7OuH79OuLi4vD1118r7h8wYAASEhKwadMmGBkZKY2vLnxsgLpWX79+XaZJZIwx9eA2gTH9dPnyZbRv3x516tRRbPJqcRs3boSfnx8aNmxY7OTs4k4eAGU/SfHmpG0jIyOMGjUKwcHByMzMxOHDh5WGRgFQqjDHJxjUhxMLpjYjR47E9u3bcezYMQiCgLy8PMTFxRXZr2nTprh9+zays7MV91lbW8PR0RHBwcH4+uuvcfbsWQCAlZUVBg8ejFmzZmHUqFFFvjgUdv78+XJP6mKMqQ+3CYzphwcPHiAlJQV+fn5FHktOTsaPP/6Iixcv4vr161i3bl2RfapUqYKsrCylz/jbyD/bSUlJiIyMROfOnZUeHzVqFE6ePIng4GC0aNGiSKnqwm1DfHw8XF1dy3xsVnacWDC1adOmDbZt24YvvvgCjo6OqF69OpYvX15kP3Nzc3Tq1Al79uwp8pi/vz+mTJmCIUOGIDU1FQDQv39/pKenKyZtlSQ4OBh9+/ZVzS/DGKs0bhMY0w+XLl2CsbExfH19izxmYmKCpKQkzJw5Ezdu3FDqhSisdevWOH/+fLmPvX37dgwZMkQxfEpOJpPhvffew9y5c9+6tsVff/2Ftm3blvvYrAwExrTAxYsXBS8vLyE3N/et+3788cfC6NGjS93n/v37gkwmE16+fKmqEBljGsRtAmPaa9asWYKPj0+Jj6elpQnBwcFC/fr1hQMHDhS7z5o1a5Q+twCEM2fOKG537NhRWLx4seL2mDFjhBEjRgi+vr5CdHR0sc+5adMmwcnJScjKylK6/83nDggIEP7444/Sf0lWIRJBKDQLjjERBQcHo0OHDqhZs2axjycnJ2POnDnYv38/oqKiSq2VffLkSUgkErzzzjvqCpcxpmbcJjCme+7du4e6desCAMaNG4euXbuif//+RfZ7+fIlvL29ceLEiSLDlipCEAS8++676NSpE7788ssS9zt16hQmT56My5cvlzp0klUMJxZMZ1y6dEmxcq680WKMGS5uExjTPsOHD8f58+dhZWWF9u3bY9WqVUWGLcmdOnUKgiAUWeiyvOLi4jB58mTcuXMHFy9eVCox/ab9+/fD09OzwutvsNJxYsEYY4wxxnTWvn37sH37dqxbtw7VqlUTOxyDxokFY4wxxhhjrNK4KhRjjDHGGGOs0jixYIwxxhhjjFUaJxaMMcYYY4yxSuPEgjHGGGOMMVZpnFgwxhhjjDHGKo0TC8YYY4wxxlilcWLBGGOMMcYYqzROLBhjjDHGGGOVxokFY4wxxhhjrNI4sWCMMcYYY4xV2v8Db3Uo1fl6n04AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "########################################################################\n",
    "## plots:\n",
    "########################################################################\n",
    "print(\"plotting...\")\n",
    "plt.figure(1,figsize=(8,4))\n",
    "\n",
    "Fs_to_m_per_year=S0*area/year\n",
    "plt.subplot(2,3,1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,DeltaS_steady[0,:],'r.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,DeltaS_steady[1,:],'g.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,DeltaS_steady[2,:],'b.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,0*Fs_range,'k--',dashes=(10, 5),linewidth=0.5)\n",
    "plt.title('(a) Steady states',loc=\"left\")\n",
    "plt.xlabel('$F_s$ (m/year)');\n",
    "plt.ylabel('$\\Delta S$ (ppt)');\n",
    "plt.xlim([min(Fs_range/Fs_to_m_per_year), max(Fs_range/Fs_to_m_per_year)])\n",
    "\n",
    "plt.subplot(2,3,2)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,q_steady[0,:]/Sv,'r.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,q_steady[1,:]/Sv,'g.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,q_steady[2,:]/Sv,'b.',markersize=1)\n",
    "plt.plot(Fs_range/Fs_to_m_per_year,0*Fs_range,'k--',dashes=(10, 5),linewidth=0.5)\n",
    "plt.title('(b) Steady states',loc=\"left\")\n",
    "plt.xlabel('$F_s$ (m/year)');\n",
    "plt.ylabel('$q$ (Sv)');\n",
    "plt.xlim([min(Fs_range/Fs_to_m_per_year), max(Fs_range/Fs_to_m_per_year)])\n",
    "\n",
    "plt.subplot(2,3,3)\n",
    "plt.plot(DeltaS_range,rhs,'k-',lw=2)\n",
    "plt.plot(DeltaS_range,rhs*0,'k--',dashes=(10, 5),lw=0.5)\n",
    "# superimpose color markers of the 3 solutions\n",
    "Fs=Fs_func(0.0,0.0,False)\n",
    "yy=steady_states(Fs,alpha*DeltaT)/beta\n",
    "plt.plot(yy[0],0,'ro',markersize=10)\n",
    "plt.plot(yy[1],0,'go',markersize=10)\n",
    "plt.plot(yy[2],0,'bo',markersize=10,fillstyle='none')\n",
    "plt.title('(c) Stability',loc=\"left\")\n",
    "plt.xlabel('$\\Delta S$ (ppt)');\n",
    "plt.ylabel('$d\\Delta S/dt$');\n",
    "\n",
    "plt.subplot(2,3,4)\n",
    "plt.plot(T[:N2]/year/1000,FWplot[:N2]/Fs_to_m_per_year,'b-',markersize=1)\n",
    "plt.plot(T[N2+1:]/year/1000,FWplot[N2+1:]/Fs_to_m_per_year,'r-',markersize=1)\n",
    "plt.plot(T/year/1000,0*T,'k--',dashes=(10, 5),linewidth=0.5)\n",
    "plt.xlabel('Time (kyr)');\n",
    "plt.ylabel('$F_s$ (m/yr)');\n",
    "plt.title('(d) $F_s$ for hysteresis run',loc=\"left\");\n",
    "\n",
    "plt.subplot(2,3,5)\n",
    "plt.plot(T[:N2]/year/1000,qplot[:N2]/Sv,'b-',markersize=1)\n",
    "plt.plot(T[N2+1:]/year/1000,qplot[N2+1:]/Sv,'r-',markersize=1)\n",
    "plt.plot(T/year/1000,T*0,'k--',dashes=(10, 5),lw=0.5)\n",
    "plt.xlabel('Time (kyr)');\n",
    "plt.ylabel('MOC, $q$ (Sv)');\n",
    "plt.title('(e) MOC transport, $q$',loc=\"left\");\n",
    "\n",
    "plt.subplot(2,3,6)\n",
    "plt.plot(FWplot[:N2]/Fs_to_m_per_year,qplot[:N2]/Sv,'b-',markersize=1)\n",
    "plt.plot(FWplot[N2+1:]/Fs_to_m_per_year,qplot[N2+1:]/Sv,'r-',markersize=1)\n",
    "plt.title('(f) $q$ hysteresis',loc=\"left\");\n",
    "plt.xlabel('$F_s$ (m/yr)');\n",
    "plt.ylabel('$q$ (Sv)');\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"stommel-model.pdf\", dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a8c60001",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation, PillowWriter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6c991f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7,5))\n",
    "\n",
    "Fs_axis = Fs_range / Fs_to_m_per_year\n",
    "FW_anim = FWplot / Fs_to_m_per_year\n",
    "q_axis = q_steady / Sv\n",
    "q_anim = qplot / Sv\n",
    "\n",
    "# choose stable vs unstable styling\n",
    "ax.plot(Fs_axis, q_axis[0,:], 'r-', lw=2, label='Stable branch')\n",
    "ax.plot(Fs_axis, q_axis[1,:], 'r-', lw=2)\n",
    "ax.plot(Fs_axis, q_axis[2,:], 'k--', lw=1.5, label='Unstable branch')\n",
    "\n",
    "traj_line, = ax.plot([], [], 'b-', lw=2, label='System trajectory')\n",
    "dot, = ax.plot([], [], 'bo', markersize=10)\n",
    "\n",
    "ax.set_xlim(np.min(Fs_axis), np.max(Fs_axis))\n",
    "ax.set_ylim(np.nanmin(q_axis)-2, np.nanmax(q_axis)+2)\n",
    "ax.set_xlabel(r'$F_s$ (m/yr)')\n",
    "ax.set_ylabel(r'$q$ (Sv)')\n",
    "ax.set_title('Stommel model tipping point and hysteresis')\n",
    "ax.legend(loc='best')\n",
    "\n",
    "def init():\n",
    "    traj_line.set_data([], [])\n",
    "    dot.set_data([], [])\n",
    "    return traj_line, dot\n",
    "\n",
    "def update(frame):\n",
    "    traj_line.set_data(FW_anim[:frame], q_anim[:frame])\n",
    "    dot.set_data([FW_anim[frame]], [q_anim[frame]])\n",
    "    return traj_line, dot\n",
    "\n",
    "anim = FuncAnimation(\n",
    "    fig, update, frames=len(FW_anim),\n",
    "    init_func=init, interval=20, blit=True\n",
    ")\n",
    "\n",
    "anim.save(\"stommel_tipping.gif\", writer=PillowWriter(fps=25), dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "601db879",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.13.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}