examples/Notebooks/flopy3_lake_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FloPy\n",
    "\n",
    "### Lake Example\n",
    "\n",
    "First set the path and import the required packages. The flopy path doesn't have to be set if you install flopy from a binary installer. If you want to run this notebook, you have to set the path to your own flopy path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:47.367638Z",
     "iopub.status.busy": "2022-07-30T20:16:47.367227Z",
     "iopub.status.idle": "2022-07-30T20:16:48.790628Z",
     "shell.execute_reply": "2022-07-30T20:16:48.788478Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.10.5 (v3.10.5:f377153967, Jun  6 2022, 12:36:10) [Clang 13.0.0 (clang-1300.0.29.30)]\n",
      "numpy version: 1.23.1\n",
      "matplotlib version: 3.5.2\n",
      "flopy version: 3.3.6\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "from tempfile import TemporaryDirectory\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# run installed version of flopy or add local path\n",
    "try:\n",
    "    import flopy\n",
    "except:\n",
    "    fpth = os.path.abspath(os.path.join(\"..\", \"..\"))\n",
    "    sys.path.append(fpth)\n",
    "    import flopy\n",
    "\n",
    "# temporary directory\n",
    "temp_dir = TemporaryDirectory()\n",
    "workspace = temp_dir.name\n",
    "\n",
    "print(sys.version)\n",
    "print(\"numpy version: {}\".format(np.__version__))\n",
    "print(\"matplotlib version: {}\".format(mpl.__version__))\n",
    "print(\"flopy version: {}\".format(flopy.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are creating a square model with a specified head equal to `h1` along all boundaries. The head at the cell in the center in the top layer is fixed to `h2`. First, set the name of the model and the parameters of the model: the number of layers `Nlay`, the number of rows and columns `N`, lengths of the sides of the model `L`, aquifer thickness `H`, hydraulic conductivity `k`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:48.847542Z",
     "iopub.status.busy": "2022-07-30T20:16:48.846627Z",
     "iopub.status.idle": "2022-07-30T20:16:48.852845Z",
     "shell.execute_reply": "2022-07-30T20:16:48.851936Z"
    }
   },
   "outputs": [],
   "source": [
    "name = \"lake_example\"\n",
    "h1 = 100\n",
    "h2 = 90\n",
    "Nlay = 10\n",
    "N = 101\n",
    "L = 400.0\n",
    "H = 50.0\n",
    "k = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a MODFLOW model and store it (in this case in the variable `ml`, but you can call it whatever you want). The modelname will be the name given to all MODFLOW files (input and output). The exe_name should be the full path to your MODFLOW executable. The version is either 'mf2k' for MODFLOW2000 or 'mf2005'for MODFLOW2005."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:48.859296Z",
     "iopub.status.busy": "2022-07-30T20:16:48.858790Z",
     "iopub.status.idle": "2022-07-30T20:16:48.867373Z",
     "shell.execute_reply": "2022-07-30T20:16:48.865591Z"
    }
   },
   "outputs": [],
   "source": [
    "ml = flopy.modflow.Modflow(\n",
    "    modelname=name, exe_name=\"mf2005\", version=\"mf2005\", model_ws=workspace\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the discretization of the model. All layers are given equal thickness. The `bot` array is build from the `Hlay` values to indicate top and bottom of each layer, and `delrow` and `delcol` are computed from model size `L` and number of cells `N`. Once these are all computed, the Discretization file is built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:48.874097Z",
     "iopub.status.busy": "2022-07-30T20:16:48.872779Z",
     "iopub.status.idle": "2022-07-30T20:16:48.930671Z",
     "shell.execute_reply": "2022-07-30T20:16:48.920291Z"
    }
   },
   "outputs": [],
   "source": [
    "bot = np.linspace(-H / Nlay, -H, Nlay)\n",
    "delrow = delcol = L / (N - 1)\n",
    "dis = flopy.modflow.ModflowDis(\n",
    "    ml,\n",
    "    nlay=Nlay,\n",
    "    nrow=N,\n",
    "    ncol=N,\n",
    "    delr=delrow,\n",
    "    delc=delcol,\n",
    "    top=0.0,\n",
    "    botm=bot,\n",
    "    laycbd=0,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we specify the boundary conditions and starting heads with the Basic package. The `ibound` array will be `1` in all cells in all layers, except for along the boundary and in the cell at the center in the top layer where it is set to `-1` to indicate fixed heads. The starting heads are used to define the heads in the fixed head cells (this is a steady simulation, so none of the other starting values matter). So we set the starting heads to `h1` everywhere, except for the head at the center of the model in the top layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:48.968777Z",
     "iopub.status.busy": "2022-07-30T20:16:48.962226Z",
     "iopub.status.idle": "2022-07-30T20:16:48.984597Z",
     "shell.execute_reply": "2022-07-30T20:16:48.983496Z"
    }
   },
   "outputs": [],
   "source": [
    "Nhalf = int((N - 1) / 2)\n",
    "ibound = np.ones((Nlay, N, N), dtype=int)\n",
    "ibound[:, 0, :] = -1\n",
    "ibound[:, -1, :] = -1\n",
    "ibound[:, :, 0] = -1\n",
    "ibound[:, :, -1] = -1\n",
    "ibound[0, Nhalf, Nhalf] = -1\n",
    "start = h1 * np.ones((N, N))\n",
    "start[Nhalf, Nhalf] = h2\n",
    "bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aquifer properties (really only the hydraulic conductivity) are defined with the LPF package. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:48.998785Z",
     "iopub.status.busy": "2022-07-30T20:16:48.998185Z",
     "iopub.status.idle": "2022-07-30T20:16:49.013455Z",
     "shell.execute_reply": "2022-07-30T20:16:49.011795Z"
    }
   },
   "outputs": [],
   "source": [
    "lpf = flopy.modflow.ModflowLpf(ml, hk=k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need to specify the solver we want to use (PCG with default values), and the output control (using the default values). Then we are ready to write all MODFLOW input files and run MODFLOW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:49.019504Z",
     "iopub.status.busy": "2022-07-30T20:16:49.018925Z",
     "iopub.status.idle": "2022-07-30T20:16:49.932211Z",
     "shell.execute_reply": "2022-07-30T20:16:49.929686Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FloPy is using the following executable to run the model: /Users/wes/dev/usgs-executables/mf2005\n",
      "\n",
      "                                  MODFLOW-2005     \n",
      "    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL\n",
      "                             Version 1.12.00 2/3/2017                        \n",
      "\n",
      " Using NAME file: lake_example.nam \n",
      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2022/07/30 16:16:49\n",
      "\n",
      " Solving:  Stress period:     1    Time step:     1    Ground-Water Flow Eqn.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2022/07/30 16:16:49\n",
      " Elapsed run time:  0.521 Seconds\n",
      "\n",
      "  Normal termination of simulation\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(True, [])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcg = flopy.modflow.ModflowPcg(ml)\n",
    "oc = flopy.modflow.ModflowOc(ml)\n",
    "ml.write_input()\n",
    "ml.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the model has terminated normally, we can read the heads file. First, a link to the heads file is created with `HeadFile`. The link can then be accessed with the `get_data` function, by specifying, in this case, the step number and period number for which we want to retrieve data. A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions are used to make contours of the layers or a cross-section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:49.941743Z",
     "iopub.status.busy": "2022-07-30T20:16:49.940708Z",
     "iopub.status.idle": "2022-07-30T20:16:50.451703Z",
     "shell.execute_reply": "2022-07-30T20:16:50.425858Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABLrklEQVR4nO2dd5xcVfn/3+fWaTuzPQkpBEgAASF0FOlSv0pTiqIgXxBUsIANUJEfqIgNRJAuTaoUKQJKR/2q9F4TWvr26TO3nd8fc7MukN1MZndndpfzfr02O/fec+95zmT2M6c853mElBKFQqFYE1qjDVAoFJMDJRYKhaIqlFgoFIqqUGKhUCiqQomFQqGoCiUWCoWiKqoWCyGELoR4Rghxd3i8nhDiP0KIhUKIm4QQVnjeDo8XhtfnjpPtCoWijqxNz+KbwCtDjs8BzpVSzgP6gWPC88cA/eH5c8NyCoViklOVWAghZgH/A1weHgtgd+CWsMjVwIHh6wPCY8Lre4TlFQrFJMaostx5wPeApvC4DRiQUnrh8RJgZvh6JrAYQErpCSHSYfmeoQ8UQhwHHAcQj8e33njjjWtsgkKhqIannnqqR0rZUev9axQLIcSngC4p5VNCiF1rrej9SCkvBS4F2GabbeSTTz45Vo9WKBSrQQjxzmjur6ZnsSOwvxBiPyACJIHfAs1CCCPsXcwClobllwKzgSVCCANIAb2jMVKhUDSeNc5ZSClPlVLOklLOBQ4HHpJSHgE8DHw2LHYUcEf4+s7wmPD6Q1LtVlMoJj2j8bP4PnCyEGIhlTmJK8LzVwBt4fmTgVNGZ6JCoZgIVDvBCYCU8hHgkfD1m8B2qylTAg4ZA9sUCsUEQnlwKhSKqlBioVAoqkKJhUKhqAolFgqFoiqUWCgUiqpQYqFQKKpCiYVCoagKJRYKhaIqlFgoFIqqUGKhUCiqQomFQqGoCiUWCoWiKpRYKBSKqlBioVAoqkKJhUKhqAolFgqFoiqUWCgUiqpQYqFQKKpijWIhhIgIIR4XQjwnhHhJCPH/wvNXCSHeEkI8G/4sCM8LIcT5YfrC54UQW41zGxQKRR2oJgZnGdhdSpkTQpjAP4QQ94bXviulvOV95fcF5oc/2wMXhb8VCsUkpppUAFJKmQsPzfBnpND+BwDXhPf9m0p+kRmjN1WhUDSSanOd6kKIZ4Eu4H4p5X/CSz8NhxrnCiHs8Nxg+sKQoakNFQrFJKUqsZBS+lLKBVQyj20nhNgMOBXYGNgWaKWSR6RqhBDHCSGeFEI82d3dvXZWKxSKurNWqyFSygEqmcj2kVIuD4caZeBK/ptDZFX6wlUMTW049FmXSim3kVJu09FRc65WhUJRJ6pZDekQQjSHr6PAnsCrq+YhhBACOBB4MbzlTuDIcFVkByAtpVw+DrYrFIo6Us1qyAzgaiGETkVcbpZS3i2EeEgI0QEI4FngK2H5e4D9gIVAATh6zK1WKBR1Z41iIaV8HthyNed3H6a8BE4YvWkKhWIioTw4FQpFVSixUCgUVaHEQqFQVIUSC4VCURXVrIYoPgTIoEDgL0H6KwiClUi/Bxn0IWUaGWRBFpGyBHggfRCCyneNiRBR0GII0YTQkgitDaG1oenTEfo6aPo6/NfBVzFZUWLxIUMGWXz3ZXz3JQLvdQLvDQLvbWSwOi9aG6GlEFoSRAQhIoAJYtX2IAmyRCDT4OWRMosM0oD7vucIhD4DTV8P3ZyPZmyEbm6KZm4cPlMxGVBiMYWRUiL9d/DK/8Z3nsB3nyHwFg5eF6IZzZyHEdkdTV8XTZ8d/lFPQ2gdCC1WU53ILEHQU+ml+MuQ/mIC7218702cwp9A5sPSOpr5EXRzSwxrO3R7ezR9+hi1XjHWKLGYYsigiFf+O175IbzyY0h/CQBCa0E3t8KMHohmboZuborQOqk44I4dQggQSXQtCcb6H7RPSqS/GN99Cd99Ad95Brd4G27hWgA0Yz6GvStGZHd0azsqUREUEwElFlMAGRTwyg/gFu/CKz8KsgQigWF/HCPxFXRrBzRj/pgLQy0IIRDGHDRjDmZ0XwCk9Ajcl/Ccf+GVHsPJX4OTvwxEEjOyJ0b00xj2JxDCarD1H26UWExSpJT4zuO4hZtwS38BWURo0zBjh2NG9kS3tp80f1xCGOjWFujWFtiJr4Ti9xhe6W+4pb/hFm9FaC2Y0QMxY4ehm5s02uQPJUosJhkyyOEW/oRTuLYy/yASlT+i6MHo1rYIMflXw4UWw4zugxndh4h08MqP4RZuw8lfh5O/Et3cGjP+RczopyaNIE4FlFhMEgJ/JU7+Dzj560Bm0MwFRJp/jRn5FEKLNtq8cUMICzPySczIJwmCftzCrbj5P1Ia+BblzM+xEv+LFTsCoTU12tQpj6js+2os22yzjXzyyScbbcaEJPC7cHIX4eT/CLgYkX2xE8ejWwsabVrDkDLAKz+Kk7sE3/k/ECnsxHFY8aMRWqLR5k1YhBBPSSm3qfV+1bOYoMggTzl3CU7+UpBlzOhnsJtORDPmNtq0hiOEhhnZDTOyG77zHOXseZSzv8TJX4nddBJm7HC1ijIOTP4B7hRDSolTuI1c1y44ufMw7N2Idz5EtOVXSihWg25tQaztSmLtd6AZG1BK/4B899545X802rQphxKLCYTvvUWh93BKA99C6DOItf+ZWOtF6MZ6jTZtwmNYWxJru5lo6+VI6VDo/TyF/hMJ/J5GmzZlUGIxAZDSp5y7jHzXXvjuS0RSPyPefgeGpfIzrQ1CCMzIXiQ6H8BKfAuveC/57k/iFu9qtGlTAiUWDSbwV1DoPYJy5iwMe6fKBz3+hSmxBNoohIgQSZ5MvOMehD6bYv8JFPu/hQxya75ZMSzqE9lA3NKDYW/iWSLNvyLaeoXaGzGG6OZGxNtvx246Cbf4Z3Ld++I7LzTarEnLaHKdrieE+E+Y0/QmEXrHCCHs8HhheH3uOLdh0iFlQDl7HsW+oxH6TOId92DFDp0Q7thTDSEM7KaTiLXdDLjkew7GKdzaaLMmJdX0LFblOt0CWADsE4b4Pwc4V0o5D+gHjgnLHwP0h+fPDcspQmRQpNj/VcrZ32BGP0u843b01Wy4Uowthr0d8fZ70K2tKQ2cRCn9U6QMGm3WpGI0uU53B1YlRb6aSu4QqOQ6vTp8fQuwh1BfmQAEfi/53sPwSvdhJ08n0vxrFc+hjmh6K7G2P2LGjsLJX0Kx/6thQB9FNdSU6xRYBAxIKb2wyNB8poO5TsPraaBtNc/8UKUvDPxlFHo+Q+C+QrTlUuzEsWrY0QCEMIikzsROno5XupdC71Fq4rNKasp1SiXH6aj4MKUvDLzF5Hs+SxB0E2u7DjO6d6NN+lAjhMBOHEu0+Xx853EKvUcgg0yjzZrw1Jrr9GNAsxBilbv40Hymg7lOw+spoHcsjJ2MBN5S8r2HIYMs8bYbMOzt1nyToi6YsQOJtlyM775IofeLqoexBmrNdfoKFdH4bFjsKOCO8PWd4THh9YfkRNit1gACv5tC7+eRQYZ423Xo1uaNNknxPszo3kRbfo/vPk+h7xg1hzEC1fQsZgAPCyGeB54A7pdS3g18HzhZCLGQypzEFWH5K4C28PzJwCljb/bERwZ5Cn1fIgiWE2u7SgnFBMaM7k20+Tf4zr8o9n9TrZIMw2hynb5JZf7i/edLwCFjYt0kRUqfYv/XCdyXiLZejmHVvCtYUSfM2EEEQQ/lzFmUM2cTSf2g0SZNONQW9XGgnP0VXvkBIqmfYEY+2WhzFFVixY8l8N7FyV+CZm6MFftMo02aUCh37zHGLd6Dk7sQM3YEVvzIRpujWAuEEERSp6NbH6c0cIpyDX8fSizGkMB7h+LAd9HNrYik/l+jzRlTpPTw/V58byme9zae9w6+v4wg6J9SY3whTKItFyL0Ngr9X1MrJENQw5AxQkqXYv83AEG05YJJGUhWygDPW4TrPo/nvornvYHnLcH3l1FZNR8Ogaa1ousz0fU5GOZ8THNjTHMBuj5r0jmfaXob0ebzKfQeSin9Q6It5zXapAmBEosxopz9Hb77DNGWC9GMWY02p2o8dyGl8sOUS4/hOE8g5SrnJBPDWB/dmI1lb4uutVdSGYrYYMg6KR2kLBAE/QR+D76/BNd9mVLpHqDS29C0Dix7B2x7ZyKR3dD1dRrT0LXEsLfDSnyjEq0ssjtmdP9Gm9RwlFiMAb77Ek7uAszoQZjRTzfanDXieUsoFm6mWLwLz3sNAF1fn2h0fyxrK0xrSwxj/ZrjWEpZxHXfwHWexnGeolz+J6XiXaQB09qaaPQAotGD0PUP7AKYUNhN38ArP0wpfTq6tSPaBLd3vFHRvUeJlD75nk8j/ZXEOx9A01oabdJqkTKgXHqQfP4KyuW/A2BZ2xONfgo7shfGOPaGpJR43uuUSn+tCJT7EmASie5DIn4slr3tuNU9Wnz3NfLd+2FE9yPW8rtGmzMqVHTvBuMWriVwXyTacsGEFAopfYrF28hlz8fzFqHpM2hqOplo7LBxFYihCCEwzY0wzY1oavoGrvsqhfwNFAo3UyrehWluTVPyJGx7twk3v6GbG2ElTsDJnYcXOwzD/kSjTWoYqmcxCgK/j1zXzujm5sTarptQH3QpJaXSfWQzZ+N5CzHNzUgkvkYk+j8TJkx+EOQpFm4ml7sY31+MZW1LMvkjLHtiObFJWSLX9UmEsIh3/I3/bomaXIy2Z6GWTkeBkzsPZJ5I6owJJRSeu5C+3s/T33cMIGhpvZz2jr8SjR04YYQCQNPixBNH0znt76Saf47nLaanZ3/6+76J73c12rxBKjE9TyfwFuIWbmi0OQ1DiUWNBN7bOPk/YsYORzc3bLQ5QGXIkc1eQFfXnjjOMyRTZ9LR+SDR6H4TSszejxAW8fiRdE77O4nEiRSLd9C1clcKhduYCD1fACOyJ7q1HeXsuR9a3wslFjVSzp4LwsRuOqnRpgDge0vp6TmIbOZnRCJ7hH94x06qLrOmxUmmTqOz8wEMcwMG+k+kv/8rBBMg1oQQAjv5A2TQg5O/stHmNAQlFjXge4twi3dgxY5E0zsbbQ6l0iN0d++F575Kc8sFtLRehq5P3oBChjmP9vY/05Q8jVLxHrq79sV1X2q0WRjWlhj27ji5y5BBttHm1B0lFjXgZC8CYWEljm+0KeRyl9PXewSaNoOOjnuJxQ6e0EOOahFCp6npRNrab0XKEj3d+1Mq/q3RZmE3nYSUAziF6xttSt1RYrGWBP4K3OLtmNHD0PT2htkhZUB64Mdk0qcTiexNe8edGOYGDbNnvLDt7ejovBfD2JC+vv8ln2vsEEC3tkC3Po6TuwIpnYbaUm+UWKwlTv4awMdOHLPGsuOFlAEDA98mn7+MePxYWlovQ9Ni41SXy0DuepZ0f5Gu/jNw3LcACGQZx12EHwyMS71D0fVO2tpvxY7sQTr9A7LZC8a9zpGwEscjgxV4xb801I56o8RiLZCyjFu4ESOyR8Mymq8SimLhJhJNJ5NqPhMh9HGrL1e8j2zhdtqS30DTEqTzNxDIMgPZK1jWewKLuw5hZf8PQtuGX7nwRznG17QYra1XEI0eSDbzs4YKhmHvgqavF35xfHhQYrEWeKW/IoMerFhj4lRIKUmnfzgoFMnkd8a1vkCWKbuvEY/sTtTelqbYAYBgIPsHcqUHmTv9PuZMuxtdayNb+MuwcyWF0r94c/n29KR/MdiOWhDCoLnl/EHByOeuqrFlo0MIDTP+RXz3KXz3lYbY0AiqCdg7WwjxsBDi5TB94TfD82cIIZYKIZ4Nf/Ybcs+pYfrC14QQUybuvVO4CaHPQrd3akj9+dwlFPJXEU8cP+5CAaAJmyDI4AeV4Oyevxw/6CNXehBTnz1YRgiLQvlfQMXXo/K7suu05DxPtngnMXt7KvGeoZKjqjZWCYYd2ZN0+ocNm/Q0Y58BLNzCTQ2pvxFU07PwgG9LKTcBdgBOEEJsEl47V0q5IPy5ByC8djiwKbAP8Hsxnv3kOhH4y/DL/8CMfrYhGc6LxfvIZM4kEv00yeSPxrUuP8iSK/6VTP52mmL743iLWN77DTL5P+EHPUStrdC0KMXyk5SdV8nkb0HTmgbvlzJACA3XW062eC+GPouotT1C2GNinxAGLS0XYZofpb//q7jOi2Py3LVB01owInviFm9HSrfu9TeCatIXLpdSPh2+zlJJAzBzhFsOAG6UUpallG8BC1lNYN/Jhlu8G5BYsYPrXrfnvc1A/zcxzS1oaTlvXMQqkCUyhTtY0v0FFi7dlKU9R9Ofu5yItTXtqVNJxg4mYm2Joc8kGT+EqL0DK/tPIVOsZIAw9VWb0uSgfen8DdjmRrQk/hfXewdDGzvfD02L0dp2NUKk6Os7viGOW2bsYGTQj1f+R93rbgRr9akLM6JvCfwnPHWiEOJ5IcQfhBCrtlwOpi8MGZracOizJlX6Qrd4F5r50bpPbErp0N93PKDT0nrZkK782OD7/fRmfsuby7Zjee9XKbuv0NL0ZWZ33MKcztsBD9vckHh0V8ruS1jGBpj6TJKxA5g7/QHakydjGrMxjTkVe6lktMwV76foPImUZRxvIV7Qg2GsCnwzNi7cut5JS+tF+P67DAx8t+6u4Ya9M4gkXvHOutbbKKr2BRZCJIBbgW9JKTNCiIuAs6j8z58F/Br432qfJ6W8FLgUKrtO18boehP4ywjc57Cb6p8CJZs9F9d9gdbWK8d0S7mUDv25q+hN/4ZAZohHdqel6Vhi9k7vWV1Z1nsCnr8MTTSha60k4wfheu/iBd3Y5nx6M78jYm1FxNoCAC1M9KxrzVjGXNL5G/H9Xlx/Cb7fxbTWX2MZY+cPYtvb05T8HtnM2RQjexOrY89PCBszsidu6QEi0p1Qm/TGg6rEQlTehVuB66SUtwFIKVcOuX4ZcHd4OJi+MGRoasNJiVe6HwAjsldd63Wc58llLyAaO5TIGOZHLZafZkXft3C8hcQiu9KR+iERa5PVlu1s/jGutwTPX0nM3gFda6boP01/7nI8fyXxyC60Nh2PpsUIggI9mV/SnvwOUXtbomFQG9dbwsr+05jW8nO6Bn6MH/QwvfVcrDHqpSUSX6NU+hvpgR9i259Ar6MLvhHZG7d4K77zOIa9Y93qbQTVrIYIKlnGXpFS/mbI+RlDih0ErJpluhM4XAhhCyHWA+YDj4+dyfXHKz2Epq+Hbs6rW50VD81T0bQ2UmMUKVxKj570r3i36wACWWRm+9XMar9uWKEAMPROovZWNMX2RdcrI81EdHdmd1zPetMfpLP59EGHMCFM4vYuaFocKSVSepUgwH4PQpjoWgfxyK6UnVd4e8WepHM3jMnQQQid5ubfIGWBTOZno37e2mDYOwMWXunhutbbCKqZs9gR+CKw+/uWSX8hhHghTGu4G3ASgJTyJeBm4GXgPuAEuWo9bRIiZQnP+RdGZJe61lss3IzrPkMy9UM0LTXq5/lBmiU9R9Kb+Q3J2EHMnf4gieieY7qPRAiTeHTX8LVACAMhNKL2Ama2X4GmmTQnjmDu9AeJWgtY0f9tVvZ/f0xWE0xzPonEcRQLN+OU6xdISWgxdGtbvPIjdauzUVSTvvAfwOo+UfeMcM9PgZ+Owq4Jg+88BbKEbu9ctzqlLJLJ/ALT3JpodPRjcM9fweLuz+O4i5jW8kuaE0eMgZW1YxozmdVxEz3pc+jL/g7XX8LMtstH7bKeaPomhcKfyGTOoq39z3XbUGfYO1POnk3gr0TTp9WlzkYweYIdNAiv/G9Aw7Dqt/qbz11DEKygpfXCUS+Tut5yFncfguevZFbHH4lH1t6hzA+KlLxFFJ3XcINe/CCDQEPXkph6J1FzIyLmemhrkStFCI2O5lMxjXVZ2f89lvR8gVntfxyVYGhanKbkSaQHTqFcfoRIZLean7U2GPbHKWfBd/6DNoVTBiixWAO+8ziauQlCS9alPilL5HIXYtk7YdsfG9WzfL+fJd2H4/tdzO64fnDCsRrK3hJ687eTLj1KrvwMMHQkqbEqL8gqhLBpsrcnFdmV9vhBGHp1wYubE59H02Is7z2Rpb3HMqv96lGtKsRih5PLXkg285u6iYVmbgoijlf+z5TOL6LEYgSk9PHd5zBjh9atzmLhdoKgh6amr4/qOVK6LOs9Dtd7h1lrIRR550WWpy+kv/hXICBmbc6M5FeIWZsSNTfE0qejiTgQ4MscjreEoruQvPMs6eJjLB44kyXpc2iPHcyM1NewjdlrqpJk7EBkUAznMH7AtJZzah5CCGGRSBxPOv1DnPKTdQn+K4SBbm6B7z4z7nU1EiUWIxB4b4AsoJsL6lKflJJc7nIMcxMsa3TLcN3psymU/8n01t8Si3x8jeVdv5slA7+kJ/8ndC3J9ORxdCa+iG0M56yrY4gUhpUiZm1KW/wAaIGC8xpd2avoKdxGT/5Wpie/zIzkCehrGF6kEp/D8d6iL3sBUWtLUonP1dDqCtHYYWQyvySXv4zWOkUK162tcHIXI2UJEfqaTDXUrtMR8N1KFm3d/Ghd6nPdZ/C8V4jHvzSqybl86R/0Zy8mFT+SVPyQNZbPlP7Ji8v3pTd/O9OajmXzdf7O7OZTRhCK4YlZGzG37Ww2n/EIrbH/YXnmQl5e8WmKzutrvLc9dQoxe0dWDvwIx3t7retehabFicUOoVS8D9/vrfk5a0PlM+IRuK/Vpb5GoMRiBAL3ZRARNGP9utRXyN+EEFGi0QNqfkYQFFjRdxKmsT6dzaevsfyKzGW81vUFDK2ZTaf/hTktP8B43/yMlAF5ZyHLsrfw1sCFvNH3M97oO4e3By5mRe4uiu6SD/hLWMYM1m8/l406r8cL0ry88gD6CveNaIsQGtNbz0NgsKLvpFH5YMTinwNcisXba37G2qCbFV8V363/prZ6oYYhI+C7r6EZG41rcJlVSOlRKt1NJLL3e3Zwri192Qvw/KXM7rx9jSsLy9K/Y2n617TE/of1Wn/5gaFC0VvK0sz1LM/dhhv0DZ7XRQyJJJDFwXMRYxYzmw5jncQhmEMmN5ORj7PpjHtY2P0VFvWcgGw7j7b48PlgTWMmHc0/YmX/d8kW7yAZO3At34HwOeZHMMxNKRbvJJE4tqZnrA1CnwMiTuCtuQc1WVFiMQKB90bd0tU55X8RBP1ERpFY2fWW0Je9mKbYgcTs7UcsuzxzEUvTv6YtfjDrtf7yPYLo+hneHDiXpdkbAeiI7UFbdFdSka2IGrMGl0j9oETBfYt0+Wm6CvexqP/XvDVwIeumjmfd1LGD5Sy9k406r+WN7qN5s/dbCGHSGttnWNtS8cMZyF1N98BZJCJ7o2m1bZ6LRj9FNnMOvr9s3LO3CyHQjA3w3akrFmoYMgwyyCODlWhjuOlpJEqlBwEb267dU7QvezFSenSkThux3EDxIZYMnENrbP8PCEVP4WH+vXQflmZvZGbT5/j4rAf5aOcFTE8cQN5dwduZ23ip97e80nchi3N348gSM5s+x1bTr2W7de6kLboLbw38lseXfppM+b9dcl2LM7/jSuLWFrzVe9KIcxhC6HQ0n47nLydT+FPN70ckUhGkUp1csXVjAwL/7brU1QhUz2IYAv9dADRj3brUVy4/hmVvW7NTku/3kc7fQDJ2EOYIu1PL3lLe7D2JmLkJ67X9YlAopJQsyVzNG/0/J2FtzIK2y2iyNx287/GV32Zl4TEABAaSgFW+FhF9GrMS+zGv+Yt8tPN8eguP8Wrv6Ty94gts2vFrOmJ7ABXBmNdxMS8t34+FPV9lk+l3DbtKErN3JGItoC97Can4ETUNBQ1jQzRtOuXyY8Tj4++1KvR1kf6fkbI8ZoF+JhKqZzEMgV8JyaHpc+pQVx+e9yr2KIY86cKtSFmkpWnkXCbv9p+JlC6zW35OT/FZugqVcHiLM1fyWt/ZtEV3Y+vpN5B43+ayDVJfYLPW0/hI6w+Zkzye9VJfZ5PWH7Fp63dJ2RuxMH0197/7P7zS93uaozuw7Tq3kDDn82LXN+guPDj4HEvvZIP28yl5i1ie+f2wdgohaGk6Htd7i0L57zW9J0IIbHtHnPK/6xLrQjNmA5LAXzbudTUC1bMYBukvB0DoM9ZQcvQ47rMAWNbWNT8jk7+JiLnFiDtIM6X/Y6D4V2amvss7uXvJue9S9FbSkpvPivydRI152NYCPFnEDucJpJS8m3uYVwZupLv0fPgkwX8D2AhmxXdki/az6C4+yusDl9JTfJztpv+GBdOv5NkVR/NS98lsNf06kvZmQGXSszV2ICsyl9KROGxYx61EdB80rZl0/mbikV1rel8sayuKxVvx/aVjGg9kdWjhvIj0l4Ox3rjW1QhUz2IYAn8FYCC08U8k5DrPAgLT3KKm+x13EWX3ZZLxz4xYbmn6XCx9JtObjmFl4R98bMaF7Djj9yzJ3YutT2e91HEU3MW83Pc73CBL0evloWUn89iK08i5K2iyNkXT1yUfxCjIJIa+PilrE7qKL/D3lWfh0cKWHT8h7bzKP5Ydiy8dNp92EabWwsvd3yYYkpRndvP3QQiWp4fvXWjCpim6P7niXwlkuab3xrS2AsANBXk8WfXFUvnsTD2UWAyDDHoRWmtdgvO63mvo+rpoWqKm+/OlRwGIR/ccvozzArnyE0xrOhonyGMb7ZT9fl7tPZMAyQ4zLmfd5AFs1fkTBsov0Vt6nbve/QIrik+RsDahy82wsryUoh/gSBtX2hQCl6Wld+j3fJLWJizM3MkzvdewefsZFL1lPL7iZHQtycZtZ1Lw3uad9OWD9ljGDNpiB9FTuA3P7x/W7kR0D6QsUiw/UdN7Y5obAgLXfbWm+9cGTWsDKp+dqYgSi2GoiEVbXery3NfDD3VtFMr/xNTXxRphMrY7dwOaiNGeOBRTizM9tgv/WHoUi/OPkjTXpbf0MkVvBb3FpwF4ZPn3MUQUKaazsrQYSTN9nqDHdRjwdPpdQbdTJu0bCC3J0tLbRM2NKfkDPNP7BzZv/yF95Wd5vf9y2mI70xHbm3fSl+AOEYZpTV9CyjK9hT8Pa3fM/jhgUCjVFhRXiCi6PhevHkuaIgmYSiw+bMgggxiDoDNrrEdKfP9ddL32VZeS8wIRe8GIdaSLj5KM7IShJdG1CPNSX6TVmkVM09mi40zSzqu80ncBr/RfSMH3kNKnEJjkvAEyvkbGlxR9i37XIOtFyHg2A65B0bcZcCtlu8rvoOvrkHEXs6zwIuvE92RR+lrKfh/rNZ9AIEssy90yaFfM2piIsT7p4iPD2q5pcWxzQ8qj8Iw0jDn4/uI1FxwlQgiElkQ2INJ4PVBiMQwyyCJG4UlZfT39SFlEr3HyzQ8yeP4SbHP4ic2y9xaOv5RUpBLAR0qJLwsMlB5jZtNBtEYWMKdpf9aJ70lrZFf6nKWY+mxyfj8ZX1AObLKuTtaLkPai9Lox+t0YaS9G1rPJeBpOEAkF412S1sa8kfkzHdFd8WWZhQPXkrA2ojmyHUuzN75nZSIV3YVM+d/vmc94P7a5KWX35ZreHwBdn4XvL6n5/rVBiCRSKrH4cCELIOLjXs2qbF9ajROprvcOANYI+1cKYYq9uF2ZQBVCMFB8EolLR6wShDhpzWdabCdez9xJ0lyPLuddHGnjE6XoG+R8m4wXod+J0leK01OOMuBGyXoR8p5N3tfwZARda2F5eQmWluLt3GNMi32CJbl7kDKgM7Y3JW8JpSFLi3FrC6QsU3bfHtZ+05iL568gkKWa3iNN7yAI+gezpI0rIgZD3OCnEqNJX9gqhLhfCPFG+LslPC+EEOeH6QufF0JsNd6NGA+kLNZlq3EQVMbwmlZdsJj343qVb8yRHLGKzmuARmSIN2rGeQ6BTsregiAMrLs49yiu9yYuFoaIk/c9Ch4UfJO8bzHgROkuJliaS7Iin6S7GCPtRcj7FkXfpOhDv5tHSklzZCu27/weM+N7U/JX0ld+jpRd+SikS/+N+xA15wNQGGG3phkurXpebUHiNa0ZCOryjS9EBBkUxr2eRjCa9IWnAA9KKecDD4bHAPtSieg9HzgOuGjMra4LLjD+eSBk+C0kRG2em0GQBkAfYTLW8bswtTb00HdioPwq/eWXsfRp6FoMLQys219eiCMNep3lBEQIpEk5MCgGFhk3Qm8pRk8uQW9/nO6BBD35OH2lKDnfphhYlAMdsLH0FgbcLqJGG+3RSjyJrLOImFnp/ZSG/NFboci5/kqGQ9dagUrQ4VoQYQ9R1uMbX5gQJlqaalQTsHc5sDx8nRVCrEpfeACwa1jsauAR4Pvh+WtkZWD6byFEsxBiRvicyYMMqM8orRLZWqxF/MqhBLIQ3j/8ZitfZtHDbecLB64m7bxBd+llEkYTy3L3s7zwCO2Rbci6i7H1DtJunkKQww1sfKlRDgxKvkG2bJPN2si0CbokC0RMl4RZJq7r+FLHDTyKfgkvyCBlgK23ARpFbyW6FkETEbxgYNA2XVSWi/0gO6z9Wti2VW1dW0Qo+vXJSaqx6v90qjGa9IXThgjACmBVWOMpkb5QIqEe0aFH7YZcxf3Sh9BfZHHuL8xLHcl2086nGGj0lZ8jbsxice5uit5bQ/ZgyMF/JSClQMpQQz3AF0hf4Ac6SG0wImdl63qAlD4SiUBDoCHDGJ5DXxOeAS3cazIcq/4fapxzEKO8f62Z0An2amY06QsHr0kppRBird6hiZ6+UAi98kc27hWFQ50av/VW9ShG6mLrWhI/yFL0upAyIGVviONnkARs1vYdAMp+Hw8uPoSiV0QQIabFKWk+ugiwNQ9bd0nYDuWkQZ6K9iQSZZrsMrbhYmk+mgiwNI2YHiFmtJN136W39AoSj4jeSSAdfFnA0JoHbQtkHgg+EHBnKKvaptU4VJPhSkutvbe1IwDGP/5JI6g5fSGwctXwIsxO1hWenyLpCw3q0Z0Uo+xirxpe+EEak9VPchp6G57fi4ZOZ+xjvNL3ewJZQhcR3CCHqSUo+31YWhJJnhZzGlk/jy48LE0jork0mWX8QEMgsQ0PIaDJLtNiF0kYDlHNxdY8JA6+LGJIh0eWnUKbvQFSQsKcQ9GtdDjtIbk1nHCuwtCHn3Pxw2GLVmOE9f8KaR12gkoXpuCOUxhF+kIqaQqPCl8fBdwx5PyR4arIDkB60s1XAIgI1LgfYW3Qwm9ZWePknaFXRnjeCH4EMXMjJC5B0Mvcps+ScV4nYa7PR9u+x3PdZ/B8z895Y+Aqpsc+QVRvI66buEEWW0BMl8R0l7hRptku0h7Ns05TlplNGTqieZqtInG9IhZRHVK6jZAlZsc/zgFzbyKqSyw9RXt0W9LlindoaogDWTFcBYmO4MHqhm0z9bWPCQoQDIpNc033rw2S8phnup8oVNOzWJW+8AUhxLPhudOAnwM3CyGOAd4BVsXLvwfYD1gIFICjx9LgeiFEtC6z51o40x8EPTXdbxqVLfRO6G+xOlb9Ia4o/I1AtDIt9gl6S8/g+e9Q9t6lIF1mJ49hTurzBFoLz/VeSos5l6yXwZFlIlpllKgjsTSfcmAgAFvziGkOUd0lqnvolNCQtFidtEU2ouR18272YabHdkITJr3FRzG1VmLmf3dkFpyXAP09y7rvx/XeRdda0LTa/F6CoBchmusSHhFZhA+rWIyQvhBgj9WUl8AJo7Sr4QitqS5uu5XNRxaeX9tIzdDb0LVOys7wHo5Rc2N00cKi9C3MTR2DptnknRdB9tEa3QcokXWewPP3ZZPmz/N29m+U/DS60GkxJP04mJpAFz524NGklWg1c+R9m7I0MXUDnRJJHZrMaWza/ClWFJ7g1f5rKAUBWT/PS31X0pV/kLnNx75nc1669HcS9laDy7qro+y+jG1uOuz1NeF7S9GN8Q2rt4rKNoHx9/xtBMqDcxjq5bYrhFZxR/Zq37sQsTal5D4/7HUhNCxrS9xggPVThzCn6dM0W+uCiDE7+RnyXpqCuxhNRHH9bjZu2p7AH2Ca3Y4hdFoNnajmkDAcZtv9TLPSlAOTuOEy3UqT0AOShqTJ7OSA2eeSstalv/w6jreYTZr3Z9OWL7Akdy+ONJjZ9LlBd++yt5SC8wKpEZJOB7JM2X1ldGLhv4uhrznZ0WiRUlbEQtQne129UWIxDEJrQw6JaD2emOY8vFFEhY7aH8NxX8Xzu4Yt0xLbB4HLqz0/I+u8hRsU0ESUpDWPlLkO5cDjnfQlvNJ7GlJmmBWJIr3X2Tj5ceKaR0LL0W5qJHSPgrSRmoatRTGEgU6GdaILOGD2eSTMDiwtQslbhq418ZHWrxPVLLLOQtpjuxMx1hnMidKVvQbQaYsfNKzdpfKTSFkiFqktlaOUDp73JoYxv6b7166yAlBGjDBZO5lRYjEMQm9D+j11CcdmGBvjeYuQNU6oxsNv5nzxoWHLrJM4hGazjRWF+3mu5yw8NAwcXuw+mZX5O7D1BMuyN7F+8zfZqO3HfLTjPNrs9RDSZ491fsaM6Po06Qk0XJr0Ei16Dpte4kYbu087mf1n/5KokaK/9CJPrvwuUU1H09open082XUaEotNWr83aI8XpOnO30hLdO8RkxnlSg8BRrhVfe3xvDcBD8PcqKb71wYZVPyF6hXaoN6osHrDILRpgIOUA4TbXsYN09oc8HHdF2sKrWebm2Hqc8gUbiOVOHy1ZYQQzG/5Dq93f5OZsSOYlfoKmjBYkrmOgruclL0A39qAhLUxAJbRSd59iw5zMzqjm7Nl65d5N/8o06NfIu/1sDT/D7Zp/wYDzps833c5GyR3Y9HAtbyRvpqI3saO61xMT2kRjy0/gZLXw5btJxG3ZiKlRAjBsvQF+EGGGamvDdsuKX2yhT8Tj+xSc2Ag16nsQzHNzWu6f20IwmVgTZu2hpKTEyUWw6CFIdKkvwxq3ORVLZa1JQCO83RNYiGEIBn/LL2Zc3G9JcNuKmuNfRrbPJ83By4gwCRhzae7cD+zmg5meuIAFvX/mld7fsDM5BG8PfB7NCF4PXMLXaXH6IzuRl/pZSwtQSBdPj7tNJrM2RiawAvy/OXt/RBkSVhbsm7TgcSNWbxb+h0xuYgFnacxO/nFQVsr+VCvpD1+KHFrs2HbVSj/Hc9fTkfzj9f6PVmF4zyNEEmMOqR0qGfc1kagxGIYtHBNP/CWoI9icq0adH0Guj4Hp/xPSHy5pmek4ofTm/kt/dnL6Ww5Y7VlhBBs3PFbXlh+EMuyVzI98Tnmpr5CS3QHAOa1fI8l2etJl54iXX6WeS1nku27iqJf5OWBW9BFwOLMQnRRJp2/Gh+TUhAQSIvWyAI+2vYNAnSQBZ5ZcSQZ5znmt36f2cmjBm3wgzyLek5A15qZ1fydEdvUn70MXWsnEd2rpvcEoFz+Pyxru7qERwxCfxCtRn+QiY6asxiGVflC6pU0xrZ3plz+v5o3O5nGLJpiBzCQv27Q43F1xKxNWa/1R5hyBZbwaI5sx7Lsn1iZvwchdKbH98cN0nTG9mLd5P7sOuNXbNZ2AvvMvooNmvamSQ9IGO3kZTuaAKG1s27T/nxind+TMOeQLv6V51ceRcF9k806zmdO6n8H65bS562+71Hy3mSD9vMx9Y5h7Sw5L5EvPUxL0zFoNYYK8LzF+P5b2GHQn/Em8N5GaB2IGnO/THRUz2IYhNaM0FoIvEV1qc+O7Eah8Eec8r+xIzvV9Iy2phPIFm6nN33esL0LgI7EF8g7L7I88zs0YROP7MRrPaezNHM9hp4iaX2UdVOVHk7SmkWzvT6O30tENzD1Gewy+04KXg/PrTyejzQdStycw+u9Z7Eyfye+LDAt/mnmtXwf2/ivGFSE4rv0F/7C7OYfkIwMP2EppaQ7/RM0kaQ5cdSw5dZEqfQ3ACL2rjU/Y20I/DfRpmAKgFUosRgBzZhP4C2sS12RyK4IEaNYvLtmsbCtj5CKf57+3B9IJY7ANle/XCiEYG7rz5DSZWn6V8yQZbZd53bS5aex9LZBD0spJZqweCf7IAsHrsH1X2PDlm9VniGLmELjjb4zgAANi87Ep5jd9EWa7PeG+AuCEm/1fY++wp3MTH2b6cmRh1r50gMUSo/S0fz/0Efhol0q3o1hbIxhzqv5GdUipcR338CMfmrc62oUSixGQDM2xC3eOTiDP54IESUS2ZNi8W5SzWfWnP6uPfV9soW7WNn3HWZ33jasi7MQOuu1/RIhTPLlp3H9bnKlR0mXHqXJ3obOxNGsyN9DW3R7ys7jtJoGdmQn1k9VVlvy7pvoWpwZic/QGduLVGRLjNV4Lpbct1jY8zWK7ivMav4+M5JfHdF+P0izsv80LGMeLYkv1fQeAHjeEhzncZqaRp4XGStksBJkGt2oPUr7REfNWYyAbm4CMoOsQ2RogGjsMKTsp1T8a83PMPR2Olt+StF5gr7shSOWFUJnbuvPmddxOd2563D95cxvvwSBSbr0MJrQ6M3fhCaX4LvPopMDIOe8TtZ5gebINnyk/SfowuPt3u+xqOcbrMxehRdkCIISy9IX8NKK/XD85czvuHKNQgGwsv8HeP4KpreehxC1RyorFCoZ4KOxQ2p+xtrguy8BoI0QOHmyo3oWI6CbHwXAd19EM8Y/56lt74SuzyKfv5ZobP+an5OMHUy+eD896V8StbYmFtlx2LJCiIqXY5AhFd0dy5hBMvJxiu4iXPc5DK2VDdouYEnmOnrz19GVvZEVhceImxsyPfFpAumgESUV3Q1DS9FXuBspXVZmr8bxl9AS3YfZLT8a0fFqFQO568gWbqMt+W2idu2hW6V0KRRuwLZ3GveUhasI3BcBUfmCmaKonsUIaObGgIXvPLPGsmOBEDrx+NE4zj9xndrzZAghmNb6SyxjfZb2HkfZfWPE8rpI4vor8PxKpPGy9y4l9w2K7kJMvR0hBOs0HUwgcyzO/IGEtRHrpo4jbm2AJiyaItvRkTiUltjeJOytKXlv0hbbn406r2dex8VVCUW+9Bgr+08lHtmNtuS3am47QLF4N4G/nHji2FE9Z23wnWfQjA2m7CYyUGIxIkLY6Oam+O7TdaszFv88QsTJ5YbPAVoNutbEzI5rEOgs6f7ciFvYhRB0NH2BgeKDvNF9LJnyv/Flnvb4Z3D9XrpzNzBQvB8hc8xJfZn1mr+BrkWG3K8P5v0ouW9hap3MbP72iCseQymWn2BpzzFY5nxmtP1+VFvJpQzI5X6PYWyAbe9e83PWrk6J7z6Dbm5Zl/oahRKLNaBb2+I7zyFrzFmxtmhainj8aIrFO0adn9My1mVWx/UEssjirs+M2MNI2jsyp/VMZqZOImnvgKl30Brbn1RkJ9LFvyOlg6G1YOnTPzDZG0gHTVh0526k7C2mObpH1X/whdL/saT7CAx9GrM7rkMfZRa4UulePPclEk3frIsjFkDgvY4M+tHt7epSX6NQYrEGdHsHwMF3nqpbnYmmryJEgmzm56N+VsTajNkdNyGlw7tdB1Io/XO15SRlIsYcYtamDBQfImFtiam30Rz7JPM6fk8y8glsYz0sY5Urc0UwpPTRhEVX9jr6CnczPfll4nZ1+zDS+T+xpPvzGPo6zO74E4Y+fVRtldIlm/kFhjGPaHT4naxjje/8GwDD2qFudTYCNcG5Bgxre8DAK/8dwx5+onAs0bQWEk0nks2cTan0CJHIrqN6XsTajDnT7mJp95Es7j6c9tR3aW064T3f/ot6vonrr0AIm4S1Jc2xvciWH6fgvIypd9BfuI+2+KeJGHMBBnsXQuj05G4lV36Sua3nVDU/EQQFugbOJJ2/hpj9CdZpv3RU/hSryOevxvPeoKX1D/WJihXilR9D6LMR+vhPgjcS1bNYA0JrQre2wis/Vtd6E4nj0PW5ZNKn17x1fSiWsS5zpt1NU3Q/etI/Z3HXwe8ZlqzX9gvWbzuX2c3fZ53Ut9CEhSYiOP5S+ot/Ixndibb4gQhh4AVpXu/6X/ygSNFdyFt93yHvPM/K7BUsGfgFBeeVYe0olP/D2yv3Ip2/lpamrzCr47oxEQrf7yab+RW2vSuRyN6jfl61SOnglf8Pw9553H1xGk01AXv/IIToEkK8OOTcGUKIpUKIZ8Of/YZcOzVMXfiaEKJ+/2vjiGHvRuC+SODXL+6wEDap5p/geQvJZs8bk2fqWhMz2i5meuvvKLuv8faKPejq/zF+MIChpYiY65Owtx4McZewt2JOy4+Z134hnYnPDz7H0FLM77gcXYti67P4yLRbWLf1rEr4Pi31vrwgFVxvCct6v8riroOQsszsjpvpbD59VL4Uq5BSkh44FSnLpFJn1vWP1i//G2QeI7Jb3epsFNUMQ64CLgCued/5c6WUvxp6IkxreDiwKbAO8IAQYkMp65GAY/wwIntRzp6DV7ofK35k3eqNRHYnGjuUXPYCIpG9BreyjwYhBKn4Z4hHdqEn/Qv6c5czkL+e5vgRtDQdO2LO1Pc+p/I9o2kREvbw2+pLzsv0Zy8hU7gdIQzakifR2nQC2hhutioWb6dUuoem5Gl1ce0eilu6H0QEw6rNRX8yscaehZTyMaDa+HIHADdKKctSyreoRPie9FPEmjEPzdgAt3hP3etOpf4fuj6d/r6vDoa0HwsMvZ3prb9g7rQHSET3pj93BW8u357FXYeSzt2IN0Lu0TXheG/Tn72ct1fswzsrP0m2eDfNiaNYb/pjtKe+O6ZC4bmLSA+cgmVtSyKxZg/RsURKH690L4a9C2KEgMNThdFMcJ4ohDgSeJJK4uR+KmkK/z2kzGpTF0IlfSGVxMnMmTOxJ4aEEBiRT+Hkfkfgd6HpnXWrW9NStLReTE/3wQz0fyucvBu7qSbb+gjrtF2AmzqFdP4mMvlbWNF/MvSDZW5MxNwc29oE05iDqc9C15rDvBgSKYv4QR+utxjHe4uy8xIl51lc/93Ks83N6Gw+k2TsYHS9dcxsXkUQFOjrOw4hTFpaRuefUQu+8wQy6MKMfrqu9TaKWsXiIuAsKkkdzwJ+DfzviHe8j4mevvD9mNH9cXK/xS3eiV1Hz0AAy9qaZOpHZNJn4DpPYdnbjnkdpjGL9tS3aUueTNl9kXzpMYrlf5EvPUKmcHN1z9BnY1ub0WIfTzyyC1aYNX08kNJnoP8EPO9VWtuuQ69iFWascYu3g4hi2J+se92NoCaxkFIO9lGFEJcBd4eHUyR14QfRzflo5gLcws1Y8WPqPvMdjx+DbX8MM4za5bkL8fzF2PYuY9rTEEIQsT5KxPooq9K/eH4vrrcYz19KEGTCVIsamhZD01KY+ixMY85gKsV6kEmfQan0V5Kps0a9tFwLMijgFu/CjHxqyga7eT81icWqHKfh4UHAqpWSO4HrhRC/oTLBOR94fNRWThCs2KGU0qcRuM+hWwvqWrcQAtPcFCklTvkfZDLnoOszSA/8kFTzWUQiu4/bVnpDbwtzkS4Y82fXQibzK/L5K4jHjyOROKYhNrilv4DMYcYOXXPhKUI1S6c3AP8CNhJCLAnTFf5CCPGCEOJ5YDfgJAAp5UvAzcDLwH3ACZN9JWQoZvQAEHGc/PsXhuqHlFlKpQeIxg6ite0yUqkzKJceDa8VcNewaWyyk8n8ilz2N0Rjh5FMnd4wO5z8NWjGPHRr0s/fV0016Qs/t5rTV4xQ/qfAT0dj1ERFaE2Y0YNxCzdjJ09D09sbYke5/DBtiRsACGQGw5xHPn8djvMErvMiuj6NpuR3serc+xlPpAzIZM4in7uEaOwwmpt/Xbe9H+/Hc54hcJ8jUmefjkajPDjXEit+NFDGyV/dkPqFMLHsHSkU/4zrvkqxeCe+v5JM+kxiscPpnPYAlr0DrlPZKVuPJEnjjZRF+vu+Qj53CfH40Q0VCgAndwmIJGb0sw2zoREosVhLdHMeRmQv3MI1yCBf9/qFiJJIfBnfX0Iu+zsEBkHQi2lujm1XNjJFIvvguq8iZfED33yTTTw8bzE93QdTKv2FZPLHJFM/aahQBN7beKX7sOJfQNSY+GiyosSiBqzE15BBP07h2rrW63mLcd1XMIz1aW4+G8PcCN2YQyLxZQxzPsXiXfh+N7ns7whkBiGiHxAHIQSet6SudtdKsXgf3V174Xlv0dJ6BYmm4xve7S9nfweYYQ/zw4XadVoDhrUVur0LTu5irFj9vmFc52kymV9iWpuiaW34/lJSqR9jGOtj2ztSyF+L675MsXgXbe03hXcFSCkQQsN1X6JUeohy6TEkJeLxI4nVKUbl2hAEA6TTZ1As3IxpfpSW1ksxwjwujcT33sIt3oYV/xKaPjVTFI6E6lnUSKTpO8igj3LukrrVGY0dQOe0B7GsbYlE9qSl5SKESCClQzT6P7S134iudWKam2GHgViE0BFCCzdb/QhNJGhtu4qmppNx3Vfx/e662b8mpAwoFP5E18pdKRZuJdH0Tdo77poQQgFQzpwDwsaqs1v5REH1LGpEt7bAiHwaJ38pVvwItFEGbqkWIWwSoQeplAHF4p8p5K8iGj0Q31+OlEWakt8OrzsIYQFQKt2NlCUMcxM0LU4kshu2PXE2P5XL/yKT/hmu+xSmuRWtbddgWeOfzLhaPOcpvNI92E0n1dXdfyKhehajIJL8HkifcubshtQvhEY8/nlSzWfj+0vRjbk0JU8d9GhcJRRSSnLZi9C0ZkrFe+nu2g+n/CRCGAjRuO8LKSXl0t/p6TmM3p7P4PtLaW4+l/aOOyeUUEjpU0qfgdA6seLHNdqchqF6FqNAM9bFShyHk7sAM/Y5DLsxYdUsa+sPZF8vlR5E12dgmpsQBMsJgn5a265F19vIZM7B9V7FsrepSwKl9xMEGYrFOyqRrdyX0bR2kqkziMe/GG5Sm1i4hRsI3OeINv8WocUbbU7DUGIxSuzE13GLd1BKn0q8415EjUl8xxopHYIgC4Drvoptfwxdb0NKH0NfD8d5koz3M2z7E1j2x8YkCM1IBEGBcvkRSsW7KZXuqwyJjE1obv4N0diBE+Z9ez+Bv4JS5ufo1scwogc22pyGosRilAgtSjR1NoW+L1DO/oZI8rRGmwRANLrv4GvL2p587gpy2d+j6R2Uy//AMOaRy51PLncBQiSw7I9h2ztimVthWh+tOX3iKoIgh+s+i+M8Tbn8D5zyE0AZTWslGj2EWPxwTHNBw5dCR0JKSWngNJBlIs0/n9C21gMlFmOAEdkZM3Y4Tu5SDPuTGBMsJLymxUmmfkw+fxWB+zzx+JFY1tYkmr5MqfQI5dIjlMuPUS7dH96hYxjrYxgbohuzMfRZaFpbJbO8iIZzIRIpHaQsEAT9BEEPvrcEz38Hz30d33+XSgQDMIyPEE98iYi9e9iLmRwfO7d4M175Aezkj9CncHb0ahETwaNvm222kU8++WSjzRgVMsiR694X8Eh03IsYgyC044GUHkIYq52r8P0VOM4zuM5zeN5ruN4b+N5SoLqAwULE0PXZGOaGmMbGmNYCLGsBmtYyDi0ZX3xvEfnu/dDNLYm1XVf3wDrjgRDiKSnlNrXePzkkfhIgtASxlgvI9xxEceD7RFsunpDd1lXf6quzTdenE43u+54hjJQBQdAb9h76K8mWpANCQ2AiRBRNa0HTWhFay4Rs89oiZYli34kIESHact6UEIqxQInFGKJbW2AnT6WcOQsndxF209cabdKoEUJD1zvQ9Y5Gm1IXKvMUpxJ4LxFtvbJu/jOTAeVnMcZY8WMxovtTzp6DW3qw0eYo1hInfzlu8Vbspm9jRvZotDkTCiUWY4wQgmjql2jmphT7T8B3a8+GrqgvbvGvlDM/wYjsi5X4eqPNmXAosRgHhBYl1voHhNZMofcoAu/tRpukWANe+XGK/SeimwuINp/X0G3wExX1jowTmj6dWOu1IF3yvUcQ+MsabZJiGHzneQp9R6PpM4m2XvmhyAFSC7WmL2wVQtwvhHgj/N0SnhdCiPPD9IXPCyG2Gk/jJzq6OZ9Y27XIYIB8z2FKMCYgvvM8+d4jEFqKWPv1aOOQ32SqUE3P4ipgn/edOwV4UEo5H3gwPAbYl0pE7/lUEghdNDZmTl50a4tQMHrJ9xyihiQTCK/8BPnezyG0JuJtN6Hp6zTapAlNrekLDwBWBaG8GjhwyPlrZIV/A81CiBljZOukxbC2It52AwRZ8j2fUZOeEwC39ACFviPQtHbi7begGbPXfNOHnFrnLKYNyRuyAlgVNmgmsHhIuRHTFwohnhRCPNndPXECsIwXurUFsfZbQZjkez6rllUbiJO/mmLfsWjGhsTab1U9iioZ9QSnrPiLr7XPuJTyUinlNlLKbTo6PhwOP7o5n3j7n9GM9Sn2/S/l7IWTLoDuZEbKMsWBUyilf4Rh7x4OPRqTzmEyUqtYrFw1vAh/d4Xnp2z6wrFC06cTb7t10HGr2P8V5BhmR1esnsBbTKHnMNzC9ViJrxNtvexDHZuiFmoVizuBo8LXRwF3DDl/ZLgqsgOQHjJcUYQILUq0+Xzs5I/wSveT694XrzxlsjxOONzineS698X33iDacjGR5HfVfo8aqDV94c+BPYUQbwCfDI8B7gHeBBYClwGTf3PEOCGEwE58mXj7bYBBofdQSumfVjZqKcaEIOin0P+NirOVMY9Ex72Y0f0abdakRW1RnwDIIEcpcyZu4UY0YwMiqbMbFqJvKiClxCvdRSl9BjIYwG76OlbixHGPBjbRGe0WdeXBOQEQWoJo8y+Itf4RKR0KvYdS7P8mgb+i0aZNOnz3dQq9n6fYfyJCn068427sppM+9EIxFiixmEAYkZ1JdDyAlfg6bvEv5Lp2oZT5NTLINdq0CU/gd1EcOIV891747gtEUmcRb78L3dyk0aZNGZRYTDCEFiWS/C6JzgcxInvi5H5LrusTlLMXNSS36kQn8HsopX9KrusTuIU/YcW/RKLzMaz4UWoSc4xRwW8mKJqxLrGWC/DjX6aU/RXl7Nk4uYsw40dhxY/60PsHBN7blHOX4xZvBulgRg+oJAAy5jbatCmLEosJjm5tQbztWjznaZzcRTi583FyF2NGP1URDXOLKRHKrhqkDPDLf8fJX41XfhAwMWMHYSW+gm5s0GjzpjxKLCYJhrUVRutl+N4inNyVuMVbcYu3ohkbYcYOw4zuP2XT6gXe2ziFW3GLtyD9pQitHStxIlb8yA9lguJGoZZOJykyyOEW/4xb+BO++wwg0K0dMKP7YUT2QtMn9/49312IV/orbukvBO6LgEC3d8KKHYoR2WcwNaOieka7dKrEYgrgu2/gFu/CK91F4C0CQDM2xYjsjGHvjG5tXXXGr3LuCtzCDYDEjH0OO3EsvvsypYHTkDKPps8i2nI+Qmtazb2Xh/cKNHNjos2/qrpeGaTxnMfxyn/HLz1K4L9VaYe5BWZ0f8zopya9ADYaJRaK91D5Rv4bXvlhfOcpwAMsdGtzdHNrdGsBurkZQp/zwbwh7msU+08g3n4XCJNC7xeJNp9Nof9EIskfYtg74BRuIvAWE0l+5z33Bv4K8j2fIdH5IEJEKPR9FSOyO1bskA/YKKVP4L2J776A7zyD7z5F4L4MBCAiGNbHMSK7YUT2VDtCxxCVN0TxHnRzHro5D7vpa8ggi+c8jl/+N77zOE7+Ssg7lYIiiW7MRzPnoxnro+lzkUE3urXlYFg5w94Bt3gvgfcWurV9eG4nCrkvwPvEAgDpgSwhMUAW0bRp+O5CAu9lAu9dAm8hvreQwHsdVrm1ixi6uQVW4usY9g5r1QtS1BclFlMYoTVhRvYYDGkvZZnAfRXffRHffZnAex2v9DdkUIltZMa+jF9+nCDoRxDBKz2MZm2ObmyIV/obZnRv3OJfCPwP7g3U9OlYiePIrtwBISLo9s4YkZ0p9p+MW7wltGcamrkhVuwINHNTdHMzNGO+8oeYJCix+BAhhI1ubYFubfGe8zJIE3jvIPR2PHM+hd4jECKGZm6CQCPS/EtKmR9Tzv0WM7Ingg+6TstgAK90P4nOfyK0JMX+r+IUbsNKfA0rcRyaPgehxerVVMU4oMRCgdBS6NbmAFjxw7HihwNQypyDps9AN+cRb7sOAN97E7f00Aee4ZX/gWbMRtPbADAi++A7T2HFDq5TKxTjjXL3VryHwO+p/PaW4pXuw4weMHhOygAnez5W/AsfuE/oM/Gdp5FBESklfvmf6Ma8utquGF9Uz0LxHor9xyODfsAkkjoLoaVwclfg5q8BwIjugxk9FKisgJQGvk+s7WoMa0uMyH7ke/YDdHRzU8z45xvXEMWYo5ZOFYoPCSqehUKhqAtKLBQKRVWMas5CCPE2kAV8wJNSbiOEaAVuAuYCbwOHSin7R2emQqFoNGPRs9hNSrlgyFhouNSGCoViEjMew5DhUhsqFIpJzGjFQgJ/E0I8JYQ4Ljw3XGpDhUIxiRmtn8UnpJRLhRCdwP1CiFeHXpRSSiHEatdmQ3E5DmDOnDmjNEOhUIw3o+pZSCmXhr+7gNuB7Rg+teH77/3Q5TpVKCYzNYuFECIuhGha9RrYC3iR4VMbKhSKScxohiHTgNvDACoGcL2U8j4hxBPAzWGaw3eAQ0dvpkKhaDQ1i4WU8k1gi9Wc7wX2GI1RCoVi4qE8OBUKRVUosVAoFFWhxEKhUFSFEguFQlEVSiwUCkVVKLFQKBRVocRCoVBUhRILhUJRFUosFApFVSixUCgUVaHEQqFQVIUSC4VCURVKLBQKRVUosVAoFFWhxEKhUFSFEguFQlEVSiwUCkVVKLFQKBRVocRCoVBUxbiJhRBiHyHEa0KIhUIIlcJQoZjkjItYCCF04EJgX2AT4HNCiE3Goy6FQlEfxqtnsR2wUEr5ppTSAW6kkgNVoVBMUkabvnA4ZgKLhxwvAbYfWmBo+kKgLIR4cZxsaTTtQE+jjRgHVLsmHxuN5ubxEos1IqW8FLgUQAjxpJRym0bZMp5M1bapdk0+hBBPjub+8RqGLAVmDzmeFZ5TKBSTlPESiyeA+UKI9YQQFnA4lRyoCoVikjIuwxAppSeEOBH4K6ADf5BSvjTCLZeOhx0ThKnaNtWuyceo2iaklGNliEKhmMIoD06FQlEVSiwUCkVVNFwsJrNbuBDiD0KIrqE+IkKIViHE/UKIN8LfLeF5IYQ4P2zn80KIrRpn+cgIIWYLIR4WQrwshHhJCPHN8PxUaFtECPG4EOK5sG3/Lzy/nhDiP2Ebbgon5hFC2OHxwvD63IY2YA0IIXQhxDNCiLvD4zFrV0PFYgq4hV8F7PO+c6cAD0op5wMPhsdQaeP88Oc44KI62VgLHvBtKeUmwA7ACeH/y1RoWxnYXUq5BbAA2EcIsQNwDnCulHIe0A8cE5Y/BugPz58blpvIfBN4Zcjx2LVLStmwH+BjwF+HHJ8KnNpIm2pow1zgxSHHrwEzwtczgNfC15cAn1tduYn+A9wB7DnV2gbEgKepeBf3AEZ4fvBzSWVF72PhayMsJxpt+zDtmUVFxHcH7gbEWLar0cOQ1bmFz2yQLWPFNCnl8vD1CmBa+HpStjXsnm4J/Icp0rawq/4s0AXcDywCBqSUXlhkqP2DbQuvp4G2uhpcPecB3wOC8LiNMWxXo8ViSiMrsj1p16aFEAngVuBbUsrM0GuTuW1SSl9KuYDKN/F2wMaNtWj0CCE+BXRJKZ8arzoaLRZT0S18pRBiBkD4uys8P6naKoQwqQjFdVLK28LTU6Jtq5BSDgAPU+meNwshVjkpDrV/sG3h9RTQW19Lq2JHYH8hxNtUdnnvDvyWMWxXo8ViKrqF3wkcFb4+isp4f9X5I8OVgx2A9JAu/YRCCCGAK4BXpJS/GXJpKrStQwjRHL6OUpmLeYWKaHw2LPb+tq1q82eBh8Je1YRCSnmqlHKWlHIulb+jh6SURzCW7ZoAkzL7Aa9TGTf+oNH2rKXtNwDLAZfKePAYKuO+B4E3gAeA1rCsoLLyswh4Adim0faP0K5PUBliPA88G/7sN0XatjnwTNi2F4HTw/PrA48DC4E/AXZ4PhIeLwyvr9/oNlTRxl2Bu8e6XcrdW6FQVEWjhyEKhWKSoMRCoVBUhRILhUJRFUosFApFVSixUCgUVaHEQqFQVIUSC4VCURX/HwG/+CSRrEaEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hds = flopy.utils.HeadFile(os.path.join(workspace, name + \".hds\"))\n",
    "h = hds.get_data(kstpkper=(0, 0))\n",
    "x = y = np.linspace(0, L, N)\n",
    "c = plt.contour(x, y, h[0], np.arange(90, 100.1, 0.2))\n",
    "plt.clabel(c, fmt=\"%2.1f\")\n",
    "plt.axis(\"scaled\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:50.484840Z",
     "iopub.status.busy": "2022-07-30T20:16:50.480024Z",
     "iopub.status.idle": "2022-07-30T20:16:50.792229Z",
     "shell.execute_reply": "2022-07-30T20:16:50.791312Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7eElEQVR4nO3dd5xcVfn48c9zy8ydsrM9oSRI7woCQlRERUEEpFcFEVGKwFcUQVERQVFRFEEUBEGKIiDlR5cuggpIEykikZqQZHuZetv5/TF34xKS7GRndmZ2c96v17525s695zx3s3n2lnPPI0opNE3TJmI0OgBN06YHnSw0TauIThaaplVEJwtN0yqik4WmaRXRyULTtIpUnCxExBSRp0Xk9uj9eiLymIjMF5HrRCQWLY9H7+dHn687RbFrmlZHq3Jk8WXgxXHvzwHOU0ptCAwCR0XLjwIGo+XnRetpmjbNVZQsRGQOsAfwm+i9ADsDN0SrXAnsE73eO3pP9PnHovU1TZvGrArX+zlwKtASve8EhpRSfvR+AbB29Hpt4E0ApZQvIsPR+n3jGxSRo4GjAVKp1LabbrrpJHdB07RKPPnkk31Kqe7Jbj9hshCRPYEepdSTIvKRyXa0LKXUJcAlANttt5164oknatW0pmnLISKvV7N9JUcWHwT2EpHdAQfIAOcDbSJiRUcXc4CF0foLgbnAAhGxgFagv5ogNU1rvAmvWSilTlNKzVFKrQscAjyglPoM8CBwQLTaEcAt0etbo/dEnz+g9NNqmjbtVTPO4uvAV0VkPuVrEpdFyy8DOqPlXwW+UV2ImqY1g0ovcAKglPoz8Ofo9SvA9stZpwgcWIPYNE1rInoEp6ZpFdHJQtO0iuhkoWlaRXSy0DStIjpZaJpWEZ0sNE2riE4WmqZVRCcLTdMqopOFpmkV0clC07SK6GShaVpFdLLQNK0iOllomlYRnSw0TauIThaaplVEJwtN0yqik4WmaRXRyULTtIpMmCxExBGRx0XknyLyvIicGS2/QkReFZFnoq+to+UiIhdE5QufFZFtpngfNE2rg0rm4CwBOyulsiJiA4+IyF3RZ6copW5YZv1PAhtFXzsAF0XfNU2bxiopBaCUUtnorR19rWxq/72Bq6LtHqVcX2TN6kPVNK2RKq11aorIM0APcK9S6rHoo7OjU43zRCQeLVtavjAyvrShpmnTVEXJQikVKKW2plx5bHsR2RI4DdgUeB/QQbmOSMVE5GgReUJEnujt7V21qDVNq7tVuhuilBqiXIlsN6XUouhUowT8lv/VEBkrXzhmfGnD8W1dopTaTim1XXf3pGu1appWJ5XcDekWkbbodQLYBfj32HUIERFgH+C5aJNbgc9Gd0XmAcNKqUVTELumaXVUyd2QNYErRcSknFyuV0rdLiIPiEg3IMAzwLHR+ncCuwPzgTxwZM2j1jSt7iZMFkqpZ4H3Lmf5zitYXwHHVx+apmnNRI/g1DStIjpZaJpWEZ0sNE2riE4WmqZVpJK7IdpqQCkPFSwiDBahwiWEQR8qHECpEQhHUKoAqoTCB+UDAmIi2CAOIknEaAHJYJidiNGJmGtgmGshxmzKN9O06Uwni9WMUkVC7yUC73lC/z8E/suE/muoYCEQLrO2AdKCGBlEnHJSwAYxQClQihAXVBGl8qhwFNTocnq1Mcy5GNZ6GPZGGNYmmPYWGNYGlJ9N1KYDnSxmuDBYTFB6FN/9B4H3NKH3IhCUP5QEhrUhZmxbDHM/DGsOhrEmYq6BmN2ItCKyameqSvmocBAV9hEGi1HBQsLgTUL/dUL/FfzSw4AbrR3HtLfEjG2DGdseK74DYrTVcO+1WtLJYoZRyiUo/R2/9CB+6SFC/7/lDySNaW9NLP2l8n9Qe3PEnLvKyWAiIhZidoPZjWlvtpz4fEL/FULvRQLvWQLvadzcVZC7FBAMe0us+IexnJ0x7W1qHp82eTpZzABKFfGLD+EVb8MvPgAqC8Qx4/OIJz+NFZuHYW/eFNcNRCxMe2NMe2Ns9gZAqRKB+wyB+yh+6S+42YtwsxciRheWsyt24lOYsXlNEf/qTCeLaUopRej9Ezd/HV7hVlCjiNGBndgTy9kVK75j+TrDNCASx4rvgBXfgXjLl1HhMH7xz3jFe/AK/w8vfw1izMJO7oedPBjT2qDRIa+WdLKYZlRYwCvcjJu7mtB/HsTBdnbHTuyHGf8AItP/n1SMVuzk3tjJvVFhAb/0AF7hJtzspbjZizFjHyCWOhzL+cSM2N/pQv+kp4kwGMDLX4mbuxIVDmBYm+G0/gA7sXf5luUMJUYCO7EHdmIPwqAHL/9H3PzvKQweh5hziae+gJ08BDESjQ51xpPyc1+Ntd1226knnnii0WE0pTAcxM3+Gjd3Bag8VvxjxNLHYMZ2oDw7wOpHqQC/eC9u9tcE3pOI0UUsfRyx1OHT5tSrEUTkSaXUdpPdXh9ZNCmlirjZyyllfwkqi5X4FPH0iZj2Jo0OreFETOzEbtiJ3fBLj1Ea/Tmlke/hZi8lnjkFO7Gfvhg6BfR9qSajlMIr3Em256OURn+EFdueVPc9JNsv1IliOaz4DqS6/kCy83rEXIPi0MnkevfEL/2j0aHNODpZNJHQf5P8wBEUBo9FJEOy81qSnb/VSaICVnweqa5bSLRfiFKD5Pv3pzB0CiocanRoM4ZOFk1AqRA3dzXZ3l0J3H8Qz5xBqvsOrPgHGh3atCIi2Im9SHc/QCx9LF7+BrI9H8cr3tPo0GYEnSwaLAz6KAwcSXH4W5ixbUh330s8fZS+JVgFMZI4mW+S6roVMTopDHyBwtA3UGGh0aFNazpZNJBfeoRc7yfwS3/Daf0+yY7fYVhzGh3WjGHG3k2q+zZi6ePw8teQ69uDwHup0WFNW9XUOl1PRB6LappeJyKxaHk8ej8/+nzdKd6HaUcpRSl7Mfn+wxCjrfwLnfrsansrdCqJxHAyp5Hs/D0qHCbXtzde4Y5GhzUtVXJkMVbrdCtga2C3aIr/c4DzlFIbAoPAUdH6RwGD0fLzovW0iFJFikMnURr5AZazO6mu2zDtTRsd1oxnxT9EqvtOTGszCoPHURz5Kc0wxmg6qabW6c7AWFHkKynXDoFyrdMro9c3AB8T/ScTABUOke8/DK9wM/GWU0i0/xIxko0Oa7VhmLNJdl2LnTgIN3s+xaGTUMqdeEMNmGStU+C/wJBSyo9WGV/PdGmt0+jzYaBzOW2uVuULw6CHXN+BBO4zJNovJN5yoj7taACROE7bT4i3nIJXuJn8wBf0hc8KTarWKeUap1VZncoXhsFi8n0HEgZvkuy8AjuxV6NDWq2JCPGWE3Faf0xQeoj8wBGoMNfosJreZGudvh9ok//d3xtfz3RprdPo81agvxbBTkdhsIR838GEYS/Jzt9hxXdsdEhaJJY6hETbBQTu4+QHjtRHGBOYbK3TFyknjQOi1Y4Abole3xq9J/r8AbWaXkkqX6P4DGG4hGTnVVixST/Do00RO7k3ibbzCdzHyA8eo69hrEQ1tU5fAK4Vke8DTwOXRetfBlwtIvOBAeCQKYi76SlVJN9/JKH/GsnOK3WiaGJ2cm+UylMc/jrFoVNw2n6uryctRzW1Tl+hfP1i2eVF4MCaRDdNKRVSGPwKgfcUifZfYcU/2OiQtAnEUoeiwl5Ko+ci5jo4mZMbHVLT0WOKp4Cb/QV+8Q7imW9iJ/ZodDhahWLpEwn9N3Cz52Pam+p/u2Xo4d415hXvpzT6M+zEfsRSxzQ6HG0ViAhO29mY9rYUhk7WQ8OXoZNFDYXBWxQGT8Kwt8Bp+5E+752GROIkOi5GJEVh8Dh9h2QcfRpSI0oFFAZPArzyyMxpPr2bUoow7CUI3iTwFxKG/YThEEoVUbiAIMQQI4FhtGMYXZjmHCxzDobZ0ejwq2KYs0m0/5x8/2EUR84k0fajRofUFHSyqBE3+2sC91Gctp9iWus1OpxVFgRvUSo9iuc+hes9i+/9p1zn9B3spY/PK+UB/jvWMIxOLGtTYrGtsGPbEo/tMO0SiBX/ELH0MbjZi7GcnbGdXRsdUsPpZFEDgfcypdGfYTmfxE4cMPEGTUAphec9Q6FwO6Xi/fj+fwAQSWDb7yGR3Bfb2gjTWgfTXBvD6MIw2t5Rm1SpImE4RBj04gcLCPw38P3/4Hkvks1eCvwKEGx7Kxzn4ziJT2HbG9V/hych3vI1/OJDFIe+iTVr+9W+tKKe3btKSoXk+/Yn9P9Latb9GGZzD10P/IXk89eRz19PELwB2MTj7yce/wjx+Aex7M1qNvGOUkVc91nc0sMUS3/Gc58CFJa9OcnkISST+2MY7SvYVjXFNZ/A/Re5vr2wkweQaPtJo8Opip7du8G8/PUE3pM4bT9t6kThuk+TzV5MsXAHEBKLf4iWlq/gJHbDMFqnpE8Rh3h8e+Lx7WnhZIJgMYXCHRTyNzAy/B1Ghn9AMnUQ6fTRWNb6S7crr3cbYThIOn3MlMVXCTP2bmLpL+JmL8ZOHIQVf1/DYmk0fWRRBRUOk+35MIa1AcnOG5riL+GyXPefjI78iFLpIUQypFKHkUx9Fstap6Fxed7z5LK/JZ+/AVDMXuMJTLMb33+V0ZHzEHEwjHZKpb/S3n4+lt24koUqzJHt/Rgi7aS6b5+2ZQb0kUUDlbK/RIWDOK1nNl2iCILFjAyfRaHw/zCMdjKZb5NMfRbDSDc6NABsewva2s+lJXMqbulRzOiozHWfRuHS1nYeIiaDA18mCJdgsQFhmMdowPwfYqRwMt+kMHgCXuEmYsnVc4CyThaTFPoLcbO/LdcYtbdsdDhLKRWSy13O6Mg5KOWTbjmJdPo4jCYtcWias0gk//fIvgpHMM21ETHx/VcxzW4Mo42RkR/je/MRI0lr6xkrvNYxVSznUxj2pZRGfoKd+NS0vzU+GXpQ1iSVsucDinjL1xodylK+/wb9fQcyMvwdYrEdmDX7QTKZU5s2USxPLLYtbulhRkZ+zMjwWSjlMTR4Mr73EpnW72AasxgcPAmlSnWNS0RwMt9ChYtxc1fXte9moZPFJIT+G3j5G7BTn8aw1p54gzooFG6nt2cXPO852trOo6Pzaixr3UaHtUrCMIcdezdd3XchWOVblWKglEem9btY1hzSLV/GtjZBKbfuc2ha8XmYsQ/iZi9eLUd26mQxCaXsJYBBPP2lRoeCUgEjw99ncOBoLGsjumfdTzJ1cNNdQ5lIoXAn/X37Uir9FaWyBMFCTGM2qdThOIld8f3ycxqe9y9c9wnAb8g+xltOQoW9eIXr6t53o+lrFqsoDAbwCtdjJ/fDMNdoaCxKFRgcOIFi8S6SycNpbfseUUWGaSeR2B0IGR76Dra9OZa1Aan05zGMDJa1ISMj5+C5/yKXuyy6UNuOUiEi9f17Z8a2x7S3pZS9FDt52GpVDEofWawiL/87UEViqS82NI4wHKG/71CKxT+RaT2LtvZzpm2iGJNI7Mms2ffT1v5TWjIn4XsvoVRIMrkf3d13AYIYbWQypy7dRqkQAN9/lWLxPoYGv06x+NCUxSgixNLHoII38Vezsog6WawCpXzc3O8x4zthNnDIcjlRfBrXfYr29l+RTn+hYbFMBZEYSrmUSo/Q27MLhcJtjI7+HM97htbM6UD530LEQMQgDEcZHPgSnvsclr0RI8OnMzx8Fv+bfL62LGcXxJyDm7tqStpvVqvPMVQN+MX7UeEiYsnvNSyGMMzT338Ynvcs7R2/JpH4ZMNimUoiMVoyX8GyN6JYvI+YvTWJzDewlynIFIZZstmLMIw2WjInAZBOf4FS6W9AOEWxmcSSn6E0eg6B/19Mq3EDxuqpkgl754rIgyLyQlS+8MvR8u+KyEIReSb62n3cNqdF5QtfEpFPTOUO1JOXvw4xZmE5Ozekf6UCBgePw3Ofor3jVzM2UYyXSOxJe/vPSaU/97ZEMXatwHOfIjt6Eaa5Nv19h5DP3wyUb8FO5WmZnTwQMPHy109ZH82mktMQHzhZKbU5MA84XkQ2jz47Tym1dfR1J0D02SHAFsBuwK9kuo6PHScMevFLD2In92/YRa2R4TMpFe+ltfV7JBJ7NiSGRhm7Tep5zzMwcGz0eDzkcleSSh9JW/u5pNLHUCreFz2EFp/SeAxzFlb8o3j5G1EqmNK+mkUl5QsXKaWeil6PUi4DsLLBBXsD1yqlSkqpV4H5LGdi3+nGL94JBNiJ/RrSfz5/I7ncb0ilvkAqfWRDYmiksduktr0FrZnTEbEJgn4goLX1OwCY5hqAIggW1CUmO7kfKuwhcB+rS3+NtkoXOKOK6O8Fxn46J4jIsyJyuYiMjb9dWr4wMr604fi2plX5Qq9wG4a1Maa9Sf379l5meOhUYrH3k2k9ve79NxszGghnmp0YRhcjI+fie/MplR5GKQ/TXLMucVjxj4Ek8Qq31qW/Rqv4eFpE0sCNwElKqRERuQj4HuUiyd8Dfgp8vtL2lFKXAJdA+anTVQm63sKgj8D9B/GWk+red3m48/8hkqC941fvmHxm6vpVuP58iqUnKbr/xPVfxQsWEIbDhCoPGBiSxDTasa25xKz1ceLvJRHbDtuaU5cYAVoypzIy8gOGhk4lHt+JdPpoRKy6zIchRgIrvjN+8R6U+kHdx3zUW0XJQsq/oTcCv1dK3QSglFoy7vNLgdujt0vLF0bGlzaclvzifYDCasDUatnRC/G8f9LecQmmOXvK+yu6zzGSu4Fs4S68oHyAaEiGmL0+jr1FNN9mAhSEKkcQ9uH5b5Av/RWVvRSAmL0ZLYndyKQOJjbFj8Kb5iza239OEPRgmrOWLq/X6E7b2RW/eDuB99SMLyQ1YbKQ8k/9MuBFpdTPxi1fUym1KHq7L/Bc9PpW4BoR+RmwFrAR8HhNo64zv/QAYq6FYW0+8cq17Nd/ldHR80kk9p7SC5pKKbKFuxgYvZii+wRCjKSzEx2Z/yMR356YtcGEfzWV8ih5/yZf+jvZwp/oHzmf/pGfk3I+RkfmeJLxHaYsfuBtiaKeLOejgIlffFAnC+CDwOHAv0TkmWjZN4FDRWRryqchrwHHACilnheR64EXKN9JOV5N48vFSnn4pb9iJ/as67MISimGh05HJEam9btT1k++9Bi9Q2dRdJ/Gttaju+27tKYOwlzF+SZFbJzYu3Fi76aj5Wg8/y2Gc79jKPt73uzZl5Tzcbrbvk3c3nhqdgTI528iCBbQ0vJ/U9bHssRoxYy9F7/0EHBK3fpthErKFz4CLO9/yZ0r2eZs4Owq4moagfcsqFGs+E517bdUeohS6QEymTOm5PQjDPP0Dp/NUPa3WOaarNH+MzKpA2s2C5RtrUVX66l0tJzAYPZyBkYu5PXFu9LZ+jU6Wo6bktmm3NLfyOevJ5HYG8t6V83bXxErvhOl0fNQ4dCMntR3Zl+RqYGg9HcAzNi8uvWplGJ05EeY5hxS6c/VvP2S+yKvLdmVoexvaU9/gfXWeJjW9CFT8h/YMJJ0Zk5gvTUfIZXYhb7hH/Bmz374weKa99WS+RoiFqMj59a87ZUxY+8HFH5pWp9tT0gniwkE7j8wrI0wzM669Vkq3ofnPUtLy1drPrgoV/gzr/fsRRjmmNt9A7Paz6rLVHWW2cVanZewZseFFL3neX3JHpTcF2vah2muQTJ1JIXCTfj+KzVte6X9xrYG4jN+vIVOFiuhlMJ3n8KMbVvXfrPZizHNtUgk969puyP5W1nQdzgx6128a/adJJ0P1LT9iYgImdR+rDPrFkDxRs8+FEpP1rSPdPoYwI5qltSHSBzT3pLAe7pufTaCThYrEQavgRrGtN9btz5d91lc9++kUl+s6ZiKbOFeFvUfTyK2LXNn3Yxt1Wfg0vI4sS1YZ9atmGYnC3o/XdMjDNOcRTK5P/ncdYThYM3anbDf2DYE7r+WDkOfiXSyWInQK98NNmPvqVuf+fw1iDgkU4fUrM2i+yxv9R9L3N6cOd2/w2yCOTltaw5zu6/HMFIs6DusptcwUunPA0Xy+Ztq1uZETPvdQInQ/2/d+qw3nSxWIvBeACwMa8O69KdUgUL+Zhxnz5oV1gnCYRb2HYVptDOn++qmKQUA5YSxdtfVBOEwb/UdU7MHsmx7C2z7PeTzf6hJe5Uw7PIYnMB7vm591ptOFisRei9hWBtM+ROMY4rFP6PUKMlk7eql9gyegR8sZq3OS7EaNHBpZZzYFqzR/mMK7j8YGL24Zu0mkgfgey/geS/XrM2VMaz1AZswqhk7E+lksRKB/zKGVb8ZsYqFWzGMDmLx2lx4zBbuYyR/PZ2ZE0nE63fdZVW1JPclndid/uGfUKrRf+5EYg9AKBZun3DdWhCxMKz1CT2dLFY7Srmo4E2MOs2CpJRPsfggjrNrTebLUMqnZ+i7xKyN6IxmkGpWIsLs9h8h4tA7VJtZyExzTWx7G4rF+2rSXiUMa/3yRfEZSieLFQiDhUCIUaeaoJ77NEqNEI9/pCbtDeeuxfNfobvtm9NiIl/L7KIjczy54n3kS7UZrxB3PoznPUMYDNSkvYkY1rsI/Tdn7GQ4OlmsgPKjJy7N+iSLkvsoAPH4jlW3pVTIwOhFOLFtSDXgSdnJak8fhWl0MTDyq5q0F49/CFCU6jRYqvy74qLCnrr0V286WaxAGLwFgGGuVZf+XPcJLGsDDLOj6rbyxYfw/FdpTx9V04ffQuVS9F5hqPAgw8VHKNX4r6hhJGlLH0aueB+u/0bV7cViWwE2nvtU9cFVQKLflbHfnZlGz+69Aios3/eXOswhAeC5/yQe/3BN2hrOXYtpdNKS3KMm7RW9V1gyehV9uRsIVfZtn1lGF7PSn6a75TBiNbjb0po6jP6RCxjJXUdXa3VPcYo42PbmuN4/q46rEmNFp9QUPPfSDHSyWIEw6Adprcv5fhgMEIY92PZm1belimSL95NJ7l917Er5LBw+n0Ujv0Qw6UjuScbZkbj9rvL8Ff5rDObv5a2RX7B49BLWaT+DrtQhVR3N2NZaJOLvI1u4u+pkAWDbm9XtIqcYXQCosL8u/dWbThYroMJ+DKM+D495fvl2oVWDuR4KxUdRKk9LYreq2gmVy3/7TmCocA9dqQOZ03Yqttm9zFrvpzt9KEXvNV4fPJ3XBk6j6L3KnLbTqkoY6cRu9A6diecvqHqKPsvahDC8ljAcxDDaJ96gCmKUTyFnarLQ1yxWQKkRpEajKCcyNhu1VYOLqQX3H4BBIl7dhOoLhs5hqHAP67R/l/U6f7KcRPE/jr0uG3dfyaz0Z1k8egl9uWur6jsZfz8ABbf6h8zM6G6W7785wZrVE7FA0qhwZMr7agSdLFYkHIU6PUMxlizGZq2uRsF9mri9KYaRmnQb2dKTLBm9jFnpzzG75XMVbSNisE77d8k4O/LG4Jm4weTvCMTtzRBxKJaqvzBpmuWfab3KA4iR0clidaNUHpGpn+cBytdHRFoQSVTdlue9QqzKOqyLRy/HlAxz2k6deOVxRAze1f49QlWkd/T3k+5fxCZmrY/rvzrpNsaY0RFRWKdTA5EkShXq0le9VVO+sENE7hWRl6Pv7dFyEZELovKFz4rINlO9E1NBqSKIU5e+anU+rVSAF7yFXcXpjBf0MJj/E93pQzAnMSmOY69Hq/NRerPXVHVb1bbWwavB7dOxn2vdHlcXB1bXZMGKyxd+A7hfKbURcH/0HuCTlGf03gg4Grio5lHXhYdQrxodeUQmf9owJlSjgI9pdk26jZz7PBDQlthl0m20J3fFC3txqxhvYBqdBOHQpLcfI+IAZt3+2gs2iqmp3t5o1ZQv3Bu4MlrtSmCf6PXewFWq7FGgTUQaN9PKZKkA6lSiVSmvJhPdhGEeAKOK06eS/zpQvmg5WfFostxSFUcGhqRQKj/p7ccTsVHKrUlbEzPKvzszUDXlC2ePqxuyGBgbvTRDyhfWs0harfoKo++Tv22plB+1MPnkNbZtWNV/UEEt3Z9qCdSsrUo0dYG9Sas4WSxbvnD8Z6pc4nqVfkJKqUuUUtsppbbr7l7xbbmGEQtUfX7Byn/5qj90HTuiqOYvcjy61ViMjjAmY+nRibXupNsIVb6qI6TxlPIQ6jMnCYTl350ZqKJksbzyhcCSsdOL6PvYvbIZUr7QAuozn6JIoiaH3IbRAghBFRfzkrHyKNKR4iOTbmO4+AiGpIlVMaAqDIcwjMyktx9TPv3w6zaBkcJDqM/pa71VcjdkueULKZcpPCJ6fQRwy7jln43uiswDhsedrkwbIvHyHZE6MIw2wppczLOxzDXw/MmPKYhbc8k4O9KT/d2kjnbKd1PuoCt1AEYV12E8/w1sa+7EK05g7Oc61aM3l6rjXbR6q+TIYqx84c4i8kz0tTvwI2AXEXkZ+Hj0HsqVyl4B5gOXAl+qfdh1IIm6XUE3jA6UGq7JRTjbWg/Xn19VG7NbjsQLFrF49DervO2CoR+jCJjd8tlJ969UiOu/Qsxcd9JtjBkbX2EY1T/NWwmlijUZL9OMqilfCPCx5ayvgOOrjKvhRFrKozjroDzKUBEEi6ouu+fEtmJo9DJCVcKY5KF3q7Mz7YlPsnDoXBL2JrQlPlrRdktGr6AvdwNrZf4Px15/Un0DuP5/CdUoTnzrSbcxJvDHRsdW94xJpVQ4gjTB7OlTQY/gXAExMig1XJe+zOhwO6jBIKREbDsULkX3mUm3ISKs2/kjEvamzO/9IotGLiZcST2MIMzz+sAZvDH4XdoSu7JWa3WFiQvRTFlOrPrxfEFQvjFnmlOfLJQKQY2C1OeZonrTyWIFxOhE1WnUnxXN8+n71U9Wm3Q+CNhkC3dXF5PRyiazf09bYhcWDP2IZ9/6EG8NX0i29Ax+MIgX9DBa/AdvDJ7NPxfOoyd7JbNbjmKDrgurnkM0W7gb21yHWA1KMHj+y4hkMIypv+OmwiFA1e2Up95m5j2eGigniwGU8msyge7KGMZsRNrwvH9X3ZZpZEg5O5LN30F367cRmfzfA8toZYOuXzFcfIAlo1ewcPhcFg6/veiwYNGe3I3ZLZ8nHa/BkUA4RL74MG3pz9Vkli/fexHL3rimM4atiArL44WkjnVx60knixUwzNmAQoW9iDm1A1BFBDv2btwa1crMJA9g0cDxFEqPk3Sqq/4uIrQlPkZb4mOU/DfJu/+m5L+OiE3cWodUbAvsGtYjGc5dh8Ilk6q+dopSPp73L5LJT9cgsgr6C5YAIEZ9ZlerN50sVmD8fIrGFCcLgFhsW7KjFxCGuaoeLwdoSe6Bba1b81ohcWsu8RrczlwRpXyGsleSiL0PJ7Zl1e153gsoVcCuwbWPStR73tZ609csVsCI5kFQVYxZWBWx2A5AiBvN8l0NkVhTFxVakZHcH/H812hvObYm7bmlvwIQj+9Qk/YmEgYLABOJ5uKcaXSyWAEjuoUZBpMf9rwq4vEdEHEoFf9c87b9oBmfvXm7MMzTN3IuTmwb0lVOCTimWPozlrUJZp3+0ofB64i5Vk0eCmxGOlmsgIiDGGvWrSq2iEMsviPF4t2Uh6pULwwL9A79kCUDp/D64k+SL1Z/1DJVeod/iB8sorvt9JpcjAzDIdzSo8SdnWsQXYV9+q9gWOvVrb9608liJQx7I4IqR0OuioSzJ0GwAK8GFzqVUozkb8b1XmKNzgvoaj2VodyV+EHzTSabL/6NoexltKU/T7JGpwzFwt2ARyLxqZq0NxGlQkJ/PmYda+PWm04WK2FaGxP6L9etHJ2T+AQQp5C/seq2RIRc8UGSzocxjQzx2OYoVUIkRr70OD1DZzXFkYbnL+Kt/i9hW+vT3frNmrWbL9yIaa6DbW9VszZXRgVvgipgWNXP0N6sdLJYCcPeHFSR0H+lPv0ZrSQSu5HP31yTh9hi1gZLx4iM5G8hGf8ASwa/zsDILzCNTnqHz6Z3+Jyq+5msIBxlYd/hhCrH2p2XYkxiGr/l8f3XcUuPkEweXJfxFQCB9zwApr15XfprBJ0sVsK0y7fvQu+5uvWZTH4apYYo5G+ZeOUJtCT3YCR/Cwv7vsho/nZc/3VyhQeZ1XYmnZnjWaPjXCT6FajVdZJKhWGWhb2HU/JeYu3OS4nHqi+wNCafuxowSCYPqlmbEwm85wATw96kbn3Wm04WK2FYG4IkCWo0WKoSsfiOWNZmZLO/nvR/4CAcwQ8W48TezTqzbiDp7Ejc3hwn9m4yyX0wo3qqgknR/RdBOFK3v8AAftDDGz37U3CfZM3OX5JKfKRmbYfhKLnc1TiJPWpSWqFSgfs0RlTCYKbSyWIlRCxMeyv8OhXWLfcppFuOxff/TWmSZfdyxQd4o+cA+kd+wXDuBkbzt5F05pFJHoAYDgv7vsBI/hbe6j8Wy+zGNDLvSExTdZ2mUHqa15fsjuvPZ+2uK8gk96pp+/nc1Sg1Sjp9XE3bXRmlAgLvn5j29Bvbsip0spiAGXsfofc8qk6PqwMkEvtgmu9iZPTH5ScZV1EmuQ9zuq6g5P0bP1jMrLbTyST3QcRgVtsZdLQci+8vwvPfpDPztWircrL4X38GftBH79DZeH71cxcp5dI7fA5v9OwFGKwz62bSiXfMcFCVMBwlm/0l8fhHiMW2rmnbK+3Xex5UFitWXRW4ZqeTxQSs+DwgwHf/Ubc+RWxaMl/D956nUPh/k2ojZm/IWp2/pDNzAnH73QyOXsHA6CXlz6z18PzXyaQOwLbWRKlg6QNnY98HRi+kf+Q8hrJXM5K/YdL7opRiNH8Xry3ehYGR88kk92fdNe7Fib1n0m2uSDZ7EWE4SEvmGxOvXEN+NOrWrNNI0UbRz4ZMwIxtB8QJSg9j13GATyKxD7nsJYwMfx/H2RXDSE+6LRGDlPNhlgyexkjuJmxrLinnw2SS+wD/mzdSKY9c8WGyhT8hmKQTu+L5C8gk9y2vp8KKn2JVymU0fzuD2csouk8TszZg7a6rSCc+Pun9WBnff43s6EUkEvsSm4JEtNK+S3/BsDbEmKHDvMfoZDEBEQczvgN+6aE692vS2vZD+nr3ZHT0p7S2nlFVezF7PebOuhbXexURB9v638NxRnRRbnD0N5S8F+hqPQ3L7KZ/5Dzi9ibY1hyUUhMmijDMki89TrbwJ7KFuwjCfmxrfWa3/4TW1MFT9qi/UorhoW8jYpNpPX1K+lhx30WC0mPEUofVtd9GmPBfT0QuB/YEepRSW0bLvgt8ERh76OCbSqk7o89OA44CAuD/lFLVzcLSBKz4RyiNnEXov4ZRxfT2qyoW24Zk8nBy2UtIOLsTi7+v+jbt/w1HVsqnZ+gM2tNHYlvrYRgZcsVHMI2LScS3p+T9h9bUIby+ZA9sc31i9jqYRvvSOSZDlSMI+vH813H9+ZS8fwMhIinSzs60pg4l6exU1ZwalSjkr6NUeoBM61mYdf7r7pceBkpY8cqmHpzOKkn1VwAXAlcts/w8pdTbZkKJyhoeAmwBrAXcJyIbq3oNgZwitrMrpZGz8Ir3Ek9/sa59Z1pPp1R6iMHBL9M9696qH18fT8SiNXUwMbs8I1Vb+jO0JPeid+hMFg18BSe2JYakAYeC+3dGCzfxzvIwFra5Nra9Hp2J3UjEtiPhzFt6tDLVfH8Bw8PfIRb7AKnU5+vS59v6L94L0oIZr27ekOmgkgl7/xJVIqvE3sC1SqkS8KqIzAe2B/4++RAbz7DWwbA2xy/cWfdkYRhp2tp/Tn/fAQwPnUpb+4U1HRMxdqGxfBdEMI0WWlOfwQsW4thbYllr8q7ZN0breIRhllDlAQPDSGJIGqlTmcdlKVVicOBYQGhrP2/Kj2De2b+HX7wby/koIrG69t0I1fx0T4iqpF8+VkGdCksXwnQoX/h2dmJPAu9JQr/+9ZLi8Xm0ZE6hULiZXO6yKelDxEBECFWRQunvJGLb0t32LWLjZhsXsTHNdmxrbWxrTUyjtWGJAmB46Dt43lO0tf8Mawon5VmRoPQ3VDiI7dR2rEizmmyyuAjYANgaWAT8dFUbaPryhcuwo6cXvcLNDek/nT4Rx/kEI8PfpVi8f8r6McShI/MlOjInAExqnEc9ZLOXks9fTTp9PInEHg2JwSvcDJLBcnZqSP/1NqlkoZRaopQKVPk36VLKpxowY0oXvpNhvQsztgNe/o91f44Cyn/529ovxLa3YHDgaFz3ySntb+yaQ70P7StRyN/CyPB3cZzdacmc1pAYVDiCV7wTO7HXjB7iPd6kfhPGapxG9gXGnrS6FThEROIish6wEfB4dSE2Dzt5EGHwKoHbmF0yjBQdnVdjGLPo7zsM1322IXE0UqFwF4ODJxKLbU97xy8alsy8wq2gisTq+LBao1VS6/QPlC9QbiIiC0TkKODHIvIvEXkW+CjwFQCl1PPA9cALwJ+A46f7nZDxbGdPkFbc3LI3hurHNGfR2XU9htFCf9/BuHV8bqXRCoXbGBw4Btt+Dx2dVzasTKBSCjd3FYa1BUad5stoBpXcDTl0OYtXeJVNKXU2cHY1QTUrMRLEkgfh5n5LGCxu2Ig9y5pLZ9eN9PcdRH/fgbR3/BrHmZqRkc0il/0tw8Pfxo5tS2fn72pSYX2yAvcxQv/fOK3n1PVp3UZrvhPSJhdLHQGEuNmpuStRKcuaS1f3LVjWhgz0f45s9jcNuZYy1ZTyGB46neHhbxF3Pk5n57UNTRQAbvZixOjAjobBry50slhFhrUOlrMHbv4aVFifWqgrUj4luQnH2YWR4e8wNHgCYZhraEy1FARL6O87mFzuMlKpL9LRcXnNZtOadEzeS/ilB4ilPrfaXNgco5PFJMRbjgM1ipv7baNDwTBStHdcRkvm6xQKt9Dbswul0vS/plwo3EZvz8543j9pa7+Q1rYzGzqmY0xp9AKQFHbqiEaHUnc6WUyCaW+J5exKKfubqBhuY4kYtLR8mc6uG4GQ/r59GR4+kzDMNjq0VRYESxgYOJbBgWMwzXfR1f0nksn9Gh0WEB1VFG8nlvochtE+8QYzjE4WkxRvORnUKKXRXzY6lKXi8R3onnUfydTh5LK/pmfJTuTzNzXtwKrxlCqRHf01PUs+RLFwNy0tp9LVfSu23TxT65dGfgjSQix9dKNDaQidLCbJtDfDTuyPm7uC0H9z4g3qxDDStLX9iK7u2zHNWQwNnkBv7ycoFu5tygugSnnkc9fSs+RDjIycSSy2PbNmP0BL5qQpr16/KvzS3/BLDxBPH79aHlWAThZViWdOARGKI813pzgW24au7jtoa78QFWYZGDiC3p6Pkc9dV5MyA9UKw2GyoxezZMk8hoa+imF00tH5Bzo6r8ZqsqpeSvkUh89EzDnE0kc2OpyG0cmiCoa5JvH0CfjFO/GL9Z0cpxIiJsnkfsya/Rfa2s8HhKGhr7B40TYMD30L132qrkcbSvkUiw8xOHAiixe9l5GRs7DM9ejovIqu7jtxnA835biF8tHjiziZ01e7OyDjSTMcmm633XbqiSeeaHQYk6JUiVzPJ1B4pLvvQWo430StKaVw3b+Sz11DoXAXUMI05+A4uxJ3PkIs9v6azpcBEAT9uKVHKJb+TKl4L2E4gEiGRHJfkslD6z4F3qoK/TfJ9u6CFduBRMcVTZnMKiUiTyqltpvs9s1zUjhNicRx2s4h338gpdFzcFrPanRIKyQixOM7Eo/vSGs4TLF4D4XCbeTz15DLXQ4YWPamxOz3YtmbYlsbY1pzMc21JpyvQakCvr+AIHgT33sJz3sR13uKIKrmJtKG43wEJ/EpHOcjDRuqvSqUUhSGvg4ITuvZ0zpR1IJOFjVgxXcglvocbu4KLGcXrPiHGh3ShAyjlWTyQJLJA1GqiFt6nJL7KJ77FIXCHaj879+2vkgrRjSlnogNKJTyUCpHGA6iVHaZ9mcTi21FMnkw8fgHsO2tm2KcxKrw8lcRuI/gtJ6NYc1pdDgNp5NFjcRbTsMvPUJh8CRS3fdgmJ2NDqliIg5xZyfi0bwMSinCsAffn0/gLyAIFhCGg1FSKKKUCyIIMUSSGEY7htmBac7BNOdi2xtN+zsGgfcixeHvY8V3xk7O/Ml4K6GTRY2IkSDRfiG53r0oDn2VRMdvm3IuiEqICKY5G9OcDfFGR1N/KsxRGDweMTI4beeu9qcfY6bnb3OTMu3NcVrPwC89SGn0vEaHo01C+TrFyYT+KyTaz8cwuxodUtPQyaLG7ORh2ImDcLPn4xXuaHQ42ipys7/AL95JPHMaVnzHRofTVHSyqDERwWn7Pqa9LYXBk+paVFmrjpe/hdLoudiJfYmlVs8h3Sujk8UUEHFIdFyKmGtQGDiSwJvf6JC0CfilhykMfRUztgNO24/1dYrl0MliihhmF6nOqwGTfP+hhP7rjQ5JWwHffYL8wBcwrA1IdvwGkdXwqm4FKpmD83IR6RGR58Yt6xCRe0Xk5eh7e7RcROQCEZkf1RTZZiqDb3aGtS7JzmtAlcj1H6ITRhPy3SfI938Ww1iDZOfvEaO10SE1rUqOLK4Adltm2TeA+5VSGwH3R+8BPkl5Ru+NgKMp1xdZrZn2piS7roEwR67vQH1K0kT80l/J9x+GYXST7LoOw2z++jWNNGGyUEr9BRhYZvHewJXR6yuBfcYtv0qVPQq0LVM2YLVk2luS7LoeCMj3768vejYBr3A7+f4jMMw5JLv+2LDJl6eTyV6zmK2UWhS9XgzMjl7P2PKF1TLtTUl13QiSId93sL6t2iBKKUrZiykMHo8Zew+prhswzFmNDmtaqPoCpyo/trrKj65Ot/KFtWBY65LquhnT3pLC4HEUR346LWaxmimUKlIcOonSyA+wnN1Jdl6DGG2NDmvamGyyWDJ2ehF974mWz9jyhbVimF0ku/6AnTwYN3s+hYEjCYP+Roc14wX+K+R698Ur/D/iLaeQaP/laj03xWRMNlncCoxNb3wEcMu45Z+N7orMA4bHna5oEREHp/XHOK1n45f+Rq73E/ilhxsd1oyklMLN/5Fc7+6oYCGJjsuJt5w4bZ/baaTJli/8EbCLiLwMfDx6D3An8Aown3LB5C9NSdQzgIgQSx1OqvtWxMiQ7/8MxeHvoMJ8o0ObMcKgh8LgsRSHTsa0tyI1625s52ONDmva0jNlNQEVFiiN/hA3dwViziXR+gMs58ONDmvaUirEK9xIafh7KFUg3vJVYumjp918GrVW7UxZ+lisCYiRwGk9i2TnHxFs8gOHkx84pqlmDZ8uAu858v0HUhw6GcPagFT3XcRbjlvtE0Ut6GTRRKz4DqRm3U285RT80oNkez5Kcfjspihk1OxCfyGFwZPI9e5B6L+C0/YTkl03YtobNjq0GUNPftNkROLEW07ETh5AaeRc3NwluPk/EE8fRSz1eT0ceRlhsIhS9iK83DWAEEsfSzz9Jf1zmgI6WTQpw1yTRPtPiaW/QGn0Z5RGz6OUvZRY8jPE0ket9iMOA+8/uNlL8Ao3Awo7eQDx9JcxrOWOAdRqQCeLJmfam5HsuJTAe45S9mLc3KW4ucuwnN2IpY7AjG2/2jxOrZSPX3wAN38lQelhEAc7+Wni6aMxrLkTN6BVRSeLacK0tyTZfiFhy9dxc1fi5q/DL96OYa6HnTwQO7HPjJ2BOvBewivciJe/ERX2IsaaxFtOwU5+BsPsaHR4qw1963SaUmEBr3g7Xv6PBO6jAJj2NliJ3bGdXTGsdRsbYBWUUoT+83iFe/CLdxD6LwMmlrMzduJgLGfnpqqDOl1Ue+tUJ4sZIPRfxyvcjle4jdB/AQDD2hArvhNmfCes2PsQo6WitkrZy/Dyf6B8HeBQ4ukvEHgvUBz6JkrlMMw5JNovWG57pexvom0Fw96URNu5FQ+pDoM+AvdR/NJf8IsPocJFgGDG3oed2AvL2V1PnlslnSy0twn9N/CK9+KXHiQoPQaUAAPD3gIzti2mvTWmvSWGtcE7xh4E3ksUBo8n1XUbiE2+/3ASbT8kP3gCTubbWPF5uPnrCP03cTJfe3u/wWJyffuTnnU/Ig75geOwnJ2JJQ98R4xKlQj9lwncfxF4TxO4TxD60TwfksGKfxAr/lEs5+M6QdSQLl+ovY1hrUM8fRTx9FEoVSRwn8QvPUrgPo6Xvx5PXVFeURwMa0NMayMMa30M812ocAgz9l7EKJcWtOLz8Ap3EfqvYsZ2iJZ9iHz2MFgmWQCgfFBFFBaoAoYxm8B7kcB7jtB/ndD/L6H/MqH/CuCXw5A2zNh7sRP7Y8bnYdpb6VOMJqX/VWYwESf6K/1BAJQKyn/RvecIvRcJ/P/gu4+hCjcDYKeOJig9ThgOIjj4xQcxYu/BtDbGL96DnfgEXuEOwuCdzwYa5hrE0kczumQeIk759MfZicLgyXiFPwIGYs7FtDfGcnbBtDfHsLfEMNddbe7mTHc6WaxGRExMe1NMe9O3LVeqSOi/gRgZfGtD8v2fKZcltDdHMHDafkJx5AxK2fOxnV0Q7He0rcIh/OK9pGf9FTEyFAaPw83fRLzlxPJTnubaUY1UbbrSyUIrHwnYGwMQSx1CLHUIAMWRczDMNTHtDUl1lgslB/4reMUH3tGGX3oEw5q7tMar5exG4D5JLLlfnfZCm2r62RDtbcKgr/zdX4hf/BN2Yu+ly5QKcUcvIJZ6Z6FgMdcmcJ9ChQWUUgSlv2Ja+rmMmUQfWWhvUxg8BhUOAjZO6/cQoxU3exle7ioArMRu2ImDgPIdkOLQ10l2XokVey+Wszu5vt0BE9PeAjv16cbtiFZz+tappq0m9HwWmqbVhU4WmqZVpKprFiLyGjAKBICvlNpORDqA64B1gdeAg5RSg9WFqWlao9XiyOKjSqmtx50Lrai0oaZp09hUnIasqLShpmnTWLXJQgH3iMiTInJ0tGxFpQ01TZvGqh1nsaNSaqGIzALuFZF/j/9QKaVEZLn3ZqPkcjTAOuusU2UYmqZNtaqOLJRSC6PvPcDNwPasuLThstuudrVONW06m3SyEJGUiLSMvQZ2BZ5jxaUNNU2bxqo5DZkN3Bw9XmwB1yil/iQi/wCuj8ocvg4cVH2YmqY12qSThVLqFWCr5SzvB3RBSU2bYfQITk3TKqKThaZpFdHJQtO0iuhkoWlaRXSy0DStIjpZaJpWEZ0sNE2riE4WmqZVRCcLTdMqopOFpmkV0clC07SK6GShaVpFdLLQNK0iOllomlYRnSw0TauIThaaplVEJwtN0yqik4WmaRXRyULTtIpMWbIQkd1E5CURmS8iuoShpk1zU5IsRMQEfgl8EtgcOFRENp+KvjRNq4+pOrLYHpivlHpFKeUC11Kugapp2jRVbfnCFVkbeHPc+wXADuNXGF++ECiJyHNTFEujdQF9jQ5iCuj9mn42qWbjqUoWE1JKXQJcAiAiTyiltmtULFNppu6b3q/pR0SeqGb7qToNWQjMHfd+TrRM07RpaqqSxT+AjURkPRGJAYdQroGqado0NSWnIUopX0ROAO4GTOBypdTzK9nkkqmIo0nM1H3T+zX9VLVvopSqVSCaps1gegSnpmkV0clC07SKNDxZTOdh4SJyuYj0jB8jIiIdInKviLwcfW+PlouIXBDt57Misk3jIl85EZkrIg+KyAsi8ryIfDlaPhP2zRGRx0Xkn9G+nRktX09EHov24browjwiEo/ez48+X7ehOzABETFF5GkRuT16X7P9amiymAHDwq8Adltm2TeA+5VSGwH3R++hvI8bRV9HAxfVKcbJ8IGTlVKbA/OA46N/l5mwbyVgZ6XUVsDWwG4iMg84BzhPKbUhMAgcFa1/FDAYLT8vWq+ZfRl4cdz72u2XUqphX8D7gbvHvT8NOK2RMU1iH9YFnhv3/iVgzej1msBL0etfA4cub71m/wJuAXaZafsGJIGnKI8u7gOsaPnS30vKd/TeH722ovWk0bGvYH/mUE7iOwO3A1LL/Wr0acjyhoWv3aBYamW2UmpR9HoxMDt6PS33NTo8fS/wGDNk36JD9WeAHuBe4L/AkFLKj1YZH//SfYs+HwY66xpw5X4OnAqE0ftOarhfjU4WM5oqp+1pe29aRNLAjcBJSqmR8Z9N531TSgVKqa0p/yXeHti0sRFVT0T2BHqUUk9OVR+NThYzcVj4EhFZEyD63hMtn1b7KiI25UTxe6XUTdHiGbFvY5RSQ8CDlA/P20RkbJDi+PiX7lv0eSvQX99IK/JBYC8ReY3yU947A+dTw/1qdLKYicPCbwWOiF4fQfl8f2z5Z6M7B/OA4XGH9E1FRAS4DHhRKfWzcR/NhH3rFpG26HWC8rWYFyknjQOi1Zbdt7F9PgB4IDqqaipKqdOUUnOUUutS/n/0gFLqM9Ryv5rgoszuwH8onzd+q9HxrGLsfwAWAR7l88GjKJ/33Q+8DNwHdETrCuU7P/8F/gVs1+j4V7JfO1I+xXgWeCb62n2G7Nt7gKejfXsO+E60fH3gcWA+8EcgHi13ovfzo8/Xb/Q+VLCPHwFur/V+6eHemqZVpNGnIZqmTRM6WWiaVhGdLDRNq4hOFpqmVUQnC03TKqKThaZpFdHJQtO0ivx/iFurfyR+Ux4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = y = np.linspace(0, L, N)\n",
    "c = plt.contour(x, y, h[-1], np.arange(90, 100.1, 0.2))\n",
    "plt.clabel(c, fmt=\"%1.1f\")\n",
    "plt.axis(\"scaled\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:50.798444Z",
     "iopub.status.busy": "2022-07-30T20:16:50.797246Z",
     "iopub.status.idle": "2022-07-30T20:16:50.960532Z",
     "shell.execute_reply": "2022-07-30T20:16:50.959638Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAABFCAYAAABHYOExAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAc6UlEQVR4nO2deZgdVZ33P6e2u2+9p9OddBKyse+gAvoCsskIPuO8w/M6wCiKivq6MLIIIoqAjA/LvKPoMICI4Cg648YoDgiOMEpIgBACZE+g00m603v33epW3fP+cas7TdI36dt9l2q6Ps/TT99bVefUt373nvur3zmnzk9IKfHw8PDwmLsotRbg4eHh4VFbPEfg4eHhMcfxHIGHh4fHHMdzBB4eHh5zHM8ReHh4eMxxPEfg4eHhMcfRai1gIg0NDbKjo6PWMiqKtHcjZQpFW1JrKUWx7d1ImUTTDqtI/VJmMa3t6Op8FCUy5XKmtYu8zODXF09WK0lzC7pah6HWA2DlR0laXUT0RUhg2HyLkN5CyhpEIMjks6jCh5k3saV4259lK9h5BWyByIGaBVsXSJ8ERaKoEk3Joyp5VCHH/zShAjkCaggrnyZqtDJsvklYb8VQwgybm/Cp9fgcjXY+SdraSVBfiCL8JV5zMfvmMK0t6GorihKbcrlSsO0upDTRtEUVqb8c5O0ukCaKizWWixdffLFXStk43fKucgQdHR2sWbOm1jIqSnrgi1jmC0Sa/6fWUooyMPBFzOyzNLdU5rMwrR1s3/1uWuq+SSz0N1Mut6P/RvpTj3N82+S6nus8jfrA6axsuB2AEXMHT++8mGMabmRB5GJ+tu0C2kLvwSTIluFnyMo6FMVPV7qflK2Tsg2Gc35SlsFA1k/vSIh0RkeO6qhdKnZzHhGxMAI56iJJ6v1pQrpJVM8SVE0iao7WQCPDuU4WBVoQQnBU7CzW9N7DRQsfw5Yj/KnrI5zQdBtt4QsA2Dn8CJv6b+E9bX/ApzUfcE0be/4OOz/K4S2/nLKdsuYGdnSfSWv9XUSCfzXlcqXQ3/cxLOstmpqfqkj95SDV/0ny1jbCTU/WWkrFEUK8OZPyXteQxwEIFCr5oKEiwgDk86MllTPUedj5Qex8atL9QX0xydy28fdhfSE+tZHe9BoUodEePoPO5LMcFnkfOZlmSeRIhnK7CKg+wqqKT8kT1Gx8Wo6okSUWShMOZtHjJvbCHGqdSTBiEg1mSPiyBDQLv2phKBZ+RSKxSVp7OCx8PH3mVpZF38+bo0+TMJYSNRbQm14NQIP/pHGNydxWNBHBUJsmvSbT2oWhtpZkp7ws2HXMzpVAYoMQFavfo7p4jsBjEhTArlztIgTs+8GaKj5tPgCmvWvS/SF9CUlzy7gTE0LQGDiJ3sxqpJQsDJ9JLj8KMkNYayJn92EoAZp8deTJoAoLQzHxqxZ+zSLhyxIPZoiG0oTrMkTDGWLBFHXBFH4th0+18Cs5NGHjV2B+oBVLZqnTYwgU2oJHszezjgWR/wVAb3o1YX0Rfm1fBJ/MbSFoLEFM8qMqpcS0d2Fo03MEQgmWVK4kZB6BWrn6PaqK5whqgsuX9RAqlXQEQvgBnXx+qKRyhtYGQNaaPAoOGyux5SjJ3KbxbY2Bd5G1++jNrKYleBI+JcYbg//G8ug57EytYXn4ePrM7QRVjYRuYCgWIdUkqGUJ6lli/gwJf4ZoIEPcnybhSxPWcgQ1k5BacBpxTUfKLKNWF+3BlWwbfYb20EnsGPk9IOgIn03G6qU3s5qmwKnj2ux8hpHsa4T15ZNeTy6/l7zM4HOue6rY+WEAVBEtqVxp2OB6R+DyduYiPEdQddwfTgs0pLQqV78QaGojlr23pHJBfSUgSJrrJ93fGDwbgcqe0d+Mb5sfOgef2sCmwQdQhc6xDZ+iO/0ScT1B3GinN/MybYGVaCJNXiap03V8So6wahLSTCKaSdRXcAJxX4awXtge1kwCqklc07DlCFFdIazG8Is0EsmRsfPZMPhTDot+kKixgK1DPyIvLRbFLhnX1pt+GlsmaQqdP+n1pMxXnes+oiQ72XYPAGqR7qZyILEQwlVDjEVwf3tzA54j8DgQoQO5ip5CU5ux7O6SyqhKCL+2ZPwHcn8MtY66wOl0J3+DlHmnjI/DYpfRm15Ff+YVlkYvojlwHC/33cvpjZ8hmx/BEIMkjBZCmkRXJGHVJqDmiGhZonqaqJ6mzkgRNdLE9AxhLUNQMYmoEkmGRl8IVSi0B+bRk3mD05s+y6v99xHQ6jmh4XNk7QG2Dz9GW/hcwvrCcb3do7/GpzaT8J886fUks+sAhaBxeEl2suweQENVEiWVKwmZA/TK1e9RVTxH4HEAlY4IADSlESvfU3K5kHE0SXNd0f0toQ+StfcwmFk9vq0j+jcYSoKNA/chhMKpTV8hL3NsGf4Z57V+g4HsW8Q1Bb8SRBXDNPsS+JUcETVHRMsQ19PEjRRxLUVEyxBRnWhADxDXIZsf4bDQEjpTq3hX46fIWjsZNLdySuO1GGqYbUOPYssMS+MfH9dk2v30pZ+lOXQhQkzexZI0X8WvL0FVQiXZyLJ70NRGhKhc85ZytkQEHlPBcwQ1weV9l05EUMmZQ6rajF1iRAAQ8h1Nzu4ha3VOur8heCaqCNE18uPxbZoSYEns7+hJ/w+7k88QNdo5tv6T7Ew+x7C5gXNab6Ivu4UmX4Bl4ZMYtbto9PmJ6T58ikVIzRHTcoS1XMFB6DrNvgim7KHBWMCK8Ao6U6s4qf7vaTTm82r/Q3REzqEtfBoj5na2Dj1Ka+gsosa+Z0d2j/4ciUVL+IOTXoeUNknzZULGUSXbyMp3o6kHTkUtL7MhInB5O3MRniOoOu7vsxQYFBpR5aICXWvDzg+QzydLKhfzF2bgDKR+P+l+VQnQFr2UntQTjGRfH9++JH4pMWMlL++9iWSui5XxS1gUOY+1fd9jILuWc1tvZsjsoi+zipMTZ2EIH5J+GnwGS8OLWRRayLLwYczzR1DFIJZMcWLibJA72ZlazakNVxLT4vz3nq+Q8C3llMZrsfJpVnd/GVX4ObL+y+NacvYwbw49QF3gdMLGikmvYyS7Gis/QDxwVkn2AchZnWhqaQPMpSKliRBGRc9RHtzf3tyA5whqgrvvVMYauJRmxc6hawsAyFlvlVTOry8koK9kIPW7oscsiF6BrsTZPHD7eFSjCoOTmr+NRLK6+x+wZIp3N9/IstiHeH3gUTYNPsyF829lXuBoNg4/TlhNc2LiTBYEljOY62Rv9g32ZjfRYLRyYuJs5vsTbB7+DQE1xgfm304y9war9t5Ok/8Yzpp/D7riZ+3erzOS28oJTbcSmPCw2I6he7HyQyyJX130GgZSTyCEj5j/fSXZR8o8lrUTXWsvqVypzA5H4O525iZm5AiEEN8WQmwQQqwTQvxCCBGfsO96IcQWIcRGIcS5M1b6jmEW3KGMNXCZrdgpdNVxBHZpjgCgLng+o+aLmPbkYwy6GmVR/P8ymHmB3tQfxreH9DZOaLqNYXMzf9l9FbbMcErTtZzecguD5lae3XMdyyKncEHrrfjVGJuH/5Pu9LP4xQBtvjhRNUt/dhWbh39Dxh7izOZrOKXhI6ztvYs3R57imPorOWv+PfiUCC/vvZmu5BOsTHyWpuC7xzWkcjvYOfwI88J/TcS3clL9UuYZSD9BzP/eaYwPdCMxxx1t5ciC8FX4HOVgFrQ3FzDT0Z4ngeullJYQ4g7geuBaIcThwCXAEUAr8JQQYpmUsnKT02cTLk8PKpwGLql8RGCWGBEAJALn0TV0F/3Jx2mJfmzSY1ojf0vXyI/ZPHAbicApaM6aRi3BMzix+R9Z030Nf959Jae2fIeOyPup8y3nuT1f4/me29GVMMtjH+a0pqsYtfbSl93OoNlJWGui3reYuNHGUPYNXh/4F5LWHoJaE2fP/2dagidgS5OXem6iK/kEKxKfZlniinFNUubZ1HcLitBZHP9C0esbNV8iZ+8hEbymZNuMOVZDrawjkNJ0uhDdjLvbmZuYUUQgpfwvuW96yfPAWMfkRcBPpJRZKeV2YAsw+Ry5ucYseCxf4DiCCkYEqlKHoiQwc1tKLhswlhEyjmHv6CPj00T3RxEaK+pvIWvtYUPvV9828N0aOpOTm+9ixNzGs12XM5BZT9RYwPntD3Ju2/eZFziR9QMP8YeuT/Fa//fJ5DbTZDQh891sGXqEp7uuYk3vPQS1Rt477w4+1PEftARPIJXbxfO7r6Ir+QSH132e5YlPvk3TW8P30595jsMS1+DTiq8P1jPyMKqIkAiUHkiP2dPQK7Ng4BhSZsdvGNyN+9ubGyjn/K+PAT91Xs+n4BjG2Ols8wDcfqcyHhFU0BEIIfBpSzFzG6dVvilyOdv7vsRw5jligTMmPSbmP57FiS+wdeBOYiPH0x69bHxfS+gM3j3vX3ih+2r+tOtSFkQu5vC6z9EUOJamwLEMZrfzVvIZBrNbGTS30pX8MyGtmbhvMe2hM2gNnkJT4DiEENj5DBsG72fz4A8QCI5v/CbtkQvfpmUgvYptA/fQFDyf1sgl+0sdx7S6GUj9lqbIZSV3CwGYuY0IEUBTK9vcZoUjcHnk7SYO6QiEEE8BLZPsukFK+SvnmBsoTDF5tFQBQogrgSsBFiyodL+mG3D/HYpwlkOWMlPR8xj6MkbSjyOlnHStnYNRF/wAnQO30T3yYFFHALAg+nGGMmvZ3H8buhKjJXzRvjr8x3BW+y/ZNHAfW4d+zK7kkyyK/m9agmeQ8B3F0b593U77a5RSMpLbSnfqWXYM/4yUtYvW0DkcWf8lAtrbm8twdh3rej5FQO9gRcMtB73WntEfIbFpCl9W9JiDkc1twqctq/AzBJLCGMGBy2a7jlkQgbuBQzoCKeXZB9svhPh74ELgLLkv/u4CJk5baHO2TVb/fcB9ACeeeOIcceEuv8zxiKCyjsCnL2MoOYid34tW4nIIivDRHLmcrqE7SZkbCBaZhimEwhGNd7Ou50re6L0eRRhvW9JBV8IcUf8lFkQu5rX+u9ky+EM2Dz6IocRpDp5OxFhMSGsjoLWQtQdIWm+RzHXSnXqWlFVY/C7uO4JjG79OY+CkA84/kn2dtd0fR1frOa75B+NjFZNh50fpGX2EeOAc/HpHSfYYw7Q2EfSdPq2yU6fwvXB9ROD2duYiZtQ1JIQ4D7gGeK+UcuLawL8GfiyEuIvCYPFS4IWZnMujeoxFBJWcNQRg6IUf72zu9ZIdAUBT+FL2DN/HzsFvs6zpgaLHqYqPo5ru5ZXuT/Da3qvJS5uW8Nu7biLGYk5t+WdMe5i96b+wO/kM3ann6JywbtEYmghT5z+GpfEraA6e9rapoRMZzr7K2u4rUEWI45ofmjTfwET2DN+HnR9kXvSqKVz9gVh2H5a9B58+uVMsF2NdhmI2RAQeU2KmYwTfAXzAk064+7yU8lNSyteEEI8Br1PoMvqMN2NoIu6+U6lW15DfWUwtY64nVOJ8eQBNjTMvdhU7B+9gJLOKiP+U4scqIY5p/lfWdX+S13uvZii7hsMS16Dut1SzoUaZHz6X+eHCQG0uP0oyt5O0tQefWkdIb8NQEgft3pHSpnP4h2wbuBtDa+S45ocJ6Ad/wCtn97Bn5H4SwQ8Q9h1TghX2kc29BoDPOHJa5afK2PfC/Y7A3e3MTczIEUgpi05NkFLeCtw6k/rfmbi/z3LcEVBZR6CqCTR1/vgP2HRoDn+U7pGH6Ry8nZXNvzjoD3TBGTzAtsF76Bz+Af3pP3N4w7eI+Y8vWkZXwsR9K4j7pnaXnc518nrvdQxl19AQOIsVDd/AUBsOWa5r6J+Q0qQt9uVDHluMrLMqq7/E1UpLZfY4ApgN7c0NeE8WexzAvoigsl1DAH7jyPEfsOmgKH7mx75I0lzLQPq3hzxeVXwsrbuW41oeRmLx4p6PsL7nCyTNrdPWAGDafWzqu5VVuy5k1NzAyoZvcVTTd6fkBNK5rewd/QmN4f8z7bEBgExuPZraiqrWTbuOqTDeNYTbxwg8pornCGqCu0PWanUNAfj0IzGtbSWvOTSRhtBfE9CXs3PgW0XTWO5Pwn8yJ7f+moWxj9OX/hMv7PorNvR+lVSutNSvpt3H9sHv8JedZ9M18ihNwfM5Zf5vmBf+0JRmQkkp6Ry4BUUEaI19rqRz70/WXI9fr2y3EEyMCDxH8E7BW0e26rg/VBVVmjUE4DeOBiSZ3HqCvuJ9/AdDCJWFiW+woedv6Rq6kwWJr06pnKaEWZK4mvboR9kxeC9dIz9h1+hjBPXFNATeR8J/Kj6tBUNtQFfiWPlRTLsX097LYPYl+lJ/ZNhcB0gag+eyOPEFQvrikrT3pX7JUOaPtMdvQleLP2R2KPL5UUxrK5HgxdOuY6qMfy+8rqF3DJ4j8DiAakYEfmeZ5az56rQdAUDEfwpN4UvpHnmQuuAFhH0nTLmsodaxrP5GFsSuYG/qSXpTf6Rz+Ee8NfzgQctFjaNZFP8cjcGzCRuTp5s8GDm7h7cGbiZsnEBz5PKSy08kk3sNkOP2rCSza4zAYyp4jqDqzIY7lOpFBJragqo0kSmSdawU2uLXMph+mu1913DEvP9EKfGHyq/Noz16Ge3Ry7Dyo4yaGzDtvZh2H6bdj6ZEMdQGDLWesLF0Sv3/B+PN/q+Rz6fpqL+jaHKaqZJ17FeIsCrMuCMIVP5cHlXBcwQ1we1jBApgVPw5gjH8xpFkctMfMB5DVcJ01N3Opr2XsXvoXubHvzTtujQlTNx/4ow1FWMg9XsG0r+jLXYNgTKsC5Qx16MqjVVISDPxOQJv0bl3Ct5gcdWZDRFBoZFXMh/BRHz6SszcFqSceZ7kWOAM6oIXs3v4e6TMDWVQV36s/BBvDtxMQF9Bc/QTZakzm3sDnz75stblZmxV2tkxWDw72lut8RxBTXD/nYoQPiTViQgKP2C5aa1EOhkLEjeiKXG29H4aOz9SljrLhZR5tvV+EcvupaPuDhQx83SPUtqY1ubqOYLxiMDdjkDOgnbmFsriCIQQVwshpBCiwXkvhBD/z0lMs04IUfyJnTnHLLlDEUZVniMA8OmHA4W72nKgqw0safguWesttvd9uaK5l0tl9/D3GMo8TXviq9N+gnh/ctZ2pMzgM6rjCPZ1GbrbERSYJe2txszYEQgh2oFzgIkZRs6nsL7QUgori35vpud5Z+GeH6ZiCHxVGyMw9CWARja3qWx1Rvwn0x6/joH0E/Qlf162emdC0lxP19Cd1AU/SFP40rLVm3WW8q70GkNjjHUZemME7xzKERHcTWHhuYlWvwh4WBZ4HogLIeaV4VzvAGbHHYoQOvtyDlX+XLrWRs7aUdZ6myNXEDZOoHPwdix7oKx1l4qUNm/234Cm1LGw7uBLUZeK6dhN1xaVrc6DISmM5YgydGtVGjFL2lutmWnO4ouALinlK/vtmg90TnjvJaaZdejAzAdvp3w2bSE5q7Sneg+FEAoL677pDM5+raZdRN0jD5I0X6E9cQOaEitr3TnrTVSlDvUgS1yXlfFBffc7Ao+pMaPENMBXKHQLTZu5l5hmdiCEVrWIAMBQOxg2f1X2eoPGSubHPk/X0F2EfcfRHPlo2c9xKIYzf6Fz8FvEA+dQX4Enf3PWDnSto+z1FkNiAWLGzz54uIdpJ6YRQhwFLAJeccLcNuAlIcTJeIlpDsFsuMxqRwQLyOcHsfNDqGW+Y54X/SxJcx2dA7cS85+JX19Y1voPRl5m2db3RfxaB4vr7yxrl9AYOasTv+/YstdbFJljdkQDs6GduYNpdw1JKV+VUjZJKTuklB0Uun+Ol1LuoZCY5jJn9tCpwJCUcnd5JM92ZkefpRAq1UwhoTrr7Nh2X9nrFkJhYeJWQNA9cvBlI8pNX/JX5Ow9LEjcXLGuGyvfi6aUnthnukisWRQNzI72Vmsq9RzBb4FtwBbgX4HppVzyqCEKkK/a2VQlAYCdr8ygrqE1Ux+6iN7kY1j2YEXOsT9SSrpH7iegryDqP60i58jLLFImx+1XFaTEewTpnUXZPk0nMuh1Xksp5WeklEuklEdJKdeU6zwe1UIFqhgRKIU19O18f8XO0Rz5BHmZpmf00YqdYyLDmf8mndtES/TKinQJAeQdx1npHARvx6bw/fB4p+CtNVRlFHUeilalB39mgKYvqdoDZQCa2ohPPwJRwb7noLGcmP+95Oyeip1jIll7Fz61nbrghYc+eJpImcOvH4OmHjwVZjlR1XnoVch7MFMUtQOUqeWnmOsINz11KYQYATbWWscUaAB6ay1iCng6y8ts0DkbNIKns9wsl1JOexDKbRHBRill5ZZ8LBNCiDWezvLh6Swfs0EjeDrLjRBiRt3v3oiPh4eHxxzHcwQeHh4ecxy3OYL7ai1ging6y4uns3zMBo3g6Sw3M9LpqsFiDw8PD4/q47aIwMPDw8OjyrjGEQghzhNCbHSS2VxXaz0TEULsEEK8KoRYOzY6L4SoE0I8KYTY7Pyv4qOd47oeFEL0CCHWT9g2qa5aJQsqovFmIUSXY8+1QogLJuy73tG4UQhxbjU0OudtF0I8I4R4XQjxmhDi8852t9mzmE5X2VQI4RdCvCCEeMXR+XVn+yIhxCpHz0+Fk9RACOFz3m9x9nfUUONDQojtE2x5rLO9Jp/5BL2qEOJlIcTjzvvy2VJKWfM/Co8pbgUWAwbwCnB4rXVN0LcDaNhv2z8C1zmvrwPuqIGuM4DjgfWH0gVcAPyOwuIrpwKraqjxZuAfJjn2cOez91FY0HAroFZJ5zwKa2UBRIBNjh632bOYTlfZ1LFL2HmtA6scOz0GXOJs/z7waef1VcD3ndeXAD+tocaHgA9PcnxNPvMJ5/8S8GPgced92WzplojgZGCLlHKbLKQ/+gmF5DZu5iLgh87rHwIXV1uAlPJPwP5rMhTTVZNkQUU0FuMi4CdSyqyUcjuFtapOrpi4CUgpd0spX3JejwBvUMih4TZ7FtNZjJrY1LHLqPNWd/4kcCYwljJuf3uO2fnnwFlCVGhdjkNrLEbNEm4JIdqADwD3O+8FZbSlWxyB2xPZSOC/hBAvikL+BIBmuW9F1T1Ac22kHUAxXW6z8Wed8PrBCd1qrtDohNLHUbhDdK0999MJLrOp05WxFugBnqQQjQzKfYkuJmoZ1+nsHwLqq61RSjlmy1sdW94thBhLzlzLz/weCpkgx1aCrKeMtnSLI3A7p0kpj6eQi/kzQogzJu6UhRjMddOv3KqLQg7rJcCxwG7gzpqqmYAQIgz8O/AFKeXwxH1usuckOl1nUymlLaU8lkI+kpOB6iRVLoH9NQohjgSup6D1JKAOuLZ2CkEIcSHQI6V8sVLncIsjmHIim1ogpexy/vcAv6Dwpe4eCwud/9VZyezQFNPlGhtLKbudBpinsEz5WFdFTTWKQhLefwcelVL+h7PZdfacTKdbbepoGwSeAd5FoTtlbGmbiVrGdTr7Y0D5k1McWuN5TveblIVVF39A7W35HuCDQogdFLrNzwT+iTLa0i2OYDWw1BkFNygMcPy6xpoAEEKEhBCRsdcUUnOup6Dvcuewy4Hy51mcHsV0uSZZ0H79qh+iYE8oaLzEmfWwCFgKvFAlTQJ4AHhDSnnXhF2usmcxnW6zqRCiUQgRd14HgPdTGM94Bviwc9j+9hyz84eBp50IrNoaN0xw/IJCv/tEW1b9M5dSXi+lbJOFBGCXULDNRyinLSs90j3VPwoj8pso9CPeUGs9E3QtpjDr4hXgtTFtFPrc/gBsBp4C6mqg7d8odAPkKPQRXlFMF4WZDt917PsqcGINNf7I0bDO+dLOm3D8DY7GjcD5VbTlaRS6fdYBa52/C1xoz2I6XWVT4GjgZUfPeuAmZ/tiCo5oC/AzwOds9zvvtzj7F9dQ49OOLdcDj7BvZlFNPvP9NL+PfbOGymZL78liDw8PjzmOW7qGPDw8PDxqhOcIPDw8POY4niPw8PDwmON4jsDDw8NjjuM5Ag8PD485jucIPDw8POY4niPw8PDwmON4jsDDw8NjjvP/AbbluWOC/gAVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = np.linspace(-H / Nlay / 2, -H + H / Nlay / 2, Nlay)\n",
    "c = plt.contour(x, z, h[:, 50, :], np.arange(90, 100.1, 0.2))\n",
    "plt.axis(\"scaled\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-30T20:16:50.965304Z",
     "iopub.status.busy": "2022-07-30T20:16:50.964925Z",
     "iopub.status.idle": "2022-07-30T20:16:50.970878Z",
     "shell.execute_reply": "2022-07-30T20:16:50.969892Z"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    # ignore PermissionError on Windows\n",
    "    temp_dir.cleanup()\n",
    "except:\n",
    "    pass"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}