MF6 with MST3D-USGS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flopy is installed in /Users/mattijnvanhoek/miniconda3/lib/python3.7/site-packages/flopy\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace = r'/Users/mattijnvanhoek/modflow6/output/mf-mt-hatari'\n",
    "if not os.path.exists(workspace):\n",
    "    os.makedirs(workspace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# name\n",
    "name = 'mf-hatari'\n",
    "\n",
    "# exe name\n",
    "exe_name = r'/Users/mattijnvanhoek/modflow6/bin/mf6'\n",
    "\n",
    "sim = flopy.mf6.MFSimulation(\n",
    "    sim_name=name,\n",
    "    exe_name=exe_name, \n",
    "    version='mf6',\n",
    "    sim_ws=workspace\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model domain and grid definition\n",
    "Lx = 310\n",
    "Ly = 310\n",
    "nlay = 1\n",
    "nrow = 31\n",
    "ncol = 31\n",
    "delr = Lx / ncol\n",
    "delc = Ly / nrow\n",
    "top = np.ones((nrow, ncol))\n",
    "botm = np.zeros((nrow, ncol))\n",
    "perlen = 27."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis = flopy.mf6.ModflowTdis(\n",
    "    simulation=sim, \n",
    "    time_units='DAYS', \n",
    "    pname='tdis',\n",
    "    perioddata=((perlen, 1, 1.),)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_nam_file = '{}.nam'.format(name)\n",
    "gwf = flopy.mf6.ModflowGwf(\n",
    "    simulation=sim, \n",
    "    modelname=name, \n",
    "    model_nam_file=model_nam_file\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ims = flopy.mf6.ModflowIms(\n",
    "    simulation=sim, \n",
    "    pname='ims', \n",
    "    complexity='SIMPLE'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "idomain = np.ones((nlay, nrow, ncol), dtype=np.int32)\n",
    "\n",
    "dis = flopy.mf6.ModflowGwfdis(\n",
    "    model=gwf, \n",
    "    pname='dis', \n",
    "    nlay=nlay, \n",
    "    nrow=nrow, \n",
    "    ncol=ncol,\n",
    "    delr=delr,\n",
    "    delc=delc,\n",
    "    top=top,\n",
    "    botm=botm, \n",
    "    idomain=idomain\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "strt = 15 #starting head\n",
    "\n",
    "ic = flopy.mf6.ModflowGwfic(\n",
    "    model=gwf, \n",
    "    pname='ic', \n",
    "    strt=strt\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "npf = flopy.mf6.ModflowGwfnpf(\n",
    "    model=gwf, \n",
    "    pname='npf', \n",
    "    icelltype=0, \n",
    "    save_flows=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAEICAYAAAB/KknhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAASdklEQVR4nO3df7BcZX3H8fcHcnODECAxQJPwS1NE0MHAXCIUaqkgYKYzgA5WUIodnKCVUSp2zKAjsUMVHQE7U6SEhpIiglRBGYs/YgQpU6RcMEJosEEaIMlNLghIohgS+PaP57l4vO69u3d/Jnk+r5mdPXt+fvfsfvY559k99yoiMLOy7NLrAsys+xx8swI5+GYFcvDNCuTgmxXIwTcrkINvRZF0vaRLO71uSSdIWtuJ7bTDThF8SWskndTrOqokzZR0u6T1kkLSwU2s42xJg5I2SxqS9F1Jx7e/2le319Y3q6QPSLqnzjx3SfpgJ+uYiPy6Lcn7e5OkRyV9VtLuvainU3aK4PeapEk1Rr8CfA94d5Pr/DjwZeBzwH7AgcBXgNOaLNPqkDQduBfYDTg2IqYC7wD2Bub0sra2i4gd/gasAU6qMX4a8B3gaeC5PLx/nnYm8MCo+S8CvpWH+4EvAU8CG4F/BnbL004A1gKfBDYAN4xT2yQggIMn8Hz2AjYDZ44zTz/pg2F9vn0Z6B9V30XAMDAE/HVl2fnA/wCbgHXAJ4DdgRdJH1ib820WMI8Uhufzev4JmFxZVwAfAlbnfXwVIOAw4LfAy3ldz4/xPO4CPjhq3AnA2lH7Y0ne/jrgUmDXPG0O8CPgl8AzwI3A3pVljwQezM/168DNwKVj1HIp8DCwyzj7/Y3AMuBZ4OfAeyrTrh9Zd43n8Mlc+6a83Im9zMzO3uLvAvwrcBCpxXyR9MYFuB14naTDKvO/H7ghD38BeAMwF/hjYDbwmcq8fwRMz+teMNHCJC2U9J0xJh8LTAFuG2cVnwKOyfW9hRTQT4+qb69c93nAVZKm5WlLgPMjtWhvBn4UEb8G3gmsj4g98m09Kbh/C8zIdZ0I/M2oWv4CODrX8R7glIhYRfpAuDeva+9xd8j4lgLbSK/DkcDJwMjpgYDPkz6kDgMOABYBSJoMfIv0mk4H/p3xj8BOAm6NiFdqTcyH+8uArwH7AmcBX5H0pvGKl3QocAFwdN7np5Aaq97p5adOu26M0eLXmG8u8Fzl8dXAP+ThN5FarH7Sm+nXwJzKvMcC/1f5NH8JmNLANptp8d8HbKgzzy+A+ZXHpwBrKvW9CEyqTB8GjsnDTwLnA3uOWucJVFqpMbZ7IXBb5XEAx1ce3wIszMMfAO6ps767gN+QjihGbptH6iCd5mwhH23lcWcBd46xvtOBn+bht5GOhlSZ/l+M3eKvBj40Tq1/CfznqHHXAJfk4eup0eKTPrCGSR8sfb3KSfW2U7f4kl4j6RpJT0h6Abgb2FvSrnmWpcDZkgScA9wSEVuAfYDXAA9Iel7S86Tz9X0qq386In7bodJ/CcwYo+9gxCzgicrjJ/K4V9cREdsqj38D7JGH30063H9C0o8lHTvWRiS9QdJ3JG3I+/BzpNa/asMY22nURyNi75Eb6QhixEFAHzBUeS2uIbW4SNpX0s2S1uX6vlqpbxawLnL6suo+G+2XwMxxph8EvHWkjlzL+0hHV2OKiMdIH5iLgOFc76zxlum0nTr4pHPcQ4G3RsSepBYAUotORPyE1HL/KXA2vzvMf4bUYr6p8obcKyKqb+hOXtZ4L+n8+PRx5llPeiOOODCPqysi7o+I00jh+RaplYbaz+lq4FHgkLwPLybvv0Y21eB843mK1OLPqLwWe0bEyOH15/N2jsj1vb9S3xAwO3+wjzhwnG39EDhD0li5eAr4cfVDKtJpzIfrPYmI+FpEHE96zYJ0KtkzO1Pw+yRNqdwmAVNJAX4+99heUmO5fyOd92+LiHsAIp3jXQtcKWmkZZkt6ZSJFCRpCunUAaA/P64rIn5F6k+4StLp+cilT9I7JX0xz3YT8GlJ+0iakef/agM1TZb0Pkl7RcRW4AXSeTykTszXStqrssjUPM9mSW8E6r7JKzYC++dz7aZExBDwA+BySXtK2kXSHEl/VqlvM+k1ng38XWXxe0l9Ax+VNEnSu0h9IWO5AtgTWCrpIHj1db9C0hGkzuE3SDonvx59ko4e1U/0ByQdKuntkvpJH+gv8rt93hM7U/DvIO3QkdsiUk/3bqQW/Cekw/XRbiB1cN0wavwngceAn+RDyB+Sjh4m4kXSmxJSq/niyARJF0v67lgLRsQVwMdJHXZPk1qbC0gtNKQe6EHgIVJP9IN5XCPOAdbk5/UhUitJRDxK+kB5PB/KziL1+J9N6o2+ltQz3qgfAY8AGyQ9M4HlRvsrYDLpm4jngG/wu0PyzwJHAb8C/gO4dWShiHgJeBepr+E50jn6q9NHi4hngT8BtgL3SdoELM/rfiwiNpE6Ft9LOrraQGq5+2uv8VX9wGWk9+EG0pHWxY088U7R75/+lEfSbqSOl6MiYnWv6zHrhp2pxW/Wh4H7HXoryXi9xjs9SWtIHUHjdaKZ7XSKP9Q3K5EP9c0K1NVD/cnqjynsVBc5mW1XfsuveSm21P2dRUvBl3Qq8I/ArsC/RMRl480/hd15q05sZZNmNo77YnlD8zV9qJ9/9noV6cKOw4GzJB3e7PrMrHtaOcefR/pRw+P5hxI342vFzXYIrQR/NunXZCPW5nG/R9KC/FdkBreypYXNmVm7tBL8Wh0If/DdYEQsjoiBiBjoq/vLRjPrhlaCv5b0Rw9G7E+DV4eZWW+10qt/P3CIpNeR/qTQe0kXczTl++tXtFCKWVlOmTW3peWbDn5EbJN0AfB90td510XEIy1VY2Zd0dL3+BFxB+lyWDPbgfgnu2YFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCtTSv9CStAbYBLwMbIuIgXYUZWad1VLwsz+PiGfasB4z6xIf6psVqNXgB/ADSQ9IWlBrBkkLJA1KGtzKlhY3Z2bt0Oqh/nERsV7SvsAySY9GxN3VGSJiMbAYYE9Njxa3Z2Zt0FKLHxHr8/0wcBswrx1FmVlnNR18SbtLmjoyDJwMrGxXYWbWOa0c6u8H3CZpZD1fi4jvtaUqM+uopoMfEY8Db2ljLWbWJf46z6xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgeoGX9J1koYlrayMmy5pmaTV+X5aZ8s0s3ZqpMW/Hjh11LiFwPKIOARYnh+b2Q6ibvAj4m7g2VGjTwOW5uGlwOltrsvMOqjZc/z9ImIIIN/vO9aMkhZIGpQ0uJUtTW7OzNqp4517EbE4IgYiYqCP/k5vzswa0GzwN0qaCZDvh9tXkpl1WrPBvx04Nw+fC3y7PeWYWTc08nXeTcC9wKGS1ko6D7gMeIek1cA78mMz20FMqjdDRJw1xqQT21yLmXWJf7lnViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xAjfzvvOskDUtaWRm3SNI6SSvybX5nyzSzdmqkxb8eOLXG+CsjYm6+3dHessysk+oGPyLuBp7tQi1m1iWtnONfIOmhfCowbayZJC2QNChpcCtbWticmbVLs8G/GpgDzAWGgMvHmjEiFkfEQEQM9NHf5ObMrJ2aCn5EbIyIlyPiFeBaYF57yzKzTmoq+JJmVh6eAawca14z2/5MqjeDpJuAE4AZktYClwAnSJoLBLAGOL+DNZpZm9UNfkScVWP0kg7UYmZd4l/umRXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrkINvViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXIwTcrUN3gSzpA0p2SVkl6RNLH8vjpkpZJWp3vp3W+XDNrh0Za/G3ARRFxGHAM8BFJhwMLgeURcQiwPD82sx1A3eBHxFBEPJiHNwGrgNnAacDSPNtS4PROFWlm7TWhc3xJBwNHAvcB+0XEEKQPB2DfMZZZIGlQ0uBWtrRWrZm1RcPBl7QH8E3gwoh4odHlImJxRAxExEAf/c3UaGZt1lDwJfWRQn9jRNyaR2+UNDNPnwkMd6ZEM2u3Rnr1BSwBVkXEFZVJtwPn5uFzgW+3vzwz64RJDcxzHHAO8LCkFXncxcBlwC2SzgOeBM7sTIlm1m51gx8R9wAaY/KJ7S3HzLrBv9wzK5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVyME3K5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVyME3K5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVyME3K5CDb1agRv5p5gGS7pS0StIjkj6Wxy+StE7Sinyb3/lyzawdGvmnmduAiyLiQUlTgQckLcvTroyIL3WuPDPrhEb+aeYQMJSHN0laBczudGFm1jkTOseXdDBwJHBfHnWBpIckXSdp2hjLLJA0KGlwK1taKtbM2qPh4EvaA/gmcGFEvABcDcwB5pKOCC6vtVxELI6IgYgY6KO/DSWbWasaCr6kPlLob4yIWwEiYmNEvBwRrwDXAvM6V6aZtVMjvfoClgCrIuKKyviZldnOAFa2vzwz64RGevWPA84BHpa0Io+7GDhL0lwggDXA+R2p0MzarpFe/XsA1Zh0R/vLMbNu8C/3zArk4JsVyME3K5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVyME3K5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVyME3K5CDb1YgB9+sQA6+WYEcfLMCOfhmBXLwzQrk4JsVqO6/0JI0Bbgb6M/zfyMiLpE0Hfg6cDDpf+e9JyKea7aQU2bNbXZRM5ugRlr8LcDbI+ItwFzgVEnHAAuB5RFxCLA8PzazHUDd4EeyOT/sy7cATgOW5vFLgdM7UqGZtV1D5/iSds3/InsYWBYR9wH7RcQQQL7fd4xlF0galDS4lS3tqtvMWtBQ8CPi5YiYC+wPzJP05kY3EBGLI2IgIgb66G+2TjNrown16kfE88BdwKnARkkzAfL9cNurM7OOqBt8SftI2jsP7wacBDwK3A6cm2c7F/h2p4o0s/ZSRIw/g3QEqfNuV9IHxS0R8feSXgvcAhwIPAmcGRHP1lnX08AT+eEM4JnWym+r7ake11Kba6mtWstBEbFPvQXqBr9TJA1GxEBPNl7D9lSPa6nNtdTWTC3+5Z5ZgRx8swL1MviLe7jtWranelxLba6ltgnX0rNzfDPrHR/qmxXIwTcrUE+CL+lUST+X9Jiknl7VJ2mNpIclrZA02OVtXydpWNLKyrjpkpZJWp3vp/WwlkWS1uV9s0LS/C7VcoCkOyWtkvSIpI/l8V3fN+PU0vV9I2mKpP+W9LNcy2fz+Invl4jo6o30Q6BfAK8HJgM/Aw7vdh2VetYAM3q07bcBRwErK+O+CCzMwwuBL/SwlkXAJ3qwX2YCR+XhqcD/Aof3Yt+MU0vX9w0gYI883AfcBxzTzH7pRYs/D3gsIh6PiJeAm0mX+BYnIu4GRv/asSeXO49RS09ExFBEPJiHNwGrgNn0YN+MU0vXRdKWS+R7EfzZwFOVx2vp0Y7MAviBpAckLehhHSMauty5iy6Q9FA+FejKaUeVpIOBI0mtW0/3zahaoAf7ppVL5Kt6EXzVGNfL7xSPi4ijgHcCH5H0th7Wsr25GphD+stLQ8Dl3dy4pD2AbwIXRsQL3dx2A7X0ZN9EC5fIV/Ui+GuBAyqP9wfW96AOACJifb4fBm4jnYr00nZzuXNEbMxvtFeAa+nivpHURwrajRFxax7dk31Tq5Ze7pu8/ZYuke9F8O8HDpH0OkmTgfeSLvHtOkm7S5o6MgycDKwcf6mO224udx55M2Vn0KV9I0nAEmBVRFxRmdT1fTNWLb3YN229RL6bvZKV3sn5pN7RXwCf6kUNuY7Xk75V+BnwSLdrAW4iHSZuJR0JnQe8lvTHS1fn++k9rOUG4GHgofzmmtmlWo4nnf49BKzIt/m92Dfj1NL1fQMcAfw0b3Ml8Jk8fsL7xT/ZNSuQf7lnViAH36xADr5ZgRx8swI5+GYFcvDNCuTgmxXo/wGBTsUEFXbpNAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "chd_data = []\n",
    "chd = 15    \n",
    "    \n",
    "for c in range(1,31):\n",
    "    dd = [(0, 0, c), chd]\n",
    "    chd_data.append(dd)\n",
    "for c in range(1,31):\n",
    "    dd = [(0, c, 30), chd]\n",
    "    chd_data.append(dd)\n",
    "for c in range(0,30):\n",
    "    dd = [(0, 30, c), chd]\n",
    "    chd_data.append(dd)\n",
    "for c in range(0,30):\n",
    "    dd = [(0, c, 0), chd]\n",
    "    chd_data.append(dd)    \n",
    "\n",
    "chd = flopy.mf6.ModflowGwfchd(\n",
    "    model=gwf, \n",
    "    pname='chd', \n",
    "    maxbound=len(chd_data), \n",
    "    stress_period_data=chd_data, \n",
    "    save_flows=True\n",
    ")\n",
    "ra = chd.stress_period_data.get_data(key=0)    \n",
    "\n",
    "ibd = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
    "for k, i, j in ra['cellid']:\n",
    "    ibd[k, i, j] = -1\n",
    "\n",
    "ilay = 0\n",
    "plt.imshow(ibd[ilay, :, :], interpolation='none')\n",
    "plt.title('Layer {}: Constant Head Cells'.format(ilay + 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAEICAYAAAB/KknhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAARCElEQVR4nO3df5BV5X3H8ffHFX9CIyAo4iYoEiumujo7aGpVjBoJTaNOa6qZGto6wfywoy2ZhrHTBNtJapyo7UwSE6xUqkbrVK3U2ihhVOIkEleDCoMWtSjICog6YJtGXb/94zxrbta7u5e9596z+HxeM3fu+XnP9z67nz3POfecvYoIzCwve1RdgJm1n4NvliEH3yxDDr5Zhhx8sww5+GYZcvCt5SSdLOmZqusYjKQNks5Iw4sk3Vx1Ta3m4A9Q+0swWkiaImmZpM2SQtK0XVy/re8p1XhE/3hE/DgijmzX9m14Dv4oI2nPOpPfAX4I/H6by7H3KQe/QZLGS7pH0jZJr6XhQ9O88yQ9NmD5BZL+LQ3vLelbkl6UtEXS9yTtm+bNlrRJ0lckvQz808BtR8SWiPgu8GgL3tfnJD0r6dXUqzikZt7RkpaneVskXZ6mz5L0U0mvS+qV9G1Je6V5K9PqT0h6Q9If9r/Hmtc9StKDaf21kj5VM+9GSd+R9B+SdkpaJWn6ILUvlbQgDU9NPY0vpvEjUt1K45+UtDpt8yeSjmmgbfaRdLOk7Wm9RyUdtMuNPAo5+I3bgyKUHwI+CPwC+Haatww4TNJRNcv/EXBTGv4m8GGgCzgCmAp8tWbZg4EJ6bXn72phkhZKumcE630M+Dvg08AU4AXgtjRvHPAjip7GIanuFWnVPuDPgQOBjwKnA18EiIhT0jLHRsTYiPiXAdscA/w7cD8wGfgz4BZJtYcCFwBXAOOBZ4GvD/IWHgJmp+FTgefTM8ApwI8jIiQdDywBLgYmAt8Hlknae5gmmgd8AOhM632e4ue++4sIP2oewAbgjAaW6wJeqxm/Dvh6Gj4aeA3YGxDwP8D0mmU/Cvx3Gp4NvAns08A29wQCmFbGewJuAK6qGR8LvAVMowjfzxt8/cuAu2rGAziiZnw2sCkNnwy8DOxRM/9WYFEavhH4x5p5c4GnB9nudOB1ij/K36MIdv92lgJ/UfOz+dsB6z4DnDqwfYBFwM1p+E+BnwDHVP17WfbDe/wGSdpP0vclvSBpB7ASOEBSR1pkKfCZ1LW8ELg9In4JTAL2Ax5L3cXXKfaik2pefltE/F/73s27DqHYywMQEW8A2yl6JJ3Ac/VWkvThdKjzcmqLb1Ds/Rvd5saIeKdm2gtpm/1erhn+X4o/SO8REc8Bb1D8ET4ZuAfYnHoPp1L0CKDoSS3ob//0M+hMtQzlJuA+4LZ0YvWq1GPZ7Tn4jVsAHAmcEBG/QdGVhGKPTkQ8QrHnPhn4DL/q5r9C0T08OiIOSI8PRETtL3NVt0hupggFAJL2p+jSvgRspNij1nMd8DQwI7XF5aR2aHCbnZJqf/c+mLY5Eg8BfwDsFREvpfHPUhwmrE7LbKTojR1Q89gvIm4d6oUj4q2IuCIiZgK/DXwyvfZuz8Gvb0w6sdP/2BMYRxHg1yVNAL5WZ71/pjjufzsiHgZIe7brgWslTYZ3T0SdtSsFSdqH4tABYO803ux7+gHwJ5K60vHuN4BVEbGBYu95sKTL0snJcZJOSK81DtgBvCHpN4EvDNjWFuDwQepYRXHo85eSxkiaDfwe6dzCCDwEXELRAwN4kOK8wcMR0ZemXQ98XtIJKuwv6XfTeYxBSTpN0m+lXt0OisOgvqHW2V04+PXdSxHy/sci4O+BfSn24I9QdNcHugn4CL/a2/f7CsVJqkdS1/hHFL2HXfELim4tFHvbd08ySbpc0n8Os/573lNErAD+GrgD6KXYw58PEBE7gTMpQvkysB44Lb3Wlyl6NTspQvVrJ/Ao2mtp6lZ/unZGRLwJfAr4BEVbfhf4bEQ8PWwL1PcQxR+i/uA/THFo1T9ORPQAn6P4o/waxc/ijxt47YOBf6UI/bq0rffFxT1KJzGsBCo+otsKHB8R66uux2ww3uOX6wvAow69jXb1rhKzEZC0geIE1zkVl2I2LHf1zTLkrr5Zhtra1T9wQkdM63xfXP9gNipt2PgWr7zaN+w1FU0FX9Ic4B+ADorLLK8cavlpnWP42X2dzWzSzIYw66yNDS034q5+uqjhOxSfx84ELpA0c6SvZ2bt08wx/izg2Yh4Pl2UcRtwdjllmVkrNRP8qRTXQPfbxK/faAGApPmSeiT1bNv+vrja0Wy310zw651AeM9ngxGxOCK6I6J70sSOOquYWbs1E/xNFLc29juU4s4rMxvlmgn+o8AMSYelf7t0PsV/ojGzUW7EH+dFxNuSLqH4RwUdwJKIWFtaZWbWMk19jh8R91Lc7mlmuxFfsmuWIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLENNfYWWpA3ATqAPeDsiussoysxaq6ngJ6dFxCslvI6ZtYm7+mYZajb4Adwv6TFJ8+stIGm+pB5JPdu29zW5OTMrQ7Nd/ZMiYrOkycBySU9HxMraBSJiMbAYoPvYfaLJ7ZlZCZra40fE5vS8FbgLmFVGUWbWWiMOvqT9JY3rHwY+DqwpqzAza51muvoHAXdJ6n+dH0TED0upysxaasTBj4jngWNLrMXM2sQf55llyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlyME3y9CwwZe0RNJWSWtqpk2QtFzS+vQ8vrVlmlmZGtnj3wjMGTBtIbAiImYAK9K4me0mhg1+RKwEXh0w+WxgaRpeCpxTcl1m1kIjPcY/KCJ6AdLz5MEWlDRfUo+knm3b+0a4OTMrU8tP7kXE4ojojojuSRM7Wr05M2vASIO/RdIUgPS8tbySzKzVRhr8ZcC8NDwPuLuccsysHRr5OO9W4KfAkZI2SboIuBI4U9J64Mw0bma7iT2HWyAiLhhk1ukl12JmbeIr98wy5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0yNOw36Vhezjqka8j5921e3aZKrJUa+e68JZK2SlpTM22RpJckrU6Pua0t08zK1EhX/0ZgTp3p10ZEV3rcW25ZZtZKwwY/IlYCr7ahFjNrk2ZO7l0i6cl0KDB+sIUkzZfUI6ln2/a+JjZnZmUZafCvA6YDXUAvcPVgC0bE4ojojojuSRM7Rrg5MyvTiIIfEVsioi8i3gGuB2aVW5aZtdKIgi9pSs3oucCawZY1s9Fn2M/xJd0KzAYOlLQJ+BowW1IXEMAG4OIW1mht5M/p8zBs8CPigjqTb2hBLWbWJr5k1yxDDr5Zhhx8sww5+GYZcvDNMuTgm2XIwTfLkINvliEH3yxDDr5Zhhx8sww5+GYZcvDNMuTgm2XIwTfLkINvliEH3yxDDr5Zhhx8sww5+GYZcvDNMuTgm2XIwTfLkINvliEH3yxDDr5ZhoYNvqROSQ9IWidpraRL0/QJkpZLWp+ex7e+XDMrQyN7/LeBBRFxFHAi8CVJM4GFwIqImAGsSONmthsYNvgR0RsRj6fhncA6YCpwNrA0LbYUOKdVRZpZuXbpGF/SNOA4YBVwUET0QvHHAZg8yDrzJfVI6tm2va+5as2sFA0HX9JY4A7gsojY0eh6EbE4IrojonvSxI6R1GhmJWso+JLGUIT+loi4M03eImlKmj8F2NqaEs2sbI2c1RdwA7AuIq6pmbUMmJeG5wF3l1+embXCng0scxJwIfCUpNVp2uXAlcDtki4CXgTOa02JZla2YYMfEQ8DGmT26eWWY2bt4Cv3zDLk4JtlyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlyME3y5CDb5YhB98sQw6+WYYcfLMMOfhmGXLwzTLk4JtlqJEvzeyU9ICkdZLWSro0TV8k6SVJq9NjbuvLNbMyNPKlmW8DCyLicUnjgMckLU/zro2Ib7WuPDNrhUa+NLMX6E3DOyWtA6a2ujAza51dOsaXNA04DliVJl0i6UlJSySNH2Sd+ZJ6JPVs297XVLFmVo6Ggy9pLHAHcFlE7ACuA6YDXRQ9gqvrrRcRiyOiOyK6J03sKKFkM2tWQ8GXNIYi9LdExJ0AEbElIvoi4h3gemBW68o0szI1clZfwA3Auoi4pmb6lJrFzgXWlF+embVCI2f1TwIuBJ6StDpNuxy4QFIXEMAG4OKWVGhmpWvkrP7DgOrMurf8csysHXzlnlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDn4Zhly8M0y5OCbZcjBN8uQg2+WIQffLEMOvlmGHHyzDDXypZn7SPqZpCckrZV0RZo+QdJySevT8/jWl2tmZWhkj/9L4GMRcSzQBcyRdCKwEFgRETOAFWnczHYDwwY/Cm+k0THpEcDZwNI0fSlwTksqNLPSNXSML6kjfUX2VmB5RKwCDoqIXoD0PHmQdedL6pHUs217X1l1m1kTGgp+RPRFRBdwKDBL0kca3UBELI6I7ojonjSxY6R1mlmJdumsfkS8DjwIzAG2SJoCkJ63ll6dmbVEI2f1J0k6IA3vC5wBPA0sA+alxeYBd7eqSDMrlyJi6AWkYyhO3nVQ/KG4PSL+RtJE4Hbgg8CLwHkR8eowr7UNeCGNHgi80lz5pRpN9biW+lxLfbW1fCgiJg23wrDBbxVJPRHRXcnG6xhN9biW+lxLfSOpxVfumWXIwTfLUJXBX1zhtusZTfW4lvpcS327XEtlx/hmVh139c0y5OCbZaiS4EuaI+kZSc9KqvSuPkkbJD0labWknjZve4mkrZLW1Eyr5HbnQWpZJOml1DarJc1tUy2dkh6QtC7dCn5pmt72thmilra3Tam3yEdEWx8UFwI9BxwO7AU8Acxsdx019WwADqxo26cAxwNraqZdBSxMwwuBb1ZYyyLgyxW0yxTg+DQ8DvgvYGYVbTNELW1vG0DA2DQ8BlgFnDiSdqlijz8LeDYino+IN4HbKG7xzU5ErAQGXu1Yye3Og9RSiYjojYjH0/BOYB0wlQraZoha2i4KpdwiX0XwpwIba8Y3UVFDJgHcL+kxSfMrrKNfQ7c7t9Elkp5MhwJt/y9LkqYBx1Hs3SptmwG1QAVt08wt8rWqCL7qTKvyM8WTIuJ44BPAlySdUmEto811wHSK/7zUC1zdzo1LGgvcAVwWETvaue0GaqmkbaKJW+RrVRH8TUBnzfihwOYK6gAgIjan563AXRSHIlUaNbc7R8SW9Iv2DnA9bWwbSWMognZLRNyZJlfSNvVqqbJt0vabukW+iuA/CsyQdJikvYDzKW7xbTtJ+0sa1z8MfBxYM/RaLTdqbnfu/2VKzqVNbSNJwA3Auoi4pmZW29tmsFqqaJtSb5Fv51nJmrOTcynOjj4H/FUVNaQ6Dqf4VOEJYG27awFupegmvkXRE7oImEjxz0vXp+cJFdZyE/AU8GT65ZrSplp+h+Lw70lgdXrMraJthqil7W0DHAP8PG1zDfDVNH2X28WX7JplyFfumWXIwTfLkINvliEH3yxDDr5Zhhx8sww5+GYZ+n9/dMpInk69VAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# build a stress_period_data template with 1 wells over 1 stress period\n",
    "wel_periodrec = flopy.mf6.ModflowGwfwel.stress_period_data.empty(\n",
    "    model=gwf, \n",
    "    maxbound=1, \n",
    "    boundnames=True, \n",
    "    aux_vars=None,\n",
    "    stress_periods=[0]\n",
    ")\n",
    "\n",
    "wcol = 15  # x index for the well\n",
    "wrow = 15  # y index for the well\n",
    "wlay = 0\n",
    "q = 100. # postive is injection, negative is extraction\n",
    "\n",
    "# [idx_stress_period][idx_well] = ((wlay, wrow, wcol), q, name_well)\n",
    "wel_periodrec[0][0] = ((wlay, wrow, wcol), q, 'well_1')\n",
    "\n",
    "# build the well package\n",
    "wel = flopy.mf6.ModflowGwfwel(\n",
    "    model=gwf,\n",
    "    pname='wel', \n",
    "    print_input=True, \n",
    "    print_flows=True,\n",
    "    maxbound=1,\n",
    "    stress_period_data=wel_periodrec, \n",
    "    boundnames=True, \n",
    "    save_flows=True\n",
    ")\n",
    "\n",
    "ra = wel.stress_period_data.get_data(key=0)\n",
    "\n",
    "ibd = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
    "for k, i, j in ra['cellid']:\n",
    "    ibd[k, i, j] = -1\n",
    "\n",
    "ilay = 0\n",
    "plt.imshow(ibd[ilay, :, :], interpolation='none')\n",
    "plt.title('Layer {}: Location wells'.format(ilay + 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "headfile = '{}.hds'.format(name)\n",
    "head_filerecord = [headfile]\n",
    "budgetfile = '{}.cbb'.format(name)\n",
    "budget_filerecord = [budgetfile]\n",
    "saverecord = [('HEAD', 'ALL'), \n",
    "              ('BUDGET', 'ALL')]\n",
    "printrecord = [('HEAD', 'LAST')]\n",
    "    \n",
    "oc = flopy.mf6.ModflowGwfoc(\n",
    "    model=gwf, \n",
    "    pname='oc', \n",
    "    saverecord=saverecord, \n",
    "    head_filerecord=head_filerecord,\n",
    "    budget_filerecord=budget_filerecord,\n",
    "    printrecord=printrecord\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing simulation...\n",
      "  writing simulation name file...\n",
      "  writing simulation tdis package...\n",
      "  writing ims package ims...\n",
      "  writing model mf-hatari...\n",
      "    writing model name file...\n",
      "    writing package dis...\n",
      "    writing package ic...\n",
      "    writing package npf...\n",
      "    writing package chd...\n",
      "    writing package wel...\n",
      "    writing package oc...\n"
     ]
    }
   ],
   "source": [
    "# Write the datasets\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FloPy is using the following  executable to run the model: /Users/mattijnvanhoek/modflow6/bin/mf6\n",
      "                                   MODFLOW 6\n",
      "                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL\n",
      "                            VERSION 6.1.0 12/12/2019\n",
      "\n",
      "  MODFLOW 6 compiled Jan  2 2020 21:28:58 with GFORTRAN compiler (ver. 8.3.0)\n",
      "\n",
      "This software has been approved for release by the U.S. Geological \n",
      "Survey (USGS). Although the software has been subjected to rigorous \n",
      "review, the USGS reserves the right to update the software as needed \n",
      "pursuant to further analysis and review. No warranty, expressed or \n",
      "implied, is made by the USGS or the U.S. Government as to the \n",
      "functionality of the software and related material nor shall the \n",
      "fact of release constitute any such warranty. Furthermore, the \n",
      "software is released on condition that neither the USGS nor the U.S. \n",
      "Government shall be held liable for any damages resulting from its \n",
      "authorized or unauthorized use. Also refer to the USGS Water \n",
      "Resources Software User Rights Notice for complete use, copyright, \n",
      "and distribution information.\n",
      "\n",
      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2020/01/10 23:23:44\n",
      "\n",
      " Writing simulation list file: mfsim.lst\n",
      " Using Simulation name file: mfsim.nam\n",
      " Solving:  Stress period:     1    Time step:     1\n",
      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2020/01/10 23:23:44\n",
      " Elapsed run time:  0.023 Seconds\n",
      "\n",
      " Normal termination of simulation.\n",
      "\n",
      "Success is:  True\n"
     ]
    }
   ],
   "source": [
    "# Run the simulation\n",
    "success, buff = sim.run_simulation()\n",
    "print('\\nSuccess is: ', success)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAEICAYAAAA3EMMNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9d5gc13Wn/Z7qHCfnGWAQBjkRAHMCI5jEIFEybUm0LMna1YqSLetb27L92Vp/pi3vJ0uWLZtaypKVkxVJSswkCEaAIAEig4iTc+yc6u4f1TMYTOhpzHRPAOp9nnpmurvq1u3uql+fc+6554pSChMTE5OFijbXHTAxMTGZCaaImZiYLGhMETMxMVnQmCJmYmKyoDFFzMTEZEFjipiJicmCxhSxBYSI7BCRj2e57zYRaRn1+JCIbMtRPz4oIs+MeqxEZHku2k63FxSRpblqb0zb/yAifzzNYz8jIl/MdZ9MZoYpYjNARM6IyM1jnvuIiLwyV32aDKXUWqXUjkz7iEh9WpCsU7T1A6XUrbno10TCrJTyKqVO5aL9MecqAx4E/s80m3gU+JCIlOeuVyYzxRQxk/NiKoGb53wE+K1SKjKdg5VSUeBJDCE0mSeYIpZnRKRaRH4uIt0iclpEPjPqtctE5HURGRCRdhH5mojYR71+i4gcFZFBEfkaIBnO4xKRb4tIv4gcBi4d8/qI1Zg+7x4RGRKRThH5cnq3nem/A2mX7sq0ZfmqiHxFRPqAL0xibd4hIqdEpEdE/n8R0dLn+oKIfH9UP0asPRF5GLgW+Fr6fF9L7zPinopIgYh8N/35NYrIX41q+yMi8oqIfCn9vk+LyO0Zvo7bgZdG9WWbiLSIyJ+KSFf6O7hXRO4QkXdFpE9E/mJMGzuAOzOcw2SWMUUsj6RvtseBd4Aa4Cbgj0Vke3qXFPBZoBS4Mv36/0gfWwr8HPir9OsngasznO5vgGXpbTvw+xn2/SrwVaWUP73/T9PPX5f+W5h26V5PP74cOAWUAw9P0uZ9wFZgM3AP8NEM5wdAKfWXwMvAQ+nzPTTBbv8KFABLgesxrKA/GPX65cAxjM/ofwPfFJHJxH59et/RVAJOjO/nr4FvAB8CtmAI7F+Pic8dATZO9d5MZg9TxGbOr9KW1ICIDAD/Puq1S4EypdTfKqXi6TjPN4AHAJRSbyml3lBKJZVSZzBiNdenj70DOKyU+plSKgH8M9CRoR8fAB5WSvUppZqBf8mwbwJYLiKlSqmgUuqNKd5jm1LqX9P9nMwV+8f0uZvSff3dKdqcEhGxAL8DfF4pFUh/Rv8EfHjUbo1KqW8opVLAd4AqoGKSJguBwJjnEhifWwL4MYYYfjV9vkPAIWDDqP0DGKJqMk8wRWzm3KuUKhzeSFtSaRYD1WNE7i9I32QiskJEnhCRDhEZAv4e4yYCqAaahxtSxkz9ZianeszrjRn2/RiwAjgqIm+KyF1TvMdM551on8Z0f2ZKKWDn3PfSiGE1DTMi7EqpcPpf7yTt9QO+Mc/1pgUQYFigO0e9HhnTng8YnLLnJrOGKWL5pRk4PVrklFI+pdQd6dcfAY4CDWnX7i84G/dqB+qGG0q7SHVMTvuY1xdNtqNS6rhS6ncx3MN/BH4mIh5gspIm2ZQ6GXvutvT/IcA96rXK82i7B8NSWjym7dYs+jMR+zHEeyasxggPmMwTTBHLL7uBIRH5s3Tg3SIi60RkOOjuA4aAoIisAj456tjfAGtF5L3pEcHPMF4ARvNT4PMiUiQitcCnJ9tRRD4kImVKKR0YSD+dAroBHSP+dL78z/S564A/An6Sfn4fcJ2ILBKRAuDzY47rnOx8aQvpp8DDIuITkcXAnwDfn2j/LPgtZ9316XI9xgilyTzBFLE8kr4J3wNsAk5jWBb/wdmYyv8D/B5GnOUbnL3xUUr1AO8Hvgj0Ag3AqxlO978wXK3TwDPA9zLsextwSESCGEH+B5RS0bQ79jDwatr9veI83u6vgbcwROs3wDfT7+PZ9Pvan379iTHHfRW4Pz26OFEc79MY1twp4BXgh8C3zqNfo/kuxiiqazoHi4gTI1b5nWme3yQPiFkU0eRiQkT+HuhSSv3zNI79NFCnlPrT3PfMZLqYImZiYrKgMd1JExOTOUFE/khEDooxr/eP088Vi8izInI8/bdoqnZMETMxMZl1RGQd8IfAZRjJw3eJSAPw58DzSqkG4Pn044yYImZiYjIXrAbeUEqFlVJJjOlg92HM9hgeOPkOcO9UDc2LybzFpX5Vu2hsYQDB57gal61hTvpkYrIQiSWbGIy+iJEpcy4H9p7sUUqVzaR9i6tKqVRsyv1Uov8QEB311KNKqUdHPT6IkTpTgpFQfAewB6hQSrUDKKXaJYuKIVOKWHpYeSfgSO//M6XU34hIMcbQeT1wBviAUqo/fcznMbLCU8BnlFJPZzpH7aJyHt/55ZHHekrx7a+28upz/dz9e+Vce1vxVN00Mbno2fvaID/8Otxy00N88i8WYbOd62jV++7JNIsjK1QqhqNy+5T7RZt/HFVKbZ20HaWOiMg/As8CQYwE4uR0+pSNJRYDblRKBUXEBrwiIk8C78XwXb8oIn+O4bv+mYiswZgbuBZj6slzIrJi1NSOjJwjYB8s554PTjYNzsTEZDQ33V2KZhG+/29tPPL3TRMK2YwRIV1EZMYopb5JOp8wnfrSAnSKSFXaCqsCuqZqZ8reKINg+qEtvSkm913vAX6slIoppU4DJzCCd1NiCpiJycy44c4SPvSpat7ZFeCRv28ikRjvVs4EQdDEOuWWVVtpV1FEFmEYRT8CHuNsBZbfx0iizkhWkpqeLrMPQxWfVUrtYozvijEPD4zJuaMnA7dw7oTd4TY/IUZNqz19PUOmgJmY5Ii8C5loU25Z8nMxat89DnwqHY76InCLiBwHbkk/zkhWkpl2BTeJSCHwy/Tw6GRMVMtpXEZtOsj3KMCGzQ3KFDATk9xxw50lAOe4lrli8nJt54dS6toJnuvFqKuXNec1OqmUGhCRHRhz7ybzXVs4t6JBLWcrGkyIz3G1KWAmJjlmrJDlBmG+ZWZlMzpZBiTSAuYCbsYo3zLsu36Rc33Xx4AfilHyuBpj4vLuTOdw2RomFTClcqP657SZ8xZhtmdvXeyTxXJ/VUxxvjycMB/vQeTcK2O0kOXuHPNLxLLpTRXwoojsB97EiIk9wSS+a7oa5k+Bw8BTGL7ulCOTN72nZHrvwMTEJCM33FmCliPdEclpTCwnTGmJKaX2A5dM8PykvqtS6mEmr8VuYmKyYJGsRx9ni/nVGxMTk3lO7vLEcoUpYiYmJueFKWImJiYLFsFIeJ1PmCJmYmJyHpjuZE6ZbppBpnSI6bap56HN6aZtXCjpF9P9vZ9uOkSmw6Z722bsSh76mXcENG1+ycb86o2Jick8ZwEmu5qYmJiMxnQnTUxMFixixsRMTEwWOmK6kyYmJgsZ0xIzMTFZuIigaZa57sU5XLAiNt00iummSkz3fNNtc7pl7vJRFWQmjK26kC2ZbIFM73C6r2XqZcbUjGk2mo+qGbnASHbNjSUmIp8FPo7xSRwA/gBwM8naHZMxv+xCExOTeY7kpIqFiNQAnwG2KqXWARaMtTnMdSdNTEzySw5L8VgBl4hYMSywNhbqupMm8w9dVyST6S0x+q9OMqlQCqxWwWoTrFYNq02wWWXkOS2jH2WycJGcuJNKqVYR+RLQhLHu5DNKqWdEJPfrTppcWOi6oq83QXdXnJ6uBEODSYKBFMFAkkAgNfJ/MJCa9rQnEfB4LHj9Frw+Kz6fBa/P+N9fYKGk1E5ZhZ2SUhtW8wpcWAhIdtOOSkVkz6jH5yyeKyJFGFbXEmAA+C8R+dB0umReQhcoibhOS1OM5qYoXZ1xujridHcl6O1OkEyeq05ujzYiNuUVdpY1uPD6LDgcGlabhtWKYW2NWF6CCCQSYyy1pCKZ0InHFYHAsDimaGuNEQykCIfOFUbNAqWlNsoqDFErr7BTW+egrt6Jw2FGOuYjRrJrVlZ2T6bFczHK3J9WSnUDiMgvgKuYxrqTpohdAOi6oqsjzplTUc6cinDmVJSW5ih6uii4zSaUVdiprLazYZOPsoq0cJTb8futWKyz4/qlUopgIEl3V4LuzjhdnXG609uJY2FiMUPhNA2qaxwsXuqkfqmL+qVOqmocpos6T8jR6GQTcIWIuDHcyZuAPUCIidfumJR5L2L5SGtIZTouD6kZmdIhMqU8TNamUoquzjiHD4Q4eijI6eMRIhHjLA6HxqIlTm7cXsLiJU7qFrsoLLaOCMDYNnVAT3cwX9UvRt6hCG6/jcV+G4uXu4GzKQhKKYYGkjQ3Rmk8FaHxdJS33wzw6kuDANgdwpJlblat9bB6vYfqGsekFkEmrcuU0pHp1pxuSkemz9SSh5yO2ZD5XCS7KqV2icjPgLeBJLAXYwlHL/BTEfkYhtC9f6q25r2ImRjE4zrHj4Y5fCDI4f1BeroTAJRX2Nl8uZ/FS1wsXuqiosq+IC0WEaGgyEZBkY11m3yAIWzdnXEaTxvCdvxomF//Vxe//i8oKrayZr2XNRu8rFztweE03c/ZQXKWxKaU+hvgb8Y8HSOf606azC7JpOLgOwF2vTrIsUMhEgmFzS6sWOXhhu3FrFnvpaTMPtfdzBsiQnmlg/JKB5deWQBAf1+CIweDHN4f4q1dQ7z60gBWq7B8pZvLry5g4xYfDrspaHlj/lXiMUVsPtLeGuO1lwfY/dogwUCKgkIrV11XyNqNXpavdGOzaRdM4cPzpajYxlXXFXH1dUUkk4pTxw3rdO+eAN95tA2XW+PSKwq48toCFtW75rq7Fya5Wv8tR5giNk+IRFK8vXuI13YOcOZUFM0C6zf5uPLaQlat9WDJGEC5OLFahRWrPaxY7eHu+8s5fizMGy8P8NrOAXa+0E/tIgdXXlvIpVcU4PHOr/l+C5r5pWGmiM01gaEkzz/Vx0sv9BGPKapqHLzvgXIuvbIAj8/8erJF04SVqz2sXO3hdz6UYs+uQV7bOch//aCTX/6ki6uuK2T7XSUUFZuf6YwQUPNsYqf5jc4RwUCS557sY8fz/STiii2X+9l2czH1S50jo26ZRjxNJsftsXDdjcVcd2MxzU1Rdj7fzysv9fPazgGu2VbIrXcWU1hkm+tuLlzml4YtbBGbdjWKaaZDZE7NyC5VIhRM8cLTvbz0XD/xuM4ll/nZ/p5SKqocACQVI29stt/fXJAxdSHDzZIxjWLU+6+sdfKB36/ipjtKePY3Pex8oZ9XXxrg6usLufmOEgoKz4pZpgIfeoZvI5Onr03zu8j4/jIcNyvMs9HvBS1iC4lEQmfn8/08/XgP0ajOJZf6ue09pZRXO+a6axcFJWV2HvhINdvvLOXpJwwxe23nADfdVsJNt5eYMwSyJncpFrliShETkTrgu0Alxo/Ho0qpr4rIF4A/BLrTu/6FUuq36WM+D3wMSAGfUUo9nYe+LxjaWqJ859E22lpirNng4Z77y6mudQKZrTuT3FNabueDH63m1jtLefwXXTz5WA9vvjHIgx+vZlmDe667N/8RpsjSnX2yscSSwOeUUm+LiA94S0SeTb/2FaXUl0bvLCJrMOoCrQWqgedEZIVSKpXLji8EdF3x/NN9PPGLblxujU98upb1l/jmulsmQFmFnY9+spbjN4T4/jfb+co/NLL9zhJuv6cM6yxNw1qwLDRLLF0WY7g0RkBEjgA1GQ65B/ixUioGnBaRE8BlwOs56O+Cobcnzne+0cbJdyNs2OzjgQcr8flN732+0bDKw5//7RJ+/sNOnnqil0P7g/z+J2qoqjHd/EmZXxp2fhkfIlIPXALsSj/1kIjsF5FvpUtrgCFwzaMOa2EC0RORT4jInjHlOhY8Silef3mAv/9/T9PaFOPDH6vi45+qMQVsHuNyWfjQx6r5w4dq6O9L8sUvnOaFZ/rQzeHh8QhGYH+qbRbJWsRExAv8HPhjpdQQ8AiwDNiEYan90/CuExw+7mpQSj2qlNo6RbmOBUVgKMmj/9rC97/VTl29k7/4/5ZyxTWF2ZYuMZljNm3x85d/t5TV6zz8/Eed/OuXmujrTcx1t+YfksU2i2RlHoiIDUPAfqCU+gWAUqpz1OvfAJ5IP2wB6kYdXotRdnZWme7iHJkC7XqGMfjmxiiP/kszgaEU9/5OOdtuKUbThJSaqs3p9SWlT96XTMHHjNU28mR4ZPphzvRaphx7S4bchUxx54xVJQCP38rHP13L6y8P8IsfdfLFL5zm4w/V0rAyU9B/8lYz/X4tSDtPBGWZXyO5U/ZGDDPim8ARpdSXRz1fNWq3+4CD6f8fAx4QEYeILAEagN256/L8490jIb7yD2dQCj77l/XcuL1kQVaSMDEQEa66rog//ZsluD0WvvalJvbtGZrrbs0f5pkllo2kXg18GLhRRPaltzuA/y0iB0RkP3AD8FkApdQh4KfAYeAp4FMX8sjku0dCPPLPzZSU2vncX9VTt9g5110yyRHllQ7+5C+N7/RbX281hWwYkam3WSSb0clXmFhbf5vhmIeBh2fQrwXBsICVltl56H8uMoP3FyAer4VP/kkdX/9yM9/6eisf/e+waat/rrs1t8wzL8O866bJaAH7zJ8uuiAma6dSiv7+BJ1dMbq744RCKWIxnVhUN/7GdaJRnXjMmITjcGjnbE6Hht2h4XFbKC2zU1Fup6TEjrbA865cLgv/43N1/Ps/mUKWK3dRRFZiLJI7zFLgrzES689r8dyFf+fNAWMFzOe3LrjJ2qFQilOnwnR0xujsitPVFaO7O0Fq1GiCzSbjBKqgwIrToYFALKaIxVIEg0l6e9NCFzMWChlG06C01BC08nIHFRV2li1141tgom8K2Shy4C4qpY5hZDYgIhagFfglZxfP/aKI/Hn68Z9lamthXUnzgIkEbCGg64qW1ihHjoY4dixEc0sUpYzrsbjYRnmZnRUrPJSXOygvs1NebsfpnF4NrlhMp7snTldnjK7uON1dcdo7Yhw6HByp519d7WDVSg+rV3pYtMi1IOqlmUKGccHk/ru6CTiplGoUkXuAbennvwPsYKGLWMaKDNOt5DDNahRNjVG+/tWzMTCP76wFljkdIsNrGc4Xy3BcMkOKRSJ9XDSa4t2jQY4fD3HieIhwOAUCNTVOrt1WyrLlXioqHdhs48d3UkBoVIrUVHM8z7muNY2icitF5W5Wju5zUqe7K8bJEyFOHg/y4kt9vPBiH06nxtJlHhpWeFi12ofbbYjnBN0awZqhP5nmclsyLBSSzTCXw2nhv/9JHY+kY2Sf+pyFVas9k+6fqfpFJinIVOEio4TMzkoh2eyVcd3JMTwA/Cj9v7l4br4IDCV59KvNuN3Gr/F8tsDaWqO8ubufA/uHSCQUHo+F5Q2GSNQv9eL2zE3frVaNqmoXVdUurrmulGgkRePpIMffNUT28KEATzzWydp1PrZeWsiyJa55mSjscln475+t4ysPn+Gb/9bKn/11PaXlF+5aB+PI7iuZat1JoykRO3A38Pnpdmf+3onzCKUUP/jPdoKBFH/yV/Xn1KCaL8TjOvv3D7F79wAtLVFsNmHdej+btxZQW+sayVvLZMHNNk6XhbXr/Kxd50cpRXt7jL1vD/DO3iH2vzNERYWDyy8vZNMm/7Rd23zhdlv4xB/V8aX/dZpvP9rGZz+/eEG4xDNFASq3o5O3A2+PSp43F8/NB6/tHODAviD3PVBB7aL5lQfW1xfntdf6efvtQaJRnfJyO3fcVcHGeXjjZ0JEqK52Ul1dyS23lnNg/xB7dvfz2GOdPPVUFxs3FnDVVUVUVMyfidll5Xbu/1Al3/tGG88+2cttd5XOdZfyj5DrPLDf5awrCUay/IW1eO5c09MV5+c/7mLFajfX31w09QGzRCyW4sUXe3n11X5AsXatj8svL6K+3kUyU4nSBYDdrrFlayFXXFZIS0uEXbsG2Lt3kD17Brj88kJuvrkMv3d+TH3ZeoWfQ/sC/OZX3axd7704kp1zdHmlV/++Bfhvo57+IubiublD1xXf+Y82NIEPf6x6Xkwl0nXF3n1D/PbJHgKBJJs3+7n11jL8/lEu7gJL98hEba2L2loXt99eznPPdbNr1wDvvDPE9ltLufyywjl34USEBx6s4uTxCN9+tJU//8KSCQdKLhwEcjR3UikVBkrGPNfLeS6eeyF/2jPmuSd7OXU8wgc+VElxydzHwZqbI/z715v4yU87KCiw8slPLub++6vPFbALFLfbwt13V/LpT9dTVeXkV7/u4l/+tZGTJ8Nz3TU8Xgsf/GgVHW1xHvt599QHLGSymTc5H6tYzFemW6ki0/D1cMpES1OUJ37ZzaatPrZc4UefohpFpjSKRIZgeqY0iljKOC4YTPLcM93sfXsQr9fCve+tZPX6QkQTghNUiklm6GemviTztIqINcNPZcY0igm66i1x8TsPLuLUsSGefrKL//ONZtas9bH9tvKRFYx0y+QfgCOjNT299AsBVq/zcu0NRbzwdB9rN3hZkU67yLjYyfR6MvfMA49kNAtaxPJFIqHz3W+04fFaeODByjkd5j98KMBjv2onHte56ppirt9WgtNpGRG4ixERYc1aPw0rvLz6ch+vvNzL8XeD3HFXBZu3FM5Zv+75QDlHD4f4/jfb+PzfLsXlXjgDK+fFPBMx052cgBef6aO9Ncbv/UE1Hu/c6LyuK55/tpuf/KiV4mI7n3xoCdtvK19QI475xmbT2HZjKQ/90VLq6lz8+pcdPPFYxzlTp2YTh0PjwY9XM9Cf5MnHeuakD3lHjKXtptpmE1PExhAYSvLsb3pZt8nL2g3eOelDNJrie99vY+dLvWzeUsBH/3AxZWXzJ7VgvlFYaOPDH6njmmuLeXP3AP/xzWaCweSc9KV+mYsrrilg5/N99HTF56QPeceiTb3NIqaIjeH5p3qJxXTuuX/K2Q55obc3zr8/0sSxY4Z7dPe9lXlffScUSpLIV0BsltA04Zbt5bzv/iqam6N87d8aaWuPzklf7ryvDE0Tnnr8ArTGJIv6+vO1xv7FQGAoycsv9LPlcj+Vc7Co7fETIb72b40EA0k++tFaLr+iKK/xuNONQb79o9P84GeN7Hqrl4FBw3JQmSaeznM2bCrgv32iDl2HRx5p4sDBwKz3oaDQxjU3FPHm64N0dV6A1piWxTbL3TFJ88xv+0gkFLffPbuZ10opXn21n//8zxb8fiufemgxy5dNPql4OnR0R3lldw9Pv9DOmeYQAO2dUdatLuBjH1qKJsLLr49PD9B1RSicpLs3SlNLiKPvDvH2O328/U4fR94dpLEpRHdPlGAoMWexqLHU1bl46FOLqap08IMftPHssz2zvnLRzbeXYLHKhRkbW2iVXRcqGdMvJog8Dg4m2flCP1uuKKCkwjFhOkXmRT2ml0YRTsCOF3rZuaOXlau83Hd/FQ6HhUgSYhnSISLJDOcbM3LZ2hZmx852FtV58Re72PlGD5GUcKo5zLZtVQRTGqvWl/Dm947THdRxOM4OHjz122aamoKTv4ExOBwaHq+NqioPNTVuqqo9uDMMRtgzuB6ODKkSU6XJ2Nw2PvwHi/jN4508/0IvA4EUd95VARnmaWsZbj7JYJ1OVHvd47dy9bYiXnq2j+13lVJROd6yV5kqasxXhHk3OnnBitj58txve0klFdvfM3tWmFJqRMAu2VLAe+6uRPJwgVRXufjgA8sAiOtw5Eg/DoeFYDBBJJLC5wOLRaOgwE7jmQArVhailEJEWLW6kLpFXhxOK06XFZfTgtNlRYBINEU0kiQaTRKLJolEUkQiSQYH4xw92s+hg30AlJU5qanxUFPjoarKPWsZ7Vabxt33VeL1WXllZy8A77u3fNZmXtx0Wwmv7ujnycd6+MgnMq03vbBQ86yyiCliQDiU4uUd/Vx6ZQFlFbNXUuXV1wbyLmDASFytqTnIG2/2YLVq+Hw2qqs9nD4doLzcBcCixT5aWkKsWHk216p+iVH0b6Ll6pwuKxQZFoY2xqpIpRRdXWHaWsO0tYbYv7+Xfft60TRYutTPunXFVFa6yHd6t4hw483GD9MrO3vxeTRu216W13MO4yuwct2NRgLse95bRknpBVCuR5g4C3kOMUUMePONQRJxxQ23FM/aOY8fD/Gb33SxarU3rwIGjFhVAEuX+lAKdu3qwuu10tISGtnP67UxNGQEomc6oGCxCFVVHqqqPFx2aRmJhE5HR5jGxiDHjg1w4sQQpaVONq0vYuWKAqyZ0vpnyLCQhcMpduzoo6rSwcaNs1OR9fqbinnh6T7eeGWQO++dHfHML7Mf85oKM7APvP7yILWLHLNWgaCnJ84Pf9RGRYWd+95XnVcBg7OCtKjOy4YNJaxbV4zDYcFi0bDbLRw40MfQUJwzZwIsXZafm9tm06ir83LNNZU8+OAKrr++Cl1XPPdiO9/89nFeea1zREDzgYhwx50V1Ne7+NnPO2htnZ30i5JSGyvXeHj9lYFZH1zIG2aKxfyiuSlKc2OUK6+dnekq0WiK736vFRHhwQ/XYM9US3kGZEqTsFiEnp4ohYV2rrmmgr6+KE8+2YzDYaGkJP9CbrNprFlTxAc+sJT771tMba2bt/f18p/fO8FzL7YRi+VnmVKLVfjQB6vxeCx893utBAKzkxB71bWF9PcmOXYkNPXOCwFzAvj84vWXB7BahUuvKMj7uXRd8eOftNPTE+djH62juNhOJA/3UVt7mBd2tHPbrTWUpkUpmdTp6Yly/HSQ7u4IVqtGSYkTn8/GVVdVYrNpU9bRzzUiQk21h5pqD4Fggr37etm3v4/GpiA331DNiiW5TTMB8HqtPPjhGh75ehPf/0Ebf/jx2ry6sgAbNnvxeCy8/vIgq9fOzSyQnCE5r+w6Yxa0iGXKycxmoZBEQufN1wfZsNmHy2NBV5lTMzLd5NlUo3jumR6OHg1xx10V1NZ7iKUyp1GEEhlSLCaYAJ5I6Oza3cU77/Th9doYjCgc6TaU0mjtiROKKRrWlVFe5SElwkAcQIP42QVGJmIqTyhTekKmwUjXcBqFw84ll1dRu6SAl3e08avHm1i9upCrr6o4J6b39RoAACAASURBVOUjq/5kqEykpaCkwsV976vipz9u4xe/6uLue41J/mMHJ0aTsRpFhtd0BRarxtYr/Lz60gCBQAqP13g/Ga/fyV+adUtnHDkSMREpBP4DWIfxlj8KHOM81528qN3J/XuDhEM6V16bfyvs0MEhXt7Zy5atBVx2ee5d17b2MD/56SneeaePVWuKuO/+ZZRXuEdeFxFWrCxi8xVVVFR75+UCHABl5W7ufu9SNmwq5ejRAX7845M0NuY+637tOj/XbSvh7bcGeXP3QM7bH8sV1xaSTCr27BrM+7nyimAsbTXVlh1fBZ5SSq0CNgJHOLvuZAPwfPpxRi5qEdv1ygBFJdaR2k/5IhBI8tivOqipdXLHXbkt7aPritdf7+SXvzyDrivuuXsxV11Tjc2+cKtdWK0aWy+r4H3vXYLNbuGJ3zSzY0c7yRzP77zhxlIaVnh4+skuerpjOW17LLWLnNQtdrLrlQUuYsOjkzPM2BcRP3Ad8E0ApVRcKTUA3IOx3iTpv/dO1daUIiYidSLyoogcEZFDIvJH6eeLReRZETme/ls06pjPi8gJETkmItunfEdzQDSS4t0jITZf6s978uOzT3eRSCjed391TidzJ5M6v/ltE2/v7WXN6kIeeGAZtbX5FeTZpKLCxQfev4RNm0o4dLifX/26MadBf00T7n1vFVar8ORvu/I+Z3TzZX6aG6P0901QyXIhkd3oZKmI7Bm1fWJMK0uBbuA/RWSviPyHiHgYs+4kMGUlhmwssSTwOaXUauAK4FMisoZJzL70aw8Aa4HbgH9PL1M+rzh2JEwqBWvW5zfQ2tQU4Z19Q1x5dVFOkx2TSZ0nn2ymqSnEtm1V3HBDNfY8Z8InEikGBqO0tQdobR+ifyBCPJ7K681vtWpcfVUFt22vpbs7wmOP51bIvF4r199QyonjIY4cze/o4XBpp8MHsp/GNe8YnnY0tYj1KKW2jtrGLpxrBTYDjyilLgFCZOE6TsSUgf20Gg4rY0BEjgA1GGbftvRu3+HscuP3AD9WSsWA0yJyArgMeH06HcwXRw4EjZWnG9xT7zxNdF3x68e68PmsXHd9ydQHZIkhYC00NYe48YZqVq/ObYwtFkvS1DrE4GCUcCRBKJwgFI6TmCTyb7VqeNw2PG4bbpcdv89OXW0Bblfuav8vW+Znu1bH008389jjjfzOvYsmDPhPh8uvKOKtPQM8/kQXKxrceRutrKy2U1Ri5fCBEFdfP39WzjpfcjTtqAVoUUrtSj/+GYaI5XfdSRGpBy4BdjH5cuM1wBtjOjtu4ljavPwEwE3X/t35dGPGKKU4tD/IyrWevNbqeuvtIZpborz3fVU5u+GSSZ1f/Sb3Aqbritb2AKfODNDWHkBXCptNw+O243HbKC1xGyLltuN22xAgHEkQDhsiF4nECYUT9PQZ1tne/R1UVXpZsriI+hoflhwUylu6xMf27YaQ/devm3j/PbkRMotFuP2OCr73nWZefqWfG7bl7gdnNCLCmvVe9rw+RDKpsC3EWUjDgf0ZopTqEJFmEVmplDqGscLR4fSWn3UnRcQL/Bz4Y6XUUIbg9EQvjPM30ublowA3X/fwpP7IdB2VTCHgttYYA/1JblvvGdd+xkoVGRfZOPe1aDTFk091s2iRi5XrCiatiZ+pGsXYNIrRFtiV11VTs6yIiZLcQxnaDCfPFZNAIEpzUz8tLf3EYklsdiu19SVU1xbj87vOnnvU5zL82Trc4CiBIs6dThcKxmhv6aOttZ+2jmZ22yzU1BZSV1dEQYHrnIGNmHXyb2oiw6+0xs91N9Wx8/lmfvKrJu5+z+JxQpYpxDmZlNYt9bJ6tZfnX+hj46ZC/P5zb41MSxpYMqRmjJ1yuma9l1d3DHDyeJg1ayf3AjJJ89zm/ec0I//TwA9ExA6cAv4A4yvK/bqTImLDELAfKKV+kX56MrOvBagbdXgt0JbNeWaLwweM2MfqdfmLh734Yi+hUIoHHyzPybSikRhY2gKrWTYzd6S3N8iRwx3094cRgfIKP5U1xZSWz3ygw+N1sHxVFctWVtLbE6CzpY+mxj7OnO7F73eyanUlFRXTn95Ut9g/YpE99njjhEI2He64o5x//ufTPP10F+9/f/WM25uIFas9WKzC4QPBjCI2r8mRiCml9gFbJ3gpt+tOivGz+U3giFLqy6NeGl5uHM41+x4DHhARh4gsARqA3efTqXxz9GCQqhoHRcX5Wa+xpyfOa6/1sXlzAbW1rqkPmIJUSp0jYDNxIePxJHv3NvPaq6eIRhOsWVvFLbeu4bLL6imvLMjpSK2IUFrmZ8vWxdxy62rWr68hpSt27zrDm7vPEI1Of5Ru2LXs6Yny2OONk8brzoeSEjvXXFPM3r1DNDdHZtzeRDicGssaXBw9uECnIGUz5WiWUxCzCVRcDXwYuFFE9qW3OzB81ltE5DjGUuRfBFBKHQJ+iuHbPgV8SimVn8lw0yCVUpw5GWHF6vz9Cr76ah8iwq235qZqwa7dXTQ1h7hhW9WMBKyjK8hLO47T2tLP8oYytt2wkmXLynA48j9xw263Ur+khG3bGli1upKurgA7drxLS+v086aWLvGx/dZaurqi7NzZnpN+bttWgsul8fLLfTlpbyIaVnloa4kRDs2b2yJrFMa0o6m22SSb0clXmFxbJzT7lFIPAw/PoF95o701RjyuWLx05hbSRMRiOvv2DbF+vQ+fb+bi0NoaYu/eXtasKWTNmum5kMmUzv6DnRw/2YvH6+CaS5dTWDQ3roymaTQ0lFNZ6Wfv2828uquZ+kUBLtlYhd12/i7h0qV+tm4tZc+eHurrfSybYRUOh0Njy5ZCXnutj6Gh5LjYWC6oT197jaejrF63APP65tlsj4suY//MKaMES75EbN++QWIxncsvn/kQeiKh88KLbfj9Nq65unJabfT1R3jmhZMcP9lLw7JirruuYc4EbDQ+n5Nrrl3GmpVlNDYN8PRzJ+jqnp6LtXVLGWVlTl7a2U4kBzPqL7usEF2HPXvyMx1p0RJjUv6ZU/lxWfNKbqcd5YSLUMQiuD0WSstyHw9TSrF79wBVVQ7q6mZe0uaNXV0MDSW48YbqaZV0Pnm6j+d2nCSZSHH91fVs3lg9oxyoZDJF32CQxtZuzrR00TsQIJ6Yvmhomsb6tRXcuG0pmia8+PJpjhwbv1jJVFgswo03VBOLpXj5lY5p92eY0lI7y5e7efPN/NQAc7ktVFTZR35QFxICaNrU22yysKtYTOOYM6eiLF7inHT+Yqbk80wRjIQOzc1R2ttj3HV3BUklIx1MZqp+McnYfWtbiP37+1i9tpjiCu+4kj2Z0igG4haaG3s5fKCNkjIfGy9ZhM1uZSAOgcTkV1hwVJvJZJJQKEQ4FCEcDhMKR4jHJp5faLPb8bhduN1u3B4XHo8Hm+3sj4TXOvkHENcFq9fH5deu4ND+FvYf6iScFJY1VJDpEx+r6Z5CF5s2l/H2nm4OLfezfOnEbmWmxUdGjw1subSIn/yolUNHgqxa7SPTVNTpTFhYvMTF0YPBc6runtPm+Tc5a8wzb3Jhi9j5Eo2kaG+NsWGzLy/tv7mrH7tdY8MMSx8nEjrPvtCOz2dj62UV5328IWAtlJb52LS1/rwSTROJBB3tnXR1dY9MJ7I77DhdbopKSnC53bhcRq5XJBIhEg4bWyTMwIARpBeBktJSqqsqsTuyy+i0Wi1suGQRInDimGFNbVl7fou2bNhURuPpAC/uaKem2o3LOf3Le+UqLz6flTd3D7Bqde6vl0VLnOx+bZD+viTFJfkZJc8L2c3vnlUuKhFraoyhVH7iYZFIikMHA2zaXDDjnKXXdnUxOBjn9rvqz9uNPHm6j8MH2s5bwJLJFK0tXXR2dqHrOiWlJZSUFON2u7BarcQnSPT12Wz4/IZg2zVFKpUiHI7Q39dPd3cPvT29lJWXsryuErtt6ktNRFi/aRFgCJnLqlizMvsRXk0Trt1Ww2O/PMmOnR3cfmtt1seOxWIRNm8tYOeOXvr743jKcys0w8H9M6ciC0vEkHlXxumiiom1NBoxiHzU0j9yOEAyqdi8ZWa1yQKBBO/s72P92iKqqs9v5Kq3L8xb+85PwJRStHf1smvfYdrbOygsLGDdujUsWbIYv9+H1Zr975zFYsHn87JocR3r1q+hpKSYrs5udu89TEt7V1bxpWEhq6ou5MChTto7zq+WWHGJk0u3lPLu8SG6u2cWc9q8pRCl4OCB3Nczq65zolmguTG/JYByzXyMiV1UItbVGcfl0vD5c19U4/i7IXw+K9XVMxPIA4eMIpZbN5/f/L1kUmfXnhacTisbNi/OSsAGh4K8deAYx0424XI6WL1mJUuXLcHpmrnIOxwO6pcsZu261fi8bk6caWXP/iP09g9NeayIsHZjHX6fg91vtxKLnd/gwaYNJVitwjsHZ5brVVhoo7LSwYl3c191wmoVSkttdHXkb3GUvCAg2tTbbHLRiVhZpT3n5nAqpTh1MsTyBs+M2k4mdQ4e6mdpvQ+///xmB79zsINAMM7lW2qxZZFv1drRzd5Dx0kkkqxuWMwlaxvweLKz/FKpFMlUdsLicrnYsHoZ61YtRSk4cPQkZ5qnTky1WDSuuLSWeCzFW/vazqvcj9NpYdWKAo69O0g0OrOE0uUrPDQ1RWbczkSUVdjp7lxgIkZOaiLmlItKxLo7E5TnYXFc4yLXaVgxs8TF4yeGiERTbNxwfutfdnQGOXGqj4ZlJVSUTz0ftLWjm+OnWygp8nPZptVUlBZPKb7ReITO/nbebTnCO6feZv+ptznWfJiOvjYisXBGkRERSosKuHTjKirKijnT0pGVkBUVuli7uozm1iEam88vs3/jhmKSScWhIxnLs09JQ4MXXYfjJ8IzamciyivsdHfF816MMZdkX05s9rhgA/tqTPmAREKnrzfBZVcVZKxUkWkG3mTHHTkWQjRYvMQ7rpoFQCJD9Yv4qBPuO9BHUZGd8io3cX3ixUCGGa5GEY8n2fVWC16vg4ZVVYSTWsY0iqPNvXS0tOD1F1C2aBm9CQ3SUxgHImfPp5ROLB4gEu0nEh0gmTLiSxaLC7ezChBiiQHaelto623BYrHjchThchbidBQgo3yKQtfwB2ehuGYJ0ZRwpqWDUFJj7ZLJR1/tmmLR0gpa2oO8ta8NX5EXl8v4EXJYJv+mHBaFv8hJdbWbdw70s3Z9ycic0EzfxUTVr6tq3TgcGkeOhli7duJRyulcMwCl5XZiMcXgQIqCwjHB/QyVMeaaeRbXv3BFbCy9PQmUMkz4XHP0WIi6OhdO1/RjbZ2dYbq6olx77fnV4D94sI1YLMmll00dyG/t6B4RsNoly9AmicBGogP0DZwkpScAwekowOmoxGErxGJxjNqzlpQeJ54YIJkaIBTpJhjuRBMrRYVL8LjGx/VEhOpF9QB0d7RxxqpTX1c1aZ81Tbhkcx0v7TjOvn0tXHHFkqw/n3XrinnmmRaamoLU108vTcJiEZYs83DsWGjSnK7pMuwVdHUmxovYPGa+idhF404Oxx5y7U4GAkna2mIsa5hZWZ+DB/ux2TRWrMh+dLOnJ0hrywDLG8opLMw8lWjYhcwkYLrS6Rs8Q3ffUTTNSmnRCmort1Jesgq3s2KMgBlYNDsuRzllxSuprdxKWfFKrFYHvf3H6ek/ga6Pj50NC1lBcUlWrqXH42D1mkp6uoO0tWXvVtbX+/B4rBycYYB/eYOHwaEkHTmOXw1fiwsqLiagWWTKbTa5aCyx7i7DZyrLcb7PiZNGrGQmIhaPpzhxYog1awqxn8cqRYcPteN222loyLyWwsBgYCQGVrZoYgFL6Sm6e48RiwfweSop8C9CS7uEQ/ISOK0YK6cKWnQDKUsT2DsRjBvRQgM2KcPlLMLpKGAw0MpQsJVEIozfuQqr5dzPfVjIXBbFmZYO3G4n5SWTzzetry+huamfI4fbqa7OTugtFmHt2iJ27+4mEEjg803vu1+23PhuT5wIUVU5XsinS3GpDc0C3V0LR8SE3FliInIGCGBMzUgqpbaKSDHmupMTExhMYrEwsnBprujoiKFpUF4+/Yu7vT2MriuWLMne5RkcjDA4GGHpstKMbqRSiuNnWnHYbaxpWDKpgJ1uNwSspKiBooL6EQEzEKyRLajwcsIhFwOptwjEuwmHnOjhpdgil2HjbFKqiEahv46y4lUkkxFOtR0hmRpfO0xEWLlsER63i5ONbaRSk0eXRITly8uIRBL0dGef8jDsRra2Tr9+l7/AhtdroSPH6RCaJvh8VgJDeVgGPl9kMTJ5niJ3g1Jqk1JquDiiue7kZAQDSbw+a87TK7q64pSV2rHMwIRubQ2jaUJFRfbVJZqb+9E0oaYmc32x9q5eQuEIyxbXTCh2wwIWjgYoLWoYF8dK6jGiMehP7SWgH0UnjltbgstagEInqL9LX2o3A7FTJPVzqzK4nIWUlawilohmFLLl9TXEYnGa2zOvCVFR6cdms9DcnP2IY3GxA6fTMiMRAygvs9OVh7UpvT4LwcDCqiuW5xSL3K87eaEQCKTw+nKf5NrZFaO8fGZxttbWEJWVrqynGKV0ndaWfioq/Njtk0cEEskkp5vaKfB5KCsZL3ajBWxRxXLcowQslhqkN3KYjvAu4nGwWJK4HBb89krcWi0Oqwe3O4nbqWGxCMFECx3hN+mJHCCaPCsyTkcB9ZUrMwpZUYGPsuJCmlo7icUmt3YsFo2amkLa2weJx7O78UWEmhoPra2hGaUyVFQ46OrMfTrEgrPEyDrFYqp1J8GY5/6MiLw16vXzXnfyoomJBQOpkSKF0158ZMyByaROb2+CDRsyT/jOtHB1KJKipyfKlq1lpMa0P1nF5ebWIPF4ioraYqITpGEMV6Noae4ikUxSU7ecUMoQyOE0Cl1P0d1nuJClRQ1gKaGz0zihbukm5TgGWNESNWh6MbFmO3FrAmk4yECTFT1WikpWAoKzuhGnP0E84idq7SSaOoAWW4oladSpr6gopKx4Fd29RznecpTy0jVYNOO7cKdXGKmsraGnf5B3mzpYsmQxAJ4JVqIqrynmzJleTjYNsnzpxPl0Yz/Hqmo3J08O0T8Qx1s2vR+c0jI7kajO0KjrKBumuta8Pgs93RdkTKxnlIs4GVcrpdrSK6U9KyJHp9Oni8gSS+bcEutJp22UT/PGAGhvN9ycmprsE2VPN/Zjd1gpKZ08hpZIJOjq6qK4uAi3Z7ybOjDUNBIDG22BDQuY6H6s4UuxJJYgKSOQrpI24gMlWD1BVNLGcMHfWE8FFnsQS2Ix1shWJFmM7jiFkrMunNNRQFnJKhLJCH0Dp8ZZNA6Hg/LyMnp7eolGJp/z6C9w4fE6OdWYvUtZnf5sW9umn7Bakba2O7ty61L6/AvMnczh6KRSqi39twv4Jcb6tJ3phYfIdt3Ji0bEgoFUzudMdqUv6JmkbbS1hrBahbLy7CprRKIJ2joCVNcWZVzUo72tA11XVNeMz8GKRPsJhjvxearOiYGNFjBLdC2CBUiBpN0dLYXN308q4kasZ60He2EvKmmM4AkaltgqLNE1iDpXmJ2OAgp9tUSifYQiPeP6VVlVgaZptLZNnnIhIlTXFdHbF2FwKLsJ3gUFdjweK20ziIsND9zkeiTR67MSjeo5WehktshFTExEPCLiG/4fuBU4yOQLEE3KReFOJhI60YiONwc170fT1RVHBMpK7Uz30m5rC1NZ6c56YKCpeRCloLp28qlJiUSC7u4eSstKcDrPncytlE7vwClsVjeF/rMr68VSgxMIGIgljq1oH1a/QkQR7yuhvO4M4gwxFHEDgh53kBhaN9KWoCGpifvn81YTiQ3QP3gGl6OQ0b+jNpuNiopy2ts7iFZXgXPilIiqmiJOHG3nTNMAG9dNXbZbRKiu8dDcNP2J3D6fBadTozPnImZ8zsGhFEUlC8OmyNHYWAXwy/RAmxX4oVLqKRF5k3ysO7nQiccM18XhyO1FEgqlcDk1bDaN+DRXHwsMxamtzd6VHBiK4XRa8XgnrzQRDBpB7NLS8Rnz0dggup6gpHDpOVODAvFmwHaOgAGolIt4z5VE2yKIpKheu4+S8jY6311DoGXZyH7O6uxSTESEQv9iOnsOEokNAOeKXWlpCe3tHQSCQSiaOG/M4bDh9zkYHMretSsqcnD83UESCX1apb5FjHSIXK9QNHxNxuILwxLLVZ6YUuoUsHGC53vJ9bqTFwLJdGTdasttekU0pmOfgTAmEjrJpMJ5HhVIw+E4HnfmpM1w2Ij9uFzjY2GR6AAiGk7H2YTRpB4lmupDS1SmBUxnbEh6WMDKKlqpcATYuGEXldVnMvRi8pvSbvNg0WxEouPjWnaHHYvFQjiUOX7lcdsJhbP/5XA6DWEOR6YvQg6HRjSWW7Gxpa/JZGL+zpU8hyxGJmd7AvhFImLGBWKdYLRrJsRiOk7n9D/C4RvKdR5zLkPhBG7XVCIWwel0jssLU0oRifWPm6AdShgxKC1pTMa2+k5gK97LiBCNErBye5BKW5Aqa4BNm16dQMgUtsL92AoOT9o/EcHpLDKsQqWPe83ldhEOZ14JyOO2EQpn79q5XMYPRXgGqyE5HRqxHIvY8A9rMtNCDPMMzTL1NptcFO7k8K9crkUsHtNH3IGxw/rZMHxDOZxWdDW+b2NTOpRShMMJqqt8GRcfCYfCeHy+cSWlI7EIqVQcq6OGaFSl29QJptqxSwmRVoXm7MVR2US8t5xoW+wcC6zcGqBSBil3DCE6SCrJpo2vsPu1A/TxnlFnsuGsaiHV60CPeIkVFaN0G4xyU61aIUp10RcK4nefO43I6XLT29NDQp98wrXLbSOR0InFU+Pqp030WTrS1m4kgyU21Xdod2gEQ7mNiQ2vPrVQRCyX045yxcUhYsOWWB7cSbd7+j87kRFLLLuvIRZLkdIVbvfko6HxRJJEIoHLNX60czBkrKNot59NfI2pHhRJnFJFBIVr0UlUykK0tX68gFkGqfIOUGUbwKKBVUtBGC77JOx75wwdbfU4PEOUlbcTQuFefBI3CpsjQTRuo2fQj9Kd2KUEp7UaEIZCA+NEzOV2o3SdUCSG1z1x7M+TtkZD4QSFBVN/B86cWWK5FRvrAnQn51uN/YtDxPJkicViOkVF059QPnxDZRsTC0eMGJAngzs5FDTcMKd7fDxsMDSA1eLBop0VwajehgUXNinEVtSEzTdIuHEZpLRzXUgxBKzW1k9FKoBF17E4065VGDZtepV9gKOyjfYjG1COKJ5lR4mkNCKBm3E4u6gtbmSofz1R1UVQvYvd5mcwNEBN6aJzbozhvg8FI5OK2LCQhyMJCgumLqedq5hYLMcVXoevyWSmjOh5xjzTsIsjJpbIkyUWG+VOTofhGyrbOmTDMSB3hsB+IGSImGuMiCVTCULRIHbbWSssqUIkCeDUqkB0XLWnSYa8xHsqqVq9n5LyNsrsQSpsQSq8Q1TZBqhMBSix6fhFoyI5RI1zgArXEJW2ABs2vYbXGSIRdZEYKCEZt5Kw6IikSMRq0SxBNLHh1urQVQy7rYhYIkYscW6+l9NprAs6/F4mYnhwI5xlXMxmM6ZHhcPTt8QcDo1YXM/p1KMREVsolhjzrzz1RWGJDSPk9tM1iuTN4Ph00EvLspHheyeTOT+8otBE+yyramAwetZ1U+nFaQUroCOWFCpiiIPVEcUiOpqmsEoKSelYHAqbSpGyJdFsYItbEBSS0rFqKTRRkLRSUn+SUF8pdmsSXYFOHAsulLKiqzhR1Y0mjpHBhbHB/eH+Z1odafj96VkaMCLGUmMz0R+R7M93Pm3C/F4sdywLzhITkW+JSJeIHBz13BdEpFVE9qW3O0a99nkROSEix0Rke746fj7ky2R3zHC0yp2O00Si2VkHw6OSw27lRPi8RiwsGjnXirFabBR4i9BkVHAdHxZcRPV2BBvRtkXYCvqxFfTRcmAL/QMl9ES9dKa8dMd8dEQK6LD4iSU0ErrQofnpCPvpSXrpSPg4/u4GDu+5DhQU1h9H1yAe8BFIHqNff4W2PhuDqYMolcSnrSKeGMBmseGyn2s1xmMxdF3H55ncTRxxradINxkmmdBJJnXcM6i+O2x55zImlK9QR74QAas29TabZGOJfRv4GvDdMc9/RSn1pdFPiMga4AFgLVANPCciK5RSczo5LF8m+4xFLD0oEI2k8GVRSsw94kIl8E9SgcfvMUQsEg7j8WYu1CgiOLUqQvopEipArKsKe2knrrpTDB3aTNM7l8HG3cbOLiDt9cVdFkhBd9hLe7yADlXEkWOXcPL4ejRLAoUiZk2goagu60Y0B4l4OcFgNTEtgRUPgoV4YpBi3/hFSiLpPDe/d/KpWNm41ue0mf6hcGc5iDIR0RmGDyYiX6GOfCFgWNzziCm/EaXUTiDb+r73AD9WSsWUUqeBExiTOueUfOXizFTEhkclo1laYi6nFU0kYxzIYbdisVpHhGAqHFIBaET1dkAj3LQMzRHDWdWCnrTT9M5l9A6U0Bny0UkhbdECmuMlNMdLRgRs/8/8nDy+HoDqNe+gO6LoQHDQx4nXbmao9yaiccFb+CaQIKw3EYifQqkUBZ7xahwJhxGRSYP6YAg5ZC9iw0uuzWQ0eaZ5gRNxNodx4YSn51uy60xiYg+JyIPAHuBz6RKyNcAbo/ZpST83jnT9oE8A3HTt382gG1MzbIklci5iFgYGjJtpOjURh12baCQ54a/buFiZCC6XlXAkwaTehwget4toJIxdG9/m6JvQiOnZiEUrCCU7cVTXI1SSCvfgqGxBLHWoVAFdndcBO42r0wV6xI6g6FFujhy7hCbnepzVhtVk9UTRrQn0WAkpHZzVLmPOYVixyGej2F2HrjTa+vYAQonXh2VMP2ORME6XE3uGirXhSAKrRcNpt4yL0Uz0WcbS7memdJapvsNYhKxYpwAAIABJREFUTMdhz7GILTB3EubfaOB0+/MIsAzYBLQD/5R+fqJvYkLlUEo9qpTamkXNoRmTP3dSZhgTS4tYlpYYgMcz9XQbt9tNNBJFnyIKLSKkUnEcUgzo6NYOABKBBlAWbMX7QIuh646zFlnET0fKR3vKP+JCnm0wwVBCQySFCtUh6CAxdAbxewZQuhdjeEVDJRROhx/LmPRuI6E3gnuCFJHRhMMJ3G5b1vGpEUtspjGxnFti+ZkSly+G3cmpttlkWpaYUqpz+H8R+QbwRPphC1A3atdaoG3avcsR+XInnQ7LjERseNg/Uxb5WNwuG13dmUvKuN0ulFJEI9FxtcSU0glH+kgkw+h6imQyis3mQikN3daClqxFdAfxvk3Yi/diL9lDvG/ziGvJxt3E/RZaD2+kr3PUV63FsRftJaUUtmQpxZXvYLEGgRZSSiMRdREJXo5FsxJPBlE6uJzjJ3gnEgmSySRud+bSRKEs5pCOJhIZjonNTMQKC3I7oD98TdoWkiU2z7o6rW9ERKqGS8gC92HUAQKjFtAPReTLGIH9BmD3jHs5Q+xpFyCe43lvLrdGJJIiNZ05RxiWkMdjY2gw+6ksfp+dM00DRKMJnJOUqvF4jaoYff3940QskYwwGGjB6ymn0L+I/9vee0fHdV15ut+pnFGFjEIGmDNpihRFBUqWrGgFy5blILuD7J5uu8P0mzXPXv3WdLtn/F63X4/nud32rJbDWFa7ZXtadkuyIkWJkihKpEgxRzAgh0IGKqDieX/cAgmAwEWhUIUqEPfjuouFW3Vvnbrhd/c+Z5+9wxE/kagfY8RFTDdE3HQefWQFRD1EBjZjKjyCuXg/MlJN2FdJ67Ft1Gz8EGuJD3qqAYmppBtzSQuIONHBjYTDJfgnpAsrK9NhBtApltbIqPJcs5qvFbHBAWVWgcMxc2aPWDTOyGiYkuLUs3+MjEQxmnRpZbAYJxCIz2meayqMZ1gxZthNzRZCMHNXRo6YVcSEEM8Cu1ByZrcDfw3sEkJsQnEVm4E/ApBSnhJC/Bo4DcSAr+V6ZBIUETOZBH5/ZptSUmwikYD+/gh2T2pJDadSXmGjpXkk5cKsNVUFHD/lo6t9gPpl01fPNpvNFBZ68PX0UlZWitF4VexMRjsOexnBUB9Ws4fQ2CCxWIjiglX4hlpJmDuJEcQQ2QRRD+G+HRhd57BWtWAq9hFqbaT5ox0A6O0jWKsvYrAHiIcLiQ2vQsbtSCSIEEJe6xIGQ/2ExgYocFZjMJiBq+ckHo/T1dWN0+nEbrdPem8iPV1DxOOS2qrUa3R2dgSoKLelHR4RCMQIBOKUzCOL73T4RxULMRv1H7KFWISjk5+TUlZIKY1Syiop5U+klE9IKddLKTdIKR+cYJUhpfy2lLJRSrlSSvlKdpufOg6nIeNpgMezffrmkSivstJGOJygvz+1LKVOh5mSYhud7QOqkePeygoSiQRdXVc8/yufdzkqMBqs9PSfJp6IUORZRiIRRUaMEHaDMUDMdEwRo7iV6OAm/E1rAIljxUlsdeex1lzEueo4OmOUwKWVRAe2XBGwhLGZmPXIpPTUANFoiIHhy5iMDlwO7zVt9vl6icViVFbNXBEcoLN9AJfTTFFhag8Ovz/K8HDkSprqdBg/x/MtCjOV0dE4RpPIeOhGtlD6xDI3OimE0Ashjgghfpf8u1AIsVsI0ZT8f+ZipEkWx5HLAE6XPuNVZcafyvPJ9un1KjdWZ0fq+d8baj0EAxGGB2fexmKxUFRcRK+vl3BYSR6oRKxLwpFRQKDXGdHrjAihwx/oQcoQxtg6dOEqMPqJWQ6TMHQjiRMbKWT09BZCHTUY3QOYinyMdVcycmoL0cESRbz0vcQtJ0iYOtDFymGCJRaNhujpP40QOoo9y66xiGKxGN1dPRQUuHCoxLcF/GGGBoPU17pTtqrG01LPpY7BVK6KWOaK58LkAjaLBV0Kyxz4c+DMhL/nXHdycR29eeBw6hkZHp9mkx5TnzBWiw6P20DvLPUI1UKAXE4jbreJzo4AmzdNzsQ6U/dNQ42LQ0e76OkYoLzkWpfNYVAsruU15RwcGKS7o4O1K+oB8NhgOBAFq5kybwPRWASjQWLRl3Mpcp6CAj8WcwOjYTN+0UlcfwFhbsaxogyLrgK9WIUM1wEJLC4bDmeYsUQXYbpJyCh6YcFpbMRu914RGYs+SGfPGXQCGryrsZjMjA9aj7f1fGsn8Xic5XXeK+ss+mstzcudSsji8jr3jMdnaqhEV2cAi0VPSbE57Wjy3t4IJpPA4zbMzdKY5X1/FgrYZBNB5kYfhRBVwP3At4G/TK5+CKX7CpS6k3uB/1NtP0tIxAx0tme++GlpqXle7iQoFsL588PE4zKlXPtGgx6vt4DOziHWrvPOGChpMZuoqSyjua2LoeFR3AXKtIACuweXzU3/iA+3XUkPbTFZcdorCI0NYDY5cZorcZi8RBLD+KOdhGKdhOIdGIUHi6hAoGNMNhOR/cr2+kIcRi9mvWeShRSNhujsUR60ioBd6wL6gyE6e/qoLC/GoTIqKaWkvW2Q0lLnrIkhJ27T0RHA6x3vD0vvBvT5wpSUmFSLs6TD6EjmC9hkmxQ79ouFEIcm/P2UlPKpKZ/5/4D/DEycrzKp7mSynJsqS8addDj1jI7EM178tLTURG9v5Mpk7nSorLQTjSbo7VXPZjqR6hoPsViCrq5h9c9VlGI2GWlq7pj024UQjIUDdPa3XlknZQKd7mrslRACs95NkWUNHv0N2HQ1xGWA0cRpRhInicphrKIKj/4Giq3rsBgKrxGwnn4lw+tMAial5EJzOwa9nrpZ+sJ6e/2MjUWprpm1m+QKIyNR/P7YvFxJULoMyjLsSoJiiS0md1Kknp66bzwONLk8NXk/4gHAJ6U8PN82LRkRc7kMxGKSsVBmwyzKykxEo5L+gfn0iykuYescqvEUFtqx2UxcvtSnKsx6vY7G2koCwRAXpghZZUk98USMlu4mugfaCQR7MRqmt4T0woxNV4tHfwNO3WqculUU6rdj19ejF9dODwqE+ujuUyJvZhIwgNbOHoaG/dRVV2A0znwzSym5dLEXo1FPWZl6seKJtLSMKr91HiIWDMYYHo5RNo/SfNORSEhGR+I4FpkllqFg153Ag0KIZuCXwB1CiH9Bqzs5M0UlivuR6bqBjQ2KAF26kH5NQ6vVQG2tgzNnhlKOORNCsGpVGcPDIVqa1ae2lhS5qSwvoaO7l56OtklCVlu2DIetAJPBjKegFpv1avWhWOxa91sIHWZdMWZdyaQ8/Ve3GaNvsIn+wQsYDVbKi9fOKGB93V1cbu2itNhDZXmx6m/o6hqht9fPihWl19QOmAkpJadODVJaasHjSd+KunRRObcNDeqzCObK8FCMWExSnOGwjWySqdFJKeU3k9EOdShJI96UUn6RNOpOLhkRGy9w29uTZm21GSgqMlFcbORiU/o1DQHWrSskGIxx6dJIytt4K90UFzs4c6ZLdVK4EIJldZVUlpcw0OubJGRC6ChylVLoKplUASkYGqDTd4TOnqMMDjcTiQ4jp8n7JaUkEh1laKSVLt8xOn1HCYb6KXBWUVa8FoNh+kncfd1d+Lo6KC32sHpZrepIYzgc4+SJDlwuC3X16mI3kY6OAENDEdatm7lGZypcbApgs+mprpo9g+xc6O1Jjnhm2MLLNhkenZzK3wF3CSGagLuSf6uyeJzxeVKcjO/x9WTWEgNYucLOgYPp1zQEqK62U1Bg4uTJAZYvTy2IUwjBxk1VvL33PEePtLHjpoYZxWBcyEJxwUCvYqGXVVZP+vzEgFuTyY6noI7Q2BCjgR6gGyH0mIwFmIxuBDrC0cGkuMUAgdnkxO0qxWbxzChecFXAXJ5CVi+rURUwKSXHj7cTicTZfmP9nDrWT5wYxGLR09iYuvt5zfcnJBeb/KxYbst4p/74tViy2EQsw8GuUsq9KKOQadWdvG5FbGpUscUiKHAb6POFVc1dNQmaabvVqxy8t3+I9pYAy1dcG+OkpmumKzsVbFzv4Z19PQz1hygtsWKdJsRgnLBBsYpsLgNbNlZw8HAHbc29rFpefE2Vo4msb/RyQS/x9Q+wskbp9B/Hlhx26urqxqzTUVlSjF5fQjwRp98/ynBgiJHAEKMRxX016A14HAUU2D0UOZwY9BMvp6uBxeMhE1JKWjt78HV1JS2wGlwmyUwjhjZDgsstg3R3jbBhbRkVRWbGy8ipHRuTDkZGIrS0jLJ1SzG2CVN6jNNk9hjHMM17nd1jBAJxVq+yz3j+07lmQOnaMBgEhUX6vIuCn4nxpIj5xHUrYtNRWmbEl2F3EqCxwYrBIGhqml7EUmXNKjfvf+Dj2IlB7roj9WlMdTVuOjpHOHGqh/JSB1hn7sQet8hqKssmCRgoItPT7WN4aIR4In6lgvhA3wBCb6S6pBZRWkcoHEQisZntV6yo6QRgIoHgGBea2xkcHp3iQs68XSAY4cixLoqLbKxckbobCXD8lFKYd/3a1Ecyp6PpvNJNsGLF/EY3p6O3J0JxqTHjFl42EeRfH1S+tSerlJaZsuJOGo066htsXDg/v34xs1nPqpVuzp0fTjllNSjCtHVLJUajjgOH2knEZ0/BM1XAxsIRfD29hEIhCos8uFxOhBC0t3UwMjLKQH8/zRcvAmCz2LFbHClFzMdiMZout/PhsTOM+IM01lXO2gcGiqAePNyBBLZvrUq5DoHynQlOnRqiod6J05l+NSqAC+cDeCstOB2Zf977eiKLrj8M8i8Vz5ISsZIyE/7ROMFg5uekL1tup78/Sl/f/ERy43oP8bjk2PHBOW1nMRu4YUslQ8NjnDnVOad4uGgsRnuXD6PJiLeygrGxMZxOBwP9g0QiEeob6qhftgyz2Yx/dDSlfcbjCbq7ejhx/BQd3b14y4rZvnk11RWlKYnfhfM9+HoDbN5QjsM+txv91JkhxsJxNm2YX4d+MBinrS3EsuWZt8ISCUmfL7ro+sMg/zK7LikRq6xWhtnbW1KbbD0X1qx1IgQc+WhoXvspKrKwfJmLQ4f7GBycWzsrK1ysWlFMe2s/Z052pCxk7V29jIwGsCZLpUXCEdxuN319fXgKPeh0ySrn8TiR5DzMmfYdjcbw9fg4eeIk7e0d2O12tm5YyYqGakzG1KyiC+e7udTUQ12Nm/raubmDfn+U/R/4qKq0UemdX0jE0SPDSAlr16VQAGGOdHeGiUYlVdWZD6DNJuPuZBZHJ+fMkuoTq21Q+plaLoVYsTqzT1eXy8iKlQ6OHB7m9jtK5pVueNct5bS3B3h3bycPPDS3EbkNa8sYiwmaL/UCsHpd5ayWT7GngEg0ysWLl4nH4+j1Ovx+P6HQGG731ZHSwOgopeXlgOKSxqJRAoEAoWCQcChIMBAkGlX6HB1OBw2NXpxOx5WO/VS4cL6bi+d78FZ5uOFj3pSstnGklLz3TieJhOTO2+e27XT7OvzhEFXVFsrLLWS6qFrLZeUBVdeQ2bCNhSDfuvCWlIg5HHpKyoxXLqBMc8M2N+fO+jl7ZpR169Mf1rfZDOy6tZxXXu/g5PE+NmwqSXlbIQQrVivTd8aFbNXaSlUhdDpsrHTUMBqF1pY24vE4vb4+nM6r/V5joRASsFivDjgMDw3R0apMW7JYLDidDmw2Gw6HHbvDrioi8XhC6bubMHQbi8Xp6hjEW+Vh3cZqdGJusyuazg3R0e5n1y3lFBTMz01rvhykry/CI4+qT4VKl5ZLIaxW3aJzJ4XgmpoIuWZRi1g6D4S6eivnz82cwkbt4a02OcSog1Ur7Hg8Rg4dHGTzxqsipmaUmWcIFVi30sn5C04+OtTL8gYnhYWT3Y6o6v2dYNvGUsx6ybmmPgIjQW7cWoXTaZ62eMiVwNdYDBEPs6a2HINez8U2H6XmGEII+rv6KCt0UmGNXYknKyx3UOVehsNuwWKYenSuNnBiNoqhoSAtLQMMDQZZt96L2zPB5TPBJ3Y1YDTp0YkEdhULzjrlvdHRKAc/6MbrtXHDJveMIQtq52JiKMyhg0NYrTo2bXBi1Kmf+3QMvpbLIWobLDM+XPLM2JlEvvVB5Vt7sk5do4XhwRhDg5kPtdDpBNu2ubl8OYTPN7+MGUIIbr21ApNJx543O1SrYc+0/ab15ezYVo3fH+G1Ny/QdLF/2r6s8erYNouJVfVeHFYLdquZSCzOkbOtnLzQTiwep9Y7OczBbDLidtkx6FOb+9fdNcyhQy0YDXpuvmUZRUXXhqOYzYY5jUSCIsJ731bcyDvm6UYCjIzEOH16lI99rGBe6axnIhJO0NUepq4hvWzAuWQ8FY82OplD6uqv9otlg499rAC9Hg4cmF8HPyhu5a23lOPzjfHRR32zbzANNVUF3H3nMkqK7Hx0rIsDH1xmbEylgrjtqnWwflklhQWKSK2orcCZLMybqkiEghHOne7ko8Ot9PaOYjQZ8Fa4qa7xoNfrCIUyE+5y+swQra0BdtxYNm83EuDQoSESCdi2bX4xZjPR1jJGIrE4+8NAG53MOVU1ZvR6aLmUnX4xh8PA2rUujhwZnlclpHGWLXOxfJmLgx/20taWXhyazWrk1p21bNlYQf9AgL1vnaetbXDGEcZxkTKbjNR5i1lVX4HhGndRnd6eYQ4fuIheJ/BWujl9WslgbjDoOHOmm337LnDsaDunz/oYGFQeKOmkSfL5Qrz7bjdVVXbWr5+/6MTjkg8/HGLZMhvFxdnpr2q5rPze2vrFZ4mBJmI5x2jSUV1npelc+lknZmPnTg9jYwnefrt/3vsSQrBrVwUej5mXX2mjoyO9dgshWN5YxG23rcDhMHP0SBvv779ET8/InF3VVCgsdrLj1pUsW1VBebkLi8WIwaCjsMiOxWxg1apyNmyoIh6XnDjdc6WNc6G3N8QLL7Zgsxm4687ZR2FT4YMPBhkejnHTTfOLMVPjwrkgRSVGCtyLr0taAEYhZ10WkiUnYgCr1tppvTxGIMPVj8aprrayaZOLffsGGJhHnrFxTCY9Dz1Yi8tl4ncvtaYtZAAOh5mdNzeyYWMlo6NjHDzQzBu7z3DmdBcBf+asU71ep4RqjI7x9t7zBAMR+vv8FBRYWb+hkuJiBza7iYZ6D/F4gmBobn2Uvb0hnn+hBZNJz8MP1WKzzV8Q/P4Ye/b0sXy5nZUrMx/gCkqdyaYzQVatzc7+s80ckiIuGEtSxFavdyAlnD01v2lCatxzTwk6Hbz88qw53VLCZjNMErLuzvSFTAhBbW0Rd31iNVu31uJ2W7l4sZf9b5/j4HtNtLf2E41mRuBNJj0NDcXcfMsyhkfGaG6evO+2jhGKC23YrMaU3cn+vskC5nJlxu3bvbuXSCTB/fenNqsgHS41BQmHE6xZn/4c21yTbyKW9/ZsusdDTZ3r6i3YHXpOHw+wdfvktDdqJ0AtPmZqNECh28Adtxfx6mt9bL3sp3HZ9E9etV6zqV6e1anns4/U8OvftvLW6y08cH/NtBlL1QbUzPrJ31hU72RNvZPQWJSm5mEuNQ9y5kQ750914HZbsNuM2KwmbDYjNqsRm82I3WZCCAgEowSDUYKhCMFglFgswdbNk0uxSbOg1uMGYGWdi4vNQ1hklLbL/bR2jGAx6ti8pgS3afzHXv3RU8MoQLHA9rzSgtmk57OfqqHAZWRqIOp02139/dO/19kxxqFDw9y804O33HTNPtXOvdo1M/Wt08f9GAyClavtaVsQuQy/EFxbiCXX5L2IZQOdTrB6nZ0zJ/0kEjJrWQRuudnDh4eGeeWlHv746/UpFQGZDbvNwGOP1PCr37byu5daZxSyuWK1GFm9ophVy4sYGAzR3DrM8MgYg0NjtHeOztpvpktWM59aBHji63hcotcJClxmunw6tm/xUlWc+gjdRBdSEbDMWGBSSl55qQebTc+dHy+afYN5cPq4n2UrbZgtOjI9C2ChyMTtIoSwAO8AZhQd+jcp5V8LIQqBXwF1KIW5H5NSqk4kXpIiBrBmvYNDH4zQ1jKWtVEio1HHA/eX8vTPOzh4YJAdGeostiddy+dfaOF3L7Vy/301VFVlpo9FCEFRoY2iwqtBqFJKQmNxAsEIwVCUQDCKlBK7zYTNakxaa4ZpXbBwOMblrhGaW4eJRONsWqdULV+1bFwsUruRfT6lE3/chVQssMxw4vgIra0hPv1oGVZr9vLd9/dF6O6KcNNt2QndWAiU9NQZEd8wcIeU0i+EMAL7hBCvAJ9CqTv5d0KIb6DUndRKtk3H6nV2hIBTx/1ZHepes9pO4zI7e9/sY/0GF44MpXSxTRCyF3/Xwo03lrFpY2FW+nKEEFgsBiwWAzPZKTN9rcGgIx6XrF1ZTGXF3CdSSyk5c2aId/d1Y7UaJvSBZcaKCYcTvP5qL16vha0fSy2jbrqcPqH0Y65Zvzg79UE5z8YMXGJS6QAd75Q2JhdJGnUnl2THPih1KOsarRw9NJrxMm4TEUJw732lRKMJXny+O6PfZbMZ+NQjddTWOtm/v4d/f76F4eHM50ubD3q9jtXLi9ISML8/yksvt/HW3i7Kyqx8+tH6jHXij/PaKz78/hj3PlCa9eSERw+NUFJqpLR8cc2XnEqKHfvFQohDE5avTt2PEEIvhDiKUtFot5TyAFPqTgLzrzsphPipEMInhDg5YV2hEGK3EKIp+b9nwnvfFEJcEEKcE0LcncIxyRnbbyqgqyOctQnh45SUmrnzE6WcPeNn3zvqlYnmitms5957qrjjdi99fWP88lcXOXVy+ulFiwUpJU3nh3j2lxfp6Ahw884yHvxkZsIoJnL40BCHDw2x85ZCamoyW8loKn2+COfPBNm20521kc+FIsVpR6p1JwGklHEp5SagCtgmhFiXVntS+MzPgHumrPsGit+6HNiT/BshxBqU8ktrk9v8UAiRt0X1tmx3YTQJPtg3/ylCs7HjJg/rN7jY80Yv585mNrRDCMHq1W4+93gjXq+dA/u7efnFZkaGM1/xPNsEAlF2v9rKu3s7KCoy89nHGti4sSjjVlJrS5CXXuymcZmdOz6eepaQdDnw3jBCwPad2XVZs8346ORsy1yQUg6huI33kEbdyVkfbVLKd4QQdVNWz+S3PgT8UkoZBi4LIS4A24D3Z/0laaD2QFM7juP3g92mZ8tWF4cPjPDpx8swmXWo9VmqnRyzyuMgkRzW/8yjZQz0hXnuf3fytT+pobTUfE0YxSRU+q6nu6ftbgOffaiKo6dGeHtfN//+3EV2bCtl/TrPlYnMM4UYAKiVvExI9StTrbNX7biN1/CIxROcOTvMvv09JBKS224uY9tmz4x5+FXDKFTCIcx6ydBQlF8924HHbeQLn6/AZgKQqudQ7TfMdq0lEpIP9g2xep2doqKrJ1X1+s1jYy1Do5MlQFRKOSSEsAJ3An/P1bqTf0eW607O5LdWAm0TPteeXDfdj/jquL+cZhsywo5b3YyFEhw5nFra5flgMul44olKjEbB0z/vIBjK/IwBIQRrVrt54nONVFXaeXd/Dz95+jzv7OvOu/4ygFF/lPcP+Pjp0028ubeL4iILX/hsI5s2FmXF7YpEEjzzTAfRqORLX6rElsXRyHHOngowNBhjxy3urH9XthmvdjTbkgIVwFtCiOPAhyh9Yr8jD+pOTnfVTftYTPrITwHceeu3c9aBs2yFlZJSI++/M8T2m7Jv6rvdRr74BS8/+nEbz/5rJ49/sRpdFqIHHQ4jD95fTWdXkOMnBjl2YoAjxwaorXWwbl0h1dXqSQuziZSSrq4gJ04McPnyKFJCfZ2DjesLqcliu6SU/Ntz3XR2hfnylyopLV2Y1ND73x3C7tCzftPijdIfR3En53+7SimPA5unWb9gdSd7hBAVUsquKX5rO1A94XNVQGea37EgCCHYcYubF57rxdcTwbMAJeXr6mw8/FAZz/2mh92v+7j73rKsfI8QgkqvnUqvHb8/yolTg5w4NchLL7VSUGBi7VoPdXVOXC7jggia3x+lpWWUkycHGRgIYzbr2LKpiPXrPBkLXFVj3zsDHD8+yj13F7Nq1cIIin80xvGPRrn144VZyU2WC/LtV6QrYjP5rS8A/yqE+C7gBZYDB+fbyGyzfWcBv/ttL+++NciDj2VHUKZyww1uurrD7N8/iMNpYOfN2Y0UdziM7NheyqYtJVy8OMLJkwPs39/D/v09OBxGKittVFQqgmd3ZCaQNBSK0dkRoKszQEdH8Io7W1RkZteuCpYtK8Cu1hGVQY58NMSbe3rZuNHJbbdlL0PFVPa/M0Q8DjtuXtwd+uMowa65bsVkZhUxIcSzKJ34xUKIduCvUcTr10KIPwRagc8ASClPCSF+DZwGYsDXpJTZSRWRQdweIzfcWMC7bw2y6+4iXAULEwN8/32lDI/EeeM1JRd+toUMQK8XrFhRwIoVBQwOhunoCNDREaC52c+5c8MAFBSY8Fba8XhMWCwGrFY9JosRq1WPxWK4MlKYSEjC4TihUIzIWIxQKEZoLM7wUITOzgCDA8roqMmko6LCxtq1Hior7RQVmRfUlT3y0RAv/Hs3jY12Pv1o+YJ9d3gswZ5XB1iz3k5l9eJMgDgdi07EpJSfm+Gtaf1WKeW3gW/Pp1G54J5PFvPhB8PseaWfRx5fGGtMrxd86tNeoJM3XuslGIjz8btKFqwitMdjxuMxs25dIVJKfH1hOjsCdHYEuNA0THSGRP5msx4EhMemfz7pDYLychvLlxfg9dopK505l3w2kVLy3rsD7Hmjl8ZGO5/9vDKoslC8++YAfn+c+x/OfgjHQiGEzEifWCa5bqcdqV2q0xWRKKswsm2Hi/17B7nr3kIK3Ne6VGqnTu3EmtVuYDN8/vEKXnhRx/59A/T3hvnc4xVYLHp0KjasmhOmFkYRTczcFkeFmYYKM2xVRG1sLE4wNL7ECIXiBIMxgqE4UkpsVgM2mx6wj8OSAAAW6UlEQVSr1YDNqr/yt8WsnyRaymjV9G0yqmUGUTtsKr/RrJdEIgme+003x46NsmGDk898uhyjUcwSRpFeporpxmXGQnHefHWAtRvs1DdOX/ItzwyalBCkPPq4YFy3IpYO9z5UzMH3R3jj5X4e/Xz5gn2vXi945OFyKiosvPBCDz/4YStfeqISV2HuCqsKIbBaDVitM8+XzFeGh6P8/JkOOjvD3H13Mbtuy86cUjXe3jNIMBDn/oeLZ//wIiPf3Mk809TcUlJqYttNBezbO8RwFqohzcaN2908+YfVBINxfvDDFprOZy9p4/VKa2uQ7/9TC319UZ54opLbd2Un3kyNUCjOm68NsG6jY1FWNFIjGxH780UTsSnc/UAxUkpeej696kLzpaHBxte+VovbbeQXz7Tz3r7FPQ9yIfno8BA/+0krZpOOP/njGtaszk1c1u6X+wkG4tz70PVnhSFSnju5YGgiNoXiUhO33FHIB+8O0dqcnbJus1HoMfInf1zD6jVOXn+1l7a2a9vhH43y/HMd7Hu7D/9oLAetzB8iESVDyPO/7aa2TnkIlJXlxhXv7Ymw9/UBtt7ooqbu+rLCxtGlsCx0ezSmcO+DxbgKDDzzo04ikfmXXUsHk0nHY497+cITVddkWJBSkkjAqjVOgoEYv/pFK20tM1c1v16RUnLyxAj/9L1LHPpwiJ03F/LFL1Vjs+Um50A8Lnnmx50YDIIHPz1rBplFyXicmJZjP8+x2fV84Q8q+OF323jxuV4e/dzChFxMRQjBipXXukRCCFwFRlwFRlaudrH3DR+tLUHcHiP73u7D7jSy8WNu7BlKwJiPdHeN8fJLPbQ0hyivMPPoZ7zU1o2LfW7c7zde6efyxRBf+ooXT2HmMs/mEwL1EeVcsKivctUwCpX31MzP8WQNa9c7uPUOD3t3D7B+k4NVa+yoPt9VbVq1ofuZW6rWt2DUJXPW6wXB0SgrVtjpbg1APEFfd4hf/3yYdetc3HRzISbT1cbF0jQs1TJcQPqduQaVG2K6WTqBQIzdu/v48MMhrFY9n3qkjG03FCRDOpR9pRtGoU8zi4VOKFW9X36+ly03ONm2w3UlC4XaZZHu9Ztr8m10clGLWLZ55LFSzp4O8MyPO/mr/9aAeQEyHszGxMImg4NRjh0bob8/gtdbxJEjI6xY4WDNehdSSrq7w+iukw6DeFxy8OAQb7zRSzic4KYdbu66szhnruNEotEEP3uqA6fTwONfqlj0SQ/VEDlwF2dDEzEVTGYdv/cVL//w7WZ+/S/dPPGVabMKLSg6nSAYjPPq63309IRpbLTz+ONeDAYdnZ1jXLgQwF1owltpoaJi8U918ftjHD48zMGDQwwORmlstHH//WVUefMnxfMLz/XS3RnhT/6yGrsj96KabfLtuaiJ2CzUNli555PFvPx8H+s2O9m81ZWztiQSkn3vDfLhoWE2bHBx992lmJO+UzwuefxxL6dOjbL7NR+fuLd00YqYlJLWlhCHPxzi5MlR4nFJfb2V++8vZfVqx4yJEnPB+TMB3np9gFvv8LB2ERfEnQv5ZmhqIpYC9zxQzMljfn71dDd1DdacdtomEhCNStrbx7h4McCaNUoBDr1eYLXq2brVzeBwjPNn/ZSXL+xE6/kSCsU5fmyEQx8O4etRUvVs2+Zm2zZ3zkIm1Aj44/zip52UlJp4+LHrczRyKpnKYiGEqAZ+DpSj1JB+Skr5Pa3uZJbQGwS/91Uvf/+3zTz1j238xTfqksVPFxadTrDrtkJ23VbIiVN+DhwYIhaTxGISr9dMeblieSUSkkAgjhACKSXtbWNcuhhg+QoHJWUWRJ50akgp6ekeo6kpQNP5AK0tQRIJ8HotPPhwOVs2uSYNSuQTsZjkJz9sZ2Q4zl9+s+qKRbwUyNAvjQH/h5TyIyGEEzgshNgN/B5a3cnsUFZh5vf/QyX//L02fvKDdp7806qc3mCrVztZvdpJNJrgo4+G+bd/68LlMuByGRkeiXHb7VdnPDY3B3lzTx9v7unD7tDTuMxB43I7jY0OrAvcMR4ei3PpUoCLTQEuNvkZGVECdcvKzOzYWci6dS68lYoY52sOwXhc8i8/7uTC2SBPfMV73U0tmo3pEijMlWRa+/EU96NCiDMoqeznXHcy70VMzRtSDZVQeU9tu4TKlus2OvjCH1Twi5928ePvt/PVP0tRyFQ+IlSmFKn99njyPatesOsmN7tucnPytJ/hkRgb7ytOVrJW9n3PxwvZuc3F+aYAZ84GOH9+lONHldxhHo+R0lITZaUmSkpMlJaaKS01zbsSdjicwOcL4/NF6O2N0JN8PTAQRUqwmHUsX25j1Uo7K1fYcU/KGqK0W60FepXQjNnCIeazXTwueeZHnXx0cISHHytl+00Fqje12vepXtt5XEQkxa8vnlI/46npyrYBJAsRbQauqTsphJjVT897Ecs3brzZjZTwr/+ri6f+cQ5CtgCsW6N0LMdRXLWJ/WFOp4GPbSlg8+YCEglJW/sYFy4E6ekJ0+OLcPFikFhMTvi8HqfTgNmsu7JYJrxGCCLhBGPhOOFwgnA4wdhYgkg4gT8QZ3j46lQovR5Kik1Uei1s3uRiWaON2lorer3Iu+F6NeJxydNPXRWwO+9dbPk95o8gZRHtk1JunXV/QjiA54C/kFKOpNOHq4lYGoxXrZkoZIY8833ULgadTlBbY6W25qoblEhI+gei+HwRxYLqjRAIKAI1MhK7IlThcIJoVBE7g0FMEjmzWYfTaaC0zExpiYnSUsWyKy4yol/o1AYZRhOwq2TqTAohjCgC9gsp5W+Sq2eq3zEjmoilyVQhy3Uf2XzR6QRFRSaKikysniX7QzwukVIRsZT2vbj16xoXcikLGBlKtSOUp+xPgDNSyu9OeGvB6k5qoAjZ53+/gnOnA/z4++1EwrmZLL7Q6PUiZQFb7MRimoBNZNydnG1JgZ3AE8AdQoijyeU+8qDu5JJjokX2j3/fwlf+rGra1NYai49gIM5PftBO09mgJmATyMTjS0q5T2VXc6o7qVliGWDHLW6e/NMqurvCfOdvLnPqWPariWtklwvngnznW5e5dCHEF5/0agI2AZHCspBct5ZYuhkC1P39mYfSN212Uvp/1fG//rmTf/5eOztvc/Opx8uwWHSq36dWz049HGDmtqglgs03h1c1y0O6IQgq+5wtjCIaTfDib3p587UBikuM/Mdv1FLfaFUNo1Db52LNVKFGvvVxXrcilgu8VRb+83+p46Xf9vHGq/2cOx3gy1/xUrfMNvvGGjmnvXWMp3/USWd7mJtvd/PIY8pDSOMqubC0ZkMTsQxjNOp4+LFS1m508MyPO/nu/9PCnfcVcd9DJUumM3yxkUhI9rzaz0u/7cXu0PPHf1HNuo1LYzJ3Oix0Dv3Z0EQsSyxfaeObf1vPc8/2sPulfs6c8PPEk168VYszs8T1Sp8vwjM/6eRSU4hNW5187kvlOJzabTEjqY8+Lhja2coiVqueL/6Bl3WbnDz7sy6+863L3PrxQu5+oHhJ5J3KZ8JjCd54pZ89r/Wj1wueeNLL9ptciyrrRy4Q5N9ooCZiC8CGzU7qG628+JyPvbsHeP/dIXbdVcjtnyhUzRbbcilEYbERp0s7TZkiHE7w7puD7HmlH78/zpZtLh76TCmFRca8szDylXw7TvO6O4QQzcAoyiBbTEq5NZ18QEsBp8vA53/fy667Cnn5+T5efaGPt3cPcNtdhdx2V+E12SSGh6Ls3ztIW/MYeoPggUdLWLlW66dJl0g4wb63Bnnz1X78o3FWr7Nz38MlSy4DRSbIMw3LiCV2u5RyYqXZbzDHfEDpoloMROVN1QwXKm+qPYHUsl9M3K66xswf/WllsrCEImZ2m4477p4ch+Ry6Hnsi+WYzTpOHhtl/1uDrF5jZ2wsoWzj0LP5BhclZXNL05xQ+X3Z6q5N9zxlYp/BYJz33x1i98v9jI4o4nX/w8U0LBsXr6u/Ws1NSjekIxthIrlmKYRYzDkf0FKkutbCH/1ZFa3NYUrLrxUik0lHNJlVwmjSoTcILjWFOHxwGKNRx0B/lOf/t48nv1610E3PexIJSdPZIO+/O8TRw6NEo5JVa+zc93Axy1Zoltd8yFRm10wyXxGTwOtCiQT852S+oJTyAQkhvgp8FeDjt/y3eTZj8VJdO/1oZSwmr2R+eG/vICvX2Hl/3xAul4GHPlOKBH79TBenT/hZs0Ryu8/G4ECUA/uG+OC9Yfp7o1itOnbc4mbHLQUTqnHnV3jAYiTPNGzeIrZTStmZFKrdQoizqW6YFLynAO689dvalTUFg0HQ3RXm2ae7Ka8wseMWN2++NsATT3qvfKbPFyWWTIszNX/YUmF4MMrpEwGOHBrh7KkAUsLKNTY++UgJGz/mXNSZRfITmZHMrplkXiImpexM/u8TQvwW2EYa+YA0JhMMxjl6aJTDB0fYss3FLbd7GB6KEgzEJ3VEt7WM0Zh0j5aKgMXjkuaLIU6f8HPmhJ/21jAAniIDdz9QzPabCygtzZ9ybtcj+XalpS1iQgg7oEvmx7YDnwD+ljTyAWlcJRpN8PMfddJ8KcTjX65gw2almlHL5bFJAnbxfBCbXY/dYSCRkITDCfR6cd1ZHrGYpKNtjJZLIZrOBjl/JkAomECng4blNh78dAlrNzioqFxclZ0WK3NItbNgzMcSKwN+m7xwDMC/SilfFUJ8CPxaCPGHQCvwmfk3c+mg1wtu2FGAq8DAnlf6eeX5Xu55sAS7Q0+B20A8rvSVHT4wwqatzivb7XtzkJf+vZfKKgu1DRZq663UNlgpLTddqRie70gp6fVFabkUUpbLIdpbw8STAxwFHgMbtzhZvcHBqjX2vKj+vRTJt6OetohJKS8BG6dZ388c8wGpoXr7qb2p4rarZR1Qy/Kg1hOgOuSvts+p/Qs62H6jk+03OpESOjvCSAmVVWbee2uAf/p/W7A79FRUmtl5awEGnQQdrFpjIxIuovliiEPvj7DvrSEAzGZBSamJkjITpeWmq6/LTDhd+gW3XqSUBPxxfD0RepOLrydCry+KryfCWEg5AyaToKbewu13eqhrsFDXaMVTaLjSXuWYTn9G1H5Sulkl0t1nNkIzcv1IytQlI4T4KfAA4JNSrkuu0+pOXm94K68Wjf2DP/Jy/myQXl+UHbcWTLKw6hut1Dcq7mYiIenpjtByaYy2ljF8PRE62sY4dmSUxITcP2aLDpdLj8NpwOHU43Qqr53JdWaLDoNByeJqMIpJrwWCWCyh1L2MKrUvo8n/I+EE/tE4o6Mx/KNx/KMxRkfj+EfijIzErggVKDdEUbGR0jIT9Y0uqmoUS7LCa170efmvTzKax+JnwD+hFNEdZ85xppqILSJ0OsGqNXZWrVG3GHU6QYXXTIXXzI03F1xZH49L+npjigXki9DnizAyoojMQH+Ulssh/P74JKGbD0KQFEdFJKtqzDhddopLjFesweJSo5bdYxGhSFhmzpeU8p1kubaJXH91JzUyh14vKE2Kx0xIKQkFE4yOxoiEZdLSSlypND5ubSG5xjobf20y63A49dhs+ln74/JtuF5jdoTI6uCRVndSY34IIbDZ9djs+dZ9q5E/pGSJpVw8d75oIqahoTEHBCK1ZDwpFc+dhjnHmV5fQUUaGhpZRwjdrMs8GI8zhRTjTK9bS0x1GFilGybdkKp0QzPUUNtOzdlT/b5F1AeVbvdxuiEAC51tI9+CRlMnMw0XQjyL0olfLIRoB/4aJUh+TnGm162IaWhoZB6R/JcJpJSfm+GtOcWZaiKmoaExJzIlYplCEzENDY05IUR+jVxrIqahoTEH8q/ypCZiGhoac0JzJzU0NBY5+RWZtahFLO3nQRYyBGTl2ZRmsZOlwELbAtkIh8gveyZ1NEtMQ0Nj0SKEyLvkk5qIaWhozAmRZ2kRNRHT0NCYI5olpqGhsWjR3EkNDY1FT36JWN6MlR7/cDTXTdDQuC45eXiUhFqGgjki0M26LCR5YYlF4l389H8omUe373JPei8bmT8XOhxCQyNXnDw8yvf/toXqBgu8m6m95tfFnheW2FDoDVastfOjf2jjwN6hXDdHQ+O6YFzAKqrN/Kf/uz4j+xQIdEI367KQ5IWIQYw//1adJmQaGhliqoA5XJl0unQpLAtHXriToJQP+/Nv1fG9v27mR//QBnCNa5lv/PR/tHP84AhOt4H/+j9XAPD8v/TwzmsDOAuUQ/upL5ex4QZXLpupscTIroBpEfuqLDYh23mnh49/sogf//e2SevveriYex4tyVGrNJYy2RawTGaxEELcA3wPJVHxj6WUf5fOfvLEnbzKuJAtBtdy5Xo7dmd+RS9rLF2yL2AK41OP1JYU9qEHfgDcC6wBPieEWJNOe/LKEhtnokX21HfaePr7HbluEgDf+sFySspnrtk4zpsv9vP+niFql1v57JMVk4QukZB8/TOns9lMjSVKJJygqs6SVQGDjE072gZckFJeAhBC/BKlcO6cbw4hZe7zIQgheoGWHDejGOjLcRsmorVndvKtTfnenlop5bz6OYQQryb3OxsWYGzC35PqTgohPg3cI6V8Mvn3E8B2KeXX59qmvLDE5ntgM4EQ4lCadfKygtae2cm3Ni2F9kgp78nQrqbzOdOyqPKuT0xDQ2NJ0A5UT/i7CuhMZ0eaiGloaOSCD4HlQoh6IYQJeBylcO6cyQt3Mk94avaPLChae2Yn39qktSdFpJQxIcTXgddQQix+KqU8lc6+8qJjX0NDQyNdNHdSQ0NjUaOJmIaGxqJmSYqYEKJZCHFCCHFUCHEoua5QCLFbCNGU/N+T5Tb8VAjhE0KcnLBuxjYIIb4phLgghDgnhLh7gdrzN0KIjuRxOiqEuG8B21MthHhLCHFGCHFKCPHnyfU5OUYq7cnJMRJCWIQQB4UQx5Lt+VZyfc6uoZwhpVxyC9AMFE9Z9x3gG8nX3wD+PsttuBXYApycrQ0o0zKOAWagHrgI6BegPX8D/KdpPrsQ7akAtiRfO4Hzye/NyTFSaU9OjhFKnJUj+doIHABuzOU1lKtlSVpiM/AQ8HTy9dPAw9n8MinlO8BAim14CPillDIspbwMXECZtpHt9szEQrSnS0r5UfL1KHAGqCRHx0ilPTOR7fZIKaU/+acxuUhyeA3liqUqYhJ4XQhxWAjx1eS6MillFygXLFCag3bN1IZKYGKqjHbUb6BM8nUhxPGkuznumixoe4QQdcBmFGsj58doSnsgR8dICKEXQhwFfMBuKWVeHJ+FZqmK2E4p5RaUGfRfE0LcmusGzULGpmjMkf8JNAKbgC7gvy90e4QQDuA54C+klCNqH12INk3TnpwdIyllXEq5CSXafZsQYp3Kx3N1DWWdJSliUsrO5P8+4LcoZnWPEKICIPm/LwdNm6kNGZuiMReklD3JGyUB/Iir7seCtEcIYUQRjF9IKX+TXJ2zYzRde3J9jJJtGAL2AveQZ9fQQrDkREwIYRdCOMdfA58ATqJMefhy8mNfBp7PQfNmasMLwONCCLMQoh5YDhzMdmPGb4Ykj6AcpwVpj1CSUv0EOCOl/O6Et3JyjGZqT66OkRCiRAjhTr62AncCZ8mza2hByPXIwkIvQAPKKM0x4BTwV8n1RcAeoCn5f2GW2/EsivsRRXlK/qFaG4C/QhlROgfcu0DteQY4ARxHuQkqFrA9N6O4O8eBo8nlvlwdI5X25OQYARuAI8nvPQn8l9mu42yfs1wt2rQjDQ2NRc2Scyc1NDSuLzQR09DQWNRoIqahobGo0URMQ0NjUaOJmIaGxqJGEzENDY1FjSZiGhoai5r/H7xie6fUGPIPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot model results\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy.utils.binaryfile as bf\n",
    "\n",
    "# create the headfile object\n",
    "fname = os.path.join(workspace, headfile)\n",
    "headobj = flopy.utils.binaryfile.HeadFile(fname)\n",
    "head = headobj.get_data(totim=27)\n",
    "times = headobj.get_times()\n",
    "\n",
    "# Setup contour parameters\n",
    "levels = np.arange(0, 90, 5)\n",
    "extent = (delr/2., Lx - delr/2., delc/2., Ly - delc/2.)\n",
    "\n",
    "# Make the plots\n",
    "plt.subplot(1, 1, 1, aspect='equal')\n",
    "plt.title('Head distribution (m)')\n",
    "plt.imshow(head[0, :, :], extent=extent, cmap='YlGnBu', vmin=0., vmax=90.)\n",
    "plt.colorbar()\n",
    "\n",
    "contours = plt.contour(np.flipud(head[0, :, :]), levels=levels, extent=extent, zorder=10)\n",
    "plt.clabel(contours, inline=1, fontsize=10, fmt='%d')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "BTN error. Required input is None, but no modflow model. If no modflow model is passed to Mt3dms, then values must be specified in the BTN constructor for: nlay, nrow, ncol, nper, laycon, delr, delc, htop, dz, perlen, nstp, and tsmult.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-2b82baaa0bf0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     24\u001b[0m     \u001b[0mthkmin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m     \u001b[0mmunit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'g'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m     \u001b[0micbund\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0micbund\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m )\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, MFStyleArr, DRYCell, Legacy99Stor, FTLPrint, NoWetDryPrint, OmitDryBud, AltWTSorb, nlay, nrow, ncol, nper, ncomp, mcomp, tunit, lunit, munit, laycon, delr, delc, htop, dz, prsity, icbund, sconc, cinact, thkmin, ifmtcn, ifmtnp, ifmtrf, ifmtdp, savucn, nprs, timprs, obs, nprobs, chkmas, nprmas, perlen, nstp, tsmult, ssflag, dt0, mxstrn, ttsmult, ttsmax, species_names, extension, unitnumber, filenames, **kwargs)\u001b[0m\n\u001b[1;32m    218\u001b[0m         \u001b[0;31m# they are specified in the constructor.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    219\u001b[0m         self.setmodflowvars(nlay, nrow, ncol, nper, laycon, delr, delc, htop,\n\u001b[0;32m--> 220\u001b[0;31m                             dz, perlen, nstp, tsmult)\n\u001b[0m\u001b[1;32m    221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    222\u001b[0m         \u001b[0;31m# Make the rest of the assignments\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36msetmodflowvars\u001b[0;34m(self, nlay, nrow, ncol, nper, laycon, delr, delc, htop, dz, perlen, nstp, tsmult)\u001b[0m\n\u001b[1;32m    334\u001b[0m                 \u001b[0ms\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m'perlen, nstp, and tsmult.'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    335\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 336\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    337\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    338\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mnlay\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: BTN error. Required input is None, but no modflow model. If no modflow model is passed to Mt3dms, then values must be specified in the BTN constructor for: nlay, nrow, ncol, nper, laycon, delr, delc, htop, dz, perlen, nstp, and tsmult."
     ]
    }
   ],
   "source": [
    "## Trial 1\n",
    "# name\n",
    "name_mt3d = 'mt-hatari'\n",
    "\n",
    "# exe name\n",
    "exe_name = r'/Users/mattijnvanhoek/modflow6/bin/mt3dusgs'\n",
    "\n",
    "#MT3D-USGS\n",
    "mt_model = flopy.mt3d.Mt3dms(\n",
    "    modelname=name_mt3d,\n",
    "    version='mt3d-usgs',\n",
    "    exe_name=exe_name,\n",
    "    modflowmodel=sim\n",
    ")\n",
    "\n",
    "#BTN file\n",
    "icbund = idomain.copy()\n",
    "icbund[0, 15, 15] = -1 #constant concentration\n",
    "\n",
    "btn = flopy.mt3d.Mt3dBtn(\n",
    "    model=mt_model, \n",
    "    sconc=0.0, \n",
    "    prsity=0.3, \n",
    "    thkmin=0.01, \n",
    "    munit='g', \n",
    "    icbund=icbund\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'ModflowGwfdis' object has no attribute 'nper'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-4e713c67edf9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     24\u001b[0m     \u001b[0mthkmin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m     \u001b[0mmunit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'g'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m     \u001b[0micbund\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0micbund\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m )\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, MFStyleArr, DRYCell, Legacy99Stor, FTLPrint, NoWetDryPrint, OmitDryBud, AltWTSorb, nlay, nrow, ncol, nper, ncomp, mcomp, tunit, lunit, munit, laycon, delr, delc, htop, dz, prsity, icbund, sconc, cinact, thkmin, ifmtcn, ifmtnp, ifmtrf, ifmtdp, savucn, nprs, timprs, obs, nprobs, chkmas, nprmas, perlen, nstp, tsmult, ssflag, dt0, mxstrn, ttsmult, ttsmax, species_names, extension, unitnumber, filenames, **kwargs)\u001b[0m\n\u001b[1;32m    218\u001b[0m         \u001b[0;31m# they are specified in the constructor.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    219\u001b[0m         self.setmodflowvars(nlay, nrow, ncol, nper, laycon, delr, delc, htop,\n\u001b[0;32m--> 220\u001b[0;31m                             dz, perlen, nstp, tsmult)\n\u001b[0m\u001b[1;32m    221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    222\u001b[0m         \u001b[0;31m# Make the rest of the assignments\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36msetmodflowvars\u001b[0;34m(self, nlay, nrow, ncol, nper, laycon, delr, delc, htop, dz, perlen, nstp, tsmult)\u001b[0m\n\u001b[1;32m    354\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnper\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    355\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 356\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnper\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    357\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    358\u001b[0m         \u001b[0mnlay\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnlay\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'ModflowGwfdis' object has no attribute 'nper'"
     ]
    }
   ],
   "source": [
    "## Trial 2\n",
    "# name\n",
    "name_mt3d = 'mt-hatari'\n",
    "\n",
    "# exe name\n",
    "exe_name = r'/Users/mattijnvanhoek/modflow6/bin/mt3dusgs'\n",
    "\n",
    "#MT3D-USGS\n",
    "mt_model = flopy.mt3d.Mt3dms(\n",
    "    modelname=name_mt3d,\n",
    "    version='mt3d-usgs',\n",
    "    exe_name=exe_name,\n",
    "    modflowmodel=gwf  # <--- changed from sim into gwf\n",
    ")\n",
    "\n",
    "#BTN file\n",
    "icbund = idomain.copy()\n",
    "icbund[0, 15, 15] = -1 #constant concentration\n",
    "\n",
    "btn = flopy.mt3d.Mt3dBtn(\n",
    "    model=mt_model, \n",
    "    sconc=0.0, \n",
    "    prsity=0.3, \n",
    "    thkmin=0.01, \n",
    "    munit='g', \n",
    "    icbund=icbund\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{internal}\n",
       "(1)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Trial 3\n",
    "tdis.nper  # time discretization package contains a nper attr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'MFArray' object has no attribute 'get_value'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-20-9e2899bfa9bd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     24\u001b[0m     \u001b[0mmunit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'g'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m     \u001b[0micbund\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0micbund\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m     \u001b[0mnper\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtdis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnper\u001b[0m  \u001b[0;31m# <--- from the time discretization package\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m )\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, MFStyleArr, DRYCell, Legacy99Stor, FTLPrint, NoWetDryPrint, OmitDryBud, AltWTSorb, nlay, nrow, ncol, nper, ncomp, mcomp, tunit, lunit, munit, laycon, delr, delc, htop, dz, prsity, icbund, sconc, cinact, thkmin, ifmtcn, ifmtnp, ifmtrf, ifmtdp, savucn, nprs, timprs, obs, nprobs, chkmas, nprmas, perlen, nstp, tsmult, ssflag, dt0, mxstrn, ttsmult, ttsmax, species_names, extension, unitnumber, filenames, **kwargs)\u001b[0m\n\u001b[1;32m    218\u001b[0m         \u001b[0;31m# they are specified in the constructor.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    219\u001b[0m         self.setmodflowvars(nlay, nrow, ncol, nper, laycon, delr, delc, htop,\n\u001b[0;32m--> 220\u001b[0;31m                             dz, perlen, nstp, tsmult)\n\u001b[0m\u001b[1;32m    221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    222\u001b[0m         \u001b[0;31m# Make the rest of the assignments\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/flopy/mt3d/mtbtn.py\u001b[0m in \u001b[0;36msetmodflowvars\u001b[0;34m(self, nlay, nrow, ncol, nper, laycon, delr, delc, htop, dz, perlen, nstp, tsmult)\u001b[0m\n\u001b[1;32m    368\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    369\u001b[0m             self.delr = Util2d(self.parent, (ncol,), np.float32,\n\u001b[0;32m--> 370\u001b[0;31m                                \u001b[0mmf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    371\u001b[0m                                \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'delr'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    372\u001b[0m                                \u001b[0mlocat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munit_number\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'MFArray' object has no attribute 'get_value'"
     ]
    }
   ],
   "source": [
    "# name\n",
    "name_mt3d = 'mt-hatari'\n",
    "\n",
    "# exe name\n",
    "exe_name = r'/Users/mattijnvanhoek/modflow6/bin/mt3dusgs'\n",
    "\n",
    "#MT3D-USGS\n",
    "mt_model = flopy.mt3d.Mt3dms(\n",
    "    modelname=name_mt3d,\n",
    "    version='mt3d-usgs',\n",
    "    exe_name=exe_name,\n",
    "    modflowmodel=gwf  # <--- changed from sim into gwf\n",
    ")\n",
    "\n",
    "#BTN file\n",
    "icbund = idomain.copy()\n",
    "icbund[0, 15, 15] = -1 #constant concentration\n",
    "\n",
    "btn = flopy.mt3d.Mt3dBtn(\n",
    "    model=mt_model, \n",
    "    sconc=0.0, \n",
    "    prsity=0.3, \n",
    "    thkmin=0.01, \n",
    "    munit='g', \n",
    "    icbund=icbund,\n",
    "    nper=tdis.nper  # <--- from the time discretization package\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'MFArray' object has no attribute 'get_value'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-bc45a0b6e600>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgwf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m: 'MFArray' object has no attribute 'get_value'"
     ]
    }
   ],
   "source": [
    "gwf.dis.delr.get_value()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10.,\n",
       "       10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10.,\n",
       "       10., 10., 10., 10., 10.])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gwf.dis.delr.get_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # advaction package\n",
    "# adv = flopy.mt3d.Mt3dAdv(mt, mixelm=-1, percel=0.75)\n",
    "# # dispersion package\n",
    "# dsp = flopy.mt3d.Mt3dDsp(mt, al=10.0, trpt=0.1, trpv=0.1, dmcoef=1e-09)\n",
    "\n",
    "# # source/sink package\n",
    "# ssm_data = {}\n",
    "# itype = flopy.mt3d.Mt3dSsm.itype_dict()\n",
    "# ssm_data[0] = [(0, wrow, wcol, 10.0, itype['WEL'])]\n",
    "\n",
    "# ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)\n",
    "\n",
    "# # matrix solver package\n",
    "# gcg = flopy.mt3d.Mt3dGcg(mt, cclose=1e-6)\n",
    "\n",
    "# # write mt3dms input\n",
    "# mt.write_input()\n",
    "# # run mt3dms\n",
    "# mt.run_model()\n",
    "\n",
    "# # post-processing and plotting\n",
    "# # plot flow\n",
    "# import matplotlib.pyplot as plt\n",
    "# import flopy.utils.binaryfile as bf\n",
    "# fig = plt.figure(figsize=(10,10))\n",
    "# ax = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "\n",
    "# wpt = ((wcol+0.5)*delr, Lx - ((wrow + 0.5)*delc)) # origin at low upper..\n",
    "# hds = bf.HeadFile(modelname+'.hds')\n",
    "# times = hds.get_times() # simulation time, steady state\n",
    "# head = hds.get_data(totim=times[-1])\n",
    "\n",
    "# cbb = bf.CellBudgetFile(modelname+'.cbc') # read budget file\n",
    "# frf = cbb.get_data(text='FLOW RIGHT FACE', totim=times[-1])[0]\n",
    "# fff = cbb.get_data(text='FLOW FRONT FACE', totim=times[-1])[0]\n",
    "\n",
    "# # create flopy plot object, plot grid and contour\n",
    "# modelmap = flopy.plot.ModelMap(model=mf, layer=0)\n",
    "# lc = modelmap.plot_grid() # grid\n",
    "# cs = modelmap.contour_array(head, levels=np.linspace(head.min(), head.max(), 21)) # head contour\n",
    "# plt.clabel(cs, fontsize=20, fmt='%1.1f', zorder=1) # contour label\n",
    "# quiver = modelmap.plot_discharge(frf, fff, head=head) # quiver\n",
    "# plt.plot(wpt[0],wpt[1],'ro') # well location\n",
    "# plt.show()\n",
    "\n",
    "# # plot conc\n",
    "# ucnobj = bf.UcnFile('MT3D001.UCN')\n",
    "# #print(ucnobj.list_records()) # get values\n",
    "\n",
    "# times = ucnobj.get_times() # simulation time\n",
    "# times1 = times[round(len(times)/4.)] # 1/4 simulation time\n",
    "# times2 = times[round(len(times)/2.)] # 1/2 simulation time\n",
    "# times3 = times[-1] # the last simulation time\n",
    "\n",
    "# conc1 = ucnobj.get_data(totim=times1)\n",
    "# conc2 = ucnobj.get_data(totim=times2)\n",
    "# conc3 = ucnobj.get_data(totim=times3)\n",
    "\n",
    "\n",
    "# # conc 100 day\n",
    "# fig = plt.figure(figsize=(10,10))\n",
    "# ax = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "# modelmap = flopy.plot.ModelMap(model=mf, layer=0)\n",
    "# lc = modelmap.plot_grid() # grid\n",
    "# cs = modelmap.plot_array(conc1) # head contour\n",
    "# plt.colorbar(cs) # colrobar\n",
    "# plt.plot(wpt[0],wpt[1],'ro')\n",
    "# plt.title('C  %g day' % times1)\n",
    "# plt.show()\n",
    "\n",
    "# # conc in 500 days\n",
    "# fig = plt.figure(figsize=(10,10))\n",
    "# ax = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "# modelmap = flopy.plot.ModelMap(model=mf, layer=0)\n",
    "# lc = modelmap.plot_grid() # grid\n",
    "# cs = modelmap.plot_array(conc2) # head contour\n",
    "# plt.colorbar(cs) # colrobar\n",
    "# plt.plot(wpt[0],wpt[1],'ro')\n",
    "# plt.title('C  %g day' % times2)\n",
    "# plt.show()\n",
    "\n",
    "\n",
    "# # conc in 1000 days\n",
    "# fig = plt.figure(figsize=(10,10))\n",
    "# ax = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "# modelmap = flopy.plot.ModelMap(model=mf, layer=0)\n",
    "# lc = modelmap.plot_grid() # grid\n",
    "# cs = modelmap.plot_array(conc3) # head contour\n",
    "# plt.colorbar(cs) # colrobar\n",
    "# plt.plot(wpt[0],wpt[1],'ro')\n",
    "# plt.title('C  %g day' % times3)\n",
    "# plt.show()\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}