main_case6.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import flopy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "# ! Custom modules\n",
    "from utils import pht3d_fsp\n",
    "from utils.tools import calculate_MAS_sums, clear_directory\n",
    "from utils.pht3d_utils import copy_pht3d_datab, run_pht3d_program\n",
    "\n",
    "clear_directory()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_folder = 'case6'\n",
    "model_name = 'case6_modflow'\n",
    "model_ws = './models_folder/' + data_folder\n",
    "\n",
    "Lx = 0.2\n",
    "Ly = 1.0\n",
    "\n",
    "nrow = 1\n",
    "ncol = 180\n",
    "nlay = 1\n",
    "\n",
    "delr = Lx / ncol\n",
    "delc = Ly / nrow\n",
    "\n",
    "top = 1.0\n",
    "botm = 0.0\n",
    "\n",
    "laytyp = 1\n",
    "\n",
    "hk = 1.0\n",
    "\n",
    "nper = 2\n",
    "perlen = [0.05944444, 0.1488889]\n",
    "nstp = [250, 600]\n",
    "steady = [True, True]\n",
    "\n",
    "model = flopy.modflow.Modflow(modelname=model_name, model_ws=model_ws, exe_name='./bin/mf2005')\n",
    "dis = flopy.modflow.ModflowDis(model, nlay=nlay, nrow=nrow, ncol=ncol, \n",
    "                               delr=delr, delc=delc, top=top, botm=botm,\n",
    "                               nper=nper, perlen=perlen, nstp=nstp, steady=steady)\n",
    "\n",
    "ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)\n",
    "ibound[0, 0, -1] = -1\n",
    "strt = np.ones((nlay, nrow, ncol), dtype=np.float32)\n",
    "bas = flopy.modflow.ModflowBas(model, ibound=ibound, strt=strt)\n",
    "\n",
    "lpf = flopy.modflow.ModflowLpf(model, hk=hk, laytyp=1, ipakcb=53, ss=0.1)\n",
    "\n",
    "wel_data = {0: [[0, 0, 0, 2.839]],}\n",
    "wel = flopy.modflow.ModflowWel(model, ipakcb=66, stress_period_data=wel_data, unitnumber=78)\n",
    "\n",
    "lmt = flopy.modflow.ModflowLmt(model, output_file_name='mt3d_link.ftl')\n",
    "\n",
    "stress_period_data = {}\n",
    "for kper in range(nper):\n",
    "    for kstp in range(nstp[kper]):\n",
    "        stress_period_data[(kper, kstp)] = [\"save head\", \"save drawdown\", \n",
    "                                            \"save budget\", \"print head\", \"print budget\",]\n",
    "oc = flopy.modflow.ModflowOc(model, stress_period_data=stress_period_data)#, compact=True)\n",
    "\n",
    "pcg = flopy.modflow.ModflowPcg(model=model)\n",
    "\n",
    "model.write_input()\n",
    "success, mfoutput = model.run_model(silent=True, pause=False, report=True)\n",
    "if not success:\n",
    "    raise Exception('MODFLOW did not terminate normally.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R2 score:  0.9999999403953552\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAADaCAYAAADZqVrgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3FUlEQVR4nO3deVyU9fbA8c8zw7AKyCIigop7buCSCoJopYmKe26VmlpqLvjLbmX3pla3LLt5A3HLXLLcyi03TC0XUNxBcd8VF8QVENnn+f1BciNNQRmGGc779ZrXK555HuZ8tTOeMzPfM4qqqipCCCGEKNM0xg5ACCGEEMYnBYEQQgghpCAQQgghhBQEQgghhEAKAiGEEEIgBYEQQgghkIJACCGEEEhBIIQQQgikIBBCCCEEUhCUeZMmTUJRFG7evGnsUID/xSOE+J/ffvuNZs2aYWdnh6IorF69GoBly5ZRv359bGxsUBSFuLi4p86hQYMGUa1ataeKr7DXtmnThgYNGjzVYxhKmzZtaNOmjbHDKBWkIBBCiFJMVVV69+6NTqdjzZo1xMTEEBQUxI0bN3j99depUaMGGzduJCYmhtq1azN06FBiYmKK/DgfffQRq1atMsAKhKmwMHYAQggh/t7Vq1e5ffs23bt358UXX8w/vnPnTrKzs3nttdcICgrKP25ra4unp2eRH6dGjRrFEq8wXfIKgQDg+vXr9OvXD0dHRypWrMjgwYNJTk7Ov19VVWbMmIGvry82NjY4OTnRq1cvzp07V+D3bN68ma5du+Lp6Ym1tTU1a9Zk2LBhj3xLYv369fj6+mJlZYW3tzf/+c9/DL5OIUqT6OhoXnzxRezt7bG1tcXf35/169fn3z9p0qT8f9zff/99FEWhWrVqDBo0iICAAAD69OmDoij5L3v/3VsGixcvxs/Pj3LlylGuXDl8fX2ZO3du/v2Petl/+vTptG7dGjc3N+zs7GjYsCFTpkwhOzv7mda9b98+AgMDsbW1pXr16nzxxRfo9foC56SkpPDuu+/i7e2NpaUllStXZuzYsaSlpT1VjKqqMmXKFKpWrYq1tTVNmjQhMjLymdZhbuQVAgFAz5496dOnD0OGDCE+Pp7x48cDMG/ePACGDRvGggULGDNmDF9++SW3b9/mk08+wd/fn0OHDlGxYkUAzp49i5+fH0OHDsXR0ZELFy4wdepUAgICiI+PR6fTAXnviXbt2hU/Pz+WLl1Kbm4uU6ZM4fr168b5AxCihG3fvp127drRqFEj5s6di5WVFTNmzCAkJIQlS5bQp08fhg4dio+PDz169GD06NH0798fKysrHBwcaN68OSNHjuTzzz+nbdu2ODg4/O1jTZgwgU8//ZQePXowbtw4HB0dOXLkCBcvXnxsjGfPnqV///75/ygfOnSIzz77jBMnTuQ/NxRVYmIir776KuPGjWPixImsWrWK8ePH4+HhwYABAwC4f/8+QUFBXL58mQ8//JBGjRpx9OhRJkyYQHx8PFu2bMkvegob48cff8zHH3/MkCFD6NWrFwkJCbz55pvk5uZSp06dp1qL2VFFmTZx4kQVUKdMmVLg+Ntvv61aW1urer1ejYmJUQH166+/LnBOQkKCamNjo7733nuP/N16vV7Nzs5WL168qALqL7/8kn9fixYtVA8PDzU9PT3/WEpKiurs7KzK/5aiLGjZsqXq5uampqam5h/LyclRGzRooHp6eqp6vV5VVVU9f/68CqhfffVVgeu3bt2qAurPP/9c4PiDnH7g3LlzqlarVV999dXHxjNw4EC1atWqf3t/bm6ump2drS5cuFDVarXq7du3C33tA0FBQSqg7tmzp8DxevXqqS+//HL+z5MnT1Y1Go26b9++AuctX75cBdQNGzYUKcY7d+6o1tbWavfu3Qucv3PnThVQg4KCnhh7WSBvGQgAunTpUuDnRo0akZGRQVJSEuvWrUNRFF577TVycnLyb+7u7vj4+LBt27b865KSkhg+fDheXl5YWFig0+moWrUqAMePHwcgLS2Nffv20aNHD6ytrfOvtbe3JyQkxPCLFcLI0tLS2LNnD7169aJcuXL5x7VaLa+//jqXL1/m5MmTxfJYmzdvJjc3l5EjRxb52tjYWLp06YKLiwtarRadTseAAQPIzc3l1KlTTxWPu7s7zZs3L3CsUaNGBV6tWLduHQ0aNMDX17fAc87LL7+MoigFnnMKE2NMTAwZGRm8+uqrBR7X398///lJyFsG4g8uLi4FfraysgIgPT2d69evo6pq/tsCf1W9enUA9Ho97du35+rVq3z00Uc0bNgQOzs79Ho9LVu2JD09HYA7d+6g1+txd3d/6Hc96pgQ5ubOnTuoqkqlSpUeus/DwwOAW7duFctj3bhxA6DIHzS8dOkSgYGB1KlTh7CwMKpVq4a1tTV79+5l5MiR+flcVH99roG855s//77r169z5syZ/LcY/+rBZ5IKG+ODP0t5znk8KQjEE7m6uqIoClFRUfmFwp89OHbkyBEOHTrEggULGDhwYP79Z86cKXC+k5MTiqKQmJj40O961DEhzI2TkxMajYZr1649dN/Vq1eBvLwrDhUqVADg8uXLeHl5Ffq61atXk5aWxsqVKwt00XFxccUS1+O4urpiY2Pzt59TePBnU9gYHxQhf/ec87TzF8yNvGUgnqhz586oqsqVK1do1qzZQ7eGDRsC5H/I569Fw+zZswv8bGdnR/PmzVm5ciUZGRn5x1NTU1m7dq2BVyOE8dnZ2dGiRQtWrlxZoDPW6/X8+OOPeHp6Urt27WJ5rPbt26PVapk5c2aRrntUPquqypw5c4olrsfp3LkzZ8+excXF5ZHPOQ/+AS9sjC1btsTa2ppFixYVOL5r164nfrCyLJFXCMQTtWrVirfeeos33niD/fv307p1a+zs7Lh27RrR0dE0bNiQESNGULduXWrUqMEHH3yAqqo4Ozuzdu1aNm/e/NDv/PTTT+nQoQPt2rVj3Lhx5Obm8uWXX2JnZ8ft27eNsEohStbkyZNp164dbdu25d1338XS0pIZM2Zw5MgRlixZUmwTO6tVq8aHH37Ip59+Snp6ev724mPHjnHz5k0+/vjjR17Xrl07LC0t6devH++99x4ZGRnMnDmTO3fuFEtcjzN27FhWrFhB69at+b//+z8aNWqEXq/n0qVLbNq0iXHjxtGiRYtCx+jk5MS7777Lv//9b4YOHcorr7xCQkICkyZNkrcM/kReIRCFMnv2bCIiItixYwd9+/alU6dOTJgwgbS0tPwPCOl0OtauXUvt2rUZNmwY/fr1IykpiS1btjz0+9q1a8fq1atJSUmhT58+vPPOO/Ts2ZPBgweX9NKEMIqgoCB+//137OzsGDRoEH379iU5OZk1a9bQp0+fYn2sTz75hIULF3Lx4kVeffVVunXrxvz58/H29v7ba+rWrcuKFSu4c+dO/rZHX19fwsPDizW2R7GzsyMqKopBgwbx7bff0qlTJ3r37k14eDienp75rxAUJcZPPvmEyZMns2nTJrp06cK0adOYNWuWbDn8E0VVVdXYQQghhBDCuOQVAiGEEEJIQSCEEEIIKQiEEEIIgRQEQgghhEAKAiGEEEIgBYEQQgghMJHBRHq9nqtXr2Jvb19swzqEKItUVSU1NRUPDw80mtLbD0jOC1F8Cpv3JlEQXL16tUgzuIUQj5eQkFDkL7spSZLzQhS/J+W9SRQE9vb2QN5iHBwcjByNEKYrJSUFLy+v/JwqrSTnhSg+hc17kygIHrxk6ODgIE8OQhSD0v4yvOS8EMXvSXlfet9EFEIIIUSJkYJACCGEEOZVEGRk5xo7BCFECcrIzkWvl+9nE6I4mE1BsP3UDdp8tY1dZ24aOxQhRAn5cGU8Q77fx817mcYORQiTZzYFwaxtZ0lMyeDVuXuYsvEE2bl6Y4ckhDCgCzfTWBd/ja0nbxAcFkX0aWkGhHgWZlMQzB3UjH7NvVBVmLHtLL1nx5Bw+76xwxJCGEg1VzvWjGpFLbdy3EjN5PV5e/hSmgEhnprZFAS2lhZM7tGI6f2bYG9tQeylu3QMj2Ld4avGDk0IYSB13R1YMyqA/i2qoKowc9tZXpklzYAQT8NsCoIHOjWqxIYxgTSpUp7UjBxGLY7lgxWHuZ+VY+zQhBAGYGOp5fPuDZnxahMcrC2IS7hLx7Ao1hySZkCIojC7ggDAy9mWn4b5MaptTRQFlu5LIGRaNMevpRg7NCGEgXRsWIkNoYE0q+pEamYOY5bE8t7yQ9IMCFFIZlkQAFhoNbz7ch0WDWmBm70VZ2+k0XX6ThbGXEBVZZuSEObI08mWpW+1ZMwLec3AT/sv03laNEevJhs7NCFKPbMtCB7wr+nKxrGtebGuG1k5eib8cpQ3Fx7gTlqWsUMTQhiAhVbDO+3rsHhoSyo6WHHuRhrdp+9i/s7z0gwI8RhmXxAAONtZ8t3AZkwMqYelVsOW49cJDoti97lbxg5NCGEgfjVciAxtzUvPuZGVq+fjtcd4c+F+bkszIMQjlYmCAPK+1OGNVt6sfNuf6q52JKZk0H/ObqZuPkWObFMSwiw521kyZ0AzJuU3A0kEh+1g11mZWSDEX5WZguCBBpUdWTs6gN7NPNGrEP7bafp+u5srd9ONHZoQwgAURWFQK29Wj2xFjQp2XE/J5NXv9vD1ppPSDAjxJ2WuIACws7JgSi8fwvr6Us7Kgv0X7xD8zQ4i468ZOzQhhIHU83Bg7egA+jTLG2A27fcz9Pl2N5fvyMwCIaCMFgQPdPWtzIYxgfh4lSclI4cRiw7y4ap4+ZIkIcyUraUFX/ZqxLR+jbG3suDAxTsEh0WxQZoBIcp2QQBQxcWW5cP9GB5UA4DFey7RJSKak4mpRo5MCGEoIT4ebAgNpPEfA8zeXnSQ8SvjSc+SZkCUXWW+IADQaTV8EFyXH4Y0p4K9Faeu36NLRDQ/7r4o25SEMFMPBpi93aYGigJL9uY1AycSZYCZKJukIPiTwFoViAwNJKh2BTJz9Pxr9RFG/HiQu/dlm5IQ5kin1fBeh7r8OKQFFeytOJ10jy4RO/lBBpiJMkgKgr9wLWfF/EHP869Oz6HTKmw8mkjHsCj2Xbht7NCEEAbSqqYrG0MDaVunAlk5ej765SjDfzwgzYAoU6QgeASNRmFoYHVWjmhFNRdbriZn0Gd2DGFbTpOrl65BCHPkUs6KeYOe56PO9dBpFX49mjfAbI8MMBNlhBQEj9HQ05F1YwLp0aQyehX+u+UU/ebs5lqyzCwQwhwpisKQAG9Wvd0Kb1c7riVn0G/Obv4rA8xEGSAFwROUs7Jgam9f/tvHBztLLXvP3yY4LIpNRxONHZoQwkAaVHZk3egAejbJG2AW9ttp+s/Zw1UZYCbMmBQEhdS9sSfrxwTSsLIjd+9n89YPB5j4yxGZWSCEmbKzsuDr3j580ydvgNneC3nNwMYj0gwI8yQFQRFUc7VjxQh/3mpdHYDvYy7SbfpOziTJzAIhzFW3xpVZPyYAH09HktOzGf7jAf61WgaYCfMjBUERWVpo+LDjcyx443lcy1lyIjGVztOiWbr3kmxTEsJMVXWx4+fh/gz7oxn4cfcluk3fyenr0gwI8yEFwVNqU8eNDaGBBNZyJSNbzwcr4xm1JJbk9GxjhyaEMABLCw3jOz7HwsHN85uBkIhoFu+RZkCYBykInoGbvTXfv9Gc8cF1sdAorD98jY5hURy4eMfYoQkhDKR17QpEhramde0KZGTr+XBVPCMXHyT5vjQDwrRJQfCMNBqFYUE1WD7CnyrOtly5m07v2TFM33pGZhYIYaYq2FuxYNDzfNgxrxnYEJ9Ix/AoDlyUAWbCdElBUEx8vcqzfkwAXXw8yNWrfPXrSV77bg/XUzKMHZoQwgA0GoW3WtdgxQh/qro8aAZ2M+03GWAmTJMUBMXI3lpHWF9fvurVCFtLLTHnbtHhmx38dvy6sUMTQhiIj1d51o0OoJtvXjPw9eZTvPrdbhKTpRkQpkUKgmKmKAqvNPNi7egA6ns4cOd+NkO+38/Ha4+SmSPblIQwR/bWOr7p25ivX/HB1lLL7nO3CQ7bwZZj0gwI0yEFgYHUqFCOlW/780aragDM33mB7tN3cfbGPeMGJoQwmJ5NPVn3p2Zg6ML9TFpzVGYWCJMgBYEBWVlomRhSn7kDm+FsZ8mxaymETIvm5/0Jsk1JCDNV/Y9mYGiANwALdl2g+4xdnEmSZkCUblIQlIAXn6tIZGgg/jVcuJ+Vyz+WHyZ0aRypGbJNSQhzZGWh5V+d6zF/0PO42Fly/I9m4Kd90gyI0ksKghJS0cGaH4a04B8v10GrUVhz6Codw6OIS7hr7NCEEAbStq4bkaGBtKrpQnp2Lu+tOMzoJbGkSDMgSiEpCEqQVqMwsm1NfhrmR+XyNiTcTqfXzF3M2n4WvWxTEsIsuTlY88PgFrzfIW9mwbrD1+gUHkXsJRlgJkoXKQiMoGlVJzaEBtKpUSVy9CpfRJ5g4Py9JKXKNiUhzJFGozCiTQ1+Gu6Hp1NeM/DKrBhmbDsjzYAoNaQgMBJHGx0R/RrzZc+GWOs0RJ2+SfA3UWw7mWTs0IQQBtKkSl4z0PmPZmDKxpO8Pm8PSTLATJQCUhAYkaIo9Hm+CutGB1DX3Z5baVkMmr+Pz9YfIytHb+zwhBAG4GCtY1q/xkzp2QgbnZadZ24RHBbFVmkGhJFJQVAK1HSzZ/XIVgz0qwrAnKjz9Jy5i/M304wcmRDCEBRFoffzeQPMnqvkwK20LN6Yv49P1x2TAWbCaIpcEOzYsYOQkBA8PDxQFIXVq1cX+tqdO3diYWGBr69vUR/W7FnrtHzctQHfvt6U8rY64q8k0zk8ipUHLxs7NFHGSc4bTk23cqx6259B/tUAmBud1wyckwFmwgiKXBCkpaXh4+NDREREka5LTk5mwIABvPjii0V9yDKlfX13IkMDaeHtTFpWLu/8dIh3lsVxLzPH2KGJMkpy3rCsdVomdanPdwOa4WSr48iVFDpPi2bFgcsys0CUKEV9hv/jFEVh1apVdOvW7Ynn9u3bl1q1aqHValm9ejVxcXGFfpyUlBQcHR1JTk7GwcHhacM1Kbl6lelbz/DNllPoVajmYsu0fk1o6Olo7NCECXvWXJKcN6zE5AzGLotl97m8r1Hu5uvBp90aYG+tM3JkwpQVNp9K5DME8+fP5+zZs0ycOLFQ52dmZpKSklLgVtZoNQpjXqzFsmF+eDhac+HWfXrM3MmcHedkm5Io9STnn467ozWLhrbk3fa10WoUVsddpfO0aA7JADNRAgxeEJw+fZoPPviARYsWYWFhUahrJk+ejKOjY/7Ny8vLwFGWXs9XcyYytDUd6ruTnavy2YbjvLFgHzdSM40dmhCPJDn/bLQahVEv1OKnYS2pXN6Gi7fu03PmLr7dIQPMhGEZtCDIzc2lf//+fPzxx9SuXbvQ140fP57k5OT8W0JCggGjLP0cbXXMfK0Jn3VvgJWFhu2nbhAcFkXU6RvGDk2IAiTni0/Tqs5sGBNIcAN3cvQqn284wSBpBoQBGfQzBHfv3sXJyQmtVpt/TK/Xo6oqWq2WTZs28cILLzzxccrq+4mPcjIxldFLDnLqet6nkIcFVefd9nXQaWUHqXgyQ3+GQHK++KmqytJ9CXy89igZ2Xpcy1kytbcvrWtXMHZowkQUNp8K93reU3JwcCA+Pr7AsRkzZvD777+zfPlyvL29DfnwZqmOuz1rRgXw6bpjLNpzidnbz7H73G2m9W1MFRdbY4cnyjjJ+eKnKAr9mlehWVUnRi2O5eT1VAbM28uw1tUZ174OlhbSDIjiUeSC4N69e5w5cyb/5/PnzxMXF4ezszNVqlRh/PjxXLlyhYULF6LRaGjQoEGB693c3LC2tn7ouCg8a52Wz7o3JLCWK+8tP8yhhLt0DI/is+4N6Opb2djhCTMjOV861Kpozy+jWvHZ+uP8sPsis3ecY/e5W4T3a0xVFztjhyfMQJFLy/3799O4cWMaN24MwDvvvEPjxo2ZMGECANeuXePSpUvFG6V4pA4NKhE5tjXPV3PiXmYOoUvj+MfPh0iTmQWiGEnOlx7WOi2fdmvA7Neb4mij49DlZDqFR7M69oqxQxNm4Jk+Q1BS5P3Ex8vJ1RP++xkifj+NXoXqFewI79uYBpVlZoEoyFRyyVTiNKard9MZuzSOvRfyZhb0bOLJJ13rY2dl0HeChQkqVXMIhGFZaDW80642i99sibuDNedupNFjxi7mRZ+XSWdCmCmP8jYsfrMFY1+qhUaBFQcv03laNEeuJBs7NGGipCAwIy2ruxAZGki7ehXJytXzybpjDPl+P7fuyTYlIcyRhVbD2Jdqs/QtPyo5WnP+ZhrdZ+xkrjQD4ilIQWBmnOws+fb1pnzStT6WFhp+P5FEcFgUu87cNHZoQggDae7tTGRoIO3rVSQ7V+XTdccYvGCfNAOiSKQgMEOKojDArxq/jGxFTbdyJKVm8urcPXz16wmyc/XGDk8IYQDlbS2Z/XpTPu3WAEsLDVtP3qBDWBQ7pRkQhSQFgRl7rpIDa0a1ol9zL1QVpm89S5/ZMSTcvm/s0IQQBqAoCq+3rMqaUa2o5VaOG6mZvDZ3D1M2SjMgnkwKAjNna2nB5B6NiOjfGHtrCw5eyptZsP7wNWOHJoQwkLruDqwZFUC/5lVQVZix7SyvzJJmQDyeFARlROdGHmwYE0jjKuVJzchh5OKDfLDiMOlZucYOTQhhADaWWib3aMiMV5vgYG1BXMJdOoZFsfbQVWOHJkopKQjKEC9nW34a5sfItjVQFFi6L4GQiGiOXyubXzUrRFnQsWElNoQG0rSqE6mZOYxeEsv7yw9zP0sGmImCpCAoY3RaDf94uS6LhrTAzd6KM0n36Dp9JwtjLsg2JSHMlKeTLcveasnoF2qiKLBsfwKdp0Vz9KrMLBD/IwVBGeVf05XI0EBeqOtGVo6eCb8c5a0fDnAnLcvYoQkhDMBCq2Fc+zosGtqCig5WnLuRRvfpu1iwU2YWiDxSEJRhLuWsmDuwGRM618NSq2Hzset0DI9i97lbxg5NCGEg/jVciQxtzUvPuZGVq2fS2mO8uXA/t6UZKPOkICjjFEVhcIA3K9/2p7qrHdeSM+g/ZzdTN58iR7YpCWGWnO0smTOgGZNC8pqBLceTCA7bQcxZaQbKMikIBAANKjuydnQArzT1RK9C+G+n6TdnN1fuphs7NCGEASiKwqBW3qwe2YrqFey4npJJ/+928/Wmk9IMlFFSEIh8dlYWfPWKD2F9fSlnZcG+C3foGBbFxiMys0AIc1XPw4F1owPo3cwTVYVpv5+hz7e7uXxHZhaUNVIQiId09a3MhjGB+HiVJzk9m+E/HuSfq+LJyJaZBUKYI1tLC6b08iG8X2PsrSw4cDGvGYiMl2agLJGCQDxSFRdblg/3Y3hQDQAW7blEl4hoTl1PNXJkQghD6eLjwYbQQHy9ypOSkcOIRQcZvzJeBpiVEVIQiL+l02r4ILguPwxpjms5K05dv0fItGgW7bko25SEMFNezrb8PNyPEW3yBpgt2XuJrtOjOZkozYC5k4JAPFFgrQpsHBtIUO0KZObo+eeqI7y96CDJ97ONHZoQwgB0Wg3vd6jLD4NbUME+rxnoEhHND7ulGTBnUhCIQnEtZ8X8Qc/zr07PodMqRB5JJDhsB/su3DZ2aEIIAwmo5crG0EDa1slrBj5afYThPx7g7n2ZWWCOpCAQhabRKAwNrM6KEf5Uc7HlanIGfWbHELblNLl66RqEMEd5A8z+1wz8evQ6wWFR7D0vzYC5kYJAFFkjz/KsGxNIj8aV0avw3y2n6D9nN9eSZWaBEOboQTOw6u1WeP8xwKzvtzF8s+WUNANmRAoC8VTKWVkwtY8vU3v7YGepZc/52wSHRbHpaKKxQxNCGMiDAWY9m+QNMPtmS94As6sywMwsSEEgnkmPJp6sGxNIw8qO3L2fzVs/HGDiL0dkZoEQZqqclQVf9/bhmz6+2Flq2ftHM/CrNAMmTwoC8cy8Xe1YMcKfNwO9Afg+5iLdpu/kTJJsUxLCXHVrXJn1YwJp5OlIcno2w344wEerpRkwZVIQiGJhaaHhn53qMf+N53Gxs+REYioh03aybN8l2aYkhJmq5mrH8uH+DGtdHYAfduc1A6dlgJlJkoJAFKu2ddyIHBtIYC1X0rNzeX9FPKOWxJKcLjMLhDBHlhYaxnd8ju8HN8e13B/NQEQ0i/dIM2BqpCAQxc7N3prv32jOB8F1sdAorD98jU7hURy4eMfYoQkhDCSodgUiQ1sTWMuVjGw9H66KZ9RiaQZMiRQEwiA0GoXhQTX4ebgfXs42XL6TTu/ZMUzfeka2KQlhpirYW/H9G80Z/6AZiL9Gx7AoDlyUmQWmQAoCYVCNqzixfkwgXXw8yNWrfPXrSV6fu4frKRnGDk0IYQAajcKwoBqsGOFPVRdbrtxNp/fs3UT8LgPMSjspCITBOVjrCOvry1e9GmGj07Lr7C2Cw6L4/cR1Y4cmhDAQH6/yrBsdQDffvGbgP5tO8dp3e0hMlmagtJKCQJQIRVF4pZkX68YEUK+SA7fTshi8YD+frD1GZo5sUxLCHNlb6/hvH1/+84oPtpZaYs7dIjhsB78dl2agNJKCQJSoGhXKsWqkP2+0qgbAvJ3n6TFjF+du3DNuYEIIg1AUhV5NPVk3OoD6Hg7cuZ/NkO/3M2nNUZlZUMpIQSBKnJWFlokh9Zk7sBnOdpYcvZpC52nR/Lw/QbYpCWGmqlcox8q3/RkSkDfAbMGuC/SYsYszSdIMlBZSEAijefG5ikSGBuJX3YX7Wbn8Y/lhQpfGkZoh25SEMEdWFlo+6lyP+YOex9nOkmPXUgiZFs1P+6QZKA2kIBBGVdHBmh+HtuAfL9dBq1FYc+gqncKjiUu4a+zQhBAG0rauGxtDA2lV04X07FzeW3GYMUvjSJFmwKikIBBGp9UojGxbk5+G+VG5vA2Xbt+n18xdzNp+Fr1sUxLCLLk5WPPD4Ba81yGvGVh76CqdwqOIvSQDzIxFCgJRajSt6sSG0EA6NaxEjl7li8gTDJy/l6RU2aYkhDnSaBTeblOTn4f74elkQ8LtdF6ZFcOMbWekGTACKQhEqeJooyOif2O+6NEQa52GqNM36RgWxfZTN4wdmhDCQJpUyWsGOjfKawambDzJgHl7SZIBZiWqyAXBjh07CAkJwcPDA0VRWL169WPPX7lyJe3ataNChQo4ODjg5+fHr7/++rTxijJAURT6Nq/CutEB1HW35+a9LAbO28vnG46TlaM3dnhljuS8KAkO1jqm9WvMlJ55A8yiz9wkOCyKrSeTjB1amVHkgiAtLQ0fHx8iIiIKdf6OHTto164dGzZs4MCBA7Rt25aQkBBiY2OLHKwoW2q62bN6ZCsG+lUF4Nsd5+g5cxcXbqYZObKyRXJelBRFUej9vBdrR7fiuUoO3ErL4o35+/h0nQwwKwmK+gx7PRRFYdWqVXTr1q1I19WvX58+ffowYcKEQp2fkpKCo6MjycnJODg4PEWkwtRtOprIeysOc/d+NnaWWv7dvQHdG3saOyyT86y5JDkvSkpGdi5fRJ5gwa4LADSo7MC0fk3wdrUzbmAmqLD5VOKfIdDr9aSmpuLs7FzSDy1MWPv67kSGBtLc25m0rFz+b9kh3vkpjnuZOcYOTTyB5Lx4GtY6LZO61Oe7Ac1wstVx5EoKncKjWHHgsrFDM1slXhB8/fXXpKWl0bt37789JzMzk5SUlAI3ISo52rDkzZa80642GgVWHrxC5/Ao4i8nGzs08RiS8+JZvFSvIpGhrWlZ3Zn7WbmM+/kQ/7dMmgFDKNGCYMmSJUyaNIlly5bh5ub2t+dNnjwZR0fH/JuXl1cJRilKM61GYcyLtVg2zA8PR2su3LpPj5k7+S7qnGxTKoUk50VxcHe0ZtHQloxrVxutRmFV7BU6hUdx+PJdY4dmVkqsIFi2bBlDhgzhp59+4qWXXnrsuePHjyc5OTn/lpCQUEJRClPxfDVnNoQG8nL9imTnqvx7/XHeWLCPm/cyjR2a+IPkvChOWo3C6BdrseytllQub8PFW/fpMWMX3+6QAWbFpUQKgiVLljBo0CAWL15Mp06dnni+lZUVDg4OBW5C/FV5W0tmvdaUf3drgJWFhu2nbhAcFkX06ZvGDq3Mk5wXhtKsmjMbxgQS3MCdHL3K5xtOMGjBPm6kSjPwrIpcENy7d4+4uDji4uIAOH/+PHFxcVy6dAnIq/QHDBiQf/6SJUsYMGAAX3/9NS1btiQxMZHExESSk+V9X/HsFEXhtZZVWTMqgNoVy3EjNZPX5+3hi8gTZOfKzILiIDkvShtHWx0zXm3C590bYmWhYccfzcAOGWD2bNQi2rp1qwo8dBs4cKCqqqo6cOBANSgoKP/8oKCgx55fGMnJySqgJicnFzVcUYbcz8xRx688rFZ9f51a9f11apeIaPXizTRjh1WqPE0uSc6L0uxkYorafur2/Lz/fP0xNTM719hhlSqFzadnmkNQUmRPsiiKyPhrvL/iMCkZOdhbWfBZj4Z08fEwdlilgqnkkqnEKUqHjOxcPlt/nB92XwTAx9OR8H6NqeoiMwugFM8hEMLQghtWInJsa5pVdSI1M4cxS2L5x8+HuJ8l25SEMEfWOi2fdmvArNea4mij49DlZDqFR/NL3BVjh2ZSpCAQZqlyeRuWvtWSMS/URFHg5wOX6TwtmqNX5X1sIcxVhwbubAgN5PlqTtzLzCF0aRzjfjpEmswsKBQpCITZstBqeKd9HRYPbYm7gzXnbqTRffou5u88jwm8UyaEeAqVy+cNMAt9sRYaBVYcvEzItGiOXJFm4EmkIBBmz6+GC5Ghgbz0XEWycvV8vPYYQ7/fz+20LGOHJoQwAAuthv9rV5slb7akkqM1526m0WPGLuZGSzPwOFIQiDLByc6SOQOa8knX+lhaaPjtRBIdvtnBrrMys0AIc9WiugsbxgTSvl5eM/DpumMMXrCPWzLA7JGkIBBlhqIoDPCrxuq3W1Gjgh1JqZm8+t0e/vPrSXJkZoEQZsnJzpLZrzfl0z+aga0n82YW7DojzcBfSUEgypx6Hg6sHR1A3+e9UFWI2HqG3rNjSLh939ihCSEMQFEUXverxppRrajpVi6vGZi7hykbZYDZn0lBIMokW0sLvujZiIj+jbG3tuDgpbt0DI9i/eFrxg5NCGEgdd0dWDsqgH7Nq6CqMGPbWWkG/kQKAlGmdW7kwYYxgTSuUp7UjBxGLj7I+JWHSc/KNXZoQggDsLHUMrlHQ6b3b4K9tQWxl+7SMSyKtYeuGjs0o5OCQJR5Xs62/DTMj5Fta6AosGRvAiER0Ry/lmLs0IQQBtKpUSUiQwNp+scAs9FLYnl/+eEyPcBMCgIhAJ1Wwz9ersuPQ1rgZm/FmaR7dJ2+k4UxF2SbkhBmytPJlmVvtWT0HwPMlu1PIGRaNMeuls1mQAoCIf6kVU1XIkMDeaGuG1k5eib8cpS3fjjAHZlZIIRZstBqGNe+DouGtqCigxVnb6TRbcZOvt9V9poBKQiE+AuXclbMHdiMCZ3rYanVsPnYdTqGR7Hn3C1jhyaEMBD/Gq5EhrbmxT+agYlrjvLmwrLVDEhBIMQjKIrC4ABvVr7tT3VXO64lZ9Bvzm7+u/mUzCwQwkw521ny3cBmTArJawa2HL9OcFgUMWfLRjMgBYEQj9GgsiNrRwfQq6knehXCfjtN/zl7uHI33dihCSEMQFEUBrXyZtVIf6pXsCMxJYP+3+1m6ibzH2AmBYEQT2BnZcF/XvEhrK8v5aws2HvhNh3Doth4JNHYoQkhDKS+hyPrRgfQu5knqgrhv5+h77e7uXzHfGcWSEEgRCF19a3M+jEB+Hg6kpyezfAfD/Cv1fFkZMvMAiHMka2lBVN65TUD9lYW7L94h45hUUTGm+cAMykIhCiCqi52/Dzcn2FB1QH4cfclukbs5NT1VCNHJoQwlLxmIBBfr/KkZOQwYtFBPlwVb3YDzKQgEKKILC00jA9+joWDm+NazoqT11MJmRbNoj0Xy9w2JSHKiioutvw83I8RbfIGmC3ec4mu06M5mWg+zYAUBEI8pda1KxAZGkjr2hXIzNHzz1VHeHvRQZLvZxs7NCGEAei0Gt7vUJcfBreggr0Vp67fo0tEND/uNo9mQAoCIZ5BBXsrFgx6nn92fA6dViHySCIdw6PYf+G2sUMTQhhIQK28AWZt6uQ1A/9afYQRPx7k7n3TnlkgBYEQz0ijUXizdXVWjPCnmostV+6m03t2DOG/nSZXb/pdgxDiYa7lrJg38Hn+1SmvGdh4NJGOYVHsPW+6zYAUBEIUk0ae5Vk3JpDujSujV2Hq5lP0n7Oba8kys0AIc6TRKAwNrM7KEa2o5mLL1eQM+n4bQ9gW02wGpCAQohiVs7Lgv318mdrbB1tLLXvO3yY4LIrNx64bOzQhhIE09HRk3ZhAejTJawb+u+UU/ebs5qqJDTCTgkAIA+jRxJP1YwJpUNmBu/ezeXPhfib+ckRmFghhpspZWTC1ty//7eODnaWWvedv0zE8ik1HTWeAmRQEQhiIt6sdK0e0YmiANwDfx1yk2/SdnEkyn21KQoiCujfOawYaeTpy9342b/1wgAkm0gxIQSCEAVlaaPhX53rMf+N5XOwsOZGYSsi0nSzbd8kstikJIR5WzdWO5cP9eat13gCzhX80A6dL+QAzKQiEKAFt67gRGRpIQE1X0rNzeX9FPKOXxJKSITMLhDBHlhYaPuz4HN8Pbo5ruT+agYholuwtvc2AFARClBA3B2sWDm7O+x3qYqFRWHf4Gh3Dojh46Y6xQxNCGEhQ7QpsCA0ksJYrGdl6xq+MZ9TiWJLTS18zIAWBECVIo1EY0aYGPw/3w8vZhst30nllVgzTt55Bb4LblIQQT+Zmb833bzRnfHBeM7A+Pq8ZOHCxdM0skIJACCNoXMWJ9WMCCfHxIFev8tWvJ3l93h6up2QYOzQhhAFoNArDgmqwfIQ/VZwfDDDbTcTvpWdmgRQEQhiJg7WO8L6+TOnVCBudlp1nbhEcFsXWE0nGDk0IYSC+XuVZPyaArr55zcB/Np3ite/2kJhs/GZACgIhjEhRFHo382LdmADqVXLgdloWbyzYxydrj5GZU/q3KQkhis7eWsc3fXz5zyt5A8xizt0iOGwHvx037gAzKQiEKAVqVCjHqpH+vNGqGgDzdp6nx4xdnLtxz7iBCSEMQlEUejX1ZN3oAOp7OHDnfjZDvt/PpDVHjdYMSEEgRClhZaFlYkh95g5shpOtjqNXU+g8LZrlBy6X2m1KQohnU71COVa+7c/gVnkDzBbsukD36bs4a4RmQAoCIUqZF5+rSGRoa/yqu3A/K5d3fz7E2GVxpMrMAiHMkpWFlgkh9Zg3qBnOdpYcu5ZC5/BoftqfUKLNgBQEQpRC7o7W/Di0Bf94uQ5ajcIvcVfpFB5NXMJdY4cmhDCQF+pWJDI0EP8aLqRn5/Le8sOELo0rsQFmUhAIUUppNQoj29bkp2EtqVzehku379Nr5i5mbT8rMwuEMFMVHaz5Ycj/moE1h67SKTyK2BIYYCYFgRClXNOqzmwIDaRTw0rk6FW+iDzBwPl7SUo1/jYlIUTxe9AM/DzcD08nGxJu5w0wm7nNsM1AkQuCHTt2EBISgoeHB4qisHr16ides337dpo2bYq1tTXVq1dn1qxZTxOrEGWWo42OiP6NmdyjIdY6DVGnb9IxLIrtp24Y/LEl54UwjiZ/DDDr1CivGfhy4wkGzNtLkoEGmBW5IEhLS8PHx4eIiIhCnX/+/Hk6duxIYGAgsbGxfPjhh4wZM4YVK1YUOVghyjJFUejXvAprRwVQ192em/eyGDhvL59vOE5Wjt5gjys5L4TxONroiOjXmC975jUD0Wdu5g0wO1n8A8wU9Rk+wqgoCqtWraJbt25/e87777/PmjVrOH78eP6x4cOHc+jQIWJiYgr1OCkpKTg6OpKcnIyDg8PThiuE2cjIzuXzDcdZGHMRgEaejoT3bUw1V7vHXvesuSQ5L4TxnElKZdTiWE4k5n2N8tAAb97rUBdLi8f39oXNJ4N/hiAmJob27dsXOPbyyy+zf/9+srMf/cnJzMxMUlJSCtyEEP9jrdPySdcGzH69KY42Og5fTqZTeBSrY68YOzTJeSEMpKabPatHtmKQfzUAvos+T8+Zuzh/M61Yfr/BC4LExEQqVqxY4FjFihXJycnh5s2bj7xm8uTJODo65t+8vLwMHaYQJunl+u5EhgbS3NuZtKxcVhw0/hAjyXkhDMdap2VSl/rMGdCM8rY64q8kc/Bi8exAKJFdBoqiFPj5wRPWX48/MH78eJKTk/NvCQkJBo9RCFPlUd6GJW+25P0Odfm6t8/f5lVJkpwXwrDa1cubWfBBcF16NvUslt9pUSy/5THc3d1JTEwscCwpKQkLCwtcXFweeY2VlRVWVlaGDk0Is6HVKIxoU8PYYQCS80KUlEqONgwPKr68N/grBH5+fmzevLnAsU2bNtGsWTN0Op2hH14IUcIk54UwTUUuCO7du0dcXBxxcXFA3hajuLg4Ll26BOS99DdgwID884cPH87Fixd55513OH78OPPmzWPu3Lm8++67xbMCIYRBSc4LUUaoRbR161YVeOg2cOBAVVVVdeDAgWpQUFCBa7Zt26Y2btxYtbS0VKtVq6bOnDmzSI+ZnJysAmpycnJRwxVC/MnT5JLkvBCmrbD59ExzCEpKcnIy5cuXJyEhQfYkC/EMUlJS8PLy4u7duzg6Oho7nL8lOS9E8Sls3hv8Q4XFITU1bwiDbEUSonikpqaW6oJAcl6I4vekvDeJVwj0ej1Xr17F3t7+b7ctPaiAzK2jkHWZDlNYk6qqpKam4uHhgUZTer/brDA5D6bxZ/40zHFd5rgmMI11FTbvTeIVAo1Gg6dn4fZZOjg4lNq/lGch6zIdpX1NpfmVgQeKkvNQ+v/Mn5Y5rssc1wSlf12FyfvS2yIIIYQQosRIQSCEEEII8ykIrKysmDhxotlNO5N1mQ5zXFNpZ65/5ua4LnNcE5jXukziQ4VCCCGEMCyzeYVACCGEEE9PCgIhhBBCSEEghBBCCCkIhBBCCIEZFQQzZszA29sba2trmjZtSlRUlLFDKrRJkyahKEqBm7u7e/79qqoyadIkPDw8sLGxoU2bNhw9etSIET/ajh07CAkJwcPDA0VRWL16dYH7C7OOzMxMRo8ejaurK3Z2dnTp0oXLly+X4Coe9qR1DRo06KG/v5YtWxY4pzSuyxxI3hufOeZ9Wc15sygIli1bxtixY/nnP/9JbGwsgYGBBAcH5389qymoX78+165dy7/Fx8fn3zdlyhSmTp1KREQE+/btw93dnXbt2uXPey8t0tLS8PHxISIi4pH3F2YdY8eOZdWqVSxdupTo6Gju3btH586dyc3NLallPORJ6wLo0KFDgb+/DRs2FLi/NK7L1Enelw7mmPdlNucN+I2LJaZ58+bq8OHDCxyrW7eu+sEHHxgpoqKZOHGi6uPj88j79Hq96u7urn7xxRf5xzIyMlRHR0d11qxZJRRh0QHqqlWr8n8uzDru3r2r6nQ6denSpfnnXLlyRdVoNOrGjRtLLPbH+eu6VDXv63+7du36t9eYwrpMkeR96WOOeV+Wct7kXyHIysriwIEDtG/fvsDx9u3bs2vXLiNFVXSnT5/Gw8MDb29v+vbty7lz5wA4f/48iYmJBdZnZWVFUFCQSa2vMOs4cOAA2dnZBc7x8PCgQYMGpX6t27Ztw83Njdq1a/Pmm2+SlJSUf58pr6u0krw3Deac9+aY8yZfENy8eZPc3FwqVqxY4HjFihVJTEw0UlRF06JFCxYuXMivv/7KnDlzSExMxN/fn1u3buWvwZTXBxRqHYmJiVhaWuLk5PS355RGwcHBLFq0iN9//52vv/6affv28cILL5CZmQmY7rpKM8l702CueW+uOW8S33ZYGH/9ilRVVR/7tamlSXBwcP5/N2zYED8/P2rUqMH333+f/0EVU17fnz3NOkr7Wvv06ZP/3w0aNKBZs2ZUrVqV9evX06NHj7+9rrSvyxSYcl5I3ptu3ptrzpv8KwSurq5otdqHqq6kpKSHqlJTYWdnR8OGDTl9+nT+p45NfX2FWYe7uztZWVncuXPnb88xBZUqVaJq1aqcPn0aMJ91lSaS96ahrOS9ueS8yRcElpaWNG3alM2bNxc4vnnzZvz9/Y0U1bPJzMzk+PHjVKpUCW9vb9zd3QusLysri+3bt5vU+gqzjqZNm6LT6Qqcc+3aNY4cOWJSa7116xYJCQlUqlQJMJ91lSaS96ahrOS92eS8sT7NWJyWLl2q6nQ6de7cueqxY8fUsWPHqnZ2duqFCxeMHVqhjBs3Tt22bZt67tw5dffu3Wrnzp1Ve3v7/Pi/+OIL1dHRUV25cqUaHx+v9uvXT61UqZKakpJi5MgLSk1NVWNjY9XY2FgVUKdOnarGxsaqFy9eVFW1cOsYPny46unpqW7ZskU9ePCg+sILL6g+Pj5qTk6OsZb12HWlpqaq48aNU3ft2qWeP39e3bp1q+rn56dWrly51K/L1Enelw7mmPdlNefNoiBQVVWdPn26WrVqVdXS0lJt0qSJun37dmOHVGh9+vRRK1WqpOp0OtXDw0Pt0aOHevTo0fz79Xq9OnHiRNXd3V21srJSW7durcbHxxsx4kfbunWrCjx0GzhwoKqqhVtHenq6OmrUKNXZ2Vm1sbFRO3furF66dMkIq/mfx63r/v37avv27dUKFSqoOp1OrVKlijpw4MCHYi6N6zIHkvfGZ455X1ZzXr7+WAghhBCm/xkCIYQQQjw7KQiEEEIIIQWBEEIIIaQgEEIIIQRSEAghhBACKQiEEEIIgRQEQgghhEAKAiGEEEIgBYEQQgghkIJACCGEEEhBIIQQQgikIBBCCCEE8P/IeZQIjJaZIgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head = flopy.utils.binaryfile.HeadFile(model_ws+'/'+model_name+\".hds\").get_alldata()\n",
    "head_official = flopy.utils.binaryfile.HeadFile(\"./official_examples/ex06_pht3d/mf2k/heads.dat\").get_alldata()\n",
    "\n",
    "plt.figure(figsize=(6,2))\n",
    "plt.subplot(121)\n",
    "plt.plot(head[-1, 0, 0])\n",
    "plt.title(\"head\")\n",
    "plt.subplot(122)\n",
    "plt.plot(head_official[-1, 0, 0])\n",
    "plt.title(\"official head\")\n",
    "\n",
    "print(\"R2 score: \", r2_score(np.ravel(head_official), np.ravel(head)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'./models_folder/case6\\\\pht3d_ph.dat'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Create PHT3D Reactive Transport Model, supported by PHT3D-FSP\n",
    "## Define \"spec\", calling PHT3D-FSP (adopt function variables as needed)\n",
    "spec = pht3d_fsp.create(\n",
    "    xlsx_path=\"./\", \n",
    "    xlsx_name=f\"./data/{data_folder}/pht3d_species.xlsx\",\n",
    "    nlay=nlay, \n",
    "    nrow=nrow, \n",
    "    ncol=ncol,\n",
    "    ph_os=2,\n",
    "    ph_temp=25,\n",
    "    ph_asbin=1,\n",
    "    ph_eps_aqu=0,\n",
    "    ph_ph=0,\n",
    "    # ph_print=0,\n",
    "    ph_cb_offset=0,\n",
    "    # ph_surf_calc_type=\"-diffuse_layer\",\n",
    "    write_ph=\"yes\"\n",
    ")\n",
    "# Move the pht3d_ph.dat file into the models_folder\n",
    "source = 'pht3d_ph.dat'\n",
    "destination = os.path.join(model_ws, 'pht3d_ph.dat')\n",
    "shutil.move(source, destination)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: [0,\n",
       "  0,\n",
       "  0,\n",
       "  0.0,\n",
       "  2,\n",
       "  1.42e-05,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.005,\n",
       "  0.005,\n",
       "  7,\n",
       "  4,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.0],\n",
       " 1: [0,\n",
       "  0,\n",
       "  0,\n",
       "  0.0,\n",
       "  2,\n",
       "  1.42e-05,\n",
       "  0.0,\n",
       "  0.005,\n",
       "  0.005,\n",
       "  0.0,\n",
       "  7,\n",
       "  4,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.0]}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define SSM and RCH data (for the SSM package) for each stress period\n",
    "ssm_data = {}\n",
    "itype = flopy.mt3d.Mt3dSsm.itype_dict() # ... to check the key words\n",
    "# print(itype)\n",
    "# {'CHD': 1, 'BAS6': 1, 'PBC': 1, 'WEL': 2, 'DRN': 3, 'RIV': 4, 'GHB': 5, 'MAS': 15, 'CC': -1}\n",
    "\n",
    "wel_spec={}\n",
    "\n",
    "## Initiate RCH data...\n",
    "# rch_spec={}\n",
    "# for key in spec.keys():\n",
    "#    rch_spec[key]={}\n",
    "\n",
    "for i in range(1):\n",
    "    # rch and wel\n",
    "    for key in spec.keys():\n",
    "        # rch_spec[key][i]=spec[key][0,:,:]\n",
    "        wel_spec[key]=spec[key][0,0,0]\n",
    "\n",
    "    wel_spec['c4'] = 0.00001420\n",
    "    wel_spec['ca'] = 0.0\n",
    "    wel_spec['cl'] = 0.0\n",
    "    wel_spec['na'] = 0.005\n",
    "    wel_spec['tenside'] = 0.005\n",
    "    \n",
    "    wel_spec['ph'] = 7\n",
    "    wel_spec['pe'] = 4\n",
    "\n",
    "    wel_spec['calcite'] = 0.0\n",
    "    wel_spec['calciumtenside'] = 0.0\n",
    "\n",
    "    wel_spec['cax2'] = 0.0\n",
    "    wel_spec['nax'] = 0.0\n",
    "\n",
    "    ssm_per=[0,0,0,0.0,itype['WEL']]\n",
    "    for key in spec.keys():\n",
    "        ssm_per.append(wel_spec[key])\n",
    "    \n",
    "    ssm_data[i] = ssm_per\n",
    "\n",
    "for i in range(1,2):\n",
    "    # rch and wel\n",
    "    for key in spec.keys():\n",
    "        # rch_spec[key][i]=spec[key][0,:,:]\n",
    "        wel_spec[key]=spec[key][0,0,0]\n",
    "\n",
    "    wel_spec['c4'] = 0.00001420\n",
    "    wel_spec['ca'] = 0.0\n",
    "    wel_spec['cl'] = 0.005\n",
    "    wel_spec['na'] = 0.005\n",
    "    wel_spec['tenside'] = 0.0\n",
    "    \n",
    "    wel_spec['ph'] = 7\n",
    "    wel_spec['pe'] = 4\n",
    "\n",
    "    wel_spec['calcite'] = 0.0\n",
    "    wel_spec['calciumtenside'] = 0.0\n",
    "\n",
    "    wel_spec['cax2'] = 0.0\n",
    "    wel_spec['nax'] = 0.0\n",
    "\n",
    "    ssm_per=[0,0,0,0.0,itype['WEL']]\n",
    "    for key in spec.keys():\n",
    "        ssm_per.append(wel_spec[key])\n",
    "    \n",
    "    ssm_data[i] = ssm_per\n",
    "\n",
    "ssm_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DSP: setting dmcoef for component 2 to zero, kwarg name dmcoef2\n",
      "DSP: setting dmcoef for component 3 to zero, kwarg name dmcoef3\n",
      "DSP: setting dmcoef for component 4 to zero, kwarg name dmcoef4\n",
      "DSP: setting dmcoef for component 5 to zero, kwarg name dmcoef5\n"
     ]
    }
   ],
   "source": [
    "icbund = np.abs(ibound)\n",
    "\n",
    "## Initiate model object\n",
    "mt = flopy.mt3d.Mt3dms(model_name, model_ws=model_ws, exe_name='./bin/mt3dms.exe',\n",
    "                      ftlfilename='mt3d_link.ftl',modflowmodel=model, namefile_ext='nam_pht3d')\n",
    "\n",
    "prsity = 0.424\n",
    "## Define FloPy packages\n",
    "save_times = np.array([0, 2.083E-03, 4.167E-03, 6.250E-03, 8.333E-03, 1.042E-02, .0125, 1.458E-02, 1.667E-02, .01875, 2.083E-02, 2.292E-02, .025, \n",
    "                       2.708E-02, 2.917E-02, .03125, 3.333E-02, 3.542E-02, 3.750E-02, 3.958E-02, 4.167E-02, 4.375E-02, 4.583E-02, 4.792E-02, 5.000E-02, \n",
    "                       5.208E-02, 5.417E-02, 5.625E-02, 5.833E-02, 6.042E-02, 6.250E-02, 6.458E-02, 6.667E-02, 6.875E-02, 7.083E-02, 7.292E-02, 7.500E-02, \n",
    "                       7.708E-02, 7.917E-02, .08125, 8.333E-02, 8.542E-02, .0875, 8.958E-02, 9.167E-02, 9.375E-02, 9.583E-02, 9.792E-02, 1.000E-01, .1020833, \n",
    "                       .1041666, .10625, .1083333, .1104166, .1125, .1145833, .1166666, .11875, .1208333, .1229166, .125, .1270833, .1291666, .13125, \n",
    "                       .1333333, .1354166, .1375, .1395833, .1416666, .14375, .1458333, .1479166, .1499999, .1520833, .1541666, .1562499, .1583333, .1604166, \n",
    "                       .1624999, .1645833, .1666666, .1687499, .1708333, .1729166, .1749999, .1770833, .1791666, .1812499, .1833332, .1854166, .1874999, .1895832, \n",
    "                       .1916666, .1937499, .1958332, .1979166, .1999999, .2020832, .2041666, .2062499, .2083332])\n",
    "exec(f'btn = flopy.mt3d.Mt3dBtn(mt, \\\n",
    "                                nper=nper, perlen=perlen, nstp=nstp,\\\n",
    "                                nlay=nlay, ncol=ncol, nrow=nrow, \\\n",
    "                                laycon=[1], prsity=prsity, icbund=icbund, nprs=1, mxstrn=50000, tsmult=1, timprs=save_times, \\\n",
    "    ncomp={pht3d_fsp.create.ncomp}, mcomp={pht3d_fsp.create.mcomp}, {pht3d_fsp.create.sconc_btn})') # call BTN package in this way to invoke PHT3D-FSP variables\n",
    "\n",
    "adv = flopy.mt3d.Mt3dAdv(mt, mixelm=-1, percel=0.75, mxpart=5000, nadvfd=0)\n",
    "dsp = flopy.mt3d.Mt3dDsp(mt, al=0.000576, trpt=0.1, trpv=0.1, multiDiff=True, dmcoef=0.0,)\n",
    "gcg = flopy.mt3d.Mt3dGcg(mt, isolve=2, cclose=1.e-10)\n",
    "# exec(f'ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data, mxss=1000, {pht3d_fsp.create.crch_ssm})') # call SSM package in this way to invoke PHT3D-FSP variables\n",
    "exec(f'ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data, mxss=1000)')\n",
    "# kwargs={'sp12':2e-4}\n",
    "# rct = flopy.mt3d.Mt3dRct(mt,rhob=1800.0,ireact=0,isothm=1,sp1=0,igetsc=0,**kwargs)\n",
    "mt.write_input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Manually add the PHC reaction package to the nam file and push out as pht3d.nam\n",
    "s='PHC               39  pht3d_ph.dat\\n'\n",
    "namfiletxt=open(model_ws + '/'+ mt.namefile, 'r').read()\n",
    "pht3d_nam = namfiletxt+s\n",
    "file = open(model_ws + '/pht3d.nam','w')\n",
    "file.write(pht3d_nam)\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "copy_pht3d_datab(source_folder=data_folder)\n",
    "\n",
    "run_pht3d_program(work_dir=data_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAIcCAYAAABrUjh1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAADklUlEQVR4nOzdd3wT9f8H8Ndld+9JS1tm2XvKKnsry4GAgKAy/MoWEKXIRhRwgchPQJkqiDIVoSAqIDjYmwKlpXTvlfH5/XEkNG3aJm0uq+/n43GPJHeXz32STy+9930WxxhjIIQQQgghhJAqEFk7A4QQQgghhBD7R4EFIYQQQgghpMoosCCEEEIIIYRUGQUWhBBCCCGEkCqjwIIQQgghhBBSZRRYEEIIIYQQQqqMAgtCCCGEEEJIlVFgQQghhBBCCKkyCiwIIYQQQgghVWY3gcXx48cxfvx4REZGwsXFBTVq1MCzzz6Lv//+29pZI4QQQgghpNrjGGPM2pkwxogRI5CamooRI0agYcOGSE5Oxocffojz58/j559/Rvfu3a2dRUIIIYQQQqotuwkskpKS4O/vr7cuJycHderUQePGjfHrr79aKWeEEEIIIYQQu2kKVTKoAABXV1c0bNgQcXFxVsgRIYQQQgghRMtuAgtDMjMz8c8//6BRo0bWzgohhBBCCCHVmsTaGaiKKVOmIDc3F++88065+xUWFqKwsFD3WqPRIC0tDT4+PuA4TuhsEkIIIYQQYhWMMWRnZyM4OBgikbB1CnYbWLz77rvYvn07PvnkE7Rq1arcfZcvX45FixZZKGeEEEIIIYTYlri4OISEhAh6DLvpvF3cokWLEB0djaVLl2L+/PkV7l+yxiIzMxM1a9bEzZs34e3tLWRWiYUolUrExMQgKioKUqnU2tkhZkBl6nioTB0PlanjoTJ1PGlpaahXrx4yMjLg4eEh6LHsrsZCG1RER0cbFVQAgFwuh1wuL7Xe29sbPj4+5s4isQKlUglnZ2f4+PjQD6GDoDJ1PFSmjofK1PFQmTouSzT/t6vO24sXL0Z0dDQWLFiAhQsXWjs7hBBCCCGEkCfspsbiww8/xHvvvYe+fftiwIABOHPmjN729u3bWylnhBBCCCGEELsJLPbv3w8AOHLkCI4cOVJqux12FSGEEEIIIcRh2E1gceLECWtngRBCCCGEEFIGuwksrEWtVkOpVFo7G6QCSqUSEokEBQUFUKvV1s6OVUmlUojFYmtngxBCCCHVDAUWZWCMITExEZmZmdTMyg4wxhAYGIi4uLhqP+khx3Hw8PBAYGBgtf8uCCGEEGI5FFiUITMzExkZGfDz84OLiwtdoNk4jUaDnJwcuLq6Cj6rpC1jjCE3NxfJyclwcnKCp6entbNECCGEkGqCAgsDGGNISkqCu7s7fH19rZ0dYgSNRoOioiIoFIpqHVgAgJOTEwoLC5GUlAQPDw8KigkhhBBiEdX7CqwMarUaarUa7u7u1s4KIZXi7u6u+zsmhBBCCLEECiwMUKlUAACJhCp0iH3S/u1q/5YJIYQQQoRGgUU5qAkJsVf0t0sIIYQQS6PAghBCCCGEEFJlFFgQQgghhBBCqowCi2rq4sWLGDduHCIiIqBQKODq6oqWLVti1apVSEtLM+uxvv76a/j5+SE7O1u3Ljw8HBzH6RZXV1e0a9cOX3/9td57u3XrhsaNGxtMNyUlBRzHITo6GgDQtGlTiMVivXQNLVu2bAEATJgwAY0bN4anpyecnJxQr149zJ49GykpKXrH2bJli977FQoFAgMDERUVheXLlyMpKalU3t599120bNkSGo2mCt8cIYQQQoj9oN7J1dCXX36JyZMno379+pg9ezYaNmwIpVKJ8+fPY8OGDTh9+jR++OEHsxwrLy8P8+fPx9tvvw03Nze9bc888wxWr14NAHj48CFWr16NV155Bbm5uZg0aZLJx9q2bRskEoluuNlNmzbh//7v/3DkyBF4eHjo9qtduzYAIDc3F6+99hrq1KkDhUKB8+fPY+nSpTh06BD+/fdfyGQyvfQ3b96MyMhIKJVKJCUl4ffff8fKlSuxevVq7N69Gz179tTtO2vWLHz66afYunUrxo0bZ/JnIYQQQgixNxRYVDOnT5/GpEmT0KtXL+zbtw9yuVy3rVevXpg5cyaOHDlituNt3boVqampmDBhQqltnp6eaN++ve51z549ERYWho8++qhSgUXTpk3h7u6uCyy0n6NVq1YG5yPZuXOn3uvu3bvDzc0NkydPxu+//47u3bvrbW/cuDFat26tez1s2DBMnz4dnTp1wtChQ3Hr1i0EBAQAADw8PDBq1CisWLECY8eOpc7UhBBCCHF41BTKBIwxqHPVNrMwxkz+DMuWLQPHcdi4caNeUKElk8kwePBg3evdu3ejd+/eCAoKgpOTExo0aIC5c+ciNzfXqOOtX78egwYNMmoGaE9PT9SvXx/37983+vOYm5+fHwDjhxquWbMmPvzwQ2RnZ+OLL77Q2zZ69GjcvHkTMTExZs8nIYQQQoitoRoLE2jyNDjlesra2dDpnNMZYhex0fur1WocP34crVq1QmhoqFHvuXXrFvr3749p06bBxcUF169fx8qVK/HXX3/h+PHj5b734cOHuHTpktG1D0qlEvfv39dd3BdnaD4Gc03+plKpUFhYiP/++w/vvvsuOnXqhGeeecbo9/fv3x9isRi//fab3vpWrVrB1dUVBw8eLFX7QQghhBDiaCiwqEZSUlKQl5eHiIgIo9+zYMEC3XPGGJ555hk0aNAAXbt2xcWLF9G0adMy3/vnn38CAFq2bGlwO2NMFzA8fPgQ0dHRSEpKwuzZs/X2u3LlCqRSqdF5NsWZM2fQoUMH3ev+/ftj165dEIuND9hcXFzg6+uLhIQEvfVisRjNmjXDH3/8Ybb8EkIIIYTYKgosTCByFqFzTmdrZ0NH5Cx8S7a7d+9iwYIFOH78OJKSkvSaX127dq3cwEJ7oe3v729w+6FDh/QCBicnJ7z55ptYsmSJ3n61a9fGrl27Sr0/MzNTr8N0ZTRp0gTnzp1DXl4e/vvvP6xYsQK9evXC8ePH4ezsbHQ6ZTVL8/f3x7lz56qUR0IIIYQQe0CBhQk4jjOp6ZGt8fX1hbOzM2JjY43aPycnB507d4ZCocCSJUtQr149ODs7Iy4uDkOHDkV+fn6579duVygUBrd36tQJa9asAcdxcHZ2Ru3atUuNxKR9f/FO01olh4WtDBcXF13aXbp0Qbt27dC+fXt88cUXmD59ulFp5ObmIjU1FU2aNCm1TaFQVPg9EUIIIYQ4AosFFhqNBgUFBSbdBSbmJRaL0aNHDxw+fBgPHz5ESEhIufsfP34cCQkJOHHiBLp27apbn5GRYdTxtCMxpaWlISgoqNR2Dw8PgwGDNbVu3RoikQg3b940+j0HDx6EWq1Gt27dSm1LS0szOCIVIYQQQoijEawtTUFBAbZs2YIRI0YgODgYMpkMbm5ucHZ2RuvWrTFnzhxcuHBBqMOTMsybNw+MMUycOBFFRUWltiuVSuzfvx8AdEOklhw9quToR2WJjIwEANy5c6cqWbaokydPQqPRoE6dOkbt/+DBA8yaNQseHh54/fXXS22/e/cuGjZsaO5sEkIIIYTYHLPXWOTn52PVqlVYt24dMjMzERkZiR49esDf3x8KhQJpaWm4e/cuvvzyS3z44Yfo2LEjVq1apdeBlginQ4cOWL9+PSZPnoxWrVph0qRJaNSoEZRKJf79919s3LgRjRs3xqBBg9CxY0d4eXnhjTfewMKFCyGVSrF9+3ajA8J27drByckJZ86c0RvC1hYcOHAAX375JQYPHoywsDDdBIFr165FnTp1DM67cfnyZahUKqhUKiQlJeHUqVPYvHkzxGIxfvjhh1KjWaWmpuLWrVt48803LfWxCCGEEEKsxuyBRd26deHi4oIFCxbg5Zdf1k0YVhJjDDExMdi8eTOioqLw6aefGryYI+Y3ceJEtG3bFmvWrMHKlSuRmJgIqVSKevXqYeTIkZg6dSoAwMfHBwcPHsTMmTMxatQouLi44Nlnn8Xu3bvLHOmpOJlMhuHDh+PHH3/EsmXLhP5YJqlTpw5kMhkWL16Mx48fAwDCw8Px6quvYu7cuXozdWtpZ9CWyWTw9PREgwYN8Pbbb2PChAkGh8j98ccfIZVK8fzzzwv7YQghhBBCbADHKjPLWjm++uorvPLKKyYN13n37l08ePDAYBt1IWRlZcHDwwMpKSnw8fEptb2goACxsbGIiIgos+MxMc758+fRpk0bnDlzBu3atRPsOBqNBllZWXozb1tb586dUbNmTWzfvt3ix3aEv2GlUolDhw6hf//+gg03TCyLytTxUJk6HipTx5OamgpfX19kZmbC3d1d0GOZ/Qps/PjxJgUVAFCrVi2LBRXEslq3bo3nn38eixcvtnZWLOq3337DuXPnqt3nJoQQQkj1ZRu3dolD+/DDD9GmTRtkZ2dbOysWk5qaiq+//hq1atWydlYIIYQQQizC7H0svv76a5P2HzNmjLmzQGxMSEgIFi5caO1sWNSQIUOsnQVCCCGEEIsye2AxduxYo/flOI4CC0IIIYQQQhyA2QMLY2d1JoQQQgghhDgOswcWrq6uBkdaIoQQQgghhDgus3feDgwMRM+ePbFhwwYkJiaaO3lCCCGEEEKIDTJ7YLFnzx6EhIRg/vz5CAkJQefOnbFu3To8ePDA3IcihBBCCCGE2AizBxaDBw/Gli1bkJSUhAMHDqBBgwZYunQpIiIi0LZtW6xatQq3b98292EJIYQQQgghViTYPBYSiQR9+/bFxo0bkZiYiKNHj6JNmzZYu3Yt6tevj2bNmmHx4sW4cuWKUFkghBBCCCGEWIhFJsgTiUTo3r07PvvsM8THx+PkyZPo3r07Nm3ahKZNm1oiC6SEixcvYty4cYiIiIBCoYCrqytatmyJVatWIS0tzazH+vrrr+Hn51epCfLCw8MxcOBAs+bHHpw4cQIcx+HEiRO6daNHj8Zzzz1ntTwRQgghhJTH7KNCVYTjOHTq1AmdOnXCmjVrcO7cOUtnwWqYmiHjVAaKHhVBFiSDZ2dPcGLO4vn48ssvMXnyZNSvXx+zZ89Gw4YNoVQqcf78eWzYsAGnT5/GDz/8YJZj5eXlYf78+Xj77bfh5uZmljSrq+joaERGRuL48ePo3r27tbNDCCGEEKLH4oFFSW3atLF2FiwieW8ybr91G4UPC3Xr5CFy1FlXB35D/SyWj9OnT2PSpEno1asX9u3bB7lcrtvWq1cvzJw5E0eOHDHb8bZu3YrU1FRMmDDBbGlag1KpBMdxkEisd8rUrl0bffv2xYoVKyiwIIQQQojNMXtTKJFIBLFYbNRizYs0S0rem4wrw6/oBRUAUBhfiCvDryB5b7LF8rJs2TJwHIeNGzfqBRVaMpkMgwcP1r3evXs3evfujaCgIDg5OaFBgwaYO3cucnNzjTre+vXrMWjQIHh6euqtLygowLx58xAREQGZTIYaNWpgypQpyMjIMJjODz/8gKZNm0KhUKBWrVr4+OOP9bZrNBqsXr0aDRo0gJOTEzw9PdG0aVOsW7dOb79bt25h5MiR8Pf3h1wuR4MGDfDZZ5/p7aNthvTNN99g5syZqFGjBuRyOa5cuQKO4/B///d/pfJ3+PBhcByHn376yaRjAcD169fRt29fODs7w9fXF2+88UaZzcZGjx6NX3/9FXfu3DG4nRBCCCHEWsx+Zf/ee++B4yzfvMdWMTXD7bduA8zQRgAccHvabfg+6yt4syi1Wo3jx4+jVatWCA0NNeo9t27dQv/+/TFt2jS4uLjg+vXrWLlyJf766y8cP3683Pc+fPgQly5dwqRJk/TWM8bw3HPP4dixY5g3bx46d+6MixcvYuHChTh9+jROnz6tF/T8999/mDZtGqKjoxEYGIjt27fjrbfeQlFREWbNmgUA+OCDD7BixQq888476Nq1K5RKJa5fv64XqFy9ehUdO3ZEzZo18eGHHyIwMBA///wz/ve//yElJQULFy7Uy+e8efPQoUMHbNiwASKRCKGhoWjRogU2b96MV199VW/fLVu2wN/fH/379zfpWI8fP0bXrl0hlUrx+eefIyAgANu3b8fUqVMNfqfdunUDYwyHDh3Cm2++We73TwghhBBiUawayszMZABYSkqKwe35+fns6tWrLD8/v8rHSotJYzGIqXBJi0mr8rEqkpiYyACwF198sVLv12g0TKlUspMnTzIA7MKFC+Xuv3v3bgaAnTlzRm/9kSNHGAC2atUqg/tv3LhRty4sLIxxHMf+++8/vX179erF3N3dWW5uLmOMsQEDBrAmTZowtVpdZn769OnDQkJCWGZmpt76qVOnMoVCwdLS+DKIiYlhAFiXLl1KpfHxxx8zAOzGjRu6dWlpaUwul7OZM2eafKy33367zM8HgMXExJTKQ40aNdgLL7xQ5udkzLx/w9ZSVFTE9u3bx4qKiqydFWImVKaOh8rU8VCZOp6UlBQGoNQ1iRAsMipUdVb0qMis+1na3bt3MXLkSAQGBkIsFkMqlaJr164AgGvXrpX73oSEBACAv7+/3nptTcfYsWP11o8YMQIuLi44duyY3vpGjRqhWbNmeutGjhyJrKws/PPPPwCAtm3b4vLly5gyZQp+/vlnZGVl6e1fUFCAY8eOYciQIXB2doZKpdIt/fv3R0FBAc6cOaP3nmHDhpX6TC+//DLkcjm2bNmiW7dz504UFhZi3LhxJh8rJiamzM9XFn9/f8THx5e5nRBCCCHEGgQPLO7cuYPRo0cjODgYcrkcNWrUwCuvvFJt2ojLgmRm3a8qfH194ezsjNjYWKP2z8nJQefOnXH27FksWbIEJ06cwLlz57B3714AQH5+frnv125XKBR661NTUyGRSODnp99pneM4BAYGIjU1VW99YGBgqbS167T7zp07F4sXL8bZs2fRr18/+Pj4oEePHjh//rxuP5VKhU8++QRSqVRv0TZfSklJ0TtGUFBQqeN6e3tj8ODB+Prrr6FWqwHwzaDatm2LRo0amXys1NTUcj+fIQqFosLvnhBCCCHE0gTtPX39+nV06NABBQUF6N69O4KDg5GQkIBvv/0WBw4cwB9//IHIyEghs2B1np09IQ+RozC+0HA/C44fHcqzs6fgeRGLxejRowcOHz6Mhw8fIiQkpNz9jx8/joSEBJw4cUJXSwGgzA7WJfn6+gIA0tLS9C7SfXx8oFKpkJycrBdcMMaQmJhYaqSwxMTEUmlr1/n4+ADgJ2ScMmUK5s2bh6ysLPz666+YP38++vTpg7i4OHh5eUEsFmP06NGYMmWKwfxGRETovS6rr9C4cePw3Xff4ejRo6hZsybOnTuH9evX67abciwfH59yP58haWlpCA8PL3M7IYTYK6ZhUOeqoc5RQ52rhiZXA00Bv6jz1fzz/KfrNIUasEIGTeGT50XFnisZNEX8uuLPNUoNmIrxi5I9fa5iELuK0eDrBnBp5GLtr4IQuyRoYDF//nz4+PjgxIkTehexDx8+RPfu3fHOO+9gz549QmbB6jgxhzrr6uDK8CsAB/3g4sl1a521dSw2n8W8efNw6NAhTJw4ET/++CNkMv2aEqVSiSNHjmDQoEG6C+uSo0d98cUXRh1LGzTeuXNHdzcfAHr06IFVq1Zh27ZtmD59um79nj17kJubix49euilc+XKFVy4cEGvudCOHTvg5uaGli1bljqup6cnhg8fjvj4eEybNg337t1Dw4YNERUVhX///RdNmzYt9blN0bt3b9SoUQObN29GzZo1oVAo8NJLL+m2Ozs7G32sqKgorFq1yuDnM0SlUiEuLk5X80EIIbZKlaNC4f1CFNwrQGF8IVTpKijTlE8f01RQpaugylDpgglNvsba2Ubc6jhEbnbsm56ECEXQwOLkyZP4+OOPS90ZDwkJwXvvvYf//e9/Qh7eZvgN9UOj7xsZnsdirWXnsejQoQPWr1+PyZMno1WrVpg0aRIaNWoEpVKJf//9Fxs3bkTjxo0xaNAgdOzYEV5eXnjjjTewcOFCSKVSbN++HRcuXDDqWO3atYOTkxPOnDmjN4Rtr1690KdPH7z99tvIysrCM888oxsVqkWLFhg9erReOsHBwRg8eDCio6MRFBSEbdu24ejRo1i5ciWcnZ0BAIMHD0bdunXRsWNHBAQE4P79+1i7di3CwsJQt25dAMC6devQqVMndO7cGZMmTUJ4eDiys7Nx+/Zt7N+/v8JRrrTEYjHGjBmDjz76CO7u7hg6dCg8PDz09jH2WNOmTcNXX32FAQMGYMmSJbpRoa5fv27w2BcvXkReXh6ioqKMyishhAiNy+CQvCsZ+RfzUXCvQLcoU5RVSBQQu4ohchZB7CSGSCGCyEmk/6gQQSTnF07OQSQz8FzKgZPxrzkZB5H0yTopB06ivxTcLcCNCTeQvCcZdT+rC7Gz2HxfEiHVhKCBRV5enq6pSkm+vr7Vqp2431A/+D7raxMzb0+cOBFt27bFmjVrsHLlSiQmJkIqlaJevXoYOXKkbqhTHx8fHDx4EDNnzsSoUaPg4uKCZ599Frt37zZYU1CSTCbD8OHD8eOPP2LZsmW69RzHYd++fYiOjsbmzZuxdOlS+Pr6YvTo0Vi2bFmpGpLmzZtj3LhxWLhwIW7duoXg4GB89NFHerUdUVFR+Pbbb7Ft2zZkZWUhMDAQvXr1wrvvvgupVAoAaNiwIf755x8sXrwYCxYsQFJSEjw9PVG3bl2TawDGjRuH5cuXIzk5WddpuzhjjxUYGIiTJ0/irbfewqRJk+Ds7IwhQ4bg008/xbPPPlsq3X379sHX1xe9e/c2Kb+EEGIuqhwVMn/LRPqv6Ug7mgb3y+64iZsG95V4SqAIV0AeIofUVwqJlwQSbwmkXlL+0VsKsYcYYtcniwv/KFKILD50PevGcH/JfRTcK0Dq/lT4v+Bf8ZsIIXo4xpihlv9m0bx5czRp0gTffPNNqW2vvPIKLl68iH///Veow5cpKysLHh4eSElJMRj4FBQUIDY2FhEREaU6HhPTnD9/Hm3atMGZM2fQrl07wY6j0WiQlZUFd3d3iESOOdiZWq1GnTp1MHLkSCxdurTcfR3hb1ipVOLQoUPo37+/Ljgk9o3K1H5pVBokbk7E428eI+t0FphK/9LBpZkLPLt4wqmOExThCn4JU0DiYV8T4ca+G4v7S+7De4A3mh5oau3sWAWdp44nNTUVvr6+yMzMhLu7u6DHEvSM/9///ocJEyYgMzMTr7zyCoKCgvDo0SNs27YNP/30EzZt2iTk4YkNaN26NZ5//nksXrwYBw4csHZ27Nq2bduQk5OD2bNnWzsrhJBqJO1oGu7MuIPcy7m6dYpwBbx6esGtmxvOq8/jmZeecYiLUP+X/XF/yX2kHUlDUXIRZH7Cj9hIiCMRNLAYP348Hj9+jCVLluDgwYMA+JF/nJycsHTpUoNNSIjj+fDDD/F///d/yM7Ohpubm7WzY7c0Gg22b98OT09Pa2eFEFIN5F7PxZ1Zd5B2MA0AIPGSoOa8mvAb5genWk4A+Lvb7JBgDR8sziXSBW6t3ZB9PhtJu5MQMrX80RMJIfoEr6OcN28eJk+ejNOnTyM1NRU+Pj7o0KFDqc6uxHGFhIRg4cKF1s6G3aNAnBBiCcpUJe4tuoeE9QlgKgZOwiF4SjDC3wuH1Nv+ayUqEjAqANnns/F422MKLAgxkUUaP3p4eKBv376WOBQhhBBCKunRlke4M+MOVOkqAIDPIB/U/qA2nOs7WzlnluP/oj9uz7yN7LPZyLuVB+e61eezE1JVFgksrly5gvv376OgoKDUtqFDh1oiC4QQQggpx4OVD3B37l0AgEsTF9T+qDa8e3pbOVeWJwuQwbuXN9KOpOHx9seIiI6o+E2EEAACBxZ37tzB8OHDcfHiRQB8/4riOI6DWq02Or3s7GwsXrwY//33H/7991+kpKRg4cKFiI6ONme2dQQcMIsQQdHfLiHEWIwxxL4biwdLHwAAas6viYj3I6wyHLqtCBgdwAcW2x4jfGG4xYe+JcReCRpYvPbaa0hMTMSaNWvQoEGDKs12DPDDZW3cuBHNmjXDc889J9ioUtqRLfLy8uDk5CTIMQgRUl5eHgA4xCgthBDhMA3D7em3Ef9xPACg1spaqDmnppVzZX2+z/pC5CJCwZ0CZJ3Ngkd76hdKiDEEDSz++usvfPnll3jxxRfNkl5YWBjS09PBcRxSUlIECyzEYjE8PT2RlJQEAHB2dqa7FTZOo9GgqKgIBQUFDjuPhTEYY8jLy9NNyCcW08yxhBDDmJrhxms3kPhVIgCg7md1UWNyDSvnyjaIXcTwG+qHx988xuNvHlNgQYiRBA0s/Pz8zDr6kyUv7gMDAwFAF1wQ28YYQ35+PpycnCgIBODp6an7GyaEkJI0Sg2ujb6G5N3JgAiI/CoSga/Qb0ZxAaMC8Pibx0janYQ6a+pAJKu+N60IMZaggcWkSZPw5Zdfol+/fkIeRhAcxyEoKAj+/v5QKpXWzg6pgFKpxG+//YYuXbpU++Y/UqmUaioIIWVSF6hx9fmrSN2fCk7KocGOBvAf7m/tbNkcz+6ekAXKUJRYhLSf0+A7yNfaWSLE5gkaWMyePRszZ85Eq1at0K9fP3h7648uwXEcpk+fLmQWAACFhYUoLCzUvc7KygLAX4waEzTQRZrt02g0UKlUEIvF1b68NBoNNBqNtbNRZdpzkwJ7x0Flan2aAg2uDb2GjF8zIFKIUH93fXj186p0mTh6mfq+4IuEdQl49PUjePStHs2hHL1MqyNLliXHBBw+5uzZsxgwYADS0tIMH9zEUaGKS0lJgZ+fn1GjQkVHR2PRokWl1u/YsQPOzjQ+NSGEkGpACTivdIb0vBRMwZD7Ti7UTSr3P7i6EN0RwW2mG5iMIWtzFuBi7RwRYrq8vDyMHDkSmZmZcHd3F/RYgtZYTJ06Fb6+vvjqq6/MMipUZc2bNw8zZszQvc7KykJoaCiioqLg4+NjlTwR81IqlTh69Ch69epV7ZtCOQoqU8dDZWo9TMVwY9QNpJ5PhUghQsP9DeHRtep34B29TBlj+PfLf5F/PR+t8lohYESAtbMkOEcv0+ooNTXVYscSNLC4cuUKdu3ahcGDBwt5mArJ5XLI5fJS66VSKZ00DobK1PFQmToeKlPLYhqG669eR+pevk9Fox8awaeneW+qOVKZMjVDxqkMFD0qgixIhoBRAbi34B5SdqYgZEKIWdP27Oxps/OFOFKZVneWLEdBA4uaNWvSRF2EEEKIlTDGcHPyTTz+5jEgBhp+2xA+fR2jpl6Ii/Tkvcm4/dZtFD582i9TFsS3tsiIyUDBwwIoQhRmS1seIkeddXXgN9SvSvm2p4CFODZBA4u5c+di9erV6NOnDxSKyp2IJR0+fBi5ubnIzs4GAFy9ehXff/89AKB///7UZ4IQQggBH1TcmXEHj754BHBAg28awO+5ql3AViofFgoAqnqRnrw3GVeGXwFK3A8tSizinzAgeXcyQmeGmi3twvhCXBl+BY2+b1SlfAsVsBBiKkEDi3/++Qfx8fGoXbs2oqKiDI4KtW7dOpPSnDRpEu7fv697/d133+G7774DAMTGxiI8PLzK+SaEEELsXey7sXi49iEAoP6m+gh4qfz+AfYeAFTlIp2pGW6/dbtUmvzGp09TDqWYHFhUmDYH3J52G77P+pr8fQsZsBBSGYIGFp9++qnu+Y4dO0ptr0xgce/evapmixBCCHFYjDHcX3ofD5Y+AADU/bQugsYHlfsehwkAKnmRnnEqQ++zlyXzt0yoclSQuBp/+VRh2gwojCtExqkMeHXzMjpdIQMWQipL0GkktePpl7VUdqhZQgghhJTGNAy3p9/GvXfvAQBqraqFGlNqlPsebQBQ8uJXGwAk7002PR9G1ADcnnYbTG1aP0xTLtJNUfSoyLgdVUDGcWHSNjoPTwj1XegloWZIP5GOxzsfI/1EusnlRaofQWssCCGEEGIZmkINro25huRv+UCg9uraFTbbsVoNQCXv0gt1ka7toG2MtCNp8B1s/CzcxqZtSh4A4b4LrdQfUhE7M5b6bhCTmL3GIjc316LvI4QQQqo7VaYKF/teRPK3yeCkHBpsb2BUXwBr1wAIFQCYepHu2dkT8hA5UFbsxAFSX37IzrTDaSaNeGlM2vJQOTw7e5qUZ6G+CwCQnJbg+ovXzVqLRaoHswcWERERWLNmDbKysoza/9y5cxg8eDA++ugjc2eFEEIIcXiF8YX4t/O/yDiRAbGbGE0PN0XASOMmcnPEAKAyF+mcmEOddXV0aZRMEwDqrKsDTsah4F4B8m7kmTfttXVM7gch1HfB1AxOm5zM3oyNVA9mDyxWr16Njz76CEFBQXjxxRfxxRdf4Ny5c7h//z4eP36Ma9eu4cCBA1iwYAGaNGmC9u3bw8XFBePHjzd3VgghhBCHlnstF/90/Ae5l3IhC5Sh+W/N4dXD+KZFDhkAVOIiHQD8hvqh0feNIK+hP6GuPESORt83QsDIAHh25fObdjjNrGlXpmmRUN9F1u9ZEKWWc3lohr4bxHGZvY/FmDFjMGLECGzZsgUbNmzAt99+C47T/6NmjMHJyQnDhw/Hli1b0KpVK3NngxBCCHE4xYeEVaYqEftuLNQZajjVc0LTI03hFOFkUnraAKAwvtDwHWqOv/itbABwZfgV/iK3eNpmCgAMjmK1tmrt//2G+sH3Wd8yh9317uuN9KPpSDuShtDppg07W1Halc2vub8LoftuEMcmSOdtJycnTJo0CZMmTUJ8fDz+/PNPJCQkID8/H76+voiMjES7du1oqnhCCCHESIaGhAUARV0FWvzRAjJf09vSO2oAUBWcmCuzQ7l3P2/cmXkHGSczoM5TQ+wsNlvalWXu70LIvhvE8Qk+KlSNGjUwYsQIoQ9DCCGEOKyy5oQAgILbBcj8LbPSF+qOGAAIxTnSGfIwOQrvFyIjJgM+A3wsevyymPO7cO/kDo2PBqI0kVlrsUj1QMPNEkIIITas3CFhn6jqRGiOFgAIheM4+PTzQcKGBKQdSbOZwMKcODGH/An5cFnlYvZaLOL4BJ0gjxBCCCFVY4mJ0ICnAUDASwHw6uZFF45l8O7rDQBIPZxq5ZwIR9VBhchdkWbtbE6qB6qxIHalML4Q6cfTIXGXwKuPF8QKw+1bi3dwNOedN0IIsTTqTGtbPLt7gpNyKLhTgLxbeXCu62ztLAnCZ4gPAoYF0P9SYhIKLIhNUxeokfl7JtKOpCH953TkXn46kaLYQwy/oX7wH+kP106uuvWGOjjSbKGEEHtFnWlti8RNAo/OHsg4noG0I2kOG1gAjtWMjVgGBRbEJiV9m4TErYnIiMmAJl/zdAMHuLV2Q9GjIhQ+LETi5kQkbk6ENFAKRWsFHl55iPsL7pdqi6ydLZSqcAkh9sazsydkITIUPSyjRoI601qcd19vPrA4nIaQN0OsnR1CbAYFFsTmpPyUgqsvXNW9lgXJ4N3XG959vOHV0wtSHymYhiHz90w83vEYyd8lQ5mohPyAHPcP3DecKAPAVb2DIzWxIoRYGifmEPBiAOJWxxnYyD9QZ1rL8u7njbtz7iIjJgPqfDXETqYNO0uIo7JIYHH58mVcu3YN+fn5pbaNGTPGElkgdkKVpcLNyTcBAAGjAxA6OxQujV1KTbLIiTh4dvGEZxdP1P24LpIPJ+PiexchvVjO3CjFOjhWpmqXmlgRQqxBla3C452PAQBidzHUWWrdNnMMCUtM59LIhZ9Y8GEhMn/LhHcfb2tnya7QTTrHJWhgkZeXh8GDB+P48ePgOA6M8e1Til8kUmBBirs77y6K4ovgVMcJ9b6oZ9RdIJFMBO/+3lD+piw/sHiiMh0cyxpDnppYEUKEdm/hPRTFF0FRW4HWF1oj+1w2XZBZGcdx8O7rjUebHiH1cCoFFiagm3SOTdDhZhcvXox79+7h5MmTYIxh7969OHr0KIYOHYq6devin3/+EfLwxM5k/pGJhM8TAAD1NhoXVBSn8dJUvBNM7+BY7hjyT9bdnnYbTF3OIPOEEFIJ2f9l4+G6hwCAep/Vg8RFQkPC2gjvfnwwkXY4zco5sR/am3Qlh0/W3qRL3ptspZwRcxE0sPjxxx/x9ttvo2PHjgCAmjVrokePHvjuu+/QsmVLrF+/XsjDEzuiKdTgxoQbAIDAVwPhFWV6UyV1QzVkNWS6NsdlST2YClWOyuh0LTWGPCGEFMcYw60ptwAN4Pe8H90VtzFePb3ASTjk38xH/t3STb2JPrpJVz0IGljcu3cPkZGREIvF4DgOeXl5um0vv/wy9u3bJ+ThiR25v/w+8q7nQRogRe0PalcuETFQ66Na/POSwUWx1w9XP8S5hueQvDdZ1zyvPDSGPCHEGrLOZCHrzyyIFCLUWVPH2tkhJUjcJXB/xh0AkHaEai0qQjfpqgdBAwtPT0/k5vLzDvj7++PWrVu6bUqlUreNVG+5V3LxYNkDAEDdT+pC6lVxP4my+AzxQaPvGxmeLXRPIzT+qTEU4QoUxhXiyrAruNT/EvJu55WRGo/GkCeEWEPCer5pqP+L/pAHyyvYm1iDdhZuag5VMbpJVz0I2nm7SZMmuHnzJvr27YuoqCgsW7YMdevWhUwmw/vvv49mzZoJeXhiB5ia4caEG2BKBp/BPvAbXvWOW35D/eD7rG+ZI0549fDCg+UP8GDVA6QdScO5xudQc25N1Hy7psF+HZ6dPfnRP+ILDVfh0hjyhBAzU6YqkfRtEgAgeFKwlXNDyuLTzwex82KRfjwdmkINRHJB79faNbpJVz0Iega8+uqryM7OBgAsXboUeXl56Nq1Kzp06ID79+/jww8/FPLwxA7Er49H1pksiN3EqPtZ3VLDylaWdrZQQx0cxc5iRCyOQJtLbeDV2wuskOH+ovs41/gcUvanlGoexYk51Fn3pBlCGU2saAx5Qog5Pdr8CKyQwbWlK9zauFk7O6QMLk1dIAuSQZOnQfqxdGtnx6Zpb9KV2Q+SA+ShdJPO3gkaWDz//PN45513AAARERG4efMm9u3bhx9//BG3bt1Cu3bthDw8sXEFcQWInRcLAKi1ohYUIQqLHt+5njOaHmmKht82hKyGDAV3C3B58GVcGngJebf0m0f5DfUru4kVDTVLCDEjpmFI2MA3gwqeFGy2Gy7E/DiOg89AHwDAlRFXEP9ZPJiGOh8bQjfpqgeLzrzt4uKCQYMGWfKQxEYxxnBr8i2oc9Rw7+iO4DesU9XPcRz8R/jDu583Hix9gLgP45B2KA3nfj2H0JmhqDm/JiSu/GlSURMrQggxh/Rf01FwpwBiDzECXgqwdnZIBSKWRiD/Tj4yjmfg1tRbSNmXgvpf1Yci1LI3y+yB9iadwXksaKJHh2D2Gov09HQMGzYMBw4cKHOfAwcOYNiwYUhNTTX34YmdSD2YitQDqeCkHOp/WR+cyLoX5xJXCWotr4U2l9vAu683WBHDg+UP8FfkX0janfR0csdymlgRQog5aDttB44JhNjFtPl8iOXJ/GRodrQZ6nxcByInEdJ/Tce5JueQ+HWiUSMPVjd+Q/3Q/l57NItphgY7GqBZTDO0j21PQYWDMHtgsWnTJly4cAF9+/Ytc5++ffvi0qVL+Oyzz8x9eGIHNEoN7sy6AwAImR4Cl4YuVs7RU871nNHkUBM0/rExFBEKFMUX4eqLV/Fft/+QeTrT2tkjhDi4gocFSPkpBQCsVpNLTMeJOIS8GYLW/7WGWzs3qDPVuP7KdVwZegVFSTTKUUlC3qRjaob0E+l4vPMx0k+k07wYFmb2wGLXrl2YOHEiJJKyW1lJJBJMnDgRP/30k7kPTwRg7pP00cZHyL+RD6mvFGHzw8yUS/PhOA6+g33R5mobhL8fDpGTCJm/ZeLfjv/i0uBLyLmUY+0sEkIc1KMvHwEawKOrh03ddCHGca7njBa/t0DEsghwUg4p+1JwrtE53Ft0D/l3aBI9oSXvTcaZ8DO4EHUB10Zew4WoCzgTfoZm9LYgswcWN2/eROvWrSvcr2XLlrh586a5D0/MzNwnqSpThXvR9wAA4YvCIfGwaDcfk4gVYoS/G462N9oiaEIQIAZS96fifLPzuPry1QrnvyCEEFNolBo+sADVVtgzkUSEsHlhaHWuFVyaukCZosS96Hs4W+cs/un4D+I/j0dRCtVimFvy3mRcGX6l1CR8hfGFuDL8CgUXFmL2wEKlUkEqrXiCM6lUCqVSae7DEzMS4iS9v+w+lClKOEc6I+i1ILuoslSEKlD/y/poe6Ut/J73AxiQtCMJ5xqcw403bvDzWxBCSBWl/pSKokdFkPpLqb25A3Bt5opW51oh8ptIePXxAkRA1uks3JpyC6eDTuPS4EtI2p0EZTpdC1UVUzPcfuu24bmmnqy7Pe22TV5jOBqzBxZBQUG4evVqhftduXIFgYGB5j48MRMhTtL82Hw8XPsQAFB7dW2k/pRqV1WWzvWd0Wh3I7T6pxW8+3uDqRgeffEIZ2qdwfXx15F7hWaSJ4RUXvz6eABA0KtBEMloojVHIJKJEDgqEM2ONEOH+A6ovaY2XFu5gqkYUven4uqLV/GHzx841/wcbv3vFpK+T6I+GZWQcSqj1E1QPQwojCtExqkMi+WpujL7L1fXrl3x+eefl1sboVQqsX79ekRFRZn78NWaOe/+C3GS3p13F6yIwbOHJ9QFarutsnRr4YamB5ui+W/N4dHFA6yIIXFzIs41PoeL/S8i/Xi6ySOB2EPNDSHEeKae03k385BxLAPggKDXgiyTSWJR8kA5QqeFovX51mhztQ1qvlMTTvWdAAbkXshF/CfxuDriKv4M+BN/NfgLN167gYSNCcj6KwvqPLW1s2/Tih4ZF4wZux+pPLM3cJ8+fTpat26NIUOGYOPGjQgO1m8nmpCQgIkTJ+LGjRvYvn27uQ9fbSXvTTY8LvS6yo0Lbe6TNPNMJpJ3JwMcUHtVbVx+9nLZtSEcXxvi+6yvTQ/n6tnZEy1OtkDm6UzEfRiHlL0pSDuchrTDaXBt4YrQmaHwe94PImn58bu5y44QYl2VOae1E+J59/eGU7iTRfJJrMelgQtqLamFWktqoehxETJOZSDzt0xknMxA7qVc5F3PQ971PF2fG4j4WnPX5q780swVzg2dIQ+R0wSKAGRBMrPuRyrP7IFF06ZN8dlnn2Hy5MmIiIhAq1atEBERAQCIjY3F33//DY1Gg/Xr16NJkybmPny1pO0LUfJCXXv3vzIzQ5vzJGWM4c4MfnjZwHGBUGWpjK4N8ermZVQ+rMmjgwc8vvdA/p18xK2JQ+JXicj5NwfXRl3DnTl3EPhKIALHBcK5rnOp9wpRdoQQ66nMOa3OVyNxSyIAoMakGpbKKrERsgAZ/If7w3+4PwBAmaZE5h+ZyPwjEzn/5SDn3xwok5TIu5aHvGt5SNqZpHuvyEUE53rOcI58stTnHxW1FbrJXasDz86ekIfI+T6Phm5acnxw79nZ09JZq3YE+aubOHEiGjdujGXLliEmJgZnzpwBADg7O6Nv376YN28e2rdvL8Sh7QJTM7PN3lxhX4hK3v0350ma/F0ysk5nQeQsQsTiCGSczDAqD/ZWZelU2wn1Pq2HiEURiF8fj/hP4lGUUIQHyx/gwfIH8OjsgcDxgfAb7geJq0SwsiOEWEdlz+nkb5OhSldBHiaHd19vi+WX2CaptxS+g3zhO8hXt64wsZAPMp4suRdykX87H5pcDXL+5YOPUun4S+FUywmKWgrdoyJCAUW4AvIa8gpr0+0JJ+ZQZ10dPqjnoH8OPjnV6qytQ/9LLUCwcLZDhw7Yv38/NBoNUlL4yX58fX0hEjnOH3JlmLvZiyl9IUy5+2+uk1RTqMHduXcBADXn1IQ8WO7wVZZSHynCF4Sj5uyaSPkpBYmbE5H2cxoyT2Ui81Qmbk29Bf8X/OHS1MWham7MRRt458flQ3xfDNaHARUPNEeI1VXm91hTqMHDj/lBLYJfD6YLH2KQPFAOeV85fPr66NZplBoU3C3gm03dyNM1n8q7kQdVmgrKJCWUSUpknckqnSAHyIJlUNRUQF5TrnuUh8ghDhCDS+P4fkF29NvrN9QPjb5vZPgaa63tNi02581mWyB4PZlIJIK/v7/QhxGEuQtbiGYvQnZYMsdJ+vCThyiILYAsWIbQWaEAqk+VpUgugv8If/iP8EfBwwI8/voxEjcnIv92PhK/SjQ6narW3NjTj1bJwNsVrji/4TzqflzXZv8pEKJl6u+xOk+NK8OuIOefHIicRQgaT522ifFEUhHf9Kl+6Wa2ygwlCmILUHC3APl381EQ++TxbgEKHhSAFTIUxRehKL4IOF06bXe448+Jf0IeJIeshgzyGvxNQXmQHLJAGWRBsqeP/jKb+Z/iN9QPvs/62u3/PMB8fSyL/+/Pc7bcvFvVpwGeicxd2EI1exH67n9VTtL0mHTcX3IfABCxNAJiFzGA6lllqQhRIGx+GGrOq4nM3zOR+FUikr5LgiZXU+F7q1JzY6kfLSED76KEIupvQuyCKb/HqmwVLg26hMyTmRA5i9D4x8aQBdhnLS2xPVJPKaQtpHBr4VZqG2MMyiQlCh4UoPBB4dPH+wUojC/kl0eF4NQcCh8WovBhIbKRXfbBRIDUVwqZvwzSAP5RFiCD1F/KP/pJ+cVXCpmfDGJ3saAdzjkxJ0gtvz3cbC6edvH//bmw3HD4FFgYIERhC9VkyRJ3/009SdW5atx5+w4SPuNHOXFr54bA0fpzlthrlWVVcRwHz86e8Ozsibqf18WZmmegTCl7aGaJpwSyQBkYYyb/EFvyRwuwzcCbEEsy9vfYpZELLvS8gOy/siF2F6PJwSbw7ORp6eySaorjOMgC+It/tCm9XalU4tD+Q+jRqgc0SRo+0HhYiKLEIhQ9Knr6+KiIn3NDA12zK1w24vgyDlLfJ8GGz9NF4iPh12tfe0sg9eYfJZ4SiCTWa0pvT//zyvrfbykUWJQgVGEL1WTJ1u7+Z5zKwPWx11FwtwAAEPxGMGqtqmXw+PZWZWluYicx6n1Rjy87wODfnCpDhXMNzkERroBXby949/aGZ3dPSL3Kb/hqjR8tWwy89ZKwoyZhxD4Z83scvigcF3peQO7FXEi8JWj2SzO4tSp9V5kQqxID8mA5pGFSg8GHFlMzFCUXQZmkRNHjp49FSUVQPlbyj8lKKJOVKEougiZXA1bEUJRQhKIE0653xB5iXaAh9ZJC4imBxOvJ4lns0VMCiUexRw8JRM6iSteS2NP/vHL/91sIBRYlCFXYQjZZsoW7/+o8NWLficXDdQ8BBshD5aj/f/Xh3av8EU6EqrK0F2WVnTRACq+eXihKKELm75kouFeARxsf4dHGRwAHuDR2gUcnD92iqKnQS9cqP1o2GHhr0VwhxFLK+z0OezcMD1Y+QP6NfMgCZWh6tClcG7taMbeEVA0n5viO5YFyo/ZX56uhTFHqgg1l6pMlhX9Upap0r1XpKijTlFBn8ZMDqjPVUGeqgdhK5FPCQewh5gMNdwn/3J0POsTu/HOxuxhityfPnzyKnEW4Ofmm3fzPq/B/vwVQYFGCUIUtdJMla9z9V+eroUzmx9a+9b9byL+ZDwAImhCE2qtrQ+JBf17GqKjsVDkqZP6WibRf0pD+czryruch91Iuci/lImE939xMHiqHxzMecH/GHW6t3VBwv8CoY5v9R8sGA29LzBUiVG0IpSt82kKka+iclteQ42Lviyi4VwB5qBzNjjUzOLcNIY5M7CSGOFQMRaii4p2f0Cg1UGWo+JGu0pT8YzofeKgyVHqPynQl1Jlq/nUmv0ADMBWDKlUFVarKvB/oyf+8v9v/DUWoAmJXMcQuYohcRLrnYucnr4s9V2Ualw9T/+fZwjD9dOVXTGFCod7EM+VJ3JKItCNpJqWvqK0o+6KMAYpaClwff92kNMuTuNn4kYe0edA+Msb0XoMBqkwVf4ch5WmVZnGyYBnqb6oPn34+IKYpr+ZG4iqBT38f+PTnv9eix0X85Em/80vOvzkojCtE0q4kJO168vdrZFNUoX60bCXwtkTfDaFqQyhd+80zYwz5sfkojCtE1uksZP6ZidxLuYAGcKrjhGa/NoMizPgLK0KqM5FUBJmfDDI/028sMcagzlHzQUaGCuosNVRZKj74yOIDD3UWv12drYY6+8n2J49FiUVGBSM553OQc770XCJVIgZuTr4JsbMYYicxRE6ip4tCxK9TiJ4uTiIUPDDupqKQKLAo5tabt5C6P9WofdN/STf78TNPZiITmWZPV0iclIPUTwqfAT6otbJWhW3/SdXJAmTwG+qnu/BR56qR9VcWMk9lIutsFrL/zobycdkdwrW0HeU0Ko3RneKEqlkQqq+Q0H03hKoNsUq6w64gYmkE3Du6gxUxMOWTnUQAJ+L4RzGn/1rEgRNzSD+Zjti3S7dP0Oa3zsd14N3HGyq1Clwih4LYAqhkKr7NMwfdomsDrS1mDkg9lIqbE0s3RdCmXf//6sNnYImbGcX+TMpqV52yPwU3xt8oM916G+vBu483P5a/mg9SmZqBqRg0BZqnFyjFLkTUWWrkXs1F1uksKJNLn4Nu7dzQ+IfGkAcZ12yEEFI1HMdB4iaBxE0ChJj+/vQT6bgQdaHC/ULnhUJRUwF1jhqaXA3UOWp+yVNDnftkXfHnuXytijpbXXaiaiD/Wr7pmbYyCiyeyDqXhZS9KQAH+L3gh+RdyWXuGzAmAC5NXCp9LKZhKIgtgDpLDbG7GIoIBf+P2gYY+seufZS4SZ4OG+dnmWHjSMXELmJ4RXnBK4q/MGaM7xj38JOHiFsZV+b7VKkqnG9yHpyMg1MdJ36owGLlqx0aUOon1XWQc23qClmIjB/73MxN+spqmy6rIUPddZWbx0LIvhtC1YaYO11VpopvPnctF7f/V066AGLfqUTj5fI8Sff2m7d1q9zhjr/xt9nSvjH+RtXTMpDuzYk3q5QMJ+Pg1soN7h3c4d7BHR4dPCCvQQEFIfbE2Nr0WosND1JTkeS9ybj11i0UPXz6P0gWKEPIrBB4dPSAJl8DTb4G6ny17rmmoNhS4rU6X42C2AJk/1XO8MACo8DiCe0/1IDRAWiwtQGSRxioIg917KFQiWPgOA7yGnLUXlEb7m3dS/0diz3FcGvlBnWOGrmXc6HJ1SDvah7yrho5gU7JGoXiGBD6diiYmlXqR7Z42/T8uHz8c/8f9J7VGzJF5cb3F7LvhlC1IeZIV5mhRMreFCTtTEL68XSg4ulSdORhcn5oR6kI4PgbIdDwAQ80/Gum5ptKMjWDOlcNZWLFNWQiZxEgAtQqNcQisa6Jpa7ZZfGmlyi23oS8m52Eb4bBiTlAzNfYcGIOIoWI7/TpJn7a8fNJZ095qBzuHdzh1tINIrn1hsckhFSd0CNvCtU/1lAzT0uhwAL8RG7pR9PBSTmELwoHQEOhEsdQ0d8x0zAUPChA/u18Xf+Z4kMDatep0vmOcZoCTYXD2N2eehu3/3cb8hpyyEPlkIc8fVSEKvjXIXJIA6QGm2Bp+5u4Kl2hPqSu0jkn5KAJQtWGVDZddZ4aqQdS8XjHY6QdTgMrevqBZTVkkHhIjAoeay2vhYCXAozO7+Odj3Ft5LUK96u/qT68h3vj0KFD6N+/P6TSiptNGpt2g+0N4P+Sf+kNZfytJu1MwrVRRqT7dQOTvgtCiOMReuRNIUbHLPm/P885D3jOrIcoU7UPLBhjuDvvLgAg+PVgOIU76bZV96FQiWMo7++YE3FwCnfS+7svj7pArQsylClKZJzIQN7NPL7NaJEahbGFKLhXAE2+BoVxhSiMK+duCQe+SV2QDPIgOWRBsqeLvwycFwfRfRGKkoogCZRUKsAQ8m6TULUhpqSrzlcj7ec0JH+fjNQfU6HOedpe16WxC/xf8of/i/5wquVkdFthIfNrKqPTDpYZbpJZRrHKagiXZ0KI47HHm83F//enphrXf9gcqn1gkbo/FdlnsyFyFqHmOzWtnR1CbJpYIYY4SKzrfOrZxbPUPowxFD0uQuF9frbWgrgCFD7kgwzdY0IhoH46W2vuhVyDx3ODG869dY4PQnyeztQq8Zboz9bq/eR58QmSvPjxyoW62yRUbYgx6Uq8JYj/LB6XBl7SG51NEaHgg4mX/EvNj2DN/GrTVWlMG+rRFvJMCCEA3Ww2VrUOLJia6fpWhLwVYvQEL4SQsnFcsQmT2hneh6kZlClKFD4qRNGjIt1S+KgQRYl8E6yix0XITciFKFsEMPDNtFIqbstfkthdDIknPyGSs4czODEHibcETrWc+OF6L+bw440bWlyejEnuLOLHH3cW6ZpvCVUbUm664F+rUlVI+T4FACCvKYffMD/4Pe8H93buZQ6mYJX8lkzXxP4SNpFnQgghRqvWgUXKnhTkXs6FxFOC0Nmh1s4OIdUGJ+YgC5BBFiADmhveR6lU4tChQ+jXux+QBRQlPenzkVZsdtaSzzOeTpSkyeOvYtVZat3MrcVlnqjc0M6cjHsaZChEkIXIoExUPh2eFYBIIYJbGzekHkpF+q/pEMlF4OQcRDIROBkHjuP4zsla2g7LKobCeL45mSxAhqJEw/0tnOo6wW+YH3yH+cKtlZvRI7MJVXsjZBtke8wzIYRUV3YVWOTk5GDBggX49ttvkZaWhsjISMydOxcvvvhipdJ7uPIhJJAgdE4ozb9AiI3iJBykAVI+CDGBpkijmxRJlV5svoGs0o+6Mcdziz0vNg65Ju9pp3VWxKAqUgEZ5Rw7X4PM3zKR+ZsZ5qURARIvCeQ15PB9zhd+I/zg0sil0sM8C9VWWMg2yPaYZ0IIqY7sKrAYOnQozp07hxUrVqBevXrYsWMHXnrpJWg0GowcOdLk9ArvF8IpwAkh/6vErCmEEJsmklV+ttaSGGPQFGqgyeMnNtLk8ZMd6cYQz386hrj2NSt68p7C0s91is8X82SCOFmwDIowBb+EKyCrITN6AkNjCdVWWMg2yPaYZ0IIqW7sJrA4dOgQjh49qgsmACAqKgr379/H7Nmz8cILL0AsFpucbvi74RC7mP4+Qkj1wXEc33FdIYbUm2o3CSGEEEPsZvaeH374Aa6urhgxYoTe+nHjxiEhIQFnz541OU15TTmCJgaZK4uEEEIIIYRUW3YTWFy+fBkNGjSARKJfydK0aVPddlOFzAmBSGY3XwEhhBBCCCE2y26aQqWmpqJWrVql1nt7e+u2l6WwsBCFhU9H/cjM5DtUinqILDppCBGOUqlEXl4eUlNTjZrRl9g+KlPHQ2XqeKhMHQ+VqeNJS0sDAP3RCAViN4EFgHJHQSlv2/Lly7Fo0aJS6+s3qG+WfBFCCCGEEGLLUlNT4eHhIegx7Caw8PHxMVi7oI3CtDUXhsybNw8zZszQvc7IyEBYWBgePHgg+BdMLCMrKwuhoaGIi4uDu7u7tbNDzIDK1PFQmToeKlPHQ2XqeDIzM1GzZs1yr5XNxW4CiyZNmmDnzp1QqVR6/SwuXboEAGjcuHGZ75XL5ZDLS8+q7eHhQSeNg3F3d6cydTBUpo6HytTxUJk6HipTxyMSCd+v2G56Lg8ZMgQ5OTnYs2eP3vqtW7ciODgY7dq1s1LOCCGEEEIIIXZTY9GvXz/06tULkyZNQlZWFurUqYOdO3fiyJEj2LZtW6XmsCCEEEIIIYSYh90EFgCwd+9evPPOO3jvvfeQlpaGyMhI7Ny5Ey+++KJJ6cjlcixcuNBg8yhin6hMHQ+VqeOhMnU8VKaOh8rU8ViyTDlmibGnCCGEEEIIIQ7NbvpYEEIIIYQQQmwXBRaEEEIIIYSQKqs2gUVOTg6mTZuG4OBgKBQKNG/eHLt27bJ2togRTpw4AY7jDC5nzpzR2/eff/5Bz5494erqCk9PTwwdOhR37961Us6JVnZ2NubMmYPevXvDz88PHMchOjra4L6mlOEnn3yCyMhIyOVyREREYNGiRVAqlQJ+EqJlbJmOHTvW4LkbGRlpMF0qU+s4fvw4xo8fj8jISLi4uKBGjRp49tln8ffff5fal85R+2BsmdI5aj/+++8/DBgwADVr1oSTkxO8vb3RoUMHbNu2rdS+VjtPWTXRq1cv5unpyTZs2MCOHz/OJkyYwACw7du3WztrpAIxMTEMAFu2bBk7ffq03pKdna3b79q1a8zNzY117tyZHTx4kO3Zs4c1atSIBQcHs6SkJCt+AhIbG8s8PDxYly5ddOfewoULS+1nShkuWbKEcRzH5s2bx2JiYtiqVauYTCZjEydOtNCnqt6MLdNXXnmFOTk5lTp3//vvv1L7Uplaz/Dhw1lUVBT7/PPP2YkTJ9h3333H2rdvzyQSCTt27JhuPzpH7YexZUrnqP2IiYlhr7/+Ovvmm2/Y8ePH2f79+9mLL77IALDFixfr9rPmeVotAouDBw8yAGzHjh1663v16sWCg4OZSqWyUs6IMbSBxXfffVfufiNGjGC+vr4sMzNTt+7evXtMKpWyOXPmCJ1NUg6NRsM0Gg1jjLHk5OQyL0KNLcOUlBSmUCjYa6+9pvf+pUuXMo7j2JUrV4T5IETH2DJ95ZVXmIuLS4XpUZla1+PHj0uty87OZgEBAaxHjx66dXSO2g9jy5TOUfvXrl07FhoaqnttzfO0WjSF+uGHH+Dq6ooRI0borR83bhwSEhJw9uxZK+WMmItKpcKBAwcwbNgwvZlCw8LCEBUVhR9++MGKuSPaqvXymFKGR44cQUFBAcaNG6eXxrhx48AYw759+8yaf1KaMWVqCipT6/L39y+1ztXVFQ0bNkRcXBwAOkftjTFlagoqU9vl6+sLiYSfQcLa52m1CCwuX76MBg0a6L50raZNm+q2E9s3ZcoUSCQSuLu7o0+fPvj999912+7cuYP8/HxdmRbXtGlT3L59GwUFBZbMLjGRKWWoPWebNGmit19QUBB8fX3pnLYx+fn5CAwMhFgsRkhICKZOnYq0tDS9fahMbU9mZib++ecfNGrUCACdo46gZJlq0TlqXzQaDVQqFZKTk/H555/j559/xttvvw3A+uepXU2QV1mpqamoVatWqfXe3t667cR2eXh44K233kK3bt3g4+OD27dv44MPPkC3bt1w8OBB9OnTR1eG2jItztvbG4wxpKenIygoyNLZJ0YypQxTU1Mhl8vh4uJicF86p21Hs2bN0KxZMzRu3BgAcPLkSaxZswbHjh3DuXPn4OrqCgBUpjZoypQpyM3NxTvvvAOAzlFHULJMATpH7dHkyZPxxRdfAABkMhk+/vhjvP766wCsf55Wi8ACQLlV9uaszifm16JFC7Ro0UL3unPnzhgyZAiaNGmCOXPmoE+fPrptVM72z9gypLK2D9OnT9d73atXL7Ro0QLDhw/Hl19+qbedytR2vPvuu9i+fTs++eQTtGrVSm8bnaP2qawypXPU/syfPx8TJkxAUlIS9u/fj6lTpyI3NxezZs3S7WOt87RaNIXy8fExGHVpq/kMRXXEtnl6emLgwIG4ePEi8vPz4ePjA8Bw7VNaWho4joOnp6eFc0lMYUoZ+vj4oKCgAHl5eQb3pXPatg0ZMgQuLi56w0VTmdqORYsWYcmSJVi6dCmmTp2qW0/nqP0qq0zLQueobatZsyZat26N/v37Y/369Xjttdcwb948JCcnW/08rRaBRZMmTXDt2jWoVCq99ZcuXQIAXfUfsS+MMQB8RF27dm04OTnpyrS4S5cuoU6dOlAoFJbOIjGBKWWobQ9act/ExESkpKTQOW0HGGMQiZ7+C6IytQ2LFi1CdHQ0oqOjMX/+fL1tdI7ap/LKtDx0jtqPtm3bQqVS4e7du1Y/T6tFYDFkyBDk5ORgz549euu3bt2K4OBgtGvXzko5I5WVnp6OAwcOoHnz5lAoFJBIJBg0aBD27t2L7Oxs3X4PHjxATEwMhg4dasXcEmOYUoZ9+/aFQqHAli1b9NLYsmULOI7Dc889Z6Fck8r4/vvvkZeXh/bt2+vWUZla3+LFixEdHY0FCxZg4cKFpbbTOWp/KirTstA5al9iYmIgEolQq1Yt65+nJg9Qa6d69erFvLy82MaNG9nx48fZxIkTGQC2bds2a2eNVOCll15ib7/9Nvvuu+9YTEwM27hxI6tfvz6TSCTs6NGjuv2uXbvGXF1dWZcuXdihQ4fY3r17WePGjWmCPBtx6NAh9t1337GvvvqKAWAjRoxg3333Hfvuu+9Ybm4uY8y0MtRO6jN//nx24sQJ9sEHHzC5XE4TNVlQRWV679491rFjR/bxxx+zQ4cOscOHD7O5c+cyhULBGjVqxHJycvTSozK1ntWrVzMArG/fvqUmSjt9+rRuPzpH7YcxZUrnqH2ZOHEimzlzJtu9ezc7ceIE+/7779kLL7zAALDZs2fr9rPmeVptAovs7Gz2v//9jwUGBjKZTMaaNm3Kdu7cae1sESMsX76cNW/enHl4eDCxWMz8/PzYkCFD2F9//VVq3/Pnz7MePXowZ2dn5u7uzp577jl2+/ZtK+SalBQWFsYAGFxiY2N1+5lShuvWrWP16tVjMpmM1axZky1cuJAVFRVZ6BORiso0LS2NDRkyhIWHhzMnJycmk8lY3bp12Zw5c1hGRobBNKlMraNr165llmXJe5B0jtoHY8qUzlH78tVXX7HOnTszX19fJpFImKenJ+vatSv75ptvSu1rrfOUY+xJQ3VCCCGEEEIIqaRq0ceCEEIIIYQQIiwKLAghhBBCCCFVRoEFIYQQQgghpMoosCCEEEIIIYRUGQUWhBBCCCGEkCqjwIIQQgghhBBSZRRYEEIIIYQQQqqMAgtCCCGEEEJIlVFgQQghhBBCCKkyCiwIIYQQQgghVUaBBSGEEEIIIaTKKLAghBBCCCGEVBkFFoQQQgghhJAqo8CCEEIIIYQQUmUUWBBCCCGEEEKqjAILQgghhBBCSJVRYEEIIYQQQgipMgosCCGEEEIIIVVGgQUhhBBCCCGkyiiwIIQQQgghhFQZBRaEEEIIIYSQKqPAghBCCCGEEFJlEmtnwBo0Gg0SEhLg5uYGjuOsnR1CCCGEEEIEwRhDdnY2goODIRIJW6dQLQOLhIQEhIaGWjsbhBBCCCGEWERcXBxCQkIEPUa1DCzc3NwAALGxsfD29rZybog5KJVK/PLLL+jduzekUqm1s0PMgMrU8VCZOh4qU8dDZep40tLSEBERobv+FVK1DCy0zZ/c3Nzg7u5u5dwQc1AqlXB2doa7uzv9EDoIKlPHQ2XqeKhMHQ+VqeNRKpUAYJHm/9R5mxBCCCGEEFJlFFgQQgghhBBCqowCC0IIIYQQQkiVUWBBCCGEEEIIqTIKLAghhBBCCCFVVi1HhSKEEEKIZanVat3oNMR2KZVKSCQSFBQUQK1WWzs7pBwSiQRisdimJnumwIIQQgghgmGMITExERkZGdbOCjECYwyBgYGIi4uzqQtWYphYLIa/vz88PDxsorwosCCEEEKIYLRBhb+/P5ydnW3i4oeUTaPRICcnB66urhCJqMW8rWKMQaVSISsrC48ePUJ+fj6CgoKsnS0KLAghhBAiDLVarQsqfHx8rJ0dYgSNRoOioiIoFAoKLOyAm5sb5HI5UlJS4O/vD7FYbNX80F8MIYQQQgSh7VPh7Oxs5ZwQ4rhcXFzAGLOJPkwUWBBCCCFEUNT8iRDh2NL5RYEFIYQQQgghpMoosCCEEEIIIYRUGQUWhBBCCCEm4DjOqOXEiRMWzVd0dLTRzWLCw8MxduxYsx5fqVQiMjISK1as0K3bsmWL3ncikUgQEhKCcePGIT4+XrffiRMnwHEcvv/+e4NpT506VffZSqZZ1hIeHg4A+PXXX9GrVy8EBwdDLpfD398f3bt3x6FDh0odJzw8XPd+kUgEDw8PNGjQAGPGjMEvv/xSav/09HR4enpi3759VfjmHAeNCkUIIYQQYoLTp0/rvV68eDFiYmJw/PhxvfUNGza0ZLYwYcIE9O3b16LHLO7zzz9Heno63nzzzVLbNm/ejMjISOTn5+O3337D8uXLcfLkSVy6dAkuLi4mHWfAgAGlyqBDhw4YPnw4Zs6cqVsnl8sBAKmpqWjUqBEmTJiAwMBApKWlYcOGDRgwYAC++eYbjBo1Si+tZ555BqtXrwYA5OTk4MaNG9i1axf69OmDYcOGYefOnZBKpQAALy8vTJ8+HbNnz0b//v0hk8lM+iyOhgILQgghhBATtG/fXu+1n58fRCJRqfWWFhISgpCQEKscW6VS4YMPPsD48eMNBgqNGzdG69atAQBRUVFQq9VYvHgx9u3bh5dfftmkY/n5+cHPz6/U+oCAAINl8MILL+CFF17QWzdw4EBERERg48aNpQILT09PvXR69uyJKVOmIDo6GosWLcKCBQuwcuVK3fY33ngDS5Yswffff4+RI0ea9FkcDTWFIoQQQggxs6KiIixZsgSRkZGQy+Xw8/PDuHHjkJycrLdfeHg4Bg4ciCNHjqBly5ZwcnJCZGQkvvrqK7398vLyMGvWLEREREChUMDb2xutW7fGzp07dfsYagqlVCoxZ84cBAYGwtnZGZ06dcJff/1lMM+JiYl444030KhRIygUCkRERGDRokVQqVQVft6ffvoJ8fHxGD16tFHfj/bC/f79+0btb25SqRSenp6QSIy/xx4dHY1GjRrh008/RUFBgW59QEAAevXqhQ0bNgiRVbti1zUWmzZtwsSJE+Hi4oKcnBxrZ4cQQgghxmAMyMuzdi6ecnYGzDhkp0ajwbPPPotTp05hzpw56NixI+7fv4+FCxeiW7duOH/+PJycnHT7X7hwATNnzsTcuXMREBCATZs24dVXX0WdOnXQpUsXAMCMGTPwzTffYMmSJWjRogVyc3Nx+fJlpKamlpuXiRMn4uuvv8asWbPQq1cvXL58GUOHDkV2drbefomJiWjbti1EIhFmz56NRo0a4ezZs1iyZAnu3buHzZs3l3ucgwcPwt/f3+jmX7dv3waAUjUPGo3GYCDDGDMq3fJoNBpoNBokJSXhiy++wM2bN/VqHowxaNAgrFixAufPn0enTp1067t164Z58+YhIyMDnp6eVc6rvbLbwCI+Ph6zZs1CcHAwMjMzrZ0dQkh1plYDp04Bjx4BQUFA586AodlPGQPi44GsLKCgACgs5B+1z0UioGdP/iKHEEeWlwe4ulo7F0/l5AAmtvMvz7fffosjR45gz549GDp0qG59s2bN0KZNG2zZsgWTJk3SrU9JScEff/yBmjVrAgC6dOmCY8eOYceOHbrA4o8//kDv3r0xffp03fsGDBhQbj6uX7+OrVu3Yvr06Vi1ahUAoFevXggICCjV/Cg6Ohrp6em4dOkSPD094e7ujl69esHJyQmzZs3C7Nmzyw0aTp8+jZYtW5a5Xa1WQ6VSoaCgACdPnsSSJUvg5uaGwYMH6+1XssmSOfXv3x8///wzAMDd3R27d++u8DssKSwsDACQkJCgt75ly5bQaDQ4c+aMVfu5WJvdNoV644030KVLF/Tq1cvaWSGEVGd79wLh4UBUFDByJP8YHs6vB/gLlp9+AiZNAmrVAkJDgUaNgFatgI4dge7dgf79gSFDgGef5dclJVnzExFCqujAgQPw9PTEoEGDoFKpdEvz5s0RGBhYarSo5s2b64IKAFAoFKhXr55eM6G2bdvi8OHDmDt3Lk6cOIH8/PwK8xETEwMApYKI559/vlQToAMHDiAqKgrBwcF6ee7Xrx8A4OTJk+UeKyEhAf7+/mVub9++PaRSKdzc3DBw4EAEBgbi8OHDCAgI0Ntv5cqVOHfuXKnl+eefr/DzVuSTTz7BX3/9hR9//BF9+vTBCy+8oNeUzBhl1ZxoP3vxka6qI7ussdi2bRtOnjyJq1evYsGCBdbODiGkutq7Fxg+nK+JKC4+Hhg2DGjSBLhxAygqerpNIgE8PACF4ukil/OPN24AFy7wNR6//soHIYQ4ImdnPui2FWauJXz8+DEyMjLKHCEoJSVF77WPj0+pfeRyuV7w8PHHHyMkJAS7d+/GypUroVAo0KdPH3zwwQeoW7euweNom0kFBgbqrZdIJKWO+fjxY+zfv183klJFeS4pPz8fCoWizO1ff/01GjRoAIlEgoCAAAQFBRncr1atWrpO3sUZ6qxtquLf0+DBg9GvXz9MmTIFL7zwAkQi4+61a4O94OBgvfXaz25MwOfI7C6wSEpKwrRp07BixQqrjXxACCFQq4G33iodVABP1126xD9GRAD9+vFLVFTZTS5u3eKbQt28CXTqBBw7BtSpI0z+CbEmjjNr0yNb4+vrCx8fHxw5csTgdjc3N5PTdHFxwaJFi7Bo0SI8fvxYV3sxaNAgXL9+3eB7tMFDYmIiatSooVuvUqlK9c3w9fVF06ZNsXjxYuTm5sLFxUXvYrvkhXRJvr6+SEtLK3N7gwYNDAYM1tS2bVscOXIEycnJpWpODGGMYf/+/XBxcSn1WbSf3dfXV5C82gu7CywmT56M+vXr67VNrEhhYSEKCwt1r7OysgDwIyUolUqz55FYnrYcqTwdh9nLVK0G9/vvun4QrFMnw/0gjMSdPAnJw4cV7qf66iuwl1/W7xha1mcKDwdiYiDp2xfcrVtgnTtDdegQ0LhxpfNpS+g8dTwVlalSqQRjTNdp1lFpm8doP2P//v2xa9cuKJVKtGvXzuB7in8f2u+oov20/Pz8MGbMGPz3339Yt24dcnJy4OzsXCof2v4Z27ZtQ4sWLXTv37VrF1Qqld5xBwwYgMOHD6NWrVqQSCRwc3MrNcJUeWUYGRmJO3fulNpH+7qiv4GK9iv52Qwp73s0tO/Jkyfh6ekJLy8vo8pj0aJFuHr1KubNmweZTKa3j7YzemRkpMX/1jUaDRhjUCqVEBv4v2bJ31y7Ciz27NmD/fv3499//zV6ZkkAWL58ORYtWlRqfUxMDJypk6RDOXr0qLWzQMzMHGUadPo0mmzaBKdid+jyfXxwacIEPOrQwaS05GlpCD15EhEHDxr1A/rf5cuIP3zYtGO88w46REfD4949sK5dcfq995BRRlMHe0TnqeMpq0wlEgkCAwORk5ODouJNAh2M9sJNe+Oyf//+6NmzJwYMGIA33ngDLVu2hFQqRUJCAk6dOoX+/ftj4MCBAJ6OgqR9r5Z2ZCTt+p49e6JPnz5o1KgRPDw8cPPmTXzzzTdo06aN7v3am6ja99SoUQPPP/881q1bB41Gg27duuHatWv49NNP4ebmBqVSqdt31qxZ+OWXX/DMM8/gtddeQ926dVFQUIAHDx7g6NGj+Oijj/RqPUpq3749Vq1ahcTERL1rK+2wrLm5uaU+Y3F5T0YJy8/PN7if9u+nvDSKiooMbh85ciQaN26MJk2awNvbG48ePcLOnTtx8uRJfPDBB7pjA3x5pKam4tixY7p83bp1C3v37sXp06cxZMgQzJgxo9RxTp06BW9vb4SFhZWbRyEUFRXpJh40NKJWngVHYLObwCInJwdTpkzBm2++ieDgYGRkZAB4+oeWkZEBqVRqcFKWefPmYcaMGbrXWVlZCA0NRVRUlMF2jcT+KJVKHD16FL169dLNhknsm7nKlPvhB4hXrSrVZEmRloY2q1ZBvWsX2JAh5SdSWAhu/36IvvkG3M8/gzPhblTzfv3QrGtX0zPerx80gwZB9tdf6PL++1Dv2wfWubPp6dgQOk8dT0VlWlBQgLi4OLi6upbb/t7eaT+7u7u7bt3Bgwfx8ccfY9u2bVizZg0kEglCQkLQpUsXtG3bVrevSCSCRCLRey8AXedq7fqePXvil19+wfr165GXl4caNWpgzJgxmD9/vm4fbf+I4mlt3boVISEh2LZtGzZu3IjmzZvrJnKTSqW6fd3d3XH+/HksWbIEn3zyCRISEuDm5oaIiAj06dMHNWvWLHeG7LFjx2L58uU4deoURowYoVuvLXcXF5dSn7E4bTDi5ORkcD9tf5Xy0pDJZAa3d+3aFXv27MGmTZuQlZUFT09PtGrVCj/99FOpUaFEIhHOnj2L3r17g+M4uLi4oEaNGmjTpg3ee+899O7du1T6jDH8/PPPGDlyJDw8PMrMn1AKCgrg5OSELl26GDzPKhqS2Jw4Zo6BgS3g3r17iIiIKHefZ599Fvv27aswraysLHh4eCAlJYUCCwehVCpx6NAh9O/fny5YHIRZylSt5psXldVkieOAkBAgNrZ0syjGgL//BrZsAXbsANLTn27r0AEYPRpYsoRvWmXoZ7S8tI2Vnc2PFBUTw3fuPniQH0XKTtF56ngqKtOCggLExsbqJnUjtk+j0SArKwvu7u5Gd2jW0o6CddjEWlp7d+zYMfTu3RtXrlxBZGSkxY9f0XmWmpoKX19fZGZmlhuYmYPd1FgEBgbqhk0rbsWKFTh58iQOHz5c7TvMEEJKOHWq7KAC4AOCuDh+v27d+HVZWcD27cCGDcDFi0/3DQkBxozhl/r1+XUBAfyoUBynH1xom2quXVulfhxwcwMOHQJGjAAOHAAmTACuXeNHkSKEEBuzfPlytGjRAufOnUObNm2snR2LWbJkCcaPH2+VoMLWCB5YZGVl4cyZM4iPj0d+fj58fX3RsGFDNDaxM6JCoUA37T/+YrZs2QKxWGxwGyGkmnv0yPj9zp8HvvgC2LkTyM3l1ysUwNChwNixfE1BySBh6FDg++/50aGKBzAhIXxQUWxirEpTKIBdu4B69fjaj08/BWbOrHq6hBBiZo0bN8bmzZuRmJho7axYTHp6Orp27YrJkydbOys2QZDAQqVS4fvvv8eGDRvwxx9/6Hqra3EcBx8fH7z88suYPHlymeMvE0JIlZQxTnop0dH8EK9aDRoAr7/ON3fy9i7/vUOH8s2VjJl5u7JcXPhmV+PH849jxwLUjJMQYoNGjRpl7SxYlJeXF6Kjo62dDZth9pm3f/rpJzRs2BBjxoyBi4sLli1bhl9++QUXLlzAjRs3cPr0aWzbtg0vvvgi9u3bh4YNG+KNN96ocOKVsmzZsgU5tjTJDiHEdnTuzNceVDSK3M2bgEzGz5x98iRw5QpfC1FRUKElFvNNqV56iX80Z1ChNWYM0LQpkJHBBxeEEEKIjTF7YPHKK69g1KhRePjwIQ4ePIjZs2ejR48eaNKkCerWrYt27drhpZdewscff4zY2FgcOXIEN2/exOeff27urBBCqjuxGFi3jn9eVnAREACsXs3Plr19O9ClS8WBiDWIxXw+AeCzz4AnY6YTQgghtsLsTaFiY2Ph6elp9P49evRAjx49dMPHEkKIWQ0YAEyeDGzcqD8xnVwOzJjB3/03ceQTq+nVC+jbFzhyBJg7l+/fQQghhNgIs/83NSWoMMf7CCEORq0GTpwAt2sXfC5d4l9XRkYGsGIFP9zsZ5/xQYVMxnfC/uorvoP2smX2E1RoffABn+c9e4A//rB2bgghhBAdO/uPSghxaHv38oFAVBQkY8ag07vvQlKnDr/eWHFx/KhJoaHAvHlAYiJQowYfZCQkAMeOAePGCdMPwhIaNwZefZV/PnOm4Tk0CCGEECswe1Oo7iZM3sRxnG7KdEJINbd3Lz8nRMkL5YQEfv3335c/fOvFi/zd/F27AJWKX9e4MTBnDvDCC3xthaN4/31+0r6zZ4Fvv+U/HyGEEGJlZq+x0A4ta8yi0WjMfXhCiD1Sq/lRmAzcfee066ZNK90sijF+mNf+/YFmzYBt2/igont34PBhPtgYPdqxggoACAwE3n6bfz53LlBYaN38EEIIIRAgsDhx4gRiYmKMXgghxKQZsgFAo+Fnou7UiR/F6fBhvt/BCy/wE90dO8Z3crbF0Z3MZcYMIDgYuHcP+OQTa+eGkGqF4zijlhMnTlg0X9HR0eCM/N0LDw/H2LFjzXp8pVKJyMhIrFixwuT3btmyBRzH4fz582bNkz3o1q2b3kTPN2/ehEwmwz///GO9TFWS2ZtCrVu3DsOGDUNISIi5kyaEOCpjZ8h++JAfEnbFCuDyZX6dXM73mZg9G6hVS7g82pqSk+aNG0eT5pHqQ60WdlLKCpw+fVrv9eLFixETE4Pjx4/rrW/YsKHF8gQAEyZMQN++fS16zOI+//xzpKen480337RaHhxBvXr18PLLL2P69Ok4efKktbNjErMHFsuWLcOMGTPQunVrDB8+HEOHDkXt2rXNfRhCiCMxdobsWbOAx4/5525uwKRJfBMpY9/vaMaMAdau5Zt8zZ0LfPmltXNEiPD27uWbThav5QwJ4eesKa8flhm1b99e77Wfnx9EIlGp9ZYWEhJitRu7KpUKH3zwAcaPHw8XFxer5MFc8vLy4OzsbNU8TJ06Fa1bt8aff/6Jjh07WjUvpjB7U6hHjx7h119/RZs2bbB27VrUq1cPzZs3x5IlS3D16lVzH44Q4giMnSH78WPAzw9YuhR48ABYubL6BhUAf4d27Vr++aZNwNatVs0OIYLTDvJQsulkfDy/3pQR5ARWVFSEJUuWIDIyEnK5HH5+fhg3bhySk5P19gsPD8fAgQNx5MgRtGzZEk5OToiMjMRXX32lt19eXh5mzZqFiIgIKBQKeHt7o3Xr1ti5c6duH0NNoZRKJebMmYPAwEA4OzujU6dO+OuvvwzmOTExEW+88QYaNWoEhUKBiIgILFq0CCrtgBjl+OmnnxAfH4/Ro0eX2vb777+jR48ecHNzg7OzMzp27IiDBw8aTCc9PR3jxo2Dt7c3XFxcMGjQINy9e1dvn3///RcDBw6Ev78/5HI5goODMWDAADws9nfBGMPnn3+O5s2bw8nJCV5eXhg+fHiptLp164bGjRvjt99+Q8eOHeHs7Izx48fjueeeQ1hYmMH+wO3atUPLli1NPhZjDKtWrUJYWBgUCgVatmyJw4cPG/weWrVqhQYNGmDDhg0Gt9ssJrDff/+dTZs2jYWFhTGRSMQiIyPZO++8w/755x+hD12mzMxMBoClpKRYLQ/EvIqKiti+fftYUVGRtbNCKmvPHsY4jl/4XhX6i48PYx9/zFhurrVzansWLuS/I4WCsb//tnZuykTnqeOpqEzz8/PZ1atXWX5+ftUPplIxFhJi+PcB4H87QkP5/SzslVdeYS4uLrrXarWa9e3bl7m4uLBFixaxo0ePsk2bNrEaNWqwhg0bsry8PN2+YWFhLCQkhDVs2JB9/fXX7Oeff2YjRoxgANjJkyd1+73++uvM2dmZffTRRywmJoYdOHCArVixgn3yySe6fRYuXMhKXtq98sorjOM4Nnv2bPbLL7+wjz76iNWoUYO5u7uzV155Rbffo0ePWGhoKAsLC2Nr1qxhv/zyC1u8eDGTy+Vs7NixFX4H48ePZ/7+/qXWnzhxgkmlUtaqVSu2e/dutm/fPta7d2/GcRzbtWuXbr/NmzczACw0NJSNHz+eHT58mG3cuJH5+/uz0NBQlp6ezhhjLCcnh/n4+LDWrVuzb7/9lp08eZLt3r2bvfHGG+zq1au69CZOnMikUimbOXMmO3LkCNuxYweLjIxkAQEBLDExUbdf165dmbe3NwsNDWWffPIJi4mJYSdPnmQ//vgjA8COHj2q93muXbvGALCPP/7Y5GNpy+fVV1/Vfb4aNWqwwMBA1rVr11Lf3aRJk5ivry/TaDTlfvcVnWcpKSkMAMvMzCw3HXMQPLAo7uzZs2zOnDmsTp06TCQSsYiICDZ79mxLZoExRoGFI6ILFgexZQtj7u56FwsasZixSZMYKyy0du5sl1rNWP/+/HcWFsaYjf620XnqeCwaWMTElB1UFF9iYqp+LBOVDCx27tzJALA9e/bo7Xfu3DkGgH3++ee6dWFhYUyhULD79+/r1uXn5zNvb2/2+uuv69Y1btyYPffcc+Xmo2Rgob0Inj59ut5+27dvZwD0AovXX3+dubq6stjYWJaens7UajVjjLHVq1czAOzKlSvlHrtBgwasb9++pda3b9+e+fv7s+zsbN06lUrFGjduzEJCQnQXzdrAYsiQIXrv/+OPPxgAtmTJEsYYY+fPn2cA2L59+8rMy+nTpxkA9uGHH+qtj4uLY05OTmzOnDm6dV27dmUA2LFjx/T2VSqVLCAggI0cOVJv/Zw5c5hMJtNdQxp7rPT0dKZQKMr8fIYCiy+//JIBYNeuXSvzszJmW4GFRSfIa9u2LVauXIlbt27h77//xssvv1xmVRghpBpJTQXeeQd4800gKwsAwIKCcGPYMKgyMoDPP3e8IWPNSSTih9qtXRu4fx8YObLyM5YTYquMHeTB2P0EdODAAXh6emLQoEFQqVS6pXnz5ggMDCw1WlTz5s1Rs2ZN3WuFQoF69erh/v37unVt27bF4cOHMXfuXJw4cQL5+fkV5kM7+ubLL7+st/7555+HRKLfzfbAgQOIiopCcHCwXp779esHABV2Ik5ISIC/v7/eutzcXJw9exbDhw+Hq6urbr1YLMbo0aPx8OFD3LhxQ+89JfPasWNHhIWF6T5LnTp14OXlhbfffhsbNmww2Mz+wIED4DgOo0aN0vssgYGBaNasWanv38vLq9Q8bBKJBKNGjcLevXuRmZkJAFCr1fjmm2/w7LPPwufJYBnGHuv06dMoKCgo8/MZov0+4+PjDW63RVabebt58+ZYvHgxrly5Yq0sEEKqQq0GTpwAdu7kHytzIZuSws+OHR4OLFsGZGcDTZsC330HVWwsro8ezY/6RCrm5cW3L3d2Bn75BXjvPWvniBDzMrY/lQ30u3r8+DEyMjIgk8kglUr1lsTERKSkpOjt72NgRDe5XK4XPHz88cd4++23sW/fPkRFRcHb2xvPPfccbt26VWY+UlNTAQCBgYF66yUSSaljPn78GPv379f1B5HL5ZBKpWjUqBEAlMpzSfn5+VAoFHrr0tPTwRhDkIEyCQ4O1sujVsm8atdp9/Pw8MDJkyfRvHlzzJ8/H40aNUJwcDAWLlwIpVKp+yyMMQQEBJT6/s+cOVPqsxjKHwCMHz8eBQUF2LVrFwDg559/xqNHjzBu3DjdPsYeq6yyKGsdAN33aUwQaSvMPirU119/bfS+HMcZ7ORDCLFxVR2VJTkZWL0a+OwzIDeXX9eiBX8xPHgwfwf+yT8IYoKmTflO3CNH8oFamzbAc89ZO1eEmId2kIf4eIOTaYLj+O2dO1s+byX4+vrCx8cHR44cMbjdzc3N5DRdXFywaNEiLFq0CI8fP9bVXgwaNAjXr183+B5t8JCYmIgaNWro1qtUqlIX9L6+vmjatCkWL16M3NxcuLi4QCR6ev9ZGwiUxdfXF2lpaXrrvLy8IBKJ8MhALVJCQoLufcUlJiaW2jcxMRF16tTRvW7SpAl27doFxhguXryILVu24P3334eTkxPmzp0LX19fcByHU6dOQW7g5lTJdWXN/dGwYUO0bdsWmzdvxuuvv47NmzcjODgYvXv31vvcxhyreFkY+nzh4eGl1mu/z5LfkS0ze2BhymQrFFgQYoe0o7KU/MeuHZXl++/LDi6SkoAPPuCbNuXl8etateIDikGDHHtCO0t56SXgr7/40aLGjAHOnQPq17d2rgipOrGYv3kxfDj/W1H8N0j727F2rUXnsyjLwIEDsWvXLqjVarRr187s6QcEBGDs2LG4cOEC1q5dW+bwqNpJ17Zv345WrVrp1n/77belRnoaOHAgDh06hNq1a0MsFsPd3V0vsKhIZGQk7ty5o7fOxcUF7dq1w969e7F69Wo4OTkBADQaDbZt24aQkBDUq1dP7z3bt2/HsGHDdK///PNP3L9/HxMmTCh1TI7j0KxZM6xZswZbtmzRTSg3cOBArFixAvHx8Xj++eeN/gyGjBs3DpMmTcLvv/+O/fv3Y8aMGRAX+xsz9ljt27eHQqEo8/MZCizu3r0LkUiE+nb0G272wCI2NtbcSRJCbIVazddUGLpbyBj/z33aNODZZ/X/uScm8gHF+vWAtkq3TRtg4UKgf38KKMxt1Srgn3+A334DhgwBzp7l5/0gxN4NHcrfvDBUY7p2rcXmsajIiy++iO3bt6N///5466230LZtW0ilUjx8+BAxMTF49tlnMWTIEJPSbNeuHQYOHIimTZvCy8sL165dwzfffIMOHTqUOedCgwYNMGrUKKxduxZSqRQ9e/bE5cuXsXr1ari7u+vt+/777+Po0aPo1KkTJkyYgGbNmqGoqAj37t3DoUOHsGHDhnLnyOjWrRvef//9UkHO8uXL0atXL0RFRWHWrFmQyWT4/PPPcfnyZezcubNUbcH58+cxYcIEjBgxAnFxcXjnnXdQo0YNTJ48GQDfp+Hzzz/Hc889h1q1aoExhr179yIjIwO9evUCADzzzDN47bXXMG7cOJw/fx5dunSBi4sLHj16hN9//x1NmjTBpEmTjPreX3rpJcyYMQMvvfQSCgsLS91AN/ZYXl5emDVrFpYsWaL3+aKjo8tsCnXmzBk0b94cXl5eRuXVJgjePdwG0ahQjodGm7EQU0dlSUhgbNo0fhhU7ba2bRk7eJCxCobPozKtosRExoKD+e88IoKxDz5gLDXVqlmiMnU8Fh0VqjiViv+d2bGDf7TCELPFlRwVijF+VKHVq1ezZs2aMYVCwVxdXVlkZCR7/fXX2a1bt3T7hYWFsQEDBpRKs2vXrnojBc2dO5e1bt2aeXl5MblczmrVqsWmT5+udy1jaLjZwsJCNnPmTObv788UCgVr3749O336NAsLC9MbFYoxxpKTk9mbb77JwsLCmFQqZd7e3qxVq1bsnXfeYTk5OeV+B7dv32Ycx7Fvv/221LZTp06x7t27MxcXF+bk5MTat2/P9u/fr7ePdlSoX375hY0ePZp5enoyJycn1r9/f73v6/r16+yll15itWvXZk5OTszDw4O1bduWbdmypdRxv/rqK9auXTvdcWvXrs3GjBnDzp8/r/c9N2rUqNzPNnLkSAaAPfPMM2XuY8yxNBoNW758OQsNDWUymYw1bdqU7d+/v1RZM8ZYdnY2c3Z2LjXalCG2NCoUx5ihW4/mlZ2djdOnTyM1NRW+vr5o3759pdoXmktWVhY8PDyQkpJisMMUsT9KpRKHDh1C//79IZVKrZ0dx7VzJ99+vyKffQbcvAl88QVQUMCva9+er6Ho08eoGgoqUzM4exYYMIAfdQsAFArg5ZeBKVP4Pi0WRmXqeCoq04KCAsTGxuomdSO2T6PRICsry+SmUAB0o2CVNekbMd7//d//4a233kJcXFyFNRYVnWfa6+/MzMxSNVXmZvamUCWtXr0aixYtQl5eHhhj4DgOzs7OWLRoEWbMmCH04Qkh5mTsaCvTpj3tfN2xIx9Q9OpFTZ5MpVYDp07xw2cGBfGdUk1pP96uHT9D+Y4dwKefAhcuAP/3f/zSsSMfYLRpA7i68ouLC99xvjwqFVBYyAeMhYX6z7WPRUV8+ZdYuMJChP37L0RxcU/TUqufLhpN+Yuh+jHtrLja1yWfl4fj+M9raJHJ+BG2tIuLy9PngYFA3bo0YhkhJSxfvhwtWrTAuXPn0KZNG2tnx26pVCqsXLkS8+bNs69mUBA4sPj6668xZ84c9OvXD2PHjkVwcDASEhKwdetWzJ49G35+ftR5mxB7UtGoLFpKJfDMM0B0NNCjBwUUlVHVkbe0nJ2BCROAV18F/vyTDzC+/55//uefpfd3cXkaaGg0pYOGKsyPIQHQvNLvtjFiMVCnDtCwYemF5lwh1VTjxo2xefNmgyMfEePFxcVh1KhRmDlzprWzYjJBm0K1aNECjRo1wrZt20ptGzVqFK5evarrwW9J1BTK8VATCwvSjgoFGA4uIiOBTz6pckBRrcu0rJG3tN9neSNvGePRI+DLL/lJ9R4/BnJynt75N4VYzN+1Vyj0H2UyQCrlF4lE91wjkeBxSgoCgoMh0m4Ti/UXkah0TYL2dfHH4ou2lkX7/ZR8LIu2xqPkolbzwXFeHj8ccl7e0yU3l68FejKRYylBQcCaNcDzz1eLgJqaQjmeqjSFItZRbZpCXb9+HcuXLze4bdSoUSaPiEAIsTKlkr9rPmYMsHv30/4TAH8xOXcuX0tRDS6oBFPZkbdMERTED/GrnUSPMb4ss7P5IEP7KBaXDhq0j3I5HxiY8tGUSvz15CJUZM/BImNAQgJw9ar+cvkyH7S9+CKweTPf16h2bWvnlhBCLEbQwMLJyanUZClaaWlpuvGMCSE2Kjubby7z++/8cvbs0+FiterVA157jb8YNvFCkxhw6pR+86eSGAPi4vj9noxRX2UcBzg58Yu/v3nSdGQcB9SowS9PhrcEwDcVW7mSn5zw55+Bxo2BBQuA2bOpeRQhpFoQ9Cqgc+fOiI6ORrdu3fRmbExMTMT777+PLl26CHl4QoiplEo+ePj1V+DYMeDMGb6DbXHe3nz/iU6d+AvbNm2ohsKcDMxQW6X9LKmqnc3tnVzO1wK99BIwaRJ/Di1YAGzfzs/h0rWrtXNoNRYYgJKQasuWzi9BA4tly5ahY8eOqFOnDnr06IGgoCA8evQIx48fh1Qqxd69e4U8PCHEGHFxwJ49fDBx8iTfBKa48HCgSxc+kOjUiZ/FmdrdCsfYkbeM3c9SzNXZvCz2FLTUrQscPcoPzzx9OnDtGh+ET57Md56vRoG4tt9FXl4etVIgRCC5ubngOM4m+iMKGlg0atQI586dw8KFCxETE4PU1FT4+Pjgueeew8KFC0tN404IsaB794Dly/m24NqhYQHAx4fveN2zJ/9Yq5bVslgtVTTyFsfx2zt3rvwxzH2RXlZn8/h4fn1VO5sLHbQIgeP4OV/69QPmzePndPn8cyAq6ungB9WAWCyGp6cnkpKSAADOzs6lZlomtkWj0aCoqAgFBQXUeduGMcagUqmQlZWFrKwseHp6QmwDN1sEbxBdr1497Ny5U+jDEELKUvIiMiSEbwe+ZcvTZk6dOvGdgXv2BJo2pRoJaxKL+Qvm4cP5i9PiF+vaC7K1aysfCJj7It2UzuaVza+QQYvQvLyADRuAgADg/feBOXOAgQP5TvDVRGBgIADoggti2xhjyM/Ph5OTEwWBdkAsFiMoKAgeHh7WzgoACwQWhBArMnQRWVyPHvzkdVW5+03Mb+hQ/oLZUACwdm3lL6SFuEg3pbP5M8+YlrYlRsiylDlzgE2bgNhYPoh7+21r58hiOI5DUFAQ/P39oSxeO0psklKpxG+//YYuXbrYRNMaUjaJRAKxWGxTAaDggcW///6LHTt24P79+ygoPjQl+B+bH3/8UegsEFI9lXURqbV0KTB/vmXzRIw3dCh/wWyuJktCXaQL2dncGiNkCcXFhW96+Mor/Lk3dixfi1GNiMVim2iqQconFouhUqmgUCgosCAmE3zm7XHjxkEkEsHf3x+yEsPt2VKERYhDUav5jqJlBRUcxzfPePtt27/TW52Jxea7YBbqIl3Izub2PEKWIaNG8ZNHnj/Pjx71xRfWzhEhhJiVoA2ply5digEDBiApKQnx8fGIjY3VW+7evSvk4Qmpni5cALp352dULkvxi0hSPQh1ka7tbF7WjSKOA0JDK9fczhIjZKnVwIkT/AhOJ07wr4UiEgEffcQ/37QJuHhRuGMRQogVCBpYxMfH43//+x+8vLyEPAwhBAD+/ht47jmgeXPgt9+Me4+93OklVSfURbq2szlQOrioamdzIYMWgG8uGB7Oj9Q0ciT/GB7OrxdK5858E0WNBpgxo+xaRUIIsUOCBhYtWrRAfHy8kIcghPz+Oz+kZevWwI8/8hdbUVHGvdfW5kIgwhHyIl3b2bxGDf31ISFVG7VJyKBF2wepZPMwbUd2IYOLVav4mbiPHQMOHBDuOIQQYmGCBhYffPABVqxYgYtU3UuIeTHGT2jXrRt/IXjkCH9xNXo0cPUqPzmXkHd6if0R8iId4IOHe/eAmBhgxw7+MTa26kPBChG0VNSRHeA7sgvVLCoigp84DwBmzQKKioQ5DiGEWJignbfbt2+PoUOHokWLFggKCoK3t7fedo7jcOHCBSGzQIhjYYy/w7lkCfDXX/w6qRQYN47viF18Mjsh50Ig9kmoYWy1zNnZvDhzj5BlC6NNzZ/PT0558yY/cd60acIchxBCLEjQwGLlypVYvnw5/Pz8EBYWVmpUKEKIkVQq/oJw+fKnHT4VCuC114DZs/kLw5KEvogk9sncF+mWYs6gxRZGm3J3BxYvBl5/HVi0iK9t9PER7niEEGIBggYW69atw/jx4/HFF19Ueezq48ePY9u2bfjzzz8RFxcHT09PtG7dGu+99x5atWplphwTYmMKC4GtW/k22Xfu8OtcXYEpU/imFBWNg2+vF5FEWELVLNgLS4w2ZYxXXwU++4y/WbBoEfDxx8IejxBCBCZoYJGVlYWRI0eaZUKc9evXIzU1FW+99RYaNmyI5ORkfPjhh2jfvj1+/vlndO/e3Qw5JsRK1Gr9i//mzfnhKD/66OldUx8fvvZhyhSgRLPCclX3i0hCStJ2ZI+PN9zPguP47UL3QRKL+XO8Z8+nzaGKN2ckhBA7I2hg0alTJ1y9etUsF/2fffYZ/P399db17dsXderUwbJlyyiwIPZr797SzZWK94uoUYPv4DlxIj97LyGkarQd2W2hD1KPHkDv3sAvvwArV9KkeYQQuyboqFDr1q3Dhg0b8OOPP6KoiqNelAwqAMDV1RUNGzZEXFxcldImxGrKGvJSe6EzeTJw9y5/J5OCCkLMR6ghcitjwQL+ccuW8juVE0KIjRM0sGjdujVu376NoUOHwtnZGe7u7nqLh4dHldLPzMzEP//8g0aNGpkpx4RYkFrNBw5lTZDFccD+/dQfghChCDVErqk6dwa6dOGHnV292rLHJoQQMxK0KdSwYcPAlTWOvhlMmTIFubm5eOedd8rdr7CwEIWFhbrXWVlZAAClUgmlUilY/ojlaMvRLsqTMXCnTkE0bx5Ejx+Xux/i4qCKiQHr2tVy+bMRdlWmxCg2W6bPPPP0uUbDLxbGzZ0LyW+/gW3cCNXs2YCBWnpbZLNlSiqNytTxWLIsOcbKul1q2959910sWbIEn3zyCaZOnVruvtHR0Vi0aFGp9Tt27ICzs7NQWSRED6dWI/jPP1Fn3z54akd4MsL5GTMQ36WLgDkjhFgdY+gyZw68bt3CzWHDcG30aGvniBDiIPLy8jBy5EhkZmbC3d1d0GPZZWCxaNEiREdHY+nSpZg/f36F+xuqsQgNDcWjR4/gQ+OGOwSlUomjR4+iV69ekEql1s6OvuxsiL76CqJPPgH34AEAgCkU0PTsCfGBAxW+XXX0aLWtsbDZMiWVQmVaPm7/fkiGDQNzdYXq9m3TRn+zEipTx0Nl6nhSU1MRFBRkkcDC7E2hvv32Wzz//PMmvSchIQGxsbF4pnh1dBm0QUV0dLRRQQUAyOVyyOXyUuulUimdNA7Gpsr04UN+XPqNG4HMTH6dnx8wdSq4SZMg9vYGwsMrHPJSEhVVrftZ2FSZErOgMi3Dc88BTZuCu3gR0g0bgIULrZ0jo1GZOh4qU8dhyXI0e+ftKVOmoHnz5ti0aZOuL0NZ/v77b0yZMgV169bFhQsXKkx78eLFiI6OxoIFC7DQjn5wSTVz9izw0kt80PDBB3xQUb8+H2Dcvw+89x4fYGiHvASeDnGpZekhLwkh1icSAdo+g+vWAdnZ1s0PIYSYyOw1Frdv30Z0dDTeeustTJ06FS1atEDLli3h7+8PhUKBtLQ03LlzB2fOnMGjR4/QuHFj7N27F3369Ck33Q8//BDvvfce+vbtiwEDBuDMmTN629u3b2/uj0KI8ZRKfujYtWuB4n+bXbvyc1D0789fNJSkHfKy5DwWISF8WpYenYYQYl3DhvE3Im7cANavB+bMsXaOCCHEaGYPLDw8PLBmzRq899572Lx5Mw4dOoStW7ciLy9Pt0+tWrXQt29fvPzyy4iKijIq3f379wMAjhw5giNHjpTaboddRYg9KjlDdqNGwFdfAZ9++jQwkMn4Gou33gJatKg4zaFDgWef1U+3c2eqqSDEnpX8rTD2nBaLgXnzgLFjgQ8/BKZOBWiQEUKInRBsuFkvLy/MmDEDM2bMAMDPOZGfnw8fH59KtfU6ceKEmXNIiIkMzZBdnL8/MGkS8MYbQGCgaWmLxUC3blXOIiHEBhj6rQgJ4Zs3GVMLOXIkEB3Nz7GxaRPwv/8JlVNCCDErQSfIK87DwwOBgYHUEYjYp7JmyNaaMoXvPxEdbXpQQQhxHGX9VsTH8+v37q04DakUmDuXf/7BB0CxUQ0JIcSWWSywIMRuXbkCvPJK+TNk//QTfzFACKm+1Gq+psLQb4V23bRp/H4VGTsWCA7mA5SvvzZnLgkhRDAUWBBiSGYmP4pThw5A48ZATk7Z+z6ZIRunTlkuf4QQ23PqVNm1moBpvxVyOTB7Nv986dLy0yWEEBtBgQUhWmo1cPQo3745MBB4/XV+hKeSQ8GW5dEjYfNHCLFtxv4GGLvfa6/xtRb37wNNmwLffVf5vBFCiAVQYEGqN40G+Osvvj1zeDjQuzewcydQUMCP+LR6NT8crDGCggTNKiHExhn7G2Dsfs7OwIkTQJs2QHo68PzzwPjxNL8FIcRmmT2w+O2335BTXrMRQqxNqQSOHeOHcaxZE2jXDli5km9q4OXFd8Q+dw64dAmYOZMfCjYkpOyaC44DQkP54SQJIdVX587m/62oWxf44w9+4jyOAzZv5oexPnvWPHkmhBAzMntgERUVhatXr5o7WUKqJjGRr3kYOxYICAB69gQ++4wfqcXVlb8T+O23QEICPydF69ZPLw5ohmxCiDGE+q2QSoElS4CTJ/mbIXfuAM88AyxeDKhUVc42IYSYi9kDC5qojlgdY8C1a8CXX/KBRJ06fNODESOArVv5JgW+vsCrrwIHDgDJycDu3fx2hcJwmtoZsmvU0F8fEsKvpxmyCSGAsL8VnTsDFy7wE3Cq1cB77/HNpFauBC5fLnvkOkIIsRDBJsgjRHC5ucDdu8CdOxDdvImmMTEQf/YZ8PffQFqa/r4cx3d+jIoChgzh7/aZeteQZsgmhBhDyN8KT09gxw5gwABg8mTgv//4Ze5cvjajf39+6d4dcHGp+vEIIcQEggQWnLGj6BBiqj//5O/O/fUX37zpCTGAiOL7KRRA/fpARAQ/ZOyECYC3d9WPTzNkE0KMIfRvxcsvAz168BPuHTwIHD8OPHgAbNjAL3I5f/y+ffmlfn3jR7gjhJBKEiSwiIqKgkhUcSsrjuOQmZkpRBaIo/njD2DRIn442OK8vIDataGJiMBtxlCrd29I0tKATz7hmwxcuADs28e/XreOmiwRQhxHYCBfazF5MpCfD8TE8EHGwYP8ELU//8wv06cDYWFPg4zu3QF3d2vnnhDigAQJLLp16wY/Pz8hkibVzR9/ANHRwK+/8q8lEr7fxKuvAvXq6Woh1Eolrh06hFqFhfz8EyXbGsfHA8OHU38IQohjcnJ62gzq00/5fmZHjvDLyZN8oPHFF/wikQDt2/NLu3b8Y0iItT8BIcQBCBJYvPfee2jbtq0QSZPq4uxZYMEC/YBi3Dhg/nx+vglD1GqIZ8ww3IGRMb4ZwLRpfNtn6hdBCHFUHAc0bMgvM2bw/dFOnnwaaNy6Bfz+O79oBQfzQUa7dvxwtmFh/NC4zs7W+xyEELtDnbeJ7fnnH6BTJ34YRamUDyjmzSs7oHjC5+pVcPHxZe/AGBAXx3eopH4ShBB7pFab3incxeVpbQbAD1d76hR/A+fsWeDiRX6o7R9+4JfivL35AEO7BAXxzajc3AA3N3AKBbyvXuVHwfL05Pt2yGT8b7dM9vQ59e8gpFqgwILYFpUKmDiRf+zdG9i4kb9zZgRFerpxx3j0qAoZJIQQK9m7F3jrLX4yT62QENP7j9WuzS9jx/Kv8/L40fS0gcbVq/xNmOxsfoS9tDS+v5oBEgCdAb42uTwSCR8AlVy060WisheOK/ux5HNTF0D/uYcH8NFHfBBFCDEZBRbEtnz6KV9j4ekJfP01P5mdkQq8vIzbMSiocnkjhBBr2buX7ycmRP8xZ2e+5qPkjOCZmXyA8eAB/xgXBzx+zAcc2dlATg5YVhZyHz+Gi0YDLjcXUCqBwsLSx1Cp7GcyPy8v/qYWIcRkZg8sNBqNuZMk1cWDB3y/CgBYtcqkoAIAUhs2BKtRA1xCguF+FhzH390r+c+TEEJsmVrN11QI3X/MUDOrxo35pQwqpRLHDh1C//79IZVKn+ZJrQaKiviloIBPNzGRn5y0Zcunx1Op+EfGAI2GX9Rq/eeMPV00mqePKhU/MWBqKh8MNGjA115o96lo0eZV+xgXx9e8bNsGDB7MB0/mnq+oMk3ZCLEjVGNBbANjwNSpfCfDTp34UZ9MJRZD/dFHkLz4Iv+Ptvg/YW0199q19CNOCLEvp07pN38qyRz9x8zVzArgf28lEn45csR86QqVXy3GgM8/59McNMi8+RUqz1oUsBAbUfFkE4RYwt69wP79fCe/L77g7zpVAhsyhG8SUKOG/oaQEBpqlhBin4ztF1bZ/mPaZlYlgxdtM6u9e6tHuj/8YDiAq2q6gHB51qYdHg5ERQEjR/KP4eFVS5OQSqLAglhfZibw5pv887ff5odIrIqhQ4F79/jJonbs4B9jYymoIITYJ2P7hVWm/1hFzawAvpmVWl090jWkKukWT9vceQaEDVgIqQQKLIj1vfMOf6etbl3+uTmIxXyTgJde4h+pSpgQYq86d+ZrXcsaspXj+FGMKtN/zJRmVpSu6ekKmbaQAUvxY5w4AezcyT9WJS1SLVBgQazrzBm+TSsAbNgAKBTWzQ8hhNgasZhvhw+UDi6q2n9MqGZWlK7waQsZDAHgfviBmlgRk1FgQaxHqQRee43/8RszBuje3do5IoQQ2zR0qDD9x4RqZkXpCp+2gMFQ0OnTEL/4IjWxIiajwIJYz5o1wKVLgI8P8OGH1s4NIYTYNiH6jwnVzIrSFT5toQIWtRpNNm0StokVcVgUWBDrePAAiI7mn3/4IT+2OSGEkPKZu/+YUM2sHCldrco2NxMqzwIFLNzvv8MpNRVlpFrlJlbEsVFgQaxj/nwgPx/o0oVvBkUIIcQ6hGpm5SjpKhRVH65ciDwLFbAIPbwxcWg0QR6xvPPnge3b+ecffVT23RZCCCGWMXQoP3O3uSdZs+d0//4bmDULKCwEWreuWrpC5VkbsBiaeG/t2soFLEL2NyEOjwILYlmMATNn8s9HjwZatbJufgghhPC0zawoXT7dbt2AQ4eA48eBjRuBJUvMl7Y5mTlgYZ06Id/HB4q0NHCG+llwHB+4VKa/CXF41BSKWNaPPwK//cZXLS9dau3cEEIIIWV74w3+cdMmfiRDW2XOvjdiMS5NmMA/N/fwxsThUWBBLEepBObM4Z/PmMF3KqPJdwghhNiq554DAgOBx4+BffusnRuLedShA9S7dpm/HwtxeBRYEMvZsAG4dQvw9wfmzuXHwabJdwghhNgqqRTQ3r1fv966ebEwNmSI+Yc31qKbig6LAgtiGRkZwKJF/PP33weOHuUn2aHJdwghhNiy114DRCL+wvr6dWvnxrLMPbwxQDcVHRwFFqRi5rizsHQpkJoKNGwIjB3Lj2BBk+8QQgixdaGhwMCB/PO5cw3/7yLG2buXbio6OAosSPnMcWchNhb4+GP++QcfAKdPl/5RKY4m3yGEEGJLlizhm0X9+CPw7bfWzo19UqvppmI1QIEFKZu57izMmwcUFQE9ewL9+tHkO4QQQuxLkybAO+/wz998E0hOtm5+7NGpU3RTsRqgwIIYZq47C2fOALt380PUrV7NP9LkO4QQQuzNvHl8gJGczP9/JKax1E1F6hhuVRRYEMPMcWdBpXo6Gd64cUCzZvzzzp35IevKmnGb4/g2rTT5DiGEEFshkwFffcV35N65E/jpJ2vnyL5Y4qYidQy3OrsKLHJycjBt2jQEBwdDoVCgefPm2LVrl7Wz5Ziqemfh3j2ga1fgzz8BZ2dg8eKn28RiYN06/jlNvkMIIcRetG4NzJrFP3/jDX7EQ2IcoW8qUsdwm2BXgcXQoUOxdetWLFy4EIcPH0abNm3w0ksvYceOHcIcUKjqNHtItyp3Fnbt4msn/vwTcHcHtm8HgoP19xk6lJ9khybfIYQQYk+io4F69fgba9paeVIxIW8qWqJjuD1cu5WRLvf77+ZJ0xjMThw8eJABYDt27NBb36tXLxYcHMxUKpXRaWVmZjIALCUlpeyd9uxhLCSEMf5Pkl9CQvj1VWEv6apU/Ps5Tj9N7cJxjIWG8vtpZWUx9sorT/fp2JGx2NiKjxMTw9iOHfyjCeVYXFFREdu3bx8rKiqq1PuJ7aEydTxUpo6n2pbpqVNP/z/+8ou1c2NWgpepoeuV0NCqXQfFxBi+Vim5xMSYL8+2eO1WRrqZAAPAMjMzq5auEeymxuKHH36Aq6srRowYobd+3LhxSEhIwNmzZ813MKGq0+wpXVPvLPz1F9CiBbB1K9/+9L33gJMn+baNFR3H3JPvEEIIIULq1AmYOpV/PnEikJNj3fzYk6FDzT+jt5Adw+3p2q28dC3EbgKLy5cvo0GDBpBIJHrrmzZtqttuFkJVp9lbukDZzZWCg4E1awAvL77qbu5c4JlngDt3+PaRJ07ws2yXKCtCCCHEYSxbxt88u3+fHzGKGM/cNxWF6hhub9du5aVrIXZz5ZeamopatWqVWu/t7a3bXpbCwkIUFhbqXmdmZgIAcrt0KRWoICcHnBGjIbF69QBXV+M/gD2kq/1DLF4pB/Dvq1kTyMoCsrLAxcfzf/AlaAYPhvqjjwBPT36WbQtSKpXIy8tDamoqpFKpRY9NhEFl6nioTB1PdS9TbvVqSIYPBz79FOzoUb5Gv6zFTjDG0CI7G7nvvgvOXvLNGH8zU6WCoRwzgN/+v/+ZVhb2cO1mRLpZumSFDzjsJrAAUO4feHnbli9fjkWLFpVaH3b1auUzc/du5d/rSOkW99NPNPweIYSQ6unGDWvngJRHpQIuXRImbTu5dktNTYWHh4dZ0yzJbgILHx8fg7USaWlpAJ7WXBgyb948zJgxQ/c6IyMDYWFhePDggeBfMLGMrKwshIaGIi4uDu7u7tbODjEDKlPHQ2XqeKhMHQ+VqePJzMxEzZo1y71WNhe7CSyaNGmCnTt3QqVS6TVfuvQk+mzcuHGZ75XL5ZDL5aXWe3h40EnjYNzd3alMHQyVqeOhMnU8VKaOh8rU8YhEwnettpvO20OGDEFOTg727Nmjt37r1q0IDg5Gu3btrJQzQgghhBBCiN3UWPTr1w+9evXCpEmTkJWVhTp16mDnzp04cuQItm3bBjENU0r+v707D4uq6uMA/r2swyIgm4Ib7rjgbpavKKikWaagJlpplpapvZop7gsumdqmlqmVoZi7aLmUZUCZqblk2au+vu4m5gKK7DBw3j9uMzDMADMww8zA9/M89xnm3jNnzszhwv3dsxERERGR2VhNYAEAcXFxmDVrFubOnYuUlBQEBgZiy5YtiIyMNCgfR0dHzJs3T2f3KLJOrNOqh3Va9bBOqx7WadXDOq16KrNOJVEZc08REREREVGVZjVjLIiIiIiIyHIxsCAiIiIiogpjYEFERERERBVWbQKL9PR0TJo0Cf7+/lAoFGjXrh22bt1q7mKRHhITEyFJks7t2LFjGmlPnz6N3r17w9XVFR4eHoiIiMCVylh1nEqVlpaGqKgoPPnkk/Dx8YEkSZg/f77OtIbU4apVqxAYGAhHR0c0bNgQ0dHRyMvLM+EnIRV96/Sll17See4GBgbqzJd1ah7x8fF4+eWXERgYCBcXF9SpUwcDBgzAqVOntNLyHLUO+tYpz1HrcebMGTz99NOoX78+nJyc4OnpiSeeeAKbNm3SSmu281RUE2FhYcLDw0OsWbNGxMfHi9GjRwsA4ssvvzR30agMCQkJAoB4++23xdGjRzW2tLQ0dbrz58+LGjVqiODgYLF//36xa9cu0apVK+Hv7y/u3r1rxk9AV69eFe7u7qJ79+7qc2/evHla6Qypw0WLFglJksSMGTNEQkKCWLZsmXBwcBBjxoyppE9VvelbpyNHjhROTk5a5+6ZM2e00rJOzWfw4MEiNDRUrF69WiQmJoodO3aIxx9/XNjZ2YkffvhBnY7nqPXQt055jlqPhIQE8dprr4nY2FgRHx8v9u7dKyIjIwUAsXDhQnU6c56n1SKw2L9/vwAgNm/erLE/LCxM+Pv7C6VSaaaSkT5UgcWOHTtKTTdkyBDh7e0tUlNT1fuuXbsm7O3tRVRUlKmLSaUoKCgQBQUFQggh7t27V+JFqL51eP/+faFQKMSrr76q8frFixcLSZLEf/7zH9N8EFLTt05HjhwpXFxcysyPdWped+7c0dqXlpYmatWqJXr16qXex3PUeuhbpzxHrV+XLl1EvXr11M/NeZ5Wi65Qu3fvhqurK4YMGaKxf9SoUUhKSsLx48fNVDIyFqVSiX379mHQoEFwc3NT72/QoAFCQ0Oxe/duM5aOVE3rpTGkDr/99ltkZ2dj1KhRGnmMGjUKQgjs2bPHqOUnbfrUqSFYp+bl6+urtc/V1RUtW7bEzZs3AfActTb61KkhWKeWy9vbG3Z28tJ05j5Pq0Vg8eeff6JFixbqL12lTZs26uNk+caPHw87Ozu4ubmhT58++Pnnn9XHLl++jKysLHWdFtWmTRtcunQJ2dnZlVlcMpAhdag6Z4OCgjTS+fn5wdvbm+e0hcnKykLt2rVha2uLunXrYsKECUhJSdFIwzq1PKmpqTh9+jRatWoFgOdoVVC8TlV4jlqXgoICKJVK3Lt3D6tXr8bBgwcxbdo0AOY/T61q5e3ySk5ORqNGjbT2e3p6qo+T5XJ3d8fEiRMREhICLy8vXLp0CcuXL0dISAj279+PPn36qOtQVadFeXp6QgiBBw8ewM/Pr7KLT3oypA6Tk5Ph6OgIFxcXnWl5TluOtm3bom3btmjdujUA4Mcff8QHH3yAH374ASdOnICrqysAsE4t0Pjx45GRkYFZs2YB4DlaFRSvU4DnqDUaN24c1q5dCwBwcHDAypUr8dprrwEw/3laLQILAKU22RuzOZ+Mr3379mjfvr36eXBwMMLDwxEUFISoqCj06dNHfYz1bP30rUPWtXV48803NZ6HhYWhffv2GDx4MD799FON46xTyzFnzhx8+eWXWLVqFTp27KhxjOeodSqpTnmOWp+ZM2di9OjRuHv3Lvbu3YsJEyYgIyMDU6ZMUacx13laLbpCeXl56Yy6VM18uqI6smweHh545pln8McffyArKwteXl4AdLc+paSkQJIkeHh4VHIpyRCG1KGXlxeys7ORmZmpMy3PacsWHh4OFxcXjemiWaeWIzo6GosWLcLixYsxYcIE9X6eo9arpDotCc9Ry1a/fn106tQJ/fr1wyeffIJXX30VM2bMwL1798x+nlaLwCIoKAjnz5+HUqnU2H/27FkAUDf/kXURQgCQI+rGjRvDyclJXadFnT17Fk2aNIFCoajsIpIBDKlDVX/Q4mn//vtv3L9/n+e0FRBCwMam8F8Q69QyREdHY/78+Zg/fz5mzpypcYznqHUqrU5Lw3PUejz22GNQKpW4cuWK2c/TahFYhIeHIz09Hbt27dLYv2HDBvj7+6NLly5mKhmV14MHD7Bv3z60a9cOCoUCdnZ26N+/P+Li4pCWlqZOd+PGDSQkJCAiIsKMpSV9GFKHffv2hUKhQExMjEYeMTExkCQJAwcOrKRSU3ns3LkTmZmZePzxx9X7WKfmt3DhQsyfPx+zZ8/GvHnztI7zHLU+ZdVpSXiOWpeEhATY2NigUaNG5j9PDZ6g1kqFhYWJmjVrinXr1on4+HgxZswYAUBs2rTJ3EWjMgwbNkxMmzZN7NixQyQkJIh169aJ5s2bCzs7O/H999+r050/f164urqK7t27iwMHDoi4uDjRunVrLpBnIQ4cOCB27Ngh1q9fLwCIIUOGiB07dogdO3aIjIwMIYRhdaha1GfmzJkiMTFRLF++XDg6OnKhpkpUVp1eu3ZNdO3aVaxcuVIcOHBAfPPNN2L69OlCoVCIVq1aifT0dI38WKfm8+677woAom/fvloLpR09elSdjueo9dCnTnmOWpcxY8aIt956S2zbtk0kJiaKnTt3iqFDhwoAYurUqep05jxPq01gkZaWJv7973+L2rVrCwcHB9GmTRuxZcsWcxeL9LBkyRLRrl074e7uLmxtbYWPj48IDw8Xv/76q1bakydPil69eglnZ2fh5uYmBg4cKC5dumSGUlNxDRo0EAB0blevXlWnM6QOV6xYIZo1ayYcHBxE/fr1xbx580Rubm4lfSIqq05TUlJEeHi4CAgIEE5OTsLBwUE0bdpUREVFiYcPH+rMk3VqHj169CixLovfg+Q5ah30qVOeo9Zl/fr1Ijg4WHh7ews7Ozvh4eEhevToIWJjY7XSmus8lYT4p6M6ERERERFROVWLMRZERERERGRaDCyIiIiIiKjCGFgQEREREVGFMbAgIiIiIqIKY2BBREREREQVxsCCiIiIiIgqjIEFERERERFVGAMLIiIiIiKqMAYWRERERERUYQwsiIiIiIiowhhYEBERERFRhTGwICIiIiKiCmNgQUREREREFcbAgoiIiIiIKoyBBRERERERVRgDCyIiIiIiqjAGFkREREREVGFmDyzS0tIQFRWFJ598Ej4+PpAkCfPnz9eZ9vTp0+jduzdcXV3h4eGBiIgIXLlypXILTEREREREWsweWCQnJ2PdunXIycnBwIEDS0x34cIFhISEIDc3F9u3b8f69etx8eJFBAcH4969e5VXYCIiIiIi0mJn7gI0aNAADx48gCRJuH//Pj777DOd6ebOnQtHR0fs27cPbm5uAICOHTuiadOmePfdd7F06dLKLDYRERERERVh9hYLSZIgSVKpaZRKJfbt24dBgwapgwpADkpCQ0Oxe/duUxeTiIiIiIhKYfbAQh+XL19GVlYW2rRpo3WsTZs2uHTpErKzs81QMiIiIiIiAiygK5Q+kpOTAQCenp5axzw9PSGEwIMHD+Dn56fz9Tk5OcjJyVE/LygoQEpKCry8vMpsLSEiIiIislZCCKSlpcHf3x82NqZtU7CKwEKltCCgtGNLlixBdHS0KYpERERERGTxbt68ibp165r0PawisPDy8gJQ2HJRVEpKCiRJgoeHR4mvnzFjBiZPnqx+npqaivr16+PixYs6W0HI+uTl5SEhIQGhoaGwt7c3d3HICFinVQ/rtOphnVY9rNOqJyUlBc2aNUONGjVM/l5WEVg0btwYTk5OOHv2rNaxs2fPokmTJlAoFCW+3tHREY6Ojlr7PT091UELWbe8vDw4OzvDy8uLfwirCNZp1cM6rXpYp1UP67Tqqozu/1YxeNvOzg79+/dHXFwc0tLS1Ptv3LiBhIQEREREmLF0RERERERkES0W33zzDTIyMtRBw7lz57Bz504AQL9+/eDs7Izo6Gh07twZzzzzDKZPn47s7GzMnTsX3t7eeOutt8xZfCIiIiKias8iAovXX38d169fVz/fsWMHduzYAQC4evUqAgICEBgYiMTEREybNg2DBw+GnZ0devbsiXfffRc+Pj7mKjoREREREcFCAotr167pla5jx444dOiQaQtDREREREQGs4jAgoiIiMic8vPzkZeXZ+5imF1eXh7s7OyQnZ2N/Px8cxeHSmFra2txA+wZWBAREVG1JYTA33//jdTUVAghzF0csxNCoHbt2rh58yYXEbYCjo6O8Pb2hpubm7mLAoCBBREREVVjqampePjwIXx8fODi4lLtL6YLCgqQnp4OV1dXk6/STOUnhEBeXh5SU1Nx69YtALCI4IKBBREREVVLQgjcvXsXbm5u8Pb2NndxLEJBQQFyc3OhUCgYWFg4Jycn1KhRA3/99Rfu379vEYEFf2OIiIioWsrPz0d+fr5FXJARlYckSXB3d0dOTo5FjBFiYEFERETVklKpBCAvxEtkrVQDuC1hsD0DCyIiIqrWqvu4CrJulvT7y8CCiIiIiIgqjIEFERERURUTExMDSZKgUChw/fp1reMhISFo3bq10d93wYIFaNmyJQoKCtT7JEnS2Nzd3RESEoL9+/drvDYgIADPPPOMznxPnjwJSZIQExOjM8+StsTERGRkZCAyMhLNmzdHjRo14OLiglatWmHRokXIyMjQeJ/58+drvN7Z2Rl169ZFnz59sGrVKqSlpWmV7cUXX8TAgQMr9sVVEexUSERERFRF5eTkYPbs2YiNjTX5eyUlJWHZsmWIiYnRmlFq8ODBeOutt1BQUIArV65g0aJF6N+/P/bu3Yunn37a4Pc6evSoxvOFCxciISEB8fHxGvtbtmyJvLw8CCEwefJkNGzYEDY2Nvjpp5+wYMECJCYm4tChQ1r5f/vtt3B3d0dubi6SkpLwww8/ICoqCsuXL8fevXvRtm1bddr58+cjMDAQ8fHx6Nmzp8GfpSphYEFERERURfXt2xebN2/GlClTNC6GTWHFihXw8PBARESE1rFatWrh8ccfBwB07doVTzzxBJo0aYIPP/ywXIGFKi8VHx8f2NjYaO1X2bZtm8bz3r17IycnB8uWLcOVK1fQqFEjjeMdO3bUmII4MjISEyZMQI8ePfDss8/i4sWLcHR0BAA0btwYffv2xTvvvFPtAwt2hSIiIqPIzwcSE4EtW+RHC5ighKjai4qKgpeXF6ZNm1Zm2o8//hghISFo2rQpatSogaCgICxbtkyvaUxzc3Px+eefY/jw4Xqtf9G4cWP4+Pjo7KZVWXx8fADoPytY27ZtMWvWLNy4cUMrUHnxxRdx6NAhXL582ejltCZssSAiogqLiwMmTgT++qtwX926wIoVgI6bl0QWSwggM9PcpdDk7AyUd+KfGjVqYPbs2Zg4cWKZXXUuX76MyMhI1KpVCx4eHjh79iwWL16MCxcuYP369aW+z/Hjx5GcnIzQ0FC9yvXgwQMkJyejadOmGvuFEOppgIsyxlSqQgjk5+cjMzMTv/zyC9577z0MGzYM9evX1zuPZ599FlFRUfjpp58wYsQI9f6QkBAIIXDgwAG88cYbFS6rtWJgQUREFRIXBwweLF+QFXXrlrx/504GF2Q9MjMBV1dzl0JTejrg4lL+148dOxYrVqzAtGnT8Ouvv5Y4Pen777+PgoICPHr0CK6urujRowe8vLwwatQovPfee6hZs2aJ76Ea89ChQwedx1UBgxACly9fxuTJk1FQUIDnn39eI92BAwfU6zIY27Zt2zBs2DD181GjRmHdunUG5dGgQQMA8niSonx9fVGnTh0cOXKkWgcW7ApFRETllp8vt1QUDyqAwn2TJrFbFJE5OTg4YNGiRTh58iS2b99eYrrffvsNAwYMQKNGjWBvbw97e3uMGDEC+fn5uHjxYqnvkZSUBEmSNMYlFLV69WrY29vDwcEBLVq0wC+//IIFCxZg3LhxGum6deuGEydOaG0bN240/IMX06dPH5w4cQLx8fFYvHgxdu3ahUGDBmnMYFUWoeuP3T98fX1x69atCpfTmrHFgoiIyu3wYc3uT8UJAdy8Kaf7178qr1xE5eXsLLcQWBJn54rnERkZiXfffRezZs3SObj6xo0bCA4ORvPmzbFkyRK0aNECzs7O+PXXXzF+/HhkZWWVmn9WVhbs7e1ha2ur8/hzzz2HqVOnQpIk1KhRA40bN9aZ1t3dHZ06dSrfhyxDzZo11XmHhoaicePGiIyMxFdffYXw8HC98lCNCfH399c6plAoyvyeqjoGFkREVG63bxs3XUny8+Xg5PZtwM8PCA4GSrh+IaoQSapYtyNLJUkSli5dirCwMJ3df/bs2YOMjAzs3LkTNWvWhJubG2xsbHDmzBm98vf29kZubi4yMjLgouML9PHxMVnAUF6PPfYYAJTZGlPU119/DUAeU1FcSkoKAgICjFE0q2VwV6jExERMnToVXbt2RYMGDeDr64uWLVti8ODBWL16Ne7du2eKchIRkQXy8zNuOl3i4oCAACA0FBg+XH4MCJD3E5H+evfujbCwMCxYsADpxZplVOMuVFOoAnK3n08//VSvvAMDAwHAqmZFSkhIAAA0adJEr/S///473n77bQQEBOC5557TOKZUKnHz5k20bNnS6OW0Jnq3WGzYsAHvvPMO/vvf/8LV1RVt27ZFp06doFAokJKSgrNnzyIuLg6TJ0/Gc889h4ULF6oHuBARkfUQQm4ZuHsXSE0FHj2SH1U/P3oEZGUBublAdrbcTaO0WXQcHIC33gKUSjukpoZgxgw7CAEUFMhbfn7hz6r9QhTOzvPokXaef/0FDBoEODkBqnGeqq7PpXSBLvHzCiG3gNjby492doCjo7wpFIU/OzoCXl5A/fpAgwaFjw0aAO7uhr0vkTksXboUHTt2xN27d9GqVSv1/rCwMDg4OOD555/HuHHjYGNjg7Vr1+LBgwd65au6g3/s2DG0adPGFEUvt7Vr1+Lw4cN48sknUa9ePWRkZODw4cNYtWoVunbtigEDBmi95tSpU3B3d0deXp56gbzY2Fj4+vpi7969cHBw0Ej/xx9/IDMzU+9ZsaoqvQKLjh074urVq3jhhRewceNGdOzYUeccxSkpKdizZw9iYmLQokULbNy4EYMHDzZ6oYmIyHjS04FTp4Bjx4Djx+XHinZdKio3Fzh9GgAkAMa9+s7KkjdL4OYGdOkCvPkm0Ldv+acHJTKl9u3bY9iwYdi8ebPG/sDAQOzatQuzZ8/GiBEj4OXlheHDh2Py5Ml46qmnysy3Xr16CA4OxldffYVXX33VVMUvl6CgIOzbtw8zZszA/fv3YWdnh6ZNm2LmzJmYPHmyznUs+vbtC0BuwfH09ERQUBCWLl2KUaNGoUaNGlrp9+zZA29vbzz55JMm/zyWTBKlDW//x9y5czFlyhS4ubnpnfFPP/2ElJQUDBw4sCLlM4lHjx7B3d0d9+/fh5eXl7mLQ0aQl5eHAwcOoF+/fiabpo4qF+vUdPLzga+/BtauBS5dAq5c0b7Lb2sLeHvLd+FVm5tb4aOzs9wS4eAg38U/fx7YtQt4+LAwDy8v4JVXgG7d5BYAIZQ4efI4nniiCxwc7GBrC9jYaG6SVPh48iQwZkzZnyc2Vr6gV13IF3/U5eBBYPx47c+tes3y5XK5s7OBnBx5y86Wt/v3gevX5e3GDfkxOVkzn9atgalTgchI+Tuqqqz9PM3OzsbVq1fRsGFDKBQKcxfHIqimm1WNsTDErl27MHToUFy/fh116tQxUQktT35+Ppo0aYLhw4dj8eLFlf7+Zf0eJycnw9vbG6mpqQZdy5eHXi0WCxYsMDjj7t27G/waIiIyrbg44PXX5W5ORXl5yWMXunSRt44dDZ+JZu3a0gdY5+UJ5OXdR0iIgD7XoOfP6/e+trZAsTW2SpWfD7z9dslT5EqSvLDfpEn6DxBPTweuXgU2bJC/hz//BEaOBGbOlFswxoyRAzKiqiwiIgKdO3fGkiVL8NFHH5m7OJVm06ZNSE9Px9SpU81dFLPjOhZERNVEXJw8LqF4UCFJQEoKMGwYMGWKHBCUZ3pLW1sgJETOJySk4rM2mWpguCFT5OrL1RUICgLefVd+7TvvALVry4sETpkC1KtXcjBDVFVIkoRPP/0U/v7+Bq0NYe0KCgrw5ZdfwsPDw9xFMTu9Wix++ukngzJlawURkWXJyABGjNB9THWXftIkYMAAy5nGNTgYqFtXvjjXdUEuSfLx4GDD8jX1FLkeHsC0afL3+eWXcrBx/jwwaxZQp47ckkFUVbVu3RqtW7c2dzEq1ahRo8xdBIuhV2AREhJS4vLvRQkhIEkS8rnEKhFRhRlr7Ybr14Enn5SDi5IUvUuvY3p2s7C1lbskDR4sBxFFgwvVv6QPPzT8O6mMKXIBeezJyy8DL70EzJ8PLFwITJgg12OjRhXLm4jIEukVWMTHx+sVWBARkXHExQETJ2p22albV77Q1rFobokOHZIHEBcfXFwSY84GZQwREcDOnbq/iw8/NOy7UDFVS0hJbGyAefOAxEQ5cHvhBeCnn+QB7UREVYlef9a6d+9u8MwARERUPnFx8l364he9t27J+3fu1O+C+tAhoE8feV2IZs0AfRaXrehdelOIiJC7aBlr5W1TtYSU9Z4bNwJt2wJHjwJLlgBz5hgvfyIiS6BXtODr64sxY8bg22+/RV5enqnLRERUbeXny3fnS5qxCJD77pfV4zQvT+52U1AAPPecvI5E3bolT8EqSfIAY2PdpTc2Yw8MV7WEFJ8Rs25d/QM3QwUEAB9/LP8cHS2vGUJEVJXoFVjMnTsX//vf//DMM8/A19cXL774Ir766itkZ2ebunxERNWKsWYs+vhj4L//BXx8gHXrABcX+S49oB1cmOouvaWLiACuXQMSEoDNm+XHq1dNE1SoPP+83DUtP1/uEpWebrr3IiKqbHoFFv/+97+RmJiIpKQkLFmyBHfu3MGQIUPg4+OD5557Dtu3b0dGaaMCiYhIL8aYsejePXmwMAAsXiwvageY5y69pTN2S0hZJAlYvVpuHbp0SV7jgoioqjBo4ISvry/Gjh2L7777Dn///TdWrFiBjIwMjBgxAj4+PhgwYABiY2PxsOjSq0REpDdjzFg0dy6Qmgq0ayfPSlSUOe7Sk6aaNeXxFpIEfPYZsGePuUtERGQc5R6R7enpiZdffhn79+/H3bt3sW7dOtjY2OC1115DrVq1jFlGIqJqQzVjUXnHQvz+u9z1CZC7Pum6A1/Zd+lJW0gIoFqkd/Roy5uNi4ioPIwy1ZObmxteeOEF7N69G/fu3cOmTZuMkS0RUbWjmrEIMHwshBDywG7VgG2uVWrZFiyQW5WSk4FRo7gqd1WSny9PL7xli/xojuW9YmJiIEkSFAoFrl+/rnU8JCTEJAvZLViwAC1btizXytuSJGHChAlGL5OlU9XVtWvX1Pu6d++OSZMmma1M5WX0OWRdXFwwZMgQY2dLRFRtlHcsRFycfBGjUADLlpm8mFRBjo5ydzSFAjh4EPjhB3OXiIwhLk6eASw0FBg+XH4MCJD3m0NOTg5mz55dKe+VlJSEZcuWYcGCBVymoIIWLlyI1atX47///a+5i2IQvdaxaNiwod4L5EmShMuXL1eoUERE1Z2hazdkZwNTpsg/T50KNGhQeWUlbfqumt6iBTBmDLBqlRwM9u5d+WUl4zHWGjTG1LdvX2zevBlTpkxB27ZtTfpeK1asgIeHByKsfNBWVlYWFAqFWReH7tGjB5o3b4733nsP61T9W62AXuFkjx499N66s+2diMgoDBkL8f778qDsOnWAadMqqYCkk6F3rCdPluv2+++B336rzJKSMRlrDRpji4qKgpeXF6bp8Yfh448/RkhICJo2bYoaNWogKCgIy5Yt02sNs9zcXHz++ecYPny4VmtFSkoKxo0bhzp16sDBwQGNGjXCrFmzkJOTozOvtWvXolmzZnB0dETLli2xdetWjeOZmZmYMmUKGjZsCIVCAU9PT3Tq1AlbtmzRSHfy5Ek8++yz8PT0hEKhQPv27bF9+3aNNKpuSN999x1efvll+Pj4wNnZGdu2bYMkSfhBR1PiJ598AkmS8Mcffxj0XgBw7Ngx/Otf/4JCoYC/vz9mzJhR4vf74osvYvPmzUhLS9N53BLp1WIRExNj4mIQEVF53boFvP22/POyZfKaFWQe5bljHRAgj4nZsgVYvlzuHkXWx5A1aEJCKq1YqFGjBmbPno2JEyciPj4ePXv2LDHt5cuXERkZiVq1asHDwwNnz57F4sWLceHCBaxfv77U9zl+/DiSk5MRGhqqsT87OxuhoaG4fPkyoqOj0aZNGxw+fBhLlizBmTNnsH//fo30X3/9NRISErBgwQK4uLhg9erVGDZsGOzs7DB48GAAwOTJkxEbG4tFixahffv2yMjIwJ9//onk5GR1PgkJCejbty+6dOmCNWvWwN3dHVu3bsXQoUORmZmJl156SeN9X375ZTz99NOIjY1FRkaGeu22L774Ar169dJIGxMTgw4dOqBNmzYGvde5c+fQq1cvBAQEICYmBs7Ozli9ejU2l3DSh4SEYNq0aUhMTET//v1L/f4thqiGUlNTBQBx//59cxeFjCQ3N1fs2bNH5ObmmrsoZCSsU/29+KIQgBBduwpRUGDu0pSsqtepUilE3bpyXejaJEmIevXkdMWdPi2nsbUV4urVSi96uVl7nWZlZYlz586JrKysCue1eXPJdV9027zZCAXXwxdffCEAiBMnToicnBzRqFEj0alTJ1Hwzx+JHj16iFatWmm9Lj8/Xzx48EDk5eWJvLw8sXHjRmFraytSUlJKfb+lS5cKAOLvv//W2L9mzRoBQGzfvl1n+u+++069D4BwcnLSyEOpVIrAwEDRpEkT9b7WrVuLgQMHllqewMBA0b59e5GXl6ex/5lnnhF+fn4iPz9fCFH4PY0YMUIrj8mTJwsnJyfx8OFD9b5z584JAGLVqlUGv9fQoUNL/HwAxNViJ39ubq6QJElMmzat1M9a1u/x/fv3BQCRmppaaj7GUK6RNSkpKZgzZw66du2K5s2b41//+hfmzZuHBw8eGCfaISIivfz6KxAbK/+8YkXJ09SS6VVk1fT27YGwMLmbzAcfmK6MZDrGWIPGVBwcHLBo0SKcPHlSZ/ccld9++w0DBgxAo0aNYG9vD3t7e4wYMQL5+fm4ePFiqe+RlJQESZLg7e2tsT8+Ph4uLi7q1gYV1V384l2NevXqpbFsga2tLYYOHYpLly7hr39OsMceewzffPMNpk+fjsTERGRlZWnkcenSJVy4cAHPP/88AECpVKq3fv364fbt21qDogcNGqT1mV5++WVkZWVh27Zt6n1ffPEFHB0dMXz4cIPfKyEhocTPp4u9vT08PDxw69YtncctkcGBxa1bt9ChQwcsXrwYqampqF+/Ph4+fIiFCxeiQ4cOSEpKMkU5iYioGCEKB2yPHAl06mTe8lR3FV01PSpKfvzsM3kKWrIuFV2DxtQiIyPRoUMHzJo1S2ef/hs3biA4OBhJSUlYsmQJfvzxR5w4cQIff/wxAGhdvBeXlZUFe3t72BYbDJacnIzatWtrDYT29fWFnZ2dRvclAKhdu7ZW3qp9qrQrV67EtGnTsGfPHoSGhsLT0xMDBw7E//73PwDAnTt3AABTpkxRB0iqbdy4cQCA+/fva7yHn46Ir1WrVujcuTO++OILAEB+fj42bdqEAQMGwNPT0+D3Un0XJX0+XRQKRZnfvSXRa4xFUTNnzkRWVhaOHz+Ozp07q/efOHEC/fv3x8yZMzkmg4ioEuzfL9/9ViiARYvMXRqq6B3rXr3klovffgNWrwbmzDFe2cj0VGvQDB4sBxFFx9mUtQZNZZAkCUuXLkVYWJjOWYb27NmDjIwM7Ny5EzVr1oSbmxtsbGxw5swZvfL39vZGbm4uMjIy4FJkoJeXlxeOHz8OIYRGcHH37l0olUqtFo6///5bK2/VPi8vLwDy0gbR0dGIjo7GnTt31K0X/fv3x4ULF9R5zpgxo8QZqpo3b67xvKQZoEaNGoVx48bh/PnzuHLlCm7fvo1Ro0ZpfG5938vLy6vUz6fLgwcPtL4jS2Zwi8W3336LRYsWaQQVANC5c2csWLAA33zzjdEKR0REuuXnA9Onyz9PnCjfKSXzqugda0kqbLVYuRKwopuU9I/yrkFTWXr37o2wsDAsWLAA6enpGsdUF9aOjo7qfUIIfPrpp3rlHRgYCABaSw706tUL6enp2LNnj8b+jRs3qo8X9cMPP6hbAQC5lWDbtm1o3Lgx6ur4Q1erVi289NJLGDZsGP773/8iMzMTzZs3R9OmTfH777+jU6dOOrcaNWro9bmGDRsGhUKBmJgYxMTEoE6dOnjyySfVxw15r9DQ0BI/ny5JSUnIzs5Gy5Yt9SqrJTC4xSI1NRUBAQE6jzVs2BCpqakVLRMREZVhwwbgP/8BatYsDDDIvIxxx3rwYGDGDHnq4JgY4PXXTVhgMglD16CpbEuXLkXHjh1x9+5dtGrVSr0/LCwMDg4OeP755zFu3DjY2Nhg7dq1eo+fDflnqqtjx46pZ0sCgBEjRuDjjz/GyJEjce3aNQQFBeHnn3/G22+/jX79+qF3scVbvL290bNnT8yZM0c9K9SFCxc0ppzt0qULnnnmGbRp0wY1a9bE+fPnERsbiyeeeALOzs4A5Clrn3rqKfTp0wcvvfQS6tSpg5SUFJw/fx6nT5/Gjh079PpcHh4eCA8PR0xMDB4+fIgpU6ZoTaer73vNnj0bX3/9NXr27Im5c+fC2dkZH3/8MTIyMnS+97FjxwBAa6Yti2boaO/AwEAxceJEncfefPNNERgYWLHh5JWAs0JVPdY+MwlpY52WLCNDiDp15Blm3n3X3KXRX3Wp0127tGeHqldP3q+PVavk1zRqpHsGKUti7XVqzFmhLE3RWaGKGz58uACgNSvU3r17Rdu2bYVCoRB16tQRU6dOFd98840AIBISEsp8z+DgYNGvXz+t/cnJyWLs2LHCz89P2NnZiQYNGogZM2aI7OxsjXQAxPjx48Xq1atF48aNhb29vQgMDBRffvmlRrrp06eLTp06iZo1awpHR0fRqFEj8eabb2pd1/3+++/iueeeE76+vsLe3l7Url1b9OzZU6xZs0av70nlu+++EwAEAHHx4kWdafR5LyGEOHLkiHj88ceFo6OjqF27tpg6dapYt26dzlmhXnzxRREUFFRiuVQsaVYogwOLpUuXCkmSxBtvvCFOnjwpbt26JU6ePCkmTZokbG1txfLly01RTqNiYFH1WPs/N9LGOi3ZkiXyhWf9+kJY0/VQdapTpVKIhAR5atGEBMMChPR0Iby85DouNkOnxbH2Oq3KgUV5qaabVU2RaoidO3cKW1tb8ddff5mgZNVLamqqcHFxEevWrSszrSUFFgZ3hZo6dSouX76Mjz76SD1TwD8tH3j11VcxRTVFCRERGV1yMvDOO/LPixbJA7fJ8qhWTS8PFxdgwgQgOhpYurSwaxWRpYuIiEDnzp2xZMkSfPTRR+YujlX74IMPUL9+fY2B4tbA4MBCkiSsXbsWkydPRkJCApKTk+Hl5YWePXuiWbNmpigjERH9Y/FiIDUVaNsW+GfadKqCxo+XV1E/dQpITASsqYs1VV+SJOHTTz/F119/jYKCAq2xCKQ/Nzc3xMTEwM7O4Et1syp3aZs3b641VRcREZnOtWuAqqF46VKA/7OrLh8f4OWX5fpetoyBBVmP1q1bo3Xr1uYuhtV78803zV2Ecil3YJGeno4bN24gOztb61iHDh0qVCgiItI2Zw6Qmyuvd1BktkOqoiZPBj75BPj2W+Cjj+TuUURElszgwOLevXsYM2YM9u7dq3VM/LP4SX5+vlEKR0REsjNngC+/lH9eupR97quDRo2AefPk7Y03AFdX4KWXzF0qIqKSGRxYvPbaa4iPj8fEiRPRokULODg4mKJcWhITE0ucx/fo0aN4/PHHK6UcRETmMG2aPHFpZCTQsaO5S0OVZc4c4MEDef2LV16RB3YPGWLuUlU9ouiiI0RWxpJ+fw0OLOLj4/Hee+9hzJgxpihPmd5++22tAIN9+YioKtu8GfjuO8DeXh68TdWHJAHvvw+kpwOffQYMGwZcvgw0aGB5C69ZI3t7ewBAZmYmnJyczFwayyOE/LuXmws4OMitZsZoLTVVvtVVRkYGJElS/z6bk8GBhYuLCxo0aGCKsuiladOmbJ0gomrj88+BV1+Vf544Ue4eQ9WLJAFr1gAXLgA//yyvzK1St6682ndEhPnKZ81sbW3h4eGBu3fvAgCcnZ0hWeEVrhBAZmbhhbqzc/kv1AsKCpCbm4u7d3Nw546EvLzCY/b2QO3agLt7+cuamgr8/TeMnq+KMb8LS85XCAGlUolHjx7h0aNH8PDwgK0F3GUwOLB48cUXsWPHDjzJkYNERCb17rvA1Knyz2PGFK5fQdXPV18BR45o7791S17nYudOBhflVbt2bQBQBxfWJjMTSEkBig5vtbUFPD3li1RDCSHw8GEuHj1y1Hn89m151rLy5J2ZCdy7p/tYRfItmr8xvwtryNfW1hZ+fn5wLxaV5ecDhw/L36uzc+UFywYHFosWLcIrr7yC8PBwPP300/D09NRKE2HCv27jx49HZGQknJ2d8cQTT2DOnDno1q2byd6PiKiyCQHMmgUsWSI/nzZN/tkKb6SSEeTny61VurpRCyH/XkyaBAwYUL26RRW9cKpItzBJkuDn5wdfX1/kFb2NbgW++07374bqb8WKFYbPIJeTk4fnn8/HgweOALT/6EiS3Lpw6JBh33d+vjyj3d9/6z5e3nxVTPFdWHq+dnZ2sLW11Wpli4uT8/7rL3VKwwtYTga/09WrV3H8+HFcvHgRX331ldZxU80K5e7ujokTJyIkJAReXl64dOkSli9fjpCQEOzfvx99+vQp8bU5OTnIyclRP3/06BEAIC8vz+r+iJBuqnpkfVYd1bVO5YtIG6xbJ/9nXbw4H1OnFkCpNHPBjKC61mlF/fijhL/+KvnftRDAzZtAQoISPXqUfxBnfj7w88+S+kK9WzdR5gWePnVannzLsnu3hMmTbXHrVuEFVZ06Au+/n4/w8PJ/B7a2tiYprynk5wPjxtnh1i1AdwAgMH488L//KQ0q/9GjwJkzNUpNc+0acPSoYb9vP/8s4fjx0i87y5MvYLrvwpz5jhsHnDqlREGB3E0qL09+VP2sVMrd1PLy8v95lLcjRyS8/75qkaPKvxslCQOHkoeGhuLChQuYMWNGibNC9ejRw2gFLM3Dhw8RFBQET09P/P777yWmmz9/PqKjo7X2b968Gc4VacMiIjKivDwJK1d2wOHDdSFJAq+//juefPK6uYtFZvbTT3Xw/vudykzXu/c1jB37B+zsDL+wPnrUD599FoTk5MIBzF5eWRg9+iyeeOK2wfmZMt+jR/2wdGnnf54VvXCSP/e0aScqlLcpvgeV/Hzg3DkvPHigQM2a2WjZMrncQcvZs16YM6fsHhsLF/6MoKBkvfPV9/dt8uST6N79ltnzBcr/XQihukC3RW6uDfLybJGXZ6P++b//rYnPPmtTZr7du99EzZo5yMuzQV6eDZRKeVP9XHxfero97t51MegzVswjAO5ITU2Fm5ubSd/J4MDC2dkZn3/+OYYNG2aqMhnk9ddfx5o1a0qd0UFXi0W9evVw+/ZteHl5VVZRyYTy8vLw/fffIywszCJmRaCKq251evEiMGWKLb791gb29gIxMfkYMsRyphA0hupWp8by448SwsL062DQuLHA3Ln5GDpU6L0y++7dEiIjbf/pklF4oS5J8u/f1q0ltwKUVqcVybck+flAkyal3+mtU8fwO8imKm/x/I3ZyrJ1q4QRI8r+vdi4UYnISP3z/+GHfDz1lKLMdN9/b1jLgr6/x/v3K9Gxo0BmJtRbVpb0zyPUj1lZkvr5qVMS9u0r+xfe11fAzg7Izi7cCgosv4+pra2Ag4M8yN3eHho/29kV/pyVJXDhgq7vofICC4O7QtWqVQseHh4mKEr5qOKi0mZxcHR0hKOj9iAke3t7/nOrYlinVU9VrdP8fOD4cXlQ7tdfyzP+AICTExAXJ6Fv38rrE1vZqmqdmkpoqDz7061busdZSJI8m469PXD5soSRI+3w7rvy1MTPPFP62Jz8fOCtt0oavyFBkoApU+wwaFDp/d6L16mx8i3uyBH8E1ToJoSEv/4Cjh2zR0iI/vmaqrwqcXHyGjTF809KkhAZaVeuwff16umbzg6GnG4hIXIrTUqKAkLoHmNRpw7Qrp0d/v4bSEuTp45NTwcyMgp/LrovI0NO5+QkBwWlefpp0/3tu3u39CDC0VHeFAp5y88v/fdNJSICaNxYvuBX5aHaVPscHAp/PndO/n0ry8GDQO/egI2NfsHPli0Shg/XK6nJGFx7r7/+OtauXYunnnrKFOUxyIMHD7Bv3z60a9cOCkXZ0TURkTndvSv3X/76a2DfPvm5ip2dfAG5cCHQpYv5ykiWx9ZWHsg5eLB8UVf04lQVNHz+uTzQc9UqYNky4OxZ4Nln5cGyH3wABAXpzvvw4aIDPLWpxm8cPgyDLtRNle9tPXsj6ZtOxVTlBUw3+D44uOyAs25dOZ1Kbq483WtqKvDwYeHPqanAo0fy9vChDerWTUNysu7rKiHk78rHR/+yloe9vTwzUtHNyUneij93dAQ2bJCDl5L4+MgBnotLYeDg5CQ/qoKA4q18+flAQEDZ3/H27YbVXViYfF6WlW+vXtplKo2fn/5pTcXgwMLGxgZ//PEHOnTogH79+mnNCiVJEt58802jFVBl+PDhqF+/Pjp16gRvb2/873//w3vvvYc7d+4gJibG6O9HRKQvXbPTPHgAnDolbydPytvNm5qvc3cH+vWTLyj69jXOHO5UNUVEyFPKas70Il98fPhh4d3uGTOAsWPl4OKDD4AffgDatZP3RUcD3t6a+ZrqQt1U+ep74WToBZapyguYJmjJygLu35cDkilTSs7XxQVo21b+e/TwodxtqGy2AHz1KodqcbsaNeT3cnXV3FxcCjfV8wsXgK1b5WlWVWrXBubOBcLDCwMGQxs1Q0Pl4Fv12VVUwfeaNYChk4jqE9R/+KHhrVimyresYLMyGBxYREVFqX8+c+aM1nFTBRZt2rTBtm3bsGbNGqSnp8PT0xPdunVDbGwsOnfuXHYGREQmEBcH/Pvfms3ltraa85KrSBLQrJl8Z/nZZ4Hu3eV/zET6iIiQg9CyplitWVOenvjVV4GoKDkgWb1aXsF9/nxg3LjCizZTXaibKt/y3KU3ZjnKc0dY32AkKUkOAG7eBG7ckIORv/8G7tzR3tLT9ctT1cWyuBo15BsZujYXl3zcvn0RnTo1g7u7Lf76Sx6LULeufGHu4VEYTJT379fKlcaZKrgofYPvqpxvaQFLZTF48Pb162XPUGLOlbn18ejRI7i7u+P+/fscvF1F5OXl4cCBA+jXrx/7blcRllynd+/KKyBv2CB3ayqJn598B7JTJ6BjR6B9e8DE4+YsmiXXaVX244/yxYtq8sTAQOD99+VWsoIC/bp6XL2q+8KvpDrVtwtJSfmWJi6u9DvT5RmvYMryJibKd9PLolDIF/D6srMDvLzkVihPT7ncqhWs27SRj9WsWbh5eMiPbm6lfwZrP0+Ntb6JNeervY6FBQ/etvSggYjI2G7fli/OVNv582W/RpLkf/yxsdVr0TKyPD16yF3y1q+XF168cEHugvevf8ktGB9+CAwZYh1dPYDKv9Nb3vIKIY91OXlS7t5TVjckVVDh7S0PzK5bV77QrFVLe/P1lVsXuGimNltbw8fBVLV8i7duOjsrMXCgcd+jJFV32hEioiIMuSt07x4QHy/3T//xR3kq2OIaNpTvXpakIoM9iYzN1hYYMwZ47jlg0SLgo4/kGZbCwuQAY84cOfCw9K4eRfPWp1tYZZc3Lw/46Se5JfPrr+UF3/Qxbx4wfLj8Xlxei4yhaMCSnFx5faL0Cixat26NhQsXIjw8XK9Mb9++jSVLlsDf3x/Tp0+vUAGJiCpKu1lY/ge+YoV8sZCZKV+gHDokb8WHj0mSPAAyJES++xscDHz3HfSa1q88gz2JTMXdHVi+HJg8WR7gvWaNHGAcOQJ07SoPBvbzA/z9jdfv3dgBgEpl3OnVp7w5OcDu3XIgceCAPMOSikIhz+zTqZM8LmLTJnmMhEq9ehUPsogsiV6BxXPPPYcRI0agZs2aeP755xESEoIOHTrA29sbkiQhKysLly9fxrFjx/DVV1/h4MGD6Ny5M8aOHWvq8hMRlUrVH7t4v+lbt4BBg4BWrYD//U+ehrGoNm3kC4LQUHnAYs2amsdNOdiTyNT8/ORZo6KigKVL5QDjl1/k7bHH5Atdb2/5/KhodxtTdSExFUPKe/EiMHSo5s0IHx+gf385QOndW7MFYulS0wRZRJZCr8Bi7ty5GDNmDD788EN89tlnWLp0KSRJgiRJsLe3R+4//5GFEAgODsbWrVsRwfCbiMysrPnjAeA//5Ef69WTLwJ695YDilq1Ss/bVLPTEFUmPz/5jrkqwFi7Fvj1V3mbPh2oX18ej9GvH9CzpzxdKMliY4HXX5cXf/PyAkaPlmd769Kl5GDB2oIsIkPpPcbCz88PS5cuxaJFi3D8+HEcPXoUSUlJyMrKgre3NwIDAxESEoK6deuasrxERHrbsqX0+eNVNm4EXnjBsDuzphycSlTZ/P3l3+fp04Fdu+QuPQkJ8rSna9bIm6Mj8MQT8pTJjRsDjRrJgUdGRvUarpmeDkyYIM8KB8iBwqZN8mrURNWdwX8N7O3t0a1bN3QzdJURIiITy88Hjh8vHDipz+xNgDx7U3m6e5hycCqROfj5yRfNEybIY48SEoD9++Xtxg156tTExKKvsAfwNN54Q6B+/cJpTVWb6rmbW+Eqx6rVkouufOzgULjZ2xc+WtqsR2fOyF2fLl6UV0SePx+YOZM3EIhUqtdtBiKyeKrZm27elHD9uhf69Cl9BdbMTOD774GvvgL27ZNndFKxsZHn6S9LRcZBmHJwKpE5OTsDTz8tb0IA584BJ04AV64UbpcvC9y9KyElRdJYSdlY7Ozkzd6+8Oeiz21tCx+L/1zSZmNT+Fh8s7eXgyEvL+3t8GFg6lR5sHadOvKCg927G/8zE1kzBhZEZDE0Z2+yA9ANa9YIrFypefc/JUUOIvbsAQ4e1Jwf3sND7g/+7LPyVJpt25p+HAT7TVNVJ0nyQO5WrTT35+UpsWvXQTRr1gd37tjj4UPo3FJT5XUasrIKH4v+nJenPYECACiV8mbIwnGm1r8/8MUXcrBBRJoYWBCRRShp9qakJHn/2rWF0zr++KPcsqHSoAEwcKDcctCtm2YLB8dBEJmWk1M+goKADh0qlo8QchChCjJUmyq4UB0r+nN+fuGmVGo+L20rKNC95ebKNy6Sk7U3AJgxA/j3vy2vixaRpWBgQURmV/rsTfJ/8Fdf1dwfFASEh8sBRbt2Jf+j5zgIIusgSfJNAXt7LhJHZK0YWBCR2R0+rN/sTa1aAaNGycFE48b6589xEERERKbHwIKIykU1yNoYF+q//aZfulmzgGHDyvceHAdBRERkWuUOLFJTU3Hx4kVkZWVpHevOaRKIqjTNQdayunXl8Qz6di26cgXYvl3e9A0suIo1ERGR5TI4sFAqlRg7diw2btyI/KKjJ4soaT8RWb+SBlnfuiXv37mz5ODi+vXCYOLkycL9qmkec3J0v46rWBMREVk+G0Nf8MEHH2Dv3r1Yv349hBD46KOPsHbtWnTq1AlNmzbFN998Y4pyEpEFKH2Qtfw4aZLmjE1JSXJLRteuQEAAEBUlBxU2NkDv3sC6dcCdO/Kc8JKkPQhbkuSMOXsTERGRZTM4sIiNjcWsWbMw7J+Ozl26dMHo0aNx/PhxNGjQAAkJCUYvJBFZhrIGWQsB3LwJ7N0LfPKJPKahbl052Dh6VA4aQkLkY7dvywvbjRkDeHsXzt5Up45mnnXqlN4KQkRERJbB4K5QV65cQdu2bWFjI8ck2UVWrRk7diwmTpyIJUuWGK+ERGQxbt/WL11EhGarRteuwNChclcpf//SX6eavenmTSWuXz+GKVO6QKEoZeltIiIisggGBxYuLi7Izc2FJEnw9PTE9evX0bVrVwCAk5MTklWryBBRlaPv4GkhgI4dgchI4LnngPr19X8P1exNeXkCBw4ks/sTERGRlTA4sAgMDMTVq1cBAF27dsX777+P4OBgODg4YNmyZWjevLnRC0lE5WesaWEfPpRfV7Mm8OBByenc3IDjx4HAwHIXmYiIiKyQwYHF0KFDcfHiRQBAdHQ0unfvjgYNGgAA7O3tERcXZ9wSElG5lWdaWKUSuHBBngL2zBngjz+Ac+fkQdilUQ26/uILBhVERETVkcGBxbhx49Q/t2/fHufOncOePXsgSRLCwsLYYkFkIfSdFvbiReDQocJA4uzZkqd9rVMHaNkScHQEjhzRbLmoW1eeuYmDrImIiKqnCq+8Xa9ePbzxxhvGKAsRGUlZ08JKEvDKK8A77wAnTmincXUF2rWTt7ZtgVat5IDC3V3zPYy18jYRERFZvwoHFkRkefSZFvbhQzmosLUFQkOBxx6TA4n27YFGjeR1JkqjGmRNREREBOgZWLRp00bvDCVJwu+//17uAhFRxek7LeyIEcCyZUCtWqYtDxEREVV9egUWnp6ekIovh1tMeno6Tp06VWY6IjK9vDz90o0axaCCiIiIjEOvwCIxMbHEY0qlEuvWrcOCBQsgSRKGDx9urLIRVRvGGK+QnAxs3Qps3Aj8+mvpaSVJHmwdHFz+MhMREREVVUYv6tLt2LEDLVu2xBtvvIG2bdvi1KlTiI2NNVbZiKqFuDggIEAe5zB8uPwYECDvL4tSCezbJ8/y5OcHTJggBxW2tvICdUDhNLAqqucffsjB1kRERGQ85QosEhMT0aVLFwwdOhRubm747rvvcPDgQbRr187IxSOq2lRTwhYfaK2aErak4OI//wGmTpVbHfr3B3btkrs/tW8vBwxJScDJk/L+OnU0X1u3buFUs0RERETGYtCsUGfPnsW0adNw8OBBNGzYEJs3b0ZkZKSpykZUpekzJeykScCAAXLLwoMHwJYt8gJ0J08WpvX1BV54ARg5Eig+z0JEhPx6TgtLREREpqZXYHHz5k3Mnj0bmzdvhqenJz788EOMHTsW9vb2pi4fUZWlz5SwN2/KLRCnT8utF9nZ8jE7O7ml4qWXgKeeAko7FTktLBEREVUGvQKLZs2aITc3F3379kVUVBRq1KiBs2fPlpi+Q4cORisgUVWl75SwU6YU/hwUJC9sN3w44ONjmnIRERERlYdegUVOTg4A4JtvvsG3335bYjohBCRJQn5+vnFKR1SF+fnpl87FRV5v4uWX5QHZnNGZiIiILJFegcUXX3xh6nIQVTvBwfJA6lu3dI+zAABPT+D6dcDVtXLLRkRERGQovQKLkSNHmrocRFbBGOtNAPIMTl9/DXh56R5noWqV+PRTBhVERERkHQyaFYqoOouLk2dxKhoI1K0LrFih/9StN28C69bJAcOdO4X7HR2Bf3ocqvP98ENOCUtERETWg4EFkR5U600U77KkWm+itHUhCgqAQ4eA1auBvXvl5wBQu7Y8bmL0aKB+fU4JS0RERNaNgQVRGQxdb0IlJUVec2LNGuDSpcL9oaHAuHFy+qLTxHJKWCIiIrJmDCyIyqDvehOHD8vBwZ9/AitXAps2AVlZcho3N3nNibFjgRYtKqPURERERJWLgQVRGfRdb2LvXmDhQiA+vnBf27bA+PHyuhMuLqYpHxEREZElYGBBVAZ915t4/3350cYGCA+Xu09168Z1J4iIiKh6YGBBVAZ91psAAA8P4NVX5RaK+vUrrXhEREREFsHG3AUgsnS2tvKUsqV57TU58Fi6lEEFERERVU8MLIjKoFTKW5Mm2sd8fOSpZtesAZydK79sRERERJaCXaGoSjLGCtlpacBnn8mtFdevy/scHIDeveXZnzp35noTRERERCoMLKjKqegK2X/9JU8Xu24dkJoq7/PxkcdOvP464OtrmnITERERWTMGFqUwxl1vqly7d0uIjCzfCtlnzgDvvQds3Sp3fQKAwEBg8mTghRcAJyeTFp2IiIjIqjGwKEFF73pT5cvPByZPtjVohWylUl5/YtUqICGhMH1ICPDWW0C/fvL0sURERERUOgYWOsTFyXe3y3PXm8zn3Dkv3LpV8qIRRVfIbtVKHj/xySfyPkAONp57Tg4oOnaspEITERERVREMLIrJz5dbKgy5603lZ8zuZg8eKPRKN3cucPw4kJsrP/f2BsaMAcaO5VSxREREROVlVZ080tPTMWnSJPj7+0OhUKBdu3bYunWrUd/j8GHN7k/FFb3rXV75+UBiIrBli/yYn1/+vKw537g4ICAACA0Fhg+XHwMC5P3lUbNmtl7pDh+Wg4pOnYANG+T6fPttBhVEREREFWFVLRYRERE4ceIE3nnnHTRr1gybN2/GsGHDUFBQgOHDhxvlPW7fNm664kw1dsMa8zV2d7OWLZNRp45AUpJU6grZw4fLn+mxxwwvNxERERHpZjUtFgcOHMD333+P1atX47XXXkNoaCg+/fRThIWFYerUqcg30m10Pz/jpitKdTFdvEVEdTFd3jv11pZvWd3NALm7mSFVmp0NXLpUE717F5QaVKxfD3z5JYMKIiIiImOzmsBi9+7dcHV1xZAhQzT2jxo1CklJSTh+/LhR3ic4WL4jL5UwBliSgHr15HSGMMXFtDXmC1S8u5lSCZw9KwcJY8fKA609Pe0wbVp3bNige4BGvXrArl3AqFGGl5eIiIiIymY1XaH+/PNPtGjRAnZ2mkVu06aN+njXrl0NynPvXglubtr7hw0Dli/X/RohgMhI4KuvDHor/PmnfhfTixcDrVtXnXyLByZCAD//rF8ZVq0Ctm8H7t2Tt7t35cfkZF0BjwR39xz861/2CAmxQUSEXD6uQUJERERUOawmsEhOTkajRo209nt6eqqPlyQnJwc5OTnq56n/LKc8alQmyvMVlBR0GMO8ecxXpbSuVi4uAu3ayVuHDgJt2uTh4sV49OwZCnt7ewCAhwcQFCSnf/jQdOUk08jLy0NmZiaSk5PVdUrWjXVa9bBOqx7WadWTkpICABCl9RU3EqsJLABAKql/UhnHlixZgujoaB1H6hmhVGQOGRnAkSPyRkRERESlS05Ohru7u0nfw2oCCy8vL52tEqooTNVyocuMGTMwefJk9fOHDx+iQYMGuHHjhsm/YKocjx49Qr169XDz5k246erfRlaHdVr1sE6rHtZp1cM6rXpSU1NRv379Uq+VjcVqAougoCBs2bIFSqVSY5zF2bNnAQCtS+no7+joCEdHR6397u7uPGmqGDc3N9ZpFcM6rXpYp1UP67TqYZ1WPTY2pp+zyWpmhQoPD0d6ejp27dqlsX/Dhg3w9/dHly5dzFQyIiIiIiKymhaLp556CmFhYXj99dfx6NEjNGnSBFu2bMG3336LTZs2wZZT/hARERERmY3VBBYAEBcXh1mzZmHu3LlISUlBYGAgtmzZgsjISIPycXR0xLx583R2jyLrxDqtelinVQ/rtOphnVY9rNOqpzLrVBKVMfcUERERERFVaVYzxoKIiIiIiCwXAwsiIiIiIqqwahNYpKenY9KkSfD394dCoUC7du2wdetWcxeL9JCYmAhJknRux44d00h7+vRp9O7dG66urvDw8EBERASuXLlippKTSlpaGqKiovDkk0/Cx8cHkiRh/vz5OtMaUoerVq1CYGAgHB0d0bBhQ0RHRyMvL8+En4RU9K3Tl156See5GxgYqDNf1ql5xMfH4+WXX0ZgYCBcXFxQp04dDBgwAKdOndJKy3PUOuhbpzxHrceZM2fw9NNPo379+nBycoKnpyeeeOIJbNq0SSut2c5TUU2EhYUJDw8PsWbNGhEfHy9Gjx4tAIgvv/zS3EWjMiQkJAgA4u233xZHjx7V2NLS0tTpzp8/L2rUqCGCg4PF/v37xa5du0SrVq2Ev7+/uHv3rhk/AV29elW4u7uL7t27q8+9efPmaaUzpA4XLVokJEkSM2bMEAkJCWLZsmXCwcFBjBkzppI+VfWmb52OHDlSODk5aZ27Z86c0UrLOjWfwYMHi9DQULF69WqRmJgoduzYIR5//HFhZ2cnfvjhB3U6nqPWQ9865TlqPRISEsRrr70mYmNjRXx8vNi7d6+IjIwUAMTChQvV6cx5nlaLwGL//v0CgNi8ebPG/rCwMOHv7y+USqWZSkb6UAUWO3bsKDXdkCFDhLe3t0hNTVXvu3btmrC3txdRUVGmLiaVoqCgQBQUFAghhLh3716JF6H61uH9+/eFQqEQr776qsbrFy9eLCRJEv/5z39M80FITd86HTlypHBxcSkzP9aped25c0drX1pamqhVq5bo1auXeh/PUeuhb53yHLV+Xbp0EfXq1VM/N+d5Wi26Qu3evRuurq4YMmSIxv5Ro0YhKSkJx48fN1PJyFiUSiX27duHQYMGaawU2qBBA4SGhmL37t1mLB2pmtZLY0gdfvvtt8jOzsaoUaM08hg1ahSEENizZ49Ry0/a9KlTQ7BOzcvX11drn6urK1q2bImbN28C4DlqbfSpU0OwTi2Xt7c37OzkFSTMfZ5Wi8Dizz//RIsWLdRfukqbNm3Ux8nyjR8/HnZ2dnBzc0OfPn3w888/q49dvnwZWVlZ6jotqk2bNrh06RKys7Mrs7hkIEPqUHXOBgUFaaTz8/ODt7c3z2kLk5WVhdq1a8PW1hZ169bFhAkTkJKSopGGdWp5UlNTcfr0abRq1QoAz9GqoHidqvActS4FBQVQKpW4d+8eVq9ejYMHD2LatGkAzH+eWtUCeeWVnJyMRo0aae339PRUHyfL5e7ujokTJyIkJAReXl64dOkSli9fjpCQEOzfvx99+vRR16GqTovy9PSEEAIPHjyAn59fZRef9GRIHSYnJ8PR0REuLi460/Kcthxt27ZF27Zt0bp1awDAjz/+iA8++AA//PADTpw4AVdXVwBgnVqg8ePHIyMjA7NmzQLAc7QqKF6nAM9RazRu3DisXbsWAODg4ICVK1fitddeA2D+87RaBBYASm2yN2ZzPhlf+/bt0b59e/Xz4OBghIeHIygoCFFRUejTp4/6GOvZ+ulbh6xr6/Dmm29qPA8LC0P79u0xePBgfPrppxrHWaeWY86cOfjyyy+xatUqdOzYUeMYz1HrVFKd8hy1PjNnzsTo0aNx9+5d7N27FxMmTEBGRgamTJmiTmOu87RadIXy8vLSGXWpmvl0RXVk2Tw8PPDMM8/gjz/+QFZWFry8vADobn1KSUmBJEnw8PCo5FKSIQypQy8vL2RnZyMzM1NnWp7Tli08PBwuLi4a00WzTi1HdHQ0Fi1ahMWLF2PChAnq/TxHrVdJdVoSnqOWrX79+ujUqRP69euHTz75BK+++ipmzJiBe/fumf08rRaBRVBQEM6fPw+lUqmx/+zZswCgbv4j6yKEACBH1I0bN4aTk5O6Tos6e/YsmjRpAoVCUdlFJAMYUoeq/qDF0/7999+4f/8+z2krIISAjU3hvyDWqWWIjo7G/PnzMX/+fMycOVPjGM9R61RanZaG56j1eOyxx6BUKnHlyhWzn6fVIrAIDw9Heno6du3apbF/w4YN8Pf3R5cuXcxUMiqvBw8eYN++fWjXrh0UCgXs7OzQv39/xMXFIS0tTZ3uxo0bSEhIQEREhBlLS/owpA779u0LhUKBmJgYjTxiYmIgSRIGDhxYSaWm8ti5cycyMzPx+OOPq/exTs1v4cKFmD9/PmbPno158+ZpHec5an3KqtOS8By1LgkJCbCxsUGjRo3Mf54aPEGtlQoLCxM1a9YU69atE/Hx8WLMmDECgNi0aZO5i0ZlGDZsmJg2bZrYsWOHSEhIEOvWrRPNmzcXdnZ24vvvv1enO3/+vHB1dRXdu3cXBw4cEHFxcaJ169ZcIM9CHDhwQOzYsUOsX79eABBDhgwRO3bsEDt27BAZGRlCCMPqULWoz8yZM0ViYqJYvny5cHR05EJNlaisOr127Zro2rWrWLlypThw4ID45ptvxPTp04VCoRCtWrUS6enpGvmxTs3n3XffFQBE3759tRZKO3r0qDodz1HroU+d8hy1LmPGjBFvvfWW2LZtm0hMTBQ7d+4UQ4cOFQDE1KlT1enMeZ5Wm8AiLS1N/Pvf/xa1a9cWDg4Ook2bNmLLli3mLhbpYcmSJaJdu3bC3d1d2NraCh8fHxEeHi5+/fVXrbQnT54UvXr1Es7OzsLNzU0MHDhQXLp0yQylpuIaNGggAOjcrl69qk5nSB2uWLFCNGvWTDg4OIj69euLefPmidzc3Er6RFRWnaakpIjw8HAREBAgnJychIODg2jatKmIiooSDx8+1Jkn69Q8evToUWJdFr8HyXPUOuhTpzxHrcv69etFcHCw8Pb2FnZ2dsLDw0P06NFDxMbGaqU113kqCfFPR3UiIiIiIqJyqhZjLIiIiIiIyLQYWBARERERUYUxsCAiIiIiogpjYEFERERERBXGwIKIiIiIiCqMgQUREREREVUYAwsiIiIiIqowBhZERERERFRhDCyIiIiIiKjCGFgQEVm5X375BfPnz8fDhw+1joWEhCAkJKTSy1SWjRs3wsfHB2lpaUbL86WXXkJAQEC5Xvv555+jTp06yMjIMFp5iIiqG0kIIcxdCCIiKr93330XU6dOxdWrV7UurM+dOwcAaNmypRlKpltmZiaaNWuGSZMmYcqUKUbL9/Lly3j06BHat29v8GuVSiVatmyJYcOGITo62mhlIiKqTthiQURUhbVs2dKiggoA2LBhA5KTkzF69Gij5tu4ceNyBRUAYGdnh9deew0rVqxAZmamUctFRFRdMLAgIrJi8+fPx9SpUwEADRs2hCRJkCQJiYmJALS7Ql27dg2SJGH58uVYunQpAgIC4OTkhJCQEFy8eBF5eXmYPn06/P394e7ujvDwcNy9e1frfbdt24YnnngCLi4ucHV1RZ8+ffDbb7/pVeZPPvkE/fv3h4eHh8Z+SZIwYcIEfPHFF2jevDmcnJzQqVMnHDt2DEIILF++HA0bNoSrqyt69uyJS5cuabxeV1coVZ6xsbFo0aIFnJ2d0bZtW+zbt0+rXM8//zwePXqErVu36vU5iIhIEwMLIiIrNnr0aLzxxhsAgLi4OBw9ehRHjx5Fhw4dSn3dxx9/jCNHjuDjjz/GZ599hgsXLqB///545ZVXcO/ePaxfvx7Lli3DoUOHtFoW3n77bQwbNgwtW7bE9u3bERsbi7S0NAQHB6u7XpXkr7/+wtmzZxEaGqrz+L59+/DZZ5/hnXfewZYtW5CWloann34ab731Fo4cOYKPPvoI69atw7lz5zBo0CDo05t3//79+Oijj7BgwQLs2rULnp6eCA8Px5UrVzTS1a5dG4GBgdi/f3+ZeRIRkTY7cxeAiIjKr27duqhfvz4AoH379noPXvbw8MCePXtgYyPfX7p//z4mTZqEwMBAfPXVV+p0Fy5cwIcffohHjx7Bzc0NN2/exLx58zBhwgSsXLlSnS4sLAxNmzZFdHQ0tm3bVuL7/vLLLwBQYuCTk5OD7777Di4uLgDkFoeBAwciISEBp0+fhiRJAIB79+5h0qRJ+PPPPxEUFFTqZ83KysKhQ4dQo0YN9Xv7+/tj+/btmD59ukbaDh064NChQ6XmR0REurHFgoioGurXr586qACAFi1aAACefvppjXSq/Tdu3AAAHDx4EEqlEiNGjIBSqVRvCoUCPXr0UHfBKklSUhIAwNfXV+fx0NBQdVBR9P2feuopdVBRdP/169fL/KyhoaHqoAIAatWqBV9fX52v9fX1xd27d6FUKsvMl4iINLHFgoioGvL09NR47uDgUOr+7OxsAMCdO3cAAJ07d9aZb9FgRZesrCwAgEKhMGq5SuPl5aW1z9HRUV2WohQKBYQQyM7Ohqura5l5ExFRIQYWRESkN29vbwDAzp070aBBg3K/PiUlBX5+fkYtmzGkpKTA0dGRQQURUTkwsCAisnKOjo4AoPMOvLH16dMHdnZ2uHz5MgYNGmTw6wMDAwHIa060atXK2MWrsCtXrljc9LxERNaCgQURkZVTDV5esWIFRo4cCXt7ezRv3lxjXIGxBAQEYMGCBZg1axauXLmCvn37ombNmrhz5w5+/fVXuLi4lLrAXJcuXeDk5IRjx47h2WefNXr5KqKgoAC//vorXnnlFXMXhYjIKnHwNhGRlQsJCcGMGTOwd+9edOvWDZ07d8apU6dM9n4zZszAzp07cfHiRYwcORJ9+vRBVFQUrl+/ju7du5f6WgcHBwwePFhj5ilLkZiYiNTUVDz//PPmLgoRkVWShD6TgBMRERnJyZMn0blzZxw7dgxdunQxd3HUXnzxRVy5cgVHjhwxd1GIiKwSAwsiIqp0Q4cORUZGhs4VsM3h8uXLaNGiBeLj49GtWzdzF4eIyCqxKxQREVW69957D507d0ZaWpq5iwJAXqfjo48+YlBBRFQBbLEgIiIiIqIKY4sFERERERFVGAMLIiIiIiKqMAYWRERERERUYQwsiIiIiIiowhhYEBERERFRhTGwICIiIiKiCmNgQUREREREFcbAgoiIiIiIKuz/V5FC+5k/uUoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x550 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Ca = np.loadtxt(f'{model_ws}/PHT3D002.ACN') * 1000\n",
    "Na = np.loadtxt(f'{model_ws}/PHT3D004.ACN') * 1000\n",
    "T  = np.loadtxt(f'{model_ws}/PHT3D005.ACN') * 1000\n",
    "\n",
    "T_dist = np.array([0, 0, 0, 0.3, 0.43, 1.4, 1.92, 2.39, 2.93, 3.6, 4.10, 4.58, 4.76, 3.3, 2.97, 3.08, 3.0, 2.72, 1.90, 1.20, 1.06, 1.06, 1.06, 1.06, 1.06, 1.0, 0.86, 0.41, 0.16, 0, 0, 0, ])\n",
    "T_meas = T_dist * 0.007 / 8.38 # cm * mol/l / cm\n",
    "T_meas = T_meas * 1000 # mmol/l\n",
    "\n",
    "Na_dist = np.array([0.42, 0.1, 0.08, 0.1, 0.15, 0.63, 1.61, 2.31, 3.0, 3.7, 4.12, 4.72, 5.12, 8.58, 8.95, 9.8, 9.3, 9.12, 7.58, 5.95, 5.56, 5.6, 5.7, 5.55, 5.48, 5.9, 5.98, 5.98, 5.9, 5.92, 5.91, 5.8])\n",
    "Na_meas = Na_dist * 0.007 / 8.38\n",
    "Na_meas = Na_meas * 1000 # mmol/l\n",
    "\n",
    "Ca_dist = np.array([0.3, 0.33, 0.4, 0.48, 0.43, 0.78, 0.72, 0.4, 0.3, 0.28, 0.32, 0.25, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.87, 1.01, 1.03, 1.05, 1.06, 1.06, 1.05, 0.97, 0.7, 0.49, 0.28, 0.31, 0.3])\n",
    "Ca_meas = Ca_dist * 0.007 / 8.38\n",
    "Ca_meas = Ca_meas * 1000 # mmol/l\n",
    "\n",
    "times_meas = np.arange(4.7498, 4.7498 + 9.4995 * (31 + 1), 9.4995)\n",
    "\n",
    "times = np.array([0, 2.083E-03, 4.167E-03, 6.250E-03, 8.333E-03, 1.042E-02, .0125, 1.458E-02, 1.667E-02, .01875, 2.083E-02, 2.292E-02, .025, 2.708E-02, 2.917E-02, .03125, 3.333E-02, 3.542E-02, 3.750E-02, 3.958E-02, 4.167E-02, 4.375E-02, 4.583E-02, 4.792E-02, 5.000E-02, 5.208E-02, 5.417E-02, 5.625E-02, 5.833E-02, 6.042E-02, 6.250E-02, 6.458E-02, 6.667E-02, 6.875E-02, 7.083E-02, 7.292E-02, 7.500E-02, 7.708E-02, 7.917E-02, .08125, 8.333E-02, 8.542E-02, .0875, 8.958E-02, 9.167E-02, 9.375E-02, 9.583E-02, 9.792E-02, 1.000E-01, .1020833, .1041666, .10625, .1083333, .1104166, .1125, .1145833, .1166666, .11875, .1208333, .1229166, .125, .1270833, .1291666, .13125, .1333333, .1354166, .1375, .1395833, .1416666, .14375, .1458333, .1479166, .1499999, .1520833, .1541666, .1562499, .1583333, .1604166, .1624999, .1645833, .1666666, .1687499, .1708333, .1729166, .1749999, .1770833, .1791666, .1812499, .1833332, .1854166, .1874999, .1895832, .1916666, .1937499, .1958332, .1979166, .1999999, .2020832, .2041666, .2062499, .2083332])\n",
    "times = times * 24. * 60.\n",
    "\n",
    "ncol = 180\n",
    "index = np.arange(ncol - 1, len(times) * ncol, ncol)\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(8, 5.5))\n",
    "\n",
    "# 第一个子图\n",
    "axes[0].plot(times, Ca[index], 'm', label='Ca (PHT3D)')\n",
    "axes[0].plot(times_meas, Ca_meas, 'mo', label='Ca (observed)')\n",
    "axes[0].grid(True)\n",
    "axes[0].set_xlim([0, 300])\n",
    "axes[0].set_ylim([0, 2])\n",
    "axes[0].tick_params(axis='x', labelsize=12)\n",
    "axes[0].tick_params(axis='y', labelsize=12)\n",
    "axes[0].set_ylabel('Ca (mol/l)', fontsize=12)\n",
    "axes[0].legend(fontsize=12)\n",
    "\n",
    "# 第二个子图\n",
    "axes[1].plot(times, T[index], 'r', label='Tenside (PHT3D)')\n",
    "axes[1].plot(times_meas, T_meas, 'ro', label='Tenside (observed)')\n",
    "axes[1].grid(True)\n",
    "axes[1].set_xlim([0, 300])\n",
    "axes[1].set_ylim([0, 5])\n",
    "axes[1].tick_params(axis='x', labelsize=12)\n",
    "axes[1].tick_params(axis='y', labelsize=12)\n",
    "axes[1].set_ylabel('T (mol/l)', fontsize=12)\n",
    "axes[1].legend(fontsize=12)\n",
    "\n",
    "# 第三个子图\n",
    "axes[2].plot(times, Na[index], 'b', label='Na (PHT3D)')\n",
    "axes[2].plot(times_meas, Na_meas, 'bo', label='Na (observed)')\n",
    "axes[2].grid(True)\n",
    "axes[2].set_xlim([0, 300])\n",
    "axes[2].set_ylim([0, 10])\n",
    "axes[2].set_xlabel('time (min)', fontsize=12)\n",
    "axes[2].tick_params(axis='x', labelsize=12)\n",
    "axes[2].tick_params(axis='y', labelsize=12)\n",
    "axes[2].set_ylabel('Na (mol/l)', fontsize=12)\n",
    "axes[2].legend(fontsize=12)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}