Front2D/IC_gen.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "ac48945a-00b5-4e5b-b3c2-f3b3935803a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "plt.style.use('/mnt/c/Ryan_Data/Python_Tool/ryan_mplstyle.mplstyle')\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import h5py\n",
    "import cmocean.cm as cmo\n",
    "import dedalus.public as d3\n",
    "\n",
    "from scipy import special"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "97e265d3-39f9-4490-88f3-11939bdf8221",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "Ly, Lz = 16, 1\n",
    "Ny, Nz = 512, 16\n",
    "\n",
    "dealias = 3/2\n",
    "dtype = np.float64\n",
    "\n",
    "Ro = 0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "07868419-285c-4847-8535-baef9d3d7fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bases\n",
    "coords = d3.CartesianCoordinates('y', 'z')\n",
    "dist = d3.Distributor(coords, dtype=dtype)\n",
    "ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)\n",
    "zbasis = d3.Chebyshev(coords['z'], size=Nz, bounds=(-Lz/2, Lz/2), dealias=dealias)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "e89404ff-63f9-4a33-bd9d-dfcc4758fa5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Substitutions\n",
    "dy = lambda A: d3.Differentiate(A, coords['y'])\n",
    "dz = lambda A: d3.Differentiate(A, coords['z'])\n",
    "\n",
    "y, z = dist.local_grids(ybasis, zbasis)\n",
    "\n",
    "lift_basisy = ybasis.derivative_basis(2)\n",
    "lifty = lambda A, n: d3.Lift(A, lift_basisy, n)\n",
    "lift_basisz = zbasis.derivative_basis(2)\n",
    "liftz = lambda A, n: d3.Lift(A, lift_basisz, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0ff79a4f-d352-4362-803b-6fd7857c0990",
   "metadata": {},
   "outputs": [],
   "source": [
    "Lz = 1\n",
    "\n",
    "dealias = 3/2\n",
    "timestepper = d3.RK111\n",
    "dtype = np.float64\n",
    "nu = 0\n",
    "\n",
    "Ro = 0.3\n",
    "\n",
    "# Bases\n",
    "coords = d3.CartesianCoordinates('y', 'z')\n",
    "dist = d3.Distributor(coords, dtype=dtype)\n",
    "ybasis = d3.Chebyshev(coords['y'], size=Ny, bounds=(-Ly/2, Ly/2), dealias=dealias)\n",
    "zbasis = d3.Chebyshev(coords['z'], size=Nz, bounds=(-Lz/2, Lz/2), dealias=dealias)\n",
    "\n",
    "# Substitutions\n",
    "dy = lambda A: d3.Differentiate(A, coords['y'])\n",
    "dz = lambda A: d3.Differentiate(A, coords['z'])\n",
    "\n",
    "y, z = dist.local_grids(ybasis, zbasis)\n",
    "\n",
    "lift_basisy = ybasis.derivative_basis(2)\n",
    "lifty = lambda A, n: d3.Lift(A, lift_basisy, n)\n",
    "lift_basisz = zbasis.derivative_basis(2)\n",
    "liftz = lambda A, n: d3.Lift(A, lift_basisz, n)\n",
    "\n",
    "# Fields\n",
    "q = dist.Field(name='q', bases=(ybasis,zbasis))\n",
    "\n",
    "P0 = dist.Field(name='P0', bases=(ybasis,zbasis))\n",
    "tau_y1 = dist.Field(name='tau_y1', bases=zbasis)\n",
    "tau_y2 = dist.Field(name='tau_y2', bases=zbasis)\n",
    "tau_z1 = dist.Field(name='tau_z1', bases=ybasis)\n",
    "tau_z2 = dist.Field(name='tau_z2', bases=ybasis)\n",
    "tau_P0 = dist.Field(name='tau_P0')\n",
    "\n",
    "P1 = dist.Field(name='P1', bases=(ybasis,zbasis))\n",
    "tau_P1y1 = dist.Field(name='tau_P1y1', bases=zbasis)\n",
    "tau_P1y2 = dist.Field(name='tau_P1y2', bases=zbasis)\n",
    "tau_P1z1 = dist.Field(name='tau_P1z1', bases=ybasis)\n",
    "tau_P1z2 = dist.Field(name='tau_P1z2', bases=ybasis)\n",
    "tau_P1 = dist.Field(name='tau_P1')\n",
    "\n",
    "G = dist.Field(name='G', bases=(ybasis,zbasis))\n",
    "tau_Gy1 = dist.Field(name='tau_Gy1', bases=zbasis)\n",
    "tau_Gy2 = dist.Field(name='tau_Gy2', bases=zbasis)\n",
    "tau_Gz1 = dist.Field(name='tau_Gz1', bases=ybasis)\n",
    "tau_Gz2 = dist.Field(name='tau_Gz2', bases=ybasis)\n",
    "\n",
    "bt = dist.Field(name='bt', bases=ybasis)\n",
    "bb = dist.Field(name='bb', bases=ybasis)\n",
    "\n",
    "u = dist.Field(name='u', bases=(ybasis,zbasis))\n",
    "b = dist.Field(name='b', bases=(ybasis,zbasis))\n",
    "zeta = dist.Field(name='zeta', bases=(ybasis,zbasis))\n",
    "\n",
    "C0 = dist.Field()\n",
    "C1 = dist.Field()\n",
    "D = dist.Field()\n",
    "\n",
    "# Substitution\n",
    "Q2 = lambda A: dy(dz(A))\n",
    "dyy = lambda A: dy(dy(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "ec7f72ca-1ab4-4172-8ea3-fd5defddd4ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'LHS': <Field 140290968281760>,\n",
       " 'RHS': Add(Diff(<Field 140291357154272>), Mul(0.3, Diff(<Field 140291222776160>))),\n",
       " 'condition': 'True',\n",
       " 'tensorsig': (),\n",
       " 'dtype': dtype('float64'),\n",
       " 'M': 0,\n",
       " 'L': Convert(<Field 140290968281760>),\n",
       " 'F': Add(Diff(<Field 140291357154272>), Mul(0.3, Diff(<Field 140291222776160>))),\n",
       " 'domain': <dedalus.core.domain.Domain at 0x7f9809546ab0>,\n",
       " 'matrix_dependence': array([False,  True]),\n",
       " 'matrix_coupling': array([False,  True])}"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Problem\n",
    "# Hacky way of using IVP to solve BVP, just don't put time derivatives\n",
    "problem = d3.IVP([P0, tau_z1, tau_z2, tau_y1, tau_y2, D, \\\n",
    "                  G, tau_Gz1, tau_Gz2, tau_Gy1, tau_Gy2, \\\n",
    "                  P1, tau_P1z1, tau_P1z2, tau_P1y1, tau_P1y2, C1, \\\n",
    "                  u, zeta, b\n",
    "                 ], namespace=locals())\n",
    "\n",
    "problem.add_equation(\"lap(P0) + liftz(tau_z1,-1) + liftz(tau_z2,-2)  + lifty(tau_y1,-1) + lifty(tau_y2,-2) = 0\")\n",
    "problem.add_equation(\"integ(P0) = 0\")\n",
    "problem.add_equation(\"dz(P0)(z=Lz/2) +D = bt\")\n",
    "problem.add_equation(\"dz(P0)(z=-Lz/2)-D = bb\")\n",
    "problem.add_equation(\"dy(P0)(y=-Ly/2) = 0\")\n",
    "problem.add_equation(\"dy(P0)(y=Ly/2) = 0\")\n",
    "\n",
    "problem.add_equation(\"lap(G) + liftz(tau_Gz1,-1) + liftz(tau_Gz2,-2)  + lifty(tau_Gy1,-1) + lifty(tau_Gy2,-2) -2*Q2(P0) = 0\")\n",
    "problem.add_equation(\"G(z=-Lz/2) = 0\")\n",
    "problem.add_equation(\"G(z=Lz/2) = 0\")\n",
    "problem.add_equation(\"G(y=-Ly/2) = 0\")\n",
    "problem.add_equation(\"G(y=Ly/2) = 0\")\n",
    "\n",
    "problem.add_equation(\"lap(P1) + liftz(tau_P1z1,-1) + liftz(tau_P1z2,-2)  + lifty(tau_P1y1,-1) + lifty(tau_P1y2,-2) - C1/Ro = -( dyy(P0)*dz(dz(P0))-dy(dz(P0))**2 )\")\n",
    "problem.add_equation(\"integ(P1) = 0\")\n",
    "problem.add_equation(\"dz(P1)(z=Lz/2) = D/Ro\")\n",
    "problem.add_equation(\"dz(P1)(z=-Lz/2) = -D/Ro\")\n",
    "problem.add_equation(\"dy(P1)(y=-Ly/2) = 0\")\n",
    "problem.add_equation(\"dy(P1)(y=Ly/2) = 0\")\n",
    "\n",
    "problem.add_equation(\"u = -(dy(P0)+Ro*dy(P1))\")\n",
    "problem.add_equation(\"zeta = -dy(u)\")\n",
    "problem.add_equation(\"b = dz(P0)+Ro*dz(P1)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "6760f02f-cd38-442a-b41a-11e793391bd5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2024-06-21 15:35:01,223 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 1s, Remaining: 0s, Rate: 1.9e+00/s\n"
     ]
    }
   ],
   "source": [
    "# Solver\n",
    "solver = problem.build_solver(timestepper)\n",
    "solver.stop_sim_time = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "489cd777-ccc1-4a02-a387-d1cd2a208bad",
   "metadata": {},
   "outputs": [],
   "source": [
    "q.change_scales(1)\n",
    "bt.change_scales(1)\n",
    "bb.change_scales(1)\n",
    "\n",
    "q['g'] = 0\n",
    "bt['g'] = -special.erf(1.0*y)\n",
    "bb['g'] = -special.erf(1.0*y)\n",
    "\n",
    "###\n",
    "snapshots = solver.evaluator.add_file_handler(\"front_IC_Ly_%d\" % Ly, sim_dt=1e-10, max_writes=10)\n",
    "# snapshots.add_task(-(-q), name='q')\n",
    "snapshots.add_task(-(-b), name='b')\n",
    "snapshots.add_task(-(-bt), name='bt')\n",
    "snapshots.add_task(-(-bb), name='bb')\n",
    "\n",
    "solver.step(1e-10); solver.step(1e-10); solver.step(1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "9ce35a18-f07f-4267-895e-62d839249964",
   "metadata": {},
   "outputs": [],
   "source": [
    "# b.change_scales(1)\n",
    "\n",
    "# fig=plt.figure(figsize=(3.2, 2.0))\n",
    "# # plt.contour(y.ravel(), z.ravel(), zeta.allgather_data('g').T, np.arange(-3.0,3.1,0.1), cmap=cmo.balance, norm=colors.CenteredNorm())\n",
    "# plt.contour(y.ravel(), z.ravel(), (b.allgather_data('g')+z).T, 13, cmap=cmo.balance)\n",
    "# # plt.gca().set_aspect('equal')\n",
    "# # plt.xlabel('$y$')\n",
    "# # plt.ylabel('$z$')\n",
    "# plt.xticks([])\n",
    "# plt.yticks([])\n",
    "# plt.title(\"$b$\")\n",
    "# plt.colorbar()\n",
    "# plt.xlim([-2,2])\n",
    "\n",
    "# figname = 'plot/figs/2DFQGpl_bt0'\n",
    "# figname = figname.replace(\".\", \"d\" ); figname = figname+'.pdf'\n",
    "# # plt.savefig(figname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "39b38f43-9573-4bb3-ba4c-5f96d5f8736c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAADpCAYAAAC0n/tpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABthElEQVR4nO2dd3hT1RvHP91770VLC2WvFpQhu+ACQSyCivxEpagoLmQ4QEApBUQFRVpwC4oURVQcLaIgIqNlb1pK924606RJ7u+P2trd3DSlw/t5Hp/HJPecvElzv7znnHcYCIIgICEhIdGJMWxrAyQkJCRaG0noJCQkOj2S0ElISHR6JKGTkJDo9EhCJyEh0emRhE5CQqLTIwmdhIREp0cSOgkJiU6PJHQdlPj4eIKDgzEwMGDt2rXIZLJmx0yYMIHo6GidX68iMTGR+Ph4Mea2mMWLFxMcHFz9WF+fReI/giDRYYmMjBTE/AljYmKEgoKC6sdxcXG1Htd9vSHi4uIEf39/YdGiRSKtbRm7du0S/P39qx/r47NI/HeQPLr/ECEhIdjb2wMgk8mYPn16o683RlBQEEFBQa1kYePUtUsfn0Xiv4MkdJ2I2NhYgoODiY6OZvr06Tg4OBAbGwtULnWnT59OVFRU9bWJiYmEh4cTHx9f73WAtWvXEh0dzbx581i7dm2z7x8fH8/atWtZu3Yt06dPr15OVz0fFRXF9OnTSUxMbNbeqnGLFy9m7dq1REZG1npe7GfR1QaJTkJbu5QSutPQ0hUQIiIiBEEQhEWLFglBQUHVrwUFBVW/VnVtzeVdzdcTEhIEe3t7QRAEoaCgoNb7hIaG1lu6FhQUCCEhIbXmioyMFAoKCqrnEYTKJWXNx83ZW0VEREStpauYz9ISGyQ6B5JH18mwt7evXloOGTKk1iGFv79/k2Nrvu7v709cXBwAJ06cAGjywCMqKooJEyZUP46LiyMsLIyoqCgGDx5c/XxISAhA9UFBY/bWHVd3uSzms+hqg0TnQRK6TkxL96iq9r6qlnlNkZCQoPXz/v7+Dc5Z094qkdUHutog0XmQhO4/Tn5+foPPJyYmMn36dHbt2kVYWFiz89jb2xMTE1PrOZlMRkBAQIOC0tyBhr29vVYCW5PGPouuNkh0HiSh68TUXYI1tCSr+VzN/6+5Ia+N4MyYMYPY2NjqA4Do6GhOnDhBWFhYrbg7mUyGTCarXj42ZsuECROIjY2tHhcfH19LyMR8Fl1tkOhEtPUmoYRuxMXFCUFBQdUb6QUFBUJMTIwACGFhYUJBQYEQGhoqAEJcXJwQFxcn2NvbCyEhIdWb9mFhYUJQUJAQExNT7/WCggLB399fCAoKEnbt2iUEBQUJoaGh1dcFBQUJCQkJtWyKjIwU7O3tBXt7eyEyMrKWraGhoUJkZGS1rYIgNGmvIPx7ABEaGiqEhYUJISEhDdra3GdpiQ0SnQMDQZBKqUtISHRupKWrhIREp0cSOgkJiU6PJHQSEhKdHknoJCQkOj2S0ElISHR6JKGTkJDo9EhCJyEh0ekxbmsD2hMajYb09HRsbGwwMDBoa3MkJJpFEASKi4vx9PTE0LB5v6W8vBylUtno66amppibm+vTxHaBJHQ1SE9Px8fHp63NkJAQTUpKCt7e3k1eU15ejoOlFeWCptFr3N3duX79eqcTO0noamBjYwPATy/9Dysz0xbNVXThtM5j887lih5zJaP2n1IQBC6rSogtz8HSwIi7LNxwN6r94z0rlGo1d5FayV/yXHxNrOlpqr23m6As4aqyiLlWftgbmjR7faCHiktyOW+lZPCslzt9rSxrve7U17nJ8Ta9+rPt9xN8F3eB8Bm308/Hvfo16x4Dmx7bc1Cz9hUaWzX+WmEhp06d4vSpeM6fO8+lixcoLy/H188PR08f3L274Orpiau7F45u7tjZO2Bg2PJVg7yslGfvv6f6t9sUSqWSckHDnVYemBjU9/4qBA0/ZWagVColoevMVN3AVmamWJu3TOjUJrp/tVaD3Mk9ky1qjLmBUfX/p6jk/KLIIlejZIK5K4NM7DBsQJyGGNgCcEooaXTePLWCw/JcAk1t6Wlmq7U9aRVlXFUW8Yilbz2BbYwCVTmb0rOY4+7CLbbW9V63auI7te4zgIgfDnLo0g0+nhdKVxeH2q9bWjT53jbW9d+vLmqTf69RqVScOhnPHwcOcPjwYS5dvICvnx/BwYOZMvVeXnltGd26d8fU1JQTKTnNzt1SxGy1mBgYNih0nRlJ6DoJPb0q+D1ZzQFFLtdUJQwzdWSWpU8tAWyMgQbWDYpdpkrOUXkefc3sCTBtXgiqyFMriC/P5z4LL3yMmxaYKjzcynklKZ3x9raEONjVe925v2ujY1UaDa/uiuVyRg4fhU3Dw762d2PdK7iRkZXY9h7c5OtVqNVqjv59hJ9+/JGYmF8xwIBRY8Yw57HHGDpsOC4uLlrNI3HzkYSuE3BVXs53eQXElZYx2NSe56wDsNFiqViTumJ3o6KUk+UFBJs74GPS+JKtLl4aE36SpzPR3JXeJs0vpwD8PRSsTM7A39ycB12cRNmtVGtYcz6LtIIiPpw7DQcr7YRVDDdSU9kR/Q1f7fkBjUbN7XfexfsfbGHgoCCMjJr+h+RmeHMSzSMJXSth128QhWdP6jzeub9rk8tXtSBworiUfQUyEuQKJjjY8V43X3Kydb/RBxpYIwgCXynSuaIsZriFM67GzS87BxpUentFmgq2lt1gsIk9t5o6avWegZ5K3k3PRoPAfE/XBpfYjXlzSrWGZYcvUWJoSuSjU7GzrG9rc95cYwiCwOFjx4n69HN+//MvRo8dy8o332TEbSMxNpZum9ZCJpO1SpVn6S/WwcirUPGbrIjfZEUIwB0Odiz09sDmH8+iJf6DStCwpzyTVFUZYVZ+eNTYWzsllFQLWkOUCWo+LUvB39iSEDPtl3Bf5eRzTV7Oaj9vzLQIj6iiSuTKjM3Y/L/J2JibaT22JnWXrYIgEPP7Qd7ZEkXijRvMnjGd8NdewapLN53m76wkJiYSERFBcHAwCQkJRERENHl9dHQ04eHh7N+/v56Q1awAHRISUq9StT6QhK4DUKpWc6y4lEOFxVwokzPA2pJH3Z0JsrbCqI4H1NOrgktp4patAMUaFV/JU1ELAk9Y+dVb+jYlckpBwxdlKTgamnCPuYfWG+Nplrn8ml3Im37e2DXiJTXkzVWJXLFSxQePheoscnX58+gxVm94l+S0NJ569BFmz5iOtVXlsl2ml3foPEyYMIFdu3YRFBREbGwsEyZMaFKgQkJC6vXehcpK1hEREdXVnlurZ4ckdO2UfLmSo1YCv13K4ExpGT6mZoy0s+ZpTzccW3Ci2xCpajlflqXiZ2zJFHMPTEWcyKkEDTvKUjHGgPstvOoJb2Mo7Ar5ODWXpT4eeDUSytOQyFWoNaz46zKFChWRT81sUuS0PYS4mnidleve4mjcSZ5+fA6PPfQgVjVCW2Qm2h/E/BeIjY0lPz+/uudGSEgIEyZMIDExsdHubI0JWGRkJEOGDCExMbFVe3hIQtdOkKvUnM8t5mR2IcczZCQUltLL0YYBFhY84uaMh6n24S7aenWCIHCiQsZP5VmMM3NhhKmjqDAFtSDwtTwdhaDmEasuWocsWDiVsPpGJmEervSpEyvX5PtpBN74+wo5ciVvjenTYk+utKyMtz+I4sMvdjBz2lTefnMlzo7a7S1qQ2c9iIiPj6/VPhIqu6rFxsZq1UipJjKZjMWLFwOVvT1qNirXJx1C6MTuB1QRHBzc4J5AW6MRBNJLyrlSUMqF3GIu5BVzVVaKk7kJA1ztuL+nJ4Pd7LE3NxEdT1dFc2KnENR8X55JgqqUhy196NpEMGxjn2G3PJ18jZJHrXwx0yKMBcDFVc4rSRlMcbJnlF3jp7J1vTm1RiD86FVSi8vZMLYP3sFDmnyf5ry5P0+d47WnX8HT3Y0fv/qC3j0CtbL/v0BRUVGtx2ZmZpiZ/fuPSkJCQr17yt7evtGWl01RtdyNiopi3rx5BAQEsGjRIvFGN0OHEDqx+wFQ+cVVdX1qK6z7DODq0aOklpRzo6iM5CI5ibIyrheWUaHR4G9nSQ9Ha+7t7kE/FxvcreqfGjZ3+toUjYldmlrO1/I07AxMeNKqK7YiQ1E0gsC35Rlkasp5zNIXSy1Frou7gmU30ulvZcE0J4dGr6srcoIg8HZcIlcKSnhnXF/szMTvQVZRUiZn9Sdf8fPf8Sx76QVmTb+vyRzRzrhs7Wdg1WB8ZTlq9kK9NMjly5fz+uuvt6pNYWFhyGQydu7c+d8UOl32A1q7ZZ0gCJQqKpCVyckvlZNfUkZOURk5xaVkFZaQVVhCuqyIzMISNBoNHlbmdLG1oIutBfd0c8PfzgpfWwtMjG5udLpaEDikzOOgIpfRZs6MNHVqMJyjKTSCwJ7yDFLUch6z7IKVoXY/oQBPJeEpGdgbGxPm4ar1ElkQBLacvsGJTBkbx/fF0dwUu35Np2s15s2dvHyNF9+JwsfNhQPf7cbHy1MrG3ShIy9bU1JSsLX9NwumpjcHDffJlclkDBnStJfdHKGhof/dpasu+wHh4eFEREQwb948nd4z8rdjqDUaypQVlCkqKFEoKSlXUqJQUiwvp1CuQKXWYGRogIOVBU7WljjbWOJsbYWHvQ2DfD3wcLDB094Wy9RrGIsIm9A3VV5dllrBHnk6CjQ8ZuWLl5H4eLsqTy71H5HTJii5p1cFgiDwXkY2JSoNK3y9MG5C5Op6c9svphFzI4dN4/riamnWrMg1aLdGQ+S3+9iy+weenXkvCxYt1arSx38VW1vbWkJXl6CgoHqCpK/DhNY6kGj3Qid2PyA2NpYZM2ZoNbdCoUChUFQ/rtqbKFEosDU3x93OHCszk+rcV2tzM+wszLCxMMPB0gIrM1MMm0nMLkxv+Q3VkuWrUqPhtGkm3+XKuNXUkXFmzjrlOaoFgWh5Olmach4VIXIAO3LyuFgm500/byya8GLritzea5l8fTmdt8f2wctGO2Gu680VFJew8N0oktKz2L5qCX0D/LQWufaybFWUl1NWUkxZSTHyslLkpaWUy+WUy+WUFBXedHuqQkGqVlWxsbEEBQVVr7Di4+Oxt7evteKqWmXl5+dX38+JiYkkJiZWzxcZGcnSpUtbxeZ2L3RiiYmJ0fqwIjw8nBUrVtR7/sU7R7Y4qb+KlmZItISTJaV8mJmLpZEhb/h5oczXLiWrLhWChp3yNIo0FTxm6avVcrVK5H7Kl7FfVsQqX28cmsgoqCtyvyXnsuV0EmtH9ybAvvKgRKw3dzEpmfkR79Graxf2rFuOjYgT3pagzbJVEARKiorIzUwnNyuT/Jxs8nOykeXnUpifT5GsgOJCGSVFRVQoK/8xNjM3x9LaBnMLS8wtLTEzt8DERPf9ypZQdZ8FBwcTFxfH/v37q18LDw9nyJAh1XttMpmMqKgooDJwOCwsDHt7exITE5k+fXr1dtSMGTNazaNr9w2s165dS0xMTK3Dh4CAACIiIggNDa13bdWXCJUVHQoKCho9dW3Io/Px8eHgq3P1JnSA3oROW68uTaHks+xcLpeV84CrEyH2trXi28QEFJcLanaUpaJC4GFLHyyaOXioEjiAw0XFbMnIZlkXL7pbNJ5KVlfkjmUUsPyvy6wY3oNbPP49tBCzN/fj4WO8uvljwu69iyfum1S9J6htAn9LvLmaQicIArlZGdy4dpW0pETSkq6TkXKDzLQU5KWl2NjZ4+TqhpOrGw7OLtg7OWPn6IStvQO29g5Y29phZWODpbVNg6lnZaUlhE0aT2FhYZPLTaj8fdvZ2fGKTWDDhxGCmjeLr2g1V0ej3Xt0YvYDdu7cSXh4eK3nunbtytKlSxs8yal7bN7RkalU7M4tYL+siHH2tswPcMPWuP4PWts4u1KNis/KUrAyMGKWpU+TgcQ1BQ7gdEkZH6Rn86K3e5MiV5cLucW8/tdlFg4OECVyVWg0GjZ9vZdPf4xhw/PzGDt4oNbvrQ+Ki4s4fewIV86eJuHieRIvX0QhL8Ojix/eXf3x9vMn+LbRuHv74OrphYWluLAeCd1o90InZj8gLi6u1lgDAwOuX7/e7uLodKWxvboilZof8mXsy5fR38qSiK4++DRTOLRKmBoTvDyNks9Kk/E0suA+C89GDxDqChxUVlN5Ky2DMA9XBlk3fSPX9OaSCstYeugic/v7Mt7333xZbUTOulcwCmUFi9/7kDNXE9n55st07+JV6xptvTkxqNVqTp86ycE//uCvw39y9uxZXNw9Cezbn8Ejx3D/3Kfw6RqAiYiAbwn90+6FDsTtB7RHWmufLq9CxY/5Mn4tKKSnpQXLfL0IFOE9QcOCl6qW80VZCgNM7LjdrOGKIg0JHECqQkl4SjozXZyaDAiG2iKXXaZg0R8XmNLNnXu7e4j6DFB56PBUxCZUKjW71ryKk53uS6/mlq0VFRX8dfhPfvn5J37757c4cuQoHpz1MJa+PXB0abx2nkTb0CGEzt/fv9H4ml27djU6rp1vP+qEc39Xjh5N5qeCQv4qKmGglaVOAleXKuE6UlTCJ+lZPODqxN2O9oBK6zmylBWsTE5jor0ddznaN3ltTZErVFTw0u8XGOrpwJy+tYNVtfHmipz9ePSVcPy93HnruTAsGtiO0Ic3d/7cOXZH72Lfjz9gYW7OxDvurFeXriPHz3VmOoTQSYC8Qs1vKbn8kJDF9YJSxtjbsq6rT6MJ8WIRBIFv8gr4Lq+A57zcGWwjbu8or0LFyuQ0htpYM8Ol6XzRmiInV6lZevAifnYWPBvkXyuQWBuRy7R2Z87Lqxk3eADLHp+FUQuDsOt6cwqFgh9/+J4dX3xOYkICd951Nxvf28zgIUPqhalIItd+kYTuJqHL8lWtEYjPLiT2Rg4HU/PwtrbgLn9XJozuTfmlfL3ZJtdo2JyeRUK5gpW+3viJTJYvVKlYlZxGH8vKAgRNZT3UFLkKtYblhy9jbmzIK0MDMRLZLOZ8ahbPfPEJD905ngUzpjT6vrp4c4WFhWz//DO2f/4Z9vYOzJo9m8n3TMFaiyY07Z1ADxWWRvU7gZWpNVDcBgbdBCSha2coVGrisws5nJbPodR8TAwNGNfFmY3j+tLd4V9vw7oFQcQ1SVMoeSstE1sjI8L9fLBr4JS2KYpValYmp+NnbsY8j4b386qoKXIaQSDi2DVk5RW8Pa4PpnU8sea8ubjraTz/5S8smDGFRyZNFGVzUxQWFvLRtq188dln9OnblzXr3uK2kSOlPr8dHEno2hhBEEgqkhOXJeNEpoxT2UU4mpswzNORVbf1pK+zjeh8VG356584twn2djzg6tRkalZDFKvVrEpOx8PUhGc83ZqsRVdT5ARB4P2TSVzKL2Hj+L71uns1J3J/XU3mpS9/4tXHZjE9ZFST12rrzcnLy9m0fRfboiLp07cvWz/6iKBg7T1BadnavpGE7iZi128QmSdPcK2glIt5xZzLK+ZsThHlKg19XWwZ4m7PvAF++NlaaOVB6JoaptBo+CQrlyNFJTzt6cYtNuKDY4vValbdSMfZxJhnvdy1FjmozF/9PSWX98b3w7FOYHZzInfwUhJLdv7Cm/MfZfLIoaLtrosgCHz30y+sWv829k5ObHx/M0OHDW/xvBLtC0noWgmNRiCzsJjE7HwSsvO5kpnH1cxcErPzsTMzoZejNX2dbZge6ElPR2udK5mIFbtEeTkb07OwNjJinb8PLjqkEBWrKj05ZxNjnvd2x0SEyO29lsnOS2m8M64vHtbiTor/uHidl3f9yrpnw7h9WPPeVnPe3NXE6yxd9SZXriXy/KLFTJl6r07J/pI31/6RhE5HNBqBgjI52UWlZBeWkCErJkNWTGp+ISn5RaTkyVBpNPg42ePv4kCguzMT+gbgXZSNs4WpXvd8tBE7tSDwXV4Bu3MLmOrkwDRnB63LntekUKVi5T/L1We9Ghe5hsqg/5GSywenk1g76t/81Zo05c0duJDIq9ExvBE6QSuRawqlsoKNW7fx/raPmTX9Pj7a+DYaR/GxexIdB0noGuDQ5euo1ALF5QqK5AqK5OUUlpUjkysoKCkjr0ROfmkZao2AnYUZrnbWeNjbVJZo8vPknqBedHG2x9vBFpN6m/v+rRI83JTYJZcr2JyRTblGw0pfLwJ0jLkrqKgUOV9z0yb35BoSuROZMtYcvcby4T3o51I/mLcpkfvtQiKv7ooh/P6J3D0ttNHratKYN3fq3HmeXfoqxsbG7PniEwb06d1uqpR0JPTVBUzX6uFikYSuAT45dBI7C3NszE2xsTDD1sKMADcn7CzNcbK2wMHKAmdrK5xsLDHXc6OallBX7JQaDd/mFbA3T8adjnbc7+yIqY512LKVFaxMTqenpTlPeriKErkLucUsO3yJFwYHMNSzfmXhZj25XTGEz5jI3ffqLnIVFRW8E7mVzR9+woKwx3n68Tl6qfzxX1226qsLmC7Vw3Wh/dyl7YgPH5+m1+olN5MqsTtdUsaHWTmYGxq2yIuDyhCUVcnpDLax4lE3Z1FNpq/LSlnyT/7qBL/6/V6bErk/Ll3nlV0xrL5/AqN7dtXZ/hupqTy1cAklpaXs3f4p/Xr3qn6tPXtz5XI5RQX5FMnyKS4srC7bVFpcRFlJMaUlJZS2QT06fXUB02UeXZGErg1ozRp1acVyNhfmcSqzkBkujtzuYKfTXlwVV+XlhKekE2JvxwMuDXcJa0zk0krKWfjHBe4L9BCdv3r4yg2W7vyVVaEhjOnl32yzmyrqenM//BrDC6++zrRJd7F80YtYmLcsVa4mLfXmigsLyUxNJis9lZyMdHIzM8jLziI/J5uCvBzkpaUYGRlhY++ArZ091nZ2lWWbrG2wsrHF3skFI2NjTh87oqdPpB366gKmz25izSEJXSdBVl7BZxdS+DExm4m+LmyfHIzqSkGL5jxZUsqGtExmulTlvdamMYEDyClTsPD384zv4szs3t4NXtOYN3csIZVFX/3M8nvHMb5PgNYiV5OKigreeOsdvvxmDxveeJ1JEyfUu+ZmeXMVSiXJiddIunKJ5IRrlXXpblynuFCGnYMjbl7euHp44ezuQbc+/XB0ccXxn9p0Vja2TZ4El5WW8FXkJr3ae7O6gOmzm1hzSELXRujLqytSVLDzcjrfXs1gkKsdkRP642f3TyXdf4RIl1i7/bIiPs7M4QkPV26rU4WkKYEDKChXsvD3CwS52vHkQL8GvcDGRO7UjQxe2L6PpZNHc3v/7qJsrvLmcnLzmPvcixSVlPDz1zvw9/Otd21LRa4pb664UMbFU3FcPnuaq+fPcuPaFcwtLPHr3oMuAd0ZdeckvPz88fDpgqXVzV86O/V1rhekDVBaoYIriW3SBay1kYSug5InV7LrSjp7r2XSx8mG9aP70Nu54TxMZxGCpxEEdubk80tBIUt9POljZVFrjuYoUlSw8PcLdHOw4oXBAaJE7kJaNgs++4Hn7hjOpEE9geb7s1ZRJXKnz19gztPPckvQIHZs/QBLC/FNgMSiUqm4cvYUZ479zZljR0i5noCXnz89+w/kjvtm0q13X1w8PDtMGtnN6gLWWt3EGkISujZEF6/uWkEpu6+ksz85lyHu9k0KXF3qilVd4VNoNLyXnsX1cgXv394fX1txPRZKlCoWH7yIh7U5S2/t1mCSfmMidyUzl6c+2UvYuCGE3tIX0F7kqvjh1xgWLH2V55+Yx9OPz2lUWPThzakqKjh74ihHf9/PySN/YmxsTP9bhzH5of/RJ2gItvaN961t79ysLmCt2U2sLpLQdQAUKjUHU/P57lom12SlTPRzYdvtA+giUojqUlP4skoVrDh8CQtzIyLHDxLdJLqsQs2SgxexMzNm2bDABls8NiZyN3JlPPXxXmYNH8isEQMBcSJn0yuYTVs/5J0tW/lg3RpuHzdWlO3aIggCX8ce5OBPP/D3gRjMLS25dcx4FoZvIKBXn/9MC0V9dQFrbh59IgldG9OYV6cRBM7lFhOTlMOBlFycLUyZFOBG+Khe2Jjq9892KruQFX9dZpS3E08P6io6HU1eoWbpoYuYGxuyckTPepVIoHGRSy8o4omPvuOeoF48Nkb8oYNl4EAWvb6KmN8P8t0Xn9C3V88mr9fFmysrK+P7vd/x5fYvSLpxg2FjJ/Dim+vp1qfff0bc6qKPLmDNzaNP2n0XsJtJVZckfXcBa44qoVNrBC7kFXMwNY+DqXmUqzSM9nHijq6u9HK01vsejyAI7LyczqfnUnhqkB+TA9xFzyGvULPk0EWMDAxYPbIn5g2UeWpM5HKKSnls2zeMCPRl0d3/lkLS1psrK1ewMOpLUtPT2R65GS+Ppu0XK3LZ2dl88dmn7PxyB15e3gy7awrDxk3E/Cbs+2mLLl3Afph2S6OHEZO+OSZ1AZPQPwWlco4J1hw4dpKjGQUIAozwcuS5YH8Gu9nrnOzfHEWKCtYeT+BaQSlvj+1DTyfxBSVbInIFpXKe/OQ7gvw8eeku8SJXUFzCU29vw9TUlL3bP8VWjwUx09PT2BYVxe5dXzNs+Ag2bd7CkFtuIS41V2/vIXFzkYTuJpNXUsaZ5EziktKJu57G1cw8urs7Mbx7F1b5u9HbyUZ0pV2xnM8tZtWRywTYWxE5sb/o/Tio2pO7gImRIW/eJk7kiuUK5n+ylwBXR16bOhZDQ3Eil56Tx+MR79OjWwDvr1uDuRYtK7Xx5rKzs4n8YDPRX+9k3PgQvt79DT16VmZR/FdTvToLHULoxCT+xsbGMm/ePBITEwkNDW2yeU5rk19SxpXMPC5n5HAhLYcLadmky4oIcHVkoK8Hj4wKYnBXL5ysKw8VWitbogq1RmD7xVR2XEzjsX5dCA300Gk5XKJUseTgRSxMjHhjRA/M6ohcU2ldZQolCz7/AWdbK94InYDRP3tc2opcQmoGc1auZ2JICOGvvVzdlKYpmhO50tJSPtq2lY8/3MaI225j17d7CAzsUf26JHIdnw4hdNom/spkMnbt2kVMTAwymYzx48czb968RjuI6YNShbK6PFNqfhE3cmUk5RSQkJ2PrKwcbwdbuns409vThanBvejr7YaNRcMeSGunhq05do1ChYqN4/oS6KhbiEWRooJFf1zA3tykwYOHpkROUaHi+e37MDEyYu3MOxqo7NI05xNv8Oiqt/jfAw+w+NmnW7xnKQgCP+z9jrURa/D29uHDTz5jUCuENki0Pe1e6MQk/sbGxtYStaVLl7Jz507R76lUqcgpqqBQXo6stBxZmZz8Ujl5xWXklpSRU1xZgy6rsIRCuQJLUxO8HG3xdrTD18meSYN6EuDmSFcXB6xEdunSt9hpBIE9VzPZevYGd/i5Mm+Ab4PLTG3IL1fy0u8X8LI259VhgaJErkKtZtFXv6CoULP5kcm1qr5o483FXbpK2Jvv8MLTT/HknP9pbXNj3ty1q1d5fdmrJCcns3jpy9w9aXKDwil5c52Ddi90YhJ/Q0Nrl/GpG8ujLSFrPgbAyswEe0sL7C3NcbS2xNHKAmcbSwLdnXG3s66uQ2dnYabXE1F9id31wjLeOp5AfrmS1bf1YpCbnc5zZZVW5q72cLRmya3d6sXJNSVyao2GV3fFklNUQuSjU7GsIf7aiNxfZy4wP2ITy5csYvaM+qV+GqMhkVMqFGz5YDMfbdvKAw/NInLbR1hZNdzasbOKnG3vAQ1GFRiVK+GbY21gUevT7oWuJYm/MTExzJs3r9HXFQoFCoWi+nFVMvN3zz+Eu70NJlrs/7QWLRE7uUrNFxdS2X0lg6nd3flfHx8sdPTiAG4UlbHojwsM83RkQVDXWmWamuvxoNEIrPz2AAnZeWx97N5ay3ZtRO7gybMsWL+Z1+c+LErkGuLC+fMsWbQQI0Mjtn/1NX369m3RfO2Bihq/X4nGafdCpyuJiYk4OjpWR183RHh4OCtWrKj3vIOVZZuKnK4IgsBvyblEnr6Bm5UZ743vRzcHcY2o63I+t5iXD11kSjd35vT1EdVgWhAEIn44yMkb6Xz4+DQcrP6NP9NG5H47forn397C6vmP8sDcp0TZXdObU6vVbI3cwpbN7/PY3DCeeGp+s0U328qbEwSB0uKi6nJNVtY2BPYbwMGff+TkXwcpKpRRXChjzF33cNf9D/H+m8vaxM6ORrsXOl0TfyMiIpo9hFi6dCkvvPBC9eOioqJ6lRvakioh0cazO5NTxAenksiVKwnr70uIb9ONpLXhSHo+q45cYW5/31r15JoTOKi8YTf8dJg/r9xg2+P34mL7r+BqI3Kxx07y4juRrH3mcaY/9oQou2uKXHp6GotefIHc3Fw+//Ir+vXr3+z4myFy5XI5yQlXSb+RREbKDYxNTZj+6BP8tOtLdnzwLpZW1ji6uHLrmPEE9huArb093fsOwNbBARtbOzx9/QB48uXXeSZ0Uqvb29Fp90KnS+JvVFQUixcvbnbuunW22itNLWOv5Jfw0blkzuYUM6OnJ/f38NT5sKEme65mEHnmBotv6cYYH+dqO7RBEATei/mbmHPX2PrYvXjY/xvMq43IxRyNZ+G7Uax/Noz75ogrwFhT5A78tp8li15i4sTbifrwYywtW5YbrCvFhTKuXTjH9SuXSE9OYv6rq7h24RwfrF6Ol68fHj5+dPWpLCU15u57GHP3PfXKNw0cOoKBQ0fUm9vMvP1kabRn2r3QiU0gjo6OZvDgwdWPExMTSUxMbHIJ2xGoK3YX8or54kIq8VmFTOnmztJbu+sU+FsXtUYg6swNfr6ezdpRvasb2WgrcgBbfjvG9/GXiHpsKj5O/x6AaOvJLXw3ireem8e0R+aK/wBUlk16d8Nb7Ni+nRVvvMGkyfdoPVYf3lxBXi7n44/Tf8hQjIyNeereO3Dz9sG/R28CevZBrVbTJ2gw7+/eV29sW9Sna0/IZLIGy663lHYvdKB9AnFsbGyDDTg6SzqvTZ+BHLqSxMe//MnVghLu6ebOwiEB9ZpA60pZhZo3/r5CanE5m0P60Xu4+EbO2w6cIPrYebY+NhU/l39LFWkjcvuPVy5XdRU5mYk1+Xl5vPDcAnJzctn17R5Rp+4tFbkfvvycw7E/k5qUiF/3Hvh0DcC3WyBb9sZgZa2/FLX2gJgg/uaurbk9FRIS8t9tjuPv79/oflvNzIeQkJBOI2o1kZWV88PJS3z991nKVSoeGNqfN6w1eq1iklZSzmt/XsLB3ITPFzyIrQ7NdLb9foIdR04T+ehU/F0dq5/XRuQOnDjFC29Hsu7ZuYTcor33WIXMxJqLFy4w/4kw+g8YyPtbohoNG2kIsSKn0Wi4dOYkf/8Wg5uXN3fPmIWVrS1TZs2hb/AQrG3/9WQ7m8iBuO5dTV0bGxtLRERE9YqrNbw56CBC919EpdZwLDGV7+MvceBiIn28XJk/cSjjevvXOxFuaczd8UwZbxxLYPKgnjx7+3CMdSgk8NEfcez46zRbHp1Cd3en6ue1DSF5/u1IIp55nIm3Bjfak7UxZCbW/PzTPl5espiweU8w78mnRB3EiBW5c3HHiIpYRYVSya1jxtOj30AAxt49RdQ8HRWxQfxNXRsZGcmQIUNareBmFZLQtSM0GoGzqVn8evYqMWevgQHcOSCQL+fPoGuNZWBdau6faSt6dv0GodEIfPjHCT45coWlk8cwaVCP5gc2wLYDlZ7cljlTCHR3rn5eG5E7fPo8C9Zv5s2n5nDHsMGiRa7A2IqtWz5gywebWbd+A+Mn1G+C0xTaiJxarSL+8CH+/HUfT726Cg8fX+Y8v5h+Q4ZibPzfu4XEBPE3d61MJqs+OAwLC2u1dM3/3l+pnVFeoSLuehoHLyXx+8XrKFQqxvX2Z2VoCEP8vaqT3rVF20ODwrJylu3ez/XsfD6aO40eHvV7rmrD1gPH+fLIGSIfnSrak/v77EWeXvseK8Jmc/eIW0SLnFJZwdJXFnHk8GG+2PEVvfv0EW1/c/wc/RU/RX8JwPjJ9yIIGpxc3XByddP7e7UX9NkFrLlrq5awUVFRzJs3j4CAgOqCnfpEErqbjFqj4UpmHicSUzmakEp8UjoOVhbc1sOXlfeNJ6irZ6sHK59JyWTJV7/Qw8OZL56artN+nCAIfLD/GLuPnyfqsal0cxMncscvXOHJNRt57bGHmDJ6mGiRKywqYvazL1FcXMyub/fg6qpd856aNObNFckKyMnMIKBnbxTlcmY99RzBI0Zh2AGDyBvCusdArC0bCEspkwO0SRewKu9u586dktB1RPJKyriYls251GzOpmZxJjkDAwMDBvl6Mrx7F164cwRdXRxuSocotUbDxwfj+eiPOJ4KuZWHhg/Q6X0FQWDTr3/zw8lLbH1M/MHDycvXCFv9Dksemcm0cbeJFrm0jExmhM3H28eHzZHiDh2qaEjkSouL+HHndn7ZvZMRE+4goGdvpsyaI3rujo4+u4CJuTY0NFRaurZ3iuUKbuTJuJ5dQGJOPtey8riamUd2USldnOzo6+3GyB6+LJg4lG5uTqKXpC0lvaCI5bv3k1NcyrbH76W3l3gPCCr3EdfvO8RvFxLZ+vi9+DrbV7+mjciduXadx994mxcfuo8ZE0aLFrmLV64yM2w+Y8aN5bXlK3TaI2tI5Arz81j0yEx8uwWyZP0muvfpJ3rezoI+u4CJDfhvrQMJSei0QKXWUFAmJ6eo9J8STaVkFhaTXlBMWkERaQVFFJTKsbc0p6uLA11dHLgt0I9HRgbTw8P5pvafqIsgCOyJu8iGn/7k9n7deffhu2tVDxGDWqPhze9+50RiGh8+Pg0vx8qbQduimReuJ/PYqg08ff8UZt05XrTIHYs/yawnn+GRRx/jyfm61aOrKXKCIBB3+CCW1tb0HhjM0rfew6+7bgcy/yXEBPE3dW3dYP7IyEiWLl3aKjZLQtcAr+yKoVShrKxDVyKnUF6OIIC9pTkutla42lrhbmdDNzdHxvTqipejLV2c7HTa62pNUvMLeWPP7yTlFrBmxu2MCKzfsV5blCo1r0XHcDUrn21z78XVtjKCX1uRu5KcypyV6wm79y7mTJ4oWuRi/zhE2IuLWLRkKTMfeFC0/VBb5DJTk/nk3XUkX7vKnOcrT/0kkdMeMV3AGrs2MTGR6dOnV4eczJgxo9U8OqkLWA2quiQ9PWEorrbW2FuaVdehc7K2FF0Rt62oUKnZ/tdptv5+nLsG9GDB7cOwMdc9p1eurGDhlz9TWFbOptmTqquQiCl/PmtZBLPuHMf86ZXpWGKEbs++n3n+1eWsXrOWO+++W/wHoLbIpScn8WrYbEbfOZnQx57o0AG9unQBi//8/QYPI0rK5AQ9PF/qAvZf4f5b+7XpcrMlHL2WwtofD2FkaMCm2ZMJ8vNs0XwFpXKe/fwHzE1M2DJnSvX3oq3IJaVn8b/X1zJjwiidRG5H9De8Gh7Bxvc2M3L0aPEfgH9F7vqVSxgaGtIloDurIj/Fy7erTvO1BVX+iFKh4PqVixQXFlJSVEhpcXEbW9YxkISuk3AjV8Y7Px8mLimdsLFDmDG0X4vDVNILinj60+/p5ubEG9MnYGpspLXAASRnZjP79bXcM2oYz868FxAncts+307ExveJ3PYxQ265RbT9UClyarWK7774hB++/IyHn34B326B7U7kKpRKcjLTyU5Po3uffpiambMlfAX5OdkU5OYgy89l064fkZeW8MGby7Gxs8fazg5L687lebUWktB1cHKKSok6cJzvT17inqCefPf8rFoFLnXlQlo2z37+AyF9u7HwrtswMjQUJXKp2bnMXr6OO4YN5qWHp2NgYCBK5N7b9hHvbfuIDz/9nP4DBujyETiRkkNxoYz1S1+gvKyMZZu2tvk+nEqlIiM5iRsJV/HtFohP1wBWLgjj6vmzGBub4OrhyZOvrMC3WyABvXoTfNsoHJxccHB2xsraGhs7O97dubd6vrLSEo79HtuGn6hjIEroTp06xcCBA1vJFAkx5BSX8snBeL49cYERgb58NX9GrWohLeGPS9d55esYnhh/S3WsnRiRS8vOZfbytYwbMpClj8wULXIbNkey7YsdfPz5Dnr17q3LR+BESg6K8nKsbGy5bcKdjL77HkxNb37twbzsLMrlZXj5dmXHB+8Ss2d39fJ56sNz8OkawOxnXsTG3h4HJxcMa4Qd3XX/Qzfd3s6KKKGbPn06S5YsYfz48fj5+bWSSRJNkZJXyOeHT/J9/CWGde/CR3On0dNTt/StugiCwPa/TvPB/qOsmDaekL7dRAkcQEZuPg8vX8uoQf147bEHRYmcIAis3bSZz7/exafbv6R7YKAuH4O/rqXw2aYNFObn8eLqt5hwb8t6TYjl+pVLxO6J5sKpOPKyspg4bTqz5j/PyNvvZuyke3Hz8q4laG3tZf4XECV0CQkJzJ07FwMDg+oYmfvvv5+xY8cCksfXWgiCQNz1dL48cpo/r9xgQt9ufPbE9Fq5pS1FqVKz5vs/+PPyDaIenUofbzfRIpeek8fDy9dy28C+LHv8IdEit+bdTXy5ew+fbv+KgG7ddPkY7Dt6indfX4qJiQnPLF+t0xxiSU++wYlDB7h++SILVqyhQqnE2taOR55bRGDf/lhYVmZu+Pjr9pn0jU3PQdhY1y/waVBS0gbW3BxECV1MTAzjx48nPj6enTt3sm7dOiIjIzEwMCAkJIT8/HyOHz/eWrb+5ygolbPv9BV2Hz9PQUkZ04b04aW7R+Jur99wiNziUl768mcq1Bo+f3I6AUPFn25Wi9yAPiyfOwtDQ0NRIhf+zkZ2fruXT3Z8pVOLSqhcrm5d9wY9+g3gwScWYNxMA5yWIi8rZeUzYWSmJNN38C0EjRiFRqMmsG9/Avs23JviSLp2p6TDPDtuyEt7RJTQjR8/HqhM0wgKCqKwsJCIiAhiYmLYuXMnsbHSpmhLUVSoOHw1mX2nLnPochJ9vN14bHQwIX0CMDPR/9nR6eRMFn35M4P9vXht6lic+98qeo607FweXr6WkQP76iRyq9/eyNd79vLx9i91EjlBEIj8ajfBI0bx4uoNmFu0Th8FRXk5xw/+xh8/fc+EqdO5ZfQ47n/8SXr0H9hgCXRtRa0hao6VRK/liLpzkpKSau3NGRgYYGdnR2hoKKGhoWzdulXf9v0nKK9QcfRaCvvPJ/D7pevYW5pzZ/9AFtw+jC5O9q3ynoIgEH3sPG//fJhnJg7l8TlzdEqpSs3O5eFlEYwO6s+yxx+66SJXXl7OU8+9wMVTcfj36IWji245vM3x+4972bFlIw5Ozoy+azI9B1SWwxo07LZa17VE3Bqjak5J8HRHlNDNmzeP6OhobGwa/sLnztWtmcl/DUEQSMkv5O9rKRy+kszxxFScrS0Z1yeAzY/cQx8v11atZlKmUPLmd79zLDGVTbMnM/quyTrNk5yZzezl6xg3ZKBOBw9VIvfJji/p2lW8yOXk5PDInEcxMDRgVeSn2Dvqb88S4NKZk5z6+y9mhs2nS7fuvLh6PYF961d8aQ1xa4gj6cWS2OmI6D26qmTd0NBQ8vPza71eVFTU6VJH9IFGI5CUW8CZ5Ezik9I5cT2NvJIyBnTxYHj3Ljx7+7CbVqrpSkYuS77+FWdrS/a+sxoXB7vmBzXA9fRMZi9fx53DB7P0kZnY9Wm6z25Nqg4eqjw5XUQOIDbuLF5+XXnk+UV6Cx0RBIEzx//muy8+JjUxgfFTpqFSqfDv0avWdTdL3OoiiZ1uiBK6iIgIgoKC+PXXX4mJiSE+Pp7o6OjqpNyEhAQ++OADvRupz45DrU2FWk1KXiFXMnK5nJnLxbQcLqbnUKFW08vThSA/T5bdO44BXdyxMG3dzfKaaDQCO4+eYdOvfzN32t08dd9kjHToDQFwNTmNR1asZ+qY4SycFSpK5ADWvbeZHbu/5dPtuh08HPnrMKeupzFk1Fi9llNSKsoxNjHlhy8/Y+Ctw1m4egOWNU4n20rc6tIexE5f9+TNul9FCd1LL70E/HsoUVhYSExMDF9//TWrV6+msLCwVYROXx2H9EWFWk1OUSnpBcWk5heSkl/IjVxZ5X95MowNDenm5kh3d2cm9OvGgtuH0d3dqdUrBzdGVmEJK775jaTcAj5avpDBvXSLTwM4n3iDOSvfYtad43jm/imiRe6t97fw2c5dfPLFl/gHBIh+/+++/ZZlr73CnBeWiB7bGDeuXeGrqPexsbXjqVdX8vKGzbW86/YicDVpa7HT1z15M+5X0HP1kiVLlrBmzRp9TQdQ3au1oKCg+jkDAwMSEhIa7Dik7bUNUVXd4d1Zd6NUqykorSzTlFdSRm7xP7XoikrJLS7DwABcba3xdrTF29EOP2d7ujjbE+DqiKe9LYaGrb8MbQ5BEPj+5CXe2neYCcMG8+qjDzZcQltL4i5dJezNd3jyvkk8PvVO0aWW3o3cSuSnX/DJ9h0EBooPkt0WFcn777/HsyvW0DdYt9zXmpSVlPDllo38GfMzE6bexz0PPVLdprA9iltDDLAzEF295OrxvxqMoysuKaH7kOHNzqWve7KqTJOu96sY9BqvoG+RA/12HNKWNT8cxMHKAgcrcxytLHG0tiDI2QsXG0vc7KxxtbXGxcZKp7aAN4u0/CLe3Ps7iXlFrH0ujPFDxPdKrcmhk+dYsP59Fs2+nwduHyta5N7/8GO2fPyZziJ3IiWHFFkJr7z9QYszCTRqNUqlAmNTEwwMDVnz0Q7cvLyBjiNwVRzPuPlBvvq6J2UymV7uV23QSugKCwsJDw/HwMCAGTNm3NTsB312HKqLQqFAoVBUP67qfrTz6ZkdtkxThUrN54dP8eHBeCaNHMp7y+7H1sqyRXPuO3yMlzd/zBtPPsKk224VLXKRn3zOpq0f8tFn2+nRs1fzA2qgUqlYvHwlE6fdz+QHZosa2xDJCVf5cH043fv2Z9b853j0nyVwRxO41uRmdQGTyWSi7teWoJXQ2dnZVXtru3fvZsuWLXTr1o2wsLAOfcoaHh7OihUr2toMvfHn5STe3n8CE2NjPlr2IsE9u7d4zs/37WfD9mjeffFJRgf1Fy1yH27/krc2b+HDTz8T3Y5QoVDw6LwnyEpL5Y7pD4gaWxdVRQXfffExP+78gtvvm8m9sx8F2q/AqSuUaDQaTMzMSTxxGFlmKmWyAuRFBQyaNAPnLv5sX/gIJfm5oucuNLZCbVJ/6VpiXLmL1RZdwFob0UvX++67j/vuu4/CwkIiIyNJTExkwoQJTJs2rTXsa7WOQwBLly7lhRdeqH5cVFRU74/cEbiSmcv7h89z5up1np05lZkTx2DcwoMPQRDYsOMbdv76O58sX8iAwADRIvfpV1+z5p1NbPvkM/r1F1dqqaysjNlzHkVeVsqr73xQvXemK4pyOYmXL1aXaqoUuIoWzdlS1BVKZBmpGBob4+DZhfgfdpIU9xeF2ZmUyfIY+b+n6T9xKjdO/o1CXoqlnSMOnl0wMa8s2T9+3iIwMGD38mf0atfN6gJW1TNCm3lais57dHZ2dtWnsCdPnmTJkiWtsrRtzY5DdV3yjkaevTfv7drLz0dO8MDEsax/NgwHm/r/UotFWaHi5c0fE3/pKl+++TIB3h6iRW77rt2sWv82UR9/wgAdfg9xKdl4+vox/bEnqpPixaLRaPg5+itk+bk8+MQCXlrzNkfSi8loAy9OKS+jtCAXB88uXD68n7jvdiDLSMHY1IyBd4Vyy33/w9HbDysHZ2xd3LF1ccPSrrKN5Ni5LzY4p3v33ijLSvVu683qAmZvby/qfm0JejmMGDRoEIMGVW5211zaLly4sMVz66vjUGcix86LqG/38f2hSCaPHMq+d97Ax00/pZoKS0p5Zt37FJfJ2bn6FVwc7ESL3I7d37I8Yj1btn1IUJC4CiiFhYWs3bSZKbPmMPuZhm9wrebJz2PLmpVkptzgyVdWtskStTA7nVM/RpNx5Tz5Kdfx7D2AqS+vw8mnK8NmPo6Tjz82zv9mwfgNFJ9n3Bbo656suv5m3K86C11jWRA1l7b6Qh8dhzo61r2COXUlgY/3/sL+E6eYPHIoP2xYhZ+nm97eIzkzm7lvvoO/lzsfLFmAlYW5aJHbuec7XlsdwZZtHzF4iLgQkIKCAmY++CCOLq6o1WoMdVx+a9Rq3nz+KboEdGPKqg/Is2y5l9scapWKzCvnSDp1lLTzp5j62gYMDAwxNDZmyL0P4d69N1YOzgA4d/HHuUvH/sdXX/fkzbpfdY6j6969OzExMbWS/Pfv309sbCzh4eH6su+mUrNLEjcutLU5AJgE9OfnIyfY/vNvXE1JY/r4kTwyaSKeLvrN6zxy9iLPrt/M1DEjWDz7foyMtE/Or2Lnnu94eVU4m6O2cevQoaLGFuTnM2PmTNy8vJn/2hs6lVgSBIHESxcI6NWHX05dwcbFvVXT6lRKJRq1ChMzcz577iHUKhW+A2+hS/8hdA0ahvFNqGisLCtl69wpouLoTpw8jXUD+eolxcUMHjRA6gJWE0EQqjcTx40bB1RmTMhkMn3Z1qZUFZ0suRjXJu99+UYqu387xJ6IT3GwsebB28cybewIbFoYKlIXQRDY/vNvrPt8F8sef4j7xo0ExDWxgZaJHEBs/Dm69ujFnBeWYGws/mcpLyslKmIVN65d5d43o7B19RA9hzYIGg0Zl89x6dCvJBw/xC3TZjPgzvuY8vJ67Fw9MDBsv7GV/2V0FjoDAwPGjRvHyZMn+eabb5g2bRrffPMNJ06c4L777tOnjW2Kda/gmyJ21r2CSc/JY9/hY+zdtpyk9CxuHxbMpoVPcUufHq3imSiUFby+9XP+iD9TKxxFrMh99e13vPJGOO9HbhUtcnl5ueyOPUjQ8JHMXfSqqLFVZKQk8/arL2Fo68ikZe9ibKr/GEiVUoGxqRlx33/F6X3RdB8+lsmL1+AW0BMAe3cvvb+nhP5o0anrqVOnGDRoEP7+/qxbt46oqCji4m6+B9TatIZ3Z90rGEEQuJKcxoG408R+uooLicnc2rcnj0yayMRbg1qUrtUc6Tl5LHhrMwDfrF2Ou1NlYx1dRe6zDzbRZ9gwUWPz8nKZMWMmvt0DCRo+UtTYmpyLO4Zbv1sYOvMxDA31l08sCAJZ1y5y+qfdpF8+x+x3t9Nvwj0Muns6RsY3ryCDRMvRWehOnDhR/f9VoSZ5eXmdbm1fk5o9FMSIXs1xmXkFHL9wmSPvf8Sfp89TWFLK8P59mDlhDONvGaSX8JDmOHTyHC++G8XtQ4N59dEHMfunioro09Xob3gtfG2lyI0YK2pslcj5+HcjbPFrosZCpQjF7ommwqc3DrfczvCWp77WQqNWs+fNF8lLTqTXmLuYtvxtjIyNMTJu/b+PhP5p97mu7RVtGscoKyq4fCOVsz8f4OSVa5y8nEBqdg69u3ZhaN9ehM9/lMG9AquFprWpUKl496s9fPHTfpbPncW9Y0ZUvyZW5D7/OpoVa9/iiy3v0Wv4GNG2bPtqN97+ATzx8nKMjMT9DNVqFZ9tfIu//jjA3S+uEv3ejSFoNFw7dhAjI2P8h9xG0OSZePbsj6mFfvdFJW4+UgNrPVChUpGancv1tEyupaZzNSWNS0kpJKRmYGVhTh9/XwYFBjDptlsJ6tFN7wcK2pCcmc1LG7dSKi8nes2rdPOp3FMSK3AA277YQcS777EjajOBt4pbchbk53M0IZXx90xj7KSptdr+aUO5XM7KpYsoyctm+sr3sHFueXiNIAjcOHmUIzu3oSwrZfiDlQnlfoPEH6p0BOLTcrG0Kq/3fFmp1AXsP4sgCBSWlJIrKyQrX0ZWvozMvHzScvJIy84lJSuH9Jw8jIwM8XV3o5u3B926eDHh1iB6+vrg4+ZyUyoHN2X/Nwf+5I2PvmTK6GEsmT0Dc7PKzXpdRO79Dz9mY9Q2dn4YiX+wuD05mUzGzAcfJLDfAP63YKFokQM4kSPHxa8bE59+GTM9xcflp1znt63rCZ7yEH3H342RSccs6CDROJLQNUBY+LsUl5ZRUFxCQVExFSo1luZmuNjb4e7kiJuTPZ7OTvQL8KOLuytd3F3xcHLUuWJva5GZV8CyyE85n3CDd198klGD/q3GK1bkBEFg/fsf8PH2r4j+eBs+A8TlIxYVFfHAQw/h7O7BQ089J2oswE/xF4nf+xWjH32Wofc/Knp8XcoKCzjy1Ta6Bg/Hf/AIZr+7o1VOayXaB5LQNcDdI4bg4eSEvY01zva2ONrZYN1KLfRaA41Gw679h1j72deMGzKQH99ZhX2NQw5dRO71tW+x58ef+Pbzj3HrJS5BH+CVlW9g5+DIM8veFB0nt/dwPN+vXUrPkRNbfKqq0ag5H/sDf+/6CJ8+QTj7VlY5lkSucyMJXQNMGTW8VUM7WpMryaksi/yMzNwC1j8Xxtjg2qIkVuRUKhUvLV/Fn0eP8t0Xn2AfIK6enFwu50yWjBlh8zE0MhKd8bD7wF/8uO4VBt87i0F33y9qbENkXDrLqZ93c/szr9Klv/6rZLQWFYpy5EUybF3ckRcXkhT/N+UlRZTJ8psfLCEJXWehsKSUjTv38HXsQWbdMZ5nXpuCpXntFCSxIleuUPDUwiUkJCWxd/unWHiL6/Egl8t5aPZs+gTdwpRZj4gaC5W14ixs7Bn5v6fpOXKi6PFVqFUq4r//Cq+e/fHqPZAH137Y7uLgSgtyKczKoCg7g6KcDPpNnIqBgQHfrHiO0oJclPIyDI2MmLvte8qLi7j8Zwzm1rYYd+DqOy2hoaKdTSEJXQdHWaFixy+/sXnX9/QN8GPPutcJ8K6f/iRW5AqLinjk6edQKBTs+fwTDFzERf4rFAoenjMHtVrN7ffNEDX2SHox147+QWlBPgPuuBcHT91rBBakJxP7wRoqFAr8BlVWB2lLkSuV5ZObdI28lETyUpLoGzIZj8A+7HzlSYyMjbF19cDW2R2VUoGVgxMjZj2BlYMzVg5OmFvZYGBoiIOnD1NfWQ9U5rpePqT/ZjL6Qp8d/GrWtgsJCRHVREcSug6KSq1m7x9HeG/XXizNzVj/XBgjB/bVywlvZnY2D8x9Ek93d77Y8h4Vdq6i53hywXOUl5WxeN1GzEXsbx5JL+bSoV/54+ON3P6MbilhVRSkJ/P1q0/RN2Qyt4bOuen7cBWKcrKuXSLz2gUcPH0IGDKSA1vXU5iVgVOXrjh5d8XcujLA/n8bv8Sogb1L3wF6joS+yeirW1hsbCwRERHVZZ/EeHMgCV2HQ1mhYu/BI0R9uw+NILBgxlQmjxza5ImvGG/u4pUrPDRvPqOGD2Pd669Raukgyj5BEIhLzeXumbNwdvcQVTTzSHoxF37/iT8/38xdL6zEp69uBRgVpSUoSoux9/Bh2rK3cfFreUl5bVCrVBTnZGLv4c2ZX/dw+IstWNja4R7YByefrgDc9eKqBg9UGhK5jk5sbCz5+fnVhTSr+j9X1Z8Tc21kZGR1VWJdCnN2vm+3k1JUWsbXMX/w+U/7MTc15Yn7JjF55K2YNHODiBG5Pw4fYe7zC3nikdk8/2SYaO9QpVLx+FPzGTZuIsEjRokaW1UY08DAgEkL38SzV39R46vITrzMLxtX4Rc0jJGz57e6yFWUy0k8cZjrcX+RcvYEdm5e3P/mB/gHj8B34C3YunjU+h71mYvb3tFnBz+ZTMbixYsBCAsLq1eZuDkkoWvnXEpKYccvB9h78Ai9unbh5TkPEDJkkFYxe2JE7pMvd7Ji3VtELHuV+6feA4CsgQYqjaFWq3n8qfmkJV0nsK84kary5Fz8utNr9B2ixtbkwoF9HPp8M4OnPkTQJHH7gmJQlpWSeOIwXYOHo1KWc2rfLroGjyDonpm4+HYDwNpJPxWf24LmuoBpiz47+FUtYaOiopg3bx4BAQHVhT21QRK6dkhhSSk//XWc6N8OceVGGpNG3soXKxfTN8BP6zm0FTmVSsXKdRuI/v5Hvtq6hVuDK5cFYkROEAReXryIlIRrvPz2Zmzs7LUeeyS9mDO/7uHo1x8xebHuudIajZrk08e5+8VVePdpWQ/bxki7cIrzv/3I9bi/sHP3xqVrd5x8ujJjtTjvoq05nlGCqWX9ertV/SfacxewKu9u586dktB1RErK5Pwed5p9fx3nj/iz9PTzIXTcbdx9262i+7JqK3KywiKeeHER6ZlZ7PvqC/y6VP7AxYgcVC43HbsGsmTW49g5OGo97kh6MWd//a5a5Ny79xb1vgAl+bmc2reL4Q+Eccdzy2u9FpeiaXZ8sE/TnrG8uJCi7AzcAnqSdPIo1o4uhK58r3rPrTPSXBewqKioJnuvTpgwgZCQkFbr4BcaGiotXTsSyZnZ/BF/lt/jT/P32Ut09XTj9mGDeenh6XT1dNdpTm1F7uKVqzz6zPN0D+jKj199jo11pbiJETmNRsM7G96ix20hjL9H+3aXVftxgiCQlXhJZ5HLSrjEvg3L6NJ/CIJGTVy6+BPnmmJYU/QK0pM5+eMurhzeT/fhY3EL6MmIh+aJnr8j0lwXsLr7a43Rmh38xB5ItHuhExOHExsby7x580hMTCQ0NJRdu3bdREubRhAE0nLyOHHhCscuXObvsxfJyi8guGd3Rgf359VHH9RZ3KrQVuT2/PgTLy5bwbxHHmbh/Cerk+vFitzy117lj0N/8sq4O7UeVyVy1+OP4N17ICFPLNZ6bE2ux/3Fr++vxivkf9jdNo1TmS0Pq6kSPUVBFmc2PEH34WOZvur9Tu29tSb66hZW1bKh6prIyEiWLl0qypZ2L3TaxuHIZDJ27dpFTEwMMpmM8ePHM2/ePNEurj4QBIGM3Hwu30jhfGIy5xOTOH01EVlxKX38fRnSO5Dlc2cxpHePetkLrUm5QsHyNevYs+8nPlgfwcSxo3WaR6PR8Pqy1ypF7p3NOLloVyqpSuQu/vEzhz57n3tf3YBLV3GnooIgYGBgQKrgTsDMl3HopZ9SSqVp10iJ+RSX4Ak49RvFoKU7uLWHvV7m/i+jj25hiYmJTJ8+vTrkZMaMGZ3LoxMbh1NT1JYuXcrOnTtbzTZBEJCVlJKalUNKVg43MrO5kZFFQmoGiWkZlJUr8PN0o3dXX4b07sHjU++kT1ff6hJJ+qY5b+7a9es8+eJijE2Mifnma7p41c50EOPNKZVK0vMLeeVt8SJ35fB+Dn76HpMWviFa5DRqNbvf34i1dyCuQ+7E0s1P1PiGKM/PJOXnDym4+Ddut07CtmtlhRdjS9tqD6+5fTyJxqmKgWuIuiuuxq4NCQmhoKCgRXa0a6ETE4cTGhpa63GVS6wLSRlZqNRqCoqKKSgqIa+omLzCInJlhWTny8jKLyAzrwC5QomDrTXeri74urvi6+7KbQP64O/tSYCXx02rHNyUyAmCwBe7drM8Yh2zZ0zn5eeexbSOXdqKnEaj4eudX+E3dCxPLF3e/IB/qBI5lVLJ8W+/4K7nV+DVe6DW4wGOJZRx7ctwyvPT8Rr7gKixDaGpUGBoYoYiLx1DE3MGLPwIM7uGQ0LiUjSS2HVw2rXQiYnDqUtMTAzz5jW9eaxQKFAoFNWPq+KHZrz8Jg62NthbW+Foa4ODrQ3O9rYEeHkwvF9vXB3tcXd2xMPJESsLc/EfTI80JXKZ2dksXLaSsxcu8smmdxk1vP4yT1uRU6vVvLxkMadOxrOodxB2jtr1la0SuaKcLGxd3Ji5JkpUrmlcigZBo+bSx6+AINBn3gaMLev3JNUWQRDIO/UbN/ZtJfDh5dh1D8Kue/PLIEnsOjbtWuh0JTExEUdHx+rNy8YIDw9nxYoV9Z4/+vHGDlGmqTGREwSBnXv2snzNWiaMGc0f33+LvZ3uTYsqKipY/NKLXLl0mYXr39da5KpIORfPvg3LuP+NzTh4dtF6XFyKpnJPztAIz9EzsAsYgKGJ7nua5blpJO7ZiDzrBr6T5mHt01PU+PYkdlV7lYXZ6W1tSoegTYSuNeJwahIREaHVIcTSpUt54YUXqh8XFRXVC5bsaCRcT2LxyjdIuJ7ExvA3uH1c4925tPXm0tJSycnO4fm1G0XHyWVcPsdPby9nzJxntRa5mqefV3e8SbcHXsahp+7J7VWiUJqRiIWLD4GzlmFsrn0Obl3bWlPs5EUyCrPSKcnPpbQgF7duvXDv1ouDn2wi8+oF5MWFlJcUEfLkYgKGjOSnDdpvIfyXaROha404nCqioqKqc+KaQ9fUlvZAXW+utLSMdyK3su3z7Txw3718vOmd6ti4htBG5ORyOcePHWXU6DEsWPOuqNzXI+nFlJcW8+NbrzH8gTB6jJyg1bgqkSvPTePC1kU49LoVM3vx1VOqkOemkRj9Fr53h+HUbyRO/XTvH1vTxpaKXUF6MqnnTyLLSEWWmYZHYB8GT32Iv3ZEkXz2BNaOLlg5OGHn5glAlwFD8OozEAtbe8ytbbFxqvxOZq7Zyta5U1r8mTo77XrpKiYOByA6OprBgwdXP64bf9NZqClyarWand9+R8TG9/H18eb7HZ/Rt1fTSzJtRK6wsJAn5j6OmZkp5n69RDWyOZJejKDRYG5lw32vvyvakyvPz+R85Is4DxpHlzvn6lR6ShAEco7/TNIPW3AdcgcWejihrWurNmJXXlpM6rmT5Fy/TM71axibm3PX8yvISrhEUvzfOHj64DfoVty6VVZuHjfvpQY/rz47kp1M02BsXj9rRFXefCZJR6VdCx1oH4cTGxvL9OnT640XhPo5fZ0BQRD4KfY3Ija+R7lCwaqXFzP59gl6qUeXnZ3N43P+h7e3Dxve3ci5nGKtxx5JL0aWmcZPby9n8qJwrUSubqqWqa1TpQc2YKzOn0d2+RgpMZ8SOOs17APFdzvThrpiJwgCsowU0i6cJvPaBcaFLUSWnkLcd9tx9e+B/5DbcPUPBKDnyIkNVk1uy45xnZl2L3TaxuGEhIR0WlGriXXPIH78NZZ3tkSRnZvHc0/M5cHQaZhpWVRSG2/ujwMH6N2nL2+sDudUhvbxS0fSiynOzeK71S/RbehorBydmx1TU+TKMq+Tun873WYuwXngOK3ftyYlKZcxMrfEvsctDHhhG8YW+mmJ2BjHrpUwyNsIQyMjtr/0KGWyfNwD++DVawDqigrcu/fucEn/nZF2L3QSlSgrKvj1WiZbFr1JcUkx8x+bw4Oh07AUUb23OZE7e+Y0VtY2TJ8xg9D77ycuNVfruY+kFyMvkvFd+CJ8B9zC8Aear2dXV+QubF2E+/ApGBqJ/1kKgkDW39+TvG8rXacuwMLFp9VETlVWRP75v8g/9yeF107i8OQiug8by8T5L+Ps203qKNYOkYSunZOek8fOmD/Ytf8gDk7OPPG/hwmdMllrD05bDvy2nxeff46XX3kVf39/0QcPAKaW1gy88z76jJuko8hNxXv8Q6Jt11QoSNi9gaKE0/R8dHV1doM+0aiUKPIzsXDtQvJPH1KakYhTv5H4Tn6SIufKLBNdChNI3BwkoWuHlCuU7D9xim8O/MnfZy8yalB/Nq2LYPTwYTp1t4emvbmvdmxn7ZpwVkes5Y477wLgREqO1nMry0r5e9dHDL3/MfqG3NPs9XX35AyMTfEa+yAet92r9XvWnsAQE2sH+i3YjKmN9uEv2iDPSSHr6I/kxsdg5dWdXo+twW/qMzp5nRJth/TXaifIFQr+PHWeX/4+wf5jJ3FxsGfauNt488lHcHdyFN3FS1vUajVH//6bbR9/QlBw5XuIEbmDidn8sP4VTMwstOp7UFPk5DmpFF6Lx33YPVjoIHJFSefIOrKXbjOW4DfpCdHjG0PQ/BPHJ8vmzLtP4th3BIGzlmPzj6fYkMi1p2Di9oSY6kNQGTkRHh7O/v37a2VFiZ2nLpLQtRGCIJCUkcXh0+c5ePIsf5+9hJujPROHBvPFqiX07tqlevnXUpFryJsrLi7io23beHL+07y9cZNO8x5KyuOnt1/HwNCQO55bjpFJ08vpmiJXnp/Bxa2LcB40Xqf3zomL4fqejXS5OwwDHb3cughqNblnDpB+YCfeEx7Gqd8ogl/dqXNwsYS4LmBQeajYUPSE2HnqIgndTaJCpeJqchonryQQf+kaR89foqCohOCe3Rg1qB8vzZpONx9PvYcXNCRyiQkJPPVEGF26dEGpVGJaY79PW2+u8vChAFMrK8bNXYiJWdM5vzVFTlGYw8Vti3Hsexs+dzyq5Sf5l4zD35Ia8xmBs18n2eBWsq/Vfr1/N5XoOWVXTnB9zyYwMMBz9P3V5Z/EiJzk1dVGTPWhKhpqY6jLPHWRhE7PCIJAXmERCakZXEtN51JSSuV/N1IwNjJiQHd/gnp2I3T8SAYFBjRbtknfS9ZDBw/y/IKneeChWTz3wosYGf3blUpbkVOrVSSeOIz/4BHcsWBZs9fX3ZMrSb6IXbcgfCc/KUrYT181xMDAEI3pLZjcNZhkg4YLYp65Vvtn3ZjwCYJAUeIZbP37Y2Rmgff4WTgPHIeBke6duiSx+xcx1Ydaex5J6ESiVmvILyomu0BGVl4B6bl5ZOTmV9aly84hKSObkjI5ns5OBPh40NPXh9l3h9DH3xc/DzdRhwn6XLJWVFRgbGyMl5cXqyPWMvF23bptadRq3ly+jNwbCfj0CxblyankJSgLc3DqNwqnftq3QzxzzRhBVY5i/0qM/cdg3L1+oG1z46G24JWkXCbphw9Q5GXQb8FmbHz7YOPbR9S8nZW26ALW2vNIQtcAW7/7CaVSRWFpKYUlpciKSygoKiG/qJiC4hIEQcDexho3h8pyTZ7OjvTt1pU7hg/Bz9MNX3e3m1o5uDluJCXx4vPP8r9H5jB5ylT8AwLqXaONN6fRaAhftZLsxCvc+9oGUSKnVpZz+ZPXMHVwpftM7cpgVwmUoChG8cvLYGiEUZfhWo1taj7HpI9J/+NrPEbeh+ej4RiZtf9KNfrkfKIxBqb1b31BWflce+4CpiuS0DXAjYwsnOxscXOwJ7CLF3b/1KVztLXB2d4OJzsbTE1at6imPrw5QRD4Znc0q1et4r7poUy8o+HeDtouWb87dJyMK+eZ8vI6LO0cmry2pshpVBVc+WIlGBoRcN8LTYyqpObSU9CoKd+3EEMrZ0zHLcPAWLd/QARBQJN5BiOPAeRYDqX/8yGYO3roNFdzdPTla1t0AWsKfcwjCV0DvPHEIx2iHl1jVC1ZExMS2PTOO2x4911Gj2m8XFNzCIJATkY67t17a1U4s+6eXNaRvahKZfSau67JenJ199aqatGZ3joPQ/f+GBjq9nPVFGeiPLQeoTAF82nbMHLtxZV86O8o/tDiZlEqy6cwI5WyIhnlxTLs3L3w6RvMmV/3kHI2DqW8FGVZKb3GaN+YSFvaogtYa88jCV0nQxAEdu/6Gjt7e0ImTOTX/b9h2sT+SnPenEaj4ZN31nLm7FnuX7VZlMgJgoCgrsB9+FRcgic2eYJZV+Q0xZko96/AdPxyjDzF3Rg1UV39FeVf72IcMB6TkBUYmP5rw5lrxjqd0GpDU16dslxOmSwPe3dvshIuceHAPopzsyjJy8HR2487nl3Gxd9/4vxvP2JpZ4+5jT3+g0cAYOvsjk/fIMwsrTGxsMTWtXW8Un0gtvoQVHpqAPn5+dX7cs3Now2S0LVDdF22Xr6WwMJVa0hJTmbV6nCAJkWuOTQaDR9vWEPc8WNMfWV9s/FqdT251JhPkeemEfjgK42WP68rcACaonQU+17EyOdWDKx1q0UnqJQYGJtiYGaD2fgVGHk3/J22pthpNGpk6ak4eHWhMCudQ5++R35aEiV5Odh7+PDQ+o8xMDTEwtYeV/9ArJ1cq+vPDZ76EIOn1k+H8wuqXa5JWVbaKrbrCzFdwGQyGVFRUUBl4HBYWFi12DU1jzYYCP+Fkh9aUlRUhJ2dHfGfv9+mS1ddhC45LY3Rk+5l5kOzmP/MAqybKLpZRXPe3IEfv+P7HZ9x55J1WDs13DimJjWFLuPwt6Tt306fJzZg4dpwqaYGRa40B8XeZzDyuw2TofN1iitUpxxFeWg9Zne9haG9drXw9CV2glpNSswnFN+4SGnaVQzR8PA7X2BobMzVv37D0dsPRy9fzG3s9BIzqSwrZevcKRQWFja53IR/f98Ws7+v5dlW264sRf7ZZK3m6mhIHl07Q4zIKZUVfL4rmkD/rowcNpQffonBq04bw8ZoSuTUahUVSiWj7rgbjf+gZg8eoLbI5Z35g9RfP6XX42tEiRyAgYUDJsFzMOo+UbQQCBoVFce3obr0PabDFmBgp31ZfF08O41aRXHSeQqvxlF4NQ7v8bNw6D0MQaXCZdB4uk6Zz4ggPwz/icvrN0GqBNxWSELXAamoqCD6+x/ZsHkLlhYWvPnKEgCtRa4pVBUVbH5zORaWlvSd9axokQOw8g4kcPbrjTafadCTK86k4sRHmI58EePA23WyXX39EOrUE5hP+UBrT66uXc2JnVopp6IoH3NnLxJ2RlCUeAa7wGDcb7sXa9/K6iW+k/7tPmdo1HFPXzsTktC1I5rz5uTl5ViYm3Po76NsivqQRc/MZ9qku2plN2hDY96corycja8voaiggEVr3+WcFts/NUWuOPkixdfP4jn6/gZDNxrz4jQl2dV7chiJLz+lzjyHgYUdRv5jMPIdgYGxfktYadQqCq+cIPfUbxRcOIJDn+F0n7mUrlOfwcjcusm9y44eatJZkP4CHYAbKam88dY7DBoTwvlLlxl72wgO/biH6VMmV4ucth29muLjtyOQl5Wx9K33OVcqTjzl2clc/uTVRsNHGhM5oSwfxU8LMfQMwmTY06KXq6pLP6D4eRGanCsYGBi0WORq2lmel4EgCJSkXOL6nk2YObjRd/5Gus2o9KCNLW31VlBAonWRPLp2zstvhPP517uYOGY0H218h949AjEwMBDtxVXRkDdXVlKChZUVoY/Ow8bOHjNzcyhqvk9EdVtCWTYXP1yC29DJuA+vvw/VmMgBCCo5Rl2GYXLLPAwMtBcNQaOm4lgkqqu/YjZxNUaeA7Ue2+S8gsDJA/FYJO2m8Nqp6vSwQYs/l0StAyMJXTvBtvdgSkvLOHD4MD/8EoOFuQVvv7mC0Hsm8WzY47i5Nn7q2RJvLu3GddYveZ7ZCxYyaNhtwL8Vg5ui5pK1NO0aDr1H4D3hf/Wua9STq5CjTjqEcfeJmN76pHjDBQ2CsgzzKe9jaNvyvckqlL+Ho06Lw3nE3fjf9wKmtv806+5ETWt8EpMxMq4fVaBWybnSBvbcDNr9P1GJiYnMmzdPVL9WgODg4Orgw/aKIAhcS63stH7k+Al6DRtJ+Nsb8fXx5rFZDwAQ1L9fkyInhrre3KXTJ1n5TBjDJ9zBwKGVAaliRE6tlFOanoBjn+F0naJ9KIigVqKIXYbq6i8IGrWozyDIC1Ae/QAMDDAbtbDFIido1Kiu/ILyr3cBMAn6HxYzd1DQ9dF/Ra6F1D2skbj5tHuPTpeCe1FRUcTHx98kC8WhUqv5aO8vnLx8jbhL11BWVHDitxj6m5ny+95v8PfzFTWfrt5cYUE+G15ZyANPPMPYSVNFj6/MX12FoZExPf63ssFrGvLmBI0a5e9rQFGC2V1vYWDY9BLc91pS9f8r5dkkn96Apa0fHteuY2j4b5bGjW5+ouwXBAH1jT+pOPERqJWYDJqNIAgY2v0rnK0ZTCxxc2nXQqdLwb229uIEQUClVmNibEzssZOcupJAQmoGl2+ksOzxWYwJ7k9SeiaDewXyxLS76dXVF0cHewDRIieGKm9Oo1aTnnwD767+rP3sa+wd//VatPXmBI2GhF3rq/NXG6LRE9bsi2gKrmN+99sYmFo2eE1NcauivCSVlDNvYeMyBLduM+vt51WN0UbwBEEDGhWqMzsx7j0V4x53YWDUukUaJNqWdi10uhTcCw8PJyIignnz5jX4eksoKZOTmV9QXbLJyc6Gwb0C+Tr2IL/+HUdmXj5p2bk8PvVO5k+/h7iLVykuK+OWPj14+K7x9O9WWShy9XzxVXUbQqw3V1JUyPurXqOspITl72+rJXLaULUEy4mPpTTtCn2e2NBg/mqj+3JyGUbufTG/N6pBYWlI4KowMDDA0Xsijj53NLlEbkrwhPJCKk58hKBRYzZqIWaTNzW73NaXVyeFmbQt7VroxBbci42NZcaMGVrPr1AoUCgU1Y+rCg7+fe4im6O/p1ReTqm8nMAu3ny07EW+/PV33vv6O+xtrHGys2HCLUEM7hWIu5MD4wYPxN3JAW83Z7z/2VNb/L/7m7WhtZre1ORESg5JVy+z8fWlePv5s3j9xnoFQLXx5qByyeoSNB6HnkMwsa4fTNyYyKku/UjFqS8wn/5ZPZFrSuBKCy6hlGfi4DkGMyvt9+N8ryVVi50gCKiv/oLy6BYMXXtjOmw+gN7L1kvcPGQyWYNl1xujXQudWKoSf7UlPDycFStW1Hve192NedMmYW1hjrWFOQ62lQnpj0+5g7lT65fFGTVI/31E9c2R/b8y6o5J3PPQIzq1TIxL0ZDx57fIrhyn16OrRYmcOvkIyr/fx+z2cFEiV5J/jtRz7+Pe/UHR9tacO6mLC6oLezAduRBjv9tEz3Oz9+oEjQZFaTHlJcWUlxTh7NsNQdBw+c9YKsrLUMrlqBTlDH9gLoU5mTfNLl3QVxcwqF2XLiQkpP03x2mNwn1r165l6VLtKtdWsXTpUl544d9CkEVFRfj4+ODh7Ej3LvW9h/bkAWizbJXJZLz91jrGzJjDA0880+h12nhzOXExpMZU5q+KQVNwHcWBNzAdvRgjjwHVzzclcAAleWdIPf8BHj0ewc7tVlHvCZVenCzjIBXlufhxH0lTPmgXfz+VUoFGo8HU3ILLh2LIT0+mND+HUlk+Q+9/FFf/HkQ9fg8qRTlGJiaYWdkwbfk7WFjbcf3EYUwsLDAxt8DU3BKNRtNslee2Rl9dwGJjY4mIiKgu2STGm4M2ErrWKNy3c+dOwsPDaz3XtWtXli5dWl0Gpi661sLvCPx56BCvLFmMZ0B3xtCyAjWxMUe4/t0mesxeISp/FcDArgtmIasw8vr3b9acyAFUlOfh2XMOtq63iLZXpSwm4/JHlBcn49FzDgB+CTdEn8y2BLWiDLVCjqmtE5lH9pJ//jDxOSkoC3O57eH5DLjjXnKSrlKhVGDn7o1nz/5YOzpjYGDAQ+s+xszKGmMz81riPHlxeL33sbS1v2mfSSz66gIGEBkZyZAhQ3Qq3AntfOkqpnBfXFxcrbEGBgZcv35dtPLfTFprf273rq9Z/cYbvLRkCQEjQpr0ZJrz5uJSNFh6+BM4axl23QY1eE2DYSTlhSj/3ozp8GdFiVxpwQWMTGxx8NKtIrIgaLhxag1mlp74D1mBkR5S40C75WvR9bNkH/2RkrQrlOem4Xbr3XSdugBTG0ec+4/B3NmLW/r7YGnvCMBtDz/V4DzalMTqCOirCxhUrk6q4mjDwsLqOUDN0a6FDsQV7vuv0NCyVRAEft63j9FjxzJ+wkSGDh9BhqZleZ9pF0+T+vdpvENmYWbX8M3XoMipFCh+fQUDKxcw+Xdp1azIyS6TcnYTHj0ewdzaW5StgqChvDgJC1t/fPo9i4m5Sz2Br3lA0RIEQaA8NxXZlRMUXTtJRVkRfZ98B0GjwczJE6eB47Dy6o6pTeU+pmPff/cFrRza/8lre+sCBlQvd6Oiopg3bx4BAQGi7vt2L3T+/v6NqveuXbsaHfdfqid6/tw51qx+k6Sk6/h17Uqv3r2xt7cno5nCmk15cxmXz/H9utfwnSwuPUsQNCj/CAcMMB29pDrerTmRkxddJ/XsRty6PSB6T06tKiP94jaUZRl0HbISU4vGqxLrKnaCoOHknxcZMDQARUEWZ959Etuu/bDrHoStf+Xeo13AAOwCBjQzU9ujyP4TwwaqxGjUSqB9dwELCwtDJpOxc+fOziV0nRV9LVt3fvUl4W+sYvYjc/ggMgprm4ZLlosh89pFflj3Ml3uehzXwY3XhmtwX668EEGlwGziG9Udu7TZk8tPi8XFbwoOntr3ewVQlGWSem4TJmZO+AW9WitbojHEiJ0m9yqqazGoEw8gqCuQB6zF0r0rg5dFY2Tavg8CdKW9dQGrS2hoaOdbukrURmZiTfKNG+za+RULnn+B0aPHcNvIUfWKbjZXJr0pb65Mlo/nhDm43Tqp0WsaLJ6Zdw0DxwDMb/9307w5katQFGBkbI1nz8dEVS+pojT/PNZOA3D1D9VpfEMIimKEslwMHbpSceYrMDLBdNQiDD0HYfVPmT1dRa4jBA63ty5gjc0thvb9jUvUIv70GV54bgF333E7WVlZlJSU4O7hoZfKwgBZCZdIv3gG/8EjGiy31BTqlGOUf78AoSit+rnmRE6lLCL51FpkGX+IFilZxiEUZZk4eo/HLeB+0eMbsk1TkITyz7eRfzkD1flvATAb9xpmo5dg5D0EA0PjJktOSdSm5mEi0OBhYkMeH1R2AasiMTGR2NjY6seRkZGiQ8kkoWvnFBYVcSWh8sew9fPtODk589OvMax9awMODg2XOde2IXVNMi6fY2/4IvJSrzdbbaNea8K8BBS/rcR05IsY2lUeIjQnchpVOSln3sbcxhcHr3Fa2ykIAtmJu8lO+Bp1RYnW4xqipo2avATK98xDUCsxn/Quprc132hbonmqDhOjoqLYtWtXvcPE6Ojo6sd1u4BViV5iYiLTp09n+vTpREVFMWPGDNEendQFrAY3qwtYc/tzpWVl/HboT7776Rdifz9I6D2TWL9yuda5rWKXrannT7Fvw2sMf3AefcdPalLo6jWZVldQHv0IxoF3YjJoFqDdnlzaxa2oFDK69H++VmPq8qw/Gh0jCBpyM/6ivCyTLgNfxsyqZT1N5UUJJGf9gsnAhzB07w/lMgwsmu+R0dIsCX0uXXXpAuYb+ECjhxE3rnwpdQGTaB0EQeBq4nX+Onac2TOmc+b8BSLefY/Jd9zOS08/RY9uAa36/oknDjHqf8/Qc9RE0bXTDIxMMJvwBgYOfoB2Igfg4jcFIxPrapFrSuCqEQQMDI3x8L0LoeQK5SVXMHcbLcpeqDzAyEmMpiT/PI7eIZQ5+FWGomghciCVb+qISELXBpSUlmJoaIiZqSlPLFzM3yfiKC4uYejgYO6543aGDg7mz317dZpbjDeXcPwQXr0GMOp/jaeH1aSmN1dZV241xr3vxci9r9b25SZ9j5GJdXVAsDYCp1FXUJh3DnuX/ji7127gXDVeG8HTaCowNDRBUZKMkakt3W5dg7GZHTfM7bS2Xx90hAOJzoYkdK1EhUpFVr4ME2Mj3Bwd+OzHGI6cvcjlG6mkZucS9fZ67rljIqOGDeXxWQ8yoG8fzJsIytRH85sqBEEgfu+XxH//FZMXr8G9e2WbPjHeXMXfm9HkJ2L4jycHzXtzhZlHyEv5Gd9BS7Xz4KhcTmWmxGJgaIwgaDAwaLhQZ3nWH42KnSAIFGb9RXbCLnwHLcbW9ZZaqWX6CiSWaL9IQtcMFSoVRaVlyBVK5OUKDAwN6ObtybWUNOIvJ1BYUoqsuISunu6Ejh/JB9Hfs/2XA+QUFGJoYMCCGVN4MnQy5mZmDO3biyeefIqegd1xcaqsBffw/aE39fNo1GoOffYeiSf+YuqrG3Dx66bVuJreXMW53agSD2B+z3sYmFUKsDZZDxlXPsW77zNQql1kvFqtJCs5BkMjM1y9x2LYTDXihsROWZZFxpXPUJZl4h74MKYW7lq9t0TnQhK6GlSdy5TI5dXPxRw7yeJN2wCwNDNlYI9uvPfSfM5cvc73h/7GzsoKOxvLyrZ4ZXJGDOhNv27+uDs54Gxvi7GRESVlcu4aXhkkadOzcplXXCLuxLDEuPkzo/i03CZfP55RQnlZCcW5WdyzdA02jq4oyyqbt55Ma9qbE5T//lQ0xZmYjn0VA3M7BGUpPonJNNf5QVmaiat/KIbl6WjrN1YoizA2tcXJdTAIajTq5vtLlKXHYOb6b8pVcd5pTM1d8Oj5GEbGFmjU5Q2O8750kRR/7Zten74Affx136dTluln6aqUlwHiMoE0mgpRz3cGpFPXGiQmJhIQ0Lob/xISrUFCQkKjFUGqKC8vp2vXrmRmNl7Dzt3dnevXr2Nu3rmyPiShq4FMJsPBwYHk5GTs7G7uBrVYqmrn1U3XaY90JFuhY9lbWFhIly5dKCgo0KpST3l5OUqlstHXTU1NO53IgbR0rUVV5V07O7t2/wOvorl0nfZER7IVOpa92laNNjc375RC1hzSGbeEhESnRxI6CQmJTo8kdDUwMzNj+fLlHaK8umRr69GR7O1ItrYl0mGEhIREp0fy6CQkJDo9ktBJSEh0eiSha4TY2FgCAgIwMDBosM9keyM6Oprg4ODqGl7thcTERObNm0dUVFR1F6f2Snv9Dhujo/1G2xRBoh4FBQVCWFiYkJCQIMTFxQn29vZCWFhYW5vVJAUFBQIgFBQUtLUptfD39xfi4uIEQRCEmJgYISQkpI0tapz2+h02REf8jbYlktA1wK5du2o9joiIEIKCgtrIGu1pbzdpTEyMYG9vX+s5QEhISGgji5qnvX2HjdFRf6NthbR0bYDQ0NoVRaqaZEuIo6kGxhItQ/qNikMSOi2IiYlh3rx5bW1Gh0OfDYwlmkb6jTaNJHTNkJiYiKOjY3VHIwmJ9ob0G22e/1RSv7aNd2sSEREhulmuPtDF1vZGazUwlqhNW/1GOxJSZkQTREVFERIS0mH2PgwMDLQu13MziI2NZd68ebUE28DAQKvaaW1Fe/sOm6Oj/UbbCmnp2gjR0dEMHjy4+gdUt4lue6Ohxr9tTXMNjNsb7fE7bIqO9httS/5TS1dtiY2NbTAAs706v3Ub/4aFhbUbj6SqgXFwcDBxcXG1Ghi3J9rzd9gQHe032tZIS1cJCYlOj7R0lZCQ6PRIQichIdHpkYROQkKi0yMJnYSERKdHEjoJCYlOjyR0EhISnR5J6CRuKh2lqKVE50ISOombSnh4eFubIPEfRBI6iTZl7dq1GBgYMGHChHoFACQk9IWUAiZx04iOjmbGjBm1ngsKCgJg165d7TrlSqJjI3l0EjeNmJiYamGr+Zy/v78kchKtiiR0Ei0mPj6e4OBg1q5dC1QmnE+YMKHWNTKZjICAgHpjY2Nj65UFl5DQN5LQSbSYoKCgWsUfQ0JC6nluUVFRhIWF1XpOJpMRHx9fTxQlJPSNJHQSeiEkJIT8/Pzq8JG64tVQ/4iq2mntvVKyRMdHOoyQ0BshISHExsZib29fS7wa89p27tzZoMjFx8fX8wglJFqC5NFJ6I0hQ4YQExNT7/mdO3c2uA/X0F5ebGwsJ06caDUbJf6bSB6dhN4ICgoiJiZGq6Xo4sWLkclk7Ny5k7y8PKCyFHh0dDQFBQWtbarEfwypwrCE3oiOjiYkJKTWXlx0dDT+/v7SUlSiTZGWrhItRiaTkZiY2GA8XEOxcxISNxtJ6CRaRHx8POPHjycxMbGeoDUWOychcbORlq4SEhKdHsmjk5CQ6PRIQichIdHpkYROQkKi0yMJnYSERKdHEjoJCYlOjyR0EhISnR5J6CQkJDo9ktBJSEh0eiShk5CQ6PRIQichIdHp+T+S8EHK8daCoAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 280x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b.change_scales(1)\n",
    "u.change_scales(1)\n",
    "\n",
    "fig=plt.figure(figsize=(2.8, 2.0))\n",
    "plt.contour(y.ravel(), z.ravel(), (Ro*b.allgather_data('g')+z).T, np.arange(-3,3.3,0.1), colors='black', linewidths=0.8)\n",
    "plt.contourf(y.ravel(), z.ravel(), Ro*u.allgather_data('g').T, 15, cmap=cmo.balance)\n",
    "\n",
    "plt.xlabel(\"$y/L$\"); \n",
    "plt.ylabel(\"$z/H$\")\n",
    "plt.colorbar()\n",
    "plt.xlim([-2,2])\n",
    "plt.title(\"Initial condition\")\n",
    "\n",
    "figname = 'plot/figs/2DFSemiu_t0'\n",
    "figname = figname.replace(\".\", \"d\" ); figname = figname+'.pdf'\n",
    "plt.savefig(figname, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a77b15-6a90-4ead-bc61-87de84bf4e9b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}