{ "cells": [ { "cell_type": "markdown", "id": "637d7016-f6fc-40d1-a28d-2afa6fb50494", "metadata": {}, "source": [ "# 2D Shallow Water Equations 🌊\n", "\n", "This notebook is a Devito reimplementation from this [notebook](https://github.com/daniel-koehn/Differential-equations-earth-system/blob/e291791b6d6b32f776af01532b37a9869c5cf569/10_Shallow_Water_Equation_2D/01_2D_Shallow_Water_Equations.ipynb). We will use an approximation to the Navier-Stokes equations - the 2D Shallow Water equations - in order to model the propagation of Tsunami events.\n", "\n", "I've only made minor modifications to the text, in order to make it more compatible to the new coding format.\n", "\n", "Author: Átila Saraiva Q. Soares." ] }, { "cell_type": "markdown", "id": "e8035395-5fb0-4d9f-b8cd-b2585c9a9aec", "metadata": {}, "source": [ "## Governing Equations\n", "\n", "Starting from the continuity and momentum conservation equations, we want to model the following problem:\n", "\n", "\n", "\n", "For a given bathymetry model, which can include a complex seafloor topography, we want to model the amplitudes, speed and interaction of waves at the seasurface. At a given point $(x,\\; y)$, the thickness of the water column between the seafloor and undisturbed water surface is defined by the variable $h$, while the wave amplitude is $\\eta$ and therefore the whole thickness of the water column $D = h + \\eta$.\n", "\n", "Using appropriate boundary conditions at the water surface/seafloor, assuming that the horizontal wavelength of the modelled waves are much larger than the water depth and integrating the conservation of mass and momentum equations over the water column, we can derive the following equations to decribe wave propagation \n", "\n", "\\begin{equation}\n", "\\begin{split}\n", "\\frac{\\partial \\eta}{\\partial t} &+ \\frac{\\partial M}{\\partial x} + \\frac{\\partial N}{\\partial y} = 0\\\\\n", "\\frac{\\partial M}{\\partial t} &+ \\frac{\\partial}{\\partial x} \\biggl(\\frac{M^2}{D}\\biggr) + \\frac{\\partial}{\\partial y} \\biggl(\\frac{MN}{D}\\biggr) + g D \\frac{\\partial \\eta}{\\partial x} + \\frac{g \\alpha^2}{D^{7/3}} M \\sqrt{M^2+N^2} = 0\\\\\n", "\\frac{\\partial N}{\\partial t} &+ \\frac{\\partial}{\\partial x} \\biggl(\\frac{MN}{D}\\biggr) + \\frac{\\partial}{\\partial y} \\biggl(\\frac{N^2}{D}\\biggr) + g D \\frac{\\partial \\eta}{\\partial y} + \\frac{g \\alpha^2}{D^{7/3}} N \\sqrt{M^2+N^2} = 0\n", "\\end{split}\n", "\\tag{1}\n", "\\end{equation}\n", "\n", "known as **2D Shallow Water Equations (SWE)**. The derivation of these equations is beyond the scope of this notebook. Therefore, I refer to the [Tsunami Modelling Handbook](http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf) and the lecture [Shallow Water Derivation and Applications by Christian Kühbacher](http://www.mathematik.tu-dortmund.de/lsiii/cms/papers/Kuehbacher2009.pdf) for further details.\n", "\n", "In Eq. (1) the discharge fluxes $M,\\; N$ in x- and y-direction, respectively are given by\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", "M & = \\int_{-h}^\\eta u dz = u(h+\\eta) = uD\\\\\n", "N & = \\int_{-h}^\\eta v dz = v(h+\\eta) = vD\\\\\n", "\\end{split}\n", "\\tag{2}\n", "\\end{equation}\n", "\n", "with the horizontal velocity components $u,\\;v$ in x- and y-direction, while $g$ denotes the gravity acceleration. The terms $\\frac{g \\alpha^2}{D^{7/3}} M \\sqrt{M^2+N^2}$ and $\\frac{g \\alpha^2}{D^{7/3}} N \\sqrt{M^2+N^2}$ describe the influence of seafloor friction on the wave amplitude. $\\alpha$ denotes the Manning's roughness which can be as small as 0.01 for neat cement or smooth metal up to 0.06 for very poor natural channels (see [Tsunami modelling handbook](http://www.tsunami.civil.tohoku.ac.jp/hokusai3/J/projects/manual-ver-3.1.pdf)).\n", "\n", "The Shallow Water Equations can be applied to \n", "\n", "- Tsunami prediction\n", "- Atmospheric flow\n", "- Storm surges\n", "- Flows around structures\n", "- Planetary flows\n", "\n", "and easily extended to incorporate the effects of Coriolis forces, tides or wind stress." ] }, { "cell_type": "markdown", "id": "f93f5cf1-7a65-49a4-82b6-0940c9ab44d0", "metadata": {}, "source": [ "## Operator implementation\n", "\n", "This is a simple function which returns the operator that solves this equation. The important part is the operator, which contains all the equations expressed above." ] }, { "cell_type": "code", "execution_count": 1, "id": "23c0a5bd-d332-46ba-a123-2f2dcc7a4a9b", "metadata": {}, "outputs": [], "source": [ "from devito import Eq, TimeFunction, sqrt, Function, Operator, Grid, solve, ConditionalDimension\n", "from matplotlib import pyplot as plt\n", "import numpy as np\n", "\n", "\n", "def ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid):\n", " \"\"\"\n", " Operator that solves the equations expressed above.\n", " It computes and returns the discharge fluxes M, N and wave height eta from\n", " the 2D Shallow water equation using the FTCS finite difference method.\n", " \n", " Parameters\n", " ----------\n", " eta : TimeFunction\n", " The wave height field as a 2D array of floats.\n", " M : TimeFunction\n", " The discharge flux field in x-direction as a 2D array of floats.\n", " N : TimeFunction\n", " The discharge flux field in y-direction as a 2D array of floats.\n", " h : Function\n", " Bathymetry model as a 2D array of floats.\n", " D : Function\n", " Total thickness of the water column.\n", " g : float\n", " gravity acceleration.\n", " alpha : float\n", " Manning's roughness coefficient.\n", " etasave : TimeFunction\n", " Function that is sampled in a different interval than the normal propagation\n", " and is responsible for saving the snapshots required for the following\n", " animations.\n", " \"\"\"\n", " \n", " eps = np.finfo(grid.dtype).eps\n", " \n", " # Friction term expresses the loss of amplitude from the friction with the seafloor\n", " frictionTerm = g * alpha**2 * sqrt(M**2 + N**2) / D**(7./3.)\n", "\n", " # System of equations\n", " pde_eta = Eq(eta.dt + M.dxc + N.dyc)\n", " pde_M = Eq(M.dt + (M**2/D).dxc + (M*N/D).dyc + g*D*eta.forward.dxc + frictionTerm*M)\n", " pde_N = Eq(N.dt + (M.forward*N/D).dxc + (N**2/D).dyc + g*D*eta.forward.dyc + g * alpha**2 * sqrt(M.forward**2 + N**2) / D**(7./3.)*N)\n", "\n", " stencil_eta = solve(pde_eta, eta.forward)\n", " stencil_M = solve(pde_M, M.forward)\n", " stencil_N = solve(pde_N, N.forward)\n", "\n", " # Equations with the forward in time term isolated\n", " update_eta = Eq(eta.forward, stencil_eta, subdomain=grid.interior)\n", " update_M = Eq(M.forward, stencil_M, subdomain=grid.interior)\n", " update_N = Eq(N.forward, stencil_N, subdomain=grid.interior)\n", " eq_D = Eq(D, eta.forward + h)\n", "\n", " return Operator([update_eta, update_M, update_N, eq_D] + [Eq(etasave, eta)])" ] }, { "cell_type": "markdown", "id": "eadd14a7-e8c5-411c-8847-00fb9f72380b", "metadata": {}, "source": [ "The following function transforms the saved wavefield snapshots and transform them into a video compatible with jupyter notebook." ] }, { "cell_type": "code", "execution_count": 2, "id": "d55ba615-2aa8-4e09-a16d-8c3eb84eade9", "metadata": {}, "outputs": [], "source": [ "from IPython.display import HTML, display\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "\n", "\n", "def snaps2video (eta, title):\n", " fig, ax = plt.subplots()\n", " matrice = ax.imshow(eta.data[0, :, :].T, vmin=-1, vmax=1, cmap=\"seismic\")\n", " plt.colorbar(matrice)\n", "\n", " plt.xlabel('x')\n", " plt.ylabel('z')\n", " plt.title(title) \n", "\n", " def update(i):\n", " matrice.set_array(eta.data[i, :, :].T)\n", " return matrice,\n", "\n", " # Animation\n", " ani = animation.FuncAnimation(fig, update, frames=nsnaps, interval=50, blit=True)\n", "\n", " plt.close(ani._fig)\n", " display(HTML(ani.to_html5_video()))" ] }, { "cell_type": "code", "execution_count": 3, "id": "7c2913ce", "metadata": {}, "outputs": [], "source": [ "def plotDepthProfile (h, title):\n", " fig, ax = plt.subplots()\n", " matrice = ax.imshow(h0)\n", " plt.colorbar(matrice)\n", " plt.xlabel('x')\n", " plt.ylabel('z')\n", " plt.title(title)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "32ba2ef5-9e17-4668-ad81-02d85a6b66dd", "metadata": {}, "source": [ "## Example I: Tsunami in ocean with constant depth\n", "\n", "After writing the `Shallow_water_2D` code and all the required functions attached to it, we can define and run our first 2D Tsunami modelling run.\n", "\n", "Let's assume that the ocean model is $ L_x = 100\\; m$ in x-direction and $L_y = 100\\; m$ in y-direction. The model is discretized with $nx=401$ gridpoints in x-direction and $ny=401$ gridpoints in y-direction, respectively.\n", "\n", "In this first modelling run, we assume a constant bathymetry $h=50\\;m$. The initial wave height field $\\eta_0$ is defined as a Gaussian at the center of the model, with a half-width of 10 m and an amplitude of 0.5 m. Regarding the initial discharge fluxes, we assume that \n", "\n", "\\begin{equation}\n", "\\begin{split}\n", "M_0(x,y) &= 100 \\eta_0(x,y)\\\\\n", "N_0(x,y) &= 0\\\\\n", "\\end{split}\\notag\n", "\\end{equation}\n", "\n", "In order to avoid the occurence of high frequency artifacts, when waves are interacting with the boundaries, boundary conditions will be avoided. Normally Dirichlet conditions are used for the M and N discharge fluxes, and Neumann for the wave height field. Those lead to significant boundary reflections which might be not realistic for a given problem.\n", "\n", "Let's assume the gravity $g = 9.81$ and the Manning's roughness coefficient $\\alpha = 0.025$ for the all the remaining examples." ] }, { "cell_type": "code", "execution_count": 4, "id": "9b704039-286f-4cb7-9e5d-465584c3744a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.00022222222222222223 13500\n" ] } ], "source": [ "Lx = 100.0 # width of the mantle in the x direction []\n", "Ly = 100.0 # thickness of the mantle in the y direction []\n", "nx = 401 # number of points in the x direction\n", "ny = 401 # number of points in the y direction\n", "dx = Lx / (nx - 1) # grid spacing in the x direction []\n", "dy = Ly / (ny - 1) # grid spacing in the y direction []\n", "g = 9.81 # gravity acceleration [m/s^2]\n", "alpha = 0.025 # friction coefficient for natural channels in good condition\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 3.\n", "dt = 1/4500.\n", "nt = (int)(Tmax/dt)\n", "print(dt, nt)\n", "\n", "x = np.linspace(0.0, Lx, num=nx)\n", "y = np.linspace(0.0, Ly, num=ny)\n", "\n", "# Define initial eta, M, N\n", "X, Y = np.meshgrid(x,y) # coordinates X,Y required to define eta, h, M, N\n", "\n", "# Define constant ocean depth profile h = 50 m\n", "h0 = 50. * np.ones_like(X)\n", "\n", "# Define initial eta Gaussian distribution [m]\n", "eta0 = 0.5 * np.exp(-((X-50)**2/10)-((Y-50)**2/10))\n", "\n", "# Define initial M and N\n", "M0 = 100. * eta0\n", "N0 = 0. * M0\n", "D0 = eta0 + 50.\n", "\n", "grid = Grid(shape=(ny, nx), extent=(Ly, Lx), dtype=np.float64)" ] }, { "cell_type": "code", "execution_count": 5, "id": "a1790290-c565-4220-a1d5-c5d32a387db3", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 34\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 17.97 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=17.90456999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.061892000000000044, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "nsnaps = 400\n", "\n", "# Defining symbolic functions\n", "eta = TimeFunction(name='eta', grid=grid, space_order=2)\n", "M = TimeFunction(name='M', grid=grid, space_order=2)\n", "N = TimeFunction(name='N', grid=grid, space_order=2)\n", "h = Function(name='h', grid=grid)\n", "D = Function(name='D', grid=grid)\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "\n", "# Compile the operator\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator\n", "op.apply(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 6, "id": "c5fff152", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling a tsunami in ocean with constant depth\")" ] }, { "cell_type": "markdown", "id": "2b719339-5442-4caf-8682-29045758c7e6", "metadata": {}, "source": [ "Let's take a look on what the code looks like" ] }, { "cell_type": "code", "execution_count": 7, "id": "d27a6e45-3009-4bd7-841b-14783312fadb", "metadata": {}, "outputs": [], "source": [ "# To look at the code, uncomment the line below\n", "#print(op.ccode)" ] }, { "cell_type": "markdown", "id": "9a447c4b-9bc0-450d-af34-ec5f5eaf6b66", "metadata": {}, "source": [ "## Example II: Two Tsunamis in ocean with constant depth\n", "\n", "In example II, we will model the propagation and interaction of two Tsunamis, where the initial conditions for the wave height field $\\eta_0$ consists of Gaussian functions with opposite sign located at $(x_1,\\;y_1) = (35\\; m,\\; 35\\; m)$ and $(x_2,\\;y_2) = (65\\; m,\\; 65\\; m)$. All other modelling parameters remain the same:" ] }, { "cell_type": "code", "execution_count": 8, "id": "8c71df64-a160-44ec-8bf1-bc6eef334919", "metadata": {}, "outputs": [], "source": [ "# Define constant ocean depth profile h = 50 m\n", "h0 = 50 * np.ones_like(X)\n", "\n", "# Define initial Gaussian eta distribution [m]\n", "eta0 = 0.5 * np.exp(-((X-35)**2/10)-((Y-35)**2/10)) # first Tsunami source\n", "eta0 -= 0.5 * np.exp(-((X-65)**2/10)-((Y-65)**2/10)) # add second Tsunami source\n", "\n", "# Define initial M and N\n", "M0 = 100. * eta0\n", "N0 = 0. * M0\n", "D0 = eta0 + h0\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 3.5\n", "dt = 1/8000.\n", "nt = (int)(Tmax/dt)" ] }, { "cell_type": "code", "execution_count": 9, "id": "599ada28-8ff4-49d2-865b-f1c37814913a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 70\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 37.93 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=37.860470999999926, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.06278300000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "# Compile the operator\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator. No need to recompile since the etasave factor is smaller\n", "# than before\n", "op(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 10, "id": "195e4fa1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling two tsunamis in ocean with constant depth\")" ] }, { "cell_type": "markdown", "id": "c8edc6c1-9f2f-40e0-b042-005ce38edc56", "metadata": {}, "source": [ "## Example III: Tsunami in an ocean with 1D $\\tanh$ depth variation\n", "\n", "So far, so good, we can achieve stable and accurate modelling results. However, a constant bathymetry model is a little too simple. In this example, we assume that the bathymetry decreases with a $\\tanh$ function from the left to the right boundary in x-direction. Let 's place a Tsunami source at $(x_1,\\;y_1) = (30\\; m,\\; 50\\; m)$ and see what 's happening:" ] }, { "cell_type": "code", "execution_count": 11, "id": "66379acd-139c-4ac9-82eb-0aed1a1611d9", "metadata": {}, "outputs": [], "source": [ "# Define depth profile h [m]\n", "h0 = 50 - 45 * np.tanh((X-70.)/8.)\n", "\n", "# Define initial eta [m]\n", "eta0 = 0.5 * np.exp(-((X-30)**2/10)-((Y-50)**2/20))\n", "\n", "# Define initial M and N\n", "M0 = 100. * eta0\n", "N0 = 0. * M0\n", "D0 = eta0 + h0\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 2.\n", "dt = 1/4000.\n", "nt = (int)(Tmax/dt)" ] }, { "cell_type": "markdown", "id": "4ac035d5", "metadata": {}, "source": [ "Let's take a look into how this depth profile $h$ looks in the plot below." ] }, { "cell_type": "code", "execution_count": 12, "id": "a9dd81a8", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEWCAYAAAAQBZBVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdyElEQVR4nO3dfZRdVZ3m8e+TIiQI2MiLDiRRohNpwSXRiQw2rYMvNJFWaXu1Thh11KGNTmMLjrMEdE2ja3W6o8uXtte0L1FoWSOCGZCeyNgiqLTDEnkVeQs0URDKxAABJIgkqbrP/HFOmUtRVfdW5Z5T59Z9PmudVffue865+1TgV3ufffb+yTYREYNm3mxXICJiNiT4RcRASvCLiIGU4BcRAynBLyIGUoJfRAykBL8BIeldkq6e7XpMRdK9kl43C98rSf8o6RFJ10l6paS7ZrteUa0EvwaS9Hjb1pL027b3b6upDveW37td0qOSfiTpfZJ68t+MpK9K+us9OP5dkkbL38ljkm6W9IYZnu4PgROAxbaPsf3/bB8x07pFf0jwayDb+41twH3AG9vKLqixKm+0vT/wPGAtcCZwbo3f38k15e/oAIp6rZd04PidJO3V4TzPA+61/ZveVzGaKsGvj0g6RtI1ZUtsi6T/KWnvts9dts7uLrtw/yBJ487xqfKzeyS9vpvvtf1r2xuA/wi8U9KLy3MtKM93n6Stkr4oaZ/ys+MlDUv6iKSHypbk28rPVgNvAz5ctty+1fZ1yyXdIunXkr4haWEX9WsB5wH7AM+X9DFJF0v6mqTHgHdJOkzSBkkPS9ok6T1lXU4FvgK8oqzLx8fqPsm/wTxJZ0n6maRtkiYMuNF8CX79ZRT4IHAw8ArgtcBfjNvnDcDLgaOBtwIntn3274G7yuM/CZw7PjhOxfZ1wDDwyrLoE8ALgeXAvwUWAX/Vdsi/Kb9rEfBOYJ2kI2yvAy4APlm2Zt/YdsxbgZXAUuAlwLs61ats2f058Dhwd1l8MnAxRavwAuDCsu6HAX8G/I2k19o+F3gfZSvS9jkdvu4DwJ8A/6E81yPAP3SqYzRPgl8fsX2j7R/bHrF9L/Aliv8J2621/ajt+4AfUASmMb+w/WXbo8D5wKHAc6ZZjc3AgWXQfA/wQdsP294O/A2watz+/8P2Dtv/AvxfiuA2lb+3vdn2w8C3xtV/vGMlPQr8CjgFeLPtX5efXWP7n8pW4cEU9/XOtP2k7ZspWnvv6O6Sn+K9wEdtD9veAXwM+LMuutbRMPkH6yOSXgh8BlgBPIPi3+/Gcbv9qu31E8B+E31m+4my0df+eTcWAQ8Dh5R1uLGt8ShgqG3fR8bdR/sFRWtpKuPrP9X+P7b9h5N8dn/b68OAsQDdXpcVHeoykecBl0pqtZWNUvwR+eUMzhezJC2//vIF4E5gme1nAh+hCDi1kPRyiuB3NfAQ8FvgKNsHlNvvlQMQY54lad+298+laDkCVL2cUPv5x1qr+4+ry0yC1f3A69uu+QDbC20n8PWZBL/+sj/wGPC4pN8H/msdXyrpmeVjJBcBX7N9a9md/DLwWUnPLvdbJOnEcYd/XNLekl5JcT/yf5flW4Hn11F/2/cDPwL+VtJCSS8BTqW4FzhdXwTWSHoegKRDJJ3cu9pGXRL8+st/B/4TsJ0i8Hyj4u/7lqTtFK2dj1J0ud/d9vmZwCbgx+Wo6pVA+/Nxv6IYENhMEWjeZ/vO8rNzgSPLket/qvQqCqcAh5d1uRQ4x/YVMzjP54ANwHfL382PKQaSos8oi5lGFSQdT9FKXDzLVYmYUFp+ETGQGhf8JK2UdFf5IOpZs12fiJibGtXtlTQE/CvFPMth4HrgFNt3zGrFImLOaVrL7xhgk+2f295JMbqYkbSI6LmmPeS8iKc+nDrMFCNpBx845MOXzK+8UjEz7vJRvvF7tR/ncfuN9VQMGJU/215btBBGjDLvd+9bLt63LEY9jxblTxeftVrCFjbQUnHSllD5BWqVPw202F1uUMtjlUAee+1yn/I9u8t/95q28qeUje03/nfjCV9WbTuPPGT7kD05x4mv3tfbHh7tuN+Nt+y43PbKPfmubjUt+E30wO5T/pnLSfGrAZ67aC+uu3xJHfWKGRh1a8Ly1rj/c1u0yv391PeYVlk2imkBO138HDXsQuzyPEYRT3qIXR5il/fiSc/nSc/nN60FPNmazxOtBTzR2pvtowt5fHQBvxlZwPaRBTy+awG/2bU3v901nyd2zmfHjvmM7ByitWMIds5j3pPzGNop5u0QQ0/C0E6YtxOGdpihHTC007/b5u0y83a2mLdzlHkjLbRrtNhGWjAyikZGoVW8ZnQUt1rQMrRGi5+jo0Vgb7WKwNcqfgdFWfn7avt9utUeCCf+PdOjW1pX+uJf7Ok5Hnp4lGsv7zzwP//Qnx28p9/VraZ1e4eB9mi2mN0zAgCwvc72CtsrDjloiIjoB2bUrY5bnZoW/K4HlklaWi7VtIrigdKI6GPF3QJ33OrUqG6v7RFJ7wcup5ggf57t22e5WhHRA2O3M5qiUcEPwPa3gW/Pdj0ioneM2VVzt7aTxgW/iJh7TDFo1SQJfhFRi7rv6XWS4BcRlTO7H2VqiqaN9kbEHNXqYuuGpNMl3SbpdklnlGUHSrqiTN51haRndTpPgl9EVM6Y0S62TsrMge+hmAp7NPAGScuAs4Dv2V4GfK98P6UEv4ionA27uti68CKK3C1P2B4B/gV4M8UaAOeX+5xPkWFvSgl+EVEDMdrFBhws6Ya2bfW4E90GvErSQZKeAZxEMSvsOba3AJQ/n92pRhnwiIjKmd1TlDt4yPakWfVsb5T0CeAKijzNPwVGZlKntPwiohZdtvw6sn2u7ZfZfhVFGtW7ga2SDgUofz7Q6TwJfhFRueIh594Ev7Zsgc8F/hS4kGINgHeWu7wT+D+dzpNub0RUzsAu96ytdYmkg4BdwGm2H5G0Flgv6VTgPuAtnU6S4BcRlRtbXLYn57JfOUHZNuC10zlPgl9E1KLl7rq1dUnwi4jKjd3za5IEv4ioQZEzpUkS/CKicsVKzgl+ETFgbLHTzcq5k+AXEbVo5Z5fRAyaYsAj3V4k3QtsB0aBEdsrJB0IfAM4HLgXeKvtR2ajfhHRa80b8JjN2rza9vK2SczTXo8rIvrD2IBHp61OTQrF016PKyL6x6jVcavTbN3zM/BdSQa+ZHsd49bjGpu8PF65vtdqgOcuyi3LiH5gxC436//X2arNcbY3lwHuCkl3dntgGSjXAaw4emGzMqJExIQy4FGyvbn8+YCkSynW498q6dCy1dfVelwR0R9M/d3aTmoPxZL2lbT/2GvgjyiWpp72elwR0T+aNuAxGy2/5wCXShr7/q/b/o6k65nmelwR0R9sGveoS+3Bz/bPKVLOjS+f9npcEdEfigGP3kxvk/RB4M8pbiXeCrwbeAbTfE64WaE4IuasUeZ13DqRtAj4ALDC9ouBIWAVydsbEU1kRMudty7tBewjaS+KFt9mkrc3IpqqFy0/278EPkUxLrAF+LXt7zKDvL0JfhFRuSJv77yOGx2Slkt6FkUrbylwGLCvpLfPpE7NeuQ6IuaorlNTTpm0HHgdcI/tBwEkfRP4A2bwnHCCX0RUrkhd2ZPR3vuAYyU9A/gtxRMiNwC/oXg+eC3J2xsRTWFrrFu7h+fxtZIuBm4CRoCfUEx33Y/k7Y2IJurVQ862zwHOGVe8g+TtjYimKdbza9bc3gS/iKhB81ZyTvCLiMoVj7qk5RcRA6aXc3t7JcEvImqRpOURMXCKJa3S7Y2IAZR7fhExcIpVXdLtjYgBU0xvS/CLiIHTvJZfZbWRdJ6kByTd1lZ2oKQrJN1d/nxW22dnS9ok6S5JJ1ZVr4iYHS3UcatTlaH4q8DKcWUTLjUt6UiKpaiPKo/5vKRmPRQUETM2NtrbaatTZcHP9g+Bh8cVT7bU9MnARbZ32L4H2ESRyzci5oguFzOtTd2d8MmWml4E3N+233BZ9jSSVo+t8vrgttFKKxsRvdHjHB490ZQ7kBNdtSfa0fY62ytsrzjkoPSMI/qBgRHP67jVqe7R3smWmh4GlrTtt5giI1NEzBEDM9o7iQ0US0zDU5ea3gCskrRA0lJgGXBdzXWLiKp00eXtptsr6QhJN7dtj0k6Y6onSSZT5aMuFwLXAEdIGi6Xl14LnCDpbuCE8j22bwfWA3cA3wFOs50behFzxNhipnv6qIvtu2wvt70c+HfAE8ClzCBpeWXdXtunTPLRhEtN214DrKmqPhExuyoY0Hgt8DPbv5B0MnB8WX4+cBVw5lQHZ4ZHRFRuGouZHizphrb362yvm2TfVcCF5eunPEkiqWPS8gS/iKicESOtru6ydcrbC4CkvYE3AWfPtE7NGn6JiDmrx9PbXg/cZHtr+X5r+QQJ3SYtT/CLiOqZXj/kfAq7u7ww+ZMkk0q3NyIq18sERpKeQfG0yHvbiteSpOUR0US9Cn62nwAOGle2jSQtj4imMWK0uwGP2iT4RUQt6l6vr5MEv4ionJ0ERhExoJzgFxGDp/71+jpJ8IuIWqTlFxEDx4bRVoJfRAygjPZGxMAx6fZGxEDKgEdEDChPmJJs9iT4RUQtmtbtrTKHx3mSHpB0W1vZxyT9si35yEltn50taZOkuySdWFW9IqJ+xWjvvI5bnar8tq8CKyco/+xYAhLb3waQdCTFktRHlcd8XlKS8kbMIXbnrU6VBT/bPwQe7nL3k4GLbO+wfQ+wCTimqrpFRP1sddzqNBtrzLxf0i1lt3gst+Yi4P62fYbLsqeRtFrSDZJueHBbsltG9APTOfB1G/wkHSDpYkl3Stoo6RWNyts7iS8ALwCWA1uAT5flE131hI1g2+tsr7C94pCD0jOO6BfuYuvS54Dv2P594GhgIzPI21tr8LO91fao7RbwZXZ3bYeBJW27LgY211m3iKiQwS113DqR9EzgVcC5ALZ32n6U4tbZ+eVu5wN/0ulctQa/sexKpTcDYyPBG4BVkhZIWgosA66rs24RUa0uu70Hj93WKrfV407zfOBB4B8l/UTSVyTty7i8vcDs5e2VdCFFBvWDJQ0D5wDHS1pO0cK9lzIBie3bJa0H7gBGgNNs54ZexBzS5Whup7y9ewEvA/7S9rWSPkcXXdzJTlQJ26dMUHzuFPuvAdZUVZ+ImD09nNs7DAzbvrZ8fzFF8Nsq6VDbW5K3NyKaw4DVeet0GvtXwP2SjiiLXkvRY0ze3ohoph4+xPyXwAWS9gZ+DryboiGXvL0R0TTdjeZ2w/bNwET3BZO3NyIaKKu6RMTAcfNWdUnwi4h6pOUXEYMpLb+IGESt2a7AUyX4RUT1xp7za5AEv4ioRXJ4RMRgSvCLiIGUbm9EDCKl5RcRA8eCHk1v65UEv4ioR1p+ETGQEvwiYiAl+EXEwGngQ86VreQsaYmkH5R5NW+XdHpZPml+TUlnS9ok6S5JJ1ZVt4ion9x56+o80r2SbpV0s6QbyrJG5e0dAT5k+0XAscBpko5kkvya5WergKOAlcDnJSUxb8Rc0cPEvcCrbS9vS3ZUTd5eSd+TdNK4snVTHWN7i+2bytfbKRILL2Ly/JonAxfZ3mH7HmATu/P6RkSf61XLbxKV5e1dCpwp6Zy2sqnSyz2FpMOBlwLXMnl+zUXA/W2HDZdl48+1eiyn54Pbkt0yom90l8CoU95eKNqI35V0Y9vnleXtfZRiffy/l/Qt4O1dHoek/YBLgDNsPyZNetNzog+e9rfA9jpgHcCKoxc2bPwoIibUfbe2U95egONsb5b0bOAKSXfOpErdtvxke8T2X1AEsqvpIrJKml/uf4Htb5bFW8u8mozLrzkMLGk7fDGwucv6RUTT9eien+3N5c8HgEspbo9NFlcm1W3w+2LbF38VeBfw3akOUNHEOxfYaPszbR9Nll9zA7BK0gJJS4FlwHVd1i8iGk6tzlvHc0j7Stp/7DXwR8BtVJW31/aXxr2/EfgvHQ47DngHcKukm8uyjwBrmSC/pu3bJa2nSEA8ApxmOzf1IuaK3tykeg5waXn7bC/g67a/I+l6mpK31/bVTL5o/4T5NW2vAdZUVaeImB09GM0FwPbPgaMnKN9G8vZGRCM1bIZHgl9E1KNhz2Yk+EVELbKYaUQMHnc3mlunBL+IqEdafhExkBL8ImIQNe2eX5VLWkVENFZafhFRj4a1/BL8IqJ6Ge2NiIGVll9EDBrRvAGPBL+IqEeCX0QMnB6t6tJLCX4RUY+GDXjkOb+IqEUvs7dJGpL0E0mXle8blbc3ImK33ubtPZ0iHe6YavL2RkTskW4CX5fBT9Ji4I+Br7QVV5a3d9okLZH0A0kbJd0u6fSy/GOSfinp5nI7qe2YsyVtknSXpBOrqltE1K/Lbm83eXv/DvgwT72LWFne3pkYAT5k+6Yy29KNkq4oP/us7U+17yzpSGAVcBRwGHClpBcmiVHEHNGDvL2S3gA8YPtGScfvSXWqTGC0BRiLxNslbQQWTXHIycBFtncA90jaRJGP85qq6hgR9enR9LbjgDeVPcaFwDMlfY0yb6/tLb3O27tHJB0OvBS4tix6v6RbJJ3XNiqzCLi/7bBhJgiWklaPNYkf3JZGYURf6NE9P9tn215s+3CKnuL3bb+dGeTtrTz4SdoPuAQ4w/ZjwBeAFwDLKVqGnx7bdYLDn/brsL3O9grbKw45aKiaSkdET6nLbQ+sBU6QdDdwQvl+SpU+5CxpPkXgu8D2NwFsb237/MvAZeXbYWBJ2+GLgc1V1i8iatTjGR62rwKuKl9PO29vlaO9As4FNtr+TFv5oW27vRm4rXy9AVglaYGkpcAy4Lqq6hcR9erlQ869UGXL7zjgHcCtkm4uyz4CnCJpOcXfgXuB9wLYvl3SeuAOipHi0zLSGzGHDMrcXttXM3E3/ttTHLMGWFNVnSJilmQx04gYWIPS8ouIaJclrSJiMCX4RcQgSssvIgaPadxipgl+EVG5JDCKiMGV4BcRg0huVvRL8IuI6k1/mfrKJfhFRC1yzy8iBlKmt0XEYGpYyy/Z2yKiel0sZ9VNt1jSQknXSfppmRjt42V58vZGREP1JnXlDuA1to+mWA1+paRjSd7eiGiisYec97Tl58Lj5dv55WaalLc3IqKdWu640UXeXklD5QLJDwBX2L6WhuXtjYgodN+tnTJvL0C5wvtySQcAl0p68UyqlJZfRNRCrc7bdNh+lCKB0UrKvL3wuzxBs5e3dyajMpLOlrRJ0l2STqyqbhExC3ow4CHpkLLFh6R9gNcBdzKDvL1VdnvHRmUeL1NYXi3pn4E/pRiVWSvpLIpRmTMlHUmRhPgo4DDgSkkvTBKjiLmhRzM8DgXOlzRE0Xhbb/sySdcA6yWdCtwHvKXTiapMYGRgslGZ48vy8ymarWeW5RfZ3gHcI2kTcAxwTVV1jIiaGOjBwga2bwFeOkF5c/L2wrRHZRYB97cdPlyWjT/n6rGRoAe3pVEY0S96fc9vT1Ua/GyP2l4OLAaO6TAqM1Gay6f9qbC9zvYK2ysOOWioRzWNiCr16jm/XqpltLfLUZlhYEnbYYuBzXXULyIqZne31ajK0d7pjspsAFZJWiBpKbAMuK6q+kVEvZrW8qtytHdaozK2b5e0HrgDGAFOy0hvxBzSsFVdqhztnfaojO01wJqq6hQRsyeLmUbE4DEw2qzol+AXEbVIyy8iBlOyt0XEIErLLyIGT1JXRsQgEqAMeETEIFLu+UXEwEm3NyIGU/1zdztJ8IuIWjRttDc5PCKiHj1Y1UXSEkk/kLSxTI9xelmepOUR0UAuRns7bV0YAT5k+0XAscBpZQqMJC2PiIbqQQIj21ts31S+3g5spFjxfdpJy3PPLyJq0eWjLgdLuqHt/Trb6yY8n3Q4xcpRT0uPISlJyyOiIboLfh2TlgNI2g+4BDjD9mPSRFkwppZub0RUz0Cri60LZSrcS4ALbH+zLG5O0vKIiDHCyJ23jucpmnjnAhttf6bto0YlLY+I2K3Vk9yUxwHvAG4t0+ICfARYS1OSlktaCPwQWFB+z8W2z5H0MeA9wIPlrh+x/e3ymLOBU4FR4AO2L6+qfhFRo7Fu756exr6aidPcwjSTllfZ8tsBvMb242Uf/WpJ/1x+9lnbn2rfuXxWZxVwFHAYcKWkFyaJUcTc0LSFDSq75+fC4+Xb+eU21dWfDFxke4fte4BNwDFV1S8iajYoeXsBJA2V/fIHgCtsX1t+9H5Jt0g6r20ayiLg/rbDh8uy8edcLekGSTc8uC2Nwoj+MEBJywFsj9peDiwGjpH0YuALwAuA5cAW4NPl7hP145/227C9zvYK2ysOOWioknpHRI+NZW/rtNWolkddbD8KXAWstL21DIot4Mvs7toOA0vaDlsMbK6jfhFRvV486tJLlQU/SYdIOqB8vQ/wOuDOsQcRS28GbitfbwBWSVogaSmwDLiuqvpFRM0a1u2tcrT3UOB8SUMUQXa97csk/S9JyykawvcC7wWwfbuk9cAdFCs3nJaR3og5wkCrWaO9lQU/27dQTDoeX/6OKY5ZA6ypqk4RMVuyknNEDKoEv4gYOAZGezK9rWcS/CKiBgYn+EXEIEq3NyIGziCN9kZEPEVafhExkBoW/LKSc0RUz4bR0c5bF8oFUR6QdFtbWfL2RkRD9W5621eBlePKkrc3IhqqR8HP9g+Bh8cVJ29vRDSRux3t7Tpv7zjJ2xsRDWRwdw85d5W3txcS/CKiHtVOb9sq6dCy1Ze8vRHREHaRurLTNnPTztub4BcR9ejRgIekC4FrgCMkDZe5etcCJ0i6GzihfD+ldHsjohbuTdJybJ8yyUeNydsbEVHKYqYRMYgauLBB5ff8yty9P5F0Wfl+0mkoks6WtEnSXZJOrLpuEVEPAx4d7bjVqY4Bj9OBjW3vJ5yGIulIYBVwFMXUlc+XyY8iot+5XMy001ajSoOfpMXAHwNfaSuebBrKycBFtnfYvgfYxO6cvhHR59xyx61OVbf8/g74MNAe0p8yDQUYm4ayCLi/bb/hsiwi5oKGtfwqG/CQ9AbgAds3Sjq+m0MmKHvanwJJq4HV5dvHhw7dtA14aKb1bLCDyXX1k7l8Xc/b05Ns55HLr/TFB3exa22/wypHe48D3iTpJGAh8ExJX2PyaSjDwJK24xcDm8eftJzk/LuJzpJuqGsuYJ1yXf1ljl/X4Xt6Htvjl6CadZV1e22fbXtx+YtbBXzf9tuZfBrKBmCVpAWSlgLLgOuqql9EDLbZeM5vLbC+nJJyH/AWANu3S1oP3AGMAKfZrnfsOyIGRi3Bz/ZVwFXl621MMg3F9hpgzTRP381aX/0o19Vfcl19Rm7YlJOIiDpkVZeIGEgJfhExkPo2+ElaWc4B3iSpY6amJplu6r1+mfMsaYmkH0jaKOl2SaeX5X19bZIWSrpO0k/L6/p4Wd7X1zVmYOff2+67DRgCfgY8H9gb+Clw5GzXaxr1fxXwMuC2trJPAmeVr88CPlG+PrK8vgXA0vK6h2b7Gia5rkOBl5Wv9wf+tax/X18bxQP4+5Wv5wPXAsf2+3W1Xd9/A74OXDZX/lvsZuvXlt8xwCbbP7e9E7iIYm5wX/D0Uu/1zZxn21ts31S+3k6xoMUi+vzaXHi8fDu/3EyfXxcM9vz7fg1+c3Ee8Jya8yzpcOClFK2kvr+2smt4M8WMpCtsz4nrYoDn3/dr8OtqHvAc0XfXKmk/4BLgDNuPTbXrBGWNvDbbo7aXU0y7PEbSi6fYvS+uq33+fbeHTFDWuOvqVr8Gv67mAfeZreVcZ2Yy57kpJM2nCHwX2P5mWTwnrg3A9qMUD+yvpP+va2z+/b0Ut45e0z7/Hvr2urrSr8HvemCZpKWS9qaYO7xhluu0p/p+zrMkAecCG21/pu2jvr42SYdIOqB8vQ/wOuBO+vy6POjz72d7xGWmG3ASxWjiz4CPznZ9pln3C4EtwC6Kv6anAgdRrGx9d/nzwLb9P1pe513A62e7/lNc1x9SdINuAW4ut5P6/dqAlwA/Ka/rNuCvyvK+vq5x13g8u0d758x1TbVleltEDKR+7fZGROyRBL+IGEgJfhExkBL8ImIgJfhFxEBK8IuIgZTgFxEDKcEvKifp5ZJuKdfF27dcE2+qubERlctDzlELSX9Nkb95H2DY9t/OcpViwCX4RS3KOdjXA08Cf+CkJY1Zlm5v1OVAYD+KFZ4XznJdItLyi3pI2kCxbNJS4FDb75/lKsWAqyVpeQw2Sf8ZGLH9dUlDwI8kvcb292e7bjG40vKLiIGUe34RMZAS/CJiICX4RcRASvCLiIGU4BcRAynBLyIGUoJfRAyk/w8E3+ZzW7cGkwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "plotDepthProfile(h, \"Tanh Depth Profile\")" ] }, { "cell_type": "code", "execution_count": 13, "id": "f9c3fc00-b277-401a-af3c-3d69d141535b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 20\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 10.70 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=10.634207999999974, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.06243400000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "# Compile the operator\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator\n", "op(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 14, "id": "b20891b6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling a tsunami in an ocean with 1D Tanh depth variation\")" ] }, { "cell_type": "markdown", "id": "94d6e979-c316-4f75-94c4-b1271ecf4770", "metadata": {}, "source": [ "## Example IV: Tsunami in an ocean with a seamount\n", "\n", "What happens if we have a constant bathymetry, except for a small seamount reaching 5 m below the water surface? Lets define the model parameters:" ] }, { "cell_type": "code", "execution_count": 15, "id": "cafa2271-d13a-4d83-8838-5137130b400d", "metadata": {}, "outputs": [], "source": [ "# Define constant ocean depth profile h = 50 m\n", "h0 = 50. * np.ones_like(X)\n", "\n", "# Adding seamount to seafloor topography\n", "h0 -= 45. * np.exp(-((X-50)**2/20)-((Y-50)**2/20))\n", "\n", "# Define initial eta [m]\n", "eta0 = 0.5 * np.exp(-((X-30)**2/5)-((Y-50)**2/5))\n", "\n", "# Define initial M and N\n", "M0 = 100. * eta0\n", "N0 = 0. * M0\n", "D0 = eta0 + h0\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 2.\n", "dt = 1/8000.\n", "nt = (int)(Tmax/dt)" ] }, { "cell_type": "markdown", "id": "e76012b0", "metadata": {}, "source": [ "Let's take a look into how this depth profile $h$ looks in the plot below." ] }, { "cell_type": "code", "execution_count": 16, "id": "52b08af0", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEWCAYAAAAQBZBVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqB0lEQVR4nO3de5BcZ3nn8e/v9HWkkSzJlo2RBHaI2cWmQGEVh42TLYebjZONIbWk5NoQs3HFZNfUklp2YztJcUmixbAQnGwWghyIFQh2lIsXxSGAYkIRAmtbBmN8xQo4WEiRLMm6jKSZ6T7n2T/et0c9o57pnlGfnu7p5+M6Nd2nzzn9nrH06H3Pe3lkZjjn3LBJFrsAzjm3GDz4OeeGkgc/59xQ8uDnnBtKHvycc0PJg59zbih58HNLiqQ7Jf3OIn33f5a0X9KYpHPjzx9a7HK51jz45UTST0j6mqSjkg5L+kdJP7rY5eoWSVdK2tPmmDslTUo6HrdHJb1f0jldKsPbJH31LM6/SJLFIDUm6RlJtyzwWiXgd4E3mNmomR2KP7+70PK5fHnwy4GklcC9wP8G1gDrgPcBE4tZrkXyQTNbAawF/hPwauAfJS1f3GJNs8rMRoHrgHdLunrmAZKKba5xAVAFHsuhfC4HHvzy8VIAM7vLzFIzO2VmXzSzRxoHSPolSU9Iel7SFyS9uOmz35P0rKRjkh6S9JNNn71X0p9L+nSsTX1b0ksl3SrpQDzvDU3Hv1DSjlj73C3pl5s+m9YUm1mbizWh/y7pkViD/TNJ1Ri4/hZ4YVOt6YVz/ULMbNzMHgR+FjiXEAg7+V2YpP8q6buSDkr6X5ISSS8D/hD4t/H7jzR93WpJfxN/P/dLekmb/1+NMn6dELxe3vhdSLpZ0r8AfyypIul2SXvjdnvc91LgqXiZI5K+1FT2H271XZJ+RtLDko7EFsIrOimj6x4Pfvn4DpBK2ibpjZJWN38o6U3ArwM/R6gR/QNwV9MhDwIbCbXGzwB/Lqna9Pm/Bz4FrAa+CXyB8P9yHfBbwMebjr0L2AO8EPgPwP+U9Np53MvPA1cDFwOvAN5mZieANwJ7Y9Nu1Mz2dnIxMzsO7AR+Ejr6XQC8GdgEvAq4FvglM3sC+BXg6/H7VzUdfx2hpr0a2A1saVcuBVcAlxF+pwAvIPw/eDFwI/AbhJrrRuCVwOXAb5rZd+J5EGqRr2nzXa8CPgm8nfAPwceBHZIq7crpuseDXw7M7BjwE4ABdwDPxdrXBfGQtwPvN7MnzKwO/E9gY6PGY2afjs+M6mb2YaAC/Kumr/gHM/tCPPfPCUHjNjOrAXcDF0laJWlDLMfNseb1MPBHwFvncTu/b2Z7zeww8NeEv/hnay8hqECb30X0ATM7bGbfB24nBLe5/JWZPRCv96cdlPkgcJjwu7nFzO6L+zPgPWY2YWangP8I/JaZHTCz5wgBdj6/y4ZfBj5uZvfHlsE2wiORVy/gWm6BPPjlJP5lfpuZrQdeTqh53R4/fjHwe7HJc4TwF0+EmhuS3hWbgUfj5+cA5zVdfn/T61PAQTNLm94DjMbvPBxrWw3/3PieDv1L0+uT8bpnax3hnqHN7yJ6tun1PxPuay7zLfN5ZrbazF5mZr/ftP85Mxtvev/C+P3zKUsrLwbe1bjneN8bFngtt0Ae/HrAzJ4E7iQEQQh/md9uZquathEz+1p8vnczobm5OjbnjhICwnztBdZIWtG070XAD+LrE8Cyps9eMI9rL2g5IEmjwOsIzVuY43fRdNqGptcvItzXgsswDzOvv5cQuFqVZT6eBbbMuOdlZjazue9y5MEvB5L+day9rY/vNxCaav8vHvKHwK2SLoufnyPpLfGzFUAdeA4oSno3sHIh5TCzZ4GvAe+PHRWvAG4gNAUBHgaukbRG0guAX53H5fcD56rDYSuxY+DfAP8XeB744/jRXL+Lhv8haXX8Pb4T+LOmMqyXVJ5Huc/GXcBvSlor6Tzg3cCnF3CdO4BfkfRj8Vnjckk/PeMfKZczD375OA78GHC/pBOEoPco8C4AM7sH+ABwt6Rj8bM3xnO/QOhJ/Q6hWTXO9GbffF0HXESoodxDeIa1M372KeBbwDPAFzkdVNqKtdm7gO/GpttsTbZfk3Sc0Jz9E+Ah4Mdjp0m730XDZ+N5DwN/A3wi7v8SoXf2XyQd7LTsZ+F3gF3AI8C3gW/EffNiZrsIz/3+gPAPwW7gbV0rpeuIfDFT188kGXCJme1e7LK4pcVrfs65odR3wU/S1ZKeUhiQu6CpRs65pUth8P234yDxXXHfGkk7JT0df65ue51+avZKKhCedb2eMDD3QeA6M3t8UQvmnOsbkp4BNpnZwaZ9HyQM67otVppWm9nNc12n32p+lwO7zey7ZjZJGLB77SKXyTnX/64FtsXX24A3tTuh3WTtXlvH9J7NPYRe05bOW1OwizaUci+Uc8PsoUcmDprZ2rO5xlU/tdwOHU7bHvfQIxOPEUY4NGw1s60zDjPgi7Ez7OPx8wvMbB+Ame2TdH677+q34NdqIO+0drmkGwnzLHnRuiIPfGFDi1Occ91SuHD3P7c/am4HD6fc/4X1bY8rXfhP42a2qc1hV5jZ3hjgdkp6ciFl6rdm7x6mj+Zfz4wR9Ga21cw2mdmmtecWelo459xCGallbbeOrhQX0TCzA4Sxq5cD+yVdCBB/Hmh3nX4Lfg8Cl0i6OI7a3wzsWOQyOefOkgEZ1nZrJ86GWdF4DbyBMDB+B3B9POx6wsD4OfVVs9fM6pLeQZjlUAA+aWa+OKRzS0BGZzW7Ni4A7pEEIX59xsw+L+lBYLukG4DvAzOnSJ6hr4IfgJl9DvjcYpfDOdc9hlHrsFk753VCWoBXtth/CJjPOpX9F/ycc0uPAWnui/DMjwc/51xPdPJMr5c8+DnncmdA2kezycCDn3OuR7rS3dFFHvycc7kzzJ/5OeeGjxnU+iv2efBzzvWCSBeUhiY/Hvycc7kzIPOan3NuGHnNzzk3dMIgZw9+zrkhY0DN+msdFQ9+zrncGSLts0WkPPg553oiM2/2OueGjD/zc84NKZH22TO//iqNc25JCis5J223TkkqSPqmpHvj+/dK+kHM5fuwpGvaXcNrfs653JmJSetqzp13Ak8AK5v2fcTMPtTpBbzm55zriQy13TohaT3w08AfnU15PPg553IXOjyStluHbgd+jTNXyXqHpEckfVLS6nYXWZTgJ+kZSd+ObfNdcd8aSTslPR1/ti28c25QhA6PdhtwnqRdTduN064i/QxwwMwemvEFHwNeAmwE9gEfbleixXzm91NmdrDp/S3AfWZ2m6Rb4vubF6dozrluanR4dOBgm6TlVwA/Gzs0qsBKSZ82s19oHCDpDuDedl/UT83ea4Ft8fU24E2LVxTnXLelprZbO2Z2q5mtN7OLCHm9v2Rmv9BIWB69mZDLd06LVfMz4IuSDPi4mW0FLjCzfQBmtk/S+a1OjNXgGwFetM47q50bBIaoWa5/Xz8oaSMhtjwDvL3dCYsVPa4ws70xwO2U9GSnJ8ZAuRVg0yurfbZCmHOulUaHR1evafZl4Mvx9Vvne/6iBD8z2xt/HpB0D3A5sF/ShbHWdyFwYDHK5pzrPqOzZm0v9fyZn6TlklY0XgNvILTPdwDXx8OuBz7b67I55/LTzRke3bAYNb8LgHskNb7/M2b2eUkPAtsl3QB8H3jLIpTNOZcDM/pubm/Pg5+ZfRd4ZYv9h4DX9ro8zrn8hQ6Prk5vO2veXeqc6wlfzNQ5N3QM+WKmzrnh5DU/59zQCXl7Pfg554aOfBl759zwCakrvbfXOTdkzOTNXufccBr6Qc7OueET1vPzZ37OuaHTf6krPfg553IXhrp4zc85N2T6cW5vf9VDnXNLVs5Jy+edAM2Dn3Mud2FJq7PP4dGkkbS8oZEA7RLgvvh+Th78nHM9kZnabp2YJWn5vBOg+TM/51zuwqouHdW1zmvk8o62xrw9zW4nJC1f0bSvowRozTz4OedyF6a3nX3e3uak5ZKuPJsyefBzzvVA16a3tUxazgISoOX2zE/SJyUdkPRo075Ze2Qk3Sppt6SnJF2VV7mcc4sjQ223dmZLWs4CEqDl2eFxJ3D1jH0te2QkXUq4kcviOR+V1F+DgpxzC5ZDb+9MtwGvl/Q08Pr4fk65NXvN7CuSLpqx+1rgyvh6GyHh8M1x/91mNgF8T9JuQi7fr+dVPudcb3V7VZcZScvnnQCt10NdpvXIAI0emXXAs03H7Yn7ziDpRkm7JO167lCaa2Gdc93RyOHRjaEu3dIv4/xa3bW1OtDMtprZJjPbtPZcbxk7NwgMqFvSduulXvf2ztYjswfY0HTcemBvj8vmnMtRvy1m2uvSzNYjswPYLKki6WLgEuCBHpfNOZeXDpq8vW725lbzk3QXoXPjPEl7gPcQemC2S7oB+D7wFgAze0zSduBxoA7cZGb+QM+5JWKoFjM1s+tm+ahlj4yZbQG25FUe59zi8vX8nHNDxxczdc4NJUPUs/7q8PDg55zriaF55uecc1PMm73OuSHkz/ycc0PLg59zbugYIvUOD+fcMPIOD+fc0LE+7PDor3qoc27JMlPbrR1JVUkPSPqWpMckvS/uf6+kH0h6OG7XtLuW1/yccz3QtYULJoDXmNmYpBLwVUl/Gz/7iJl9qNMLefBzzvVEJzW79tcwA8bi21LcWq792Y43e51zuTODNFPbjZi3t2m7cea1JBUkPUxYD3Snmd0fP3qHpEdi8rTVM8+byYOfc64nOszedrCxUnvcZiYsx8xSM9tIWPT4ckkvBz4GvATYCOwDPtyuPB78nHO5M7rT4THtmmZHCAmMrjaz/TEoZsAdhARoc/Lg55zrge6s5CxpraRV8fUI8DrgyZgWo+HNwKMtTp/GOzyccz1hC+qWOMOFwLaY1zsBtpvZvZI+JWkjoZL5DPD2dhfy4Oec64ku9fY+AvxIi/1vne+1cmv2xh6XA5Iebdo360BESbdK2i3pKUlX5VUu51zvhd7epO3WS3l+253A1S32f8TMNsbtcwCSLgU2A5fFcz4aq7XOuSXCrP3WS7kFPzP7CnC4w8OvBe42swkz+x6wmw56a5xzg6Pbvb1nazF6e1sNRFwHPNt0zJ647wySbmwMgHzukGe3dG4QGO0D31IPfrMNRGx11y0rwWa2tTEAcu253jJ2blBYB1sv9bS318z2N15LugO4N77dA2xoOnQ9sLeHRXPO5cnAsiFe0mqOgYg7gM2SKpIuBi4BHuhl2Zxz+eq3Zm9uNT9JdwFXEiYq7wHeA1zZaiCimT0maTvwOFAHbjIzf6C3BKSWzfl5QT7JaFj0uje3ndyCn5ld12L3J+Y4fguwJa/yuN5oF+zaHe/BcGlqzO3tJz7Dw3XFzCCWtXh8nRGOSVo8bUlin1fjOh4ElxgDPPi5paQ56DUHvEagS1u0dVJOP9EoqPEX4nSwS5AHwSVoaJq9bmlrFfSaA97U6zYDGBoreSTYGYHQg+BSor7r7fXg5+ZtZuCbGfRSjMxsKvBlU+eF940glwA1oIBIZGSm0CRW44xkWnPYA+CA85qfG2SNwNdc22sOejULoTA1I4X4+nQABMCMBCgoBMACFgIgUDojCE6vBXoAHFDmHR5ugLUKfDVLp2p6NYxJMzKgZjBpCRkijQtZprEWV8BIZBTMSDDKyigRAmJGRqmpJlhSqP9lmAfAQec1PzeImgNfq9pezYwaMG6iZgnjVqBmBSYJP1ObHrDKSkmUUSalpJSqUkrKyAjN5RICJWBp0/PAxAPgQDv7mp+kKvAVoEKIX39hZu+RtAb4M+Aiwhjinzez5+e6lgc/19bMGl8j8NUsm6rtTRiMWyFuRU5kFcatRM2KIQhagWyqCZtRVgh6JdWpqsbyZIKq6tRiIKzEIFsi1P6whER4DXCQzW8I6Gxmy9v7c8B9ZnabpFuAW4Cb57qQBz83L1n8rxH4xs0YN3EyKzJuRY5lVU5ahRNZ2E5mZSay0rTaX0EZJaVUkhrLkkmWJxOctAmWaYKVyThpUiMlDcFWgGWUFLJ/tRoj6AZAl8b5zZG391rCjDKAbYTERh783MLNbO42nvHNDHzHrMLxrMqRdDnH0ypH02UcT6ucSkucyspMpMVGakKKyiglKSPJJCOFGisK45xTOMmKwjg1CkwyzkpNQFKHGAATDCwNtUBv/g6kDsf5nSdpV9P7rTPTV8aFjh8Cfhj4P2Z2v6QLzGxf+B7bJ+n8dl/kwc+11fycr/GMb7Ip8B3JRjiWVTmcjnK4Psrh+nKO1JZxrF5hrFbhZL1MLSucHtMno1KoUy3UGC1NsLI4wfFSlTXFE6GGWEjIkgQ4BUmdBCMhAyUkZlO9wEkXniG5Huos+B00s01zXibM+98Ys7jdE/P2zpsHPzerM6esZaFXNz7ja9T4jmVVnquv5GB9BQdroxycHOX5iWUcnaxyYrLMeK1IvV4gi4Nck8QoFlOqpTrLy5OcUx7nRKXMRLlIikLzuAiJMhIzCtTjcBgjU0aBwrQyeu1vQHR5qIuZHZH0ZULqi/2SLoy1vguBA+3O9+Dn5nRGrY9Gr27o2Dgea3wH6ys4MLmC/RMree7UKM+PjzB2qsLkRIl0ogA1nf7DL4OScbKSMlapcWKkzKl66XTtsBKeCxZib3AhMUpKwxAZzGt/A0pdGOoiaS1Qi4Gvkbf3A4Rl8a4Hbos/P9vuWh78XEutan3Nzd1xK3AsPuM7XB/lYG2U/RMr2X9qBQfHljN2okr9RAmNJxTGE5IaKE7ptQJkJciqBU6NFKhNFqnVC1PPBBOFWl5JKaU4FKZkGQlGyWt/g8kE3ZneNlve3q8D2yXdAHwfeEu7C3nwc201an0ZxAHMCeNW5KRVOJ5WOVxfzsHJUZ47NcrBseUcHxshPVaieLxA8aQonILCJKgermdFSMuQjhSoL0uo1xOON83xLSUplaTOsmSC5ckEy22CqqVkSqk11f4Sr/QNli7U/ObI23sIeO18ruXBz82q0eRtPOtLzajF8XyNoSxH02UcqS3j+Ylloal7okp6rETpaIHSUVEag+JJozhuofMWyIpQr4r6MpFMCGUFasCYjFIxZaRYY3lhkhWFcVYUxlmmClXVqVoaan6x9hdqit70HRg+w8MNopQwV3fSEmpWYNxKnMhCze9YvcLRySpjpyrUT4QaX+moqBwxyseM8lhG8VRGUovDZkoJ9ZGEydEEpQkgLClQLxYZK1c4Wq6yshyufSKrsCIpxYHSCRVlcQaIGzge/Fy/a37el85YnSVDccpakZNZmVNpibFahROTZSYnwjO+4slQ4ysfM6rPp5SP1iicmES18NDPSgWKy8sktRDCrJCQVUQ6UmBypMSJyTJjtQqnyiVOZmVqVmSS8ExwaoWYGU1ff+7X5/pwMdPc/rRI2iDp7yU9IekxSe+M+9dI2inp6fhzddM5t0raLekpSVflVTY3f43VWVJELc7bncjCAOaT9TCcJZ0okIwnFE6Fpm55LKN8tEbxyCmS54+jQ0fQoSMkzx+neOQU5aO1UCs8aRTGIRlPSCcKjNeKjMfB0Y3vyiyJw2C6NUvK9Zqs/dZLef5TWQfeZWYvA14N3CTpUsKcu/vM7BLgvvie+Nlm4DLCuJ2Pxh4d1ycywuKjqSVMNgJgWqSWFajXw3CWpBY6N4rjRvFUFmp8YyexY2Nkx46THTuOHRtDYycpnJikeCqjOG4UJiCpATVRr5++bi0rMtmYG2zywDfI+ixxb0fBT9J9kq6ZsW/rbMdDmGJiZt+Ir48DTwDrCHPwtsXDtgFviq+vBe42swkz+x6wG7i8w/twOciY3vwFppalykhI45JVmSkMYDahNPTqJnVIallo6k5MYhMTZONhs4kJmAjN4KSWkdTDMBilQLxWZqLeuH78Y9r47ual8TMPhwNjUGt+FwM3S3pP0745p6A0k3QRoXv6fmDaHDygMQdvHfBs02l74r6Z17pR0i5Ju5475NktnRsYpvZbD3Ua/I4QxtBcIOmvJZ3T6RdIGgX+EvhVMzs216Et9p3xb4GZbTWzTWa2ae253irOU/MKKo2l5wvxf0lCmIFRVEYiI0nCP91WCOP4smLo1bVSASplVKmQVMOmSgUqZbJykayUkBXDwGcrAEm4VjGJ127M62367tO5PlpngnN9qJMmb49rfp329srM6sB/kfQ24KvA6rlPgbje1l8Cf2pmfxV3zzYHbw+woen09cDeDsvneiAhzL4oKz29Hl+SUinUKRZTKBlZKQxgrldFfSShuLyMassQoIlyuFCljI0uI11eoj6SUK+KtBJmfVAM835LSWOrT31fIvNQN8j6bKhLp3+W/rDxwszuBN4GfHGuEySJkKT8CTP73aaPGnPwYPocvB3AZkkVSRcDlwAPdFg+l7NGzo2wBP3p9fhGkkmqhRrVUp1CJSWrZqQjUF8mJkcTJs8pUV81QrZ6BXbuKuzcVaRrVlJbPUJtZYnJ0YT6MpFWIatmFCphwYNqIVy7FINsoiwsbqB8e+lcfpS133qpo5qfmX18xvuHgF9qc9oVwFuBb0t6OO77dcLE4zPm4JnZY5K2A48TeopvikvXuB4rKGlKGRk6H2rEmh8Wl56vsyyuxzdammB5eZKxSo1TccpaMqE4gBmykiguL7Yc5DxxTkJtBdSXGTaSUq7UWF6eZLQ0wUghLHZaUp0yaWwCxzKiac1fH+M3APqs5pfbIGcz+yqzL9rfcg6emW0BtuRVJrdwBUShkWxI6dTS8ysK46wsToRlqUbK1CaL1OsJygqAsEJo1hbHk5bT22orYHKlUV+ZUlxWZ3QkXGtlMVw7LG9fo6SUsjIKsSxusCxGb247PsPDzSqZmjsbs60hShhVpVNLz59TOMnxUpUTlbgsVb3A8VhTtKRAVhH15aIwoWmruqQVSKuhxldfmVJYUWN0+Tirq6dYXTnJqtJJzimcDAsbJBMxwVGoiSYKy9knPq93sPTZDA8Pfq6tRtM3NHuhpIyq6ixTqJ2tKZ5gohwGJTeWpRqTUS8WSUfCrI/WS1pl2Eis8S0f57zRE6wdGeO88hhriidYVTjJMoXERiWF0X4ldEaT1w0Ir/m5QdD83A9C/a+kxiovRk0pK5OQc6NmhancvOHYsDrLWLnC5EhYzDRtXsw0MSgahUp4xjc6MsHq6inWjoxxQeUY55ePsaY4xorCKVYm41SVUpVRbqr1zSyr63/e7HUDpdH0RVlIIo5IY9M3TWpMMk5aCLM9qIShMKUkLEt1tNz5MvbnVk6wpnyC88vHWFs8zprCGCuScZYltdDkZXqtz5u8A8a605sraQPwJ8ALCDMut5rZ70l6L/DLwHPx0F83s8/NdS0Pfm5WrWp/iUJC8QpGSspKTYRkQ8W49DxGJamzvDDJynJIYDSelkL2tljzKyYhe1tzAqNVpZOsKZ5gTXGMNYUxViYhg1tVKRVByWt9g687Nb/GmgHfkLQCeEjSzvjZR8zsQ51eyIOfa2tm7Q8lU83f0IV7KozDiz3By2Iv8PG0yqlyXJ0lK1CPeXtnS125qnCSFYVTrIiBb1lSn2rulpR4rW/QdWcl531AY3rscUmNNQPmzYOfm1Oj9tcIgCUR8ueiOJApBMDE4vg/0qkhMI2k5TULq7NM5ejAKCX1GCgnp3p0G0nLG03dqoyqNNXcLakwLfB5rW+wdPjMr23e3qnrTV8z4ArgHZJ+EdhFqB0+P9cXefBz85LEoS8hCGZTCcUL1KeyrC23CZapsQJzMa7CXJhanSUhOz09TvWpMYNV1WPQC03dcgx8JZ0e2OKWvLZ5e+HMNQMkfQz4bUL98reBD9NmIoYHP9dWc+0vw8IwEwu1wMbCA1PZ1iyjailV1UPQIyxEOmnTF6Eoxylr5ZidLYzjy6gqLFFfamrqhp7mcL7X+gZYl3p7W60ZYGb7mz6/A7i33XU8+LmOzGz+JgpL2mMpKJlKKp5gZEqpWspkXI+vMQymsR5fYWrQdDg+zBoJYwgbtb1E05u64IFvoHWvt7flmgGNxVLi2zcDj7a7lgc/17GZNcAkPgNMLGRTa+TVrcWfFWVTy9/PTNnaWCghgTBlbUbQS0imOjfC8R74Bl53an6zrRlwnaSN8VueAd7e7kIe/Ny8TK8BQqMXuEBhWhDMaEp8NJVkaPp6fM2LFMwMeuHz6b26HvgGl+jOIOc51gyYc0xfKx783Lw1j/+bGgYDZwTBEkwFwPB5i2txOtA1B73T1z79nW7A+QwPtxQ0gtEZtUA4/TyQ0CSe+zrTA1547UFvyfFVXdxS0zoIQnMgbGgkG2o1ZGXmoGUPektQn+Wa8uDnuqI5CMKZwSzsa59zxYPe0uU1P7ekzQxezXODOzneLWEe/Nww8eDmgEXJztZObn8yJW2Q9PeSnpD0mKR3xv3vlfQDSQ/H7Zqmc26VtFvSU5Kuyqtszrne67ek5XnW/Oa19IykS4HNwGXAC4G/k/RST2Lk3BIxLDU/M9tnZt+Ir48D7ZaeuRa428wmzOx7wG7g8rzK55zrrX5LXdmTBzIzlp6BsPTMI5I+KamR/Hwd8GzTaXtoESwl3Shpl6Rdzx3ySqFzA8E63Hoo9+A3c+kZ4GPAS4CNhEUJP9w4tMXpZ/w6zGyrmW0ys01rz20/dMI5t/jU4dZLuQa/2ZaeMbPUzDLgDk43bfcAG5pOXw/szbN8zrkeGpaa31xLzzQd1rz0zA5gs6SKpIuBS4AH8iqfc663hqm3d15Lz5jZY5K2A48Teopv8p5e55aQPuvtzS34LWTpGTPbAmzJq0zOuUXSpcVMu8mH3zvneqMLz/zmmDyxRtJOSU/Hn6vbXcuDn3OuJ7r0zK8xeeJlwKuBm+IEiVuA+8zsEuC++H5OHvycc73RhZrfHJMnrgW2xcO2AW9qdy1f2MA51xM55+29oJHAyMz2STq/3Rd58HPO5c/odDHThebtnXeRvNnrnMtdI4FRN8b5tZo8AexvjCGOPw+0u44HP+dcb3Snt7fl5AnCJInr4+vrgc+2u5Y3e51zPaE2yaw6NNvkiduA7ZJuAL4PvKXdhTz4Oefy16W5u3NMngB47Xyu5cHPOdcTnsDIOTeU+m16mwc/51xveM3POTd0FmHJqnY8+DnnesODn3Nu2DQGOfcTD37OuZ5Q1l/Rz4Ofcy5/i5Cjox0Pfs65nui3oS55JjCqSnpA0rfiiqvvi/tnXXFV0q2Sdkt6StJVeZXNObcIhiV7GzABvMbMXknI0Xu1pFczy4qrcTXWzcBlwNXARyV5Yl7nloh+y96WW/CzYCy+LcXNmH3F1WuBu81swsy+B+zmdE5f59wgM8Cs/dZDeSctL8SVFw4AO83sjBVXgcaKq+uAZ5tO3xP3zbzmjZJ2Sdr13CHPbOncoFDWfuulXIOfmaVmthFYD1wu6eVzHN5qpYYz/ikws61mtsnMNq0911vFzg2Cbi5m2i09WczUzI4AXyY8y5ttxdU9wIam09YDe3tRPudczjpp8i6VZq+ktZJWxdcjwOuAJ5l9xdUdwGZJFUkXA5cAD+RVPudcb3VxGftPSjog6dGmfe+V9ANJD8ftmnbXyXOc34XAtthjmwDbzexeSV+nxYqrZvaYpO3A44TcnDeZmT/Uc26p6F7F7k7gD4A/mbH/I2b2oU4vklvwM7NHCGnlZu4/xCwrrprZFmBLXmVyzi2ebj3TM7OvxLSVZ8UTGDnn8mdAau23mLe3abtxHt/yDkmPxGbx6nYHe/BzzvVEh8/8DjZGc8StZcLyFj4GvIQwoWIf8OF2J/jcXudcb+TYm2tm+xuvJd0B3NvuHK/5Oed6Is9xfo3hc9GbgUdnO7bBa37Oufx1ceECSXcBVxKeD+4B3gNcKWlj/JZngLe3u44HP+dc7gQo7U70M7PrWuz+xHyv48HPOdcT6vEMjnY8+Dnn8ucrOTvnhlPv5+6248HPOdcTnr3NOTecvObnnBs61r3e3m7x4Oec643+in0e/JxzveFDXZxzw8mDn3Nu6BjQZ0nLPfg553InzJu9zrkhlfVX1S/PBEZVSQ9I+pakxyS9L+6fNdGIpFsl7Zb0lKSr8iqbc67HGs3edlsP5VnzmwBeY2ZjkkrAVyX9bfzsjEQjki4FNgOXAS8E/k7SSz2JkXNLQ781e3Or+VkwFt+W4jbX3V8L3G1mE2b2PWA3cHle5XPO9diw5O0FkFSQ9DAhMflOM7s/ftQq0cg64Nmm0/fEfTOveWMjuclzh7xS6Nxg6F7S8lny9q6RtFPS0/Hn4iYwMrPUzDYC64HLJb2c2RONqNUlWlxzayO5ydpzC7mU2znXZZ1nb+vEncDVM/bdAtxnZpcA98X3c+pJDg8zOwJ8GbjazPbHoJgBd3C6absH2NB02npgby/K55zLn8zabp0ws68Ah2fsvhbYFl9vA97U7jp59vaulbQqvh4BXgc8OUeikR3AZkkVSRcDlwAP5FU+51yP5fvM7wIz2xe+xvYB57c7Ic/e3guBbZIKhCC73czulfSpVolGzOwxSduBx4E6cJP39Dq3RBiQdRTczpO0q+n91nnk7p2X3IKfmT0C/EiL/W+d45wtwJa8yuScWywd1+wOmtmmBXzBfkkXmtm+2Lo80O4Ez9vrnOuNfJu9O4Dr4+vrgc+2O8Gntznn8mdA2p0pHLPk7b0N2C7pBuD7wFvaXceDn3OuBwysO8Fvlry9AK+dz3U8+DnneqPPprd58HPO5a/z3t6e8eDnnOsNr/k554aSBz/n3NAxg7S/5ix48HPO9YbX/JxzQ8mDn3Nu+Jj39jrnhpCBdWmQc7d48HPO9UaXprd1iwc/51z+zPoudaUHP+dcb3iHh3NuGJnX/Jxzw6f3qSnb8eDnnMtfHy5skPtKzjF37zcl3Rvfz5pfU9KtknZLekrSVXmXzTnXGwZYmrbdeqkXy9i/E3ii6X3L/JqSLgU2A5cRcnJ+NCY/cs4NOouLmbbbOiDpGUnflvTwjGRH85Jr8JO0Hvhp4I+ads+WX/Na4G4zmzCz7wG7OZ3T1zk34Cyztts8/JSZbVxgsiMg/5rf7cCvAc0hfbb8muuAZ5uO2xP3OeeWgi7V/Loltw4PST8DHDCzhyRd2ckpLfad8U+BpBuBG+PbscKFuw8BBxdazj52Hn5fg2Qp39eLz/Yix3n+C39nf3FeB4dWO8jba8AXJRnw8YXm9c2zt/cK4GclXQNUgZWSPs3s+TX3ABuazl8P7J150XijUzcradfZVH37ld/XYFni93XR2V7HzK7uQnEarjCzvZLOB3ZKetLMvjLfi+TW7DWzW81sffzFbQa+ZGa/wOz5NXcAmyVVJF0MXAI8kFf5nHODycz2xp8HgHtYYN/AYiQtvw14vaSngdfH95jZY8B24HHg88BNZtZfS7865xaVpOWSVjReA28AHl3ItXoyyNnMvgx8Ob4+xCz5Nc1sC7BlnpdfUHt/APh9DRa/r964ALhHEoT49Rkz+/xCLiTrsyknzjnXC4vR7HXOuUXnwc85N5QGNvhJujrOAd4t6ZbFLs98SPqkpAOSHm3aN/BzniVtkPT3kp6Q9Jikd8b9A31vkqqSHpD0rXhf74v7B/q+GoZ2/r2ZDdwGFIB/An4IKAPfAi5d7HLNo/z/DngV8GjTvg8Ct8TXtwAfiK8vjfdXAS6O911Y7HuY5b4uBF4VX68AvhPLP9D3RhiAPxpfl4D7gVcP+n013d9/Az4D3LtU/ix2sg1qze9yYLeZfdfMJoG7CXODB4KFAZmHZ+we+DnPZrbPzL4RXx8nLGixjgG/NwvG4ttS3IwBvy8Y7vn3gxr8luI84CU151nSRcCPEGpJA39vsWn4MGFG0k4zWxL3xRDPvx/U4NfRPOAlYuDuVdIo8JfAr5rZsbkObbGvL+/NzFIz20iYdnm5pJfPcfhA3Ffz/PtOT2mxr+/uq1ODGvw6mgc8YPbHuc4sZM5zv5BUIgS+PzWzv4q7l8S9AZjZEcKA/asZ/PtqzL9/hvDo6DXN8+9hYO+rI4Ma/B4ELpF0saQyYe7wjkUu09ka+DnPCsPuPwE8YWa/2/TRQN+bpLWSVsXXI8DrgCcZ8PuyYZ9/v9g9LgvdgGsIvYn/BPzGYpdnnmW/C9gH1Aj/mt4AnEtY2frp+HNN0/G/Ee/zKeCNi13+Oe7rJwjNoEeAh+N2zaDfG/AK4Jvxvh4F3h33D/R9zbjHKznd27tk7muuzae3OeeG0qA2e51z7qx48HPODSUPfs65oeTBzzk3lDz4OeeGkgc/59xQ8uDnnBtKHvxc7iT9qKRH4rp4y+OaeHPNjXUudz7I2fWEpN8h5G8eAfaY2fsXuUhuyHnwcz0R52A/CIwDP26eltQtMm/2ul5ZA4wSVniuLnJZnPOan+sNSTsIyyZdDFxoZu9Y5CK5IdeTpOVuuEn6RaBuZp+RVAC+Juk1ZvalxS6bG15e83PODSV/5uecG0oe/JxzQ8mDn3NuKHnwc84NJQ9+zrmh5MHPOTeUPPg554bS/wf5lJYd6GR/fQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "plotDepthProfile(h, \"Seamount Depth Profile\")" ] }, { "cell_type": "code", "execution_count": 17, "id": "149fb8d8-5c92-4730-a39d-74f8ce340242", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 40\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 21.54 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=21.467746999999985, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.06220900000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "# Compile the operator\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator\n", "op(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 18, "id": "3cce7e29", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling a tsunami in an ocean with a seamount\")" ] }, { "cell_type": "markdown", "id": "b08916a5-9a5e-4231-9313-cb5d7a0ccd1e", "metadata": {}, "source": [ "## Example V: Tsunami in an ocean with random seafloor topography variations\n", "\n", "Another modelling idea: what is the influence of the roughness of the seafloor topography. First, we add some random perturbations on the constant bathymetry model $h_0$ by using the `random.rand` function from the `NumPy` library. To smooth the random perturbations, we apply the `gaussian_filter` from the `SciPy`library:" ] }, { "cell_type": "code", "execution_count": 19, "id": "97a8e17b-eb3d-4e74-b962-b9eaa3af98c2", "metadata": {}, "outputs": [], "source": [ "from scipy.ndimage import gaussian_filter\n", "\n", "# Define constant ocean depth profile h = 30 m\n", "h0 = 30. * np.ones_like(X)\n", "\n", "# Add random seafloor perturbation of +- 5m\n", "pert = 5. # perturbation amplitude\n", "\n", "np.random.seed(102034)\n", "r = 2.0 * (np.random.rand(ny, nx) - 0.5) * pert # create random number perturbations\n", "r = gaussian_filter(r, sigma=16) # smooth random number perturbation\n", "h0 = h0 * (1 + r) # add perturbations to constant seafloor\n", "\n", "# Define initial eta [m]\n", "eta0 = 0.2 * np.exp(-((X-30)**2/5)-((Y-50)**2/5))\n", "\n", "# Define initial M and N\n", "M0 = 100. * eta0\n", "N0 = 0. * M0\n", "D0 = eta0 + h0\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 3.\n", "dt = 1/4000.\n", "nt = (int)(Tmax/dt)" ] }, { "cell_type": "markdown", "id": "d6ef2c5a", "metadata": {}, "source": [ "Let's take a look into how this depth profile $h$ looks in the plot below." ] }, { "cell_type": "code", "execution_count": 20, "id": "bf708c5a", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAEWCAYAAAAQBZBVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAADQpElEQVR4nOz9e7Rty1kfBv6+qrn2OVcPELIEXCTZsi1BjOi2GFbUTuO0MeCA07YVZxgb2jYQZMvdhqbBpI1EOgbsKIE2iJCAsS/BjXDAoBgwsgY0r7RClOhhxMCAAMeyEVhICASS9bj37L1m1dd/fI/6qmattfe5Z517z+HuGmPttfZ81KxZj1/9vkd9RcyM63SdrtN1eqKl9HgX4Dpdp+t0nR6PdA1+1+k6XacnZLoGv+t0na7TEzJdg991uk7X6QmZrsHvOl2n6/SETNfgd52u03V6QqZr8HuMExF9NRH99493OR5tIqJPIaJ/SUQfIqL/iIheT0R/+fEu191KRMRE9LzH4bkPENE/JaJ/S0T/AxH9BSL60ce7XL+T0jX4ASCidxDRIzqgf52IvoOInvJ4l+seTX8LwDcz81OY+Z88lg/WAf9hbaffIqKfIKI/f8L87wjIdWLba/neT0T/KxH9e48yuz8L4GMA/C5m/mxm/i5m/g8ebdmu0zZdg19Lf4qZnwLghQA+GcArHt/i3LPp9wB42918ABEtR07/QW2nTwDwHQC+mYi+6m6W5zbT92r5ngngDQC+n4hovIiI8iX5/B4A/xszr3ehjNcJ1+C3Scz86wB+BAKCAAAiejkR/Ssi+iAR/QIR/Zlw7guI6A1E9PVE9D4i+mUi+hPh/O8lov9J7/0xAM+IzyOiP01Eb1Om8Hoi+gPh3DuI6P9JRD+rjOfbiehjiOiHNb8fJ6KPmr0HET2DiF6n+f42Ef3PRJT03McR0fcR0W9qeb8k3PdiInqj3vduIvpmIjrTc/8KwO8D8E+V3dwYnpmI6P9FRL9CRL9BRN9JRB95G+/6FUT0swA+fAkAgpnfy8z/EMD/DcAriOh3aT4fqfX0biL6NSL6LwxotK3+FyL6b1Wc/CUi+nQ990oA/z4ETD9ERN8cHvcZKuq/j4i+ZQZmk/LtAbwawMcC+F0qTXwrEf0QEX0YwB8joj+g9fB+rZc/rWX5GgB/E8Cf17K81PrZ7FlEdEP7368S0XuI6O8R0QOXlfEJn5j5Cf8B8A4An6G/nw3g5wB8Uzj/2QA+DjJZ/HkAHwbwoJ77AgB7AH8FQIYMxncBID3/RgCvAnADwP8JwAcB/Pd67uM1rz8OYAfgbwB4O4CzUK43QcSfZwH4DQA/DWGmNwD8jwC+6sA7/VcA/p7mu4MMbNJ3eCtkcJ1BwOxfA/hMve8PAfjDABYAzwXwiwC+dFZX+v/rAfxl/f2FWv7fB+ApAL4fwD+8jXf9GQDPAfDAgXdiAM8bju0ArAD+hP7/TwD8fQBPBvDRAN4C4K+GtloBfJne9+cB/FsATx/fZXjm6wA8DcDvBvCbAD7rQPm+OrTtDQB/B8C/0f+/Q5/1KdoGT9X3/0pth0+D9I1PGPMKZX/DrC4A/NcAXgvg6ZrvPwXwXz3e4+pe/zzuBbgXPjrwPqSdjwH8BICnHbn+ZwC8RH9/AYC3h3NP0jw+VgfLCuDJ4fx3hwHynwN4TTiXAPwagE8N5foL4fz3AfjW8P//HcA/OVDGvwXgBydg8X8A8KvDsVcA+P8cyOdLAfzAUFeHwO8nAPy1cO4TIBPDcsV3/cJL2mkDfnr81wH8BcgkcY4AngA+F8D/L7SVT0x67C0A/tL4LsMz/0j4/zUAXn6gfF8N4ALA+yET1f8I4A/pue8A8J3h2n9fy53CsX8E4KtDXpeCH2RC+zCA3x/O/XsAfvnxHlf3+ueoaPEES/8RM/84Ef1RCEA9A9KJQUSfB+CvQ5gQIKwmiq+/bj+Y+WGViuya9zHzh8O1vwJhN4CwyV8J91Yi+jcQlmfpPeH3I5P/Dxlm/g5kAP2oluchZv5aiC7p44jo/eHaDOB/1nf9eAhTfREEyBcIU7xK6t5Hfy8QULrKu/6bKz7HExHtIPq134a82w7Au4NkmoZ8f40VIUIZP+6Sx/x6+P0wDtc5IAD/Fw+ci+X4OAgrrENZnoXbS8+EtNNbwzsTpE2v05F0rfMbEjP/T5BZ+usBgIh+D4BvA/DFEMvb0wD8PKSDXZbeDeCjiOjJ4djvDr/fBRmw0GcRBBh/7dG/gSRm/iAzfzkz/z4AfwrAX1f91r+BsIKnhc9Tmfk/1Fu/FcAvAXg+M38ERCy7yrtu3geN+b5nPHfgXR9NiKGX6DPeAnm3cwDPCO/2Ecz8gnD9swad3e/Wsj3a599Oivm/C8BzTA8bynK7bf9eyCT4gvDOH8lidLlOR9I1+M3Tfw3gjxPRCyG6I4boekBE/wmAT7pKJsz8KwB+CsDXENEZEf0RCBBZeg2A/zMRfboymC+HDN7/9U5fgIj+JBE9Twf6BwAU/bwFwAfUuPAAEWUi+iQi+nf11qfq9R8ion8HosO8avpHAL6MxMjzFAD/JcT6uZ76XYno6UT0FwB8C4CvY+bfYuZ3A/hRAN9ARB+hBpjfr2ze0kcD+BIi2hHRZwP4AwB+SM+9B6KvfCzSmyHi6t/QsnwqpG98z+1koszx2wB8IxF9NAAQ0bOI6DNPW9zfeeka/CaJmX8TwHcC+M+Z+RcAfAPEcPEeAP87AP/LbWT3f4Ho2X4bwFdpvvacfwHgLwL4byEz+J+CuNxcnOA1ng/gxyG6zDcC+LvM/HpmLvqcFwL4ZX3ufwfArLL/qZb5g5BB9b238cx/AOAfAvhJzfsWRC95ynf950T0IYix4C8D+DJm/pvh/OdBDAi/AOB9AP4xgAfD+TdD6ua9AF4J4M8y82/puW8C8GfVqvvf3Ga5bivpe/9pAH9Cy/J3AXweM//So8juKyD18SYi+gCk3T/hVGX9nZrMInmdrtPv+EREXwAxaPyRx7ss1+nxT9fM7zpdp+v0hEz3HPgR0WcR0b8gorcT0csf7/Jcp+t0nX5npntK7FVP/P8N4gj7TgD/DMDnqt7tOl2n63SdTpbuNeb3YojD8L9WhfD3QFwZrtN1uk7X6aTpXnNyfhZ6R9B3Qiyl07TcfDKfPfXp4oVGANt30t8JQGKAAEoMIkai+A0k6G/ocdjxKtnaMbA6uxlTJjCACkJlubpw0g9hrUmOVwJXAioAJlCVLOwbAIi7bOUpk3dC0osTIyVGThVLqlhIPtk+qF25reRSVkJFQmVCsW8tf+X2YRCYAWZ9a/+OLUCt3PE4T85vzvXZ9L+5HSOrIID0N+lvazP7nbRdU/idUZGoIoP1u3ZtLtWqNUS0cWhkZljzVFD7aP0VTv69WvtX6QOlJjCH9q8kr1Kl/f23PqB995UkVU+tn1vVEjbH4gvwVb0zWxVP061ff+d7mfmZV89tmz7zjz2Zf+u3y6XXvfVnz3+EmT/rTp511XSvgd+sufrhRvQyAC8DgN1TPgr/zn/8ZagZqDtCPQPKGVBuAuUGozzA4BsVdLNgubHi7MaKG7s9bu5W3FxWPLDscZZW3MwrbuQVZ2nFjfDZUcHNtMeNtMeOig6kVpw9Z5zXHR6uZ/hguYkPrjfx/osH8P6LB/CB85v44K0beOTWDuutHfhWRnokIZ8T8i1COgfSCqQ9QCWAoYH3AtQdUHZAucmoN4DyQAU/UJCftOKBBy7wEQ/cwkfdfARPv/FhPP3sYXzU8jA+cnkYT0238OR0LmUmWUBwwRl7XnCr7vDBehMPlxv4kH4+vN7Ah8sZHl7P8Mi6w/m64KJk7EvCvmSUklCrAHktScDAQV1AEvob9psBmg36CJYR6BPAiYGs34mBLJ+0VFBmpFSRMyPnil0uWHLF2bJilypu5BU3lhU38x5PXi7w5HyBJy/neFK6wFPzLTwl38JHpEfwpHSOJ6dzPJkucJNW3KCCM6rYkcwtGUBWH+jCjAJgz8AehFuccYszHq438GE+wwfLA3h/eRL+bXkS/u36AN6/Pgm/ffEkvP/8AXzw4gY+FNv/EW3/W4TlEUK+BfmcM/I5kC8Yec+gAqSVHQRZkB11IXAGaibURfv7AtQF4IVQs/QZTgDnYdKcjazQDl2bTNIv/Zd//VfmZ66e3vvbBW/+kWdfet3uwX/1jEsvOlG618DvnWhLvwAJMvCueAEzPwTgIQB48jOfw9Z4kT2REgeqcogZwmgqgQO7WWvCQjJjp5qRwMjESCwfAMhcgQpUSthRQYKASeWEPWfcqjvcqjuc1wXnZZHZX2f+ygSuKTC/8EH7HWd8hgAFh3MCIuwAw1XeZV8y9iXjoi64qIuUoe6wowJU4GbaS/kRwI932HP2Tx1Y35icZRnDShWoyRkZJykfG0M7RDdGZsgKlMrk2I5VbTw7r+yT9Lysy9T29G/q3qEwYc9J2Hcif9dbvEPiiszCBC1VEAqUHcIKY+eAC07YI2HPCbd453V5wXlgf30dygpcwjG1+iHGNbI/Ynb23epQq7wClKSwTDqRJq3egsYSjz37rqv+GaVbyff4p3sN/P4ZgOcT0e+FLPP5HIjD7Typ+Mja8KQf1Pa/MQ8DjLUmLFUBKhnwJRGTFAQ9ab6FEgol7Ckj0Rb8Hintc6vscFEzVmVMzKRloVa+ANidqMMivjG3a6k2ZkgF4ELgNWFdM/ZrxnlZ8PB6hrO0YqHizHSfZMDb/5VJwU8A8rwuylwX7GtWMX0LgIkY1YCPyEVPPQlUDuKXiWaB1oXk71p72YzJlAgQ8CeAiaXtkrYfCQtibhMYq2iZiTuRfa0ZK1Xsk7RPAjfmrpMXkkwIF2mPMy56njtmb/VWIAB6AcnvYb6BD9cbuMVnOK8ymaw1Ya25qRG0fOP7dyqOSZqCXliNR4yub6DKxIRC9lpt4hxE48czCfG/d4yrwD0Gfsy8EtEXQ+LpZQD/gJmPBs5sHYE7gCETtwoJYCj41ZqwloR9SkhF2N6YDAAqE1ZK2KfsIm+m6p3bGMUjZYcPrzdcbLy1LjhfM9Y1oa4kHbM2MXCj94u/XedFAnpZwa8AaSXwSuB9QtllXKwZD+93ovujgmzAxxlPShd4mAp2adV3SY0BKfg9UnY4r1nAmnMHIIBiG4z1EVISNKrSWMIAEwEOXPYOwtqYGARqA96AbwAAuZwaACrQMgEoNoBJRLpKqERISYEvsevXCiWsKWPlij0nnNfF9X5ZQhKiIGHPC/a84CbtcYt3OKOCHa3OBm2Cs3oT8Fs69nxL1R3O+utO2aYCoOn7mBqLHftulFRmLMz6Q+gjnBjE5P0dxeZo05kASMIEpR7btzwMXsfT/9EfP2WqqJdf9Bimewr8AICZfwhtreUlFwdmlMLvAlAhUGEBvypsqeaMdWWklJBLVqNHY3qVE2omYQ2csKOKJRXsqijNI4sSMEm4UBC5te7w4fUMD+93ON8v2O8X1JLBawIKgVYDZupnbQXuKP4aihARaAVSgoBeZtCewAuhXCRc0A4pNeV+5YR9zXjyco5H0g430uqD34wbexYx+byImHxRMy7KgrWK+MZM3TglksGWUkWtyQGQWTvzCIDHBrqem7KfYudlELMbPYxdErgyKonoXyuBiFAqgShh1bak0gweffbyfnvOwghZ1ANntAZ9bhOHbZKTe6W9DfxMhH643MDD9UzBT1QPpvZgFXdtQqAgxo/gH/vz5rdKA6gkKoIq7JpMvFWQEsYnaghO1ETd0TASmWA4ZpoGoOV5OwaTyxKDsb8We0+bpLMLKzCwSwXgVRXAK0ArgTOhroSaRGSM7cqqM1qpooKwUMHCCRfEWBT4chhMRcFv5YSLknGr7Brj2y+4dbHDus+o+wTskwBfESCLImxjrIGtAu6ARCS4wkkBMJFafROqEsrzJE1YmbCeJVzUjFtlwVkquJFXLB2LkffcG9urAt5rFdAsNU3F36RiZUpVRHKqqCq6MlXUIgDIqsNrYpaNIrSJSuSfrYXbH6dASqSsT0U6ArjIDMdIqFUmh5USiBgriepCVBMscV4We29lhCnhPC3YpTPsSETdRDWIvFJXecJQZPLYsufzuuCReiYMWnWwa00oKmVwZLphYuiMDAPouei70ZMa4CkTLyw+CZaHMz4ewI+8jjkywiQf+72xFJ8U/IByLfaeLhEDqQjlT+KfIpauAgGalURsXAHOBE4JJbF2nqDUzgVrTdjlgsIJWQfQkswtoh8MBnwGGufrIhbSNeNiv2C/z6gXGXwhwJf2hLQKG00RAEe93jAgEqRzJ3Vzsc6KRAKEyFgBQI05a0m4OMt4ZN0J8KkbTNRTmmFgVXZiTGWtogAotemrKsPryYyG4rEjI42IZYCzMEBKUm7udHr2Tab46fSerTHDtSCgqBhng9jE3yLicC1J9ZCCc4mAfWxTm9BSwpoK1ixi8C4VLA580R2G3TI+U4XIu1On8nDwc13v4pPKWqObU2O8nYi7Ab35/zKn2CQgLJBUr5oC2zMA2zK+ppftrtVn1NwmI+7uOy37u9b5nTiJW4B0Cl4JlBmUCSkP7C8TmBJqMnEKLuKVSlhyRWESay+JD50NDkvGhgwsChPWIrq3tSTs9xllVcZ3kUD7BLog0F7F13UQzQfW5wYBlf48YC9FMUYHABKK+qDtTZe5ZlysC3a5YLcU5FT9fXqRXX0UB+NBDX1zo6wPIBhZIFDBnFBRQZznCvY48P19g3hsAw0cbiNnfWSXjuJvITXCCPDHZO9k4HdRdUJT9cUuFX0fqx9lff6erTLi5AHA2fOqagZh2zucl0Ut8GnL+vz3VvyNRo7o/6kNIdrQClBSAIRKO2oBZga4HmB79kgDRnUp4iyEoGqX6pp7wgTvNDHEdeheSvc3+LH6RbEABHvjAkgyo6W9+Ee1xk9NqLFBnytKrQIWgzN0e1S0MJIzpHXNKIVQi4AerwSsCbQnJAW+tJdyCBsN4Bc/pvsLZRP2B6BZQsI5VsaXUCukDLuEdZ9xkSvyUpFSRdL3Gfvx2A1FnA2DEHDA7I6TuAOJqkDq2QahucOMedM42M3wMzA/YhJ2oIUlZTpsVl8Tf1XvxUgog99sZRnU2SzBSay+mdpklgaAi8zez43fw1tVkFt4V06iQ11FhWCW/sokojrThvWRgl78vzMCMWz26wEwWkG0H1ANEyMMyAYwTAKINhZEZyvi88Yf8K4xv3sr3dfgRwzki4q6tBbiREBicBbQcX0ZESoBoGTSl/QfFRlrJpSUdADPwUJEPPkuRXU6hVDVqIFV9Xurgp6CXTIAVPaXCpAKD+IvO/sTIiuDvjGwBoBmPJDBQ6CSUM8YvCdwTqgLY80MJJaVLUlBzOqNxlGG3pePZAVJNbAAvF5iSlYqqzNzcYmuLoHpUAS+geEYaJIxW4K7LLEaiZRoNvBTIKxEsGki6QSVkwAfqRuTgfloCOn6UzgXV4vYuakbjKo/ShX978UqutOiXgb2DtHPc5wM0AFd+z8+Tj1AGwDqZNNEXQXEid5OSAEJ46s94zN9H3ErgrywMsQTgR+Dr3V+J00M0MpIQXzKMlKREysQQnVkwg6rdxTTbYn4VksCOVMaHsNQ/Y0tVwK4CPBBPw56BaLnK4HtlQZ8wvQGcTeKv9p/TSldVcTxFw4KbgE+gHbyTF5YxPzMqv8UsdiW+DWaAf/fl4wlFos5sddDSiwDIFXfEIKogQGzLAezQeR5h+RADTjbc/E+UAEbgDALprEODgMzipAEMJn4qw6ZSFJeAKUCidpkloZyXSbi2/UGejlVf7VRFWKuNquuiFlXFXuLWvrL4Lc58/WEMcFQsFBISsr2kwKgdhSZFJqbULuh1WsTd3WFiDI+SvAJND7WvXNG9cUdJGZZKXMvpfsa/IiBtGdwYeft3nBJ9H4u/iYVoYiQIDoqYX5qOcwVqFl0gnDBAuah70u5DPDUhxAKeMlcWUbDRmmgJ988GD164PMZX1HFXBgM9IhN3AmDaSUkW+qU4aAvwGV5WWfm9r+f1++sIJgJNSwnAwAMADgCxrx9Ig3pGc9IwHqAk9onKDCZz0j4tBUhOqirOUCzusE0FhtTXKPsfnjoWY8BuTFgA9DcMdw+z2KO8yWJ3nfVyZFDXwkMMPp4jvUxiscARKcHNCN6EHEj8HVs0SYlBb+6kDBZZdec1CBY+3sQ+omLw3ecxFn8Xkr3NfiBGemigJcEMsWQdgxOQE4yU2bXZ4gSvUJ8w2Db2FUW0SqxWjLDI6p2YO+8NLiuUOeE3Onwioi4TbfHSAXdahQKbMj8/QDtxGSDXjut6tkyBCBr1ectACvTjDpPD/BgndjqBujAkTOAzMIUs6gMaKnAQkD0zVIAjPrBzWRuoMTt/06fFUHsttravinkwTDVhSy5YxSGGkG2DzDjgOl64+S2VVQacKi460AIpFQ76cCAtJSEWhJqaZOk9w9G74Q/TgCTMky0E+CqKhlzh4o31D4DuY587XRiaN9nUKK+XYYbo5HkFIm3xXvc0/0NfhUCfpXBO+sNCUyM7CZ9dtM+J4hLjM1qNutrXpyCvskaylwVqgEduRuNA5uxvwH42ndgd8Pxjc4HoaND+574nKj+qQ1Wqoya5dm8ogc+A7ogjnZAaOCY4QvmRWxm8E4mg8oVrD2EGcgadCCZYhINUGIEGDkRQSocuywF8ZxjZcxudeDSeoGyTWp2466NTWWhbMzZY8c4QyNofVUNtECpsUmK6GVgqvnXNQGuB26szxnWwPg2ILh5z1Y3fm3h6TWbiDBVQLwuev/gaE+hHF09U/icKF0zvxMmYgbtq4AAA8mtXgCnNAXA1qjULQkyvyfvYTYQIuhFoDMAjOJtAD8MIm0npla4pa/r/BPxpwGgHHPfP9aOXbhZ8MzXMXRaHjqzz+iq/7EIMrSwRAdZCFwY9UxcWAxom3goekBnfmpxNibVrfCYMJd5Q/ZMI4J4z1qxRQsrHwDTfG1EWYYz+BhezL8dANrg5AB+TKz1yupIbA2Hrh5sDXmnB55a9YeJb6iLg+lAHdq6cP8dyDpp308rUG3JpPX1A8yvib63T9CPFf0a/E6ZmJEuVrHsMqNiERE3ib+eiL4KgJ27i4i45BSP1DAizAYg+HIkA78AemkNfnudIUM6dSpwsItiX/vmXudjSR85BYqg60uFm84qAVzs3Q6zpGj5A0GdvlVPuKgi/AygnejMUBP4rEVQYZapwvSAPfhRF+Kq+bS1d+7SwCgc3BJQMzv4mT8a1PDSu2G0TMXwwr3S3ievIeRWAD4K/3cWWBIg9fpK3PoGIZRlsIbGEF4lfCqQgtHjkN5z1mau/51cOwW9YCVmAlD0diJEjwIH+tg+oW59Djuh2Ls/nQLxJOm+Bz/sV1BVNwcAFQsSyQqIvIcbPnjf2F9S5gfV/wEMZBkkvl6S2+CwJXPmtmKAF+PxpbV3XRnBTsorXxtRA1IeE9a6fj50PmODJvaKywNkbeymfuweBargLM0KMh4PbicWUhmsYgGUkFEJLYwUmh5QDQomOlYOAOhAMin8xK/MdFI1BxZrv1UP6YFpLTits8DYHeKADqA3gB8pO+tZ+YTmhLLZyhqoHllY4PB+Q5+htUkD0dAxLBhq7UzWRvqOSRy6HQAnbeq/I/DpObmvOUizAX+87xC60fB9h4khwV/vpXSfgx9A+xVckzsIJ10QK0CXkDILCGZyZT5bpJSVQjggRZUUBxC2Bo0B+NKeB+YHBQ9cTWYIY46gOrrREz7q7uyQsxTuj1m9BOtxO2fLnJT1ZbEA1oXdeFMKXCdkBZOf6uPIAkoUmabqunrGF+oRjXH6OlIbdxGMI/AtA/Cp36Lp3mxQbgwbBnwd6Mm3RUEhNRTBdF/G8MesHGSbA30HhAP42rO9H7guuB2b6vn8edofk4AFOJC+aO1lvZZjAwXg6/Lm7UR0AOQ7F5dw/FRpFi/y8Uz3Ofg15scAaC+9lZKEPOKFkFZ1eF7hyn0bZJQAUncY8YDmBkQ6cDpDRTRWRJHXZ/jBbSUWdZhJx04W1Yyz6Xa836+MgOfMggPL4DDo2JmFgEyLBlxWQi0SNaQ4cDcwqwzVAaplPKsm3VLUd5krzmyQpfCiJjYmdMBXF4j1eVHACw7brmcbq8jrIVhvvTxoRgefzKiBnjNA2riKWBnJrefk4aJcvxrbIrK/2ljfJnhFrBOCGNt0YtFpUOct8uuI0ay9V5lYDyXvgzRndhQ+J0qMa53faRMzsAqqSfw3AqUEWiqoVNCakPYszGYIc9VATcEACDpAyEAw/UjtO3Qq/SqN1Im9A+uLoGVuJjMgPPaex/pMYHmd36CtICncyqTfwrRIVQKEeiaAUGLE6AAm9l2ZwRrZmj2acytHr1cb6sBYXkI7of93Yu6ibC98bJWKsc2RbMHdVkKZQzmmwBct7x52P2aqLNUDuKLzm+Q0aU9rD2uLEp4Z6yOwYWOUBtgVYtXnpOw0uDs5W2ZpV4+BaJP3VdbOHgO72fnLFJNXTuQhwu6VdH+DHxi8rtJeusCdUwKtCZQrUqmoRUL/pEKos2gqVfoMsX5rSHbvrLV9+8zdgUtkfE2hHAc+6wQbcMSPu55Hr0W4ztMBlmO+gRLOq52j2oAvFQatCoCF2wAhUhWA+KeJXlNEw2aFVhaItpSuchKAsnDztrwmDFAyEIxFVqYU34mBoNcDMAAfLRWUgqNx6gciaxnrWGlW/0y9Dm4Evqj3G9vNihnaCYnacjBjhCNDCuwvAuxo/LF7PcqK6d+06qjqczojyTAx2URFEKZIYQI3oLTlN1p+eTZ1Ft2x+rxtTsz86rXO74SJAZQKpgqiIqbdWgFlOyjcQGoCeq2jIri6wMUKlXo6Hc1B15VO5G2dz23KBB2t1gFV1OlAMHx7QcLvwBwcrEnXvXoYfL2cFfBWRioVtIrsSqXJXZzF+kOFUYusTmiD1B4Y6anQxloBZBLRNDrCRpFvwn4dV4w5eXQRFXEzA4ustkmLrDAx8JMgqggRbzQoRRUXp2rifBDV26RlQDQA3wEDFbTdXTq3CawY8MnB+O5Bku/6VWTlIxO2PsdZ3VDcNcp0ruyTkedX5Lgz0tBGVABuPQ6uK04TwDPWOfa3yybeR5mYCRecL7/wMUz3N/iBIV6lRfRQRcRd1AoO4p50KhP9qLE4Az5uyuPoG9aJUWFwUPc7+OtFHRt0oED1bAxVkkM7PnUi1MFYbGjMMFqNvewqhrs/XoJbfknLJ5NABa06Mbj1twJZAiNQkSjIxAnUhfgIg0sBkLjtFoaJY/hoOeVxUPnKEpbd2jIDSwUtwvZyrkhZv0OkaksekquKBUJC/0k5ObSZl8OAjyfAF0TSCFCdhVWbDQ3/5XBoJzvfTVCDGG15sbEwLawxveZ8zI09erkUCJUNprU5chMbC5fMGG2ShaqDXO1ixpsBCDdgaDP/CVM9FZKeKN3f4McQ6j9OKLZAPDKADrjYAQ1oLCrqVOIzolGhe7Z983BOE1XAJ7tIoCbOvM23bcIGW469m0wVvWPVkZcgxghb1mZL9UzfJyBYRPy1V80ElIxUuRWWgSajmhjVZKHKDCyQZ2X2Mo914/UW/A+bgSOIuAuDdhVJgW9ZinxriLGcGluNYcVWGAtMsnLDn9+oprOoTsfXA98o9s7cJR0QzM5jbHDTRqG1Yn7Wt6z97QnETbw1oOZmbW8TagPGVCRAB61aJ4sY3jzL6PYUJlXR9YZyjGA3SacSfRm4dnUBACJ6B4APQjjKyswvIqKnA/heAM8F8A4Af46Z33c8I2x9rcZLgqvHOJFNdTwIXGcCdqNoNHvWoRR1fE3Rb9/DutwJA7RyxAEreij2B1BlcXrOANZ2o4E+FVbVgIqRq1huueatsSbG09daMZGYKlCXALbR9cUqMnw3fzl28DPgS7uClBnLrmBZiu/Ju8vFg7Fa8q0piwC1ib5kPoTWVgokLvYyehG49vXY6WqH9jWdLSg8IgBf1PuNYDH2MbccKwga0PmSRVNfdOoVc5QWYx0nNbJZTyXtCmv/LKt/cfIPUseE+XVsL547Wbo2eMT0x5j5veH/lwP4CWb+WiJ6uf7/Fcez0J6kxo4N1QfCt94yw6fY6QcSc4jVzVLne4U2w8fO5ta9yPjU99CBMHzsNbvyoDEXJMi7q5LbBpPop9Sx299dKa6qBsw3kmvSQcebuZlpzgBrhYhYy7C8rhtM3LM925TcgU/Z3lKx7CT69NlSfCPyXS5IGFaToAe+Uod1tl27UQ9ubvyAA0zv6jKf3Lz6IqDrpzOIxHNjsutUIHFgHSULjqKusT5rT9nSgEoDeqLG+ir03FAAKZ/6vtp3MvGXMAPBqd75DpLMPdfgdyi9BMCn6u9XA3g9LgU/qM5JtdBJgnlC17hGHRtwnN57mom+mo6qQEK+/eoFK4M5yhJM9Km5AV8LR4UpA+yNCmgbVZsCXgtteq2U44xPjb0AALPrRgGAjAlWFlE6vhMBfiTo/Kx+uJIuj+PgC6fsRsW6DQCmOfDd2K04ywU3lhW7VHCm4Jeo+v4jq67mqUy61UDqHZ2Dxbnz3RtBbgZ8o6FmIhXYjyjuuuP20Fabtgt9Y9aVOp9NK1Mw0qQirD6tWq8w1k89rCRM3Jrg/c7LNIrEBnyhKk+ZyqkzvMP0eIEfA/hRkpb7+8z8EICPYeZ3AwAzv5uIPnp2IxG9DMDLAOBmegqwLKCUgGUBUgJSkgGfUmjUkfLPG8EtfGHmp/B78xJj5x47OgX9WzRqBP1e53jtIrD9bk7AXTLxrcAHopSRJDahuvWkhXXLS9KJwEalMsCq34C4lnD7fztHN0Q3nSkxUHfQoBLkobGQBZwR9Z22PledlimzrL5ZqgPfzWXFjWXFjbzKJuy6AZO8sgBfQnbgo2rBSq3C7UP+mwLrG8XdznA1iL1x3Wxs/8iISNu4hbS2d91OtK7rRTs+nYwDswfgIdBkR8Lm+2eGC9H6DABY0YOYTbiq72s64VDmUC5/4auQhSsmhux9fC+lx6s0n8LM71KA+zEi+qWr3qhA+RAAfOTZRzMZ6C0ZvGT5zllmuyU1kTIFqxfQU3so8AE984sD6VB5NEagd6QUrvfO1dwMNszOGSDaWltz/7BrLK9Y1mpOv2TLlF0UpcKoK1AzIWVb06xO4K0i5YaqCqdCQO639EwYZ38ZXquzTGVZS8PRtlYasECjXhfGVB34RMd3thRnfDfyipt5j7NUsKS2EXthQlI0XThhX5suMMbo64xboxEjHN8aweL/3N+DMAl6u1vbmhxL/n7Oyo15Bca1+dY8nEUG2kUVshOhfWdZUgkijTHYlQgOgMS9n6VNtFHcjRNxBMEZez1BYlwbPAAAzPwu/f4NIvoBAC8G8B4ielBZ34MAfuPSjIiAJQvLWzKwW8BLApYEzsnFvs6YEGfmTcEa7s3cWjb6vwikSRiRY2dkmpHphc41ireNAXInAiP1g8OcsZEg+2xAd15TZsO2VaeBvlmRVQSOy9K4qribmr7SzyYgB51pDvcVFYFhz4ZmAwVAG9DjxKHszzZXWlLFksXIsUsFZ2md7jkMJN9TZJYa8CkgWxDaUcQ91K7W/gH4OleV4bGGOy7+puZh5xOWXRjb29uV2zrhyfxKYNdP+trgFW0SsRlPryaGh2hrlmz2CZ/zlu1tdNHxPML3CRKDrsVeInoygMTMH9Tf/wGAvwXgtQA+H8DX6vcPXiEz8NnOQZB38qm7pJ8GfDXq06bAZ4M5zPIzwNuUAc3yZ75bJiK5UrldNwPC3sgxLPDvdH+6b4Yxv0Kq3hFXCNaIzhKctAUuSEsCdC1ux/6qOSRW3Rc3QRxIbHDLlTn8BnIYFD5VeHVUrwthP7zRI8Dj4ZnzciYFwVRxloXxxf2GASBNGs233NS13WBylcB0ze4AdM7u0B8fJ7sZA4xt2uYE9b0bQN/bPU5suR33yS2+ojHXopFdVuhWDOE667NaxgSJUZlsyabPmHES3E7CM9G39dtjnf/20ikMHkR0E8BPArgBwa9/zMxfFc7/pwD+DoBnDgbVTXo8mN/HAPgBsVBiAfDdzPz/JaJ/BuA1RPRSAL8K4LMvzYlI2F5Strck1F0G75IE5cwSnLNm6hp/0+ChY/tc6uwB2wEBdLNlm/3lHo4dKHa4TtTp8xj1RdEVpm1ABEDDS5n7CTOpBVWAr2YgWQSbNDw/E2hkfy6vCnVkXWjajUUi1S0tqmSKnbgBIIesUyLUxM2NY9J0tlGQ7ZJmnxwGXCbuGIPtO+x7ZrA4NrdwWtSvMBnAK7b3pkyDE/nBb3vrAPQhE4Cbs3XUEUZ9bjWrt05yrT80dQIbgK8Qi75uwtW/hzD2qpOMASBRAHxC1web8Ylaf7SyW5lPTNLE1/4kYu85gE9j5g8R0Q7AG4joh5n5TUT0HAB/HIIfl6bHHPyY+V8D+IOT478F4NNvKzMi1Bs77UACfgJ8jfk1HRqC7iM0uGXFgC8UD0xh4waxKQNU6c3qbBXeKXa6YZadKZt7IFQZ0sQj1/EoKGs5Ocl+DObWMjpMewBO9/QPwJcSgKIAqCzLHGTdlYKARB4qrNk9UntHW1oWDEs1QUii1qHp5SLyxM2A7LuqeGS7kRYmrDVjzwkryz65vml81X2LS2oRmpXpoU7Y3iQxRTZHHQDady8uh8yiDtnacPK8kdXXLoCDgV8AQHfQV9aXZGWH+G0Gpu3vZw9uAGj5uBHPfUtD/59MxBvQOxEIisHjzpe3sSyC/pD+u9OP1fg3AvgbuIrUiHvL1eW2EycC38hNkbsI46u7hHpGCn7G/pro2zU20PRoQDsYnWBxAAAD47OOx5afP4MGBgaMhpeWH3d5+3sOljfWKNTMMsMz2nt1Yrbqhnprt4T86hgEgLa7TAG4tiGmoEfUPsnddhKyvwdpOZUBJtFVlcwhwooNan0PE1vtG6TMLmFlRtJjhQkXdZF9cat8us3Ba9tG1K26h8Raq9th4hlFVWuO3iLMQ/tbpyFtp8YEoxThz3DdLvfrmS1mYbBceyBWdWoW/S7aBBKMO8Sh+fS86PwCAx374ACAHRk45uJwB+lUBg8iygDeCuB5AL6Fmd9MRH8awK8x8z+nA94cY7qvwQ8E1AB+LupG0FPg8y0dhxmuGRB6MBythuNAcl1PBBobBKF8G93KmDpxrBlNutdkBZbZ7RMc64EzloOCspwc2Pxa03lWSMAIrAJ4ailOdr89IziVu/tFEkLJGSKuaeBTMqfo8KlVlqkxgJUTck1YlVHGwJfC+BJulR3Oy+Kbg68loayyW5rHETSmXrGZqFq526Tl7R7bs6tg9MAX+4FfL1MFWZ4cAix4/whqiAxRpHogB0bcIEnaArp6RgFQrSjVmN6iel4W1S1VYYnVRGHqwVpAtwGpi8JjveDuJNaJ7QrpGUT0U+H/h9TDo+XFXAC8kIieBlGh/e8B/GcQ+8GV030NfkyEciOHhiU1ckCBT0FPv6N1rc9Ivjon00Hs7dwj/PloYNCkjpYGFraxpg1lcCDWji+7FG2fG8WdqWgXCeRYXgDm7e/FSNTE3VrBKYFK0XuKxEg0Y4mBoLFa0o3CSbcItTDvyrLTKvpHKgAXyEBWcbWUhDXLRt8EuCV34YolLFBdOeOiZOxrxq11wa39opuDC/PrtonUcFrTYWaDO2lTqT9cN4lRWC7odch9f9iAn1Q66YbiG8NKeD7U2GP6PsoSxYZy2yLTGqvqpFMhelgPXlHhy/pQIa5GKu7LBlfqH0gh0IP1waASOTop34V0Reb3XmZ+0VUuZOb3E9HrIQskfi8AY33PBvDTRPRiZv71Q/ff1+CHBJQz6cm+/WImt3LayomqrM8tvpH1IUzelgwMA/BtlOZAz/pm5zFhfhPgJXNKNVcVO1YhHVi1kX4vIzAcWw7VWM9RHWVXfnKxWQBQbzAAtIg5axN5jT3KPiklvFfq9IzddqFa/5QhcQALoRZCTQKAe2q6oLoQSq24oLB2F7Kk7VwZ38W6YL9X4LNtIjVE/Wjs8GpOrU07ABzaw1QWGzYdgW/YI8PamdS/rxOX+yz6CdPD8gvwk2+Qriy8CvAlVGFNKgqbccu2ZagVMpL1AQkS7MIj1uhzuyWVo/pn2jfH2fzRJwZQT2DwIKJnAtgr8D0A4DMAfB0zf3S45h0AXnQvWntPljgRyg1lIRHcchN1hf1RF83iEPMD0Is2gfkhdugoIqR2joe8XEcXO1l8ZnhOAzzStZsIe4voDB56cotLFz41fFcVe0I4r14MpLbqgwjQZWKcUnOBqXpTKX4dhftsIyjW+s+DocV9DVd1tl7F2sxJgLIkBvYN+JhFvxeDGVj4qn0RMfdiXbCuCeuaUfcJvPaszyeNmKwdIks3kIvXcmOBHOaavq+MkbqDyEvcl+HYBGRuP9QCtVrQVmsee15NFnwiBSMJNKqOMkH7MDwadGeA0fcyIsAjCN71RKcKY/8ggFer3i8BeA0zv+7RZHR/gx8B5SyuWww+fTFayiEfv0nn7FwkRhbRiZC6z4LhE450dEQADAXQZ3T6w4IWqGDRG7O6jJgFVhlAC6FPvqmSR5aOewh7+ScFTEnoiu2tUVitwFBroRSQSwXRKsC3D+KvWZizrKvOcWBFpm2W50SoavSp1JiAAF/BWpI4QGs1VQaqWXdLwrom1JJRVwLvhfWRA6Aw4Mjk2URNVRE4mVH2vGknZU9UeQOAnfHD8tebSKPbtElNH+Z1r54EB5KDYOhEZI7MakjywoS+zJnBahTh3MBO9ujtDTQcxkjMYzMxo70jnwgZGTiVtfdnAXzyJdc89yp53dfgBwW/rT4DE7cPbGa4jSVw7NxDR4/PBQYgs2wm/bvbOjI+OzK+wEjlMh2lCWrwCMW38WA7yvn2mbL8yXeXcyDUyNbjgHeXF6M6LUiAWZKd/VGVgLFr6fR/lMUgktXoBKrK7MzXMIi/CnqyLC+pKCQAwQpuOdduvS4r+NVKqEWMG9VEXd1PGUXAv9s/pGsnycgnKAMlGj4MB7BuQtq0Z3uGA22YLMcIMaP1eeb2w8ruu03QNw9GA8IIYMb+uK2yQdT5hfqIk9IG2b0wbUY/FSlkteLfS+m+Bj8moJ6F2cvcSAbAc30H0FpzBLz4ewS9sV/EfGa9I1wvnX1yURgYY0ormrVPnZU7nLWF+qWxvXxhwMftWz8OgKVO/NQIlJLo4hLDohL04qDS08ouAvt+KXszgAgAghKyWzZ7/Z+zP9J9N1Q9XzmDlyrW35LUkNxAgFkMCWzGjZWE+a7kUZl9p71Y7Bl4KQNmFsLrukLoq9ukZGVXtr7ZInNoR4t243spR/1rvC5uNFTR4iFyA3oidjDcUrL+/cS1JtRxbu6m7sNoQI8ZSWiPiI/qWPKJ0nU8v1OmyPwife/AsF3bpctY3qFGpyHPQ33TxKLZc+x2JWOj3glZvxPaipHYMcNa1bZ/MJD2jHwB5D0j7xlpr3t4rAp8tppjBMCUQNlcbIp4KKcBlSsDKFJeagDY9H/qa6miVcpQHWCvX0omykEBsOp2mIXUAsquP5P6oObAXBT0qgBg03HK+SjujpNTbFa/rgJM6hSsxzhLvTqbsrKbemLS1JHFd7rhaDCzsmbIbnIan4/NRUXWAuomTeR9h/39TYQeEMrEYZvsG64quIeuFcdFEHm7vmX1oBPz1bxTLk9SHafikadJ9zX4MQF1h168iQAIbHvrAEQbEeXQBB/zHnUlB9ifg11t390AtWuo/5+5n5kJ22t8MCn4CdgB+UJAL35oVeCz3dsM/Ay8lEJQ1jWy1uvrpDJYxd9ahQXumxU4ZXL9X9bAEu7+0jl2G6UQRlcryc5tauiJq1l8wNcAch7jjjaTV9fmihcxHKFl6Vbd1IJDgBUrqmC/qE64n0zNPSa2Y0ixzT30VwBpC07KBRJoYk3K1Cqqrg3vXr2STAxxR7woiWh1uhU5lol7QIzW3XGZp7FEjnV+Mo2fPPya+Z0yGfjpb+AA07OxHr8j8I2uLJPOEmfMze9JD3GfPRsstXUsN0KE8rneL5Stm52Ha83YQS7mCvAZ40sXjHRRkfYVaa0KgFXD2A/MLxs6SGdnAz7zuxnDtpTSVGYT/V/KYpXMeSL++vvIQK9MAjQVajxRoIxtxwH4IggCzmyiS4e1XRvk3LUjIuPWrR8rGgAa60OJk5AaQYboPZs2t8dZO0VvATNCrbZczbpBQmXyjeBN7JU6p7YZvO4/3C3jiw83lU+out4yHeqIWv/aeCFY3RCdlvmdKrMTpfsf/K7wBiPbmiqjfRDJhZ0xI4oKqmiuttdPnEEtDcDnMytB9EspAGBkgWFwgtAGeGR9wY0lrcr8VjjoOfMz4NsXAT7duc3j9/m7aeY5STTnnIXVkUQNhoV2ngBgdH9hBT8XfxUwNu4vkQEqALbVDGis2q5wVjcsWxsnkFhPEfi6jcYjZSJ9D63sDFRmCcqqYm/NUDGbXQSWttVjmCTuP3Gi4qSsNbEYauwW7X+2EXyna3MjDrXN1+PeyuF1jOF6F+Ihr1A31qen6iBqXeRSX9ErplOt7T1luq/BT0z9mLI8+33QkdlApBNR2mlxYyGfGTfuG2Yxc3bGQcFsYad0pjfgs2QRmEOZ7Pkx5BQN7yV78YYBtZphw/R9sj9vvqggY3wKfrZzG4LeT9hMeEZOyvia+HUMALlUEFYHQOzFMpOWBOxp7v4SRC1rRAEb3bjcnK4j4Nt3mDCabpa87hHvIzjwGQh6gAjA28eMGc76WDdlCgzQVlKACclWWihwd07E1NrS91KxndcsGKnF5COoBRQyCjO1jeBjCqzXjTpBzzmbAJjQ9MYz8KKhHYLfDkH77V0gadd7eJwyEdBFxo4sqrXnxpobPxH4+o5MrXPYwF0C8GWgWkiiMIMyANs/ARUSg63AAc8ti9Bnx5mZG/P0V5wxVd3PwcFvH9jeykgXBWlfgbXKVpWrsr1S0en8gAaANt2bvK5A1wEgFDEAMYq4AaT3/6OUOvGXM7c9RdTXL4JbU+zzVMfZ3p+2BikSVuE+l+REsbG/bCH0Y1eR/2Wjb9JIKj3oNQdiLZuWU0Rke0hsrFBua89hspKq1lUYDCQD/wSPztMle2cXn2mYtMO3peCuc5AdW93NRF5rlxMCoIS0uhZ7T5pqbq3bsa5DotEolgzA1106YXx1QduxLAakNBnEyqAdsiaNrUZoe6sOKZbzoNuNWXfVgVmcm6NhoyJdiHib9qUxPtucaLTy+sOFockzzR+wrxAHwMQCeracQIGcV7T1vymBloS0T+BcwQuB9wZ+0DXYoU6gGgcGYKGbEFnU0G7hmF8TfvsBA7sUgM+MFyE/tvvZRHB1LzLxN37YSK/GzBv61wgk3g+LhqUqcOdt1j7i23+6KNoDRNcHQrk3fUYf7tLHyISHNOJQ77Stx06MVdc6v1MmFXu94Vj/qKjgerYxjcASxQPqPzPg6/bZcIdR7vN3kVfFKlKv+xUHyL/5iIUyKYA2UYcd9MS3r6rBQ9meibkGfGuBb1E56Pk4hrUK35yAuK0lIP+7VdgqKobAAtr635SANYOzAbEEk5UAB5CNlUz/pYyZIiAFPakXN7ZPaCdidKy7azNCAD34OtoODDQuog947TPO+MqW/bkvoBUgqg2AJjWEfgBS5ofmg9esyqTSA2Oqg7OyTYDuqGfCmMaL40Xh96l0fJsy4drJ+fQpGegBnVV1BL4AKs3qG6c7uBgwbipeLTLMwqg7PbZYPDYdZN2gYhlISZTbiUjDu1ukXTjr6ZJO21HUjWJuE3ebD1++UMa3L6C9fpcCrKVtTzkzcFiBPVgBddcwBa8szmLhtGVwxvhUBwliMYAktQStBbQkUEmgNcmWi6tMGqkAXAgpLHWjqArQNnDRNxa9n1+mQMF2oeljNYqKM0Fil/Cb1Vmv1cALzAQuujZZRVIzfpi4Ku2nE+yMUccCKfsDZNsBM3AgqT7RAHNgvTF1oBQnajsUJoyOCR9CM7eGDJfdNfAD9tfgd7rkE6IBXaT7V7kZk74xYXzG+ozx1Z0AHy9hUEVdjbMEkggmCV0o8kQEqP4nWactrQDkkVWB3tChABiAzwFvX5EuVmF7pYDW0sRdz1oHWGR85urSWbcFmLw6M8MXjo5MJ67/rVWeu+S2oqS25XXd9pHBZcOdcRWHXGybDXhro7E5KQDA8NuA0PYOsVdno2EmnppVOqo7tB+44zPDt+SsgOvhvIwzcdLrjP0lkz6bCnrf0QPvN3tfu9YlHGrHIwAeEnG731HisP9Pmp5AzI+I/gGAPwngN5j5k/TY0wF8L4DnAngHgD/HzO/Tc68A8FIIDHwJM//I1R4EdLob13lse2LnDBtFKc0nsr5uXfACF3HrjsEKfliqBB1Q5uf+WQz3zbJAqx6KfE9gYoknsJIDd2SCFWJV7CzCygLTKgCUVlu90QAQawHtV2Fhpcp3THkwjQfgc70fIFZFZmAxMBKQp5TAKftqj3a9jhoVr6lUAUL1KSTdDH3j96ZWUd/iMjbZbQy+3nJp/7M6S4+f/vXHe7vJL3xqlvmt35BJrzXQGPpTdDXxiOCsFmjzyRvBGsAIfht9ph1LDfgMsEm//VoOxiBEcOtjQW5+A7fVBldJT6QVHt8B4JsBfGc49nIAP8HMX0tEL9f/v4KIPhHA5wB4AYCPA/DjRPTxGrH1cIodZsYW4qVXbcix85tuzwwdSwM+2kkQShFt2wM8PHshXYuamn8gKeglAUuY4p/kXFOIYxvwwEEwrNU1i24pCnrF2R8H1kcm65mfXlia1jZ5lwdRYmF0K3oA5CxyWhajh/hA6zPMKuD7AFcHQhGbLQ/qJ57hu5PGZkAYz/mH54wvfBrZZW+jzQQZ8rAldsYEKQDhxow6TKadQYZbt3QArICtYomgNwbKGC2znqcBXOinXYRtKBA6mw4RZaK+MIJd8KHszp8oPaGsvcz8k0T03OHwSwB8qv5+NYDXA/gKPf49zHwO4JeJ6O2QvXzf+OgLMBlk48xmaRg0o8gzGjwM+NLCSLnoFowcFqSLy0RNCTUn1MTgNSmrFGdgcwvhzOIOs7YBKnohda+oNihUXGMFPrPMatACrEHUNeCLRgtjdl0MvwB8tqsbZLCgCPOTfTEYvGR9tgZByBlAaX6CITyVi9rcygxb6nVoQMXBHf4ff4+ibVylIGyoAWH3GfM8lnRS2oiTDn6WmbJz0y5wz/ZaoXucdZ0mgB7vhsqZvLe8MzWGl+A+fcTCUFuboLkOdQYOPW0MvLb7RxZ4yvSEEXsPpI9h5ncDgG5ObtFXnwXgTeG6d+qxTSKilwF4GQAsT/uo1pBA33B2fQC7qbEjJHNqnoMg60Yz8m3AtywVOdfNZtqVCTVX2WAnJZTMssPcKu4gNUEWua/wQZUG/zdxmZGVFl2QBkD1ZQKCzraYG/BFkZf6TmeWXk6pAZ99W1LDBpEyPlbAW+SZ7W2Pk/NtHQ8HCFvAmx2b3WNWdgMBBwee3wf0EXbGwkwmRbuMQp9o1wWUM+AzAJnk6cxvXrQ52EwZX9wPJPgmalX4SpTUQDe+qnfVwPbMuDYeO1W6jT08HrN0rxg8ZrUyrXrdzOQhALjx7OfwOMNGi+7MR+xgzhPQiy4IJv4iy4YzKVXkzMhZNttOqSKnJlJZEM6iAEhrRk2MmgnVgDBbZGNoYFDuqQDDVwa0dcTUW2aNXdXA9Oybh/wC40NOiEFIOac+X3dGrl49jMbiCGtfjUQCjmkwnox1HOp6M9nEAR6vt1cd77e2MZIfgxBMHTdDngGoth9q14dnUyhrfAfY/sR99p3hZtSj9SImd+dm9WZO9+5Xqq5ClbWdGG1Vh8U49nuHJpmwvE4XG5ngiRJDNqK6l9JjDX7vIaIHlfU9COA39Pg7ATwnXPdsAO+6SoZxhqXYec0JNALh8BmdZb1zdWJUOw79n3SnLQO+JVcsucg+DK04bWPtJOxQIhFn1MyoS0XdZwUfBaAL0kHbBp9FAkkr6Ybk3Mo1SzO3iyjqGttLCaxgxYtEYEHWNbrMqpyUJX5J8yXOat1VwMUq91UZgBTBNbjQ2H7BcQD3ABLAfazvSWO7QSPqu1Kfp+sCQ7WQIZGqJWCRUsLm9DPpoatefS0rHyVMN5tylyvEiUP/9+fwFoj8QeM7c5uYLXw962Mzo4Jc7LX7N36QsSrjM80QZb8NAGM5TpCe6GLvawF8PoCv1e8fDMe/m4heBTF4PB/AWy7LjIBNh3XfuPB9VN83y3QEQ3eQZQc+0fMZ46vIxFhSRQx6yUk2nilMKDVhXxLWXGXT7VVE4bIm8F4MIjVRWPpFbUYuhLrIcjbbpU4iJU9YljFDjeQr0ZmVlRnDy0l0gJlQl+Qs1Lf29GeLPtE8Qqx6GHCgo6pGDmDzDIvoPAO1MVBEdC/pmCDQg5gdMzVAam20UQ34TVJyVuADEIBP+8WwB4hHj5n0mU4MRgNAfYLUkYmcMyAdgG/DtCYTGBlLZqhuGbB9a8CQNcdZl8txX7cbN5hhPERxt18/fUrq9wQSe4noH0GMG88goncC+CoI6L2GiF4K4FcBfDYAMPPbiOg1AH4BwArgiy619NpzRoZns3ig8BRnNm3U2wLBySlLiYCchAHmJLq/hGb8AITuVyac5YR9LVhLli0blyqb8SwZNWfVBUbDgUbvKEBd4SJzikBFFlVZ3GmQkujkcrZ2EFAaWd8SP+Q730UwparuNMks0eTvb8/15XPGLlMCcgYvArAO1ENQ09FaOQKhsM6hsjvmh441OmhGg4fX45gB+mgptt2lxdwbGdABdmaGCwdAl797Mdgsr87+KJSpY3uhX479U5mfbWlgz0kgd78U5qcAb0aQ2uqm687xnUZxN+xRfKokVfoEAT9m/twDpz79wPWvBPDK23sIgoK2XwA+ZX0BKP1+2GC+vceOSZawMTKpDjCMkp2BYE0oLJvx7IvuRZszLnLFPlWUnFFTRlUvWmMftAJ1JdSVZeOehV1Uhu6GRgpsclOLDON6uCXL9UsGlO1VA75dQrVApNRGNFUIoO4VRAEkLIhhrFCrbqDNKjpnfVYDQNtGtAdB6tlg1LHa/7Fdhu/OP04Brx2btFDX7gH4NKpyFwrfWF/Fpq8Ac+khWnAZvRhsAORApEt8CNwDYSjr5hkGmu6r2J4GJhV/mxhszNAnmFh/1u834DdszH5C8AOu1/aePPlMHRqyE3sPsb5H2bCDF0eXEgQAEzEWKv7bUlXqv3LGWhPOy4LzvGC3ZNzKCy4SY09myM3KTiT2W12BuofvT1yzRE7mJYnrTE6yt0UITCqFSk0UXQSQeJflsyQBvsX0jjpgvG6hKxBkiZocJFXnafBSjRHIeo5TApbUANZE6rCd6Mj2OneVQfeH2FRhAEe9n19/kPVR/79GRG7A10JFQb+7vjM4Z3f9oQM9KDvrWSAT5jo01as6KyMRnKe6PwSWKZ0JgAXxCBFnjPF1EwltVjVGULcVRL3efA7yjzYxrsHvtGkmlkQKP7FgjQ17sD26GZBgfmpg3XeBSay5VfR6trm2GNoUAFPFQhVJe34mRmFZ5rNywkXJuJV3OC+L6w0BYM+Q3coqQCWhnqnYuyOUHZD2kGgpJlYusobWfPG6RAResoiiysac8dln18DJ/dvs/ZU4YhUQ5CT79VJOSHtZA0um8wMc/OrOPlbOAHQDAJohKS7u70Teju01uj66f2xYX3S4i6yPG7PzWHsBANsxbIxls4kzsqro59eBoLK8aAhpL8e+1le2wOQm0QwiMtm72I8QRAHMovMb9H0cmfGQZr59vSV6e8+jSQzCWp/YBo+Tp07nFzzoZ8A30/V5Z40pijpsLFLX61ZohA+S3cbUorsSI4dNf4T9VexScQDM4eH7Knq/m2XFw+uZXiPnmQl7e8bKAnwLqZO1rTMW4Eo5ucWWagJz1r0tGB69xUTRnYi+BoAGTsb8bNP3VjfkopqAnlUGkFbZc0KWsqVWj2qIMXBlFXtNp2h5db55Cnwby3pgeFKgHmi68wesmu4HE9q22wvE99aQfpL8WJMeOmtozGvyPC8fo+n4GNIWDPi2nN1HUdP7rwKg9b/YLxX4TMUjFS/13rFAi4wdxN4N+4t5hrIe1DveYXrC6Pweq+QRQQb214m9YVbbzN5H2sOv8wEgy9W4EGpOKEWsvmtJyEmsuDWRN3Ki6gB4llZkNICrScJ6n6UFSypYUnFLcVF2uS/iFF33hLqDfBYooKjebxGLbVqThqAfwqubvm+nur5ddj1fx/4WdHo4r9MqWfSszXwUGcn2lTDGSeQGlAauCtgzsddY3wh6HZsb2iqCXgeUfLg9vW8cYHwKhm1/jCPA58xr24dGZ2KO1ysIxnIbKFmAU8uDTOVht8U82lwQwFBFXJ+wLPS+/h8ni1jm0FnuBuPzxNdi70nTSNk78NuIuVew8AJd49uCdQNSrhA9TSHUlUCUUBJjTQm5ZBVrEyrXrqEXKtjZJxXdrVZ2szqvC3aBFRYWY0gpCWU1hqahtByoOOjRdLncknSlR9MVQq2/ZtWtQRyVIA0N+JpRAhvw46QhqKgxN07i7M1FB6nWGxNa2XbUMUvTV86Az8EA9jsAXxCBRyZ40Co/sj1jgNGY4QwQUz0fQvtH0OnSAChxfPv8yugMHo1pSpQee4cK1af6skQT0cOG656nIKv4ZcqlBLRo2BabcGMMok05W54H6u8EiXENfidPLVoGpu4JY8ftGvhIJ7AZu203qJtj++CXNbuFEogE+HIS4MokIqE4dQo1TcTYpYIbtCJroQsnB8VEFYUJa83Y7zL2a8Z+qag7jYa8BIOBxhisi+6Xkau4qSxJETtLcAIisQZHt5ZsrKzl58xvBn6MFuopQXwNFbwsNp8HYbU6TRTypT7/3D9nzuD67x5gBuAb02TQzoAvMj5nggeMZZcC3yGQjmUyANRysD1TwS5paKvk/plaJ1Ue5AAY31En9VgeAUC5Xxyi0SYtgrhgBRF4E0zB6+zU1O8a/E6eporajZjbGtKV+YOo4p1gGPjsg0bYAQdH5ErJWUpKWfU5ovu7qAvOuHQNniCuMAZ2AHxHq4KElTMu6oIb6w67pSAvFetSdQVGD3wCMNLBa5ZNx6mIoytRBUN6eFy+5no4E0tdL9dAcLN7WoU6dsNZH9la49Imhg78lP3VAVgj64sWXRd5Y/0fbGz0wDedvAb6tRF5x48xwUn/uSxFAIyAHlJkgHD2h7Cvi4bFJ/eC0Y3UdVLR/P31PVr58AAtD4Eb0yRl6Sr6elAECteOdXgXMIoh+vF7Kd3f4KedCWjA1zpIuMyWbPn/+uMQuxieQawztW47aBZRUBKXMQL2ezTwU2PHRco4SwJqldsGHokqslPRgptpjwrCeVpwllbscsEuF+QsIbMssEJ0RxHgC35zCjjEBDYnuQQNWCDAJ/eNYii2Ym/oo2zsNwkI2ppSTmbppY0LSFx0b1GwI/MbRd3YnLedwkTmDC+eY8D0YE0lQlMJIX42ecQU+kuns0ztPA/v5uwvTL7O9ABBvBL6qJYbSVUKxgat2wQj3pTtKvCRMUADPXO9GdZV+3tNXvdU6drgcRfSRpc3zIQ2QNj+R3/enEE3IGgDQUVfibxLfedHUsmEsd9nZX6MJS3Y5YKLWrCvxQ0cN7EHAGd+4IREVVghKnZqJMkky+dIl9aZRTQGwRTAgy4zESOJrOiwdzOQE+CrSwPLxh4b4EVjRKxbVpcJMtHXmGAmdwcZ69yegbQF19G378pjwpEjPI/C7w1oNR/QBn5tCVtkgx3IRQCkIe/Y9mnye+aqY1loOSzeRGR6VGRpmq3y2ITIoknX3TDAKAJrXRGCxV7rhCDGoWTXtXrkyXNOkSTC2jX43f1EjRAw0C89ml07uATMZmwTU8hu6k4nVGKsaH0pp4q87rBQxVlacV53WFLFrbprfn9DgcwpOqnbixeF+g9P/vdlbgpIANp63W49cBsI5toS3VwcDMP7UdVn6XfVvUmoApTRfNLsHmNDSUXdqOvLY/nbjQS0vTW4xzWyfKPvHrBtUx1gnTqkwtnfaBTr9MIhv4izXZNPwG7z2/vRUDiWics3aYrZs05kTLqNZs/euvdlwKNj++8tA6Bk7I7EKBV9K801xh3DqQEfezWeNPE1+J0+daIEhsnaxsogDs9E33FhvWeiDIHNR8tutEsYytugAMi4RTv39VtSESZHFRlVDB4J0HURqJxQQb4C5LJ3te/2DhTeo43aGE2lAdJ8rW0Xtn9UzWjd6Mo52LI3NxaMopeCsec5AGtXz5sXVCDQhusA0IGRNoOz802LgGbsr+rjan/ukCEMs/5kdX8Z8B3aiU198MhmyP6UlFuNFKaX6y5TMbcT1zX4xMb4gZiPAqExS2OCCrYeFzGCYHzxk6QnUGCDxyQFkOLQYUewY7T/ebzf7jXwmIhj7vZS21AUliJx+0wvzchq211krW9a1ArcAh6Yn59ZegExeuxrlm82HSG115iJdceSLVGLwK7hqkYWifHdZ0p7K4hdu3EHob5skwmlsyTPDB0Iz3DW1/46G9FrZm4ZPQAOoNexQPS6ykuAr7OfTCbK/h35oE6z1WN4gPVLN1BgGg4xvpsYmRrwmbWdhg2mzOJP2mYSiFcfYJG1SbdfpQEEgXlB7iBdM78Tp1E/BQxgx5PzsQ1mDPAQIzG3Aw0fBG4hhGSJW3UAvMAOgLraUVvzC4j+a89Zrb6MyoTzuuC8LrioiwQ+qAm1JnFbqNSJaFPFfHwfRrBeU2fJ5vAxhmZiXGSV3cTiTEkHp4puXMPEMGPUgzohMsxO5A31TaHNTNFvImRjftu2sYd34m5st44RUseexvwi45pZQqdibwrANwM/eyebJ2zzJgVoy3eKN6zXd643cVc8lu0MvB4YTaEov2VfFg13RQqCLOUGkwTp1X7j1mCvxDtPsq3LNfidNI36KUszYJhhxSEzfxzEcSBBB6R0HNZourLbGnPy5UUFwJ4WPBLyr6z+dSCsee/BDwz8Hik7XJSM87JgVUdn9hBd1A3ajg1GhD/EDqMrTwT5OEDjuWC5dFF0AoIGfBs8GvLYGARs0b0VK+RBBuBo56z8l4q78TePvw/sWDaU2/E1MMyjapIECSUVgW+m92MBJEr6fho0diriwiY4DvpJbqxvBL7I/rxCtN5IgI8qgByIX4Zv5NSAj93AdcgH8NGmU1h7iegmgJ8EcAOCX/+Ymb+KiP4OgD8F4ALAvwLwnzDz+4/ldV+Dn1s7/cD05/F01fbo2Ah85QMnINlG11X1ziwLldZhpHpUl5pwnhfcyKuv9tjXjFtlwcPrGc7XBRdrRikEXtM26kgNAyIq8TuWy7fXeQdGOAJiZ30MABJBapZfBLx2bOJbFuvZ6jge3pLE+KgrAx4wgN5sktSyuqUV6MrbAfmM8dmxIEVI22jYr7Cq4xDwRTHdmX5o/2TgtzLaXi6TVUzE8DD39kLmEcBQsZfabnbWt2PIrRMkmaNPktk5gE9j5g8R0Q7AG4johwH8GIBXMPNKRF8H4BWQzdEOpvsa/IAB/Cw9mjrm7e+pEj0mUygnAT5eZCYXEExgBlYsbTc3hkZ2TrioK87K6qLwquB3a93hfF2wloRaMlBId1DDZgkWzNBgnR/wQdDtx/FoUxi8DgYb1nnkvg1b4i6/zW2zySsQmGm6XYY3lHdq0R3ZXjjXsbpg3BiBrzN4xMmjMjr3kvG9JuU2huefdQC+UlW/qQq/Tu8HZZeqPsjJQZBMD2wgSIBtYkWEwxPUo0qnMXiw7JX6If13px9m5h8Nl70JwJ+9LK/7G/wIqAt3/4/nj6ZRQhgGjk2W7Vz4HbJnguywVgi0iBgLNgAkrNxmvlqTM79dlogvgER73peMW+uCW/sF+31G3SdgJaQVoBUDAAalN6MFF4ibGEXEcB+Sy1On+wzoIAYH8jqYroCKIOFW5/DoY4wvZjOC0visQ2036vJmE9mkvOOxqTpkYMaXAl9kfkM+h+rBWd7I8D3WYBN1Hfhc52eoHWcRs4yzxF/k2kAw7r2ShO1xbQU5tX3iiivmnkFEPxX+f0g3LfNERBnAWwE8D8C3MPObhzy+EMD3Xvag+x78eGm/LTG2x7o0ilBxoJjFbNSvdTMxukFkIYpqhTgah2VTRbd5LEw4rySbGNWEi1yw5OphsJglisvFuuDiYkHZZ/A+gfYEWglpJaQiKwK2i/AFiWxBfBTX4BsO0WGmNtRp+z1b+sR9/rPqnbG+Q487BGpHjvWW3jnLG9so3rexMM/KOYB1B3oKfFH0jWDYRP1+4rBgA67nHOplfI8GgOyibirB0DEC34T5Iam1XKkck67W0cjblJoXQFv1Qd0qn1OlK4q972XmFx3PhwuAFxLR0wD8ABF9EjP/PAAQ0X8G8Tj7rssedDf38PgHAP4kgN9g5k/SY18N4K8A+E297CuZ+Yf03CsAvBSiEfkSZv6Ry57BJKsH2kNHVNsOvMgAeNQHRaYwMie9J8Z668AkyTnOQNWQ6OZOUDnB9vmxOH3rkrEsEsbKiI3t7rbuM+pFlo3NI/NbISKPgaCCnotTHoacQSpr8WwGiEAxS4O4G4+bsePQfV2dH2J7sQ1oy8ra//MHjeoIAFvGh8lvhOeNZR3eYVMHo+g+GDY6ttf5zQ3v6oy4f/C03AaAA+uLH1gwWZ20N+xPX0LKqr2hakEU5GwzKjFiBe+AWB93mMTae1pEZeb3E9HrAXwWgJ8nos+HYM6nq3h8NN1N5vcdAL4ZwHcOx7+Rmb8+HiCiTwTwOQBeANm97ceJ6OMv3cSIIHo2/Q3EDsvdda7H0YsMJOIyJzafsAlgtpm4dcgOALXTcwJoacxPAIqAmtQoQtgXQtlVlDXJSgm1xtWSUAuh7jOwT6C9REsmZX0GeqbsTquyAHdyxUa2GF1RbidtLJyTuukf1n9vGLg11Yx5hfJ6+xy4ZgoSE/CYWXOju8zGmTuUd2R6DnrxWOfCE7bSJG7MaWDSnb5vliLodUx2FHd74JM+MACf/ba17bqWzsVg1Qc2TwAFQy/vaeXeUwSKIaJnAtgr8D0A4DMAfB0RfRbEwPFHmfnhq+R1Nzcw+kkieu4VL38JgO9h5nMAv0xEbwfwYgBvvPQ5A/hJx9wyQLkYEN8qbop784eCsqXaOr72kb5DGvAV65BoYqABoEUDVlC1gWgbmXMl1CKhpqDgB0Ci765JgY+Q9hJCKu2F9aUVncjbolWr2MPc9zB1bbBE5lt2LPkA58BihvOH2mK8ZsaqDIwnN2+MFSHTKRheBnoT0O9AXMsxZXjBVWULhO23gR3HenMn7vZw8g41KUtXB+zfndgbJjgqIkb4BvIGfNb2k6VuANTqG8RgW1gMqHVX010QeQGcytr7IIBXq94vAXgNM79OMeMGgB/TzbbexMz/12MZPR46vy8mos8D8FMAvpyZ3wfgWRALjaV36rFNIqKXAXgZAOSnP20LfsCW9QGN+XVgRmADPW6d0bOK2QzAJ7oXAx/dNcvALzN8I+xxFUQl1KpRdgvrILNOS2KZ28snrSTMrygIqsgrH4b7gTGCoUML7JtXX2G6HZjIxuWFBgJ1qA+P7Hs8zdi0R/d7pr/rzod8hjxGXd8hA0cQAjrg674vY3t6swNfdHnxa3uEZ3AXSYXC83ySHd7FrfkGegZ23q+4fUdj16E2L/EhfSM14CMLQXnSxKCTgB8z/yyAT54cf97t5vVYg9+3AvjbkCb92wC+AWKZmdXKtAXV8vMQANx47rPZd+xCyGXoeJ6bsj77bcBn9oANfka2F5hWKixMLFheXfQlWcgv+hkB2KjDMvZVa9ItKNk7IlXpeLQS0gWBRsanoNuA15xbm67vjlIY9JH1dQYMHAa3Dc2KSAPl1uMlAdCmYuwRwPPjh+4disRWnmPgHRhczwDRJgHX94VrrB/aZDbUWcv/WAH693NVi71PXN3hQDhS20v+H4+P5+8C8Pkj717Wjyo9puDHzO+x30T0bQBep/++E8BzwqXPBvCuK2WaA/OL4Icwu8aBo5ZcdzT161uHHMUv64TJGJfp3tbe2djZQoG6FaBZWlvuMJppeyyY7klEGriOL+1JGd/A+sK6zm5v2Ziu4NoSWU5ke1GBH0U7L/4lyRn0Iaa3uUFviuBVcRzwwvGN8Wm87kg5ZyJt3Dt4ZHt2vmOCI/DNxMY4Qeu9HfmKv3ko+wh2Lt6281tRF3PgM3/QeP0hIDxlYshqpXsoPabgR0QPMvO79d8/A+Dn9fdrAXw3Eb0KYvB4PoC3XJ4hgmI5RFwZBqqBGYOF1ldjhrQZzCaCdIyito+7mgQDBIV1lQYWPmO7cilGIhEGymyRVdrNDqwKfAKG4VmRaQ4soHNxiEpsWPahDFNW0kBvo8yfMZljiQcQHM75t00uM+CbsLeNeHjo2PjIjXTQfm/0eA5q7VgnCneWXe5cXabSh/U/4hbRxWaG22SCB0Vb69cxjbP/CHwRSEcAvQtA+IQJbEBE/wjAp0KcFt8J4KsAfCoRvRDSpO8A8FcBgJnfRkSvAfALEB+dL7rU0usPCiwvdLgYA61FB6HWplEk00bpGETnRxdApzaLa2N/3GLoJWVzyrxsMCULROqDjlCZZa2lrcYwMWdt2yemzrm5PYt0VcdRdnNolUcAgH5gwwc3A50V06+xipp05GlZJsA3TjDG+kbgmzK8Q3nPSM6EWXVW2hH0DgQniK4rXV24rx+3dj0k8trEq1VHoUwHse+q+GOLde2338/9+XjdaBE2t6l4rsbZ9M7T3SSWjybdTWvv504Of/uR618J4JWP5lnNd8r+5+G8+ruNg2XCPMwQEvUskeVF0dfcTAyMALT9LTomSI1JmDsMAYlk7wYPbmngF9he8+9qIDzqtXzGTnK9R+jFAG5ADwLoz0XWAwM+Cz5q7NTHFndfzvRMXRAnF4TjXsf6bsEpfAN8EwC8LG2MVsO7bsT6rl0G4Jvp9qyvRYuuSR8OggPrIzQDGOw8Oe5f7cUmxxJ6HV2c7KZ+qvrcZP1oNmNwA70TohXjCcT8HsskrhPkI4WZNgB4+Ob26RTMB1if/d8ZHmzJUcwTEq0jJYALq/sLt42yteMmtJBN4LCGl8fnhzIe65Qd8PXfI/CNoGCgJ6yPOybkA31IbsV2sOPOsLEBs1DnFOq8A75B7D0KfAfY00HWF9luAL5pVJbZio0Z6PnxgfUZQ7b6USMKIURNRl+mKydlcZ2LSkzZXtpYXLsHQAPAeM3dTLeF9I9Nur/BL1Sou6zoSNmsbDBWb/cEBbuxvY3IOwzQ2XrLzg8rJor30AbIqKLN2gZYE8Aby9HKGJ5nUTuA1rldDDPgaxsdWfnidZ3Rw0BgYETTJnC/NMtwaB97iJW/9vXSNoTXYoXJBwjMb3zubBzR5PyM6Q6sj5XhXuq0HEVcYMv0lB27JMIyIXOl5gvVgeP83abJn0nCImHM7ggAemGsPIrC4nDaADCpyGBgmLSxr8XeezxZ+2x8qIaaHgEPxlpGcSy4YwRDxwiAUXxzYPJnKRut7ZpRlPO+GD/xfcZz8ZpjyZctaacf9XVAP/i639yudQAYWE14x5gZQx2obf4Z3nerR6WmUhhY3yjWR2OSPbJTOw5A3l2L8P4b4DvC9rrzgdFFdhfAL+qc3Wmd/fG6ly76T2yLK6TIFn1HQsWpzpI6m6hq2MHQK5SFwbrjs3ZK0//lWcikR5voiW3tvSup0NAx+8Xjnhz4BoYXgS+yk8j6fN1kvC74WA2bds/SIYulnVOS0K71Z/Hm+jGxBqy0+V8wndDt32Flo+EabNlf1HGhYz4cHwoTc2H61Ea8/T0d5LrJYwt804liqCuystqz0Mo96jU7IIzWWWNvg17P9kA+CHqR5WldbADPJ14rhjqzp1jmIwBwVWZk0V7jw0acmkRNEL03e4grXXALE59ltUiV9xtXC50iXTO/EyaGrk3U/8eZdexnLtoGANT/nSxG1sHbY5tB2U3vB4p5iHmF87edkjzQ2UAisBs+AqgdeOax4+7SYS87Y34RgSLwdfVHbbK4DPgGEPR8Yzm5AWA3wY0AOLI8e5+0Pb4RcZM5nqOBnvvxKeDZkkQDO2KdV/oOwEzBtSrodidrZq+sorYgo0nY3naOP8QAqakSEou6onADQQVCA0EAvR77ThPj2uBx0sTk++iydarYwTY9Aw6AUUe1ATvLIrLBMfnAI+2INkCoWXfHfWqDiBUHqeXnGBKKzjbLT5+tz9dZnVK4x88hgGH/Pc9v+O2VgG19gnuEDaDV6SqDZfcg8E1YuGdrVTA+n/q6PBgqf3RWjlbcAHq9wzL7MdmonvVb/newMwZo/1uZ2VhWkj6hjGtT56HPIVT3MZbEBI3Bx3Nxd+Lj6edtXXvVkPXBR9SBUEHP/U9Pla6Z34lTkQamzvEUYRAbXQj3RGYyS7cxC7vYE/ZIdV+/A9tEjkDYzdBHgUjfj8RhVmKzsYfPcj0QsO34xybdq7yvAaG8+DyPjlVTrzO9XeCzMmkVb1hsZG6EAGID6HXsD1sL7mWgl1gj7wi4EbFH4onHNtXBhFoJFRVck4bTMwbYVBSx/kapgobvvj30/jS0d1RzTP08hYnK/isB+JKI6PZbX2Jy/52ka+Z3usTiBMwEoJJ3ZmIKVjnaAl/4PsrqJsfdbpKo2/HK+4kCX13sg7Bhd9vAuw2+NrCtuLPP1j2FvPDCBKBELNDGNGF5AGzJnegTB+Z2rL9HxaT97/VIjfFFK7VZu8v2eAd8FZvB779peGSsGzsWgS/1dbsJQOAgdwXQS7LrngFemgDfWMMMid4joJgAZhSKu6JNbhjqvwFfpIVa/4kkSMHA8CLodf1gdHxO0l99348EcK0A5SYS27WnTKc1Ht9xur/BD3CxV1wJACjwWQghp+5R4gg9sPO3Gq6L7h8N9NBWcLCAm63ucAaSCTVLoFUDwLrIsbiBt/uXGTMNjKkT5zbAF8toG8309HY2yOI64wg0oz6zt0Qbk7NyUihrK2/vsEyDng8b/8UrAd/YLsNkEEHPgc/ql9CB4Qb0ohHHGaCCXO5BLyVGSlX9NqtqErhjfXF1GTOhALqdqYCJrP3evtyM7V1aF4gTJs1BrzN0DRMcQhv76o+0YYNy3YkA8Jik9Til+xv8bOAQROdFrETB1vHKQGQDxpgGkKPJwIrHmu+b5MemiKIWmKCJX9QD3w7gRT6RCbpjrU7osrE0pBMO4lo38NNYRrp8lvYBxR2A9fpO3ZPYf0PqkaGGpfCM4BdJGopLQJC6VTHu1H1M1B3LOLRTBP9NfQy+elGlsNHp2WRo7TVheymbaFuREiNn+fa9l1NFUtaXaCvyykb20hfXksGJUatdZ3VoL6SvyAPoHaoLwEHN9bt2+SHQi3pgtGdwsKp3IFiVHJhUcULyd+3nd8Jk/d986jj1g7exFQxgx62Bo8HAACUaJgLg+UDL5jCvRgb3mm+6vQh8dSH53inwLUBdGiMBWlGj1/7ogrIZ/Cb+HKsjDqHsuR9onZ6t87WTCNRcrd5CXVkyw5GB3ejAPTh32zOiG9HMuNE1bqyDTqcX22b4OBscluZ1ris965sBX84CdDlXZAW9nBg5VQE9oNuIHoBvTVqUdsom9gBRcutwfLdNihPReMr1u1ofYyYR+MwIN6zssXzAcF1xXBxgRjOuQSc5UZs86nQNfidObMwvzGKHEg2/baAwfANuaAAC6liFfMiBT877zmmasbM+BzjIJuX2OwKfDcxYpgRnsqOur7EdQ/PjKTJJaCh9dyMZAM+XOSlQsX/rhABGt4rE6t3iD5pOr1AAQGOC7VkbpnmkfXiog849ZQN83IAvIzA+DuwwAF8AwKjfQ2B8Efhso6mcFAiJfeOpEfxKTUjaaMyEUkVUFp3fZEMopqs0Z3e5r2SKk/oM+HwC6R/qkHcZAJ4arK7F3ruXNmAyAkc47npBdddg7QgGcFwFrFKF72xfmcVbgMIAj3nq4BORllzcFfYXgG+JA1qBSUX4bkXGZUxhWgnhem5528SwYX0bvRyBMstewdSGhVSTjIYWlToCXliXbKLuZp9h9EzPRLADbdbpPDfMDg3Uwv9m2OhDcw2sa+wPAILU2D6AMr4GfDlVLKkigV30BQT4EggpM9aaBAhdNL6N9psk00v7pGwRWlzXiynwbTYiugzMGM0IchdY2t3I807SfQ9+vQ6IO7HHdTxxxge0EzTDCDPLQnCdTFFlF3ty4BPylECoxC6+dSxTn1vNomtsT3V8AoRBL5UbCyAO/lrUOrv9f/XKiN8Wnp/d54xUHdDWKFP71mALlAW0QAJixgbEn408/07cjYwv6vgGEddF60PvNjDdKPJ2DNAmKBNns52Lsfa0/TvgCw8fVSFBqjdDhoi18n8Evky16QAD+FUmoCY/nlNFSQRSV5eDbXmVNh5YsagOt/q9zsp7lfy9X3APfKde4cHU+zneA+m+Bj8Xc6yjR+X26MZgA8JvhltVPc5fVhDKECud6rVYlw7JRldxY2zWclB7vjG/JXwvLN9JATDo+gAdeODb7xyToJReJpvqq4g3EmrLdpMLTG10TSkQn8VkWiXVo0afMQO0mMco9o6W3XEcDYN51PF1YKfXb3z5gmV3XMI2Y3cwI5VnyM6GbeKTquxvjEBnwGfsL16zIgGpolba6APNPaYrl/dJCv/3H2N7bPcnuKgaMh9WcwysD+214wZJsCC8FhsyboRUW186WbpmfqdNvHDr8KNF7zLwIwgA6imGLhFjYW+VXQ3n4q5ugKUdiLqBNrqydGKusT1jLIGF2mbWU/GIh9+21ncQHaPIAgiYck7wUFzpgPtJtL5G4IqDydxZVEQHAoubMb4Z8A0df8bwNqLujP0FUdddWrSNWz42KejE5iNfH2wNro3gRSOAkcCpiIOyNoixOls+a6CWXFkKVBbGJ+Jv7wIzdYIewH/8379NEkjkBih3bN/kOWF94QUPAp9tfj4BwJOma/A7YUoA7wbxxsAuRCPxZUluvtdZPyqFFx0byg5sEbixgAoBhGpoOJYjKNyF8XFwbmb/Np2eJzbmN0kR4JShkHdk+2hIrdhxAUHSUsGUxCiR9DrTVTIOi74VSEUCrTZXCHjBSctk4JYKHV69MQDfZWLtBgw7kbcxfKR4P/f3DnXYducLIAiSkGLG+lhmQiagFhFTa2KUSsiJHADNlUWALnUAeNuJtvUxA3v3QrBzoLYK41DW44ofe3XtI7751bj5eVHx1wHwhIh1DX4nTMTgXQ1iwgTwbDE6NecANiNHFdGWKQGrAJUMBgqsj5GMIZgldoig0blgKNDVAfiQOegfvSACRIMbSacn4/5/sICMdWDofq4ehih451NK8Ii+FnMwUweCHdvL7TdIdJwer28A7FF07mIgHmF8M3CbDfh2HR8Ex27lRqxXK+Pwu1tSZnkwhEUlBiOhogJL0kX/4qJCRXWANWElAYfqVgihyBUNICsaSFYH1tC4/k2h74YyDUAoH4Itj9sAG9DE3EZoO+DzTcsD8MkGWHXzbZFeTqvzw2RmenzT3dzD4zkAvhPAx0Kq/SFm/iYiejqA7wXwXMg+Hn9O9+4FEb0CwEshc/KXMPOPHH8IQDeCQmkAu6SoYUuRLDGTGkKTSg8CcUAYM9ZpFzlLCaBiy4JaEcwS6Q7LKvYa83PQM8tuAD4gbGpjrxRBzoBkA4Q8gGOYvaMeEhWEJIOmQPYLqTy4pNhHjyvT8ECrnZjYjnUib3BnmcUu9HoagGuzLtd/t0miA0qE38YCR9Dzh7XGHK3L1lecTaU28TErACJjBXz5WiLV/dWmD5j5+TXQI+1nOnF2smxfJ43ZogNy9xu1ydt0LqRlHQCwZdz/6yHZTLQ14FtrY3ulCOjVKmGtTg1+6LvQvZDuJvNbIZuS/zQRPRXAW4noxwB8AYCfYOavJaKXA3g5gK8gok8E8DkAXgDZwe3Hiejjj25kRIx8VmDRNQDoEqS29OhQxI2qm4dTSaiUUCkAoA7cCgG8BGV9SZcphSJsnG+DiOtibgyLBLRZMO6zwMMnHOsBRplYsZlbv01sMcOHykqMKkwuZRWPde/gjejbg2FfVvTUz6p7ZHvxezYmB+Dr6m3Q5W3Azpg9wv92PnzTWH/x2MA8SJkUmyoj23I0yOZSAPbDO9RcsFNwG40exv72JaNEMbnSHEuU/fUgqCuIkkgEpkqhGsFSlzQeA0AMbI8xBT5aK7CWBnqlwCI9nzqS830p9hLRTwD4Bmb+oXDsIWZ+2aF7dIvKd+vvDxLRLwJ4FoCXAPhUvezVAF4P4Cv0+Pcw8zmAXyaitwN4MYA3HnpGSowbN/cObHGR+bj2MqZaE2oSh9SaGOvKACXtJwlm+fPrEZye4+CyQXgM9HTtbtTRsALHVAwIjK+5pAR217FBc0+oDnwWTt8MQIQkg0g7/Sj6NtcXuFOygR9TeNYgUpqqa1yf6yqwsY4mwNetb44GjMB2Nnq8Q5LTDPDiRlQ293TtR411WcCJLH2AK4tGYWdtRk27kAmZ1aE59DFjfUUBcC0JxdifMeihO0ZGbO9+iAm6VdeMXrFaZgB4VeArRYBPQRCliLqjnHYH8/uV+f1eCDv7d5n5a/TYi676ECJ6LoBPBvBmAB9je/cy87uJ6KP1smcBeFO47Z16bMzrZQBeBgBnH/0ReNKNi41LAdAscsOYBQCUmvRTsaoYY+0s+DKIwCoGMofB5J2VnfHBjRrswDcztogRgVqnDgXsxMU4cKNuznQ2EdAGS50zP7JzFUgU7u3ZXmSBKPJuCWHwje1wiO0Nk4O/2gh+bhlvLivRPSk+dzZHbETZ7vnUAXLPnts2BVGnyEVFzAWiB94ZAxTQ4pqUxSXUWrrlbgDc0ssKiqUS1ppQinyqg+D4IuE9g9tOE3fbkknbb5qgjG+cpS35BDQBvtpEXWN3DnzrKtFdTATmUzO/QzPX45OuCn7vB/DpAP4bIvqnAP7iVR9ARE8B8H0AvpSZP0B0sAJmJzZNy8wPAXgIAJ76CR/LT71x7uciCKbh1qrZ2xKkVT/GlCSwI4BFVnLY+mBzpTGjAbqBY4M2sL2soJerhkRCN+URA7WQsxp7bf9Xmdh832B0gEc1dmD7rQBIBGRT3CuzclcHNBE46uv0m9Xi64QvgvQAcqNe8mDLRkAzoMsN+HwC8esYx3oEgzoA7JjeaHiZbY9pt5pzsD6/VgIrAUJNYGbUCvCuCrDlgrqQBDwg7iQM0+sVXeFRK6GUJGBqwR9mABjrmGLfsriNcKuvGT7Gd2hbKtjB6KyM1i/c6CGGDVI9H2oV4FtXYX52/FTpWP94nNJVwY+YeQXw14joCwC8AcBHXXoT0Q4CfN/FzN+vh99DRA8q63sQwG/o8XcCeE64/dkA3nUs/0wVT92db46P/lf+GzobU0JOCankjS5wZQIvgCxMl+HvfmIJveUuOlObiJu5hURS3WOkGhUQQKw0Qfae+fXirzGzwPxKYH7RSicvA0BZQoUOvOR6QhuEDoQOEMr+gipyM18dYKgbQ5B+O5sJrKZm7kDPxWAF6qNTYfz2OgrvY/sdxziC0fdwYKdNFJfr6yITg2CC+NdxJZRaxXezJl3720JddcVkVa0o6NWahFmqS5DrH80Y0gHe/NvAUKz48s04wv5CosgAAwBqQaXflOIfdt3ftc4PAP6e/WDm7yCinwPwRcduIKF43w7gF5n5VeHUawF8PoCv1e8fDMe/m4heBTF4PB/AW449IxPjKRPwi6lqB3NXhERYa8ZFzTJzl+a3YtY5ZkLNVfFDep/vThbdC6IVNwDfGBap5S9rhWtwW2gn0Q1K18E5OPXA1ymlmZulzkOQk5xL1HR5pcqyPQcrY4LUPb8zeABbIDoEenzg+lHnF3V8nf6Pt2tx/Xm9DN2Ls9QzZbNmj+uLSwO/DqijKK4rVbgAtJPra03uFlULgRZZ/lgLvJ3dbmVzTxXZgysp8CVnfjbRdHUWXzcYOxoDhEeGseNX0qGN14RlbA6AHMPZt37FJ/bzuxOXyLuRrgR+zPz3h//fCuALL7ntUwD8JQA/R0Q/o8e+EgJ6ryGilwL4VQCfrXm+jYheA+AXIJbiLzpq6YUwvCcvc/DrGJ/6XhUW4BNxpXbrMjkTaq36TWDWqLbZZtbm8wb5Vz8KfAke+ZeIkTQOXM8sNQLVCHwxdazKjB1wEISt2CjVQc0Yn7kyAE18lUi91sm3jK1jmrUXez1NwG8s6zQNolxvFUcn+m6WJE7qpXtmBLrA9tIYVKG0T3JGOIqNCjC6KocKUAsErJQJ1wLUHQGFwCuDFwo63SFkFcP1ey7uVrmXAvtzVQO3ckQ9Z9T9eWShCkTDh7E/UQNwW5IU2+Aq+KV9iyMYVsZJ9X73KfO77cTMb8C8GwOiP5zd80oAr7zqMwjAjbQevaZycuCrnLCjivO69AvSs1rpckVhQs4irnAW867jaNwpDuhENAe+ZExgbm0mjdVngVflxeFspmNUQXfV6wFZj0crb20zuT2rKpiMx80oEkTFXodnWxmid7EI5d3o2+xfu476QdvEyxB+aoivNwW+ULZ+L5DJuuIAeqkAtDbQE/VAqMsBdNzwsgK0kLO+jgUWoO4go6ZkFde5X1q3KbeBnQJfWAN9zPobXV9c/zn5X7BtBoAmrbR7SDogppuc69ptook65gSJhn5yL6T7eoUHgbGjo+QQJYIcJ+wNyeqCmghnXAT4UkLJBaUSakoCYmydW8VC6UmzgrQOaYdUTOkkW8bU9WaT3QBIzgK5AdfRTyxE/N+vsf8HlsnBidu2W9TX24pPrax+KLC8Y6wPHRiiibqzxBCRdwS+wOyc7a3h2BoY38oNBIPudCy3ReQRoNQ8dyIFUiHQTiYMXgl1YYkENFu5M5Y/RLnuwHqcdCIYD3pAc3j2fWPsYYV7AByDlSZ5X+nDCVyr6JxzUrWI0Eki3XVQ5Gv9H0Bw6r7jdJ9ae+/blNtmpQAqMiQs1S4J6K0q/uZUkXTTGduoxnQ6bd2wyhBjI9oMH/WBKvo0F5crFHYDeCMQcjjPOqAC4N2GR37EmhFoO6aVDksrU8YHHGZ9NpgxsL4jddH0eUPA1ALQOrA9Y3r+zQEAgVS4scVho3lOJG2egbQSys6uJdX7wcGrLgqCpic0YJroKgE0hh0B0HWQgQV2lRvqMU4grIDGIdKOTVYc/tdZi7jPmjQjTgClBM4Mqkl2QUwSzYNylttLEYA/VbpmfqdLDMJe402lYRDloPgoId5Pgur7uIUdWlLFvlrMtt59AR5SyB4aqRw5IJlrDOl3I2IRECfLnTQdJIQKdE0/FCx2V0lHImnG3du2QCvPZW4iUizyQVF3NGR0rG9Yh9vlp1Z1bv9P2VIVhpcKuVi7BT7eMD8XfauAoIORvYOuqkgahowqUGLQBwUsY4G1QMKTaUSe6fpiYAuANqmET3NPGfpBx/paXEbp8tTUExWwvYGdPWp/8Qm49JkLm9dglTkDizHIdp5TOqnF9xRiLxHdBPCTAG5A8OsfM/NXHVs2eyjd5+AHnNf2CllrN4GxIm8ibpj+L6YU7rn0YUFcBCA9JAQhBVk/TrZEHgjAO0qdLb9hxo4gdDvJxFwNRSxWw9nniAQyY4D6fseYHjAAwAQEo7i7eSYAcloIZ32H2F4yEXc1nZ795u5YciBkOadO4s0BWupLyicA2BgfBwAMLFABTPTCcZOkub7M1QYGuNzychecoa5tb2aKE4peoyu23TUnyss211IAQ1maSUga5YeIvFca4JmfoAOgx348EfgNk80dpHMAn8bMH1JXujcQ0Q8D+I8xWTZ7LKP7GvwKJ3x4veH/J9I9FrQXmVV3TGvNXRSOoymwOveNC7HURFdlbA5AbgAomyn1bhBckwdJtfwfVbLem0gNMQbC4ZsItp0h6wDvietEfB/K1G2AZO8xEXE737Tcfsdw8/H5PZMOM4Kx54nLSgpi7ijiCrhFMOTt98pqEWa0iMX6wgTZV1l1ZBQckpsrELdji3oZWeRuFXs3wG7Zd5JBBMHwPSatUzFckfudAoyqQXVFBG5tZvmY7tYlksrglMT1ilgC2Wj7OwBqf6GUgJLE3++U6QTMj6USPqT/7vTDOLxs9mC6r8GvMuED+5sAgjhLYZtBsANgDrTFLL9jGKKYos9fXJfpymvrbCp22MxvbgKVq4SI8oBYLV8wFACDKHSQMoR7nRGQdlQ48JllWjq5AV+STyZZtJ/TwAi3zzAdkZVpo96MIm4QXzeuLN2HHQQ7PR8bK+qZ7+ir56A3sjxndAEgSzBumJi79g7hzWLOTeRMyuQ0dp8zI28XR2sAFFgYHIQ4kcfcaxXWvkcDh+t1ByCM9WvMD4CuMZeli6T6WDA3cB7LzNR8OXVvFi7qfF8kc5kjtT/lBOwV/NYEWgw9T5SultUziOinwv8P6aouT0SUAbwVwPMAfAszv5mIDi2bPZjua/ArNeED5zfdqBA3lll0i8GFyma/BUCAc60JK+e24YwBHtDEU2UhonsiYJ0oqAkCfFmXxVVjeM39pfMBC/vcjh3X8osgA4304YAVwI2TsBDz9RMGIjM4UgIv+knKAnPTUblObgS0WYpghwnbC2AX2V+3hG3MPzIWtjqG++11Iq59GwCuc9AznZ4AJzsrpCrMz53ELZKxLwsTuZ6yvExFI9fNim0TBqvoqecTHLymLi/2fiNAGQgO/n5dnaf+PqjbSmN8tL2vey7BdcaVkFbxXaWiOs5EoFyRMjnw8VqAJTf3qROlK+r83svMR+MGqP/vC4noaQB+gIg+6dGU5/4GPyZ84NYNsc4SQNT2VR03nDEwNJHYWN9aEwonCUOkC9Lj0iRW51RUCPBFPy1LJOyGCzX3B3MaVSW/RXVxEfrQOs8BiBpQtVBHrCAGFlbJqzJQV/gIO0QmAb5sIBjuDaxtA3rx2TE0+sj2gi5vY9wYghaMoBD9BN3HUZX/AnYUgE6ttWsAvXVgdx7wQb/DOV8K2EUtVkakdeYsC9KP0iorcYxAcwKS6QSFQAN6j4GYi/pjk44gP2GDnaHJ61qADrnVvUTaDhIDwj2TZEYtN9gsYixK+9Y3UyYVdwm0VNCaWuCDUzK/Eydmfj8RvR7AZ+HwstmD6b4Gv1oTPvzIDfenS7bnqgFfblsOLrr94AwA9yVjr5FeLBIH63pMAT4CrdQGZVgmBcA7qg38uqCB4MI6U1PrrWbJ5OGDmF/PpgwA6yIMlKP/lQHfDPxSA8CaFQCz5ONAGgGwY5yhPBHoJoyvY3/jzmpo90WRt9PvsbG7XsRNewM/bsxvRWe86NY+FyAG75R8jwBf2PPEHH3FMhru9/yDmwqzM1V7HUKos/Cqnahr/4/n0V/TRN44qck1GzY4po550vAOsnolZQFBzrpiJRPSQqBVHOZpteVtR55zu+kEOEpEzwSwV+B7AMBnAPg6HF42ezDd1+DHlbB/ZCfiikZQsbWWObNE3lD2Z0CYiFVyVABUcXe1+GsloRRd4laSL0milQZLI/mgdfBQy5/5gNUFAp4W428MCDrO4BGEOrYH6aTqWlEzAbvU9DWpwmV1QMAvq5FD9X0GfOLIq7qtGFElRA1uIip1rO7Q7w3bG0VgNQZ0arOB7TnwjaC3b8CX1uiwHMAvAp7Wg4AqO1B1Fl4MwGdt6G2jUFbHE7H8oR31f45tiXbNQXaHcH7Wv1W0Zda1vcxNFTNL1mbjYdP7qS41FaAuIgIn9WvMGagraf0m0NrWiZ8kWVvfeXoQwKtV75cAvIaZX0dEb8Rk2eyxdF+DHyrAj2RtdIYFEKXMKLlKkIFUsSzCCH3zaRVBLQxRqdRi/JUkC9lXAT4UY33K/JwBYjNTC+gAvDZHWNvLA5ngmyqF8ncbn6NnfAYoNYuymhaVy8T1HomSAKuvFojMD878HPByYH3BObc56SIYJdAYaAS7CbC1sFS2TjeyQ26Ajgg+8LXJY72mvX1YwQ/I+8FlZW2srwO7DbviAFTzwbzBNwOQNAE+YCra+jMDwM9A71IXJuq/ufu91e8dMkBtMg3suqqRqGq/znvZ6S+tDN7bpKIRgE4p9Z4gL2b+WUhs0PH4b+HAstlD6T4HP0J6JLWB6vtoiFkfC6PmKo2dqgcbkE+bhizU+Lpm1JJQ1wReG/DRXhnJnkD7poeKQQciENQMd4SlhQV8FtXXWdgmANGKHP3jLK+aAXKnW+oHDiWwWS8zXHcl2VIDrsD0JICoAbKy0+CqwTMmGI6NbM+BbtxDV99xqgPj9t5u2CjUmF0EvosGemnPygAj4wvi65GBNXMKZ6J+e04gWMnhrkFN90Y9uASw4Qg89o76vXFtmYGfgRuFvoAhT72OZ/cMwNcdQw/CVt9Vnb9loibkvTDAtEwY9QlS1HjcK+m+Bj9iYHmEAjNRMdHARnUZZWHUnESZmywGW4vlZxsa1arAt1dT/56Q9jowFfiyDs62rCqAjgIMLZDF8KtGAtkJw+QECY2UaDtQYKCi5a/Q2HHUmA2ipZbcerkN0Ik2EFLP8mqmBnAOgOiZXbwmAF81xXsExc7JFy72uo4vtlcAIBN540oNB74LYXz5grW+2cGP1qrWXAW0uKrBWZAxt7HD+PJ/AG0dLMJ9rgf1/hTrYvs7Wrw3Vl4Ds0vEXruWwjuMLG5mad8AXSjLzOocgzqkFbJGWVUJvJD2a+tX6tx90Pz/KNI1+J0wVSDfiuDH7nRaV/W+X8S3iTMLm8sKhGGXNzB5wEpek1h1jfHtyQdl3gPpoumhOsOHit4CfiShkXY2yCUwJi+65aGxv2FWjwBUfRtNgFijS5tobdcYIzQFfEjO/pypTQZxx/SU+XQgOJ7H3JIbQLCLzhIBnuGg5CBgDswu+jZRN1+wgh8LGK6MtK+ikFdLuuvtTDwNvotMEIu6AaB963K0qI/z8Z0D+C0pqAgobEXa6sYnBAN6zdebdTbY+UBbRRAOfeIQszPfvw6AbQxYXwr5RQZKRcfIMoDgolZuN+rRyZhfp9u+R9J9DX5UgXweOgO1RqUss1ktJOC3qOXVgFD1Ue4hwtTp+EbWlwz4TBQz9hf0TM0iy6g7uEHDPm1DdAOgYRAAai1mJG7rgwtkgCTVGdJKvqOXLY+K/l6jHqgNkJ7RdRbacG4DgIOYWz0G34z9qeoh9HRbp8tGeH0gkq64QLPo7qV+s36ncwG9bMBXQsh+SwVovo26SkNZstetsUFXEVCrcz3fJgXy2H5tD2YBwJrhOtMawLC9bHs/Z3QBGPvGtltmLDS0YWBzHbPrJibeGqSirNlZ1tkNTDMQdINSjZ3pBOnEC0buNN3X4AcG0jn6TqCAUBcCF4jIqfsycFE2aPttELe2Nctr2QIfBV1UvmgDNKnSPbq8iL4PLjZEUJLFGA0ACQ0AOVlQIn01liVMQBs8NjjFUx99hJBRnBpEJB8UAfBinj0QDnV6QMy13eqc7WU0w1NspsoY927p3C90RUaz6BoQBuDbV6R9he9XMoTuMmBlq7AR+GZAOJxngjI+6T9VmZ98FAwj8C09IFmfdB1fUjx7lIwn6hKn7efSzth2Td/aAFD7ote5jAVSa6+1JykQ1uBwfqp0zfxOmIgFkLqBa9ZVViBkYWKsom1doOJQEJcg19t60hQMHA30AvCZOLaq1dEcZRMhZUYthKKzpoWrt/JR/A6gKfezD1wBSvmfVrkn6btZgE33ZZuBX8x3ox+iwSKLOQDSBPSMCfo2ncAmCvOogIwicKxrs0AGdxdbgpb2LA7Nq37va7/JtjsnS1tSSupqQj3BcmAz0bivItexZQqSQ7OKG7urC/x/B8JQHx34DbrNKLKaU3RfiMb+Nm43BsoD8EXVQ5uoGhBOtwTgpuJxnXUSEExZmB8VuFQx6pLvOF2D3wkTC0NonaKJiuZ3VRkiQlYWel+pAeTQaX1NabTsOusbdFEXOihtQyAAULcSj69KQFuO1AAwJUK12Rnj4OFwH1DNf7GQrOTIwpKqiou2Y9lmyd2oQ+rAr737lA0aGB7S7xkYKvBFJr3R9U2UWZGMtAjVzVk5BiIQMFTGZ8AXRF4CZBUOzZRmVr+N2Vk7tXNw5igqh94dyFmfibgKfO7CZKAfM4zvn/Q9gxh8bKPxY8nbzeo6MnRfRqgM3CaiUbViwRnCJEyr/E+65pc96CtOB1izyflxTncN/IjoOQC+E8DHQubCh5j5m4joqwH8FQC/qZd+pW2GTkSvAPBSiBbnS5j5R44+A6qrIGNU7PusJoavzxQQ1DBAuuKA1RHZ8wpsxAEvWCBNGZ/3rAxQFfC2exqEPVAm3z/EQKRmAuneEGQrhmJniEAEKTTbesxE4spiLi2lzdyV+07a+ZUdSDOr4WghvEysuhT4Oq9tYAPK8TejA+8u0CjDjRvNshtE3ii+6u8O4E0c9f+DHnC8fnAFqjm4AQUQrJHxLb146b6WFN5dq4NrK8tUB2j1McfvdkmYnJoxy1QR3PaODn6vM9Yt6onGyk0lQ9mWBtJcmriD9EQSe1cAX87MP01ETwXwViL6MT33jcz89fFiIvpEAJ8D4AWQ3dt+nIg+/ugmRjZoqPWbBFl0DlKGBAAKPH4hS8ewe2Ne4nJB7nphivhsLNDEsIuq4Bd2TCsEia4CMCVXKKdFHEvjvhOjfs4jn9gpHeOkOkRks8BxY6iWVxB7o5PvVUBw4yR7SIkegM8HVdR1RZbhg22U4dph14VGAAzllf/Zf3dpXMWQjN0RPOhDdFshATW7N76zMbGN03ew7rpTuE4CIvYO65bju8V3N6xOE/YHnoPLoclilqj/OAv1Nhq8GjhcD8g2rQbWxGL003xODX5PGOan4WUsxMwHiegXATzryC0vAfA9zHwO4JeJ6O0AXgzgjccfBLcgsg8kdpcGIgg4dOOF2gxtMzADLXgmPGJIW1rFnRJe3C4q0lrd8kgpiZidAN4ri1jQlhQdAj4HnQaA9m6oMujabmVoxpSwXtOU2Q56Khke7cCHGODEesiBRZhI1cL7B3YRRT4d/N5OdszKdaxZDZSMVSW0NvVrSAd5QvPzbGDn95s+bwA9y8NUIB3wpcYAoy5584mTiNV71ufkxvooTdjfIeCjVl+M9n1ZfbU2YG8Tiu1DDdBi3GZvItKQWYVaXNxLnns76UTL206WHhOdHxE9F7Ik5c2QLS2/mIg+D8BPQdjh+yDA+KZw2zsxAUsiehmAlwHA2ZM/qp0IDILDb4TOh9j5eHA14QZ8MYBmv/UhOyBSUeAzJTxkkOgYFb1eYXHirey6uVH3wcqYWsy70ONYgYa1J9Ym8jbgo+a+EERHA8NuvelsBEUADGLU6DZhlnG3bMbBFtOMtRz6jpcagJCxOBUpNQ4hsbAUWZJoDFMAzv3zsobucn89OBPswcGe2cQ+czyPYmUD1DAJmA40smPLT8uEAte3uO/kyP6wZX8WOHZcK+xtGD481OVRgKTxt4FhfH/2yYntfJy87jSdmkWeIKXLL7mzRERPAfB9AL6UmT8A4FsB/H4AL4Qww2+wSye3b6qLmR9i5hcx84uWm08OD0KbgWcpMhAHCQof9J8AhLaigFbIgm9bYbBWUCkSBaSw/NaNw/vPAQYWQcSATzc/x1KBhYEdAzv5zbsK3jHqjQo+Y9T42YlvYV2AuoM4VO/Et7FT2I9ibBzoHasb6vQQ2MX65YAe7uNIjREZITlSD028NudiCcdVY2iuXZbPkoAlxCzUEF41p5ZHiFrTibnBCdxfYXjnaGDYuAHRAILBymqGELuv06fGfGf90+uy76vjZwTD/tMGAo/9blb3oQ28/e1d0njDo0t0xc9jme4q89MY+98H4LuY+fsBgJnfE85/G4DX6b/vBPCccPuzAbzrsmeMnQuHOvDRTMIHCAO1Z1O+faTHh+t9zhiyGJwzb/Vw8XdXSfbhpqNxNmEUAXB9jQZCNUbIGnKdS9trgoOPFqsBgQ3sR/YVmV9ntDAmA2d90wFrDGHUX9ksZOyUDyjQB7ZloGWrV4RTJ1ElGIuOZR+COFi8Qg5srmN6sd7Rjkerb88IexB0RprasU6chjI4zYOo7QvCQc/sKpkBmHwBnk4YTO0eEFx6Ydtjw9QBulWl7LbHbcJJKi1Y+1o/sknpmDSAbZPfUbrHmN/dtPYSgG8H8IvM/Kpw/EELNw3gzwD4ef39WgDfTUSvghg8ng/gLccfEmbubnZujruxE1/KDtHAaQpYfoy3HwBUq+yWZdtJ1gCeBzqZDDb2gUQalcY7qyXttBYFmhkea5Ar1LzdQBDqBmPPoRoGhSrOOp+y4PB9lD3bgGWAIOArI4S21yrgmoi+mWTsNoIyPVbLOGudkY7t5uazCWJgTE5F3JHpHfLvu0oajUG9aB6BkVu9WFm0zTkJ8BjYuc5yKI8vxIACoAKZ+wzGvmuqD4LvL+Sg6kY8ewCjE/u9kdEA0P6/y+mJZO39FAB/CcDPEdHP6LGvBPC5RPRCSHW/A8BfBQBmfhsRvQbAL0AsxV901NKrGcR1lRH4ulUJnSIf/WyN0A9D59zO5uPDwwHz88vUgHACeHFZVX9CPhbynjQg6yYD1sHNAnicqEWaLtQGmlm2SVeNVMtfmSCUAVwVEJhgseREt0TqZ4O2ciPmFQdUp1LomefI+GqGWOatYCSMhgoj5RTYdKgXmoCdje+JUudKSnfCBvj8HQ0vorEnis5JgwGQMjjaAuAxIPZJNykADucdDwmyFNPKUXQ2K3KHrSGSvCzjMLv5RKSTkx9r504OVk8U8GPmN2DezD905J5XAnjllR9CIdJIDLwZrXNTIBzEODUQuI7FwVJWaHAcAP7sKyLHZZ2os5iKQ3MiBUHiPv4ps0egYVK9VmX5TQCKMQdVWjOJGO2IjvZ9SZH8Fh6+XYQ6Dp4dcw6GmRnz8/apBGb2tcxsIqMuV/QADrMHH2J5zuTlh4uet5NiX+mAdjhuzxuuGZnelP2NE2VFA9DuHUPXIyAVQo0AGIVVZZ5T+TWqJRQAo1HuqJHs0SR+FPV+l9P9vcIDgO1HERX3Fmmj5nZstNJ1nYGgUUICa4xictg/w/RM3edQ0MuuoMP3mHzAKAiS7rIVAFAkbGF+lapEdTZfHgvOb6KWRTQJg1EJSXukiVokfxQypwBJ7jrE2Az4IT9/Tx4GkzFAbvezMdJsBLqxGU79agO3Yo/sD9hORtFHMDAaYj4OgEO/4Nn/EzYYf3fgpm0qlt12rV8T6m1Tf6Q/uN0YJWEUzcdYH4z52jpx69M8qZ/A9mK9jsa5K4sHV0hPFOb3mCRC87UarZgj+4shl2JnNZaTgAqNplLR+3KZOJ25syCSWg198yADw8uK3XXyORUjBZo0RJ0WAJRAAQUJqQLV9nQdAW/23AhKQ3l8yLBXTYeDIjm1jDePGADexm40+mwihdikxZA9jwmy3jSJlZ2NkUTgHNn05F3AFPSMLIPaCmWM6FiaANoIhA1wDGgbokVw2zC9Y20zFsF1eK1R1cYBELqllO2Z2h9BQQyfvKO3CW3ayV72lKLvE0nnd9cTkwYqsI6Zghhs3vlj+CWC6FO0M3hjq3GgQq2nBoAbQ4paE3MC57odzEcLPMzyESziAKbxm5VwsoiCLPEBU2IUhjJFG4STskQRlCfPtudAAbAqO4qnD4nMI+CNDLd738HXMTAp36FMfTJJrb5xuVvnMhQH5yCmuXispJjcGiqV2q2sCNU1GnxMjG46RT7M+iD1f1nwz0M41C4Ik5HNpRROQideNXj4edP/+ZX2nhTqOTTeUH8j8N0VoLoGv9Omzjvfvp218Zz5hWVY5gLClWHbi1dbMrYAtEikljHEexR/m5GE5uyPZVD0VmOCRSH2PqHHNxLd8NsImG9gjSP9isfPMMvrNT4ZIADgZvDN8+7Y1gh+s3sQ7gnMzxi4rWqJou4W/OCW355VwnWMzp6h7MbAy4AhViqFstj3qP6IIOLHOFA8yy+A1CG2d4CJxYnJ5rKRRTrzK33+vf6vvWdrx7FfDs8age/EYHXN/E6ZjPmh76wOeh5wkmF7e2zC/ACwjXREaiTfJY01vE9d4iJ36lggUgJxEavrZaLUwL5MrCHbH5jJo26IcUP6a+Rzsf+wIhZPnmG6HLOwjlsvjsyvG2gpMofAdKzaJyxuwyYPtNfmd2RbOhGQicGYAN7A/vz/uMqlxGdRWxFi7wp5P3u+TJrNxy/6jo6O4D7RPgrn30tZn5URgf3pq7CK7bZ23ZlfdHkhNJ9AXRIpFufAWGn+nGmbXVrY20jaRvdSuq/BL3ZKJmzDL9ni88x9tItRXNHlZ4C2D5MYQFaJxpLcBw0OgL7sqhhqMny5VUikrM9Bz4DH2CVrx2WIawPDjRqm49OMWnld9wd3WDX/P9fh2J6y48qV8NyW4cD8BovuRiE/Mq2RQaD934Gb/t+JjcBG/8Yhn3Gp4kEDShXmI8yon9x8lUNShK3t9Ahwze2Ghlh57dNY4rYvHdPrxfcnqxMaWCgm9WnP0/tJQ/fHZZskyKjrn9mQUiWLBveXQrZamE+djDDfS+m+Bj8AHpWi7W8b2N9iLFAYHzL3kS6AsGqiDUrfz2MBeAWi8aNmEjBMpB0uHW7VjrFwYyo1AKG5gFQBXAu11Sy5UPGtdUjbbEn8/ZKXF8oe294YwYI3gOGY4sDw3zxIShHouOWzFUfDLTpwzc+NbUAq4HXO6AhtYGUZQTCItjDRWMW/qkyoQsKX2XUe8JSlMmV1THtet8lTF6cPYcKzvTGGFR5XwYnZNVe4b2R/Pmkq03MVgU0eExZ+yLXUUpz0ugffDaC6Br8TJkLT5TnzQxN7FwM+lnWywYnYpR5lXFykB4mhg4KekNBZfH3tY1xNMAlO6bHn0MTcqJMq/betzuCUUKm20QwBPhqYX60k22zWwPpsg/UQ7koCMlAfDCGyswA2EZSccQwMcQN68Z0CAFr7kLWTAQYCAEZ2ptSAB0Dh8INsogrM2Z6D0l/HBnwpHLN3s02MMDC+MMFtojWHDzIH1hdQ9BigHWKHdIVbDYu6CYZ1WVv/zldKk4ed0qPl4GNPtRnSidJ9D36jtXcDfPqhzKBcxZJIjG7bygTxm0NSZbuJyw34ZMtH3oAeAHUvodb7go6JlPU1fRu7a4HtjkUFsiojk0YEESCsFg6EeMP8uCrwFQXuogzSo9CQfyMwvs2+DIxODKPIJNBfF0XOKD4fjClo1ZLQsRGggZK1I2jcd2L7fDDEdUOttn5NtTawyUzOm/gnrFPldxMP0Z7DCS1YqUVrzuH3wmEyHICPQgFZnjnuV+JpI/KGyUHLeUw2NRA8CnSz46bzC/6tvdvOYwBKQ/vfC+n+Bj+0Bo0RNZqoOwBf5m71BKDgV2VEVK7gkhvjG/Q8/bMM+KixDsjsxqqgEQW8DIoxWkxSwLN9FKjAV2hwU+CgJhm01leFAVATzUsCbIN123nOgS98OkfhUIHUxpsDoLGycN2oc3O2N7A/ZyX6HoY5DqgGOCFfRg9E8XxfuFgmXfvA4fwIIuGYAJ5mEsiaMfu2NweFyDiQ/V98MoUYErSvOQCOsvplKYI7WT0FZk+H8WgjkfJwLDLa8Kx+YuG+DBh+A+C7AIjXOr8TJjNyOGtw5ncA+FLVVRMI4KdqFGYPhQ8Lkhk6jivFNwPKChNYn4q84qDLfQTmqHszhqahw3nVXl9EBc6s+wunvsPzIOr6dpsrIa0R+AgtLBe23vvU/WxsIBAaH4gczkUGGUHP/PEistngtvwHEJ2NhzhoW+EU46zMAfDs43rcyGRHIEzUAYXr9Qz4DPQcDAUA6wIJK2bSRGR9xKpw5FYIY7Ogfl2v1QPBN7Fy/Azs1drmUsCwvhjrr3tfDhNLY76dx8NIVO8SSl0vbztxEgtuZH7ssdTMwGHrZVPaLhmrql6zndIiW/BBeOXCaB4VuucGd6JuY34S5HTcRYvMYuz5UROzYjLgq+gZnwHfCgHBdQC+QeSNhKpjWCHF60132Rty7DsAn97TshTHbAe8AyJQHLTjoPafPB44kGIe9hlE7ajrq7u2R0ddAN5B4iEqMHbAF4xmHh6KQuVNwEgc1KkFNA1lo6GMDLhFugPqCdAdqjuMeW522MMA4GiNfQgU7zRdM78TpsD8ujDwpudLTcSlwPjsm5ncUfhRpRq+jUImZX3GAGsAwEKd6Eu6Vy3UkJKIPAqVs5tRv+XAI6AHY3grhdD7PeNzo8rA+iLemVWwE3njM2dpOD66aHTgGsreKfBPlEadYn8STccYgcdZXxN1/VuBzwLCdvrj1Jhf1DFGtjo6iMf/R0bKUWccr2u3+3tYXhtgR3+MgeYBEYEvoYHeAIRdGLW7AHzXYu+pU9xXwnQxwarboqNwB3zTNLCceMwttlFsA0Ss9UX03LZRNOCzb42zJ3t5BNE0NxBsm5rbu2HjCtIYFynrU7CrgfEN+5DMfPsigXosLH0zoBwB8FIPi1j3PLQJ+mMtU/SgocfcLzT3jE8AUHw6HfgW9ujanbcAQ8TaJBPpyJhmYncPeOGdCOI+pQDoovtYBw58W1rsFvPYZ/x5YXwY8IXNjeyd2ooRjlmfJl2D3wkTNbHXaftgiSOgzWqT1ByGxam465AONOiPqy7PlqcBcBZIujIDdk2xqM9b0dd1cgku9iYSPzUwmvg+gl9c6F8G4FsPAN+gb/GBdah3X6XXd8AioE9axt6aeEl+B9ggjdfYeQP/oU1GK7aUa/t8X6KYetBzy+4O4BnwmTRh/UvZmk94G2ZHPePTgLOjGO6uOYzmjgP4qg4v94TxdQ7jw//cgR7a+Ehw9uqgl9r68WPj5dGmSCrvlXR/gx+4GR2CMrcxPmAbEw+wLmWgV2sAvsmeHhEMR9YRDRxu9KiqvwvXx93XzNJrH3G/gYi9YN37Rn0Kh1DsnZ6tE6ObqNsB3sy9BZcA35iM5Wg5OjaD/puHbDmyFOqPX/K4BoQB+My1xr9H/0Ju99nm4J1foTIs1/WZP18OzM+Pc+cqldTQ4UU35Z3qi0cgshEfw+LTIH775kZeboKLv/GcnRrr0EBuBoL+WzPwibQBn7t+pbpRCx2UkB5lsi1e75V0f4NfYHhNh8Hh2Lyyu7h4NaySMH85CwGvgyt11trmvHwwnH3Q+TW9X9T/Sf7O/pIOigIVe1nXZuqeDPFdzHBgAND59GED2p2RItbbHda7Dy53HFaxrfIwQAkm5nXMZdowrZyEbswHxi114E7OE6Y+lhUIIBzcl9yaG/351K0FHfBVpMwOEuHlFK8aAHaiZpwoFPRs7baXexR/WVkpuFuG5+8yYbK3lbRMI/CZF0QKICivdiLAmrXN45zub/ADNgrcjvUBHQByGHUN+CwsPPlgsv17UxnXx3JjcXEvidGE3wGiXOtAFQCVlQWaWGK+utWGfnSlcR8Pfa8o9irwufgXnY9vs9NtmIVVoaFR/B5YBsEGbrgG82unDx6i3PhlBmzBVWh8z80YddALbW461LBszYDP9XxLb+Ag+wzsiLkFC23vpQDo76wBcA3wBqADELa5RCf+en2FejCGt2HStwuEobIa09uyvjudIzePPYGrCxE9B8B3AvhYyJB5iJm/SbfG+HsAbkK2wfhrzHx0D6CZL/9JEhHdJKK3ENE/J6K3EdHX6PGnE9GPEdG/1O+PCve8gojeTkT/gog+80oPCorcjTgWUnPD06VhCnyVCTWskojiY7cqIrIqn6W5AZ9uWLSpB7vGjB8RoMwlZtgvOA2Gi81AHxjeaIzpGNQh4LuNnj0C4kbcCmzKle3BITyyvkvZy/AevRg/B74N+4vvaGUIm4+39br6/8I98+v8RM3IUYOrVA8UCMY0kz66vY+jn6haXcdYkdsACkN0maG+OvYc33tkwQx4jEE/Nm/8zjAIdQ87peg79s/Z5/K0Qvb6/gMA/jCALyKiTwTw/wbwNcz8QgB/U/8/mu4a+AE4B/BpzPwHIXv0fhYR/WEALwfwE8z8fAA/of9DX+BzALwAwGcB+LtElI8+IXQ263jeCSc6i0bISLFqcBbu1sdiKkqabm+zTlGUdoeLOjRyFIE79xcT52pbkzuWpduwXBnRoc6zYVpR7Dwmgh5LA6ObAeAmEkpwEj/kr7ax2kZ9awzHNQG8UeT158Ry5Mln4Q4Y4X6iTSeWNuxoAD6vC5NAtF7Uotq8ELQeMtAbIgaA83rtN2fa1BUPv6cTB20AJo6BPkIQddLRqVOcpA99LkvM/G5m/mn9/UEAvwjgWfp2H6GXfSSusO3t3dzAiAF8SP/d6YcBvATAp+rxVwN4PYCv0OPfw8znAH6ZiN4O4MUA3njsORQ6ne8vMaHsvcgLcJXR2FZLwMVd2rAM3nawQxT+UEBT2L2sg5oGVhN+d5IjtYHV5UVbpuNlgA8o8nzabYcU55tsGM1/d3wVy/vQM4byGBBt2ItJcPqsaOSw44fq6tiSvS7/EQSjsaMLhsE9YHnf6t2l4oPaRGtAxf5MJEKn52st2h9XfQcx2taTWq+d8nOWIgiaLjnu51tZozxrEUL+IHM+jwDIqDUhpYrCdDrRl4ErRl54BhH9VPj/IWZ+aHYhET0XwCcDeDOALwXwI0T09ZC3/z9e9qC7vWl5BvBWAM8D8C3M/GYi+hjbt5eZ301EH62XPwvAm8Lt79RjY54vA/AyAMi/62kb4HN1iHfa/n6Pfwf4/hCuRB+ZFaMfbMZGYrLZm9vDY0RnU7K3Amj5rNOaLlHBlgP4uc6H0HRpCPdvKgcN2LQIBkpTsDsAel7OA8A3PtOfge09LnGNvmeTZ4GH2HZm3LF6H8T8GdsdXWziemzb72UT6TtsdNWzPnP/mM9lAHTlBqSvuVFDa6NIvraZ0GaasHez+5Q1UgSoyQQkBjXy6zjcw6oHthtId3STSarVrnyLhrkiIaH6BvEC8PJ9Jbi6Yrqizu+9zPyiS/MiegqA7wPwpcz8ASL6LwB8GTN/HxH9Ocie4Z9xLI+7Cn667+4LiehpAH6AiD7pyOWz7rWpe50FHgKAG7/v2UwhJL3PwpeWSxCiGToaAM70StMekATkSExkYNRmsSDyjY5m4kvTGbZPF68u6RDRzjILAeXv4u8dwMcGPULeVpsRWANQzisqnLtkFBgD3FwanxGfOXmWDNN28iDwjRPU2EbxWVGfFoCvaig0AZywEiJIDwA6ptcsoPEtSTUejGoGjeawFOomAKAxvgyA1Shi7M+Ab1rJYZIJEkgqmG5o1L50UgYfBUAB1YqU4KufTpVCld55XkQ7CPB9FzN/vx7+fAD/D/39PwD47y7L5zGx9jLz+4no9RBd3nuI6EFlfQ8C+A297J0AnhNuezYuldt7tmfAF1nfGAdvvN2zsZmUt58IgC3AAekuYwRkAiE5w+OsgJiTA6D76yEM/tgLI/DVHsDCqzYmFcGmczeB1ElFE3MOsMQRlAyYR2V6d9shxjmCj90+DmTaXjO+oz83DPBTAp9vZ9DpJNnvdeCYlJG6cwEME3wnPSild6tsbRNTJ8uyujKZq5OxN2sLe1Z8R6uj2vb05QB6lo3NhkkhtxoA2mSbtQgMCaCR1fXlQBi1O07RD/YOEkmhvh3ALzLzq8KpdwH4oxA12qcB+JeX5XXXwI+Inglgr8D3AISCfh2A10JQ+mv1+wf1ltcC+G4iehWAjwPwfABHTdXyHBwFPdefMdxFoUuH2mM8Hjqlx/RLBORAsUKh6pJkze6S2vUmBg/PkWK3UFiz842lBTDXnk8cBpeCg20A1AHT5J1uR6kzGiRiGY/fiFZ/Y4psxn4YmVfgvjLwWZYR1CPwTYwwHpl5AL5YL6Our/2WAifFC5m7qkTgjmXwPUWU/dnLGh4k3ddX62gKuvFdvZzSZxrgSUpo5TGmKcdsbxhjgq3CmFnrQj5XlaJuJ52I+X0KgL8E4OeI6Gf02FcC+CsAvomIFgC3oKqxY+luMr8HAbxa9X4JwGuY+XVE9EYAryGilwL4VQCfDQDM/DYieg2AX4CYs79IxeajqXPGDMA3sj6JwNEU1L4W81DrhkHQu3eQdJKlbTokm/6aiEMOkFDRty5tRQGrWDzTfxmIWUd3rAnXRUutHLCbFDiVhTiLGJjrTB/X5YcDnXS8nw8cf5Spe3fL9xjw1e0zpwaO8NsW98fwZJ2RIhYm5ttZYsJl2s/MUECkWw/IIm3ZT7nIb4sUDpDsJ8OQNcEk7M+s4V13nEghdtz8qm2jqxSMFpUBZCGzvooJrGK2FtjAN7e8/bfvA3Ja0fcUCkRmfgMOj9o/dDt53U1r789CLDHj8d8C8OkH7nklgFde9Rm9n1UPejPmBwApocXvI3gjewTeCEo2QEZ3iUSyP/BCSAp8XLgVyv20FPhy0wFuYgNi+939jueCJdHOO+tj2+0tgKB16gCATh6u0hEjmR0GX5fHobxoOB+7bGR8Y5lCeTfAV/vzm2fF/+3TuZ6gX+vq9Xt5hRxyALZ7ZXIlpCR+pKgERkIXKRwN5NlW8NTA+I6wvuk7OzA2QHPQy8r+GEggVFYADGI3FCjZCpYU+DSI7im53ylx9BTpvl/hQcG3bxa5pYvYrDMlkMBcQZxl9tU1tJzaOsxeZyTnaiaJwrKzDpZQkwQuQBZA9fsj2GUDQQQAbM8B0InDUz0ZcWMwkQn6TfpJ7JvbuEVwYA+u1xs640yU3eibZr/jPbT9dwPU4/vNynEA+KbgR32+cTLh0JY+6dg1E/bdMmiFGkVegvn+teOsJtlK4jhvAFgroxbpSNUMIWx9S/oVKdBsXRNCXdS+XqLe2HS9zNpH2QWRTiROxgCt0rSFbKWKR5ohvemkOj8AZewwj2+678GvBSfdLs3pvNOJZTWHWgGYZcrllBRQFAAzD3t3ADSuCKgAdgJeFojUZ14EYPPdwHSfkByYXyeKYTsQJwx09ttTBAPTWde2qB8TZ9cRzNylYjg+MrHN82IKY6sNr3CpiWt6cDO8DHSDeLsJxlqtSdknjRHLeVOPvKnjg5bnA6lb+ZDq4NMuT5dgP+JJUBKh1oRC0IGfZauEKvpg2J4w2Jajm2Si6M/y3jKpUD+RMZBY3HcSKLC+AIB0AABDG4/+mKdK18zvhImIkXP130SMTA0MxwUX1jFXSuEYgbRDcmXV06EHwQzUhRA3H4IqXiweX9d7SYGSGvOLgNeB4ABqPPzujrlrhh0b0EPfz1mRzeQV8I1/jProgHIXG72/SxPWsQFAe2Xuq8DewazWPsx4YISHWN/o+O0f7sugdSCx8HC5QdFBkbfHjt3mfUqALyfu2F/SyRWQNig1gSqF7V10ss3azwpjE7EnJgO5rj7CuzPcx9D02WbkEObHDQApkDhSYPS+o42gCMzc1isDuJI64MrpBNbeU6b7GvwAIGuUDeuIObVOat+WqnXKoLMxn7+6sILfsHObrgSgKoyPvDfJoE4FKvK2MvVWYQSxqzd8dOs4Z2JYYC1xzWzbPjNcB8B1PzaDq9hrjE6cYbmJNMGa6sAU2F8HfIPYNer/vBhD2c0B1wGQB0YY54wDwNeYH/dlQHtGZKKjfrObWMYyztIgMbjPOuDAl1NF0skWaFKGSRc5VZSasA+ini2pJDe6BHeSkcXbdwD92WoWX7mhvoNu61K2SUQaLUjrQQ1hnSN1lRaRNiEVYli7+eko4DXzO2EiYiy5yAQ6AF4OwBddXFZilEQdK/SlPbocqOqa2rqwr7wQ4FOlseQq+SbqwUBODWxtvhNcp/sbrJOdyJtanp0u0tifPtPK4OwvLHOygdQ6vHbuKmLXxi1Gs47A1621DVFtrAgRXNyKacBnPohDcb0Z7PgB4EuFu+e3PhAyoFaXXbrKoHMmPVQABtZHMuFmYiyphuM9+DHTBs9SImGI/qwJsHArxriapYHfUAFdJQbXFwoTj/Zx9yU1q69PdBr70PTFrkg9URrHyD2Q7nPwA3a5dUACvEMuqSKhdV5AGjynhFJTEFmWtuJjUdeEhVGLAoOB4aK6Pu00TNKZozjSCgaMcewi6EVQ3Oj8BjeMDSM01hd+byZniyZtwOeyJg0ACAB9WKZudg4sy8VPbiA0zuTO5Iytgpq13XSPoaoi1ka2GRlPB3xR74eWbytAm4y6VS3YPvtYct0xmk65ibsN+HKqWKh2EoYBXw2NUlkMH0WRzJfExTSy6TDhGOszZt8tE+vEV60DDnXAAfjs/cMxjNegXXtKvR8BYhi8h9L9DX5gnOXiQGaA57+JHQAB8XLfMWGtCbnmfmbWGbsUA0Cg6jcVBb4FPkvaDN5FzmgF63R3G/YW2N6G9UWRzECE/LGB+bEC5QCAYQBR8NdyuZYgC91TGwgeT+6Q6OVg1EJybfR/Wj531yC5OTLjqF+MFuD4nG3QAhX3wiZMDpTU7gMCwwzX+TuOyes51F+ocxgAkgBgDuLukip2uSBTddYX1StrTSAWPz8TgVdKLj5v0sDuOmbdTUC8XXIZ6tomt7glpq0Dnr77oRTr9pQAeK3zO11KxHhg2QOA62ASMRYqXac0AKyQzrjWhIu6+HkJTAmJ8bcL0Z0Lg4voA6moq0yFKpZJ8jUwGtt1An4xnPmMDXb7dcT7PT/uvsUfKz4zIIudN+yrQFvXyQ6AriecjEoXv3QQEaMPmhqBi9qjOaETR9lPordOxhSf5ayHu71Ipizb9Bc6ERnDdKboID95qNcxt//jnhbJrLrswLcLwLco80vjiv0kALhQRVHjWtwn2ussTAYIhqpuZUv8RJBka0maTnyb9xx+jxLF1Op9QuC7FntPnBIxbiwrEhhJZ2Fjf0sqovcbanzPCSvliZ4GqLsker+FwCtJrLcCPabAtzQGU0lD0c8atmN/PehNGeEg7sZ8Gks59JuHDm4AoSAIC7kUFrYrm4vRhTsleEyjj10UPwfGuymfgavS1479oV3jedX23Ym6cSsBY3rGePQ9HRzrROwbXmmrUuAumkvc5zmriLtzsbcB35LKxqiWWO4vQKdy8Wfrw+N68ulSvjDpRP2cs8LInuM7xgk2BNXojW8hsEM3oR7oA3ecrEHunXTfg9+TlgsFPwHAnYsi9rvNypUTboCwrxlLzUilgV9h9claEkpJCnwcfPQYlGWj8apLhySR+tCEgvnsehj0Dh2zNAXB8XfHBhGomGWg5UrqwoD27YAU7rXydv0+DrYNE5m8t4X1crMjOpazSeFYxzQD0I3WTn9Nkt9V6a2d5wiiwzOsajZ1p6oEUsBLuSJnBTn9zg54FWdp9cl2TMmtOyPoQXz8LIpQhQD1uPfKMNl0rHisc6AD8k5fbCtZzMNA3zMej0v+pixwVOzeQTphVidJ9zf4gfHkfCG/ibFLpbE+NDE4o6KofCjiSEai9uprFSPImjOydnqJyAL397PYbNF3zVmFBq2MTreXAdyGAd7mrNv5X8WB3K6Q/6uhjzK/ZIvrtwC4kXy5ddhOJB3cTjblCGKuMxLu85y5uIzX9Dqw4bn+vBBQVBmfXENe5lhE9pdBcxvykPWsu7RJH8gKfLtcsMsFZ8m+14M65ZhMzVK0f9Wa4GHUmNruewPgbaJ8H5s8vM+0lUNxZVKNblUxjqGzQN56D5yc9Wm6Zn6nS4kYD+S9i7wGeDvV+WUSEThTRdGZeE8ZC9cg8iacpYJ9LlhyQU6y5A3+oY4VzRhaE8N64LPfG7Y3Ob7pcBMg6lmSxYebpHiv7RULeZe46U7HeoLI64aL8dmBfXSGj/jcBMCU7AfEzql+afKsJo5zB3xO5nXwUg1Ls8Yyjc+lKPZF4JMd2lKqWJaCJVecLQVny4pdqg58xvxEt9yzPnNyrpzc5UXArxnTaqHNXjHdFqRRvRBF3ph8AqWO5UWAi/uU1KWdq+bH6hM72gqn25h8bzsxrq29p0wJjAfyRQd8maqAn4LeKPbuuOCcF/2fsKaMi5xxUbPrdJKKP+ZMPDI2w47L0sFBfggQr5IGhgQTYQOo+jcPxzg8Myn90uqZhlMayhXZ2MzdxRneBPAO5Tm9RAd7FIMj+7OMaXheJ1bNGOkgCiKh7cm7MFIuWBZhfRH4buS1Y3xLUK1Yqjq5VnVzMba3KuMrRZmfgh9sy9Fxk/nCB0VeIEywCIDnDC840QcnfU4KgFk3akoIe5dwC/xwN1nf2Cb3QLq/wY8qnpQvGtAFwMtqBMmhxgsR9hrqolDCSglL6mfzuE54lmJHnFnfAMDWjxDmwObANx6bXd+xoCCqDkrvzpq5ATATb4NI6h8WkYmg/mdX9+nvWMmRRfCHdJoG0JE993mHV+rYpxx0R9xgUBnz6XRghBazzramXBi0VAe+3VIc+M4moq4ZOQAgE6M42zNPgoy1JqycsK8J+yIAuK4N+AT0wvaos42yJj6NPqF1bK+tH7fN1msm+daP7UxXR+DzoK4T0LsLCrprV5cTJgLwpGQ6vx7w5DuKt4RkfiEJ2HHGPmXsatPd2HrgaQyzOBjjzDxR/Ktw6fcdAsGOAU6eZ+zJixPKIH5cGoSyklpgIt2bVJazB7P+ovPLc0Y7MN0OS6P/WVdengPgZWL+kSLHujiUjjHMzk0oKP5hu7QtFbRU5KVgt+tF3bNccCMbAK6dW4t5EVQQMjH2NTc3Kk66rC0L8ysZ65pRa0JdE7AmYKW2NekaWV8EPe7wxyaI0V9URFkK+w+rZ4J9cgDB3EAQvlkTmsgb2+tupGvwO11q+r3G8CLo5RDfNhOavx8nZL3O7onJFom7d3Hc1yOIfdH9AkBT5Cuj8tUGQM/SbiN19wfwbdEsMQDIFdBkwsBGEOzE++76KwS4DEA3NeyYZfFgma6CiKH8x57b6cLYP1gYSUVcA74bOwW9ZRXgSwVneXW2tyPrLzaLiOsUgI71XdSMwoSLVQCwrA34qBCSibpxb+aiK1kGvV8EjL4OA/AtBnztWxgg+v2II/BlBTxX7UykhlNilfXbeyjd1+BHaOAHoAO9PACaGTy6Y7COm1Chy5IYnbJ+s61k2FzcnHA78axjNrrES5XxFvVqtHYeTJ2oK/+TLVerBk5q+FCLc2OAk7z6ynORFxSWos0+en3HAkd86tx6yIEHOlCbrrHlZXkYM47BPEdW2pV78y40AB91Lh1u3VQAaIyvYtmJmLtbCm4uawC+xvqi+5SlGvqTi7ycUFhY38W6uLhbawKvSXR9KzVd3xr60jFH7lD/zaqr0YGCgaMuKu7uFAh3QN2puGvAtxjwsbv3eL2O6epz0KWJwNdi76lT6hhcz/ZmqXJCAaGYRQ7UWeZqTbLCoyT1wTJLXLPIyczN8m3AqMkHtY0Nmlg+DwFfON5dZ7osA9CqrspkArYCIKGt5b0KuA4g52xvYH/sDtLtOgcsfVCnV+uAr//EDYW8DLU9M7LQWE8OhmzMcJu/Pzc+M7p3LE3PZ8C3W4ozvgeWPW4sCnxJxN6FKnaq4+tXCkkRKot+z3V9NaEwYS1q5CgJdVXQWwPrGzekt2g7E5VCq6st8DW9HnWirjBAVqvvyPikHjq975i41fHJUr23qN+WDp0oEdFNInoLEf1zInobEX2NHv9qIvo1IvoZ/fyH4Z5XENHbiehfENFnXvoMcAd8h1LhhIKEPWcUECrL7z1n7GvWzisK6rVKh+Ua3BFMLFkH4HMgjKILNkrrmSI+pmMsMBpVOounO8eSslNqa4yZMA9eun1QB0I6ELbH2mcDbgGAjOX1gEc9CHWrClpeMY9DYvI2bwRAaD5tnZuHsj5jP1gYaSd6vmUAvgeWPW7mPZ60XOCBvMdZWnEjrVio6JLJwbUlTJyR9a1FgXDN4jCvE2nn4qL9ZNNfjIhZf/F4Wtt3b+9JnahbFfi8zmfAF5y6N2uOu1nncN+8rWRi72WfxzDdTeZ3DuDTmPlDus/mG4joh/XcNzLz18eLiegTAXwOgBdAdm/7cSL6+GObGDEIBQkZ1X2soGspo5hbmRz8zusOt+oOe85Ya8Z5zbgoC87LgrVkdUmQmRoFOltTA7o1giA6nZ8P2th5CM2KaszNPqNIt31B/SZYNA+2V6xwzgfTUSYLR27M73DmMZZcZFy+LnfGCpMMUDZ9o/nWYc706uBu0YFWSBR9DFNfBk6Q5YVJRUGPmg2AeideEXHJxUAOOi8o60u74gaOGyrqOvAtBngFN9La+Y8CpiapqJy1j4nkYKxvb6AXWB9HUTf69Vmghui/GCa41k5939qKu6GOp5+2I5tLJDrJRcDjrq+FvnfC9IQRe1k2y/iQ/rvTz7G3fwmA72HmcwC/TERvB/BiAG889hzR1yWflUc/ShNzje3dqjuc1wWPlB0eKTvcKjvcKovqabK4JKh+hlaSmXmNwMft2zqzAYACSA3WDVtoz3fQoUhFXYvGPAdAKQArAE4jvhzKfKL368TeCEY2gECwmIDyngPDm33sfCxTrBM2cLX8uImwI4O1sgVrZ2fVzI391J2wHtqJI/OiFt0b+jHgM7a3cwNHP+9a3yquJknu1lL0t/n22STaWN8AetxLBs2lJ1SIivk9wx4Y4Kx+TceZuFm3Tb8XVCIdFkXQu1sA+EQBPwDQbSvfCuB5AL6Fmd9MRH8CwBcT0ecB+CkAX87M7wPwLABvCre/U4+Neb4Muifn0x68KX57gLT6JI1i7q26U+A7c/A7LwvO1wVrSaglg00/sw/+WAp4ae3F3w78VKST8OECKsZqKAWGpumQpbJLUR89AqANDibp2A6AEKbomyQNOkCGAOSmclu5OuZFcL0fFMzdkBHK3otjPQOz700Q1mhFYQXaMLj9N+vL1lB3NFo7qbN68k4BcGFgV5F2zcBxthQHvrO8dmKurRAaw1QBtnojSYCMoOsbWZ8ExiUX50xnPO47PKo1xM+SO7fMKfB1e8OgB8Sg9xQ1RuxE1qeo/Y6AV4fjJ0vGAO6ddFfBT0XWFxLR0wD8ABF9EoBvBfC3IVX7twF8A4AvxJyfbGqLmR8C8BAAPOsFT+NbdTd9djWjhn6LmJtw7uC3w8PrGR5Zd7i1LrhYM/b7jLpXl4Q9dXq+tA+/V27ML/r4KePqRMiqnXB8kwA0m3SgJgwA3cUlhWsV9JwF2hpeAmxLzoP7MQTQZRXTm+NzO9+9D8I7RXE3ONy64r1jJ9xrmsOWi2AFSgZStcHOQDbvnhCGy0XAOfDVnSr7dyLyCusr2O3WZtlVq+7NvBfDhq0OGoCvLVujDvQuiqwM2isAthUd1K3mMF3sTH97qE/YvNZNiuF3Z1jrAJLDBDP0JQe9oc9WOPiZDtkA8WS+zown5u5tzPx+Ino9gM+Kuj4i+jYAr9N/3wngOeG2ZwN419F8AV+qNur4ANHRSIcV/d6eJY7fI2WHW+sOH17P8PB+h/P9gov9grJm8F6AL5llzkXeJu4a8KWVxUpn5dHB7wO1iq+fW3kPtT1tOyuPnVb7LNlJglh6ac4C5RYNZ6TuDPL/kQ44guB4nAKQx4FF4d2HvU+ii0lbR9qykEEY1ASsbC/D97zw4AEU6jA+cyFx71AXD3fv2LGA364iLeLPt8vFV27YZ6cWXVu5EZP1n8qE87pIP6oCehd1wb7mjvXVSmowgxqdqHeKN4YVm9cmmVDX3uQOdjToVdv+MKN+VvKZtHMEPCsL4Mw0+rOeHPyAJ47Oj4ieCWCvwPcAgM8A8HVE9CAzv1sv+zMAfl5/vxbAdxPRqyAGj+cDeMuxZ1ROeLicNWMHGuABwMribGrAt2qnNVH3kf0Ot/YLzvcLVmN9tvRoBdI+AGCw7CYDvsIN2EhZF4Sd2E5vNABftyIhAtwh4LPE4duAaWSBCR60VKBRgcRXgNi548mBbyyTOTgbC+QwOEfluwJf9DFzcSxIVm7s0BeLIcM8biLkuRsVgxo4om9bPdPPjlHPKrCouBuCFdxYdK1uXmV5o0YDismDEiB1E+d5WRrwFQG+uIytFOl/3Fnbadqo0bWon1ia6CvXGXMPYDfkM81cGZ7s0yJs3imlqT4c+NB0koyNLvIk6YkCfgAeBP7/7V1drCxZVf6+Xd19zww/kXEUJ0BkSOABDUEihIgagiYiEvVFwwPIAxEfQDHG4CCJiYkTUROjr0RJ8AdxIkQn+EBAnBhjIv8gMCKgREcnIKhxRnJOd9VePqy19l67uvqcc+/te8/te2olfbq6TnV17apdX33rH+80u18C8ICIvI/kH5J8PvS0fhnAzwCAiHyW5AMAPgegB/CG0zy9gNrz/m+4ostFNanxew58HsriE3adOxz3CnrrfoHNpsPQJ8BZ34Yjz66xvkEsCl9CAnpgJYZKzGpvG+en7ppI43i2KWk8gPaHBkiFBToI+gFFAJz67XPOxYaZ+K6BwEoCsJX8UmznlHo+qe/L6U2mzcRacl1ZX/WtYGDrJXb7YsxqcNBbCGRV1V337q4WA650A5bJUtZM1XXJQoAJg6A4NPyh6eCnqq7Oo5NhgRMLbdmESIFq79tmeeX8RTbt9lSG6xcZYsTOifkxlQXk9kV9JLftKN22VzOVbLlvP+8V/ARAviTgJyKfBvBdE+tfc8p37gdw/3l/IwvxeL/aqqjRJJm7F85j+Qa11Wz6rtj5hk0HWScNSdh49D2bnMsm/9K9vH0O8Vj1vQ1qNjZ4iso7tue0J2X7sxvC4YzTaZr9SD0kA0Da0/+0aP7zSDjGch8FNayx7bnaG0MupsDPmkSlnmrbE0ESNo3nqJjUnoto7+sAWQLDSpDtJcsMrjK6ZdaCBZ1VaTG11wtZuGSjXlnQPDijxuDg56ruJieNEDDgK5VbAvB5XcF6DvWp1XjS3bnjYBPUXj/HW/Y7O38xkoDG8jBIeWA0NRv9EDxONGSXpFBhZhyjuh85ze5zMXLQGR6DJDy+uVLBD7VHh0hNMq9R9x2GTLXR9Mr2hj4hrztzcqTC+sbBy6rGjl7hyeiam7cHHMdtCRQIZQfyTALf1kbh3dXfeJe44c+B0e022QpWuioU99eEONTjbX5vLE5Kyg3MhvlFlpejvW8R4s18/8ZQMtRLrtexNtx24GtuxAK4qM4NT+VaGetbVe9uqdRirK/2eFHWlw30dFkfoJuRxrAeuhLEvHHvblB1+76rcX0x2LwBPrSqbuOdtSIVPp/yjtMfH6oWPVCAyvvz+gTooCAX55aDZJnbrMUVwvvONLvrkRn89idZiMc2R2XZ2wb2WS04XlNtyLW+mj6hiTx0yD0h7t1dE2ljE6JJOB9V1S15vhKYHRTUTN3diukbqb47K5EALQAGhtesm1iuXmWW+LsmJa78brtDjo5N39kA4CRrjUw1hLjkyPRGSfWy8Duv7lssXhFwlqc9lTMFyRpuyziFsFGxq3NDGZ8CX7fKWCx7LJd9yeJYpWGrAnOptmxqrmdrRDNJH0DPH6A1lk/nk+QJlbc5XyGOciIkpZgmRteX4+vt19VzuwdbZ++pN7InyqjdSVK+b+wwZidtmXiidrNPtXe4tdLbDh78vrFZ1l6p9j5ke5eaq+uxVzknZG9P6cHMxc7HCcYXXue0gbCEb0gbT3e1E2kMfGX/7WcJC03MWFm2htQFTO0OGhkbiYlxjljhWEqDHLuRi/Nj7Om1DAuMmZ+I2vMCO8nulqGyF5qK5udEiKBGWxDzUlXdCHyrVY+jZV8KFmgl5tpwyBmfsz3Nza3hKw54gzs0hs4qMns4i3X5E3XUyJBqbF84X2OnhliFcEkAilcbSMLKjiqBb1PePMYTKGW1gOog8lAhDoGRo/2+z+cCehuLYvBwrlJdRvYHfhB9it1CcuDgl/CNE/X26n1kIJg5PTk94t6r6WZjeeVln0eAdy0ToAAgWjA5lfWZODadKYEJ+E1SS3HJJGAVzJPqkGkAL7zGwbjjY0RzQ1ewq/m00gLfImuwN6Epf0IL3AbQt/tEArI1hY/qYPzNsn+P5VsNSEsNaVkZ4xtXanHm19qGE3qLA3XQc0dGVG0n55SPQVDVXQ9xGZ8v95AHOx/i586vi1TP/YiVE6gAKAqAlenZd4Z6PcaaRFOco6+g122kjWHtZ7X3lpacieP1sl58QZ2YgqJ+lKfxYJOyuPVZ1FxlfF5dF7u9XW6zsRzTYsMjUHrIRgnfn/LK7ZpcZwJgvClGak38f7RFNutsPT3A1VOu/AYZA59UFb8eI1sAtPdi+/MSUgZ8aZHVflccM4rEeaDuayAESbNY7FqUGLlyEh1kA6guPWc3Y7nq1avrwDcqQd+2K+UW03PQU0dGqMwiqi1486EmlMUvWGDazTXyFMIk+r3g2EFga0lQmCCkssX44GweeG7iTSgqrgxsgW80N5KbcjyA34FvIwaEDn5u194TYAkuj7f3ZogIsD5Z2M1tEykEliJ63XJdF138yNWT6yWqGqNy/D278ZryS+F/ZRsGxgffVwCgACgS/+84GpavRRomCP9NU32bAdXfLipvyEZog3PDcU+ck+18U6lqrgFfss5o3ibAH1opEbkzc0RCCWuRzpoTbbE+369lbiw0jk/DWYKqawVJF27rG9n5HPg8vXE9dBrzOcSKLKottI4MnH6B4jkKDwcI9JyABnJ6oWUEepLcZNJeU5+Tdc5x4sEjdXmK9QXml3qpwLcRpHU2AMzgkC2OdY+ANTO/PUom8nG3/dSdALvIcgqriXFOzjJ8koTr1N7goVMY0Ez0KY/tluNgBIRjp4R+6dpU7TNFKpFgPF8h5zQCXqsWT9h/OHo5Iwu2Piy8HaSgWwzoOmsQZWqvmymGQUAm5CSVYQ2sYUPl9yQ0HaqVmFcLVXOvdINWZzG2t0qtZ7fP3RbwHfcLHDcxnxoJUAqQuonEz9VYWN8nM2gIy7oJjFDULqpqa2WA7rjyAg8tSx9HGEjZQAGQRTOJsYFx3nGQEqRfmJ4D3yaDmwHcZDBnYNhzAdI97IvkMwD8AYBvg16Nt4vI79r/fhbAG6Fxwn8pIm8+bV+HDX5C8KRrmUlhemiBDyjgFrcdG4GjkdpLN7XAp17MQvtyczj1SR+igmmOD7EJ7MdVQDSCYvnSGWMP6tDkvyPri2oSWuBjBPzRg2Kr6siIgdXzZK/g6XXWR2N9pU9Gykgpq0fXmJ8I0XfVMZVDloQErykpYIJ+v6vxe+NipEfdRtPWrDmVSy8J2U6Mh7OcGPCdbJbYbLqa6ePZPiU/t50ffh7clicGzE0xCaAEL9N01Opt13FFlTcDqolknWvM0gBgzLzQaxipoe1/dG2ah6+VX0uDAh4d/NYDUu/AN4B9BvoByLlJ37wuEQGGU3MWzis9tBjKx0k+CcDHSH4AwFOhlaGeJyInJL/1rB0dNvhlIB2z3sgOZIGx1ImynXHRAGFc5xLCEXJnk1FoMWmijCAiloFek4kg7ctZX7mZgl3HO1CeNd3KpB4DJNttmmOgr6zo3gIfmodA+Sw7WJ//XgjXEFe5HACN9Tk7W3aaYtYlbQLkhzbkhIWHJHXe4LuyQgBFVU5J2d4i5RK7F0HvqNvgSqrxfC69JV1v0BUnhwe8O+Pr153GfG6SOr42HvIU5k4YuwDmuBA7D95SIFzXcj0MBAMzM+4GiBQbonQojoviJQbgBTQ8/q7O7YBuzRxh/R6gmo0xP2YDvSGDm8j4BmDTgwZ8GPJe2Fo9qOvfl6XGPmrLj5F8GFr96acBvM1K4kFEvnrWvg4a/ChAd1xViajWRruVb3sWqrDiQom7yh2RLEwkG3h4lRRymzE2amDYb3lq2+SN5a6KE8VVnrMPtZFG3Q7qTruR73f0ECjHtIP1jdVg1P0XxhcBMDA/dgZUCysoYGDVJQUvlxKUbsHoHpuZc3yKGPAlwbIbSj/dCHp3dBtcsbJUsQjpRjokadtMDla5e9NrwHu/CcC3tnJmm9hMvD2pkiTE5jmrQ2FgReLXjBHS2g8AykATqgNEMhpNIz51GNibOy3G3fTq5n6x7MFVGr6bM2PImqG0yeBgwNcP4KZX0MsZGAa1S+5Lzrevu0l+NHx+u1Vy2hKSz4Rmkf09gN8C8H0k7wdwDOAXReQjp/3QQYMfsoFfZG8+CTJ2A164ecfsyZ+21UBtIBi+7CWfNNneGZSEfUSVF0XdZWIJ2HWwK8BnvxkB8DSZCtXbJeV2NABsHB3u5JhifQ7Ugi3mt2Vsb9LYFPyS9crw+nnuhHDwiw6IPicLNk4NGLokSvne0qqxeAHSO7oN7kjrAnzew9n3fZyXQAbWXCAxl74bXpggqrpcE8nBr4dW8bbzUQ/GGH4Hc2AAsZteOVWmCkd2DwnzzuehUOdZri8MdW7Uk47C3lIv5qwTBbY4BwsQ1omr19E8uIOoStvnouIq6A0KfH2voLcfNbUe/PlU6K+JyHeftRHJJwJ4D4CfF5H/JbkA8BQALwbwQgAPkHyWnILeBw1+FKA7wRb4TamzQDvpijcshQkW/4e6TquKADCbSrKJTmr+L8x4H0MR4qRtVF3C0pCkBPc2XcuCsfssBNw67tO2H2lHTUhLZHfjz6OHSmNQT/VdS6WjFCytjcBrZ7SjxQZHXY9V1zdZFkALgA5OgLI0r6+3oNoMV2koLO+OTkHvzrTGMvU4Yt+krW2k02BgYQOIJfNnUC8zelV1HfjSuha32Bq/2j3qCkoN9A6rLRCvVNQJg22ulwc4SxaLRa3sb2ouuXMubbLZ8Nw5AWN4UnZcbIMGhMyqyrLPheEVNbcflO0NjsLnBqyzRQDZU5CztcV4D4A/FpH32upHALzXwO7DJDOAuwH85679HDb4ZaA7xgTgxcevb9yyFXh0PaDqS5xk4eZ2wGICQL0RM2HMDwqCA2tA6Ph3LZxBWV5VmdNgKVwDtxyIW6rvOdhdHGe7M1s9YhCV3YWudAOa2LpmTDLadwTAmK7lwcydoOukFBK40vU46vrC1rwBuMu4Gk9sDemNwr2TmldcvjOtcZQ2uDOtcSVtsGKPBN1W+2skeLHbE1k09fqyWK63hbJ4UYvUa5qjZzvEsvM+bq/ZWD7L6Hr5+XHgsyZBTtUlCby0vTPyJiXNS4ONApVjuEupMBRVWGd3Q9mh2W3jZ3/XcBZ1RAQ190YAn8se0tuojWd+H8DDIvLb4V9/DuBlAB4i+RwAKwBfO21fBw1+qva27v+x+gugPj2NneSOZb1P3tpqEtUvIDr5km1XQK9kKdhnAzENFB6FI0TwcOYH1LxMWBL/WP22bAA3rO+UyPoKeE98Q0bLghAOhC2v97ipTvlq+L0mBGik8rqtz+vnXVn0pSvaHd1mskeG18/zZRe33y1SxlHa4Ap7fU8bHHGDI3tfsi/9mjeywNpaG2xkgSWHYgOsubxWcbn00p2u4+jnJ57jhHK5lamPzo+fPLrX29Vf6HmVTIg12yr25AwwlP1Prl5T4I3ci+NjcBVWHRbozXZnKm0BOanqLqZA0IAQIvu1743Ff+f65SUAXgPgH0h+0tb9MoB3AHgHyc8AWAN47WkqL3Dg4EcBFsf2wRmfLdeNKuPLppKVnMgJFdeftF4nFPaQZnL2B+t6JuUG8J4dCQooJaLL2ZJIMZj7b6cBXm4tACEaduV2R1+1pdaelxGW4wjsIZ/+atihTP+2n7to90PyuL7csD4HvicsToptbslhZ3N5dQTohXQbnqu1CnrrAnpH3GBlrSU75NKp71iU9a1lgWVeKYjag8HzwEsRgiF4dmPQuzM/oNqC/Vp0aEOKynWThvmlTsxbrRuJWO8N2qwxddULwrITeGe65FpK1Ez8gTSY+uphKoM5LZzJjUEOqADkBTmCiqwrwmA71u/tQ/bj7f1b7J75r76afR08+HUnUpYBbJ9gstijINZUhyhBtj55I4txO6DYft1J4ba5UkQ0AqCwTKSG0DjghNiU6jABMBgghlSxKMXxwt0g1Gw/Vm99XWRxzvocBAPr03UxlKKe6+rlrYUIir3PgA8e0NxlLAwAj6w5kLM+V1eXBbDaazbYoHx9ovbWWFEBU9neugDfkgNWaFPXjrFAB8HAVL4X/w+gVGGpPZDRvvvLzSjUE+7OJrfllutQGB/sIZCtL27NbAFQQRcJWTJEks6dDA3wjs2exiYZ/wl3XmQxxhfCVMbe2iwV1KaEaXrZx7UnkVusaflBgx8y0K0d/KafKk3YgLEupjB5gaAuojbZiSqxaLFH0JleUFVEakUNBxbP3UIFDwEsYNRCG6K6S1VjStMgoMbyOeicAoBj7+FYYqDrGPgiC2zWOVCWg6/nqHk5M0mm8qYai7f0qsndYMC3xp1pjTu7Exyxx5LT7A8Akqmonam8iRkrA74lewO+HkfUHhxdeBAN1PAkZOCYywCy9Xc8uLqekwnQC+q/nwZkTJefqgde3pkEKSnwpSRIFt6jzI/W0KdDFoEMAnQ000Fr89vK0w3XElZVHP2gwNcP6q3trbaVe2yj/a503XObTh59RunpvD+JT4pbQw4a/CiCbu00fmoDhEojekFjeMqWehwCVpvcyKzOCZLlBoNYGaFc66ZNReX7sRX8MgCMzpIS6mLsssT8MXw3AOBVSVB5o/e23OzO+iTc7FMAuOvcjhgzk4BJLCQllzQzV3XdQeEgtppgfxEMExT0/N2ZnoPeERX0VnazDiLYKCpgw2QA21cWOR6QeHe1UdB3YMDl4WFPPdmi31V7KKzP4gCZFPg0rS84eLKrGLm2ocwGgDQ27amUrOzPbX+Asb/ovHDg22yq88LV3+baEUjJmGhnavUI9GybvYlg/w6U65QbDn7Ww+OjAP5dRF5J8i4AfwrgmdAeHj9pfXtB8i0AXgeNcvo5EXn/qTsXIG12MD5XTc3bmgxFvJva2IjfhMF04TNQGQ5CHqZpElrHLU7OClapgIirTaGnbtYNPVAaRAl+jqxzi51eiwTgQ2A7zXq/ySPo+bKfT5cw1iabgfpESEnQJY3LU+AbCvDtclJUcHLGZ+/QvNwV1EGigDdgCcERgSWJpalqHYiBghPRE7/hYOqyAuy4HaWzdQc8jJjgJANGfYBsCYHaKF43isDXseY0s7BAPRZ2xv4C6NUHC7d+kK7OeiaGq7m9xuuJxe0VtdclsWgm0nWgZGjJZ/1fAT2yAOI+RADIXuMGr19uBvN7E4CHATzZPt8H4K9E5G0k77PPv0TyuQBeBeA7oN3bPkjyOac2MRKAfb2wpaqy2dOE+iQFoKloXZjo9v365TDZ6Mwv3Pk2R7IIUuZoWwJBVZ4EKgHo3hP/PQNGWk27YkMKwCPG3Cbtfdcgjfc5fJ56b+63KRAM//KHBU2961LGgjUoecn46hvm58wsFcBrga+DaHtJZKyYsSRwRGIJBb4la8P6JIIMwQCNI1xpPz90yJPqdRwQ47hl+lxNaW6N2aEAn3qpFUPU0VJzmi3XG0BKFj7llWrchJDYPFim4j4bO7e/hqyML4StSAA/ChX0cq4ARwEWXcv6yKoe70OcLdxCskdeuy0knw7gRwD8Xlj9YwDeacvvBPDjYf27ReRERP4FwBcBvOjU/cNVE8tT3OqzgVKGftdTvN2ZNBO4qnNSgC7mshZGGJnQeLcx5tDBToAYEtPecNLedLtYxnlkPFY7wMmbfOqc7FgXb0QZ3fDlfjKWkwyEvLKKM7HObHhu91uyxwpDYHhDcWQo01PguzIBfAvU15IdlkzoSAUcA7w0cRKL+joCuZ1z5Lxi4werl5cUI1bKikmNCNBiDVLCpzw/OgLemHU34niSLVRFlAX6sgwDkO3lQJhPCWtx1hdZ4J5Eg7hPf91MuaHgB+B3ALwZbXLQU71vr7179YWnAfi3sN0jtm63mM1DY5sQYtakqpsmdV27Xv9puxurus0klBbk4qS8hgfkFvAA5wO5qwXDyPDC513bbQHBeWXiHHhmhgYpK4NzVpeQjZHpOmd5zvRW0LaSSwtzccbnfVATiSU7JCR0rC/dd7IYdh1EY0M86+RNnKdTvxLGvTV/gBbk7HzocVRmuLWvqfm04+GqPyyuP9fPDiZxrkdbjX/eOa79qr36e/ns102UG9m0/JUAvioiHyP50vN8ZWLd1tUh+XoAr7ePj3/oobd+HWdEch+o3I15XIckhzEugZKEzbm/cTeAb7/en30M//3+D8qf3X2OTW/aObyRNr+XAPhRkq8AcATgyST/CMBXSN4jIo+SvAeAl555BMAzwvefDuA/xju1Cg+lygPJj54nEfrQZB7XYcltPq5nXu9+ROTlezicvcoNU3tF5C0i8nQ7ca8C8CEReTWABwG81jZ7LYC/sOUHAbyK5BWS9wJ4NoAP36jjm2WWWS63XESc39ug5WZeB+BfAfwEAIjIZ0k+AOBz0GqtbzjV0zvLLLPMch1yU8BPRB4C8JAtfx3AD+zY7n4A91/l7icLHd4GMo/rsGQe14EJb2glh1lmmWWWW1RudKjLLLPMMsstKTP4zTLLLJdSDhb8SL6c5OdJftHS5A5GSL6D5Fet8KKvu4vkB0h+wd6fEv73Fhvn50n+0MUc9dlC8hkk/5rkwyQ/S/JNtv6gx0byiOSHSX7KxvWrtv6gx+VCsiP5CZLvs8+3xbjOFLEKrof0ggb6fwnAs6Dlqj8F4LkXfVxXcfzfD+AFAD4T1v0mgPts+T4Av2HLz7XxXQFwr427u+gx7BjXPQBeYMtPAvBPdvwHPTZoAP4TbXkJ7Rb24kMfVxjfLwB4F4D33S5z8TyvQ2V+LwLwRRH5ZxFZA3g3NDf4IERE/gbAf41W7y3n+aJERB4VkY/b8mPQghZPw4GPTVQet49LewkOfFzAjc+/v5XlUMHv6vOAb33ZX87zLSCjnqoHPzZTDT8JzUj6gIjcFuPCjc6/v4XlUMHvXHnAt4kc3FjHPVVP23Ri3S05NhEZROT50LTLF5H8zlM2P4hxxfz7835lYt0tN67zyqGC37nygA9MvmK5zriWnOdbRXb0VL0txgYAIvI/0ID9l+Pwx+X591+Gmo5eFvPvgYMd17nkUMHvIwCeTfJekito7vCDF3xM1ysHn/N8Sk/Vgx4byW8h+U22fAeAHwTwjzjwccllz7+/aI/Ltb4AvALqTfwSgLde9PFc5bH/CYBHoYWFHoGW7v9mAH8F4Av2flfY/q02zs8D+OGLPv5TxvW9UDXo0wA+aa9XHPrYADwPwCdsXJ8B8Cu2/qDHNRrjS1G9vbfNuE57zelts8wyy6WUQ1V7Z5llllmuS2bwm2WWWS6lzOA3yyyzXEqZwW+WWWa5lDKD3yyzzHIpZQa/WWaZ5VLKDH6zzDLLpZQZ/Ga54ULyhSQ/bXXxnmA18U7LjZ1llhsuc5DzLDdFSP4atH/zHQAeEZFfv+BDmuWSywx+s9wUsRzsjwA4BvA9MrclneWCZVZ7Z7lZcheAJ0IrPB9d8LHMMsvM/Ga5OULyQWjZpHsB3CMib7zgQ5rlkstNaVo+y+UWkj8FoBeRd5HsAPwdyZeJyIcu+thmubwyM79ZZpnlUsps85tlllkupczgN8sss1xKmcFvlllmuZQyg98ss8xyKWUGv1lmmeVSygx+s8wyy6WUGfxmmWWWSyn/D/sNDkAEIKhZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "plotDepthProfile(h, \"Random seafloor Depth Profile\")" ] }, { "cell_type": "code", "execution_count": 21, "id": "f1471ddf-9c13-4014-b6f5-cc57a1d568dc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 30\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 16.40 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=16.326802999999938, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.06265700000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "\n", "# Compile the operator again\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator\n", "op(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 22, "id": "b06f9ed5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling a tsunami in an ocean with random seafloor topography variations\")" ] }, { "cell_type": "markdown", "id": "275eb854-965a-4a85-934d-113d88c640fc", "metadata": {}, "source": [ "## Example VI: 2D circular dam break problem\n", "\n", "As a final modelling example, let's take a look at an (academic) engineering problem: a tsunami induced by the collapse of a circular dam in a lake with a constant bathymetry of 30 m. We only need to set the wave height in a circle with radius $r_0 = 5\\; m$ to $\\eta_0 = 0.5 \\; m$ and to zero everywhere else. To avoid the occurence of high frequency artifacts in the wavefield, known as numerical grid dispersion, we apply a Gaussian filter to the initial wave height. To achieve a symmetric dam collapse, the initial discharge fluxes $M_0,N_0$ are set to equal values." ] }, { "cell_type": "code", "execution_count": 23, "id": "2c3c9bb7-f37d-4c18-925b-eac7bfb2ef11", "metadata": {}, "outputs": [], "source": [ "# Define constant ocean depth profile h = 30 m\n", "h0 = 30. * np.ones_like(X)\n", "\n", "# Define initial eta [m]\n", "eta0 = np.zeros_like(X)\n", "\n", "# Define mask for circular dam location with radius r0\n", "r0 = 5.\n", "mask = np.where(np.sqrt((X-50)**2 + (Y-50)**2) <= r0)\n", "\n", "# Set wave height in dam to 0.5 m\n", "eta0[mask] = 0.5\n", "\n", "# Smooth dam boundaries with gaussian filter\n", "eta0 = gaussian_filter(eta0, sigma=8) # smooth random number perturbation\n", "\n", "# Define initial M and N\n", "M0 = 1. * eta0\n", "N0 = 1. * M0\n", "D0 = eta0 + h0\n", "\n", "# Maximum wave propagation time [s]\n", "Tmax = 3.\n", "dt = 1/4000.\n", "nt = (int)(Tmax/dt)" ] }, { "cell_type": "code", "execution_count": 24, "id": "d1fa80fd-f9bd-457b-a53f-78d81f45e1d3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The nt/nsnaps factor is 30\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Operator `Kernel` ran in 16.72 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=16.65377599999995, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n", " (PerfKey(name='section1', rank=None),\n", " PerfEntry(time=0.06264100000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "\n", "# Inserting initial conditions\n", "eta.data[0] = eta0.copy()\n", "M.data[0] = M0.copy()\n", "N.data[0] = N0.copy()\n", "D.data[:] = eta0 + h0\n", "h.data[:] = h0.copy()\n", "\n", "# Setting up function to save the snapshots\n", "factor = round(nt / nsnaps)\n", "print(f\"The nt/nsnaps factor is {factor}\")\n", "time_subsampled = ConditionalDimension('t_sub', parent=grid.time_dim, factor=factor)\n", "etasave = TimeFunction(name='etasave', grid=grid, space_order=2,\n", " save=nsnaps, time_dim=time_subsampled)\n", "\n", "# Compile the operator again\n", "op = ForwardOperator(etasave, eta, M, N, h, D, g, alpha, grid)\n", "\n", "# Use the operator\n", "op(eta=eta, etasave=etasave, M=M, N=N, D=D, h=h, time=nt-2, dt=dt)" ] }, { "cell_type": "code", "execution_count": 25, "id": "649a4227", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# NBVAL_SKIP\n", "snaps2video(etasave, \"Modeling a 2D circular dam break problem\")" ] }, { "cell_type": "markdown", "id": "9ce6aee8-b425-435c-9427-856797293673", "metadata": {}, "source": [ "## What we learned:\n", "\n", "* How to solve the 2D Shallow Water Equation using the FTCS finite difference scheme.\n", "\n", "* Propagation of (multiple) Tsunamis in an ocean with constant bathymetry.\n", "\n", "* Tsunamis reaching shallow waters near a coast will slow down, their wavelength decrease, while their wave height increases.\n", "\n", "* The influence of a seamount and roughness of the seafloor topoghraphy on the Tsunami wavefield.\n", "\n", "* Tsunami modelling for engineering problems, like the collapse of dams." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }