examples/groundwater_paper/Notebooks/example_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example for FloPy methods note"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the `modflow` and `utils` subpackages of FloPy and give them the aliases `fpm` and `fpu`, respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flopy is installed in /Users/jdhughes/Documents/Development/flopy_git/flopy_us/flopy\n",
      "3.7.3 (default, Mar 27 2019, 16:54:48) \n",
      "[Clang 4.0.1 (tags/RELEASE_401/final)]\n",
      "numpy version: 1.16.2\n",
      "matplotlib version: 3.0.3\n",
      "flopy version: 3.2.12\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import flopy\n",
    "import flopy.modflow as fpm\n",
    "import flopy.utils as fpu\n",
    "import matplotlib as mpl\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": [
    "Create a MODFLOW model object. Here, the MODFLOW model object is stored in a Python variable called {\\tt model}, but this can be an arbitrary name. This object name is important as it will be used as a reference to the model in the remainder of the FloPy script. In addition, a {\\tt modelname} is specified when the MODFLOW model object is created. This {\\tt modelname} is used for all the files that are created by FloPy for this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "exe = 'mf2005'\n",
    "ws = os.path.join('temp')\n",
    "model = fpm.Modflow(modelname='gwexample', exe_name=exe, model_ws=ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The discretization of the model is specified with the discretization file (DIS) of MODFLOW. The aquifer is divided into 201 cells of length 10 m and width 1 m.  The first input of the discretization package is the name of the model object. All other input arguments are self explanatory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Discretization Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.Modflow`) to which\n",
       "        this package will be added.\n",
       "    nlay : int\n",
       "        Number of model layers (the default is 1).\n",
       "    nrow : int\n",
       "        Number of model rows (the default is 2).\n",
       "    ncol : int\n",
       "        Number of model columns (the default is 2).\n",
       "    nper : int\n",
       "        Number of model stress periods (the default is 1).\n",
       "    delr : float or array of floats (ncol), optional\n",
       "        An array of spacings along a row (the default is 1.0).\n",
       "    delc : float or array of floats (nrow), optional\n",
       "        An array of spacings along a column (the default is 0.0).\n",
       "    laycbd : int or array of ints (nlay), optional\n",
       "        An array of flags indicating whether or not a layer has a Quasi-3D\n",
       "        confining bed below it. 0 indicates no confining bed, and not zero\n",
       "        indicates a confining bed. LAYCBD for the bottom layer must be 0. (the\n",
       "        default is 0)\n",
       "    top : float or array of floats (nrow, ncol), optional\n",
       "        An array of the top elevation of layer 1. For the common situation in\n",
       "        which the top layer represents a water-table aquifer, it may be\n",
       "        reasonable to set Top equal to land-surface elevation (the default is\n",
       "        1.0)\n",
       "    botm : float or array of floats (nlay, nrow, ncol), optional\n",
       "        An array of the bottom elevation for each model cell (the default is\n",
       "        0.)\n",
       "    perlen : float or array of floats (nper)\n",
       "        An array of the stress period lengths.\n",
       "    nstp : int or array of ints (nper)\n",
       "        Number of time steps in each stress period (default is 1).\n",
       "    tsmult : float or array of floats (nper)\n",
       "        Time step multiplier (default is 1.0).\n",
       "    steady : boolean or array of boolean (nper)\n",
       "        true or False indicating whether or not stress period is steady state\n",
       "        (default is True).\n",
       "    itmuni : int\n",
       "        Time units, default is days (4)\n",
       "    lenuni : int\n",
       "        Length units, default is meters (2)\n",
       "    extension : string\n",
       "        Filename extension (default is 'dis')\n",
       "    unitnumber : int\n",
       "        File unit number (default is None).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package. If filenames=None the package name\n",
       "        will be created using the model name and package extension. If a\n",
       "        single string is passed the package will be set to the string.\n",
       "        Default is None.\n",
       "    xul : float\n",
       "        x coordinate of upper left corner of the grid, default is None\n",
       "    yul : float\n",
       "        y coordinate of upper left corner of the grid, default is None\n",
       "    rotation : float\n",
       "        clockwise rotation (in degrees) of the grid about the upper left\n",
       "        corner. default is 0.0\n",
       "    proj4_str : str\n",
       "        PROJ4 string that defines the projected coordinate system\n",
       "        (e.g. '+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs ').\n",
       "        Can be an EPSG code (e.g. 'EPSG:32614'). Default is None.\n",
       "    start_datetime : str\n",
       "        starting datetime of the simulation. default is '1/1/1970'\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "    heading : str\n",
       "        Text string written to top of package input file.\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> dis = flopy.modflow.ModflowDis(m)\n",
       "\n",
       "     _ModflowDis__thickness = <flopy.utils.util_array.Util3d object at 0x114760b38> ('flopy.utils.util_array.Util3d)\n",
       " _name = DIS\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " _sr = xul:0; yul:1; rotation:0; proj4_str:None; units:meters; lenuni:2; length_multiplier:1.0 ('flopy.utils.reference.SpatialReference)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " botm = <flopy.utils.util_array.Util3d object at 0x114760668> ('flopy.utils.util_array.Util3d)\n",
       " delc = <flopy.utils.util_array.Util2d object at 0x114760550> ('flopy.utils.util_array.Util2d)\n",
       " delr = <flopy.utils.util_array.Util2d object at 0x1147604a8> ('flopy.utils.util_array.Util2d)\n",
       " extra = \n",
       " file_name = gwexample.dis\n",
       " fn_path = temp/gwexample.dis ('str)\n",
       " itmuni = 4 ('int)\n",
       " itmuni_dict = {0: 'undefined', 1: 'seconds', 2: 'minutes', 3: 'hours', 4: 'days', 5: 'years'} ('dict)\n",
       " laycbd = <flopy.utils.util_array.Util2d object at 0x114760438> ('flopy.utils.util_array.Util2d)\n",
       " lenuni = 2 ('int)\n",
       " ncol = 201 ('int)\n",
       " nlay = 1 ('int)\n",
       " nper = 1 ('int)\n",
       " nrow = 1 ('int)\n",
       " nstp = <flopy.utils.util_array.Util2d object at 0x1147607b8> ('flopy.utils.util_array.Util2d)\n",
       " perlen = <flopy.utils.util_array.Util2d object at 0x1147606a0> ('flopy.utils.util_array.Util2d)\n",
       " start_datetime = 1-1-1970 ('str)\n",
       " steady = <flopy.utils.util_array.Util2d object at 0x114760908> ('flopy.utils.util_array.Util2d)\n",
       " top = <flopy.utils.util_array.Util2d object at 0x1147605f8> ('flopy.utils.util_array.Util2d)\n",
       " tr = <flopy.utils.reference.TemporalReference object at 0x103cc3630> ('flopy.utils.reference.TemporalReference)\n",
       " tsmult = <flopy.utils.util_array.Util2d object at 0x114760860> ('flopy.utils.util_array.Util2d)\n",
       " unit_number = 11"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpm.ModflowDis(model, nlay=1, nrow=1, ncol=201,\n",
    "               delr=10, delc=1, top=50, botm=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Active cells and the like are defined with the Basic package (BAS), which is required for every MODFLOW model. It contains the {\\tt ibound} array, which is used to specify which cells are active (value is positive), inactive (value is 0), or fixed head (value is negative). The {\\tt numpy} package (aliased as {\\tt np}) can be used to quickly initialize the {\\tt ibound} array with values of 1, and then set the {\\tt ibound} value for the first and last columns to -1. The {\\tt numpy} package (and Python, in general) uses zero-based indexing and supports negative indexing so that row 1 and column 1, and row 1 and column 201, can be referenced as [0, 0], and [0, -1], respectively.  Although this simulation is for steady flow, starting heads still need to be specified. They are used as the head for fixed-head cells (where {\\tt ibound} is negative), and as a starting point to compute the saturated thickness for cases of unconfined flow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Basic Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which\n",
       "        this package will be added.\n",
       "    ibound : array of ints, optional\n",
       "        The ibound array (the default is 1).\n",
       "    strt : array of floats, optional\n",
       "        An array of starting heads (the default is 1.0).\n",
       "    ifrefm : bool, optional\n",
       "        Indication if data should be read using free format (the default is\n",
       "        True).\n",
       "    ixsec : bool, optional\n",
       "        Indication of whether model is cross sectional or not (the default is\n",
       "        False).\n",
       "    ichflg : bool, optional\n",
       "        Flag indicating that flows between constant head cells should be\n",
       "        calculated (the default is False).\n",
       "    stoper : float\n",
       "        percent discrepancy that is compared to the budget percent discrepancy\n",
       "        continue when the solver convergence criteria are not met.  Execution\n",
       "        will unless the budget percent discrepancy is greater than stoper\n",
       "        (default is None). MODFLOW-2005 only\n",
       "    hnoflo : float\n",
       "        Head value assigned to inactive cells (default is -999.99).\n",
       "    extension : str, optional\n",
       "        File extension (default is 'bas').\n",
       "    unitnumber : int, optional\n",
       "        FORTRAN unit number for this package (default is None).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package. If filenames=None the package name\n",
       "        will be created using the model name and package extension. If a single\n",
       "        string is passed the package name will be set to the string.\n",
       "        Default is None.\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "    heading : str\n",
       "        Text string written to top of package input file.\n",
       "    options : list of str\n",
       "        Can be either or a combination of XSECTION, CHTOCH or FREE.\n",
       "    ifrefm : bool\n",
       "        Indicates whether or not packages will be written as free format.\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> bas = flopy.modflow.ModflowBas(m)\n",
       "\n",
       "     _name = BAS6\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " extra = \n",
       " file_name = gwexample.bas\n",
       " fn_path = temp/gwexample.bas ('str)\n",
       " hnoflo = -999.99 ('float)\n",
       " ibound = <flopy.utils.util_array.Util3d object at 0x114671eb8> ('flopy.utils.util_array.Util3d)\n",
       " ichflg = False ('bool)\n",
       " ixsec = False ('bool)\n",
       " options =  ('str)\n",
       " stoper = None ('NoneType)\n",
       " strt = <flopy.utils.util_array.Util3d object at 0x103cc32b0> ('flopy.utils.util_array.Util3d)\n",
       " unit_number = 13"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ibound = np.ones((1, 201))\n",
    "ibound[0, 0] = ibound[0, -1] = -1\n",
    "fpm.ModflowBas(model, ibound=ibound, strt=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hydraulic properties of the aquifer are specified with the Layer Properties Flow (LPF) package (alternatively, the Block Centered Flow (BCF) package may be used). Only the hydraulic conductivity of the aquifer and the layer type ({\\tt laytyp}) need to be specified. The latter is set to 1, which means that MODFLOW will calculate the saturated thickness differently depending on whether or not the head is above the top of the aquifer.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Layer Property Flow Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which\n",
       "        this package will be added.\n",
       "    ipakcb : int\n",
       "        A flag that is used to determine if cell-by-cell budget data should be\n",
       "        saved. If ipakcb is non-zero cell-by-cell budget data will be saved.\n",
       "        (default is 53)\n",
       "    hdry : float\n",
       "        Is the head that is assigned to cells that are converted to dry during\n",
       "        a simulation. Although this value plays no role in the model\n",
       "        calculations, it is useful as an indicator when looking at the\n",
       "        resulting heads that are output from the model. HDRY is thus similar\n",
       "        to HNOFLO in the Basic Package, which is the value assigned to cells\n",
       "        that are no-flow cells at the start of a model simulation.\n",
       "        (default is -1.e30).\n",
       "    laytyp : int or array of ints (nlay)\n",
       "        Layer type, contains a flag for each layer that specifies the layer\n",
       "        type.\n",
       "        0 confined\n",
       "        >0 convertible\n",
       "        <0 convertible unless the THICKSTRT option is in effect.\n",
       "        (default is 0).\n",
       "    layavg : int or array of ints (nlay)\n",
       "        Layer average\n",
       "        0 is harmonic mean\n",
       "        1 is logarithmic mean\n",
       "        2 is arithmetic mean of saturated thickness and logarithmic mean of\n",
       "        of hydraulic conductivity\n",
       "        (default is 0).\n",
       "    chani : float or array of floats (nlay)\n",
       "        contains a value for each layer that is a flag or the horizontal\n",
       "        anisotropy. If CHANI is less than or equal to 0, then variable HANI\n",
       "        defines horizontal anisotropy. If CHANI is greater than 0, then CHANI\n",
       "        is the horizontal anisotropy for the entire layer, and HANI is not\n",
       "        read. If any HANI parameters are used, CHANI for all layers must be\n",
       "        less than or equal to 0. Use as many records as needed to enter a\n",
       "        value of CHANI for each layer. The horizontal anisotropy is the ratio\n",
       "        of the hydraulic conductivity along columns (the Y direction) to the\n",
       "        hydraulic conductivity along rows (the X direction).\n",
       "        (default is 1).\n",
       "    layvka : int or array of ints (nlay)\n",
       "        a flag for each layer that indicates whether variable VKA is vertical\n",
       "        hydraulic conductivity or the ratio of horizontal to vertical\n",
       "        hydraulic conductivity.\n",
       "        0: VKA is vertical hydraulic conductivity\n",
       "        not 0: VKA is the ratio of horizontal to vertical hydraulic conductivity\n",
       "        (default is 0).\n",
       "    laywet : int or array of ints (nlay)\n",
       "        contains a flag for each layer that indicates if wetting is active.\n",
       "        0 wetting is inactive\n",
       "        not 0 wetting is active\n",
       "        (default is 0).\n",
       "    wetfct : float\n",
       "        is a factor that is included in the calculation of the head that is\n",
       "        initially established at a cell when it is converted from dry to wet.\n",
       "        (default is 0.1).\n",
       "    iwetit : int\n",
       "        is the iteration interval for attempting to wet cells. Wetting is\n",
       "        attempted every IWETIT iteration. If using the PCG solver\n",
       "        (Hill, 1990), this applies to outer iterations, not inner iterations.\n",
       "        If IWETIT  less than or equal to 0, it is changed to 1.\n",
       "        (default is 1).\n",
       "    ihdwet : int\n",
       "        is a flag that determines which equation is used to define the\n",
       "        initial head at cells that become wet.\n",
       "        (default is 0)\n",
       "    hk : float or array of floats (nlay, nrow, ncol)\n",
       "        is the hydraulic conductivity along rows. HK is multiplied by\n",
       "        horizontal anisotropy (see CHANI and HANI) to obtain hydraulic\n",
       "        conductivity along columns.\n",
       "        (default is 1.0).\n",
       "    hani : float or array of floats (nlay, nrow, ncol)\n",
       "        is the ratio of hydraulic conductivity along columns to hydraulic\n",
       "        conductivity along rows, where HK of item 10 specifies the hydraulic\n",
       "        conductivity along rows. Thus, the hydraulic conductivity along\n",
       "        columns is the product of the values in HK and HANI.\n",
       "        (default is 1.0).\n",
       "    vka : float or array of floats (nlay, nrow, ncol)\n",
       "        is either vertical hydraulic conductivity or the ratio of horizontal\n",
       "        to vertical hydraulic conductivity depending on the value of LAYVKA.\n",
       "        (default is 1.0).\n",
       "    ss : float or array of floats (nlay, nrow, ncol)\n",
       "        is specific storage unless the STORAGECOEFFICIENT option is used.\n",
       "        When STORAGECOEFFICIENT is used, Ss is confined storage coefficient.\n",
       "        (default is 1.e-5).\n",
       "    sy : float or array of floats (nlay, nrow, ncol)\n",
       "        is specific yield.\n",
       "        (default is 0.15).\n",
       "    vkcb : float or array of floats (nlay, nrow, ncol)\n",
       "        is the vertical hydraulic conductivity of a Quasi-three-dimensional\n",
       "        confining bed below a layer. (default is 0.0).  Note that if an array\n",
       "        is passed for vkcb it must be of size (nlay, nrow, ncol) even though\n",
       "        the information for the bottom layer is not needed.\n",
       "    wetdry : float or array of floats (nlay, nrow, ncol)\n",
       "        is a combination of the wetting threshold and a flag to indicate\n",
       "        which neighboring cells can cause a cell to become wet.\n",
       "        (default is -0.01).\n",
       "    storagecoefficient : boolean\n",
       "        indicates that variable Ss and SS parameters are read as storage\n",
       "        coefficient rather than specific storage.\n",
       "        (default is False).\n",
       "    constantcv : boolean\n",
       "         indicates that vertical conductance for an unconfined cell is\n",
       "         computed from the cell thickness rather than the saturated thickness.\n",
       "         The CONSTANTCV option automatically invokes the NOCVCORRECTION\n",
       "         option. (default is False).\n",
       "    thickstrt : boolean\n",
       "        indicates that layers having a negative LAYTYP are confined, and their\n",
       "        cell thickness for conductance calculations will be computed as\n",
       "        STRT-BOT rather than TOP-BOT. (default is False).\n",
       "    nocvcorrection : boolean\n",
       "        indicates that vertical conductance is not corrected when the vertical\n",
       "        flow correction is applied. (default is False).\n",
       "    novfc : boolean\n",
       "         turns off the vertical flow correction under dewatered conditions.\n",
       "         This option turns off the vertical flow calculation described on p.\n",
       "         5-8 of USGS Techniques and Methods Report 6-A16 and the vertical\n",
       "         conductance correction described on p. 5-18 of that report.\n",
       "         (default is False).\n",
       "    extension : string\n",
       "        Filename extension (default is 'lpf')\n",
       "    unitnumber : int\n",
       "        File unit number (default is None).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package and the output files. If\n",
       "        filenames=None the package name will be created using the model name\n",
       "        and package extension and the cbc output name will be created using\n",
       "        the model name and .cbc extension (for example, modflowtest.cbc),\n",
       "        if ipakcbc is a number greater than zero. If a single string is passed\n",
       "        the package will be set to the string and cbc output name will be\n",
       "        created using the model name and .cbc extension, if ipakcbc is a\n",
       "        number greater than zero. To define the names for all package files\n",
       "        (input and output) the length of the list of strings should be 2.\n",
       "        Default is None.\n",
       "\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> lpf = flopy.modflow.ModflowLpf(m)\n",
       "\n",
       "     _name = LPF\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " chani = <flopy.utils.util_array.Util2d object at 0x103c9db38> ('flopy.utils.util_array.Util2d)\n",
       " extra = \n",
       " file_name = gwexample.lpf\n",
       " fn_path = temp/gwexample.lpf ('str)\n",
       " hani = <flopy.utils.util_array.Util3d object at 0x114786198> ('flopy.utils.util_array.Util3d)\n",
       " hdry = -1e+30 ('float)\n",
       " hk = <flopy.utils.util_array.Util3d object at 0x114786160> ('flopy.utils.util_array.Util3d)\n",
       " ihdwet = 0 ('int)\n",
       " ipakcb = 0 ('int)\n",
       " iwetit = 1 ('int)\n",
       " layavg = <flopy.utils.util_array.Util2d object at 0x103cc3278> ('flopy.utils.util_array.Util2d)\n",
       " laytyp = <flopy.utils.util_array.Util2d object at 0x103cc3550> ('flopy.utils.util_array.Util2d)\n",
       " layvka = <flopy.utils.util_array.Util2d object at 0x103c9dba8> ('flopy.utils.util_array.Util2d)\n",
       " laywet = <flopy.utils.util_array.Util2d object at 0x114786080> ('flopy.utils.util_array.Util2d)\n",
       " nplpf = 0 ('int)\n",
       " options =   ('str)\n",
       " ss = <flopy.utils.util_array.Util3d object at 0x114786390> ('flopy.utils.util_array.Util3d)\n",
       " sy = <flopy.utils.util_array.Util3d object at 0x114786320> ('flopy.utils.util_array.Util3d)\n",
       " unit_number = 15\n",
       " vka = <flopy.utils.util_array.Util3d object at 0x1147862e8> ('flopy.utils.util_array.Util3d)\n",
       " vkcb = <flopy.utils.util_array.Util3d object at 0x1147864a8> ('flopy.utils.util_array.Util3d)\n",
       " wetdry = <flopy.utils.util_array.Util3d object at 0x114786550> ('flopy.utils.util_array.Util3d)\n",
       " wetfct = 0.1 ('float)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpm.ModflowLpf(model, hk=10, laytyp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aquifer recharge is simulated with the Recharge package (RCH) and the extraction of water at the two ditches is simulated with the Well package (WEL). The latter requires specification of the layer, row, column, and injection rate of the well for each stress period. The layers, rows, columns, and the stress period are numbered (consistent with Python's zero-based numbering convention) starting at 0. The required data are stored in a Python dictionary ({\\tt lrcQ} in the code below), which is used in FloPy to store data that can vary by stress period. The {\\tt lrcQ} dictionary specifies that two wells (one in cell 1, 1, 51 and one in cell 1, 1, 151), each with a rate of -1 m$^3$/m/d, will be active for the first stress period.  Because this is a steady-state model, there is only one stress period and therefore only one entry in the dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Well Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which\n",
       "        this package will be added.\n",
       "    ipakcb : int\n",
       "        A flag that is used to determine if cell-by-cell budget data should be\n",
       "        saved. If ipakcb is non-zero cell-by-cell budget data will be saved.\n",
       "        (default is 0).\n",
       "    stress_period_data : list of boundaries, or recarray of boundaries, or\n",
       "        dictionary of boundaries\n",
       "        Each well is defined through definition of\n",
       "        layer (int), row (int), column (int), flux (float).\n",
       "        The simplest form is a dictionary with a lists of boundaries for each\n",
       "        stress period, where each list of boundaries itself is a list of\n",
       "        boundaries. Indices of the dictionary are the numbers of the stress\n",
       "        period. This gives the form of:\n",
       "\n",
       "            stress_period_data =\n",
       "            {0: [\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux]\n",
       "                ],\n",
       "            1:  [\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux]\n",
       "                ], ...\n",
       "            kper:\n",
       "                [\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux],\n",
       "                [lay, row, col, flux]\n",
       "                ]\n",
       "            }\n",
       "\n",
       "        Note that if the number of lists is smaller than the number of stress\n",
       "        periods, then the last list of wells will apply until the end of the\n",
       "        simulation. Full details of all options to specify stress_period_data\n",
       "        can be found in the flopy3 boundaries Notebook in the basic\n",
       "        subdirectory of the examples directory\n",
       "    dtype : custom datatype of stress_period_data.\n",
       "        If None the default well datatype will be applied (default is None).\n",
       "    extension : string\n",
       "        Filename extension (default is 'wel')\n",
       "    options : list of strings\n",
       "        Package options (default is None).\n",
       "    unitnumber : int\n",
       "        File unit number (default is None).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package and the output files. If\n",
       "        filenames=None the package name will be created using the model name\n",
       "        and package extension and the cbc output name will be created using\n",
       "        the model name and .cbc extension (for example, modflowtest.cbc),\n",
       "        if ipakcbc is a number greater than zero. If a single string is passed\n",
       "        the package will be set to the string and cbc output names will be\n",
       "        created using the model name and .cbc extension, if ipakcbc is a\n",
       "        number greater than zero. To define the names for all package files\n",
       "        (input and output) the length of the list of strings should be 2.\n",
       "        Default is None.\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "    mxactw : int\n",
       "        Maximum number of wells for a stress period.  This is calculated\n",
       "        automatically by FloPy based on the information in\n",
       "        stress_period_data.\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "    Parameters are not supported in FloPy.\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> lrcq = {0:[[2, 3, 4, -100.]], 1:[[2, 3, 4, -100.]]}\n",
       "    >>> wel = flopy.modflow.ModflowWel(m, stress_period_data=lrcq)\n",
       "\n",
       "     _name = WEL\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " dtype = [('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4')] ('numpy.dtype)\n",
       " extra = \n",
       " file_name = gwexample.wel\n",
       " fn_path = temp/gwexample.wel ('str)\n",
       " ipakcb = 0 ('int)\n",
       " iunitramp = None ('NoneType)\n",
       " np = 0 ('int)\n",
       " options (list, items = 0)\n",
       " phiramp = None ('NoneType)\n",
       " specify = False ('bool)\n",
       " stress_period_data = <flopy.utils.util_list.MfList object at 0x11478e208> ('flopy.utils.util_list.MfList)\n",
       " unit_number = 20"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpm.ModflowRch(model, rech=0.001)\n",
    "lrcQ = {0: [[0, 0, 50, -1], [0, 0, 150, -1]]}\n",
    "fpm.ModflowWel(model, stress_period_data=lrcQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Preconditioned Conjugate-Gradient (PCG) solver, using the default settings, is specified to solve the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Pcg Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which\n",
       "        this package will be added.\n",
       "    mxiter : int\n",
       "        maximum number of outer iterations. (default is 50)\n",
       "    iter1 : int\n",
       "        maximum number of inner iterations. (default is 30)\n",
       "    npcond : int\n",
       "        flag used to select the matrix conditioning method. (default is 1).\n",
       "        specify npcond = 1 for Modified Incomplete Cholesky.\n",
       "        specify npcond = 2 for Polynomial.\n",
       "    hclose : float\n",
       "        is the head change criterion for convergence. (default is 1e-5).\n",
       "    rclose : float\n",
       "        is the residual criterion for convergence. (default is 1e-5)\n",
       "    relax : float\n",
       "        is the relaxation parameter used with npcond = 1. (default is 1.0)\n",
       "    nbpol : int\n",
       "        is only used when npcond = 2 to indicate whether the estimate of the\n",
       "        upper bound on the maximum eigenvalue is 2.0, or whether the estimate\n",
       "        will be calculated. nbpol = 2 is used to specify the value is 2.0;\n",
       "        for any other value of nbpol, the estimate is calculated. Convergence\n",
       "        is generally insensitive to this parameter. (default is 0).\n",
       "    iprpcg : int\n",
       "        solver print out interval. (default is 0).\n",
       "    mutpcg : int\n",
       "        If mutpcg = 0, tables of maximum head change and residual will be\n",
       "        printed each iteration.\n",
       "        If mutpcg = 1, only the total number of iterations will be printed.\n",
       "        If mutpcg = 2, no information will be printed.\n",
       "        If mutpcg = 3, information will only be printed if convergence fails.\n",
       "        (default is 3).\n",
       "    damp : float\n",
       "        is the steady-state damping factor. (default is 1.)\n",
       "    dampt : float\n",
       "        is the transient damping factor. (default is 1.)\n",
       "    ihcofadd : int\n",
       "        is a flag that determines what happens to an active cell that is\n",
       "        surrounded by dry cells.  (default is 0). If ihcofadd=0, cell\n",
       "        converts to dry regardless of HCOF value. This is the default, which\n",
       "        is the way PCG2 worked prior to the addition of this option. If\n",
       "        ihcofadd<>0, cell converts to dry only if HCOF has no head-dependent\n",
       "        stresses or storage terms.\n",
       "    extension : list string\n",
       "        Filename extension (default is 'pcg')\n",
       "    unitnumber : int\n",
       "        File unit number (default is None).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package. If filenames=None the package name\n",
       "        will be created using the model name and package extension. If a\n",
       "        single string is passed the package will be set to the string.\n",
       "        Default is None.\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> pcg = flopy.modflow.ModflowPcg(m)\n",
       "\n",
       "     _name = PCG\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " damp = 1.0 ('float)\n",
       " dampt = 1.0 ('float)\n",
       " extra = \n",
       " file_name = gwexample.pcg\n",
       " fn_path = temp/gwexample.pcg ('str)\n",
       " hclose = 1e-05 ('float)\n",
       " ihcofadd = 0 ('int)\n",
       " iprpcg = 0 ('int)\n",
       " iter1 = 30 ('int)\n",
       " mutpcg = 3 ('int)\n",
       " mxiter = 50 ('int)\n",
       " nbpol = 0 ('int)\n",
       " npcond = 1 ('int)\n",
       " rclose = 1e-05 ('float)\n",
       " relax = 1.0 ('float)\n",
       " unit_number = 27"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpm.ModflowPcg(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency and type of output that MODFLOW writes to an output file is specified with the Output Control (OC) package.  In this case, the budget is printed and heads are saved (the default), so no arguments are needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "    MODFLOW Output Control Package Class.\n",
       "\n",
       "    Parameters\n",
       "    ----------\n",
       "    model : model object\n",
       "        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which\n",
       "        this package will be added.\n",
       "    ihedfm : int\n",
       "        is a code for the format in which heads will be printed.\n",
       "        (default is 0).\n",
       "    iddnfm : int\n",
       "        is a code for the format in which drawdown will be printed.\n",
       "        (default is 0).\n",
       "    chedfm : string\n",
       "        is a character value that specifies the format for saving heads.\n",
       "        The format must contain 20 characters or less and must be a valid\n",
       "        Fortran format that is enclosed in parentheses. The format must be\n",
       "        enclosed in apostrophes if it contains one or more blanks or commas.\n",
       "        The optional word LABEL after the format is used to indicate that\n",
       "        each layer of output should be preceded with a line that defines the\n",
       "        output (simulation time, the layer being output, and so forth). If\n",
       "        there is no record specifying CHEDFM, then heads are written to a\n",
       "        binary (unformatted) file. Binary files are usually more compact than\n",
       "        text files, but they are not generally transportable among different\n",
       "        computer operating systems or different Fortran compilers.\n",
       "        (default is None)\n",
       "    cddnfm : string\n",
       "        is a character value that specifies the format for saving drawdown.\n",
       "        The format must contain 20 characters or less and must be a valid\n",
       "        Fortran format that is enclosed in parentheses. The format must be\n",
       "        enclosed in apostrophes if it contains one or more blanks or commas.\n",
       "        The optional word LABEL after the format is used to indicate that\n",
       "        each layer of output should be preceded with a line that defines the\n",
       "        output (simulation time, the layer being output, and so forth). If\n",
       "        there is no record specifying CDDNFM, then drawdowns are written to a\n",
       "        binary (unformatted) file. Binary files are usually more compact than\n",
       "        text files, but they are not generally transportable among different\n",
       "        computer operating systems or different Fortran compilers.\n",
       "        (default is None)\n",
       "    cboufm : string\n",
       "        is a character value that specifies the format for saving ibound.\n",
       "        The format must contain 20 characters or less and must be a valid\n",
       "        Fortran format that is enclosed in parentheses. The format must be\n",
       "        enclosed in apostrophes if it contains one or more blanks or commas.\n",
       "        The optional word LABEL after the format is used to indicate that\n",
       "        each layer of output should be preceded with a line that defines the\n",
       "        output (simulation time, the layer being output, and so forth). If\n",
       "        there is no record specifying CBOUFM, then ibounds are written to a\n",
       "        binary (unformatted) file. Binary files are usually more compact than\n",
       "        text files, but they are not generally transportable among different\n",
       "        computer operating systems or different Fortran compilers.\n",
       "        (default is None)\n",
       "    stress_period_data : dictionary of of lists\n",
       "        Dictionary key is a tuple with the zero-based period and step\n",
       "        (IPEROC, ITSOC) for each print/save option list. If stress_period_data\n",
       "        is None, then heads are saved for the last time step of each stress\n",
       "        period. (default is None)\n",
       "\n",
       "        The list can have any valid MODFLOW OC print/save option:\n",
       "            PRINT HEAD\n",
       "            PRINT DRAWDOWN\n",
       "            PRINT BUDGET\n",
       "            SAVE HEAD\n",
       "            SAVE DRAWDOWN\n",
       "            SAVE BUDGET\n",
       "            SAVE IBOUND\n",
       "\n",
       "            The lists can also include (1) DDREFERENCE in the list to reset\n",
       "            drawdown reference to the period and step and (2) a list of layers\n",
       "            for PRINT HEAD, SAVE HEAD, PRINT DRAWDOWN, SAVE DRAWDOWN, and\n",
       "            SAVE IBOUND.\n",
       "\n",
       "        stress_period_data = {(0,1):['save head']}) would save the head for\n",
       "        the second timestep in the first stress period.\n",
       "\n",
       "    compact : boolean\n",
       "        Save results in compact budget form. (default is True).\n",
       "    extension : list of strings\n",
       "        (default is ['oc', 'hds', 'ddn', 'cbc', 'ibo']).\n",
       "    unitnumber : list of ints\n",
       "        (default is [14, 51, 52, 53, 0]).\n",
       "    filenames : str or list of str\n",
       "        Filenames to use for the package and the head, drawdown, budget (not\n",
       "        used), and ibound output files. If filenames=None the package name\n",
       "        will be created using the model name and package extension and the\n",
       "        output file names will be created using the model name and extensions.\n",
       "        If a single string is passed the package will be set to the string and\n",
       "        output names will be created using the model name and head, drawdown,\n",
       "        budget, and ibound extensions. To define the names for all package\n",
       "        files (input and output) the length of the list of strings should be 5.\n",
       "        Default is None.\n",
       "\n",
       "    Attributes\n",
       "    ----------\n",
       "\n",
       "    Methods\n",
       "    -------\n",
       "\n",
       "    See Also\n",
       "    --------\n",
       "\n",
       "    Notes\n",
       "    -----\n",
       "    The \"words\" method for specifying output control is the only option\n",
       "    available.  Also, the \"compact\" budget should normally be used as it\n",
       "    produces files that are typically much smaller.  The compact budget form is\n",
       "    also a requirement for using the MODPATH particle tracking program.\n",
       "\n",
       "    Examples\n",
       "    --------\n",
       "\n",
       "    >>> import flopy\n",
       "    >>> m = flopy.modflow.Modflow()\n",
       "    >>> spd = {(0, 0): ['print head'],\n",
       "    ...   (0, 1): [],\n",
       "    ...   (0, 249): ['print head'],\n",
       "    ...   (0, 250): [],\n",
       "    ...   (0, 499): ['print head', 'save ibound'],\n",
       "    ...   (0, 500): [],\n",
       "    ...   (0, 749): ['print head', 'ddreference'],\n",
       "    ...   (0, 750): [],\n",
       "    ...   (0, 999): ['print head']}\n",
       "    >>> oc = flopy.modflow.ModflowOc(m, stress_period_data=spd, cboufm='(20i5)')\n",
       "\n",
       "     _name = OC\n",
       " _parent = MODFLOW 1 layer(s) 1 row(s) 201 column(s) 1 stress period(s) ('flopy.modflow.mf.Modflow)\n",
       " acceptable_dtypes (list, items = 3)\n",
       " allowDuplicates = False ('bool)\n",
       " cboufm = None ('NoneType)\n",
       " cddnfm = None ('NoneType)\n",
       " chedfm = None ('NoneType)\n",
       " compact = True ('bool)\n",
       " extra = \n",
       " file_name = gwexample.oc\n",
       " fn_path = temp/gwexample.oc ('str)\n",
       " ibouun = 0 ('int)\n",
       " iddnfm = 0 ('int)\n",
       " ihedfm = 0 ('int)\n",
       " iubud = 0 ('int)\n",
       " iuddn = 0 ('int)\n",
       " iuhead = 51 ('int)\n",
       " iuibnd = 0 ('int)\n",
       " label = LABEL ('str)\n",
       " savebud = False ('bool)\n",
       " saveddn = False ('bool)\n",
       " savehead = True ('bool)\n",
       " saveibnd = False ('bool)\n",
       " stress_period_data = {(0, 0): ['save head']} ('dict)\n",
       " unit_number (list, items = 5)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpm.ModflowOc(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally the MODFLOW input files are written (eight files for this model) and the model is run. This requires, of course, that MODFLOW is installed on your computer and FloPy can find the executable in your path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FloPy is using the following  executable to run the model: /Users/jdhughes/.local/bin/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: gwexample.nam \n",
      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/05/31 19:40:04\n",
      "\n",
      " Solving:  Stress period:     1    Time step:     1    Ground-Water Flow Eqn.\n",
      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/05/31 19:40:04\n",
      " Elapsed run time:  0.006 Seconds\n",
      "\n",
      "  Normal termination of simulation\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(True, [])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.write_input()\n",
    "model.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After MODFLOW has responded with the positive {\\tt Normal termination of simulation}, the calculated heads can be read from the binary output file. First a file object is created. As the modelname used for all MODFLOW files was specified as {\\tt gwexample} in step 1, the file with the heads is called {\\tt gwexample.hds}. FloPy includes functions to read data from the file object, including heads for specified layers or time steps, or head time series at individual cells. For this simple mode, all computed heads are read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpth = os.path.join(ws, 'gwexample.hds')\n",
    "hfile = fpu.HeadFile(fpth)\n",
    "h = hfile.get_data(totim=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The heads are now stored in the Python variable {\\tt h}. FloPy includes powerful plotting functions to plot the grid, boundary conditions, head, etc. This functionality is demonstrated later. For this simple one-dimensional example, a plot is created with the matplotlib package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n",
      "   could not parse units from None\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4XOWZ8P/vPapW77Ilq9qW3Iss94JpxvTQNjiNBBKS98fukmUhIW8Km32z2WRJI8nuEhICpJk0SDDFxlQD7t2yLRc1S1bvsnp5fn/MyAhbsiV5Zs6U+3NdujQ6Omfm1tEzc5/zVDHGoJRSSl0um9UBKKWU8g2aUJRSSjmFJhSllFJOoQlFKaWUU2hCUUop5RSaUJRSSjmFJhSllFJOoQlFKaWUU2hCUUop5RSBVgfgTAkJCSYzM9PqMJRSymvs3bu33hiT6Izn8qmEkpmZyZ49e6wOQymlvIaIlDnrubTKSymllFNoQlFKKeUUmlCUUko5hSYUpZRSTqEJRSmllFO4LKGISJqIvC0ix0TkiIg86NgeJyJbROSk43vsCMff49jnpIjc46o4lVJKOYcr71D6gH81xswAlgIPiMhM4FHgTWPMNOBNx88fISJxwGPAEmAx8NhIiUcppZRncNk4FGNMFVDleNwmIseAVOBWYI1jt+eAd4Cvnnf4dcAWY0wjgIhsAdYBGy72mjWtXTz5bhFx4cFkxoczY1IkkaFBTvqLlBq9urZujle3UdbYTktnLz19A4QGBWjZVJarbeuyl82GDlq7ep363G4Z2CgimcACYCeQ7Eg2GGOqRCRpmENSgfIhP1c4tg333PcD9wMET5zK914rPPc7m8DMlCiumZHMzfNSmJIY4YS/RqkLGWPYd7qZVw5V8VZhDaUNHRfd3yYwKyWaa2fay2ZWQribIlX+xhjD3rImXj5UxdvHaym7RNm8HGKMcdmTA4hIBPAu8B/GmBdEpNkYEzPk903GmNjzjnkECDHGfMfx8zeBDmPMDy/2Wvn5+Wbrth3UtnZT0tDOgdPNbCuqZ09ZE8bAiqnx3L96CqunJSAizv9jld/p6x/gxf1nePr9Egqr2wgOtLFyagJLs+OYnRpNZnw4ceHBBAXY6Oztp76tm5L6dvaXN/PBqXr2nbaXzZVTE/jiFdmsnKplUzlHb/8AL+47w6/eL+ZEzVmCA22smprA0ux4ZqdGkxEfRlx4MBOCA/caY/Kd8ZouTSgiEgS8DGw2xvzIse04sMZxdzIJeMcYk3vecesd+3zR8fMvHPtdtMorPz/fDDf1Sm1rF3/eW8HvdpRR1dLF4sw4vnXzTGanRjvl71T+acvRGr776jFK6tuZPjGSe5ZnctPcSWOqyqpp7eLPe8r57Y4yalq7WZJlL5uzUrRsqvExxrD5SA3/+doxyho6mDkpinuWZ3Dj3BQiQi6slBIRz08oYr/Meg5oNMZ8ecj2x4EGY8z3RORRIM4Y85Xzjo0D9gJ5jk37gIWDbSojGSmhDOru6+dPeyp44o2TNLZ385llmXx13XQmBAeM629U/qm2tYv/++Jh3jhWy7SkCB65LpdrZyZf1p1Fd18/f9xdzhNvnKSpo4fPrcjiketyCQ3SsqlGr7rFXjbfKqwlNzmSR67L5eoZSRctm96SUFYC7wGHgQHH5v+LvR3lT0A6cBq4yxjTKCL5wJeMMZ93HH+vY3+wV5c9c6nXvFRCGdTS2csPXz/Ob7aXMSUxnJ+uX6BXhGpUthyt4St/OUhnbz8PXZvD51ZkERTgvM6SLR29/NfmQn6/8zTTkiL46foFzJgU5bTnV75rU0EVX/3rYXr6BvjXtTl8dnkmgaMom16RUKww2oQy6INT9fzrnw7S1NHD9+6Yw20LJrswOuXN+gcMP95ygp+/fYo5qdH8+OPzmZrkuk4eW0/U8fCfD9LW1cf37pjDrfOH7ZOiFP0Dhsc3H+fJd4uYN9leNrPH0AFJE8oIxppQAOrPdvPA7/exs6SRf7pqKg9dm6ONouojunr7efD5/Ww+UsPdi9L49q2zCAl0fVVUbVsXD/x+H7tLm/jyNdN48OppWjbVR3T29PNPG/bxxrFaPrkknW/dPHPMZdOZCcXvp15JiAjhd59fwt2L0vjZW6f46l8PMTDgO0lWXZ62rl4+9audvH60hm/dNJPv3THXLckEICkylD98YSl35E3mJ2+c5P++WKBlU53T0tnLJ361gzcLa/n3W2fxH7fNcVvZHIlPLbA1XkEBNv7z9jkkRYbw07dOYRPhu7fNwWbTq0F/1t7dx2ef2c3B8mZ+vj6PG+dOcnsMQQE2fnDXXJKjQvifd4qwCXznY7P1TsXPtXX18tlndlFwpoX//eRC1s2eaHVIgCaUc0SEh9bmYoCfvXWKwADh/92qb1x/1dnTz73P7uZAeTM/X7+A6+e4P5kMEhEeuc5eNv/3nSKCAmw8dvNMLZt+qr27j3uf3c3hihb++5N5XDfLM5IJaEK5wEPX5tDTP8Av3i0mKMDGt27SN66/6ert5wu/2cPu0kZ+/PH5liaTQSLCV67LpadvgKffLyE40MbXrp+uZdPPdPb0c99zu9lb1sTP1ntWMgFNKBcQER5dN53ePsOvPyhhUnQo96+eYnVYyk2MMTzyl0N8UFTP43fO86jeVSLCN26cQW//AE9tLWZiVCj3rsyyOizlJgMDhof+dICdJY38+B/mW1IFeymaUIYhInzzphnUtHbxn68VMi0pkiunDzflmPI1//NOERsPVvKVdbncudDzupGLCP928yxqWrv4zitHmZoUweqcRKvDUm7w07dO8lpBNV+/YQYfW+A5FzpD+X0vr5GICI/fNZcZE6P45w37OVXbZnVIysVeP1LN45uP87H5KfyfKzz3rtRmE370D/PJSY7kH/+wj+K6s1aHpFzstcNV/OSNk9yel8rnV3nuXakmlIsICw7kl/fkExJk4/PP7aG5o8fqkJSLFFa38uU/HmDe5Gi+d8dcj2+bCA8J5JefyScwwF42WzqdOw258hxHKlt46E8HWZAew3dvm+PRZVMTyiWkxkzgF59eSGVzFw//+SC+NBBU2bV39/Gl3+4lIiSQX3w632vmz0qLC+N/P5lHeVMHX/mLlk1f1NbVy//53T6iJwTxi08v9PiyqQllFBZmxPHV66fzxrFafrfztNXhKCf79sYjlDV28NP1C5gYHWp1OGOyJDueR67LZfORGp7fXX7pA5RXeeylI1Q0dfDzTywgKdLzy6YmlFH63PJMVuck8p2Xj3KyRttTfMWrh6v4054K/r81U1iaHW91OOPy+ZX2dVT+feNRirQ9xWe8dLCSF/ad4R+vmkZ+ZpzV4YyKJpRRstmEH9w1l4iQQP75+QN09/VbHZK6TJXNnTz610PMS4vhy9fkWB3OuNlswg//YR6hQTYefH4/PX0Dlz5IebSKpg6+/uJh8tJj+OerplodzqhpQhmDpMhQ/uvOuRyrauXxTcetDkddhsE+/X0Dhic+Pt+pU9BbITkqlO/dMZeCM638cIuWTW/WP2B46I8HMQZ+8vEFo5qC3lN4T6Qe4uoZyXxqaTpPf1DC/tNNVoejxukPu06zo7iRx26eSaaPrOd+3ayJrF+cxi+3FnO4osXqcNQ4/W5HGbtKG/m3W2aRHh9mdThjogllHL66bjrJkaF87YXD9PZr9YK3qW3t4vubClkxNZ5/yE+zOhyn+toNM0iICOHRFw7Rp2XT61S3dPH45uOsmpbAHXmeOXjxYjShjENkaBDfvnUWhdVtPP1+idXhqDH69sajdPcN8J2PeXaf/vGICg3i326ZxZHKVp7dVmp1OGqMHnupgL6BAf7DS8umJpRxum7WRK6dmcxP3jhBeWOH1eGoUXqrsIZXDlfxT1dOJctHqrrOd/3siVw1PYkfvn6CiiYtm97i9SPVbD5Sw4NX53hdVdcgTSiX4du3zCJAhK//rUAHlXmBjp4+vvm3I0xLiuCLHjy1yuUSEf791lkAfOvvR7RseoGz3X089tIRpk+M9OipVS5FE8plSImZwMPX5bL1RB2bCqqtDkddwn+/fYozzZ189/Y5BAf6dtGfHBvGv67N4a3CWt44Vmt1OOoSfvbmSapbu/ju7XO8useh90buIT69NIOc5Ai+t6lQx6Z4sIqmDn75Xgm3LUhlkZcMErtc9yzPZEpiOP/56jHtPOLBTjd08MwHpdyZN5m89Firw7ksmlAuU2CAja/fOJOyhg5+u73M6nDUCB7ffBwBHrku1+pQ3CYowMbXb5xBcX07v9+hZdNTfX9TIQE24WEfKJsuSygi8msRqRWRgiHb5onIdhE5LCIbRSRqhGMfFJECETkiIl92VYzOckVOIlfkJPLEmydpbNcZiT3N/tNN/P1AJfevziYlZoLV4bjVlblJrJyawE/ePElLh85I7Gn2lDbyyuEqvnTFFJKjPH+urktx5R3Ks8C687b9CnjUGDMHeBF45PyDRGQ28AVgMTAPuElEprkwTqf4+o0zaO/u46dvnrQ6FDWEMYbvvHKMxMgQvuTDDfEjERG+fuMMWjp7+dlbWjY9ycCA4f+9coyJUaF8YbX3NsQP5bKEYozZCjSetzkX2Op4vAW4Y5hDZwA7jDEdxpg+4F3gNlfF6Sw5yZGsX5zOb3eUcapWJ+jzFK8ermZvWRMPr80hPMQ/FyidMSmKj+en8dz2Ukrr260ORzlsPFTJwfJmHr4ul7Bg3yib7m5DKQBucTy+CxhumHIBsFpE4kUkDLhhhP0AEJH7RWSPiOypq6tzesBj8S/X5jAhKIAfbNa5lDxBX/8Aj28uZPrESO5c6Fsj4sfqobU5BAXY+MHrWjY9QU/fAI9vPs6slChu99DlfMfD3QnlXuABEdkLRAIXNDgYY44B38d+B7MJOAj0jfSExpinjDH5xpj8xERr19ZOiAjh3pVZbDpSzZFKnUvJai/sO0NpQwcPr80lwOZ9o46dKSkylM+tyOSVw1Ucr9blF6z2l70VVDR18vB1udh8qGy6NaEYYwqNMWuNMQuBDUDRCPs9bYzJM8asxl5t5jWVv/etzCIqNJCfvOE1Ifuknr4BfvrWSeZOjubqGUlWh+MRvrAqm/DgQH7yxgmrQ/Fr3X39/Pytk+Slx7Amx9qLYGdza0IRkSTHdxvwDeDJS+yXDtyOPfl4hegJQXxhVTZbjtbojK8WGrwC/Jdrc7xyTiRXiAkL5t6VWbxWoHfQVvrT7nIqW7p46Npcnyubruw2vAHYDuSKSIWI3AesF5ETQCFQCTzj2DdFRF4dcvhfReQosBF4wBjjVfPEf3ZFJjFhQfxI16WwxOAV4AIfvAK8XPetzCJS76At09Xbz8/fPsXizDhWTPXOFUIvxmVdC4wx60f41RPD7FuJvfF98OdVrorLHSJDg7h/dTb/tek4+043ef3oV28zeAX4/Tvn+twV4OUavIP+0ZYTHK5oYc7kaKtD8it/2HmamtZufvLxBT5ZNnWkvIvcsyyTuPBgfrxF66vdafAKcFFmLCunJlgdjkf63IpMoifoHbS7dfX28z/vFLEsO55lU3zv7gQ0obhMeEgg96/O5r2T9dqW4kYv7DtDTWs3D16tbScjiQwN4gursnj7eB1HK1utDsdv/HlPOfVnu3nwGo8fpz1umlBc6BNL0okMCeQXW4ftzKacrH/A8Mv3ipmTGu2T9dPO9OmlmYQFB/CUlk236Osf4JfvlbAgPYYlWb47OakmFBeKCg3iE0vTefVwFWUNOkLZ1V4/Uk1JfTtfvCJb704uITosiPWL09l4qEoXiHOD1wqqOd3YwRdXT/HpsqkJxcXuW5FFoM3GL98rtjoUn2aM4cl3i8iID+P62ZOsDscr3LcyCwFdxtrFBstmdmI4a2cmWx2OS2lCcbGkqFBuz0vlz3sqqD/bbXU4PmtHcSMHK1r4wqpsvx8VP1opMRO4dX4qf9xdTpPOku0yH5xq4EhlK19cne1To+KHownFDb6wOpue/gGe21ZqdSg+6xdbi0iICObOhZOtDsWrfPGKbDp7+/mNruXjMk++W0RSZAgf86E5u0aiCcUNpiRGsHZmMr/ZXkZ794jTkqlxOlbVyjvH6/js8kxCgwKsDser5CRHcvX0JJ7bXkpnj6446mwFZ1p4/1Q9967MIiTQ98umJhQ3uX/1FFo6e/nrvgqrQ/E5v36/hAlBAXxqaYbVoXil+1dn09jew4v7z1gdis95+v0SIkIC+cSSdKtDcQtNKG6yMCOWeZOjeW5bKcYYq8PxGY3tPfz9YCW356USExZsdTheaXFWHDMnRWnZdLK6tm5ePlTJnQsnExUaZHU4bqEJxY3uWZ5JUV0775+qtzoUn/H87tP09A1wz/JMq0PxWiLCZ5dncrymjR3F56+Jp8Zrw67T9PYbPrPMf+6cNaG40Y1zJ5EQEayN807S1z/A77aXsXxKPDnJkVaH49VumZ9CbFiQlk0n6e0f4Pc7y1idk0h2YoTV4biNJhQ3CgkMYP3idN4srOV0gw4mu1xvHKuhsqVL706cIDQogI8vSuf1o9Wcae60Ohyvt6mgmprWbj673H/uTkATitt9ckkGASL8dkep1aF4vWc+KCU1ZgLXzPDtwWLu8qml9obj32oX4sv23LZSMuLDWJPjX4u7aUJxs4nRoVw3eyJ/3F1OR492IR6vY1Wt7Cxp5NPLMnQgo5NMjg3j2pnJPL/7NF292oV4vArOtLCnrIlPL83w+YGM59OEYoHPLs+ktauPv+2vtDoUr/Wb7WWEBNr4eH6a1aH4lHuWZ9Lc0ctLB7RsjtdvtpcyISiAu/ywbGpCsUB+RizTJ0byh11atTAeZ7v7+PuBM9wyL4XYcO0q7EzLsuOZlhTB73edtjoUr9TS2ctLByv52IJUoif4R1fhoTShWEBEWL84nYIzrRSc0bVSxuqlA5V09PRz92L/GCzmTiLCxxelcbC8mWNVulbKWL104AxdvQOsX+x/dyegCcUyH5ufSkigjQ16JThmf9x9mpzkCPLSY6wOxSfdnjeZ4AAbf9xdbnUoXsUYw4Zd5cycFMWcVP9cWlkTikWiw4K4cc4kx9W2Ns6P1tHKVg5WtHD3onSfXlfCSnHhwVw3eyIv7KvQxvkxKDjTytGqVtYvTvPbsqkJxUJ3L06nrbuPVw5VWR2K13h+92mCA23cnuf7M7daaf2iNFq7+nitQMvmaG3YfZrQIBu3+sGswiPRhGKhRZmxTEkM53mtWhiVrt5+Xtx/hutnT9R5u1xsaXY8GfFhPL9Ly+ZotHf38dKBSm6ck+I383YNx2UJRUR+LSK1IlIwZNs8EdkuIodFZKOIRI1w7L+IyBERKRCRDSIS6qo4rSQi3L0onb1lTZyoabM6HI/36uEq2rr6uHuRNsa7ms1mb5zfWdJIcd1Zq8PxeK8cquJsd5/fNsYPcuUdyrPAuvO2/Qp41BgzB3gReOT8g0QkFfhnIN8YMxsIAO52YZyWuj0vlaAA0SvBUXh+VzmZ8WEszY6zOhS/cGfeZAJsoo3zo7Bh92mmJkWwMCPW6lAs5bKEYozZCpw/dWkusNXxeAtwxwiHBwITRCQQCAN8dpRVfEQIa2dO5MX9FfT0DVgdjscqrjvLrtJGPq6N8W6TFBXK1dOT+Ou+Cnr7tWyO5GRNG/tPN3P3Iv9tjB/k7jaUAuAWx+O7gAvuD40xZ4AfAKeBKqDFGPP6SE8oIveLyB4R2VNXV+eCkF3vzoWTaero5Z3jtVaH4rFe3H8Gm6CN8W5258LJ1J/t4b2T3vnecocX9p8hwCZ+scTvpbg7odwLPCAie4FIoOf8HUQkFrgVyAJSgHAR+dRIT2iMecoYk2+MyU9MTHRR2K61aloCCRHBvLBPV8wbzsCA4cX9Z1gxNYHkKJ9sTvNYa3KTiA0L0rI5gv4Bw9/2n+GKnEQSIkKsDsdybk0oxphCY8xaY8xCYANQNMxu1wAlxpg6Y0wv8AKw3J1xultggI1b56fyZmENzR0X5Fi/t7u0kYqmTu7Im2x1KH4nONDGLfNSeP1oDS2dvVaH43F2FDdQ1dKld84Obk0oIpLk+G4DvgE8Ocxup4GlIhIm9grJq4Fj7ovSGrctSKW33/Cyjkm5wAv7zhAeHMDaWTpNvRVuz5tMT98Arx3Wsnm+F/adITIkUJdQcHBlt+ENwHYgV0QqROQ+YL2InAAKsTe0P+PYN0VEXgUwxuwE/gLsAw47YnzKVXF6ilkpUeQmR/LCvgqrQ/EoXb39vHq4inWzJxEWHGh1OH5p7uRoshPDtdrrPB099oGfN86dRGhQgNXheASXvUONMetH+NUTw+xbCdww5OfHgMdcFJpHEhFuz0vlP18rpKS+nayEcKtD8ghbjtbQ1t3HHVqlYBkR4Y68yTy++TjljR2kxYVZHZJH2Hykmo6efm7XqthzdKS8B7l1fioi9h5Nyu6FfRWkRIeyNDve6lD82mAPJi2bH3ph3xkmx04g38/HngylCcWDTIwOZeXUBF7YV8HAgLE6HMvVtnWx9WQ9H1uQ6ncr33ma1JgJLMuO54V9FRijZbO6pYsPTtVzu5bNj9CE4mFuW5BKRVMne083WR2K5TYerKJ/wGgPGg9xW14qpQ0d7C9vtjoUy7108AwDBm7T6q6P0ITiYa6bNZHQIBsbD/rs5ACjtvFgJTMnRTE1KdLqUBSwbvZEggO0bIL9Ymfe5Ght6zyPJhQPEx4SyFXTk3j1cBV9fjzdRXljBwfKm7l5XorVoSiHqNAg1uQm8soh+52jvyqtb+fwmRZumqtl83yaUDzQzXNTqD/bw86S86dC8x+D43FumjvJ4kjUUDfNS6G2rZvdpf5cNu13aDdq2byAJhQPdOX0JMKDA84VXH/08qFK5qfFaBdVD3PNjCQmBPl72awiPyOWlJgJVoficTSheKDQoACunZnMawXVfjnLa3HdWY5UturdiQcKCw7kqhlJvHa42i+rZE/WtFFY3aZVsSPQhOKhbp6XQnNHL++fqrc6FLd7+VAVImgdtYe6eW4KDe09bC9usDoUt9t4qAqbwPVzJlodikfShOKhVk1LJCo00C971Gw8WMmijDgmRuvMwp5oTW4iESH+VzaNMbx8sJKl2fEkRWrZHI4mFA8VHGhj3eyJbDlSQ1dvv9XhuM3x6jZO1p7l5nla3eWpQoMCWDszmU0F1X61KNzRqlaK69u1uusiNKF4sJvmptDW3ce7J/xncaONByuxCaybrQnFk900bxKtXX1+tfDWxoNVBNqEdbO0umskmlA82PIp8cSFB/OKn0xpb4zhlcNVLJsST2KkLlbkyVZOTSR6QpCflc1KVkxNIDY82OpwPNaoZxsWkXxgFfZVFDuxL+f7hjHGfzuku1hggI21M5N5+VAV3X39hAT69hTZJ2rOUlLfzn0rs6wORV1CcKCNa2Yks+WovdorONC3r02PVrVS3tjJP1451epQPNolS4GIfFZE9gFfAyYAx4FaYCWwRUSeE5F014bpv66bPZGz3X1sO+X7PWo2FVQjAmtn6mJF3mDd7Im0dvWxww96e20uqMYm6EJalzCaO5RwYIUxpnO4X4rIfGAa9pUWlZMtnxJPZEggmwqquXJ6ktXhuNSmI9UsTI8lSdeN9wqrpiUQFhzApiPVrM5JtDocl9p0pJrFWXHE67rxF3XJOxRjzH+PlEwcvz9gjHnTuWGpQSGBAVw1I4ktx2p8eiBZWUM7x6paWTdbGzy9RWhQAFfmJvH6kRqfnturqO4sJ2rOamP8KIy64lNEskTkRyLygoi8NPjlyuCU3fWzJ9LY3sPuUt+d0n7zkWrAPtuy8h7rZk+k/mw3+3x4uYVzZVMvdi5pLEsA/w14GtgI+O6lsgdanZNIaJCNzUeqWTbFN1cu3FRQzezUKJ27y8tcOT2J4AAbmwqqWZQZZ3U4LrG5oJr5aTFMita5uy5lLF0zuowxPzXGvG2MeXfwy2WRqXPCggO5IieRTQXVPrmSY01rF/tON2uVgheKCAlk1bQENhVU++RKjmeaOzlY0aJVsaM0loTyhIg8JiLLRCRv8MtlkamPWDd7ItWtXRw602J1KE73uqNKQd+03um62RM509zJkcpWq0Nxute1KnZMxlLlNQf4NHAVH1Z5GcfPysWump5MoE14raCK+WkxVofjVK8VVDMlMVxXZvRS18xIJsBRNmenRlsdjlO9VlDN9ImRujLjKI3lDuU2INsYc4Ux5krH14jJRER+LSK1IlIwZNs8EdkuIodFZKOIRA1zXK6IHBjy1SoiXx7bn+V7oicEsXxqApt9rGqhqd2+kJjenXivuPBglmTFsamg2upQnKr+rH0hMb07Gb2xJJSDwFgujZ8F1p237VfAo8aYOcCLwCPnH2SMOW6MmW+MmQ8sBDoc+/q9tTOTKW3ooKjurNWhOM3bx2vpHzBcO1PftN5s7cxkiuraKalvtzoUp3mrsBZj4FodaDtqY0koyUChiGweTbdhY8xW4PxpWXKBrY7HW4A7LvGaVwNFxpiyMcTps66eYR/Y+MaxWosjcZ43j9WSGBnCXB+rKvE3VztGkL95rMbiSJznzWM1pESHMivlgooUNYKxtKE85oTXKwBuAf4O3AWkXWL/u4ENTnhdnzApegKzUqJ481gNX7piitXhXLaevgHePVHHTXMnYbOJ1eGoy5AWF0ZuciRvHKvh86uyrQ7nsnX19vPeyXpuz0tFRMvmaI1mLi8BGNpV+PxuwzL6M34v8ICI7AUigZ6LvG4w9uTz50vEd7+I7BGRPXV1vj+V9tUzktlb1kRj+4inzmvsLGngbHffuatb5d2unpHE7tImWjp6rQ7lsm0vbqCjp1/L5hiNpsrrbRH5p/MngBSRYBG5SkSeA+4ZzYsZYwqNMWuNMQux33kUXWT364F9xpiL3kMbY54yxuQbY/ITE317PiGAa2ckM2Dg7ULvr/Z681gtIYE2Vk5NsDoU5QTXzEymf8DwzglfKJs1hAUHsCzbNwcSu8poEso6oB/YICKVInJUREqAk8B64MfGmGdH82IikuT4bgO+ATx5kd3Xo9VdF5idGkVyVAhvFnp3XbUxhjeO1bByagITgn17Wn5/MX9yDAkRwV7fxme1zW8NAAAej0lEQVSM4a1jtayalkBokJbNsRjN5JBdxpj/McasADKwN5QvMMZkGGO+YIw5MNxxIrIB2A7kikiFiNwHrBeRE0AhUAk849g3RUReHXJsGHAt8MJl/n0+R0S4anoyW0/U093nvUsDn6g5S0VTp1Yp+BCbTbgyN4l3jtfS68UTmR6taqWypUvL5jiMpg0lbvALe7tHN2Absm1Yxpj1xphJxpggY8xkY8zTxpgnjDE5jq9HjWNAhTGm0hhzw5BjO4wx8cYY3xsW7gTXzEjibHcfO4u9d22zNxy9gQZ7rinfcPWMZNq6+thd4sVl82gtInCVjy8X4QqjqfLaC+xxfK8DTmCv7qpzbFNutmJqAqFBNq/uovnGsRrmTo4mWdc+8SmrpiUQHGjz6mqvNwtrWJAWQ4KufTJmo6nyyjLGZAObgZuNMQnGmHjgJrRKyhKhQQGsnJrIG8dqvXLUfF1bNwfKm7l6ulYp+JrwkECWT4nnzcIaryybNa1dHKpo0equcRrLwMZFxphz7RzGmNeAK5wfkhqNa2Ykcaa5k8LqNqtDGbO3j9tHIGt1l2+6ekYyZQ0dnKr1vhkd3nL0ntSlfsdnLAmlXkS+ISKZIpIhIl8HfH8xaQ81WL/79nHvq1p4u7CW5KgQHYHso6724rL5VmEtqTETyEmOsDoUrzSWhLIeSMQ+r9bfgCTHNmWBpKhQZk6K4p3j3jWYs7d/gPdP1rMmJ0lHIPuoFMcHsreVzZ6+AbadqmdNbqKWzXEa9dQrxphG4EEXxqLGaE1uIr/YWkxrVy9RoUFWhzMq+8qaaOvuY02u7w9C9WdrcpN45oMSznb3EREylhmerLOntJH2nn7W5GpV7HiNZU35RBF5XEReFZG3Br9cGZy6uDW5SfQPGD44WW91KKP2zok6Am3Cimk6Ot6XrclJpLffsO2Ud5XN4AAby310mW13GEuV1++xD0jMAr4NlAK7XRCTGqW89BgiQwO9qmrhneN15GXEes0dlRqf/Mw4woMDeOeEN5XNWhZlxRLuJXdUnmgsCSXeGPM00OuYGPJeYKmL4lKjEBhgY9W0BN49UecVXTRrWrs4VtWq1V1+IDjQxvKpCbx73DvKZmVzJydqzrImR6u7LsdYEsrgFKJVInKjiCwAJrsgJjUGV+QkUt3a5RXdh9913Enpm9Y/rMlN5Exzp1d0Hx68y79CL3Yuy1gSyndEJBr4V+Bh7Ksv/otLolKjdoXjw9kbqr3eOWHvLjxjkq4d7w8GG7e9omweryUlOpRpSdpd+HKMOqEYY142xrQYYwoc68kvNMaMuGKjco+J0aFMnxjJOx7e57+vf4D3TtZzRY52yfQXqTETmJYU4fHT2ff0DfDBqXquyNWu7JdrLL28ckTkTREpcPw8V0S+4brQ1GityU1ib1kTbV2eu7DRvtPNtHX1aZdMP7MmN5HdJU20d/dZHcqI9pQNdhfW6q7LNZYqr18CX8PRlmKMOYR9iV5lsTW5ifQNGD7w4C6a7xyvJcAmrNDFtPzKmtwkevoH2FbkuZNqvHu8jqAALZvOMJaEEmaM2XXeNs+97PAjCzNiiQwJ5F0P7qL57ok6FqbHEj1Buwv7k/zMWMKCA3jXg6u93j1RR35GnNcMwPRkY53LawpgAETkTqDKJVGpMQkKsLFsSjzvnaz3yC6aDWe7OVLZyuocvQL0NyGB9mV03/fQwbe1bfYekqtztLrLGcaSUB4AfgFMF5EzwJeBL7kkKjVmq6YlUNHUSVlDh9WhXOADR3XHymn6pvVHK6clUNrQQXmjB5ZNRzXxKp25wSnGklDOYF+y9z+A54EtwD2uCEqN3eCH9XsnPa/a6/2TdURPCGJOarTVoSgLDH5Yv++BbXzvnagnLjyYmZN05mtnGEtC+TtwM/ZG+UrgLNDuiqDU2GXGh5EaM4H3PKxqwRjD+yfrWT4lngCbdsn0R1MSI5gYFepx1V7GGN4/ZS+bNi2bTjGWVqjJxph1LotEXRYRYdW0BF45VEVf/wCBAWO5VnCd4vp2Klu6eOAqrVLwVyLCymkJvHGshv4B4zEXFidqzlLb1q3VXU40lk+dbSIyx2WRqMu2cloCbd19HKxosTqUcwavSldN1fYTf7ZqWgLNHb0cqfScsjlYPaxte85zyYQiIodF5BCwEtgnIsdF5NCQ7cpDrJiSgAgeVbXw3sl60uPCSI8PszoUZaHBMR6eVCX7/ql6shPCSY2ZYHUoPmM0dyg3YW87uR6YCqx1/Dy4fVgi8msRqR0cWe/YNk9EtjuS0UYRGbYlTERiROQvIlIoIsdEZNlY/ih/FRsezOyUaN4/5RkN8739A+wobmClVin4vYSIEGZMivKYTiPdff3sLG7Usulkl0woxpiyi31d5NBngfPbXH4FPGqMmYN9KeFHRjj2CWCTMWY6MA84dsm/RAH2aq/9p5s56wFTXRwot8exSkcgK+zVXnvLmujosb5s7itrprO3n5VaNp3KZS23xpitQON5m3OBrY7HW4A7zj/OcdeyGnja8Tw9xphmV8Xpa1ZNTaBvwLDDA6a6eO9kPTaB5VP0Tatg5dQEevsNO0vO/1hwv/dP1RFgE5bq6oxO5e6uQAXALY7HdwFpw+yTDdQBz4jIfhH5lYiEuytAb7cwM5bQIJtH9Pl//2QdcybHEB2m060oWJwVR3CgzSPa+N4/Wc/8tBhdOdTJ3J1Q7gUeEJG9QCTQM8w+gUAe8L/GmAXYx7o8OtITisj9IrJHRPbU1XlG/ayVQgIDWJwVz1aL66pbu3o5WNGi1V3qnNCgABZlxlqeUJo7ejh0pkWru1zArQnFGFNojFlrjFkIbACKhtmtAqgwxux0/PwX7AlmpOd8yhiTb4zJT0zU7n9gr/YqrmunuqXLshh2FTfSP2B0Blf1ESumJnC8po26tm7LYthR3IAxaIO8C7g1oYhIkuO7DfgG8OT5+xhjqoFyEcl1bLoaOOq2IH3AMke98PZi664EtxU1EBJoY0F6jGUxKM8z2J62o9i6Nr5tRQ1MCApg3mQtm87msoQiIhuA7UCuiFSIyH3AehE5ARRin77lGce+KSLy6pDD/wn4vWOcy3zgu66K0xfNmBRF9IQgtlvYML+9uIH8zFhCgwIsi0F5ntkpUUSGBLLdwoSyvaiBRY72HOVcLlsAwBizfoRfPTHMvpXADUN+PgDkuyg0nxdgE5ZkxVm2qFFjew/Hqlp5eG2OJa+vPFdggI3FWXGWXezUtnVxsvYst+dNtuT1fZ2maB+1fEo8FU2dlkwZvtNx9blMuwurYSybEk9JfTtVLZ1uf+0dxfYuy8u1u7BLaELxUcsdjeFWXAluK2ogLDiAuZN1unp1oXNtfBaUze1F9USGBjIrRaerdwVNKD5qWlIECRHBltRVbyuqZ3FWHEEeMuOx8iwzJkYRExZkSZXs9qIGlmTFecxs3L5Gz6qPEhGWZsezrci9ywLXtnZRVNfOsmytUlDDs9mEpVnxbC9qcGvZrGzupLShQ6tiXUgTig9bNiWemtZuSurdtw7a4B2RTreiLmb51HjONHdS3ui+dpTBKja92HEdTSg+bPBD3Z1VC9tONRAVGshMraNWFzH4ob6tyH1jpbYVNRAbFsT0iZFue01/ownFh2XGhzEpOtSt7SjbixtYkq3L/aqLm5oUQUJEiNvKpjGGHcUNLNPlfl1KE4oPExGWZcezw0111RVNHZxu7NAqBXVJIsKyKfFsc1PZPN3YwZnmTi2bLqYJxcctmxJPQ3sPx2vaXP5ag3XUy6fqm1Zd2vIp8dS1dVNUd9blr3Wu/UTb9lxKE4qPW+q4IttZ7Po1KHYUNxIbFkROktZRq0sbLJs73FI2G0iICGFKoq6E4UqaUHxcWlwYqTET2OWGRY12lTawOCtO66jVqGTGh5EUGeLysmmMfVGvJdlxiGjZdCVNKH5gcVYcO0tcW1dd6egCuiRLq7vU6IiIW8pmRVMnVS1dLMmKc9lrKDtNKH5gSVYc9Wd7KHbheJTBq8zF+qZVY7Ak2z5W6rQL55wbXHJYL3ZcTxOKHxj8kHdlO8rOkkYiQwOZMUnHn6jRG7xrcOU68zuLG4gJC2JaUoTLXkPZaULxA1kJ4SREhLCrxHV9/neWNLAoM07Hn6gxmZYUQVx4sEsvdnaVNrI4U9v23EETih8QEZZkx7GzpNElddV1bd0U17VrdZcaMxFhUWYsu0pdc7FT3dJFWUOHlk030YTiJ5ZkxVHV0kVFk/PnTtpdOlhHrW9aNXZLsuIpb+ykstn5ZXOn465c20/cQxOKn1jswrrqncX29U9mp+r6J2rsBsumK7oP7yppJCJE55ZzF00ofiInKZKYsKBzqyk6086SRhZmxOr6J2pcZkyKIjI00DUXOyWN5GfGatuem+gngJ+w2YRFmXHsKnXum7a5wz6ty+JMre5S4xMwWDad3Gmk/mw3p2rPavuJG2lC8SNLsuIoa+iguqXLac+5u7QJY+zjCZQar8VZcRTVtVPX1u2059yt40/cThOKHxl8Y+104pXgrpIGggNtun68uiyDHTp2O/EOemdJI6FBNuZo257buCyhiMivRaRWRAqGbJsnIttF5LCIbBSRYVvKRKTUsc8BEdnjqhj9zYxJkUSEOLeueldJI/PTYggNCnDacyr/Mzs1mrDgAKe28e0qaSQvPZbgQL1udhdXnulngXXnbfsV8KgxZg7wIvDIRY6/0hgz3xiT76L4/E5ggI28jFj2ljY55fk6evooqGzV9hN12YICbCxIj2FPmXPKZltXL4XVrSzSsulWLksoxpitwPmXwrnAVsfjLcAdrnp9Nbz8jFiO17TR0tF72c91sLyF/gHDwoxYJ0Sm/N3CjDiOVbVytrvvsp/rQHkzAwYtm27m7nvBAuAWx+O7gLQR9jPA6yKyV0Tud0tkfiI/0/4G23f68q8E95bZrxfy0vVNqy7fosxYBgzsd0LZ3FPahE1gQXqMEyJTo+XuhHIv8ICI7AUigZ4R9lthjMkDrnfsv3qkJxSR+0Vkj4jsqaurc37EPmZ+WgwBNmFP2eW3o+wpayInOYLosCAnRKb83YL0WGxi7zl4ufaWNZE7MYrIUC2b7uTWhGKMKTTGrDXGLAQ2AEUj7Ffp+F6Lva1l8UWe8yljTL4xJj8xMdEVYfuUsOBAZqdEXfabdmDAsK+siYUZWketnCMixD5b9d7LvNjp6x9g/+km8rW6y+3cmlBEJMnx3QZ8A3hymH3CRSRy8DGwFntVmXKS/Mw4DpY309M3MO7nOFl7ltauPq2jVk6VnxHL/tPN9PWPv2wWVrfR3tN/rnpXuY8ruw1vALYDuSJSISL3AetF5ARQCFQCzzj2TRGRVx2HJgPvi8hBYBfwijFmk6vi9Ef5GbF09w1QUNky7ufY6+iNo1eBypnyM+Po6OnnWFXbuJ9jsH1Q2/bcL9BVT2yMWT/Cr54YZt9K4AbH42JgnqviUrDQceW2t7Rp3G+6PWWNJEQEkxEf5szQlJ8bvKvYXdrInHEOlt1T2kRyVAiTYyc4MzQ1Cjrixw8lRYaSER92WaOS95Y1sTAjFhGddE85z6ToCaTGTDh3Bzwee8uayM+I07JpAU0ofio/I469ZU3jWnCrrq2bsoYObT9RLpGfGcvu0vEtBlfV0smZ5k4tmxbRhOKn8jNjaWjvoaS+fczHDl49ag8v5Qr5GbHUtnWPazG4D8umJhQraELxU4scddXjmepib1kjwYE2ZqfqokXK+fIzxz9R5J7SJiYEBeiCWhbRhOKnshMiiAkLGte8XnvKmpibGk1IoE4IqZwvJzmSyJDAcV7sNDEvLVoXe7OInnU/ZbMJ+Rmx7B7jILKu3n4KzrSc6ymmlLMF2IS8jFj2jPEOpb27j6NVreRrVaxlNKH4sbyMWIrr2mlqH2kGnAsdPtNCb7/RN61yqYUZsZyoOUtL5+gnMT1Y0WyfrFQvdiyjCcWPDY5B2V8++qqFPaXa6Klcb7BsHihvHvUxg9W3eWlaNq2iCcWPzZ0cTYBN2Fc2hjdtWSPZCeHEhQe7MDLl7+alRSMC+8bQjqKTlVpPE4ofCwsOZMakyFHfoRhjzg1oVMqVIkODyE2OZP8o71AGBgz7TutkpVbThOLn8tJjOXDaXvd8KUV17TR19Oqke8otFqTHsv90EwOjKJsnatto6+rTueUspgnFzy1Ij6G9p58TNZeejG+fDmhUbpSXHkNbVx9FdWcvue+5yUr1YsdSmlD83GDj52hWcNxT1khMWBDZCeGuDkspFoyhbO4tbSIhIpj0OJ2s1EqaUPxcelwY8eHBo2qY31vWxML0WGw2nXRPuV52QjjRE4JGVzZP62SlnkATip8TEXtd9SUa5ls6eimqa9c1upXb2GzCgvSYS5bNhrP2yUoX6PonltOEoliQHkNxXTvNHSMPcDxYYb9KnK99/JUb5aXHOlYHHXmA44dlUy92rKYJRX04wPH0yFULB8qbEYG5aeNb9Eip8chLj8UYOHCxsnm6GZvAnFQtm1bThKKYl+YY4HiRxs+D5c1MSYwgKlQHjSn3OTfA8SJl80BFCznJkYSHuGwBWjVKmlAUYcGBTJ8YOeIdijGGA+XNWqWg3C4yNIicpIuXzYPlzdq25yE0oSjA3o5yoHz4AY4VTZ00tPdoQlGWyMuIGXGAY0l9Oy2dvVo2PYQmFAXY66rPdvdxsvbCAY6D01/om1ZZYUF6LK1dfRTXXzjAcXDyyHlaNj2CJhQFXLxh/sDpZkKDbEyfGOnusJQaMvh2mLJZ3kx4cADTkrRsegKXJRQR+bWI1IpIwZBt80Rku4gcFpGNIjLiOp0iEiAi+0XkZVfFqD6UER9GTFgQB4eZjO9AeRNzUqMJ1FXwlAWyE8KJDAkcoWw2M8cxa7aynis/IZ4F1p237VfAo8aYOcCLwCMXOf5B4JhrQlPnExHmTY65YP2Jnr4BCipbtbpLWcZmE+amRZ8bbzKoq7efY1WtOjbKg7gsoRhjtgLnr+GZC2x1PN4C3DHcsSIyGbgRewJSbjIvLYYTNW20d/ed21ZY3UpP34C+aZWl5qfFUFjVRldv/7ltRypb6e03erHjQdxdh1EA3OJ4fBeQNsJ+PwG+Agy4IyhlNz8tmgEDBWdazm0bvGOZr90ylYXmTY6hb8BwpPLCsqldhj2HuxPKvcADIrIXiAQumOtDRG4Cao0xe0fzhCJyv4jsEZE9dXV1zo3Wz8ybbH9jDq1aOHC6mcTIEFKiQ60KS6lzdyEHyj+aUCZFh5IcpWXTU7g1oRhjCo0xa40xC4ENQNEwu60AbhGRUuB54CoR+d1FnvMpY0y+MSY/MTHRJXH7i/iIENLiJnykHWVwQKPO4qqslBQVSkp06Eca5g+UN2l1l4dxa0IRkSTHdxvwDeDJ8/cxxnzNGDPZGJMJ3A28ZYz5lDvj9GfzJsdw0HEV2NLRS3F9u75plUeYl/Zhp5GGs92UN3Zq2fQwruw2vAHYDuSKSIWI3AesF5ETQCFQCTzj2DdFRF51VSxq9OanxXCmuZPati4OOKq+FuibVnmAeWkxnG7soLG958O2PS2bHsVls6kZY9aP8Ksnhtm3ErhhmO3vAO84NTB1UYNv0IPlLRytbEUE5kzWWVyV9c6VzYpmDpQ3E2ATLZseRqfnVB8xK8U+SOxgeTNHKluYlhRBpM4wrDzAnNRobGLvKHKgvJmc5EjCgvUjzJPo0Gf1EROCA5g+MZID5c06w7DyKOEhgUxLimS/lk2PpQlFXWBeWgw7ihto6ujVAY3Ko8xLi2bbqXrauvq0bc8DaUJRF5ifZh9EBvY3sFKeYn5a7LmyqYNtPY8mFHWBwaqECUEB5CbrLK7Kcwxe4ESEBDIlMcLiaNT5tEVLXWBKYgThwQHMStEZhpVnyU2OJDTIxpxUnWHYE2lCURcIsAnfvGkmKTETrA5FqY8IDLDxrZtmkREfZnUoahiaUNSw7l6cbnUISg3rE0u0bHoqrc9QSinlFJpQlFJKOYUmFKWUUk6hCUUppZRTaEJRSinlFJpQlFJKOYUmFKWUUk6hCUUppZRTiDHG6hicRkTagONWx3EJCUC91UGMgsbpXN4QpzfECBqns+UaY5wyaZ+vjZQ/bozJtzqIixGRPZ4eI2iczuYNcXpDjKBxOpuI7HHWc2mVl1JKKafQhKKUUsopfC2hPGV1AKPgDTGCxuls3hCnN8QIGqezOS1On2qUV0opZR1fu0NRSillEZ9IKCKyTkSOi8gpEXnU4ljSRORtETkmIkdE5EHH9n8TkTMicsDxdcOQY77miP24iFznpjhLReSwI5Y9jm1xIrJFRE46vsc6touI/NQR4yERyXNTjLlDztcBEWkVkS97wrkUkV+LSK2IFAzZNubzJyL3OPY/KSL3uCnOx0Wk0BHLiyIS49ieKSKdQ87rk0OOWegoL6ccf4tTl0scIc4x/59d+VkwQox/HBJfqYgccGy38lyO9Bnk+vJpjPHqLyAAKAKygWDgIDDTwngmAXmOx5HACWAm8G/Aw8PsP9MRcwiQ5fhbAtwQZymQcN62/wIedTx+FPi+4/ENwGuAAEuBnRb9n6uBDE84l8BqIA8oGO/5A+KAYsf3WMfjWDfEuRYIdDz+/pA4M4fud97z7AKWOf6G14Dr3RDnmP7Prv4sGC7G837/Q+BbHnAuR/oMcnn59IU7lMXAKWNMsTGmB3geuNWqYIwxVcaYfY7HbcAxIPUih9wKPG+M6TbGlACnsP9NVrgVeM7x+DngY0O2/8bY7QBiRGSSm2O7GigyxpRdZB+3nUtjzFagcZjXH8v5uw7YYoxpNMY0AVuAda6O0xjzujGmz/HjDmDyxZ7DEWuUMWa7sX/S/IYP/zaXxXkRI/2fXfpZcLEYHXcZ/wBsuNhzuOlcjvQZ5PLy6QsJJRUoH/JzBRf/AHcbEckEFgA7HZv+0XFL+evB202si98Ar4vIXhG537Et2RhTBfZCCSRZHONQd/PRN6snnctBYz1/VscLcC/2q9NBWSKyX0TeFZFVjm2pjtgGuTPOsfyfrTyfq4AaY8zJIdssP5fnfQa5vHz6QkIZrv7R8q5rIhIB/BX4sjGmFfhfYAowH6jCfnsM1sW/whiTB1wPPCAiqy+yr6XnWESCgVuAPzs2edq5vJSR4rL6vH4d6AN+79hUBaQbYxYADwF/EJEorItzrP9nK8/nej56wWP5uRzmM2jEXUeIacyx+kJCqQDShvw8Gai0KBYARCQI+z/y98aYFwCMMTXGmH5jzADwSz6sirEkfmNMpeN7LfCiI56awaosx/daK2Mc4npgnzGmBjzvXA4x1vNnWbyOBtabgE86ql5wVCE1OB7vxd4ekeOIc2i1mLvK6Fj/z5acTxEJBG4H/ji4zepzOdxnEG4on76QUHYD00Qky3ElezfwklXBOOpSnwaOGWN+NGT70DaH24DBniIvAXeLSIiIZAHTsDfauTLGcBGJHHyMvZG2wBHLYE+Oe4C/D4nxM47eIEuBlsFbZzf5yNWfJ53L84z1/G0G1opIrKM6Z61jm0uJyDrgq8AtxpiOIdsTRSTA8Tgb+/krdsTaJiJLHeX7M0P+NlfGOdb/s1WfBdcAhcaYc1VZVp7LkT6DcEf5dGbvAqu+sPdSOIH9KuDrFseyEvtt4SHggOPrBuC3wGHH9peASUOO+boj9uM4ucfHCDFmY+8BcxA4MnjOgHjgTeCk43ucY7sA/+2I8TCQ78bzGQY0ANFDtll+LrEnuCqgF/uV3H3jOX/Y2zBOOb4+56Y4T2GvGx8sn0869r3DUR4OAvuAm4c8Tz72D/Qi4Oc4BkW7OM4x/59d+VkwXIyO7c8CXzpvXyvP5UifQS4vnzpSXimllFP4QpWXUkopD6AJRSmllFNoQlFKKeUUmlCUUko5hSYUpZRSTqEJRSmllFNoQlHKyURkkoi8PMZjfiAiV7kqJqXcQROKUs73EPapQsbiZ9inFFfKa2lCUWqcRGSRYybcUMd0NkdEZDb2UdKbHPt8VkT+JiIbRaRERP5RRB5yzEK7Q0TiAIx9Wv54EZlo4Z+k1GXRhKLUOBljdmOfEuQ72Bcv+h3QDjQZY7qH7Dob+AT2yQ3/A+gw9llot2Ofy2nQPmCFG0JXyiUCrQ5AKS/379gnJewC/hlYAtSdt8/bxr7QUZuItAAbHdsPA3OH7FcLpLg2XKVcRxOKUpcnDogAgoBQoNPxfaihdysDQ34e4KPvwcHjlfJKWuWl1OV5Cvgm9kWqvo99ptvMcT5XDh9O0a6U19E7FKXGSUQ+A/QZY/7gWPtiG/YqryIRmWqMOTWG5woCpgJ7XBOtUq6n09cr5WQichuw0BjzjTEek2eM+abrIlPKtfQORSknM8a8KCLxYzwskA/XTFfKK+kdilJKKafQRnmllFJOoQlFKaWUU2hCUUop5RSaUJRSSjmFJhSllFJO8f8DmzuVveBhAWAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "ax = plt.subplot(111)\n",
    "x = model.dis.sr.xcentergrid[0]\n",
    "ax.plot(x,h[0,0,:])\n",
    "ax.set_xlim(0,x.max())\n",
    "ax.set_xlabel(\"x(m)\")\n",
    "ax.set_ylabel(\"head(m)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 1
}