notebooks/model_swan_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of generating a SWAN simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.core.options.set_options at 0x7fecfb6965e0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import rompy\n",
    "import intake\n",
    "import numpy as np\n",
    "from rompy.swan import SwanModel, SwanGrid\n",
    "import xarray as xr\n",
    "xr.set_options(display_style='text')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiate a model simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Default context: (from cookiecutter.json)\n",
      "OrderedDict([('run_id', 'run_0001'),\n",
      "             ('compute_start', '20200221.040000'),\n",
      "             ('compute_interval', '0.25 HR'),\n",
      "             ('compute_stop', '20200224.040000'),\n",
      "             ('out_start', '20200221.040000'),\n",
      "             ('out_intvl', '1.0 HR'),\n",
      "             ('output_locs',\n",
      "              OrderedDict([('coords',\n",
      "                            [['115.61', '-32.618'],\n",
      "                             ['115.686067', '-32.532381']])])),\n",
      "             ('cgrid', 'REG 115.68 -32.76 77 0.39 0.15 389 149'),\n",
      "             ('cgrid_read', ''),\n",
      "             ('wind_grid',\n",
      "              'REG 115.3 -32.8 0.0 2 3 0.3515625 0.234375  NONSTATION '\n",
      "              '20200221.040000  10800.0 S'),\n",
      "             ('wind_read', \"SERIES 'extracted.wind' 1 FORMAT '(3F8.1)'\"),\n",
      "             ('bottom_grid',\n",
      "              'REG 115.68 -32.76 77 390 150 0.001 0.001 EXC -99.0'),\n",
      "             ('bottom_file', 'bathy.bot'),\n",
      "             ('friction', 'MAD'),\n",
      "             ('friction_coeff', '0.1'),\n",
      "             ('spectra_file', 'boundary.spec'),\n",
      "             ('_generated_at', ''),\n",
      "             ('_generated_on', ''),\n",
      "             ('_generated_by', ''),\n",
      "             ('_template', '')])\n"
     ]
    }
   ],
   "source": [
    "new_simulation = SwanModel(run_id='test_swan',\n",
    "                           template='../rompy/templates/swan',\n",
    "#                            checkout='88c529557db', #allows for the optional specification of a particular commit for reproducability\n",
    "                           output_dir='simulations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'_template': '../rompy/templates/swan', '_checkout': None, 'model': 'SWAN'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting start out empty, only need to add values from the default context we want to replace\n",
    "new_simulation.settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse the parameters of the grid defined in the template and create a SwanGrid object to conveniently get the bounding box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid=new_simulation.grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 720x509.002 with 1 Axes>, <GeoAxesSubplot:>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAGsCAYAAAAbq4Z0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABbNElEQVR4nO3deVyVVf7A8c9BTRFwFwSXwgVcEEhQcakoJzPLpqZ1JqeapilnaZqa+s1U0960jK3aTJvZ6rSqZU77QqWGCoqAyiK4pOBFcWMRN87vD+7Fq3LhIvdZ7r3f9+t1Xy94eO7zfI/Ll3PPc873KK01Qggh7C3E6gCEEEK0TJK1EEL4AUnWQgjhByRZCyGEH5BkLYQQfkCStRBC+AFJ1sdRSs1VSlUopfLdjl2ulFqrlKpXSqW6HT9NKbVfKZXjfL3g4Zr3K6W2uZ031e1nM5VSWUqpszxcM0cpdY2RbXaLxedtd557s1Kq0Hmdf7kdD+i2K6XedTtnk1Iqx+1ntmi7Qe1OVkplOs/JUkqNcfuZLdrtl7TW8nJ7AWcCo4B8t2PDgHggA0h1O36a+3nNXPN+4PYmjg8FZgKdgfdac00/avvZwFdAR+f3kcHS9uOu/yRwr93abtDf+RfA+c6vpwIZdmu3P77aN5G/g5rW+nul1GnHHVsPoJTy9e3aAfWABnx+8dYyqO2/Bx7TWh9wXq/CeTwY2o7z/Qq4AjjHecg2bTeo3Rro4vy6K1Dm/No27fZHMgzSdrFKqdVKqe+UUmc0c96flFK5zo+d3QG01mtp6GUsAZ53O3fQcR8Lm7uulbxpexxwhlJqufO80RA0bXc5A3BorYvB79vuTbv/AsxUSv0EPAHcCX7fbstJz7ptyoEBWutKpVQK8KFSaoTWet9x5z0PPERDj+IhGj4SXw+gtb65ieuWaK2TjQvbJ7xte3ugO5AGjAbeU0oN1A0Cve0uvwTedj/gp233tt2/B27VWs9XSl0BvAL8DPy23bYgPes20Fof0FpXOr/OBkpo6Ekef55Da31Ea10PvAyMOf4cf+Nt24GtwAJncl5Bw8fgXuZF6nutaDtKqfbAL4B3zYvQGK1o97XAAufX7xMA/97tQJJ1Gyileiul2jm/HggMAUqbOC/a7dtLgPzjz/E33rYd+BDnWK1SKg44BdhpUpiGaEXboaFHWaC13mpWfEZpRbvLgLOcX58DFJsTYYCz+gmn3V40fFwtBw7R0Cv8LQ0JditwAHAAnzvPvRRYC6wBVgHT3K4zB+eTdOBNIA/IBRYB0c3c/zRgP5Dj9vqzH7f9FOAtGn5BrQLOCZa2O79/DZjhxf0tabtBf+cTgWznecuBFLu12x9fyvkHJoQQwsZkGEQIIfyAJGshhPADkqyFEMIPSLIWQgg/IMlaCCH8gCRrIYTwA3673HzKlCl65077rq2oqqoCICIiwuJIzBfIbT906BAOh4OOHTvSu3fvE35uVdvr6+spLi4mLCyU6OgoU+8NUFOzH4CwsFDT730ycnLyP9daT7E6jtbw22S9c+dOsrKyrA7Do4yMDADS09MtjcMK0nbz237JJRczaFA/XnhhJu3atTP13gCZmbkApKUlmn7vk9G9+xC/K3ngt8na7uLimiwVERSk7eapq6vjN7+5js8++5yCgmWWJGqA2Ni+ltw3mEiyNkhMTIzVIVhG2m6Ouro6fvGLS+jQAXJzM+ja1bphp6ionpbdO1hIshbCD7kSdceOIbz00hN06NDB6pCEwWQ2iEHKysooKytr+cQAJG03tu1aay699Be2StQORyUOR6XVYQQ06VkbpKioCAjOIQFpu7Ft//LLL1myZAkbNiy3RaIG2LhxGyDDIUaSnrUQfuTbb7/l6qt/xZtv/ts2iVqYQ5K1EB4cOHCAu+66i3nz5lkdCjt37uSOO+7giisu59VXn+XMM8dZHZIwmSRrITz4xz/+waJFC7nzzr/x73//27I4Pv74Y/r168fmzUUsWDCXiRPHWhaLsI6MWQvhQV1dHf379+Wxx+7moouuoaCggHvvvbfJlYtGqKmp4cEHH+TVV+fyv//NIyUlyZT7CnuSnrUQHkydOpXVq/Po1asnH3/8Fjt2bKVv377cdNNNvPvuuxi9y9I11/yaNWuy+PDD1yVRC0nWQnhy3nnnMXToUBYv/pIBA/rywgszWbPmW7p0ac8DD9zLoEEDG2d/+NJ3333HeedN5osvvuA//3mM4cODd0WoOMpv92BMTU3Vdq4NIgLDr389ncTEwfzud9OPOX7kyBHefPN9/vWvf5Obm0evXr4pNfHf//6Xv/zlFv7wh+u48cZr6NzZPwoj+Zvu3Ydka61TrY6jNSRZC9GMb7/9lssuu5SePXsAoJTi73+/mUsumQrApZdez7ZtDv7xj3vo0aMHU6ZMISTk6AfWnJwcNmzYwOTJk+nSpYvH+yxZsoQHH3yAnJwcPvzwdelNG0yStYkkWQuzVFRUsHv3bgDKy8uZOnUqJSUrCA3txP79dbz//iIyMpZRULCBkJB23H33P5g37y02b97M9u0Ohg0bQnFxKVFRUVx++RUkJiailEIpxbp163jzzTfYvn07v/vddG688RpLa3wEC0nWJrJ7snbFlprqV/8efCLQ2x4VFcmAAX3p378vt902g4SEoQDU1R1g3rwPWbNmPX379mTq1ElERvYmMrIXxcWlVFbu5qWX3mT//jq0blg23rFjB2666RqGD4+jR4/uFrfs5OXlFQMwcuQQiyPxjj8ma5m6Z5Dq6mqrQ7BMILf9rrvu5LTTTuPqq6ezZk0O06ZN5+GH/87VV19Gp04dGTFiGCNGDDuhrnNc3CAAxo3zq/zgNdfmA8I4kqyFaIXHHnscrTUDBgzggw/mc8UVl/Pgg0+RnDySESPirQ5PBDBJ1kK0wu23305+fj5aayIjezNp0iQqKytZvjxbkrUwlCRrIbxw4MABQkJCmDlzJunp4/nxxx+pqKigQ4cOREX14cknH+f6639ldZgigMmiGCFasHbtWvr2jWHo0HieffZZ8vMLueeeexqr3iUlJVFQUMy6db5fICOEi/SshWjB999/zymndGDXrl0MGTKEHTt2HvPz0047jdTUVFauXE18/DCLohSBTpK1QaKjo60OwTKB1vYZM2YQERFB165dOf/885s856abZvDww/fz4otPmxydPURG9rA6hIAnydog8fHB+7Ap0NqulGL69OnNnvPDD98THd2HgQP7mRSVvQRru80kyVqINlq9ejVz5sxh8WLrNykQgUseMBqkqqqKqqoqq8OwRCC2fd++fVx44QV8/fXXJ/zsj3/8Q+NKxurqWqqray2I0FrB2m4zSbI2SHZ2NtnZ2VaHYYlAbPvcuXP53/8+Ye3atSf87Ne/voaXX36L7dsryM/fQH7+BgsitFawtttMkqyF8MKf/vQn1q1bx5/+9KcTfvb73/+e008fxaxZLxu+IYEIXi0ma6VUJ6XUCqXUGqXUWqXUA87jM5VSBUqpXKXUQqVUN2/f6/xZjFLqG6XUR0qpcOexeKVUhlIqRym1Xin1kg/bKsRJa9++PcOGDTum/Km7OXNe4bPPvmPLlm0mRyZ8xe65zpue9QHgHK11EpAMTFFKpQFfAgla60SgCLizFe8F+DNwMzAHcD1qnwU8rbVO1loPA2Z7EZ8QlouMjCQubgiHDh2yOhRx8myd61pM1rqBq4xaB+dLa62/0Fofdh7PBE6Yu+Ppvc7v2wH1zpdyHosGtrq9P6+l+ISwg/fff581a9YQGembHWOE+eye67yauqeUagdkA4OBf2utlx93yvXAu61873PAm8BewFVU4WngG6XUMuAL4FWt9Z6mrltVVUVGRkaT8cbFxRETEwNAWVlZs/vkpaenN36dlZXlsbxndHR04/zhqqqqZh+gpaSkNH5dWFhIeXl5k+eFh4cfU/PZU3vAHm2KiGgoit9Sm9wFSptSUlJ44okn6NChA8nJycf8/ODBg7z99n+5667badfu6H8ph6OSjRs9D4u4l1HNyyv2WGY0MrJH4zzm6uraZh/kJSQMJjy8MwClpVupqNjV5HlhYaHH1J7OzMz1eM3Y2L5ERfUEPLdp61bHCcfs3CZP7JjrXLx6wKi1PqK1TqbhN8oYpVSCW4B3A4eBJieZenqv1nqz1vpMrfU0rXWV89irwDDgfSAdyFRKdfQmRiGMNHfuXJ566gluvfXWEx4ifvLJ/+jbN5rBg2Mtik74ip1zXat3ilFK3QfUaK2fUEpdC8wAJmmtW5xk6f5eL++VD1yrtT6he2T3nWJc84xdPbdgEohtnzVrFq+/Ppe6ujquvPJX3HvvvY0/Gzt2DA89dAdjxpzeONc4PLwz5eUOXnvtHY4cOYJSiiuvvDhgE7p7u/2BNzvF2CXXuXgzG6S36+mnUioU+BlQoJSaAvwNuMhT8J7e28y9piilOji/7gP0BPzy8XpERERAJavWCMS2b9myhQkTRrN79x5mzpx5TO96yJAhPPXUC8ya9TLjx08lLi6NmJhEkpPPobJyP926xXDoUAcmTLiQuLhx5OTkW9gSY4SHd/abRO2J3XNdiz1rpVQi8DoNg+QhwHta6weVUhuAjkCl89RMrfUMpVQMMEdrPdXTe5u511PABUCd89BMrfVbTZ1r9561CByHDx9m6NB42rdXjB07nldemUv79kfHpuvq6pg5cya7du0iPT2dc889F4CQkBA6derUeN6hQ4f43//+x/Tp0wkJUYSHh/O7302nc+dOtG/fniuu+DkREeEn3F/4XlM9a7vmusb3+Oskfrsn68LCQiDwihp5IxDbPmPGDF588UXefPPNZos6edP2uro6Dh48yKpVq5g/fz4AGzZsYMWK5Zx++kheeeUZv9vhvLS0YWKDvxR08scNc2UFo0HKy8s9zi4IdIHY9ueff54dO3a0WH3Pm7Z36tSJLl26kJ6ezuzZs5k9ezaffPIJa9bkEh+fwKWXXs/evft8Gb7hKip2eZylIXxDkrUQXlBK0auXcXOolVL069eP5557jnHjJvKLX/hfwhbGkmQthI0opZg1axYjRiTyj388ZnU4wkYkWQthM0opHnjgAb7++gc+/vhzq8MRNiHJWggbio2N5emnn+HBB5+iqqrp1ZoiuEiyFsKmLr/8cs444yz+/veHrQ5F2IAka4OEh4efUCcjWEjbfdP2kJAQ7r//fr766ns+++wbn1zTKGFhoYSFhVodRkCTedZC2NwHH3zAc889w4cfvmZ1KAFD5lkLIXxu9OjRrF1bQEnJJqtDERaSZC2EzZ166qmMGDEiIGuKCO9JsjZIRkZGs7WcA5m0PcPn173xxpt49NFZHD58uOWTLZCZmXtS9aOF9yRZC+EHrrjiCjp3DueDDz62OhRhEUnWQviB9u3bc++99/Kf/7xmdSjCIpKshfATw4YNo6am9oSdakRwkGQthJ8YNGgQHTuG8vnn31odirCAJGsh/ESnTp1ISkqU5edBSpK1EH4kJCSE+noZBglG7Vs+RZyMuLg4q0OwjLTdOCEhIbYcs46N7Wt1CAFPkrVBYmJirA7BMtJ24yil0Lre0HucjKionlaHEPBkGEQIP6KUor7efslaGE+StUHKysooKyuzOgxLSNuNafuOHTv49NNPGTbMfsNMDkclDkdlyyeKkybDIAYpKioCgnNIQNpuTNsdDge9evVk1KhEn1+7rTZu3AbIcIiRpGcthJd27NgRtDVPhPUkWQtbWLFiBaNGjeIXv/gFTz75JJmZmRw8eNDqsBrV19dzzjnncOGFFzBuXJolBZV2796NUsr0+wp7kGEQYQvPP/88q1evZvXq1SxcuBBoWAQyZswYJkyYwIQJExg3bhw9evQwPbbly5cTGhpKeXkZ06dfxty5b3P99dfz009bePvtd+jTp48pcTzzzNNcccU0U+4l7EeStbCFpUuXAnD/nXeytayMpZmZrC8s5Pvvv+f7779vPG/48OGNyXvChAkMGjTI0N7m3r17SUtLo2fPHjz11NN8+uknPPbYY9x9993ExETxyCOPMGvWLMPu727Xrl0MHTrElHsJ+5FkLSxXUVFBcXExnTt35q7bb6dDhw4AVFZWsmzFCpZmZrI0M5OVq1axbt061q1bx8svvwxAVFQU48ePZ+LEiUyYMIHTTz+dU045xafxtWvXjtraWvbv38/bb7+D1poDBw5QXFzEww+bs5ltZmYm+fn5pKYmmXI/YT+SrIXlli1bBsDY1NTGRA3Qs2dPpp1/PtPOPx+AAwcOsConh6XLlzcmcIfDwcKFCz0OnYwfP57u3bufdGxdu3Zl7dq13H33XWzevInDhw/z4IMPsHv3HubMeYWQEHMe+9xzzz/47W9/Rc+e5g8DCXuQDXOF5e644w6eeOIJ/nHHHTx0zz1ev09rzYaSEpYuX86SH39kaWYmBc6pc+5cQyeu3vfAgQNbNXRy11138e23X7B+fTHTp/+a1157lZqaWsrLy00Zr162bBmXXHIxmZmf0L17N8PvFwz8ccNc6VkLy7nGqyekpbXqfUophgwezJDBg7nu6qsB74dO3Me9Wxo6SUpK4sknn+Sss87koYceYvfuXezYsYOePY2fU7xnzx5uuukm/vCH30iiDnLSsxaWqquro2vXrhw6dIjdW7bQtWtXn17ffejE1fveWXnsSjtvhk6ysrJISEigU6dOPo2vJY8++iiffbaYDz98Tabt+ZA/9qwlWRvEFVtqql/9e/CJ1rR9yZIlnHHGGYwcMYLcH380OjS01hRv2HDMuHdTQycjRow4pvft7dCJL//e33zzTf72t//jww9fIy5uUJuvZ6S8vGIARo70j9kq/pisZRjEINXVwVsgvjVtX7JkCQATWzkEcrKUUsQNGULckCH8Zvp0AHZWVvLjcUMna9euZe3atbz00kuA90Mnvvh7r6ys5J///Cdvv/1fv0jUADU1+60OIeBJshaWOtnxal/q1cSsk+ycnMbk7Zp1smDBAhYsWABAaGjoCQt22jLrZPv27bzzzjs4HA7+9a9/cfHFU1m06A2GDBnokzYK/yfJWlimvr6+cdqelcn6eB07dmT82LGMHzuWO265xePQyXfffcd3333X+L4RI0YQGxtLQkIC/fv3b3HoZNOmTTz00EO8//77HDp0iIsumkL37l1ZsmQxw4b5x3CCMI8ka2GZwsJCdu3aRUx0NKcOGGB1OB55GjpZ5pa83YdOFi9ezGOPPUafPn1OGDrZuHEjK1asIDMzk1deeYXLL7+IH3/8H126RBAREW5xS4WdSbIWlnEfAvG3mQ69evbkoqlTuWjqVODo0Mkb77xD/rp1FBYVsX37dubPn8/8+fOBhqQfEhLC4MGx9O/fly+/fI+EhGFWNkP4EUnWwjKNyXrsWIsjaTvX0MnBI0cAOGvcOIo2bGBpZia33v5XVPv27K2q5siRIxQWbqCwcANfffUdQ4cOIS0thbFjU0hLG8Wpp/b3u19cwhySrA0SHR1tdQiW8bbtdni46GvRUVFAQy86fsgQ4ocMYdaTT3DLLb/lzLPGszxrjfOVw+rctRQUFFNQUMxrr70DQFRUb8aOHcXYsSmMHTuKxMThxyzBt6vISFkGbzSZZy0sUVFRQVRUFJ07d2bPTz/5RUI6WXl5eVx88c954IH/48ILz208Xld3gJy89SzPynG+1rBr955j3hsa2olRoxKdPe8Uxow5na5du5jcgsAj86yF8JKn4k2BaOTIkXz44UdcfPHPARoTdqdOHUkbnUza6GTAWeukdPMxybu4ZBNLl65g6dIVQEOPfejQIY29bxk6CR6SrA1SVVUFQEREhMWRmM+btgfSeLW7KueimIjwY2d2uCdspeCCC8494b1KKYYMOo0hg05j+pUXA1C5azfLs9aQubIhgefkrWP9+iLWry+y1dBJdXUtAOHhnU29bzCRZG2Q7OxsANLT060NxALetL1x5eK4cWaEZJrs3FwA0sePP+FnI0eOZOHCD7nkkotRSjF16s9avF7PHt2ZOjmdqZPTgaaHThyOHSxa9DmLFn0ONAydpKQkNSZwM4ZO8vM3AJCWZr/NfAOFJGthuv3795OdnY1SirTRo60Ox1SJiYmNCVtrzdSpP2vVEEZzQyeu3veG0s0sWbKcJUuWAzJ0EigkWQvTZWVlcejQIRITEnxeZc8fJCYmsmDBQqZNm8brr7/Pq68+Q1jYyQ0fNDV0srNyFyuy15C5co3HoZM+fSK55poruPPOW3zVLGEw2d1cmC5Qx6tbIykpiU2bNjFkyFCmT/8DS5cu99m1e/XswdTJZ/Pg3X/h84WvsTn/Bz6d/yr333kL5597Fj26d2P79gpeeOF1n91TGE+StTBdIM6vPhkhISHMnj2byy67iltuuYdnn32J6uoan9/HNXRyy++v47+vPMPsmfcCMHKkrJ70J5KshansWrzJKiEhIfzxj3/k448Xs2pVAddc8yeqqowtr5u5cg0AY8emGHof4VuSrIWp/KV4k9mGDRvGwoULSUwcxbXX3kx9fb1h91qelQPA2LGjDLuH8D15wGiQlJTg7bU013Z/Lt7kjZTEk5+6FhISwtNPP83AgbEsXvwlF110ng8ja9Aw9W8dAGPGnO6z6yYkDPbZtUTTJFkbJBgXw7gE42IYl+MXw7SWUooPP/yIX/ziEtq1C2ly8Uxb5OSt5+DBQwwdOoRu3Xw3E0cWwxhPkrUwlWsxjIxXe5aUlNS42vHIkXqf9rBlCMR/SbI2SGFhIQDx8fEWR2I+T213OBxs2LCBsLAwktswXGBnhSUlAMQPatu+iSNHjuSjjxZx8cU/R+t6fv7z830RXmOyTkvz7TBdaelWAAYO7OfT64qjWnzAqJTqpJRaoZRao5Raq5R6wHl8plKqQCmVq5RaqJTq1sw12imlViulFrsdi1FKfaOU+kgpFe48Fq+UylBK5Sil1iulXvJBGy1RXl5OeXm51WFYwlPb3Ys3tW8fmP2EcoeDcofDJ9dKSEhg0aKPefjhZxqXkreF1prlWQ0zQXydrCsqdlFRscun1zSb3XOdN7NBDgDnaK2TgGRgilIqDfgSSNBaJwJFwJ3NXOMWYP1xx/4M3AzMAaY7j80CntZaJ2uthwGzvYhP+IlAH682wvDhw3nmmWd57LFZ1NTUtulaxSWb2LV7D1FRvTn11P4+ijCg2DrXtZisdQPXxM8OzpfWWn+htT7sPJ4JNPn5RynVD7jAGai7dkC98+WaFhANbHW7d15L8Qn/IYthTs7555/P2WdP4q9/va9N13Efrw7EmThtZfdc59U8a2fXPgeoAL7UWh+/NvZ64FMPb38G+D9noO6eA14EZgBvOY89DXyjlPpUKXVrcx83hH8J5uJNbRUSEsJ9993HihWr+eyzb076Oq4hEFkM45mdc51XA4da6yNAsvOCC5VSCVrrfACl1N3AYWDe8e9TSl0IVGits5VS6cddczNw5nHHXlVKfQ5MAX4O3KSUStJaHzj+2lVVVWRkZDQZb1xcHDExMQCUlZVRVFTksW3uZTyzsrKorm569Vh0dHTjA7OqqqrGMqBNcZ9nXFhY6HHsOjw8nNTUo5tVeGoP2KNNril5LbXJnatNubm5HDp0iIGxsaxeu7ahTYMGEePcBqvM4aDI+XCuyTa5lRzNys313KaoqMaHe1XV1Y0lS5tsU2Ji41S7wpISj2PN4eHhpLo9EM1wjr03paq6uvGavm7TggUL+eXVV1NX354+fSKbPDdh6KmEdw4FoHTLdip27jka95KGDQzCwrqQl1fMyJFDGn+Wmen5zyk2ti9RUT0BcDgq2bhx2wnnbN164p9dXl4xNTX7m7xmZGSPxoeR1dW1jSVWm2xTwuDGqYGlpVs9jo2HhYV63SZP7JjrXFq1glFrvQfIcN4ApdS1wIXA1brp/cEmABcppTYB7wDnKKXeauI893uUaa3naq1/TsMfTEJrYhT2lJ+fD0DCiBEWR+K/kpOTueqqq1j240qOODfm9daevfvYVlZOx46nMGhQrEERBg475roW92BUSvUGDmmt9yilQoEvgMedF38KOEtrvaPZizRcJx24XWt9YTPnTAG+1lofUkr1AVYDp2uttx9/rt33YHTF5t5zDhZNtX3atGksXryYt15+mauvvNKq0AyX5ezJpxo0NfHw4cNMmXIeU6eexfXXX+31+z754luuvuE2JkwYw+LFJ3QM2ywvrxjgmJ6tnTW1B6Ndc52LN8Mg0cDrSql2NPTE39NaL1ZKbQA6Al86H1Zkaq1nKKVigDla66leXPt4k4FnlVJ1zu/vaC54OwvGJO1yfNuDqXiTUUnapX379tx7731cffWvmD79ck455RSv3md08SZ/SdItsHWuk93NheHWrVvHiBEjiImOZmtBgcxEaCOtNRdccAEXXJDO1Vdf6tV7zrvkOlZkr+Hdd19msnOLsGDmj7ubS9U9YTjXlL2J48ZJovYBpRQJCQnU1TX98O54RhVvEuaSZG2QjIyMZmd3BLLj2x5Mi2Eyli1rdraIL2ityc1dQ+fOYV6d7yreNGxYnE+LN7nLzMw9qdkXwnuSrIXhZDGMb+3evZs1a9Zw+eXTvDo/c+VqQIo3+TtJ1sJQ7sWbkkaOtDqcgNGhQwev66scXbkoi2H8mSRrYahgKN5kBW8nBmitWZHdMDyRliY9a38myVoYKpjGq80SERFBr169+c9/XmvxXCneFDgkWQtDyXi173Xo0IGPPvqIZ555scVzpXhT4JBkLQwjxZuMExMTQ2Rkbz7//Ntmz5PiTYFDBhENEhcXZ3UIlnG1feXKlRw6dIjEhAS6djVmypjdxLVxhxhvtW/fniuvvIrMzGzOO+9sj+eZtY1XbGxfQ68vpGdtmJiYmMYqecHG1fZgHK+OiYpqrCRotKuuuoq33vrA4893Vu5iQ+lmQkM7kZg43NBYoqJ6NlbmE8aQZC0M475yUfje4MGD6dIlgq+++r7Jn6/IbhgCSUlJokOHDmaGJgwgydogZWVllJWVWR2GJcrKyti6dWvQFG9yV+ZwUOajPRhbopTij3/8E5mZTdch/3FlDmDOYhiHoxKHo9Lw+wQzGbM2iGtzgGAcCikqKmLTpk3s3r2bvjExDOgfPFPGXJsNmDUU0q1bNxYvzqW+vp6QkGP7XmYuhnFtSCBDIcaRnrUwhGuzgQlpaTJlzEDXXnstRUUl7Nmz95jjdXUHWJO3HqWUFG8KEJKshSEak3UQPVy0Qrt27QgJCTnhF+Lq3HUcPHiIoUOH0LVrF4uiE74kyVoYwr1nLYx3fLI2a8qeMI8ka+Fzu3btYtu2bVK8ySRN1QmR4k2BR5K18Lm1zt3LpXiTObTWx/SspXhTYJJkLXwuLy8PkPFqq0jxpsAk3R6DpKenWx2CZbZs2QIE53h1+vjxpt5v9+7dHDhw8Jhpe1YUb0pLM3ajYCE9a+Fj+/fvZ9WqVYSEhDBuzBirwwl4b731Fj/72VlERIQ3HnMVb0pLk/HqQCLJWviUq3jTyBEj6NJFpowZ7eDBgwwYcOzCq6PbeEmyDiSSrA2SlZVFVlaW1WGYzlUPZMjgwRZHYo2s3Fyycs3bOHbv3r3HDIHs2LmLko1b6Nw5lJEjh5kWR15eMXl5xabdLxjJmLVBqqurrQ7BEq5kHR+kJWLN/HvXWvP888+zaNEbjcdcxZtGjUo0tXhTTc1+0+4VrKRnLXymvr6+sXhTwnBjS3KKhmRdd+AAw4Yd/cWYKYthApYka+EzBQUF7N69m169ehEVGWl1OAFPa83xcz1kMUzgkmQtfMY1BJIwfLgUbzLBkSNHjvlzluJNgU2StfAZ92QtjDd37lzS0lIbv5fiTYFNkrXwmSVLlgCQMGKExZEEh+3btzNhwtGNiKV4U2CT2SAGiY6OtjoEUzkcDkpKSggLC2NiWlrQ1gSJNmnTAfBcac+9t22WyMgept8z2ATn/ygTxMfHWx2CqVxDIGmjRzMiyNruLt6k3c2hIVm7Ku7V19e7rVw0v2c9cGA/0+8ZbGQYRPhEMO5kbjX3ZF1csonde/bSp08kAwZI4gxEkqwNUlVVRVVVldVhmKYxWaelUVVdTVWQLgoyq+319fV8+ukn9OvXsNTciuJN7qqra6murjX9vsFEkrVBsrOzyc5uetfpQONevClt9Giyc3PJNnHJtZ2Y1fbDhw+zfn0Bl156IXC0eJNVDxfz8zeQn7/BknsHC0nWos2keJP5QkJCjtl0QBbDBD5J1qLNZLzafN988w3du3cDrCveJMwlyVq0mft4tTDHV199xU03/RqwrniTMJcka9Em9fX1kqwt4prLLsWbgoMka9Em69evZ8+ePfSNiWFAf9nvzyxKHd3VXMarg4Mka9Em7r1qKd5kHqUaHjDur6sjJ3edFG8KArKC0SApKcHRy3El64luQyApicG7eapZbS8tLSUqaiSr16zj0KHDDB8eb2nxpoSE4NwZyEzSszZIREQEERERVodhuKbGqyPCw4kID/f0loBmRtv37NlDRkYGl1/+c9sUbwoP70x4eGdLYwh0kqzFSXMv3pSYkGB1OEGjrKyMTp06Eh4e5rYYJjg+yQUzSdYGKSwspLCw0OowDOVevMm9yl5hSQmFJSVWhWUpM9r+0ksvMX36ZdTX1zdO27OieJO70tKtlJZutTSGQCfJ2iDl5eWUl5dbHYahPC2GKXc4KHc4rAjJcka3va6uji+++JyUlCRbFW+qqNhFRcUuS2MIdJKsxUmT+dXmmzdvHl27RjBp0hmWF28S5pJkLU5KbW3tMcWbhPFqa2t56qknufXWGYD1xZuEuSRZi5MixZvMt23bNkJCFJMnnwXIYphgI8lanBQp3mS+ZcuWSfGmICbJWpwUGa823+LFH3PDDVcDR4s3paQkSfGmICHJ2iDh4eGEB+jCkPr6epYtWwbAxHHjTvh5ILe9JUa1ffPmzaxYsYLU1GQAMleuBuwzXh0WFkpYWKjVYQQ0WW5ukNRU83eYNoureFO/vn2bLN6UGsTLzY1q+5IlS4iNPZWYmIbd0+02Xj1y5BCrQwh40rMWrSbj1eaqra3lsccebaxfvb+ujpy89SilGD1aijcFC0nWotVkvNpcK1eupEuXcM4772yAxuJNw4bF0bVr4NefEQ1aTNZKqU5KqRVKqTVKqbVKqQecx2cqpQqUUrlKqYVKqW4e3r9JKZWnlMpRSmW5HY9RSn2jlPpIKRXuPBavlMpwnrteKfWSj9ppuoyMDDIyMqwOwxAtJeuMZcvIcI5pBxsj2j5v3luMGTOKkJCG/652Kd7kLjMzl8xM/94k2e65zpue9QHgHK11EpAMTFFKpQFfAgla60SgCLizmWucrbVO1lq7D+T+GbgZmANMdx6bBTztPHcYMNuL+ISJpHiTuQ4dOsSqVas577z0xmNSvMkwts51LSZr3aDa+W0H50trrb/QWh92Hs8EWlucoB1Q73y51spGA43VYLTWea28pjCYp+JNwhjz58/nyJFDjBo1EsBWxZsCjd1znVf/25RS7YBsYDDwb6318uNOuR5418PbNfCFUkoDL2qtXd3954A3gb3Ar5zHnga+UUotA74AXtVa72nqolVVVR6HGeLi4oiJiQEaykkWFRV5bFt6enrj11lZWVRXVzd5XnR0NPHx8Y33zs7O9nhN940HCgsLPRZ0Cg8PP2bWSHPDJnZoU0REBEuWLAGgb9++TX7cP37aWnNDAnGDBhET1TC7oczhoKiZanXp48c3fp2Vm+u5TVFRxA8aBEBVdTXZuZ4/mqckJjbWni4sKfFYgCk8PPyYWR7Ntamqurrxmm1tk9aaxV98wU1/+C2nnHIKAGvWFrJ7z1569uxBWdkuyst3H/OehITBjXWlS0u3eiyuFBYWeswMjuaGMGJj+xIV1RMAh6OSjRu3nXDO1q0n/tnl5RVTU7O/yWtGRvZg4MCGnFddXUt+/gaP9zeiTZ7YMde5ePWAUWt9RGudTMNvlDFKqcbPv0qpu4HDwDwPb5+gtR4FnA/8USl1pvOam7XWZ2qtp2mtq5zHXgWGAe8D6UCmUqqjNzEKc7h61gnDh1scSWCrr6/nhyVL0LqeUSlHf1FkrWpIQMOHx0vxJgPYOdcp16ab3lJK3QfUaK2fUEpdC8wAJmmta7147/1Atdb6CS/vlQ9cq7U+ocuXmpqqs7KymniXPbh6ye69XH9XW1tL165dqa+vZ/eWLR5rgrh6nu49x2Dhi7YfOXKE3/3uBtavX8vf/vYn0tMnNP7sD7fdy9sffMyjj97NjBnXtTVcn3H1YtPS/GOOfffuQ7KPG1c+gV1ynYs3s0F6u55+KqVCgZ8BBUqpKcDfgIs8Ba+UClNKRbi+BiYD+c3ca4pSqoPz6z5AT+DEz1zCEitXruTw4cMkJiRI8SaDZGVlccYZZ7Bt208sXPjaMYka7LcYJpDYPdd5M2YdDbzuHMsJAd7TWi9WSm0AOgJfOj+OZWqtZyilYoA5WuupQBSw0Pnz9sB/tdafNXOvycCzSqk65/d3aK23exGj7cTFxVkdgs95uxgmzjlmHIza0vbly5dz1VVX8te//p7LLptGp07Hfiqu2FFJ6aaf6Nw5lISEoW0N1adiY/taHYIv2DrXtZistda5wAnLpLTWTW5nrLUuA6Y6vy4Fklq6h9t7bwNu8/Z8O3M9DAwk3i6GcT00DEYn2/aDBw9y+eWX8eyz/+Tssyc0eY6dize5HkD6M7vnOlnBKLziXrxJVi761uHDh7nxxhtJSUnymKjBnothhHlkoqxBysrKgMDpYbdUvMldmXMKXDD2sE+m7Y888ghlZVt4443nmj3PzuPVDkclEBg9bLuSZG0Q1zzoQEnWrSne5JpbHIzJurVtX7VqFXPnvsLcuc+cMEbtzu7Fm1xzryVZG0eGQYRXXIthZAjEd+rr6/n73//OH/5wHaNGNT/lTYo3CUnWwitSaa9p33zzDf369aNHjx707NmTd959h5+2bm3xfevXr+dXv/oV9fUHufrqy1o83zUEkpZmvyEQYQ5J1qJF27dvp7S0VIo3HScjI4Prr7+eOXOepLBwGXl539GzR1fmz/+Aurq6Jt+Tn5/P9OnTOfvsdKKju/P667MJC+vc4r1kJ3MhyVq0SIo3Hauuro7rrruWa675Nc8//zhpaal06NCBzp1DmThxLD16dOfWW2/l4ot/Ts+ePfn00085dOgQN9zwWyZNmkS/fj1ZseJz7rvvdq8SdX19vfSshTxgFC1zJeuJMgQCwDXX/Jr27eGHHxbRteuxKzmVUpx5ZhrffvY5u3bt4OChQ/Tv358//emP7Ny5neXLP6Vbt66tul/Rho3s2buP6Ogo+vcPiMUn4iRIshYtkvHqow4dOsTXX39Dbm6Gx15xn6hIZs68j5Ej01nyww/U1NTw/vsfkJn5CXV1B/juu2UkJY3wOmm7D4FI8abgJcnaIIFSwKm2tpZVq1YREhJC2ujRXr0nkAs4zZ07lzFjRtG5c9M7eaeNalgGXlvbUBo0KSmJU08dwKxZj1BVVcNVV93EdkcFl1w8lVmz/unVPf1hMYy/FHDyZ5KsRbNcxZuSExOJiJApY/v27SMpaXiLPdz27dsBkJ2dzeTJk3n66ZfYvXs33bp1JyIinLvuusXre9p5MYwwjzxgFM2SncyPKigo4IUXnmf8+GYrawJwyimnMHv2I1x22aVccskveOKJJ3nxxZfYvn07s2b9kz59Ir26p52LNwlzSc/aIK5a2+47wfijkxmvznLuzuK+u0ogePzxx7nssgs580zPwzx5BZsAGDn0NCZPTqd9+/bccsufmTjxDCZNmsT48amtSrp2Lt7kLi+vGOCYnVqEb0nP2iDV1dUet57yFydbvCkQ2n68Q4cOsXDhQv7wh980e15NbR01tUfnWJ9zzkR++GER5eVbefrppwkNbXmqnjt/mbJXU7Pf4xZewjckWQuP1q1b53XxpkA3d+5ckpNHnjBVzxudOnXksssuYNeuXfzwQybLl6/y+r0yXi1cJFkLj2S8usGsWbO4/Y47+OMfm+9VN+f88yehFNx999387ne3sWVLy0vS3Ys3paYmn/S9RWCQZC08kvnVUFJSwqxZz/L2f5/n3HPPOunrhIV1Jjw8jHPOOYcLLriAl1/2tOfqUavWrOXQocMMHx4vxZuEJGvhWePKxXHjLI7EGlprbrvtNq644iImTmzbL6z8/AIqK3fx17/+ldraWhYu/B9bt5Y1+57lK3MAe8+vFuaR2SCiSa7iTeHh4YwcMcLqcCyxcuVKsrJW8vTT97T5Ws8++zLR0TF8/c03hHbqRO/evSgp2US/fp7rnUvxJuFOkrVBoqOjrQ6hTdpSvCk6ADYd0Frz0EMP8rvfTadHj+5evy+yV7cTjv34YxZLliznoYce4pNPPmLv3n1s3vwTp5/ueWpjfX1947Q9u88EAYiM7GF1CAFPkrVB4uPjrQ6hTdrycDE+AHY3f/TRRykvL+P66x9u1fsGDuhzzPcN+yv+lUceeYQBAwbw9dffs3Tp/+jduycdO57i8Tr+Vrxp4MB+VocQ8CRZiyYF88PFmTNn8sEH7/Huuy/RpUvbHuz9739fcuDAQf58yy106tiRCy6YTL9+LX/qkuJN4niSrA1SVVUF4Jf1NE6meJO7KueCmIjwcF+HZrjvvvuOl19+icWL36R379bvJ1jtLOAU7iz0dOGFk8nJWUdMTCQPPvRUY82QlvhD8SZ31dW1AISHt27Rj/CeJGuDZGdnA/5ZfW/FihVtKt6U7Vxu7m/V98rKyrj55j9x6603ERnZ+6SukV+wGThafa9du3bcd99fAXjnnY+or9deXefoykX/KFeQn78BkOp7RpKpe+IEwboY5te/ns6ll07ll7+8xJDrn332BOrr61s8z1W8KSyssxRvEo2kZy1OEIzj1bNnz2b79nJuuGG6Yff41a9+wbp1RS2e5+pVp6QkyTZqopH8SxDHqK+v58cffwSCJ1nPnj2bl156gffem+PVnogna+DAUxk48NQWz/O38WphDhkGEcdwFW/q369fUBRvciXqd955ib59+7T8BhMcnQli//nVwjzSsxbHCKbx6vnz5/P88//m/fdfsU2i3l9Xx5p8Kd4kTiQ9a3GMYBmv3r59Ow899BB//esM2yRqkOJNwjPpWRskJcU/P8L6Ilmn2HyHmMrKSn77298yZkwil146zafXThja8ph0c/y1eFNCwmCrQwh4kqwN4o+LYcrLy31SvMnui2H+9re/sXfvLh577G8+v3a4h13PveWvxZtkMYzxJFmLRm0p3uQP9u7dy4svvsjSpT+wYMGrREfbq+DUscWb/GMxjDBP4P2PtInCwkLAvwo6+erhYmFJCWCvgk7V1dVcd911lJdv5Z13XjQsUZdu2Q6cWNDJG67iTTExUfTv77l0qh2VljbsfCMFnYwjDxgNUl5eTnl5udVhtIqvHi6WOxyUOxy+CMknKisrmTz5XPbs2cnrr8/m1FONm5JYsXMPFTv3nNR73afs+VvxpoqKXVRU7LI6jIAmyVoADcWbVq9efdLFm+yqsrKSadMuZMKEVBYseJWoqJOr+WGGzJWrAf8brxbmkGQtgKPFmxITEvzy4WhT3BP1XXfdYvveqiyGEc2RZC0At/0WA2R+tb8l6oodlWzcLMWbhGeSrAUQWIthdu3a5VeJGqR4k2iZJGsRcMWbZs+eTZ8+vfwmUYMUbxItk1/hBgm3+cIQd+7Fm/r3a/vUKyvb/sYbbzBv3lu8/fYLliTqsM6dTup9/j5eHRbWtsVAomWSrA2Smuo/ixp8Xbwp1cTl5ps3b+bbb7/lyiuv5N133+XRRx/hnXde9KoUqRFGDj2t1e8JhOJNI0cOsTqEgCfJWrBkyRLA/4ZANm/ezJQp5xEREcb9999Hp04deeedF4mNHWB1aK3iKt6UkDBUijcJjyRZC799uHjLLX/msssu5LbbZlBeXkH37l0JDT25YQgrHS3e5J9DIMIc8oDRIBkZGWRkZFgdRovKy8vZuHFjm4s3uctYtoyMZct8cq2mVFZWcs8997B58yauv/6XKKWIiYmyRaLOXFVA5qqCVr3HX4s3ucvMzCUzM9fqMAKaJOsg50/Fm7TWPPXUU4wfP45Vq1bwzjsv0r17N6vDapP6+nq3mSDSsxae2ft/pzCcP+0M88AD9/PZZ5/wn/88TkqKvWtme6uweCN791X5ZfEmYS5J1kGuceXiuHEWR9K88vJy/v3v//D55+9aNtPDCO69an+ZEy6sIcMgQcy9eNNYm081nD17NuefPymgEjXIYhjhPelZBzFX8abTk5JsW7xJa80rr7zC888/z5dfvm91OD7n74thhHkkWQcxu49X19fX889//pMPP1zACy/MZPDgWKtD8ikp3iRaQ5K1QeLi4qwOoUVGLYaJ88EOMe+88w733XcvXbtG8O67L9GrVw8fRGa82FbsEOMaAklNTbb9TJyWxMb2tTqEgOff/0JsLCbG3k/2jSzeFBN18ltmaa355z//yfvvv8s999zGhRee68PIjBfVq5vX5wbSeHVUVE+rQwh4kqyD1Nq1a9m7d6/Pijf5ynPPPceHHy5g/vy5ftObPlmBsBhGmKfF2SBKqU5KqRVKqTVKqbVKqQecx2cqpQqUUrlKqYVKqW4e3t9NKfWB89z1SqlxzuMxSqlvlFIfKaXCncfilVIZSqkc57kv+bCtpiorK6OsrMzqMDwycry6zOGg7CT2YNyzZw/PPvsMDz/8N79N1I6de3B4sQdj7f795OStJyQkxG+LN7lzOCpxOCqtDqNN7J7rvJm6dwA4R2udBCQDU5RSacCXQILWOhEoAu708P5ngc+01kOBJGC98/ifgZuBOcB057FZwNNa62St9TBgthfx2VJRURFFRUVWh+GRkfVAikpKKHLucN4aBQUF9OrVg7Q0e08jbM7GLdvZ6NzhvDmr16zl8OHDDB8eR5cu9pyJ0xobN25j48ZtVofRVrbOdS0ma92g2vltB+dLa62/0Fofdh7PBE74LK2U6gKcCbzivNZBrfUe54/bAfXOl2s1QDSw1e3eeS3FJ06OHYs3vfzyS0yenG51GKbIlCl7tmP3XOfVmLVSqh2QDQwG/q21Xn7cKdcD7zbx1oHADuBVpVSS8xq3aK1rgOeAN4G9wK+c5z8NfKOUWgZ8Abzq1uBjVFVVeSyUFBcX1/iAr6ysrNkebnp6euPXWVlZVFdXN3ledHQ08fHxjffOzs72eM2UlKP/AQsLCykvL2/yvPDw8GPqXjdX+MmXbaqsrGTjxo2EhYU1Fm+qqq4mO9dzIZ6UxEQinJsKFJaUUO5hmOP4jQeaK+gUN2hQ48PIkk2bqDl8hL6xg5oshJQ26ujUtryCTdTU1jV5zche3RjonJFRXbuf/ILNHu+fMPRUwjs3FM0v3bKdCg/DF2GdOx1Tp7q5Qk01tXWNGxA4du7x2Mv+/OuGmTiu8eq8vGJqavY33abIHgwc2JAfqqtryc/f4LlNCYMJD+/c0KbSrVRU7Gq6TWGhx9Sgbq4IU2xs38YHiA5HZZM96K1bT/z3YOc2eWLHXOfi1QpGrfURrXUyDb9RxiilEtwadzdwGJjXxFvbA6OA57XWpwM1wN+d19ystT5Taz1Na13lPPYqMAx4H0gHMpVSHb2JUXgvPz8fgKSRI20zZaykpIR2IYoBp9rnYadR6uvrWV9QDEjP2m7snOuU1rpVjVFK3QfUaK2fUEpdC8wAJmmta5s4tw+QqbU+zfn9GcDftdYXeHmvfOBarfUJ3djU1FSdlZXVqtjN5Oolu/dy7eLWW2/lmWee4f477+S+Oz0Nv508V286ffx4r99zzTW/pl+/Xvzf//3J5/GYydXrdv8kcLz1hSWMP/cyYmKiyM//ISBqgrh6sWlp/lFgq3v3Idla62Yfjtgl17l4Mxukt+vpp1IqFPgZUKCUmgL8DbioqeABtNbbgZ+UUvHOQ5OAdc3ca4pSqoPz6z5AT8Dvn1rYjd3Gq4uKivjhhx+49torrQ7FFFK8yZ7snuu8+QwcDbzuHMsJAd7TWi9WSm0AOgJfOv/BZWqtZyilYoA5WuupzvffDMxTSp0ClAK/aeZek4FnlVKuAck7nH8IwkdqampYtWqVrYo3LV26lBEjhhIV1dvqUEzhStZpaTIEYjO2znWtHgaxC7sPg9jVt99+yznnnMPpSUms+uEHq8OhpqaGlJRR/POfdzJp0hlWh2OKUWdcxMbNP5GR8SFJSb7ZnUe0jjfDIHYjJVKDjN2KN7311lvExERx9tkTrA7FFI6KnY3Fm0aMiG/5DUI42WMqgDCNncarc3JyePzxx3jjjecICQmOfkMgFW8S5gqO/yEWyMrKwm7DNEYWb3KXlZtLVjNztl1+/PFHzjhjbECVB80r2ERewSaPPw+k4k3u8vKKycsrtjqMgCa/2g3iaXGNlcwq3uRN21evXs3jjz/Oyy8/aVgcVvC0WMclUIs3eVr8InxHetZBpHG/RZsMgUyefBajRydbHYppavfvZ01+QcAUbxLmkmQdROw0Xh2M84sDrXiTMJck6yBit2Ttr9NGT5YUbxJtIck6SJSVlbFx40YiIiIaizdZRWvNihUrCA0NtTQOs8liGNEW8oAxSLh61WmpqbRr187SWBYsWMCKFZm8//7LlsZhpvr6elZIz1q0gSRrg0RHR1sdwjHMHAKJbmEPxkcffZSbb/4N3bp1NTwWs0V62IOxsHgje/dVERPTh/797b0/58mIjPTPnX38iSRrg7hqX9uFmck6voXdzcvKtnH22YG5tHygh93NM1euBgK3V+2qTS2MI2PWQaCmpobVq1fbpnjTkSP1tG9v7VCM2Y6OVwfW/GphHknWBqmqqqKqqsrqMABYsWIFR44cIWnkSCIijJ8yVlVdTZWHhTEff/wxffpE0rHjKYbHYYXq2v1U1564QGR5gI9XV1fXUl3dZPVQ4SOSrA2SnZ3d7NZfZjK7eFN2bq7HLcIyMzM577yzLX/IaZT8gs0nbCXmqNjJpi1bA7p4U37+hma35hJtJ2PWQcAu86t3797Nf/7zH7766gNL4zCbFG8SviA96wDnXrxp4rhxlsby4osvcsEFkxk06DRL4zBboBZvEuaSZB3gXMWbBvTvT7++fS2LQ2vNV199ybhxgTlm2xzXeLUshhFtIck6wC1ZsgSwfrOB+vp6Vq7M4pe/vMTSOMwmxZuEr0iyDnB2Ga8+fPgw9VoHzSYDLq7iTSNGxBMREW51OMKPBdf/nCBkl2T9yCOPMOmcwFwI0xwp3iR8RR5NGyQlxfr/nGVlZWzatMn04k0piYnHfL9v3z7mzJnDe++9ZFoMVkkYeuox3x9duRjYDxcTEgZbHULAk2RtEDMWn7TEquJNEeHHftx/4omZTJgwhhEjAmf7Lk/COx+tJFhfX8/K7Ib55oHesw4P72x1CAFPhkECmB2GQFzlUKdMOduyGKxSUFwa0MWbhLkkWRuksLCQwsJCS2OwKlkXlpRQWFICwPbt21m67EfGjbO+JokZSrdsp3TLdgCWr8wBAr9XDVBaupXS0q1WhxHQJFkbpLy8nPLycsvub2XxpnKHg3KHA4A+ffpw/pTz+PbbpabGYJWKnXuo2LkHCK7iTRUVu6io2GV1GAFNknWAchVvSk5MtHT8XClFVFQU9fX1lsVglUAv3iTMJck6QJldvKk5e/fuC7oNcl3Fm8LDwwK2eJMwlyTrANW4ctHi+dU7duzgq6++5Nxzz7Q0DrMdLd6UJMWbhE9Isg5AR44caSzeZHWy/u9//8vAgacRGdnb0jjMdrR4kwyBCN+QX/kBaO3atezbt8/y4k1z5rzMnDkv8fbbL1oWg1WOjlcH/sNFYQ5J1gYJD7euDoTV49Xh4eEcqKtj5syZPPbYP+jXz16bBxsprHMn6uoOBF3xprCw0JZPEm0iwyAGSU1NJdWi/Q6tXgyTmphIB6Br1y6kp4+3JAarjBx6Ggf2VwVd8aaRI4cwcuQQq8MIaJKsA5DVyRrghRee5+yzJwTdLBCATBmvFgaQZB1grCredLwdO3YwceJoy+5vJdkZRhhBkrVBMjIyyMjIMP2+VhVvcpexbBnRp51myb2ttixrHZkrcoDg6llnZuaSmdn0JsnCNyRZBxhXsrZyv8UdO3awo2IHsbGnWRaDVbb8tI3qmhr69o2W4k3CpyRZBxg7jFfvq6oiMrJXUM0CcVlX0FC8S4ZAhK9Jsg4g1dXVrF69mnbt2plevMldiFJorS27v5XWrS8CgmsIRJhDknUAcRVvSho50tJ53iqIk/Xa9Q0962CotCfMJck6gFi9GMalIVlbGoIltjt2sN1RQWhoJ4YPl+JNwrckWQcQO4xXA6iQkKDsWbum7A0dGifFm4TPyb8og8TFxZl6PzsVb/pp40YO1OyzNAYruOqBjB8ffPPLY2Otq0ETLCRZGyQmxtxpW3Yp3gTw6eKPueeev1gagxVcPetJkyZaG4gFoqJ6Wh1CwJNhkABhl/HqzMxMKit3Mny4uZ8srFa7fz+5awuDqniTMJcka4OUlZVRVlZm2v3sMl69a9cuRqWOovbAYUvjMNuqnLUcPnyY+PjB1NYesDoc0zkclTgclVaHEdAkWRukqKiIoqIi0+5nh5WLAFpruvbqzUbnDt/BwlW8afDgQWzcuM3aYCywceO2oGy3mSRZB4Bt27axadMmunTpQsLw4ZbG0r59e/btq+JI/RFL4zCba7x6+PCh1gYiApYk6wBgh+JNLpMmTaL+iMbh2GlpHGaqr69nZXZDEaMRIyRZC2NIsg4AdhmvhoaedWhoKLq+3upQTFNQXMrefVX07RtN7969rA5HBChJ1gHATskaoGfPnhQVlVodhmkyV6wGpHiTMJYkaz9XXV1NTk6O5cWb3A0fPpzKXbutDsM0RzfHleJNwjiSrP2cXYo3uQsNDWXnzl2UlzusDsUUroeLUrxJGElWMBokPT3dlPvYZTGMuysvuYSP3n+Pl1+ex7333mZ1OIba7tjB5p+2ER4exvDh8UFbEyQtLdHqEAJeiz1rpVQnpdQKpdQapdRapdQDzuMzlVIFSqlcpdRCpVS3Jt4br5TKcXvtU0r9xfmzGKXUN0qpj5RS4W7nZzjPXa+Uesm3zQ08dhuvdvnLX25l6dLlVodhOFevOjU1KWgTdaCwe67zZhjkAHCO1joJSAamKKXSgC+BBK11IlAE3Hn8G7XWhVrrZK11MpAC1AILnT/+M3AzMAeY7jw2C3ja+Z5hwGwv4gtadiredLzhw4dTXb2fLOd4bqByjVenpdnjeYFoE1vnuhaTtW5Q7fy2g/OltdZfaK1da4ozgX4tXGoSUKK13uz8vh1Q73wp57FoYKvbvfNais+usrKyyMrKMvQeruJNpw4YYHnxJndZubkUlJYyZsxovv/+R6vDMdTxO5nn5RWTl1dsYUTWCIR22z3XefWAUSnVTimVA1QAX2qtj/98ez3waQuXuQp42+3754AXgRnAW85jTwPfKKU+VUrd2tTHDX9RXV1NdXV1yye2wZIlSwB7jVfD0bbPmPF7PvvsW6vDMYx78aaUlCQAamr2U1Oz3+LIzBco7bZzrvNqkE1rfQRIdl5woVIqQWudD6CUuhs4DMzz9H6l1CnARbh9fHD+1jnzuPu8qpT6HJgC/By4SSmVpLU+oTJOVVUVGRkZTd4vLi6usURpWVlZszU63B8EZmVleUyw0dHRxMfHN947Ozvb4zVTUo5O4SosLKS8vLzJ88LDw0l1m27nqT3QdJsWLFgAQO/ISDKWLTvapvHjj7YpN9dzm6KiiB80qKFN1dVk5+Z6blNiIhHO2SaFJSWUO5qe6eE+IyUuLo7+g4fw9oIviD1twAnnxg7oQ1SvbgA4du5ptp5I2qijKwPzCjZRU1vX5HmRvboxcEAfAKpr95NfsLnJ8wAShp5KeOdQAEq3bKdi554mzwvr3ImRQ09r/D5zVQEAuXnrOHz4MIMGxbJ2bcO88pqa/YSFNVzT4ahstl6G+0O5vLxij8kuMrIHAwc2dOaqq2vJz9/guU0JgwkP79zQptKtVFTsarpNYaGMHDnkaJsyPf/dx8b2bSyB6qlNW7ee+O/Bzm3yxI65zqVVU/e01nuADOcNUEpdC1wIXK2b3xrkfGCV1rrFuVxa6zKt9Vyt9c9p+INJaE2MwSQ/Px+AhBEjLI6kad27d2fihIl8992P1NTUWB2Oz6117mQu9UACjx1znWpp+yWlVG/gkNZ6j1IqFPgCeNx58aeAs7TWO1q4xjvA51rrV1s4bwrwtdb6kFKqD7AaOF1rfUKXKzU1VRs9JtwWrl6yUVP4tm3bRr9+/ejSpQu7Nm+2vCaIO1cv39XDf+yxx/jkk49ZtOgNK8Pyucuu+SNfZyxjzpynufTSC4Gjvblgm8rmb+3u3n1Ittb6mKfCds11Lt70rKOBb5VSucBKGsZxFtMwDhMBfOmcfvKCM4gYpdQnbkF1Bs4FFnhxr8lAvlJqDfA5cEdzwQczOxVvaskNN9xAcXEJ+/ZVWR2Kz9TX17NyVUOCkpWLAcPWua7FMWutdS5wehPHB3s4vwyY6vZ9LeDVnj9a69uAwF5F4SN2nV/dlF69enHWWeksWLCY6677pdXh+MT6ohL27aumb99o+vWLtjoc4QN2z3Uyi98g0dHG/ge2c7KOjoo64diAAQM4fDhwalwvX5kDnFi8KTKyhwXRWC9Y220mSdYGcc0cMYIdize5c80wcae1RqkmTvZTnoo3uWY4BJtgbbeZpJCTH7Jj8aZg41oMM26cjFcLc0iyNkhVVRVVVcY8UGvcb9GGQyDQMGe76ri53Q0968DoWruKN0VENBRvclddXUt1da1FkVknWNttJknWBsnOzm524UxbNK5ctGmyzs7NPWGBTSAl66PFm5JPmImTn7+h2QUegSpY220mSdZ+xs7Fm5rTuXNnVq9ea3UYPiGbDQgrSLL2M/n5+VRVVXHqgAH0dS4/9we33347Cz/8hJYWYfmDzJWyjZcwnyRrP2PHzQa80aVLF1JTUpg3b77VobRJTe2JxZuEMIMkaz9j5/nVzVFKkZ6eTmVl00V4/MWqnHyOHDlCQsJQIiJkJo4wjyRrP+OvyRogJiaGDz5Y7NdDIUfrV8t4tTCXJGs/sm3bNjZv3kyXLl1IGD7c6nBaraKigsGDT/PrWSGZx202IIRZZAWjQdxrWvvKMmc1u5CQEO55+GEmjB3LuDFj6NHDXkt9UxKbrrw2ffp0nnrqKZOj8R334k1paU3//SYkNFlGIuAFa7vNJMnaIBERET6/Zt++fTnllFPYs2cPjz75ZOPx4UOHMiEtjQljxzIhLY1BAwda2nuN8LCqsn///vTpE8WiRZ9z0UXnmRxV27mKN/XrF0Pfvk3XfnEVyQ82wdpuM0my9iPjx4+nrKyMZcuWsXTpUpYuXcrKlStZV1DAuoICXn7tNQCiIiMZ70zcE9PSOD0piVNOOcXa4Gn4RHDjjTeRm7vWL5O1p+JNQphBkrVBCgsbdhHxdUGnnj17Mm3aNKZNmwbAgQMHWLVqVWPyXrp0KY6KChZ+/DELP/4YgE6dOjEmJaWx9z1+7Fi6d+/u07jcFZaUAJ4LOvkrbxbDlJY27IEabIWNgrXdZpJkbRDXvotGVt8D6NixI+PGjWPcuHHcfvvtaK3ZsGEDS5cuZcmSJSxdupSCggK+X7qU750zSeDo0MnEtDQmpKUxMDbWZ0Mnrv0Zm0rWDfzzAaNrJkhamueetWt/wGBLWsHabjNJsg4wSimGDBnCkCFDuO666wCorKz0auhkgjNxTxg71rChE3+tEVK+vcJj8SYhzCDJOgh4Gjpx9bxdQycLFi1iwaJFgHFDJ/6arF1DIE0VbxLCDJKsg5D70Mkdd9yB1pri4uJjxr2bGjoZMWzYMbNOTnboxB/HrWUxjLCaJGuBUoq4uDji4uL4zW9+A8DOnTv58ccfjxk6Wbt+PWvXr+elVxs2bj6ZoZOxY8dy9dXP8ve/3+xXPdTlshhGWEyStWhSr169Thg6yc7OPmHWifvQSWhoKGNSUujXrx8Jw4eTNGzYCUMnEyZM4MCBgxQXb2ToUP9YSCHFm4QdSLI2SKBtt9WxY0fGjx/P+PHjmx06+c65MQLAnffe2+TQyYwZM3j77YU88MAdFrbIe67iTUlJI1os3hQWFmpSVPYSrO02kyRrg6TacCNbX/I0dOI+6yQrK+uEoZM+UVFE9+5NWGgHsnPySRwRT4cOHaxsSotaMwQycuQQg6Oxp2Btt5kkWQuf6dWrFxdddBEXXXQR0PTQyXaHg+3Oedg/u+jXhHbqxKjkEYxNTWZsahJjRiXRrVsXK5txgkx5uChsQPnjk3mA1NRUnZWVZXUYohVcQyf33XcfOTnZHDlyhOLi0mPOUUoxNG4QY1OTGJuaTFpqMqcO6GvZdL/6+npiE89i375q8vO/91gTRPiX7t2HZGut/erjr/SsDZKRkQFAenq6pXFYwVPbXUMnY8eOpVu3Tjz66N1UVu5i+fLVLF+ezfLlq1i9Opf1hRtYX7iB15y7ykT17tWQvEcnMzY12dShE2+KN7nLzHRV5Wu68mCgCtZ2m0mStTCd+6e5nj17MHXqJKZOnQQ0DJ3k5KxtTN7Ll2fj2LGTRZ9+zaJPvwY4ZugkbXQyY0Yl0bWr76scguy3KOxDkrWwxIEDB5s83rFjR8aOHdWYHBtqnWx0S96rKC4uZWlmNkszswFjh05kJ3NhF5KshenOO+887rrrLv7v//5Inz6RzZ7bUOtkIEOGDGT69MsBWjV0kjb6dMaOTmbk8LiTGjrxpniTEGaQZC1MN3z4cK688gpmzZrDI4/c1er3Hz90Uld3gJyc/GN638cPnXQO7cSo5IRjZp20NHRSvr2CLT+VSfEmYQuSrIUl/v73O0lLG8ttt91Er14923StTp06kpaW0rjVlqehkyU/ZrHkx4YZRN4MnbiGQEaPPt2vlsaLwCTJWlhi6NChRER0Yf/+Az6/dnNDJ5mZWSxfvoqcnLwThk76RPZy9ryTGTs6maXLG8bE5eGisANJ1gaJi4uzOgTLeNP2zz//nHbtFL16mbPZrzdDJ9srdvLRJ1/x0SdfHfPePn2ivL5PbGxfn8btL4K13WaSZG2QmJgYq0OwjDdt37JlC2eeOY7Q0E4mRHSiloZOMjKWsXVrGQDnnXe219eNimrbkI6/CtZ2m0mStbDEiBEj+Mc/7mbPnr1069bV6nCaHDopLNzABx98TFRUb4ujEwJCrA4gUJWVlVFWVmZ1GJbwpu3jx4+na9euFBaWmBRV68XHD+buu29t1XscjkocjkqDIrKvYG23mSRZG6SoqIiioiKrw7CEt22/5ppree21d02IyDwbN25j48ZtVodhumBtt5kkWQvLnHXWWSxZksm+fVVWhyKE7UmyFpY544wzOPXU03jzzfetDkUI25NkLSz10EMPM3v2K+zevcfqUISwNUnWwlJnn302AwYMoKBgg9WhCGFrkqyF5Tp37sy6dYVWhyGErUmyFpa79977mD37Faqra6wORQjbkm29hOW01owdO4a7776Zs84ab3U4Igj447Ze0rMWllNKMWnSJN54Q2aFCOGJJGthCxdccCFbt5ZbHYYQtiXJ2iBZWVkE6zDNybRda01IiP//c8zLKyYvr9jqMEwXrO02kxRyMkh1dbXVIVjmZNqutfbJnolWq6nZb3UIlgjWdpvJ/7syIiD06dOH4uISVq3KtToUIWxJkrWwhbi4OJ588inuuONBq0MRwpYkWQvbOOOMM9i0aQslJZusDkUI25FkLWwjNjaWyZMns3LlaqtDEcJ2JFkL2wmEWSFC+JrMBjFIdHS01SFYpi1tDw0NJT+/kCuu8GFAJoqMNGcDYLsJ1nabqcXl5kqpTsD3QEcakvsHWuv7lFIzgWnAQaAE+I3Wek8T778VuAHQQJ7zvDqlVAzwFlAFXK21rlZKxQMvAt2c9/tBa31jU3HJcvPAlJ+fz9SpU8jNzbA6FBHAmlpubtdc5+LN580DwDla6yQgGZiilEoDvgQStNaJQBFwZxPB9wX+DKRqrROAdsBVzh//GbgZmANMdx6bBTyttU7WWg8DZnsRnwggr732KuPG+VXJBhE4bJ3rWkzWuoFrlUMH50trrb/QWh92Hs8E+nm4RHsgVCnVHugMuHZSbQfUO1+u1RDRwFa3e+e1FJ9dVVVVUVUVnNtVtaXt27ZtY9KkM3wckXmqq2uprq61OgzTBUK77Z7rvBqzVkq1A7KBwcC/tdbLjzvleuCEnU+11tuUUk8AW4D9wBda6y+cP34OeBPYC/zKeexp4Bul1DLgC+DVpj5uQENCyMjIaDLeuLg4YmJigIadtpvbvDU9Pb3x66ysLI+r76Kjo4mPj2+8d3Z2tsdrpqSkNP48Ojqa8vKma16Eh4eTmnq0F+mpPWCPNkVERABQWFjYbJtc90tPT291m3r16kVV1QEyM49dHJOWltj4dV5esccVc5GRPRg4sOH/UnV1Lfn5njc1SEgYTHh4ZwBKS7dSUbGryfPCwkIZOXJI4/fHx+aupmY/YWGhpKUl4nBUNruJrL+0KTa2L1FRPQE8tmnrVgcAl112rl+0yRM75joXrx67a62PaK2TafiNMkYpleDWuLuBw8C849+nlOoO/ByIBWKAMKXUdOc1N2utz9RaT9NaVzmPvQoMA94H0oFMpVRHb2IUgaHhEYr/LzsX/snOua7V9ayVUvcBNVrrJ5RS1wIzgEla6xM+AymlLgemaK1/6/z+GiBNa/0HL++VD1yrtT6hy2f3B4yuHqV7LzdYtKXtV155BeeeO47LLpvm26BM4urNufeag4G/tdubetZ2yXUuLfaslVK9lVLdnF+HAj8DCpRSU4C/ARc1FbzTFiBNKdVZNVTpmQSsb+ZeU5RSHZxf9wF6Ap4/R4qA46+bYQj/Z/dc582YdTTwunMsJwR4T2u9WCm1gYYpJ186q6Vlaq1nOKepzNFaT9VaL1dKfQCsouHjw2rgpWbuNRl4VilV5/z+Dq31di9iFAFEtvcSFrF1rmsxWWutc4HTmzg+2MP5ZcBUt+/vA+5r6T7Oc28DbvPmXBGYrr56Onfc8Veuu+6qlk8WwofsnutkXa+wleTkZOrq6lo+UYggI8vNDZKSkmJ1CJYJ5rYnJDTZCQt4wdpuM0myNohrTnIwCua2u+YDB5tgbbeZZBhECCH8gCRrgxQWFlJYWGh1GJYI5raXlm6ltHRryycGmGBtt5kkWRukvLzc45LsQNeWtnfo0IF9+6rYsaPSx1GZo6Jil8fl0IEsWNttJknWwlZiYmI455xJLF78RcsnCxFEJFkL2+nevRvt2rWzOgwhbEWStbCdI0eOsGvXbqvDEMJWJFkL27nxxpt4/vnXrQ5DCFuRZC1sZ/jw4dTU1HDgwAGrQxHCNiRZGyQ8PJzw8HCrw7BEW9vetWtX0tLSmDdvvg+jMkdYWChhYaFWh2G6YG23mWQFo0Hcd4AJNm1tu1KKkSMT2b/f/2qEuO9UEkyCtd1mkp61EEL4AUnWQgjhByRZGyQjI6PZzWIDWTC3PTMz96Q2avV3wdpuM0myFrbUsWNH1q8vtjoMIWxDkrWwpTvuuIN58z5gyxYpDiQESLIWNtWzZ09uvPFGHn/8OatDEcIWJFkL27rvvvv473/ns3VrcFYvFMKdJGthWzExMdxwww3Mm/eB1aEIYTlJ1sLWJk6cyKJFn3PkyBGrQxHCUrKC0SBxcXFWh2AZX7b9mmuu4c9/vhmHYwcxMX18dl2jxMb2tToESwRru80kPWuDxMTEEBMTY3UYlvBl25VSXHHFlbzwwhs+uZ7RoqJ6EhXV0+owTBes7TaTJGthe2eeeSYZGUs5ePCg1aEIYRlJ1gYpKyujrKzM6jAs4eu2T58+ndraOubPX+yzaxrF4ajE4fDP/SPbIljbbSYZszZIUVERQFAOhfi67Uop5s9fwPjx4znrrPG2HrveuHEbQNANCQRru80kPWvhF5KSkrjsskt59tmXrQ5FCEtIshZ+Y/DgIezZs9fqMISwhCRr4TfWrVtL9+7drA5DCEtIshZ+48477+Lddz/E4dhhdShCmE6StfAbiYmJ/OxnP+O5516xOhQhTCfJWviV6677DUVFG60OQwjTKa211TGcFKVUltUxCCH81k6t9RSrg2gNv03WQggRTGQYRAgh/IAkayGE8AOSrIUQwg9Isj6OUmquUqpCKZXvduxypdRapVS9UirV7fhpSqn9Sqkc5+sFD9e8Xym1ze28qW4/m6mUylJKneXhmjlKqWuMbLNbLD5vu/Pcm5VShc7r/MvteEC3XSn1rts5m5RSOW4/s0XbDWp3slIq03lOllJqjNvPbNFuv6S1lpfbCzgTGAXkux0bBsQDGUCq2/HT3M9r5pr3A7c3cXwoMBPoDLzXmmv6UdvPBr4COjq/jwyWth93/SeBe+3WdoP+zr8Aznd+PRXIsFu7/fElVfeOo7X+Xil12nHH1kND9TcfawfUAxrw+cVby6C2/x54TGt9wHm9CufxYGg7zvcr4ArgHOch27TdoHZroIvz666Aq16ubdrtj2QYpO1ilVKrlVLfKaXOaOa8Pymlcp0fO7sDaK3X0tDLWAI873buoOM+FjZ3XSt50/Y44Ayl1HLneaMhaNrucgbg0FoXg9+33Zt2/wWYqZT6CXgCuBP8vt2Wk55125QDA7TWlUqpFOBDpdQIrfW+4857HniIhh7FQzR8JL4eQGt9cxPXLdFaJxsXtk942/b2QHcgDRgNvKeUGqgbBHrbXX4JvO1+wE/b7m27fw/cqrWer5S6AngF+Bn4bbttQXrWbaC1PqC1rnR+nQ2U0NCTPP48h9b6iNa6HngZGHP8Of7G27YDW4EFzuS8goaPwb3Mi9T3WtF2lFLtgV8A75oXoTFa0e5rgQXOr98nAP6924Ek6zZQSvVWSrVzfj0QGAKUNnFetNu3lwD5x5/jb7xtO/AhzrFapVQccAqw06QwDdGKtkNDj7JAa73VrPiM0op2lwFnOb8+Byg2J8IAZ/UTTru9aPi4Wg4coqFX+FsaEuxW4ADgAD53nnspsBZYA6wCprldZw7OJ+nAm0AekAssAqKbuf9pwH4gx+31Zz9u+ynAWzT8gloFnBMsbXd+/xoww4v7W9J2g/7OJwLZzvOWAyl2a7c/vqQ2iBBC+AEZBhFCCD8gyVoIIfyAJGshhPADkqyFEMIPSLIWQgg/IMlaCCH8gCRrIYTwA5KshRDCD/w/n8O+F0kv6gwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x509.002 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_lon, min_lat, max_lon, max_lat = grid.bbox(buffer=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set simulation timeframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'_template': '../rompy/templates/swan',\n",
       " '_checkout': None,\n",
       " 'model': 'SWAN',\n",
       " 'compute_start': '20201214.000000',\n",
       " 'compute_stop': '20201215.000000',\n",
       " 'out_start': '20201214.000000'}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "startdt = pd.to_datetime('2020-12-14 00:00')\n",
    "enddt = startdt + pd.to_timedelta(1.,'d')\n",
    "\n",
    "new_simulation.settings['compute_start'] = startdt.strftime('%Y%m%d.%H%M%S')\n",
    "new_simulation.settings['compute_stop'] = enddt.strftime('%Y%m%d.%H%M%S')\n",
    "new_simulation.settings['out_start'] = new_simulation.settings['compute_start']\n",
    "new_simulation.settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gather boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Winds\n",
    "\n",
    "Use the winds in the BoM Wavewatch model available on NCI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_wind = rompy.cat.bom.wavewatch3_nci(startdt=startdt,\n",
    "                                       enddt=enddt,\n",
    "                                       fmt_fields={'fcdate':[startdt.strftime('%Y%m%d')],\n",
    "                                                   'hr':['0000'],\n",
    "                                                   'grid':['PER'],\n",
    "                                                   'output':['msh']},\n",
    "                                       ds_filters={'subset':{'data_vars':['longitude','latitude','uwnd','vwnd']}},\n",
    "                                       hindcast=True\n",
    "                                      ).to_dask()\n",
    "\n",
    "# Crop to the region of the model domain and load data to memory for later interpolation\n",
    "ds_t0 = ds_wind.isel(time=0,drop=True)\n",
    "inds = (ds_t0.longitude > min_lon) & (ds_t0.longitude < max_lon) & (ds_t0.latitude > min_lat) & (ds_t0.latitude < max_lat)\n",
    "ds_wind = ds_wind.isel(node=inds).load()\n",
    "\n",
    "ds_wind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup a regular swan grid to interpolate the unstructured wind mesh onto"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "uwnd: 100%|██████████| 25/25 [00:07<00:00,  3.54it/s]\n",
      "vwnd: 100%|██████████| 25/25 [00:06<00:00,  3.64it/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (m: 257, n: 167, time: 25)\n",
       "Coordinates:\n",
       "    init     (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead     (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23 24\n",
       "  * time     (time) datetime64[ns] 2020-12-14 2020-12-14T01:00:00 ... 2020-12-15\n",
       "Dimensions without coordinates: m, n\n",
       "Data variables:\n",
       "    lon      (m, n) float64 115.5 115.5 115.5 115.5 ... 115.8 115.8 115.8 115.8\n",
       "    lat      (m, n) float64 -32.81 -32.81 -32.81 -32.81 ... -32.3 -32.3 -32.3\n",
       "    u10      (time, m, n) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
       "    v10      (time, m, n) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0</pre>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (m: 257, n: 167, time: 25)\n",
       "Coordinates:\n",
       "    init     (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead     (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23 24\n",
       "  * time     (time) datetime64[ns] 2020-12-14 2020-12-14T01:00:00 ... 2020-12-15\n",
       "Dimensions without coordinates: m, n\n",
       "Data variables:\n",
       "    lon      (m, n) float64 115.5 115.5 115.5 115.5 ... 115.8 115.8 115.8 115.8\n",
       "    lat      (m, n) float64 -32.81 -32.81 -32.81 -32.81 ... -32.3 -32.3 -32.3\n",
       "    u10      (time, m, n) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
       "    v10      (time, m, n) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tqdm import trange\n",
    "from scipy.interpolate import LinearNDInterpolator\n",
    "import xarray as xr\n",
    "\n",
    "dx = 0.002\n",
    "wind_grid = SwanGrid(gridtype='REG',\n",
    "                     x0=min_lon,\n",
    "                     y0=min_lat,\n",
    "                     dx=dx,\n",
    "                     dy=dx,\n",
    "                     nx=int((max_lon-min_lon)/dx),\n",
    "                     ny=int((max_lat-min_lat)/dx),\n",
    "                     rot=0\n",
    "                    )\n",
    "\n",
    "interp_vars = {'uwnd':[],'vwnd':[]}\n",
    "for var,values in interp_vars.items():\n",
    "    for t in trange(len(ds_wind.time),desc=var):\n",
    "        ds_t = ds_wind.isel(time=t)\n",
    "        interp = LinearNDInterpolator(list(zip(ds_t.longitude, ds_t.latitude)), ds_t[var], fill_value=0)\n",
    "        values.append(interp(wind_grid.x,wind_grid.y))\n",
    "\n",
    "ds_wind_interpolated = xr.Dataset({'time':ds_wind.time,\n",
    "                                   'lon':(('m','n'),wind_grid.x),\n",
    "                                   'lat':(('m','n'),wind_grid.y),\n",
    "                                   'u10':(('time','m','n'),np.array(interp_vars['uwnd'])),\n",
    "                                   'v10':(('time','m','n'),np.array(interp_vars['vwnd'])),\n",
    "                                  })\n",
    "ds_wind_interpolated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the interpolated winds to swan format, this function returns the corresponding swan commands for including in the simulation INPUT file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "REG 115.48384449028222 -32.809999999999995 0 166 256 0.002 0.002 NONSTATION 20201214.000000 1.0 HR\n",
      "1 'wind.inp' 3 0 1 0 FREE\n"
     ]
    }
   ],
   "source": [
    "inpwind, readwind = ds_wind_interpolated.swan.to_inpgrid(new_simulation.staging_dir + '/wind.inp',\n",
    "                                                          var='WIND',\n",
    "                                                          grid=wind_grid,\n",
    "                                                          z1='u10',z2='v10')\n",
    "print(inpwind)\n",
    "print(readwind)\n",
    "new_simulation.settings['wind_grid'] = inpwind\n",
    "new_simulation.settings['wind_read'] = readwind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wave Spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.Dataset&gt;\n",
       "Dimensions:       (direction: 24, frequency: 29, station: 167, time: 25)\n",
       "Coordinates:\n",
       "  * direction     (direction) float32 7.5 22.5 37.5 52.5 ... 322.5 337.5 352.5\n",
       "  * frequency     (frequency) float32 0.035 0.0385 0.04235 ... 0.4589 0.5047\n",
       "  * station       (station) int32 1 2 3 4 5 6 7 ... 161 162 163 164 165 166 167\n",
       "  * time          (time) datetime64[ns] 2020-12-14 ... 2020-12-15\n",
       "    init          (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead          (time) int64 0 1 2 3 4 5 6 7 8 ... 16 17 18 19 20 21 22 23 24\n",
       "Data variables:\n",
       "    frequency1    (frequency, time) float32 dask.array&lt;chunksize=(29, 25), meta=np.ndarray&gt;\n",
       "    frequency2    (frequency, time) float32 dask.array&lt;chunksize=(29, 25), meta=np.ndarray&gt;\n",
       "    station_name  (station, time) |S64 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    crs           (time) int32 -2147483647 -2147483647 ... -2147483647\n",
       "    cur           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    curdir        (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    mean_dir      (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    dp            (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    dp_sea        (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    dp_sw         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    dpt           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    efth          (station, frequency, direction, time) float32 dask.array&lt;chunksize=(167, 29, 24, 25), meta=np.ndarray&gt;\n",
       "    hs            (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    hs_sea        (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    hs_sw         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    latitude      (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    longitude     (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdir0         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdir1         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdir2         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdp0          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdp1          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    pdp2          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    phs0          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    phs1          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    phs2          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ppe0          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ppe1          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ppe2          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    psig1         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    psig2         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ptp0          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ptp1          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    ptp2          (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    swe           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    t01           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    t02           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    tp            (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    tp4           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    tp_sea        (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    tp_sw         (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    cd            (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    wnd           (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    wnddir        (station, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "Attributes:\n",
       "    product_name:                    per.2020_src.nc\n",
       "    history:                         Mon Dec 14 14:01:41 2020: ncks --overwri...\n",
       "    NCO:                             netCDF Operators version 4.7.7 (Homepage...\n",
       "    date_created:                    2020-12-14T14:02:20Z\n",
       "    program:                         WAVEWATCH_III(R)_version_number_6.07\n",
       "    Conventions:                     CF-1.6, ACDD-1.3\n",
       "    id:                              f1762_5247_9737_2398\n",
       "    naming_authority:                au.org.nci\n",
       "    title:                           AUSWAVE\n",
       "    date_modified:                   2020-12-14T14:02:20Z\n",
       "    date_issued:                     2020-12-14T14:02:20Z\n",
       "    references:                      http://www.bom.gov.au/research/\n",
       "    source:                          Perth triangular\n",
       "    comment:                         WAVEWATCH_III_switches: O0 O1 O2 IS0 IC0...\n",
       "    summary:                         (C) 2020. Research version of marine for...\n",
       "    standard_name_vocabulary:        CF Standard Name Table v47\n",
       "    keywords_vocabulary:             NASA Global Change Master Directory (GCM...\n",
       "    keywords:                        Earth Science &gt; Oceans &gt; Ocean Waves\n",
       "    institution:                     Australian Bureau of Meteorology\n",
       "    project:                         ep4\n",
       "    creator_name:                    Australian Bureau of Meteorology\n",
       "    creator_url:                     http://www.bom.gov.au\n",
       "    creator_email:                   oceanservices@bom.gov.au\n",
       "    contributor_name:                Stefan Zieger\n",
       "    contributor_role:                Research Scientist at Bureau of Meteorology\n",
       "    publisher_url:                   http://www.bom.gov.au\n",
       "    publisher_name:                  Australian Bureau of Meteorology\n",
       "    publisher_email:                 oceanservices@bom.gov.au\n",
       "    acknowledgement:                 Any users of this data are required to c...\n",
       "    license:                         Data is not to be copied or reproduced w...\n",
       "    disclaimer:                      You accept all risks and responsibility ...\n",
       "    distribution_statement:          see license\n",
       "    notes:                           [1] This dataset is a preliminary foreca...\n",
       "    geospatial_lat_max:              -31.4953\n",
       "    geospatial_lat_min:              -33.9909\n",
       "    geospatial_lat_units:            degrees_north\n",
       "    geospatial_lat_resolution:       n/a\n",
       "    geospatial_lon_max:              115.7349\n",
       "    geospatial_lon_min:              114.7644\n",
       "    geospatial_lon_units:            degrees_east\n",
       "    geospatial_lon_resolution:       n/a\n",
       "    geospatial_bounds:               POLYGON((114.764397 -33.990898,114.76439...\n",
       "    time_coverage_start:             2020-12-14T00:00:00Z\n",
       "    time_coverage_end:               2020-12-21T00:00:00Z\n",
       "    time_coverage_resolution:        PT1H28S\n",
       "    time_coverage_duration:          P7D\n",
       "    processing_level:                0 (raw model output)\n",
       "    DODS.strlen:                     16\n",
       "    DODS.dimName:                    string16\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time</pre>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:       (direction: 24, frequency: 29, station: 167, time: 25)\n",
       "Coordinates:\n",
       "  * direction     (direction) float32 7.5 22.5 37.5 52.5 ... 322.5 337.5 352.5\n",
       "  * frequency     (frequency) float32 0.035 0.0385 0.04235 ... 0.4589 0.5047\n",
       "  * station       (station) int32 1 2 3 4 5 6 7 ... 161 162 163 164 165 166 167\n",
       "  * time          (time) datetime64[ns] 2020-12-14 ... 2020-12-15\n",
       "    init          (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead          (time) int64 0 1 2 3 4 5 6 7 8 ... 16 17 18 19 20 21 22 23 24\n",
       "Data variables:\n",
       "    frequency1    (frequency, time) float32 dask.array<chunksize=(29, 25), meta=np.ndarray>\n",
       "    frequency2    (frequency, time) float32 dask.array<chunksize=(29, 25), meta=np.ndarray>\n",
       "    station_name  (station, time) |S64 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    crs           (time) int32 -2147483647 -2147483647 ... -2147483647\n",
       "    cur           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    curdir        (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    mean_dir      (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    dp            (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    dp_sea        (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    dp_sw         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    dpt           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    efth          (station, frequency, direction, time) float32 dask.array<chunksize=(167, 29, 24, 25), meta=np.ndarray>\n",
       "    hs            (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    hs_sea        (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    hs_sw         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    latitude      (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    longitude     (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdir0         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdir1         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdir2         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdp0          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdp1          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    pdp2          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    phs0          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    phs1          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    phs2          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ppe0          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ppe1          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ppe2          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    psig1         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    psig2         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ptp0          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ptp1          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    ptp2          (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    swe           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    t01           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    t02           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    tp            (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    tp4           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    tp_sea        (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    tp_sw         (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    cd            (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    wnd           (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    wnddir        (station, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "Attributes:\n",
       "    product_name:                    per.2020_src.nc\n",
       "    history:                         Mon Dec 14 14:01:41 2020: ncks --overwri...\n",
       "    NCO:                             netCDF Operators version 4.7.7 (Homepage...\n",
       "    date_created:                    2020-12-14T14:02:20Z\n",
       "    program:                         WAVEWATCH_III(R)_version_number_6.07\n",
       "    Conventions:                     CF-1.6, ACDD-1.3\n",
       "    id:                              f1762_5247_9737_2398\n",
       "    naming_authority:                au.org.nci\n",
       "    title:                           AUSWAVE\n",
       "    date_modified:                   2020-12-14T14:02:20Z\n",
       "    date_issued:                     2020-12-14T14:02:20Z\n",
       "    references:                      http://www.bom.gov.au/research/\n",
       "    source:                          Perth triangular\n",
       "    comment:                         WAVEWATCH_III_switches: O0 O1 O2 IS0 IC0...\n",
       "    summary:                         (C) 2020. Research version of marine for...\n",
       "    standard_name_vocabulary:        CF Standard Name Table v47\n",
       "    keywords_vocabulary:             NASA Global Change Master Directory (GCM...\n",
       "    keywords:                        Earth Science > Oceans > Ocean Waves\n",
       "    institution:                     Australian Bureau of Meteorology\n",
       "    project:                         ep4\n",
       "    creator_name:                    Australian Bureau of Meteorology\n",
       "    creator_url:                     http://www.bom.gov.au\n",
       "    creator_email:                   oceanservices@bom.gov.au\n",
       "    contributor_name:                Stefan Zieger\n",
       "    contributor_role:                Research Scientist at Bureau of Meteorology\n",
       "    publisher_url:                   http://www.bom.gov.au\n",
       "    publisher_name:                  Australian Bureau of Meteorology\n",
       "    publisher_email:                 oceanservices@bom.gov.au\n",
       "    acknowledgement:                 Any users of this data are required to c...\n",
       "    license:                         Data is not to be copied or reproduced w...\n",
       "    disclaimer:                      You accept all risks and responsibility ...\n",
       "    distribution_statement:          see license\n",
       "    notes:                           [1] This dataset is a preliminary foreca...\n",
       "    geospatial_lat_max:              -31.4953\n",
       "    geospatial_lat_min:              -33.9909\n",
       "    geospatial_lat_units:            degrees_north\n",
       "    geospatial_lat_resolution:       n/a\n",
       "    geospatial_lon_max:              115.7349\n",
       "    geospatial_lon_min:              114.7644\n",
       "    geospatial_lon_units:            degrees_east\n",
       "    geospatial_lon_resolution:       n/a\n",
       "    geospatial_bounds:               POLYGON((114.764397 -33.990898,114.76439...\n",
       "    time_coverage_start:             2020-12-14T00:00:00Z\n",
       "    time_coverage_end:               2020-12-21T00:00:00Z\n",
       "    time_coverage_resolution:        PT1H28S\n",
       "    time_coverage_duration:          P7D\n",
       "    processing_level:                0 (raw model output)\n",
       "    DODS.strlen:                     16\n",
       "    DODS.dimName:                    string16\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds_spec = rompy.cat.bom.wavewatch3_nci(startdt=startdt,\n",
    "                                       enddt=enddt,\n",
    "                                       fmt_fields={'fcdate':[startdt.strftime('%Y%m%d')],\n",
    "                                                   'hr':['0000'],\n",
    "                                                   'grid':['PER'],\n",
    "                                                   'output':['spec']},\n",
    "                                       ds_filters={'rename':{'dir':'mean_dir'}},\n",
    "                                       hindcast=True\n",
    "                                      ).to_dask()\n",
    "ds_spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (dir: 24, freq: 29, site: 167, time: 25)\n",
       "Coordinates:\n",
       "  * dir      (dir) float32 7.5 22.5 37.5 52.5 67.5 ... 307.5 322.5 337.5 352.5\n",
       "  * freq     (freq) float32 0.035 0.0385 0.04235 ... 0.4171 0.4589 0.5047\n",
       "  * site     (site) int32 1 2 3 4 5 6 7 8 9 ... 160 161 162 163 164 165 166 167\n",
       "  * time     (time) datetime64[ns] 2020-12-14 2020-12-14T01:00:00 ... 2020-12-15\n",
       "    init     (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead     (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23 24\n",
       "Data variables:\n",
       "    dpt      (site, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    efth     (site, freq, dir, time) float32 dask.array&lt;chunksize=(167, 29, 24, 25), meta=np.ndarray&gt;\n",
       "    lat      (site) float32 dask.array&lt;chunksize=(167,), meta=np.ndarray&gt;\n",
       "    lon      (site) float32 dask.array&lt;chunksize=(167,), meta=np.ndarray&gt;\n",
       "    wspd     (site, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "    wdir     (site, time) float32 dask.array&lt;chunksize=(167, 25), meta=np.ndarray&gt;\n",
       "Attributes:\n",
       "    product_name:                    per.2020_src.nc\n",
       "    history:                         Mon Dec 14 14:01:41 2020: ncks --overwri...\n",
       "    NCO:                             netCDF Operators version 4.7.7 (Homepage...\n",
       "    date_created:                    2020-12-14T14:02:20Z\n",
       "    program:                         WAVEWATCH_III(R)_version_number_6.07\n",
       "    Conventions:                     CF-1.6, ACDD-1.3\n",
       "    id:                              f1762_5247_9737_2398\n",
       "    naming_authority:                au.org.nci\n",
       "    title:                           AUSWAVE\n",
       "    date_modified:                   2020-12-14T14:02:20Z\n",
       "    date_issued:                     2020-12-14T14:02:20Z\n",
       "    references:                      http://www.bom.gov.au/research/\n",
       "    source:                          Perth triangular\n",
       "    comment:                         WAVEWATCH_III_switches: O0 O1 O2 IS0 IC0...\n",
       "    summary:                         (C) 2020. Research version of marine for...\n",
       "    standard_name_vocabulary:        CF Standard Name Table v47\n",
       "    keywords_vocabulary:             NASA Global Change Master Directory (GCM...\n",
       "    keywords:                        Earth Science &gt; Oceans &gt; Ocean Waves\n",
       "    institution:                     Australian Bureau of Meteorology\n",
       "    project:                         ep4\n",
       "    creator_name:                    Australian Bureau of Meteorology\n",
       "    creator_url:                     http://www.bom.gov.au\n",
       "    creator_email:                   oceanservices@bom.gov.au\n",
       "    contributor_name:                Stefan Zieger\n",
       "    contributor_role:                Research Scientist at Bureau of Meteorology\n",
       "    publisher_url:                   http://www.bom.gov.au\n",
       "    publisher_name:                  Australian Bureau of Meteorology\n",
       "    publisher_email:                 oceanservices@bom.gov.au\n",
       "    acknowledgement:                 Any users of this data are required to c...\n",
       "    license:                         Data is not to be copied or reproduced w...\n",
       "    disclaimer:                      You accept all risks and responsibility ...\n",
       "    distribution_statement:          see license\n",
       "    notes:                           [1] This dataset is a preliminary foreca...\n",
       "    geospatial_lat_max:              -31.4953\n",
       "    geospatial_lat_min:              -33.9909\n",
       "    geospatial_lat_units:            degrees_north\n",
       "    geospatial_lat_resolution:       n/a\n",
       "    geospatial_lon_max:              115.7349\n",
       "    geospatial_lon_min:              114.7644\n",
       "    geospatial_lon_units:            degrees_east\n",
       "    geospatial_lon_resolution:       n/a\n",
       "    geospatial_bounds:               POLYGON((114.764397 -33.990898,114.76439...\n",
       "    time_coverage_start:             2020-12-14T00:00:00Z\n",
       "    time_coverage_end:               2020-12-21T00:00:00Z\n",
       "    time_coverage_resolution:        PT1H28S\n",
       "    time_coverage_duration:          P7D\n",
       "    processing_level:                0 (raw model output)\n",
       "    DODS.strlen:                     16\n",
       "    DODS.dimName:                    string16\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time</pre>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (dir: 24, freq: 29, site: 167, time: 25)\n",
       "Coordinates:\n",
       "  * dir      (dir) float32 7.5 22.5 37.5 52.5 67.5 ... 307.5 322.5 337.5 352.5\n",
       "  * freq     (freq) float32 0.035 0.0385 0.04235 ... 0.4171 0.4589 0.5047\n",
       "  * site     (site) int32 1 2 3 4 5 6 7 8 9 ... 160 161 162 163 164 165 166 167\n",
       "  * time     (time) datetime64[ns] 2020-12-14 2020-12-14T01:00:00 ... 2020-12-15\n",
       "    init     (time) datetime64[ns] 2020-12-14 2020-12-14 ... 2020-12-14\n",
       "    lead     (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23 24\n",
       "Data variables:\n",
       "    dpt      (site, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    efth     (site, freq, dir, time) float32 dask.array<chunksize=(167, 29, 24, 25), meta=np.ndarray>\n",
       "    lat      (site) float32 dask.array<chunksize=(167,), meta=np.ndarray>\n",
       "    lon      (site) float32 dask.array<chunksize=(167,), meta=np.ndarray>\n",
       "    wspd     (site, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "    wdir     (site, time) float32 dask.array<chunksize=(167, 25), meta=np.ndarray>\n",
       "Attributes:\n",
       "    product_name:                    per.2020_src.nc\n",
       "    history:                         Mon Dec 14 14:01:41 2020: ncks --overwri...\n",
       "    NCO:                             netCDF Operators version 4.7.7 (Homepage...\n",
       "    date_created:                    2020-12-14T14:02:20Z\n",
       "    program:                         WAVEWATCH_III(R)_version_number_6.07\n",
       "    Conventions:                     CF-1.6, ACDD-1.3\n",
       "    id:                              f1762_5247_9737_2398\n",
       "    naming_authority:                au.org.nci\n",
       "    title:                           AUSWAVE\n",
       "    date_modified:                   2020-12-14T14:02:20Z\n",
       "    date_issued:                     2020-12-14T14:02:20Z\n",
       "    references:                      http://www.bom.gov.au/research/\n",
       "    source:                          Perth triangular\n",
       "    comment:                         WAVEWATCH_III_switches: O0 O1 O2 IS0 IC0...\n",
       "    summary:                         (C) 2020. Research version of marine for...\n",
       "    standard_name_vocabulary:        CF Standard Name Table v47\n",
       "    keywords_vocabulary:             NASA Global Change Master Directory (GCM...\n",
       "    keywords:                        Earth Science > Oceans > Ocean Waves\n",
       "    institution:                     Australian Bureau of Meteorology\n",
       "    project:                         ep4\n",
       "    creator_name:                    Australian Bureau of Meteorology\n",
       "    creator_url:                     http://www.bom.gov.au\n",
       "    creator_email:                   oceanservices@bom.gov.au\n",
       "    contributor_name:                Stefan Zieger\n",
       "    contributor_role:                Research Scientist at Bureau of Meteorology\n",
       "    publisher_url:                   http://www.bom.gov.au\n",
       "    publisher_name:                  Australian Bureau of Meteorology\n",
       "    publisher_email:                 oceanservices@bom.gov.au\n",
       "    acknowledgement:                 Any users of this data are required to c...\n",
       "    license:                         Data is not to be copied or reproduced w...\n",
       "    disclaimer:                      You accept all risks and responsibility ...\n",
       "    distribution_statement:          see license\n",
       "    notes:                           [1] This dataset is a preliminary foreca...\n",
       "    geospatial_lat_max:              -31.4953\n",
       "    geospatial_lat_min:              -33.9909\n",
       "    geospatial_lat_units:            degrees_north\n",
       "    geospatial_lat_resolution:       n/a\n",
       "    geospatial_lon_max:              115.7349\n",
       "    geospatial_lon_min:              114.7644\n",
       "    geospatial_lon_units:            degrees_east\n",
       "    geospatial_lon_resolution:       n/a\n",
       "    geospatial_bounds:               POLYGON((114.764397 -33.990898,114.76439...\n",
       "    time_coverage_start:             2020-12-14T00:00:00Z\n",
       "    time_coverage_end:               2020-12-21T00:00:00Z\n",
       "    time_coverage_resolution:        PT1H28S\n",
       "    time_coverage_duration:          P7D\n",
       "    processing_level:                0 (raw model output)\n",
       "    DODS.strlen:                     16\n",
       "    DODS.dimName:                    string16\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import wavespectra as ws\n",
    "ds_spec = ws.read_dataset(ds_spec) #Convert it to wavespectra format\n",
    "ds_spec = ds_spec.sortby('dir')\n",
    "ds_spec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the spectra within some distance of the boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAGsCAYAAAAbq4Z0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABxRklEQVR4nO3deVxUVf/A8c8BERAQXBFcUlPcWFRQUHJJS832vadsfVpsNftlqfWUrVr2ZGpppW1PtrpmtmhllBsgKAIuoOKSgqAoyC4w5/fHzCAqA4PMzJ3lvF8vXi+4c+fe73H5cubcc75HSClRFEVR7Jub1gEoiqIoDVPJWlEUxQGoZK0oiuIAVLJWFEVxACpZK4qiOACVrBVFURyAStbnEUJ8KoTIE0Kk1zp2qxBipxBCJ4SIqnW8qxCiTAiRYvj60MQ1ZwghjtY6b3yt12YLIZKEECNMXDNFCHGPNdtcKxaLt91w7pNCiAzDdd6uddyp2y6E+K7WOQeFECm1XrOLtlup3f2FEPGGc5KEEINrvWYX7XZIUkr1VesLGA4MBNJrHesD9ALigKhax7vWPq+ea84Anq3jeG9gNtAC+L4x13Sgtl8O/A54Gn5u7yptP+/6/wVesre2W+nvfB1wleH78UCcvbXbEb+a1ZG/XZqU8m8hRNfzju0GEEJY+nbugA6QgMUv3lhWavujwCwpZYXhenmG467QdgzvF8BtwCjDIbtpu5XaLYGWhu/9gWzD93bTbkekhkGarpsQYrsQ4i8hxLB6zntCCJFq+NjZCkBKuRN9L2MjsLDWuZee97GwvutqyZy2hwDDhBAJhvMGgcu03WgYkCul3AsO33Zz2v00MFsI8Q/wDjANHL7dmlM966bJAbpIKfOFEJHAKiFEPynl6fPOWwi8hr5H8Rr6j8QPAEgpn6zjuvullP2tF7ZFmNv2ZkArIAYYBHwvhOgu9Zy97Ub/Ar6pfcBB225uux8FJksplwshbgM+Aa4Ah223XVA96yaQUlZIKfMN3ycD+9H3JM8/L1dKWS2l1AGLgMHnn+NozG07cARYYUjOieg/Bre1XaSW14i2I4RoBtwEfGe7CK2jEe2+F1hh+H4pTvDv3R6oZN0EQoh2Qgh3w/fdgZ5AVh3nBdX68UYg/fxzHI25bQdWYRirFUKEAM2BEzYK0yoa0XbQ9yj3SCmP2Co+a2lEu7OBEYbvRwF7bROhk9P6Cae9faH/uJoDVKLvFf4bfYI9AlQAucBaw7k3AzuBHcA24Npa11mM4Uk68CWQBqQCq4Ggeu7fFSgDUmp9PeXAbW8OLEH/C2obMMpV2m74+XNgohn316TtVvo7vwxINpyXAETaW7sd8UsY/sAURVEUO6aGQRRFURyAStaKoigOQCVrRVEUB6CStaIoigNQyVpRFMUBqGStKIriABx2ufm4cePkiRP2u7aiqKgIAD8/P40jsT1nbntlZSW5ubl4enrSrl27C17Xqu06nY69e/fi4+NNUFCgTe8NUFJSBoCPj7fN730xUlLS10opx2kdR2M4bLI+ceIESUlJWodhUlxcHAAjR47UNA4tqLbbvu033ngDl17aiQ8/nI27u7tN7w0QH58KQExMuM3vfTFaterpcCUPHDZZ27uQkDpLRbgE1XbbKS8v5/777+PXX9eyZ89mTRI1QLduHTW5rytRydpKgoODtQ5BM6rttlFeXs5NN92Ihwekpsbh76/dsFNgYBvN7u0qVLJWFAdkTNSenm58/PE7eHh4aB2SYmVqNoiVZGdnk52d3fCJTki13bptl1Jy88032VWizs3NJzc3X+swnJrqWVtJZmYm4JpDAqrt1m37b7/9xsaNG9m3L8EuEjXAgQNHATUcYk2qZ60oDuTPP//krrvu5MsvP7CbRK3YhkrWimJCRUUF06dP56uvvtI6FE6cOMGUKVO47bZb+eyzuQwfPkTrkBQbU8laUUx48cUXWb16JdOmPc8HH3ygWRw//vgjnTp14tChTFas+JTLLovWLBZFO2rMWlFMKC8vp3Pnjsya9QLXXXcPe/bs4aWXXqpz5aI1lJSU8Oqrr/LZZ5/y009fERkZYZP7KvZJ9awVxYTx48ezfXsabdu24ccfl3D8+BE6duzII488wnfffYe1d1m655672bEjiVWrvlCJWlHJWlFMGTt2LL1792bNmt/o0qUjH344mx07/qRly2a88spLXHpp95rZH5b0119/MXbsGNatW8eCBbPo29d1V4QqZznsHoxRUVHSnmuDKM7h7rsnEB7eg4cemnDO8erqar78cilvv/0BqalptG1rmVITX3/9NU8/PYnHHruPhx++hxYtHKMwkqNp1apnspQySus4GkMla0Wpx59//sktt9xMmzatARBCMHXqk9x443gAbr75AY4ezeXFF/9D69atGTduHG5uZz+wpqSksG/fPsaMGUPLli1N3mfjxo28+uorpKSksGrVF6o3bWUqWduQStaKreTl5XHq1CkAcnJyGD9+PPv3J+Lt7UVZWTlLl64mLm4ze/bsw83NnRdeeJGvvlrCoUOHOHYslz59erJ3bxaBgYHceutthIeHI4RACMGuXbv48sv/cezYMR56aAIPP3yPpjU+XIVK1jZk78naGFtUlEP9e7AIZ297YGB7unTpSOfOHXnmmYmEhvYGoLy8gq++WsWOHbvp2LEN48ePpn37drRv35a9e7PIzz/Fxx9/SVlZOVLql417enrwyCP30LdvCK1bt9K4ZRcvLW0vAGFhPTWOxDyOmKzV1D0rKS4u1joEzThz26dPn0bXrl25664J7NiRwrXXTuD116dy11234OXlSb9+fejXr88FdZ1DQi4FYMgQh8oPZjNuPqBYj0rWitIIs2a9hZSSLl26sGzZcm677VZeffVd+vcPo1+/XlqHpzgxlawVpRGeffZZ0tPTkVLSvn07Ro8eTX5+PgkJySpZK1alkrWimKGiogI3Nzdmz57NyJFD2bJlC3l5eXh4eBAY2IH//vctHnjgTq3DVJyYWhSjKA3YuXMnHTsG07t3L+bOnUt6egb/+c9/aqreRUREsGfPXnbtsvwCGUUxUj1rRWnA33//TfPmHpw8eZKePXty/PiJc17v2rUrUVFRbN26nV69+mgUpeLsVLK2kqCgIK1D0IyztX3ixIn4+fnh7+/PVVddVec5jzwykddfn8FHH82xcXT2oX371lqH4PRUsraSXr1c92GTs7VdCMGECRPqPWfDhr8JCupA9+6dbBSVfXHVdtuSStaK0kTbt29n8eLFrFmj/SYFivNSDxitpKioiKKiIq3D0IQztv306dNcc83V/PHHHxe89vjjj9WsZCwuLqW4uFSDCLXlqu22JZWsrSQ5OZnk5GStw9CEM7b9008/5aeffmbnzp0XvHb33fewaNESjh3LIz19H+np+zSIUFuu2m5bUslaUczwxBNPsGvXLp544okLXnv00UcZMGAg8+YtsvqGBIrrajBZCyG8hBCJQogdQoidQohXDMdnCyH2CCFShRArhRAB5r7X8FqwEGK9EOIHIYSv4VgvIUScECJFCLFbCPGxBduqKBetWbNm9OnT55zyp7UtXvwJv/76F4cPH7VxZIql2HuuM6dnXQGMklJGAP2BcUKIGOA3IFRKGQ5kAtMa8V6Ap4AngcWA8VH7PGCOlLK/lLIPMN+M+BRFc+3btyckpCeVlZVah6JcPLvOdQ0ma6lnLKPmYfiSUsp1Usoqw/F44IK5O6bea/jZHdAZvoThWBBwpNb70xqKT1HswdKlS9mxYwft21tmxxjF9uw915k1dU8I4Q4kAz2AD6SUCeed8gDwXSPf+z7wJVAIGIsqzAHWCyE2A+uAz6SUBXVdt6ioiLi4uDrjDQkJITg4GIDs7Ox698kbOXJkzfdJSUkmy3sGBQXVzB8uKiqq9wFaZGRkzfcZGRnk5OTUeZ6vr+85NZ9NtQfso01+fvqi+A21qTZnaVNkZCTvvPMOHh4e9O/f/5zXz5w5wzfffM306c/i7n72v1Rubj4HDpgeFqldRjUtba/JMqPt27eumcdcXFxa74O80NAe+Pq2ACAr6wh5eSfrPM/Hx/uc2tPx8akmr9mtW0cCA9sAptt05EjuBcfsuU2m2GOuMzLrAaOUslpK2R/9b5TBQojQWgG+AFQBdU4yNfVeKeUhKeVwKeW1Usoiw7HPgD7AUmAkEC+E8DQnRkWxpk8//ZR3332HyZMnX/AQ8eeff6JjxyB69OimUXSKpdhzrmv0TjFCiJeBEinlO0KIe4GJwGgpZYOTLGu/18x7pQP3Sikv6B7Z+04xxnnGxp6bK3HGts+bN48vvviU8vJybr/9Tl566aWa16KjB/Paa1MYPHhAzVxjX98W5OTk8vnn31JdXY0Qgttvv8FpE3rtdjsCc3aKsZdcZ2TObJB2xqefQghv4ApgjxBiHPA8cJ2p4E29t557jRNCeBi+7wC0ARzy8bqfn59TJavGcMa2Hz58mNjYQZw6VcDs2bPP6V337NmTd9/9kHnzFjF06HhCQmIIDg6nf/9R5OeXERAQTGWlB7Gx1xASMoSUlHQNW2Idvr4tHCZRm2Lvua7BnrUQIhz4Av0guRvwvZTyVSHEPsATyDecGi+lnCiECAYWSynHm3pvPfd6F7gaKDccmi2lXFLXufbes1acR1VVFb1796JZM0F09FA++eRTmjU7OzZdXl7O7NmzOXnyJCNHjuTKK68EwM3NDS8vr5rzKisr+emnn5gwYQJubgJfX18eemgCLVp40axZM2677Xr8/HwvuL9ieXX1rO0119W8x1En8dt7ss7IyACcr6iROZyx7RMnTuSjjz7iyy+/rLeokzltLy8v58yZM2zbto3ly5cDsG/fPhITExgwIIxPPnnP4XY4z8rST2xwlIJOjrhhrlrBaCU5OTkmZxc4O2ds+8KFCzl+/HiD1ffMabuXlxctW7Zk5MiRzJ8/n/nz5/Pzzz+zY0cqvXqFcvPND1BYeNqS4VtdXt5Jk7M0FMtQyVpRzCCEoG1b682hFkLQqVMn3n//fYYMuYybbnK8hK1Yl0rWimJHhBDMmzePfv3CefHFWVqHo9gRlawVxc4IIXjllVf4448N/PjjWq3DUeyEStaKYoe6devGnDnv8eqr71JUVPdqTcW1qGStKHbq1ltvZdiwEUyd+rrWoSh2QCVrK/H19b2gToarUG23TNvd3NyYMWMGv//+N7/+ut4i17QWHx9vfHy8tQ7Dqal51opi55YtW8b777/HqlWfax2K01DzrBVFsbhBgwaxc+ce9u8/qHUoioZUslYUO3fJJZfQr18/p6wpophPJWsriYuLq7eWszNTbY+z+HUffvgRZs6cR1VVVcMnayA+PvWi6kcr5lPJWlEcwG233UaLFr4sW/aj1qEoGlHJWlEcQLNmzXjppZdYsOBzrUNRNKKStWJTV91ajHDTIYREuOm46la14MNcffr0oaSk9IKdahTXoJK1YjNX3VrMr8t8aCarAAHSjV+X+dQk7POTkEpK57r00kvx9PRm7do/tQ5F0YBK1orN/Lq8BSA4RFfK8GI5NzKBJWQsO8aMP19m8tqz+xtKKZm8djIz4mZoGrM98fLyIiIiXC0/d1Fm7W6uKBYhBaAjkFzc0XETq7iJVQAU3tyCdUGlrFq5nRv+PZv/O76EuUnzmRQ9CSklQghNQ7cXbm5u6HTqE4crUsnaSkJCQrQOQTMm2y4kvWQm7ug4QRveZgqxbCaWTbQ9lc+tp4Bdf6P7OJpPn4dJIycxZ+wch0rU1v57d3Nzs8vhoW7dOmodgtNTydpKgoODtQ5BM6baPu7mUoKXbQTgTy5nNs8zG8m4m4v55fWjyI0b+eyDh2hdBoXeOFyiBuv/vQshkFJn1XtcjMDANlqH4PTUmLViMz9/78PtXf8CYBNDQVQz7pYSflnmhy4khMmd0vn3DXDjv/Tn1x7DVvSEEOh09pesFetTPWsryc7OBlyzh11X22fEzaCgvIA5zRMB6PxeAk+FPEMr71bodC8RuSiSlGMpTIrWD31MXjuZuQlzAcfqYVvz7/348eP88ssv3HvvDRa/dlPl5uo3/lY9bOtRydpKMjMzAddM1ue3XUpJQXkBX62fy3uZUOHZjKn531GVCE8NforJ6yaTciyF/h368+6YdxFCMGfsHAACvAIcJlGDdf/ec3Nzadu2DQMHhlv82k114MBRQCVra1LJWrE6Y/IN27wfWMPmoCqaV0PyR7BpzTweuxomxUzi3THv4ubmds577ClRHz9+nJ07dzJy5EitQ1FckErWik0IIXjgtpnMWr+GvW1g8FEIz4NKd0DAv1r9i6ioKLp27UpsbCyxsbEMHDiQ5s2bax06ADqdjlGjRnHgQBZhYWFs2LCRZs1s+9/n1KlTdvXLS7EtlawVm5BSMvnIYuZeCaTewYvxGcB2NnXWv37/f+5n9/bdbN++nZUrVwL6RSCDBw+uSd5DhgyhdevWNo89ISEBb29vcnKymTDhFj799BseeOAB/vnnMN988y0dOnSwSRzvvTeH22671ib3UuyPmg2iWJ1xNeLchLl0Tl0Epe0ZVrkdgE2J30D8E+xO3g3Ay1OnEn5FH9p1bkN5eTl///03M2fO5JprrqFNmzb069ePhx9+mC+++IJ9+/ZZfbZIYWEhMTExjBp1Oe++O4fCwgpmzZrF0qVLOXz4IG+++aZV71/byZMn6d27p83up9gX1bNWrE4IQYBXAJOiJzH3lQdoNuwlrsjSvzaOX3H7NYZNwDFPD4oHnCRL7OaxRx9lysDn2JKYyKb4eDbFx7N12zZ27drFrl27WLRoEQCBgYEMHTqUyy67jNjYWAYMGGDxoRN3d3dKS0spKyvjm2++RUpJRUUFe/dm8vrrttnMNj4+nvT0dKKiImxyP8X+qGSt2MSMkTOQUjJXwvC/L0fwBgD38wX38wUAS90reSDlYx6LepRZV8xCCMG1V13FtVddBUBFRQXbUlLYlJBQk8Bzc3NZuXKlyaGToUOH0qpVq4uO29/fn507d/LCC9M5dOggVVVVvPrqK5w6VcDixZ/UPBC1tv/850X+/e87adPG9sNAin1QG+YqNiXcdCDdaM0JruZnOnGEWBYTywEWdIKZd8PpqYVmPUiTUrJv/342JSSwccsWNsXHs8cwda62vn37EhsbW9P77t69e6Me1E2fPp0//1zH7t17mTDhbj7//DNKSkrJycmxyXj15s2bufHGG4iP/5lWrQKsfj9X4Igb5qpkrdiUsUwq1E6WQxDEE3ATVPbinJ51Y+Xn57P5vKGTioqKc84JDAys6XmbM3Ty3Xffcc899zBixHC+++57nnjicY4fP85PP/2Mh4dHo2NsjIKCAoYNG8Ytt4xn0qSHrHovV6KStQ2pZO24rrq1WF8uVQrwyoYzXUDqOHXoMDO3zmRB0sImJezaag+dGHvfJ/LzzznHnKGTpKQkQkND8fLyalI8jTVz5kx+/XUNq1Z9rqbtWZBK1jZk78naGFtUlEP9e7CIxrT9gfkP8NlTnxHWrx+pW7YgpWTq71Px9/Jn+rDpFo9NSsneffvOGfeua+ikX79+5/S+zR06seTf+5dffsnzzz/HqlWfExJyaZOvZ01paXsBCAtzjNkqjpis1QNGKykudt0C8Y1pe0iJvqToZTExgH7miCV61KYIIQjp2ZOQnj25f8IEAE7k518w62Tnzp3s3LmTjz/+GDB/6MQSf+/5+fm88cYbfPPN1w6RqAFKSsq0DsHpqWStaGrTpk0AxBqSNWDzj/tt27S5YNZJckpKTfI2zjpZsWIFK1asAMDb2/uCBTtNmXVy7Ngxvv32W3Jzc3n77be54YbxrF79P3r27G6RNiqOTyVrRTM6nY7NmzcD5yZrrXl6ejI0Opqh0dFMmTTJ5NDJX3/9xV9//VXzvn79+tGtWzdCQ0Pp3Llzg0MnBw8e5LXXXmPp0qVUVlZy3XXjaNXKn40b19Cnj2MMJyi2o5K1opmMjAxOnjxJcFAQl3TponU4JpkaOtlcK3nXHjpZs2YNs2bNokOHDhcMnRw4cIDExETi4+P55JNPuPXW69iy5SdatvTDz89X45Yq9kwla0UztYdAHG2mQ9s2bbhu/HiuGz8eODt08r9vvyV91y4yMjM5duwYy5cvZ/ny5YA+6bu5udGjRzc6d+7Ib799T2hoHy2boTgQlawVzdQk6+hojSNpOuPQyZnqagBGDBlC5r59bIqPZ/Kz/4do1ozComKqq6vJyNhHRsY+fv/9L3r37klMTCTR0ZHExAzkkks6O9wvLsU2VLK2kqCgIK1D0Iy5ba/r4aKjCwoMBPS96F49e9KrZ0/m/fcdJk36N8NHDCUhaYfhK4XtqTvZs2cve/bs5fPPvwUgMLAd0dEDiY6OJDp6IOHhfa2+8MYS2rdXy+CtTc2zVjSRl5dHYGAgLVq0oOCffxwiIV2stLQ0brjhel555TmuuebKmuPl5RWkpO0mISnF8LWDk6cKznmvt7cXAweGG3rekQwePAB//5Y2boHzUfOsFcVMxlkg0VFRTp2oAcLCwli16gduuOF6gJqE7eXlScyg/sQM6g8Yap1kHTonee/df5BNmxLZtEm/d6UQgt69e9b0vtXQietQydpKioqKAPDz89M4EttrqO1SynPGq6WUTpNsigyLYvx8z53ZUTthCwFXX33lBe8VQtDz0q70vLQrE26/AYD8k6dISNpB/FZ9Ak9J28Xu3Zns3p1pV0MnxcWlAPj6trDpfV2JGgaxkri4OACX3K+vvrYbdzlPeCWB+Ph4flm+nD/53WrLy20tzvCJYeTQoXW+npqayo033sBrrz3P+PFXNPr65g6dREZG1CRwWwydxMenAhATY3+b+dZFDYMoSj2Mu5zP3TAXsdUNENwyfw1iyGeEi/uZdpnz9LBNCQ8PZ+XKVdx44w1IKRk//opGtbm+oRNj73tf1iE2bkxg48YEQA2dOAuVrBWbMe5Ynjs/iberN7HZE7Z6fMYvf99FasI7THM/xawpzj+rIDw8nBUrVnLttdfyxRdL+eyz9/Dxubjhg7qGTk7knyQxeQfxW3eYHDrp0KE999xzG9OmTbJUsxQrU3swKjYlhKDzOh86A7dXwDu/wZBNI6DKg8XzvLUOz2YiIiI4ePAgPXv2ZsKEx9i0KcFi127bpjXjx1zOqy88zdqVn3MofQO/LP+MGdMmcdWVI2jdKoBjx/L48MMvLHZPxfpUslZsSkrJOwEZ9AMeD4evQ2Fz7N+ApKLQtrWitebm5sb8+fO55ZY7mDTpP8yd+zHFxSUWv49x6GTSo/fx9SfvMX/2SwCEhanVk45EJWvFZqSUPP3L08iSQ+wCPk9N5ZEOj/LP8CVwxVSa+7temU03Nzcef/xxfvxxDdu27eGee56gqMi65XXjt+4AIDo60qr3USxLJWvFZoQQVB2vgjIQdEBwCfw+C7Y+Cmf8eOipcq1D1EyfPn1YuXIl4eEDuffeJ9HpdFa7V0JSCgDR0QOtdg/F8tQDRiuJjHTdXouptkspGVA5AICQPpEczy7nTKEXzbe+zENPlTvFw8XI8Iufuubm5sacOXPo3r0ba9b8xnXXjbVgZHr6qX+7ABg8eIDFrhsa2sNi11LqppK1lbjiYhijutpunF9duKkQgI6dIznW7UU43RqR9oytQ7Sa8xfDNJYQglWrfuCmm27E3d2tzsUzTZGStpszZyrp3bsnAQH+FruuWgxjfWoYRLG6mvnVCXMRP3xHb2BLQRrVEYvA6zQVhV4smNmSqbNPah2qXYiIiGDVqh946aW3Wb16rUWvrYZAHJfqWVtJRkYGAL169dI4Ets7v+3G+dVuB0/w7qmvAMhPXUli1iXEnejCAgSyqhmL53kza4pmYVtExv79APS6tGn7JoaFhfHDD6u54YbrkVLH9ddfZYnwapJ1TIxlh+myso4A0L17J4teVzmrwZ61EMJLCJEohNghhNgphHjFcHy2EGKPECJVCLFSCBFQzzXchRDbhRBrah0LFkKsF0L8IITwNRzrJYSIE0KkCCF2CyE+tkAbNZGTk0NOTo7WYWiirrYLIRijG8r3wFF3aFMOV504xHX8WHOOM0zdy8nNJSc31yLXCg0NZfXqH3n99fcs0sOWUpKQpJ8JYulknZd3krw8x/5kZO+5zpxhkApglJQyAugPjBNCxAC/AaFSynAgE5hWzzUmAbvPO/YU8CSwGJhgODYPmCOl7C+l7APMNyM+xQFIKXnhj0+4HegWDf0ehX/3G818nqg5x9PfdWeDmJJVGUD7++Yza283blyUxbo9py/6Wnv3H+TkqQICA9txySWdLRil07DrXNdgspZ6xomfHoYvKaVcJ6WsMhyPB+r8/COE6ARcbQi0NndAZ/gyFikIAo7UundaQ/Ep9k9KyeS1k9mWuA2AGy75nMP7H+H76/7gpys2ABLRrIoHn3K9edb1WbMjm5dX76SgQiCEIK+omrfW5V10wq49Xq3qglzI3nOdWWPWQgh3IBnoAXwgpTx/bewDwHcm3v4e8Bxw/hSB94EvgULgTsOxOcB6IcRmYB3wmZSywJwYFfslhMBH+OB2zA0pJIveHs2bH43go206qstb4ulfzoNPlTnF1D1LWLMjm7l/7CW74MJfXuVVko825jOmd+Or6BmHQNRiGNPsOdeZlayllNVAf8NYzUohRKiUMh1ACPECUAV8df77hBDXAHlSymQhxMjzrnkIGH7esc+EEGuBccD1wCNCiAgpZcX51y4qKqopxXm+kJAQgoODAcjOziYzM9Nk22qX8UxKSqK4uO7VY0FBQTUPzIqKikhOTjZ5zdrzjDMyMkyOXfv6+hIVdbZKo6n2gH20yTglr6E21WZsU2BGILoqHd27dWP7zp1cdRk8dc8UOnbooG9T7uma8qJ1tqlWydGk1FTTbQoMrHm4V1RcTHJqquk2hYfXTLXL2L/f5Fizr68vUbXmT9cXZ1Fxcc01s3NzyTQ8cKxLXW3KzC0iLuM4rXQ6WhmG8E9WeZNTpf+z9xKVtDmTR/y2C/5LABDa+xJ8W+hrrGQdPkbeiYKzcW/Ub2Dg49OStLS9hIX1rHnNWOK0Lt26dSQwsA0Aubn5HDhw9IJzjhy58M8uLW0vJSV1f1pq3751zcPI4uJS0tP3mbx/aGiPmqmBWVlHTI6N+/h4m90mU+wx1xk1auqeIfPHGW6AEOJe4BrgLll3YexY4DohxEHgW2CUEGJJA/fIllJ+KqW8Hv0fTGhjYlTsU3p6OgCh/frVHFMfxS+UcOAkVQ2sXnSrrqDasDGvuQoKT3M0OwdPz+Zcemm3poToEuwx1zW4+YAQoh1QKaUsEEJ4o++yv2W4+LvACCnl8Xovor/OSOBZKeU19ZwzDvhDSlkphOgAbAcGSCmPnX+uvW8+YIytds/ZVZzfdikl1113HWvWrGHJokXcedttTpuokww9+aiLXMkY+tKv1Pc/0svDDc+dP3DL4EAeeOAus6/787o/uevBZ4iNHcyaNRd0DJssLW0vwDk9W3tW1+YD9prrjMwZBgkCvjCM5bgB30sp1wgh9gGewG+G/3jxUsqJQohgYLGUcrwZ1z7fGGCuEMI4LWBKfcHbM1dM0ka12z4jbganSk/V7Lk4NDqaqb9PdZqdYc53sUnaKCjAu86xaoDgAG8mje5JyxFe3HXXnUyYcCvNmzc367rWLt7kKEm6AXad6xpM1lLKVOCCIgJSyjqLAUgps4ELgpdSxqH/WFHfvZ4BnGftsYszrlz8cPU8qk5CcFAQd33zLml8Dlsn8t87S52mJoilTBrdk5dX76S88uwwh5eHO69c149rIvTPLKQMIiwsnKVLf+Suu24267pq5WLD7D3XqRWMitUYVy52mb2WB9lD/OkcEjd+TuuSa0na9h9KaMGCmc2BkyphGxgT8tw/9pJTUEaQoTdtPA76P9fQ0FDKy82b6mit4k2KbalkbSVqw1x924UQ9D3mRktgTAmM2QDwI0tpzgN85jTLzI0a2jDXHNdEBJ+TnM8npSQ1dQfXX29ekSdj8aY+fUIsWrypNkfbMNcRqWStWJWUkjsKjtICuHwkDC2GmMx2bDwdW3OOMywzt6VTp06xY8cOliyZZ9b58Vu3A2oIxNGpZK1YjZSSh799mMLsQkq9mvHLtoOsHvwGPL4Qtu6F3yUgDMvMXWf/RUvw8PCgWTPz/vueHa9Wi2EcmUrWitUIITi9T780etjgoQwaUcUHM1/Xv1juD4hay8xVsm6Mhqbc1j4vMdk4RKF61o5M1bNWrKpzob5gUGxMDLOmtObxaUU03/oybJyGp38Zj007rR4uNpKfnx9t27ZjwYLPGzxXFW9yHqpnrVjVpk2bAH2yBpg1pXWth4neqB5143l4ePDDDz8QFRXFY4/dV++5qniT81A9a8VqysrKSE5ORghBzKBBWofjVIKDg2nfvh1r1/5Z73mqeJPzUD1rKwkJCdE6BM0Y275161YqKysJDw3F3986U8bsTUgTd4gxV7Nmzbj99juIj09m7NjLTZ5nq8Uw3bp1tOr1FdWztprg4OCaKnmuxtj2miGQ6GiNI7Kd4MBAggMDbXKvO+64gyVLlpl8/UT+SfZlHcLb24vw8L5WjSUwsE1NZT7FOlSyVqzGmKwvGzJE40icU48ePWjZ0o/ff/+7ztcTk/VDIJGREXh4eNgyNMUKVLK2kuzsbLKzs7UOQxPZ2dkcOXKkpniT8eGiK8jOzSXbQnswNkQIweOPP0F8fN11yLdsTQFssxgmNzef3Nx8q9/Hlakxaysxbg7gikMhmZmZHDx4kFOnTtExOJgunV1nyphxswFbDYUEBASwZk0qOp0ON7dz+162XAxj3JBADYVYj+pZK1Zh3GwgNiZGTRmzonvvvZfMzP0UFBSec7y8vIIdabsRQqjiTU5CJWvFKmqStQs9XNSCu7s7bm5uF/xC3J66izNnKunduyf+/o3fr1GxPypZKxZlXAZtTNZDVbK2ifOTtapf7XxUslYsZkbcDCavnUx+fj5Hjx7Fx8eH7/K+4c0Nb2odmlOrq06IKt7kfFSyVizCuCvM3IS5fLfufwC06R7Ah9s/4u0P8vFrqaNt5zKmzq57Z2rl4kkpz+lZq+JNzknNBlEswrgrDEDfBz7hPWBL9lESVwwhPuNfpFJNRaE3C2Z6oHaGsS5VvMk5qWRtJa64Q4wxYa87OZdg4OZ8uDl/CzCKnfQlhnin2xnmfE3ZIeZinDp1ioqKM+dM29OieJPaIcb6VLJWLEZKyZOrn+SDM9ALGHUlDEnvTUxOFSn0rzlP7QxjOUuWLOGKK0bg5+dbc8xYvCkmRo1XOxOVrBWLkFIyee1kPlj1AUjwCgvja48YvrxvEWx9FPffX685V+0MYzlnzpyhS5dzF16d3cZLJWtnoh4wWklSUhJJSUlah2EzQggCvAIYin4YoGePHjwycDokPwLl/lSjr01xdmcY55SUmkpSaqrN7ldYWHjOEMjxEyfZf+AwLVp4ExbWx2ZxpKXtJS1tr83u54pUsraS4uJiiouLtQ7DZqSUzBg5g1bHWwHQKySEt55rw2ODp+KZNhmQLrEzjC3/3qWULFy4kOuuG1tzzFi8aeDAcJsWbyopKaOkxHl/CdsDNQyiNNmMuBkUlBfw3yv/W1O8qdD3JG9ueJO3npvOW88Zz1Q7w1iSlJLyigr69DlbOz1eLYZxWqpnrTRJ7fnV9y2+j1OnTuEX4EtGeaaaX21lUkrOn+uhFsM4L9WzVmqs2n6U2WszyC4oIzjAmylje3HDgPp3AKk9v/rktLnEASnu1ZxIbM3u316gADc1v9pKqqurz5map4o3OTfVs1YAfaKetiKNowVlSOBoQRnTVqSxavvRBt9rTNijD8EIYFJ+Ga/9tpJ/6Eoig2nLiZr51YrlfPrpp8TERNX8rIo3OTeVrBUAZq/NoKyy+pxjZZXVzF6b0eB7jdP2ntHBtcD/wgNJb9+Rcjxpx3FOoK9xrOZXW9axY8eIjT27EbEq3uTc1DCIlQQFBWkdQqNkF9T9JN/UcSNjop77x1wohD99fDjS9xZ01/7GkS1X03X9Q2AYWXWF+dVBNtp0AExX2qvd27aV9u3V8Ja1qWRtJb169dI6hEYJDvDmaB2JOTig/uRqnF99tdfV/MRPxAwaRGSPZ1mQWAVnWpGJ/s/h7Pxq507WvWy0uzno/+yNFfd0Ol2tlYu271l3797J5vd0NWoYRAFgytheeHu4n3PM28OdKWMb/qUzY+QMQkr008dio6Ndcn61Fmon6737D3KqoJAOHdrTpYtKnM5IJWsrKSoqoqioSOswzHbDgI7MvCmMjgHeCKBjgDczbwprcDaIUe3NcYuKi3nxMU9O/ONN0WnBiX+8XSZRFxUXU2SDRTE6nY5ffvmZTp30S821KN5UW3FxKcXFpTa/rytRwyBWkpys33Hakarv3TCgo9nJ2UhKSXl5Odu2bcPNzY3oqCiSDcutbV2Bzh7Yqu1VVVXs3r2Hm2++BjhbvEmrh4vp6fsAVX3PmlSyVi6aceXijd43UllZSURYGG8mvolvfhuGdblM6/Ccmpub2zmbDqjFMM5PJWvlotReubj1QCIAbh11LEhaSOyJZ9j4teSW7WU8+FSZywyBXIw1O7KZ+8decgrKCArwZtLonlwTEdzg+9avX0+rVgGAdsWbFNtSyVq5KLVXLg6fPZd3gfjMnSTmXEl1iyhOeQq1crEBa3Zk8/LqnZQb5rdnF5Tx8uqdAA0m7N9//51HHrkb0K54k2JbKlkrF00IwX+veIeTJXNpB0RnAfwG/Ea2VxARpFBe5e3UO8M0xdw/9tYkaqPyymrm/rHXrN51s2b6/76qeJNrULNBlIsmpeS+T+7nUmCMN7wZC+tbd6LMzYtq4U65YU61WrlYtxwTC45MHa9NiLO7mqvxategetbKRTGuXFyyZgkArcbdyNsVbal+ZBGX5TxDmz0DIF5/riusXLwYQQHeda4QDWpgIRKAEPoHjGXl5aSk7lLFm1yAStZWEhnp3L0c48rFPmV92M1uLhsyhK6nb2VBohvbq5rDztH681xk5aJRZLj5U9cmje55zpg1gJeHO5NG92zwvVlZWQQGhrF9xy4qK6vo27eXpsWbQkN7aHZvV6GGQazEz88PPz8/rcOwqhkjZ3Dm4BlAvxjGuHKxKnkyJUU+Lrly0c/XFz9f34ZPRP8Q8ZXr+hFsWIgUHODNK9f1a3C8uqCggLi4OG699Xq7Kd7k69sCX98Wmsbg7FTPWrloubm57N+/Hx8fH8JDQwF467k2ameYRrgmItish4m1ZWdn4+Xlia+vT63FMM79SU5RPWurycjIICOj4fKijmzTpk0AxAwaVDMzASBj/34y9u/XKixN2aLtH3/8MRMm3IJOp6uZtqdF8abasrKOkJV1RNMYnJ1K1laSk5NDTk6O1mFYlTFZx0ZHn3M8JzeXnNxcLULSnLXbXl5ezrp1a4mMjLCr4k15eSfJy1Nbt1mTStbKRatJ1jExGkfiOr766iv8/f0YPXqY5sWbFNtSyVq5KKWlpTXFm2IGDWr4DUqTlZaW8u67/2Xy5ImA9sWbFNtSyVq5KFu3bqWyspKwfv1o2VLt92cLR48exc1NMGbMCEAthnE1KlkrF8XUeLViPZs3b1bFm1yYStbKRVHj1ba3Zs2PPPjgXcDZ4k2RkRGqeJOLUPOsrcTXzIURjkin09XsDHPZkCEXvO7MbW/I+W2/2BKo5zt06BCJiYm88sozAMRv3Q7Yz3i1j4+aT29tKllbSVSU7XeYtgUpJbt376agoIBOHTvSudOFU8aiGrHk2tnUbntTSqCeb+PGjXTrdgnBwfrd0+1tvDosrOEl8krTqGStmM24M0yfQ/ox0tjoaKb+PhV/L3+mD5uucXT2p6klUI1KS0uZNWsm06c/BaAv3pS2GyEEgwap4k2uQo1ZK2apvTPM3KVzATjRKo8FSQt5+4N8/FrqaNu5jKmz1cIIo6aUQK1t69attGzpy9ixlwPUFG/q0ycEf3/nrj+jnNVgz1oI4QX8DXgazl8mpXxZCDEbuBY4A+wH7pdSFtTx/oNAEVANVEkpowzHg4ElhtfuklIWCyF6AR8BAYb7bZBSPtzENmoiLi4OcKwNc+tTe2eYm1+fyzEgMW0jleJGUhNmAm41O8M0E39wxVBvl9wwN84wlj9y6NAmlUCt7auvljB48EDc3PR9K3sp3lRbfLx+o2BH3jDX3nOdOT3rCmCUlDIC6A+ME0LEoN8SJFRKGQ5kAtPqucblUsr+xuANngKeBBYDEwzH5gFzDOf2AeabEZ9iI0IIpgc/wLAquBWYnQh/JazkCF34iavpwV5kVTO2JajZCaAvgerl4X7OMXNLoBpVVlaybdt2xo4dWXNMFW+yGrvOdQ0ma6lXbPjRw/AlpZTrpJRVhuPxQGOLE7gDOsOXca1sEFBTDUZKmdbIaypWJKXk399Pow/wcBtYEgb7WvjjQynD2UA+bQCorlCPQuDiS6DWtnz5cqqrKxk4MAzAroo3ORt7z3Vm/a8SQrgDyUAP4AMpZcJ5pzwAfGfi7RJYJ4SQwEdSyo8Nx98HvgQKgTsNx+YA64UQm4F1wGd1fdwAKCoqqhlqOF9ISAjBwfr/ENnZ2WRmZppsW+1hiqSkJIqLi+s8LygoiF69etXcOzk52eQ1a288kJGRYbKgk6+v7zmzRky1B7Rv08CBA3lp80us+eNnAAaPu4slLQL58vpt+O/vw6U7O9Gn9R6Ki3xx96zC+E/LOCxQZ5suvZTgQP3shuzcXDLrqVZXe0glKTXVdJsCA+l16aX6NhUXk5yaavKakeHhNbWnM/bvN1mAydfX95xZHvW1qai4uOaa2bm5+JYc5IUY4/97oOQgcZsPmtUmKSVr1q3jkcf+TfPmzQHYsTODUwWFtGnTmuzsk+TknDrnPaGhPWrqSmdlHTFZXMnHx/ucGRzGYYy6dOvWkcBA/S/i3Nx8Dhw4esE5R45c+GeXlraXkpK6x+fbt29N9+76nFdcXEp6+j6T97dGm0yxx1xnZNYDRilltZSyP/rfKIOFEKG1GvcCUAV8ZeLtsVLKgcBVwONCiOGGax6SUg6XUl4rpSwyHPsM6AMsBUYC8UIIT3NiVKzLuDNM4Cl9cg3t25eooGGQE0mhRxu2tdb/ghLu1QyMrtQyVKeg0+nYsHEjUuoYGHn2F0XSNn0C6tu3lyreZAX2nOuEcdNNcwkhXgZKpJTvCCHuBSYCo6WUpWa8dwZQLKV8x8x7pQP3Sikv6PJFRUXJpKSkRsVuS872gBH0U8j8/f3R6XScOnyYli1b8vzb+XwyvwUVhV54+pfz4FNljIvdA+DyDxgvVnV1NQ899CC7d+/k+eefYOTI2JrXHnvmJb5Z9iMzZ77AxIn3NTVci3G0B4ytWvVMPm9c+QL2kuuMGuxZCyHaCSECDN97A1cAe4QQ44DngetMBS+E8BFC+Bm/B8YA6fXca5wQwsPwfQegDXDhZy7F5qSUbN26laqqKsJDQ2u2LHvruTac+MebotOCE/94u9QWXtaQlJTEsGHDOHr0H1au/PycRA32txjGmdh7rjNnzDoI+MIwluMGfC+lXCOE2Id+yslvho9j8VLKiYZpKoullOOBQGCl4fVmwNdSyl/rudcYYK4Qotzw8xQp5TEzYrQ7ISEhWodgMcbFMO2S2wEwtIHFMCGGMWNX1JS2JyQkcMcdt/N///cot9xyLV5e534qzjueT9bBf2jRwpvQ0N5NDdWiunXrqHUIlmDXua7RwyD2wt6HQZyFlJLJayczN2Eu//6qPcl78+jy1CjW+6zHfcdDVP88G0//Ch58qkz1qpvgzJkz9OhxKXPnvsHll8fWec6aX9dz98P/x7BhMaxe/aWNI3Qu5gyD2Bs1x0qpV81imCodL8+YTyugaOF6tvp1JP5kG+KJI7FwMAtmtgROqoR9Eaqqqnj44YeJjIwwmajBPhfDKLajlptbSXZ2NtnZ2VqHYRFCCB5ueyc/AvsF+FXCqJNHmc4sVnMD1/IjsqoZi+fpV+Zl5+aS7aJ7MF5M2998802ysw+zcOHb9Z5nz+PVubn55Obmax2GU1M9aysxzoM2zo12ZFJKnvz6BdYDzXrDJVdAdNwYhqT1JJp44tHXtK4o9AKomS9tnEPtShrb9m3btvHpp5/w6afvXTBGXZu9F28yzr02zsdWLE8la6VexjHr9X+tB+C/97/NS3F7WX3NIlYHXgq/r8e4KMvTvxxQdY3NpdPpmDp1Ko89dh8DB9Y/5c1YvKlv316qeJOLUsMgSr2Mi2H88/wBuCwmhkcGTofkR6DcH2OiFs2qePCpxlWTcwbr16+nU6dOtG7dmjZt2vDtd9/yz5EjDb5v9+7d3Hnnneh0Z7jrrlsaPN84BBITY39DIIptqGStNGhi74kUHivEx8eH8NBQ3nquDY8Nnopn2mRA4ulfxmPTTrvcw8W4uDgeeOABFi/+LxkZm0lL+4s2rf1ZvnwZ5eXldb4nPT2dCRMmcPnlIwkKasUXX8zHx6dFg/dSO5krahhEqZeUsma/xZhBg3B311eRe+u5Nrz1nPEsb1xp+KO8vJyJEx9h/fr1fPTRO8TE6GeAeXh4cNll0eQcy2Py5Mnk5GSzYcNGlixZwhVXXMGjj07kxx/XcP/9d5CYuJaWLc0bztDpdKpnrahkrZhmXAwjNuqHOmJjYtTOMMA999xNs2awYcNq/P1bnvOaEILhw2P489e1nDx5nDOVlXTu3JknnnicEyeOkZDwCwEB/o26X+a+AxQUniYoKJDOnZ1i8YlyEVSyVupUe2eYwF/0Mxv2NE/n56Q1uO94iJlX61xyMUxlZSV//LGe1NQ4k8MXHQLbM3v2y4SFjWTjhg2UlJSwdOky4uN/pry8gr/+2kxERD+zk3btIRBVvMl1qWRtJY5ewKlmMUxhKbdkLCIBSNq+hran7+HEX+8AomZnmPMXwzhzAadPP/2UwYMH0qJF3cM+MQP1y8BLS/UPWyMiIrjkki7Mm/cmRUUl3HHHIxzLzePGG8Yzb94bZt3TERbDOEoBJ0emkrVikhCCe3P6MgC4DOBHgP+xj40kEMNHPMz2qoEsnufNrCmahmozp0+fJiKib4M93GbN9GP7ycnJjBkzhjlzPubUqVMEBLTCz8+X6dMnmX1Pe14Mo9iOStaKSVJKnt36Dc2A4R0g1hOijjSjR3UWPchiKfopZ8bFMM5uz549fPjhQubMeaXBc5s3b878+W9yyy0388EHC7jzzruoqqri3nvv5fvvF9GhQ3uz7mnPxZsU21LJ2kqMRaZq7wTjSGoWw2xNBODeNz7hhuXxcMciwv68hZjEQSQyGLhwMUySYXeW2rurOIO33nqLW265huHDTQ/zpO05CEBY766MGTOSZs2aMWnSU1x22TBGjx7N0KFRjUq6xi28IiMj8PCw370t09L2ApyzU4tiWSpZW4mpbacchRAC/+b+eOZ4UkGFfjHMgVEsSHQjpXkrUnhUf17NYpizydrR216XyspKVq5cyY4df9Z7XknpufOrR426jA0bVnPffZOYM2cOffs2Lpk5ypQ9U9t3KZajFsUoJt3a/lYqiivo1LEjXTp3dunFMJ9++in9+4ddMFXPHF5entxyy9WcPHmSDRviSUjYZvZ71Xi1YqSStWKScTFMbHR0zTFX3Blm3rx5PDtlCo8/fv9FX+Oqq0YjBLzwwgs89NAzHD7c8JL02sWboqL6X/S9FeegkrViUk2yjonROBLt7N+/n3nz5vLN1wu58soRF30dH58W+Pr6MGrUKK6++moWLTK15+pZ23bsVMWblBoqWSsmGZP1ZUOGaByJNqSUPPPMM9x223VcdlnTfmGlp+8hP/8k//d//0dpaSkrV/7EkSP11ztP2JoC2Pf8asV21ANGpU7Hjh0jKysLX19fwvr10zocTWzdupWkpK3MmfOfJl9r7txFBAUF88f69Xh7edGuXVv27z9Ip06m652r4k1KbSpZW0lQUJDWITRJ7eJNzZo17p9JkBNsOiCl5LXXXuWhhybQunUrs9/Xvm3ABce2bEli48YEXnvtNX7++QcKC09z6NA/DBhgemqjTqermbZn7zNBANq3d/5nF1pTydpKevXqpXUITVLXw0Vz9XKC3c1nzpxJTk42DzzweqPe171Lh3N+1u+v+H+8+eabdOnShT/++JtNm36iXbs2eHo2N3kdRyve1L17J61DcHoqWSt1cuWHi7Nnz2bZsu/57ruPzS5jaspPP/1GRcUZnpo0CS9PT66+egydOjX8qUsVb1LOp5K1lRQVFQHg5+d4T/FLS0vZtm0bbm5uxAwa1Oj3FxkWxfj5+lo6NKv766+/WLToY9as+ZJ27Rq/n2CxoYCTr6HQ0zXXjCElZRfBwe159bV3a2qGNMQRijfVVlxcCoCvb8MbKSgXRyVrK0lOTgYcs/peYmIiVVVV9A8Pv6hfNsmG5eaOVn0vOzubJ598gsmTH6F9+3YXdY30PYeAs9X33N3defnl/wPg229/QKeTZl3n7MpFxyhXkJ6+D1DV96xJTd1TzlF7Z5jY6GikNC+5OIO7757AzTeP51//utEq17/88lh0Ol2D5xmLN/n4tFDFm5Qaqmet1DDuDJO5KROAoS60M8z8+fM5diyHBx+cYLV73HnnTezaldngecZedWRkRKNn4ijOS/1LUIBaO8NsmYvnBk8ANlT+yZKkJU6/M8z8+fP5+OMP+f77xWZtXnuxune/hO7dL2nwPEcbr1ZsQyVrBTi7M0zzLfu5ongNmz0hOW4Jvnn3U7y+/p1hHJkxUX/77cd07Nih4TfYwNmZIPY/v1qxHTVmrdQQQnD1vpaMAWZUwI/fwpH1X7CB4czmWWLZhKxqxuJ5zrGT+fLly1m48AO7StRl5eXsSFfFm5QLqZ61UkNKybR/EgkGhneBoZUQcQz6yx30Zwe5BLKJWKfYGebYsWO89tpr/N//TbSbRA1nizf169dbFW9SzqGStZVERjrWR1jjzjBbDFOwpi/5m1GzvsSzzyIi/7yOmG3h/Mo44MKdYc4Xaec7xOTn5/Pvf/+bwYPDufnmay167dDeDY9J18dRizeFhvbQOgSnp5K1lTjaYhghBO4l7nAKfH19CQ8N5ZGB01mQ6MaGFq3YwHP68+rYGeZ89r4Y5vnnn6ew8CSzZj1v8Wv7mtj13FyOWrxJLYaxPpWslRpDpL4UqrF401vPtYG3p/LJ/BZUIPH0L3fo2SCFhYV89NFHbNq0gRUrPiMoyL4KTp1bvMkxFsMotqOStZVkZGQAjlPQqa7FMEII3nquDW89ZzzLm/p61EYZ+/cD9lXQqbi4mPvuu4+cnCN8++1HVkvUWYePARcWdDKHsXhTcHAgnTubLp1qj7Ky9DvfqIJO1qOStZXk5OQAjpGsjYthNm/aDMDQ6OgmLYbJyc0F7CdZ5+fnc+211+Dp2YwvvphPYODFLSU3R96JAuDiknXtKXuOVrwpL+8koJK1Nampey6uZjHM33NJ3paMm5sbP5esYUHSQn7alEubzqX4tZS07VzG1NkntQ630YyJOjY2ihUrPrNqom6q+K3bAccbr1ZsQyVrF2dcDPNS/nAyq3V83lyHbsknXH/8BtLeepszhS04uyCmpUMl7NqJevr0SXbfW1WLYZT6qGStIIRgRFZzLgXuLod5v8KSxas4VN2TZdzMXeg3d3WkBTGOlqjzjudz4JAq3qSYppK1gpSSh05lMgB4OgyW9oF/PH1pRQFj+Y3+bK851xEWxJw8edKhEjWo4k1Kw9S/ChcnpeTpX54ma+dhAH74dicf7HmfB5IW0mnDBGI2Xs5+zj4obGhBjD2YP38+HTq0dZhEDap4k9IwlaytxNfOF4YYCSGozK2EcujcqRNdOndmVqdZbEwuJ9WtHcu49ey5ZiyIAW3b/r///Y+vvlrCN998qEmi9mlxcZ88HH282sfHvn+BOwOVrK0kKspxFjVEnIkAzm6OK4Rg4/PvMc39FIt3lFFR6NWoBTFRNlxufujQIf78809uv/12vvvuO2bOfJNvv/3IrFKk1hDWu2uj3+MMxZvCwnpqHYLTU8nahRkXvmzcuBHQbzZgJIRg1pTWzJpiPGLeghhbOnToEOPGjcXPz4cZM17Gy8uTb7/9iG7dumgdWqMYizeFhqriTYppKlm7KONCmDlj59SsXEwmgTc35DvMrjCTJj3FLbdcwzPPTCQnJ49Wrfzx9rb/B6DnO1u8yTGHQBTbUMnaSuLi4gD73DC3ZiFMwlxKTpZw4MABPLw8WHZiGY+XP1rT475YcZv1KyGttWFufn4+7733HocOHWTOnJcQQhAcbB91PuK37QHObphrDkct3lRbfLx+k2S1Ya71qKl7Lsi4EGZS9CQOfLSYRcDdPpW83OlfzBo9k2nvnKJt5zK7W7kopeTdd99l6NAhbNuWyLfffkSrVgFah9UkOp2u1kwQ1bNWTFM9axdlTNjv3T+XB4EH84Env6H4+XVcVhZNMzmEvxnG9sKBdrOV1yuvzODXX39mwYK3iIx0jh5cxt4DFJ4ucsjiTYptqWTtooybDfxRDf8AwzvB6CJv/ArzuYafuYaf+YVx3Mb3NSsXzz5stL2cnBw++GABa9d+p9lMD2uo3at2lDnhijZUsnZBxkQ9d8Nc3E66scsNdLPu5+70T+j+97+I2TycGBJIZHDNe7ReuTh//nyuumq0UyVqUIthFPOpZO2ChBAEeAVws//NLK9ezoCICOZc+y6e3s35KKkFWdzF19x1znu0WrkopeSTTz5h4cKF/PbbUpvf39ocfTGMYjvqAaOLmjFyBv3P9Af0i2GEEMy6YhYThzyJaFZ1zrlnVy7alk6n4/XXX2fhwg/48MPZ9OjRzeYxWJMq3qQ0hupZW0lISIjWITSoZmcYw2IY40IYOMnied6NXrloFGKBTQe+/fZbXn75Jfz9/fjuu49p29YxthLr1ohNB4xDIFFR/R2+eFO3bh21DsHpOfa/EDsWHGzfT/Z1Oh1btmwBziZro6auXAwOvPg5z1JK3njjDZYu/Y7//OcZrrnmyou+lhYC2waYfa4zjVcHBrbROgSnp5K1i9q5cyeFhYV07tSJzp3sZyum999/n1WrVrB8+acO05u+WM6wGEaxnQbHrIUQXkKIRCHEDiHETiHEK4bjs4UQe4QQqUKIlUKIABPvDxBCLDOcu1sIMcRwPFgIsV4I8YMQwtdwrJcQIk4IkWI492MLttWmsrOzyc7O1joMk2pvjmtp2bm5ZBv2YWyMgoIC5s59j9dff95hE3XuiQJyDfsw1qe0rIyUtN24ubk5bPGm2nJz88nNzdc6jCax91xnzgPGCmCUlDIC6A+ME0LEAL8BoVLKcCATmGbi/XOBX6WUvYEIYLfh+FPAk8BiYILh2DxgjpSyv5SyDzDfjPjsUmZmJpmZmVqHYdL549WWlLl/P5mGHc4bY8+ePbRt25qYGMepWHi+A4ePccCww3l9tu/YSVVVFX37htCypeMXbzpw4CgHDhzVOoymsutc12CylnrFhh89DF9SSrlOSmmcNhAPXPBZWgjREhgOfGK41hkpZYHhZXdAZ/gyrgYIAo7UundaQ/EpF8eayfpiLVr0MWPGjNQ6DJuIV1P27I695zqzxqyFEO5AMtAD+EBKmXDeKQ8A39Xx1u7AceAzIUSE4RqTpJQlwPvAl0AhcKfh/DnAeiHEZmAd8FmtBp+jqKiopljS+UJCQmoe8GVnZ9fbw61daCkpKYni4uI6zwsKCqJXr141905OTjZ5zcjIs/8BMzIyyMnJqfM8X1/fc+pem2oPWLZN+fn5HDhwAB8fH8L69dO3qbiY5NRU020KD8fPsKlAxv795JgY5jh/4wFjUac623TppTUPI/cfPEhJVTUdu11aUwypttqFkdL2HKSktLzOa7ZvG0B3w4yM4tIy0vccMnn/0N6X4NtC//A06/Ax8kwMX/i08DqnTnVd8RmVlJbXbECQe6LAZC977R/6srTG8eq0tL2UlNQ9PbJ9+9Z0767PD8XFpaSn7zPdptAe+Pq20Lcp6wh5eXXXdfHx8T6nBrWxEFNdunXrWPMAMTc3v84e9JEjF/57sOc2mWKPuc7IrHnWUspqKWV/9L9RBgshQms17gWgCgy7qp6rGTAQWCilHACUAFMN1zwkpRwupbxWSllkOPYZ0AdYCowE4oUQnubEqJgvPT0dgIiwMLuZMrZ//37c3QRdLrGfh53WotPp2L1nL6B61vbGnnOdkFI2qjFCiJeBEinlO0KIe4GJwGgpZWkd53YA4qWUXQ0/DwOmSimvNvNe6cC9UsoLurFRUVEyKSmpUbHbkj2XSJ08eTLvvfceM6ZN4+VppobfLt7FlEi955676dSpLc8994TF47Elc0qk7s7Yz9ArbyE4OJD09A1OURPE0UqktmrVM1lKWe/DEXvJdUbmzAZpZ3z6KYTwBq4A9gghxgHPA9fVFTyAlPIY8I8Qopfh0GhgVz33GieE8DB83wFoAzj8Uwt7YfzFbByvHmqFmSAXIzMzkw0bNnDvvbdrHYpNqOJN9snec505n4GDgC8MYzluwPdSyjVCiH2AJ/Cb4R9cvJRyohAiGFgspRxveP+TwFdCiOZAFnB/PfcaA8wVQhgHJKcY/hCUJjLuDPP6Za+zbds23Nzc+Kn4R7ZuiNd8Z5hNmzbRr19vAgPbaRqHrRiTdUyMGgKxM3ad6xo9DGIv7H0YxJ7UVNlLmMtNXjexYuoK2ndrS+kdJwgX97Png5mcKfS+qKXlTVVSUkJk5EDeeGMao0cPs9l9tTRw2HUcOPQPcXGriIjop3U4LsmcYRB7Yx9PlxSrMm40AFA2aS5fAPHVJzh9/AZWf/Y2slr/XKOi0NvmGw0sWbKE4OBALr881ib301pu3oma4k39+vVq+A2KYqCq7rkIY8K+4QjcAyw4DEsWr+JQdU+WcTPPMpu2nKjZaMAWUlJSeOutWbz66vO4ubnGP0VnKt6k2JZr/A/RQFJSEvY0TCOl5Olfnua5SngCWNoT/vH0pRUFjOU3XuY1BDqg6RsNJKWmklTPnG2jLVu2MGxYtFOVB03bc5C0PQdNvu5MxZtqS0vbS1raXq3DcGrqV7uVmFpcowXjmPW8NfOgEgo7dcJn1rU8kLSQThsmELPxcnqRwXHaA03faMCctm/fvp233nqLRYv+e9H3sUemFusYOWvxJlOLXxTLUcnaBRh3hrnc/XL+5E8ui4lh1hWz2JhcTqpbO5Zx69lzazYasO5QSEpKCmPGjGDQoP5WvY89KS0rY0f6Hqcp3qTYlhoGcREzRs6gY4G+QHxsTAxCCDY+/x6Pj3gMT/8yQOLpX8Zj007b5OGiK84vdrbiTYptqZ61C5BSIoS4YDGMcWeYpmw0cLGEEDjqtNGLpYo3KU2hetZObkbcDCavnczRo0c5cOAAfn5+fJ2zhDc3vKlZTFJKEhMT8fa2/Qa8WlKLYZSmUD1rJyalpKC8gLkJc8mKywCgzLMlH27/CPcdD1G4JZ+3nrP9dkwrVqwgMTGepUsX2fzeWtHpdCSqnrXSBCpZW0lQUJDWIZyzGKbr/XOZB2wRR0lcMYxVGVNZuM4fISy/ACaogT0YZ86cyZNP3k9AgL9F72sP2pvYgzFj7wEKTxcRHNyBzp3te3/Oi9G+vWPu7ONIVLK2EmPta60ZE/ZfJ+fSFeh6HP51fAMbOcGJqnYsnudda8zaMno1sLt5dvZRLr/cOZeWdzexu3n81u2A8/aqjbWpFetRydrJSSl54ocn+PAM9ANGjYYBu3uTkR0CNH0BzMWortbRrJm7ze+rpbPj1c41v1qxHfWA0UqKioooKirSNAbjYpgFqxagA3YSzifiUZ64dw/yiunop+vVv4jjYhQVF1NkYmHMjz/+SIcO7fH0bG7x+9qD4tIyiksvXCCS4OTj1cXFpRQX11k9VLEQlaytJDk5ud6tv2zBuBhmiBwCgLuIht9nwdZHodwf0azasADGspJTU01uERYfH8/YsZfj7u6cPev0PYcu2EosN+8EBw8fceriTenp++rdmktpOpWsndyMkTPwz9M/yBt3Q5i+J/37TDzTJttsAYzRqVOnWLBgAbfffr3N7mkPVPEmxRLUvxwnNfmNkyx4x4szBZ7AFgA+eOMKOncyzm223QIYo48++oirrx7DpZd2tel9teasxZsU21I9ayc0+Y2TvLf+Pc7ETAd2AoUIOnHHktmaLYaRUvL7778xZIhzjtnWxzherRbDKE2hkrUTWvCOFzQ/zT1ec/nVZyQvAuFtBKnyMwrLCzVZ5q3T6di6NYl//etGm99bS6p4k2IpKlk7oTMF3vDrHK5M6MXYklO8Boxo9w9sncisK2ZpUkSpqqoKnZQus8mAkbF4U79+vfDz89U6HMWBudb/HBfhGVAOCJ7L+YPbWsBc4K++0HzrDM2q3b355puMHuWcC2Hqo4o3KZaiHjBaSWSkdv85H322jPdeaUZO7MssjYOlzcGnJ0T0noaU71k9YUeGh5/z8+nTp1m8eDHff/+xVe9rD0J7X3LOz2dXLjr3w8XQ0B5ah+D0VLK2Ej8/7eoVvzu9FXHVj5KS/gkAzXQjiHDvTqr8jKm/e1l9KMTP99yP+++8M5vY2MH06+c823eZ4tvi7AwbnU7H1mT9fHNn71n7+rbQOgSnp4ZBnJAQguuHd6D/mf4AvDhlKBuff4/Hoh7F38vfpkMhxnKo48ZdbrN72os9e7OcuniTYlsqWVtJRkYGGRkZmt1/xsgZNDuq/+Bk3Blm1hWzmD5sutXvnbF/Pxn79wNw7NgxNm3ewpAhUVa/rz3IOnyMrMPHAEjYmgI4f68aICvrCFlZR7QOw6mpZG0lOTk55OTkaHb/kpIStm/fjpubG9FR+kRpqx51Tm4uObm5AHTo0IGrxo3lzz832eTeWss7UUDeiQLAtYo35eWdJC/vpNZhODWVrJ1UYmIi1dXV9A8P13T8XAhBYGAgOp1Osxi04uzFmxTbUsnaSRn3W4w17LeopcLC0y63Qa6xeJOvr4/TFm9SbEslaye1ceNGQD9eraXjx4/z+++/ceWVwzWNw9bOFm+KUMWbFItQydoJVVdXs2WLvniT1sn666+/pnv3rrRv307TOGztbPEmNQSiWIb6le+Edu7cyenTp+nSuTOdOnbULI7FixexePHHfPPNR5rFoJWz49XO/3BRsQ2VrK3E11e7OhBaj1f7+vpSUV7O7NmzmTXrRTp10n7zYFvxaeFFeXmFyxVv8vGxbbldV6SGQawkKiqKqCjrzC2e/MZJvFqVIYTEq1UZk984d8pUTbLWaAgkKjwcD8DfvyUjRw7VJAathPXuSkVZkcsVbwoL60lYWE+tw3BqKlk7mMlvnGTuKy1xK5CAoKLAm7mvtGTyGydrSp8ak/VQDWeCfPjhQi6/PNblZoEAxKvxasUK1DCIg1n4jjfBlcc4QDdS6M8mYtlUGcuy3zcgI92YEjaFgwcP4ufnx9c5S/jl9E82WbV4vuPHj3PjjaNtfl97oHaGUaxBJWsriYuLA2DkyJEWvW5FgRehpCOQDCKJQSTxNHMhDrJSIKPTUh4Bdnf04cPtH/FY1KNIKW3aw43bvJmgrl1tdj97sjlpF/GJKYBr9azj4/UFq2Jiwhs4U7lYKlk7GM+ActYWjCOAAqJJIFbft2YIW+heUET3gmxGAa+UHiMm6lFNNhs4fvw4x/OO061bV5ve1x4c/ucoxSUldOwYpIo3KRalxqwdzKPPliE8qijBl/WM5jVe4iqPn5jx6n7ktm087g9fA3EhaLYrzOmiItq3b+tSs0CMdu3RF+9SQyCKpalk7WDmvNCaSS+fxjOgDJB4BpQx6eXT/PfFtjx+cBELTsNdArZGwNTfp2qy36KbEJrc1x7s2p0JuNYQiGIbahjEAc15oTVzXjD+5I2UXkxeO5mFqxaChIEREVwxdCgLkhYCtu9hCxdO1jt363vWrlBpT7EtlaydgBCCAK8AYmQM8cQTGxPDrCtmAdh8swFjPK6Yq4/lHudYbh7e3l707auKNymWpYZB7FxDC2CMZoycgX+eP2D7zQbOJ9zcXLJnbZyy17t3iCrepFic+hdlJSEhIU2+xuQ3TrJiRgGvVs3Sz/koiGXuKwHASea80LrmPCklOp2upniTcTGMVgtS/jlwgIqS05rcW0vGeiBDhw7SOBLb69ZNuxo0rkIlaysJDm76tK2F73hzW9UanmM2MBuAPZW9iH81GoJGQGwsM7K/pqCikPuD7q8p3vT+7vn4H/DXpFcN8MuaH/nPf57W5N5aMvasR4++TNtANBAY2EbrEJyeStZ2rKLAi+0M4HVeIJZNRJNAbzLofSYD/v0/AJ5t6c0fQWXEt/6RoUBx5+YsSFqoyWIYgPj4ePLzT9C3b9M/WTiS0rIyUndmuFTxJsW2VLK2kuzsbKBpPWzPgHLSC8JIJwyAZlTSnxRGeMXxzjWJsGkTvjk5XH8aIItHgKdO7mecRothAE6ePMnAqIGUVlSh4W5iNrctZSdVVVX06RNCaWmFyxRwMsrNzQdUD9ua1ANGK8nMzCQzM7NJ1zAugDGqwoNkjwFUv/hvWLoUjh6FrCzkF1/wYXNIA+J7a7cYBvTj5/5t23HAsMO3qzAWb+rR41IOHDiqbTAaOHDgqEu225ZUsrZjphbA1DxcFALZtSsPuf3No2cg3BN29dJuMQxAs2bNOH26iGpdtSb314pxvLpv397aBqI4LTUMYkdWbT/K7LUZZBeUERzgzZSxvZjzQsdzFsDov/SklExeO5lPfvgEgDGxo+g/uJdmi2EARo8ezRfffENu7gmb3ldLOp2Orcn6Qkb9+qlkrViHStZ2YtX2o0xbkUZZpb5HerSgjGkr0gC4YUDd06KMi2H6V/YnhRRiY2L4zxXPA9oshgF9z9rb2xup09n83lrZszeLwtNFdOwYRLt2bbUOR3FSahjETsxem1GTqI3KKquZvTaj3vfNGDkD9yPugPaLYYzatGlDZmaWZve3tfjE7YAq3qRYl0rWdiK7oKxRx42Ki4tJSUnB3d2daMM2YlrvztK3b1/yT57SNAZbOrs5rirepFiPStZ2Ijig7g1HTR03SkxMpLq6moiwME036a3N29ubEydOkpOTq3UoNmF8uKiKNynWpJK1lYwcObJRu8RMGdsLbw/3c455e7gzZWz9BYG03sm8LrffeCMt3N1YtOgrrUOxumO5xzn0z1F8fX3o27cXMTHhLrlbiqu225YaTNZCCC8hRKIQYocQYqcQ4hXD8dlCiD1CiFQhxEohREAd7+0lhEip9XVaCPG04bVgIcR6IcQPQgjfWufHGc7dLYT42LLNtV83DOjIzJvC6BjgjQA6Bngz86Ywkw8XjbTeydyUp5+ezKZNCVqHYXXGXnVUVIQq3uTg7D3XmfOvqwIYJaUsFkJ4ABuFEL8AvwHTpJRVQoi3gGnA87XfKKXMAPobgnMHjgIrDS8/BTwJdAcmAB8C84A5UsofDO8JMyM+p3HDgI4NJufaqqura4o32Vuy7tu3L8XFZSQl7SAqKkLrcKzGOF4dExOlcSSKBdh1rmuwZy31ig0/ehi+pJRynZTSuLwuHujUwKVGA/ullIcMP7sDOsOX8YlYEHCk1r3TGorPXiUlJZGUlGTVe+zcuZPTp09zSZcudOpoP1XPklJT2ZOVxeDBg/j77y1ah2NV5+9knpa2l7S0vRpGpA1naLe95zqzxqyFEO5CiBQgD/hNSnn+59sHgF8auMwdwDe1fn4f+AiYCCwxHJsDrBdC/CKEmFzXxw1HUVxcTHFxccMnNsHGjRsB+xqvhrNtnzjxUX799U+tw7Ga2sWbIiP1nx5KSsooKal/Bo8zcpZ223OuM2uQTUpZDfQ3XHClECJUSpkOIIR4AagCTD5NEkI0B65D//HBeM1DwPDz7vOZEGItMA64HnhECBEhpaw4/5pFRUXExcXVeb+QkJCaAkrZ2dn11uio/RAwKSnJZIINCgqiV69eNfdOTk42ec3IyLNTuDIyMsjJyanzPF9fX6Kizn58NtUeqLtNK1asAKBd+/bEbd58tk1Dh55tU2qq6TYFBtLr0kv1bSouJjk11XSbwsPxM8w2ydi/n5zcumd61J6REhISQucePflmxTq6de1ywbndunQgsG0AALknCuqtJxIz8OzKwLQ9BykpLa/zvPZtA+jepQMAxaVlpO85VOd5AKG9L8G3hX62TdbhY+SdKKjzPJ8WXoT17lrzc/y2PQCkpu2iqqqKSy/txs6d+nnlJSVl+Pjor5mbm19vvYzaD+TS0vaaTHbt27eme3d9Z664uJT09H2m2xTaA1/fFvo2ZR0hL6/uzSp8fLwJC+t5tk3xpv/uu3XrWFOgyVSbjhy58N+DPbfJFHvMdUaNmg0ipSwA4gw3QAhxL3ANcJesvxjFVcA2KWWDc7mklNlSyk+llNej/4MJbUyMriQ9PR2A0H79NI6kbq1ateKy2Mv4668tlJSUaB2Oxe007GSu6oE4H3vMdaKhgj9CiHZApZSyQAjhDawD3jJc/F1ghJTyeAPX+BZYK6X8rIHzxgF/SCkrhRAdgO3AACnlBV2uqKgoae0x4aYw9pIbM32vMY4ePUqnTp1o2bIlJw8dwt3dveE32Yixl2/s4c+aNYuff/6R1av/p2VYFnfLPY/zR9xmFi+ew803XwOc7c252jQ2R2t3q1Y9k6WU5zwVttdcZ2ROzzoI+FMIkQpsRT+Oswb9OIwf8Jth+smHhiCChRA/1wqqBXAlsMKMe40B0oUQO4C1wJT6gndlxil7MVFRdpWo6/Lggw+yd+9+Tp8u0joUi9HpdGzdpk9QauWi07DrXNfgmLWUMhUYUMfxHibOzwbG1/q5FDCrIrmU8hngGXPOdXX2Or+6Lm3btmXEiJGsWLGG++77l9bhWMTuzP2cPl1Mx45BdOoUpHU4igXYe65Ts/itJCjIOv+BjVt1GZP1UDubCQL6B5fn69KlC1VVzlPjOmFrCnBh8ab27VvXcbbzc9V225JK1lZinDliKVfdWsyvJ2aDVwH8+ipC6Is3/VT0I4kbtmhaZe98xhkmtel/yWgQjJWYKt5knOHgaly13bakaoM4gKtuLebXZS0YWZJG+CXzYNADSFmNR6sAFu9cTGF5oWY7w7gq42KYIUPUeLViG6pnbSVFRfqHaX4W2DX21+UtADc+35rMJVshz2MFa4HFbfKJ0XBzXFOKDPO6/WrNudZip3VrMRZv8vPTF2+qrbi4FKBmXrCrcNV225LqWVtJcnJyvQtnGkXqk1wb9Ftlta+Eu4FOLbXdHNeU5NTUCxbYOFOyPlu8qf8FM3HS0/fVu8DDWblqu21JJWtHICSdOYwvpZxyb86AZvA4sClM281xG6NFixZs375T6zAsQm02oGhBJWsHMO7mUmLZAMCW4DOkVMFCOuPT7EEWJC10iIT97LPPsnLVz3Yfpznit6ptvBTbU8naAfyy1Jd/Xfo3ABvLuwHQsVMUaW+/w2NRj2q2OW5jtGzZkqjISL76arnWoTRJSemFxZsUxRZUsnYAUkoGidUAeIV3B2Dq07FM+2MaLT1b2tW0PVOEEIwcOZL8/LqL8DiKbSnpVFdXExraGz8/+9hGTXENKlnbuRlxM3h+xWMEZuVS6QbxyZvxBL44/ikLkhby1eoT+LWsxq+lxK+ljhvvtt8l3cHBwSxbtsahh0LO1q9W49WKbalkbceklBSUF7D9hw9x00myO7Xi54IyTgPzv9zF7DVd6P/xMNpxAn1Nczd+X+1JxF1LNY68bnl5efTo0dXuh2zqE3/eZgOKYisNVt2zV/Zedc9S86yllKydEMO4rxNZ0hnC/9HXUTz/t+wZPPiIh5nOTPDKp/BYK9zc3M65jq2SZF3zrAEOHTrEwIED2b8/0SZxWJpOp6Nb+AhOny4mPf1vOna8sKSAq843drR211V1z96pnrWV+Pn5WWRBjBCCsfe8yqKB8HlfiHCH1uiL7L4GrAeqgeZUUsx+dOxHlrdi2h/TaoYbpJRM/X0qb254s8nxmMPP1/eCRA3QuXNnOnQIZPXqtTaJw9KMxZs6dQquM1GDPlk5SsKyJFdtty2pFYx2TkpJzOGXSLzOcCAChnldQ2JCImt35+HxjxuFUoc38BG/UsqvCNry32knqPj4Lx685lG+bJXCgp2f8FjUo5ouTnFzc+Phhx8hNXUn1103VpMYmsJU8SZFsQWVrK0kI0O/i0hTCjpJKXn616dJzE4ECf07RDCsy3Dmb50PQ8CvfxglW8Yx6u8ehLGOQ1QDm5AcpzoT7sjcRejaJ0lxg4DeQbhVePNz9VqGRkfTqlUrC7X0Qhn79wOmCzo5KnMWw2Rl6fdAdbXCRq7abltSydpKjPsuNiVZCyFo5d2KpwY/RUjKP9zw3Eo2dd6B6AKFkf0IiB3O3L1HiPd+nviyh8D7JIx7gnuubs6oyiEMfOxRdBLidXB6Vw6zdr1bc+2+vXsTGxPDZTExxMbE0L1bN4v1uI37M9aVrA0ts8h9bM04EyQmxnTP2rg/oKslLVdtty2pZG3nZoycoe+Nbp+JKFrJbbvgtl3Arzs545XJY70C6fnMHMSwYcjoaCbHB7Lmvzdz4O8C7pWww603p6/qy7VRlfQp68PmzZvZunUru/bsYdeePSz6/HMAAtu3J9aQuGOjoxkQEUHz5s0t3h5HrRGScyzPZPEmRbEFlawdxDP9c/n1cYg9DLH/wNXHA2h/tICQHUdgxxuAvif+VMt+hBSW0YYCADbpRsIvc/nR8z66P1TB37P+5syZM2zbto2NGzeyadMmNm3aRG5eHitWr2bFasPiGy8vBkdG1iRvSw2dOGqyNg6B1FW8SVFsQSVrOyelZPLayczdOo9J10xiztg5TF47mQcS5vJCyIO85nkVYvNm2LQJkpPpXpjOY6TXvL+S5qBrDr/MJ+DJ+Qgh8PT0ZMiQIQwZMoQpU6YgpWTv3r01iXvTpk3s2bOHvzdt4m/DjjQA/fr0qUneTRk6ccRxa7UYRtGaStZ2TghBgFcAk6L1iVoIwZyxcwBo5hWAGHkT3HST/uSyMoa12Eosm3mC9wkmm++4zfBaa2aMnGHyHiEhIYSEhHD//fcDcOLECbZs2VKTvLdu3crO3bvZuXs3H3+m37j5YoZOoqOjueuuuUyd+qRD9VAT1GIYRWNqUYyVxMXFATBy5EiLXO/84QNTwwnCTQdSP33ejSp0uAFuIHRI3cVPq6+oqCA5Ofmc3veJEyfOOcfb25vBkZF06tSJ0L59eeTee+scOgkKCmLlys/p3bvOfUjtTklpGZf0G4aUkoMHk+utCRIfr6/jHRMTbqvw7IKjtdsRF8WonrWV+NaxKKQpzk/MpoYfxt1cyq/LfACBruavVzLu5lLg4mPy9PRk6NChDB06tN6hk782bqx5z7SXXqpz6GTixIl8881KXnllykXHY0vG4k0REf0aLN7k4+Nto6jsi6u225ZUsraSqChtfmn/stRXv2fj8hb6HWaEPlH/stTyvzzqGjrZvHlzTfJOSkq6YOikQ2AgQe3a4ePtQXJKOuH9euHh4WHR2CytMUMgYWE9rRyNfXLVdtuSStZO6NzELGhKj7ox2rZty3XXXcd11+mXW9Y1dHIsN5djhnnYV1x3N95eXgzs34/oqP5ER0UweGAEAQEtbRKvueLVw0XFDqgxa8VmjEMnL7/8MikpyVRXV7N3b9Y55wgh6B1yKdFREURH9Scmqj+XdOmo2XQ/c4o3KY5HjVkrNSz9gNGRmGq7cegkOjqagAAvZs58gfz8kyQkbCchIZmEhG1s357K7ox97M7Yx+eGXWUC27XVJ+9B/YmO6m/ToRNzijfV5mgP2izFVdttSypZKzZX+9NcmzatGT9+NOPHjwb0QycpKTtrkndCQjK5x0+w+pc/WP3LHwDnDJ3EDOrP4IER+Ps3vcJhXdR+i4q9UMla0URFxZk6j3t6ehIdPbAmOUop2bfvQK3kvY29e7PYFJ/MpvhkwLpDJ2onc8VeqGSt2NzYsWOZPn06zz33OB06tK/3XCEEPXt2p2fP7kyYcCtAo4ZOYgYNIHpQf8L6hlzU0Ik5xZsUxRZUslZsrm/fvtx++23Mm7eYN99s/Ga/5w+dlJdXkJKSfk7v+/yhkxbeXgzsH3rOrJOGhk5yjuVx+J9sVbxJsQsqWSuamDp1GjEx0TzzzCO0bdumSdfy8vIkJiaSmBj9UIWpoZONW5LYuEU/g8icoRPjEMigQQMcamm84pxUslY00bt3b/z8WlJWVmHxa9c3dBIfn0RCwjZSUtIuGDrp0L6toefdn+hB/dmUoB8TVw8XFXugkrWVhISEaB3COVZtP8rstRlkF5QRHODNlLG9uGFAR6vcy5y2r127Fnd3Qdu2ra0Sw/nMGTo5lneCH37+nR9+/v2c93boEGj2fbp1s86fqb1z1XbbkkrWVhIcHKx1CDVWbT/KtBVplFVWA3C0oIxpK9IArJKwzWn74cOHGT58CN7eXha/vzkaGjqJi9vMkSPZAIwde7nZ1w0MbNqQjqNy1XbbkkrWLmD22oyaRG1UVlnN7LUZVutdN6Rfv368+OILFBQUEhDgr0kMtdU1dJKRsY9ly34kMLCdxtEpClx8zUylXtnZ2WRnZ2sdBgDZBWWNOt7k+5nR9qFDh+Lv709Gxn6rxGAJvXr14IUXJjfqPbm5+eTm5lspIvvlqu22JZWsrSQzM5PMzEytwwAgOKDu8pWmjjeVuW2/5557+fzz76wSg1YOHDjKgQNHtQ7D5ly13bakkrULmDK2F94e50498/ZwZ8pYbecOjxgxgo0b4zl9ukjTOBTFEahk7QJuGNCRmTeF0THAGwF0DPBm5k1hmo1XGw0bNoxLLunKl18u1TQORXEE6gGji7hhQEfNk3NdXnvtdf71rzu4886baNUqQOtwFMVuqZ61oqnLL7+cLl26sGfPPq1DURS7ppK1orkWLVqwa1eG1mEoil1TyVrR3Esvvcz8+Z9QXFyidSiKYrfUtl6K5qSUREcP5oUXnmTEiKFah6O4AEfc1kv1rBXNCSEYPXo0//ufmhWiKKaoZK3YhauvvoYjR3K0DkNR7JZK1laSlJSEqw7TXEzbpZS4uTn+P8e0tL2kpe3VOgybc9V225KaZ20lxcXFWoegmYtpu5TSInsmaq2kxDr1Vuydq7bblhy/K6M4hQ4dOrB37362bUvVOhRFsUsqWSt2ISQkhP/+912mTHlV61AUxS6pZK3YjWHDhnHw4GH27z+odSiKYndUslbsRrdu3RgzZgxbt27XOhRFsTsqWSt2xxlmhSiKpanZIFYSFBSkdQiaaUrbvb29SU/P4LbbLBiQDbVvb5sNgO2Nq7bblhpcbi6E8AL+BjzRJ/dlUsqXhRCzgWuBM8B+4H4pZUEd758MPAhIIM1wXrkQIhhYAhQBd0kpi4UQvYCPgADD/TZIKR+uKy613Nw5paenM378OFJT47QORXFidS03t9dcZ2TO580KYJSUMgLoD4wTQsQAvwGhUspwIBOYVkfwHYGngCgpZSjgDtxhePkp4ElgMTDBcGweMEdK2V9K2QeYb0Z8ihP5/PPPGDLEoUo2KM7DrnNdg8la6hlXOXgYvqSUcp2UsspwPB7oZOISzQBvIUQzoAVg3EnVHdAZvoyrIYKAI7XundZQfPaqqKiIoiLX3K6qKW0/evQoo0cPs3BEtlNcXEpxcanWYdicM7Tb3nOdWWPWQgh3IBnoAXwgpUw475QHgAt2PpVSHhVCvAMcBsqAdVLKdYaX3we+BAqBOw3H5gDrhRCbgXXAZ3V93AB9QoiLi6sz3pCQEIKDgwH9Ttv1bd46cuTImu+TkpJMrr4LCgqiV69eNfdOTk42ec3IyMia14OCgsjJqbvmha+vL1FRZ3uRptoD9tEmPz8/ADIyMuptk/F+I0eObHSb2rZtS1FRBfHx5y6OiYkJr/k+LW2vyRVz7du3pnt3/f+l4uJS0tNNb2oQGtoDX98WAGRlHSEv72Sd5/n4eBMW1rPm5/Njq62kpAwfH29iYsLJzc2vdxNZR2lTt24dCQxsA2CyTUeO5AJwyy1XOkSbTLHHXGdk1mN3KWW1lLI/+t8og4UQobUa9wJQBXx1/vuEEK2A64FuQDDgI4SYYLjmISnlcCnltVLKIsOxz4A+wFJgJBAvhPA0J0bFOegfoTj+snPFMdlzrmt0PWshxMtAiZTyHSHEvcBEYLSU8oLPQEKIW4FxUsp/G36+B4iRUj5m5r3SgXullBd0+ez9AaOxR1m7l+sqmtL222+/jSuvHMItt1xr2aBsxNibq91rdgWO1m5z6lnbS64zarBnLYRoJ4QIMHzvDVwB7BFCjAOeB66rK3iDw0CMEKKF0FfpGQ3srude44QQHobvOwBtANOfIxWn46ibYSiOz95znTlj1kHAF4axHDfgeynlGiHEPvRTTn4zVEuLl1JONExTWSylHC+lTBBCLAO2of/4sB34uJ57jQHmCiHKDT9PkVIeMyNGxYmo7b0Ujdh1rmswWUspU4EBdRzvYeL8bGB8rZ9fBl5u6D6Gc58BnjHnXMU53XXXBKZM+T/uu++Ohk9WFAuy91yn1vUqdqV///6Ul5c3fKKiuBi13NxKIiMjtQ5BM67c9tDQOjthTs9V221LKllbiXFOsity5bYb5wO7Gldtty2pYRBFURQHoJK1lWRkZJCRkaF1GJpw5bZnZR0hK+tIwyc6GVdtty2pZG0lOTk5JpdkO7umtN3Dw4PTp4s4fjzfwlHZRl7eSZPLoZ2Zq7bbllSyVuxKcHAwo0aNZs2adQ2frCguRCVrxe60ahWAu7u71mEoil1RyVqxO9XV1Zw8eUrrMBTFrqhkrdidhx9+hIULv9A6DEWxKypZK3anb9++lJSUUFFRoXUoimI3VLK2El9fX3x9fbUOQxNNbbu/vz8xMTF89dVyC0ZlGz4+3vj4eGsdhs25arttSa1gtJLaO8C4mqa2XQhBWFg4ZWWOVyOk9k4lrsRV221LqmetKIriAFSyVhRFcQAqWVtJXFxcvZvFOjNXbnt8fOpFbdTq6Fy13bakkrVilzw9Pdm9e6/WYSiK3VDJWrFLU6ZM4auvlnH4sCoOpCigkrVip9q0acPDDz/MW2+9r3UoimIXVLJW7NbLL7/M118v58gR16xeqCi1qWSt2K3g4GAefPBBvvpqmdahKIrmVLJW7Npll13G6tVrqa6u1joURdGUWsFoJSEhIVqHoBlLtv2ee+7hqaeeJDf3OMHBHSx2XWvp1q2j1iFowlXbbUuqZ20lwcHBBAcHax2GJizZdiEEt912Ox9++D+LXM/aAgPbEBjYRuswbM5V221LKlkrdm/48OHExW3izJkzWoeiKJpRydpKsrOzyc7O1joMTVi67RMmTKC0tJzly9dY7JrWkpubT26uY+4f2RSu2m5bUmPWVpKZmQngkkMhlm67EILly1cwdOhQRowYatdj1wcOHAVwuSEBV223LameteIQIiIiuOWWm5k7d5HWoSiKJlSyVhxGjx49KSgo1DoMRdGEStaKw9i1ayetWgVoHYaiaEIla8VhTJs2ne++W0Vu7nGtQ1EUm1PJWnEY4eHhXHHFFbz//idah6IoNqeSteJQ7rvvfjIzD2gdhqLYnJBSah3DRRFCJGkdg6IoDuuElHKc1kE0hsMma0VRFFeihkEURVEcgErWiqIoDkAla0VRFAegkvV5hBCfCiHyhBDptY7dKoTYKYTQCSGiah3vKoQoE0KkGL4+NHHNGUKIo7XOG1/rtdlCiCQhxAgT10wRQtxjzTbXisXibTec+6QQIsNwnbdrHXfqtgshvqt1zkEhREqt1+yi7VZqd38hRLzhnCQhxOBar9lFux2SlFJ91foChgMDgfRax/oAvYA4IKrW8a61z6vnmjOAZ+s43huYDbQAvm/MNR2o7ZcDvwOehp/bu0rbz7v+f4GX7K3tVvo7XwdcZfh+PBBnb+12xC9Vde88Usq/hRBdzzu2G/TV3yzMHdABErD4xRvLSm1/FJglpawwXC/PcNwV2o7h/QK4DRhlOGQ3bbdSuyXQ0vC9P2Csl2s37XZEahik6boJIbYLIf4SQgyr57wnhBCpho+drQCklDvR9zI2AgtrnXvpeR8L67uulsxpewgwTAiRYDhvELhM242GAblSyr3g8G03p91PA7OFEP8A7wDTwOHbrTnVs26aHKCLlDJfCBEJrBJC9JNSnj7vvIXAa+h7FK+h/0j8AICU8sk6rrtfStnfemFbhLltbwa0AmKAQcD3QojuUs/Z2270L+Cb2gcctO3mtvtRYLKUcrkQ4jbgE+AKcNh22wXVs24CKWWFlDLf8H0ysB99T/L883KllNVSSh2wCBh8/jmOxty2A0eAFYbknIj+Y3Bb20VqeY1oO0KIZsBNwHe2i9A6GtHue4EVhu+X4gT/3u2BStZNIIRoJ4RwN3zfHegJZNVxXlCtH28E0s8/x9GY23ZgFYaxWiFECNAcOGGjMK2iEW0HfY9yj5TyiK3is5ZGtDsbGGH4fhSw1zYROjmtn3Da2xf6j6s5QCX6XuG/0SfYI0AFkAusNZx7M7AT2AFsA66tdZ3FGJ6kA18CaUAqsBoIquf+XYEyIKXW11MO3PbmwBL0v6C2AaNcpe2Gnz8HJppxf03abqW/88uAZMN5CUCkvbXbEb9UbRBFURQHoIZBFEVRHIBK1oqiKA5AJWtFURQHoJK1oiiKA1DJWlEUxQGoZK0oiuIAVLJWFEVxACpZK4qiOID/B7zTeNjl61PEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x509.002 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds_boundary = grid.nearby_spectra(ds_spec,dist_thres=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Times in swan controlled by spectra, make sure that they are consistent with the wind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check=[st == wt for st,wt in zip(ds_spec.time.values,ds_wind.time.values)]\n",
    "any(check)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save spectral forcing to swan format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_boundary.spec.to_swan(new_simulation.staging_dir + '/' + new_simulation.default_context['spectra_file'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup spectral output points - a grid and specific points at buoys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_locs = []\n",
    "\n",
    "bbox = grid.bbox()\n",
    "\n",
    "from owslib.wfs import WebFeatureService\n",
    "wfs_url='http://geoserver-123.aodn.org.au/geoserver/wfs'\n",
    "wfs_endpoint = 'aodn:aodn_wave_nrt_timeseries_map'\n",
    "wfs = WebFeatureService(wfs_url,version='2.0.0')\n",
    "response = wfs.getfeature(typename=wfs_endpoint,bbox=tuple(bbox),outputFormat='csv')\n",
    "df = pd.read_csv(response)\n",
    "df = df.groupby(['wmo_id','latitude','longitude']).mean()\n",
    "\n",
    "for (wmo_id,yp,xp), data in df.iterrows():\n",
    "    output_locs.append([str(round(xp,2)),str(round(yp,2))])\n",
    "\n",
    "spec_int = 0.05\n",
    "xx = np.arange(round(bbox[0]+spec_int,2),round(bbox[2]-spec_int,2),spec_int)\n",
    "yy = np.arange(round(bbox[1]+spec_int,2),round(bbox[3]-spec_int,2),spec_int)\n",
    "specout_x, specout_y = np.meshgrid(xx,yy)\n",
    "for i, (xp,yp) in enumerate(zip(specout_x.ravel(),specout_y.ravel())):\n",
    "    output_locs.append([str(round(xp,2)),str(round(yp,2))])\n",
    "\n",
    "new_simulation.settings[\"output_locs\"] = {'coords':output_locs}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate the simulation files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'_template': '../rompy/templates/swan',\n",
       " '_checkout': None,\n",
       " 'model': 'SWAN',\n",
       " 'compute_start': '20201214.000000',\n",
       " 'compute_stop': '20201215.000000',\n",
       " 'out_start': '20201214.000000',\n",
       " 'wind_grid': 'REG 115.48384449028222 -32.809999999999995 0 166 256 0.002 0.002 NONSTATION 20201214.000000 1.0 HR',\n",
       " 'wind_read': \"1 'wind.inp' 3 0 1 0 FREE\",\n",
       " 'output_locs': {'coords': [['115.57', '-32.45'],\n",
       "   ['115.58', '-32.71'],\n",
       "   ['115.63', '-32.71'],\n",
       "   ['115.68', '-32.71'],\n",
       "   ['115.58', '-32.66'],\n",
       "   ['115.63', '-32.66'],\n",
       "   ['115.68', '-32.66'],\n",
       "   ['115.58', '-32.61'],\n",
       "   ['115.63', '-32.61'],\n",
       "   ['115.68', '-32.61'],\n",
       "   ['115.58', '-32.56'],\n",
       "   ['115.63', '-32.56'],\n",
       "   ['115.68', '-32.56'],\n",
       "   ['115.58', '-32.51'],\n",
       "   ['115.63', '-32.51'],\n",
       "   ['115.68', '-32.51'],\n",
       "   ['115.58', '-32.46'],\n",
       "   ['115.63', '-32.46'],\n",
       "   ['115.68', '-32.46'],\n",
       "   ['115.58', '-32.41'],\n",
       "   ['115.63', '-32.41'],\n",
       "   ['115.68', '-32.41']]}}"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_simulation.settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/bra467/repos/rompy/notebooks/simulations/test_swan'"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_simulation.generate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### retain the extracted data in netcdf format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_spec.to_netcdf(new_simulation.staging_dir + '/downloads/wavewatch.nc')\n",
    "ds_wind.to_netcdf(new_simulation.staging_dir + '/downloads/access.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (optionally) zip up the files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# new_simulation.zip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/bra467/repos/rompy/notebooks/simulations/test_swan'"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_simulation.staging_dir "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### If a compiled version of the swan model is available this should run..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " SWAN is preparing computation\n",
      "\n",
      "+SWAN is processing output request    1\n",
      "+SWAN is processing output request    2\n",
      "+SWAN is processing output request    3\n",
      " Number of threads during execution of parallel region = 12\n",
      "\n",
      "+time 20201214.001500   , step      1; iteration    1; sweep 1\n",
      "+time 20201214.001500   , step      1; iteration    1; sweep 2\n",
      "+time 20201214.001500   , step      1; iteration    1; sweep 3\n",
      "+time 20201214.001500   , step      1; iteration    1; sweep 4\n",
      "+time 20201214.001500   , step      1; iteration    2; sweep 1\n",
      "+time 20201214.001500   , step      1; iteration    2; sweep 2\n",
      "+time 20201214.001500   , step      1; iteration    2; sweep 3\n",
      "+time 20201214.001500   , step      1; iteration    2; sweep 4\n",
      "+time 20201214.003000   , step      2; iteration    1; sweep 1\n",
      "+time 20201214.003000   , step      2; iteration    1; sweep 2\n",
      "+time 20201214.003000   , step      2; iteration    1; sweep 3\n",
      "+time 20201214.003000   , step      2; iteration    1; sweep 4\n",
      "+time 20201214.004500   , step      3; iteration    1; sweep 1\n",
      "+time 20201214.004500   , step      3; iteration    1; sweep 2\n",
      "+time 20201214.004500   , step      3; iteration    1; sweep 3\n",
      "+time 20201214.004500   , step      3; iteration    1; sweep 4\n",
      "+time 20201214.010000   , step      4; iteration    1; sweep 1\n",
      "+time 20201214.010000   , step      4; iteration    1; sweep 2\n",
      "+time 20201214.010000   , step      4; iteration    1; sweep 3\n",
      "+time 20201214.010000   , step      4; iteration    1; sweep 4\n",
      "+SWAN is processing output request    1\n",
      "+SWAN is processing output request    2\n",
      "+SWAN is processing output request    3\n",
      "+time 20201214.011500   , step      5; iteration    1; sweep 1\n",
      "+time 20201214.011500   , step      5; iteration    1; sweep 2\n",
      "+time 20201214.011500   , step      5; iteration    1; sweep 3\n",
      "+time 20201214.011500   , step      5; iteration    1; sweep 4\n",
      "+time 20201214.013000   , step      6; iteration    1; sweep 1\n",
      "+time 20201214.013000   , step      6; iteration    1; sweep 2\n",
      "+time 20201214.013000   , step      6; iteration    1; sweep 3\n",
      "+time 20201214.013000   , step      6; iteration    1; sweep 4\n",
      "+time 20201214.014500   , step      7; iteration    1; sweep 1\n",
      "+time 20201214.014500   , step      7; iteration    1; sweep 2\n",
      "+time 20201214.014500   , step      7; iteration    1; sweep 3\n",
      "+time 20201214.014500   , step      7; iteration    1; sweep 4\n",
      "+time 20201214.020000   , step      8; iteration    1; sweep 1\n",
      "+time 20201214.020000   , step      8; iteration    1; sweep 2\n",
      "^C\n"
     ]
    }
   ],
   "source": [
    "!cd simulations/test_swan && ~/swan4131/swan.exe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}