{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8d50e013",
   "metadata": {},
   "source": [
    "# **Subset CONUS and run ParFlow-CLM**\n",
    "\n",
    "To launch this notebook interactively in a Jupyter notebook-like browser interface, please click the \"Launch Binder\" button below. Note that Binder may take several minutes to launch.\n",
    "\n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/hydroframe/subsettools-binder/HEAD?labpath=subsettools%2Fconus1_subsetting_transient.ipynb)\n",
    "\n",
    "This notebook walks through an example of subsetting a HUC8 from the CONUS1 domain. This example will subset everything needed to do a transient run with ParFlow-CLM. \n",
    "This is includes all hydrogeologic datasets and climate forcing data from NLDAS2 (see additional DOI: [10.1002/2016GL069964](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016GL069964)). All of the data is written to a folder for the specified days to run. This example uses the template runscript conus1_pfclm_transient_solid.yaml and edits it to correspond with the domain subset. It also sets-up and performs the designed simulation. The result will be model output pressure and saturation pfbs according to the days specified.\n",
    "\n",
    "### This notebook has two principal sections: \n",
    "1. Subset all static inputs and climate forcings from a CONUS run stored in Hydrodata \n",
    "2. Load and alter a reference run to set up and perform your ParFlow-CLM subset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "134427aa",
   "metadata": {},
   "source": [
    "### Import the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "84c7594b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "from parflow import Run\n",
    "from parflow.tools.io import read_pfb, read_clm\n",
    "from parflow.tools.fs import mkdir\n",
    "from parflow.tools.settings import set_working_directory\n",
    "import subsettools as st\n",
    "import hf_hydrodata as hf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "93b0c616-c6c3-4b0e-ab55-b8caff0594cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You need to register on https://hydrogen.princeton.edu/pin before you can use the hydrodata utilities\n",
    "hf.register_api_pin(\"your_email\", \"your_pin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec636a8a",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1. Define variables to access datasets in Hydrodata to subset and define write paths"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbb358b8",
   "metadata": {},
   "source": [
    "#### We will be testing with the Upper Verde watershed for this example\n",
    "- HUC: 15060202\n",
    "- Size: 6496 km^2 (ni = 112, nj = 90)   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "897c0d79",
   "metadata": {},
   "source": [
    "#### Set your variables to specify which static and climate forcing data you would like to subset in Hydrodata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2fb131aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "runname = \"conus1_upper_verde\"\n",
    "\n",
    "# provide a way to create a subset from the conus domain (huc, lat/lon bbox currently supported)\n",
    "hucs = [\"15060202\"]\n",
    "\n",
    "# provide information about the datasets you want to access for run inputs using the data catalog\n",
    "start = \"2005-10-01\"\n",
    "end = \"2005-10-03\"\n",
    "grid = \"conus1\"\n",
    "run_ds = \"conus1_baseline_mod\"\n",
    "var_ds = \"conus1_domain\"\n",
    "forcing_ds = \"NLDAS2\"\n",
    "# cluster topology\n",
    "P = 1\n",
    "Q = 1\n",
    "\n",
    "# set the directory paths where you want to write your subset files\n",
    "home = os.path.expanduser(\"~\")\n",
    "base_dir = os.path.join(home, \"subsettools_tutorial\")\n",
    "input_dir = os.path.join(base_dir, \"inputs\", f\"{runname}_{grid}_{end[:4]}WY\")\n",
    "output_dir = os.path.join(base_dir, \"outputs\")\n",
    "static_write_dir = os.path.join(input_dir, \"static\")\n",
    "mkdir(static_write_dir)\n",
    "forcing_dir = os.path.join(input_dir, \"forcing\")\n",
    "mkdir(forcing_dir)\n",
    "pf_out_dir = os.path.join(output_dir, f\"{runname}_{grid}_{end[:4]}WY\")\n",
    "mkdir(pf_out_dir)\n",
    "\n",
    "# Set the PARFLOW_DIR path to your local installation of ParFlow.\n",
    "# This is only necessary if this environment variable is not already set.\n",
    "# os.environ[\"PARFLOW_DIR\"] = \"/path/to/your/parflow/installation\"\n",
    "\n",
    "# load your preferred template runscript\n",
    "reference_run = st.get_template_runscript(grid, \"transient\", \"solid\", pf_out_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "672812cb",
   "metadata": {},
   "source": [
    "### 2. Get the desired ParFlow i/j bbox from user provided geospatial information "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9d50bbca",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ij_bound returns [imin, jmin, imax, jmax]\n",
      "bounding box: (375, 239, 487, 329)\n",
      "nj: 90\n",
      "ni: 112\n"
     ]
    }
   ],
   "source": [
    "ij_bounds, mask = st.define_huc_domain(hucs=hucs, grid=grid)\n",
    "print(\"ij_bound returns [imin, jmin, imax, jmax]\")\n",
    "print(f\"bounding box: {ij_bounds}\")\n",
    "\n",
    "nj = ij_bounds[3] - ij_bounds[1]\n",
    "ni = ij_bounds[2] - ij_bounds[0]\n",
    "print(f\"nj: {nj}\")\n",
    "print(f\"ni: {ni}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83807022",
   "metadata": {},
   "source": [
    "### 3. Make the mask and solid file\n",
    "You only do this if you providin a huc or list of hucs. Otherwise, the reference run provided is for a box domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9c3aec38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote mask.pfb\n",
      "Wrote solidfile and mask_vtk with total z of 500 meters\n"
     ]
    }
   ],
   "source": [
    "mask_solid_paths = st.write_mask_solid(mask=mask, grid=grid, write_dir=static_write_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2da721b9",
   "metadata": {},
   "source": [
    "### 4. Subset the static ParFlow inputs\n",
    "Two options to subset static inputs. \n",
    "1. subset_static(): This function when provided with a variable dataset hosted on hydrodata will subset all static inputs required to do a baseline run. Note that the function will raise an error if any of the requested variables do not exist in the dataset, so we need to modify the default variable list and remove \"mannings\" and \"pf_flowbarrier\". Pressure is the steady state pressure.\n",
    "\n",
    "3. subset_press_init(): This function will write the subset pressure of the last hour in the last day before your start date in the given time zone. If no such pressure file exists in the hydrodata run dataset specifed, no file will be written. The function assumes UTC0 as the default and will return 11PM UTC0. You can override this by providing a timezone. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ec16ed87",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote slope_x.pfb in specified directory.\n",
      "Wrote slope_y.pfb in specified directory.\n",
      "Wrote pf_indicator.pfb in specified directory.\n",
      "Wrote pme.pfb in specified directory.\n",
      "Wrote ss_pressure_head.pfb in specified directory.\n"
     ]
    }
   ],
   "source": [
    "static_paths = st.subset_static(ij_bounds, dataset=var_ds, write_dir=static_write_dir,\n",
    "                                var_list=(\"slope_x\", \"slope_y\", \"pf_indicator\", \"pme\", \"ss_pressure_head\",)\n",
    "                               )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "33d3e618",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "UTC Date: 2005-10-01 00:00:00\n",
      "Wrote /home/ga6/subsettools_tutorial/inputs/conus1_upper_verde_conus1_2005WY/static/conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb in specified directory.\n"
     ]
    }
   ],
   "source": [
    "press_init_filepath = st.subset_press_init(\n",
    "    ij_bounds, dataset=run_ds, date=start, write_dir=static_write_dir, time_zone=\"UTC\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afefaac7",
   "metadata": {},
   "source": [
    "### 5. Configure CLM drivers\n",
    "This function will get the clm drivers that are associated with your run dataset (same dataset as where you got your initial pressure file). Vegm, vegp and drv_clmin will be written into your specified static input directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6f4133f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "copied vegp\n",
      "subset vegm\n",
      "copied drv_clmin\n",
      "edited drv_clmin\n"
     ]
    }
   ],
   "source": [
    "clm_paths = st.config_clm(ij_bounds, start=start, end=end, dataset=run_ds, write_dir=static_write_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "002c117a",
   "metadata": {},
   "source": [
    "### 6. Subset the climate forcing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f9e57c",
   "metadata": {},
   "source": [
    "This function will write all variables needed to run CLM for your specified forcing dataset, on your specified grid, subset to the i/j boundary that was returned previously within the specified start and end date. This function assumes UTC0 by default, but you can override it by providing a timezone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f34c4be7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading precipitation pfb sequence\n",
      "Reading downward_shortwave pfb sequence\n",
      "Reading downward_longwave pfb sequence\n",
      "Reading specific_humidity pfb sequence\n",
      "Reading air_temp pfb sequence\n",
      "Reading atmospheric_pressure pfb sequence\n",
      "Reading east_windspeed pfb sequence\n",
      "Reading north_windspeed pfb sequence\n",
      "Finished writing atmospheric_pressure to folder\n",
      "Finished writing air_temp to folder\n",
      "Finished writing downward_shortwave to folder\n",
      "Finished writing precipitation to folder\n",
      "Finished writing east_windspeed to folder\n",
      "Finished writing north_windspeed to folder\n",
      "Finished writing downward_longwave to folder\n",
      "Finished writing specific_humidity to folder\n"
     ]
    }
   ],
   "source": [
    "forcing_paths = st.subset_forcing(\n",
    "    ij_bounds,\n",
    "    grid=grid,\n",
    "    start=start,\n",
    "    end=end,\n",
    "    dataset=forcing_ds,\n",
    "    write_dir=forcing_dir,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c65a6b",
   "metadata": {},
   "source": [
    "### 7. Spot check subset static and climate forcing with plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2c235d3",
   "metadata": {},
   "source": [
    "#### Check a static input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "94ce138b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 90, 112)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1540b12aa160>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGTCAYAAABal3q3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACFgUlEQVR4nO29e3xU9bX+/0y4hEgy4SIk8AU0FQsoys0UAcuxSkFLWynYA4gWW1qrDVSgN/FUoV4atZ6jVQFrfxxoK1GrBS/0SItYoJwCYhQrB4mAtGAh0VbJAEoIZP/+QCbrs4ZZ2XtmMhmS5+0rr9ck+/bZn733uFnPetYKeZ7ngRBCCCEtkqymHgAhhBBCmg6+CBBCCCEtGL4IEEIIIS0YvggQQgghLRi+CBBCCCEtGL4IEEIIIS0YvggQQgghLRi+CBBCCCEtGL4IEEIIIS0YvggQQgghLRi+CBBCCCEZSGlpKYqLi5GXl4euXbti3LhxqKioiC7/29/+hlAodMqfp59+2vdxQuw1QAghhGQeV1xxBSZNmoTi4mIcO3YMt956K7Zu3Ypt27ahffv2OH78ON5//31nm8ceeww/+9nPsH//fuTm5vo6Dl8ECCGEkNOA999/H127dsXatWsxcuTIU64zaNAgDB48GIsWLfK939apGiAhhBDSXDly5AiOHj2a9H48z0MoFHL+lp2djezs7Aa3ra6uBgB06tTplMvLy8uxZcsWzJ8/P9CYGBEghBBCDI4cOYLOOWfgIyT/v8vc3FwcOnTI+dvcuXMxb948c7u6ujp8+ctfxoEDB7B+/fpTrvOd73wHa9aswbZt2wKNiREBQgghxODo0aP4CB6uQ3u0RajhDeLtBx5+c+gQ9u7di3A4HP27n2hASUkJtm7dGvcl4OOPP0ZZWRluu+22wOPiiwAhhBDig7YIJfUicJJwOOy8CDTE9OnTsWLFCqxbtw49evQ45TrPPPMMPvroI3zta18LPB6+CBBCCCE+yEIIWaHEXwSyAioLnudhxowZWL58OdasWYOioqK46y5atAhf/vKX0aVLl8Dj4osAIYQQ4oMsJFd8J+i2JSUlKCsrw3PPPYe8vDxUVlYCAPLz85GTkxNdb+fOnVi3bh3+53/+Jy3jIoQQQkgaWLhwIaqrq3HppZeiW7du0Z+nnnrKWe+///u/0aNHD4wePTqh49A1QAghhBhEIhHk5+fjxlAespOQBmo8D496B1FdXR0oR6CxoTRACCGE+CDd0kC6yNRxEUIIISQNMCJACCGE+CArlKRrAEAKahKlHL4IEEIIIT5ortIAXwQIIYQQH2SFTvwkvH3qhpJSMnVchBBCCEkDjAgQQgghPqA0QAghhLRgQqFQTAvhQNuncCypJFNfUAghhBCSBhgRIIQQQnxAaYAQQghpwdA1QAghhJBmByMChBBCiA9CSO5fz5maLMgXAUIIIcQHKSkxnIFk6rgIIYQQkgYYESCEEEJ80FxdA402rvnz5+Pss89Gu3btMHToULzyyiuNdShCCCGk0TnpGkjmJxNplIjAU089hdmzZ+PRRx/F0KFD8eCDD2LMmDGoqKhA165dzW3r6uqwb98+5OXlJVXBiRBCSPPH8zwcPHgQ3bt3R1ZW4/6bu7lGBEKe56W8O/LQoUNRXFyMRx55BMCJ/7n37NkTM2bMwC233GJu++6776Jnz56pHhIhhJBmzN69e9GjR49G2XckEkF+fj7mZXdAuyT+gXrE8zCv5gCqq6sRDodTOMLkSHlE4OjRoygvL8ecOXOif8vKysKoUaOwYcOGmPVrampQU1MT/f3ke8nPczsi55MJn/DjSb6O7b39tvP7oTffjX7+41vvxd2ud7t2zu8DH7st+vmZr/2HeUxrbG/c/2z0884jRxLa7sIvn+css85J71POx1+f3+Ysk+MZ/++DnWVt7vrv6OfyCy92lg34/rjo56wrpzjLXr/8q3HHHfr0pxEP67q9szviLPtUUf3DE55zs7OsbsXvTrkPAOj43B+jn4/9+h5nWeuvuS+ncvnB373sLJNzPrpf/OiWdW1+d9eTzjK5H32+1nVqdfNPnN+P/3xu9LOeb3keci4A4MOrRkc/5/9/i5xlWV3iv5Q/9emB0c8T394Sd72GkMfPm3CZs0xeG30vDvnrxrj7rP3xN6Kf9b0v72GNfPY0+p62nhN53ay5Ob7x987vrS4e6/wu96vP1zpHa27SgTwv+Vxq9D1s3W/xiBw8iJ6fPh95eXmBtw1KFkLISsIEmKkRgZS/CPzzn//E8ePHUVBQ4Py9oKAA27dvj1m/tLQUP/nJT2L+nhMKISd0YtrC7bJ9Hdtr28b5Pat1K7G/+JcgV4WTwu1zfG3X0NjkfvV+/G4XDnBOep9yPvQ5ym31MdqIN9WYuRHHyMrL9T3ukHG+1nWLOb5YJq8TANSJ/ch9AHDevo+psbRWb+ZyeUjtx5k3tSzeeoA7bzHLjPO1rlMrNf/HxXI93/I89L9Ejstlap9Zxr9anLEl8a8befw849rE3AvGMWuNe9/vs6cJ8pz4nZvj6h5updZ1nim1zDzHJv7XpjyvOjVvEn0PW/dbQ6RDSmZlwUZizpw5qK6ujv7s3bu3qYdECCGEtBhSHhE488wz0apVK1RVVTl/r6qqQmFhYcz62dnZyM6OfUOf+PaW6Ftt3Xt/j/79+H1uGPeppZujnyc/c7+zLBdP1e9vQK+42w0cf4GzzHvrjejnPjlKNlDrSg4+5YZcKz6uDw9OuvM639sVv7M1+vnDz49wluVNrA+j9tn1tLPsydt+4/yuxx5vWdZXJjrLar8/+ZRjAYDj//ts9HPd80ucZYNu/Sr80vqb9WHszZ/qH3c/g/sNcJY9cfX3o5/7XPdjZ5kca0ffI4md446r/jf6WQcbJw6oj2rJewgAxvavj4LpufdENEwvyxX35kD3dDGob9/oZ32fYMo1zq/heT+Ifq5b/pSzbOeu6ujnim69nWVyPMVdz0I8HlfbXbt/Z9x1NfK+aTVinLNMzre89wD32ZfPEwAM+v/qI4l6bpzzVdtVqOdEnkdvYz9blr3pLCsWXzf6e6H4/icQD3m/D/7NXc4yOU9A7PMnaSOOUXx/3NWahNC5g6Kf29w/zlkmr2mWut+s+yQTaK7JgikfV9u2bTFkyBCsXr06+re6ujqsXr0aw4YNS/XhCCGEkLRA+2AAZs+ejalTp+Kiiy7CZz7zGTz44IM4fPgwvv71rzfG4QghhBCSII3yIjBx4kS8//77uP3221FZWYmBAwdi5cqVMQmEFnXv70XdkROJJDIErcNzUg6IzPtZQuPVIV6g/vf1ETes+LBaU4bjrVD86z91w/gy/C3Dj4AbutWSQkx4WGDJDyEVYreQUoEOVUo8lfx56I090c8y3NsQvc/Jd36XY9XXVIbfLbTcIEOsWV++3lmmw/9yWy13yHPUY/n91qq4y2QYd8gP/+4sK1eZ4i5vGstcBhvL5PkPVtdUyggybAsA1UJ+0FKAlAoakgleEzJO8Tvj4q7X6oeuo6NswOVx15XPlL6HpBxwzRurnWV6vuX11vuR6GXHhDQRBCvcH7PPDAyP+8Hb8Xr08+PqGsp7RZ+vfjYzjRNNhxL/Z30IKXfrp4RGKzE8ffp0TJ8+vbF2TwghhKQVugYIIYQQ0uxg0yFCCCHEB83VNXBavAhIq9kQpSFJDVNr7RJtIbKwtH5rXX0Mv/vR2uPgefWWImmXA1ztuSG73vSbfxn9/PBNn3WWSa07V1krZR6GHps1x5a10tLs9T57i7yAmGViPDoPQerbQ1TFNm1Lk7RRVi95j0mt88Qx6+9Fvc8+xtzIdXWei3VvOtY+Q1sGXBukzFcAgEl967VYna8ir5vWz7W9TTJxSrE5HklDYz+JzgmQ56+vqbQSW1ZObXPV5yTzUPKXlrkDEvOh78VB4rPObZD3YrWyecr7Vufg6LybTELnj2jrn0Ra/yY/4y6zcivk/ae/T+RzKr9PDtXVxd1fqmmu0sBp8SJACCGENDXJlxjOTP9gpr6gEEIIISQNZGxEwNv1V3if1KuuFSFAHVaVoWIdKteV9uKhQ5zyGNo+OElUegOAio/rQ5JaCpAhX71Mhmet8KuuliiRFRBPxfw/PhL9bFWa6409zjIZgtVhTUlMSFtKCm+41fqsc4zZz664qzphVS03WPMtz0mHv/VY5TGOqbDy40Kq0XZNeS9UbFUykVHpTo5VSzFyLLqyn3Vv6HtaPgva2mhJQdL2N2TjIGeZYxdVx7dsiFq2sCoySvR1s663lLu05TZ3gBt+d0L+xjH0vEmrm35OZJVHLb1J9HOpZapMwpICLGIqBIrf9T0trZ4xttpl7vPeFFAaIIQQQlowoU9+ktk+E8nUFxRCCCGEpAG+CBBCCCE+SHevgdLSUhQXFyMvLw9du3bFuHHjUFFREbPehg0bcNlll6F9+/YIh8MYOXIkPv74Y9/HyVhpYMsNd0Z7bEvd0rSzKVuUJIgl0DlG1QFnmbQyAsAk8VnbsvyiLYISK3+hIUvk2P712qjWKeU5alvUQJGToZdJG1hv1bVPrqu3s8rf6jwMqX1rfdGyAcq5CqlcDm0hk8iOjhrrmlrLrHK8OrdA6veDJrp5LnI7fQ9b9801ygbXZ1m93qo1euvZsKyFDipfQ1u/LNuptUyin/2BA+qPIfMVANdOGJO7s9W1GsrroZ/vgeJ+08ewNHO3pLJrnaw19vn7JDo8pgKrM2CqkLk91vlZOV9yu0gkAnSLn4eRStLtGli7di1KSkpQXFyMY8eO4dZbb8Xo0aOxbds2tG/fHsCJl4ArrrgCc+bMwcMPP4zWrVvjjTfeQFaW/3/nZ+yLACGEENKSWblypfP7kiVL0LVrV5SXl2PkyJEAgFmzZuG73/0ubrml/h9wffr0CXQcSgOEEEKID1IlDUQiEeenpqbG1/Grq09Ezzp16gQAeO+997Bp0yZ07doVw4cPR0FBAf7t3/4N69evD3ReGRsRuPDL5yHctg0AN8yrw7FWiNuq+ietQDo8J7kk7F9SsI5hISsJArYtcNCt9XOxSFQOBFy7IOCGJ7Ervu0ySMc1GdbT4d9BE+NfJ22TkiH/aQUdnGUy5H3NG/Etazrk6toJ3fOVYUZ9TtpmOll0P4zpaidsgfp6W90fpWyhjyfv05ixWN01lc1Sjkfb2SyrnSQmTC/uG329tZVXoiv9yfMfpGQbiSW3xEga0gKrrpM8nlV1UKMtqRLrObH40JDQdPj72lWuNCHRslgqrIZ6bEG6hkp0hcQYy2AcrHOy7gUpmX3spa+y4Inug8ltDwA9e/Z0/j537lzMmzfP3Laurg4zZ87EiBEj0L//ifv0nXfeAQDMmzcP999/PwYOHIhf//rXuPzyy7F161ace+65vsaVsS8ChBBCSHNk7969CIfD0d+zs7Mb3KakpARbt251/rVf90l55W9/+9v4+te/DgAYNGgQVq9ejf/+7/9GaWmpr/HwRYAQQgjxQarqCITDYedFoCGmT5+OFStWYN26dejRo0f07926dQMAnHfeec76/fr1w5498SPdGuYIEEIIIT7ICoWS/gmC53mYPn06li9fjpdffhlFRUXO8rPPPhvdu3ePsRS+/fbbOOss/66PjI0IZH1xArI+KTEs9aaBqjuXZafzq4Xq8rdSo588bwAsZJnRsaqUqURb1JySr0r7nfOXeguPzlGQmrHOCdCY9sJdxnpCi9X2xeJV9Rqe1hfzDO3X0jMtW5i2rElNVXcxs/JF5HUK93OvqTYkPmGUEZYarj5/bT2TmHY2n50x9TlZOSnaLjpR/G6V6bbK/1q5NDp/YIjqxic7BWqsUuByrqx50yWd5T2st9PPoi5BLPGbFxDTmVDcp1rr9tuJEXDv8URzAvRzIu/TRHMCNH5zAgD3XrHOKWaeRP6GnN/IkRrc8KMFvo+fDOmuLFhSUoKysjI899xzyMvLQ2VlJQAgPz8fOTk5CIVC+MEPfoC5c+diwIABGDhwIH71q19h+/bteOaZZxrYez0Z+yJACCGEtGQWLlwIALj00kudvy9evBjXX389AGDmzJk4cuQIZs2ahQ8++AADBgzAqlWrcM455/g+Dl8ECCGEEB+kOyLgeZ6v9W655RanjkBQQp7fI6WJSCSC/Px8vNyle7SyoAzdO5Y4uOFB2bkKsDueyXWtClraFqMJnVtvb9P7kVW6vB2vO8vkeVjVEi25Q4djs74y0RyrpGT09LjLpByh51Sehw4HyhCkDMUDsXNjrSuvm1V1TluvZIhXh39DQg7Q9kxPyU0yXKnDqnKs+pw2xwldAm441rKTWTKBDnFbYWvZ/Q5w79OyAZfH3a9lEdRjkzKVfi71+cs51lUf/Vbk1Ocr5Q+9Tyvcr9lp2BDjHQ8AWgn5Q8+pVR1Tkir7XlNj2QD1MyTvBUum0s+w3E7uPxKJIL9bL1RXVwdKwAvCyf8v/Sr/TJwRSjy17iOvDlOr/9moY00EJgsSQgghLRhKA4QQQogPQqEQQgEz/53tM7QRMV8ECCGEEB+kO0cgXWRsjkD1/j2+NBSp4WstTmr0uvuc1Gwte83jqhuY1sylTqzXlda7VspOZXZy80lDtjOrxLKFVYJV6qSW9cfSBfW2uqyr1Je1tUza+axy01rPtcaaaHlUjVV+WdoHpX4MuOWXrbK9DZWw/r1h+/S7X2s7rZ9LPVvf+zqfQecMSOTzprVma9xW7pBcpkszWxZJjTWnfkufa90/SC6NhdyPZV0NgsxZSCZfwcqPskqoy9wO/QzL+ZbfJ5GjtThz0UtpyRH4TYcuSecIXHfg/YzLEWBEgBBCCPFBFpJLrMvUpLzA41q3bh2+9KUvoXv37giFQnj22Wed5Z7n4fbbb0e3bt2Qk5ODUaNGYceOHakaLyGEENIkhELJ/2QigSMChw8fxoABA/CNb3wD48ePj1l+33334aGHHsKvfvUrFBUV4bbbbsOYMWOwbds2tGvnv5OfX6QtSiPDbJ6qJuc3rKbDuJYNTq8rQ8w6VG51g5PL1kfccLC09lmVEwE3PL1o4Z/j7scipprd1voQ97X367Xr0aFSHSqX6ApiMjys51SGgHWoXFaXuzZAWDNRKUAjw8G6Ct+1+3fG3c65h5e5kpHc55Clv3eWaXlJ3g/aTtf7nPhVAeW86XHLkH7IeIYauhclMfeGCPPq48sQvw4py3thyEZ3bizrsK4WaUmB1nnJ8LS+h+V+JjZC10AgdXKARMoByUhmWg6QyHFrWXBR1YHo50tiJJz6e0N+Lx2qS1/3weZK4BeBK6+8EldeeeUpl3mehwcffBA//vGPcdVVVwEAfv3rX6OgoADPPvssJk2aFLNNTU2N04s5EokEHRIhhBDS6IQ++S+Z7TORlEoWu3fvRmVlJUaNGhX9W35+PoYOHYoNGzaccpvS0lLk5+dHf3SfZkIIISQTCKXgJxNJ6YvAyYYIBQVuFm5BQUF0mWbOnDmorq6O/uzduzeVQyKEEEJSQnN9EWhy10B2djays7PNdW5sX99/WXfck/qqZQOM0bd86l0xVhe1nVW6tk6WH1Y2IRgd13RegERanRrqRtdb6K06J0DuR2uo0rYzUDVfLBb6ptYQ5fkH6bCmLWPSaqeLJstz1Lavpi7PapUKvkbYqXSeidSzrRK3ejttyZPau2Vn0/jVwXXZXut8LUueHovMZ8hV68o8gJgOhyIvQM+NHJseS8x5iLnS+QNWSWlpkdPInBD9nFiWZ43fXKZUWWCD7MMqv25tq593v8QrBR05dhx4f19C+yQnSGlEoLCwEABQVeX+j6Wqqiq6jBBCCDkdyQKQFUrip6lPIA4pHVdRUREKCwuxenV90Z1IJIJNmzZh2LBhqTwUIYQQklZCKfgvEwksDRw6dAg7d9aHvXbv3o0tW7agU6dO6NWrF2bOnIm77roL5557btQ+2L17d4wbNy6V4yaEEEJICgj8IvDqq6/ic5/7XPT32bNnAwCmTp2KJUuW4Ic//CEOHz6MG264AQcOHMAll1yClStXBq4hcOzX9+BYuxO5A48efjf69xhd2vCrJorU/rR3WuNoYUaJY4sg5X+lnq+13ZjcArGuXrb+L/Vjm1bQwVkmNXrtf7f0PV1G2cLxWavSrRKrxGxMqVhRHlW34U1VrQCJzkmR6HPy21pZ3wtyP/p+isk7EVp3kBbNFSJfRdYUAABZqUCPW+q0ejtdHlbmGljLrDoGuky1rKOg8w5kPoHep1Vi2KpxkKX8//KaWvkpuvyxHJssLw3E3u+yjoKFLqPccdW4uOvK79Anrv6+s8wqix6k/LGVByGvh77fH/n5t6KfdXlxmXckv4fqjtYCr7j3RmOSmf+mT47ALwKXXnoprPYEoVAId9xxB+64446kBkYIIYRkEslWB8zUyoKZmrtACCGEkDTQ5PbBoOiSwrqrYCpwQmBGSWHAthdJtIRhdQaUVr+YkL5hLfRbNliva4WjLSkg6yuuuS9I6FDatGLkFyFNWGhbmA55x0OfU6IlX3WoWqJDvDJ0rcPYjoyhwshSppnUd4k5Hhly1iFup6T1LneZvBe1pCDPUVtJpe1vkvZ5KuR+nlLXV15HHWKX3fC0NCFDxbo0sNXbU0tqcr/aoqhLNUukxCHHCQD5S8vibidD41qm0GOTIf/cAa6MYI3NQobqr90/zlnmdB8M8DwnahcMUppaIp+vdJYYbq5tiE+7FwFCCCGkKchCCFlJ/O88mW0bE0oDhBBCSAuGEQFCCCHEB5QG0kzrr92C1uFwzN91ToBVdlPq+dpqJTV7S9/SdkWtZ0td3Mof0McYsrE+12Gwyh+QeRDri1zt0coDsNoZW2itXepvQTTTIEgNe+D4+Fq71tOlTS3GziYsch1XufeFvDbJtIGVFj5tS3PGosYtx6p1ePfecHMEJNpOZeVI6LbDg+fFb+crr4VVRthqJW3lJCSDPKf4Dcdjba6mPXVXfEuqzl+YOCX+MaW9Lve+W5xlx8Xv1r1wzRurnWU6l0jmCOi8k0G3Jv8s6jbAkvim3hNY35vWc2KVpi4W3+frb/6ls+zaBsaTDpqrayBjXwQIIYSQTKK5RgSYI0AIIYS0YEKeVR2oCYhEIsjPz8fLXbojN+vEe0qQTnbxiKnKFsAaI7G6fOljSBlDSxh+ZQtd+UuiQ55WBbUgsoFcV4eDZVjdsgUlE37XVjCJPGd9vtJ6pmUaKekkU2XQ6rgmq8TFhpjrx61DvLKLnqxWB9hVCC2rn55/q5qcvN5yLBpZyQ9w740ZC//sLLMkLG2BlR1FE7022r5nzZvGr4yh5Qd9TIm2BVr7sZDXTVchbGN0ApXVDPX3p3xuk3lOLeR4Yjq4CrQFVz7T2koar5Jn5EgNOv9oAaqrqxE+hZycCk7+f2lF525on5X4v58P19Xhi//a36hjTQRKA4QQQogPTnYRTGb7TITSACGEENKCYUSAEEII8UFzTRZsMS8CqcoJ0HpXnbD6BbE2ymWW7qy1fanLah168jP3u9sKnTDG+mV09ZN2J92BTNqNtEZtna+FzjXwW3ZU68BDRI6GtsFpO6FE5yRYer5fPVmfg5zHiqWXu/scUP+7vobFQjPXmrRV4thClncG3OuWqNauO1guqjrg/C5zBvS6VhldS8O2yg/Le1gvs66htmRayBLDVk5IEB1e2/lk90HLuqtzK4rfGXfK9YKOJxXo7wV53cxcCvW71UEzXTTXFwFKA4QQQkgLpsVEBAghhJBkCH3yXzLbZyIZ+yKw88gR5IROBCx6i1CS06kNyVnB4iFD9TH7V7/LdbUNUNoJtTShrYYSKT/okL4MVcdUTDPQYcWKj+tDmXo/8hgybK3X1VUWpb1HhwN1+F92LtShysFfeTb62bq+xW4U3ZRYrLHpcLy0XumqcAMH1F8PPadWhz05j3q+5bqvXacrC9b/rqUYLX/IZ8O6v/R9KrG64enjS+T91BB6TiVaipJYVfCs6n1aCtDSiMSqrDhI2decsSq5x3mGlrrSk7QP6nOKkS0MG2CiSClMVzZsDFu1ZbMMYqW07MHporlWFqQ0QAghhLRgMjYiQAghhGQSWUjuX8+Z+i/vTB0XIYQQklGEUvAThNLSUhQXFyMvLw9du3bFuHHjUFFR4axz6aWXIhQKOT833nhjoONkbETg6l/fjXD7HACuvUiX2ZR6eqxlzC1R6ZcgeQeJ5ig4WpzS5Y6Lz1qHrxB68sAGjiHzC/TcWBY9qVNadip97pulnql0eJkToLfVerZcpq19lqYotwvPc5dJ7V13sbM65w1Z6r/krtSprTwAq2uh1oGlvqqvoUY+CzrXQNrQjqtOeXWGLc3KC9DzJtElhuX9Zln9gmjU8jot+ssBZ1npcLGfXfZ+posud9a49fnKe9wq4a2XyWsaky9hjNXS2q3vOv0MyfFoy3OWYQHWx/d7TMuuqrFKr+uxNgmf/I82me2DsHbtWpSUlKC4uBjHjh3DrbfeitGjR2Pbtm1o3759dL1vfetbuOOOO6K/n3HGGYGOk7EvAoQQQkhLZuXKlc7vS5YsQdeuXVFeXo6RI0dG/37GGWegsLAw4eNQGiCEEEJ8kCppIBKJOD81NTW+jl9dfSKa1qlTJ+fvS5cuxZlnnon+/ftjzpw5+OijjwKd12kXEYjpsCdCrtoK0xgE6WLo1z5o7UNXMpRhPR1u18hwtA5PylCxDvGPfSN+xbbiVfG7D8rOddreY4UDdcc9q2qYtFtpa6W8N2TVt4bQHfeKxfUIcr3zl5ZFP8sOkgBwjbCa6RDnxAH1tkMrjGt1kQP8VxqMrZZYf59oK6W8/5687Tdx96nvL8siqO8p+bs2xN7Yvkfc/egKhRJ5D2kpyNqPfC4A9/x1qNqvhU3PhbRoaqlLn7/VfdCSmCQ6NC/PQ9/P8vnS95uWAqRUoJddu3+cr7FpzO+JZlRZsGfPns7f586di3nz5pnb1tXVYebMmRgxYgT696+/Ttdccw3OOussdO/eHX/961/xox/9CBUVFVi2bJnvcZ12LwKEEELI6czevXudNsTZ2dkNblNSUoKtW7di/fr1zt9vuOGG6OcLLrgA3bp1w+WXX45du3bhnHPO8TWeQNKAnwzGI0eOoKSkBJ07d0Zubi4mTJiAqqr49aQJIYSQ0wGdnZ/IDwCEw2Hnp6EXgenTp2PFihX405/+hB494kfJAGDo0KEAgJ07/RdrCvQicDKDcePGjVi1ahVqa2sxevRoHD58OLrOrFmz8MILL+Dpp5/G2rVrsW/fPowfPz7IYQghhJCMIyuU/E8QPM/D9OnTsXz5crz88ssoKipqcJstW7YAALp16+b7OCHP87xgQ6vn/fffR9euXbF27VqMHDkS1dXV6NKlC8rKynD11VcDALZv345+/fphw4YNuPjiixvcZyQSQX5+Pj5YvjBqH5Q8cfX3425rWeK0TpeotdAs3RrA+mSWMTaOJ49hdSoDUjNXuqSz1Ey1Ju1XswTc+dflSUOio6O26Ek7m9aspWVPl/iVNsiYroUqR0DOsdaBpYXPb2c+wNVbrbKyVrnnhq6vPEddKlhq/5bWrLeT+rKeb3kPNdSV0W9HScs+aXURXB9xj//Iz78V/WyNG3BzTfQ9LHNNdE6OdZ/Ka6XzfKx7KFVlhCWW7VBjlWbWWFZeK5dH3tMxXVEFuvR4vHvhUF0dLnt/H6qrq51weyo5+f+ldYU9kJuVeI79obo6jKx81/dYv/Od76CsrAzPPfcc+vTpE/17fn4+cnJysGvXLpSVleELX/gCOnfujL/+9a+YNWsWevTogbVr1/oeV1I5AjqDsby8HLW1tRg1alR0nb59+6JXr15xXwRqamqcjMlIJJLMkAghhJBGIZQVQijoP+vl9gFTDRcuXAjgRNEgyeLFi3H99dejbdu2eOmll/Dggw/i8OHD6NmzJyZMmIAf/1j3LLFJ+EXgVBmMlZWVaNu2LTp06OCsW1BQgMrKylPup7S0FD/5SdM0kCCEEEL8ku6mQw0F7Hv27BnoX/7xSPhFIF4GY1DmzJmD2bNnR3+PRCLo2bMnttxwZzQEo0Oi8TBDWapi1yCjG51Eh+a1LUyGB60OXBq5zAqNa7lBhsesqm9ArBUsHroKnZQDdFjT6rAnw8jWnGqsELfGsshZ0oQMeev9W/KDRkoaFcrqJ+2r2hIpr5sO/8r96HAoltWHWPV1smQEaw51iDsXQirYGj+x1294/1TIZ9O6L/WzIG1iFSrEL+UAKQVotGyg5012kQyruXGkMPUsyPtWz7eUX3SXSkduU9fUso9a3Se1pGOF7aXcoZ9huZ3/3qaxNk9pyYy1q4rroZZJWTJX2WMli6oORD8fRcLqNvmEhF4ETmYwrlu3zslgLCwsxNGjR3HgwAEnKlBVVRW36lF2drYv6wQhhBDSlLANMRrOYBwyZAjatGmD1avr/2VUUVGBPXv2YNiwYakZMSGEENIEpMo+mGkEigiUlJREMxjz8vKiuv/JDMb8/HxMmzYNs2fPRqdOnRAOhzFjxgwMGzbMl2OAEEIIyVSaa0QgkH0w3tvMyQxG4ERBoe9973t44oknUFNTgzFjxmDBggW+GyKctGk8ltcJOaETAQtLU0zUIpeoTceyD+r8AUmiXQot+6BlbQNse5vMu9B2MqnpadulZaeSBDlfbaeTliJ9DNmJ0upip0sMtxIlfhvKCZAlWK3yx7rkbNmAy6OfLf3euvd02WCZ96CtnHqO5f2gLXN+80W0Di6xug3q/BytmctrpXN+/NpF9fH9dkbUY0sm10Eiz0Pel4CbS2LZRYOUZtb3tLw35vzF/Z549PC7cffjF52voPHbwdRazyoTrpF5EHKf6bQPbujVK2n74LA9exp1rIkQKCLg552hXbt2mD9/PubPn5/woAghhJBMI9nwfrOQBgghhJCWSnOVBtiGmBBCCGnBZGxEYHS/rgi3bgXAzQPQWqf0CF8SdrUoqVNpX+0xn3UEGqoNYJUKlvq+3o/UYoN47uU+9XY7DQ1Va8aDxWft87XKfjrrGtvp3A3L12y1fi0ZPd1ZVjq8PkdC66lSQ7dyFPR2eh6lpqvXDfUbEP2s9WyZF2DVA9D3gpwr3TI2T3jQ9TlZuRV6vuW6WjOXz5SeC3es7v0l92PlBACJ1wKRLYR1ToC8Fvp+c+bRyG3Q+9V5NtqfLxkichtifPzG8eRcWbo/4M7jwAHOopi8AIl1vSU6l2XGwj/HXVcjz8Mai5UHYdUY0GPrI+aiqeoIZIVCyErin/XJbNuYZOyLACGEEJJJUBoghBBCSLMjYyMC7+yO+LJpaDlAIkNSugOXDMlpy5YM6zZkg5PLtZ3PCbPqzoTGfv12JtToLnp1zy+JftZhNokVctXWQhnWs+yJlhSgiTlH8fujKlTtV9KxroXuPGmF2K2xDrrVtdrJY+jwf7x9AMDkZ+Ku6pQVLlbna5Wc7Rh/lzFhfL+lsIH4jUyCSAHSLgi4Fk0dRrY69Q26tf6zllTc9dx50jJZb3EMHarvY9zj8vnS97u03llWUh3u18hQ+SIjbC9L+urt9JzK+bBKdjd0DCmH6O9hszOlKPfeUIdDiVxXjuVQXR0Wv3/Y936SIYQkXQMBmw6li4x9ESCEEEIyiVDWiZ+Et8/QtgiUBgghhJAWDCMChBBCiB+S7ReQodmCGfsi8KmicNQ+OHBAfM3WKp0qrWeX6DKjRslTqe81pHVbWrC0+ul2wpbWrTVUv+hjQJTAzfqmu0yOW5du7X1OfF1W6ueDJrrtVaXeqHV3PceWLm3Z97S+L5F5IFrrjbd/ILbkr84vkMh50/qq3K/eZ/E746KfddlomQeg9Wy5H53nYlnbxqp15fzrMsIyRybG9hhnLAAwUGxnlZsGbAuwLNvcW9n35HXsDddaJ7V+nfdgXX+tmct7vGLrZr16FCsPQSPzAnSrX3l8PW87l7l5GJY9WiLtdA0x0MgLePimz0Y/63FrPV+OTecPyO+JCrUfa5lTOtgoNyz3ETlaCyzaF3fdVNJcXQMZ+yJACCGEZBInXgSSKTGcwsGkEOYIEEIIIS2YQN0H08HJLk9/+0zfqDRgIStazf/jI3HXk+FXwA2V6up1C3bXh59jwu2nEZY0ocPjEqu6V6hvX716FCkxBOnuGKTqolxXX1MZLrRC3KbNE3aoXG6rq9DlLy2LftbzLfep59DqdiiRIXTAtXkC9rxZneRkyDtIlUtLFtPylq7CGA9tgZXbaVunDF3rcL/cz/H7bnGW6ZCzlFi0bCHts/q6SWnIut+sZ01LX1aFPgstG8iwvazGCdiSknze9TzpsVlShZRU9b0nr5VVVXT6zb90lknZQhI5WoszF72Ulu6Db/T9FPJaNfz/pXgcPH4cA7a/c3p3HySEEEJaKs21xDClAUIIIaQFw4gAIYQQ4gO6BjIEbWeT+pfW+mXOgNaipJ1Ja11Sp9W6bBAN1UJqzVIjBgBvx+txtwtSctjKb5DzYXVu0+VYJ06Jfzy5T122Weurch6tEsMaqdlqPdk6X8uuaZWY1vqu1HS1lVFq5puVJi5zJnSOgnW/yfm/Zun1cbfT2+pjSM1Wn5PUuq250Fj3os4fkOdv5YToe1/mqwyMe7RYrNwCrZFLzVo/C3JZ73NczVxaFLUlVN6L2hIq7z99nabtcq280hZoafIyJwCIzQuIh1V6PAg6t+LG9j2inx/5+becZTK3JjLvZ84yOaf6HOLZCQ8dOx5ssEkQSrKOQFI1CBoRSgOEEEJIC+a0iwgQQgghTUFzlQYy1j74z2mjEG7bJma5VW1KywYyrGhVINThQMv6EqSrnkSHQ2VY1QorpgPLBqYtejI0rsPY0j7YUEVAGTp89PC7cY9vWf30vEm0Xc2qrKelAnlMbdGToUwdYpaSipZ75DXV8y3lDsuuqCtAarlLjk0/CzI0b9m5NH5toA1Va/S7bkOSkkTOh1VZz+qMB9h2WYm+NnKssqOgxuoSqSsg6u83eR31Ocrz0udkSX8S63x1ZcFr3ljt/F424PK428rvWy0hWTKdvDf0vEkceeHYcZz9yva02AffuqB30vbBfm/uzDj7IKUBQgghpAVDaYAQQgjxQSgrhFBWEsmCXmZqA3wRIIQQQnzAHIE0cVKLqd6/p1E1FNnJTeupUt+2dGAgdXbCVGDlGliWLUt7tnIErHnT2qfWUP3a8KzyzxqpWerrJG14Os/D0qz1nMpjBOl4J7Vu3bVQYunQulOdzC0AXE1Xj81CXkdttbPGZqE7JcprbHVR1NdNzpWlw1ud8TRWzoDejxyPvm7ymlrfA/p89XMj0XkA0kJn5UDp85XnqPMA5H2k7yF5jXUuiZVPoHMb5PW28qM0Vt6FPF+nS2FdHS57f19acgR2DPp00jkC577+9umdI7Bw4UJceOGFCIfDCIfDGDZsGF588cXo8iNHjqCkpASdO3dGbm4uJkyYgKqq+HWtCSGEENK0BHoR6NGjB+655x6Ul5fj1VdfxWWXXYarrroK//d//wcAmDVrFl544QU8/fTTWLt2Lfbt24fx48c3ysAJIYSQdHJSGkjmJxMJlCPwpS99yfn97rvvxsKFC7Fx40b06NEDixYtQllZGS677DIAwOLFi9GvXz9s3LgRF198caCBfXjVaBw/RffBIGF8iVV5S1vdrNCdZTVLt0wQxHZoWeYqVPXAwcYxrbmRaNuXrsoWnveD6GdtEZRShSUF6K5yMhyubX/XvHF99LOWAqxwuJ5TGZLUkobV4Q+oD9tbHe40MlSrQ8xaipD71ddbyz9+jpcMltyjkXOqr5u0ielKklKmsSyQ1vEA9/pXKMucNW/F78SvQCmtb0G+s7RsIZdruadia/13mlV1cMbCPzu/Pxx3TTc0r78zxxpdG7VsIJ8Fy3ItbcSAW4Vw8G9cKazCsEGmC1YWVBw/fhxPPvkkDh8+jGHDhqG8vBy1tbUYNWpUdJ2+ffuiV69e2LBhQ9z91NTUIBKJOD+EEEIISQ+BXwTefPNN5ObmIjs7GzfeeCOWL1+O8847D5WVlWjbti06dOjgrF9QUIDKysq4+ystLUV+fn70p2fPnoFPghBCCGlsQkhSGmjqE4hD4BeBPn36YMuWLdi0aRNuuukmTJ06Fdu2bUt4AHPmzEF1dXX0Z+/evQnvixBCCGksTkoDyfxkIoHrCLRt2xa9e5/Qf4YMGYLNmzfj5z//OSZOnIijR4/iwIEDTlSgqqoKhYWFcfeXnZ2N7OzsmL/nXtADuZ+UGJbWFK3nWpYiiaWhWd3ItC6ntUCpW2qNWFqPtPbpt8SwXia7s4XOHQQLua0ueSt1ectqpsvYajuhxNFCG9CkLb1PlzKV1D2/JPpZ5yHsFPeCznOQc6rL38Z0GHyjXovXWrfOJ5FYXRvlPWVZMn9v3EMaXQJWari6/LNl+5RY1r5Ey2sDbh6Eft6kvqznRj7vutuizInQ21n2QU1vMW9Bcg0k+p7S8yiR12bigPiWRMD93rCePQv93SfvGytfQncNnH7zL53fpxV0iHtMuV+rTLg+hvye0vk51+4fF/0sv2s/9urijuN0p7S0FMuWLcP27duRk5OD4cOH495770WfPn1i1vU8D1/4whewcuVKLF++HOPGjfN9nKRLDNfV1aGmpgZDhgxBmzZtsHp1/Zd4RUUF9uzZg2HDhiV7GEIIIaRpSdYxEDAgsHbtWpSUlGDjxo1YtWoVamtrMXr0aBw+fDhm3QcffDDhiEOgiMCcOXNw5ZVXolevXjh48CDKysqwZs0a/OEPf0B+fj6mTZuG2bNno1OnTgiHw5gxYwaGDRsW2DFACCGEZBrpdg2sXLnS+X3JkiXo2rUrysvLMXLkyOjft2zZgv/8z//Eq6++im7dugUeV6AXgffeew9f+9rXsH//fuTn5+PCCy/EH/7wB3z+858HADzwwAPIysrChAkTUFNTgzFjxmDBggWBB2VhdR/UyPCgJQ3E2GREqNySDQA3XKmXSesT8EdnWcdV4+KOR6LDY3XGMm3RktKBlBQAN+S4c5cbxpdhe6tro67CJs+/IUumDDPr0LxEj1uG3PV8y2Nq2cSyeeqxOuFwNR4pFeiQp5zTsf3d6y2xujZa1du0XfDaVfGtlZrfG/NmdQrU5+8XPTfyudXSj7zG2r43UHyO6ago5LWKj12pyXreLax7OlG0vGZ1A1wfca1+0uacaPXEIB0GB4rPupLiwzd91vldWk31M+yMVcmExeL508+ilFQse66uLIiDOK3Q7rh4ErmmuvrEM9CpU6fo3z766CNcc801mD9/vinDWwR6EVi0aJG5vF27dpg/fz7mz5+f0GAIIYSQTCWUdeInme0BxLjj5s6di3nz5pnb1tXVYebMmRgxYgT6969/8Zo1axaGDx+Oq666KuFxsekQIYQQ4oNUSQN79+51eg34iQaUlJRg69atWL9+ffRvzz//PF5++WW8/vrrxpYNk3SyICGEENIiyAol/wNE+/Wc/GnoRWD69OlYsWIF/vSnP6FHj/pqjC+//DJ27dqFDh06oHXr1mjd+sS/7SdMmIBLL73U92llbETg0JvvIusUJYaDIO0tuqvXWPFZ6nCAXapYI3VLreH1RnyrmURa4gAgyyhVLPMCtF3R1H5VPoHUcHXpVlw8FvGQemfveT9zllm2NK2TylKmshwroHIf1Lgn3VlfOtfNwXDR+ROe0O/lsYFYLRJmqeB6rJLSO3e5Vi95b4xV+SLy/LUuK39vqPyv1Ffzl5Y5y2Tei17WGOjSvDIPQ9/vcl1tZZX2YP18VYjrZJUQt7r2nWp5PGKeEwOZB2LlBOhnVuv5fnOidE6EMx8qDyDuempdPTb9nFpYeVa9jTwAnT8jkbkl8vmORCJAt/idP09nPM/DjBkzsHz5cqxZswZFRUXO8ltuuQXf/OY3nb9dcMEFeOCBB2JaAlhk7IsAIYQQklEk2zko4LYlJSUoKyvDc889h7y8vGiV3vz8fOTk5KCwsPCUCYK9evWKeWmwoDRACCGE+CDdlQUXLlyI6upqXHrppejWrVv056mn4jcqSwRGBAghhJAMxPO8tGyTsS8C8UoMWy09B010NSxdElNitfe06g9ofTHedoCbh2BpX0HaF8uywdqPLduyAq4WJ9v+Aq6+nq/KD8syo7qsqdRJtdY8WPjBdWtdq9Wu1vo7GvPhzJXaTt4bg4yyprK886mQNQ6CtHq2kKWCrVbWlibcMaEjn0BeK+scdC2KhuYqHlpPl3kBVvllXSsBu+pzBPT9/p2i+vt7/h8fcZbJ+1aXJY+pVSCWD1Ilfp2xqtwZq/6CrCkxcUr8MsL63tclhn8vxmrVDdDL5PeWXiZzp4KUYtb5KzInRc+Fcx+pHAnnu1eVYo7J1xEk+uylFJHwl/D2GUjGvggQQgghGUWacwTSBXMECCGEkBZMxkYEfrBoA9p+0qFBhrJiLERb639f/xdXCki0zGiQ7WRYU4fZZKh6kZIpZCAziDRgoSWOOX+pD2tfosJzVsdBGZ4c/JuJzjJZUlnLDRa6S6QMJVohb231k2PTYcQhIox//L5bfI9NH0Pa9KxwpFVWVYexrf3IMHJ4ntulUkoq+ng6HGuFVaXEYckdiUoBmpjzFfLTZNWJU9rrehvSiO4+KL8XtIQl50ZfX23Rc+9N9z61OkrGL9zrosv4Trqz/h62JEON7vYn7zcrxK+3k1Zqa58NIceuS0r7RVsLJwtJQV9TaUFO1X0alFBWCKEkwvvJbNuYZOyLACGEEJJRNFNpgC8ChBBCiA9CoSQjAhn6IsAcAUIIIaQFk7ERgeu65CM368R7ShDLnsTSwa2yooluZ6HLGEvbkGWX03quLs8qCVIaWa7bJ4AtSpbOHaJa/cp2skHaRWsNV2q2uqyplU8hdeksVY5XHqOhUr03tq+v5a01VDkePU9SM9bnJMej9VS/OSID1e9BrI2psF4FsRY2VP5aIjX7Vj+8x1k2UOR6aGvdQGFn0zo8RD6F/v7Q+RvF98cv2y1b6Mp8AY2+FtbxysXzpq2N+jzk3Oj7VtqDJyrbpZwrfX/JXBptOba+T7Xtsk/Om3HWdJ8TnZMRc60EEVG2XLcHl8jr9LFXF3e9lENpgBBCCGnBZCHJOgIpG0lKydBhEUIIISQdZGxEYOeRI8gJnXhPscJV0uqnQ/qySpkOVcuQnA5VyePpkL4VrrpEWxvl77vcdeXxrZCrlgLMMLKyN8m50aFZ3Y1RIkOO2t4kw+E6xB1StjBnn8a8aWuhDCvq85VzpasVyrByTEfHr9TbID9soJrZo/fHHaoZAraQ49bV5H4vrpusaggA3lv13RZnLPyzs+zR+9Nbac2SAvS86HCwvDaWTCFlGUBXyHPP/+GbPhv9PPkZ96LJELO0GAOx4XC5rZ5/adkMYvOV103b4OT3grZL6uPLEPikvu7zJu9b/R0ij69tp7Lqo/7OlBKDvt8StWNrSWfyM/XPYiRIB1PxWd5fkaO1uGHRSwmNLSiJ9AvQ22ciGfsiQAghhGQUzbTEMKUBQgghpAXDiAAhhBDih2bqGgh5ifQsbEQikQjy8/PxWF6naI6AxMoXkJoh4GrPWqN2yv8qvTyIFmZZDWWJX6vMp85DkNqfZTXTuqCmZPT06GfrnHT+gMxfkNo64L+0px7368KGBdi5BhKtp8ucBd0NTWqRepnMA9Hnq+1l1jW1yvgmyodGTkaqyk83Nla5ZY3OH7B4yui+J+9pnSOQKkLCIhuka6PMC9DfL4/8/Fv1+1d5Nfr5kjZUneckv8O0DVFaDfXY5Dnpss0WWr+XuTz6+Zbft0HKKEv0PuMdO3KkBp1/tADV1dUIh8MJHashTv5/af+XPoNwm8T//RypPYZuL7zSqGNNBEoDhBBCSAuG0gAhhBDih2YqDWTsi0Dvdu2ilQVlaM0KcVt2l4ED3NDZQBGRezjAuPQxxsZZTx9fh0qlVKBDboMbCPnH4wnVYVCGIK0wmw4rSquXDh3KyoJWOFKHtN0ahDZZolOdrrooQ9D6+BYydCytVQDQW8kPMvyvKwQmirTX6RBzY8gN6UBa27ScMmiiWxFSSjW6eqDsFKnD31LGsbprrhcyGAA8evhdc+zx0GH07xTVyzZ6n/JelJY8wH3e5XMIaLnLvfdyB7iWWPlM5SkJTd63ek7l2GTXPsB9pnNVRUI5/3rZ4HmuTCjlD10BVH7fVChbs5SG9PW2ZAQ5p8Xie6F1JAL8aEHc7VIJuw8SQgghLZlmGhFIKkfgnnvuQSgUwsyZM6N/O3LkCEpKStC5c2fk5uZiwoQJqKryXwOfEEIIIekj4ReBzZs34xe/+AUuvPBC5++zZs3CCy+8gKeffhpr167Fvn37MH78+KQHSgghhDQpJwsKJfOTgSRkHzx06BAGDx6MBQsW4K677sLAgQPx4IMPorq6Gl26dEFZWRmuvvpqAMD27dvRr18/bNiwARdffHGD+z5p0/jXvd9BuF12oHFpq1l43g/iriv1Lav8rdbJEsXquKXtbFJvs+xzTyrtTSO7npVfHD+bQWuIUkPXNjy/pWKDEKT7oJwP69oE6X5o6aRav5dav9S2gYa7GsZD6tJ+7ZlNge7MJ7VenTujLbGy26XW+q17Wq4rcwIAO19IWg31erqMr0TnCMgy1vocZZ6PzhGwnjepdevOhLr8sfxu0N0uLWuh1YlUHjM2B8h/d1WZo6PLhMtztOy61r2gkXkHci4ikQjyu/VKi32w6uoRSdsHC5753+ZhHywpKcHYsWMxatQo5+/l5eWora11/t63b1/06tULGzZsOOW+ampqEIlEnB9CCCGEpIfArzZPPvkkXnvtNWzeHPsv3MrKSrRt2xYdOnRw/l5QUIDKyspT7q+0tBQ/+Un8gjKEEEJIRtBMew0EehHYu3cvbr75ZqxatQrt2iXWiUozZ84czJ49O/p7JBJBz549cfB3LyPUulWD28vwf3he/O53OuRa/E7977pzmuxc11D43UKGTq2KiLoqmu5WJtFWHIlVEU+HFaW9SHfx09UEJd6O16OfHze6uOl9amQYXR9vsHH8nbvqw4O9Ed96FFb7kNXUNDESh6j0p8Phk+68Lvo5USlA0xhyQEMdFhPZzgqpWx0bAWCs2K8eiwzHT7rT3c7v86fDz9bzFmRu5L0537gvNVIq0OF+iZaXNPK51dUbJdqeLOcjiM1WrtuQvCblTv3dozs+xhtbjBQhwv9aFpRjk7JI7dFac5ypJUnXADLzRSCQNFBeXo733nsPgwcPRuvWrdG6dWusXbsWDz30EFq3bo2CggIcPXoUBw4ccLarqqpCYWHhKfeZnZ2NcDjs/BBCCCEkPQSKCFx++eV48003KeTrX/86+vbtix/96Efo2bMn2rRpg9WrV2PChAkAgIqKCuzZswfDhg1L3agJIYSQNBMKhRBKIiKQzLaNSaAXgby8PPTv74an2rdvj86dO0f/Pm3aNMyePRudOnVCOBzGjBkzMGzYMF+OAUIIISRjYY6APx544AFkZWVhwoQJqKmpwZgxY7BgQfDyj7kX9EBu2zYAXK3KsgQG0VqlLqn1bKn99jFK8wZBd0aU2qPOCZB2RsvOo3VRPTdSz4fSwaWGqTVyOTdW97vJz7i/y/nX18LqlGhdN52/IedGn6+8jrqrm7RJaU1Ya8ZyzvUcSzul7pooj9nUNsBEyxbr7eQ5BumEqK+37+dWzdtEMd/rlX3PyllY1L5H/XYqX0Ab+2QeiLVPjTxH51mDe05aP5f3l84z6bPM/YeWnH+tmct7cZDovAm4NtsYrV88C/q6ROb97JTjBGLzEKwukhOFJVcfX9p1ta1aP2+SeB1ED9XVxd2G+CPpF4E1a9Y4v7dr1w7z58/H/Pnzk901IYQQkjFQGiCEEEJaMpQG0suhN99F1if2QSusaNnU/IZnLRuYtt5YFbx0WEuGsnQYzbIIWlXYpC1Hh+qg9qktgxK/VQF12Ny6FqnCqrQnQ9c6NC+rmxWra9pxVf1+9Hb5S8uc3wcJ+6hGdkbU1q8goXOJ1ZnQCT+rrolBjpdoiD/Rc9LX7UPHPuguswiJkPeC3fcYa7oE6T64XsgI16pl1rXx+/1ihdgHqqqa+pmWpld9LaSFTncX1TKCRFcrtcYq0d9h8rvQ2k5/h8rn1Ko6qA2/cl253sdeGqUBNh0ihBBCSHMjYyMChBBCSCYRygohlER4P5ltGxO+CBBCCCF+aKbSQELdBxuTU3UflLYsaW8BXM1YW5ZkWVlLE9fbSe3PKusJuHY23TlLlkfV+ppVxjdRUmVZ05Y9id/cAr0Pazu9rtTe9TzJcwxyDGf/hpVRHyNRrHuqJSLno2T0dGeZXz1fP4tWDoxEl4kOYhG0sPI3pA5vfS9ojdzqMKhzmaxcGrnsiau/7yyTZbJ1d1GZ1yTXA+zcAsvKqzsTyu9M3fkzpGyQEpkHoXMEbjj4QVq6D/7z+lEIt02i++DRYzhzyUu+x1paWoply5Zh+/btyMnJwfDhw3HvvfeiT58+0XW+/e1v46WXXsK+ffuQm5sbXaevMZca5ggQQgghfshCvXMgoZ9gh1u7di1KSkqwceNGrFq1CrW1tRg9ejQOHz4cXWfIkCFYvHgx3nrrLfzhD3+A53kYPXo0jh8/7vs4lAYIIYQQH6S7jsDKlSud35csWYKuXbuivLwcI0eOBADccMMN0eVnn3027rrrLgwYMAB/+9vfcM455/g6Tsa+CHhvvw3vk8qCUhrQISgrzOs3VGx1ppNhrFMhx6NDnqXD648/Q1VFe1h81qEzGfbS4TldMc9Chs515TOJDivKdfUyy04l0ct0WHfwb+6qP54KqzrHU/ZQOR59DL9jaywJpc6wHR4T56htYFaI93ShobC9fMYe+fm3EjqGXykAcEPqDUkBVvjdQsqU0vKrj6nnRj7T2vanv8/keLTtVW4ru6kC+j5ypQEpB1gdBi0pAHAlDy3ZSkvuQLg221Y/rLeBNtR9MR5Sao0crcUNi15KaD9NRSQScX7Pzs5GdnZ2g9tVV5+waXbq1OmUyw8fPozFixejqKgIPXv29D0eSgOEEEKIH5KSBeqLEfXs2RP5+fnRn9LS0gYPXVdXh5kzZ2LEiBExPX8WLFiA3Nxc5Obm4sUXX8SqVavQtm1b36eVsREBQgghJKNIkWtg7969TrKgn2hASUkJtm7divXr18csmzJlCj7/+c9j//79uP/++/Hv//7v+N///V+0axe/d4OELwKEEEJIGgmHw4EcDtOnT8eKFSuwbt069OjRI2b5ycjCueeei4svvhgdO3bE8uXLMXny5FPsLZbT7kXA0pOtcqhae5OWmiB2Il2OtmzA5dHPMicAcDW0abtcLVBqczoPYdDE+vKgWqfLE246fb66HLDTSc7Q060uekEsepbWbem7tUaZ6CCabaIEsfpZpXplodOYZWIe9ZxmUl6AdQ9Z90JD+r3Mn0i0bHEQgtw3ftetVeWApfXN6mioy4TL56v4HXsu5P2mrX4yz8aySGqLorTo7VT5STJ/QX/3aKufHI/V7TJLza++xyS/N+yLEmlzbM4lhj3Pw4wZM7B8+XKsWbMGRUVFvrbxPA81NTW+j3PavQgQQgghTUOSLwIItm1JSQnKysrw3HPPIS8vD5WVlQBORABycnLwzjvv4KmnnsLo0aPRpUsXvPvuu7jnnnuQk5ODL3zhC76Pw2RBQgghxA9ZWcn/BGDhwoWorq7GpZdeim7dukV/nnrqRBS1Xbt2+POf/4wvfOEL6N27NyZOnIi8vDz85S9/QdeuXX0fhxEBQgghJANpqPBv9+7d8T//8z9JHydjSwy/3KU7cj95e5L6o9bppCdV+7ilFqm9vFKX17UJrBLDur2no6drz7sYm8wlAFzd0Dp+kPKoloYbpPyyPGaqyrFqpPapvdTyeidaRrg5kkzZYr9tiPXzZZGq/A0rt8RvbYimwKphInOQtEYvv3u0tq71c12fQGLVKpDoNsDWdXPbRbtjS9WzKL9f9NxI9HetzFmQ85LOEsP/+s4XEc5uk/h+amrRecGKRh1rIjAiQAghhPihmTYdYo4AIYQQ0oLJ2IjAp4rCCLduFfN33Y1OlqjU5TJ3ipCztNporDK2lhQA2DZEGa7TFqJ4+wCAiVPqx6NDZ1YYVYfqZDjYsgJpe8/kZ+6PO1ZJMqFCGZ4uNqx2VqliHfKUpZqDlKMNgpzTrC9fH3c961okap9LxmZoWUKrp1wT/aztsVJui5Fw/N0mDWKdV6bJARI5bh3S9yvhWZIl4JbSlVIj4Mo41rOgvzOt7qLS8txYHTSlLVB2YgRsqWDnruqUHD8pmmlEIGNfBAghhJCMIoHM/5jtM5DMHBUhhBBC0gIjAoQQQogfKA1kBrrdpUSXwByytF5Tkzoo4Gphlp1G67naXiX1dMvCo5F6l1UCtPgdd2zyGLr1qLYaXvPG6ujn3ystbqKYq0QtipZ+a5UtBmy90bGBKh1eaqFa+xwi25s2kr6ZqL6fjrK6FrI0t36GHF1YtYWVczxwvFviNlFrn9W+Och9kkloC7CFnCvZchwAile5uS1OToqa43KRB6BzC+Rzoq+3/J6U3zWAe5/Glpseh0TQz6I8Rh+Vd7Ko6kD08zS1TJ5jxVZ33tJGM30RCCQNzJs3D6FQyPnpK26kI0eOoKSkBJ07d0Zubi4mTJiAqqr4PlhCCCGENC2BcwTOP/987N+/P/ojWyLOmjULL7zwAp5++mmsXbsW+/btw/jx41M6YEIIIaRJOBkRSOYnAwksDbRu3RqFhYUxf6+ursaiRYtQVlaGyy67DACwePFi9OvXDxs3bsTFF18c6Djv7I7UVxYUf9fhfxm61J0JdaVBieycZdlpdDhOywgyVK/DfJLeyI+7TFsULVuaXFeH36XVCAC8Ha9HP2tLoAy5Tv7Ks84yK+TrNwSsQ7w6PJm/NH5lQ7ltzPH9dopTIe7GqMJ3OiHto/oZsrCkuPh3tE3MNU3BHKeq6l2i+0lUwmiocqe0141VFmBpkdVhfBnylx39ALdTYljJa9b32ST1feP32Xjtuh87v0tJU0oBGr3sYSF39BHS6qG6OuCgr6EkD10DJ9ixYwe6d++OT33qU5gyZQr27Dnh3S8vL0dtbS1GjRoVXbdv377o1asXNmzYEHd/NTU1iEQizg8hhBCScTTTiECgF4GhQ4diyZIlWLlyJRYuXIjdu3fjs5/9LA4ePIjKykq0bdsWHTp0cLYpKCiItk48FaWlpcjPz4/+9OzZM6ETIYQQQkhwAkkDV155ZfTzhRdeiKFDh+Kss87Cb3/7W+Tk5CQ0gDlz5mD27NnR3yORCF8GCCGEZB4hJOkaSNlIUkpS9sEOHTrg05/+NHbu3InPf/7zOHr0KA4cOOBEBaqqqk6ZU3CS7OxsZGdnx/z9wi+fh3Db2C5P2u4itblao/tf7nbX+pRo5zRtH7TyAqQtUHcRk8s8NTZ8Of7xZOlQjc41kNqcHufkZ069HuDadHTJWYmln8ruhkAwe5W8plbXxPKLxzrLhmz8ff0+VDlWC8velMlY3e+0Zi3tqgPd9I2Y8tMSmU8QZE51/orfOQ2i0btW2tSUlG6MksbaViyfRWnxBdy8Hr3cynnS5XeLDYuefBZ1eXOZszC4AQuuzEvQeSfyXtHX5sb2PfTwfSFLyMs5/NirS2h/CUH7YCyHDh3Crl270K1bNwwZMgRt2rTB6tX1N25FRQX27NmDYcOGJT1QQgghhKSeQBGB73//+/jSl76Es846C/v27cPcuXPRqlUrTJ48Gfn5+Zg2bRpmz56NTp06IRwOY8aMGRg2bFhgxwAhhBCSaYSyshBKIvM/mW0bk0AvAu+++y4mT56Mf/3rX+jSpQsuueQSbNy4EV26dAEAPPDAA8jKysKECRNQU1ODMWPGYMGCBQkN7K/Pb6u3Dwrnmw4xyvColg2cimVqWaLoDocWMnyluw9KqWCg2s7qqCjD/7obnEaG+AfPc7svyo6LuiqZrDSXaKg0ZrsA+5HXVFtCpZ1NhxwTtf3pYyDN1ewSrdCn5RdZPTNvontO8r6JkaIEpj1Xha3l8WMkHGUfleeoQ9zyWgU5f91xTyLvhXRUK9T2PfkMSckKAIqNc3x8wOXO79L2q+9peY76GZbVQrV12OqYKgkyT5bUqmWiRMmI7oNINvM/M6WBQC8CTz75pLm8Xbt2mD9/PubPn5/UoAghhBCSHk67XgOEEEJIk9BMkwX5IkAIIYT4gS8CTYel/Up9V2uB1jJpy7P0La39JYplH9RYeQiy5KjW9ywtTlsELX3V0tcT1bM1UlO2tEhdblnbQCX6GvvFKumcDuQ8Wp0gG+qoKPV9bQnsuOp/o5/1fZJrjM26NlYuh6VnB8nfsO6TIF1DGxs5v4B7vjonokyWDVa5Qxpr/mWuhz5+H2FZ1NdG5wzEo6EOon5LVcc8Xzf/0td2GpkHMWhifc5L5EgNbvhRYrlo5ASnxYsAIYQQ0uQ0014DfBEghBBC/EBpIL18qiiMcOtWAFyb3GAd/jdCxTIkdVx1o/NrJ9ThryD2QQun0p+qFmiF7WXXLV2xTNvpZBhdV5OTSKsT4Iantfwgw5y6a5wVxm0orB1vmZZmZKXDIPtMRUfFIOiKkDKMbYX/LRqaU4muCCnX1bKBXFeHsR17rrIryvHUqWWpqtYojxGkWqFcV4emG+N6W+j5nnTndfVjUWMbKyyggHv/x1QPFM+7nhvLLipD/JbtUcqQJ8bt/OpWOV3mfhfJ4+vrVDq8fv7n/CV+51eNc/7i/wmH6lhZMFkyM05BCCGEkLSQsREBQgghJKNophEBvggQQgghfmCyYHrJvaAHck92H9xVr0VpndSyCEq01UhqY/lKpzNLE6scgYlTiqOfrfwBre/JHIFFVQecZZcs3XzK/evtpNYInKLMqxi71imlnXFsf3fc0tqo9WzZDc3S6IPo9xqZ+zD4N25pZFnm1tLM9bJEdWGrG56lWVvWNisnQF9vC20Lk5pt8f3u+X6otGeJLE2sc2Kkvtzm/nFx95FMSWmJnm+JlROg77emtoTK503n4Mjz0Lkk2gYonz9tNbQs0Dp/SDLQOJ7cp/5+0SXNZelkqzOixjmm0YnwkrBrsQ5SXp0EI2NfBAghhJCMgtIAIYQQ0oLhi0B6CX360wi1y/7ktzfjrmeF2YD6Dnt1KsQvw1N6O9lxraGKfH7thLo7WG+I31XVQRkCTMauKO0/upKh/F2Hgyu2bj7leoAbAmys6m2WLcpBhf8bo6ucDnlKy6QOPydadVHKVDpUK2nIPtdbyT/xsKpcymsPuPdinnF8S0JJBrkf6xiNce2DoMcmnylt35PrPrXUne9Jfd05duZfSQzW8eX3lP4OkTKClhBk+F3LiRop0+lQvWWBljLOIz//lrNM7sfah/w+jRw7Dry/zxwrscnYFwFCCCEko2CyICGEENKCCSFJaSBlI0kpmfl6QgghhJC0kLERgdCnL0CofQ4AoPc58bUqxzL2w3ucZZbVTOqtWkMLfyX+uPS60hao8wAkWoeXVi9tGZM5CpF5P3OWSfueRp6vxilpDFcX3rIsfg6G3k6W9ixWGrVl/UqUpi4Pq3V4v10rgyC1X6tscEM5GX67wWnkNdY5IXJsMZYt8bu0kjWEletgXV+dr+EJu3Ci9lTALofrF/3s+b03tJVU5yuF5/0g+lmXeC6/eGz080BdRljkRIXF98mpjiGR569t1RUDLndXFvlT1vU3c1vUdSv2mXcSllbxwx8DX7kp7vFTSpqTBUtLS7Fs2TJs374dOTk5GD58OO6991706dMHAPDBBx9g7ty5+OMf/4g9e/agS5cuGDduHO68807k58f//5GGEQFCCCHEDydfBJL5CcDatWtRUlKCjRs3YtWqVaitrcXo0aNx+PBhAMC+ffuwb98+3H///di6dSuWLFmClStXYtq0aYGOk7ERAUIIISSjCCWZLBgKtu3KlSud35csWYKuXbuivLwcI0eORP/+/fG73/0uuvycc87B3XffjWuvvRbHjh1D69b+/hfPFwFCCCEkjUQiEef37OxsZGdnx1m7nurqE1J0p06dzHXC4bDvlwAgg18E6lb8DnWflBiW3mqtG1llbaWmZmmv2p9bt7y+/oDW/XWpYCsvwMkfgJtbID2yWUrD03kB8ahWZWO1B/3a/eOin7UuKHMdtK9cojVjy9vr6LtK69U6odT+LV043TkBDSHzUPS9KFtda4040RoDFvqelte0o1pX5g+MVcukDq3LFsu8AH2vy302dE5yrFqH9zs31nZB0O26zVoVPtH1RhZVTY9+fvTwu773EyTvRJ6Hfr7lXOmcCFm2u9goTazzg3Quk7zfylT+gFxX10qYKPIZdF6X3I9Vilt+R0aOHY+7XspJUY5Az549nT/PnTsX8+bNMzetq6vDzJkzMWLECPTvf+oS0v/85z9x55134oYbbgg0rIx9ESCEEEIyihS9COzduxfhcDj6Zz/RgJKSEmzduhXr168/5fJIJIKxY8fivPPOa/ClQsMXAUIIISSNhMNh50WgIaZPn44VK1Zg3bp16NEjtlHTwYMHccUVVyAvLw/Lly9HmzZtAo3ntHsRiLHIGWVGZThSh9tl2F6HCi17jSUVxJQRNmQDae/RYUWJZcvRdjFZqhawy9VKdFczKRXoc5Bhvmvv97V7ALGdIWUYXZd/bozStakKzcttdchV30d+scpky1CxJa8A/q+3Rt5/MecgpIEYCUl2sHzD/73XUKlkv/i9jg1JCKkola3nzeohKZ9T/QxbXVK1hClLAOcvLXOWyXOW9kgAeOLq70c/S/kQsK2k2mZsrSu/J/T3i/zukxbIU+1HEq8Ud6tIBOiWmHU2MKGswAl/MdsHwPM8zJgxA8uXL8eaNWtQVFQUs04kEsGYMWOQnZ2N559/Hu3axZ/DeJx2LwKEEEJIk5AVOvGTzPYBKCkpQVlZGZ577jnk5eWhsrISAJCfn4+cnBxEIhGMHj0aH330ER5//HFEIpFoImKXLl3QqlUrf8MKdhbAP/7xD1x77bXo3LkzcnJycMEFF+DVV1+NLvc8D7fffju6deuGnJwcjBo1Cjt27Ah6GEIIIaRFs3DhQlRXV+PSSy9Ft27doj9PPXUiof21117Dpk2b8Oabb6J3797OOnv37vV9nEARgQ8//BAjRozA5z73Obz44ovo0qULduzYgY4d63OU77vvPjz00EP41a9+haKiItx2220YM2YMtm3bllDIghBCCMkImkAasLj00ksbXMcPgV4E7r33XvTs2ROLFy+O/k1qFp7n4cEHH8SPf/xjXHXVVQCAX//61ygoKMCzzz6LSZMm+T5W1hcnIOuTEsMSKw9A2v40WosbPG9inDVdGrILWnkA5n6EvteQRVEiNTW/xwaCtTOWmp5lLbRyErQuq3M0LA3ZKRutrneiWn9j2BB1/obMGdCasbxuetmN7euTf+b/8RFnmZWv4vWLny9jobXmCtGu2tK2tX4bJCfCuY4pahks7xNtewwlWH7YOkai+9E5EfL5HjSxr149LrotsHw2J+nyy8Kip+3J1+6PnxNh2f60na+iW+/oZ52HIHNLrO8+bUe2yp3L/QwW1+X44Y/jbpNy0lxiOF0Eej15/vnncdFFF+GrX/0qunbtikGDBuGXv/xldPnu3btRWVmJUaNGRf+Wn5+PoUOHYsOGDafcZ01NTVTXkPoGIYQQQhqfQC8C77zzDhYuXIhzzz0Xf/jDH3DTTTfhu9/9Ln71q18BQDSRoaDAzRItKCiILtOUlpYiPz8/+qMLLRBCCCEZQVZW8j8ZSCBpoK6uDhdddBF++tOfAgAGDRqErVu34tFHH8XUqVMTGsCcOXMwe/bs6O+RSAQ9e/ZEq4vHotUpfJa6mpoM11nhQW3R6+2Eyt2wtQy5WdW0GsKyFsplOjw20HXaOcj9JNptTmPux5AGdIjZCqNqe5OF3DZVVrPGQMsNMlSurYVWRUar8pysjqmts4nKJp7qVGdVcJMWsUl3Xhd3Pc3jImwMuHKTvhessVr3lHOfqGc/SBjf7zESRds8ZWhcd3SUMo1GSzPyd70f+T1hdVPV0pMMzVtWvoaQ1Qv1d6/8ftVVLi2kFCerqp6OlQUzjUCvJ926dcN5553n/K1fv37Ys+fE/xwLCwsBAFVV7v88qqqqoss02dnZ0eIKQYssEEIIIWnjZLJgMj8ZSKBRjRgxAhUVFc7f3n77bZx11ok3+qKiIhQWFmL16tXR5ZFIBJs2bcKwYcNSMFxCCCGEpJJA0sCsWbMwfPhw/PSnP8W///u/45VXXsFjjz2Gxx57DAAQCoUwc+ZM3HXXXTj33HOj9sHu3btj3LhxjTF+QgghJD2EkKQ0kLKRpJSQF9CEuGLFCsyZMwc7duxAUVERZs+ejW9961vR5Z7nYe7cuXjsscdw4MABXHLJJViwYAE+/elP+9p/JBJBfn4+qvfvicoEUseyunNZ3eC0LcV310Cj4xrg5gwE6Uwo0Vq71Pu0RU1qY9ZYANuiJ+c0SN6D7FQn9WvA1bB1vkYmafvpwOr2aJWK1dcsXllVwO6oqbVtmWuhcwQk+vnSWr9E5hYkk8sht03HfZKqstUSqzS0Rs6pzs/Q11Tq69Z3gc67kN99lg3QOp7GynOxbH+WRVB3OJR5CdqSqHMtThI5eAgdLxwRbb/bGJz8/9IHC29FOIncicjHR9Dppp826lgTIXCJ4S9+8Yv44he/GHd5KBTCHXfcgTvuuCOpgRFCCCGk8WGvAUIIIcQPzdQ1cFq8CMgwm2Uf1CE+T1bUUqErK/xvYYXRg+xHoiuGWTjjxp64ywCgtxFWludhSQw6jC0rBMaEqoU0oLsNWhXyrDCqJlVdBBsbXc0Nb9TPmw5x/l5Yxq5V+5FSjA5p6zCutGzp50TeY9b17oj4WDbDhkL6Mhx+zRurnWWpkAOChPsb474Jcg/rbnwSLekUv1P/u/UM6S5+Em2BlnOl5T2ro2BIdQmV95S89wD33tTbDRxf/3mQWvakrHKp7gspkwWxI6eUNJcYTheZOSpCCCGEpIXTIiJACCGENDmhJNsQUxoghBBCTmOaqTRwWrwISE0rSMlP2Y3QKvGrtXVpd9E5AVpflVaYIDkCVo6CtgnFG5u27Gibjl9boF5P5gVY+QtBLGOpshNKW5QuLCp1+VSUhk0lUuvXeq5fzViXLdb2KjnH8axWQOz1lvei7j5o5QUEQe5Hn4cszRwklyTdtkON33wVnb+QaGlwrbXnis9Dlt4TdztpMwSAnUYnTMtaqJ+pnbKssS7hLr7TrO+QnYbt0OraWCzmO6sdG9Uly2nxIkAIIYQ0OXQNEEIIIS0YSgNNh9N1LYD1R4aK8/q6oenePu1UGr3MkhEkWn6QYV0rdFb3/BLnd0uK0KFDGXbTIV+JPn9rPFI20OFYGUqUsgxwCjtdgvi1kjY1eizS+qSRc6or+clQrbZo6TmWZH3TfU7kNdb3qZYYUoF1bXQ4WoeAJfJ+L77fXSbvv2QqG8rQvdXh0bo2Gnm9k+kSKsemz0lKLMX3u9dbSiza5iu/C/S8SUlJWoUBYOcuN/wvJUw9b09c/f3oZ21ftKoQSvT5DjQqYqaNrCSTBZPZthHJzNcTQgghhKSF0yIiQAghhDQ5zBEghBBCWjDMEWg6EtV+5Xa1Sk+1LHJS02vIomdp9pZFUHaA02V8JUHKD+t15TG1ZctvboPG0l6lZq3zFbR9EAleU0cXNbTmpu52qLVXeU8NdKsvO2PtIy1ZCv0caA1XoksFy+uRq5ZZVsNESzpbz6y2CMr8ET2WgbgFfrDOocFt5XkZ52jlBOjnS+Zz6Osky+MG6Vqokc9wjNZv5ORYORkyJ0nnFgzU103M1YfKoqjLE0tkToruhFkc4PxJ6jgtXgQIIYSQJqeZJgvyRYAQQgjxQyiUpDTAF4EmRYfcZEgupmuf6OqnpQCNFWK3Kg0666rt8o0wp9xnzLiN4+ll8vjWfizZIqbjmwhHaikg0dCtDmPKsKK25FkdyTKpa6EO21rV9D40OkhqpPyg7XsSfZ+GhfWr9j43FB8kVG0hxxNkn1nGun6vaTpsptraKeUAfV/K8LuW0Kwuipa8p7t9WtZSCylTWfcl4D6b+nuiQnQR1J1fZffBp5ZudpZdqyyiklTdiySWFvMiQAghhCQFXQOEEEJIC6aZugYyc1SEEEIISQsZGxGoe38v6o6cMDo1hqYr9abeb8Qv/6rR+q5fLc7qcNjqh/E7h+ntKj4+Ev1sWXQa2o8kSNdEiS5/PP3mX0Y/Tyvo4CyLKQ8qztm6vtoGKDXUjgEsgmUDLo9+1iVPG0N71DkRsuSvpd/rzpPH74tvn9M5A04+wf3j3P2Iz+F58fVkfX/LfVr5Iqmya1oaucbv90I6Sk/LkrqAazXUuSzWPOpz0vMhkdblCq21W+WPRf6CHovMC9D3gj7HSXfW5zdo67L8btIlrF83LLIZD10DhBBCSAummUoDfBEghBBC/NBMkwUz8/WEEEIIIWkh5Hme19SDkEQiEeTn56N6/x6Ew+FA22rtVXrZLQ1TbxfEgyt9wNrLa5WAtTzvsqVnbCvQeq2/oRoHfksHW3UEdAtVq3Tpa9fVtym1yi0DbnngVPm8LV+5VZq4MbC0XZ1bIdH3qdSXrTnUaA+4vE91WVd5TbUOPLZ/QfSzvmctjT6ZtsCpJkj7YI1138hrXH7xWGeZ1dpZ6uk6z0PfN3K/ep9yTnWNAYmucSCfLytfpaHn0Go7Lu83q4R7TClycW/q3B15jvL7JXLsOM5+ZTuqq6sD/z/DLyf/v/TBMw8hfEZO4vv56GN0uvq7jTrWRAgUETj77LMRCoVifkpKSgAAR44cQUlJCTp37ozc3FxMmDABVVVVjTJwQgghJL2E6uWBRH7QDKSBzZs3Y//+/dGfVatWAQC++tUTb6qzZs3CCy+8gKeffhpr167Fvn37MH78eGuXhBBCCGlCkpIGZs6ciRUrVmDHjh2IRCLo0qULysrKcPXVVwMAtm/fjn79+mHDhg24+OKLfe3zVNKADEHpcJVVBrOVYVGT2+kQlAxr6jCqDmVJOcCSAjQyzKrHpkOZftF2Qqv8sdUZ0ULajbQNSO4nPO8HzrJY69F10c+Jho1jOq4ZZYzTXVZYh3ilDVBfC3kvaNlA338W8j62ZDKNPIaWfpwOnobckKrQv76mVqjYQo5VdxDVyJC/Pr5ldZMh9yDdFhO1q1olhjXynC1JQZ9vENnCrzT0pCw3rNAyjdxOP8/y2ZDnEIlEkN+tV3qkgd/NR7h9EtLA4Y/RaUKJ77GWlpZi2bJl2L59O3JycjB8+HDce++96NOnT3Sdxx57DGVlZXjttddw8OBBfPjhh+jQoUOgcSWcLHj06FE8/vjj+MY3voFQKITy8nLU1tZi1KhR0XX69u2LXr16YcOGDXH3U1NTg0gk4vwQQgghGUcyskACjoO1a9eipKQEGzduxKpVq1BbW4vRo0fj8OHD0XU++ugjXHHFFbj11lsTPq2E7YPPPvssDhw4gOuvvx4AUFlZibZt28a8iRQUFKCysjLufkpLS/GTn8Tvj00IIYS0RFauXOn8vmTJEnTt2hXl5eUYOXIkgBOReQBYs2ZNwsdJOCKwaNEiXHnllejevXvCBweAOXPmoLq6Ovqzd+/epPZHCCGENApZWcn/ADFR8JqaGl+Hr64+Iel26tQppaeVUETg73//O1566SUsW7Ys+rfCwkIcPXoUBw4ccKICVVVVKCwsjLuv7OxsZGdnJzIMAInrbXI7y+qkdblBt7o5AjIvQFvtpIZq6dc6J0Bq/bKksF7WkEVP6oQxeuIuWzc9iT4nqffpnASZFxA6d5CzzNICgyB1yhgbnCjja9n3gpCotVFaQAFXew9/xV3XKn8sr6G2gVm6vx7nh0b+irxuejvLXubsP0AZXX3d5Dnqe1ref/pZtOyTcr6HqBLelrZt5b1YbaCtXKUg31H6u0Dm0uhnWM6bLk2NZfW2wxj7nph/Xd68o5E/oOdGzr+VP6K/J/R3Wjx0ee2MaEOcooJCPXv2dP48d+5czJs3z9y0rq4OM2fOxIgRI9C/f3y7aCIk9CKwePFidO3aFWPH1t9sQ4YMQZs2bbB69WpMmDABAFBRUYE9e/Zg2LBhqRktIYQQcpqzd+9eJ1nQzz+GS0pKsHXrVqxfvz7l4wn8IlBXV4fFixdj6tSpaN26fvP8/HxMmzYNs2fPRqdOnRAOhzFjxgwMGzbMt2OAEEIIyVhCoSR7DZyICITD4UAOh+nTp2PFihVYt24devTokfjx4xD4ReCll17Cnj178I1vfCNm2QMPPICsrCxMmDABNTU1GDNmDBYsWJD0IGXVKh26tKrJWcsklvVJh+N0mM2p9OcWFnTkAH18GXaT1dsA4Pdb64swmWG1Xe7xrLCqZaGyOiPq7SzZAqKyoJ63LBXWS9RuJi1E2upmWUIlQareJVrl0NpOh9Hl9X9KdZGToWFtT5X3CQCMNaxf2s4ZDy0FyGPqfQSpwKmfm3jE3IsD4qyoSMaiJyUHy2qncWzNAeQHuV1DEpolOcix6uPJ74KYDqZb6++xa++Pb6vWz75fKQawZSuJlkLks2BJMfK7ve7gIV/HSglp7jXgeR5mzJiB5cuXY82aNSgqKkr82AaBXwRGjx6NeKUH2rVrh/nz52P+/PlJD4wQQgjJKNLcfbCkpARlZWV47rnnkJeXF3Xg5efnIyfnRD2DyspKVFZWYufOEy+Rb775JvLy8tCrVy/fSYVsOkQIIYRkIAsXLkR1dTUuvfRSdOvWLfrz1FP10bhHH30UgwYNwre+9S0AwMiRIzFo0CA8//zzvo/DNsSEEEKIH7JCJ36S2T4Afgr/zps3r0HHQUOcFi8CVsc7qY0FsYxZOp3cp9Zstb1r0MTEyqy6nQpd/VTqhFZXsUVVB5zfL9FWQ6UNSqzypDMW/jn6eVpBB2eZzAsw8xeUvjj4K886vyeqvcs51nq2X124KTrhybwEba2TtqzJz9wfdx87d7n2LX0v+u02qZFavy6hLceqcwLk8Sy7IGCXtJb3jcyJANzcAm2flM97ouWHAf/dKPX3SyryR/Q+9dhk7oHucDhkY/0ynfci80diui2K7xTrO7Oh7qZWKXZ5PYr18yaOr22P0kp7rTqezCeQ9/7xo7XmOFNKmqWBdJGZoyKEEEJIWjgtIgKEEEJIk5Nm10C6SKr7YGNwqu6DFpZFUIaOtZ0lYfuaCqXpCnJ+8RtW1KHChjqpxUNb/XRYOd4xrGpmWpqQMkJDVQ/9hmODIEOVTRH+t5DXUYfwrY6OlpVTI6/Ho4ffdZZZz4llUZPoceuqk363taQRHf63kFKF1TVR2zUtW5p5PKP7nmW71PPkt6ppQ/i1y1pjs9DzpOdRnpfVlVVLr/I7U38vu5Kpixy3/D45VFeHy97fl57ug3/4DcLtz0h8P4c/Qqcx1zXqWBOB0gAhhBDSgqE0QAghhPggFAohlER4P5ltGxO+CBBCCCF+aKaugdP+RcCyD1plLq3ud5ZO9Zooowu4WreloVkWSEt71NrfQNTrgtraqJH2Pm31kxbBh2/6rLNM6rRaT5RlTS8RnQgboiErUiJYHc+sroGJdhRMBnn983/oaqbSFqZzK1zt1c0R0HkfpcPjl7WVumyt6uom0XkATgltdQ21LmyxU3ZRVM/X4N/Ef95kHoDOkZBzldfXfdZlt8UgOQH6npLPmM6rqTbyHvKXlkU/y7LYgGsJbignQo7HsjLrUr3ShhqkFLQ8D/2c6FwHtzNm/BwN/V0rrYaxeSf1y/T9JUtcDxLf7ZEjNcCPki9l35I57V8ECCGEkLTAiAAhhBDSggklWVmQOQKNg9MBTIV4s3xacXQ4UNqZdDcybXuTYT4d1pVY4WcdOrTGLUP6l4TbxV0POEV3QIEVRpZ2n5hucCKstz7i7l92DrPCiABQHL+Anu+ukUGqyUkakgLk/aDDk4naEuUx9f0mQ+5WdcCGLKByjrU51LKMyXBwrlo2eF59yFdLbbLqX0OVBa0QuLzfqqdc4yyTYX1dnVKOR9surWfRQt9T14r7VIfKO94/ztc+y5SEds0bq6Of9f2tK4nKynuT+i6Je4yx/d37Rt5vVtdKXclRyqJPXP19Z5let/c59cfUz5SUKnR3VUkQC6o8xmYh0R6qq/O9j6RpphGBzBwVIYQQQtLCaR8RIIQQQtJCM60syBcBQgghxA+hUJLSAF8EEsYqHWtq70Jr1qWAre2knmlpX4Crm+mxyXFr7c/ShXNVJy+JVcZX622W9UnmGpSqY3Q0dHkLJ59gV2KlkAE7L8Bv/kAQYuyEouObdQyt9fvNWbDKsYa/En+7yc+4Nju/pWIBV2u37Kra6iWtZ5a1EEqj17k1ln3TmmO5rj5fmZdglcKOX0w7GLpUrl909z9535QpS2BMp0CBmfP0lNvBVH5Paeuw/M7SXQMPqv1IrPLX2jotj6m/l0wrtSwb/eXrnWXyGDJf5mMvjTkCzZTT4kWAEEIIaXIoDRBCCCEtmGbqGsjY7oMfLF+IcPscAG5Izm8YEbDD/3pdiVWRUIfSZBhZVogDXOvPcVXNzQojy9Ct3k6GZ7W1T0sFiXZZk+FhXWXR6uomkRYlIDY86bf7YDLd2eJhhcb1MYMcr7FlC13VUtsJrZCvdb/JMLKWkOT9pu8nqzqnFUaWzwzg2mf1uOXxZbU+wH02rM6IDUk2iV5vOW/SEqjR90Jj3CcW1nOq57RswOXRz5a1EHDvRy3NSFlS2wdlhUD9XSvvG/0dIZ8FKSFEjtbizEUvpaX74IfrnkU4t33i+zl0GB1Hjsu47oOMCBBCCCF+yEqyoFAy2zYifBEghBBC/NBMpYHMHBUhhBBC0kLG5ghU79+TUg1F68JSX7Q0y4aQuplfKyPgaoNWvoK2LFlWL43UW3WJX5m/kKhOqe1z8hiNka8ApCZHoCHkeWl92VqWCvS9IK+/zhfRWOWJ/Wrm+ppaVkpLe5Y6MJCeDo+NTUO5JYkQ5P4OYle1rNNS29fPqcw70hq9ZZ3W+QO6PLHEbylybUGVyO+2tOYI/GVF8jkCw7/IHAFCCCHktITSAHD8+HHcdtttKCoqQk5ODs455xzceeedkEEFz/Nw++23o1u3bsjJycGoUaOwY8eOlA+cEEIISSsn6wgk85OBBHoRuPfee7Fw4UI88sgjeOutt3Dvvffivvvuw8MPPxxd57777sNDDz2ERx99FJs2bUL79u0xZswYHDkSvxMeIYQQQpqGQDkCX/ziF1FQUIBFixZF/zZhwgTk5OTg8ccfh+d56N69O773ve/h+98/oRFVV1ejoKAAS5YswaRJk2L2WVNTg5qamujvkUgEPXv2THmOgEZqaDonQHpZpZYOBCtVLPVe7ZeVWmCqdHAr10CP09J3pW6n8yekLqr3IUu+as2wofFIGsNnLfVV65w0QWpTJEqiPnY9NqnLag+4LB2r609ILVjfi1arYTlWfS9YOSF6XXm/WfUAdGlav9eioVwav3Ouy4RLLD3b2n9DeQdWm3W/9T70uGVegFXvQ+cP6LwPWfLZaif8lCqjLI9htUrXNS0k8j6JHDuOs1/Znp4cgU1/SD5HYOiYjMsRCBQRGD58OFavXo23334bAPDGG29g/fr1uPLKKwEAu3fvRmVlJUaNGhXdJj8/H0OHDsWGDRtOuc/S0lLk5+dHf3r27JnouRBCCCGNR1ZW8j8ZSKBkwVtuuQWRSAR9+/ZFq1atcPz4cdx9992YMmUKAKCyshIAUFDgVpIqKCiILtPMmTMHs2fPjv5+MiJACCGEkMYn0IvAb3/7WyxduhRlZWU4//zzsWXLFsycORPdu3fH1KlTExpAdnY2srOzzXUaI1QsQ/w6rFb8jhEqDHB8p+OaWiZDfg2F0SXWXASxL+owXzyqp1zj/C5Dvjr867c0MgDkL623xWlpJlFpxLJMBbH6WeFYiRVyDnIOftfVUoAsswq4odTYeyp+VzlrLGGxn5juf+JzQ/ZQK6zeG/X3hg4xSxknxtZrlLuWMkaQcssW2k4nj9lYtlZ9jSWyO59+Tnfu+jHiIedmrLovpFTQG640oJ8FeYyB6naTHTaveSO+PbtYzZucU0tS0LJFugiFQgglkfCXzLaNSaAXgR/84Ae45ZZbolr/BRdcgL///e8oLS3F1KlTUVhYCACoqqpCt27dottVVVVh4MCBqRs1IYQQkm5CoSTtg5n5IhDojD766CNkKY2jVatWqKs70Q+6qKgIhYWFWL26vgFHJBLBpk2bMGzYsBQMlxBCCCGpJFBE4Etf+hLuvvtu9OrVC+effz5ef/11/Nd//Re+8Y1vADgR9pg5cybuuusunHvuuSgqKsJtt92G7t27Y9y4cY0xfkIIISQ9JFsLIEMjAoFeBB5++GHcdttt+M53voP33nsP3bt3x7e//W3cfvvt0XV++MMf4vDhw7jhhhtw4MABXHLJJVi5ciXatWtn7NlGakpZKdCPAVfviik/HMCG5zd/QW+XaKnaIDkS1tiOCw3dGpvWbK1xWyVwtfVMjsdTerZfjV6j8wJSgaXLBynVm6guLbV1rVFr/Vgew7qHtdXtxvY9op8fPfxu3GPolrXWvaifKcng39zl/O7k0qhxy2uqdXg9x+lG2jArApTfte6hmBLTQmvXZXuv3b/zlGMBXFuebgEuLaExuR3iPPQzq/Mw5DEGiXEC7jnr+0R+h+vzleOZaJQwl7kMH3t1SB9JVhbM0PY+gV4E8vLy8OCDD+LBBx+Mu04oFMIdd9yBO+64I9mxEUIIIaSRYa8BQgghxA+UBpoOGWayKnFZFfqsMKYOKcuObw2F4v3KFjHjVlXSGgO/UoVlg9Pj1BUSJVbnMGl10ugQt2MZU/Pm2D5VZ8ZUVf2z9uP3GKnqTKjlAL/H0OPMmxj/ul0Sji/byWs6EEr6EddJP0P6PrHOw2/1PKhl0pKqQ+O9RfhbW2Ut2SJRe7IOv8vvLNntD7Dnwnq+rnljtfO7VR3UQobftbwiz0NKCEBsFUKJtmg6qOsm5S4tE0n0/V2xtF62kHMROXgIN1yY2FwEJtmiQAG3LS0txbJly7B9+3bk5ORg+PDhuPfee9GnT5/oOkeOHMH3vvc9PPnkk6ipqcGYMWOwYMGCmHo+5rACjYoQQghpqaS56dDatWtRUlKCjRs3YtWqVaitrcXo0aNx+PDh6DqzZs3CCy+8gKeffhpr167Fvn37MH78+EDHOS0iAoQQQkhLY+XKlc7vS5YsQdeuXVFeXo6RI0eiuroaixYtQllZGS677DIAwOLFi9GvXz9s3LgRF198sa/jMCJACCGE+CGUlfwPTtTXkT+y8Z5FdfUJCaxTp04AgPLyctTW1jr9ffr27YtevXrF7e9zKk6LiIBfXTZImc9E7XuJHlNrj+k0vARF5gXEWAuNkqcSrRnWYYl7DDFvQboBpnveUtUZ0i8NdaNLFJmjoUv1an1bYunZltZuWf1Cymrm7FPlfcTbB+Bq/71V+WO5LIhGb6HLJMv9aL1elmMeOP4CZ5lle9RzI6+bvhekva/itt84y6wOf0OM6yZze3SOgGZs/3oNWo9bXn99n8hjxF6bcdHPTW0PPSUpShbU/XTmzp2LefPmmZvW1dVh5syZGDFiBPr3P3EvVlZWom3btujQoYOzrtXf51ScFi8ChBBCSHNh7969ThvihvrtAEBJSQm2bt2K9evXp3w8fBEghBBCfBH65CeZ7YFwOOy8CDTE9OnTsWLFCqxbtw49etQXACssLMTRo0dx4MABJypQVVUV7f3jhxbzIpBopTdd+SpVFjUZ5ktVFbqEx6KsfdbYrLCuDIE2dA5+pZlEQ/Opum6JSgGPq0pzk5+539dYEq0cGXRbibxuetyyep1GHu+4sg9asoGWjZzOlCocLcPf+t6TtlNtHxwolqWqU52WFOQ9pitwyrHqELsca8zYVJfOjuL+szosBqFswOXRz9qSKKuD6kqSE1XlUOu7QD63Md8v4nfdfdCSA9ZH6iWOa+Ku1cikuY6A53mYMWMGli9fjjVr1qCoqMhZPmTIELRp0warV6/GhAkTAAAVFRXYs2dPoP4+LeZFgBBCCDmdKCkpQVlZGZ577jnk5eVFdf/8/Hzk5OQgPz8f06ZNw+zZs9GpUyeEw2HMmDEDw4YN8+0YAPgiQAghhPgjzRGBhQsXAgAuvfRS5++LFy/G9ddfDwB44IEHkJWVhQkTJjgFhYLAFwFCCCHEF6nJEfCL53kNrtOuXTvMnz8f8+fPT3RQCHl+jpRGIpEI8vPzUb1/zymTKVJl55Iapu6aJ+1sDemucjyNbS1LJVZnQqlFWjpkMudrHV/qhNpaWH7x2OhnXcb4SWGhsrTtVGHp4Nqil+jYrO6DGute9Gvfm37zL+Mu050JJVq/1pq5LFWsdXG5rh6b1Nf1vSiX6RwBiZ43PdaYDnwJoK2FVl6C07VQWfu0lVMut+4b6/i6NLDMC7C+34J81+p15bXRx5f5MrqksizNHFGW0Dl/qX/eFuyuv2aRg4fQ8cIRqK6uDpSAF4ST/186sG0zwnm5ie/n4CF0OK+4UceaCIwIEEIIIX5g0yFCCCGkBZNeZSBtnBYvAjKsGaRrnxV+dn4PUNlO09hyQKqshUHCfDJUqm14ulNgvGM0NC/WHMtz1OF3Gea1uiYGOV99DBnWN7dT4X+57nFlkUtUqhiy8ffRzw2dkzVWS+6SofKHb/qss8ySEeS9ocPr+r61KtbJY2hroZQNrEp3OhQvZQR9fS25K4jtVFotJ915nbPs9Z8+Hf2sw/+yIh92xd19zLqWtdOSybTcIKXQOnV9ZWi+oWfYkpusToWymqAet1y2qOqAefymoXm+CbDXACGEENKCOS0iAoQQQkiTwxwBQgghpAUTQpIvAikbSUo5LV4EEtXFpYabZehdqerwlio9X+qU2j4XBKmNWnYya5w6J0CWg+2o1m2M7o+eKl17TFxTrSdLW1JIdbGTc6H3aWHlIUg9VfPU0s3O79fWDy2QDl09pb6YqrbkBUHe4/qcpGYutW0A6H2OuN7q+j5x9fejn6/dP849nsrlkdrv4N/cFXec+vhy3d8v/b6zbOKU4uhnq/udPl993ZznzbgW2qInNXptSZQ5CxVb3RwBqZ/LcwBin3cnD2Wr22FQHlPfG5ZFT+YPjO3vljTOk5etgbLc0hI66FZ3/q2OlrKktX6GSdNwWrwIEEIIIU1P80wW5IsAIYQQ4gfmCJzeNFYXQUmqugY2RodDjd+x6lBlR2OffkOsDWHZEDeLELMMMQJuJTJtEbMq+2lJR4b1tS3MCmU+LkLluqublCZkmBwABo6vl18O6e5zwpanQ9PF98M3lkzkPhtuaF6Gsa9V+5ThXz02bQuzKu3J0LXeTl5TGe4GXNnKkqVSZfHVkoacU6uyoUbOm5Y0ZGdAjdUpUD/P0mqo523ilPp7OJnvLHlNtaRjoWUziSUpkMajxbwIEEIIIclBaYAQQghpuVAaSA8neyBFDh5Mel/HjtREP4cOf+wsaxWJJL3/lkDdwUPO71nt4s/bcTHHDc1v7dHa6Oc2al153VqrZYfq6qKfI2IfAHDo2PHoZ0/sAwA+9sR2ap+1aj/Oumo/nlg3ZB1DzZszTnEOgHse8hwAdx5jtgtwD8vrqK+hvG76GNa86XWdsam5OSjPSz2L+lr53a7OuIcag+Pq+CExp9ZcyDnU6+p50utK9D113Dh/57o10rxFxLWxzl9jnaPcz1HE74cn5yJy6DAAf536yKnJuO6D7777Lnr27NnUwyCEEHIasXfvXvTo0aNR9h3tPrjrTYTz8hLfz8GD6HDOBew+2BDdu3fH3r174XkeevXqhb1792bUhGUCkUgEPXv25NycAs5NfDg38eHcnJrTYV48z8PBgwfRvXv3NByNOQJpISsrCz169IiGIcPhcMbegE0N5yY+nJv4cG7iw7k5NZk+L/n58V0pqSQUCiGUhM6fzLaNCZsOEUIIIS2YjIsIEEIIIRkJXQPpJTs7G3PnzkV2dnZTDyXj4NzEh3MTH85NfDg3p4bzommeOQIZ5xoghBBCMomTroHq3W8l7RrIL+pH1wAhhBByepKkNJChEQG+CBBCCCF+aKY5AnQNEEIIIS0YRgQIIYQQXzTPZMGMjQjMnz8fZ599Ntq1a4ehQ4filVdeaeohpZXS0lIUFxcjLy8PXbt2xbhx41BRUeGsc+TIEZSUlKBz587Izc3FhAkTUFVVFWePzZd77rkHoVAIM2fOjP6tJc/NP/7xD1x77bXo3LkzcnJycMEFF+DVV1+NLvc8D7fffju6deuGnJwcjBo1Cjt27GjCEaeH48eP47bbbkNRURFycnJwzjnn4M4773Rq1LeUuVm3bh2+9KUvoXv37giFQnj22Wed5X7m4YMPPsCUKVMQDofRoUMHTJs2DYcOxe+x0Sw4KQ0k85OBZOSLwFNPPYXZs2dj7ty5eO211zBgwACMGTMG7733XlMPLW2sXbsWJSUl2LhxI1atWoXa2lqMHj0ahw8fjq4za9YsvPDCC3j66aexdu1a7Nu3D+PHj2/CUaefzZs34xe/+AUuvPBC5+8tdW4+/PBDjBgxAm3atMGLL76Ibdu24T//8z/RsWPH6Dr33XcfHnroITz66KPYtGkT2rdvjzFjxuDIkSNNOPLG595778XChQvxyCOP4K233sK9996L++67Dw8//HB0nZYyN4cPH8aAAQMwf/78Uy73Mw9TpkzB//3f/2HVqlVYsWIF1q1bhxtuuCFdp0BSiZeBfOYzn/FKSkqivx8/ftzr3r27V1pa2oSjalree+89D4C3du1az/M878CBA16bNm28p59+OrrOW2+95QHwNmzY0FTDTCsHDx70zj33XG/VqlXev/3bv3k333yz53kte25+9KMfeZdccknc5XV1dV5hYaH3s5/9LPq3AwcOeNnZ2d4TTzyRjiE2GWPHjvW+8Y1vOH8bP368N2XKFM/zWu7cAPCWL18e/d3PPGzbts0D4G3evDm6zosvvuiFQiHvH//4R9rGni6qq6s9AF71nh2ed6Ay4Z/qPTtO7Ke6uqlPySHjIgJHjx5FeXk5Ro0aFf1bVlYWRo0ahQ0bNjThyJqW6upqAECnTp0AAOXl5aitrXXmqW/fvujVq1eLmaeSkhKMHTvWmQOgZc/N888/j4suughf/epX0bVrVwwaNAi//OUvo8t3796NyspKZ27y8/MxdOjQZj83w4cPx+rVq/H2228DAN544w2sX78eV155JYCWPTcSP/OwYcMGdOjQARdddFF0nVGjRiErKwubNm1K+5jTRygFP5lHxiUL/vOf/8Tx48dRUFDg/L2goADbt29volE1LXV1dZg5cyZGjBiB/v37AwAqKyvRtm1bdOjQwVm3oKAAlZWVTTDK9PLkk0/itddew+bNm2OWteS5eeedd7Bw4ULMnj0bt956KzZv3ozvfve7aNu2LaZOnRo9/1M9X819bm655RZEIhH07dsXrVq1wvHjx3H33XdjypQpANCi50biZx4qKyvRtWtXZ3nr1q3RqVOn5j1XzdQ+mHEvAiSWkpISbN26FevXr2/qoWQEe/fuxc0334xVq1ahXbt2TT2cjKKurg4XXXQRfvrTnwIABg0ahK1bt+LRRx/F1KlTm3h0Tctvf/tbLF26FGVlZTj//POxZcsWzJw5E927d2/xc0NaNhknDZx55plo1apVTIZ3VVUVCgsLm2hUTcf06dOxYsUK/OlPf0KPHj2ify8sLMTRo0dx4MABZ/2WME/l5eV47733MHjwYLRu3RqtW7fG2rVr8dBDD6F169YoKChosXPTrVs3nHfeec7f+vXrhz179gBA9Pxb4vP1gx/8ALfccgsmTZqECy64ANdddx1mzZqF0tJSAC17biR+5qGwsDAmefvYsWP44IMPmvdc0TWQHtq2bYshQ4Zg9erV0b/V1dVh9erVGDZsWBOOLL14nofp06dj+fLlePnll1FUVOQsHzJkCNq0aePMU0VFBfbs2dPs5+nyyy/Hm2++iS1btkR/LrroIkyZMiX6uaXOzYgRI2Jspm+//TbOOussAEBRUREKCwuduYlEIti0aVOzn5uPPvoIWVnuV16rVq1QV1cHoGXPjcTPPAwbNgwHDhxAeXl5dJ2XX34ZdXV1GDp0aNrHnD6aZ45ARroGnnzySS87O9tbsmSJt23bNu+GG27wOnTo4FVWVjb10NLGTTfd5OXn53tr1qzx9u/fH/356KOPouvceOONXq9evbyXX37Ze/XVV71hw4Z5w4YNa8JRNx3SNeB5LXduXnnlFa9169be3Xff7e3YscNbunSpd8YZZ3iPP/54dJ177rnH69Chg/fcc895f/3rX72rrrrKKyoq8j7++OMmHHnjM3XqVO///b//561YscLbvXu3t2zZMu/MM8/0fvjDH0bXaSlzc/DgQe/111/3Xn/9dQ+A91//9V/e66+/7v3973/3PM/fPFxxxRXeoEGDvE2bNnnr16/3zj33XG/y5MlNdUqNStQ18O5uz4v8M+Gf6nd3Z6RrICNfBDzP8x5++GGvV69eXtu2bb3PfOYz3saNG5t6SGkFwCl/Fi9eHF3n448/9r7zne94HTt29M444wzvK1/5ird///6mG3QTol8EWvLcvPDCC17//v297Oxsr2/fvt5jjz3mLK+rq/Nuu+02r6CgwMvOzvYuv/xyr6KioolGmz4ikYh38803e7169fLatWvnfepTn/L+4z/+w6upqYmu01Lm5k9/+tMpv1+mTp3qeZ6/efjXv/7lTZ482cvNzfXC4bD39a9/3Tt48GATnE3jE30R+Mduzzv4r4R/qv+RmS8CbENMCCGEGETbEO/7W1LtgyORCPK7n51xbYgzLkeAEEIIISdoqBx0VVUVrr/+enTv3h1nnHEGrrjiisBlsfkiQAghhPgi/cmCVjloz/Mwbtw4vPPOO3juuefw+uuv46yzzsKoUaOccvQNwToChBBCiA8ihw4lZQGMJNCU6corr4xWv9Ts2LEDGzduxNatW3H++ecDABYuXIjCwkI88cQT+OY3v+nrGHwRIIQQQgzatm2LwsJC9Pz0+Unvq7CwMKaJVXZ2NrKzswPvq6amBgCcwmpZWVnIzs7G+vXr+SJACCGEpIJ27dph9+7dOHr0aNL7uu+++2LKN8+dOxfz5s0LvK+TPVTmzJmDX/ziF2jfvj0eeOABvPvuu9i/f7/v/dA1QAghhKSJmpqa6L/kT+I3IhAKhbB8+XKMGzcu+rfy8nJMmzYNb7zxBlq1ahVt/uR5Hl588UVfY2JEgBBCCEkTicoA8RgyZAi2bNmC6upqHD16FF26dMHQoUOdzpANQdcAIYQQcpqTn5+PLl26YMeOHXj11Vdx1VVX+d6WEQFCCCEkQzl06BB27twZ/X337t3YsmULOnXqhF69euHpp59Gly5d0KtXL7z55pu4+eabMW7cOIwePdr3MZgjQAghhGQoa9aswec+97mYv0+dOhVLlizBQw89hJ/97GeoqqpCt27d8LWvfQ233XYb2rZt6/sYfBEghBBCWjDMESCEEEJaMHwRIIQQQlowfBEghBBCWjB8ESCEEEJaMHwRIIQQQlowfBEghBBCWjB8ESCEEEJaMHwRIIQQQlowfBEghBBCWjB8ESCEEEJaMHwRIIQQQlow/z9t+pbvIzxXCAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "os.chdir(static_write_dir)\n",
    "test = read_pfb(\"pf_indicator.pfb\")\n",
    "print(test.shape)\n",
    "\n",
    "cmap1 = \"Reds\"\n",
    "sub_flip = np.zeros((1, nj, ni))\n",
    "sub_flip[0, :, :] = test[0, :, :]\n",
    "sub_flip = np.flip(sub_flip, 1)\n",
    "\n",
    "plt.imshow(sub_flip[0, :, :], cmap=cmap1)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99b69c80",
   "metadata": {},
   "source": [
    "### 8. Set up a baseline run from a reference yaml\n",
    "This function will return the correct template yaml file to do your run based on the grid, if you're doing spin-up and if you're using a solid file with the necessary keys changed to run your subset with selected climate forcing at baseline for your specified start and end dates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ea1f0e8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New runname: conus1_upper_verde provided, a new yaml file will be created\n",
      "Climate forcing directory has been changed to /home/ga6/subsettools_tutorial/inputs/conus1_upper_verde_conus1_2005WY/forcing  in runscript.\n",
      "ComputationalGrid.NY set to 90 and NX to 112\n",
      "GeomInput.domaininput.InputType detected as SolidFile, no additional keys to change for subset\n",
      "Updated runscript written to /home/ga6/subsettools_tutorial/outputs/conus1_upper_verde_conus1_2005WY\n"
     ]
    }
   ],
   "source": [
    "runscript_path = st.edit_runscript_for_subset(\n",
    "    ij_bounds,\n",
    "    runscript_path=reference_run,\n",
    "    runname=runname,\n",
    "    forcing_dir=forcing_dir,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42573426",
   "metadata": {},
   "source": [
    "### 9. Copy over your static files to your run directory\n",
    "You may only need to do this once, or you may want to copy subset static files to different run directories for different runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "71ae6dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "st.copy_files(read_dir=static_write_dir, write_dir=pf_out_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57ca532c",
   "metadata": {},
   "source": [
    "### 10. Change the file names in your runscript if desired\n",
    "If you have changed the name of a static input file either from those used in the reference yamls provided, or have changed the name of an individual file for an ensemble or other experiment, you can change it with this function by providing the target runscript (yaml or pfidb) and the new file name(s) as an arguments. Only those arguments with a specified file name will be updated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e0288e27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial pressure filename changed to conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb\n",
      "Updated runscript written to /home/ga6/subsettools_tutorial/outputs/conus1_upper_verde_conus1_2005WY\n"
     ]
    }
   ],
   "source": [
    "init_press = os.path.basename(press_init_filepath)\n",
    "runscript_path = st.change_filename_values(\n",
    "    runscript_path=runscript_path,\n",
    "    init_press=init_press,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab74c48",
   "metadata": {},
   "source": [
    "### 11. Change processor topology if desired and then distribute your inputs and forcings to match the new topology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "24c155b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distributing your climate forcing\n",
      "Distributed mask.pfb with NZ 1\n",
      "Distributed slope_x.pfb with NZ 1\n",
      "Distributed slope_y.pfb with NZ 1\n",
      "Distributed pf_indicator.pfb with NZ 5\n",
      "Distributed pme.pfb with NZ 5\n",
      "Distributed ss_pressure_head.pfb with NZ 5\n",
      "Distributed conus1_baseline_mod_2005.10.01:00.00.00_UTC0_press.pfb with NZ 5\n"
     ]
    }
   ],
   "source": [
    "runscript_path = st.dist_run(\n",
    "    topo_p=P,\n",
    "    topo_q=Q,\n",
    "    runscript_path=runscript_path,\n",
    "    dist_clim_forcing=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6081f81",
   "metadata": {},
   "source": [
    "### 12. Do a baseline run.\n",
    "Load in the yaml run file you've created which is in the same folder as your static forcings and points to your climate forcing data. This assumes you do not want to make any changes from the parent model you used (Ex. conus1 baseline) and will run your subset at baseline conditions. Outputs should be almost identical to the parent model at your subset location for the same time period if you make no additional changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "abbc70fa",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/ga6/subsettools_tutorial/outputs/conus1_upper_verde_conus1_2005WY\n",
      "Loaded run with runname: conus1_upper_verde\n",
      "\n",
      "# ==============================================================================\n",
      "# ParFlow directory\n",
      "#  - /home/SHARED/software/parflow/3.10.0\n",
      "# ParFlow version\n",
      "#  - 3.10.0\n",
      "# Working directory\n",
      "#  - /home/ga6/subsettools_tutorial/outputs/conus1_upper_verde_conus1_2005WY\n",
      "# ParFlow database\n",
      "#  - conus1_upper_verde.pfidb\n",
      "# ==============================================================================\n",
      "\n",
      "\n",
      "# ==============================================================================\n",
      "# ParFlow ran successfully 💦 💦 💦 \n",
      "# ==============================================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "set_working_directory(f\"{pf_out_dir}\")\n",
    "print(pf_out_dir)\n",
    "\n",
    "# load the specified run script\n",
    "run = Run.from_definition(runscript_path)\n",
    "print(f\"Loaded run with runname: {run.get_name()}\")\n",
    "\n",
    "# The following line is setting the run just for 10 hours for testing purposes\n",
    "run.TimingInfo.StopTime = 10\n",
    "\n",
    "run.run(working_directory=pf_out_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9354ce3-0336-47d7-b9d7-5f10f81ae6b6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}