{
 "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%2Fconus2_subsetting_transient.ipynb)\n",
    "\n",
    "This notebook walks through an example of subsetting a HUC8 from the CONUS2 domain. This example will subset everything needed to do a transient run with ParFlow-CLM. \n",
    "This is includes all hydrogeologic datasets as well as climate forcing data from CW3E. All of the data is written to a folder for the specified days to run. This example uses the template runscript conus2_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 Middle James-Buffalo watershed in Virginia for this example\n",
    "\n",
    "REMEMBER: CONUS 1 and 2 have different domains, a HUC that will run in CONUS2 may not be in the CONUS1 domain\n",
    "\n",
    "- HUC ID: 02080203\n",
    "- Size: 5240 km^2 (ni = 91, nj = 89)   "
   ]
  },
  {
   "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 = \"conus2_mjb\"\n",
    "\n",
    "# provide a way to create a subset from the conus domain (huc, lat/lon bbox currently supported)\n",
    "hucs = [\"02080203\"]\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 = \"conus2\"\n",
    "var_ds = \"conus2_domain\"\n",
    "forcing_ds = \"CW3E\"\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: (3701, 1544, 3792, 1633)\n",
      "nj: 89\n",
      "ni: 91\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": 10,
   "id": "9c3aec38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote mask.pfb\n",
      "Wrote solidfile and mask_vtk with total z of 2000 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 from the default argument var_list without the user specifying specific files. Pressure is the steady state pressure. If a user would like the override this, they may pass in their own value for var_list and their specifed variables in the target dataset will be subset. \n",
    "\n",
    "2. 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. \n",
    "\n",
    "**Note: In when working with the CONUS2 domain, there are no currently no available transient runs, so we will start runs from the steady state pressure file returned by subset_static().** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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 mannings.pfb in specified directory.\n",
      "Wrote pf_flowbarrier.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)"
   ]
  },
  {
   "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": 12,
   "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=var_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": 13,
   "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 downward_longwave to folder\n",
      "Finished writing north_windspeed 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 specific_humidity to folder\n",
      "Finished writing atmospheric_pressure to folder\n",
      "Finished writing air_temp 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": 14,
   "id": "94ce138b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(89, 91)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'pf_indicator.pfb')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAG0CAYAAAAIF+8NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABA2klEQVR4nO3deVxUZf8//tcZkAHBGYKUJUHJNAx3LSGtXCg+ZqbJp7KscHnUXeGCWCZ3uRai3fdtZqGWH0P9GbdpZWV+0zspsQVcKC3tziU3MgezglHMweD8/jAnZxiYOTNnZs4583rej/O4mzNnuTjM8Pa63tciiKIogoiIiBRD5+8CEBERkS0GZyIiIoVhcCYiIlIYBmciIiKFYXAmIiJSGAZnIiIihWFwJiIiUhgGZyIiIoVhcCYiIlIYBmdSnP/85z/o168frrrqKgiCgBEjRkg6v3379mjfvr1XygYAx44dgyAIGDNmjM3+AQMGQBAEr91XKxYvXoyUlBS0bNkSgiBg0aJFTT5TokAV7O8CEF3p2LFjGD58OCIjIzF27FgYDAYkJyf7u1iKNGDAAJSWlkJNM/CuXbsWkydPRs+ePTF58mTo9Xqkpqb6u1hEisPgTIqydetWXLhwAf/617/w4IMPunWNkpISmUvlmtWrV+P8+fN+ubdafPjhh9b/j4+Pt+4/duyYn0pEpEwMzqQoP/30EwDY/OGWqkOHDnIVR5LExES/3FdN5Pj9EgUC5pzJq7Zt2wZBEDB79mx8/vnnGDBgAFq1aoXIyEhkZmbi8OHDAP7K486aNQsAMHDgQAiCAEEQsG3bNkn3dJRznj17tvVaxcXF6NGjB8LCwhAXF4fJkyfj999/b3Sd+vp6LFiwANdddx1CQ0Nx3XXXoaCgAA0NDQ7v21zO+f3338cdd9yB6OhohIaGon379nj44Yexb98+6zEHDx7EtGnT0KtXL+txnTp1wvTp03Hu3Dmb6wmCgNLSUut/X97sc7YbN27EwIEDYTQaERYWhu7du2PhwoX4448/bI67Muf73//+F/fccw+io6MhCEKTtdorn+mKFSvQtWtXhIaG4pprrsGUKVNw9uxZ67ErV66EIAj49NNPG5XZ3v79+zF06FBERkYiIiICd9xxByoqKhyWgUirWHMmnygvL0dBQQH+53/+BxMnTsT+/fuxYcMGfPbZZygvL0dUVBRmzZqFbdu2obS0FFlZWdYAK2fnrldffRWbN2/G8OHDMWjQIGzevBmLFy/GmTNn8Oabb9oc+9hjj+GNN95AUlISsrOzceHCBSxcuBBffvmlpHtOnToVCxcuRFRUFEaMGIE2bdqgsrISW7duRe/evdGlSxcAwLvvvosVK1Zg4MCBGDBgABoaGlBeXo4FCxagtLQU27dvR4sWLQAAs2bNwsqVK3H8+HHrP2gAoEePHtb/XrhwIaZOnYqoqCg8+OCDCA8PxwcffICpU6fis88+w7vvvtsoOB4+fBipqano2rUrxowZg19++QUhISHN/nwLFy5ESUkJ7r//fgwdOhRbt27FokWLUF5ebi1zjx49mizzlY4cOYJ+/fqhV69eeOKJJ3D8+HGsX78et956Kz755BP07dtX0rMnUi2RyIs+/fRTEYAIQFy2bJnNe8uWLRMBiHfddZd136xZs0QA4qeffur2Pdu1aye2a9fOZt/l6xqNRvH777+37j9//rzYqVMnUafTiSdPnmxU7u7du4vnzp2z7v/xxx/Fq6++WgQgZmVl2dzjtttuE+2/Uhs3bhQBiF27dhXPnDlj897FixdFk8lkc22LxdLo55kzZ44IQFyzZo3T+112+PBhMTg4WGzTpo144sQJ6/4LFy6I/fv3FwGIq1evtu4/evSo9fc0c+ZMh9e0d/mZhoSEiHv37rXub2hoEB988EERgPjPf/7TpTJfef/p06fbvLd582brMyQKFGzWJp/o1KkTHn30UZt9jz76KDp27IhNmzbh559/9kk5Jk+ejOuvv976OiwsDA888AAaGhpsmk5Xr14NAJg5cybCw8Ot+6+55hpMnjzZ5fstWbIEAPDyyy8jOjra5r3g4GDExMTYXNtRLXXChAkALnWWc1VxcTH++OMPTJ06FQkJCdb9er0eCxYsAHCpqdlebGwsnn32WZfvAwCPPPIIunXrZn0tCALmzZuHoKAgh/doTmRkZKP7Z2RkYPDgwfj222/ZvE0Bg8GZfKJfv37Q6Ww/bjqdDv369YMoiti7d69PytG7d+9G+9q2bQsAqK6utu67XJ5bbrml0fGO9jVl586d0Ov1uO2225weK4oi3njjDdx6662IiopCUFAQBEGwBvXLnalc8fXXXwO4lAe3l5aWhtDQUOzZs6fRe927d3fajG3P0fNo164dEhISsH//ftTV1bl8rZ49eyIiIqLJe1z+uYi0jjln8okra4iO9tfU1PikHAaDodG+4OBLX4P6+nrrvpqaGuh0Olx99dWNjm/qZ3GkpqYG11xzTaN/mDgyadIkvPrqq0hISMDdd9+NuLg46PV6AMCcOXNgsVhcvq/ZbG6yrIIgICYmBidPnmz0npSfzdk5MTExOHbsGM6ePduo1cCdawG++5wQ+RuDM/lEVVVVs/uNRqMvi+OU0WhEQ0MDzpw5g9atW9u819TP4khkZCRMJhMaGhqaDdCnT59GYWEhunXrhrKyMrRs2dL6nslkwpw5cySV//I/QqqqqtCuXTub90RRRFVVlcN/qLgzw1lzv1tBENCqVStZrgUo73NC5C1s1iaf+OKLLxoNQWpoaMCXX34JQRDQvXt3P5XMscvl+eyzzxq952hfU2666SZYLBbrsKemHDlyBKIoIj093SYwN3e/oKAgALY1/st69uwJAA6Hoe3YsQMXLlyw6dntCUflO378OCorK5GSkiKpmfzrr79uNGzsyntc/rmItI7BmXzi4MGDWL58uc2+5cuX4+DBgxg6dGij2qm/PfzwwwCAuXPnora21rr/5MmTePnll12+TnZ2NoBLHdF+/fVXm/f++OMPa43wcu32yy+/tPlHzI8//oi8vDyH146KigIAVFZWNnrvwQcfRHBwMBYuXGiTq66rq8MzzzwDAC7PY11TU4Pvv/8ep06dcvj+6tWr8c0331hfi6KIv//976ivr5c8V3Z1dTXy8/Nt9m3ZsgUlJSXo0qWLwz4DRFrEZm3yiYyMDEyaNAn/7//9P6SkpGD//v3YuHEjrr76aknBzlcGDhyIsWPHoqioCF27dsU999wDi8WCt956C6mpqdZpKJ2588478dRTT+Gf//wnOnbsiHvuuQdt2rTByZMnUVJSgqeeego5OTmIi4tDZmYm3nnnHfTp0weDBw9GVVUVPvzwQwwePBg//PBDo2sPGjQIb7/9NjIzMzFkyBCEhoaie/fuGDZsGDp06IAFCxZg6tSp6NatG+677z6Eh4dj48aNOHDgAIYPH46HHnrIpZ9hw4YNGDt2LLKyshz2vs7IyEBaWhpGjRqF1q1bo6SkBLt370ZqaiomTpzo0j0uu+WWW7B06VLs2LEDqampOHbsGNavX4+wsDD83//9n6RrEakZa87kE6mpqSgpKUFNTQ0WL16Mbdu2YcSIESgrK8O1117r7+I5tHz5chQUFEAQBLz66qv46KOPkJubi0WLFkm6zj/+8Q+888476N69O95++20sXLgQ27dvx6BBg3D77bdbj1u5ciWmTp2K3377Da+88grKy8uRm5uL4uJih9d99NFHMW3aNJw5cwYLFizAjBkz8M4771jfz83Nxfvvv48uXbpgzZo1eOWVVxASEoJ//etfePvtt2VbQSs3Nxcvv/wyduzYgUWLFuHkyZOYPHky/vOf/0ju+X3ttdfiiy++QFhYGAoLC/HBBx9gwIAB+OyzzzgBCQUUQRRVtKQNqc62bdswcOBAzJo1C7Nnz/Z3cUhGs2fPxpw5c/Dpp586HLJFRO5jzZmIiEhhGJyJiIgUhh3CSPFWrlzp0nq/I0aMkG14EBGRPzHnTIo3YMAAp+OEAaCoqEjy0B0iIjnV19dj9uzZWLNmDUwmE+Lj4zFmzBg899xzkjphsuZMiid1PWciIn9ZsGABli5dilWrViElJQW7d+/G2LFjYTQaMWnSJJevw5ozERGRTO666y7ExMRgxYoV1n2ZmZkICwvDmjVrXL6O4mrODQ0N+Omnn9CqVSvZxmESEZFyiKKIs2fPIj4+3qVFYdx14cIFSauiNUcUxUYxSa/XWxenuezmm2/G66+/joMHD6JTp07Yu3cvPv/8cyxcuFDyDRWlsrLSuug6N27cuHHT7lZZWem1WPL777+LLSHIVtaIiIhG+2bNmtXovvX19eIzzzwjCoIgBgcHi4IgiPPmzZNcfsXVnC+vYHPsg5UwhLd0cjRJNXXwWMnn/KukyONrNHc9fwm6IbXZ9+u/K/f6PeRgfmCY1++hRoZ/b/T5PXPiOtu8XnTqv07PkeNzpjbm2vNof/cYSSuWSVVXV4fzEPEwwhECz1ph6yDi/zt3DpWVlTarudnXmgFg3bp1ePPNN1FcXIyUlBTs2bMHOTk5iI+PR1ZWlsv3VFxwvtxsYAhvCUMEg7Pc3PmQ2v8ePP2gK+X3GuRgycQr1ctQTmf3kEVwkPfvoUKOlsT0NvvvhitlkONzpla+SF2GQPD4b9ZlBoPB6e/06aefxvTp0zFq1CgAQNeuXXH8+HEUFBSoOziT8mT3HSXp+MIda71UEnnV7/vc30VoxDx8kL+LoFqG9z/x+T0fD2/r83uSNDoI0Hn4jwCd6Pqx58+fb5RHDwoKarRkrjMMzkREpFk6eD4VppTzhw0bhvz8fCQmJiIlJQVff/01Fi5ciHHjxkm6J4MzERFpliAAOg9btQXgUhcwF7zyyiuYMWMGnnzySZw+fRrx8fH429/+hpkzZ0q6J4MzERGRTFq1aoVFixZJXlrWHoMzSaaWnLK31T4r7V/CQOO8KHPM8rF/lv7IQZPy+LpZWy4MzkREpFk6QYYOYYDLzdpy4ZKRRERECsOaMxERaRabtUmRpI5RdoQ5ZsfC8+favHYlB80cs/f4IsfsbFzzstofnV5DiePrtUwnQ29tfwRnNmsTEREpDGvORESkWWzWJiIiUhhBEDyew9sfixezWZuIiEhhWHMOMOzcJR93JiEhbXGlAxj5F5u1iYiIFEatvbUZnImISLMEeB5cmXMmIiIiddac7XN99pNBAEBQl/42rwN14D9zzKRVSph0hJRPtrm1fUyVwZmIiMgVau0QxmZtIiIihWHNmYiINIu9tb3InfGkgZpjJt9xZ+ELUjep45r5d8j/2KxNREREslBFzZmIiMgdOgjQeThSmc3aREREMmLOWWa1z+cjKDjI4Xv2uT77Mc0Acz1E5DnOnU3+otjgTERE5Cm1dghjcCYiIs1Sa7O2pHvW19djxowZSEpKQlhYGDp06IDnn38eoihajxFFETNnzkRcXBzCwsKQnp6OQ4cOyV5wIiIiZy4tfCF4tPlj4QtJNecFCxZg6dKlWLVqFVJSUrB7926MHTsWRqMRkyZNAgC8+OKLWLx4MVatWoWkpCTMmDEDGRkZ+O677xAaGuqVH4L5ZVICR3O8c+yzfHwxl7an+LeI5CIpOH/55ZcYPnw4hg4dCgBo3749/v3vf2Pnzp0ALtWaFy1ahOeeew7Dhw8HAKxevRoxMTF47733MGrUKJmLT0RE1LSAaNa++eabUVJSgoMHDwIA9u7di88//xxDhgwBABw9ehQmkwnp6enWc4xGI/r27YuysjKH17RYLDCbzTYbERGRHHQybb4mqeY8ffp0mM1mJCcnIygoCPX19cjPz8fo0aMBACaTCQAQExNjc15MTIz1PXsFBQWYM2eOO2UnIiLSJEnBed26dXjzzTdRXFyMlJQU7NmzBzk5OYiPj0dWVpZbBcjLy0Nubq71tdlsRkJCAsJnPIvwiJZuXZNIKTj/tnzMwwfZvFZDDpr8T63N2pKC89NPP43p06dbc8ddu3bF8ePHUVBQgKysLMTGxgIAqqqqEBcXZz2vqqoKPXr0cHhNvV4PvV7vZvGJiIiaJs/0nb7vry3pHwTnz5+HTmd7SlBQEBoaGgAASUlJiI2NRUlJifV9s9mMHTt2IC0tTYbiEhERaZ+kmvOwYcOQn5+PxMREpKSk4Ouvv8bChQsxbtw4AIAgCMjJycELL7yAjh07WodSxcfHY8SIEd4oPxERUZMColn7lVdewYwZM/Dkk0/i9OnTiI+Px9/+9jfMnPlXHm3atGmora3FY489hurqavTv3x+bN2/22hhnIiKipgh/bp5ew9cE8crpvRTAbDbDaDTi15J1MPzZIcx+YQsO9Ce1YAcw96mhwxf/FrnHfO48ogbfh5qaGhgMBu/c489Y8mpEFMIEz+q+v4sNmHDuV6+W1x7n1iYiIs0KiGZtIiIiNVFrb20GZyIi0izWnGUWdEMqgv5s22dex7/sc/7ewN8xEdFfFBuciYiIPHVpyUjPr+FrDM5ERKRZah1K5Y+mdCIiImqGYmvO9d+Vo54LX/iFL3LMgYILX7hGqWOa2RdC/XSCAJ3A3tpERESKwWZtIiIikgVrzkREpFlqrTkzOAc4peSXOX86KQE/d9qj1uDMZm0iIiKFYc2ZiIg0SxAECB721hbYW5uIiEg+am3WZnAOMErJMQcKjmt2TAnjmplfDgw6eJ6/9Uf+lzlnIiIihWHNmYiINEsQLm0eXUOeokjC4ExERJol/Pk/T6/hawzORG5iPlndHPW/YB6alILBmYiINIu9tYmIiBRGrcGZvbWJiIgUhjVnIiLSLB0AnYdVX50oS1EkYXDWOLVMOqKGjjjsAKYtavjMkefU2lubzdpEREQKw5ozERFpmj86dHmKwZmIiDRLlhnC/BDdGZzJ59SS62OOWT5KWOiCApOvh1K1b98ex48fb7T/ySefRGFhocvXYXAmIiKSya5du1BfX299vW/fPtx+++249957JV2HwZmIiDRLBwE6D+vOUs5v3bq1zev58+ejQ4cOuO222yTeU4L27dtDEIRGW3Z2NgDgwoULyM7ORnR0NCIiIpCZmYmqqipJBSIiIpKLINMGAGaz2WazWCzN3ruurg5r1qzBuHHjIEhMXEuqOTurrk+ZMgWbNm3C+vXrYTQaMWHCBIwcORJffPGFpEKRtqglx0zaxs8heSohIcHm9axZszB79uwmj3/vvfdQXV2NMWPGSL6XpODcXHW9pqYGK1asQHFxMQYNGgQAKCoqQufOnVFeXo7U1FTJhSMiIvKEnL21KysrYTAYrPv1en2z561YsQJDhgxBfHy85Hu6nXO+XF3Pzc2FIAioqKjAxYsXkZ6ebj0mOTkZiYmJKCsrazI4WywWm6YBs9nsbpGIiIhsyNlb22Aw2ATn5hw/fhxbt27Fu+++69Y93Z4hzL66bjKZEBISgsjISJvjYmJiYDKZmrxOQUEBjEajdbNvNiAiIlKboqIitGnTBkOHDnXrfLdrzp5U16+Ul5eH3Nxc62uz2cwA3Qy558p2JQ+nlvm5neG4ZW15PLxts+8vq/3RRyUhJfPH3NoNDQ0oKipCVlYWgoPdC7NuneWouh4bG4u6ujpUV1fb1J6rqqoQGxvb5LX0er3TdnsiIiJ36AQZVqWSeP7WrVtx4sQJjBs3zv17unOSo+p679690aJFC5SUlFj3HThwACdOnEBaWprbBSQiIlKTO+64A6IoolOnTm5fQ3LNuanqutFoxPjx45Gbm4uoqCgYDAZMnDgRaWlp7KlNRER+4evpO+UiOTg3V11/6aWXoNPpkJmZCYvFgoyMDCxZssStgtU+n4+g4CAAQHj+XLeuQY25M9bT2TlayUmT+7wxd7Z5+CCP7sFxzQQEUHC+XF13JDQ0FIWFhZIm9yYiIvIWf3QIk4PbQ6mIiIjIO7jwBRERaRbXc6aA4SzHrNRcn33fBY57dp/cOWb7/LIjzsY1Ezmig+dNxP5oYmazNhERkcKw5kxERJoVML21iYiIVEMQJK+l7OgavqbY4Bw+41mER7T0dzF8yj4HKkdezxv5X6XmlKViDtox++fijXHsruSYpSrcsVb2axL5i2KDMxERkafYrE1ERKQwag3O7K1NRESkMKw5ExGRZgkydAjzuEOZGxicFUQpHcCkdgBSawcxdgDzH/vPuqMOYo2+D3aTkLADGLnCH+s5y4HBmYiINEvQCRA8jK5c+IKIiIhYcyYiIu3iwhckmTuTjsid33VnggnmmLXFftIRe45+355OTOKNSUiIHFFrcGazNhERkcKw5kxERJrFoVREREQKo9ZmbQZnH5JjYQv7XJ/U/K83FjHwB+aPfccXnxlXvgsc10yBhMGZiIg0i83aRERECqPWZm321iYiIlIY1pwVzF/jidU6jpnk4Y0cszvjmh/nXNoOOetv4WzceqDRCQJ0HlZ9PT3fHQzORESkWWpt1mZwJiIizRIgQ4cwLnxBREREiq051z6fj6DgIADqzaF4OhbXUe7P03ww88nki++TsxyzWsY1ezqvgDs8/bvhzvlq/RvrCkF3afPoGqI8ZZFCscGZiIjIYzKMc/ZH0pnN2kRERAojOTifPHkSDz30EKKjoxEWFoauXbti9+7d1vdFUcTMmTMRFxeHsLAwpKen49ChQ7IWmoiIyBWXe2t7uvmapGbt3377Df369cPAgQPx0UcfoXXr1jh06BCuuuoq6zEvvvgiFi9ejFWrViEpKQkzZsxARkYGvvvuO4SGhrpVSPscilLzI85yPVLn0mZ++C+cS9t9/vi+SP2s249pBryfc1bKPPNK+GxLHTvtSpmV8nf6UnD1dPpOmQojgaTgvGDBAiQkJKCoqMi6Lykpyfrfoihi0aJFeO655zB8+HAAwOrVqxETE4P33nsPo0aNkqnYRERE2iWpWfuDDz5Anz59cO+996JNmzbo2bMnli9fbn3/6NGjMJlMSE9Pt+4zGo3o27cvysrKHF7TYrHAbDbbbERERHJQa7O2pOB85MgRLF26FB07dsSWLVvwxBNPYNKkSVi1ahUAwGQyAQBiYmJszouJibG+Z6+goABGo9G6JSQkuPNzEBERNXJ5+k5PN5+XW8rBDQ0N6NWrF+bNm4eePXvisccew6OPPoply5a5XYC8vDzU1NRYt8rKSrevRUREpAWScs5xcXG44YYbbPZ17twZ77zzDgAgNjYWAFBVVYW4uDjrMVVVVejRo4fDa+r1euj1+mbvq5SOBVJJ7RRDf1FCJxlyn/0kJPbfBUcdwLzNnQ5g7JR5iTvfR6V05FXr3NqSas79+vXDgQMHbPYdPHgQ7dq1A3Cpc1hsbCxKSkqs75vNZuzYsQNpaWkyFJeIiMh1wp+TkHi6+ZqkmvOUKVNw8803Y968ebjvvvuwc+dOvP7663j99dcBXHoIOTk5eOGFF9CxY0frUKr4+HiMGDHCG+UnIiJqklprzpKC84033ogNGzYgLy8Pc+fORVJSEhYtWoTRo0dbj5k2bRpqa2vx2GOPobq6Gv3798fmzZvdHuNMREQUaARRFP0wpXfTzGYzjEYjKgd0h0HhC1/Y51TkzjEz3/UX5qBd5+n3xReTc9jnnH2xyIU/cs783Dpm/qMeCdv2oqamBgaDwTv3+DOW/LfrdWgVFOTRtc7W16Pzt4e9Wl57XPiCiIg0S9AJEHQezhAmKnwoFREREXkfa85ERKRZAdEhzJfCZzyL8IiW/i6GFXNHvsXn7To15piJfEWOGb4UP0MYEREReZ9ia85ERESeYrM2ERGRwsgxw5fiZwgjW5w7m9TAFzllqXwxrtkdnFuAlILBmYiINEuADM3aspREGgZnIiLSLDZrExERKY0MHcL8UXVmcG6CEsbZOsoVajUnpoTnrQaujGlWQo7ZH3NnO+OL7w4/xyQXBmciItIsNmsTEREpjKC7tHl6DV/jDGFEREQKw5rzn9SSK3KWT9RqTjpQKTXHbB4+yOb1tK0HfV4GJVDL341AxmZtIiIipdEJlzZPr+FjbNYmIiJSGNaciYhIu1S68gVrzkREpFmXc86eblKcPHkSDz30EKKjoxEWFoauXbti9+7dkq4RsDVnqR05lLrIhT86gLETjO/YP2tXOojJzb7zlyuUMOkIEQCf55x/++039OvXDwMHDsRHH32E1q1b49ChQ7jqqqsk3TJggzMREZHcFixYgISEBBQVFVn3JSUlSb4Om7WJiEi7LuecPd0AmM1mm81isTS63QcffIA+ffrg3nvvRZs2bdCzZ08sX75ccrEZnImISLMEnSDLBgAJCQkwGo3WraCgoNH9jhw5gqVLl6Jjx47YsmULnnjiCUyaNAmrVq2SVG7FNmtPHTwWIU0sBcJ8lrq5MmHFi+mdfFAS9XGU75e7P4R9jtml69stdKFV7G8R2CorK2EwGKyv9Xp9o2MaGhrQp08fzJs3DwDQs2dP7Nu3D8uWLUNWVpbL92LNmYiItEvGZm2DwWCzOQrOcXFxuOGGG2z2de7cGSdOnJBUbMXWnImIiDwlCH81S3tyDVf169cPBw4csNl38OBBtGvXTtI9WXMmIiKSyZQpU1BeXo558+bh8OHDKC4uxuuvv47s7GxJ11FFzVmOHLNWckVKWNjCfqytVp6tGnhjvL2zccyPu5BP1ko/EH6WNcjHM4TdeOON2LBhA/Ly8jB37lwkJSVh0aJFGD16tKRbqiI4ExERuUUHGSYhkXb4XXfdhbvuusuXtyQiIiJvkxScZ8+e3Wi+0eTkZOv7Fy5cQHZ2NqKjoxEREYHMzExUVVXJXmgiIiJX+GNubTlIbtZOSUnB1q1b/7pA8F+XmDJlCjZt2oT169fDaDRiwoQJGDlyJL744guPCpndd5TNa63kt9RCah7OlXHMUq/Bcc++Y5/XXubgGCX0fSByiUrXc5YcnIODgxEbG9tof01NDVasWIHi4mIMGnSpg0lRURE6d+6M8vJypKamel5aIiIiKQJlychDhw4hPj4e1157LUaPHm0dWF1RUYGLFy8iPT3demxycjISExNRVlbW5PUsFkuj+UqJiIgCmaTg3LdvX6xcuRKbN2/G0qVLcfToUdxyyy04e/YsTCYTQkJCEBkZaXNOTEwMTCZTk9csKCiwmas0ISHBrR+EiIjInqCTZ/M1Sc3aQ4YMsf53t27d0LdvX7Rr1w7r1q1DWFiYWwXIy8tDbm6u9bXZbG4UoJ3lmLU6NlEpeT2p45rt88Ny5KADFcc1e49W/26QnUBp1r5SZGQkOnXqhMOHDyM2NhZ1dXWorq62OaaqqsphjvoyvV7faL5SIiKiQOZRcD537hx++OEHxMXFoXfv3mjRogVKSkqs7x84cAAnTpxAWlqaxwUlIiKSSs4lI31JUrP2U089hWHDhqFdu3b46aefMGvWLAQFBeGBBx6A0WjE+PHjkZubi6ioKBgMBkycOBFpaWnsqU1ERP6h0mZtScH5xx9/xAMPPIBffvkFrVu3Rv/+/VFeXo7WrVsDAF566SXodDpkZmbCYrEgIyMDS5Ys8UrBvc0buT61cJaLc5aD1kqO2Rdjrf3xOXN6T7ucsxbyy0RqIyk4r13b/Jc0NDQUhYWFKCws9KhQREREsgiUSUiIiIjUQo7pN/0xfScXviAiIlIY1pwVLKhL/0b75B777M5YT3+MD5Uj/+tpLtyd8+3L6Y8cs/24Zmc/x7LaH21eK2W8vTMct0wOsVmbiIhIaWTorQ0GZyIiItkw50xERESyYM2ZiIi0izln/7GfFMMV7DxyiaNnJ/XZ+GLSEakdwJQyEYp9OZb54J7OFrawxw5gpGVs1iYiIiJZaKLmTERE5BCbtYmIiBQmEBa+IN/yRe7PnbydUvK5amCfz/WU1Hwy4CDvrdIcM1EgYXAmIiLNkmM9ZsWv50xERKQqKm3WZm9tIiIihVFszflfJUUwRLSU7XrOcqv+WJDAGV8sfOEKJeSYlVAGd9jniH3xOXP2rNSaY+a4ZnKLDjL01palJJIoNjgTERF5Sq2TkDA4ExGRdql0nDNzzkRERAoTsDVnJeaYXWGfh1ZC/tB+3mu15oc9JXX+b1c8Ht7W6TH245bt5+9WwmeEyG9U2ls7YIMzEREFAJUGZzZrExERKQxrzkREpGEy1JzBZm3ZBMqYSEdjoa/kzlzMngqUHLQrOWapfRuc/b4czdWt1ZxyoHyHyct0ukubp9fwMTZrExERKYxma85ERERq7RDG4ExERNrF4Kxsah3XrASBkkO2J3XcsjufMfscM9deJiIggIIzEREFINaciYiIFCYQe2vPnz8fgiAgJyfHuu/ChQvIzs5GdHQ0IiIikJmZiaqqKk/LSUREJN3lmrOnm4+5HZx37dqF1157Dd26dbPZP2XKFGzcuBHr169HaWkpfvrpJ4wcOdLjghIREQUKt5q1z507h9GjR2P58uV44YUXrPtramqwYsUKFBcXY9CgSx1dioqK0LlzZ5SXlyM1NVWeUtvhZAVNs++kJMekJM46SrnTYUyrnc6kPm92APtLeP5cm9f8npNbVJpzdqvmnJ2djaFDhyI9Pd1mf0VFBS5evGizPzk5GYmJiSgrK3N4LYvFArPZbLMRERHJQqXN2pJrzmvXrsVXX32FXbt2NXrPZDIhJCQEkZGRNvtjYmJgMpkcXq+goABz5syRWgwiIiLNklRzrqysxOTJk/Hmm28iNDRUlgLk5eWhpqbGulVWVspyXSIiImtvbU83H5NUc66oqMDp06fRq1cv6776+nps374dr776KrZs2YK6ujpUV1fb1J6rqqoQGxvr8Jp6vR56vV5SoV3JPQXqpCP+WOjCnrP8sSuTe/gjB+2NSUecHWOfUw7kHDORV6g05ywpOA8ePBjffvutzb6xY8ciOTkZzzzzDBISEtCiRQuUlJQgMzMTAHDgwAGcOHECaWlp8pWaiIhIwyQF51atWqFLly42+8LDwxEdHW3dP378eOTm5iIqKgoGgwETJ05EWlqa13pqExERNUmADDVnWUoiiewzhL300kvQ6XTIzMyExWJBRkYGlixZIvdtiIiInAuEZm1Htm3bZvM6NDQUhYWFKCws9PTSTQqU8Y+O8seBkktX4jhnqfljIiJ3cW5tIiLSLEGng+Bhb2tPz3cHgzMREWmYHJOIqLBZm4iISLECNefsC1rNKbtDCeOYpZI6ftjROUrMQRNR05rto2E2A3GJviuMCqkiOBMREbmFNWciIiKFkWP6TT90CPP9HYmIiKhZqqg5OxvXrJWxv2rMJwcyjmumQOXK39zm/p6Z/6iXszjNY7M2ERGRwqg0OLNZm4iISGEYnImISLsu15w93Vw0e/ZsCIJgsyUnJ0suNpu1/Yg5Zv9aVvujpOOZY/YezmWgbPY5ZlX97fJDb+2UlBRs3brV+jo4WHqoZXAmIiJygdlstnmt1+uh1+sbHRccHIzY2FiP7sVmbSIi0i4Zm7UTEhJgNBqtW0FBgcNbHjp0CPHx8bj22msxevRonDhxQnKxWXMmIiLtkrG3dmVlJQwGg3W3o1pz3759sXLlSlx//fU4deoU5syZg1tuuQX79u1Dq1atXL4lgzMREWmXjDlng8FgE5wdGTJkiPW/u3Xrhr59+6Jdu3ZYt24dxo8f7/ItGZxlpMROEvYLRrizCIUS7umPhTDYAcx72AFM+7jwxSWRkZHo1KkTDh8+LOk85pyJiEi7BMiQc3b/9ufOncMPP/yAuLg4SecxOBMRkXb5eJzzU089hdLSUhw7dgxffvkl7rnnHgQFBeGBBx6QVGw2axMREcnkxx9/xAMPPIBffvkFrVu3Rv/+/VFeXo7WrVtLuo4qgjPzU02Tmnv1Rj5Y6jVdKbOnOWapE4wQkevcWWzoyj45Wl74Yu3atZ7d60+qCM5ERERuEWTorS1wPWciIqKAx5ozERFpl0qXjFRscK59Ph9BwUH+Lobi+WP8rzP+KIMcOWaOayZyzH4OB3dyzjbn+HKcs0qDM5u1iYiIFEaxNWciIiKPCTrPO3T5oUMYgzMREWmXTri0eXoNH1NlcHYn30GOOcoP+2L+bU95YxxzUJf+Nq8DNQdtP69AeP5c2a8ZyJz9/VLiHP32PC2jb8c5q7PmzJwzERGRwqiy5kxEROSSQOitvXTpUnTr1s26pmVaWho++ugj6/sXLlxAdnY2oqOjERERgczMTFRVVcleaCIiIpdcXs/Z083HJNWc27Zti/nz56Njx44QRRGrVq3C8OHD8fXXXyMlJQVTpkzBpk2bsH79ehiNRkyYMAEjR47EF1984a3y+5UackP2vJFP9sZYa3+MlQ7UHLMzjvLFzvLQzDFf4vF4YKjz7wx5TlJwHjZsmM3r/Px8LF26FOXl5Wjbti1WrFiB4uJiDBp06cNUVFSEzp07o7y8HKmpqfKVmoiIyBWB0Kx9pfr6eqxduxa1tbVIS0tDRUUFLl68iPT0dOsxycnJSExMRFlZWZPXsVgsMJvNNhsREZEsLvfW9nTzMcl3/PbbbxEREQG9Xo/HH38cGzZswA033ACTyYSQkBBERkbaHB8TEwOTydTk9QoKCmA0Gq1bQkKC5B+CiIhISyT31r7++uuxZ88e1NTU4O2330ZWVhZKS0vdLkBeXh5yc3Otr81ms18CNPM6rlPr3Nn2mGN2nzfGQpP3eDq22pXcuaT5t306tzZkaNaWpSSSSA7OISEhuO666wAAvXv3xq5du/Dyyy/j/vvvR11dHaqrq21qz1VVVYiNjW3yenq9Hnq9XnrJiYiInJGjt7Ufemt7fMeGhgZYLBb07t0bLVq0QElJifW9AwcO4MSJE0hLS/P0NkRERAFDUs05Ly8PQ4YMQWJiIs6ePYvi4mJs27YNW7ZsgdFoxPjx45Gbm4uoqCgYDAZMnDgRaWlp7KlNRET+odLe2pKC8+nTp/HII4/g1KlTMBqN6NatG7Zs2YLbb78dAPDSSy9Bp9MhMzMTFosFGRkZWLJkieyFlmNtUSXmmO1zud4Yk+zOPZSwRrQcmGN2TI4xyRzXrBxKGVvd3DU4t7ZzkoLzihUrmn0/NDQUhYWFKCws9KhQREREshBkWJVKTeOciYiIyDu48AUREWlXIDRrExERqUogdAhTKiV27nKFs45WvuiIpdTOXt6YdIRIiTz9++VOBzBvXJMLdshLE8GZiIjIITZrExERKYxOht7anp7vzi19fkciIiJqFmvOfmQ/AYhS87+kbfaLVnBCEfl4I+/qjRyzprFDGBERkcKoNOfMZm0iIiKFYc2ZiIi0S6Udwhicyed8MYaZi1y4jjlm0jRBkKFZm8GZiIhIPirtEMacMxERkcKw5kxERNql0t7aDM6kSUFd+jfaxzw0aYH959jRZ52uoNIOYWzWJiIiUhjWnImISLvYrE1ERKQwKu2tzeBMXse1mYnkYz8uXS1zbduU02wG4hL9VxgVYHAmIiLt0ukubZ5ew8cYnImISMNkaNYGe2sTEREFPNacVc5+TWh7ztaIdiUf/Hh4W0llImXh3NnaZr9mtP363IAyxkJfWU7zH/W+uzF7axMRESkMe2sTEREpjEo7hDHnTEREpDCsORMRkXaxWZs85axzF+C8g5eza0o9n5SPHb7oSo4+D/adxHzRQcxRx7TL6s+dB7bd5/UyAPgzOHvaIYxDqYiIiAKepOBcUFCAG2+8Ea1atUKbNm0wYsQIHDhwwOaYCxcuIDs7G9HR0YiIiEBmZiaqqqpkLTQREZFLLjdre7r5mKTgXFpaiuzsbJSXl+Pjjz/GxYsXcccdd6C2ttZ6zJQpU7Bx40asX78epaWl+OmnnzBy5EjZC05EROTU5XHOnm6+LrYoiqK7J//8889o06YNSktLceutt6KmpgatW7dGcXEx/vd//xcA8P3336Nz584oKytDamqq02uazWYYjUZUDugOQ3CQu0WjJsgxSb79pAf27PPaSl34wn7Rem+Qmg9uLk/nzvWIHHH2OXPGlZx1c98v87nziBp8H2pqamAwGDwqS5P3+DOW/PrhChjCW3p2rdrziLprvFfLa8+jfw7U1NQAAKKiogAAFRUVuHjxItLT063HJCcnIzExEWVlZQ6vYbFYYDabbTYiIiJZ6AR5Nl8X290TGxoakJOTg379+qFLly4AAJPJhJCQEERGRtocGxMTA5PJ5PA6BQUFMBqN1i0hIcHdIhEREdlSabO223fMzs7Gvn37sHbtWo8KkJeXh5qaGutWWVnp0fWIiIjUzq1xzhMmTMCHH36I7du3o23bvxZFiI2NRV1dHaqrq21qz1VVVYiNjXV4Lb1eD71e704xyAVqWYidHGOOmZx9h+37gDg63v4YTz9XruSsm7tHrU8XvlDnJCSSas6iKGLChAnYsGEDPvnkEyQlJdm837t3b7Ro0QIlJSXWfQcOHMCJEyeQlpYmT4mJiIhcpdJmbUk15+zsbBQXF+P9999Hq1atrHlko9GIsLAwGI1GjB8/Hrm5uYiKioLBYMDEiRORlpbmUk9tIiIiOQmCAMHDmq+n57tDUnBeunQpAGDAgAE2+4uKijBmzBgAwEsvvQSdTofMzExYLBZkZGRgyZIlshSWiIgoEEgKzq4MiQ4NDUVhYSEKCwvdLhS5Twk5Zvv5vF3JiamRO3k7T8eX+oMrOUySj9Qcsz+48tlv9ucwm4G4RBlL1Aw5mqWV3qxNRESkKioNzlz4goiIyEvmz58PQRCQk5Mj6TzWnImISLsEGWb4crND2K5du/Daa6+hW7duks9lcFY5f+Rv7e+phByYO+znB3Y21zbHHJM3+OI77Ol31pUyKvbvgIzN2vbTSzc3T8e5c+cwevRoLF++HC+88ILkW7JZm4iIyAUJCQk2000XFBQ0eWx2djaGDh1qs9aEFKw5ExGRdsk4Q1hlZaXNqlRN1ZrXrl2Lr776Crt27XL7lgzORESkXYIgQ7P2peBsMBicLhlZWVmJyZMn4+OPP0ZoaKjbt2RwJq9T6rhmX6znrEZK/X1plbN5AJzlcr3x+3JnrLWUc8y+nFvbxyoqKnD69Gn06tXLuq++vh7bt2/Hq6++CovFgqCgIKfXYXAmIiLt8vHCF4MHD8a3335rs2/s2LFITk7GM88841JgBhiciYhIy3w8CUmrVq3QpUsXm33h4eGIjo5utL85DM5ERKRdOhnGOXt6vhsYnFWG+UDXyJFPDtRxzVqdC90fXJmXXOrz9UeO2Z3jFTvu2Q+2bdsm+RwGZyIi0i6Vzq3N4ExERNrl4w5hcuEMYURERArDmjMREWkXm7XJGSVMLuAPvuhgpNQJRcLz5/q7CJJp5XPnC+5M1iH3PZRyTcViszYRERHJgTVnIiLSLjZrExERKYxOd2nz9Bo+xuDsRVrN63ByAaJLAnXCFv4N8D4GZyIi0ixBECB42KHL0/PdweBMRETaJeN6zr7E4ExERNql0qFUDM5/kjoG2dE5nt6TXKfEhS3cGdOshsU1mF9smrPvsFa/464s6NHsOWYzEJcod7E0hcGZiIg0TIahVH6YEoTBmYiItEulzdqcIYyIiEhhNFFz5tyzvuVOft5TSswxqxU/2+QLivmccRISIiIihWGzNhEREclBcnDevn07hg0bhvj4eAiCgPfee8/mfVEUMXPmTMTFxSEsLAzp6ek4dOiQXOUlIiJy3eWFLzzdfExys3ZtbS26d++OcePGYeTIkY3ef/HFF7F48WKsWrUKSUlJmDFjBjIyMvDdd98hNDRUlkKTuvgi96SWfLI/yumN8fiBOvZZMXlUBVLss1Fps7bk4DxkyBAMGTLE4XuiKGLRokV47rnnMHz4cADA6tWrERMTg/feew+jRo3yrLREREQBQNa6+tGjR2EymZCenm7dZzQa0bdvX5SVlTk8x2KxwGw222xERETyEGTafEvW4GwymQAAMTExNvtjYmKs79krKCiA0Wi0bgkJCXIWiYiIAtnlZm1PNx/z+1CqvLw85ObmWl+bzWYkJCTA8O+NMBgMfiwZucobuSY5xjV7mztzafuCYnN/RP6g0pyzrDXn2NhYAEBVVZXN/qqqKut79vR6PQwGg81GREQUyGQNzklJSYiNjUVJSYl1n9lsxo4dO5CWlibnrYiIiFygzpyz5Gbtc+fO4fDhw9bXR48exZ49exAVFYXExETk5OTghRdeQMeOHa1DqeLj4zFixAg5y01EROScSpu1JQfn3bt3Y+DAgdbXl/PFWVlZWLlyJaZNm4ba2lo89thjqK6uRv/+/bF582aOcSYiInKR5OA8YMAAiKLY5PuCIGDu3LmYO1eZnWVImaR2APPFZB5ydPhSy+Qo5Bg718nnyu94/bnzvruxHK3Svq84+7+3NhERkfeoMzpz4QsiIiKFYc2ZiIi0K1A6hBGRsjAvKh8+S/nY9yO5sv9F7R/1viuIABmCsywlkYTN2kRERArDmjMREWmYOjuEMTgTEZF2MedM5BqlLmrh6bhmX4xp9kdO1Dx8kM/v6Qv+yi87e57Me8tNnTVn5pyJiIgUhjVnIiLSLjZrExERKQyDM5FjSs0xe4rzZpMzruTr1ZhjdvSd5vdBXgzORESkYersEMbgTEREmiUIAgQPm6U9Pd8d7K1NRESkMKw5E/3JPmcmx3rOnvJHPlKrfQS88Sy1OgZcUz8XO4QREREpjTpzzmzWJiIiUhjWnImISMNkaNZmb23ylH2uSI1jKB3xxRhKJeSYST5q+exr9Tvb7M9hNgNxib4pCHPORERESsOcMxEREcmANWciItIuNmsTEREpjDpbtRmcSRvsO3O504HM2Tne6DCmhM4/Wp10xBukTs6hhN+vHOT4fpE0DM5ERKRh6qw6MzgTEZF2qTTnzN7aRERECsOaM3nMGzlLf+R/ndFKns3+9xXUpX+z76uFL/K7WskhS6Xqz75Ka84MzkREpGHqzDl7rVm7sLAQ7du3R2hoKPr27YudO3d661ZERESa4pXg/NZbbyE3NxezZs3CV199he7duyMjIwOnT5/2xu2IiIgcE/BX07bbm++L7ZVm7YULF+LRRx/F2LFjAQDLli3Dpk2b8MYbb2D69OneuCXJyFnO0T5HKQe5c1qOctKqzptJ4GnOWK3PSasLSDgbWy3Hzyl1/Lan1zT/US/7/Zrk45zz0qVLsXTpUhw7dgwAkJKSgpkzZ2LIkCGSbil7zbmurg4VFRVIT0//6yY6HdLT01FWVtboeIvFArPZbLMRERHJQ5Bpc03btm0xf/58VFRUYPfu3Rg0aBCGDx+O/fv3Syq17MH5zJkzqK+vR0xMjM3+mJgYmEymRscXFBTAaDRat4SEBLmLRERE5BPDhg3DnXfeiY4dO6JTp07Iz89HREQEysvLJV3H77218/LykJuba31dU1ODxMREmM+e9WOp1KtRc5EbLRH15843+36Q3TWdHe+KWonNXO7cU+o9fMILLUVSn43979OnTY7epJFWOKe/Dxl+Tl//zs/+eT9RFL1+L/O5cx43a5vPnbv0/3bPWq/XQ6/XN3lefX091q9fj9raWqSlpUm7qSgzi8UiBgUFiRs2bLDZ/8gjj4h333230/MrKytFANy4cePGTeNbZWWl3CHI6vfffxdjY2NlK2tERESjfbNmzXJ472+++UYMDw8Xg4KCRKPRKG7atEly+WWvOYeEhKB3794oKSnBiBEjAAANDQ0oKSnBhAkTnJ4fHx+PyspKiKKIxMREVFZWwmAwyF3MgGI2m5GQkMBnKRM+T/nwWcpHTc9SFEWcPXsW8fHxXrtHaGgojh49irq6OlmuJ4oiBLsaeFO15uuvvx579uxBTU0N3n77bWRlZaG0tBQ33HCDy/cTRFH+doW33noLWVlZeO2113DTTTdh0aJFWLduHb7//vtGueimmM1mGI1G1NTUKP6DpnR8lvLi85QPn6V8+CyVKz09HR06dMBrr73m8jleyTnff//9+PnnnzFz5kyYTCb06NEDmzdvdjkwExERaUVDQwMsFoukc7zWIWzChAkuNWMTERFpRV5eHoYMGYLExEScPXsWxcXF2LZtG7Zs2SLpOn7vrd0UvV6PWbNmNdsTjlzDZykvPk/58FnKh89SGU6fPo1HHnkEp06dgtFoRLdu3bBlyxbcfvvtkq7jlZwzERERuY/rORMRESkMgzMREZHCMDgTEREpDIMzERGRwjA4ExERKYxig3NhYSHat2+P0NBQ9O3bFzt37vR3kRSvoKAAN954I1q1aoU2bdpgxIgROHDggM0xFy5cQHZ2NqKjoxEREYHMzExUVVX5qcTqMH/+fAiCgJycHOs+PkdpTp48iYceegjR0dEICwtD165dsXv3buv7oihi5syZiIuLQ1hYGNLT03Ho0CE/lliZ6uvrMWPGDCQlJSEsLAwdOnTA888/b7OABJ+lRrg1o7iXrV27VgwJCRHfeOMNcf/+/eKjjz4qRkZGilVVVf4umqJlZGSIRUVF4r59+8Q9e/aId955p5iYmCieO3fOeszjjz8uJiQkiCUlJeLu3bvF1NRU8eabb/ZjqZVt586dYvv27cVu3bqJkydPtu7nc3Tdr7/+KrZr104cM2aMuGPHDvHIkSPili1bxMOHD1uPmT9/vmg0GsX33ntP3Lt3r3j33XeLSUlJ4u+//+7HkitPfn6+GB0dLX744Yfi0aNHxfXr14sRERHiyy+/bD2Gz1IbFBmcb7rpJjE7O9v6ur6+XoyPjxcLCgr8WCr1OX36tAhALC0tFUVRFKurq8UWLVqI69evtx7z3//+VwQglpWV+auYinX27FmxY8eO4scffyzedttt1uDM5yjNM888I/bv37/J9xsaGsTY2FjxH//4h3VfdXW1qNfrxX//+9++KKJqDB06VBw3bpzNvpEjR4qjR48WRZHPUksU16xdV1eHiooKpKenW/fpdDqkp6ejrKzMjyVTn5qaGgBAVFQUAKCiogIXL160ebbJyclITEzks3UgOzsbQ4cOtXleAJ+jVB988AH69OmDe++9F23atEHPnj2xfPly6/tHjx6FyWSyeZ5GoxF9+/bl87Rz8803o6SkBAcPHgQA7N27F59//jmGDBkCgM9SSxQ3feeZM2dQX1/faJGMmJgYfP/9934qlfo0NDQgJycH/fr1Q5cuXQAAJpMJISEhiIyMtDk2JiYGJpPJD6VUrrVr1+Krr77Crl27Gr3H5yjNkSNHsHTpUuTm5uLvf/87du3ahUmTJiEkJARZWVnWZ+boO8/naWv69Okwm81ITk5GUFAQ6uvrkZ+fj9GjRwMAn6WGKC44kzyys7Oxb98+fP755/4uiupUVlZi8uTJ+PjjjxEaGurv4qheQ0MD+vTpg3nz5gEAevbsiX379mHZsmXIysryc+nUZd26dXjzzTdRXFyMlJQU7NmzBzk5OYiPj+ez1BjFNWtfffXVCAoKatTztaqqCrGxsX4qlbpMmDABH374IT799FO0bdvWuj82NhZ1dXWorq62OZ7P1lZFRQVOnz6NXr16ITg4GMHBwSgtLcXixYsRHByMmJgYPkcJ4uLiGi0y37lzZ5w4cQIArM+M33nnnn76aUyfPh2jRo1C165d8fDDD2PKlCkoKCgAwGepJYoLziEhIejduzdKSkqs+xoaGlBSUoK0tDQ/lkz5RFHEhAkTsGHDBnzyySdISkqyeb93795o0aKFzbM9cOAATpw4wWd7hcGDB+Pbb7/Fnj17rFufPn0wevRo63/zObquX79+jYb0HTx4EO3atQMAJCUlITY21uZ5ms1m7Nixg8/Tzvnz56HT2f7ZDgoKQkNDAwA+S03xd480R9auXSvq9Xpx5cqV4nfffSc+9thjYmRkpGgymfxdNEV74oknRKPRKG7btk08deqUdTt//rz1mMcff1xMTEwUP/nkE3H37t1iWlqamJaW5sdSq8OVvbVFkc9Rip07d4rBwcFifn6+eOjQIfHNN98UW7ZsKa5Zs8Z6zPz588XIyEjx/fffF7/55htx+PDhHP7jQFZWlnjNNddYh1K9++674tVXXy1OmzbNegyfpTYoMjiLoii+8sorYmJiohgSEiLedNNNYnl5ub+LpHgAHG5FRUXWY37//XfxySefFK+66iqxZcuW4j333COeOnXKf4VWCfvgzOcozcaNG8UuXbqIer1eTE5OFl9//XWb9xsaGsQZM2aIMTExol6vFwcPHiweOHDAT6VVLrPZLE6ePFlMTEwUQ0NDxWuvvVZ89tlnRYvFYj2Gz1IbuJ4zERGRwigu50xERBToGJyJiIgUhsGZiIhIYRiciYiIFIbBmYiISGEYnImIiBSGwZmIiEhhGJyJiIgUhsGZiIhIYRiciYiIFIbBmYiISGH+f1VKj73tMbBWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "os.chdir(static_write_dir)\n",
    "file_name = \"pf_indicator.pfb\"\n",
    "data = read_pfb(file_name)[7] \n",
    "print(data.shape)\n",
    "\n",
    "plt.imshow(data, cmap=\"Reds\", origin=\"lower\")\n",
    "plt.colorbar()\n",
    "plt.title(file_name, fontsize = 14)"
   ]
  },
  {
   "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": 15,
   "id": "ea1f0e8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver: Field BinaryOutDir is not part of the expected schema <class 'parflow.tools.database.generated.Solver'>\n",
      "New runname: conus2_mjb provided, a new yaml file will be created\n",
      "Climate forcing directory has been changed to /home/ga6/subsettools_tutorial/inputs/conus2_mjb_conus2_2005WY/forcing  in runscript.\n",
      "ComputationalGrid.NY set to 89 and NX to 91\n",
      "GeomInput.domaininput.InputType detected as SolidFile, no additional keys to change for subset\n",
      "Updated runscript written to /home/ga6/subsettools_tutorial/outputs/conus2_mjb_conus2_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": 16,
   "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": 17,
   "id": "e0288e27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver: Field BinaryOutDir is not part of the expected schema <class 'parflow.tools.database.generated.Solver'>\n",
      "Initial pressure filename changed to ss_pressure_head.pfb\n",
      "Depth to bedrock filename changed to pf_flowbarrier.pfb\n",
      "Updated runscript written to /home/ga6/subsettools_tutorial/outputs/conus2_mjb_conus2_2005WY\n"
     ]
    }
   ],
   "source": [
    "init_press_path = os.path.basename(static_paths[\"ss_pressure_head\"])\n",
    "depth_to_bedrock_path = os.path.basename(static_paths[\"pf_flowbarrier\"])\n",
    "\n",
    "runscript_path = st.change_filename_values(\n",
    "    runscript_path=runscript_path,\n",
    "    init_press=init_press_path,\n",
    "    depth_to_bedrock = depth_to_bedrock_path\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": 18,
   "id": "24c155b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver: Field BinaryOutDir is not part of the expected schema <class 'parflow.tools.database.generated.Solver'>\n",
      "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 10\n",
      "Distributed mannings.pfb with NZ 1\n",
      "Distributed pf_flowbarrier.pfb with NZ 1\n",
      "Distributed pme.pfb with NZ 10\n",
      "Distributed ss_pressure_head.pfb with NZ 10\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.\n",
    "\n",
    "**Notes about run speed**\n",
    "\n",
    "There are several things that may cause your run to take longer than expect. First, look at the file ending with your_runname.out.kinsol.log. If it is updating, it means your system is trying to solve. \n",
    "\n",
    "1. You may be trying to run too large a domain for the number of cores \n",
    "2. There is heavy rain on the day or days you selected, making the system harder to solve\n",
    "3. The initial (time 0) step can often take longer to run\n",
    "4. You have made a large change to one of your inputs. The change you made may mean the initial pressure file doesn't \"match up\" with your input. Try to run your subset model in spinup mode first, and then return to transient mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "abbc70fa",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/ga6/subsettools_tutorial/outputs/conus2_mjb_conus2_2005WY\n",
      "Solver: Field BinaryOutDir is not part of the expected schema <class 'parflow.tools.database.generated.Solver'>\n",
      "Loaded run with runname: conus2_mjb\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/conus2_mjb_conus2_2005WY\n",
      "# ParFlow database\n",
      "#  - conus2_mjb.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",
    "# The following line updates the root name of your climate forcing. You can comment it out if using NLDAS datasets. \n",
    "run.Solver.CLM.MetFileName = 'CW3E'\n",
    "\n",
    "run.run(working_directory=pf_out_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13b62bf-c43b-42df-a884-db8d30289d38",
   "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
}