sill_velocities.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "14f62e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate velocities at the sill\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os, sys\n",
    "sys.path.append(os.path.abspath('./TEF'))\n",
    "import TEF_Variables as tef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c669d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "case_id = 1\n",
    "file1 = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/weak_Tide_minhs/state_' + str(format(case_id,'03d')) + '.nc')\n",
    "State = file1.isel(T=~file1.get_index(\"T\").duplicated())\n",
    "Grid = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/weak_Tide_minhs/grid_' + str(format(case_id,'03d')) + '.nc')\n",
    "\n",
    "\n",
    "\n",
    "# Confine to the range of fjord\n",
    "# state = State.isel(X=range(260), Xp1=range(261), Y=range(35,45), T=range(1,len(State.T)))\n",
    "# grid = Grid.isel(X=range(260), Xp1=range(261), Y=range(35,45))\n",
    "\n",
    "HFacC = Grid.HFacC\n",
    "areaC = Grid.HFacC * Grid.drF * Grid.dyF # Z Y X\n",
    "areaW = Grid.HFacW * Grid.drF * Grid.dyG # Z Y Xp1\n",
    "sW = np.where(HFacC[-1,40,1:] == 0)[-1][0] # Sill west edge index\n",
    "sE = np.where(HFacC[-1,40,1:] == 0)[-1][-1] # Sill east edge index\n",
    "sC = np.ceil((sW+sE)/2).astype(int) # Sill crest index (178)\n",
    "print(sC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db839d27",
   "metadata": {},
   "outputs": [],
   "source": [
    "areaW.data[:,:,-2].sum(), areaC.data[:,:,-1].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "66bc3198",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sill_vel(datapath, case_id):\n",
    "    \n",
    "    file1 = xr.open_dataset(datapath+'/state_' + str(format(case_id,'03d')) + '.nc')\n",
    "    State = file1.isel(T=~file1.get_index(\"T\").duplicated())\n",
    "    Grid = xr.open_dataset(datapath+'/grid_' + str(format(case_id,'03d')) + '.nc')\n",
    "    \n",
    "    HFacC = Grid.HFacC\n",
    "    areaC = Grid.HFacC * Grid.drF * Grid.dyF # Z Y X\n",
    "    areaW = Grid.HFacW * Grid.drF * Grid.dyG # Z Y Xp1\n",
    "    sW = np.where(HFacC[-1,40,1:] == 0)[-1][0] # Sill west edge index\n",
    "    sE = np.where(HFacC[-1,40,1:] == 0)[-1][-1] # Sill east edge index\n",
    "    sC = np.ceil((sW+sE)/2).astype(int) # Sill crest index (178)\n",
    "    \n",
    "    \n",
    "    state = State.isel(Y=range(35,45), T=range(1,len(State.T)))\n",
    "    grid = Grid.isel(Y=range(35,45))\n",
    "    s = state.S.data.mean(2) # Salinity in T,Z,X\n",
    "    u = (state.U.data[:,:,:,1:].mean(2) + state.U.data[:,:,:,:-1].mean(2)) / 2 # Along channel velocity\n",
    "    ot = state.T.data # Time in seconds\n",
    "\n",
    "\n",
    "    HFacC1 = grid.HFacC.data.mean(1)\n",
    "    dyF1 = grid.dyF.data.mean(0)\n",
    "    drF1 = grid.drF.data\n",
    "    gridA = np.broadcast_to(drF1[:, np.newaxis], HFacC1.shape) * np.broadcast_to(dyF1[np.newaxis, :], HFacC1.shape) * HFacC1 # Grid Area\n",
    "    da = np.broadcast_to(gridA[np.newaxis,:,:], u.shape)\n",
    "\n",
    "    S = state.S.data\n",
    "    U = (state.U.data[:,:,:,1:] + state.U.data[:,:,:,:-1]) / 2\n",
    "    drF = np.broadcast_to(grid.drF.data[np.newaxis, :, np.newaxis, np.newaxis], U.shape)\n",
    "    dyF = np.broadcast_to(grid.dyF.data[np.newaxis, np.newaxis, :, :], U.shape)\n",
    "    HFacC2 = np.broadcast_to(grid.HFacC.data[np.newaxis, :, :, :], U.shape)\n",
    "    DA = drF * dyF * HFacC2\n",
    "\n",
    "\n",
    "    tef_q1, tef_vel1, tef_da1, tef_qs1, tef_qs21, sbins1 = tef.process_section(U,S,DA,ot,sC,23,testing=False)\n",
    "    qin1, qout1, sin1, sout1 = tef.bulk_calc(tef_q1, tef_vel1, tef_da1, tef_qs1, tef_qs21, sbins1, ot)\n",
    "    \n",
    "    \n",
    "    U0 = 1e-4 * case_id\n",
    "    Ut = U0 * areaC.data[:,:,-1].sum() / areaC.data[:,:,sC].sum() # Maximum tidal velocity at the sill\n",
    "    t0 = 89\n",
    "    Usil = U[t0:,:,:,sC].mean(axis=(0,2))\n",
    "    Qin = qin1[t0:].mean()\n",
    "    Qout = -qout1[t0:].mean()\n",
    "    Qe = (Qin + Qout) / 2\n",
    "    Ue = Qe / (areaC.data[:,:,sC].sum()) # Exchange velocity at the sill\n",
    "    Ur = round(Ut/Ue,2)\n",
    "    \n",
    "    return Ut, Ue, Ur, Usil, Qe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77d45e94",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange(5) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1361d1ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  time 0 out of 120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:269: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sm = QSm/Qm\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:273: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2m = QS2m/Qm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  time 0 out of 120\n",
      "  time 0 out of 120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:268: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sp = QSp/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:269: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sm = QSm/Qm\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:272: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2p = QS2p/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:273: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2m = QS2m/Qm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  time 0 out of 120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:268: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sp = QSp/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:269: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sm = QSm/Qm\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:272: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2p = QS2p/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:273: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2m = QS2m/Qm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  time 0 out of 120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:268: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sp = QSp/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:269: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Sm = QSm/Qm\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:272: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2p = QS2p/Qp\n",
      "/home/1959/Scripts/Jupyter-Iceplume/TEF/TEF_Variables.py:273: RuntimeWarning: invalid value encountered in true_divide\n",
      "  S2m = QS2m/Qm\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: '/work/oceans/wbao/MITgcm_results/iceplume/6_sTide_minhs/state_004.nc'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    198\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 199\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_key\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    200\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/lru_cache.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m     52\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmove_to_end\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: [<function _open_scipy_netcdf at 0x7fe3906b79e0>, ('/work/oceans/wbao/MITgcm_results/iceplume/6_sTide_minhs/state_004.nc',), 'r', (('mmap', None), ('version', 2))]",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_29553/220920369.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcase2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 25\u001b[0;31m     \u001b[0mut2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mue2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mur2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mus2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQe2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msill_vel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcase2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/tmp/ipykernel_29553/1224564468.py\u001b[0m in \u001b[0;36msill_vel\u001b[0;34m(datapath, case_id)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0msill_vel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdatapath\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcase_id\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0mfile1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdatapath\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'/state_'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcase_id\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'03d'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'.nc'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0mState\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfile1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m~\u001b[0m\u001b[0mfile1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_index\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"T\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mduplicated\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mGrid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdatapath\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'/grid_'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcase_id\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'03d'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'.nc'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/api.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, backend_kwargs, *args, **kwargs)\u001b[0m\n\u001b[1;32m    497\u001b[0m         \u001b[0mdrop_variables\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdrop_variables\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    498\u001b[0m         \u001b[0;34m**\u001b[0m\u001b[0mdecoders\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 499\u001b[0;31m         \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    500\u001b[0m     )\n\u001b[1;32m    501\u001b[0m     ds = _dataset_from_backend_dataset(\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, mode, format, group, mmap, lock)\u001b[0m\n\u001b[1;32m    288\u001b[0m                 \u001b[0mdrop_variables\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdrop_variables\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    289\u001b[0m                 \u001b[0muse_cftime\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0muse_cftime\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 290\u001b[0;31m                 \u001b[0mdecode_timedelta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdecode_timedelta\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    291\u001b[0m             )\n\u001b[1;32m    292\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/store.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(self, store, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta)\u001b[0m\n\u001b[1;32m     22\u001b[0m         \u001b[0mdecode_timedelta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m     ):\n\u001b[0;32m---> 24\u001b[0;31m         \u001b[0mvars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mattrs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     25\u001b[0m         \u001b[0mencoding\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_encoding\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/common.py\u001b[0m in \u001b[0;36mload\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    121\u001b[0m         \"\"\"\n\u001b[1;32m    122\u001b[0m         variables = FrozenDict(\n\u001b[0;32m--> 123\u001b[0;31m             \u001b[0;34m(\u001b[0m\u001b[0m_decode_variable_name\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    124\u001b[0m         )\n\u001b[1;32m    125\u001b[0m         \u001b[0mattributes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFrozenDict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_attrs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mget_variables\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    169\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    170\u001b[0m         return FrozenDict(\n\u001b[0;32m--> 171\u001b[0;31m             \u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_store_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    172\u001b[0m         )\n\u001b[1;32m    173\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mds\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    158\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    159\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 160\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_manager\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    161\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mopen_store_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvar\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36macquire\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    179\u001b[0m             \u001b[0mAn\u001b[0m \u001b[0mopen\u001b[0m \u001b[0mfile\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mreturned\u001b[0m \u001b[0mby\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0mopener\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m         \"\"\"\n\u001b[0;32m--> 181\u001b[0;31m         \u001b[0mfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_acquire_with_cache_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneeds_lock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    182\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mfile\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    183\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    203\u001b[0m                     \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    204\u001b[0m                     \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"mode\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 205\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_opener\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    206\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"w\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    207\u001b[0m                     \u001b[0;31m# ensure file doesn't get overriden when opened again\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36m_open_scipy_netcdf\u001b[0;34m(filename, mode, mmap, version)\u001b[0m\n\u001b[1;32m     97\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnetcdf_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmmap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mversion\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mversion\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# netcdf3 message is obscure in this case\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m         \u001b[0merrmsg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/local/bin/miniconda3/envs/py37/lib/python3.7/site-packages/scipy/io/netcdf.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, filename, mode, mmap, version, maskandscale)\u001b[0m\n\u001b[1;32m    245\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfilename\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    246\u001b[0m             \u001b[0momode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'r+'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mmode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'a'\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 247\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'%sb'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0momode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    248\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmmap\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    249\u001b[0m                 \u001b[0;31m# Mmapped files on PyPy cannot be usually closed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/work/oceans/wbao/MITgcm_results/iceplume/6_sTide_minhs/state_004.nc'"
     ]
    }
   ],
   "source": [
    "#case = np.arange(5) + 1\n",
    "path1 = '/work/oceans/wbao/MITgcm_results/iceplume/1_BaseCase'\n",
    "case1 = np.array([1,2])\n",
    "ut1 = np.empty(len(case1))\n",
    "ue1 = np.empty(len(case1))\n",
    "ur1 = np.empty(len(case1))\n",
    "us1 = np.empty(shape=(len(case1), 90))\n",
    "Qe1 = np.empty(len(case1))\n",
    "\n",
    "for i in range(len(case1)):\n",
    "    \n",
    "    ut1[i], ue1[i], ur1[i], us1[i,:], Qe1[i] = sill_vel(path1, case1[i])\n",
    "\n",
    "\n",
    "path2 = '/work/oceans/wbao/MITgcm_results/iceplume/6_sTide_minhs'\n",
    "case2 = np.array([1,2,3,4,5])\n",
    "ut2 = np.empty(len(case2))\n",
    "ue2 = np.empty(len(case2))\n",
    "ur2 = np.empty(len(case2))\n",
    "us2 = np.empty(shape=(len(case2), 90))\n",
    "Qe2 = np.empty(len(case2))\n",
    "\n",
    "for j in range(len(case2)):\n",
    "    \n",
    "    ut2[j], ue2[j], ur2[j], us2[j,:], Qe2[j] = sill_vel(path2, case2[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "48da0c18",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.14927527, 0.15236062]),\n",
       " array([0.15686093, 0.2028279 , 0.2780752 , 0.528     , 0.66      ]),\n",
       " array([0.88, 1.16]),\n",
       " array([0.84      , 1.3       , 1.42      , 0.07447517, 0.07479217]),\n",
       " array([2388.40434111, 3656.6548045 ]),\n",
       " array([2509.77482834, 3245.2463326 , 4449.20315977, 2383.20547266,\n",
       "        2393.34934488]))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ue1*2, ue2*2, ur1, ur2, Qe1, Qe2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d620af20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.066, 0.132, 0.66 , 1.32 , 1.98 , 2.64 , 3.3  ]),\n",
       " array([0.14925309, 0.14924685, 0.14946376, 0.14891398, 0.14873957,\n",
       "        0.14895034, 0.14958433]),\n",
       " array([ 0.44220189,  0.88444077,  4.41578619,  8.86417793, 13.31185744,\n",
       "        17.72402778, 22.06113375]))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Qe=[2388,2388,2391,2393,2510,2393]\n",
    "Ut = np.concatenate([ut1, 10*ut2])\n",
    "Ue = np.concatenate([2*ue1, 2*ue2])\n",
    "Usill = np.concatenate([us1, us2])\n",
    "Ut, Ue, Ut/Ue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c94474a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Velocities at the sill\n",
    "U0 = 1e-3 # Tidal amplitude at OBC\n",
    "Ut = U0 * areaC.data[:,:,-1].sum() / areaC.data[:,:,sC].sum() # Maximum tidal velocity at the sill\n",
    "\n",
    "t0 = 89\n",
    "Ue = (qin1[t0:].mean()-qout1[t0:].mean()) / (2*areaC.data[:,:,sC].sum()) # Exchange velocity at the sill\n",
    "\n",
    "print(Ut,Ue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f3503b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save outputs\n",
    "Vel_ds = xr.Dataset(\n",
    "    data_vars={'Ut' : ut,\n",
    "    'Ue' : ue,\n",
    "    'Ur' : ur})\n",
    "\n",
    "outdir = \"/home/1959/Parameters/\"\n",
    "if not os.path.exists(outdir):\n",
    "    os.makedirs(outdir)\n",
    "\n",
    "Vel_ds.to_netcdf(outdir + 'wTide_maxhs_vel.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "934433be",
   "metadata": {},
   "outputs": [],
   "source": [
    "(qin1-qout1) / (2*areaC.data[:,:,sC].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "c52e8d27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAHpCAYAAAAvepG7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABzTUlEQVR4nO2deXgURfrHP5VjMrlDiElIQgKogIoKijfeNyp4iyuCJ67n4oq7Hqy64rqoeJ/rfeCNB6DgAYriLQgqSEBAEpKQhCTkIpm7fn9Mkl8Ik0kn6U73lP15nnky00f1+031vFNVXfW+QkqJjY2NjZWIMtsAGxsbm47YjsnGxsZy2I7JxsbGctiOycbGxnLYjsnGxsZy2I7JxsbGctiOycbGxnLYjsnGxsZy2I5JEYQQS4UQH3Wyb7MQ4qlO9i0RQsiW11Qd7TlXCFEuhBAajr2jnQ1zdLr2RT09TgjxohBiWW/t6Ckdr9/VZxWxHZM67A380nGjEKIfkAf8Gubcz4FDgDd0tOcUYIHUtrTg2Zbrr9Dp2ucCF+l4XF8zA2va1WfEmG2ATe8RQhQAqcDPIXbv0/J3J6fVjhop5Xc62hMFnARcpeV4KWUJUCKEqNfLhkhGSrnBbBvMxm4xqUE456PFMe2AECJKCLFdCPE3IcQDQohKIcQ2IcSNLfsvFEL8JoRoFEK8K4SI71DEAUA/4NN2ZZ4khPhKCFErhKgXQiwTQhyvXSIIIQ4RQswTQpS12LdSCHFBh2NeBM4CjmzXPbwjRFldHieEOF4I8UvLtb4SQuwVopwxQogvhBBNQohqIcQzQojkLnTsJYT4SAhR01L2GiHE1e1tU72r1hV2i0kN9gE8QGEn+4qllHXdKG8IkABcD7wO/IVg1+LeltbZYODGlr+PABcDT7Q7/xRgqZSyHkAIcRTwPvBf4LaWsscAvm7YBFAAfA08BbiAw4AXhBABKeXrLcfMAPKBNP6/xVYSoqyujssH7gP+AzQDs4C3hBAjWrunQojDgMUt2s4G+gMzCTrls8PomEewriYCbmAYkNKF9j8VtmNSg32ANVJKbyf7NLeW2p0D8ICU8hEAIcTvwAXAHsBx7b6cUwh+sdpzCvBqu88XAB9KKf/dbtsH3bQJKWXbGFjLoPqXBMfPLifoQJFSbhBC1ABR4bqnGo5LBw6TUv7ecr0o4D2CWlt/AGYC30gpz2tnVymwuMWBrepYqBAig6DjP11K2Trut1jTP+BPhN2VU4POBr4FsFeofRrKqwWebLctseXvzA4D2olATbtrDgBGAR+2O6YZOFEIMVUIkd1NW9oQQvQTQjwihCgCvC2vKcDQnpYZhk2tTqmF31r+5rXYkkBwwP4tIURM6wv4qsWu/TsptwbYDDwlhDhPCJFpgO0Rj+2YIhwhRCzBL+ZvIXYPI+g4uvu0a2/gqw4tsH0Idr2+bHftBGAQ0L5lMBbYKKVc227bv4AXCXbjSoUQH4car9HAi8B5BLtYJxAcy3oecPagrK6o7fDZ0/K39Vr9gGiCXVhvu5cbiAUGhipUShkgaHs5QdvLW6Z6jNLT+EjH7spFPrsQ/ILUhNg3FvADX3SzzH2ANzts2xcolFK6223bm+CPW/sW2Sns2FqiZXzrGiHE34BjCY4RPQccrNUgIYSzpexrpJRPtdtu1o9rLSCBO4AFIfaXdXailLIQOKvlR+Vw4B7gQyFEXovj+tNjO6bIp4Lgr/R+7TcKIdIJDl6/L6XcqrWwlidsu7Lz1IN9Otm2HdjQcq4DOI5OBn6llH7gEyHEx8BorTa1EEfQAbc5xpanX+MIOoj2eNDWitJ63E5IKbcLIb4Dhkkp7+xhGV7gMyHEA8BrBAfiQ/3A/OmwHVOEI6X0CyGeBy4XQjQC3xN8ovQ3gl2Kv3ezyL0ItoJCOaFHQmxb3e5X/oiWc9taaEKI/xHsAn4OVBNsIVxKcGxIM1LKOiHEj8BtLfOdAsBNQB07P9EqBMYLIU4n+KStTEoZqgWj9bjO+AfBge4AMAdoIPi/PwW4VUq5ruMJQoh9CD7hexPYSLBL+E/gZyml7ZRasB2TGvwdqCL4WP/alvcfAbe3TF7sDnsDTbS0gmCH2eMdB9E7PvE7BVjUobu3FphA8MlcFLAaOEdKObebdkFQ39PAywSd3GMEpx5c0+G4JwgOwD9P8Iv/b4Jdro5oPS4kUsqvhBBHtJz3CsEWXRHB/31FJ6eVt+y7Fcgh2CX8nKBzsmlB2MkI/twIIZYQ/JKfB/g1LiHprKx1wH1Syme6eV4UQae1GNgqpQw3B8jmT4D9VM4G4EyCT5T+1ptCpJRDu+uUWrit5fpH9Ob6Nupgt5j+5AghhgGtSyiKpZSVJtiQQ7BbA8F1exv72gYba2G3mDpBCJEuhPhUCPF7y99+IY4ZKIT4vGWt0+qWx+GWp702guM0G6SUyzo6JSHE8y3r5HaawawnUsqylusv64lTalmHt1YIsV4IcVOI/aJlYub6lrVv+4Uqx4po0DZcCPGtEMIthJhmho1GYDumzrkJWCyl3J3g2MdONwXBp003SCn3IDgn52ohxJ59aGNP0aINghMaT+oro3qCECIaeBw4GdgTOD9EHZwM7N7ymsKOM9oti0ZtNcB1BJ/0KYPtmDpnPPBSy/uXgNM7HiCl3CKl/KnlfQOwBsjtKwN7QZfaAKSUX2L9eTUHAuullBullB6CMaXGdzhmPPCyDPIdkNaydMbqdKlNSlkppfyR4BidMtiOqXOypJRbIOiAgLBrmoQQgwg+ev7eeNN6Tbe0WZxcgmvPWilh5x8HLcdYkUi1u9f8qecxCSEWAaEWld7azXKSgHeAqa2hPsxGL20RQKjQvR2f6Gg5xopEqt295k/tmKSUx3W2TwhRIYQYIKXc0tLsD/m0qmW90zvAq1LKdw0ytdvooS1CKGHHBbN57LxOTcsxViRS7e41dleuc+YBk1veTwZ2mqncElbkOYKxkB7oQ9t6S5faIogfgd2FEINb1upNIKivPfOASS1P5w4G6lq7shZHizY1kVLarxAvgtEIFwO/t/xNb9meQzDIPgSjMEqCyzJWtrzGmm27HtpaPr8ObCE4sFoCXGq27Z3oGQusI7iM5taWbX8F/tryXhB8urWBYFKG0WbbrKO27Ja6qSe4vKUESDHb7t6+7AmWNjY2lsPuytnY2FgO2zHZ2NhYDtsx2djYWA7bMdnY2FiOiHJMKi/WtLGx+X8ixjGZvVhTBPOnKYmtLTJRWVvEOCbMX6yp7E2ArS1SUVZbJDkmlRdr2tjYtCOS1srptlizpQk8BSAhIWH/wYMHExsbi8/nQ0qJw+HA4/EQHR2NEAKfz8fAgQPZe++9Zcf9AH6/v22bEILY2Fg8Hg8xMTFIKXfaHxMTg9frJSYmhkAgQCAQaNsfFRVFdHR02P1RUVH4fD5iY2Px+/1h94fT5HA48Hq9FBQUsN9++0mVNEkpqaioIDk5mX333Veqoqn9/vz8fEaMGCEjVdOKFSuqpJS7hPqyR5Jj0m2xppTyaYLZNhg9erRctmxZlxdvaGggOTm5y+MiEVW1TZs2DY/HwyOPdMw6pQaRXm8imOo9JJHUlTN1sWZpaakexVgSVbXNmjWLq666ymwzDEPVeoMIajFJKX1CiGuAjwnm73peSrlaCPHXlv1PEUzVPBZYTzA32sV6XT89PV2voiyHrS0yUVlbxDgmACnlAjrkiZft8tjL4Irkq424ts/nM6JYS6CqtilTptDU1MTs2bPNNsUQVK03iKyunKnU11siMKUhqKqtf//+JCQkmG2GYahab2Dnlet08Nvr9VJSUoLL5QKCcauCceHUw9ZmHaKjo0lLSyMjI4OoqPDthubmZuLj4/vIMv0RQiyXUo4OtS+iunJ9SUlJCcnJyQwaNAghRMTfBOGwtVkDKSVer5eKigpKSkrIz88Pe3xRURHDhw/vI+v6Frsr1wkul4v+/fu3/dp29esVyaiq7Y8//qCsLHJCZAshcDgc5Obmsn379i6PdzgcfWCVOdgtpjC07wLExKj7r1JVm8PhIBKHKrT+UGRkZBhsiXmo+VNpAB6Px2wTDENVbbm5ufTv399sMwwjklqD3cV2TBpRtVUBtrZIxW4x2ViuS7B06VJSU1N3suvYY4/l3nvv7VZZXWmrqanhjDPOIDExkYKCAl577bUeH5uUlLTDKzo6mmuvvbZLG3uid+PGjRQVdbrqoUu6o/uxxx5j9OjRxMXFcdFFF+20/6ijjsLpdLbpHjZsWI/taqX1ibGK2I5JI36/32wTdmDlypWMHDlyp0fhK1euZNSoUd0qqyttV199NQ6Hg4qKCl599VWuvPJKVq9e3aNjGxsb214VFRXEx8dzzjnndGljT/Q6nc5eDRB3R3dOTg7Tp0/nkksu6bS8xx57rE372rVre2xXK42Njb0uw6rYjkkjcXFxZpuwAytWrNjpC7lp0yZqamq67ZjCadu+fTvvvPMOM2bMICkpiTFjxjBu3DheeeWVXh0LMGfOHDIzMzn88MO7tLEnenNycsjLy+uy7FB0V8uZZ57J6aef3qdjWgUFBX12rb7Gdkwacbvdmo99f0Uph838jME3fchhMz/j/RX6L7ZcuXIl++23Y+TgFStWkJeXF3bs4dRTTyUtLW2HV3p6OmlpaZx66qk7Hb9u3Tqio6MZOnRo27Z99903ZMuhO8cCvPTSS0yaNEnTBMie6m2tt1C6W1+91a2Vm2++mYyMDA477DCWLFnS43Ja6U031eqoOzKoM1of4b6/opSb3/2VZm+we1Ra28zN7/4KwOmj9IlZ5/V6Wb169U4thZ9++qnty/v5558zePBgBg0atMMxH3zwwU7lud3uTltNjY2NpKam7rAtNTWVhoaGXh1bXFzMF198wXPPPRfyuu3RojcUGzZsIBAIsPvuu4fUHY7uaNHCPffcw5577onD4eCNN97gtNNOY+XKley66649Kg+CXVVVsR2TRv770TrWlHfdp19RXIvHH9hhW7PXzz/m/MLrPxSHPXfPnBRuP22vLq+xZs0aAPbYY48dtv/www9t3aLnn3+eO++8s8uygLbgcKFISkraaU1WfX19yDhA3Tn25ZdfZsyYMQwePLhL+7ToDUVCQgKBQKDT/eHojhYtHHTQQW3vJ0+ezOuvv86CBQs0Dfx3RlpaWo/PtTp2V04jfo03eEen1NX2ntA6aNz+UXhVVRVffPEF48aNY968ecyfP5+LL76Yl19+eYdzTz755J2ejKWlpZGUlMTJJ5+807WGDh2Kz+fj999/b9v2888/s9deOzvQ7hz78ssvM3nyZF30Arzyyiscc8wxjB49ms8//xyAAQMGtIUGCaW79dVb3T1BCNHrJ73l5eW62GJJpJR/6tf+++8vQ/Hbb7/t8Nnr9YY8riOH/nexLPjnBzu9Dv3vYk3na2HLli3S6XTKp59+WjY1Ncm1a9fK4447Tp511llSSin9fr886qijNJfXlbbzzjtPTpgwQTY2NsqvvvpKpqSkyFWrVvX42K+//lomJCTI+vr6HbZPnjxZTp48udt6f/31V3n++efLQCAg6+rq5NixYzVr00u31+uVzc3N8qabbpITJ06Uzc3Nbdfetm2b/Oijj9q2zZ49WyYkJMjCwsJOr93x/gtFdXV1z4QZwMNvXS/vuGai/N9Fb8nHpiyST1/8tlzx5IdhzwGWyU6+l3aLSSNapwvceOIw4mN37BrFx0Zz44m9n7fSSnZ2Nm+99RaPPPII6enpnHTSSey33368+OKLAKxfv75b82S60vbEE0/Q3NxMZmYm559/Pk8++WRby+Hkk0/m7rvv1nRsKy+99BJnnnnmTt2izZs3c9hhh3Vb73vvvcdvv/3G0Ucfzbhx49rGhtavX8/GjRs1/x96o/uuu+4iPj6emTNnMnv2bOLj47nrrruA4BjZ9OnT2WWXXcjIyODRRx/l/fff7/VcJqtMF3jk7b9T8dVWMl1/wRvXH4TA40jn2+WClU8t6LqAENhhTzoJe7JmzZodxjS6s0r9/RWl3PfxWspqm8lJi+fGE4fpNvCt6frvv8+mTZuYOnWqpuOtsALf4/Gw77778ssvvxAbG9utc2+44QbOOOMMxowZAwQDqMXExFBeXo7P5+vxlAEz6Xj/haKwsNAS0QWOfXYvzvv6jqBT6oDTW8ulz50Z8jw77IkOdGce0+mjcvvUEXVk6NChTJ8+nU2bNvHQQw91ebwV5mg5HI62Qe7ucsUVV3DJJZcQGxtLbGwsL7/8MtnZ2WRnZ/d48DsSsMo8pq0xAq8jdJhfV0xqyO1dYTsmjbjdbtNbFVrZc889WbVqlebjI0lbKIYOHcpXX30Vcl+kawuHVeIx7eKTxHpqQreYfHU9KtMeY9KIqjGLQF1tv//+O5s3b+76wAjFKmGDx6eeSFHKPKL8O0apiPK72f+Ans21UvOONIBwc30iHVW1paSkRHTeta5ISkoy2wQArjvnAbLG7MK26PnBDVLi8NRwyP6SkX8d26My7a6cRlqzk6qIqtqysrJobm422wzDqKystEwKp8tPvp2bl08mtRZOPMnJbmec3avy7BaTRrr7pCiSsLVFJtnZ2Wab0MbW8k3Eu9IASBkYMut3t7Adk0asFvZET1TVtm7dOjZs2GC2GYZRW1trtgltbK5YS0qrYxrce4epXvvdIFR+7KyqtrS0NLxer9lmGIaVAsWV12wkwZNKdFQzzv4pvS7PdkwascJcH6NQVVtmZqayThesM48JoLp+M47ArsRJfWaj2105jXQnHlOkYWuLTKwUj6nWVUm0TMUZrc//23ZMGlH1kTqoq23t2rWW+vLqjVWmCwA0eGsgOpVEnUJE2Y5JI1ZLM61nMoKutOkZlH/ixIkMGDCAlJQUhg4dyrPPPqvJxp7oTU9Pp1+/fprKD4XVkxFYKVBco68RX0wqicn6jA7ZjkkjPp/PbBN2QM9kBF1p0zMo/80338ymTZuor69n3rx5TJ8+neXLl3dpY0/07rLLLqSk9Hwg1urJCKqqqnpdhm5sj0ZGRZOYpo+ztB2TRqyWjlnPZAThtOkdlH+vvfZqG2wXQiCE0PRIv6d6e1pvkZCMICcnp8+u1RWx24PdyqRMfWba245JI91qMf3yFjw4Au5IC/795S3d7dEzGUH//v37NCj/VVddRUJCAsOHD2fAgAGMHdv1soWe6F27di3r1q0D1ExGYJUWkwz4iW0KRhFIzknTpUx7uoBGND92/uUtmH8deFuWQtRtDn4G2OdcXWzROxlBuHhMegflh2AAtkcffZRvv/2WJUuWdDldoafJCPr37982j0nFZARWSe1eV12Bo3VyZUGmLmXajkkjziV3QLmGUCIlP4K/wyNTbzPMvQaWvxT+3Oy94eSZXV5C72QE4RyD3kH5W4mOjmbMmDHMnj2bJ598kuuuu67TY3uajCAjI0PpZARWmce0pXw9CZ40EAFSCvRZJmN35TTi92lcttHRKXW1vQfonYwgJSXFtKD8Pp+vyzGmniYjCAQCbbOjVUxGYJWpECWVv+P0phLrbSQ6XqfJup0FA/+zvLQmI3C73Z0GVd+BB/aS8vaUnV8P7KXtfA3onYygK216BeWvqKiQr7/+umxoaJA+n09+9NFHMiEhQb7//vtSSv2TERQWFmoK6m+0bqOSEZSWlvZMmM68+s50+dSkR+QLF73SrfMIk4zAdMdg9kurY/J4POH/y638/KaUd2Xt6JTuygpu15F58+bJESNGSKfTKQcPHiz/8Y9/yIaGBimllGvXrpVXXHGF5rK60lZdXS3Hjx8vExIS5MCBA+Wrr77atu+kk06S//nPf9o+33777RLY4XX77bdLKaWsrKyURxxxhExNTZXJyclyxIgR8umnn24795hjjtnhs1a9d955p9x3333lkUceKY888kh5/vnnSymlrKqqkhUVFZr/D0bqHj16tExKSpKpqanyoIMOkp988knYa2txTL3RpiePvnSp/N/kl+Trlz3frfNsx6SDY2pqagr/X27Pz2+2tJxSg391dkpd8d5778kHH3xQ8/Hd0mYQbrdbDh8+XPsPQDv+/ve/y6VLl7Z9bp+yyQraeoIWx7RmzZo+sKRr7nr6NPnUJe/Lede90K3zwjkme/BbI92aD7PPubo9gesJ3U1GYIU5WkYkI/D7/UoGwGslN9e8hBftaXTVkxabTHKqfkH51K01nfF6vRGzpqy7yQgiSVsoOktGsH79egKBQJdpkCKViooKS4QOjqoLPkNL6q9fDHLbMdkoyy677KJ0PCarhHSJ3R50SMlZvY/D1IrtmDSicohWVbWlp6crG50TsEwiz5jmoENKztNvOY49j0kjVpllawSqavP5fJaK8qg3xcXFZpuA191MrDsNgJRBWbqVazsmjag8iKqqtg0bNljiy2sUaWlpZptA9dbNOLxpiICXhJzO12h2FzXvSBsbgqF1rRauRjVKtqwlzpeKQ9TpmjjVbjFpROUbXFVt/fr1s1SUR72xQpaUsqoNxMg0YmWTruVGhGMSQtwnhCgUQvwihHhPCJHWyXGbhBC/CiFWCiGW6WmDFeb6GIWq2rxer7LpzwHy8/PNNoHq+s0QlUJcrL7jlJFSa58CI6SU+wDrgJvDHHu0lHKklHK0ngao/NhZVW0bN25k48aNZpthGCUlJWabwLamcvwxaSQm6TsqFBGOSUr5iZSytb/xHWCN56Q2liYrKyts0LxIxwqtwab6egLRcaSm6ze5EiJz8PsS4M1O9kngEyGEBP4npXw61EFCiCnAFAjOBSksLCQnJ4eqqio8Hg8FBQV4vV48Hg9CCHw+HzExMW2PnmNjY/F4PG1Ps3w+Hw6Ho63l0X6/lBK/309cXBxut5uoqChiYmLC7o+Ojsbr9RIbG4vf7ycQCGjeHx0d3Wazw+HA5/OF3e/1egkEAvj9fqU0ASQnJ+N2u9s+R5Imv99PZWUlNTU15ObmUlFRQSAQIC8vj+LiYtLS0oiOjqawsJD8/HxKSkqIiooiKyuL0tJS0tPT8fl81NfXU1BQQFFREQ6Hg4yMDMrKysjIyMDlctHY2Ni23+l0kpaWRnl5OZmZmTQ2NtLU1NS2PyEhgaSkJCorK8nOzqa2thZZE5zk6Y/xU1JSgtPppKqqaqfvU1FRESkpKcTExLRpCktni+j6+gUsAlaFeI1vd8ytwHuA6KSMnJa/mcDPwBFdXdeQRbx9wJdffilTUlJkIBDYYfsxxxwj77nnnm6V1ZW26upqefrpp8uEhASZn5+/wyr7Vlwul7zkkktkfn6+TEpKkiNHjpQLFizodjmh6KlWj8cj6+vrNV2jI92x9dFHH5X777+/dDgcIcO2HHnkkTIuLk4mJibKxMREOXTo0C6vHymLeO+edpZ87IrFcv17X3f7XMIs4jW/LdiClPI4KeWIEK+5AEKIycCpwAUtokKVUdbyt5KgAztQL/usNtdHzywpXWnTki3E5/MxcOBAvvjiC+rq6pgxYwbnnnsumzZt6lY5oeip1o0bN1JaWtpl+aGweoYUCM5sN5tYV3CtXkq+vl1myzimcAghTgL+CYyTMvRzSSFEohAiufU9cALBFpcudOILQ/Lhxg85Yc4J7PPSPpww5wQ+3PihXma0oWeWlHDatGYLSUxM5I477mDQoEFERUVx6qmnMnjw4LbUTN3NOqKH1uzsbHbZZZcuy++p5lbMyJACFpjmISXR7uBylNQhA3QtOiIcE/AYkAx82jIV4CkAIUSOEGJByzFZwFdCiJ+BH4APpZQf6WWA1jVXH278kDu+uYMt27cgkWzZvoU7vrlDd+ekZ5aUzMxM3bOFVFRUsG7durZQtL3JOtJTrampqSQkJHSqu7MsKZGQIQXYKSZ5X9PUWEuML40Y33YcafpGObBW/6QTpJS7dbK9DBjb8n4jsK9RNjz8y8Os3dZ1E/yXrb/gCew4p8Pld3Hb17cxZ92csOcOTx/OPw/8Z5fX0DtLSiAQ6PQJT0+yhXi9Xi644AImT57M8OHDe1xOa1k9yZACtD28gO5lSYmEDClgfjKCLVs2EBNIJSbQ84w5nREpLSbT0dpi6uiUutreE8JlDTnggAOAYJYUrd1Pt7vzRAndzRYSCAS48MILcTgcPPbYYz0upxUtWjvjjz/+4I8//gh7TCiMyJCSnJxMXFwckydP5rDDDmPBggVdn9gFZicjKClfSzRpxNCoe9kR0WKyAtP2n9Zl/jOAE+acwJbtW3baPiBxAC+c9IIutoTLGnLfffe1ZUnZvHkzl1xyCZMmTWo77uSTT2bp0qUhyz388MNZuHDhDtvaZwvZfffdgc6zhUgpufTSS6moqGDBggU7hFPpTjnd0drKK6+8wgsvvEB9fT333XcfRx99NAMGDGibGtAd3T21VSt6ZEgB82fsV9T+gYzaA0d0uf6Fd/a47s/y0jpdwOfzhTyuIx9s+ECOfmW0HPHiiLbX6FdGyw82fKDpfC3onSWlK21as4VcccUV8qCDDmpLEtDdckJlSelKq5SdZ0nRoq23mqXUP0OKlNqmC9TV1XVfmE488vB98sG/viUfm7JIPnnpB/KjEy+WtfPmdasM7GQEvXdM3ZnH9MGGD+Txbx8v935xb3n828fr6pRa0TNLipZ5TF1lC9m0aZMEdpivk5iYKGfPnq2pHCk7z5ISTquUnWdJcblcPf7ydmVr+ywpemdIkdLa85geefg++cgVC+RjVyxuez1x2QK56KgLu+WcwjkmEdz/52X06NFy2bKd1/uuWbNmh3GN1hm8kcD777/Ppk2bmDp1qqbjraDN4/Gw77778ssvv3TblhtuuIEzzjiDMWPGALTN1F+7di1SyrYB+Eii4/0XiqqqKlOW3Nx73RwSPTvPoYpzVXNk8VPs/tliTeUIIZbLTta02mNMGokkB97dLClW0GZElpScnBxlFygDpkXnTPD0C7ndHZeOb8vO46s9wXZMGomk2NHdzZISSdpC0VmWlOTkZJqb9UspZDUaG/V/GqaFJse20C0mdw0xA/SZaGlPF9CIlidykYqq2lwulyVag0Zh1jym+N02EWDHKSZRfjeDi+aTef1UXa5hOyaNhJvrE+moqq2oqMj0uT5GYpa2a6+bRm2/RcEPUhLnqmbI5rcZ/fdzSD3tNF2uYXflNGKF2DdGoaq23NxcpceYnE6naddObQj+mJ1z0S5kHnIscI6u5duOSSORnKm2K1TVlpSUZP5CVwMxK0uK3+cltmkIsYEqdjnoKEOuoeZPpQGo/Murqrbm5mbTBoj7gvJyA2Zca2DZ8o8gaijOQBHCoNa27Zg0YvY8HyNRVVtxcbFpX96+IDMz05Trrlj6GYGYBPrnG3ff2I5JI5H+SD0cqmrLy8sjOzvbbDMMw6zWoO+PYNd/9LijDbuG7Zg0EggEzDbBMFTVlpiYqOxUCICmJn1zuWlBBgLENA8mvqmMjP1HGHYd2zFpROUbXFVtTU1NyrYGwZx5TL+t/pZAzK7EBTYhDHxoYjsmjVhtrs/SpUtJTU3daQLhsccey7333tutsrrSVlNTwxlnnEFiYiIFBQW89tprYY9/44032GOPPUhMTGTXXXfdKdxIV/tD0RO9mzdvZvPmzV2W3Rnd0f3YY48xevRo4uLiuOiii3baf9RRR+F0OklKSiIpKYlhw4b12K5WzJjH9P2nHyOjHKQOFF0f3Atsx6QRq8310TMZQVfauhOY/9NPP+Wf//wnL7zwAg0NDXz55ZcMGTJE8/7O6InegQMHkpOT02XZnWH1hAStYYP7kuYNAZABDjj9REOvY61vm4Xpzlyfuvnz+f2YY1mzx578fsyx1M2fr7s9eiYjCKetu4H5b7/9dm677TYOPvhgoqKiyM3N3SGHWFf7O6MnehMSEkhMTOyy7FBEQkKCpKSkPrtWK9GuApIai8kc3b17rLvYjkkjWuf61M2fz5Z/3YavrAykxFdWxpZ/3aa7c9IzGUFGRoYuyQj8fj/Lli1j69at7LbbbuTl5XHNNde0LaTtar/eerdv394WIrc7yQi6q1sreickqKys7HUZ3WHThkL8MYOIlRsNHV8Ce+a3Zmrvf4BKDc3v5p9/Rnp2jO8tXS623Dqd2rfeDntu3B7Dyb7lli6voXcygtb4RaHoTmD+iooKvF4vc+bMYenSpcTGxjJ+/Hjuuusu/vOf/3S5vzd6Q1FSUoKUkvT09G4lI+iubi0YkZCgr6dCLP3gQxCjSMwz/imu3WLSSEBqq4yOTqmr7T1B72QE4Z5cdScwf3x8PADXXnstAwYMICMjg7///e9tgfe72t8ZPU1IkJ+fz4AehuGIhIQEtbW1vTq/uzT+3kxUwMvo008y/Fp2i0kj6dOmtX2xwvH7MccGu3EdiMnJoeCVl3WxxarJCPr160deXt5OA9Ra9/dUL4RORtC+vozU3RP0SEjQ14HihGsgids3knPQFOMv1lnM3T/LS2vMb7/fH/K4jtTOmyfX7DtS/jZseNtrzb4jux2oPRx6JyPoSlt3AvP/61//kqNHj5YVFRWypqZGjhkzRk6fPl3T/lDJCLTo7SwZQUNDQ68C9puZkEBLzO/uxKHvLaXFJfKxKxbLV8+bpluZ2MkIeu+YunMT1M6bJ9cdfYz8bfgect3Rx+jqlFqxWjKCVjwej7zyyitlamqqzMrKktdee61sbm7WtL+zZARd6e0sGUFhYaGmL7geuvVOSGC1ZARvPPqUfOyKxXLO9dfpVmY4x2QnI9CYjMDj8Ziex0sr3U1GYAVtRiQjcLlceL3eHo8LmYmWZAQlJSXk5eX1iT1PXTsLmvfguAsEux09VpcywyUjsAe/NdLdcREzGTp0KM8++6xmx2QFba3JCHoS6eCKK67gpptu4uijj+aEE06gqqoKCAZSMzOYmtH0pTbpyiOlfj27Hn5Cn1zPHvzWiM/ni5jwIN1NRhBJ2kLRWTKChoYG3G63KSmO+oK+St+0tbSGQHQmiG8QnUwr0Ru7xaQRs7s6RqKqtrKysrbWk4r0ZrlNd/hyXvCJZUx230UzsB2TRlQO0aqqtkGDBmla7hKp9JXTrf6thlhPA8NOOqZPrge2Y9KMqjGLQF1tcXFxnc5oVwGPjpN2O0NKid+dS0rdWvY+Tt+EA+GwHZNGVI1ZBOpqq6+vt1y4Gj3pi3hM1WWNBKLSkFHrierDcUjbMWlE5RtcVW1btmxhi04pq61IX8Rj+vaDzwEQWT1bI9hT1G3n6oyqKY5AXW2DBw9WNgMMQEpKiuHXKF9VgdPlZuBphxt+rfbYLSaNWGGuj1Goqs3hcCj7xBEwfPwsEJB43QNIrlvLgadcZOi1OmI7Jo2o+uQK1NVWV1fX5yvw+5KamhpDy99aXI+MSsAf9Tuxzr6Nlmk7Jo2o/Murqrby8nKqq6vNNsMwjJ4K8dOi4FKtQGatodcJhe2YNGK1sQo9kxF0pa2rQPvtmThxIgMGDCAlJaVtaUwrbrebSy+9lIKCApKTkxk1atRO4UY6oyd6hwwZ0qu1ZFZPRlBRUdHrMsJR+kspidu3kHnYQYZeJxS2Y4pQ9ExG0BVaAu23cvPNN7Np0ybq6+uZN28e06dPZ/ny5UCwyzhw4EC++OIL6urqmDFjBueeey6bNm3qstye6I2Nje3VOIzVkxEYOf/M7wvgdmeSXLeWw06/yrDrdIb9VE4j3VlLtu77cr6du4HGGjdJ6XEcMn5Xhh6kbxhUPZMRdKXtzDPPBGDZsmWUlJSEPbZ9IDUhBEIINmzYwP77709iYiJ33HFH2/5TTz2VwYMHs3z58p1CAHekJ3pra2sJBAI9Wuzamoxg1apVOyUjmDlz5k7Hd+d/pBdGRhao+KMOhANf9DoSktIMu05n2C0mjWidZbvu+3I+f7WQxprg3KDGGjefv1rIuu/LdbVHz2QE/fv37zQof0+46qqrSEhIYPjw4QwYMICxY0OHyaioqGDdunWaokL2RG9FRUVbd0fFZATFxcW9LqMzVi9dAzKAO9PYAfbOsFtMGvlhbhE1pV0vYqz4ow6/b8dxEJ8nwGevrGH1VzuH3G1PxsAkDj93aNhjQP9kBF6vV9foAk888QSPPvoo3377LUuWLAk5s9zr9XLBBRcwefJkhg8fHra8niYjGDJkSNsTRxWTEaSlpfX43K4oWrmZ5IbtpB7T+f/XSOwWk850dEpdbe8JeicjMILo6GjGjBlDSUkJTz755A77AoEAF154IQ6Hg8cee6zLsnqajKA3Y0yRkIzAKDwuH253Oin1aznsrGtNscFuMWnkwPEFmpIRvHTL123duPYkpcdxxg36/Pr0ZTKC3uLz+diwYUPbZykll156KRUVFSxYsEBTS62nyQi2bduGx+MhKytLyWQEtbW1hqRw2rK+DkQ0ruh1ZPQzJzqD3WLSiNa5PoeM35UYx47/1hhHFIeM73mTvSN77703brebZ555hubmZtatW8f555/Pqaeeyj777MOpp57KqFGjWLJkyQ5OCWDhwoVtT4ZaX3V1dTQ2NnbqlHw+Hy6XC7/fj9/vx+VyhZyUWVlZyRtvvEFjYyN+v5+PP/6Y119/nWOO+f9wGVdeeSVr1qxh/vz5Ozn6iy66KOSj9q70rlq1ioULF7J48WI+++wzZs2a1WbPtm3bOtXd+gqlOzExkTPPPJPbbruN7du38/XXXzN37lwuvPDCbv+Pamtr+fjjj9u2vfrqq3z55ZeceGLv0mzn5+f36vzOWPf9RkTAS3OmiXPAOgsGbqUXcAdQCqxseY3t5LiTgLXAeuAmLWVrTUbQPqB+V6z9bot88eav5GNXLJYv3vyVXPvdFs3nakXPZARdaQsXaL99UP7Kykp5xBFHyNTUVJmcnCxHjBixQ3KBTZs2SUDGxcXJxMTEttfs2bOllPonI/B6vbKxsVHz/6EjVk9GsHbt2p4J64Jnrn1Xzj73Sfn6o1cbUn4rRHoyAiHEHUCjlHJWmGOigXXA8UAJ8CNwvpTyt3Bla01G4HK5IiZ+dHeTEVhBmxHJCMAa2nqClmQE69evZ7fddtPtmuu+L2fxSz8Q8CcS7Wsmd/ObNB/o4dy73tTtGu35syQjOBBYL6XcKKX0AG8A4/UqPJJiYnc3GYEVtBmRjKCmpobGxka9TbUMWVlZupW17vtyPn3hZwKBJBACf2wCJQV/Ie0rB29NP0+362glkga/rxFCTAKWATdIKbd12J8LbG73uQTQbS69x+PRNPhtBbqbjCCStIWis2QEW7duJRAIKJuMoLS0tMupFlr5du4GYMcfhUB0HMWDxjH8k9vgLl0uoxnLOCYhxCIg1COGW4EngRkE++4zgPuBjnP/Q8XuCNlPFUJMAaZAcPZsYWEhOTk5VFVV4fF4KCgowOv14vF4EELg8/mIiopqS8kcGxuLx+Np6y74fD4cDkfbmrP2+6WU+P1+4uLicLvdREVFERMTE3Z/dHR029wiv99PIBDQvD86OrrNZofDgc/nC7vf6/USCATw+/1KaYLgPCa32932OZI0+f1+KisrqampITc3l4qKCgKBAHl5eRQXF7fNYSosLCQ/P5+SkhKioqLIysqitLSU9PR0fD4f9fX1FBQUUFRUhMPhICMjg7KyMjIyMnC5XDQ2NlJQUBDySTKAOy6dfvXB6yQkJJCUlERlZSXZ2dnU1tbicrnayk9KSsLpdFJVVbXT96moqIiUlBRiYmLaNIUjIsaY2iOEGAR8IKUc0WH7IcAdUsoTWz7fDCCl/G+48lRMeNldbG3WQ8sYU1lZmW6ZUl66+Wsat+3snOJc1Qz/7TbG/LBGl+u0J+LHmIQQA9p9PAMI1U/5EdhdCDFYCOEAJgDz9LLB7/frVZTlUFVbdXW14TGLzKTjBNDekJiy8z0Q5XeTv2keZSfso9t1tBIRjgm4VwjxqxDiF+Bo4HoAIUSOEGIBgJTSB1wDfAysAd6SUvZ8YVMHVA3YD+pqq6qqUjpQnF7JCLYW11OxyU1y7VoEtSAlca5qcotfo3aMcU/lwmGZMaZwSClDzmqTUpYBY9t9XgDoNs9fStkWZsPtdkf0AHE4VNW2++67t40LRhJaw5kUFRX1evA7EJDMu38Rsd4o3KlvctVTS9rt7bt0TR2JlBZTn+N0Oqmurm5bNhAVpe6/SlVtrQPQkYKUEo/HQ2lpKYmJiV0er8fY2bL5q3G500ivnMP4u17udXl6EREtJjPIy8ujpKSErVu3AsFfMVW/wKpqa2xsRErZ44W3ZhATE0NqaqqmKQ69nQbRUONi2YLN9Nu2Ad9pmWTuYswSl55gO6ZOiI2NZfDgwW2fCwsLdZszYjVU1XbUUUfR1NTEDz/8YLYphlBWVtarFE5zZy1EBOJpjnuXS6Ys1tGy3mM7Jo2oOkkP1NW2ZMmStlngKtKbelv7bTF1Nalklb/Hvg/eY7kUXrZj0kgkDqJqxdYWmfRUm7vZx+cvriBpex0Nh9aw+26dx7QyC/UGFgxC5TVXqmp75plneP755802wzB6Wm8LHv0EP4kIzxtM+sdLOlulD7Zj0ohec0asiKra3nzzTV1ia1uVntRb2bpqyjY4yKz4gj1uuILoaGt2mmzHpJGioiKzTTAMVbUtWrSIJ554wmwzDKO79eb3B/jgwc+Ic9dRN2wlBxykW/AN3bGmu7QgkRjTRyu2tsiku9qWvPQVXtmffvX/44yn3jDIKn2wW0waMTIjhdmoqu2JJ55gzpw5ZpthGN2pt9rKJtZ+10R69UoyLzuaeGfXEzjNxHZMGikv1zcvnJVQVdv8+fOZN0+3ddyWQ2u9SSl57+65RAX8NO7yCceN/avBlvUeuyunkczMTLNNMAxVtS1cuFDp6AJa623Fwl9ocmXRv/pNjo+Qp5R2i0kjqj5SB1tbpKJFm6vRy/fvF5Fc/weO8bvQv58+8ZuMxnZMGmlq6joLb6SiqraHH36Yp556ymwzDENLvb078x0k8Xgd73LmRX0cH7cX2F05jag61wfU1bZ48WJlg+BB1/W2cXkx26oy6b/1Yw5+6J4+skof7BaTRlSd6wPqaps3bx7333+/2WYYRrh683n9fPq/73E2b8V7aDWDCkZ0eqwVsR2TRhISEsw2wTBsbZFJOG0fPvIBvqj+RHnfZuINT/ehVfpgOyaNJCUlmW2CYaiqbdasWbz44otmm2EYndXb1k01lKyNp1/1D+z/j78iIjDWVuRZbBKVlZVmm2AYqmr79ttvQ+abU4VQ9SYDkvfv+ZAYv4vm3X9in5HHmWBZ77EHvzWSnR0q5Z0aqKrtnXfeUToZQah6++KVz/HIXFIaZnPeU6+YYJU+2C0mjah8g9vaIpOO2rbXuljzVRMpdevIv+wIHLGRm/3GdkwasQOORR4zZ87k0UcfNdsMw+hYb2/f8TqIWFwZH3PkCReZY5RO2F05jag61wfU1bZy5co/zTymnz9dwXZXAanb5nPac8+aaJU+2C0mjag61wfU1fbGG28wY8YMs80wjNZ687r9fPvmBuKbtpB8egapKf1Ntqz32C0mjaj6SB1sbZFKq7Y5M2bjjxlIQD7J+AveNtkqfbBbTBqxA45FHjNmzODxxx832wzDcDqdFP1SRM3WXFJrvuLkGf812yTdsB2TRlROA6SqtrVr17J69WqzzTCMysqtfPzol8R6G4k6tJrcnN3MNkk37K6cRnJyIiNcRE9QVdvs2bOpr6832wzDWPvBT3ijc4l3vcj5f7NmtpOeYreYNKJqqwJsbZHItrI6yn9PJ6VuFQf/c4rlElb2FtsxacTj8ZhtgmGoqu22225j1qxZZpuhO1JK3r3zLQC8u61gzxFjTLZIf+yunEZUnesD6mrbvHkzPp/PbDN058vXPsbFriRuf48L73jGbHMMwXZMGikqKmL48OFmm2EIqmp74YUXKCwsNNsMXXFt97Dm82YS3OWkjR9q2YSVvcXuymkkJSXFbBMMw9YWObx5yzP4o5Px9/+UAw8/y2xzDMN2TBqJiVHzlwnU1XbzzTczc+ZMs83QjV8W/0ijew+S6pdw/t1PKVtvYDsmzaicBkhVbdXV1crkzPP7Anz32jriXDXscnoGiQnJytYb2GNMmsnNzTXbBMNQVdvTTz9NQ0OD2Wb0ipen3YKrehReRzrEDsDpWsTY8+4G1K03sFtMmqmoqDDbBMOwtVmTl6fdwva6w/HG9YeWeUqehMN5edotQGRr6wq7xaSRQCBgtgmGoaq2adOmsW3bNp577jmzTekRrupRBOJ2DPYWiI7DVT0q+F7RegO7xaSZvLw8s00wDFW1NTc3R/QAsdeRHna7qvUGtmPSTHFxsdkmGIaq2h5//HGuv/56s83oEX+s/A0hQwe5i/UEB71VrTewHZNm0tLSzDbBMGxt1uLH+Z/wyaPrIOBHBLw77Ivyu3H2XwFEpjatRG4718amC6ZOncr27dt55pnIWbax4IkXKVqRRYzfQ+pe3+OpEG1P5WI9NTj7r2DSrLvNNtNwbMekkdraWmXTHKmsze12m22CZt7410NUV44g3r2Z3JOaOfEvoQK/ndP2TuV6sx2TRvLz8802wTBU1fbQQw/R1NRkthldIgOSl/42i+3e/Ulo/IWRVw5l1JixXZ6nar2BPcakmZKSErNNMAxbm3m4m9w8f/kjQafU8CXH3X2iJqcE1tfWGyKixSSEeBMY1vIxDaiVUo4McdwmoAHwAz4p5Wi9bIiKwPzvWlFV29VXX01dXR2zZ88225SQbC0p4/1/LcATtzfxzfM556nbSErup/l8VesNIsQxSSnPa30vhLgfqAtz+NFSSt3DFmZlZeldpGVQVVt8fDzR0dFmmxGSNd8vZ+n/1uGLHUi8mMNFzz9KVDdtVbXeIMK6ciIYP/Rc4PW+vnZpaWlfX7LPUFXbrFmzuOqqq8w2Yye+fHsuXzxdghSJJGUs4JInn+i2UwJ16w0izDEBhwMVUsrfO9kvgU+EEMuFEFP0vHB6euhZuCpga+s73nvwKVZ/EkeM30Xmvj8z6e6He1yW1bTpiWW6ckKIRUCoZ5+3Sinntrw/n/CtpcOklGVCiEzgUyFEoZTyyxDXmgJMgeC0/sLCQnJycqiqqsLj8VBQUEBRUREpKSnExMRQU1NDbGws9fX1BAIB8vLyKC4ubpvgVltbS35+PiUlJURFRZGVlUVpaSnp6en4fD7q6+vbynQ4HGRkZFBWVkZGRgYul4vGxsa2/U6nk7S0NMrLy8nMzKSxsZGmpqa2/QkJCSQlJVFZWUl2dja1tbW4XK62/UlJSTidTqqqqrrUlJubS0VFBbW1tSQlJSmlKRAIMHPmTLZt28aTTz5pCU2fP/oqTY0HEe8qIvP4Jg4ee6Wme6+9pvb3Xk1NDTU1NRFbT2H9gZQy7AFWQQgRA5QC+0spu3wcIYS4A2iUUoaNRj969Gi5bNmyLq9fWFioZPhZUFfbzTffTHV1NU8//bSpdvh9fl66dhbN8gASGlZy0A37s+d+R/a63EivNyHE8s4eUFmmxaSB44DCzpySECIRiJJSNrS8PwG4U6+LqxqwH9TV9t///pfm5mZTbWiqb+S1qc/hdhxAQuMSTpt1KRkD9Pl/q1pvEFljTBPo0I0TQuQIIRa0fMwCvhJC/Az8AHwopfxIr4sXFRXpVZTlsLUZQ+nGP3j1b2/ijt2LBPc8Lnj2Zt2cEqhdbxHTYpJSXhRiWxkwtuX9RmBfo67vcDiMKtp0VNV28cUX09DQwJw5c/r82iuWfMmPL2/BH5NLkuN9Jj31qO5JKVWtN4ggx2Q2GRkZZptgGKpqGzhwoClr5T55+TU2fpmEEE765X7GhDseM+Q6qtYb2I5JM2VlZcqlAmpFVW133nlnn+eVe2vmQ1Rt3AOHr4bcw8o5+fJ7DbuWqvUGtmPSjMq/Tra23iOl5KVp/2V744HEN//BiAvTOfAEY4PUqVxvtmPSiMvlMtsEw1BV28SJE2lqauLdd9819Dpej5eXrr0ftziYxMafOGr6sQzaY5Sh1wR16w1sx6SZxsZGs00wDFW1DRs2jK1btxp6jW3V1cy58VU8jgNJbPyMsx69juT0TEOv2Yqq9QYRNMHSKLROsGxubiY+Pr4PLOp7bG09Y8Pq1Xw26wc8joEk+j5k4lP3EBPbd0/KIr3ewk2wjKR5TKai8pwRW1v3+XbhQhbfvxpfzABSEz7kouce7FOnBGrXm92V04jT6TTbBMNQVduECRPYvn078+fP17Xc+f97htIfs4nCQeaQbznr5gd1LV8rqtYb2I5JMypnpFBV28iRI3VfkjL73zOpLx2Fw1vNkOMbOWbiv3UtvzuoWm9gOybNlJeXK3sjqKrtpptu0m0ek5SS56fOwOUeQ0LT74y+YhB7H/4XXcruKarWG9iOSTOZmX3zpMUMbG3hcbncvHLtg3iix5DYsIwTZ5zOgCFDdbCud6hcb/bgt0ZUfjSrqrazzjqLCy64oFdlVJSX8cqVT+OJPpCk7Ys476krLOGUQN16A7vFpJlISAPUU1TVdsghh1BZWdnj81ct/5FvHluN17EHyYEPmfjivT0KgWsUqtYbhJnHJIT4oQflSWCSlHJtr6zqQ+x5TLa2UHz23jv8Pl8SiE6kX7/PmXCvcWveekqk11tP5zGNBiqB1Rpfv7Wck6ib5RZC5TkjtrYdmfPII6z7MA5BNAP3+NmSTgnUrreuunJ3Sik1tZxaQt9O6r1J1iQhIcFsEwxDVW3jxo2jubmZTz/9VPM5L06/k6bKg4nzbmWPcYJDz7rJQAt7h6r1BuEd08XABq0FSSl9QoiLgT96bZUFSUpKMtsEw1BV27HHHqt5HMbv8/H81Bl4fEeS0LSWMVP3YffRhxlsYe9Qtd4gjGOSUr7U3cJ6ck6kUFlZqWy6HFW1/e1vf9M0j6mhsYHXpz6ON+ZIkhp+4NR7/kL/vPw+sLB3qFpvYD+V00x2dqjMUmrwZ9ZW9McGPr7zI7xxB5LcvJjzn72F2PjIGCZVud40OyYhxLnAGUAusNMiHSnlgTraZTlqa2uVnWWrqraTTz6Z5uZmlixZEnL/sq8+56fnNuNzDCVNLOQvL96re1xuI1G13kCjYxJCzAT+AfwIrAc8RhplRVQOyqWqttNOO43y8vKQ+xa89hKbFyUho/uTmbmUs2fc18fW9R5V6w00xmMSQlQCD0op/2u8SX2LPY/pz6ft1fvuoX7dCGJ8LgYdWMrxV15nknW9I9LrTY94TF5guX4mRR4qzxn5M2l75h//ou73/XB4ahh9TnTEOiVQu960jjE9DFwmhPhU/klDXqr8aFZVbccddxxut5ulS5fi9fp47ro78cujSWhaw/E3HUreCOPjchuJqvUGGh2TlPJeIcQsoFAI8QVQu/Mh8p96G2clVA7Kpaq28847j8bGRqprqphz47P4Y48iqfF7znroUpIyIn9lvqr1BtoHvy8ApgIBIImdB78loLRjqqqqUjZdjoraXp52C7J6FA7HEN7853Jk7IGkuBbzl+duIzouzmzzdEHFemtFa1duJvAm8FcpZYOB9liWnJwcs00wDNW0vTztFrbXHU6gxQHJ6FhEwIfsL5VxSqBevbVH6+B3CvD8n9UpQfDXSVVU0+aqHkUgOo6H5v2dh+b9HQAZFYOrJrLHlDqiWr21R2uL6R3gaGCxgbZYGo9H3albKmnbtq0GryO4TOPgYSfssK91uyqoVG8d0eqYPgZmCiGygc/YefAbKeUCHe2yHAUFBWabYBiqaFuy4B3Wv1UNzt0AOHjYSTvsj/XUmGGWYahSb6HQ2pV7HSgALgFmAx90eOmbH8eCqDxnJNK1SSl5YfrNFL4Tizd2IPHNXxHld+P3+/D7fQBE+d04+68w2VJ9ifR6C4fWFtNgQ62IAFJSUsw2wTAiWVtZWREfTn8Wj+N4nL4S9jkjmQNOv42Xp93Cf2bPJRAVw40nTsfZfwWTZt1ttrm6Esn11hVa5zGp65o1EhOjbiCGSNW28I3nKF0o8cQfTaL7e85+YApJ/fsDMGnW3USN3JP6+nqmXHUOcI65xhpApNabFjrtygkhUkQ3l1r35JxIoaZGrfGJ9kSaNr/Pz3M3XE/Rp5n4YgeQlfUDF71wc5tTamXixIkcc8wxJllpPJFWb90hnMvdBhwCaA2tG91yzgHAT703zVrk5uaabYJhRJK2jb+v4vO73sYVfxpx/mIOmpTP3seEDn/b1NREv379+tjCviOS6q27hHNMAjhUCKF1aqnSOeoqKipITk422wxDiBRt7z1zPzVL03ElHEmC/3v+8ti1xIVZLzZ27Fiam5v5/vvv+9DKviNS6q0ndNVJfaBPrIgAAoGA2SYYhtW1uT1uZl//d7zuseCQ5O26mvE33tzleVdeeSVbtmzpAwvNwer11hvCOaaePokr6+F5liYvL89sEwzDytpW/fQVPz6wCFfSWTj8mzj66pHsNvpUTeeed955SieFtHK99ZZwyQj+9E/i2lNcXMzw4cPNNsMQrKrt7fv/RcPPu9GcdASJUT9ywVPXExvn0Hx+XV0d69at44ADDjDQSvOwar3pgbrPG3VG1djKYD1tjY11vDX1n3gYBw4fQ0YWcfJfux+8Yvz48Xg8Hr755hsDrDQfq9WbntiOycZS/PD5fNY8vZLm1HNxeDZy0j+OZODw3XtU1nXXXce2bdt0ttCmL1D6SZqe1NbWmm2CYVhF26t3TmXVC7U0ph5GUvwKLnn64h47JYAzzzyTww6zdtLK3mCVejMCu8Wkkfx86ydA7Clma6up3sK8v99Os+MMhMPDnmNqOXriDb0ut6qqSuk02mbXm5HYLSaNlJSUmG2CYZip7cu5L/Lh1S+wPX4C0bKEs247kKMnnqlL2WeffTbnnKPeUpRWVL4nu5PwMgc4Fchj54SXysf8jopS14eboS3g9/PqTX/FW3E0zakHkpq2mr/cfQ1RUfqtaLrhhhsoK1Ny9gqg9j2pNeb3GQRDn0QDlfwJY35nZWWZbYJh9LW20k1r+eyW+2lMPAvhaGbUyZJDT79W9+ucdtppNDSoG3RV5XtSq8u9G/gEyJJS5kopB3d4DemtIUKIc4QQq4UQASHE6A77bhZCrBdCrBVCnNjJ+elCiE+FEL+3/NV1kVRpaamexVmKvtT2ycsP8Pk/36c+ZQLRYjMT7j6aQ08/3pBrlZeXs3LlSkPKtgIq35Nau3IDgWullEYuZ14FnAn8r/1GIcSewARgLyAHWCSEGCql9Hc4/yZgsZRyphDippbPurXi0tPVCsvanr7Q5vP6eO36i/E2nIIrLYP0nPVMmH45QseuW0cmTJiA1+vl66+/NuwaZqLyPanVMX0DDAMWGWWIlHINQIioKeOBN6SUbuAPIcR64EDg2xDHHdXy/iVgCTo6Jp/Pp1dRlsNobetX/8D3dzxPY7+JRMVu59Bzkxh17BRDrwlw0003UV1dbfh1zELle7JTxySEaP+c9e/Aq0KIRuBTQsf8NmpRUi7wXbvPJS3bOpIlpdzSYssWIUSnGQ2FEFOAKRBcb1RYWEhOTg5VVVV4PB4KCgooKioiJSWFmJgYampq8Hg8NDU1EQgEyMvLo7i4uG3mbW1tLfn5+ZSUlBAVFUVWVhalpaWkp6fj8/mor69vK9PhcJCRkUFZWRkZGRm4XC4aGxvb9judTtLS0igvLyczM5PGxkaampra9ickJJCUlERlZSXZ2dnU1tbicrna9iclJeF0OqmqqupSU25uLhUVFVRUVJCWlmaIpvceuhW5Iofa/ucSw3pOuflUAtE+CgsLDdUUCAQ44ogj+OmnnygvL1einjree2VlZdTX10espnCIzjJ+CyECBAe12za1/A15gpQyOuyVgmUuArJD7LpVSjm35ZglwDQp5bKWz48D30opZ7d8fg5YIKV8p0PZtVLKtHaft0kpuxxnGj16tFy2bFlXh9Hc3Ex8fHyXx0UiRmhzbW9kzrWX4PafjcvZj6zdyjnrhomGdt06snnzZlwuF7vv3vNJmlYm0u9JIcRyKeXoUPvCdeUuoRMn1FOklMf14LQSgmNcreQROoJBhRBiQEtraQDBp4e6UVRUpOyCSb21/frNR6yeNY/6jMsQ0fUcc3EOexxszAB3OC688EKampr44QdNsQ4jDpXvyXDRBV7sQzvCMQ94TQjxAMHB790JHVVzHjCZYNbgycBcPY1wOLSvao809NT27n+uo/nnIdRmno3DsZGJ/7mQ+GRzst9Onz5d6XhMKt+TmqYLCCE2CiH27WTfCCHExt4aIoQ4QwhRQjCc74dCiI8BpJSrgbeA34CPgKtbn8gJIZ5tN7VgJnC8EOJ34PiWz7qhao540Edbw7atvHbBWdT8fgS1/fZk4N61XPbwpaY5JYDjjjuO8ePHm3Z9o1H5ntT6VG4Q0NkdlkCwe9UrpJTvAe91su8/wH9CbL+s3ftq4Nje2tEZZWVlyqbL6a22Hxa8StHT31Gb9VeErOXka3ZnyN7mZ/zauHEjGzZs4Pjj+74b2ReofE+GeyqXAqS125QthOi4atBJcI6RujO9WlD516mn2qSUzLnpIrx/jGZb9hnEJRZz4V0TiYu3xtrwSy65BK/Xq6xjUvmeDHcHXQ/cTnAAXNJJa4bg07reLwW3OC6Xy2wTDKMn2ipLNvD59TdSnzQRT1oiux3s4YTJk0PNQzONf//731RW6voMxFKofE+Gc0yvAcsIOp55wDRgbYdjPMBaKWWxMeZZh8bGRrNNMIzuavvi9Ueoev0PqgZcRRQ1jL9hX/KGhpoFYi5HHnkkhYWFZpthGCrfk+Geyv0O/A4ghDga+ElKqe6KyC4oKCgw2wTD0KrN7/UyZ+r5+KqOpjbnNBL6lXHB7RNwOK3RdevI2rVrcbvdZpthGCrfk5qeykkpv5BSNgghhgkhJgohbmz5q+YkihAUFambm0GLtuLCFbx39gTqmydSm7Yrex4Tw0V3X2BZpwRwxRVXcPnll5tthmGofE9qDXuSDDwLnEXQmTUCSUBACPEucJmUst4wKy2A09kxBJU6dKXt0ydvp2HBdipyriQqqpqzbzqSrALrZ7i9++67lZ7HpPI9qfXn7kngBGAS8K6U0iWEcBJ0VI8BTwATjTHRGqickaIzbd7mJt676hy8TadRmzuU5KytnH/L2cTGdbn6yBIceuihSsfFVvme1BqPaTxwo5TyNSmlC0BK6ZJSvgr8o2W/0rQuBFWRUNrWff8ZH559IbViCrVpBYw6LZVJ/z4vYpwSwKpVq/jyyy/NNsMwVL4ntbaYGoHO2sRlwHZ9zLEumZmdBiuIeDpq+/C/1+H5Op6ygVcTFV3FhFuOpn9u5E3ku+aaa/D5fIwbN85sUwxB5XtSq2N6HJgmhPhMStncurElNMo0gl05pWlsbFQ2MFertu211Sz86/m4xDnU5e1Kv4J6zr3hLGIckdNKas99992ndKtC5XtSq2NKJbh4drMQ4lOCK/czCa5JawaWCSHubTlWycQETU1GhZsyn6amJn755G2KHnib6txr8MXEcMg52ex37DFmm9YrDjjgAKXnMal8T2p1TGcD3pbXwe22N7Tb34qSiQlUnTMipaTwhbuQvxVQNuivRDuqueDW40jLivx8bCtXrlR6drSq9yRodExSSvNXZJqMirFvtpUV8cVVF9OQMJH63EFkDXVx+rVnEhMbmV23jkydOtWOxxShWHd2nMVQIaPrW9PPw/mjg8oB43DHpRPja8SfcSMyOsBRFw5mr8PU+v156KGHlJ7HpMI92RmaM+YJIfYRQrwphNgghHALIfZr2f4fIcTJxploDZKSksw2oVe8Nf080r5yUJr/F9zO/iAEvthkZHQso07ZVTmnBDBy5EgOOuggs80wjEi/J8OhNVDcycBygvG6XwZi2+12A/pnK7QYkb5KPeeTXygeNI5AdIewWiKKNUsrzDHKYH788UcWL15sthmGEen3ZDi0tpj+C7wopTySnQO2rQRG6miTJcnOtt7q+e6QXg/uuNCPll0NHVP0qcGNN97Iww8/bLYZhhHp92Q4tI4xDSc4Xwl2TlBQD6g5maIdtbW1EbsEYOOyL6hLzidYdTvHS4pPUXOo8bHHHlM6W20k35NdofWOrAQ6SwO+F6B8PKZIfey8+tN3KbvrJdaP/BvR3kZktJNAdPsg9l4GH5pmlnmGMmLECGJi1HS6ELn3pBa0duXeAO4UQoxpt00KIYYSnLP0qu6WWYxInDPy/VtPUnr32/w+7Cqik32kZs0lt/hV4lzVICVRUQ0cf/G+HHziHmabagjffPON0k/lIvGe1IrWn5N/AXsCXwCtc/znEhwM/wS4W3/TrEWkzRn54n8zaH5rExuGXkp8uou/TD8ZZ2LoNWOFhYURpU0rt9xyiz2PKULROsHSDZwqhDiWYCaSDKAGWCyl/NRA+yxDJD2a/ejua/F/Idi0+wWkDPAw4aaTw0YFiCRt3eF///uf0mvlVK036OYESynlYkDd569hiJSgXPP+eQGsGsTmwceSuTucOfUEoqPD99gjRVt3GTZsGP379zfbDMNQtd5Ag2MSwbQXxxNcI5fVsrkC+IZgi0nXNOJWpaqqytLpcqSUzL1yHP7yQygfeDCD9nMy9rJDEFFdZy2xurae8sUXX1BUVMSkSZPMNsUQVK036MIxCSFGAW8CuwJ+oIrg8+b+LeeuE0JMkFKuNNhO08nJyTHbhE6RgQBzJx2H2zWOqux92OvodI48d1/NqZSsrK033H777fj9fmUdk6r1BuETXmYBHxMMEDcW+FxK6WnZFwccA9wDfCyE2FtKqe40VIK/TlbMeup3ufjgguNpjL2Quv5DOOiMPEafOLRbZVhVW295/vnnKS5WdyaLqvUG4acLXEsw1tLhUsqPW50SBAfDpZQLgSNajrnGWDPNx+PxdH1QH9NcX8uHZx9Hbfzl1KUO4uiLhnXbKYE1tenBkCFDlJ4drWq9QXjHdALwRLjsJ1LKWoKJCk7S2S7LYbU5I7UVm/nsnPFs7X8tjUmZnHLNKPY8ZGCPyrKaNr1YtGgRGzZsMNsMw1C13iC8Y9oN+ElDGctbjlUaK+Xw2rJ+Fd9ecDFlOVNxJSRx1j8OYdDeu/S4PCtp05O77rqLf//732abYRiq1huEH/xOBeo0lNEAqNnRbYdV+vIbln/J7/+4h6LBUwnESc6/9SjScxJ7VaZVtOnNK6+8ovQ8JlXrDcI7JsHOC3bDHas0Vlhz9cuidyn772v8MeRaRLybibefQEr/+F6XawVtRjBw4EDi4uK6PjBCUbXeoOu1ch8LISrDvYCFfWGo2dTU1Jh6/e/efoqyme+zYbcriElxM+muk3VxSmC+NqP46KOPmDt3rtlmGIaq9QbhW0zqds57QG5urmnX/uzpO3HN2cIfu11EYv/t/GX6KTji9fu1NFObkcycORO/38/ll19utimGoGq9QRjHJKW0HVM7KioqSE5O7vPrfnjPdQSWxFE85Bz65bk5959jdU8WYJY2o3njjTf4448/zDbDMFStN7CTEWgmEAj0+TXfv3kigdW7s6XgcAbsAadfcyJRXax76wlmaOsLsrOzqa/vdLZLxKNqvYHtmDSTl5fXZ9eSUvLu1eMJlB9BZc5+7HZQPCdcdLDmJSbdpS+19SXz58/H7XYzdGj3J51GAqrWG3QjS8qfnb5a2hAIBHjnouNxV51IZeZ+jDwxgxMvPsQwpwR9p62vuf/++7n33nu7PjBCUbXewG4xaaYvYiv73C7mTjyR5tjJ1PcbyGHnFTDy6F0Nv66qcaPnzJlDRYWaGWBA3XoD2zFZhqb6Wj6+cDz1KVfQlNCf4y/fi6Gj1V093hdkZGTg8/nMNsOmB9hdOY3U1tYaVnZNeQmfTjiT6n7X0pSYxvjrR/epUzJSm5m8++67vPXWW2abYRiq1hvYjkkz+fn5hpRbuv5Xvpl8GRXZf8MTH8e5txxO3vC+jbpolDazeeSRR5gzZ47ZZhiGqvUGtmPSTElJie5l/r78C1Ze+S9KBl6HPx7+csex7JLf9/NSjNBmBebOncuDDz5othmGoWq9gT3GpJmoKH19+MrP3qXsv29TNPhKouObuPDfY0lMM2ddl97arEJqaiqpqalmm2EYqtYb2I5JM1lZWV0fpJGv33mS2v8tY9OQS3GmNHDBHafiTIzVrfzuoqc2K/Hmm2/S3NzMbrupGZVH1XoDuyunGb1STX/yzAxqnilk05ALSM5sYtJd40x1SqCfNqvx5JNP8sQTT5hthmGoWm9gIcckhDhHCLFaCBEQQoxut/14IcRyIcSvLX+P6eT8O4QQpUKIlS2vsXral56e3usy5t33N5rmuigeNJ6MAjcX3HFq2HxvfYUe2qzIggULlH4qp2q9gbW6cquAM4H/ddheBZwmpSwTQowgmCChs2XVD0opZxlhXG/nw7xz64XI1XtQkXcweSME4646SVNqpb5A1bk+CQkJOBwOs80wDFXrDSzUYpJSrpFSrg2xfYWUsqzl42rA2ZKlpU/pzWLQt649He/a0VRkH8zwwxIZd/VRlnFK0DttVmb27NnMnj3bbDMMQ9V6Aws5Jo2cBaxoSVkeimuEEL8IIZ4XQvTT88I9Cfzu9/t58+IT8GwdS3X/vRh96i4ce+FBhq576wmqBrV/9tln+eCDD8w2wzBUrTfo466cEGIRECqfzq1SyrChBoUQexHMY3dCJ4c8CcwgGA54BnA/cEknZU0BpkBwhXZhYSE5OTlUVVXh8XgoKCigqKiIlJQUYmJiqKmpwePxkJCQQCAQIC8vj+Li4ra1SrW1teTn51NSUkJUVBRZWVkU/bGR1f/5G67YS2hMzebQ8wYSnyPZuHEjGRkZlJWVkZGRgcvlorGxse2aTqeTtLQ0ysvLyczMpLGxkaamprb9CQkJJCUlUVlZSXZ2NrW1tbhcrrb9SUlJOJ1OqqqqutSUm5tLRUUFFRUVjBo1qktNpaWlpKen4/P5qK+vbyvT4XBYTlMgEGDu3LmsWLGiLe63Cpra33slJSUkJSVFrKaw33erZfgWQiwBpkkpl7Xblgd8BlwspfxaQxmDgA+klCO6Onb06NFy2bJlXR3Gxo0bGTJkSJfHATQ11LJg0pnUp1xBszOVk68ayZB9rZvfrDvaIg1bm3URQiyXUo4Otc/yXTkhRBrwIXBzOKckhBjQ7uMZBAfTdUNrjviaihI+Pv9ctqVdiys+kTNuPMjSTgm0a4s0XnzxRT7++GOzzTAMVesNLOSYhBBnCCFKgEOAD4UQrXfUNQTz1v2r3VSAzJZznm03teDelikFvwBHA9fraV9ZWVmXx2zesIqlF/2Viuzr8DqjmHDb0eTsputQlyFo0RaJvPjii7zwwgtmm2EYqtYbWLAr19do7cpVVVWF/YUq/OkL1t3yKKUDL0M4mvjLnfplMTGarrRFMrY26xKuK2eleUyWxuVydbrvp8/eofTeuWwuuILY+AYm3nkqCSmRM38mnLZIx9YWmVimK2d1GhsbQ27/8t2nKLlvEcWDLiI+tYGL/jsuopwSdK4t0nnmmWd4/vnnzTbDMFStN7Adk2ZCzRlZ+Oyd1Dy3gc0F55Ca1cCk/4zXNd9bX6HqfJg333yTJUuWmG2GYahab2B35TRTVFTE8OHDeXnaLbiqR+F1pBPtH41/oJPMIS7OuuE0Q1Ir9QWt2lRj0aJFFBYWmm2GYahab2C3mDTjdDp5edotbK87HG9cfxACf4wTEfDTVPFVxDolCGpTFVtbZBK536Y+Ji0tDVf1KALROy7Tk1HRuGtGmWSVPqiabeOJJ55QOrSuqvUGtmPSTHl5OV5H6DATnW2PFFqXbKjG/PnzmTdvntlmGIaq9Qb2GJNmMjMzifX8GuzGdSDWU2OCRfqRmZlptgmGsHDhQmpqIrtuwqFqvYHdYtJMY2Mjzv4riPLvGNggyu/G2X+FSVbpg8qPnW1tkYntmDTS1NTEpFl3k1T3JgT8ICWx7moSU5cyadbdZpvXK5qamsw2wRAefvhhnnrqKbPNMAxV6w3srpxmWueMxMYWIgiw7/EFHHb2scA55hqmA6rOh1m8eDF+v99sMwxD1XoDu8WkmaKiIgCiXAIZFUtCqjqPalu1qca8efO4//77zTbDMFStN7Adk2YSEhIAiHYHF+ZG2rKTcLRqUxFbW2RiOyaNJCUlARDtD/5VyTG1alONWbNm8eKLL5pthmGoWm9gOybNVFZWAhAdCKbwVskxtWpTjW+//ZavvvrKbDMMQ9V6A3vwWzPZ2cEolFEyBVDLMbVqU4133nmH2tpas80wDFXrDewWk2Zqa2uprd2KjE4G6Tc9e66eqPzltbVFJrZj0ojL5WLz+lV4HClERTVZKi9cb1E14NjMmTN59NFHzTbDMFStN7Adk2YKCgrYuuk3PI5komM8ZpujK6rOh1m5ciXFxcVmm2EYqtYb2I5JM0VFRdRv+QOPI4VYZ8Bsc3RF1fkwb7zxBjNmzDDbDMNQtd7AdkyaSUpKwl1djic2hfiUPs9QbigqP3a2tUUmtmPSiNPpxLutGo8jmbTMNLPN0RVVA47NmDGDxx9/3GwzDEPVegPbMWmmqqoKGn3IqBh2GZhjtjm6UlVVZbYJhrB27VpWr15tthmGoWq9gT2PSTM5OTnQHAOxkLJLqtnm6EpOjlqOtpXZs2dTX19vthmGoWq9gd1i0kxVVRVRnuDaJJUmV4Lav7y2tsjEdkwa8Xg8RPsTAYhXzDF5PGpNf2jltttuY9asWWabYRiq1hvYjkkzBQUFCKneAl5Qdz7M5s2baWhoMNsMw1C13sAeY9JMUVERQqaA9BGXoNa/TdX8ZC+88IKdVy5CsVtMGolPcEJUClGyESHUWY4CkJKSYrYJhmFri0xsx6SRqq3F+GKSiYpSL85yTIxaLcBWbr75ZmbOnGm2GYahar2B7Zg0s3ntCjyOFKJj1RtwVDXFUXV1tdK511StN7DHmDQjtlfhcQzE4dxmtim6k5uba7YJhvD0008rPfitar2B3WLSTG1ZEV5HMvGpaj2RA6ioqDDbBMOwtUUmdotJI4G6ZqSIJjUzstOBhyIQUCtaQivTpk1j27ZtPPfcc2abYgiq1hvYLSbtNAX/VdmD8002RH/y8vLMNsEQmpublR4gVrXewHZM2mmOBiCjQL1+varB1B5//HGuv/56s80wDFXrDWzHpJkoXzCfXKJCiS5bSUtLM9sEw7C1RSbqtnN1RrSsk1NtOYrKTJ06le3bt/PMM8+YbYpNN7FbTBqJkkmIgIdYZ7TZpuiOytk23G632SYYhsr1ZreYNCJJISrQoNxyFID8fPUG9AEeeughmprUm6nfiqr1BnaLSRNSSohKJopGs00xhJKSErNNMAxbW2Rit5g0UF1XiT86hagoNWcRR0Wp+ft09dVXU1dXx+zZs802xRBUrTewHZMmNq//GY8jmZjYarNNMYSsrCyzTTCE+Ph4oqPVGxNsRdV6A7srp4mtm9bijU0iJl6abYohlJaWmm2CIcyaNYurrrrKbDMMQ9V6A9sxaaK2dAuIKOVC6raSnq7eMptWbG2Rid2V04CrKphpIyWzn8mWGIPP5zPbBEOYMmUKTU1Nyo4xqVpvYKEWkxDiHCHEaiFEQAgxut32QUKIZiHEypbXU52cny6E+FQI8XvLX128yMvTbqG59DgANv0Qz8vTbtGjWEuhaoqj/v37k5CQYLYZhqFqvQEIKa0xbiKE2AMIAP8Dpkkpl7VsHwR8IKUc0cX59wI1UsqZQoibgH5Syn92dd3Ro0fLZcuWhdz38rRb2F53OIHo/08JHuV3k5i6lEmz7taozPo0NzcTHx9vthmGYGuzLkKI5VLK0aH2WabFJKVcI6Vc24sixgMvtbx/CTi9tza5qkft4JQAAtFxuKpH9bZoS1FUVGS2CYZha4tMImWMabAQYgVQD0yXUi4NcUyWlHILgJRyixAis7PChBBTgCkQDB1RWFhITk4OVVVVeDweCgoKKCoqwusIPbjodaTz9aJF7Doi2Iirra0lPz+fkpISoqKiyMrKorS0lPT0dHw+H/X19W1lOhwOMjIyKCsrIyMjA5fLRWNjY9t+p9NJWloa5eXlZGZm0tjYSFNTU9v+hIQEkpKSqKysJDs7m9raWlwuV9v+pKQknE4nVVVVITWlpKQQExNDTU0Nubm5VFRUUFtbS1NTE8XFxW0LQyNdUyAQYMaMGdTU1LTFY1JBU15eXls9uVwuCgsLI1ZTOPq0KyeEWARkh9h1q5RybssxS9ixKxcHJEkpq4UQ+wPvA3tJKXfoYAshaqWUae0+b5NSdjnOFK4r9/TFb+ON6x9aS8BLnOcXskYnMfaqKyJ6slt9fb2SGTduu+023G4399xzj9mmGEKk11u4rlyftpiklMf14Bw34G55v1wIsQEYCnT0JhVCiAEtraUBQGVv7XX2X4E/xBhTrPwC4U3AG7s/Raviee7SN4hNLOTIK89j8F579fayfU5ZWVlE3+CdceeddyqdV07VegMLjTF1hhBiFyFEdMv7IcDuwMYQh84DJre8nwzM7e21J826m8TUpcS6q0FKYt3VJKYu5bJnZ3LpS7dxwrQCEhwLifJvZ7vvKBY+XMJzF81i3kOP4PN6e3v5PiMjI8NsEwzD1haZWOmp3BnAo8AuQC2wUkp5ohDiLOBOwAf4gdullPNbznkWeEpKuUwI0R94C8gHioFzpJRd5rcJ15VrT0lJSaehTAN+Px8++yhbv63FE70//thEHO5KYhN+4/ApZ7Hrvvt2/Q8wkXDaIpmJEyfS1NTEu+++a7YphhDp9WaZrlw4pJTvAe+F2P4O8E4n51zW7n01cKxR9jU2dh5ZICo6mtOumApXwOYNhSx69DkCTbuz3X8UHz9eSZznAdJHwGnXXUuMI9YoE3tMOG2RzLBhw9i6davZZhiGqvUGFmoxmYXWFlN354xIKVnw/BOUL63EF70/vtgkYt1VxDpXceglpzDsgIN6Y7auRPp8mHDY2qxLRMxjsjrdnTMihOCUS6/m0hf/zWm370ViwkKi/dU0yaNY9Gw9z130EO/MvA+v2/zMvirPh7G1RSaW6cpZHaez50kIcvJ35aIH7gPg41eepvSzzXijR1O+KZkXrpxHrONXDpx0PHsdOkYvc7tFb7RZmQkTJrB9+3bmz59vtimGoGq9ge2YNKNXRooTL5wCF0JlWTELH3gcX0U+TY4jWfKSi+/+9zApuzUzbtpU4uL77qZTNdvGyJEjaW5uNtsMw1C13sAeY9I8xlRYWMjw4cMNsWHxmy9Q9PEGfOyP15FKrGcbMbG/sP8FR7LvEccYcs32GKnNbGxt1iXcGJPtmDQ6ppqaGsPj31RXlvHh/Y/i25JDc/weADhdhSQOruP0G6fiTEw05Lp9oc0sbG3WxXZMYdDqmIqLi/s0K8WS915l4we/4ZP74XX0I8ZbR2zUz+x73iHsf9yJul6rr7X1FWeddRZNTU0sXLjQbFMMIdLrLSLmMVmdvk4DdNQZF3DUGVBbXckH9z+IuySb5vgxfDcnipWvPEFCfjXjpl1HYmpqr6+laoqjQw45hMrKXq9Msiyq1hvYLSbD5jEZwdcfvM3a91bgC+yHNy6dGG89seJn9jpjfw465dQel2sFbUZha7Mu9jwmHbDCnJHDTj2HS567m/MfOpa0XRYT4/2D5uhDWTY/gecmP8lrt/6bhm3bul2uFbQZha0tMrG7chqxUojW5NR+XDDjPwB8//FcVs35nkDUSLZVD+PVaUtxsILhp43g0NPP0lSelbTpybhx42hububTTz812xRDULXewHZMmklKSjLbhJAcdOJ4DjpxPE2N9cydNYumjWm4nIey4qNofnvvfzgHlHHatGtJDbMS3araesuxxx6r9DiMqvUG9hiTJeYx6c3yzxfy8+tf4veMxOPchWjfdhyBlex60u4cee6EnY6PJG3dxdZmXezpAmHQ6phqa2sjbqatq6mJuQ/cR+O6RNzOfZBRMcQ1byAuczOn3PBX0rNzgMjUphVbm3WxHVMYtDqmTZs2MWjQIOMNMoifl37G8tmf4Hfvi8eZRbSvidjASgYdk8+uY46IaG2dcfLJJ9Pc3MySJUvMNsUQIv2etOcx6YDL5TLbhF6x7+HHsO/hx+Bxu3j/gftoKIzDHXcAhV/F8seni3FkbOKk6y8nMy9yJ+x15LTTTqO8vNxsMwwj0u/JcNgtpgiax6Q3v36/lB9f/JBA8964nQOI9jUT6/+ZgUdkccLFl5ptni6oWG+tRLo2ex6TDqg4Z2Tvgw7nkidncvANe5O163fEen/FE7Mfv38/mGcnv8hLU6ez5Y8NZpvZK1Sst1ZU1mZ35TSi8qPZtH7pnH1jMPX52p++55tn38PPXjS6juH9u9cR532TAYekc9KUKxBCmGytdo477jjcbjdLl4ZKQxj5qHxP2o5JIyoH5Wqvbdh+BzHsiYPweb3Me/xBtq3w4XaMYuOKOJ67+BViUtZy7DUXMnCo9R9Tn3feeUrHxVb5nrQdk0aqqqqUTZcTSltMbCxnTv0HABtWreDLp97EL/dgu+dYPrhvEw7v+2QekMTJU64gJtZ6CRYALr/8cqXzyql8T9qD3xoHvyM962k4tGrzeb18+PRjVP3QhMcxikC0kzhXKTFJhRx51QQG77l3H1jbPex6sy72dAEdqKqqiuibIBxatcXExjL+6uvhati0dhWfPz4bf/NQtvuO5aMHS3F4F5IxKpZTrrrGEq2oo446CpfLxXfffWe2KYag8j1pOyaNeDzmZzMxip5oGzRsBBc/MpOA38+CZ5+g8ttaPDH7UbImnhcuf5uYxN84/PIz2W3kfgZYrI2LLrqILVu2mHZ9o1H5nrS7cn/ieUyt6KWtZMM6Fj36PP5tu+OKH4wIeIlz/0K/EYJxU/9mSivKrjfrYs9j0gGV54zopS1v16Fc9NBMLn3pUnY9cB1O7/d4Y/dgy4b9eGHKHF648lbW/Nh33Sqv18v69ev77Hp9jcr3pN1i0thiKisrIycnpw8s6nuM1Fa+eSMfP/gMvqrBuBJ2QwR8xLl/JXUPD6df/3di4uIMuS4Ex5g8Hg/ffPONYdcwk0i/J+3Bbx2IiVH3X2WktuyBQ5j8wH8B+HT28xQvLsYXsx8Vm5J44a/vExP3K6MnncDehx6h+7Uvu+wy6uvrdS/XKqh8T9otJgXjMXWXvtZWVbaZBQ8+ibciH1fCUETAT5x7Fcm7NzF+2g3E6Thx0K4362KHPQmDVsfU0NBAcnJyH1jU95ip7bM3XmbTJxvwiv3wxSYT664mJvYX9rvgKEYeeWyvym5qaqKhoYGsrCydrLUWkX5P2l05HaioqIjomyAcZmo7ZsIkmAA1W7fwwazH8G7JpTnuaL5+LcDy5x8mYXAdp994A/E9SPY5duxYmpub+f777w2w3HxUvidtx6SRQCBgtgmGYQVt6bsMYNI9wQQLX777Br9/+Bu+6P2oKU/hlWs/IiZ6JSPOOYQDTxirucwrr7xS6XlMVqg3o7C7chq7ck1NTcpmpbCqtvqaaubPehBXSTau+OBYirN5DfGDahg/bSqJKV0n+7SqNj2IdG32PCYdKC4uNtsEw7CqtpT0/lxw911c+vI1jDy1nnj/F/hjctlWeTivTl3M85fdxrcfvt/p+XV1daxevbrvDO5jrFpvemC3mDS2mMrLy8nOzu4Di/qeSNLWWF/LvPseoLkoA1f8ngA4m9fiHFjBuBv/TnK74Pyqz2OKpHoLhT34baMMSSlp/GXGnQD8+MkH/Pr2d/hiRlFbvQev/X0JsWIFQ08bwcavVjA8ajT+pCSevvhtnP1XMGnW3SZbb6MVu8Vkz2OKeG3N27cz975ZbN+Yiit+LxDRIAMg/n+kIsrvJjF1qVLOKdLrzW4x6UB+vjrZQzoS6driExOZcMftAPz02Sf88FoT/pgEGpvrAEiKTyUQHYerepSZZupOpNdbOOzBb42UlJSYbYJhqKRtv2NOwB8dXHH/7Kf/5tlP/922z+tIN8ssQ1Cp3jpit5g0EhWlrg9XTVuspwZvXH+O3efsnbarhGr11h51lemMqssaQD1tzv4riPK72XvQoew96FAgOMbk7L/CZMv0RbV6a4/tmDRSWlpqtgmGoZq2SbPuJj75S+qbaqjfXk2su1q5gW9Qr97aY3flNJKertb4RHtU1DbxnjvZo+AAouQ21m4pAs4x2yTdUbHeWrEdk0Z8Pp/ZJhiGqtpOGDkBR+Bns80wDFXrDeyunGZUDjimqrY98w9kRN4Qs80wDFXrDSzkmIQQ5wghVgshAkKI0e22XyCEWNnuFRBCjAxx/h1CiNJ2x2lfhq6BgoICPYuzFKpq29ZYybbGBrPNMAxV6w0s5JiAVcCZwJftN0opX5VSjpRSjgQuBDZJKVd2UsaDrcdKKRfoaZzKgd9V1fbSZzN55ou5ZpthGKrWG1hojElKuQZACBHusPOB1/vEoA44HA4zLtsnqKrtpP0uwOFXa4pAe1StN7CQY9LIecD4MPuvEUJMApYBN0gpt4U6SAgxBZgCkJeXR2FhITk5OVRVVeHxeCgoKKCoqIiUlBRiYmKoqakhLS2N9evXEwgEyMvLo7i4mLSWley1tbXk5+dTUlJCVFQUWVlZlJaWkp6ejs/no76+vq1Mh8NBRkYGZWVlZGRk4HK5aGxsbNvvdDpJS0ujvLyczMxMGhsbaWpqatufkJBAUlISlZWVZGdnU1tbi8vlatuflJSE0+mkqqqqS025ublUVFTQ1NREU1OTUpo8bjfD8/bH6a2lvLxcCU2h7r3CwsKI1RSOPl3EK4RYBISK03CrlHJuyzFLgGlSymUdzj0IeFZKuXcnZWcBVYAEZgADpJSXdGWTvYhXTW0+r5f/TH6TON833PTWE2abYwiRXm+WWcQrpTyuF6dPIEw3TkpZ0fpeCPEM8EEvrrUTGRkZehZnKVTVNnvJLKLkNm4y2xCDULXeIEK6ckKIKIIz5DpNPiaEGCClbA3wfAbBwXTdcLlcehZnKVTVdsroyTj8P5lthmGoWm9goadyQogzhBAlwCHAh0KIj9vtPgIokVJu7HDOs+2mFtwrhPhVCPELcDRwvZ72NTY26lmcpVBV2+45+zJsgLqhQVStN7BQi0lK+R7wXif7lgAHh9h+Wbv3FxpmHGrPGVFVW0XtZpw+tSIKtEfVegMLtZisjspzRlTV9vqXD/LS17pOZ7MUqtYbWKjFZHWcOqatthqqaht34KU4/MvNNsMwVK03sB2TZtLaZd9QDVW1Dcnei3hfRdcHRiiq1hvYXTnNtE7SUxFVtZXV/EFJjbqOSdV6A7vFpJnMzEyzTTAMVbW99dWjRMtt3Ga2IQahar2B3WLSjMqPZlXVdsbBUzj7wGPNNsMwVK03sFtMmmlqajLbBMNQVVtB5nCcPnXDz6pab2A7Js2oPGdEVW0lVeuJ86s7xqRqvYHdldOMynNGVNU255sneP27T802wzBUrTewW0yaSUhIMNsEw1BV29mHXoXD/4PZZhiGqvUGtmPSTFJSktkmGIaq2vIydiPeV2y2GYahar2B3ZXTTGVlpdkmGIaq2ooqC9m4Vd3Bb1XrDewWk2ays0PFt1MDVbW9993TRMtt3GW2IQahar2B7Zg0U1tbq+wSAFW1nTvmWuJ835tthmGoWm9gOybNqByUS1VtOemDcfo2dn1ghKJqvYHtmDSj8pwRVbVtLF9NnF/dMSZV6w3swW/NqDxnRFVt8354jneWfW62GYahar2B3WLSjMqPZlXVdv4R1xPn+85sMwxD1XoD2zFpRuWgXKpqy0obiNP3u9lmGIaq9QZ2V04zVVVVZptgGKpq+73sZ9Zu2WS2GYahar2B3WLSTE5OjtkmGIaq2j5c9hLRchv3mm2IQahab2A7Js1UVVWRkpJithmGoKq2iUdNI873jdlmGIaq9Qa2Y9KMx+Mx2wTDUFVbRkoOTl8/s80wDFXrDWzHpBmV54yoqq2wZDlx/k1mm2EYqtYb2IPfmlF5zoiK2mb/8zY++ulV5v68iqcvfpuXp91itkm6o2K9tWK3mDSial8e1NP28rRbaG44gsnH7AuAN64//rrDeXnaLUyadbfJ1umHavXWHtsxaSQmRt1/lWraXNWjCMTF0S/p/7OIBKLjcFWPMtEq/VGt3tpjd+U0UlNTY7YJhqGaNq8jHYDfin/gt+IfdtquCqrVW3vUdbk6k5uba7YJhqGatlhPDd64/nyy8g0A9sw/sG27SqhWb+2xHZNGKioqSE5ONtsMQ1BNm7P/Cvx1h3PJcdPbtkX53Tj7rwDOMc8wnVGt3tpjd+U0EggEzDbBMFTTNmnW3cQnf0lKQjop8f2IdVeTmLpUqYFvUK/e2iOklGbbYCqjR4+Wy5Yt6/K4pqYmZbNSqKjN5/VyzWn34Aj8wiOfvGW2OYYQ6fUmhFgupRwdap/dYtJIcbG62TZU1bb4lzl8/Ku6oXVVrTewx5g0o2psZVBX22XH306cf6nZZhiGqvUGtmOyUZik+FScvsjt6vyZsbtyGqmtrTXbBMNQVdvKjUtZvmmN2WYYhqr1BnaLSTP5+flmm2AYqmpbsuo9ouU2s80wDFXrDewWk2ZKSkrMNsEwVNV2xYl3cvXxZ5tthmGoWm9gt5g0ExWlrg9XVVt8XBJOn8JxsRWtN7BbTJrJysoy2wTDUFXb8vWf8+NGdceYVK03sB2TZkpL1U2cqKq2pb/N5/M1y802wzBUrTewu3KaSU9Xa2V6e1TVdtXJdxPn/8JsMwxD1XoD2zFpxufzmW2CYaiqzRHrJE7Emm2GYahab2B35TRTX19vtgmGoaq2H9Z9yrfrfzXbDMNQtd7AbjFpRuXA76pq+6ZwodLzmFStN7BQi0kIcZ8QolAI8YsQ4j0hRFq7fTcLIdYLIdYKIU7s5Px0IcSnQojfW/7qmrdH5cDvqmq79pR7uf7k8802wzBUrTewkGMCPgVGSCn3AdYBNwMIIfYEJgB7AScBTwghokOcfxOwWEq5O7C45bNuOBwOPYuzFKpqi46OISYq1K2iBqrWG1jIMUkpP5FSto7mfQfktbwfD7whpXRLKf8A1gMHhihiPPBSy/uXgNP1tC8jI0PP4iyFqtq+W/sRX6/7xWwzDEPVegMLOaYOXAIsbHmfC2xut6+kZVtHsqSUWwBa/maGOKbHlJWV6VmcpVBV23drP+Hr39V1TKrWG/Tx4LcQYhGQHWLXrVLKuS3H3Ar4gFdbTwtxfK/CbgohpgBTWj42CiHW9qY8G2sjRKhbyMYCdDp636eOSUp5XLj9QojJwKnAsfL/Y/6WAAPbHZYHhPqpqBBCDJBSbhFCDAAqw9jxNPB0t4y3sbHpMyzTlRNCnAT8ExgnpWxqt2seMEEIESeEGAzsDvwQooh5wOSW95OBuUbaa2NjYxyWSUYghFgPxAHVLZu+k1L+tWXfrQTHnXzAVCnlwpbtzwJPSSmXCSH6A28B+UAxcI6UUq1EYjY2fxIs45hsbGxsWrFMV87GxsamFdsx2djYWA7bMdnY2FgO2zHZ2NhYDtsx2djYWA7bMdnY2FgO2zHZ2NhYDtsx2djYWI7/A4jUrBzPZcxaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot velocity profiles at the sill\n",
    "\n",
    "fig1 = plt.figure(figsize=(4,8))\n",
    "fig1.tight_layout()  \n",
    "ax1 = fig1.add_subplot(111)\n",
    "\n",
    "depth = xr.open_dataset(path1+'/state_001.nc').Z.data\n",
    "\n",
    "ax1.plot(Usill[0,:],depth,'o-',lw=1.5,label=r'$U_t=0.07,U_e=0.15$')\n",
    "ax1.plot(Usill[1,:],depth,'o-',lw=1.5,label=r'$U_t=0.13,U_e=0.15$')\n",
    "ax1.plot(Usill[2,:],depth,'o-',lw=1.5,label=r'$U_t=0.20,U_e=0.15$')\n",
    "ax1.plot(Usill[3,:],depth,'o-',lw=1.5,label=r'$U_t=0.66,U_e=0.15$')\n",
    "ax1.plot(Usill[4,:],depth,'o-',lw=1.5,label=r'$U_t=1.32,U_e=0.15$')\n",
    "ax1.axvline(x=0,color='k',linestyle=':')\n",
    "ax1.legend(loc='best',fontsize=12)\n",
    "ax1.set_xlabel(r'$U\\;[m/s]$ at the sill', size = 15)\n",
    "ax1.set_ylabel(r'Depth [m]',size = 15)\n",
    "#ax4.set_xlim(-4,4)\n",
    "ax1.set_ylim(-20,0)    \n",
    "ax1.xaxis.set_label_position('top') # this moves the label to the top\n",
    "ax1.xaxis.set_ticks_position('top') # this moves the ticks to the top\n",
    "#ax4.yaxis.set_visible(False) # This erases the y ticks\n",
    "ax1.grid(ls='--', lw=.5)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0ad9a66",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}