{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba358251-dcb1-4696-ae7b-65156b3f5797",
   "metadata": {},
   "source": [
    "# Surface and bottom boundary layer dynamics\n",
    "\n",
    "In particle tracking simulations, it may be desirable to apply corrections near to model boundaries to account for physical processes that aren't resolved in the gridded input data that is used to drive the model. Examples include near to the surface and bottom boundaries, where the gridded field often ends a half-cell width away from the physical boundary. If a particle sits between one of these points and the physical boundary, one can extrapolate the field to the particle's position. However, this may be inconsistent with the model's boundary conditions. For example, in the bottom boundary layer, empicial studies have shown the velocity field follows a logarithmic profile near to the boundary, which can be modelled using the expression:\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\bf{u}(z) = \\frac{\\bf{u}_{_{\\text{*}}}}{\\kappa} ln\\left(\\frac{z}{z_{\\text{0}}}\\right)\n",
    "\\end{equation}\n",
    "$$\n",
    "where $\\bf{u}_{\\text{z}}$ is the velocity at height $z$ above the sea floor, $\\bf{u}_{\\text{*}}$ is the friction velocity, $\\kappa$ is the von Karman constant, and $z_{\\text{0}}$ is the roughness height. By extrapolating the velocity field down to the bed, one will tend to overestimate the near bed velocity field.\n",
    "\n",
    "In PyLag, particle objects have two properties which are used to signify whether the particle resides in the surface or bottom boundary layers. These are:\n",
    "\n",
    "* `in_surface_boundary_layer` - If this attribute is set to `True`, the particle is in the surface boundary layer.\n",
    "* `in_bottom_boundary_layer` - If this attribute is set to `True`, the particle is in the bottom boundary layer.\n",
    "\n",
    "These two properties are continuously updated as particles move in and out of the two layers, and they can be used to signify whether a particle needs special treatment near to the boundary.\n",
    "\n",
    "PyLag support for handling boundary layer dynamics is still under development. Presently, there is support for applying a logarithmic correction to the near bed velocity field in particle tracking simulations that use FVCOM driving data. This option is turned on and off using the run configuration file and the following parameters, which should be set in the configuration file section titled `OCEAN_DATA`:\n",
    "\n",
    "* `use_near_bottom_log_velocity_profile` - If this attribute is set to `True`, the Law of the Wall is applied to the velocity field near the bottom boundary.\n",
    "* `z0` - The roughness length in meters.\n",
    "\n",
    "In the run configuration file, this appears as:\n",
    "\n",
    "```python\n",
    "[OCEAN_DATA]\n",
    "\n",
    "# Apply the Law of the Wall to the velocity field near the bottom boundary\n",
    "use_near_bottom_log_velocity_profile = False\n",
    "\n",
    "# The roughness height in meters (will vary depending on the bottom type, and from simulation to simulation)\n",
    "z0 = 4.e-4 # Example value that approximates unrippled sand \n",
    "```\n",
    "\n",
    "For now, a bed roughness height that is constant across the full domain must be used. In the future, we may include the option to specify a spatially varying roughness length scale.\n",
    "\n",
    "## Example: Correcting the near bottom velocity field using FVCOM data for the Tamar Esturary\n",
    "\n",
    "This capability is demonstrated using the [FVCOM Tamara Estuary setup](../../examples/fvcom_forward_tracking.ipynb), as described in the [examples](../../examples/index.rst) section of PyLag's documentation. Here, we don't perform a full simulation using PyLag. Rather, we create an object of type `FVCOMDataReader`, which is configured to read in FVCOM data, and use it to compute the velocity profile at a point in space corresponding to the L4 time series site in the Western English Channel. The result of applying the logarithmic profile correction is contrasted with the result of extrapolating the field down to the sea floor.\n",
    "\n",
    "### Creating a grid metrics file\n",
    "\n",
    "In this example, we will use the grid metrics file that is supplied with PyLag's example [datasets](https://drive.google.com/open?id=15UX7Y9JnuLpnPAz700mzmzd917nTClxR). As in the tutorial examples, we assume this file has been saved in a directory called `data/pylag_doc`, which sits in the home directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f82e210b-6ab7-40fe-ae0f-90d243445cf1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating FVCOM grid metrics file /users/modellers/jcl/code/git/PyLag/PyLag/doc/source/documentation/cookbook/bdy_dynamics/grid_metrics.nc\n",
      "\n",
      "Calculating element areas ... done\n",
      "Grid has 45 nodes on the open boundary\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import pathlib\n",
    "\n",
    "from pylag.grid_metrics import create_fvcom_grid_metrics_file\n",
    "\n",
    "\n",
    "# Root directory for PyLag example input files\n",
    "home_dir = os.environ['HOME']\n",
    "data_dir=f'{home_dir}/data/pylag_doc'\n",
    "\n",
    "# An example FVCOM output file\n",
    "fvcom_file_name = f'{data_dir}/fvcom_tamar_estuary_0001.nc'\n",
    "\n",
    "# The file listing the location of open boundary nodes\n",
    "# NB - if this is omitted, the model treats open boundaries\n",
    "# as land boundaries.\n",
    "obc_file_name = f'{data_dir}/fvcom_tamar_estuary_obc.dat'\n",
    "\n",
    "# Create example directory\n",
    "example_dir = f\"{pathlib.Path.cwd()}/bdy_dynamics\"\n",
    "pathlib.Path(example_dir).mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "# The name of the output file\n",
    "grid_metrics_file_name = f'{example_dir}/grid_metrics.nc'\n",
    "\n",
    "# Generate the file\n",
    "create_fvcom_grid_metrics_file(fvcom_file_name, obc_file_name=obc_file_name,\n",
    "                               grid_metrics_file_name=grid_metrics_file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddbfca05-0308-4fb6-86f9-6b4c9ae395a9",
   "metadata": {},
   "source": [
    "We now create a run configuration object which specifies that we are working with the FVCOM tamar estuary data, and that we wish to correct the velocity near to the bed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bdc15524-38da-43ed-80e9-8015bb44ead7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import configparser\n",
    "import datetime\n",
    "\n",
    "from pylag.fvcom_data_reader import FVCOMDataReader\n",
    "from pylag.mediator import SerialMediator\n",
    "\n",
    "\n",
    "# Create config\n",
    "config = configparser.ConfigParser()\n",
    "config.add_section('SIMULATION')\n",
    "config.set('SIMULATION', 'coordinate_system', 'geographic')\n",
    "config.add_section('OCEAN_DATA')\n",
    "config.set('OCEAN_DATA', 'name', 'FVCOM')\n",
    "config.set('OCEAN_DATA', 'data_dir', f'{data_dir}')\n",
    "config.set('OCEAN_DATA', 'grid_metrics_file', f'{grid_metrics_file_name}')\n",
    "config.set('OCEAN_DATA', 'data_file_stem', 'fvcom_tamar_estuary_0')\n",
    "config.set('OCEAN_DATA', 'rounding_interval', '3600')\n",
    "config.set('OCEAN_DATA', 'Kz_method', 'none')\n",
    "config.set('OCEAN_DATA', 'Ah_method', 'none')\n",
    "config.set('OCEAN_DATA', 'has_is_wet', 'True')\n",
    "\n",
    "# Extra options to help the config pass PyLag's internal checks\n",
    "config.add_section('NUMERICS')\n",
    "config.set('NUMERICS', 'num_method', 'test')\n",
    "config.set('NUMERICS', 'time_step_adv', '1')\n",
    "\n",
    "# Set date limits. These carry no special meaning here. PyLag uses these\n",
    "# limits to ensure the supplied input data covers the period we are interested\n",
    "# in. We must simply specify a time between these limits when creating the\n",
    "# vertical profile.\n",
    "datetime_start = datetime.datetime(2020,5,2,3)\n",
    "datetime_end = datetime.datetime(2020,5,2,18)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9cca837-046c-4130-8c6b-22f64a53547d",
   "metadata": {},
   "source": [
    "To reconstruct a velocity profile, we use PyLag's Particle class to create a particle object. We then vary the position of this particle, and use PyLag's FVCOMDataReader to compute the velocity at the particle's position. We use PyLag's Python API to accomplish this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0ed91588-6b89-4e49-8c21-d94d6659edec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from pylag.particle_cpp_wrapper import ParticleSmartPtr\n",
    "from pylag.data_types_python import DTYPE_FLOAT\n",
    "\n",
    "\n",
    "def compute_current_profile(apply_log_correction):\n",
    "    \"\"\" Generate the velocity profile\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    apply_log_correction: bool\n",
    "     Boolean indicating whether the log correction should be applied or not. True or False\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    depth_arr: NetCDF array\n",
    "        Depth array\n",
    "        \n",
    "    speed_arr: NetCDF array\n",
    "        Current speed array\n",
    "    \"\"\"\n",
    "    if apply_log_correction:\n",
    "        config.set('OCEAN_DATA', 'apply_near_bottom_log_velocity_profile', 'True')\n",
    "        config.set('OCEAN_DATA', 'z0', '1.e-4')\n",
    "    else:\n",
    "        config.set('OCEAN_DATA', 'apply_near_bottom_log_velocity_profile', 'False')\n",
    "\n",
    "    # Create a mediator object, which controls file access\n",
    "    mediator = SerialMediator(config, 'ocean', datetime_start, datetime_end)\n",
    "    \n",
    "    # Create the FVCOM data reader\n",
    "    data_reader = FVCOMDataReader(config, mediator)\n",
    "    \n",
    "    # Set location to Station L4. Internally, PyLag stored particle positions\n",
    "    # in radians, so we convert these here. NB - this is different to how positions\n",
    "    # are speficied in the initial positions file.\n",
    "    x1 = -4.17 * np.radians(1)\n",
    "    x2 = 50.25 * np.radians(1)\n",
    "    \n",
    "    # Create a particle object\n",
    "    particle = ParticleSmartPtr(x1=x1, x2=x2)\n",
    "    \n",
    "    # Determine the host hoizontal element - needed for interpolating the velocity field\n",
    "    flag = data_reader.find_host_using_global_search_wrapper(particle)\n",
    "    \n",
    "    if flag != 0:\n",
    "        raise RuntimeError('Failed to locate the particle')\n",
    "    \n",
    "    # Define the time at which we want to compute the velocity profile\n",
    "    datetime_now = datetime.datetime(2020,5,2,12)\n",
    "    time = (datetime_now - datetime_start).total_seconds()\n",
    "    \n",
    "    # Read data for the current time\n",
    "    data_reader.read_data(time)\n",
    "    \n",
    "    # Compute the water depth and the sea surface height at the particle's spatiotemporal coordinates\n",
    "    h = data_reader.get_zmin_wrapper(time, particle)\n",
    "    zeta = data_reader.get_zmax_wrapper(time, particle)\n",
    "    \n",
    "    # Define a vertical grid of points between h and zeta\n",
    "    n_depths = 1000\n",
    "    depth_arr = np.linspace(h, zeta, n_depths)\n",
    "    \n",
    "    # Create a profile of the current speed\n",
    "    current_speed_arr = np.empty_like(depth_arr)\n",
    "    for i in range(n_depths):\n",
    "    \n",
    "        # Update the depth of the particle\n",
    "        particle.x3 = depth_arr[i]\n",
    "    \n",
    "        # Set vertical interpolation coefficients\n",
    "        if data_reader.set_vertical_grid_vars_wrapper(time, particle) != 0:\n",
    "            raise RuntimeError('Failed to set vertical grid variables')\n",
    "\n",
    "        # Compute the velocity at the particle's spatiotemporal coordinates\n",
    "        vel = np.empty(3, dtype=DTYPE_FLOAT)\n",
    "        data_reader.get_velocity_wrapper(time, particle, vel)\n",
    "    \n",
    "        # Compute speed from the horizontal component only\n",
    "        current_speed_arr[i] = np.sqrt(vel[0]**2 + vel[1]**2)\n",
    "\n",
    "    return depth_arr, current_speed_arr, h, zeta\n",
    "\n",
    "# No correction\n",
    "z_not_corrected, curr_not_corrected, h, zeta = compute_current_profile(False)\n",
    "\n",
    "# With correction\n",
    "z_corrected, curr_corrected, _, _ = compute_current_profile(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c241f14-6ac4-4549-939d-a5a767a524e0",
   "metadata": {},
   "source": [
    "We can now plot the current profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "188f4d48-b8b7-4394-ab14-d5572156bad4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAFuCAYAAAB5gpTAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPiUlEQVR4nO3deVhUdf/G8fewDZsgCgIqgruioiguiIk+KaZmmpmVZVLqU7lrZdrmUmalpumTVlbaqmYumVZiliuuiJlC4ILiAqKCIKhs8/39wc+pEVRQhsPA53Vdc+mc850ztyeam7PMOTqllEIIIYQoASutAwghhLA8Uh5CCCFKTMpDCCFEiUl5CCGEKDEpDyGEECUm5SGEEKLEpDyEEEKUmJSHEEKIErPROkB5YzAYOHfuHFWqVEGn02kdRwghSkQpxZUrV6hZsyZWVubbPpDyuMm5c+fw8fHROoYQQtyT06dPU7t2bbMtv0KWx8KFC5k1axZJSUk0a9aMefPmcd999xXrtVWqVAEKVryLi4s5YwohKpC9r66l3UdD+KtKMC3O/KpZjoyMDHx8fIyfZeZS4cpjxYoVjBs3joULFxISEsInn3xCz549iYmJoU6dOnd8/Y1dVS4uLlIeQohic9I74gI462zKxWeHuXe7V7gD5h988AFDhw5l2LBhNG3alHnz5uHj48OiRYu0jiaEEBVGhSqPnJwcoqKiCAsLM5keFhZGZGSkRqmEEKLiqVC7rS5evEh+fj6enp4m0z09PUlOTi7yNdnZ2WRnZxufZ2RkmDWjEKLiyL2ay+GPd5D+3XraRH0KQL61rcapykaFKo8bbt7Xp5S65f6/mTNnMm3atBItPz8/n9zc3LvOJyona2trbGxs5BRwC3fhSApx837B+tcN+J/ZSCD//MKZjR3XB/9Xw3Rlp0KVh7u7O9bW1oW2MlJSUgptjdwwefJkJkyYYHx+40yFW8nMzOTMmTPIPbTE3XB0dMTb2xs7Ozuto4hiUgZF3IqDJH++Affd6/HP2ksn/vn//4LOg7h6vbB+qDf+48LoWMdVw7Rlp0KVh52dHW3atGHTpk08/PDDxumbNm2ib9++Rb5Gr9ej1+uLtfz8/HzOnDmDo6MjHh4e8hukKDalFDk5OVy4cIGEhAQaNmxo1i9wiXuTlZLF4Q83k716PQ3jf6aJ4SxN/jU/1iGQ82164x7+IE0HB9HJzlqzrFqpUOUBMGHCBAYPHkxQUBDBwcF8+umnJCYm8vzzz9/zsnNzc1FK4eHhgYODQymkFZWJg4MDtra2nDp1ipycHOzt7bWOJP7lzI6THP9wPY5bNtDi4h+0559joVk4ctirG9ndH6ThmF40DapFUw2zlgcVrjwee+wxLl26xPTp00lKSqJ58+b8/PPP+Pr6ltp7yBaHuFuytVF+5F3P48jiSNK+2UDtP9fTIDuGf38f+7SNHyeaPojTwN40H9WF9lWl7P+twpUHwIgRIxgxYoTWMYQQ5Uzq0UvEzv0Vft6Af+KvtFRpxnl5WHPYNYTLHXtT+/kHqf9gU3ys5BfFW6mQ5SEqvqlTp7J27VoOHjyoyfvrdDrWrFlDv379NHl/UTzKoDi65jDnPl1PtV0baHZlFyEYjPNTddWI9e2Jrs+D+I/vQau6bhqmtSyyDV0JhIeHo9PpCj0eeOCBYr2+S5cujBs3zrwhy4BOp2Pt2rVaxxBmdi31GvumbmBb8xGcs/Ol0YAAukS8SsCVnVhjIN6+BVuCJ3No4Q5cr6cQkvANHec/TlUpjhKRLY9K4oEHHmDJkiUm04p7lllxKKXIz8/HxkZ+pETZO7fnNMc+3ID95vU0T/mdtlwzzruGPX/VuJ/r/+lN/TG9aRRch0YaZq0oZMujktDr9Xh5eZk83Nzc2LJlC3Z2dmzfvt04ds6cObi7u5OUlER4eDhbt27lww8/NG6xnDx5ki1btqDT6di4cSNBQUHo9Xq2b9/O8ePH6du3L56enjg7O9O2bVt+++03kyx+fn689dZbDBo0CGdnZ2rWrMmCBQtMxiQmJtK3b1+cnZ1xcXFh4MCBnD9//pb/vn379tG9e3fc3d1xdXUlNDSUAwcOmLwnwMMPP4xOpzM+B/jpp59o06YN9vb21KtXj2nTppGXl2ecf/ToUTp37oy9vT3+/v5s2rTpbv4TiFKUn5PPX59EsqXjq8Q5tKRmhzp0XvYC7VI24Mg1zlr7sM3/efZNWY+6cIl259fTedkL1Aq+88VRRTEpYSI9PV0BKj09vdC8a9euqZiYGHXt2jWllFIGg1KZmdo8DIbi/5uGDBmi+vbte8v5L7/8svL19VWXL19WBw8eVHq9Xq1evVoppdTly5dVcHCwGj58uEpKSlJJSUkqLy9P/fHHHwpQAQEBKiIiQh07dkxdvHhRHTx4UH388cfq0KFDKj4+Xr322mvK3t5enTp1yvh+vr6+qkqVKmrmzJkqLi5OzZ8/X1lbW6uIiIj/X68GFRgYqDp16qT279+vdu/erVq3bq1CQ0ONy5gyZYpq2bKl8fnmzZvV119/rWJiYlRMTIwaOnSo8vT0VBkZGUoppVJSUhSglixZopKSklRKSopSSqlff/1Vubi4qKVLl6rjx4+riIgI5efnp6ZOnaqUUio/P181b95cdenSRUVHR6utW7eqwMBABag1a9YU/z/C/7v5Z0gUX9qJVLVz9DK1ve5T6qKuulJgfORhpf6sEqL+CHtHxf1wSBnyS/A/SAVzu8+w0iTlcZOSlEdmpsnPb5k+MjOL/28aMmSIsra2Vk5OTiaP6dOnK6WUys7OVoGBgWrgwIGqWbNmatiwYSavDw0NVWPHjjWZdqM81q5de8f39/f3VwsWLDA+9/X1VQ888IDJmMcee0z17NlTKaVURESEsra2VomJicb5R44cUYDau3evUqpwedwsLy9PValSRf3000/GaUV94N93333qnXfeMZn29ddfK29vb6WUUhs3blTW1tbq9OnTxvm//PKLlEcZyb6SrbY984WKdu2scrE2+Z8gTVdV7azzuNrxwjfqUvxFraOWG2VVHrKDupLo2rVrocvSV6tWDSj4Zv4333xDQEAAvr6+zJs3r9jLDQoKMnmelZXFtGnTWL9+PefOnSMvL49r166RmJhoMi44OLjQ8xvvGxsbi4+Pj8llYvz9/alatSqxsbG0bdu2UI6UlBTefPNNfv/9d86fP09+fj5Xr14t9L43i4qKYt++fcyYMcM4LT8/n+vXr3P16lViY2OpU6eOyR3Zbs4uSl9OZg67nl1Mo1Uzuc9w1jj9mN6fMy0fxO2p3jQb3pGO9vIRphVZ8/fA0REyM7V775JwcnKiQYMGt5x/45L1qamppKam4uTkVOzl/tvLL7/Mxo0bmT17Ng0aNMDBwYEBAwaQk5Nzx2Xd+PKlusWFLG81HQrOKLtw4QLz5s3D19cXvV5PcHDwHd/XYDAwbdo0+vfvX2ievb19kdcwky+Jmo8hz8Du8Suo9fHrhOadACDZypu/u4+h/quP0aBzXW79UyzKkpTHPdDpoJifseXa8ePHGT9+PIsXL+b777/n6aefZvPmzcZvQ9vZ2ZGfn1+sZW3fvp3w8HDjtcUyMzM5efJkoXG7d+8u9LxJk4KrB/n7+5OYmMjp06eNWx8xMTGkp6fTtGnRF4XYvn07CxcupFevXkDBbYQvXrxoMsbW1rbQv6N169bExcXdslhvZDl37hw1a9YEYNeuXbdbBeIuKIMiamYETjMm0/FaNAAXdDWIGfAmHT4bRheX0jszUJQOKY9KIjs7u9DVhm1sbHBzc2Pw4MGEhYXxzDPP0LNnT1q0aMGcOXN4+eWXgYIzlfbs2cPJkydxdnY27u4qSoMGDVi9ejV9+vRBp9PxxhtvYDAYCo3buXMn77//Pv369WPTpk2sXLmSDRs2ANCtWzcCAgJ48sknmTdvHnl5eYwYMYLQ0NBCu8n+/b5ff/01QUFBZGRk8PLLLxe6/pifnx+bN28mJCQEvV6Pm5sbb775Jg8++CA+Pj48+uijWFlZcejQIf766y/efvttunXrRuPGjXn66aeZM2cOGRkZvPbaayVa9+L2jizZS86ESQRd/gOADKpw4P6JBH0zjlAvZ43TiVsy6xEVC1SSA+aWYsiQIQoo9GjcuLGaNm2a8vb2Vhcv/nPAce3atcrOzk5FR0crpZSKi4tTHTp0UA4ODgpQCQkJxgPmaWlpJu+VkJCgunbtqhwcHJSPj4/63//+V+iAu6+vr5o2bZoaOHCgcnR0VJ6enmrevHkmyzl16pR66KGHlJOTk6pSpYp69NFHVXJysnH+zQfMDxw4oIKCgpRer1cNGzZUK1euVL6+vmru3LnGMevWrVMNGjRQNjY2ytfX1zj9119/VR07dlQODg7KxcVFtWvXTn366afG+XFxcapTp07Kzs5ONWrUSP36669ywLwUJETEq8hajxgPgF/HTv3Rery6+PcFraNZtLI6YK5TSm5M8W8ZGRm4urqSnp5e6Cb2169fJyEhgbp168oVUe+Bn58f48aNqxDfWi8p+RmC9FOXiX7kbTpGzceOXAzoiKw3GL+vplM7pPQuYFpZ3e4zrDTJbishRJnIu55H5LOf4b/8DbqoguNR+6v3wOWTWXR6pIXG6URJSXkIIczuwKzNVHljHJ2zDwNw3K4Jqa99QNs3e2qcTNwtKQ9R5oo6+0pUTKd+P07ykxNon7wOKLiK7V+PTKPjl89R39FW43TiXsi1rYQQpS47I5st/5mO5/3NaJ+8jlxs2NpqLLqjRwldOQpbKQ6LJ1seQohSFfXuJqq9OZIuuUcLnlfrhtvXCwjt1eQOrxSWRMpDCFEqkg+c40S/CXQ8vaLguZU3J0bNJXjuQHRyR74KR8pDCHHP4n84hNejnejIFfKxYker0QT+NJ2Otc13qqjQlpSHEOKeJb39OY24Qpx9AOqLpYQ+Eah1JGFmcsBcCHHPqh3fC8DFZ1+hiRRHpSDlIUQpunGHxcuXL2sdpcwY8gzUyzwEgHfv1hqnEWVFyqMSCA8Pp1+/flrHKLcq4wd+aTobeQonrpKNHXX+IxdMryykPIRFyM/PL/LqvMW5T4gwr6TfjgBwyr4xNnJzpkpDykOwdetW2rVrh16vx9vbm0mTJpGXl2ecf+XKFZ588kmcnJzw9vZm7ty5dOnS5Y4XNly3bh1BQUHY29vj7u5ucsOltLQ0nn76adzc3HB0dKRnz54cPXrUOH/p0qVUrVqV9evX4+/vj16v59SpU/j5+fH2228THh6Oq6srw4cPBwpuZtW5c2ccHBzw8fFhzJgxZGVlGZeXnZ3NxIkT8fHxQa/X07BhQz7//HNOnjxJ165dAXBzc0On0xEeHg4U3Hzq/fffp169ejg4ONCyZUt++OEHk3/jzz//TKNGjXBwcKBr166V8tvzV/cVlMfFGs00TiLKkpTHvVAKsrK0eZTSxZDPnj1Lr169aNu2LX/++SeLFi3i888/5+233zaOmTBhAjt37mTdunVs2rSJ7du3c+DAgdsud8OGDfTv35/evXsTHR3N5s2bTe7FER4ezv79+1m3bh27du1CKUWvXr3Izc01jrl69SozZ87ks88+48iRI9SoUQOAWbNm0bx5c6KionjjjTf466+/6NGjB/379+fQoUOsWLGCHTt2MGrUKOOynn76aZYvX878+fOJjY3l448/xtnZGR8fH1atWgVAXFwcSUlJfPjhhwC8/vrrLFmyhEWLFnHkyBHGjx/PU089xdatW4GCG07179+fXr16cfDgQYYNG8akSZPu8b+I5bGJjwEgp6GUR6Vi1gu+W6AS3c8jM9N4L4Iyf2RmFvvfNGTIENW3b98i57366quqcePGymAwGKd99NFHytnZWeXn56uMjAxla2urVq5caZx/+fJl5ejoaHKPjpsFBwerJ598ssh58fHxClA7d+40Trt48aJycHBQ33//vVJKqSVLlihAHTx40OS1vr6+ql+/fibTBg8erP773/+aTNu+fbuysrJS165dU3FxcQpQmzZtKjJPUfcmyczMVPb29ioyMtJk7NChQ9UTTzyhlFJq8uTJqmnTpibr7pVXXinyPic3VMT7eRx06aQUqJ1jlmsdRaiyu5+H7KCs5GJjYwkODja5L3dISAiZmZmcOXOGtLQ0cnNzadeunXG+q6srjRs3vu1yDx48aNylVNR72tjY0L59e+O06tWr07hxY2JjY43T7OzsCAgIKPT6m+8mGBUVxbFjx/j222+N05RSGAwGEhIS+Ouvv7C2tiY0NPS2mf8tJiaG69ev0717d5PpOTk5BAYGGv8dHTp0MFl3wcHBxX6PisIj6yQALgF+muYQZUvK4144OkJmpnbvXQqUUiYffjemAeh0OpO/FzXmVm6+BWxxXntzFgcHh0LvC+B0043jDQYDzz33HGPGjCk0tk6dOhw7duy2WYty4+D8hg0bqFWrlsk8vV5vzFvZ5WTm4Jl/DoAabeVGTpWJlMe90Ongpg8yS+Pv78+qVatMPrgjIyOpUqUKtWrVomrVqtja2rJ37158fHyAgjuVHT169La/yQcEBLB582aeeeaZIt8zLy+PPXv20LFjRwAuXbpEfHw8TZs2LfG/oXXr1hw5coQGDYo+TbRFixYYDAa2bt1Kt27dCs23s7MDCs7o+ndGvV5PYmLiLf+d/v7+rF271mTa7t27S5zfkiXvP0MdDFzDHo/mnlrHEWVIDphXEunp6Rw8eNDkkZiYyIgRIzh9+jSjR4/m77//5scff2TKlClMmDABKysrqlSpwpAhQ3j55Zf5448/OHLkCM8++yxWVlZFbhXcMGXKFJYtW8aUKVOIjY3lr7/+4v333wegYcOG9O3bl+HDh7Njxw7+/PNPnnrqKWrVqkXfvn1L/G975ZVX2LVrFyNHjuTgwYMcPXqUdevWMXr0aKDgtrdDhgzh2WefZe3atSQkJLBlyxa+//57AHx9fdHpdKxfv54LFy6QmZlJlSpVeOmllxg/fjxffvklx48fJzo6mo8++ogvv/wSgOeff57jx48zYcIE4uLi+O6771i6dGmJ81uy1D9PA5Bs6yMXP6xszHpExQKV6IC5hRgyZIgCCj2GDBmilFJqy5Ytqm3btsrOzk55eXmpV155ReXm5hpfn5GRoQYNGqQcHR2Vl5eX+uCDD1S7du3UpEmTbvu+q1atUq1atVJ2dnbK3d1d9e/f3zgvNTVVDR48WLm6uioHBwfVo0cPFR8fb5y/ZMkS5erqWmiZvr6+au7cuYWm7927V3Xv3l05OzsrJycnFRAQoGbMmGGcf+3aNTV+/Hjl7e2t7OzsVIMGDdQXX3xhnD99+nTl5eWldDqdcb0YDAb14YcfqsaNGytbW1vl4eGhevToobZu3Wp83U8//aQaNGig9Hq9uu+++9QXX3xRqQ6Y73l9nVKg/nJqp3UU8f/K6oC5TinZcftvt7t5/PXr10lISKBu3brY29trlFB7WVlZ1KpVizlz5jB06FCt41iUivYzFDlmOR0XPEF01a4Epv2udRzB7T/DSpMc8xB3FB0dzd9//027du1IT09n+vTpAHe1i0lULPmZ1wDIs731CRKiYpLyEMUye/Zs4uLisLOzo02bNmzfvh13d3etYwmt/f9ZaQYra42DiLIm5SHuKDAwkKioKK1jiHJMJ3u/Kx0520oIcfeMZ9xJeVQ2Uh5CiLt3m9O1RcUm5XEX5AQ1cbcq6s+O7LaqfKQ8SsDauuCgoNxDQtytq1evAmBra6txktLxzxcDpTwqGzlgXgI2NjY4Ojpy4cIFbG1tsbKS7hXFo5Ti6tWrpKSkULVqVeMvIhZPdltVWlIeJaDT6fD29iYhIYFTp05pHUdYoKpVq+Ll5aV1jNInu60qHSmPErKzs6Nhw4ay60qUmK2tbcXZ4rjh/7c85JhH5SPlcResrKwqxKUlhLhncjHESkt22gsh7pqNS8EtCapkJWmcRJQ1KQ8hxF2rP6QTAI2vH+LC4fMapxFlScpDCHHXPJrVINah4La88Qt/0ziNKEtSHkKIe3I+IAwAtTFC4ySiLFlMecyYMYOOHTvi6OhI1apVixyTmJhInz59cHJywt3dnTFjxshZUUKYmeujBeXR8GQEyiBnXVUWFlMeOTk5PProo7zwwgtFzs/Pz6d3795kZWWxY8cOli9fzqpVq3jxxRfLOKkQlYv/8BCu4oCnIZmjaw5rHUeUEYspj2nTpjF+/HhatGhR5PyIiAhiYmL45ptvCAwMpFu3bsyZM4fFixeTkZFRxmmFqDz0LnqOeHQB4NxS2XVVWVhMedzJrl27aN68OTVr1jRO69GjB9nZ2be9F0V2djYZGRkmDyFEyVwNKdh15bxLyqOyqDDlkZycjKenp8k0Nzc37OzsSE5OvuXrZs6ciaurq/Hh4+Nj7qhCVDi1nikoj2aXtnEt9ZrGaURZ0LQ8pk6dik6nu+1j//79xV6eroiLtCmlipx+w+TJk0lPTzc+Tp8+fVf/FiEqs/oPNiXJqhYOXCfm0x1axxFlQNPLk4waNYrHH3/8tmP8/PyKtSwvLy/27NljMi0tLY3c3NxCWyT/ptfr0ev1xXoPIUTRdFY6jtUPw/voEq6sioBJ3bWOJMxM0/Jwd3fH3d29VJYVHBzMjBkzSEpKwtvbGyg4iK7X62nTpk2pvIcQ4tasHwiDo0vwPhwBzNI6jjAziznmkZiYyMGDB0lMTCQ/P5+DBw9y8OBBMjMzAQgLC8Pf35/BgwcTHR3N5s2beemllxg+fDguLi4apxei4ms8shsGdDS+fogTv8RpHUeYmcWUx5tvvklgYCBTpkwhMzOTwMBAAgMDjcdErK2t2bBhA/b29oSEhDBw4ED69evH7NmzNU4uROVQvbE7+2v0BuD88Nc0TiPMTacq6k2V71JGRgaurq6kp6fLFosQJXR0zWHq9W+JNQb++iSSFv8N1jpSpVNWn2EWs+UhhCj/Gj7cnMhGzwCgXnpZLldSgUl5CCFKVYNvp3EVBwKu7GTvaz9qHUeYiZSHEKJUeQfVYm/IBABqzHmF3Ku5GicS5iDlIYQoda2XT+Sizp26ufHsGva51nGEGUh5CCFKnUttF448MgWApsuncOXcFY0TidIm5SGEMIvgJf/lpG0DPFQKUYPmaB1HlDIpDyGEWdg525E0eiYAQVtnc/5gksaJRGmS8hBCmE2HWY9w2Kk9zmQR98RUreOIUiTlIYQwG52Vjvx3C65z1fHvzzn1+3GNE4nSIuUhhDCrlqPuY59HT2zI5/R/39I6jiglUh5CCLNzfG8qAMHHvyZhY7y2YUSpkPIQQphds2fasbfGg1hj4Nzz07WOI0qBlIcQokxU+WAaAB1OLuP4+liN04h7JeUhhCgTTZ9szW7vflhj4PxI2fqwdFIeQogy4zZvKgAdEldwdM1hbcOIeyLlIYQoM40HtmRXrUewQnFx9DSt44h7IOUhhChT7gumYkBH8NkfiPv+T63jiLsk5SGEKFMNH27Obp+BAKSNm6ptGHHXpDyEEGXOc+EUDOjokLSW2G8PaB1H3AUpDyFEmav/YFN2+Q0CIOPFqdqGEXdFykMIoYmaH79JPla0P/8TMV/u0zqOKCEpDyGEJur2aMSuek8BkPXyFI3TiJKS8hBCaMbn0zfIw5q2F36RYx8WRspDCKEZ3/sbsN+7DwDnv9qocRpRElIeQghNXW/XGQDHg5EaJxElIeUhhNCUR98QABpciMSQZ9A4jSguKQ8hhKYaPRbIVRyoplJJ+DVO6ziimKQ8hBCasnW0Jd61LQDnfpBdV5ZCykMIobnL/gW7rojcqW0QUWxSHkIIzTmFFZRHrVOy5WEppDyEEJprMDgYgHo5cVyKu6hxGlEcUh5CCM251a/GcbumABz7SrY+LIGUhxCiXDhbt2DX1bXNUh6WQMpDCFEu6Dp2BMAtVg6aWwIpDyFEuVD7sYItj8YZ+8jJzNE4jbgTKQ8hRLng170hqbpq2JPNiZ+OaB1H3IGUhxCi3LBVBVscdlUdNU4i7kTKQwhRLmScTqcKmQB4tfXROI24EykPIUS5cH5fIgCXdNVxdJctj/JOykMIUS6kHz4NQIp9HY2TiOKQ8hBClAvX4gq2PDJcZZeVJZDyEEKUC4ZTBVse1z1ky8MSSHkIIcoF26SCLQ9VW7Y8LIGUhxCiXHBKK9jysKsv5WEJpDyEEOVC9ayCLQ9nf9ltZQmkPIQQmsvPyccz7ywA1VvJloclkPIQQmjuwuHz2JFLPlZ4BtbUOo4ohhKVR1xcHFOnTuX++++nfv36eHt7ExAQwJAhQ/juu+/Izs42S8iTJ08ydOhQ6tati4ODA/Xr12fKlCnk5JhePC0xMZE+ffrg5OSEu7s7Y8aMKTRGCFH+XDxQsMvqvHVNbOxtNE4jiqNY/5Wio6OZOHEi27dvp2PHjrRr145+/frh4OBAamoqhw8f5rXXXmP06NFMnDiRcePGodfrSy3k33//jcFg4JNPPqFBgwYcPnyY4cOHk5WVxezZswHIz8+nd+/eeHh4sGPHDi5dusSQIUNQSrFgwYJSyyKEKH0Xl/8GwLmq/sh2h4VQxVCnTh21YMECdenSpduOi4yMVI8++qiaMWNGcRZ7T95//31Vt25d4/Off/5ZWVlZqbNnzxqnLVu2TOn1epWenl7s5aanpyugRK8RQtw9Q75BnbBtqBSo7cOWah3H4pXVZ1ixtjyOHj2KnZ3dHccFBwcTHBxcJruK0tPTqVatmvH5rl27aN68OTVr/vN7S48ePcjOziYqKoquXbsWuZzs7GyT3W0ZGRnmCy2EKOTIF3tonnuULBxp9dYjWscRxVSsYx7FKY57GV9Sx48fZ8GCBTz//PPGacnJyXh6epqMc3Nzw87OjuTk5Fsua+bMmbi6uhofPj5ypocQZSl17pcARNd7BGcvZ43TiOK6qyNTe/fuZcuWLaSkpGAwGEzmffDBB8VeztSpU5k2bdptx+zbt4+goCDj83PnzvHAAw/w6KOPMmzYMJOxOp2u0OuVUkVOv2Hy5MlMmDDB+DwjI0MKRIgykp2RTUDscgAcn3ta4zSiJEpcHu+88w6vv/46jRs3xtPT0+SD+XYf0kUZNWoUjz/++G3H+Pn5Gf9+7tw5unbtSnBwMJ9++qnJOC8vL/bs2WMyLS0tjdzc3EJbJP+m1+tL9eC+EKL4Dry1nmB1mSSrWrQcV/SuZVE+lbg8PvzwQ7744gvCw8Pv+c3d3d1xd3cv1tizZ8/StWtX2rRpw5IlS7CyMt3jFhwczIwZM0hKSsLb2xuAiIgI9Ho9bdq0ueesQojSZ/XNVwDEtX0KbztrjdOIkihxeVhZWRESEmKOLLd07tw5unTpQp06dZg9ezYXLlwwzvPy8gIgLCwMf39/Bg8ezKxZs0hNTeWll15i+PDhuLi4lGleIcSdXYy9QOvknwGo/arssrI0Jf6G+fjx4/noo4/MkeWWIiIiOHbsGL///ju1a9fG29vb+LjB2tqaDRs2YG9vT0hICAMHDqRfv37G74EIIcqXI28sx5Y8YhyDaPCQv9ZxRAnplFKqJC8wGAz07t2b+Ph4/P39sbW1NZm/evXqUg1Y1jIyMnB1dSU9PV22WIQwoxinIPyvRrH1kfmE/jBa6zgVRll9hpV4t9Xo0aP5448/6Nq1K9WrVy/xQXIhhDj24xH8r0aRiw3N3rr9STOifCpxeXz11VesWrWK3r17myOPEKISODPzaxoAB7x60b6ph9ZxxF0o8TGPatWqUb9+fXNkEUJUAumnLtN87+cAGJ6SA+WWqsTlMXXqVKZMmcLVq1fNkUcIUcFFP/IW7uoix+2aEDTtIa3jiLtU4t1W8+fP5/jx43h6euLn51fogPmBAwdKLZwQomJJ2BhPSNR8ANLemEt9R9s7vEKUVyUuj379+pkhhhCiMrjw9IvUJY99Hj1p+/oDWscR96DE5TFlyhRz5BBCVHBRMyNol7KeXGyovrT418AT5ZNZbkNbwq+OCCEquLzrebhMK7gAaWSrkdTr1UTjROJeFas8mjZtynfffXfH+3QcPXqUF154gffee69UwgkhKoadT39Cw+wjXNJVp+Vq2XtRERRrt9VHH33EK6+8wsiRIwkLCyMoKIiaNWtib29PWloaMTEx7Nixg5iYGEaNGsWIESPMnVsIYSHSjqfS4oc3ATg8cDqhdd00TiRKQ4kuTxIZGcmKFSvYtm0bJ0+e5Nq1a7i7uxMYGEiPHj146qmnqFq1qhnjmp9cnkSI0rW11VhC/5zPUX1z6l6Oxsb+rm4jJIqpXF6epGPHjnTs2NFcWYQQFczx9bGE/FlwIdUr0+dKcVQgZjlgLoQQAKnPTMCGfPZ4PUTrid20jiNKkZSHEMIs9k37mbYXfyUHWzy/klsjVDRSHkKIUpdxJgPvt14AILLNGPy6N9Q4kShtUh5CiFJ3sNuL1M5PJNGmLm3Wyam5FZGUhxCiVO2b9jOd4z4DIHX2EqrUrKJxImEOd3Xqg8Fg4NixY6SkpGAwGEzmde7cuVSCCSEsT9rxVHymDwNga+A4QseGapxImEuJy2P37t0MGjSIU6dOFboMiU6nIz8/v9TCCSEsS0y3MYQYkkiwbUS7397ROo4woxKXx/PPP09QUBAbNmzA29tbbkMrhABg98TVhJz8lnysyFr4JQ7VHLSOJMyoxOVx9OhRfvjhBxo0aGCOPEIIC3Qx9gL1Zz8PwPbgV+gyrIPGiYS5lfiAefv27Tl27Jg5sgghLJAyKI51ex4PdYF4+xYE/ypnV1UGxdryOHTokPHvo0eP5sUXXyQ5OZkWLVoUupNgQEBA6SYUQpRrkaOXEXJuNbnYYPjiS/Queq0jiTJQrAsjWllZodPpbnmfjhvzKsIBc7kwohDFl3zgHPqg5ripNLZ0nUaX39/UOlKlV64ujJiQkGC2AEIIy6QMitMPDKetSiPGsQ0h6ydrHUmUoWKVh6+vr/Hv27Zto2PHjtjYmL40Ly+PyMhIk7FCiIprx9Al3HfhZ7Kxw+67L7F1tL3zi0SFUeID5l27diU1NbXQ9PT0dLp27VoqoYQQ5VvSvjO0XDoOgF293qZB32baBhJlrsTlcePYxs0uXbqEk5NTqYQSQpRfyqA4/eDzuHCFw07tuW/NBK0jCQ0U+3se/fv3BwoOjoeHh6PX/3NGRX5+PocOHZIbRQlRCUSO+o6QlA1kY4f9d19gbWetdSShgWKXh6urK1Cw5VGlShUcHP759qidnR0dOnRg+PDhpZ9QCFFuXDiSQpOPxwKw6/436PKQv8aJhFaKXR5LliwBwM/Pj5deekl2UQlRCR3tOZqO6hJ/O7QiZN0rWscRGirx5UmmTCn49mhKSgpxcXHodDoaNWpEjRo1Sj2cEKL82P3KGjqe/p48rOHzL+TsqkquxAfMMzIyGDx4MLVq1SI0NJTOnTtTq1YtnnrqKdLT082RUQihscsJadSdPQKAHcETafJEoMaJhNZKXB7Dhg1jz549rF+/nsuXL5Oens769evZv3+/HPMQooL6K2wCnoZkjts1ocPP8i1ycRe7rTZs2MDGjRvp1KmTcVqPHj1YvHgxDzzwQKmGE0Job/+Mjdx3bCkGdFyd/zn2Ve21jiTKgRJveVSvXt145tW/ubq64ubmViqhhBDlw5VzV/Ca8l8AtrcaQ4vn5HR8UaDE5fH6668zYcIEkpKSjNOSk5N5+eWXeeONN0o1nBBCWwd6TKJ2fiKJNnUJ2jhD6ziiHCnxbqtFixZx7NgxfH19qVOnDgCJiYno9XouXLjAJ598Yhx74MCB0ksqhChTfy7YRujhhQBcfGcxdWrI6fniHyUuj379+pkhhhCiPLmWeg2XF4cBsK3JcDq/fL/GiUR5U6z7eVQmcj8PIWBLu4l02TeLJKtaOCYcwbVO4eOconwqq8+wEh/zALh8+TKfffYZkydPNl5h98CBA5w9e7ZUwwkhyt7uiau5b98cAE6/9rEUhyhSiXdbHTp0iG7duuHq6srJkycZPnw41apVY82aNZw6dYqvvvrKHDmFEGVg7xs/0WbWY1hjYFuT/9J5+oNaRxLlVIm3PCZMmEB4eDhHjx7F3v6f87179uzJtm3bSjWcEKLs7H/7V1q+PQBb8tjpO4iQPxdqHUmUYyUuj3379vHcc88Vml6rVi2Sk5NLJZQQomwdmLWZ5m/0Q08Ou2oNoP3fX8ql1sVtlbg87O3tycjIKDQ9Li4ODw+PUgklhCg7Bz/cSpOJfbAnmz1efQmK/w4b+xLv0RaVTInLo2/fvkyfPp3c3Fyg4OZQiYmJTJo0iUceeaTUAwohzOevTyJpMK43jlxjn0cvWsWtkKvlimIpcXnMnj2bCxcuUKNGDa5du0ZoaCgNGjSgSpUqzJgh30AVwlIcWbKXOs/3xJksoqp1o0X8KvQu+ju/UAjuojxcXFzYsWMHq1at4t1332XUqFH8/PPPbN261aw3iHrooYeoU6cO9vb2eHt7M3jwYM6dO2cyJjExkT59+uDk5IS7uztjxowhJyfHbJmEsFSx3x6g1tAeuJJBdNUuNI37US54KErEYr4kOHfuXIKDg/H29ubs2bO89NJLAERGRgIF91Fv1aoVHh4ezJkzh0uXLjFkyBD69+/PggULiv0+8iVBUdHF/3AI94FdqaZSOVQlhHrxv+Ls5ax1LFFKyuozrETlYTAYWLp0KatXr+bkyZPodDrq1q3LgAEDGDx4MDqdzmxBb7Zu3Tr69etHdnY2tra2/PLLLzz44IOcPn2amjVrArB8+XLCw8NJSUkp9kqU8hAV2bF1Mbj264KHusBhp/bU+TsCl9ryc16RlLtvmCuleOihhxg2bBhnz56lRYsWNGvWjFOnThEeHs7DDz9stpA3S01N5dtvv6Vjx47Y2hYc3Nu1axfNmzc3FgcU3GckOzubqKioWy4rOzubjIwMk4cQFVHCxniqPHw/HuoCsQ6t8TnyqxSHuGvFLo+lS5eybds2Nm/eTHR0NMuWLWP58uX8+eef/Pbbb/z+++9m/3b5K6+8gpOTE9WrVycxMZEff/zROC85ORlPT0+T8W5ubtjZ2d32+yczZ87E1dXV+PDx8TFbfiG0cmbHSex7/QdPQzJx9gF4/hmBq29VrWMJC1bs8li2bBmvvvoqXbt2LTTvP//5D5MmTeLbb78t0ZtPnToVnU5328f+/fuN419++WWio6OJiIjA2tqap59+mn/vdStqt5lS6ra70yZPnkx6errxcfr06RL9G4Qo75IPnCO/6/14G85y3K4p1fZvolrD6lrHEhau2N8EOnToEO+///4t5/fs2ZP58+eX6M1HjRrF448/ftsxfn5+xr+7u7vj7u5Oo0aNaNq0KT4+PuzevZvg4GC8vLzYs2ePyWvT0tLIzc0ttEXyb3q9Hr1eTk8UFVPq0Utc6RhGw7wTnLKph9Ou3/BoVkPrWKICKHZ5pKam3vZD2NPTk7S0tBK9+Y0yuBs3tjiys7MBCA4OZsaMGSQlJeHt7Q1AREQEer2eNm3a3NV7CGHJMs5kkNzqAfyzj5BkVROrzb/h1brmnV8oRDEUuzzy8/Oxsbn1cGtra/Ly8kol1M327t3L3r176dSpE25ubpw4cYI333yT+vXrExwcDEBYWBj+/v4MHjyYWbNmkZqayksvvcTw4cPlrClR6Vy9eJUTzfvQ6up+LurcyVqziQad62odS1QgxS4PpRTh4eG33MVzYwvAHBwcHFi9ejVTpkwhKysLb29vHnjgAZYvX27MY21tzYYNGxgxYgQhISE4ODgwaNAgZs+ebbZcQpRHOZk5HPEfQNv0baTjwoWvN9L0IX+tY4kKptjf83jmmWeKtcAlS5bcUyCtyfc8hCXLz8lnb/0nCD6zkqs4cOyjCAJGdNI6lihDZfUZVuwtD0svBSEqOmVQRLZ4jvvOrCQHW2LeXkOQFIcwk7u6Da0QonxRBsW2oAncF/85+VgR9dIygl7roXUsUYFJeQhRAWz9zzRCo+cBsGvo5wTPktsjCPOS8hDCwm3pN5cuW6cBsPWR+XT6LFzbQKJSkPIQwoJtH/IZXX6cAMCW+98i9IfRGicSlYWUhxAWKnLsCkK++i8AW4JeIjTiNY0TicpEykMIC7Rv6gbazn8KKxTbmvyX0D3vo7Mqu1siCCHlIYSFif0miubTBmBLHjt9nyDkz4VSHKLMSXkIYWHOz1uGA9eJt29Bu5gvsbaz1jqSqISkPISwMFZNGgGQ6eCBraOtxmlEZSXlIYSF8Xq44GKgDdP2knfdPBcjFeJOpDyEsDD1+/iTQRWqkMnxHw9rHUdUUlIeQlgYaztrjlbrAMD5tbs0TiMqKykPISzQlRYdAbDetV3jJKKykvIQwgJV7dcFgPqnt6AMxbqrghClSspDCAvUJLwD19HjZUji5KajWscRlZCUhxAWyL6qPbFVC866Ov3VHxqnEZWRlIcQFiq9VRcAbHZs0TSHqJykPISwUMbjHmfkuIcoe1IeQlioJkPacx09noZkEn6N0zqOqGSkPISwUPZV7Tnh2AKACzukPETZkvIQwoJddagOQO6Fy9oGEZWOlIcQFizHyQ2A/ItpGicRlY2UhxAWLM+5KgAqVcpDlC0pDyEsmMGtYLeVLu2SxklEZSPlIYQF03m4A2B7+aLGSURlI+UhhAWz8SooD4crKRonEZWNlIcQFsy5mS8AHldOaJxEVDZSHkJYsFr3NwGgdv4prl68qnEaUZlIeQhhwdybenBJV3DQ/NTGvzVOIyoTKQ8hLNxJt0AALqzfo3ESUZlIeQhh4a4EdALAetcOjZOIykTKQwgL59KroDzqntmGIc+gcRpRWUh5CGHhmg7tSDou1Mw/w8G5cmMoUTakPISwcA7VHPiz2ZMAXF+wWOM0orKQ8hCiAqjx2nAA2p5eRdK+MxqnEZWBlIcQFUCTJwI56BqKLXnEjZyvdRxRCUh5CFFB5Ix5CYDAfZ+QcSZD4zSiopPyEKKCCHqzF8ftmuJKBgdekGMfwrykPISoIKxsrDj3xIsANPplHrlXczVOJCoyKQ8hKpB2858ixcqTmvln2PviCq3jiApMykOICkTvoifm/jEAuC+dhTIojROJikrKQ4gKpuWi58nEicbXD3Hg/d+0jiMqKCkPISoYt/rViGo5tODJrFnahhEVlpSHEBVQ/f+NJx8r2qRuIu77P7WOIyogKQ8hKqDanfzY4/MoAGnjp8uxD1HqpDyEqKCqz55MHtZ0OLearb3e0zqOqGAsrjyys7Np1aoVOp2OgwcPmsxLTEykT58+ODk54e7uzpgxY8jJydEmqBAaazywJTv7zQagy8bJbB/ymcaJREViceUxceJEatasWWh6fn4+vXv3Jisrix07drB8+XJWrVrFiy++qEFKIcqH0DXj2NJhEgAdv3qO3a+s0TiRqCgsqjx++eUXIiIimD17dqF5ERERxMTE8M033xAYGEi3bt2YM2cOixcvJiNDrvMjKq/Qne+wrfEwrDHQ6v0n5PRdUSospjzOnz/P8OHD+frrr3F0dCw0f9euXTRv3txkq6RHjx5kZ2cTFRV1y+VmZ2eTkZFh8hCiItFZ6eh4cBG7vR/GnmyavtJHCkTcM4soD6UU4eHhPP/88wQFBRU5Jjk5GU9PT5Npbm5u2NnZkZycfMtlz5w5E1dXV+PDx8enVLMLUR7Y2NsQ+Pcy9tZ4EAeu0/SVPkS9u0nrWMKCaVoeU6dORafT3faxf/9+FixYQEZGBpMnT77t8nQ6XaFpSqkip98wefJk0tPTjY/Tp0/f879LiPJI76Kn5dEf2OPZBweu4z/5IaJmRmgdS1goGy3ffNSoUTz++OO3HePn58fbb7/N7t270ev1JvOCgoJ48skn+fLLL/Hy8mLPnj0m89PS0sjNzS20RfJver2+0HKFqKj0Lnpaxa9kT6NHaX/+J5q9+hD7slfRdmpvraMJC6NTSpX7bw8lJiaaHIs4d+4cPXr04IcffqB9+/bUrl2bX375hQcffJAzZ87g7e0NwIoVKxgyZAgpKSm4uLgU670yMjJwdXUlPT292K8RwtLkZOYQ3fBR2ievIw9r9jy/lJBFT2kdS5SCsvoM03TLo7jq1Klj8tzZ2RmA+vXrU7t2bQDCwsLw9/dn8ODBzJo1i9TUVF566SWGDx8uJSDETeyc7Wh9/Ad2NnuGkJPfEvLxYLampBK6aozW0YSFsIgD5sVhbW3Nhg0bsLe3JyQkhIEDB9KvX78iT+sVQoCtoy3BR79ia8uCwghdPZYt970hlzIRxWIRu63Kkuy2EpWNMii2hs2gy+Y3ANjW7AU6HfwfVjYV5nfLSqWsPsPkp0OISk5npaPLb6+z7fGFGNDR+cgitg34UOtYopyT8hBCANB52QtsCy64lInjrs0apxHlnZSHEMLIuXcoAFUzEjVOIso7KQ8hhJFby4IzGz2vn5ID5+K2pDyEEEbeHXwBcCWDywlpGqcR5ZmUhxDCyNHdkWSrgi/ZJu08oXEaUZ5JeQghTKQ41QPgcrSUh7g1KQ8hhIkMj4LyyImV8hC3JuUhhDCR51NQHlanpDzErUl5CCFM2DQqKA/n81Ie4takPIQQJlxaFZSH+xUpD3FrUh5CCBM1OhSUR828RHKv5mqcRpRXUh5CCBM1Ary4jh4b8knaK3fWFEWT8hBCGF1OSGN3vSewJxuAa+cz7vAKUVlJeQghAIj+4A+uNgyg4+nvycOaLd3eptGjLbWOJcopi7iToBDCfHIyc4i8/w06752FFYqTtg3IXPQNXYa21zqaKMekPISoxI6vjyVn4JN0uRYNwPZGQwncOg8/L2eNk4nyTnZbCVEJKYNi6+ML8e7ThqbXoknVVWP3xNXcF/cZzlIcohhky0OISubC4fMk3D+U0JQNAERV606tTUvp0LqmxsmEJZEtDyEqkX3TfoaAANqlbOA6erb2m0vg+V/xkuIQJSRbHkJUAtdSr7G3y8uE/vURAEf1zVHffkfoIy00TiYslZSHEBVc3IqDWA95ktDsGAC2thpL+z/exb6qvcbJhCWT3VZCVFCGPANb+n6A3+PtaZAdw3krL/a//Suh0fOkOMQ9ky0PISqg5APnONt9CF1SfwNgj1df6v++mKCmHhonExWFbHkIUcHsfmUNdkEtaJP6G1dxYNuTn9Du7BrcpThEKZItDyEqiPycfHa0HU/ooQUAxDq0Rv/Dt3Tu1UTjZKIikvIQogLIzsjmQLPBhJ5ZiQEd29pPpONv07FzttM6mqigpDyEsHBXzl0hvkV/glN/Iwdb9o/7hi5zB2odS1RwUh5CWLBLcRdJbtOLNln7yMSJ+PfW0nFiN61jiUpAykMIC3V2VyLZXcJolhPHJV11zi/5hdZD2modS1QSUh5CWKDj62Nx6BdGvfwznLX2IXtdBP5yYFyUITlVVwgLc/jzPbg91Ima+Wc4btcU3c6d1JPiEGVMykMICxI1MwK/YfdTTaVy2Kkdboe3U7O9j9axRCUk5SGEhYgcu4IWrz6IM1nsrx6G37HNVGtYXetYopKS8hDCAmx9fCEd5j+BHblE+jxGwMmf5KZNQlNSHkKUY8qg2NJ1GqErRmKFYmvzEbQ/9q18+U9oTs62EqKcMuQZ2N56LF3++h8AW7pMJXTzm+isdBonE0LKQ4hyKSsli4PthhN6ahkGdGx/dAFdvh+pdSwhjGS3lRDliDIoIkcvI71mE0JOLSMHW3aP/o5QKQ5RzsiWhxDlROw3UeSOGEvHKzsBOG3tx4WZn9Hx5fs1TiZEYVIeQmgs5VAycQNeI+ToEqxQZOHIvm6v0n7FBHyqOWgdT4giyW4rITSSnZHNlt6zsG/ZiPuOfoEVih11B5OxL54um17DQYpDlGOy5SFEGVMGxb4p66nx3gS65B4D4LBTO9TcD+k0vIPG6YQoHikPIcrQsXUxXH52PO0uRQBw3sqL+KHvEbLwKaxsZEeAsBxSHkKUgbTjqRzqP5WQQwuxIZ9s7IgMfpGgHyZzX80qWscTosTkVx0hzCjveh5bH1+IatiQ0EMLsCGf3d4Pk/x7LF0j36GKFIewULLlIYSZRM/5HefXxhKafRiAeH1zrkyfR4eJcuqtsHwWs+Xh5+eHTqczeUyaNMlkTGJiIn369MHJyQl3d3fGjBlDTk6ORolFZZW45QS7a/Yn8KX7aZh9mFRdNbY+9hH1LkfTRopDVBAWteUxffp0hg8fbnzu7PzPVUXz8/Pp3bs3Hh4e7Nixg0uXLjFkyBCUUixYsECLuKKSuXLuClEDZhK8aw51yCEPa3YGjCBg9VRC61fTOp4QpcqiyqNKlSp4eXkVOS8iIoKYmBhOnz5NzZo1AZgzZw7h4eHMmDEDFxeXsowqKhFDnoHIEd/Q4PNJdDEkARBVrTuuX8wltG8zjdMJYR4Ws9sK4L333qN69eq0atWKGTNmmOyS2rVrF82bNzcWB0CPHj3Izs4mKirqlsvMzs4mIyPD5CFEseTkkPbpSmKrBtNp8RC8DEmcsqnPnld/pPWFjTSQ4hAVmMVseYwdO5bWrVvj5ubG3r17mTx5MgkJCXz22WcAJCcn4+npafIaNzc37OzsSE5OvuVyZ86cybRp08yaXVQgeXmwbRusXg3ff4/bhQu4ARlU4UCvNwheNgZfF73WKYUwO03LY+rUqXf84N63bx9BQUGMHz/eOC0gIAA3NzcGDBhg3BoB0OkK3+dAKVXk9BsmT57MhAkTjM8zMjLw8ZF7Qot/OXsWtm6FTZtg3TpITTXOynTxZl7GsyQ/Mor//VD0LlUhKiJNy2PUqFE8/vjjtx3j5+dX5PQOHQou43Ds2DGqV6+Ol5cXe/bsMRmTlpZGbm5uoS2Sf9Pr9ej18pui+H8ZGXDoEBw8CAcOwI4dcPSo6Rh3d3joIejfn0WHevDGqzaEy9c1RCWjaXm4u7vj7u5+V6+Njo4GwNvbG4Dg4GBmzJhBUlKScVpERAR6vZ42bdqUTmBRMVy7BufOwYkTcOwYHD9e8OfhwwV/v5mVFQQGQmgo9OkDnTqBTcH/Orl/FgyxsZgdwEKUDov4kd+1axe7d++ma9euuLq6sm/fPsaPH89DDz1EnTp1AAgLC8Pf35/Bgwcza9YsUlNTeemllxg+fLicaVVR5eYWbCmkp5s+/j3t8mVISvrnce5cwfTbqV0bWrUqKIx27eC++8DVtciheXkFf0p5iMrGIn7k9Xo9K1asYNq0aWRnZ+Pr68vw4cOZOHGicYy1tTUbNmxgxIgRhISE4ODgwKBBg5g9e/bdvenateDgAEoVPOCfv99umrnGVob3MhggJweuX4fs7MJ//vvvmZkFWxB3y94e6taF+vWhQYOCR+PGBaVRgq1hKQ9RWVnEj3zr1q3ZvXv3HcfVqVOH9evX39N7qf//QMsYMuSeliPKkIMDuLjc+uHpCV5eBX96exf86eoKtzqRogSna2dlFfyZl1eilwlhNje+bnDjs8xcdMrc72BhTpw4Qf369bWOIYQQ9+T48ePUq1fPbMu3iC2PslStWsFlJBITE3G9xX7u8uLGacWnT5+2iOM6lpTXkrKC5DUnS8oKkJ6eTp06dYyfZeYi5XETK6uCL927urpaxA8KgIuLi8VkBcvKa0lZQfKakyVlhX8+y8y2fLMuXQghRIUk5SGEEKLEpDxuotfrmTJlikV869ySsoJl5bWkrCB5zcmSskLZ5ZWzrYQQQpSYbHkIIYQoMSkPIYQQJSblIYQQosSkPIQQQpRYhS+PhQsXUrduXezt7WnTpg3bt2+/7fitW7fSpk0b7O3tqVevHh9//HGhMatWrcLf3x+9Xo+/vz9r1qwpt3mXLl2KTqcr9Lh+/XqZZk1KSmLQoEE0btwYKysrxo0bV+S48rJui5O3vKzb1atX0717dzw8PHBxcSE4OJiNGzcWGlde1m1x8ppz3ZY0744dOwgJCaF69eo4ODjQpEkT5s6dW2icudZvaWcttXWrKrDly5crW1tbtXjxYhUTE6PGjh2rnJyc1KlTp4ocf+LECeXo6KjGjh2rYmJi1OLFi5Wtra364YcfjGMiIyOVtbW1euedd1RsbKx65513lI2Njdq9e3e5zLtkyRLl4uKikpKSTB5lnTUhIUGNGTNGffnll6pVq1Zq7NixhcaUp3VbnLzlZd2OHTtWvffee2rv3r0qPj5eTZ48Wdna2qoDBw4Yx5SndVucvOZat3eT98CBA+q7775Thw8fVgkJCerrr79Wjo6O6pNPPjGOMdf6NUfW0lq3Fbo82rVrp55//nmTaU2aNFGTJk0qcvzEiRNVkyZNTKY999xzqkOHDsbnAwcOVA888IDJmB49eqjHH3+8XOZdsmSJcnV1veds95r130JDQ4v8MC5P6/bfbpW3PK7bG/z9/dW0adOMz8vrur3h5rzmWrdKlU7ehx9+WD311FPG5+Zav+bIWlrrtsLutsrJySEqKoqwsDCT6WFhYURGRhb5ml27dhUa36NHD/bv309ubu5tx9xqmVrnBcjMzMTX15fatWvz4IMPGu/CWJZZi6M8rdviKo/r1mAwcOXKFZML45XndVtUXij9dVtaeaOjo4mMjCQ0NNQ4zRzr11xZoXTWbYUtj4sXL5Kfn1/o/uWenp4kJycX+Zrk5OQix+fl5XHx4sXbjrnVMrXO26RJE5YuXcq6detYtmwZ9vb2hISEcPTm+3KbOWtxlKd1Wxzldd3OmTOHrKwsBg4caJxWntdtUXnNsW7vNW/t2rXR6/UEBQUxcuRIhg0bZpxnjvVrrqyltW4r/FV1dTfd8EcpVWjancbfPL2kyyyJ0s7boUMHOnToYJwfEhJC69atWbBgAfPnzy/TrFot01zLLo/rdtmyZUydOpUff/yRGjVqlMoyi6O085pz3d5t3u3bt5OZmcnu3buZNGkSDRo04IknnrinZWqRtbTWbYUtD3d3d6ytrQs1dEpKSqEmv8HLy6vI8TY2NlSvXv22Y261TK3z3szKyoq2bdve029wd5O1OMrTur0bWq/bFStWMHToUFauXEm3bt1M5pXHdXu7vDcrjXV7r3nr1q0LQIsWLTh//jxTp041fiCbY/2aK+vN7nbdVtjdVnZ2drRp04ZNmzaZTN+0aRMdO3Ys8jXBwcGFxkdERBAUFIStre1tx9xqmVrnvZlSioMHD+Lt7V2mWYujPK3bu6Hlul22bBnh4eF899139O7du9D88rZu75T3ZqWxbu8lb1F5srOzjc/NsX7NlbWo+Xe1bu/5kHs5duM0t88//1zFxMSocePGKScnJ3Xy5EmllFKTJk1SgwcPNo6/cerr+PHjVUxMjPr8888Lnfq6c+dOZW1trd59910VGxur3n333VI/5bE0806dOlX9+uuv6vjx4yo6Olo988wzysbGRu3Zs6dMsyqlVHR0tIqOjlZt2rRRgwYNUtHR0erIkSPG+eVp3RYnb3lZt999952ysbFRH330kcmpl5cvXzaOKU/rtjh5zbVu7ybv//73P7Vu3ToVHx+v4uPj1RdffKFcXFzUa6+9ZhxjrvVrjqyltW4rdHkopdRHH32kfH19lZ2dnWrdurXaunWrcd6QIUNUaGioyfgtW7aowMBAZWdnp/z8/NSiRYsKLXPlypWqcePGytbWVjVp0kStWrWq3OYdN26cqlOnjrKzs1MeHh4qLCxMRUZGapIVKPTw9fU1GVOe1u2d8paXdRsaGlpk1iFDhpgss7ys2+LkNee6LWne+fPnq2bNmilHR0fl4uKiAgMD1cKFC1V+fr7JMs21fks7a2mtW7kkuxBCiBKrsMc8hBBCmI+UhxBCiBKT8hBCCFFiUh5CCCFKTMpDCCFEiUl5CCGEKDEpDyGEECUm5SFEBXXy5El0Oh0HDx687bi4uDi8vLy4cuVK2QSj4PpMHh4enD17tszeU5QuKQ9RJpKTkxk9ejT16tVDr9fj4+NDnz592Lx5s9bRbkun07F27VqtY5jVa6+9xsiRI6lSpUqpLK9Lly5F3r7532rUqMHgwYOZMmVKqbynKHtSHsLsTp48SZs2bfj99995//33+euvv/j111/p2rUrI0eOvOvlKqXIy8srND0nJ+de4lYqZ86cYd26dTzzzDOlsrzU1FQiIyPp06fPHcc+88wzfPvtt6SlpZXKe4uyJeUhzG7EiBHodDr27t3LgAEDaNSoEc2aNWPChAns3r0bKHoXy+XLl9HpdGzZsgWALVu2oNPp2LhxI0FBQej1erZv306XLl0YNWoUEyZMwN3dne7duwMQExNDr169cHZ2xtPTk8GDBxtvkgUFvyGPGTOGiRMnUq1aNby8vJg6dapxvp+fHwAPP/wwOp3O+PxmOTk5jBo1Cm9vb+zt7fHz82PmzJnG+TqdjkWLFtGzZ08cHByoW7cuK1euNFnG2bNneeyxx3Bzc6N69er07duXkydPmoxZsmQJTZs2xd7eniZNmrBw4UKT+Xv37iUwMBB7e3uCgoKKdXe477//npYtW1K7dm3jtKVLl1K1alXWr19P48aNcXR0ZMCAAWRlZfHll1/i5+eHm5sbo0ePJj8/32R5GzZsoGXLltSqVYu0tDSefPJJPDw8cHBwoGHDhixZssQ4tkWLFnh5ebFmzZo75hTl0N1frkuIO7t06ZLS6XTqnXfeue24hIQEBajo6GjjtLS0NAWoP/74Qyml1B9//KEAFRAQoCIiItSxY8fUxYsXVWhoqHJ2dlYvv/yy+vvvv1VsbKw6d+6ccnd3V5MnT1axsbHqwIEDqnv37qpr167G5YeGhioXFxc1depUFR8fr7788kul0+lURESEUkqplJQUBaglS5aopKQklZKSUmT2WbNmKR8fH7Vt2zZ18uRJtX37dvXdd98Z5wOqevXqavHixSouLk69/vrrytraWsXExCillMrKylINGzZUzz77rDp06JCKiYlRgwYNUo0bN1bZ2dlKKaU+/fRT5e3trVatWqVOnDihVq1apapVq6aWLl2qlFIqMzNTeXh4qMcee0wdPnxY/fTTT6pevXqF1unN+vbtW+ge2UuWLFG2traqe/fu6sCBA2rr1q2qevXqKiwsTA0cOFAdOXJE/fTTT8rOzk4tX77c5LUDBgxQb731llJKqZEjR6pWrVqpffv2qYSEBLVp0ya1bt06k/EDBw5U4eHht8wnyi8pD2FWe/bsUYBavXr1bceVpDzWrl1r8trQ0FDVqlUrk2lvvPGGCgsLM5l2+vRpBai4uDjj6zp16mQypm3btuqVV14xPgfUmjVrbpt99OjR6j//+Y8yGAxFzgcKfUC3b99evfDCC0oppT7//HPVuHFjk9dnZ2crBwcHtXHjRqWUUj4+PiaFpJRSb731lgoODlZKKfXJJ5+oatWqqaysLOP8RYsW3bE8WrZsqaZPn24ybcmSJQpQx44dM0577rnnlKOjo7py5YpxWo8ePdRzzz1nfH79+nVVpUoVdejQIaWUUn369FHPPPPMLd9bKaXGjx+vunTpctsxonyqsHcSFOWDKuI2vvcqKCjojtOioqL4448/cHZ2LjT2+PHjNGrUCICAgACTed7e3qSkpJQoT3h4ON27d6dx48Y88MADPPjgg4SFhZmMCQ4OLvT8xi66qKgojh07VuiA9fXr1zl+/DgXLlzg9OnTDB06lOHDhxvn5+Xl4erqCkBsbCwtW7bE0dHxlu9ZlGvXrmFvb19ouqOjI/Xr1zc+9/T0xM/Pz2R9enp6mqyr33//nerVq9OiRQsAXnjhBR555BEOHDhAWFgY/fr1K3QTIwcHB65evXrHnKL8kfIQZtWwYUN0Oh2xsbH069fvluOsrAoOv6l/3SEgNze3yLFOTk53nGYwGOjTpw/vvfdeobH/vmPazXdc1Ol0GAyGW+YsSuvWrUlISOCXX37ht99+Y+DAgXTr1o0ffvjhtq+7UagGg4E2bdrw7bffFhrj4eHB9evXAVi8eDHt27c3mW9tbQ2YrreScHd3L/KAdVHr5U7rat26dfTt29f4vGfPnpw6dYoNGzbw22+/cf/99zNy5Ehmz55tHJOamoqHh8ddZRfakgPmwqyqVatGjx49+Oijj8jKyio0//LlywDGD5CkpCTjvDt9P+F2WrduzZEjR/Dz86NBgwYmj6LK51ZsbW0LHRQuiouLC4899hiLFy9mxYoVrFq1itTUVOP8GycG/Pt5kyZNjFmPHj1KjRo1CmV1dXXF09OTWrVqceLEiULzb9yr2t/fnz///JNr167d8j2LEhgYSExMTLHWxe0opfjpp5946KGHTKZ7eHgQHh7ON998w7x58/j0009N5h8+fJjAwMB7fn9R9qQ8hNktXLiQ/Px82rVrx6pVqzh69CixsbHMnz/fuGvFwcGBDh068O677xITE8O2bdt4/fXX7/o9R44cSWpqKk888QR79+7lxIkTRERE8OyzzxarDG7w8/Nj8+bNJCcn3/KU0rlz57J8+XL+/vtv4uPjWblyJV5eXlStWtU4ZuXKlXzxxRfEx8czZcoU9u7dy6hRowB48skncXd3p2/fvmzfvp2EhAS2bt3K2LFjOXPmDABTp05l5syZfPjhh8THx/PXX3+xZMkSPvjgAwAGDRqElZUVQ4cOJSYmhp9//tnkN/xb6dGjB7t27SrROilKVFQUWVlZdO7c2TjtzTff5Mcff+TYsWMcOXKE9evX07RpU+P8q1evEhUVVWgXn7AMUh7C7OrWrcuBAwfo2rUrL774Is2bN6d79+5s3ryZRYsWGcd98cUX5ObmEhQUxNixY3n77bfv+j1r1qzJzp07yc/Pp0ePHjRv3pyxY8fi6upq3EVWHHPmzGHTpk34+Pjc8jdkZ2dn3nvvPYKCgmjbti0nT57k559/NnmfadOmsXz5cgICAvjyyy/59ttv8ff3BwqOL2zbto06derQv39/mjZtyrPPPsu1a9dwcXEBYNiwYXz22WcsXbqUFi1aEBoaytKlS41bHs7Ozvz000/ExMQQGBjIa6+9VuQuu5v16tULW1tbfvvtt2Kvk6L8+OOP9O7dGxubf/aE29nZMXnyZAICAujcuTPW1tYsX77c5DV16tThvvvuu6f3FtqQ29AKYWY6nY41a9bc9piPlhYuXMiPP/7Ixo0b73oZAQEBvP766wwcOLDYr2nXrh3jxo1j0KBBd/2+QjtywFyISu6///0vaWlpXLly5a4uUZKTk8MjjzxCz549i/2alJQUBgwYwBNPPFHi9xPlg2x5CGFm5X3LQ4i7IVseQpiZ/H4mKiI5YC6EEKLEpDyEEEKUmJSHEEKIEpPyEEIIUWJSHkIIIUpMykMIIUSJSXkIIYQoMSkPIYQQJSblIYQQosT+D5Au+4Ebzt5JAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 393.701x393.701 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from pylag.processing.plot import create_figure\n",
    "\n",
    "fs = 10\n",
    "\n",
    "fig, ax = create_figure(font_size=fs) \n",
    "ax.plot(curr_not_corrected, z_not_corrected, c='b', label=\"Extrapolated\")\n",
    "ax.plot(curr_corrected, z_corrected, c='r', label=\"Log corrected\")\n",
    "ax.set_xlabel('Current speed (m/s)')\n",
    "ax.set_ylabel('Depth (m)')\n",
    "ax.legend(prop={'size': fs})\n",
    "_ = ax.set_xlim([0, 0.35])\n",
    "_ = ax.set_ylim([h, zeta])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4de4f02-aaf5-4541-9a56-b82d84a33a8e",
   "metadata": {},
   "source": [
    "The impact of the correction is just visible toward the bottom of the water column. The extrapolated surface current near to the surface is evident in both plots."
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}