{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8a97d223",
   "metadata": {
    "tags": [
     "papermill-error-cell-tag"
    ]
   },
   "source": [
    "<span style=\"color:red; font-family:Helvetica Neue, Helvetica, Arial, sans-serif; font-size:2em;\">An Exception was encountered at '<a href=\"#papermill-error-cell\">In [14]</a>'.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a15547b",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.004325,
     "end_time": "2025-08-21T16:28:04.041027",
     "exception": false,
     "start_time": "2025-08-21T16:28:04.036702",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## ROF monthly, annual, seasonal discharge at ocean outlets <a id='top'></a>\n",
    "\n",
    "Use the following datasets\n",
    "\n",
    "1. reach-D19 gauge link ascii\n",
    "2. D19 flow site geopackage\n",
    "3. D19 discharge netCDF\n",
    "4. monthly and yearly flow netCD (history file)\n",
    "\n",
    "[1. Setupt](#setup)\n",
    "\n",
    "\n",
    "[2. Loading discharge data](#load_discharge_data)\n",
    "\n",
    "- Read monthly history files from archive. \n",
    "- Reference data: monthly discharge estimates at 922 big river mouths from Dai et al. 2019 data (D19)\n",
    "\n",
    "[3. Read river, catchment, gauge information](#read_ancillary)\n",
    "\n",
    "- catchment polygon (geopackage)\n",
    "- gauge point (geopackage)\n",
    "- gauge-catchment link (csv)\n",
    "- outlet reach information (netCDF) including discharging ocean names\n",
    "\n",
    "[4. Ocean discharge line plots](#24_large_rivers)\n",
    "\n",
    "- total seasonal flow for oceans. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4ba1cfa4",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:04.047807Z",
     "iopub.status.busy": "2025-08-21T16:28:04.047449Z",
     "iopub.status.idle": "2025-08-21T16:28:06.767142Z",
     "shell.execute_reply": "2025-08-21T16:28:06.766652Z"
    },
    "papermill": {
     "duration": 2.723649,
     "end_time": "2025-08-21T16:28:06.767789",
     "exception": false,
     "start_time": "2025-08-21T16:28:04.044140",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The Python version: 3.11.4\n",
      "xarray 2025.4.0\n",
      "pandas 2.2.3\n",
      "geopandas 1.0.1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR 1: PROJ: proj_create_from_database: Open of /glade/work/hannay/miniconda3/envs/cupid-analysis/share/proj failed\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os, sys\n",
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import xarray as xr\n",
    "import cartopy.feature as cfeature\n",
    "\n",
    "from scripts.utility import load_yaml\n",
    "from scripts.utility import no_time_variable\n",
    "from scripts.utility import read_shps\n",
    "from scripts.utility import get_index_array\n",
    "\n",
    "rivers_50m = cfeature.NaturalEarthFeature(\"physical\", \"rivers_lake_centerlines\", \"50m\")\n",
    "land = cfeature.LAND\n",
    "\n",
    "print(\"\\nThe Python version: %s.%s.%s\" % sys.version_info[:3])\n",
    "print(xr.__name__, xr.__version__)\n",
    "print(pd.__name__, pd.__version__)\n",
    "print(gpd.__name__, gpd.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "413667d2",
   "metadata": {
    "papermill": {
     "duration": 0.003078,
     "end_time": "2025-08-21T16:28:06.774579",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.771501",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "-------------------------\n",
    "## 1. Setup <a id='setup'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3aaeb9bc",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:06.782049Z",
     "iopub.status.busy": "2025-08-21T16:28:06.781320Z",
     "iopub.status.idle": "2025-08-21T16:28:06.784936Z",
     "shell.execute_reply": "2025-08-21T16:28:06.784460Z"
    },
    "papermill": {
     "duration": 0.007842,
     "end_time": "2025-08-21T16:28:06.785430",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.777588",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "parameters",
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# Parameter Defaults\n",
    "# parameters are set in CUPiD's config.yml file\n",
    "# when running interactively, manually set the parameters below\n",
    "\n",
    "CESM_output_dir = \"\"\n",
    "case_name = None  # case name: e.g., \"b.e30_beta02.BLT1850.ne30_t232.104\"\n",
    "base_case_name = None  # base case name: e.g., \"b.e23_alpha17f.BLT1850.ne30_t232.092\"\n",
    "start_date = \"\"\n",
    "end_date = \"\"\n",
    "base_start_date = \"\"  # base simulation starting date: \"0001-01-01\"\n",
    "base_end_date = \"\"  # base simulation ending date: \"0100-01-01\"\n",
    "serial = True  # use dask LocalCluster\n",
    "lc_kwargs = {}\n",
    "\n",
    "analysis_name = \"\"  # Used for Figure png names\n",
    "climo_nyears = 10  # number of years to compute the climatology\n",
    "grid_name = \"f09_f09_mosart\"  # ROF grid name used in case\n",
    "base_grid_name = (\n",
    "    grid_name  # spcify ROF grid name for base_case in config.yml if different than case\n",
    ")\n",
    "figureSave = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "36376c1f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:06.792618Z",
     "iopub.status.busy": "2025-08-21T16:28:06.792149Z",
     "iopub.status.idle": "2025-08-21T16:28:06.797848Z",
     "shell.execute_reply": "2025-08-21T16:28:06.797403Z"
    },
    "papermill": {
     "duration": 0.010084,
     "end_time": "2025-08-21T16:28:06.798665",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.788581",
     "status": "completed"
    },
    "tags": [
     "injected-parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "case_name = \"b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.192\"\n",
    "base_case_name = \"b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.188\"\n",
    "CESM_output_dir = \"/glade/derecho/scratch/hannay/archive\"\n",
    "base_case_output_dir = \"/glade/derecho/scratch/gmarques/archive\"\n",
    "start_date = \"0002-01-01\"\n",
    "end_date = \"0021-12-01\"\n",
    "base_start_date = \"0002-01-01\"\n",
    "base_end_date = \"0021-12-01\"\n",
    "obs_data_dir = (\n",
    "    \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data\"\n",
    ")\n",
    "ts_dir = None\n",
    "lc_kwargs = {\"threads_per_worker\": 1}\n",
    "serial = True\n",
    "analysis_name = \"\"\n",
    "grid_name = \"f09_f09_mosart\"\n",
    "climo_nyears = 10\n",
    "figureSave = False\n",
    "subset_kwargs = {}\n",
    "product = \"/glade/work/hannay/CUPiD/examples/key_metrics/computed_notebooks//rof/global_discharge_ocean_compare_obs.ipynb\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3be8e264",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:06.805976Z",
     "iopub.status.busy": "2025-08-21T16:28:06.805704Z",
     "iopub.status.idle": "2025-08-21T16:28:06.816558Z",
     "shell.execute_reply": "2025-08-21T16:28:06.816209Z"
    },
    "papermill": {
     "duration": 0.015695,
     "end_time": "2025-08-21T16:28:06.817458",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.801763",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# ROF additional setup\n",
    "setup = load_yaml(\"./setup/setup.yaml\")\n",
    "\n",
    "domain_dir = setup[\n",
    "    \"ancillary_dir\"\n",
    "]  # ancillary directory including such as ROF domain, river network data\n",
    "geospatial_dir = setup[\"geospatial_dir\"]  # including shapefiles or geopackages\n",
    "ref_flow_dir = setup[\"ref_flow_dir\"]  # including observed or reference flow data\n",
    "case_meta = setup[\"case_meta\"]  # RO grid meta\n",
    "catch_gpkg = setup[\"catch_gpkg\"]  # catchment geopackage meta\n",
    "reach_gpkg = setup[\"reach_gpkg\"]  # reach geopackage meta\n",
    "network_nc = setup[\"river_network\"]  # river network meta\n",
    "\n",
    "if not analysis_name:\n",
    "    analysis_name = case_name\n",
    "if base_grid_name:\n",
    "    base_grid_name = grid_name\n",
    "\n",
    "case_dic = {\n",
    "    case_name: {\n",
    "        \"grid\": grid_name,\n",
    "        \"sim_period\": slice(f\"{start_date}\", f\"{end_date}\"),\n",
    "        \"climo_nyrs\": min(climo_nyears, int(end_date[:4]) - int(start_date[:4]) + 1),\n",
    "    },\n",
    "    base_case_name: {\n",
    "        \"grid\": grid_name,\n",
    "        \"sim_period\": slice(f\"{base_start_date}\", f\"{base_end_date}\"),\n",
    "        \"climo_nyrs\": min(\n",
    "            climo_nyears, int(base_end_date[:4]) - int(base_start_date[:4]) + 1\n",
    "        ),\n",
    "    },\n",
    "}\n",
    "oceans_list = [\n",
    "    \"arctic\",\n",
    "    \"atlantic\",\n",
    "    \"indian\",\n",
    "    \"mediterranean\",\n",
    "    \"pacific\",\n",
    "    \"south_china\",\n",
    "    \"global\",\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0ba6aea",
   "metadata": {
    "papermill": {
     "duration": 0.003083,
     "end_time": "2025-08-21T16:28:06.823860",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.820777",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "-----\n",
    "### dasks (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "919faf72",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:06.831130Z",
     "iopub.status.busy": "2025-08-21T16:28:06.830649Z",
     "iopub.status.idle": "2025-08-21T16:28:06.834344Z",
     "shell.execute_reply": "2025-08-21T16:28:06.833995Z"
    },
    "papermill": {
     "duration": 0.007885,
     "end_time": "2025-08-21T16:28:06.834889",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.827004",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# When running interactively, cupid_run should be set to False\n",
    "cupid_run = True\n",
    "\n",
    "client = None\n",
    "if cupid_run:\n",
    "    # Spin up cluster (if running in parallel)\n",
    "    if not serial:\n",
    "        from dask.distributed import Client, LocalCluster\n",
    "\n",
    "        cluster = LocalCluster(**lc_kwargs)\n",
    "        client = Client(cluster)\n",
    "else:\n",
    "    if not serial:\n",
    "        from dask.distributed import Client\n",
    "        from dask_jobqueue import PBSCluster\n",
    "\n",
    "        cluster = PBSCluster(\n",
    "            cores=1,\n",
    "            processes=1,\n",
    "            memory=\"50GB\",\n",
    "            queue=\"casper\",\n",
    "            walltime=\"01:00:00\",\n",
    "        )\n",
    "        cluster.scale(jobs=10)\n",
    "        client = Client(cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caa147bf",
   "metadata": {
    "papermill": {
     "duration": 0.003053,
     "end_time": "2025-08-21T16:28:06.841164",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.838111",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## 2. Loading discharge data <a id='load_discharge_data'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e770b4",
   "metadata": {
    "papermill": {
     "duration": 0.003103,
     "end_time": "2025-08-21T16:28:06.847400",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.844297",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### 2.1. Monthly/annual flow netCDFs\n",
    "- month_data (xr dataset)\n",
    "- year_data (xr dataset)\n",
    "- seas_data (xr dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "59599f11",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:06.854834Z",
     "iopub.status.busy": "2025-08-21T16:28:06.854379Z",
     "iopub.status.idle": "2025-08-21T16:28:23.403150Z",
     "shell.execute_reply": "2025-08-21T16:28:23.402584Z"
    },
    "papermill": {
     "duration": 16.553186,
     "end_time": "2025-08-21T16:28:23.403800",
     "exception": false,
     "start_time": "2025-08-21T16:28:06.850614",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished loading b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.192\n"
     ]
    },
    {
     "ename": "OSError",
     "evalue": "no files to open",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mOSError\u001b[39m                                   Traceback (most recent call last)",
      "\u001b[36mFile \u001b[39m\u001b[32m<timed exec>:28\u001b[39m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m/glade/work/hannay/miniconda3/envs/cupid-analysis/lib/python3.11/site-packages/xarray/backends/api.py:1597\u001b[39m, in \u001b[36mopen_mfdataset\u001b[39m\u001b[34m(paths, chunks, concat_dim, compat, preprocess, engine, data_vars, coords, combine, parallel, join, attrs_file, combine_attrs, **kwargs)\u001b[39m\n\u001b[32m   1594\u001b[39m paths = _find_absolute_paths(paths, engine=engine, **kwargs)\n\u001b[32m   1596\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m paths:\n\u001b[32m-> \u001b[39m\u001b[32m1597\u001b[39m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mOSError\u001b[39;00m(\u001b[33m\"\u001b[39m\u001b[33mno files to open\u001b[39m\u001b[33m\"\u001b[39m)\n\u001b[32m   1599\u001b[39m paths1d: \u001b[38;5;28mlist\u001b[39m[\u001b[38;5;28mstr\u001b[39m | ReadBuffer]\n\u001b[32m   1600\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m combine == \u001b[33m\"\u001b[39m\u001b[33mnested\u001b[39m\u001b[33m\"\u001b[39m:\n",
      "\u001b[31mOSError\u001b[39m: no files to open"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "reachID = {}\n",
    "month_data = {}\n",
    "year_data = {}\n",
    "seas_data = {}\n",
    "for case, meta in case_dic.items():\n",
    "    in_dire = os.path.join(CESM_output_dir, case, \"rof/hist\")\n",
    "    model = case_meta[meta[\"grid\"]][\"model\"]\n",
    "    domain = case_meta[meta[\"grid\"]][\"domain_nc\"]\n",
    "    var_list = case_meta[meta[\"grid\"]][\"vars_read\"]\n",
    "\n",
    "    def preprocess(ds):\n",
    "        return ds[var_list]\n",
    "\n",
    "    year_list = [\n",
    "        \"{:04d}\".format(yr)\n",
    "        for yr in np.arange(\n",
    "            int(meta[\"sim_period\"].start[0:4]), int(meta[\"sim_period\"].stop[0:4]) + 1\n",
    "        )\n",
    "    ]\n",
    "\n",
    "    nc_list = []\n",
    "    for nc_path in sorted(glob.glob(f\"{in_dire}/{case}.{model}.h*.????-*.nc\")):\n",
    "        for yr in year_list:\n",
    "            if yr in os.path.basename(nc_path):\n",
    "                nc_list.append(nc_path)\n",
    "\n",
    "    # load data\n",
    "    ds = xr.open_mfdataset(\n",
    "        nc_list,\n",
    "        data_vars=\"minimal\",\n",
    "        parallel=True,\n",
    "        preprocess=preprocess,\n",
    "    ).sel(time=meta[\"sim_period\"])\n",
    "\n",
    "    # monthly\n",
    "    month_data[case] = ds.isel(time=slice(-meta[\"climo_nyrs\"] * 12, None))\n",
    "    # annual\n",
    "    year_data[case] = (\n",
    "        ds.isel(time=slice(-meta[\"climo_nyrs\"] * 12, None))\n",
    "        .resample(time=\"YS\")\n",
    "        .mean(dim=\"time\")\n",
    "        .load()\n",
    "    )\n",
    "    # seasonal (compute here instead of reading for conisistent analysis period)\n",
    "    seas_data[case] = (\n",
    "        ds.isel(time=slice(-meta[\"climo_nyrs\"] * 12, None))\n",
    "        .groupby(\"time.month\")\n",
    "        .mean(\"time\")\n",
    "        .load()\n",
    "    )\n",
    "    vars_no_time = no_time_variable(month_data[case])\n",
    "    if vars_no_time:\n",
    "        seas_data[case][vars_no_time] = seas_data[case][vars_no_time].isel(\n",
    "            month=0, drop=True\n",
    "        )\n",
    "    mon_time = month_data[case].time.values\n",
    "    if domain == \"None\":\n",
    "        reachID[case] = month_data[case][\"reachID\"].values\n",
    "    else:\n",
    "        reachID[case] = (\n",
    "            xr.open_dataset(f\"{domain_dir}/{domain}\")[\"reachID\"]\n",
    "            .stack(seg=(\"lat\", \"lon\"))\n",
    "            .values\n",
    "        )\n",
    "    print(f\"Finished loading {case}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "590c0a25",
   "metadata": {
    "papermill": {
     "duration": 0.003591,
     "end_time": "2025-08-21T16:28:23.414687",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.411096",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### 2.2 D19 discharge data\n",
    "- ds_q_obs_mon (xr datasets)\n",
    "- ds_q_obs_yr (xr datasets)\n",
    "- dr_q_obs_seasonal (xr datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ab451e27",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:23.422453Z",
     "iopub.status.busy": "2025-08-21T16:28:23.422047Z",
     "iopub.status.idle": "2025-08-21T16:28:23.587123Z",
     "shell.execute_reply": "2025-08-21T16:28:23.586626Z"
    },
    "papermill": {
     "duration": 0.169602,
     "end_time": "2025-08-21T16:28:23.587706",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.418104",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 155 ms, sys: 3.63 ms, total: 158 ms\n",
      "Wall time: 160 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:6: DeprecationWarning: cftime_range() is deprecated, please use xarray.date_range(..., use_cftime=True) instead.\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# read monthly data\n",
    "ds_q = xr.open_dataset(\n",
    "    \"%s/D09/coastal-stns-Vol-monthly.updated-May2019.mod.nc\" % (ref_flow_dir),\n",
    "    decode_times=False,\n",
    ")\n",
    "ds_q[\"time\"] = xr.cftime_range(\n",
    "    start=\"1900-01-01\", end=\"2018-12-01\", freq=\"MS\", calendar=\"standard\"\n",
    ")\n",
    "\n",
    "# monthly- if time_period is outside observation period, use the entire obs period\n",
    "obs_available = True\n",
    "if ds_q[\"time\"].sel(time=slice(str(mon_time[0]), str(mon_time[-1]))).values.size == 0:\n",
    "    obs_available = False\n",
    "    ds_q_obs_mon = ds_q[\"FLOW\"]\n",
    "else:\n",
    "    ds_q_obs_mon = ds_q[\"FLOW\"].sel(time=slice(str(mon_time[0]), str(mon_time[-1])))\n",
    "# compute annual flow from monthly\n",
    "ds_q_obs_yr = ds_q_obs_mon.resample(time=\"YE\").mean(dim=\"time\")\n",
    "# compute annual cycle at monthly scale\n",
    "dr_q_obs_seasonal = ds_q_obs_mon.groupby(\"time.month\").mean(\"time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77fb914b",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.003394,
     "end_time": "2025-08-21T16:28:23.594685",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.591291",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## 3. Reading river, catchment, gauge infomation  <a id='read_ancillary'></a>\n",
    "\n",
    "- gauge-catchment (or grid box) link (csv)\n",
    "- gauge point (geopackage)\n",
    "- ocean polygon (geopackage)\n",
    "- catchment polygon (geopackage)\n",
    "- outlet reach information (netCDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0156464",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.003319,
     "end_time": "2025-08-21T16:28:23.601452",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.598133",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 3.1. reach-D19 gauge link csv\n",
    "- gauge_reach_lnk (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c62c8652",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:23.609284Z",
     "iopub.status.busy": "2025-08-21T16:28:23.608768Z",
     "iopub.status.idle": "2025-08-21T16:28:23.615887Z",
     "shell.execute_reply": "2025-08-21T16:28:23.615470Z"
    },
    "papermill": {
     "duration": 0.012014,
     "end_time": "2025-08-21T16:28:23.616793",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.604779",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "gauge_reach_lnk = {}\n",
    "for case, meta in case_dic.items():\n",
    "    gauge_reach_lnk[case] = pd.read_csv(\n",
    "        \"%s/D09/D09_925.%s.asc\" % (ref_flow_dir, case_meta[meta[\"grid\"]][\"network\"])\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5707ddd8",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.003318,
     "end_time": "2025-08-21T16:28:23.623766",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.620448",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 3.2 D19 flow site geopackage\n",
    "- gauge_shp (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b5cc3b41",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:23.631158Z",
     "iopub.status.busy": "2025-08-21T16:28:23.630795Z",
     "iopub.status.idle": "2025-08-21T16:28:23.693669Z",
     "shell.execute_reply": "2025-08-21T16:28:23.693208Z"
    },
    "papermill": {
     "duration": 0.067141,
     "end_time": "2025-08-21T16:28:23.694236",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.627095",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 13 ms, sys: 4.53 ms, total: 17.6 ms\n",
      "Wall time: 41.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "gauge_shp = gpd.read_file(\n",
    "    os.path.join(ref_flow_dir, \"D09\", \"geospatial\", \"D09_925.gpkg\")\n",
    ")\n",
    "gauge_shp = gauge_shp[gauge_shp[\"id\"] != 9999999]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14d20670",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.003405,
     "end_time": "2025-08-21T16:28:23.701259",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.697854",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 3.3 Ocean polygon geopackage\n",
    "- ocean_shp (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e2b5b872",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:23.709188Z",
     "iopub.status.busy": "2025-08-21T16:28:23.708685Z",
     "iopub.status.idle": "2025-08-21T16:28:34.701385Z",
     "shell.execute_reply": "2025-08-21T16:28:34.700922Z"
    },
    "papermill": {
     "duration": 10.997272,
     "end_time": "2025-08-21T16:28:34.701976",
     "exception": false,
     "start_time": "2025-08-21T16:28:23.704704",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 182 ms, sys: 10.3 s, total: 10.5 s"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Wall time: 11 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "ocean_shp = gpd.read_file(os.path.join(geospatial_dir, \"oceans.gpkg\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec053bbf",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.004014,
     "end_time": "2025-08-21T16:28:34.712083",
     "exception": false,
     "start_time": "2025-08-21T16:28:34.708069",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 3.3 Read river network information\n",
    "- gdf_cat (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "de40a200",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:34.720369Z",
     "iopub.status.busy": "2025-08-21T16:28:34.720124Z",
     "iopub.status.idle": "2025-08-21T16:28:35.383382Z",
     "shell.execute_reply": "2025-08-21T16:28:35.382920Z"
    },
    "papermill": {
     "duration": 0.668215,
     "end_time": "2025-08-21T16:28:35.384027",
     "exception": false,
     "start_time": "2025-08-21T16:28:34.715812",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished reading /glade/campaign/cesm/development/cross-wg/diagnostic_framework/rof_data/geospatial/MOSART_routing_Global_0.5x0.5_c170601_hru.gpkg\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished reading /glade/campaign/cesm/development/cross-wg/diagnostic_framework/rof_data/geospatial/MOSART_routing_Global_0.5x0.5_c170601_hru.gpkg\n",
      "CPU times: user 574 ms, sys: 7.76 ms, total: 581 ms\n",
      "Wall time: 660 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "## read catchment geopackage\n",
    "gdf_cat = {}\n",
    "for case, meta in case_dic.items():\n",
    "    cat_gpkg = os.path.join(\n",
    "        geospatial_dir, catch_gpkg[meta[\"grid\"]][\"file_name\"]\n",
    "    )  # geopackage name\n",
    "    id_name_cat = catch_gpkg[meta[\"grid\"]][\"id_name\"]  # reach ID in geopackage\n",
    "    var_list = [id_name_cat]\n",
    "    if \"lk\" in grid_name:\n",
    "        var_list.append(\"lake\")\n",
    "    gdf_cat[case] = read_shps([cat_gpkg], var_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "146ff048",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.004067,
     "end_time": "2025-08-21T16:28:35.395693",
     "exception": false,
     "start_time": "2025-08-21T16:28:35.391626",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### 3.4 Read river outlet information\n",
    "- Apppend into gdf_cat (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2a420a6f",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:35.403983Z",
     "iopub.status.busy": "2025-08-21T16:28:35.403574Z",
     "iopub.status.idle": "2025-08-21T16:28:36.241267Z",
     "shell.execute_reply": "2025-08-21T16:28:36.240764Z"
    },
    "papermill": {
     "duration": 0.842493,
     "end_time": "2025-08-21T16:28:36.241957",
     "exception": false,
     "start_time": "2025-08-21T16:28:35.399464",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 204 ms, sys: 558 ms, total: 762 ms\n",
      "Wall time: 833 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# read river outlet netcdf\n",
    "riv_ocean = {}\n",
    "for case, meta in case_dic.items():\n",
    "    riv_ocean_file = os.path.join(\n",
    "        domain_dir,\n",
    "        network_nc[meta[\"grid\"]][\"file_name\"].replace(\".aug.nc\", \".outlet.nc\"),\n",
    "    )  # network netcdf name\n",
    "    ds_rn_ocean = xr.open_dataset(riv_ocean_file).set_index(seg=\"seg_id\")\n",
    "    df_tmp = ds_rn_ocean.to_dataframe()\n",
    "    riv_ocean[case] = pd.merge(\n",
    "        gdf_cat[case],\n",
    "        df_tmp,\n",
    "        left_on=catch_gpkg[meta[\"grid\"]][\"id_name\"],\n",
    "        right_index=True,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c0af5da",
   "metadata": {
    "papermill": {
     "duration": 0.003859,
     "end_time": "2025-08-21T16:28:36.250239",
     "exception": false,
     "start_time": "2025-08-21T16:28:36.246380",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### 2.6 Merge gauge, outlet catchment dataframe\n",
    "\n",
    "- gauge_shp1 (dataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "102c6249",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:36.258611Z",
     "iopub.status.busy": "2025-08-21T16:28:36.258305Z",
     "iopub.status.idle": "2025-08-21T16:28:36.274062Z",
     "shell.execute_reply": "2025-08-21T16:28:36.273668Z"
    },
    "papermill": {
     "duration": 0.020591,
     "end_time": "2025-08-21T16:28:36.274605",
     "exception": false,
     "start_time": "2025-08-21T16:28:36.254014",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.6 ms, sys: 1.87 ms, total: 12.5 ms\n",
      "Wall time: 11.9 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Merge gauge_reach lnk (dataframe) into gauge shapefile\n",
    "gauge_shp1 = {}\n",
    "for case, meta in case_dic.items():\n",
    "    df = gauge_reach_lnk[case]\n",
    "\n",
    "    # df = df.loc[(df['flag'] == 0)]\n",
    "    df1 = df.drop(columns=[\"riv_name\"])\n",
    "    df2 = pd.merge(gauge_shp, df1, how=\"inner\", left_on=\"id\", right_on=\"gauge_id\")\n",
    "    gauge_shp1[case] = pd.merge(\n",
    "        df2,\n",
    "        riv_ocean[case],\n",
    "        how=\"inner\",\n",
    "        left_on=\"route_id\",\n",
    "        right_on=catch_gpkg[meta[\"grid\"]][\"id_name\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e1dba0b",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.003766,
     "end_time": "2025-08-21T16:28:36.282348",
     "exception": false,
     "start_time": "2025-08-21T16:28:36.278582",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "------\n",
    "## 3. Plot annual cycle for global oceans <a id='24_large_rivers'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bfdab53",
   "metadata": {
    "tags": [
     "papermill-error-cell-tag"
    ]
   },
   "source": [
    "<span id=\"papermill-error-cell\" style=\"color:red; font-family:Helvetica Neue, Helvetica, Arial, sans-serif; font-size:2em;\">Execution using papermill encountered an exception here and stopped:</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b232c13e",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2025-08-21T16:28:36.290987Z",
     "iopub.status.busy": "2025-08-21T16:28:36.290697Z",
     "iopub.status.idle": "2025-08-21T16:28:36.849666Z",
     "shell.execute_reply": "2025-08-21T16:28:36.849076Z"
    },
    "papermill": {
     "duration": 0.5641,
     "end_time": "2025-08-21T16:28:36.850241",
     "exception": true,
     "start_time": "2025-08-21T16:28:36.286141",
     "status": "failed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2 μs, sys: 2 μs, total: 4 μs\n",
      "Wall time: 8.11 μs\n"
     ]
    },
    {
     "ename": "KeyError",
     "evalue": "'b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.188'",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mKeyError\u001b[39m                                  Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[14]\u001b[39m\u001b[32m, line 36\u001b[39m\n\u001b[32m     31\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m     32\u001b[39m     id_list = gauge_shp1[case][gauge_shp1[case][\u001b[33m\"\u001b[39m\u001b[33mocean\u001b[39m\u001b[33m\"\u001b[39m] == ocean_name][\n\u001b[32m     33\u001b[39m         \u001b[33m\"\u001b[39m\u001b[33mroute_id\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m     34\u001b[39m     ].values\n\u001b[32m---> \u001b[39m\u001b[32m36\u001b[39m reach_index = get_index_array(\u001b[43mreachID\u001b[49m\u001b[43m[\u001b[49m\u001b[43mcase\u001b[49m\u001b[43m]\u001b[49m, id_list)\n\u001b[32m     37\u001b[39m seas_data_vector = seas_data[case][q_name].stack(seg=(\u001b[33m\"\u001b[39m\u001b[33mlat\u001b[39m\u001b[33m\"\u001b[39m, \u001b[33m\"\u001b[39m\u001b[33mlon\u001b[39m\u001b[33m\"\u001b[39m))\n\u001b[32m     38\u001b[39m dr_flow = seas_data_vector.isel(seg=reach_index).sum(dim=\u001b[33m\"\u001b[39m\u001b[33mseg\u001b[39m\u001b[33m\"\u001b[39m)\n",
      "\u001b[31mKeyError\u001b[39m: 'b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.188'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAu8AAAJwCAYAAAA5s4/HAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAm/FJREFUeJzs3XtclGX6P/DPAMOAJqOIcpCDYAoiZQrKwUU6KJ7L2pJSSVuz5buZItWmWXnYXdG2zLN2IFl3E6mQsl9oYikHxWNgKlYeMBBBgmQGVM737w/X2UYOMjinZ/i8X6951XPPNQ/XPSMXFzfP3CMTQgi0kxACW7duxVtvvYWGhgYsXrwYs2bNgrW1dXtPQUREREREHSTTpXm///77cf78ebz00kuIjY1Fly5dWoxzcHDQW4JERERERHSTTs27lZXV/x4okzW7XwgBmUyGxsZG/WRHREREREQaNroE79u3z1B5EBERERHRHei08t4ev/76K3r16qXPUxIREREREQCrO4fcmRACaWlpeOKJJ+Du7q6PUxIRERER0W3uqnm/cOEC3njjDXh6emLatGno0qULtm/frq/ciIiIiIjod3S65h0Aampq8Pnnn+Ojjz7CoUOHMHr0aJSUlCAvLw8BAQGGyJGIiIiIiKDjyvtf/vIXuLm5YcOGDXjqqadQXFyMr776CjKZTGsnGiIiIiIi0j+d3rBqY2OD1157DQsWLEC3bt0043K5HCdOnIC/v79BkiQiIiIiIh1X3rdu3YojR47A1dUVUVFR+H//7/+hoaHBULkREREREdHv6NS8T506Fenp6Th16hT8/Pzw4osvwtXVFU1NTcjPzzdUjkREJEGZmZmYNGkS3NzcIJPJ8MUXX9zxMRkZGQgMDISdnR18fHywefPmZjEpKSnw9/eHQqGAv78/UlNTDZA9EZF56tCF6n379sXSpUtx8eJF/Pvf/8Yf//hHTJ8+He7u7pg7d66+cyQiIgm6du0aBg8ejPXr17crvqCgAOPHj0d4eDhyc3Px+uuvY+7cuUhJSdHE5OTkICoqCtHR0Thx4gSio6MxZcoUHD582FDTICIyK3r7kKbffvsNW7duRWJiIvLy8vRxSiIishAymQypqamYPHlyqzGvvfYadu7ciTNnzmjGYmJicOLECeTk5AAAoqKioFarsWvXLk3M2LFj0aNHDyQlJRksfyIic6HzVpGtcXR0xIwZM+Do6KivU1qkpqYmXL58Gd26dYNMJjN1OkRkoYQQqKqqgpubm2R2A8vJyUFkZKTW2JgxY5CQkID6+nrI5XLk5ORg/vz5zWJWr17d6nlra2tRW1urOW5qasJvv/2Gnj17sg4TkUEZohbrrXkHgMLCQjz33HN49tln9Xlai3L58mV4eHiYOg0i6iSKiook88nXpaWlcHZ21hpzdnZGQ0MDysvL4erq2mpMaWlpq+eNj4/H0qVLDZIzEVF76LMW67V5pzu7tcVmUVERHBwcTJwNEVkqtVoNDw8PrW19peD2lfBbV3b+frylmLZW0BcuXIi4uDjNsUqlgqenJ+swERmcIWoxm3cju/UDxsHBgT80iMjgpHRZiIuLS7MV9LKyMtjY2KBnz55txty+Gv97CoUCCoWi2TjrMBEZiz5rsTQuhCQiIosXGhqK9PR0rbE9e/YgKCgIcrm8zZiwsDCj5UlEZEo6rbyvXbu2zfuLi4vvKhkiIrIc1dXVOHfunOa4oKAAeXl5cHR0hKenJxYuXIji4mJs3boVwM2dZdavX4+4uDjMnj0bOTk5SEhI0NpFZt68eRg5ciRWrlyJxx57DF9++SX27t2L7Oxso8+PiMgUdGre33vvvTvGeHp6djgZos6gpr4Rs7cew9tP3g9Xpb2p0yEymGPHjuGhhx7SHN+67nzGjBlITExESUkJCgsLNfd7e3sjLS0N8+fPx4YNG+Dm5oa1a9fij3/8oyYmLCwM27dvxxtvvIE333wT/fr1Q3JyMoKDg403MSIiE9LbPu/UPmq1GkqlEiqVitdadlIfZV3A378+g/mjBmDeqP6mTocsFGtN6/jcEJGxGKLeGPSa9/vuuw9FRUWG/BJEknK9rgGbMy5gwTg/fHqsCI1N/N2ZiIiI2s+gzfvFixdRX19vyC9BJCn/OfQLBrk54M8jfaCwscKBc+WmTomIiIgkhLvNEBnJtdoGvJ9xAfNHD4BMJkPUMA8kH+VfpoiIiKj92LwTGcnWnF8w2KM7HvDoDgD4Y6A79v1Uhorq2rYfSERERPRfbN6JjKC6tgEfZl3A/FEDNGNO9ygQMaAXdnzPLVaJiIiofdi8ExnBvw5exFDPHrjPXak1HjXMA9uPFoKbPhEREVF7sHknMrCqmnp8mHUBsS1sCxnevxdu1DXi+C9XTZAZERERSY3em/fff8rq+++/D2dnZ31/CSJJ2XLgIoK9HRHQR9nsPmsrGZ4K8sB2vnGViIiI2kFvzXtpaSleeukl3HvvvZqxqVOnomvXrvr6EkSSo7pRj4TsAsT+7lr32z0V5I5dJ0ugruG2qkRERNQ2nZr3yspKTJs2Db169dJ8bHVTUxPeeust+Pj44NChQ/j444/bfb5Nmzbh/vvvh4ODAxwcHBAaGopdu3Zp7hdCYMmSJXBzc4O9vT0efPBBnD59WusctbW1eOmll+Dk5ISuXbvi0UcfxaVLl7Rirl69iujoaCiVSiiVSkRHR6OyslIrprCwEJMmTULXrl3h5OSEuXPnoq6uTivm5MmTiIiIgL29Pfr06YNly5bxWmVq08fZBRhxb08MdG39U9Xce3RBYF9HfHXishEzIyIiIinSqXl//fXXkZmZiRkzZsDR0RHz58/HxIkTkZ2djV27duHo0aN45pln2n0+d3d3rFixAseOHcOxY8fw8MMP47HHHtM06G+//TZWrVqF9evX4+jRo3BxccHo0aNRVVWlOUdsbCxSU1Oxfft2ZGdno7q6GhMnTkRjY6MmZurUqcjLy8Pu3buxe/du5OXlITo6WnN/Y2MjJkyYgGvXriE7Oxvbt29HSkoKXn75ZU2MWq3G6NGj4ebmhqNHj2LdunV45513sGrVKl2eQupEVNfrkXjwIuY90vqq+y3PcM93IiIiag+hA09PT5Geni6EEOL8+fNCJpOJefPm6XKKO+rRo4f46KOPRFNTk3BxcRErVqzQ3FdTUyOUSqXYvHmzEEKIyspKIZfLxfbt2zUxxcXFwsrKSuzevVsIIUR+fr4AIA4dOqSJycnJEQDEjz/+KIQQIi0tTVhZWYni4mJNTFJSklAoFEKlUgkhhNi4caNQKpWipqZGExMfHy/c3NxEU1NTu+enUqkEAM15yXK9882P4sVPjrcrtra+UQxdtkecKq40cFbUWbDWtI7PDREZiyHqjU4r75cvX4a/vz8AwMfHB3Z2dnj++ef18ktEY2Mjtm/fjmvXriE0NBQFBQUoLS1FZGSkJkahUCAiIgIHDx4EABw/fhz19fVaMW5ubggICNDE5OTkQKlUIjg4WBMTEhICpVKpFRMQEAA3NzdNzJgxY1BbW4vjx49rYiIiIqBQKLRiLl++jIsXL7Y6r9raWqjVaq0bWb6r1+rwr4MXW9xhpiW2Nlb4Y6A7V9+JiIioTTo1701NTZDL5Zpja2vru35D6smTJ3HPPfdAoVAgJiYGqamp8Pf3R2lpKQA0263G2dlZc19paSlsbW3Ro0ePNmN69+7d7Ov27t1bK+b2r9OjRw/Y2tq2GXPr+FZMS+Lj4zXX2iuVSnh4eLT9hJBF+DDrAh7y6417e3dr92OmBHngi9xi1NQ33jmYiIiIOiUbXYKFEJg5c6Zm9bmmpgYxMTHNGvgdO3a0+5y+vr7Iy8tDZWUlUlJSMGPGDGRkZGjul8lkzXK4faylPH8f01K8PmLEf9+s2lY+CxcuRFxcnOZYrVazgbdwv12rw78P/YIvXhyh0+Pu7X0PfF26YdepEjw+xN1A2REREZGU6dS8z5gxQ+t4+vTpd52Ara2tZnvJoKAgHD16FGvWrMFrr70G4Oaqtqurqya+rKxMs+Lt4uKCuro6XL16VWv1vaysDGFhYZqYK1euNPu6v/76q9Z5Dh8+rHX/1atXUV9frxVz+wp7WVkZgOZ/Hfg9hUKhdakNWb73M89j9EBn9Ot1j86PjRrmie1Hiti8ExERUYt0at63bNliqDw0hBCora2Ft7c3XFxckJ6ejiFDhgAA6urqkJGRgZUrVwIAAgMDIZfLkZ6ejilTpgAASkpKcOrUKbz99tsAgNDQUKhUKhw5cgTDhw8HABw+fBgqlUrT4IeGhuIf//gHSkpKNL8o7NmzBwqFAoGBgZqY119/HXV1dbC1tdXEuLm5oW/fvgZ/Xkgayqtrse1QIXa+9IcOPX7Cfa5Y+tVpXPi1Gj4daP6JiIjIsun9E1ZvrUa3x+uvv46srCxcvHgRJ0+exKJFi7B//35MmzYNMpkMsbGxWL58OVJTU3Hq1CnMnDkTXbp0wdSpUwEASqUSs2bNwssvv4xvv/0Wubm5mD59Ou677z6MGjUKADBw4ECMHTsWs2fPxqFDh3Do0CHMnj0bEydOhK+vLwAgMjIS/v7+iI6ORm5uLr799lu88sormD17Nhwcbu7PPXXqVCgUCsycOROnTp1Camoqli9fjri4uDtexkOdx/sZ5zEmwAXeTh17L4i9rTUee8ANnx67dOdgIiIi6nx02ZrG3t5elJWVaY7HjBkjLl++rDkuLS0VVlZW7T7fn/70J+Hl5SVsbW1Fr169xCOPPCL27Nmjub+pqUksXrxYuLi4CIVCIUaOHClOnjypdY4bN26IOXPmCEdHR2Fvby8mTpwoCgsLtWIqKirEtGnTRLdu3US3bt3EtGnTxNWrV7VifvnlFzFhwgRhb28vHB0dxZw5c7S2hRRCiB9++EGEh4cLhUIhXFxcxJIlS3TaJlIIblFmya6ob4iAt3aLi+XVd3Wek5cqReDf0kVdQ6OeMqPOiLWmdXxuiMhYDFFvZEK0/yNCraystHZv6datG06cOAEfHx8AwJUrV+Dq6oqmpib9/5ZhIdRqNZRKJVQqlWZVnyzDsq/yca22ASufvP+uzzVhbRZeerg/xga46CEz6oxYa1rH54aIjMUQ9Ubvl83wEhLqjK6oa/DpsSLMefhevZzv6WEeSD5aqJdzEZnSxo0b4e3tDTs7OwQGBiIrK6vV2JkzZ0ImkzW7DRo0SBOTmJjYYkxNTY0xpkNEZHJ6b96JOqNN+89j0mA3eDh20cv5Hn2gD44U/IbLlTf0cj4iU0hOTkZsbCwWLVqE3NxchIeHY9y4cSgsbPkX0zVr1qCkpERzKyoqgqOjI5566imtOAcHB624kpIS2NnZGWNKREQmp1PzfmuFo7Vjos6oRHUDnx+/pLdVdwBQ2ssxJsAFnx/nG1dJulatWoVZs2bh+eefx8CBA7F69Wp4eHhg06ZNLcYrlUq4uLhobseOHcPVq1fx3HPPacXJZDKtOBcXXl5GRJ2HTs27EAIDBgyAo6MjHB0dUV1djSFDhmiO/fz8DJUnkdnauO88Jg9xQ5/u9no979PDPJF8tAhNTe1+WwqR2airq8Px48cRGRmpNR4ZGYmDBw+26xwJCQkYNWoUvLy8tMarq6vh5eUFd3d3TJw4Ebm5uW2ep7a2Fmq1WutGRCRVZrfPO5GUFFfeQGpuMdLjRur93MP69oDCxgrZ58oxckAvvZ+fyJDKy8vR2NjY7EPsnJ2dm33gXUtKSkqwa9cubNu2TWvcz88PiYmJuO+++6BWq7FmzRqMGDECJ06cQP/+/Vs8V3x8PJYuXdrxyRARmRGdmvdp06bBxkanhxBZtA37zuGPQ/vAVanfVXfg5qUBUcM8kHy0iM07Sdbtl1YKIdp1uWViYiK6d++OyZMna42HhIQgJCREczxixAgMHToU69atw9q1a1s818KFCxEXF6c5VqvV8PDw0GEWRETmQ6fLZtzc3PDKK6/gzJkzhsqHSDKKfruOnXmX8ZeH9Het++2eGOqO734sQ0V1rcG+BpEhODk5wdrautkqe1lZWbPV+NsJIfDxxx8jOjpa84nWrbGyssKwYcNw9uzZVmMUCgUcHBy0bkREUqVT8z5//nx89dVXCAgIQGhoKBISElBdXW2o3IjM2oZ95/BUkDucHQy3y0WvbgpEDOiF1Nxig30NIkOwtbVFYGAg0tPTtcbT09MRFhbW5mMzMjJw7tw5zJo1645fRwiBvLw8uLq63lW+RERSoVPzvnDhQvz000/Yv38//Pz8EBsbC1dXVzz33HM4cOCAoXIkMjuFFdfx/34owf9F9DP413p6uAe2Hy2CDp+nRmQW4uLi8NFHH+Hjjz/GmTNnMH/+fBQWFiImJgbAzZ8pzz77bLPHJSQkIDg4GAEBAc3uW7p0Kb755htcuHABeXl5mDVrFvLy8jTnJCKydB3a5z08PBxbtmxBaWkpVq9ejXPnziE8PBy+vr54++239Z0jkdlZ991ZPD3MA70NuOp+S3j/Xrhe24DvC68a/GsR6VNUVBRWr16NZcuW4YEHHkBmZibS0tI0u8eUlJQ02/NdpVIhJSWl1VX3yspKvPDCCxg4cCAiIyNRXFyMzMxMDB8+3ODzISIyBzKhp+W8r7/+Gs8++ywqKyvR2Nioj1NaJH4st/RdLL+GSeuy8d0rD6JXN4VRvuZ76T/jcuUN/POpwUb5eiR9rDWt43NDRMZiiHpzV5+wev36dWzZsgUjR47Eo48+ip49e+If//iHXhIjMldrvzuLqcGeRmvcAeCpIHeknSxBVU290b4mERERmZ8ONe9ZWVn405/+BBcXF8yZMwfe3t7Yt28ffv75ZyxYsEDfORKZjQu/ViP99BW8MNLHqF/XvUcXBPZ1xM4Tl436dYmIiMi86NS8L1++HAMGDMCDDz6I06dP45///CdKSkrwr3/9CyNH6v9DaojMzdpvz2J6qBd63mO8Vfdbnv7vnu9ERETUeenUvL/33nuYMGECTpw4gcOHD+PPf/4zrxekTuNcWRW+/bEML4Qbd9X9llEDnVF89QZOX1aZ5OsTERGR6en0camXL1+GXC4HcPMd/+fOnYNMJkO/fv3QvXt3Q+RHZDbWfHsOM0L7okfXtj80xlBsbazwx0B3fHq0CEsfU5okByIiIjItnVbe5XI5Ll68iAkTJsDJyQnBwcEYPnw4nJycMHHiRFy8eNFAaRKZ1s9XqpDxUxmeD/c2aR5TgjzwRd5l1NRzRyciIqLOSKeV96KiIoSEhEAul+Nvf/sbBg4cCCEEzpw5g02bNiE0NBRHjx6Fu7u7ofIlMok1e89i5ghvdO9imlX3W+7tfQ8GON+DXadK8PgQfp8RERF1Njo174sXL4avry+++eYb2Nn978NpHn/8ccyfPx9jx47F4sWLkZCQoPdEiUzlx1I1ss+VY/kT95k6FQBA1DBPbD9SxOadiIioE9Lpspndu3fjH//4h1bjfou9vT3+9re/YdeuXXpLjsgcrE4/iz+N8IbSXm7qVAAA4+9zQX6JGhd+rTZ1KkRERGRkOjXvFRUV6Nu3b6v3+/j4oKKi4m5zIjIbpy+rcKigAs/9oa+pU9HoYmuDRwe74dNjl0ydChERERmZTs27m5sbTp8+3er9p06dgqur610nRWQuVu89i+f/4A0HO/NYdb/lmeGe+Pz4JdQ3Npk6FSIiIjIinZr3xx57DK+++ip+/fXXZveVlZXhtddew+TJk/WVG5FJnbykwrGLv2FGWF9Tp9JMQB8lnB0U+O7HMlOnQkREREak8xtW09LS0K9fP0yfPh1+fn4AgPz8fGzbtg0uLi546623DJIokbGt3vszZo/0QTczW3W/5dYnro4Z5GLqVIiIiMhIdFp579GjBw4fPoxp06Zh+/btiI2NRWxsLD799FNMnToVOTk5cHR0bPf54uPjMWzYMHTr1g29e/fG5MmT8dNPP2nFCCGwZMkSuLm5wd7eHg8++GCzS3dqa2vx0ksvwcnJCV27dsWjjz6KS5e0rwe+evUqoqOjoVQqoVQqER0djcrKSq2YwsJCTJo0CV27doWTkxPmzp2Luro6rZiTJ08iIiIC9vb26NOnD5YtWwYhRLvnTNLww6VK5BZVYkZoX1On0qpHH+iDwxcqUKK6YepUiIiIyEh0at6Bmw38pk2bUFFRgdLSUpSWlqKiogKbN29Gz549tWIPHDiA2traVs+VkZGBF198EYcOHUJ6ejoaGhoQGRmJa9euaWLefvttrFq1CuvXr8fRo0fh4uKC0aNHo6qqShMTGxuL1NRUbN++HdnZ2aiursbEiRPR2Pi/D7KZOnUq8vLysHv3buzevRt5eXmIjo7W3N/Y2IgJEybg2rVryM7Oxvbt25GSkoKXX35ZE6NWqzF69Gi4ubnh6NGjWLduHd555x2sWrVK16eRzFzigYt4NtQLXRU6/XHKqJT2cowZ5ILP+MZVIiKizkMYULdu3cT58+fbHV9WViYAiIyMDCGEEE1NTcLFxUWsWLFCE1NTUyOUSqXYvHmzEEKIyspKIZfLxfbt2zUxxcXFwsrKSuzevVsIIUR+fr4AIA4dOqSJycnJEQDEjz/+KIQQIi0tTVhZWYni4mJNTFJSklAoFEKlUgkhhNi4caNQKpWipqZGExMfHy/c3NxEU1NTu+aoUqkEAM05yfxcvVYrBr65S1yuvG7qVO7o0PlyERb/rWhsbN+/P+o8WGtax+eGiIzFEPVG55V3HX8x0ClepVIBgObSm4KCApSWliIyMlITo1AoEBERgYMHDwIAjh8/jvr6eq0YNzc3BAQEaGJycnKgVCoRHBysiQkJCYFSqdSKCQgIgJubmyZmzJgxqK2txfHjxzUxERERUCgUWjGXL1/GxYsXW5xTbW0t1Gq11o3M247vixHWzwmuSntTp3JHw70dobCxwoHz5aZOhYiIiIzAoM27LoQQiIuLwx/+8AcEBAQAAEpLSwEAzs7OWrHOzs6a+0pLS2Fra4sePXq0GdO7d+9mX7N3795aMbd/nR49esDW1rbNmFvHt2JuFx8fr7nOXqlUwsPD4w7PBJmSEALbjhRiWrCnqVNpF5lMhqhhHth+tMjUqRAREZERmE3zPmfOHPzwww9ISkpqdp9MJtM6FkI0G7vd7TEtxesj5tZfF1rLZ+HChVCpVJpbURGbLHN29OJV3KhrxMgBvUydSrs9MdQd+34sw2/X6u4cTERERJJmFs37Sy+9hJ07d2Lfvn1wd3fXjLu43NwC7/ZV7bKyMs2Kt4uLC+rq6nD16tU2Y65cudLs6/76669aMbd/natXr6K+vr7NmLKym/ts374if4tCoYCDg4PWjcxX0pFCRA3zgLVV278cmpNe3RQY2b8XdnzPN64SERFZOoM27+1ZHZ8zZw527NiB7777Dt7e3lr3e3t7w8XFBenp6Zqxuro6ZGRkICwsDAAQGBgIuVyuFVNSUoJTp05pYkJDQ6FSqXDkyBFNzOHDh6FSqbRiTp06hZKSEk3Mnj17oFAoEBgYqInJzMzU2j5yz549cHNzQ9++fXV5asgMXb1Wh29OlyJqmPQubYoafvPSGV3fZ0JkaBs3boS3tzfs7OwQGBiIrKysVmP3798PmUzW7Pbjjz9qxaWkpMDf3x8KhQL+/v5ITU019DSIiMyH3t762oJ77rmnzd1m/u///k8olUqxf/9+UVJSorldv/6/XT5WrFghlEql2LFjhzh58qR45plnhKurq1Cr1ZqYmJgY4e7uLvbu3Su+//578fDDD4vBgweLhoYGTczYsWPF/fffL3JyckROTo647777xMSJEzX3NzQ0iICAAPHII4+I77//Xuzdu1e4u7uLOXPmaGIqKyuFs7OzeOaZZ8TJkyfFjh07hIODg3jnnXfa/ZxwlwPz9WHmeTH7X0dNnUaHNDQ2iZDle8WxixWmToXMhDnUmu3btwu5XC4+/PBDkZ+fL+bNmye6du0qfvnllxbj9+3bJwCIn376Setnwu9r+cGDB4W1tbVYvny5OHPmjFi+fLmwsbHR2k3sTszhuSGizsEQ9Uan5v3KlStt3l9fXy8OHz7c/i8OtHjbsmWLJqapqUksXrxYuLi4CIVCIUaOHClOnjypdZ4bN26IOXPmCEdHR2Fvby8mTpwoCgsLtWIqKirEtGnTRLdu3US3bt3EtGnTxNWrV7VifvnlFzFhwgRhb28vHB0dxZw5c7S2hRRCiB9++EGEh4cLhUIhXFxcxJIlS9q9TaQQ/KFhrpqamsTD7+wT+35s+9+4OXt3z0/ilU/zTJ0GmQlzqDXDhw8XMTExWmN+fn5iwYIFLcbfat5vr82/N2XKFDF27FitsTFjxoinn3663XmZw3NDRJ2DIeqNTIj2/53d2toaJSUlmp1bBg4ciG+++Qaenjd35rhy5Qrc3Ny0PhyJtKnVaiiVSqhUKl7/bkYOX6hA3KcnkPnXhyR1vfvvXbp6HWPey8Sh1x9BNzu5qdMhEzN1ramrq0OXLl3w2Wef4fHHH9eMz5s3D3l5ecjIyGj2mP379+Ohhx5C3759UVNTA39/f7zxxht46KGHNDGenp6YP38+5s+frxl77733sHr1avzyyy8t5lJbW6v1gYFqtRoeHh6sw0RkcIaoxTpd8357n3/p0iU0NDS0GUMkBduOFOKZ4dJ6o+rt3Ht0QWBfR3x1ouTOwUQGVl5ejsbGxja3+r2dq6srPvjgA6SkpGDHjh3w9fXFI488gszMTE1Ma1v2tnZOgFv2EpFl0ftnv9/pTapE5ua3a3VIz7+CReMHmjqVu/b0MA+8n3EeUyWyTz1ZPl22+vX19YWvr6/mODQ0FEVFRXjnnXcwcuTIDp0TuLllb1xcnOb41so7EZEUmcVWkUSmtOP7Swjv74TeDnamTuWujRrojKKrN5B/mZ/kS6bl5OQEa2vrNrf6bY+QkBCcPXtWc9zalr1tnZNb9hKRJdGpeZfJZKiqqoJarYZKpYJMJkN1dTXUarXmRiQl4r+fqDo12MvUqeiFrY0Vnhnugdlbj2FV+s8oKL9m6pSok7K1tUVgYKDWNr4AkJ6ertmitz1yc3Ph6uqqOQ4NDW12zj179uh0TiIiKdPpshkhBAYMGKB1PGTIEK1jXjZDUnLowm+ob2xC+L1Opk5Fb14e7YuR/Xvhi7xiPLY+Gz697sHjQ/pg4v2u6HmPwtTpUScSFxeH6OhoBAUFITQ0FB988AEKCwsRExMD4OblLMXFxdi6dSsAYPXq1ejbty8GDRqEuro6/Oc//0FKSgpSUlI055w3bx5GjhyJlStX4rHHHsOXX36JvXv3Ijs72yRzJCIyNp2a93379hkqDyKT2HakEE8P84SVhN+oejsrKxmCfXoi2KcnFk8ahH0/liE1txgrdv2I0H49MXlIH4we6Ax7W2tTp0oWLioqChUVFVi2bBlKSkoQEBCAtLQ0eHnd/EtXSUkJCgsLNfF1dXV45ZVXUFxcDHt7ewwaNAhff/01xo8fr4kJCwvD9u3b8cYbb+DNN99Ev379kJycjODgYKPPj4jIFHTaKpLunqm3b6P/qaiuRfjb+7D/1QfRu5v0r3e/k8rrdfj6ZAm+yC3GjyVVGBPggseH9EGIT09J77JDLWOtaR2fGyIyFkPUm7vabeb06dNae7pbW1tj0KBBd50UkTGkfH8JEQN6dYrGHQC6d7HFtGAvTAv2QmHFdXyZV4w3vziFa3UNeOyBPnh8SB8MdGUjQ0REZM50at6zsrIQFxeHo0ePAri5C8D169c1e7vLZDJ88803GDVqlP4zJdIjIQSSjhRh2WOd85dNz55d8NIj/THn4XtxsliFHd8XY/pHh9GrmwKTh/TBYw+4wVVpb+o0iYiI6DY6Ne8bN25EdHS01ti+ffvg5eUFIQTWrl2LTZs2sXkns5dzvgKNTQIj+lnOG1U7QiaT4X737rjfvTsWTRiI7HPl+CK3GOu+PYv73bvj8aF9MDbABQ78xFYiIiKzoFPzfvToUcybN09rzN3dXfPmo+joaEyYMEF/2REZyLYjhXh6uIdFvVH1bsmtrfCQb2885Nsb1bUN+OZUKb7IK8aSnafxkF9vPP5AH0T49oLcmh8PQUREZCo6Ne/FxcVa++3+61//gouLi+bY0dERFRUV+suOyADKq2vx3Y9lWDypc14y0x73KGzwx0B3/DHQHVfUNfjqxGWsSv8Zr35+AmMDXNC3Z1e4KO3g4mAHF6UdnB3sYCfn7jVERESGplPz3q1bNxQUFGhW2p944gmt+wsKCvjOfTJ7nx+/hId8e6NXN+553h7ODnZ4PtwHz4f74KfSKuw9cwW//HYdRwp+Q6m6BlfUNai4VgelvRwuDjcbedf/NvS3GvxbY927yPlZEERERHdBp+Y9ODgYW7duxYMPPtji/YmJidxrl8xaU5NA0pFC/GPyfaZORZJ8XbrB16Vbs/HahkaUqWtxRV2DUnUNSlU3bwfOlWvGrqhqARlurtY72MFZaQcXB8V/G3t7uChv/r+b0p6XMxEREbVCp+Y9Li4Oo0aNQs+ePfHqq6+id+/eAICysjKsXLkS//nPf7Bnzx6DJEqkDzkXbl7WFdavp4kzsSwKG2t4OHaBh2OXVmOEEPjtWp1mtb5EVYMrqhr8fKUKmWfLcUVVg8uqG+jeRY6nh3niqUB39HboHNt4EhERtZdOzftDDz2EdevWYf78+Vi1ahUcHBwgk8mgUqlgY2OD1atX4+GHHzZUrkR3bdvhQjwz3LI+UVUqZDIZet6jQM97FBjkpmwxRgiBQxd+w/ajhXjwnf0Yca8TnhnugYgBvflBUkREROjgJ6wWFRXh888/x9mzZwEA/fv3x5NPPgkPDw+9J2hp+Ml+pvNrVS0i/rkPmX99CE738Hp3c3f1Wh125BZj+5FCVNc2YEqQB6YM80Cf7tx/vj1Ya1rH54aIjMUQ9aZDzXtbrl+/ji5dWv/TeWfHHxqms3H/OZy+rMaGqUNNnQrpQAiB7wuvIulIEXadLMEwb0c8PcwTjwzszW0r28Ba0zo+N0RkLIaoN3r7yVdTU4N3330XPj4++jolkd40NQlsP1KEacM9TZ0K6UgmkyHQyxHvPDUYOa8/gkcGOmPtt2cRtuI7vL37R/xScc3UKRIRERmNTs17XV0dFi1ahGHDhiEsLAxffPEFAGDLli3w8fHBqlWrmn2IE5E5OHC+HNZWMoTyjaqS5mAnR3SIF76e+wckzAjC1ev1mLg2G9M+OoSvTlxGbUOjqVMkIiIyKJ3esLpkyRJs2LABo0ePxoEDB/DUU0/hT3/6E/bv34/4+HhMnToVcjk/Rp3Mz803qnpwj3ELIZPJcL97d9zv3h1vTBiI//fDZXyUXYDFO0/jiSF98PRwT9zb+x5Tp0lERKR3OjXvn376KRITE/H444/jxIkTGDJkCNRqNU6fPg0bG51ORWQ0ZVU1yPj5V/x9coCpUyED6KqwQdQwT0QN88SZEjW2HynEExsPwM/FAU8P98D4+1z56a9ERGQxdHrDqkKhwPnz5+Hu7g4AsLOzw6FDh/DAAw8YKj+LwzdKGd+GfefwY2kV1j0zxNSpkJHU1Dci7WQJth8pwk9XqjD5ATc8PdwTA107z/cca03r+NwQkbGY/A2r9fX1sLW11RzL5XIolS3v10xkDpqaBLYfLcRUvlG1U7GTW+OJoe74NCYUKf8XCrm1FaZ+eAhPbT6I3MKrpk6PiIiow3Tebeatt95CXFwc4uLiUFdXh7///e+a41u39srMzMSkSZPg5uYGmUymeQPsLUIILFmyBG5ubrC3t8eDDz6I06dPa8XU1tbipZdegpOTE7p27YpHH30Uly5d0oq5evUqoqOjoVQqoVQqER0djcrKSq2YwsJCTJo0CV27doWTkxPmzp2Luro6rZiTJ08iIiIC9vb26NOnD5YtWwY977RJepZ1rhxyKyuE+DiaOhUykXt7d8MbE/1x6PVH8OhgN/wp8She/vQEytQ1pk6NiIhIZzo17yNHjsRPP/2E3Nxc5ObmIiwsDBcuXNAc37q117Vr1zB48GCsX7++xfvffvttrFq1CuvXr8fRo0fh4uKC0aNHo6qqShMTGxuL1NRUbN++HdnZ2aiursbEiRPR2Pi/XSemTp2KvLw87N69G7t370ZeXh6io6M19zc2NmLChAm4du0asrOzsX37dqSkpODll1/WxKjVaowePRpubm44evQo1q1bh3feeQerVq3S5SkkI9t2+Bc8M9yTb1QlKGysER3aF/teeRBdFdYYtSoDm/af5w41REQkLcJMABCpqama46amJuHi4iJWrFihGaupqRFKpVJs3rxZCCFEZWWlkMvlYvv27ZqY4uJiYWVlJXbv3i2EECI/P18AEIcOHdLE5OTkCADixx9/FEIIkZaWJqysrERxcbEmJikpSSgUCqFSqYQQQmzcuFEolUpRU1OjiYmPjxdubm6iqamp3fNUqVQCgOa8ZDhXVDfEwDd3iYrqWlOnQmboTIlKPP1+joh4+zuRfrpUp+9jKTCXWrNhwwbRt29foVAoxNChQ0VmZmarsSkpKWLUqFHCyclJdOvWTYSEhGhq+S1btmwRAJrdbty40e6czOW5ISLLZ4h6o9ePJzx58iRiY2P1cq6CggKUlpYiMjJSM6ZQKBAREYGDBw8CAI4fP476+nqtGDc3NwQEBGhicnJyoFQqERwcrIkJCQmBUqnUigkICICbm5smZsyYMaitrcXx48c1MREREVAoFFoxly9fxsWLF1udR21tLdRqtdaNjOPTY0UY7e8Mx662dw6mTsfPxQHbZgfjtbF+WLzzNGZsOYpzZdWmTsuiJCcnIzY2FosWLUJubi7Cw8Mxbtw4FBYWthifmZmJ0aNHIy0tDcePH8dDDz2ESZMmNfuLroODA0pKSrRudnZ2xpgSEZHJ3XXzrlar8f7772P48OEYPHgw9u/fr4e0gNLSUgCAs7Oz1rizs7PmvtLSUtja2qJHjx5txvTu3bvZ+Xv37q0Vc/vX6dGjB2xtbduMuXV8K6Yl8fHxmmvtlUolPDw82p446UVjk0DSkSI8wzeqUhtkMhnG3eeKb1+OwDCvHnh84wH87f/lQ3Wj3tSpWYRVq1Zh1qxZeP755zFw4ECsXr0aHh4e2LRpU4vxq1evxl//+lcMGzYM/fv3x/Lly9G/f3989dVXWnEymQwuLi5aNyKizqLDzXtGRgaeffZZuLq64i9/+Qsefvhh/Pzzz8jLy9Njemh2rbIQ4o7XL98e01K8PmLEf9+s2lY+CxcuhEql0tyKiorazJ30I/Psr7CTWyHYm29UpTuzk1vjpUf645vYkfi1qhaPvLsf248UorGJb0jvqLq6Ohw/flzrL6MAEBkZqfmr5500NTWhqqoKjo7a38fV1dXw8vKCu7s7Jk6ceMf3WvEvoERkSXRq3ktKSrB8+XLce++9ePrpp+Hk5ISMjAxYWVnh2Wefxb333qu3xG6tpNy+ql1WVqZZ8XZxcUFdXR2uXr3aZsyVK1eanf/XX3/Virn961y9ehX19fVtxpSVlQFo/teB31MoFHBwcNC6keElHS7kG1VJZ27d7bH2mSHYND0Q/z70Cx7bkI2jF38zdVqSVF5ejsbGxjb/enon7777Lq5du4YpU6Zoxvz8/JCYmIidO3ciKSkJdnZ2GDFiBM6ePdvqefgXUCKyJDo1797e3jhz5gw2bNiA4uJirFq1CkFBQQZJzNvbGy4uLkhPT9eM1dXVISMjA2FhYQCAwMBAyOVyrZiSkhKcOnVKExMaGgqVSoUjR45oYg4fPgyVSqUVc+rUKZSUlGhi9uzZA4VCgcDAQE1MZmam1vaRe/bsgZubG/r27av/J4A67Iq6BtnnyvHHoe6mToUkalhfR+yc8wdMC/ZCzL+PY25SLkpUN0ydliR15K+nAJCUlIQlS5YgOTlZ69LHkJAQTJ8+HYMHD0Z4eDg+/fRTDBgwAOvWrWv1XPwLKBFZEp2ady8vL2RnZyMzMxM///zzXX/x6upq5OXlaS61KSgoQF5eHgoLCyGTyRAbG4vly5cjNTUVp06dwsyZM9GlSxdMnToVAKBUKjFr1iy8/PLL+Pbbb5Gbm4vp06fjvvvuw6hRowAAAwcOxNixYzF79mwcOnQIhw4dwuzZszFx4kT4+voCuPlnXH9/f0RHRyM3NxfffvstXnnlFcyePVuzUj516lQoFArMnDkTp06dQmpqKpYvX464uDiu7pqZ5KNFiPR3Rg++UZXugrWVDM8M98R3rzwIp3sUiHwvE+u+PYuaem4t2R5OTk6wtrZu86+nrUlOTsasWbPw6aefamp5a6ysrDBs2LA2V975F1Aisii6bk+TnZ0tnnvuOXHPPfeIoUOHilWrVgkbGxuRn5+v81Y3+/bta3HLrxkzZgghbm4XuXjxYuHi4iIUCoUYOXKkOHnypNY5bty4IebMmSMcHR2Fvb29mDhxoigsLNSKqaioENOmTRPdunUT3bp1E9OmTRNXr17Vivnll1/EhAkThL29vXB0dBRz5szR2hZSCCF++OEHER4eLhQKhXBxcRFLlizReXs5blFmWA2NTSJ0+V5x+EKFqVMhC3P2ilpM/+iQGLHiW7Hr5GWz31rSHGrN8OHDxf/93/9pjQ0cOFAsWLCg1cds27ZN2NnZaW0d3JampiYRFBQknnvuuXbnZQ7PDRF1DoaoNzIhOvYRodXV1UhKSsLHH3+Mw4cPIyIiAlOnTsXkyZPRq1cvPf1qYXnUajWUSiVUKhVXfwxg349l+EfaGaTPH8m/iJDeCSGw90wZ/v51Pvp0t8dbk/zh52Ke38fmUGuSk5MRHR2NzZs3IzQ0FB988AE+/PBDnD59Gl5eXli4cCGKi4uxdetWADcvlXn22WexZs0aPPHEE5rz2NvbQ6lUAgCWLl2KkJAQ9O/fH2q1GmvXrsW///1vHDhwAMOHD29XXubw3BBR52CIetPh3WbuuecezJ49Gzk5OTh9+jQCAwPxxhtvaO2VTmRsnxwuxFS+UZUMRCaTYbS/M/bMH4k/9HfCU5tz8NaXp1B5ve7OD+6EoqKisHr1aixbtgwPPPAAMjMzkZaWBi8vLwA336P0+z3f33//fTQ0NODFF1+Eq6ur5jZv3jxNTGVlJV544QUMHDgQkZGRKC4uRmZmZrsbdyIiqevwyntL6uvr8dVXX2mtmJA2rvgYTonqBh55NwMHFzyM7l14vTsZ3hV1DVbu/hH7f/oV80f1x9PDPSG31utn33UYa03r+NwQkbEYot7Y6JrAndzpzUVEhvLp0UsYO8iFjTsZjbODHVZNeQDfF17F3/5fPtbvO4dpwV54ergHenfjJ34SEZH+6dS8d+/evc3LEcR/twBrbORuDGRcjU0CyUcLsfaZIaZOhTqhoZ49sOP/wpBbVImtBy/ioX/uxyh/Zzwb2hdDPduum0RERLrQqXnft2+f5v+FEBg/fjw++ugj9OnTR++JEeli/09luMfOBoFePUydCnVSMpkMQz17YKhnD/xaVYvtRwrx4iffw6mbLZ4N7YtHB7vBTm5t6jSJiEjidGreIyIitI6tra0REhICHx8fvSZFpKtt/ERVMiO9uinw0iP9EfNgP6TnX8G/Dl5EfNoZTBnmgenBXvBw7GLqFImISKLM451VRHfhcuUNHLpQgSeG8BNVybzIra0w/j5XJP85FEkvhKCqpgHj1mTh+X8dQ9bZX9HUpLf9AoiIqJPQaeWdyBwlHy3C2ABXKLvITZ0KUav8XByw/PH78NpYP3x+/BLe/OIUrKxkeDbEC38MdEc3O/77JSKiO7vrlXdepkCm1NDYhOSjRZga7GHqVIjaRWkvx6w/eOO7lx/EWxP9kXm2HGHx3+HNL07h7JUqU6dHRERmTqeV99v3b6+pqUFMTAy6du2qNb5jx467z4yoHfb99CuU9nIM9eQbVUlarKxkeNC3Nx707Y1fKq7hP4d+wZObczDIzQHPhvbFqIG9YWMme8YTEZH50Kl5v/Xx1LdMnz5dr8kQ6SrpSCGmBvONqiRtXj27YtEEf8SN9sWXecVYvfdnLPvqNKaFeOHpYR7oeY/C1CkSEZGZ0Kl537Jli6HyINJZceUNHL5QgfeiHjB1KkR6YW9rjaeHeyJqmAeO/XIViQcvIvztfRgb4IKZYX1xv3t3U6dIREQmptPfZBsbG/HDDz/gxo0bze67fv06fvjhBzQ1NektOaK2JB0uxLj7XKG05xv9yLLIZDIM6+uIDVOH4ruXH4R7jy74IveyqdMiIiIzoNPK+7///W+sX78ehw8fbnafQqHAn/70J8TGxvJyGjKIxiaBvKJK7D1zBXvzr6CsqhbbZgebOi0ig3JR2iFu9ABTp0FERGZCp+Y9ISEBr7zyCqytm39KoLW1Nf76179i/fr1bN5Jb67XNSD7bDn2nrmC734sg72tNUYPdMHSxwZhWF9HyPmGPiIiIupEdGref/rpJ4SEhLR6/7Bhw3DmzJm7Too6tzJ1Db79sQx786/g4PkK+Ll2w6iBztg2OwT9e9/DN6cSERFRp6VT837t2jWo1epW76+qqsL169fvOin6n2u1Ddi4/xycHezg7GAHl//+1+keW4vZRk4IgZ+uVGFv/hWknynDz6VVGHGvEyIHOSP+j/ehdzc7U6dIREREZBZ0at779++PgwcP4v7772/x/uzsbPTv318vidFNDY0CjU1AbmElrqhrUKquwRVVDW7UN8LpHgVclHb/bewVmsbe2cHu5ng3OzjY25jlSnVdQxOOFPx28/r1M1dQU9+EUQN746WH7sWIe51gb9v80iwiIiKizk6n5n3q1Kl44403EBYW1qyBP3HiBN566y389a9/1WuCnZ2yixwLxvk1G6+ubUCpqgZltxp6dS2uqGtwqliNK1U3G/yyqlrYWMvg4mCH3ppVe8X/mvv/jjnYy9HF1trg14+rrtdj/89lSM+/goyffoVbd3uM8u+Ndc8MwWD37rCyMr9fMoiIiIjMiUwIIdobXF9fj8jISGRnZ2PUqFHw8/ODTCbDmTNnsHfvXowYMQLp6emQy7l1X2vUajWUSiVUKhUcHBwM+rWamgQqrtXhirrmf6v26lpcUdXgSlUNSlU3x9U1DWhsErC1tkIXhTW6yK1hb2uNrgob2Mv/+19ba3S1tUYXWxt0sbVGF1tr2NvaoKvtf2NvjSv+d38XWxuob9Tjux/LsPfMFRz/5SoCvXpg1EBnjBroDM+eXQw6f6LOzJi1Rmr43BCRsRii3ujUvAM3G/j33nsP27Ztw9mzZyGEwIABAzB16lTExsbC1tZWL4lZKnP8oSGEQF1jE67XNuJ6fSOu1zbgel0jrtU14EZdI67XNeJ6XcN//3vz/6/VNuLG72L+99//jV+va4TC2goRvr0w2t8ZDw7oDWUX/mJHZAzmWGvMBZ8bIjIWs2je6e50th8aQgizvOaeyNJ1tlqjCz43RGQshqg3d32Rc01NDf71r39h48aNOHv2rD5yIgvCxp2oc9u4cSO8vb1hZ2eHwMBAZGVltRmfkZGBwMBA2NnZwcfHB5s3b24Wk5KSAn9/fygUCvj7+yM1NdVQ6RMRmR2dmvdXX30V8+bN0xzX1dUhJCQEs2fPxuuvv44hQ4YgJydH70kSEZH0JCcnIzY2FosWLUJubi7Cw8Mxbtw4FBYWthhfUFCA8ePHIzw8HLm5uXj99dcxd+5cpKSkaGJycnIQFRWF6OhonDhxAtHR0ZgyZUqLn/xNRGSJdLpsJiAgAMuXL8ejjz4KANiyZQtefvll5ObmwtPTE3/6059QVlaGr7/+2mAJSx3/XEtExmAOtSY4OBhDhw7Fpk2bNGMDBw7E5MmTER8f3yz+tddew86dO7U+7C8mJgYnTpzQLAxFRUVBrVZj165dmpixY8eiR48eSEpKalde5vDcEFHnYIh6o9NWkYWFhfD399cc79mzB08++SS8vLwAAPPmzcP48eP1kpiluvW7UlsfdkVEdLdu1RhTva2prq4Ox48fx4IFC7TGIyMjcfDgwRYfk5OTg8jISK2xMWPGICEhAfX19ZDL5cjJycH8+fObxaxevbrVXGpra1FbW6s5VqlUAFiHicjwDFGLdWreraystL74oUOH8Oabb2qOu3fvjqtXr+otOUtUVVUFAPDw8DBxJkTUGVRVVUGpVBr965aXl6OxsRHOzs5a487OzigtLW3xMaWlpS3GNzQ0oLy8HK6urq3GtHZOAIiPj8fSpUubjbMOE5GxVFRU6K0W69S8+/n54auvvkJcXBxOnz6NwsJCPPTQQ5r7f/nll2ZFlbS5ubmhqKgI3bp1k+SbOdVqNTw8PFBUVGTRf27uLPMEOs9cO8s8gf/NNT8/H25ubibN5fY6d6cdqFqKv31c13MuXLgQcXFxmuPKykp4eXmhsLDQJL/YGIOl/3u39PkBnKOlUKlU8PT0hKOjo97OqVPz/uqrr+KZZ57B119/jdOnT2P8+PHw9vbW3J+Wlobhw4frLTlLZGVlBXd3d1OncdccHBws9hvt9zrLPIHOM9fOMk8A6NOnD6ysDPvJya1xcnKCtbV1sxXxsrKyVhd5XFxcWoy3sbFBz54924xpa+FIoVBAoVA0G1cqlRb/b8HS/71b+vwAztFS6LMW63SmP/7xj0hLS8P999+P+fPnIzk5Wev+Ll264C9/+YvekiMiImmytbVFYGAg0tPTtcbT09MRFhbW4mNCQ0Obxe/ZswdBQUGaT+5uLaa1cxIRWRqdVt4BYNSoURg1alSL9y1evPiuEyIiIssQFxeH6OhoBAUFITQ0FB988AEKCwsRExMD4OblLMXFxdi6dSuAmzvLrF+/HnFxcZg9ezZycnKQkJCgtYvMvHnzMHLkSKxcuRKPPfYYvvzyS+zduxfZ2dkmmSMRkbHp1Lz/8MMP7Yq7//77O5QMmT+FQoHFixe3+CdoS9JZ5gl0nrl2lnkC5jPXqKgoVFRUYNmyZSgpKUFAQADS0tI0O5SVlJRo7fnu7e2NtLQ0zJ8/Hxs2bICbmxvWrl2LP/7xj5qYsLAwbN++HW+88QbefPNN9OvXD8nJyQgODm53Xuby/BiSpc/R0ucHcI6WwhBz1GmfdysrK8hksha3u7k1LpPJ0NjYqLcEiYiIiIjoJp1W3gsKCgyVBxERERER3YFOK+9ERERERGQ6Or9hFQDOnj2LL7/8EhcvXoRMJoO3tzcmT54MHx8ffedHRERERET/pfPKe3x8PN566y00NTWhd+/eEELg119/hbW1NZYvX45XXnnFULkSEREREXVqOu3zvm/fPrzxxhtYtGgRysvLUVJSgtLSUvz6669YsGABFixYgMzMTEPlSgYWHx+PYcOGoVu3bujduzcmT56Mn376qc3H7N+/HzKZrNntxx9/NFLWuluyZEmzfF1cXNp8TEZGBgIDA2FnZwcfHx9s3rzZSNnenb59+7b4+rz44ostxkvl9czMzMSkSZPg5uYGmUyGL774Qut+IQSWLFkCNzc32Nvb48EHH8Tp06fveN6UlBT4+/tDoVDA398fqampBppB+7U11/r6erz22mu477770LVrV7i5ueHZZ5/F5cuX2zxnYmJii69zTU2NgWdjPBs3boS3tzfs7OwQGBiIrKysNuOl9j2uy/x27NiB0aNHo1evXnBwcEBoaCi++eYbI2bbMbq+hrccOHAANjY2eOCBBwyboB7oOsfa2losWrQIXl5eUCgU6NevHz7++GMjZdsxus7xk08+weDBg9GlSxe4urriueeeQ0VFhZGy1c2dfha1RC+1RuhgypQp4oUXXmj1/tmzZ4unn35al1OSGRkzZozYsmWLOHXqlMjLyxMTJkwQnp6eorq6utXH7Nu3TwAQP/30kygpKdHcGhoajJi5bhYvXiwGDRqklW9ZWVmr8RcuXBBdunQR8+bNE/n5+eLDDz8UcrlcfP7550bMumPKysq05pmeni4AiH379rUYL5XXMy0tTSxatEikpKQIACI1NVXr/hUrVohu3bqJlJQUcfLkSREVFSVcXV2FWq1u9ZwHDx4U1tbWYvny5eLMmTNi+fLlwsbGRhw6dMjAs2lbW3OtrKwUo0aNEsnJyeLHH38UOTk5Ijg4WAQGBrZ5zi1btggHBwet17ikpMTAMzGe7du3C7lcLj788EORn58v5s2bJ7p27Sp++eWXFuOl9j2u6/zmzZsnVq5cKY4cOSJ+/vlnsXDhQiGXy8X3339v5MzbT9c53lJZWSl8fHxEZGSkGDx4sHGS7aCOzPHRRx8VwcHBIj09XRQUFIjDhw+LAwcOGDFr3eg6x6ysLGFlZSXWrFkjLly4ILKyssSgQYPE5MmTjZx5+9zpZ9Ht9FVrdGre+/btK7Kyslq9PzMzU/Tt21enBMh8lZWVCQAiIyOj1Zhbzd7Vq1eNl9hdWrx4sU5F/a9//avw8/PTGvvzn/8sQkJC9JyZ4c2bN0/069dPNDU1tXi/FF/P2wtmU1OTcHFxEStWrNCM1dTUCKVSKTZv3tzqeaZMmSLGjh2rNTZmzBizWpBozw+HI0eOCABtNgBbtmwRSqVSv8mZkeHDh4uYmBitMT8/P7FgwYIW46X2Pa7r/Fri7+8vli5dqu/U9Kajc4yKihJvvPGGznXeFHSd465du4RSqRQVFRXGSE8vdJ3jP//5T+Hj46M1tnbtWuHu7m6wHPWlPfVZX7VGp8tmrly5gr59+7Z6v7e3N0pLS3Vf/iezpFKpAACOjo53jB0yZAhcXV3xyCOPYN++fYZO7a6dPXsWbm5u8Pb2xtNPP40LFy60GpuTk4PIyEitsTFjxuDYsWOor683dKp6U1dXh//85z/405/+BJlM1mas1F7P3ysoKEBpaanWa6ZQKBAREYGDBw+2+rjWXue2HmOOVCoVZDIZunfv3mZcdXU1vLy84O7ujokTJyI3N9c4CRpYXV0djh8/3uy1jIyMbPW1lNL3eEfmd7umpiZUVVW1q7abQkfnuGXLFpw/f14Sn/bekTnu3LkTQUFBePvtt9GnTx8MGDAAr7zyCm7cuGGMlHXWkTmGhYXh0qVLSEtLgxACV65cweeff44JEyYYI2WD01et0al5r6mpga2tbav3y+Vy1NXV6XJKMlNCCMTFxeEPf/gDAgICWo1zdXXFBx98gJSUFOzYsQO+vr545JFHzPq9D8HBwdi6dSu++eYbfPjhhygtLUVYWFir19SVlpbC2dlZa8zZ2RkNDQ0oLy83Rsp68cUXX6CyshIzZ85sNUaKr+ftbi0gtPSatbW40NrrLKUFiZqaGixYsABTp06Fg4NDq3F+fn5ITEzEzp07kZSUBDs7O4wYMQJnz541YraGUV5ejsbGRp1eSyl9j3dkfrd79913ce3aNUyZMsUQKd61jszx7NmzWLBgAT755BPY2HRoIz2j6sgcL1y4gOzsbJw6dQqpqalYvXo1Pv/881bfw2RqHZljWFgYPvnkE0RFRcHW1hYuLi7o3r071q1bZ4yUDU5ftUbnf+EfffQR7rnnnhbvq6qq0vV0ZKbmzJmDH374AdnZ2W3G+fr6wtfXV3McGhqKoqIivPPOOxg5cqSh0+yQcePGaf7/vvvuQ2hoKPr164d//etfiIuLa/Ext69Ui/9u0nSnFWxzkpCQgHHjxsHNza3VGCm+nq1p6TW70+vVkceYi/r6ejz99NNoamrCxo0b24wNCQlBSEiI5njEiBEYOnQo1q1bh7Vr1xo6VaPQ9bWU2vd4R/+tJiUlYcmSJfjyyy/Ru3dvQ6WnF+2dY2NjI6ZOnYqlS5diwIABxkpPL3R5HZuamiCTyfDJJ59AqVQCAFatWoUnn3wSGzZsgL29vcHz7Qhd5pifn4+5c+firbfewpgxY1BSUoJXX30VMTExSEhIMEa6BqePWqNT8+7p6YkPP/zwjjEkbS+99BJ27tyJzMxMuLu76/z4kJAQ/Oc//zFAZobRtWtX3Hfffa2uOrq4uDRbJSgrK4ONjQ169uxpjBTv2i+//IK9e/dix44dOj9Waq/nrZ2DSktL4erqqhkvKytrtuJx++Naep3beoy5qK+vx5QpU1BQUIDvvvuuzVX3llhZWWHYsGEWsfLu5OQEa2trnV5LKX2Pd2R+tyQnJ2PWrFn47LPPMGrUKEOmeVd0nWNVVRWOHTuG3NxczJkzB8DNRlcIARsbG+zZswcPP/ywUXJvr468jq6urujTp4+mcQeAgQMHQgiBS5cuoX///gbNWVcdmWN8fDxGjBiBV199FQBw//33o2vXrggPD8ff//53rZouRfqqNTpdNnPx4kUUFBTc8UbSJITAnDlzsGPHDnz33Xfw9vbu0Hlyc3Ml9Q1WW1uLM2fOtJpzaGgo0tPTtcb27NmDoKAgyOVyY6R417Zs2YLevXt36LpBqb2e3t7ecHFx0XrN6urqkJGRgbCwsFYf19rr3NZjzMGtxv3s2bPYu3dvh5pNIQTy8vIk9Tq3xtbWFoGBgc1ey/T09FZfSyl9j3dkfsDNFfeZM2di27ZtZn/9sK5zdHBwwMmTJ5GXl6e5xcTEwNfXF3l5eQgODjZW6u3WkddxxIgRuHz5MqqrqzVjP//8M6ysrDq00GZoHZnj9evXYWWl3ZpaW1sD+N8KtZTprdbo9PZWHQUEBIjCwkJDfgnSo//7v/8TSqVS7N+/X2v7uOvXr2tiFixYIKKjozXH7733nkhNTRU///yzOHXqlFiwYIEAIFJSUkwxhXZ5+eWXxf79+8WFCxfEoUOHxMSJE0W3bt3ExYsXhRDN53hra6f58+eL/Px8kZCQYNbbyN2usbFReHp6itdee63ZfVJ9PauqqkRubq7Izc0VAMSqVatEbm6uZoeVFStWCKVSKXbs2CFOnjwpnnnmmWZbRUZHR2vteHDgwAFhbW0tVqxYIc6cOSNWrFhhFltFtjXX+vp68eijjwp3d3eRl5en9X1bW1urOcftc12yZInYvXu3OH/+vMjNzRXPPfecsLGxEYcPHzbFFPXu1vZ0CQkJIj8/X8TGxoquXbtazPe4rvPbtm2bsLGxERs2bND6N1JZWWmqKdyRrnO8nRR2m9F1jlVVVcLd3V08+eST4vTp0yIjI0P0799fPP/886aawh3pOsctW7YIGxsbsXHjRnH+/HmRnZ0tgoKCxPDhw001hTbd6WeRoWqNQZv3e+65R5w/f96QX4L0CECLty1btmhiZsyYISIiIjTHK1euFP369RN2dnaiR48e4g9/+IP4+uuvjZ+8Dm7t+S2Xy4Wbm5t44oknxOnTpzX33z5HIYTYv3+/GDJkiLC1tRV9+/YVmzZtMnLWHffNN99o9m6/nVRfz1tbWt5+mzFjhhDi5naRixcvFi4uLkKhUIiRI0eKkydPap0jIiJCE3/LZ599Jnx9fYVcLhd+fn5m8UtLW3MtKCho9fv293v53z7X2NhY4enpKWxtbUWvXr1EZGSkOHjwoPEnZ0AbNmwQXl5ewtbWVgwdOlRry1tL+B7XZX4RERFtfr+YK11fw9+TQvMuhO5zPHPmjBg1apSwt7cX7u7uIi4uTmuBzRzpOse1a9cKf39/YW9vL1xdXcW0adPEpUuXjJx1+9zpZ5Ghao1MCMP9HaJbt244ceIEfHx8DPUliIiIiIg6DZ2ueSciIiIiItNh805EREREJBEW2bxnZmZi0qRJcHNzg0wmwxdffHHHx2RkZCAwMBB2dnbw8fHB5s2bm8WkpKTA398fCoUC/v7+SE1NNUD2RESWgbWYiEj/LLJ5v3btGgYPHoz169e3K76goADjx49HeHg4cnNz8frrr2Pu3LlISUnRxOTk5CAqKgrR0dE4ceIEoqOjMWXKFBw+fNhQ0yAikjTWYiIi/dP7G1aLi4vRp08fAMC2bdvw2GOPoWvXrvr8EjqRyWRITU3F5MmTW4157bXXsHPnTpw5c0YzFhMTgxMnTiAnJwcAEBUVBbVajV27dmlixo4dix49eiApKclg+RMRWQLWYiIi/dDpE1bbUlpain/84x/46KOPcOPGDQDA1KlT9XV6g8rJyUFkZKTW2JgxY5CQkID6+nrI5XLk5ORg/vz5zWJWr17d5rlra2tRW1urOW5qasJvv/2Gnj17mu3HbhOR9AkhUFVVBTc3t2YfemKuDFWLWYeJyFQMUYt1at4rKyvx4osvYs+ePZDL5ViwYAHmzJmDJUuW4J133sGgQYPw8ccf6yUxYyotLW32Ub3Ozs5oaGhAeXk5XF1dW425/WNubxcfH4+lS5fqPWciovYoKioyy09fbImhajHrMBGZmj5rsU7N++uvv47MzEzMmDEDu3fvxvz587F7927U1NRg165diIiI0EtSpnD76sutq4l+P95SzJ1WbRYuXIi4uDjNsUqlgqenJ4qKiuDg4HC3aRMRtUitVsPDwwPdunUzdSo6MUQtZh0mIlMxRC3WqXn/+uuvsWXLFowaNQp/+ctfcO+992LAgAF3vHTE3Lm4uDRbtSkrK4ONjQ169uzZZsztK0C3UygUUCgUzcYdHBz4Q4OIDE5Kl4UYqhazDhORqemzFut08c3ly5fh7+8PAPDx8YGdnR2ef/55vSVjKqGhoUhPT9ca27NnD4KCgiCXy9uMCQsLM1qeRESWjLWYiOjOdFp5b2pq0hRQALC2tjbpTjKtqa6uxrlz5zTHBQUFyMvLg6OjIzw9PbFw4UIUFxdj69atAG7uZrB+/XrExcVh9uzZyMnJQUJCgtbOBfPmzcPIkSOxcuVKPPbYY/jyyy+xd+9eZGdnG31+RERSwFpMRGQAQgcymUyMHz9ePP744+Lxxx8XNjY2IjIyUnN862Zq+/btEwCa3WbMmCGEEGLGjBkiIiJC6zH79+8XQ4YMEba2tqJv375i06ZNzc772WefCV9fXyGXy4Wfn59ISUnROTeVSiUACJVK1ZGpERG1iznUGnOtxebw3BBR52CIeqPTPu/PPfdcu+K2bNmiw68PnYtarYZSqYRKpeK1lkRkMKw1reNzQ0TGYoh6o9NlM2zKiYiIiIhMR++f3FFWVqbvUxIREREREXRs3rt06YJff/1Vczx27FiUlJRojq9cuQJXV1f9ZUdERERERBo6Ne81NTX4/SXyBw4cwI0bN7RidLiEnoiIiIiIdKD3y2ak9IEgRERERERSovfmnYiIiIiIDEOn5l0mk2mtrN9+TEREREREhqPTVpFCCAwYMEDTsFdXV2PIkCGwsrLS3E9ERERERIbBfd6JiIiIiCRCp+Z92rRpsLHR6SFERERERKQnOl3z7ubmhldeeQVnzpwxVD5ERERERNQKnZr3+fPn46uvvkJAQABCQ0ORkJCA6upqQ+VGRERERES/o1PzvnDhQvz000/Yv38//Pz8EBsbC1dXVzz33HM4cOCAoXIkIiIiIiJ0cJ/38PBwbNmyBaWlpVi9ejXOnTuH8PBw+Pr64u2339Z3jkREREREhLv8kKauXbti1qxZyMrKwldffYXy8nIsXLhQX7kREREREdHv3FXzfv36dWzZsgUjR47Eo48+ip49e+If//iHvnIjIiIiIqLf6dC+j1lZWdiyZQs+//xzNDY24sknn8Tf//53jBw5Ut/5ERERERHRf+nUvC9fvhyJiYk4f/48goKC8M9//hPPPPMMHBwcDJUfERERERH9l07N+3vvvYfp06dj1qxZCAgIMFRORERERETUAp2ueb98+TLee+89yTTuGzduhLe3N+zs7BAYGIisrKxWY2fOnAmZTNbsNmjQIE1MYmJiizE1NTXGmA4RkeSwDhMR6ZdOK++bNm1qV9zcuXM7lIw+JScnIzY2Fhs3bsSIESPw/vvvY9y4ccjPz4enp2ez+DVr1mDFihWa44aGBgwePBhPPfWUVpyDgwN++uknrTE7OzvDTIKISMJYh4mI9E8mhBDtDfb29r7zCWUyXLhw4a6S0ofg4GAMHTpU6xeOgQMHYvLkyYiPj7/j47/44gs88cQTKCgogJeXF4CbKz6xsbGorKzscF5qtRpKpRIqlYrvFSAigzGHWsM6TESdnSHqjU4r7wUFBXeMKS4u7nAy+lJXV4fjx49jwYIFWuORkZE4ePBgu86RkJCAUaNGaX5g3FJdXQ0vLy80NjbigQcewN/+9jcMGTJEb7kTEVkC1mEiIsO4q33ef6+0tBRz587Fvffeq69Tdlh5eTkaGxvh7OysNe7s7IzS0tI7Pr6kpAS7du3C888/rzXu5+eHxMRE7Ny5E0lJSbCzs8OIESNw9uzZVs9VW1sLtVqtdSMisnSsw0REhqFT815ZWYlp06ahV69ecHNzw9q1a9HU1IS33noLPj4+yMnJwccff2yoXHUmk8m0joUQzcZakpiYiO7du2Py5Mla4yEhIZg+fToGDx6M8PBwfPrppxgwYADWrVvX6rni4+OhVCo1Nw8Pjw7NhYhIiliHiYj0S6fm/fXXX0dmZiZmzJgBR0dHzJ8/HxMnTkR2djZ27dqFo0eP4plnnjFUru3m5OQEa2vrZqs7ZWVlzVaBbieEwMcff4zo6GjY2tq2GWtlZYVhw4a1ueKzcOFCqFQqza2oqKj9EyEikijWYSIiw9Cpef/666+xZcsWvPPOO9i5cyeEEBgwYAC+++47REREGCpHndna2iIwMBDp6ela4+np6QgLC2vzsRkZGTh37hxmzZp1x68jhEBeXh5cXV1bjVEoFHBwcNC6ERFZOtZhIiLD0OkNq5cvX4a/vz8AwMfHB3Z2ds2uRzQXcXFxiI6ORlBQEEJDQ/HBBx+gsLAQMTExAG6uxBQXF2Pr1q1aj0tISEBwcHCLe9kvXboUISEh6N+/P9RqNdauXYu8vDxs2LDBKHMiIpIS1mEiIv3TqXlvamqCXC7XHFtbW6Nr1656T0ofoqKiUFFRgWXLlqGkpAQBAQFIS0vT7FpQUlKCwsJCrceoVCqkpKRgzZo1LZ6zsrISL7zwAkpLS6FUKjFkyBBkZmZi+PDhBp8PEZHUsA4TEemfTvu8W1lZYdy4cVAoFACAr776Cg8//HCzBn7Hjh36zdKCcH9hIjIG1prW8bkhImMx+T7vM2bM0DqePn26XpIgIiIiIqI706l537Jli6HyICIiIiKiO9DbhzQREREREZFhsXknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIIi27eN27cCG9vb9jZ2SEwMBBZWVmtxu7fvx8ymazZ7ccff9SKS0lJgb+/PxQKBfz9/ZGammroaRARSRbrMBGRflls856cnIzY2FgsWrQIubm5CA8Px7hx41BYWNjm43766SeUlJRobv3799fcl5OTg6ioKERHR+PEiROIjo7GlClTcPjwYUNPh4hIcliHiYj0TyaEEKZOwhCCg4MxdOhQbNq0STM2cOBATJ48GfHx8c3i9+/fj4ceeghXr15F9+7dWzxnVFQU1Go1du3apRkbO3YsevTogaSkpHblpVaroVQqoVKp4ODgoNukiIjayRxqDeswEXV2hqg3FrnyXldXh+PHjyMyMlJrPDIyEgcPHmzzsUOGDIGrqyseeeQR7Nu3T+u+nJycZuccM2ZMm+esra2FWq3WuhERWTrWYSIiw7DI5r28vByNjY1wdnbWGnd2dkZpaWmLj3F1dcUHH3yAlJQU7NixA76+vnjkkUeQmZmpiSktLdXpnAAQHx8PpVKpuXl4eNzFzIiIpIF1mIjIMGxMnYAhyWQyrWMhRLOxW3x9feHr66s5Dg0NRVFREd555x2MHDmyQ+cEgIULFyIuLk5zrFar+YODiDoN1mEiIv2yyJV3JycnWFtbN1uJKSsra7Zi05aQkBCcPXtWc+zi4qLzORUKBRwcHLRuRESWjnWYiMgwLLJ5t7W1RWBgINLT07XG09PTERYW1u7z5ObmwtXVVXMcGhra7Jx79uzR6ZxERJ0B6zARkWFY7GUzcXFxiI6ORlBQEEJDQ/HBBx+gsLAQMTExAG7+GbW4uBhbt24FAKxevRp9+/bFoEGDUFdXh//85z9ISUlBSkqK5pzz5s3DyJEjsXLlSjz22GP48ssvsXfvXmRnZ5tkjkRE5ox1mIhI/yy2eY+KikJFRQWWLVuGkpISBAQEIC0tDV5eXgCAkpISrb2G6+rq8Morr6C4uBj29vYYNGgQvv76a4wfP14TExYWhu3bt+ONN97Am2++iX79+iE5ORnBwcFGnx8RkbljHSYi0j+L3efdXHF/YSIyBtaa1vG5ISJj4T7vRERERESdGJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCTCopv3jRs3wtvbG3Z2dggMDERWVlarsTt27MDo0aPRq1cvODg4IDQ0FN98841WTGJiImQyWbNbTU2NoadCRCRJrMNERPplsc17cnIyYmNjsWjRIuTm5iI8PBzjxo1DYWFhi/GZmZkYPXo00tLScPz4cTz00EOYNGkScnNzteIcHBxQUlKidbOzszPGlIiIJIV1mIhI/2RCCGHqJAwhODgYQ4cOxaZNmzRjAwcOxOTJkxEfH9+ucwwaNAhRUVF46623ANxc8YmNjUVlZWWH81Kr1VAqlVCpVHBwcOjweYiI2mIOtYZ1mIg6O0PUG4tcea+rq8Px48cRGRmpNR4ZGYmDBw+26xxNTU2oqqqCo6Oj1nh1dTW8vLzg7u6OiRMnNlsRul1tbS3UarXWjYjI0rEOExEZhkU27+Xl5WhsbISzs7PWuLOzM0pLS9t1jnfffRfXrl3DlClTNGN+fn5ITEzEzp07kZSUBDs7O4wYMQJnz55t9Tzx8fFQKpWam4eHR8cmRUQkIazDRESGYZHN+y0ymUzrWAjRbKwlSUlJWLJkCZKTk9G7d2/NeEhICKZPn47BgwcjPDwcn376KQYMGIB169a1eq6FCxdCpVJpbkVFRR2fEBGRxLAOExHpl42pEzAEJycnWFtbN1vdKSsra7YKdLvk5GTMmjULn332GUaNGtVmrJWVFYYNG9bmio9CoYBCoWh/8kREFoB1mIjIMCxy5d3W1haBgYFIT0/XGk9PT0dYWFirj0tKSsLMmTOxbds2TJgw4Y5fRwiBvLw8uLq63nXORESWhHWYiMgwLHLlHQDi4uIQHR2NoKAghIaG4oMPPkBhYSFiYmIA3PwzanFxMbZu3Qrg5g+MZ599FmvWrEFISIhmtcje3h5KpRIAsHTpUoSEhKB///5Qq9VYu3Yt8vLysGHDBtNMkojIjLEOExHpn8U271FRUaioqMCyZctQUlKCgIAApKWlwcvLCwBQUlKitdfw+++/j4aGBrz44ot48cUXNeMzZsxAYmIiAKCyshIvvPACSktLoVQqMWTIEGRmZmL48OFGnRsRkRSwDhMR6Z/F7vNurri/MBEZA2tN6/jcEJGxcJ93IiIiIqJOjM07EREREZFEsHknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIINu9ERERERBJh0c37xo0b4e3tDTs7OwQGBiIrK6vN+IyMDAQGBsLOzg4+Pj7YvHlzs5iUlBT4+/tDoVDA398fqamphkqfiEjyWIeJiPTLYpv35ORkxMbGYtGiRcjNzUV4eDjGjRuHwsLCFuMLCgowfvx4hIeHIzc3F6+//jrmzp2LlJQUTUxOTg6ioqIQHR2NEydOIDo6GlOmTMHhw4eNNS0iIslgHSYi0j+ZEEKYOglDCA4OxtChQ7Fp0ybN2MCBAzF58mTEx8c3i3/ttdewc+dOnDlzRjMWExODEydOICcnBwAQFRUFtVqNXbt2aWLGjh2LHj16ICkpqV15qdVqKJVKqFQqODg4dHR6RERtModawzpMRJ2dIeqNjV7OYmbq6upw/PhxLFiwQGs8MjISBw8ebPExOTk5iIyM1BobM2YMEhISUF9fD7lcjpycHMyfP79ZzOrVq1vNpba2FrW1tZpjlUoF4OaLSURkKLdqjKnWZ1iHiYgMU4stsnkvLy9HY2MjnJ2dtcadnZ1RWlra4mNKS0tbjG9oaEB5eTlcXV1bjWntnAAQHx+PpUuXNhv38PBo73SIiDqsoqICSqXS6F+XdZiI6H/0WYstsnm/RSaTaR0LIZqN3Sn+9nFdz7lw4ULExcVpjisrK+Hl5YXCwkKT/EA1NLVaDQ8PDxQVFVnsn6M5R+mz9PkBN1eXPT094ejoaNI8WIdNw9L/jVv6/ADO0VIYohZbZPPu5OQEa2vrZisxZWVlzVZsbnFxcWkx3sbGBj179mwzprVzAoBCoYBCoWg2rlQqLfYfKgA4ODhY9PwAztESWPr8AMDKyjT7ErAOmwdL/zdu6fMDOEdLoc9abJG7zdja2iIwMBDp6ela4+np6QgLC2vxMaGhoc3i9+zZg6CgIMjl8jZjWjsnEVFnxTpMRGQgwkJt375dyOVykZCQIPLz80VsbKzo2rWruHjxohBCiAULFojo6GhN/IULF0SXLl3E/PnzRX5+vkhISBByuVx8/vnnmpgDBw4Ia2trsWLFCnHmzBmxYsUKYWNjIw4dOtTuvFQqlQAgVCqV/iZrRix9fkJwjpbA0ucnhHnMkXXYdCx9jpY+PyE4R0thiDlabPMuhBAbNmwQXl5ewtbWVgwdOlRkZGRo7psxY4aIiIjQit+/f78YMmSIsLW1FX379hWbNm1qds7PPvtM+Pr6CrlcLvz8/ERKSopOOdXU1IjFixeLmpqaDs3J3Fn6/ITgHC2Bpc9PCPOZI+uwaVj6HC19fkJwjpbCEHO02H3eiYiIiIgsjUVe805EREREZInYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWze9Wzjxo3w9vaGnZ0dAgMDkZWV1WZ8RkYGAgMDYWdnBx8fH2zevNlImXacLnPcsWMHRo8ejV69esHBwQGhoaH45ptvjJhtx+j6Ot5y4MAB2NjY4IEHHjBsgndJ1/nV1tZi0aJF8PLygkKhQL9+/fDxxx8bKduO0XWOn3zyCQYPHowuXbrA1dUVzz33HCoqKoyUre4yMzMxadIkuLm5QSaT4YsvvrjjY6RYbzrK0msx63DrpFKHAdbilkipFpusDutt3xrS7Gn84Ycfivz8fDFv3jzRtWtX8csvv7QYf2tP43nz5on8/Hzx4YcfNtvT2NzoOsd58+aJlStXiiNHjoiff/5ZLFy4UMjlcvH9998bOfP203WOt1RWVgofHx8RGRkpBg8ebJxkO6Aj83v00UdFcHCwSE9PFwUFBeLw4cPiwIEDRsxaN7rOMSsrS1hZWYk1a9aICxcuiKysLDFo0CAxefJkI2fefmlpaWLRokUiJSVFABCpqaltxkux3nSUpddi1uHWSaUOC8Fa3BKp1WJT1WE273o0fPhwERMTozXm5+cnFixY0GL8X//6V+Hn56c19uc//1mEhIQYLMe7pescW+Lv7y+WLl2q79T0pqNzjIqKEm+88YZYvHixWf/Q0HV+u3btEkqlUlRUVBgjPb3QdY7//Oc/hY+Pj9bY2rVrhbu7u8Fy1Kf2/NCQYr3pKEuvxazDrZNKHRaCtbglUq7FxqzDvGxGT+rq6nD8+HFERkZqjUdGRuLgwYMtPiYnJ6dZ/JgxY3Ds2DHU19cbLNeO6sgcb9fU1ISqqio4OjoaIsW71tE5btmyBefPn8fixYsNneJd6cj8du7ciaCgILz99tvo06cPBgwYgFdeeQU3btwwRso668gcw8LCcOnSJaSlpUEIgStXruDzzz/HhAkTjJGyUUit3nSUpddi1mHp12GAtbiz1mJ91RobfSfWWZWXl6OxsRHOzs5a487OzigtLW3xMaWlpS3GNzQ0oLy8HK6urgbLtyM6Msfbvfvuu7h27RqmTJliiBTvWkfmePbsWSxYsABZWVmwsTHvb6mOzO/ChQvIzs6GnZ0dUlNTUV5ejr/85S/47bffzPJay47MMSwsDJ988gmioqJQU1ODhoYGPProo1i3bp0xUjYKqdWbjrL0Wsw6LP06DLAWd9ZarK9aw5V3PZPJZFrHQohmY3eKb2ncnOg6x1uSkpKwZMkSJCcno3fv3oZKTy/aO8fGxkZMnToVS5cuxYABA4yV3l3T5TVsamqCTCbDJ598guHDh2P8+PFYtWoVEhMTzXbFB9Btjvn5+Zg7dy7eeustHD9+HLt370ZBQQFiYmKMkarRSLHedJSl12LW4f+Rah0GWItv1xlqsT5qjfn/eioRTk5OsLa2bvbbZFlZWbPfsm5xcXFpMd7GxgY9e/Y0WK4d1ZE53pKcnIxZs2bhs88+w6hRowyZ5l3RdY5VVVU4duwYcnNzMWfOHAA3C6wQAjY2NtizZw8efvhho+TeHh15DV1dXdGnTx8olUrN2MCBAyGEwKVLl9C/f3+D5qyrjswxPj4eI0aMwKuvvgoAuP/++9G1a1eEh4fj73//u1mtvHaU1OpNR1l6LWYdln4dBliLO2st1let4cq7ntja2iIwMBDp6ela4+np6QgLC2vxMaGhoc3i9+zZg6CgIMjlcoPl2lEdmSNwc6Vn5syZ2LZtm9lft6brHB0cHHDy5Enk5eVpbjExMfD19UVeXh6Cg4ONlXq7dOQ1HDFiBC5fvozq6mrN2M8//wwrKyu4u7sbNN+O6Mgcr1+/Disr7XJobW0N4H+rIlIntXrTUZZei1mHpV+HAdbizlqL9VZrdHp7K7Xp1pZICQkJIj8/X8TGxoquXbuKixcvCiGEWLBggYiOjtbE39oyaP78+SI/P18kJCSY9fZkQug+x23btgkbGxuxYcMGUVJSorlVVlaaagp3pOscb2fuuxzoOr+qqirh7u4unnzySXH69GmRkZEh+vfvL55//nlTTeGOdJ3jli1bhI2Njdi4caM4f/68yM7OFkFBQWL48OGmmsIdVVVVidzcXJGbmysAiFWrVonc3FzNFmyWUG86ytJrMeuw9OuwEKzFQki/FpuqDrN517MNGzYILy8vYWtrK4YOHSoyMjI0982YMUNERERoxe/fv18MGTJE2Nrair59+4pNmzYZOWPd6TLHiIgIAaDZbcaMGcZPXAe6vo6/J4UfGrrO78yZM2LUqFHC3t5euLu7i7i4OHH9+nUjZ60bXee4du1a4e/vL+zt7YWrq6uYNm2auHTpkpGzbr99+/a1+b1lKfWmoyy9FrMOS78OC8FaLPVabKo6LBPCAv4OQURERETUCfCadyIiIiIiibDI5j0zMxOTJk2Cm5sbZDIZvvjiizs+JiMjA4GBgbCzs4OPjw82b97cLCYlJQX+/v5QKBTw9/dHamqqAbInIrIMrMVERPpnkc37tWvXMHjwYKxfv75d8QUFBRg/fjzCw8ORm5uL119/HXPnzkVKSoomJicnB1FRUYiOjsaJEycQHR2NKVOm4PDhw4aaBhGRpLEWExHpn8Vf8y6TyZCamorJkye3GvPaa69h586dOHPmjGYsJiYGJ06cQE5ODgAgKioKarUau3bt0sSMHTsWPXr0QFJSksHyJyKyBKzFRET6wQ9pws2VnMjISK2xMWPGICEhAfX19ZDL5cjJycH8+fObxaxevbrNc9fW1qK2tlZz3NTUhN9++w09e/Y020/uIyLpE0KgqqoKbm5uzfZNNleGqsWsw0RkKoaoxWzeAZSWljb7tC9nZ2c0NDSgvLwcrq6urcbc/klZt4uPj8fSpUv1njMRUXsUFRWZ5Qe4tMRQtZh1mIhMTZ+1mM37f92++nLraqLfj7cUc6dVm4ULFyIuLk5zrFKp4OnpiaKiIjg4ONxt2kRELVKr1fDw8EC3bt1MnYpODFGLWYeJyFQMUYvZvANwcXFptmpTVlYGGxsb9OzZs82Y21eAbqdQKKBQKJqNOzg48IcGERmclC4LMVQtZh0mIlPTZy2WxoWQBhYaGor09HStsT179iAoKAhyubzNmLCwMKPlSURkyViLiYjuzCJX3qurq3Hu3DnNcUFBAfLy8uDo6AhPT08sXLgQxcXF2Lp1K4CbuxmsX78ecXFxmD17NnJycpCQkKC1c8G8efMwcuRIrFy5Eo899hi+/PJL7N27F9nZ2UafHxGRFLAWExEZgLBA+/btEwCa3WbMmCGEEGLGjBkiIiJC6zH79+8XQ4YMEba2tqJv375i06ZNzc772WefCV9fXyGXy4Wfn59ISUnROTeVSiUACJVK1ZGpERG1iznUGnOtxebw3BBR52CIemPx+7ybG7VaDaVSCZVKxWstichgWGtax+eGiIzFEPWG17wTEREREUkEm3ciIiIiIolg805EREREJBFs3omIiIiIJILNOxERERGRRLB5JyIiIiKSCDbvREREREQSweadiIiIiEgi2LwTEREREUkEm3ciIiIiIolg805EREREJBFs3omIiIiIJILNOxERERGRRLB5JyIiIiKSCDbvREREREQSweadiIiIiEgi2LwTEREREUkEm3ciIiIiIolg805EREREJBFs3omIiIiIJMKim/eNGzfC29sbdnZ2CAwMRFZWVquxM2fOhEwma3YbNGiQJiYxMbHFmJqaGmNMh4hIcliHiYj0y2Kb9+TkZMTGxmLRokXIzc1FeHg4xo0bh8LCwhbj16xZg5KSEs2tqKgIjo6OeOqpp7TiHBwctOJKSkpgZ2dnjCkREUkK6zARkf5ZbPO+atUqzJo1C88//zwGDhyI1atXw8PDA5s2bWoxXqlUwsXFRXM7duwYrl69iueee04rTiaTacW5uLgYYzpERJLDOkxEpH8W2bzX1dXh+PHjiIyM1BqPjIzEwYMH23WOhIQEjBo1Cl5eXlrj1dXV8PLygru7OyZOnIjc3Nw2z1NbWwu1Wq11IyKydKzDRESGYZHNe3l5ORobG+Hs7Kw17uzsjNLS0js+vqSkBLt27cLzzz+vNe7n54fExETs3LkTSUlJsLOzw4gRI3D27NlWzxUfHw+lUqm5eXh4dGxSREQSwjpMRGQYFtm83yKTybSOhRDNxlqSmJiI7t27Y/LkyVrjISEhmD59OgYPHozw8HB8+umnGDBgANatW9fquRYuXAiVSqW5FRUVdWguRERSxDpMRKRfNqZOwBCcnJxgbW3dbHWnrKys2SrQ7YQQ+PjjjxEdHQ1bW9s2Y62srDBs2LA2V3wUCgUUCkX7kycisgCsw0REhmGRK++2trYIDAxEenq61nh6ejrCwsLafGxGRgbOnTuHWbNm3fHrCCGQl5cHV1fXu8qXiMjSsA4TERmGRa68A0BcXByio6MRFBSE0NBQfPDBBygsLERMTAyAm39GLS4uxtatW7Uel5CQgODgYAQEBDQ759KlSxESEoL+/ftDrVZj7dq1yMvLw4YNG4wyJyIiKWEdJiLSP4tt3qOiolBRUYFly5ahpKQEAQEBSEtL0+xaUFJS0myvYZVKhZSUFKxZs6bFc1ZWVuKFF15AaWkplEolhgwZgszMTAwfPtzg8yEikhrWYSIi/ZMJIYSpk+hM1Go1lEolVCoVHBwcTJ0OEVko1prW8bkhImMxRL2xyGveiYiIiIgsEZt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCTCopv3jRs3wtvbG3Z2dggMDERWVlarsfv374dMJmt2+/HHH7XiUlJS4O/vD4VCAX9/f6Smphp6GkREksU6TESkXxbbvCcnJyM2NhaLFi1Cbm4uwsPDMW7cOBQWFrb5uJ9++gklJSWaW//+/TX35eTkICoqCtHR0Thx4gSio6MxZcoUHD582NDTISKSHNZhIiL9kwkhhKmTMITg4GAMHToUmzZt0owNHDgQkydPRnx8fLP4/fv346GHHsLVq1fRvXv3Fs8ZFRUFtVqNXbt2acbGjh2LHj16ICkpqV15qdVqKJVKqFQqODg46DYpIqJ2ModawzpMRJ2dIeqNRa6819XV4fjx44iMjNQaj4yMxMGDB9t87JAhQ+Dq6opHHnkE+/bt07ovJyen2TnHjBnT5jlra2uhVqu1bkRElo51mIjIMCyyeS8vL0djYyOcnZ21xp2dnVFaWtriY1xdXfHBBx8gJSUFO3bsgK+vLx555BFkZmZqYkpLS3U6JwDEx8dDqVRqbh4eHncxMyIiaWAdJiIyDBtTJ2BIMplM61gI0WzsFl9fX/j6+mqOQ0NDUVRUhHfeeQcjR47s0DkBYOHChYiLi9Mcq9Vq/uAgok6DdZiISL8scuXdyckJ1tbWzVZiysrKmq3YtCUkJARnz57VHLu4uOh8ToVCAQcHB60bEZGlYx0mIjIMi2zebW1tERgYiPT0dK3x9PR0hIWFtfs8ubm5cHV11RyHhoY2O+eePXt0OicRUWfAOkxEZBgWe9lMXFwcoqOjERQUhNDQUHzwwQcoLCxETEwMgJt/Ri0uLsbWrVsBAKtXr0bfvn0xaNAg1NXV4T//+Q9SUlKQkpKiOee8efMwcuRIrFy5Eo899hi+/PJL7N27F9nZ2SaZIxGROWMdJiLSP4tt3qOiolBRUYFly5ahpKQEAQEBSEtLg5eXFwCgpKREa6/huro6vPLKKyguLoa9vT0GDRqEr7/+GuPHj9fEhIWFYfv27XjjjTfw5ptvol+/fkhOTkZwcLDR50dEZO5Yh4mI9M9i93k3V9xfmIiMgbWmdXxuiMhYuM87EREREVEnxuadiIiIiEgi2LwTEREREUkEm3ciIiIiIolg805EREREJBFs3omIiIiIJILNOxERERGRRLB5JyIiIiKSCDbvREREREQSweadiIiIiEgi2LwTEREREUkEm3ciIiIiIolg805EREREJBFs3omIiIiIJILNOxERERGRRLB5JyIiIiKSCDbvREREREQSweadiIiIiEgi2LwTEREREUkEm3ciIiIiIomw6OZ948aN8Pb2hp2dHQIDA5GVldVq7I4dOzB69Gj06tULDg4OCA0NxTfffKMVk5iYCJlM1uxWU1Nj6KkQEUkS6zARkX5ZbPOenJyM2NhYLFq0CLm5uQgPD8e4ceNQWFjYYnxmZiZGjx6NtLQ0HD9+HA899BAmTZqE3NxcrTgHBweUlJRo3ezs7IwxJSIiSWEdJiLSP5kQQpg6CUMIDg7G0KFDsWnTJs3YwIEDMXnyZMTHx7frHIMGDUJUVBTeeustADdXfGJjY1FZWdnhvNRqNZRKJVQqFRwcHDp8HiKitphDrWEdJqLOzhD1xiJX3uvq6nD8+HFERkZqjUdGRuLgwYPtOkdTUxOqqqrg6OioNV5dXQ0vLy+4u7tj4sSJzVaEbldbWwu1Wq11IyKydKzDRESGYZHNe3l5ORobG+Hs7Kw17uzsjNLS0nad491338W1a9cwZcoUzZifnx8SExOxc+dOJCUlwc7ODiNGjMDZs2dbPU98fDyUSqXm5uHh0bFJERFJCOswEZFhWGTzfotMJtM6FkI0G2tJUlISlixZguTkZPTu3VszHhISgunTp2Pw4MEIDw/Hp59+igEDBmDdunWtnmvhwoVQqVSaW1FRUccnREQkMazDRET6ZWPqBAzByckJ1tbWzVZ3ysrKmq0C3S45ORmzZs3CZ599hlGjRrUZa2VlhWHDhrW54qNQKKBQKNqfPBGRBWAdJiIyDItcebe1tUVgYCDS09O1xtPT0xEWFtbq45KSkjBz5kxs27YNEyZMuOPXEUIgLy8Prq6ud50zEZElYR0mIjIMi1x5B4C4uDhER0cjKCgIoaGh+OCDD1BYWIiYmBgAN/+MWlxcjK1btwK4+QPj2WefxZo1axASEqJZLbK3t4dSqQQALF26FCEhIejfvz/UajXWrl2LvLw8bNiwwTSTJCIyY6zDRET6Z7HNe1RUFCoqKrBs2TKUlJQgICAAaWlp8PLyAgCUlJRo7TX8/vvvo6GhAS+++CJefPFFzfiMGTOQmJgIAKisrMQLL7yA0tJSKJVKDBkyBJmZmRg+fLhR50ZEJAWsw0RE+mex+7ybK+4vTETGwFrTOj43RGQs3OediIiIiKgTY/NORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkUSweSciIiIikgg270REREREEsHmnYiIiIhIIti8ExERERFJBJt3IiIiIiKJYPNORERERCQRbN6JiIiIiCSCzTsRERERkURYdPO+ceNGeHt7w87ODoGBgcjKymozPiMjA4GBgbCzs4OPjw82b97cLCYlJQX+/v5QKBTw9/dHamqqodInIpI81mEiIv2y2OY9OTkZsbGxWLRoEXJzcxEeHo5x48ahsLCwxfiCggKMHz8e4eHhyM3Nxeuvv465c+ciJSVFE5OTk4OoqChER0fjxIkTiI6OxpQpU3D48GFjTYuISDJYh4mI9E8mhBCmTsIQgoODMXToUGzatEkzNnDgQEyePBnx8fHN4l977TXs3LkTZ86c0YzFxMTgxIkTyMnJAQBERUVBrVZj165dmpixY8eiR48eSEpKaldearUaSqUSKpUKDg4OHZ0eEVGbzKHWsA4TUWdniHpjo5ezmJm6ujocP34cCxYs0BqPjIzEwYMHW3xMTk4OIiMjtcbGjBmDhIQE1NfXQy6XIycnB/Pnz28Ws3r16lZzqa2tRW1treZYpVIBuPliEhEZyq0aY6r1GdZhIiLD1GKLbN7Ly8vR2NgIZ2dnrXFnZ2eUlpa2+JjS0tIW4xsaGlBeXg5XV9dWY1o7JwDEx8dj6dKlzcY9PDzaOx0iog6rqKiAUqk0+tdlHSYi+h991mKLbN5vkclkWsdCiGZjd4q/fVzXcy5cuBBxcXGa48rKSnh5eaGwsNAkP1ANTa1Ww8PDA0VFRRb752jOUfosfX7AzdVlT09PODo6mjQP1mHTsPR/45Y+P4BztBSGqMUW2bw7OTnB2tq62UpMWVlZsxWbW1xcXFqMt7GxQc+ePduMae2cAKBQKKBQKJqNK5VKi/2HCgAODg4WPT+Ac7QElj4/ALCyMs2+BKzD5sHS/41b+vwAztFS6LMWW+RuM7a2tggMDER6errWeHp6OsLCwlp8TGhoaLP4PXv2ICgoCHK5vM2Y1s5JRNRZsQ4TERmIsFDbt28XcrlcJCQkiPz8fBEbGyu6du0qLl68KIQQYsGCBSI6OloTf+HCBdGlSxcxf/58kZ+fLxISEoRcLheff/65JubAgQPC2tparFixQpw5c0asWLFC2NjYiEOHDrU7L5VKJQAIlUqlv8maEUufnxCcoyWw9PkJYR5zZB02HUufo6XPTwjO0VIYYo4W27wLIcSGDRuEl5eXsLW1FUOHDhUZGRma+2bMmCEiIiK04vfv3y+GDBkibG1tRd++fcWmTZuanfOzzz4Tvr6+Qi6XCz8/P5GSkqJTTjU1NWLx4sWipqamQ3Myd5Y+PyE4R0tg6fMTwnzmyDpsGpY+R0ufnxCco6UwxBwtdp93IiIiIiJLY5HXvBMRERERWSI270REREREEsHmnYiIiIhIIti8ExERERFJBJt3Pdu4cSO8vb1hZ2eHwMBAZGVltRmfkZGBwMBA2NnZwcfHB5s3bzZSph2nyxx37NiB0aNHo1evXnBwcEBoaCi++eYbI2bbMbq+jrccOHAANjY2eOCBBwyb4F3SdX61tbVYtGgRvLy8oFAo0K9fP3z88cdGyrZjdJ3jJ598gsGDB6NLly5wdXXFc889h4qKCiNlq7vMzExMmjQJbm5ukMlk+OKLL+74GCnWm46y9FrMOtw6qdRhgLW4JVKqxSarw3rbt4Y0exp/+OGHIj8/X8ybN0907dpV/PLLLy3G39rTeN68eSI/P198+OGHzfY0Nje6znHevHli5cqV4siRI+Lnn38WCxcuFHK5XHz//fdGzrz9dJ3jLZWVlcLHx0dERkaKwYMHGyfZDujI/B599FERHBws0tPTRUFBgTh8+LA4cOCAEbPWja5zzMrKElZWVmLNmjXiwoULIisrSwwaNEhMnjzZyJm3X1pamli0aJFISUkRAERqamqb8VKsNx1l6bWYdbh1UqnDQrAWt0RqtdhUdZjNux4NHz5cxMTEaI35+fmJBQsWtBj/17/+Vfj5+WmN/fnPfxYhISEGy/Fu6TrHlvj7+4ulS5fqOzW96egco6KixBtvvCEWL15s1j80dJ3frl27hFKpFBUVFcZITy90neM///lP4ePjozW2du1a4e7ubrAc9ak9PzSkWG86ytJrMetw66RSh4VgLW6JlGuxMeswL5vRk7q6Ohw/fhyRkZFa45GRkTh48GCLj8nJyWkWP2bMGBw7dgz19fUGy7WjOjLH2zU1NaGqqgqOjo6GSPGudXSOW7Zswfnz57F48WJDp3hXOjK/nTt3IigoCG+//Tb69OmDAQMG4JVXXsGNGzeMkbLOOjLHsLAwXLp0CWlpaRBC4MqVK/j8888xYcIEY6RsFFKrNx1l6bWYdVj6dRhgLe6stVhftcZG34l1VuXl5WhsbISzs7PWuLOzM0pLS1t8TGlpaYvxDQ0NKC8vh6urq8Hy7YiOzPF27777Lq5du4YpU6YYIsW71pE5nj17FgsWLEBWVhZsbMz7W6oj87tw4QKys7NhZ2eH1NRUlJeX4y9/+Qt+++03s7zWsiNzDAsLwyeffIKoqCjU1NSgoaEBjz76KNatW2eMlI1CavWmoyy9FrMOS78OA6zFnbUW66vWcOVdz2QymdaxEKLZ2J3iWxo3J7rO8ZakpCQsWbIEycnJ6N27t6HS04v2zrGxsRFTp07F0qVLMWDAAGOld9d0eQ2bmpogk8nwySefYPjw4Rg/fjxWrVqFxMREs13xAXSbY35+PubOnYu33noLx48fx+7du1FQUICYmBhjpGo0Uqw3HWXptZh1+H+kWocB1uLbdYZarI9aY/6/nkqEk5MTrK2tm/02WVZW1uy3rFtcXFxajLexsUHPnj0NlmtHdWSOtyQnJ2PWrFn47LPPMGrUKEOmeVd0nWNVVRWOHTuG3NxczJkzB8DNAiuEgI2NDfbs2YOHH37YKLm3R0deQ1dXV/Tp0wdKpVIzNnDgQAghcOnSJfTv39+gOeuqI3OMj4/HiBEj8OqrrwIA7r//fnTt2hXh4eH4+9//blYrrx0ltXrTUZZei1mHpV+HAdbizlqL9VVruPKuJ7a2tggMDER6errWeHp6OsLCwlp8TGhoaLP4PXv2ICgoCHK53GC5dlRH5gjcXOmZOXMmtm3bZvbXrek6RwcHB5w8eRJ5eXmaW0xMDHx9fZGXl4fg4GBjpd4uHXkNR4wYgcuXL6O6uloz9vPPP8PKygru7u4GzbcjOjLH69evw8pKuxxaW1sD+N+qiNRJrd50lKXXYtZh6ddhgLW4s9ZivdUand7eSm26tSVSQkKCyM/PF7GxsaJr167i4sWLQgghFixYIKKjozXxt7YMmj9/vsjPzxcJCQlmvT2ZELrPcdu2bcLGxkZs2LBBlJSUaG6VlZWmmsId6TrH25n7Lge6zq+qqkq4u7uLJ598Upw+fVpkZGSI/v37i+eff95UU7gjXee4ZcsWYWNjIzZu3CjOnz8vsrOzRVBQkBg+fLippnBHVVVVIjc3V+Tm5goAYtWqVSI3N1ezBZsl1JuOsvRazDos/TosBGuxENKvxaaqw2ze9WzDhg3Cy8tL2NraiqFDh4qMjAzNfTNmzBARERFa8fv37xdDhgwRtra2om/fvmLTpk1Gzlh3uswxIiJCAGh2mzFjhvET14Gur+PvSeGHhq7zO3PmjBg1apSwt7cX7u7uIi4uTly/ft3IWetG1zmuXbtW+Pv7C3t7e+Hq6iqmTZsmLl26ZOSs22/fvn1tfm9ZSr3pKEuvxazD0q/DQrAWS70Wm6oOy4SwgL9DEBERERF1ArzmnYiIiIhIIiyyec/MzMSkSZPg5uYGmUyGL7744o6PycjIQGBgIOzs7ODj44PNmzc3i0lJSYG/vz8UCgX8/f2RmppqgOyJiCwDazERkf5ZZPN+7do1DB48GOvXr29XfEFBAcaPH4/w8HDk5ubi9ddfx9y5c5GSkqKJycnJQVRUFKKjo3HixAlER0djypQpOHz4sKGmQUQkaazFRET6Z/HXvMtkMqSmpmLy5Mmtxrz22mvYuXMnzpw5oxmLiYnBiRMnkJOTAwCIioqCWq3Grl27NDFjx45Fjx49kJSUZLD8iYgsAWsxEZF+8EOacHMlJzIyUmtszJgxSEhIQH19PeRyOXJycjB//vxmMatXr27z3LW1taitrdUcNzU14bfffkPPnj3N9pP7iEj6hBCoqqqCm5tbs32TzZWhajHrMBGZiiFqMZt3AKWlpc0+7cvZ2RkNDQ0oLy+Hq6trqzG3f1LW7eLj47F06VK950xE1B5FRUVm+QEuLTFULWYdJiJT02ctZvP+X7evvty6muj34y3F3GnVZuHChYiLi9Mcq1QqeHp6oqioCA4ODnebNhFRi9RqNTw8PNCtWzdTp6ITQ9Ri1mEiMhVD1GI27wBcXFyardqUlZXBxsYGPXv2bDPm9hWg2ykUCigUimbjDg4O/KFBRAYnpctCDFWLWYeJyNT0WYulcSGkgYWGhiI9PV1rbM+ePQgKCoJcLm8zJiwszGh5EhFZMtZiIqI7s8iV9+rqapw7d05zXFBQgLy8PDg6OsLT0xMLFy5EcXExtm7dCuDmbgbr169HXFwcZs+ejZycHCQkJGjtXDBv3jyMHDkSK1euxGOPPYYvv/wSe/fuRXZ2ttHnR0QkBazFREQGICzQvn37BIBmtxkzZgghhJgxY4aIiIjQesz+/fvFkCFDhK2trejbt6/YtGlTs/N+9tlnwtfXV8jlcuHn5ydSUlJ0zk2lUgkAQqVSdWRqRETtYg61xlxrsTk8N0TUORii3lj8Pu/mRq1WQ6lUQqVS8VpLIjIY1prW8bkhImMxRL3hNe9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiiWDzTkREREQkEWzeiYiIiIgkgs07EREREZFEsHknIiIiIpIINu9ERERERBLB5p2IiIiISCLYvBMRERERSQSbdyIiIiIiibDo5n3jxo3w9vaGnZ0dAgMDkZWV1WrszJkzIZPJmt0GDRqkiUlMTGwxpqamxhjTISKSHNZhIiL9stjmPTk5GbGxsVi0aBFyc3MRHh6OcePGobCwsMX4NWvWoKSkRHMrKiqCo6MjnnrqKa04BwcHrbiSkhLY2dkZY0pERJLCOkxEpH8W27yvWrUKs2bNwvPPP4+BAwdi9erV8PDwwKZNm1qMVyqVcHFx0dyOHTuGq1ev4rnnntOKk8lkWnEuLi7GmA4RkeSwDhMR6Z9FNu91dXU4fvw4IiMjtcYjIyNx8ODBdp0jISEBo0aNgpeXl9Z4dXU1vLy84O7ujokTJyI3N7fN89TW1kKtVmvdiIgsHeswEZFhWGTzXl5ejsbGRjg7O2uNOzs7o7S09I6PLykpwa5du/D8889rjfv5+SExMRE7d+5EUlIS7OzsMGLECJw9e7bVc8XHx0OpVGpuHh4eHZsUEZGEsA4TERmGRTbvt8hkMq1jIUSzsZYkJiaie/fumDx5stZ4SEgIpk+fjsGDByM8PByffvopBgwYgHXr1rV6roULF0KlUmluRUVFHZoLEZEUsQ4TEemXjakTMAQnJydYW1s3W90pKytrtgp0OyEEPv74Y0RHR8PW1rbNWCsrKwwbNqzNFR+FQgGFQtH+5ImILADrMBGRYVjkyrutrS0CAwORnp6uNZ6eno6wsLA2H5uRkYFz585h1qxZd/w6Qgjk5eXB1dX1rvIlIrI0rMNERIZhkSvvAPD/27n/mKru+4/jL+SnNbl3tbbIikNtKIpmHVzGL0NNV8Vq184/Fm+yhdmlS0fWpSLpNhxdLWYJ7bY2lQbobFiNmVLWoqvJcJUlFVCZS9nFLMOtptqhDmZgg0vbCRU/3z+Md9/bC8i93nvhHJ6P5PxxP3zu8f3m6uu+OVxPeXm5SkpKlJOTo4KCAu3Zs0e9vb0qLS2VdP3XqJcuXdK+ffv8ntfQ0KC8vDytXr064JxVVVXKz89Xenq6vF6vampq1N3drdra2qj0BABWQg4DQPjZdnh3u90aHBzUrl271NfXp9WrV6ulpcV314K+vr6Aew0PDw+rublZu3fvnvCcQ0NDeuKJJ9Tf3y+n06msrCy1t7crNzc34v0AgNWQwwAQfjHGGDPTRcwlXq9XTqdTw8PDcjgcM10OAJsiaybH9wZAtEQib2z5mXcAAADAjhjeAQAAAItgeAcAAAAsguEdAAAAsAiGdwAAAMAiGN4BAAAAi2B4BwAAACyC4R0AAACwCIZ3AAAAwCIY3gEAAACLYHgHAAAALILhHQAAALAIhncAAADAIhjeAQAAAItgeAcAAAAsguEdAAAAsAiGdwAAAMAiGN4BAAAAi2B4BwAAACyC4R0AAACwCFsP73V1dVq2bJmSkpLkcrnU0dEx6d5jx44pJiYm4Pjb3/7mt6+5uVmZmZlKTExUZmamDh06FOk2AMCyyGEACC/bDu9NTU0qKytTZWWlPB6PioqKtHHjRvX29k75vL///e/q6+vzHenp6b6vdXZ2yu12q6SkRKdPn1ZJSYm2bNmiU6dORbodALAcchgAwi/GGGNmuohIyMvLU3Z2turr631rK1eu1ObNm1VdXR2w/9ixY3rggQf0n//8R5/73OcmPKfb7ZbX69WRI0d8aw899JBuv/12NTY2Tqsur9crp9Op4eFhORyO4JoCgGmaDVlDDgOY6yKRN7a88j42Nqauri4VFxf7rRcXF+vkyZNTPjcrK0spKSl68MEH9e677/p9rbOzM+CcGzZsmPKco6Oj8nq9fgcA2B05DACRYcvhfWBgQOPj40pOTvZbT05OVn9//4TPSUlJ0Z49e9Tc3KyDBw8qIyNDDz74oNrb2317+vv7gzqnJFVXV8vpdPqOJUuW3EJnAGAN5DAAREbcTBcQSTExMX6PjTEBazdkZGQoIyPD97igoEAXLlzQL37xC91///0hnVOSduzYofLyct9jr9fLGweAOYMcBoDwsuWV90WLFik2NjbgSszly5cDrthMJT8/X2fPnvU9Xrx4cdDnTExMlMPh8DsAwO7IYQCIDFsO7wkJCXK5XGptbfVbb21tVWFh4bTP4/F4lJKS4ntcUFAQcM6jR48GdU4AmAvIYQCIDNt+bKa8vFwlJSXKyclRQUGB9uzZo97eXpWWlkq6/mvUS5cuad++fZKkl19+WUuXLtWqVas0NjamX//612publZzc7PvnNu2bdP999+vF154QV/72tf09ttv6w9/+IOOHz8+Iz0CwGxGDgNA+Nl2eHe73RocHNSuXbvU19en1atXq6WlRWlpaZKkvr4+v3sNj42N6emnn9alS5c0f/58rVq1Sr/73e+0adMm357CwkK98cYbeuaZZ/STn/xE99xzj5qampSXlxf1/gBgtiOHASD8bHuf99mK+wsDiAayZnJ8bwBEC/d5BwAAAOYwhncAAADAIhjeAQAAAItgeAcAAAAsguEdAAAAsAiGdwAAAMAiGN4BAAAAi2B4BwAAACyC4R0AAACwCIZ3AAAAwCIY3gEAAACLYHgHAAAALILhHQAAALAIhncAAADAIhjeAQAAAItgeAcAAAAsguEdAAAAsAiGdwAAAMAiGN4BAAAAi7D18F5XV6dly5YpKSlJLpdLHR0dk+49ePCg1q9frzvvvFMOh0MFBQV65513/Pbs3btXMTExAceVK1ci3QoAWBI5DADhZdvhvampSWVlZaqsrJTH41FRUZE2btyo3t7eCfe3t7dr/fr1amlpUVdXlx544AE98sgj8ng8fvscDof6+vr8jqSkpGi0BACWQg4DQPjFGGPMTBcRCXl5ecrOzlZ9fb1vbeXKldq8ebOqq6undY5Vq1bJ7Xbr2WeflXT9ik9ZWZmGhoZCrsvr9crpdGp4eFgOhyPk8wDAVGZD1pDDAOa6SOSNLa+8j42NqaurS8XFxX7rxcXFOnny5LTOce3aNY2MjGjhwoV+6x999JHS0tKUmpqqr371qwFXhAAA5DAARIoth/eBgQGNj48rOTnZbz05OVn9/f3TOseLL76ojz/+WFu2bPGtrVixQnv37tXhw4fV2NiopKQkrVmzRmfPnp30PKOjo/J6vX4HANgdOQwAkRE30wVEUkxMjN9jY0zA2kQaGxv13HPP6e2339Zdd93lW8/Pz1d+fr7v8Zo1a5Sdna1XXnlFNTU1E56rurpaVVVVIXYAANZGDgNAeNnyyvuiRYsUGxsbcHXn8uXLAVeBPqupqUmPP/64fvOb32jdunVT7p03b56+/OUvT3nFZ8eOHRoeHvYdFy5cmH4jAGBR5DAARIYth/eEhAS5XC61trb6rbe2tqqwsHDS5zU2Nuqxxx7TgQMH9PDDD9/0zzHGqLu7WykpKZPuSUxMlMPh8DsAwO7IYQCIDNt+bKa8vFwlJSXKyclRQUGB9uzZo97eXpWWlkq6fiXm0qVL2rdvn6Trbxjf+ta3tHv3buXn5/uuFs2fP19Op1OSVFVVpfz8fKWnp8vr9aqmpkbd3d2qra2dmSYBYBYjhwEg/Gw7vLvdbg0ODmrXrl3q6+vT6tWr1dLSorS0NElSX1+f372Gf/nLX+rq1at68skn9eSTT/rWt27dqr1790qShoaG9MQTT6i/v19Op1NZWVlqb29Xbm5uVHsDACsghwEg/Gx7n/fZivsLA4gGsmZyfG8ARAv3eQcAAADmMIZ3AAAAwCIY3gEAAACLYHgHAAAALILhHQAAALAIhncAAADAIhjeAQAAAItgeAcAAAAsguEdAAAAsAiGdwAAAMAiGN4BAAAAi2B4BwAAACyC4R0AAACwCIZ3AAAAwCIY3gEAAACLYHgHAAAALILhHQAAALAIhncAAADAIhjeAQAAAItgeAcAAAAswtbDe11dnZYtW6akpCS5XC51dHRMub+trU0ul0tJSUlavny5Xn311YA9zc3NyszMVGJiojIzM3Xo0KFIlQ8AlkcOA0B42XZ4b2pqUllZmSorK+XxeFRUVKSNGzeqt7d3wv3nz5/Xpk2bVFRUJI/Hox//+Md66qmn1Nzc7NvT2dkpt9utkpISnT59WiUlJdqyZYtOnToVrbYAwDLIYQAIvxhjjJnpIiIhLy9P2dnZqq+v962tXLlSmzdvVnV1dcD+H/3oRzp8+LDOnDnjWystLdXp06fV2dkpSXK73fJ6vTpy5Ihvz0MPPaTbb79djY2N06rL6/XK6XRqeHhYDocj1PYAYEqzIWvIYQBzXSTyJi4sZ5llxsbG1NXVpYqKCr/14uJinTx5csLndHZ2qri42G9tw4YNamho0Keffqr4+Hh1dnZq+/btAXtefvnlSWsZHR3V6Oio7/Hw8LCk6y8mAETKjYyZqesz5DAARCaLbTm8DwwMaHx8XMnJyX7rycnJ6u/vn/A5/f39E+6/evWqBgYGlJKSMumeyc4pSdXV1aqqqgpYX7JkyXTbAYCQDQ4Oyul0Rv3PJYcB4H/CmcW2HN5viImJ8XtsjAlYu9n+z64He84dO3aovLzc93hoaEhpaWnq7e2dkTfUSPN6vVqyZIkuXLhg219H06P12b0/6frV5S984QtauHDhjNZBDs8Mu/8dt3t/Ej3aRSSy2JbD+6JFixQbGxtwJeby5csBV2xuWLx48YT74+LidMcdd0y5Z7JzSlJiYqISExMD1p1Op23/okqSw+GwdX8SPdqB3fuTpHnzZua+BOTw7GD3v+N270+iR7sIZxbb8m4zCQkJcrlcam1t9VtvbW1VYWHhhM8pKCgI2H/06FHl5OQoPj5+yj2TnRMA5ipyGAAixNjUG2+8YeLj401DQ4Pp6ekxZWVlZsGCBebDDz80xhhTUVFhSkpKfPvPnTtnbrvtNrN9+3bT09NjGhoaTHx8vHnrrbd8e06cOGFiY2PN888/b86cOWOef/55ExcXZ/74xz9Ou67h4WEjyQwPD4ev2VnE7v0ZQ492YPf+jJkdPZLDM8fuPdq9P2Po0S4i0aNth3djjKmtrTVpaWkmISHBZGdnm7a2Nt/Xtm7datauXeu3/9ixYyYrK8skJCSYpUuXmvr6+oBzvvnmmyYjI8PEx8ebFStWmObm5qBqunLlitm5c6e5cuVKSD3Ndnbvzxh6tAO792fM7OmRHJ4Zdu/R7v0ZQ492EYkebXufdwAAAMBubPmZdwAAAMCOGN4BAAAAi2B4BwAAACyC4R0AAACwCIb3MKurq9OyZcuUlJQkl8uljo6OKfe3tbXJ5XIpKSlJy5cv16uvvhqlSkMXTI8HDx7U+vXrdeedd8rhcKigoEDvvPNOFKsNTbCv4w0nTpxQXFycvvSlL0W2wFsUbH+jo6OqrKxUWlqaEhMTdc899+hXv/pVlKoNTbA97t+/X/fdd59uu+02paSk6Nvf/rYGBwejVG3w2tvb9cgjj+jzn/+8YmJi9Nvf/vamz7Fi3oTK7llMDk/OKjkskcUTsVIWz1gOh+2+NfDd0/i1114zPT09Ztu2bWbBggXmH//4x4T7b9zTeNu2baanp8e89tprAfc0nm2C7XHbtm3mhRdeMH/605/M+++/b3bs2GHi4+PNn//85yhXPn3B9njD0NCQWb58uSkuLjb33XdfdIoNQSj9PfrooyYvL8+0traa8+fPm1OnTpkTJ05EsergBNtjR0eHmTdvntm9e7c5d+6c6ejoMKtWrTKbN2+OcuXT19LSYiorK01zc7ORZA4dOjTlfivmTajsnsXk8OSsksPGkMUTsVoWz1QOM7yHUW5uriktLfVbW7FihamoqJhw/w9/+EOzYsUKv7Xvfve7Jj8/P2I13qpge5xIZmamqaqqCndpYRNqj2632zzzzDNm586ds/pNI9j+jhw5YpxOpxkcHIxGeWERbI8///nPzfLly/3WampqTGpqasRqDKfpvGlYMW9CZfcsJocnZ5UcNoYsnoiVsziaOczHZsJkbGxMXV1dKi4u9lsvLi7WyZMnJ3xOZ2dnwP4NGzbovffe06effhqxWkMVSo+fde3aNY2MjGjhwoWRKPGWhdrj66+/rg8++EA7d+6MdIm3JJT+Dh8+rJycHP3sZz/T3XffrXvvvVdPP/20/vvf/0aj5KCF0mNhYaEuXryolpYWGWP0r3/9S2+99ZYefvjhaJQcFVbLm1DZPYvJYevnsEQWz9UsDlfWxIW7sLlqYGBA4+PjSk5O9ltPTk5Wf3//hM/p7++fcP/Vq1c1MDCglJSUiNUbilB6/KwXX3xRH3/8sbZs2RKJEm9ZKD2ePXtWFRUV6ujoUFzc7P4nFUp/586d0/Hjx5WUlKRDhw5pYGBA3/ve9/Tvf/97Vn7WMpQeCwsLtX//frndbl25ckVXr17Vo48+qldeeSUaJUeF1fImVHbPYnLY+jkskcVzNYvDlTVceQ+zmJgYv8fGmIC1m+2faH02CbbHGxobG/Xcc8+pqalJd911V6TKC4vp9jg+Pq5vfOMbqqqq0r333hut8m5ZMK/htWvXFBMTo/379ys3N1ebNm3SSy+9pL17987aKz5ScD329PToqaee0rPPPquuri79/ve/1/nz51VaWhqNUqPGinkTKrtnMTn8P1bNYYks/qy5kMXhyJrZ/+OpRSxatEixsbEBP01evnw54KesGxYvXjzh/ri4ON1xxx0RqzVUofR4Q1NTkx5//HG9+eabWrduXSTLvCXB9jgyMqL33ntPHo9H3//+9yVdD1hjjOLi4nT06FF95StfiUrt0xHKa5iSkqK7775bTqfTt7Zy5UoZY3Tx4kWlp6dHtOZghdJjdXW11qxZox/84AeSpC9+8YtasGCBioqK9NOf/nRWXXkNldXyJlR2z2Jy2Po5LJHFczWLw5U1XHkPk4SEBLlcLrW2tvqtt7a2qrCwcMLnFBQUBOw/evSocnJyFB8fH7FaQxVKj9L1Kz2PPfaYDhw4MOs/txZsjw6HQ3/5y1/U3d3tO0pLS5WRkaHu7m7l5eVFq/RpCeU1XLNmjf75z3/qo48+8q29//77mjdvnlJTUyNabyhC6fGTTz7RvHn+cRgbGyvpf1dFrM5qeRMqu2cxOWz9HJbI4rmaxWHLmqD+eyumdOOWSA0NDaanp8eUlZWZBQsWmA8//NAYY0xFRYUpKSnx7b9xy6Dt27ebnp4e09DQMKtvT2ZM8D0eOHDAxMXFmdraWtPX1+c7hoaGZqqFmwq2x8+a7Xc5CLa/kZERk5qaar7+9a+bv/71r6atrc2kp6eb73znOzPVwk0F2+Prr79u4uLiTF1dnfnggw/M8ePHTU5OjsnNzZ2pFm5qZGTEeDwe4/F4jCTz0ksvGY/H47sFmx3yJlR2z2Jy2Po5bAxZbIz1s3imcpjhPcxqa2tNWlqaSUhIMNnZ2aatrc33ta1bt5q1a9f67T927JjJysoyCQkJZunSpaa+vj7KFQcvmB7Xrl1rJAUcW7dujX7hQQj2dfz/rPCmEWx/Z86cMevWrTPz5883qamppry83HzyySdRrjo4wfZYU1NjMjMzzfz5801KSor55je/aS5evBjlqqfv3XffnfLfll3yJlR2z2Jy2Po5bAxZbPUsnqkcjjHGBr+HAAAAAOYAPvMOAAAAWATDOwAAAGARDO8AAACARTC8AwAAABbB8A4AAABYBMM7AAAAYBEM7wAAAIBFMLwDAAAAFsHwDgAAAFgEwzsAAABgEQzvAAAAgEUwvAMAAAAW8X9wODBB0m6wcgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 725x650 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%time\n",
    "\n",
    "nrows = 4\n",
    "ncols = 2\n",
    "fig, axes = plt.subplots(nrows, ncols, figsize=(7.25, 6.5))\n",
    "plt.subplots_adjust(\n",
    "    top=0.95, bottom=0.065, right=0.98, left=0.10, hspace=0.225, wspace=0.250\n",
    ")  # create some space below the plots by increasing the bottom-value\n",
    "\n",
    "for ix, ocean_name in enumerate(oceans_list):\n",
    "    row = ix // 2\n",
    "    col = ix % 2\n",
    "    for case, meta in case_dic.items():\n",
    "\n",
    "        q_name = case_meta[meta[\"grid\"]][\"flow_name\"]\n",
    "\n",
    "        if case_meta[meta[\"grid\"]][\"network_type\"] == \"vector\":\n",
    "            if ocean_name == \"global\":\n",
    "                id_list = gauge_shp1[case][\"route_id\"].values\n",
    "            else:\n",
    "                id_list = gauge_shp1[case][gauge_shp1[case][\"ocean\"] == ocean_name][\n",
    "                    \"route_id\"\n",
    "                ].values\n",
    "            reach_index = get_index_array(reachID[case], id_list)\n",
    "            dr_flow = seas_data[case][q_name].isel(seg=reach_index).sum(dim=\"seg\")\n",
    "            dr_flow.plot(ax=axes[row, col], linestyle=\"-\", lw=0.75, label=case)\n",
    "\n",
    "        elif case_meta[grid_name][\"network_type\"] == \"grid\":  # means 2d grid\n",
    "            if ocean_name == \"global\":\n",
    "                id_list = gauge_shp1[case][\"route_id\"].values\n",
    "            else:\n",
    "                id_list = gauge_shp1[case][gauge_shp1[case][\"ocean\"] == ocean_name][\n",
    "                    \"route_id\"\n",
    "                ].values\n",
    "\n",
    "            reach_index = get_index_array(reachID[case], id_list)\n",
    "            seas_data_vector = seas_data[case][q_name].stack(seg=(\"lat\", \"lon\"))\n",
    "            dr_flow = seas_data_vector.isel(seg=reach_index).sum(dim=\"seg\")\n",
    "            dr_flow.plot(ax=axes[row, col], linestyle=\"-\", lw=0.75, label=case)\n",
    "\n",
    "    # reference data\n",
    "    if obs_available:\n",
    "        if ocean_name == \"global\":\n",
    "            id_list = gauge_shp1[case][\"id\"].values\n",
    "        else:\n",
    "            id_list = gauge_shp1[case][gauge_shp1[case][\"ocean\"] == ocean_name][\n",
    "                \"id\"\n",
    "            ].values\n",
    "        gauge_index = get_index_array(ds_q[\"id\"].values, id_list)\n",
    "        dr_obs = dr_q_obs_seasonal.isel(station=gauge_index).sum(dim=\"station\")\n",
    "        dr_obs.plot(\n",
    "            ax=axes[row, col],\n",
    "            linestyle=\"None\",\n",
    "            marker=\"o\",\n",
    "            markersize=2,\n",
    "            c=\"k\",\n",
    "            label=\"D19\",\n",
    "        )\n",
    "\n",
    "    axes[row, col].set_title(\"%d %s\" % (ix + 1, ocean_name), fontsize=9)\n",
    "    axes[row, col].set_xlabel(\"\")\n",
    "    if row < 7:\n",
    "        axes[row, col].set_xticklabels(\"\")\n",
    "    if col == 0:\n",
    "        axes[row, col].set_ylabel(\"Mon. flow [m$^3$/s]\", fontsize=9)\n",
    "    else:\n",
    "        axes[row, col].set_ylabel(\"\")\n",
    "    axes[row, col].tick_params(\"both\", labelsize=\"x-small\")\n",
    "\n",
    "# Legend- make space below the plot-raise bottom. there will be an label below the second last (bottom middle) ax, thanks to the bbox_to_anchor=(x, y) with a negative y-value.\n",
    "axes[row, col].legend(\n",
    "    loc=\"center left\", bbox_to_anchor=(1.10, 0.40, 0.75, 0.1), ncol=1, fontsize=\"small\"\n",
    ")\n",
    "\n",
    "for jx in range(ix + 1, nrows * ncols):\n",
    "    row = jx // 2\n",
    "    col = jx % 2\n",
    "    fig.delaxes(axes[row][col])\n",
    "\n",
    "if figureSave:\n",
    "    plt.savefig(f\"./NB2_Fig1_ocean_discharge_season_{analysis_name}.png\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe26f035",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "pending"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if client:\n",
    "    client.shutdown()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cupid-analysis",
   "language": "python",
   "name": "cupid-analysis"
  },
  "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.11.4"
  },
  "papermill": {
   "duration": 35.200739,
   "end_time": "2025-08-21T16:28:37.687834",
   "exception": true,
   "input_path": "/glade/derecho/scratch/hannay/tmp/tmpcuy9s_ef.ipynb",
   "output_path": "/glade/work/hannay/CUPiD/examples/key_metrics/computed_notebooks/rof/global_discharge_ocean_compare_obs.ipynb",
   "parameters": {
    "CESM_output_dir": "/glade/derecho/scratch/hannay/archive",
    "analysis_name": "",
    "base_case_name": "b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.188",
    "base_case_output_dir": "/glade/derecho/scratch/gmarques/archive",
    "base_end_date": "0021-12-01",
    "base_start_date": "0002-01-01",
    "case_name": "b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.192",
    "climo_nyears": 10,
    "end_date": "0021-12-01",
    "figureSave": false,
    "grid_name": "f09_f09_mosart",
    "lc_kwargs": {
     "threads_per_worker": 1
    },
    "obs_data_dir": "/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data",
    "product": "/glade/work/hannay/CUPiD/examples/key_metrics/computed_notebooks//rof/global_discharge_ocean_compare_obs.ipynb",
    "serial": true,
    "start_date": "0002-01-01",
    "subset_kwargs": {},
    "ts_dir": null
   },
   "start_time": "2025-08-21T16:28:02.487095"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}