Skip to content

Instantly share code, notes, and snippets.

@matt-long
Last active December 21, 2020 18:06
Show Gist options
  • Select an option

  • Save matt-long/50433da346da8ac17cde926eec90a87c to your computer and use it in GitHub Desktop.

Select an option

Save matt-long/50433da346da8ac17cde926eec90a87c to your computer and use it in GitHub Desktop.
Plotting the POP model
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import pop_tools # https://github.com/NCAR/pop-tools\n",
"import util"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (nlat: 384, nlon: 320, z_t: 60, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * z_t (z_t) float64 500.0 1.5e+03 2.5e+03 ... 5.125e+05 5.375e+05\n",
" * z_w (z_w) float64 0.0 1e+03 2e+03 3e+03 ... 4.75e+05 5e+05 5.25e+05\n",
" * z_w_bot (z_w_bot) float64 1e+03 2e+03 3e+03 ... 5e+05 5.25e+05 5.5e+05\n",
"Dimensions without coordinates: nlat, nlon\n",
"Data variables:\n",
" TLAT (nlat, nlon) float64 -79.22 -79.22 -79.22 ... 72.2 72.19 72.19\n",
" TLONG (nlat, nlon) float64 320.6 321.7 322.8 ... 318.9 319.4 319.8\n",
" ULAT (nlat, nlon) float64 -78.95 -78.95 -78.95 ... 72.42 72.41 72.41\n",
" ULONG (nlat, nlon) float64 321.1 322.3 323.4 ... 319.2 319.6 320.0\n",
" DXT (nlat, nlon) float64 2.339e+06 2.339e+06 ... 1.473e+06\n",
" DYT (nlat, nlon) float64 5.94e+06 5.94e+06 ... 5.046e+06 5.046e+06\n",
" TAREA (nlat, nlon) float64 1.39e+13 1.39e+13 ... 7.431e+12 7.432e+12\n",
" KMT (nlat, nlon) int32 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0\n",
" REGION_MASK (nlat, nlon) int32 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0\n",
" dz (z_t) float64 1e+03 1e+03 1e+03 ... 2.499e+04 2.5e+04 2.5e+04\n",
"Attributes:\n",
" lateral_dims: [384, 320]\n",
" vertical_dims: 60\n",
" vert_grid_file: /glade/work/mclong/miniconda3/envs/analysis/lib/pyt...\n",
" horiz_grid_fname: /glade/p/cesmdata/cseg/inputdata/ocn/pop/gx1v7/grid...\n",
" topography_fname: /glade/p/cesmdata/cseg/inputdata/ocn/pop/gx1v7/grid...\n",
" region_mask_fname: /glade/p/cesmdata/cseg/inputdata/ocn/pop/gx1v7/grid...\n",
" type: dipole\n",
" region_mask_regions: {'Black Sea': -13, 'Baltic Sea': -12, 'Red Sea': -5...\n",
" title: POP_gx1v7 grid"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ds = pop_tools.get_grid('POP_gx1v7')\n",
"ds"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (nlat: 384, nlon: 321, z_t: 60, z_w: 60, z_w_bot: 60)\n",
"Coordinates:\n",
" * z_t (z_t) float64 500.0 1.5e+03 2.5e+03 ... 5.125e+05 5.375e+05\n",
" * z_w (z_w) float64 0.0 1e+03 2e+03 3e+03 ... 4.75e+05 5e+05 5.25e+05\n",
" * z_w_bot (z_w_bot) float64 1e+03 2e+03 3e+03 ... 5e+05 5.25e+05 5.5e+05\n",
"Dimensions without coordinates: nlat, nlon\n",
"Data variables:\n",
" TLAT (nlat, nlon) float64 -79.22 -79.22 -79.22 ... 80.31 80.31 80.31\n",
" TLONG (nlat, nlon) float64 -220.6 -219.4 -218.3 ... -39.57 -39.86\n",
" ULAT (nlat, nlon) float64 -78.95 -78.95 -78.95 ... 80.03 80.03 80.03\n",
" ULONG (nlat, nlon) float64 140.0 141.1 142.3 ... 320.5 320.3 320.0\n",
" DXT (nlat, nlon) float64 2.339e+06 2.339e+06 ... 5.335e+05\n",
" DYT (nlat, nlon) float64 5.94e+06 5.94e+06 ... 6.263e+06 6.263e+06\n",
" TAREA (nlat, nlon) float64 1.39e+13 1.39e+13 ... 3.342e+12 3.341e+12\n",
" KMT (nlat, nlon) int32 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0\n",
" REGION_MASK (nlat, nlon) int32 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0\n",
" dz (z_t) float64 1e+03 1e+03 1e+03 ... 2.499e+04 2.5e+04 2.5e+04"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dsp = util.pop_add_cyclic(ds)\n",
"dsp"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2b18650dad68>"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEGCAYAAABPdROvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2defzcVNX/3yfTFvBRQQS1LFqWUrYulAIVUDah/fKIiIqCihWRyvpgoZRNdpC9X3EBLQgUZBOEn+hDWxAQ5EGWtnQvhQpVyiIisilLOzm/P3Izk8kkM5k1mZn7fr3SSW5ukvtNZ/LJOefec0VVsVgsFoulXpy0G2CxWCyWzsYKicVisVgawgqJxWKxWBrCConFYrFYGsIKicVisVgaYkDaDWiE9dZbT4cMGZJ2M5rKe+8viCxfY9CIus6lwJp1HGuxJOWVd59GKPb+XH/NYSm2Jp45c+a8qqrrN3KOcXv8l/7ztXyy6y14b5aqjm/kep1CRwvJkCFDmD17dtrNSMQzKzcorA/d6EWeXTk4sl6e9SLLB/JCTddzzblchYFSfuymG71U0/ksvUX/0n1i9znilmznGIMjRSFxVUy96KEFRw17oAktrB0R+Wuj5/jna3ken/XJRHVzg5+J/jF3IR0tJGmw7PmiIDhSLM+ZN7I8xcKBKMWfXLF8RYyIBM/TDJ9jDiUn0fuef2Ew+QpDiILHbbxh94rOTct3Kqx/ffPHCuvTn9mZvHr/C9/Z4uG2t6uZnL/484X107b5fdX6Fy/pY2DM9yaKPA6oW72i4Yple5CjtP73hj1YWL/66c/goJm87wq4JP9bewUrJAkIikc1cmiJCDh4D+WBlD61BxH/S21X4GqgEPuTyAfa+/ILpX9/TkrbnldlFVoiTMHzBsXV1VIBdrW43//032hHfvL5snY98tdNyauDa+5SHik88F2EPE5g2yGv/rlNmTol20F88XCRkgfpL5btFji+uCOPw3vuwJi/s3j+vJZeyyV8/yr/jwevWQ8fyBXXT194QGSdnLg45v98oCPktfTRkAt+f9U7YdAyWUWutA7gBL4FQeFwRAv30y//xbLdCu3IMoqySpO5tiohImsCDwFr4D2Hb1fVM0XkOmA34A1T9duqOk9EBLgc2Bf4jymf23BDmoQVkhBLnt8QCFgG5jccZX1A6UM/V1JHIut49UofDE6EqOQqCk1jD5YoXKLNk3xEeVTdgebh6z8Ggsc5gewJrpQ+cCsx/28bm3N59VdpDhr424Mi4guMWyI8RnTU4X3zIHVjhCjYrqzzn/wgIN7V5OP//b44OoH/w1Xms+RBr7kSMQl/L/zvaQ4tCL8jLnktCohLrtCuHG5RVMXlF8t2K7FUskKTLJL3gD1V9W0RGQg8LCIzzL4TVfX2UP0+YKhZdgKuNJ+ZwApJAhyJFg9fOKJEIygWwQd/ad3SB1FYYLw6Tmg7+uGVk9bYMfkIl0X4geGG3iLzRjj8evmg/xwN7AcE8rgFa8YF8iIlloqDeA90yePgkEPNMQ7ggLgRb/2eWORxcNUpiEfBkjGWR97sc01dwNSXgPAE1iNEB+A9d0Dh3JeOvLXKXW09k+d/rbzQ3GMn4kEYFpmwBeGzSnMlYuKYvz8nxcdrQVwK9yqqzAmUmU8pWin+fb366c+YNnttSNvdpYHvb0Pn8XJTvW02B5ql0on3B643xz0qIuuIyGBVzYTfuauF5E8rNi986XNoyZtPuMzfDn46Uloetj6SCkiceJRbJk5kPa9uWFDihWPcBiNj9zXKrBfnl5WF39DyhbdNLdufV8UVP57kfa5St8yaWYXnDswj5FTJi5iHkAvi4KDkUN5X8MXE3x/ltvJFJCwgq4ybJigeYeHw63v1whZKubAfNfebpfdHhffc0p/atTtcW3ZcM/HF7LgnDy6UrdIcaziryeP9zbmgeIReGIKWQtCC8Ajc35BAlBD47bkBSyNY5oj3/5XDxVUpuL3CMZQsEWe9R7CeiAR7A01T1Wn+hojkgDnA5sDPVPUxETkSOF9EzgDuA05W1feADYGgr3elKbNC0ir+tGLzmo8JikicgIQtkDjxiBKO0jpOyT5vv1Oyz6dEFAJiIk7oARbcN4D4eqG6tdK3SYw1HXoQqRv9kFJXC9u+KLm45NXFRQuis0pd8mjBgnFR3tXV5iEP7+MUrAJfTHzPtYvD+5rDVfOJw7vuQFZpLmChCKs0xyp3QEEo3IIFE7RcPGsoaJ1AUUDcMksoHP+ItiC/8djhZWWr3Ry3fvrnkfXjmPD4YSXupbBAvbl6TQAGissAJ18iaFeM/lVh3RecqBcrKP27HIovZG4hXqKFlzZXcyVusSiLJChAYTHx9jsFy8SzdrLhRlSi3b0xvKqqY2LPpZoHRonIOsCdIrItcArwMjAImAacBJxDtE83Mxl3u0pI7l8xrCTYndQaqSYi1QQkbHmExSNOOIKiURCMkFhIrrwssFFcLykv/c6JU6w3883SB03f+keg778PTtQbZf0/3rIjg+4A1/s/2HfYrqCKum6xjhEgzedB3YLwzHjxSd7TVSaw77JKXVahvOW6vKu5gmXxrg7gfc0Z4RjAu+oJyNv5NT0hUc8KWa058sZS8IXCf1itcn0rRQplQcsjKB6NBsJ9/O/ilx85qnxf6Hlx285Xllzff5hXwkV4PyAiV4+5jkOfOJScKA7KAMecQ4vf/f7RN9b8d8RR6DkWJSiFde/D/4364uKLCXiuru9u8aemtasearBIEqGqr4vIH4HxqnqpKX5PRK4FJpvtlcDGgcM2Al5sakMaoKOF5I33FjPrua09N0eC3h5BESmWaWxZiSsrQkSiBCQsHsXyUoujTDjEiERYNKLEIvCAL4hE8KEfFIWcg4QskL7BR5c82GXNNbyVQV5QllWrys8ZPm/U/riyKpQc4bfL/3RdUGW/Ufugq1ebbl55T2hWr+au5/7MO7qKt/RdXsuv5nV3LZ7Pf4DX8x/gjdUf4B13EP/JD+StVWuyWh1WuzlcfCtDWO06ZWLRiDhEBbQjy8KB6UCdcP3w9n5/OpYBjmv2DWJAhe/+++4AHJTVgbLw270vMlePuS72PI3idzs+f/HnjQUiRZeXuGWuLtQtcXMFxeSap3dNLVaiwKomxEhEZH3vVPq6iKwFfA64yI97mF5aXwQWmUPuAo4RkVvwguxvZCU+Ah0uJJUIWiNR+D/OrTf2Buste36DQlA9ypUVZYVEWSBh6yNoeYTdVAXrwjzow9s4Ulk0yrYDdYPnCNTp2/BY84fkip/hh//AgZHnLJTFnLu0TlDYSvdryb744wv1/Dp+mQPqOCCw79e3xR3o4K7hsOoDDu+s5/D6li7uB/LIQBdngIvklNwAFxFFwPv0rVFHzeWKn8U/PSC25s093MyyOqHtcB3/HBIhHHFiUthPsd5q18ER9cQQJ7ZH1mq3tKtusN61O1zLhMcPizwuKcc9eXDhnP2jbqlav9o4lv6l+xTaGhUz8e9BrWJy93PbJq5bCUVrcW1VYjAw3cRJHODXqvp7EbnfiIwA84AjTP278br+Lsfr/ntoMxrRLDpeSMJvdtWIF5ZK1zDHhlxZwfXIIHvAEghbIBWtj9ADOomASCXLJXycLyLhB3hYNJKcM+o4KDzwy+oG90OkSEBRKBBBc4IKaM7x1h28zwHC6jUd8oOE/CDIrymsXgMG/csh/x8Hd5CiAxR3oLJ6kCcoiCKOtyDg5FxzeW87Smy8PzG4XVoeLHMoF6KgoPhxnKDo5NU7zo0RF/+B6iLFbYrlQMHC2u9Px5YcO6Dwf+x3DqCEoPV1yGPfLazfsNPVVOPYud8o+a/3e4r5gf6TF3ylIADnj7ij6vmKbXJiPQwuUreYNAWl4kDexKdRXQBsF1G+Z0x9BY5u/MqtoeOFJAmOaGR3xhzKkuc3ZOuNX2DoRvHuxpdf2CAyHhLlyop1Y0VYHZUskBLxSCIccaIhUfvCYuBEC0VE/RIBCNYLC0Qu+piCQPj7HFC//WIEIiAe+CLigDvAlOfANZ/5QYI7wBsj5w7wrpF737je3/fquzlwBzmQUyNI6rnlHSWfU+/2G2EBRXIaEBZAFMcpCkktglOL2EQJT5CcuGXWi/85a7cf8d8P/U+ZW251QZnhzl1+Rpgbdrqagx+diBOyuIJB/Fyg3T7TxkwH/HiRcUOZa0+adxA5UQZUeDkLc/GSPnM9Md27k/Xkmv7MzkwY+kjFc//u2RGxGR5qRYkfxNvLdLWQhK2PuG3freWLyfMvDC6JidQqIlFxkDI3lm8VBB7SkZZH6GEucQ/9OOGoRzTirIo4i8Jcu1DfFw/fogicwxONYt0S0RDv+up4QoF4IuNZIEXx8MTH2+cOBNc4B1S8BRckj/erV8jlQdS7tici5tqOmvYbETGfbs486APi4jpFi8UXGDFxijiLBkrFpJrA5ELiEBQUEfV6QwXE5aG9LgFg9/sms8f9JzAo538PzXFa6mrb70/HFs5XWFAw3YFv2/nKYm8y9f9zAVxyvkUU+A39ZPSNXlfngJgE/w5v7E/yHoJ+jMQnLB75kKWSJC7622dHVegoXw/SMQNR20lXCUmc26pS3aBpvmLl4Ng3l7CIFMud+HhIhBursB7lvqpmfRSG2Tuln3FWR9A9FbxehLWhhXNSeq4qsQofHZDzAuRBK6NwjoBwOL5oBCwOXwQKIkPB8vDLfTHxz+2LSUFEHAqWjoRfGx2Q1d6tUNc/VhHHu4Z/Ta/tirjincRvo6MFUVEjLH7nCF9ACn+qE3CV4Z3bd2mFYzPhMg3Fa4rxlKKY+K4xR5Rd/zAFgAECajoO+PsA7t29n76HjiucK2zheA96CvsBbtzpKuI49IlDCT9D/S7Ex879BuF3dUeSd1yYsvUMLl7SB5gu1gFRCbrwgDJBu+GZsRwy9NFE12kUL9huhSRMVwlJEuKCklEiFH6TCactiRsUGAyql4lIkLAbq5qIlASpQyIS5baqZH2E3VX1CEjQ8sB7+BcfvqVWjV+vKCR+/WKZFqwZXzSC5UZEJGC91PB7FjUdwRTTA1VQIwieNagoIK4pN6ImoqBeknTfzVWwytxiO0TUdDRzUGNxiHjXcx1/vyDqiU/eNf8NeCIQLgMvnQz4g/ckcNtDgkAxZhJmxmcvB/AEJfQVD8dfqlFtEGXYonBxiBpFH8UFS/YlF3GOSoMT8+owULz+aDc8M5acaEnizVagdE5qnHbSMiFJMylZSQqHGOEoqV/m8ip+UYLdfCFgmYTdXcYaiXRnUSoovihEdd2d+a/qQU6Avo3+p/TBXxInCQhIVM+uCEFIKh6xMY7CeYruqpI3/bB4BMqDrq2g26pUQKLLCi4tKBUWNZaJFh/4ha+FX0/Nk178hzme5aECrrEQzN8kStGy8t1YQYH0rRT1RERFiv8d6iLiiYWa+1QIkhcsEiFn4jAuYiwV///Ca7jf6/SJvguqfj/C+IISZv+HjymsH/B/R0fGUZJQ0p2XoiUx0AEHlynzDyyJVQYHMJ617W9LzhGFq0JOyt1bPr477dblOxSEaJA0nlwxri2WUlppkWQqKVnkOJEK34coa6SaiBQIxUSSiMjMN66p6e+ZsfLHhfW+IZPACcRc6rU+QuUaDpj7D/+A5eBfs8zyCAkIAesi6MIC+L/bTuDTB19WbpEExUUiRKQgaKGbYx76foxETJkaMSl4dMT/+4tuKFVBcsY6yZe2sxgjkYJ4lVkpeV9EvMC+JxLGEWNEQ8RYG/4tLVzAU7xwTMX/NqrZHvW/pzPvv88t/Lk7zjy1IEyPjqtNZPzeXv64lHq5YvSvOGLOIQwUt2BVFMQyZJmEX+68wYpBV5ZDTvJl1klk+2uwepqBtUiiaZmQdFpSMl8GwtZIXJkvHGVdfBOIiE94lHn9jc9FCwhALhQIhxIrpURAAm2b9cSZ7L3zeaXigS8C/pPYlOWkXEAcCuIRLvfXvfNR/AzsD4sIgfNoeB3KxER84dCi4eELjO/mEjCWhyCu/+APWBPGghAF9Q8098uzLsR7uPsuSqFoufliJn7oyLdupNT1FfwbXKfgIvPaYtxOrnd8XotWzui7f8Dcfc9jhxmnlng8P33PyQD8eZ8Ly78nEfzuMz9JVC8Kf2Q8UBjMWHjIqoMjefIq5cF3f7ChsS6CVoZvlQStk8KgxICbKxwnCdOqXF1KMbmnpUhL74iI5ERkHvAKcK+q+g7M80VkgYj0i4gZVh2blCx8zokiMltEZr+RcMpLn0oJGqHYY7XWmxIZKwn3VgmZP00TkcL1IkQkFKcI14sq98ZnOJ6IVEHNOSKtkICIlB5TrOOLyNhvXMajN55QKgohkYgUjKQvhr7BocV1ytalYJgU9gXKSl6ig/t80yQQgyn42dSzcMB8FuoIquIJRsEvZw4N1Pcyxvh1/fJiM7afcVph3TVxFh9fUFpFeCDjxNkTyuoUElyqcPHI27h45G1cOOL2msaUZBHfiqu29BItDba3IimZyZ45DWDo8LXK9vtmbqkrq/63E999Fe6pBdEDDsviIoZZ/76+6rX6Pn5kmWUw46XqPusZz3rpefq2OKmyK6tgBQTdcMKs2Wexz47nlLq2oNAVt3guIoPoZa6ukNVQ4uIKxksCPa28dhWPKR4fsmoIrVPBzPWtgvA6vpVBwfVl+kchfq6vwi0SPyzvBeKNpeL7xoKuroL/DA34vUDdogXjmzji+LOcayA+U2xmodmeZ6wkKA9+V1zjIlPPeZY3o93DwfhW4cdDglmEwz2sfCbNO6hk5Pu5w+8srJ+z6AuF8/npUwodDiLcW3FxknagCO8nyG3Wa7Sl11bWk5JFWSBRc4NEjV4PB9iDn7PeTZb0rm/9IyoHbBIy4+mL6NvqlOIDTKRkYGDRUqFEMIIiEhsLCZ4z6G4iWIeSWEhJOf4zVosCZcofvfEEc55A7EPC1wjsD1LttgWFIrju3xYtPrx9MUG1JJdg0R/m1wuIhRGHgqB4G+Z2FR+oagIzRUEp+K+KppIrxf1+Ly9AtejWcgrng3n7FeMk7eQbjx2OQ1Ew8n7Pr8BD3xeDoLCExSRIo+IQdKG1EoVCZmhLkVb22mpbUrJavoDVxpo4Jevx1kh4PSwi97x/U6L2jF9vIoUuvoGHfRJLJBJfDALbZeVB68QXhkLchFgRue+BUwDYY5+Lyq0N4MH/nVJfmw2PX3d8YX37w/tLhSVA3d1/JaIsIC6Fz2B8JXx8UExcDbTDFxfv0ztf4MxGRNW3cNxiLy6vhEJdT46KSue7tZZ88czkf3AbcE23OcekeXcCvbbCwlKJM7a9q8Qq8fNxXbp0XEmcpDhRQLrYYHs5rbRIOiopWbUUCnEzHjbC+PUmNuU8ccxcfD7jR55e6q4CZs0+i3E7nF3i9hLVYiDd4D94fQHxeeCek1ra7jiiRCVM8OFffgK/UpKLFauqubigpWIScZ5CYD4oJhK4sC8mfplSUk+NlSOF7mbeA3vZl85I0Oj2US3YXVI3xt0V5Ixt7yorm7zVLPqX7lNSlmTSq0q5uhpFtTjJmaVIK3ttpZ6UrNqXN846CYpGeBCii8tag//K6pe9ybPKBh/WSEnKE1VmvHxFzecIMmPJD0u2Z86Pdn/MeqLy2+3ndvPOc98fT22oPXVT5uZKdliS33hBbIr+rFKrRItxkcK+wAO+2AsscILAWYL/+hYOCiu+nY4AN5twskgw3XUD4uIaV1PYKqnk3opi0lb3cPlTn4tuB0Ia0YqsTLKVJXpuZDtUEJAEX5D/GvzX8sKIfELj1jqkpEvurLenlx8WOm7G368sq5MWf3gwJQExzJk2qbA++sh+AMSFJ6/wykce1192TEFsQoH1ZlPq+gqJSUnsxYhJXvjrYY25/bJIUDhKu/kWg+a1pC2KvU5CF1k78ILtPfnYrEjP3JE46yTuBTaJ+yrKGimbjCqGvvWP8BtmTiCFuUJmvFB/3/5uZO6Vk8rK5l8+iRGT+ouWRTAspLDogtJjtjqjXHj8uiWWBxTMlaBVUlwJuqoixCTo+zLbf/3uidX/yA4jbJX4X+M4l5dfv56A+KSt7gHgp0/tGfk7buegRBtsj6ZnhCRMEi+Ub6FE9eCKJGbqWwnNLDh+3cO97L/hRIqWmljQXy4wcSw9J77uFudFiwx4QrF88vElZZtePrUoMuH/toDPzE+rMuSqS1hxePeIyc1jpwFe763KVkk+sN74g/6YLe/nimV7mGul1wU4n9TP2kP0rJDUwip1WWeD56N3xlkfrjLrnRtqu5Dr1t9by9I4Cs+c6gnO5hfHi8uzx3nCsslPLgNCVknwGdPlz5twpuBDn6jeP6bRh/BRwx7gF8t2S1S3FV2C7cj2aLpeSJIMRgz22ApaH3mUj26wMvKYAZ9YDhwIxAfZx31wQtEaCU5pGzV9rSU1nv5BubWyfEp1a+e5Y08oKxtypT9MqujuWjGxe6yRSkRlBz5q7jcLcZJmjfb+3rAHE4tJK3Btr60yul5IgtQb+Pv3S58qjBlZa/BzlSsnncgn2FMrQ0F2S3J8i8THyy9WXu9TV19SiKt0Y9C9l1CwFkkE9o4Ygj22goH2N1/8ZLIThOMjjhR7aoXSjgSxItK5lHhNgv+t6i/W2rxi9K/4yegbW/YWHx7T0eqHvCKs0lyipZfoKYuk6YSsD3EE1GXWv4upUWa+eS3j1/5OsZKrzHztKi+vlqX7iOh6LAorvlubJbLpzd44HhEKsy5628WBL09/+fQGG9tZ/PQpbwiaN69JOqiWi5fFCklFosaVrH5580BqlMB4S3WZ9W70QKuouUasJdKZbDZ1aiCvWGBH3Gh30x04Sc+tTX51QXEcSqhLs5nF2Oz0PobdcQ4CPJWxUe+Ncv7iz+OISw4tfAIMzISBJ00ZkFhh4r9NgFuAdYG5wCGq+r7Jkn49sD3wT+Brqrqi4YY0iZ4TkmCcJPxe4ZTUC7q6HPMZ+gKpy73ubU1uoaWdDDu32Dtr2enxAfbNLp0amd8rKocXUDo4xewPxkrEobBebcR7Metv4LpSyMbFVneeDcDSA7KViyvIT0YnS2CahKOGPVBYv+bpXZt23iQoTbNI4ib+Ox7oV9VbROTnwGF4k/wdBvxLVTcXkYOAi4CvNaMhzaAnhKQZo2st3Y8/nsTvxbX5Rf1F4Yh6dvjD28Nl4W7ACfFTs8RaIjF1fLb57Vks3v+s2i+cMfykjRcv6Uu5JdE0Iw5TYeK/PYGvm/LpwFl4QrK/WQe4HfipiIg5T+p0pZC0WjhcXAZ8Yjn3ZiNrg6VOtjy7P1Ighl7QX3N24YqE4yah8w65/qLI8vBxZaIS2O1bLN0iJkGKsyZ6CRuvWLYHOVxy0v4x5krzJq0yCW3nAJsDPwP+AryuqqtNleDkfoWJ/1R1tYi8AXwUeLUpjWmQrhSSZhGVJiWvLmsMfjaF1liSss1JRSEomUSL0rIoN9Wy0yd5QuITl7crqrzWMhVwol96wu6sattBht/lubgWfuHsyP2dSqumz60FBVYlz7W1nojMDmxPMxPzeecKTfwHbBVzSYj/FmaCrhGSeq2QaunjfZqVOt7SWrad0l+XJREVH4mMf1RIIV9GpIgE3GHhRF9xx0W50MKnNcISzEU14ndnsGC/cxI01JIcqWU+kldVdUy1SoGJ/8YC64jIAGOVBCf38yf+WykiA4C1gddqbX2rsP3YaiAnDjlxCinkLdli2xMr58yKREtFZIvz+wvlUXUTnztcphH1wmqkFCaxqkS4TnDbd7u0a7rdVjJl6xlM3moWk7ealXZTCijeyPYkSyVEZH1jiRCY+G8p8ADwFVNtAvBbs36X2cbsvz8r8RHoIoukFmqZNiRxwkZLqgyfXKclckZCSySOuIB7sCyiTjjeEX9+EsdJohj5+9OZ//l0puRtBhcv6TNdgLNDk2ZIjJv4bwlwi4icBzwJ/NLU/yVwg4gsx7NEDmpGI5pFTwpJJawLq3fZ4rx4MYpzcyXuDlwHtcZJwiLj7xv1v6fjiDJ33/Oa07AW4k+5mxOXgVLb1LrtmLddVZoySr/CxH/PAjtGlL+Ln9wvg1ghqRGnBm/g+LW/AyZp48x/XZ3omL4hkwpPgxnPTa29gT3KwkuLlsU2J8W7uMJseXY/T505yRtP0q53iEgrhtriJNVMEbyHnuNo1ZlCs8JZi/bPvK/dC7ZnyT7KBlZIiJ8Z0UX58OC/8c5LmwTKXJP5F/YZeFBJmpTgHCOz/n19ze3oG1LqZunbdDIznr00prYljsUXefdxm1OSCUpwUGLTSfDAL3dfVbY0Ik+RoFfX6Lt/kHmrxJ/fJGr+kiTztbceO2d7FC0Tkk5PARBMIV81428E4z44gVlvT2f8R75bXwOSZhG2lFGwSKrdQqV8jveEXXgTubCSHFeHpVErwdkMLY3hBdut+ztMKy2SrkoBUA0JznYYwfh1DwdHmPnqtMj9cfQNnVI8pwgzll3YSDMtXUC1gHunCcdpC75kpuFNuyXJsGnky2nZHVGPuBQAt5vy6cAXzfr+Zhuzfy+RdLpM5ZBYd5fP3k71uNf4Dx/qrYgUXF6FudrjiEs5b8r7tjql6nUtLSbiGV21a2/EceWdNzvkSdommtQ7qqn4I9uTLL1ES6VVRHIiMg94BbiXGlIAAH4KgJbQSBqVsIhEzZBYmBnRLfp0pR53VYSw9G19au3n6RG2ndLceEeiF/skD406HyyVxoxEb4ePr+uyHUFasQoXJ9HSS7Q02N6KFAAiMhGYCLD+BgNjr52maV8QERGv15ZGO777Njy22ENr5Y9L9vlB9r4tQplh7biWWLY9sbGeV83qultv/KSuOh3OlPkHZiQ9fDJUYZXbWyKRhLb02mpmCgCTq2YawNDha3Xs+1bf4KMLXYN9ZqyIeZu24lGVsgGJrbplaT7cqwTmK+XfyiJT5pda9nl1yEX01soSnmvLCkmYVvbaWh9YZUTETwFwEcUUALcQnQLgz6SQAiDqq+FPs/vhDf5WKHNf3gIYWVuvqpAQ9H38yDIRiWPG0xclv06PMvyEUhFpmkWQBGtZ9BxZjN2kTSstkq5JAfDeS5sCXq6t2AGJlYTFKY9zWBpjxCRPPKMRgpsAACAASURBVBQiH9KLL5yUeBxJLPUIQL2iUeNYkiTWR9bHkrg4OKmPC6kN2/03mpYJSaenAPBTpVTKtaWucm/+ppKycR80edX841wXctWtj76Nj2PG85fX19geY+T3q8dCtjm5+SPVm5n+xAKXjry1zL2VfaxrKwo7sr0OZr04H+cTT0fve3t6sdsvFFxYJT22wuLkd+216VGaQtYf9mkL0vYzTmNO3/npNcAwef7XOrJvUzPmbO82ekZI6u3um1flvzZIPpHVuA9OKPbaqkS4TjiOsvmJhfUZyy9JfP1uYPQR/eUTUQGFFEdZ+R03Lc5SZXR7krxaoSrhQYn+dlaC8ZPnVx5r7KoTmSYlbbxeWzbXVphOfCHoHURsj606yMizMrM4KDvOTG8s0qR5rQ9/tmr0uR2QGE3PWCS14qKss8Hz6Vzc5tmKRNyI3FhhovJn1XOtZrmf6rAuWo3/bjJ2lpcl4dFxF7Tt2sc9eXBkf5NOSutiXVvl2CdWk5n19nRmvnktM9+4JnEXX0spc39ePtmU/9udf3nEPr+KwqJL4vdnkTTSpISv+el7Tm75NbsFv9eWtUhKsU+6CPJNskbK5iAxr2IzXr4i2QlUmfHMxQ23o2OJ+C2OPC6mS696c5I0O0VK8fzZeTA0M01KWjGTfIbuZ600Y6rdbqO3/toqtGJa3bgJrWa88JPY3luWhCgsnDqJhZdNAqkwZ3tneExSZZd7T6peqQYmPH4Y3539bb47+9tNPW/aqAqr1Um09BI2RtIGZr46LTLrbzi/lqUywd+mmFjI8OPNeJEMaHDa3XrrpRmxiYMfnYgjaialKj3fEXMOKayv0QXP115zWyXBCkkIF2X9DV5o+nln/OPnTT9nz6JerGTE8S2c2bAZRE2p22M4XWYO2pHt0Vgh6QLGjzoDFQEHZs05O+3mNIUnr5jEdkfHC4WfIqWjyLiw7PqHKTz8udpicgc+cqRniXTa/0UDWCEpxwpJgLwqn9jwxeoVs0gXuAwS0y2/4yRdg9vQXdhVIed4g/92v28yUAzCP7DnZa29eJNoV94ufxyJpZReevxUxIVMiUjf5ifSN3SKN9VuDYwbc1ZrGmTJBuGHWFkPrco9uiyN4yKJlkqIyMYi8oCILBWRxSJynCk/S0ReEJF5Ztk3cMwpIrJcRJaJyLgW/5k1YS2SLmDmvHMYt/2ZhV5f43Y4u+DqArjn0TNSbF19bHdU61xXnRoUbzWOKC5SiGs0e4Bgpww4rIQqrG7OxFargRNUda6IfAiYIyL3mn39qnppsLKIbI2XEX0bYAPgDyKyhZk8MHWsRYI3bmSDDFkjDWG7EFuayN5/7KwBnu2gGQMSVfUlVZ1r1t8CllKcdjyK/YFbVPU9VX0OWE5EFvW0sBZJ1hGhb9jJzFh2YcVqs+acHe3WEmHvXc4rrKsDaiKj9/3Rzv1uKSXYyyrLFkQeh1wKc5nUGCNZT0RmB7anmRleSxCRIXhTbjwG7AIcIyLfAmbjWS3/whOZRwOHraSy8LSVnrdI8gobb/hS2s0oo5Dxt04L457HPXeWVjh+r91/yF57tC/PUlJGH9Hmbr3ZfV6mRjC2Mu7B76fYkuyhKokW4FVVHRNYokTkg8BvgO+r6pvAlcBmwCjgJcDv7RD1Q87MN9daJF3ErNlnlWzf8+gZ7L1zhdnwMugGG/09LzZi4xjp4vfa8q2ScQ9+n1m7/ajh83rp4Wt36+eRzLz1Nitpo4gMxBORG1X1DgBV/Xtg/1XA783mSmDjwOEbAZnxx2fl/8YSQS/l2dp+Yr8nIlkgM+953Uknd59VbU6MREQEb3rxpao6NVA+OFDtAGCRWb8LOEhE1hCRTYChwONN/eMawFokHUKSOEm97LXnBdx3/ymJ63/mgEuL1oL5fPg3k+u69pjDpnoxm3Cv1s591lgq4IhbsHLyKh04kFHIN6fX1i7AIcBCEZlnyk4FDhaRUXivMyuA7wGo6mIR+TWwBK/H19HVemyJyFrAJ1V1WTMaXAkrJB1E31anMGNpk2Ma5oe8594Xcv+9rU8nvuO3pxZnPLT2cFfiqnDDTlcnTtiYVuC8XpoxNkdVHyY67nF3hWPOBxLNkSwi+wGXAoOATYw4naOqX6ijuVVp2U+52wbcpE4LX90qBeTDfOaL5dP+qsAuB3bGCOim0YPur76Hjisru23nK2s+j+/2yQV6heVxOHf4nZy17W/rb2Ab6KD5SM7C6x78OoCqzgOGtOpirbRIumrAjaUCGQzat51m5dHKwAyKUSTtCuwi3LjTVaVlCedfP33hAeTEzXbgVivP75IhVqvqG9Km32bL/s+6bcBNVujbOvnYj33GnlPTuffcu3IM5rNfKLdG6kEUnrj2+KacK6vUrSkZFJEg//3Q/7Dfn44tKbv1015ma1fLRSSO4Bu7q1Lm2spn2O/ZjBQpbWCRiHwdyInIUBH5CfBIVEURebPK8paIPF3pYm353woNuAFvwM0CEblGRD5iyjYEgtMSRg64EZGJIjJbRGa/8Vr3GyslAfYa3i722bE2EWkGO381uXtrh+9MrV6pmwgLS/quj5qJmmvE5+ax07h5bNkwifhzSeWBj7kEFkwaqAm2J1lS5lg87857wM3Am0DcgKC/qOqHKywfAv5d6WIt/2ubPeBGVaf5A3zWXjfXolZnm/HbnBa7b9wOZ7dERKpZI7U8F5PUzfCg6rZRfQrd2hM0Vjsm/CYd5+vf/+FjOOD/jq56vW5ENdmSbhv1P6p6mqruYJ6Xp6nquzHVv5zglBXrtLTXVjcNuMkKlQLj48acZeMVlpbjzT/iPSm//MhR/GbnK1JuUXvJckZlEfkdFbqCRPXaUtVnzbEXqWrJnMt+mV8njpYJSaUBN6rq5yQJD7i5SUSm4gXbMzXgplfZ7b8v9nqMVfnt+CPRd/6aZ2A+cusJbWidJU18MTnwkSMZ4Hhu5rB7K5hNuBvwrI3sCglel1+ALwGfAH5ltg/GG5dSib2Bk0JlfRFlZbTSteUPuNkz1NX3YhFZKCILgD2ASeANuAH8ATczSTDgphkM2Sh7ebbCJBmIOG77M9vQknjCv61PHxwdL3n8uoRB9u559nQ1lUTi6jHXcfWY65g2ZnobW9R6stz9V1UfVNUHge1U9Wuq+juzfB3YNeoYETlSRBYCw0zs2l+eAxYkuW4re209rKqiqiNUdZRZ7lbVQ1R1uCn/QsA6QVXPV9XNVHWYqs5oVduCrFg5uHqlDFBpIOK47ZonInvsc1HD58hwhxtLi/nGY4en3YSW0wkxEmB9EdnU3zBpVdaPqXsTsB+eV2i/wLK9qn4zycXsyPYOZ/yoM5oWF6llYGLLyVBTLMnwXV2HPPZdAG7Y6eo0m9MSFMFNv0dWEiYBfxQRP7YxBJNuJYyqvgG8ISI/AF5W1fdEZHdghIhcr6qvV7uYFZIOYsaSHxbWx4/4gScgWXr4J+SJa4/vve6/PULQ1XXoE4dy7Q7Xptia1pC+sVEdVZ0pIkOBLU3RU6r6XpXDfgOMEZHN8eLbd+FZK/tWPAorJB3J+OGndaSA1EK245mWamR5UqyGyH6wHQAzMVaQkSKCql5f4TBXVVeLyJeAH6nqT0TkySTX6wgbzdKZjP3GZYz9Rn05uILPoW59JvUCSRM3dhSacEmXHQLLZ/Byb1VL2LhKRA4GvkVxWMbAJBezQgI8/0JnBNzbxR7jvID7bvvWPx9KMOA+9pu1u7E64KXPUoVcDW8A+Q56FNUwQ2KKbdRjA8vheJlFBlU57FDg08D5qvqcCdD/qsoxgHVtkbMPrFKafT+63AVnaR7nDr8TgHMWtSTTeVNQwHU78jv9H7yxebGo6hLgfwLbzwGJJkHqnNeAFvPiCxuk3YSuZqdDilbJDofaQLullAtH3F5YP2Pbu1JsSRUUz1xOsqSIiPxORO4yy++BZXjB80rHDBWR20VkiYg86y9JrtfzFglArlP7mqoyc/65QPoDEisR/E3t+O2pya0ek1Ldurm6m2qD95KmoW8XGRgjkoTgNB2rgb+q6soqx1wLnAn04w0WP5SEv1YrJAFefmEDPrFh9tN7zVyYaJK0zLHTt6bWZgNbAekKJs6ewLQxN3DU3PKxbS4Ol468taz8tG1+X7J96dIMzXPXGUKyb1zerArHrKWq94mIqOpfgbNE5E944lIR69oy2BvRYhoVhc748VpC+ONKJs6e0NB5Jm81qxnNaQLJAu1pB9vx8maF6atyzLsi4gDPiMgxInIA8LEkF7MWiSWS3cdf1NLpfROR0dkCLT1Ohl9qRORI4ChgU5PP0OdDwP9VOfz7wAfwAu7n4rm3Er0BWCGxlNGql6maz2tFpKtwtZhF2lWhf9Qt6TaoHhQ02722bgJmABcAJwfK31LV1yodqKpPAHieLT20lotaj46lMhl4+xKFeT+dlHYzLA3QjGD5cVv+oQktaQaScKlwBpGNReQBEVkqIotF5DhTvq6I3Csiz5jPj5hyEZEfi8hyk5l3dMypVVVXAEcDbwUWRGTdKm36tIgswZsWHREZKSKJJpuxFoklmgbHf/jzkzSEwpNXWAHpFrombUpz/ozVwAmqOldEPgTMEZF7gW8D96nqhSJyMp5VcRJefGOoWXbCm2l2p4jz3gR8HphDuXNYgU0jjvH5ETAO001YVeeLyGeT/DHWIgnxjxfLponPPM1MI581tju6n+2O7k+7GZYWMGneQWk3oT6akCJFVV9S1blm/S08K2BDYH/An8BlOvBFs74/cL16PAqsIyJlKTlU9fPmcxNV3dR8+kslEfGPfz5UlGhOKCskAXJ2FHa2CPx3jDrWikk3Ua91csyW9ze5JTVS24DE9URkdmCZGHVKERmCl8LkMeDj/hxN5tPvNbUhEHzIrzRlkYjI6IhlMxGp5IV6XkR2BlREBonIZIybqxrWtRXBP1/ciI9uED92Z2/nQBAHyeUAuOf9m9rVtK6mKe4wi6XF1DAg8VVVHVOpgoh8EC99+/dV9U2Jf5mN2lGpJVcAo/FmOBRgODAf+KiIHKGq90QccwRwOZ5ArQTuwYu1VMVaJCEc21WolDpvhy8Ij/4q4dS6FktCjhr2AEcNe4DvDXswnQa4kmypgogMxBORG1X1DlP8d99lZT5fMeUrgY0Dh28EVBo9vQJvut0xqro9MApYBHwOiMzGqqqvquo3VPXjqvoxVf2mqv4z0N5T4i5mhaRG3Je3SLsJXc3cX5QH162VYskSosmWiufwTI9fAktVNZh87i6KYzcmAL8NlH/L9N4aC7wRnKY8gi1VdbG/YRIybqeqiXJnxXBg3I6WCUkLu7e1hddf3Jh/v/Qp/v3Sp3jvpU1Z/fLmsSKyz6Cvt7l1RcaPOiO1azcVhTnTbA8tS218Z4uH23vBpIH26u6vXYBDgD1FZJ5Z9sXLtru3iDyDNzrdz757N/AssBy4Cm/QYSWWiciVIrKbWa4AnhaRNYBVNfzFQWJf6RLFSERkjfA0jVFlIVrVva3l1JPEcdxahxRGgs/6d6VJyLochT/ffEJNE1qJKrOv9lxgo7/XbwciWjJMc7KIqurDxH/T94qorySMVxi+jSc23zfXeRiYjCcie9TS1mAz4nYktUj+nLCseMUWdW9rF7XESiSlVCLjR56eynWrsfPX6psV0WKphwlDH2nvBTtghkRVfUdVL1PVA1T1i6p6qar+R1VdVX27ztPWZ5GIyCfwHv5rich2gRN9GC8nS7KrV+jeJiLVureV+AFNF7qJAOtvkGgWSEs78YPsN56Q/JguGadm6RGyk9E+FhHZBW963U8ReM4nGUtSgdvidlRzbY3DM5E2AoIBobeAU5Ncudnd21R1GjANYOjwtewjyGKxtA9/HEn2+SUwCW+Ee6JBhSJyMXAe8A4wExiJ99z+FYCq/jDu2IquLVWdrqp7AN9W1T0CyxcC3dUqNayV3dt6mvEjfpB2Eyry6YOL7q3HbqjcBXj2L20XYUvn0IxeW23gDVWdoaqvqOo//aXKMfuo6pt4KVZWAlsAJya5WKJgu6r+RkT+G9gGWDNQfk7cMQm6t11Iefe2Y0TkFrwge7XubT3L+OGnZXoudE0YeROFJ66xItKLuCpcvt3NaTejPtIXiSQ8ICKXAHcAhU5Rftw6Bj9WsC9ws6q+VsGDVELSXls/x4uJ7AFcDXwFeLzKYX73toUiMs+UnYonIL8WkcOAv1Hsm3y3+QOW401UX1MaY0t22emQGqbXtVjqYMLQR7jhmbFpNyNL+D1egyPrFdizwjG/E5Gn8FxbR4nI+sC7SS6WNEXKzqo6QkQWqOrZInIZntLF0obubZZORuGJa8utkbm/mMToI2LyatmJriwZIANuq6qYkEStx5wsIhcBb6pqXkT+g9ebtipJu/++Yz7/IyIb4PVF3qTWhlp6iAoPfIkRETDjSOiMH6ulfjrardWkFCmtRETWFpGpgYSRl4nI2lWO+QDey/yVpmgDSi2aWJIKye9FZB3gEmAuXh6XDpzezJIVdjh0asX9kR1jrDXS8bjqdK6I+HTAOBLgGrzetV81y5vAtVWOuRZ4H9jZbK/E68VVlURCoqrnqurrqvobvH7JW6pqNkfDdTnjtzkt7SbUTLjXVmf0nrR0GocMfbQt1+mQXlubqeqZqvqsWc6m8qRW/jEXY1KoqOo7JHx9qzYg8UsV9pGkC7Cldxn7jctqG5hosXQC6YtEEt4RkV1NrNofoPhOlWPeF5G1MH+hiGxGoMdXJaoF2/cLbfu3UMy6FZIgrhbybfU69VgdfnwEQFzQXHxdUZjfP4nhx9sJryxFvr75YwDc9pftW3eRzhCSI4HpJi4iwGt4g8srcSbeQMSNReRGvJ631Y4BqgiJqh4KICJrAl8GhgSO6Yzb2U5aLCJ9W5+a6fEj7WbEJJvg0dJeMuK2qoqqzgNGisiHzfabCY65V0TmAmPxflnHqeqrSa6XtPvv/wNexwu0+/2KO+B2WrLCY9d7cZIdv105yJ6UKIvHzrBo8Tlwsznc8ZftWnPylHtkVUJEIrtD+gMLQ4PDw3UELwv7pqp6joh8UkR2VNVqYwYT99raSFUPUtWLTUbJyyo1yNJ8+raKnZws01RKJ7/Dd2K+Qgpzf55gbpLg71lh4aWTWHSxndMk6xz35MFpN6EhMh5s/1CVpRJXAJ8G/P+gt4CfJbloUovkEREZrqoLE9bvXazrKRnmNo05zBv17s9H4jP6iH6oECMBQL3bveAyKx6dgiPtS537pc2eBOB3z45o7okz7IsxvbPqZSdVHS0iT5pz/UtEBiU5MKlFsivexFTLzOyFC0VkQb2tbRZrr7FN2k2wJGDsN6Mtj6AbavuJgaB5NS0OdvmIYPFFk1h8oRWXLJLGOJL9Nl3Avpssas7JElojnRBHiWCViOQo9tpan4RJ85MKiT974T54Pbk+T3mPLoth3AcnVK/UKySw0LRWKy5QXSG259biC6yYWFpAZwxIrIcfA3cCHxOR8/FmVYxNHR8k6YDEv0Yt9bfXUgt9w05OuwmZZ/gJ0WKy5IdWTLJEvgu62YmbbOk0VPVGYApwAd6Egl9U1djJrIIkjZFYLA2x0yFTq85LEiRpKnqgYKEMn9zPwkutcFh6l7heWz5xnaRExAEWqOq2wFO1XrfjhWTvTZYC8OCKLVJuiSUOPxay44QWppPv/BddS6eQbbeV3zNrGLAD3jxP4IUiHoo7SFVdEZkvIp9U1b/VetFa3vssNTD+w3Y6lRLa9KDf9sRyF5d1b2WHXMafwlVpYrBdRK4RkVdEZFGg7CwReUFE5pll38C+U0Rkuen0NC6yeapnm55b6wGjVfUEVT0B2B5v1tlKDAYWi8h9InKXv1T/S7rAIski4lh9jqWFc4rYwYiWttA8LbwO+Clwfai8X1UvDRaIyNbAQXiz1G4A/EFEtlDVuPnYP4mXydfnfbzMJJWou+uwFRKLxdJWjnvy4M5OJd8kIVHVh0RkSMLq+wO3qOp7wHMishzYEfhzTP0bgMdF5E68Fh9AuWCF2/NgwraU0TVCstuQp7s2TjJj2YUdO7I9ivCkVtsfbhMv9grtHJDYCoSaemStJyKzA9vTVHVaguOOEZFvAbOBE1T1X8CGQDBP/kpTFomqni8iM4DPmKJDVfXJShcVkbcol8k3Au14Nu7YrhESS+dQmNRK4IlrkvfkqoaYNCk+25zSX3R3WW+jpRnUNtjwVVVNNMNggCuBc70rcS5wGfAdoh3CZS0RkQ+r6psisi7eBIQrAvvWVdXXKlx7KvAicJO53kHAJ4BleBNl7R53YMt+Xq0IJFWiW62RrkNKP8ccNhXR5jmdt53Sz7ZTrIVjaSEtHJCoqn9X1byqusBVeO4r8CyQjQNVN8J76Ie5yXzOwbMk/MXfrsR4Vf2Fqr6lqm8a62lfVb0V+EilA1v5nnYdMD6ivF9VR5nlbigLJI0HrjBD9S1dTLOD4zbY3jl0dOLGFgqJiAwObB4A+C/idwEHicgaIrIJXqaRsqy8qvp5k8V3N1XdNLBsoqrVZkh0ReSrIuKY5avBU1c6sGWurRYHkixdRKIUKTX+MLc5qb/0NcmlehJIiyUBzcqjJSI347mL1hORlXgTS+0uIqPwvvErgO8BqOpiEfk1sARYDRwd12NLVdUE2Wud3esbwOV4WYAVLybzTTNr4jGVDkwjRtJwICnM/SuG2WdENxH6oYrWNtLdWiaWltK8XltRZtkvK9Q/Hzg/4ekfFZEdVPWJGtrzLPE5FB+udGy7Q5BXApsBo/ByufiTVSQKJAGIyEQRmS0is//xj3+0ppVNYvza30m7CZ1J6NvQFGGQjs3IaskS2jG5tvYA/iwif0masV1EtjCDEReZ7REi8oMkF2urkDQhkISqTlPVMao6Zv31129tgxuhyfOSzFh6QVPP13OItVQsTaIzsv/24b2070nyjO1XAacAqwBUdQFe7LoqbRWSRgNJUdz73FbNbWRG6aZxJBYLwKR5iZ5RmaND5iN5K2KJfDkP8IGIaXVXJ7lYy2IkrQokWSyWzsfJwJO2bjqj6XPxvDz/wnMWrwO8JCKvAIer6pyIY14Vkc0oTmz1FbwQRFVaZpGo6sGqOlhVB6rqRqr6S1U9RFWHq+oIVf2Cqr4UqH++qm6mqsNUdUar2tWpWNdWKaKwYKpNxmhpM0ndWumLzUy8MSDrqepH8VxdvwaOwuuVFcXRwC+ALUXkBeD7wJFJLmbH+1o6moVTJ7HQztluaRNCx7i2xqjqLH9DVe8BPquqjwJrRB2gqs+q6ueA9YEtVXVXVV2R5GI2RUqrUGXmm9e25VKiivbwhBzDT+i385F0KJPnf41LR96adjNqIgMikYTXROQk4Baz/TXgX2agd0mfsrjJsMR0GIqbDCuItUgsnYfC/B/FWyGisPgia6VknVyHPJHL6AzX1tfxer/+P+C3ePGSr+MNy/1qqO6HzDIGz5W1oVmOALZOcjFrkXQ6qsycfy7jtj8z7ZZYLL1B+iJRFVV9FTgWvN6ywXg0sDxU92xT7x68ybDeMttnAYnmbLcWiaVrWXyBtUosTaaJMyS2kf9NWK+eybAAKySWTqHOH6adZjf7TJl/YNpNqI0Mu7ZEJMrLlDSC6E+GdZaInAk8BkxPcqB1bbUAdV1mvZ3o/iemb6tTykfLqzJzwXlNvU4WEYUnf2YFoRtxSD+XSK1kIP1JJR4HRofKrkpyoJkMayawqymqOhmWj7VILJmnmojY7r+WdpJx11aZ9aGqceNGvANE5gbqzlHVy83yZFSdKKxF0qH0giVSCwsvncS2J9oJrToRR5STF3yFC0fcnnZTqpONHlmVWD+uOy/EduXdqkpCRwHWrnRRKyQdQN+wk8GJd3OO2+5Ma1taOp7TFnyJ80fckXYzqpNtIckBH6S2kVVbJqhTMWWVFZKM0zfs5LSbkB4KT16R3G216JJJdprdDiXXIbESf2R7hnlJVc+p5QBV/WujF+1oIXnjvcUk7J1m6REWXdzmeIliR9U3kdMXHsC5w+9MuxkVETfTSpLKt9E6RFrE+A8f2vA5+oZOaUJLLC3tAmxFpGl0REbg7Cdt3CuNi3a0RZJVxLH6bLF0K1nWO1V9LY3rWiHJKH2bn9j0WRY7CVGYe6Xt1mvJIBkWkrSwr84txM7ZXictEJEl51lR6mXufm7bpp0r4+NIUsEKSato0JqYsfySRPXGjzqjoetkCVGY+4tJzP1Fax76S8+dxNJzrKBYGqRJMRIRuUZEXhGRRYGydUXkXhF5xnx+xJSLiPxYRJaLyAIRCY9eTxUrJJ1IZ/SUrA2FOdNa/5Df6gzbPdjSAOqlSEmyJOA6YHyo7GTgPlUdCtxntsGb4XCoWSYCVzbjz2kWNkZiSR1RZfbVsYNxLZa6+e2zo3BwyTUp3NjMcSSq+pCIDAkV7w/sbtanA38ETjLl16uqAo+KyDoR6eFTw1oklp7CurZ6g9v+sj13/GW71pxcNdkC64nI7MAyMcHZP+6Lg/n8mCnfEHg+UG+lKcsELROSbvL/NcL4j3w37SZYLD3Frct3aOn5awi2v6qqYwLLtEYuG1GWmZB+Ky2S6+gS/1/dNDCepG/zEyuc1/sYP/L0us/fy1irxFI3rR+Q+HcRGQzezIbAK6Z8Jd50uT4bAS/WfZUm0zIhUdWHgPDgmP0pTpQyHfhioPx69XgUWMe/mb1I36aTq9YZP+IHbWiJxdJZ3LR8p5Zfo4nB9ijuAiaY9Ql486375d8y3puxwBtZiY9A+2MkDfv/RGSi73N847WKCSktllieOttaJZb6aJaQiMjNwJ+BYSKyUkQOAy4E9haRZ4C9zTbA3cCzePOtXwUc1YI/rW6y0msrsf/P+BmnAQwdvlZmfISVGL/u4cx8LdEkZYmsEYvFkhKKH0hv/FSqB8fsKsuXZXprHd2UC7eAdlskHen/a4geTnNisdSCLUuX+QAADcZJREFUq5L5zL9gR7ZH0W4h6Uj/X7vo26R3x1Jsf7gdKGhpjBueGdueC2U7+28qtMy1Zfx/u+P1pV4JnInn7/u18QX+DTjQVL8b2BfP//cfoPEc7JaOQa3VZmmQ6c/s3Ja34g6Y2CoVWiYk3eT/awbj15vIzFcb6UZuARh5XD8IqNGeBf31Bc23PLPfziXSJUx/Zuf2XUw16xNbpYId2d4OKsy37tP3ye+3oSGdy6hj+z0RCTHieOsSs7QZ69oqIyu9tnqavo2Ps0H5Cow6JsZ6kPp+r1udYa0RS/1Y11Y5VkjaSN/6RzDjHz+n7xNHQS7niYcVkEi2OyrwsK9yi4Yf38/CqXZcSC9yzdO7Ni0hYyIUsK6tMqyQtAkRz4vY94lMjSPKFNtP7PdiH0LNTtfhk/sLcZNFl8SLik0j3x1csWwPcrgMTOM9zOpIGTZG0k4SxEp6lkAAvZ5jg2w7xYpFp5HH4fwRd6TdjETYcSTlWIvE0l1Yrba0GNtrqxwrJJaeYOsfBNxmlsxx4YjbE9e9YtkeLWxJFXqwR1YSrJBYupJtTi4db9IKJ664oLnmn9dSzqVLxzFQ8uRwU/UQewMSrZKEsUJi6TrqjrVYLEmoP0V812KFxGKxWGrAWiTl2F5bFkudqP31NI2TF3wlsvziJX1cunRcm1tTgdbPkNiR2J+CxVIDovCXE3o3S7PFy7WVZOklrGvLYkmKCuoom102Ne2WWNLEurbKsBaJxZIEheVTbBqWdnLBkn3TbkI52vI52zsSa5FYLFUQhWdOnsTml0y141BaRHgcyfmLP5/dRBDWIinDConFkoDNL7YZgy0GqyNlWCGxWKpgx6W0FrfDbrC4Pea3SoAVEovFkipOJ2U4VJo2IFFEVgBvAXlgtaqOEZF1gVuBIcAK4Kuq+q/mXLF12GC7xWJJlU6ySARFNNmSkD1UdZSqjjHbJwP3qepQ4D6znXlSERIRWSEiC0VknojMNmXrisi9IvKM+fxIGm3rOTroZdDSG5yz6AtpN6EyqsmW+tgfmG7WpwNfbEqbW0yaFklXKHHH0zkvg5YupaNcW9BMIVHgHhGZIyITTdnHVfUl7zL6EvCxFv0VTSVLMZL9gd3N+nTgj8BJaTXGYrFkC1cdcmkP0KgtRrKe73ExTFPVaYHtXVT1RRH5GHCviDzVpFa2nbQskrqVWEQmishsEZn9xmv5NjXXYrG0krhcW0GctEXEIK6baAFeVdUxgSUoIqjqi+bzFeBOYEfg7yIyGMB8vtLev64+0hKSXVR1NNAHHC0in016oKpO8/9j1l7XTgZhsXQLpy88gNMXHsAZ296VdlMqkNCtVcW1JSL/JSIf8teBfYBFwF3ABFNtAvDbFv4xTSMV11ZQiUWkRIlV9aVOUmKLxdI4uYC/6KxF+zMwq7E7pVkj2z8O3Cki4D2Hb1LVmSLyBPBrETkM+BtwYDMu1mraLiRGfR1VfSugxOdQVOIL6SAltlgszSP1GEgSmtBEVX0WGBlR/k9gr8av0F7SsEi6SoktFktzcTLeJ91ObFVO24Wk25TYYrH0GFZIyshS91+LxWLJNqqQ7wD3W5uxQmKxWCy1YC2SMqyQWCwWSy1YISnDConFYrEkRYEem489CVZILBaLJTEKamMkYayQWCwWS1IUG2yPwAqJxWLpCDKRtBFsjCQCO7GVxWLJNK46TN5qVmaSNrZ4PpKOxAqJpak8cusJZWWdNt2EJTvkEaZsPSPtZgRoTtLGbsO6troFF2Y9eVZJ0T5jz2n6ZcSFh+6aDMCuX7q0WK7wf7d5IvLnm8vFxGfHb0/1VhRmX3N8yb7R3+uPv7CWri/40aSyKtueWOF4S0eRV4fTtvl92s0oRwE3I5ZRhrBC0qmoMnPBeRWr3PPoGZHle+3+w9qu5SoP/e+UsuKH75hc23mAx687Pnbf3F+Ui0MtLLok/vitf2BFxtIkeszaSIIVkk5ClZlLahSBCO7746kl23vudWFkPVH4493lAtKJLDkvWmSGnWMFxlILNkVKFFZIskTeZcbfppYU9Q0zU9e7yoxl0Q/8Rrn/Pu8ae+xzUbFQlT/O7P6Zjped4QnMsHOtoKRJHofzh9+edjOqo6B2HEkZVkgSkleXAYOXA3BvG79HrRKPKB64p/uFI45lpxctli3Ot6LSTlwVLh55W9rNSI4d2V6GFZIQLspHN1iZdjMsKfL0aeVusM0vtuJiMdgYSRk9LyR5hY03ejHtZlgyzvIpk9j8kqnVK3Y5rgq/+8xPysoPfnRiCq1JAVXbayuCnhWSPMKmVkAsVRh6YT8qgL9YIrl57LTI8omzJ0SW94+6pZXNaS3WIimjZ4QkjzB8Y+uysiRjix8GBKTDURXmfz7ZmKLP3ndiU689bcz0pp4vfRTN59NuROboGSFxtQueCBZLFArLvhw9ZqhWHtrrkqacpxUct+Uf+OlTe6bbCJtGPpLMCYmIjAcuB3LA1araULelPA6jP7miGU2zWFJBVXju66ek3YxMcMyW95dsX/P0ru1vhO3+W0amhEREcsDPgL2BlcATInKXqi5Jeg4Xh50+tbxVTbT0MiqtTRymsGJC73bBrofvbPFwyfZNy3dq6fUUUGuRlJEpIQF2BJar6rMAInILsD+QWEh2G/J0i5pm6SWePjV5upZNL7e9ubLC1zd/DIDb/rJ9ay6gdmKrKLImJBsCzwe2VwKxrxhrr7ENe28yu+WNslgq8exx8fnDLOlw4GZzIkqbEye1wfZysiYkUf/TJXakiEwE/E7rb4vIP4FXW92wBlmP7LcRbDubjW1nc2m0nZ9qtAFv8a9Zf9Db10tYvRPuaVPImpCsBDYObG8ElAz2UNVpQKHTuojMVtUx7WlefXRCG8G2s9nYdjaXLLRTVcenef2skrWJrZ4AhorIJiIyCDgIuCvlNlksFoulApmySFR1tYgcA8zC6/57jaouTrlZFovFYqlApoQEQFXvBu6u4ZDo3AzZohPaCLadzca2s7l0Sjt7DlGbN8ZisVgsDZC1GInFYrFYOgwrJBaLxWJpiI4VEhEZLyLLRGS5iJycdnuCiMgKEVkoIvNEZLYpW1dE7hWRZ8znR1Jo1zUi8oqILAqURbZLPH5s7u8CERmdcjvPEpEXzD2dJyL7BvadYtq5TETGtamNG4vIAyKyVEQWi8hxpjxT97NCO7N2P9cUkcdFZL5p59mmfBMReczcz1tNb05EZA2zvdzsH9KOdlpiUNWOW/B6dP0F2BQYBMwHtk67XYH2rQDWC5VdDJxs1k8GLkqhXZ8FRgOLqrUL2BeYgTdIdCzwWMrtPAuYHFF3a/P/vwawifle5NrQxsHAaLP+IeBp05ZM3c8K7cza/RTgg2Z9IPCYuU+/Bg4y5T8HjjTrRwE/N+sHAbe26/tpl/KlUy2SQk4uVX0f8HNyZZn9AX9yhunAF9vdAFV9CHgtVBzXrv2B69XjUWAdERmcYjvj2B+4RVXfU9XngOV434+Woqovqepcs/4WsBQvxU+m7meFdsaR1v1UVX3bbA40iwJ7Areb8vD99O/z7cBeImLnikiJThWSqJxclX4c7UaBe0RkjknpAvBxVX0JvB838LHUWldKXLuyeI+PMW6hawKuwdTbadwq2+G9RWf2fobaCRm7nyKSE5F5wCvAvXjW0OuqujqiLYV2mv1vAB9tRzst5XSqkFTNyZUyu6jqaKAPOFpEPpt2g+oga/f4SmAzYBTwEnCZKU+1nSLyQeA3wPdV9c1KVSPK0mxn5u6nquZVdRReaqQdga0qtCVr38+eplOFpGpOrjRR1RfN5yvAnXg/ir/7rgzz+Up6LSwhrl2Zuseq+nfzoHGBqyi6W1Jrp4gMxHs436iqd5jizN3PqHZm8X76qOrrwB/xYiTriIg/cDrYlkI7zf61Se4OtTSZThWSzObkEpH/EpEP+evAPsAivPZNMNUmAL9Np4VlxLXrLuBbprfRWOAN32WTBqF4wgF49xS8dh5kevFsAgwFHm9DewT4JbBUVYMTkmTqfsa1M4P3c30RWcesrwV8Di+e8wDwFVMtfD/9+/wV4H5VtRZJWqQd7a93wesF8zSeH/W0tNsTaNemeL1e5gOL/bbh+W/vA54xn+um0Lab8dwYq/De6A6Laxee6+Bn5v4uBMak3M4bTDsW4D1EBgfqn2bauQzoa1Mbd8VzpSwA5pll36zdzwrtzNr9HAE8adqzCDjDlG+KJ2TLgduANUz5mmZ7udm/abu+n3YpX2yKFIvFYrE0RKe6tiwWi8WSEayQWCwWi6UhrJBYLBaLpSGskFgsFoulIayQWCwWi6UhrJBYOh4RuU5EvlK9psViaQVWSCwWi8XSEFZILB2DiAwx82pcZeasuMeMgg7W2UtEnhRvPphrRGQNU75CRM4Wkblm35bp/BUWS/dhhcTSaQwFfqaq2wCvA1/2d4jImsB1wNdUdTgwADgycOyr6iXTvBKY3LYWWyxdjhUSS6fxnKrOM+tzgCGBfcPM/qfN9nS8SbJ87og5zmKxNIAVEkun8V5gPY9ndfhUm9jIPzZ8nMViaQArJJZu4ilgiIhsbrYPAR5MsT0WS09ghcTSNajqu8ChwG0ishBw8eb5tlgsLcRm/7VYLBZLQ1iLxGKxWCwNYYXEYrFYLA1hhcRisVgsDWGFxGKxWCwNYYXEYrFYLA1hhcRisVgsDWGFxGKxWCwN8f8BneFN+uFXL08AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ds.TLONG.where(ds.KMT > 0).plot()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x2b1865175cc0>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de7wcVZXvv6v65IQAgQTCIySBBIgiDw2YQZTxisAocFVEFJEZBMWJcskMgzozMDIiDoyMD7g4OIxRkMcVQ4aIRAcFxOeoAQKEhDcBQsiDhEBeEJJzTte6f1RVd3V1dXd1n+qurj7r+/nU51Tt2rVrnTp96tdr77XXFlXFMAzDMJLiZG2AYRiGkS9MOAzDMIymMOEwDMMwmsKEwzAMw2gKEw7DMAyjKfqyNmA4TJgwQadOnZq1GSOG17YPeTtSLtu5P9cfodwwWHRjy0cVmv/uF7TVyrV55sEHH1yvqnsMp40pMka3Ef+3iLKegbtU9YTh3K9byfV//dSpU1m0aFHWZuSaXy17GYCCSIOa4Ph1wu+b4H0W9w5619Tdh21fr/PSptdL+3vvulPFcZhw1Hz4T9XKq7+o1e2EbehVROSF4baxDZdTmZio7nd5YcJw79et5Fo4mmXxqo2l/RmTxpX2H1ixobT/Z/uO76hNafPzJ9eW9k88aK9U23ZVcUSo8eW3ij8ufyW2PBCUXnruazZWvvDr6bADsd9Z19UQDYBCY11PTEHi79/IBgAJ/WJ77LJjekblBKGJv0UPT5HLtXBsGyzy+EubgcpvZK7/F6soi0x0fPDFjVVlUPmyC16QQb1BV0PnKq+Na6tY44MTVzctRjnlT/WCx1+qOl/yGqTyeNB/ZgUn/r8iKC6GbA97KW6xsr1SnZivxMEzHuXkv6tk1YbXqfHIqnBi9iX6vKLXDEMwwpdGP3HR+9TDrfFxfWXL1pr3C+4Z/T+BStGK86SCsuDS8CRlF9h3t52r2nxu/ZbINeX3QPQ+w0GA/qR/lGI69+xGci0c9Sh/+LydqAiA9xKsOK5RJ08M1vovjxCIWqUQgBujdo5Ixf9AIC7Bc6r3f+QW47u4ovde+MKrHLXfbolsz4o4kQiOox5G+FeNikPc86r1COs922i7nSIu20Stj130i4irZdFS1dIvntCJBWDlq6+V6qu294tYFM/jyOa5dxM9IRwlkUArPkhFN7QfEomwQAQvr+CDH3xDqhQPKsoGiy5FhdPeuk/bfqdmmLdkdVVZPXc66hWA9w8e/ecP6hVVK9oLxKV0nkoxCb8rCiLlZ+mWTxScctdXwMIXXi39U+al6yowP+o7hV/q4ecRffJRYUgiMnHtlOt34qVWfY+4l3fc67zSi5Dq8hr1KoXCE5+y8JQ/ewXx2nABR/3yNPuMJN1uw7zSE8IBZdGAamFwVWuKRVQoipFvMHFCEnDzQytjbdk+VOmjfubI/Vr7pRISCNjcR1aVylyt9xKJ+UeKfOXzhCRcL/TSD74lhl78wXMKBCjatRUISLl+vHh0gqBbI7hrcP+4Qeda3+prDVAnEYyK8pj2q8Qk5v5xz6zRY5Q2fjOXGgoXvWVB/P8lkYpPoUS6pbx6oMHfpiQsiojfrogvIuoJRUhA8E6n1kVVshPzOCDvwiHlvszA0yi6sL1YLHkbnnB4LzbvuNJzgEqPwg2PY0Q8jXJ5/KfRreGvX7tweVXZQNHl/KP3b+a35fv3l4NCigqffUelIAXhsqMKjv9iKXsKZxw+GYBbHl6J64RfbmWvAkKiEPEqwm/HSlGK/M6BADmVL1DPaymLUTDIHohHxb3ayLKXt7Q0bhDXJRUnLo3Eop5HUVm/sl700VSJgLpVfwrRGh1AtcpL55t/28Y+0hp/z4KUP0wa7AsgDiAlsQh3NwcWlb8A+p6Flv9PA2+k6J8Iex6uNP87xf5KmMcBeRcOPNEouspgURl0XQaLyrYhl8FAKFxl0FWG/O4lCAmKqxUD2OFurTBxA3yNztUaZA5z5e+fbXhdWFwG3WTfzoshr+rMmfsCnugEglIM/dIF0ZKopMHPn1zri0K1BxLnfQTiAZUCHXRbpdFl9cy6LdUv3kideoIQ96esV6fsyYTbj28r6vVEbSkJhFL1wq8QhqgYxAlLrWuB/l0rI0cH1q+Evv7a7adE6VcNiUmwH4hKSWhEQByvvFRHGHK9L4fqfzEsujAo4v0fuOqNz6XmeYh5HORcOATh9YEirw0U2TpYZNuQy/Yhl62DRQZd1/vQ+NvAkO9d1ImMqicQ0FgMap2PK4+WBcdRG6749TOlcwVH6toQDIw7GhxX/rMHgvJpX0zaQRACHBWQ4I0Z9T6i4hEQ/HM+sGJD0+Lx5NrNobGH+mMGUSGOjlVUvvzju6VK3kTMdfXrVHeRQVgs3FBZjEjUEYe6ogLgDkHk8zG4bnn5BS4ODA14u0PbvbLCqNrtxXkpCcUmzjNSif7lgsohEfGPC4U+r1wc1OmDQoFi/2gGii4DRWWgqGxLKcJJgFEmHPkWjseeXsE5V/+BHXcdzegdRtHX71AoOIzu9+I2ghdtwRH6Qi/fgPIL2Ykpq64Hlf2btV7+QOl+dduqsR+1I3j5Fl0tlf/9ew6sqH/l758t/a5xXWmfOXI/vnuf19X13fteqOrmakQwfhK87BoFBiSZQ/LbZ9fXFI9AXKA58Xhy7eaGdWr920dFY9L4nUrzM6pe7jUEoZZgVNap9C4qup38l2hDoWgoKjWEw42/bnDt8/6ggAuqlOLowvcJBCTaZp37SY171xK9moIX7Ptx3+qGjt0iWiyW7qFDA8jANgrbtzF66xaKm15hwxPlbt7hIGJdVZBz4ZBCHxtWreT1zbsxesxoRo3uo6/fYfQOo3D6PBERBxxHcAoOjiNISEzC3+CDn32Rb/VenRrCUoyUhaKGiqF6Q/6Lv8+p7MKJ7ofvGxxHywEufO90Lr/36SpPJOytNPKOAhGB6rGSJASRXKe9dR/mL11d8TI85dBkM2vjyGKwPA6p0XUUPldLNIg5jopG6QtzPdHQmBdt3Eu54cu20TWhQYJG7YeP4+rWEodo3SZEASgLw9Ag6oZEwi2i27d54jE0CEMD6NAgxS0bGXp9GwObX2f7xtfYvnELaWFdVTkXDoCBrZtwhwYYfH0MhdFjcPr6GT1mdEg4hEKf0DeqgDiCI4LT5/hiIjj+f7c4gohUiAvQksBUl1d6LHGeRuChDLlaajvO+wG45K4n6e9zqrrdwuIR9UjA67oKf+iDl9h373uhPCEwcq+gW8t1FceRqoHseUtWN/0N7JfPvOy3UXu8Ayq9jgdf3Mjbp4yLa67E4y9tbmrgu2ZXkV+wasPrTBpfOwVHeAJcdDyjVtdU1fhFvRd1ElFJIhTDEaA4T8CN2NfItoggaA1BKJ3z76F+Vxmu69ULxCH0E1843KFB3IEh3MEhigODDGzZSnHbAEPbBihu244mTXfQAME8Dsi7cIigxSJF/wOmbhGnr7/0s9DnICIU+hz6RrmecDiCOFAoOJ6AiCcUQThh2EsBT1AK/ttMGohJXJdYnKBEvZaCCAPANae+lfNvXwqUBQTg6x88pPQrX/Tfj9Pf55TGbOK8FPDGRgLibD3/6P25duHyUFSVHwLpC0SUIFIFKAlI2DNI01OItuU0IR5xJDGrVvjtS5ter8jf9PLmrVWeRLR7KnwOGngZ1BGMJN/SI2JR06OoK1ANRCKJQES9BSiNkYS7kbyfRV88IkISEgivjlu1r4MDJa9DXRd3YIjitgHcYhF30BMOt+gy9Po2ioNDuAODaNHFTVU4TDnyLRzgf4DK4gEgTgF1i7hDBZy+ftyii7paEghHhCHHpa/fqfAyAFzfSwlEBqBYcCuEpBi5pp534r383bpCEt5/Y7BIvy9UA375+bcvLZ0fHRGNwMO47MS38NV7niqJR9RLievyOveoqbHP9Pv3v1AlHmceEQrnLQ14+2G1gegkDOM8fvoeJa8jOlDeiFbFA5pL3+FI+UUf5G+qJQzRc44kEI2aL+E6L/mocNTrggrdp+YYSPQ+UU8ien09oYjzIgAd2FZuy+9aIiwaDYSiNH7hHxe3e96Duq73s+gy5HsUxcEhTyQGh0r7Qb20EGki5UjddmQKcBOwN17U8BxVvVpEdgNuBaYCy4HTVHWDeB+wq4GTgK3A2ar60LANaZHcC0cY1y0iwbcZH3WLuBQoFl0K/vdK1wHHheKQ4jheGJ+4nhi4jhcI7riChsRBXS9wwwG8CHEH/G/nQ74oBVR0G/XF5WMKPshlASm6yjlzH64Y3I77ZlNvXOTLf/FmAL56z1Olsmi9RpFjUH/CYiASScvjuOupdYnClb12q+d4BMkqw4kqO0Ujs+udj4uWqjiu54kkGXzWau+g+h41vIuwZxE9H9P1VCEYQVmxLAiE/w8DwQjqhz2LoG7Q9QSV3kWw798zLAbFgcHyfiAUxbJQhAUsTVLqqhoCvqCqD4nIWOBBEbkHOBu4V1WvEJELgQuBfwROBKb72zuAa/2fmdA24RCRHYDfAaP9+9ymqpeIyA3Ae4BNftWzVXXxcBQ1+ICUPY0BxCng0I+Ltw/9Za/D9xgc/1tu2MNwfDFQvxsp6KoS/9jxr3fdoic0fh3xhSY8KA4hEYntziqLQHgAvnQNypxTZzT97AMBifKN3y6r2I8bB0lCWCSis8sLArc/ugao9KgcET7wlvhIq4rJga7UnVEeHgNZusb7CB02cVceW7M5UZdU2B6IzJloooHwgHjtwfKYMZR6noZ/HNvdJA794/ZsaNfgS8/WaD9GMGLGKlr2LuLGKNywOLihugPVUVExHoaWuqDLguD6AlEcHCyVhUXDDQtIyt4GpDfGoaprgDX+/hYReQKYBJwMHONXuxH4DZ5wnAzcpN7U+oUiMk5EJvrtdJx2ehzbgWNV9TURGQX8j4j83D/396p6W6T+sBXV8y5ASh/0as/DAc9ZcEFcxXWo8DDAKXkgUB7rcPzjwAtxVEr/W+L4YqNSui7srYS/8QciEkRbBZFWBSfwUrzw4Fqewbm3PVLyRK459a3NPJ4Kb6VVgi6rmx9aWfK4gNjxjqhHFKR8L83vcKSp7qMw4fstXbOpar7GcAjbFDf2ITH1atlX0UUVJxo1Io+C/SRiEWbU3geU9odWP1XhrVRFPsV1m0UEo6ZnEbQT/l+r5VmUjv26Md5FtEsKYKdPXMyWm75S08uIdlkF+274vikjbZgAKCJTgcOB+4C9AjFQ1TUiEnwAJgEvhi5b6Zf1lnD4yviafzjK3+r1k7SkqOoWfY8i/lxp3/H2XQqVLzz/BV/uvgoG0b2uKdcFhtyyEES6sRxHwPVehIGAiOPNzpZQvWbCbYOurILjcPYtntN1wxlHMGve4oq6s+cvAZILyIXvnZ6oXhzXL1oBUDF5MBgwjxss7yRpikbLNkS8DYgfDG8UQistCkYsvhDECkasiFWKhIZe4sHPwqHHMfTgf1d4FkBFBFTp9oMDpWvDAqS+dxInGBq+N1R4EG6xWCUUbhu9i1o04XFMEJHwSnNzVHVOuIKI7AzMB/5OVTfX8XrjTqQ2H75Z2roggogURGQxsA64R1Xv809dLiJLROQqERntl9VS1Gibs0RkkYgs0gEvHFL9D18wUN5oC+O6ivpb0O2irvrl3vhHuA5QXU+1lJzNDdXTUD2gYiZ7cBz8HIopD3PO3IdL+9EJfoGAtIvwnI/rF63gxge9P1P0H8hVL73LqYftwymHTuSUQyfyoYP3Lp2Py/1VykpckVgy5V8gR/TvOqEq/UfL1AuhpYaIxUVH4XVJ6eAgQw//ItmtQ55GXBdXxX6pXuW9N//gy+x6zmWRdmuLRuBttBNvAqAk2oD1qjoztEVFYxSeaPxQVX/sF68VkYn++Yl4707w3odTQpdPBqrTYneItgqHqhZVdQbeL3mkiBwKXAQcBPwZsBte/x0kVFRVnRP8IaQ/+QpkUWEJv/DDwhDe4sTDHfIitIp+pFZQr1TX9cSjGKoXbr9YZxuqKnMp+v8IVfXbmOk0jmgCxyBhZNy5KB86eO/Y8Y1GA/WVAhOIahJrW6ccWhufjyraTRU3ZyM8X0Pqjh9450bvsltpq8XgyysYXLe8tCWhb/Ih9E0+pPZYRnQLBqD9b/+Fg96NDg6UB6mhMpQ2mIA3NFA+LgaD2p5nUfJahgbKIhLs+2JUasctj1G4A0MlAanyNkKiEUeVmKQZVYUXVZVkq9uO96G6DnhCVa8MnVoAnOXvnwXcESr/pHgcBWzKanwDOhRVpaobReQ3wAmq+k2/eLuI/AD4on+cqqJGPYugO6ty3MPrRvKEwXtRuP6AtRP8n+B3hTha+t8JXihB91Yw6F7ELXVdeddqxZhKcG2zva83/9URTV6RHuG5HsEcj/DEvLj5H/OWrE68Vsnx0/corXueFQ7NDYwnpVYU1eidd010/cD6lX47w3vxFfZ7G8Xli5NFSoX+P4YeubtCMLzqyccyKsZFwsITM0YiTqEkGgGBtxEVjTDRSKpaRK9rlbQGx4GjgTOBpX6vDMA/AVcA80TkHGAF8DH/3J14gUPL8IKHPpWKFS3SzqiqPYBBXzTGAMcD/xaMW/iK+2HgUf+SBcBsEZmLNyjesqJ6EU+V4x7BB10KhVA978UeFg/xRcBVjRUP7zpvgFn9wXVHpCQe3gvUexGpWykewbVx4hFNGTL/05lF2pW45k/Ph1KhS4VAFGuMadQb5wi8jp8/ubbCmzj2wNbEQzXZ5L5aNraqFWFvIyAaSVWL0TuNTXSPgVcj35n8F/yoPacmur6KWtFcNcYzKgaX40JrQwLT/+7TAdj+q5tiB9R3OOnc1mz22X32N0r7a7/+N+UuqohYuB0a50hjcFxV/4f4XhaA42LqK3DesG+cEu30OCYCN4pIAe9dOU9VfyYiv/JFRYDFwOf8+k0ramzfVsxgebRMQ99yvOPyPIxgPCJ8HJ3LES5vJB7le2pJoEAq2oXKwfKfznpno1+9KyhFT4W9jgRfx+ISIAbi0Q15quJIYc5XYqpEIyUC0XD2n4m7bGGpPHjRFw4+huKj95YviH5LDx87BfrfeWrF6dHHfjJdgxMQ50nERVRJSuvbi3Rm3Zhup51RVUvwQsyi5cfWqN8WRQ2LRtTrKIXoBsLgv9S9c5WC4OJFYKlLqduq5K0EdUOf4T/8U9WXhlwS1y1Vj8CTaKa7KuDYA/fgt8+ub+s/ZrMtJ+3CUmDsjmPY+sa2yoy3cWk8khBJIdKyp+FTmFbZ3ekceFR8vUPrf24H/jQf3CL9R582LHtaJW5sI4mXkVZXFQhiyaryP3O8QhiKxYquqKhoRL2O4DjcZVXRfeULQoFq8ZACsXUfuPgvOvSbZ0d4PCPqdUBlAsRW12WPW6Ojpj1oqiG5zU0k9H6O3XFM48rqsn3LxtJ+QNWgeFQ09tiXbiHqZXSaiRd9p7S/8pK/BsDp72Pypd8D4PkvntnW+4tAoT8+/H8kkXvhAHj1F1+uON7zlCtjxzkCwgJTKS5a+jYR7b4qiQO+eKiWvo2Gz/cSs985jWv+9HxF2o+kXkerXsN7DpjA7597paKtNPqUk3Q1NVMnkSfSgrcR7qZScRB3yIuiUmXUXtMStzMSCMQizLRv3sxzF5zRvrEOwTwOci4cM940iUW/vKyqfN3tn695zeS/vC7WCwlCSp/9z1Mq6h90/h1QqByTcP3xiFI0llM9ZtErzH6n97KKrpvuugpOezKFvnv/3QH443JPQIopicdwSWJBXTMbdF+FRWO4kVQjmf2vuiX+xHUpfIZEEo3j9Tq5Fo7hoMUiK27w3Npps6LZT8o8efXJvOWCBaVv0HECUSsVeS8RzaT7/fvLkwLjPAxXldPfVjV/syneNXX3knhkSb0/baIuKqjfNZUA8za6h7QG2vPMiBOOlT88p6rs+TkfrXvNE1d9KLZ8xj/dWdVlNVKIy6Ab5K+C6gWhWuVdU3dn4QuvptJWmuy6U7xgjNlhB97YVk4lHuc5bN/8au3Mt+FrzOvoOkQwj4MRKByt8JYLFgCUMukCFYtCQbmra+Yld5XqPHDJ+zpvbBfQLu+rHfmw6n13TDKOsWXrGwiwc8jziK7s1xSh6J/hRlIZ7cHGONqccqQXifMqgpxUvd5d1Yggc24v0OgfY9Prb7Bl6xsttV3RVRXjXag4XRVJZZQREQr9hURbL2MeRwKaGfhWV1l06fubav+Yq35bWiMEqFrz/M7Pvas5g3NOeC5HwcnfhKsdxoxh29bXywW+KIweW154KkhkOLDhJa9AHPrH711KM2J0KUJPBsE0iwlHHaafd3vlsrIxXSXh2eQBb//nX/Dgv5xQs90//7dfV61rHkcgHh+c86dSWV5mlSfh50+uxZGySDqSfpRWOzQnrgsrmPgH5cl/wVhFvdxU/eP3rjye0DteW28iOEknGPUwI0o4Jn2inNV41Y9m1ax3wOfml1YKrEc4mio86xwqxzoKoQ/awouPj70+elxr0aUPf38hP/lM/KzfbuCMw9N98b1rqhea+8CKDam2W+tPW7O8VjvU8IjEYfSOO7VimtHN2DwOYAQJxz4fv7Zi5viUM28A4MWbzwZg6mduBaDQV/lISnM2YiYERgle/OFzSbpZwhMP49oP8lgFInLq9d6yJt2QCHE4hHNW/fKZ5EkOiy7UmNtpGG1FTDiAESQcYcICsu/ZNyNOoeYqgkmI666qxTv/9V4KBafqmnBYbz3CHsjHb7ifW88+sllzc0HQLfjH5a943Vk5GOdQhR3H7JC1GUabsa6qHhCO8e+7pPTSF6eAFAoVx0HaEaevvyp3lTOqv+n7BS+0egPmtc7V81KSEF5mts//+fEb7qfgCLd88s+S/gq5IOqpZTV7vNafJlwefCZMNHofEaEwyoQj109g8dOrmr5G3SIrf3hO7ETAqrqRpV/jqLf6XbAcbfT6oN0k7UVXyYseV69Xnl/qzZ9st2gM5zHmLeprJPLkZz6cTkMCUnASbb1Mz/12cbn4w6y+NX5RmSDleiMhaET0+vBLpUpA6swRi9aNCsZQj81SP376Hhw/fY+szWhIve7E8Ixxo3t4atZHUm3PKUiirZfJfVdVUtQtsua/ZleWRdKwD699L9V6lGZmO0fHOcJdU7UI6pz5/xYBcPNfzUxudJfxy2deZlROPShzOrqTZ86tn06oacTW44ARJBxRwllyW6WVcY443AQCUY/wtWff8hAAN5yR3TrlSQnmcUBvdbkZvYv4XVUjnVwLx4w3TeLZJurv8/FrWX3ruUz6xJyGolHvxe+qUoiJ7K8XGdWskDgx4biNjgs5ytr5syfWNlzXo1cZ2LS+1E/ZP27PRNcMrX6qdE3fpLe0zbZeY/q1oczXqaRVxwbH6YExjg13X1q1kFMcgVCEJwE2Q6sD5DXb0+h4R7pjFufMfTjV9jpFsU7gQNaEdX/sjmMqUqqreplxA7a/vqW8bdnobZurs/yWUo40wdCqJ5q+xkgH8WeOJ9l6mbZ5HCKyA/A7YLR/n9tU9RIRmQbMBXYDHgLOVNUBERkN3AS8HXgF+LiqLm90n/Hvu8S7X0MPojjsrqmk1AuxVVcr3kDNzAGJehpDrpbCcr3zLgXHyVW3Tzsy3rYTVWXczjuWjndOuh6Hz/bNrzJ61wkMbFxnadPziE0ABNrbVbUdOFZVXxORUcD/iMjPgc8DV6nqXBH5T+Ac4Fr/5wZVPVBETgf+Dfh4G+0rUU9UoqlEyte0Z8W/pBMB88rtj66h4EiuxGJYJBCHgVdX079b7bXZh9Y8U1VWfHFpRQr2wn5va80+ozlsjANoY1eVerzmH47yNwWOBYKOxxuBIMD6ZP8Y//xx0iNvUHcY3S+N5nHUuiZJvU5z+6NrsjYhljQ/ZtveeKOcGbeGaAxsWl8+8Os0lRU3pt3i8sXJrzeGgSCOk2jrZdr624lIQUQWA+uAe4BngY2qOuRXWQkE64tOAl4E8M9vAnav134rEwCHSxpjEa2MiRiVDEeMG9HqP8W2N+qvzyHqVq4IWF4FrMU7+vjPovj8Q8Nrp4dJawKgtwKgjXG0NapKVYvADBEZB9wOxIWDBG+AuK99VW8HEZkFzAKQHXZpaIPrFktpR5KS1hri0e6s4YToNoqs6nbmL109crqnGtGg+2pw7fOl/Yq1xsWhb583A35XFZREw8ZLapPqBEARnP5cB6OmQkeegKpuFJHfAEcB40Skz/cqJgOr/WorgSnAShHpA3YFqsJQVHUOMAegb5eJmX91b9dYR7PkTUjyRNZPtW/i9KqywpTDKL7wSAbW5IvUJwD6XVUjnXZGVe0BDPqiMQY4Hm/A+9fAR/Eiq84C7vAvWeAf/8k//yutl9TJpxPRUs14IM0MbkfrdosIdSNZJTmMI1gyNhyO+8a2bc0LTMRLGFi/sqmFnGxAvD7piwb+4Ljl9G+nxzERuFFECnjdxvNU9Wci8jgwV0QuAx4GrvPrXwfcLCLL8DyN05PeaMPdl7LbCV8dtsHtEKHhiIEJSfcRFvrXfAFJ5OnFdSWJY11MKfPcBWegrosW2/NcBbGoKtooHKq6BDg8pvw5oGoRCVXdBnyslXuFU6uPNIaqxj7yN5cDqrva8ja/o4Iaq/9VRFNBhWgEg+aDL6+oODf40rOM2vuA9tjZYzx3wRntv4mAY11V+Z85PlIFIwndPHu8nVFRWaEi7DBmDDuMSTApsAlPY2jNMwytesJmjHcJllY957mqWqGZ7qhOJSpUl7ZJ+Kx5i5lz2oz2NJ6QeUtW92ReqmZnjSciIigSOR5a+Vhpv2/yIenfv4tZdelnUT90WYvl7ii3jV1TUUQEZ9SIe21WkXtZDNbR6DaGM1cjGhPQ6mS+QMhmzctuctjcRzo/1yb35LWLbiQgYh4HOfc4ZrxpEs92sKtqOIPVafbZtxJ6O3v+EgCuOfWtqdiQhF4UDVVl7E47Nq44vJvUPR31QkY6juNQbOBxpBZCaylHgB7wOIzGHkkWA+W3PNxECo0RQv+uE+jfdUJlYZJZ4yYUAEy65Ls1z0375s01zzkFhzfP+XFqdtjM8ZwLRxopRxotNdsuOp12pJvTlY804ruLGacAACAASURBVNbgUHEYtce+ia5XcShMOSxts3LN8188s+a56dfelmLKEctVBTnvqsqaNOdZpJXmpBHn376Uq09J76Xz/ftfACjZ/umZyV5+I53+8XvHrsUxas+pFSlHjOYRx+GAq+cC3iRASzmSPvYERhBpdFld8yfvpVYQLzV6NFrqxgdfjC03qukfv3dsVtyK/FRGQ4qDQxXH6rosm31a+yYBpuRNiMj1wAeAdap6qF+2G3ArMBVYDpymqhv8TOFXAycBW4GzVTWzrJa97U91mLx0RZ9/+9K2tBuISa/iKhWLOKVB/4TJpc0YHlIosP9Vt7T3HiI4hUKiLQE3ACdEyi4E7lXV6cC9/jHAicB0f5uFt4ZRZphwROjW8F5oPSw3js/f8WhT9a/8/bNc/YfnEtc/4/D8vwijDtpuY9scTZUy7tN/oPjk7yk++fusTUmNyZd+r+75TsweTyscV1V/R3Ui1/C6RNH1im7y1zlaiJcsdmJKv1LT9GRXVSeXic2SVjPiBtf8w08fK+1/7X8fnKptvYQCu+dMNEYiHQl0aS4cd4KILAodz/Gze9djL1VdA6Cqa0QkiKQorVfkE6xllMnqaOZxGABc/PPhp7MIR27d/FD+wnFr9bJ1k2gMrXyM4otLy+txJKT4+G/aY1AGdGqWeDxNRVWtV9WZoa2RaNS/cTWZhUr2pMfRzbQz422zHkiQILFd8zw6FSlmxOO86Wivm8rvfi0+em/Ft/K+t70vK9NaZuUlf53p/cVpe1TVWhGZ6HsbE/FWT4XyekUB4bWMOo55HEaJNLyOgF6IqnK1u7yNKM16HUY6tHkeR7AuEVSvV/RJ8TgK2BR0aWWBeRw9TCMPpM+8gZ6ncNC747upXJehB/+7HAziuqX9/nee2jkD84ZIauOnIvIj4Bi8sZCVwCXAFcA8ETkHWEF5qYk78UJxl+GF434qFSNaxIQjY5pZMbATXHLXk1z6/oPa0nbe1thIsABlx+mbfEjZ00iSriSMU6Dv0OMYeuTuUvdVHAN/mAdA/9GntWpmW1h16WezNsEjJeFQ1U/UOHVcTF0FzkvlxilgXVUjiKQ5rb56z1Mt36MX19nIO4WDjylt0GBso0vD0Vdfdi6QYrLClhFwnGRbD2Meh1HBcAfKox6F6ypOlw545PVfuzDlsBEzvrHma13zJdvD1hwHTDhaQl1t2+hv3tYZLyp89h37cf2iFbHnHUdwVSnkqIsqTxSXL6YwtT0LdQ38fi797z49cf2t87/p7bhuyXPZ8eMXtXTvdd88v7RoU5i4so4iAn392drQBZhwGC3jqnLuUVMb1nNEKLZPa0cufrdgW8TDf0Fv/80PGX3MX6bbdgzrv/2F0n50nka28zYqEX8ex0inbU9ARKaIyK9F5AkReUxEzvfLvyIiq0Rksb+dFLrmIhFZJiJPicj722Wb0Tyff/cBNc8FGXLrERaNoqt86OC90zBrxNKJtOrNzMR+4ydXxZa//qPL0jKnOxC8wfEkWw/TTo9jCPiCqj4kImOBB0XkHv/cVar6zXBlETkYOB04BNgH+KWIvElVu3O0LgVG8gS5BY+/ZFl0e5xuzvvWOtLzopCEtnkcqromSPurqluAJ/Byq9TiZGCuqm5X1efx4pWPbJd93UA3i8bl9z7N5fc+XVF2/tH7l/aLCbupoDrSKjwAX7QgrFQoPp88w/bQI3c31fb2X91U9/wbC77dVHv1mPC330qtrXZhCzl1KLBERKYChwP3+UWzRWSJiFwvIuP9slpJvKJtzRKRRSKyyB3c2karuwPt8EqBUPliv+LXz1Scm/3OaaUtKfXmbpjHkQJN5PNvVjTSCM9tprvq5asuGPb92oo43uB4kq2HabtwiMjOwHzg71R1M14e+QOAGXiZHYOvGImSeKnqnCBpmDOqe9NB5JXwXI9AQL7x22V847fLsjLJSEg9r6P46L0UH7039Xu+8bPvpN5mPTL/Ju+H4ybZepm2RlWJyCg80fihqv4YQFXXhs5/D/iZf9hVSbyaIW8htM3QrgSIRufopcy42SM9P7kvCW0TDn+pw+uAJ1T1ylD5xFByrlOAYEWhBcAtInIl3uD4dOD+dtlnNEfBEa7+w3MVQtJMd5XRW2y789rEcyrUddly01cAGPvJr7TPqE4QRFWNcNopnUcDZwLHRkJvvy4iS0VkCfBe4AIAVX0MmAc8DvwCOK+XI6p6ic8cuR+fnrkvn565b9ampMqeu+6UtQl1STJ3o9tW/wsEJEoeBsU9vCSHSbZepm0eh6r+D/HjFnfWueZy4PJ22WSky7ULlyeOrMoj6za9ngvxqDW24T79h9Tus/1XNzH62E+m1l6usa6q3KbrMTImSCHy3fsaT/4zOk+aomGEEAfp60+09TKWcsRomSDMNpg5/pkj98vSnBFLYdoRpX332fubCs/tJib87bdyEI6LeRyYcBjDxOZhGJ0k6ySHgvR8qG0STDiM1Lh+0QpG2bexzHCfW5S1CYnYdN3FAOx6TvN5rMRxSkkPpeB0PgGiRVUBNsaROxotxpQVzaRP6dbfwajP9l/+APBCcdNg45x/avqarD2OUq6qEZ7k0ITDyIQ8ZMftpiV9WyW1mdYpJyzsplTpTSGC9I1KtPUyJhxGR3FVOeXQiQB84C17ZWxNY17Z0vv50LIk7HWE1+ToasRJtvUwNsZhdJRTD9snaxMSk7dsK87+M72oKnFw3nRUqTwPKUeaEY1gbCMbr0V6XhSSYMJhpMqg68YOkBcVTntrtWiceFCl1/GrZS+3zbZWeWXLVnYfm4+Ems4B+VuJYMO1FzZVP+tuLjXhsK4qo7s49sA9sjbB6CBZi0DTCNZVhXkcRsoUIgPKrqs4NtkjM9LOVbXt7utSba8VMgnDLd8deiBoYrj0tiwamXHG4ZMBL0zXVY3tpqrFew6Y0C6zjGGQfSisR5ZeigJa6Eu09TImHEaqBGlIbn5o5bDaMfEwGiEFh32/9oMO31Ry3VUlIoeIyIdCx1f5K7FeLyJH1Ls2TG/LomEYXUUa3oIUHCZf+r0UrGnVgO4UhYRcAXwtdPx+4J+BHYEvAx9O0kiun4DRvUSHNeYtycVijjV5NYfzOXo5Q+6LF5/Dixefk8Gd8+1xABNV9Y+h482qOl9VbwYSu/ld+9sZ3UEa6UGcFgcT373/7sO+dxoE1udRPHqRsNey4qJPdf7+4iTaupSx4QNVPSp0uGfSRrr2tzOMbmTja7XFY9sbb7Bt6+ts2/o621/f0kGr8kWzaVCk0GWvqXx7HKtF5B3RQhE5CkjcLWBjHEZX866pntex8IVXM7YkfzPJu5lxs/61pSSHmSOS9wSG/wjcKiI3AMHSkW8HzgI+nrSRrpVFwwhz1H67ZW1CXd7Yti1rE3qafS6uzsibVVhunruqVPV+4CigAJztbw5wlH8uEW3zOERkCnATsDfgAnNU9WoR2Q24FZgKLAdOU9UN4qUivRo4CdgKnK2q8YspGyOSP9t3PA++uDFrM9iy9Q3AX5pBpOZ8sO2vb2H0TmPjT7aZbhoYD6+h0SxOwWGvf/j3lC0aDpL7FQBVdS1eBFXLJBIOERmtqtsblUUYAr6gqg+JyFjgQRG5B0/h7lXVK0TkQuBCPPfpRGC6v70DuNb/aRhdQ7Op1rdvKQvd6LHj0jYnF6jrIgWHnf/yy2z+QXPvq0A0Vl92bjtMa54g5UhOEZFf481jjENV9bgk7ST1OP4ERCeHxJWFLVgDrPH3t4jIE8Ak4GTgGL/ajcBv8ITjZOAmVVVgoYiME5GJfjuGAcDbp4zrCq8jERmt/e0uW5jJfRux5aavZG1CCuQ+O+4XY8qOAv4BWJe0kbrCISJ7473sx4jI4ZQjE3fBmzCSCBGZChwO3AfsFYiBqq4RkSAEbBLwYuiylX5ZhXCIyCxgFoDssEtSEwzDyJCgq2qXT3018TWpLUKVMurkN6ZIVR8M9kXkPXiT/0YDn1PVnydtp9ETeD9e19Jk4MpQ+RYgUUiEiOwMzAf+TlU313H1405UuVSqOgeYA9C3y0Rbg9QwjM4hufc4EJFgtvg24HJV/XWzbdR9Aqp6o6q+F2+g+r2h7UOq+uMEBo7CE40fhuqvFZGJ/vmJlN2jlcCU0OWTaSKu2DAMvIWcuphwd9W4Wf9as54UHPb84tUdsKgFRJJtDZuRE0TkKRFZ5o/3th0ReQD4LvAjvO6pTSJyRLAlbSeRz6Wq80XkfwOHADuEymv6nX6U1HXAE6oa9lYW4MUMX+H/vCNUPltE5uINim+y8Y3ewFXl9LdNytoMI0fsccFVWZtQg3Q8DhEpAN8B/gLvS/MDIrJAVR8fduP1eR14Dfiov4VR4NgkjSSNqvpPvDGN9wLf92/Y6KvN0cCZwFIRWeyX/ROeYMwTkXOAFcDH/HN34oXiLsMLx+18LgEjF7x9yjgWr8rJALlRk1xOACS1FQCPBJap6nMA/hfmk4G2CoeqHpNGO0lHed6lqm8VkSWqeqmIfAuo21Wlqv9D/LgFQFXIlx9NdV5CewzDiOA+tyhrE4aNFBwm/O23qsr3ufhaVl362ZrXdHQyYHLhmCAi4T/KHH+MFuKDgToy/cAPSDoPrwdJ8cTqO6qaTlRViDf8n1tFZB/gFWBaE7YaIxTX1dKiTobR6AUfJxpQnscxnMmEaaAIbs3vw1WsV9WZNc4lCgZKGxE5GrgFuAFvgrbgTau4X0T+UlUTzRxNKhw/E5FxwDfw8psoXpeVYdTFRMNohvXf/kJN8YD4VQg7KySKq6m837MKBvoW8GFVfThUdoeI3I43aJ7I60nkc6nqv6jqRlWdD+wHHKSq/9ysxYZhtIc8dlPVi6rqZjTh1oAHgOkiMk1E+oHT8QKE2s0uEdEAQFUXE0m5Xo9GEwA/UuccSUJyjd7FVaXQ4lobxshl8w++3NREwG5CgRSWqEFVh0RkNnAXXsLB61X1seG33BARkfGquiFSuBtNJL1t1FX1wchx8MjE3zfhMDJhxiQv79PSNZsytsRoN+E8VUnGOA64ei7LZp/WNns0na4qVPVOvGjSTnIVcLeIfJHKtOr/BvzfpI3UFQ5V/RSAiOwAnIqX0Ta4xmZtj3BaXdkvTQ6buOuIF4/i8w8lH64dATx3wRltazstjyMrVHWOiKwG/oXKqKrLVPWnSdtJOjj+E2AjnkIFCw/k+PEZvcRhE3flsTWbszbDaIHx514BeIPi7cIpOBRjBtVbQqGY8zefqv4M+Fm0XET+TlUTeR1J+7Qmq+rpqvp1Vf2Wv13Z+LL2MuNNNhu527nl4ZVZm9A1bN+c/iqGxeWLG1fqQjZdd3HNcy9fdUFsuRQcJl3y3YZtRyOvpl97Gwd9/yfNGVivfdVEWw75fNKKST2OP4rIYaq6tEWDDKOtHDJxFx5/qYu9joxSrHcjcSG1pXP++MXLV12AOzhUtYhTrUmAYYKsugdeM28YVtawD29Vuh4lcY9nUo/jz/EWYnpKRJaIyFIRWdKabemy4e5L2XD3pVmbYdSgk/M4Dt57Fw7ay9uM7qdeypFAQNZ87byqslpIwVsLvJ4wpYFqsi2HJLY6qcdxYouGdIxXf/Fldj/p8qzNMAwjAe2YtKfFYkX7y2af1havI8+D4yLyOlCMOwWMSdpO0gmAL8RtSW/SKV6580tZm2D4FDXbWeNv3rP7vA5Jubsqr+MbneSZc6MJYIeHKhRVE21dytOqukvMNlZVE69Qle8VSQzDyDXNZshtphsq8GraIR457qpKxbL8roFodD23PLwyc69j2ctbMru/kYwN13ZkDaNU8OZxdK8qJGBPEakZPZU0WrbnPA7rruoOCjYjLZaBTeuzNqGr6FSCwqdm1cye1DQp5arKigKwM15eqrgtEeZxGIaROYVRoygODralbTdlccrz4Diwpt7KrUnpOY/D6C5sAqDRa+R8jCOVvgDzOIye5sA9xvbsOEdh6oyeiqyKrsOx9ut/k5EltdHujphKQtXqq63QM8Kx2wlfRRxvAlAwEcgwjHzgFoulfFVadNnjgqtSa9spOEy/9jb4XjoDb3nuqlLVVPLetK2rSkSuF5F1IvJoqOwrIrJKRBb720mhcxeJyDJ/dvr7k9xj8dOrGHd87Zw3htGr3kavEYTZBoPl6755fqrtP/mZD6fSjpL7rqpUaOcYxw3ACTHlV6nqDH+7E0BEDsZbAesQ/5r/EBFzG3qEuY+sytoEYwST+uA4mmjrZdomHKr6OyCpW3QyMFdVt6vq88Ay4Mik9xr/vksq712Mm1FvZIHjWFyu0VuYx5FNVNVsP1Hi9SIy3i+bBLwYqrPSL2tIMK5hGEZvIwWn7nEnCCYAJtl6mU4/+WuBA4AZwBogCKOI+1oa++RFZJaILBKRRe7g1tibSKGAuuZ1GN3JwMZ1WZuQS6KTBTs1ebDingqDRU209TIdFQ5VXauqRVV1ge9R7o5aCUwJVZ0MrK7RxhxVnamqM51RO9a8l3kiRleScqLDwtQZqbZnNCJZgsOch+w2pKPCISITQ4enAEHE1QLgdBEZLSLTgOnA/Z20zWgv85bEfg9oK0+t6+KFnVKil+Zx5AHrqvJo2zwOEfkRcAwwQURWApcAx4jIDLznvxz4LICqPiYi8/AWTR8CzlNV62vqERyxAXKjR1DIoIes62ibcKjqJ2KKr6tT/3LAVmIyjCYoTJ1B8fmHsjajqzjg6rksm31aW9rugey4qdAzM8cNI8xI6KYyahOs/Jf6WhzAYJ6njqeECYdhdBJx6N91QtZWjAjSFg3A76oy4ch9dlx1i2y4+9KszTAMI2OkUGD/q25p6z2UZAPjvd6dlXvhMAwjBqeA86ajs7aiJylqsq2X6RnhePUXX87aBMMwehwLx/XoGeEAWza225m/tPNzOQwjVfwxjiRbL2OD40ZHsDXIPdIeGI+dACgOzgGJc4TmFik4TLmsZoR/W7CoKo+e8jiM7se8DiMtGonG9GtvS/2e1lXlYcJhdJzbH12TtQm9ywjxNpKSunio4rrJtl7GhMPoKEH6kQWPv5SxJfmnUZ4q9+k/dMiSkYNiUVVgwmFkQMEWdzKGQbNjG2+e8+NU729dVTkXjhlvmsTGX16WtRlGi7TL63hybe+nGym+uDRrE3LFQd//SSrteOtxuIm2XibXwhGHheQa3czApvXDbsNEIzusq8rDwnENwzCaoNe7oZLQcx6HYfQyQysfy9qETBHHYfKl38vs/tqhFQBF5GMi8piIuCIyM3LuIhFZJiJPicj7Q+Un+GXLROTCYRnQgJ4UjvU/beszM7qYPIxv2JrjrSEFh0mXfDdbIzo3c/xR4CPA78KFInIwcDpwCHAC8B8iUhCRAvAd4ETgYOATft220LNdVS/f8fcA7HnKlRlbYhghhrnmeN/kQxKNcfRSKK4UHCZe9J2szQD8MY4OzNFQ1ScApHr1zJOBuaq6HXheRJYBwcSdZar6nH/dXL/u4+2wr2eFA0w0jJGHs//MxpVyhFNw2Osf/j1rM0qowsBQYvGfICKLQsdzVHXOME2YBCwMHa/0ywBejJS/Y5j3qklPC4cxsnj8pc3YFJHeoptEA/wxjuQex3pVrankIvJLYO+YU19S1TtqXRZrVvywQ9tco54WjnW3f968jhHA0jWbcBCqvXrDSJkUVwBU1eNbuGwlMCV0PBkIEsDVKk+dtg2Oi8j1IrJORB4Nle0mIveIyDP+z/F+uYjIt/1ogCUickS77DJ6hyWrN7F0zaaszWiJdg2Qu88talzJaJlgjCPDtOoLgNNFZLSITAOmA/cDDwDTRWSaiPTjDaAvaJcR7YyqugFv1D/MhcC9qjoduNc/Bi8SYLq/zQKuTcuIdbd/Pq2mDCNzkoTjussWNqxjtIZ2KKpKRE4RkZXAO4H/FpG7vPvrY8A8vEHvXwDnqWpRVYeA2cBdwBPAPL9uW2hbV5Wq/k5EpkaKTwaO8fdvBH4D/KNffpOqKrBQRMaJyERVtTSqRiyLV20sJUwcKYz0ORzdQoeiqm4Hbq9x7nLg8pjyO4E722wa0Pl5HHsFYuD/3NMvn0R1RMAkYhCRWSKySEQWvfzyy4lu+tL881u32DDyxDDDfY36uKpsH3ITbb1Mt0wArBUpUF2oOkdVZ6rqzD322KPNZhlG+xjYkDzJ49CqJ9poidEMtnRs56Oq1gZdUCIyEQhGCOtFChg9yocOjotEHCGYZ5BLNMWoqjzTaY9jAXCWv38WcEeo/JN+dNVRwCYb3zBGAgOv2vejeqz9+t9kbUIVnchV1e20zeMQkR/hDYRP8KMDLgGuAOaJyDnACuBjfvU7gZOAZcBW4FPtssvIPw++uJFCt3SyDgfzOnJHkxMAe5Z2RlV9osap42LqKnBeu2zZ+9SrEafQruaNDvLgixuzNqHjDK1+KmsTup7nv3gmWiyWjve/6pa23KfJlCM9S0/PHDeMvDO05pmsTehqVlz0KYqDQ1Xlz55/OgdcPTf1+3kTAE04el449jzlShzzNowcYqJRnxUXxfdoq+ui7Vq6Va2rCkaAcBi9wwMrNvTkpL+B9SvpnzCZwXXLvQJVUDc2Rn2ks/KSv04sCstmn8aB18xL9f6dSqve7fS0cFiCw95g4QuvAlDoQdHA7/YYfHlFxoZ0L2u+dh5abM2LeObcj+Km6H2owpAJR9dMAGyJxU+vYtzxF2dthtEkRVdH9hwOIzFui4IBtKW7qguSHHYFPe1xpIUX9NWD33YNY4SgKQ1oq6pFVdGjwjHhg1cgTsFCcLsQV5VTD9sncf0/Ln+lJ8c1jPzS695EEnpCOMa/75IKoZBC+oLx1L9/mDf/zU9Sb9eo5FfLXi6NZTgivTHRz+gZLOWIR+6Fw7wKwzA6iZpw5F84OoV5G4ZhqIJrwpHvqCojn9z+aHz+yl8+k2x9FcPIDkU12dbLmHA0YNl/fASxwdnUsIFuI9coFIfcRFsvY11VNXBdxXGEA//Pj004DMMAvHkcltTYPI5Y1C3y/JyPZm3GiOKup9Y1rmQYXYB1VZnHUYW6RV64/gwO+Nz8rE3pSeLmcZhoGLnBBscBE45Yps26Dcex7inDMKKoheNiwlGFzQsxDKMWqlBsV8r2HGHCYRiG0QTmcZhwGB0mOr7x8yfXWoiukStMODKKqhKR5SKyVEQWi8giv2w3EblHRJ7xf47Pwjajvcxfurq0/7Mn1mZoiVEXxwIu41BVXDfZ1stk+el4r6rOUNWZ/vGFwL2qOh241z82DCMLbF3tmlg4bnfN4zgZuNHfvxH4cIa2GIZhxKJusq2XyUo4FLhbRB4UkVl+2V6qugbA/7ln3IUiMktEFonIIndwa4fMNdKkVq4qw+h21FKOANkNjh+tqqtFZE/gHhF5MumFqjoHmAPQt8vE3vYHe5BgIHzB4y8BPbqOuNG7qA2OQ0bCoaqr/Z/rROR24EhgrYhMVNU1IjIRsOnEPUrBJlcauUVxe3z8Igkd76oSkZ1EZGywD7wPeBRYAJzlVzsLuKPTthmdxcJwjbzhJTnURFsvk4XHsRdwu59xtg+4RVV/ISIPAPNE5BxgBfCxDGwzDMOojXVVARkIh6o+B7wtpvwV4LhO22MYhtEMvT5HIwk2c9wwDCMhqopruapMOAzDMJrBPI7umgBoGIbR9ahbTLQNBxH5hog8KSJLROR2ERkXOneRiCwTkadE5P2h8hP8smUi0tbMGyYchmEYSVHtiHAA9wCHqupbgaeBiwBE5GDgdOAQ4ATgP0SkICIF4DvAicDBwCf8um3BuqoMwzASomgaotD4Pqp3hw4XAsFa1icDc1V1O/C8iCzDmwcHsMwPPkJE5vp1H2+HfSYchmEYSVHFHRzo9F0/Ddzq70/CE5KAlX4ZwIuR8ne0yyATDsMwjKRoUx7HhGDZCJ85fsokAETkl8DeMdd9SVXv8Ot8CRgCfhhcFmcV8cMObRvFN+EwugpX4dgDJ/DbZ9dnbYphxNKEcKwPLRtR3Y7q8fUuFpGzgA8Ax2k5T/tKYEqo2mQgWOSmVnnq2OC40TUUXeX46XtkbYaRMuI4jP3kV7I2IxWCMY4ORFWdAPwj8CFVDacBXwCcLiKjRWQaMB24H3gAmC4i00SkH28AfcGwjKiDeRxG5riqnHjQXlmb0ZOI4+Ac9O6KsqFH7q5Rexg4BXb88AUAbP2vr1ec2ukTFwPUFY/13/4CAFJwmPC336o4t/qycxOZIIUC0755c1X5U7M+kuj6RGhTHsdwuAYYjZc9HGChqn5OVR8TkXl4g95DwHmqWgQQkdnAXUABuF5VH2uXcZLnlapmzpypixYtqiqf8MErEKeAOAUAHH8/OJbocaH6nONncBVHEJGKY4BCwak4FserIw74f2j/WEp1gqR+Tp9TOg9efXGoe88go2zB3w8f91Wdd6rqg5fCPHzcF9MWQL9vX9w9w205IhT8HtfA1lH+kqPh8op6Um4jeB4feEu1aPz22fV+3XK94J6eTV55kCfRwdsPEu8G5UHa9lI55esCd1tEcIRQW1T8/kGnsldHYtsCENXSCj5SWtHH//9St1zmH6NuaaW9ynNarh+uG24boFisbNvfdw4IgmyaY+BP8yF4KbouWiyGjotosCrgUHlweIeTkr3UuwERebBe11ESRo2fouOPuSBR3Zd/8oVh369b6UmPY/1Pq+e+TPzYNRlYYrSCjW/UQJz4peXEwZl2xLCb73/nqRXH23/zw/iKff3s8L5zhn2/PKKquEMdj6rqOnpSOOJY81+zS/uTPjGnTk2j3RRd5ZRDJ2ZtRn4Qh7593lxRVHxxaWm/MHVGW247+pi/BGD7r26qKB+pogF4UVXFjnRVdTUjRjjCrPrRrNL+lDNvyM6QEUgvCoYq7Dhmh47eszDlsI7da/Sxn+zYvfJAh8Y4upoRKRxhXrz57KqyabNu67gdRpn3HNCd4bjjdt4xaxOMrGluHkfPKZCRZAAABv1JREFUMuKFI47n53yUAz43P2szMufC906vKrvmT89nYEk2FF1lz113ytoMo6sw4QATjgr2+/QtVVFVRiWz3zmt5rnrF62oKnNd5fTDJ8XU7l5MLIxaeEvH2nocJhzAvmffXArHzTtFV7n5r5JF2PzDT9MN8/70zH1Tba9TqCoTx++ctRlGHrCoKsCEI3cUXeXOz70rlba+/sFDUmmnHbzngAn8/rlXOnKvybuZaBgJUcW1rqruEw5/qv3VeLMfv6+qV2RsUsf5zQXvydqEruDd++9ecfzAig2pta0Kk3c3wTCaQ8HCceky4QgtRvIXeMm8HhCRBaralpzy9XBdbds4h7rKwovr5jczYvizfceX9hev2tjUta7CgXuMTdskY6RhUVVAlwkH3oIkHVuMJGDFDWcmrnvQ+Xe00RIjKTMmeStpLl2zKVF9Ew0jHUw4oPuEYxIdXIykFZ68+uSsTTBCHDZx16xNMEYSNjgOdJ9w1FqkpFxBZBYQTP1+TUReAbpvtlglE+h+G8HsTJM82Agjy879hmuEvvHKXYOLfzAhYfU8PNeW6DbhqLdICQD+ClrhVbQWdXsGyjzYCGZnmuTBRjA7m0VVT8jahm6g2xZy6uhiJIZhGEbzdJXHoapDnVyMxDAMw2ierhIOAFW9E7iziUvykCM9DzaC2ZkmebARzE6jBXK9AqBhGIbRebptjMMwDMPockw4DMMwjKbIrXCIyAki8pSILBOR6kXGM0RElovIUhFZLCKL/LLdROQeEXnG/zm+UTttsOt6EVknIo+GymLtEo9v+893iYgMf1Hr1m38iois8p/nYhE5KXTuIt/Gp0Tk/Z2w0b/vFBH5tYg8ISKPicj5fnnXPM86NnbV8xSRHUTkfhF5xLfzUr98mojc5z/LW/1IS0RktH+8zD8/tRN2GiFUNXcbXsTVs8D+QD/wCHBw1naF7FsOTIiUfR240N+/EPi3DOz6X8ARwKON7AJOAn6ONynzKOC+DG38CvDFmLoH+3/70cA0/zNR6JCdE4Ej/P2xwNO+PV3zPOvY2FXP038mO/v7o4D7/Gc0DzjdL/9P4Fx///8A/+nvnw7c2om/uW3lLa8eRymnlaoOAEFOq27mZOBGf/9G4MOdNkBVfwe8GimuZdfJwE3qsRAYJyJtXzC8ho21OBmYq6rbVfV5YBneZ6PtqOoaVX3I398CPIGXMqdrnmcdG2uRyfP0n8lr/uEof1PgWCBYxzn6LINnfBtwnIjYymsdJK/CEZfTqpuWmVPgbhF50E+RArCXqq4B7x8a2DMz6yqpZVe3PePZfhfP9aFuvq6w0e8qORzvm3JXPs+IjdBlz1NECiKyGFgH3IPn7WxU1aEYW0p2+uc3AZU5+I22klfhaJjTKmOOVtUjgBOB80Tkf2VtUAt00zO+FjgAmAGsAb7ll2duo4jsDMwH/k5VN9erGlPWEVtjbOy656mqRVWdgZdm6EjgLXVsyfzvPtLJq3A0zGmVJaq62v+5Drgd7x9hbdA14f9cl52FFdSyq2uesaqu9V8sLvA9yt0nmdooIqPwXsg/VNUf+8Vd9TzjbOzW5+nbthH4Dd4YxzgRCSYph20p2emf35Xk3ZtGCuRVOLo2p5WI7CQiY4N94H3Ao3j2neVXOwvoloU9atm1APikHw10FLAp6ILpNJGxgFPwnid4Np7uR9lMA6YD93fIJgGuA55Q1StDp7rmedaysduep4jsISLj/P0xwPF44zG/Bj7qV4s+y+AZfxT4laqax9FJsh6db3XDi1J5Gq8v9EtZ2xOya3+8yJRHgMcC2/D6YO8FnvF/7paBbT/C65oYxPvWdk4tu/C6A77jP9+lwMwMbbzZt2EJ3ktjYqj+l3wbnwJO7OCz/HO87pElwGJ/O6mbnmcdG7vqeQJvBR727XkU+LJfvj+ecC0D/gsY7Zfv4B8v88/v36m/u23eZilHDMMwjKbIa1eVYRiGkREmHIZhGEZTmHAYhmEYTWHCYRiGYTSFCYdhGIbRFCYcRu4RkRtE5KONaxqGkQYmHIZhGEZTmHAYuUFEpvprS3zPX7fhbn+mcbjOcSLysHjroVwvIqP98uUicqmIPOSfOyib38Iw8o8Jh5E3pgPfUdVDgI3AqcEJEdkBuAH4uKoeBvQB54auXa9e8slrgS92zGLD6DFMOIy88byqLvb3HwSmhs692T//tH98I97CUAE/rnGdYRhNYMJh5I3tof0inlcR0Ggxn+Da6HWGYTSBCYfRSzwJTBWRA/3jM4HfZmiPYfQkJhxGz6Cq24BPAf8lIksBF2+tasMwUsSy4xqGYRhNYR6HYRiG0RQmHIZhGEZTmHAYhmEYTWHCYRiGYTSFCYdhGIbRFCYchmEYRlOYcBiGYRhN8f8BY0djm/AJ/K0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"dsp.TLONG.where(dsp.KMT > 0).plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:miniconda3-analysis]",
"language": "python",
"name": "conda-env-miniconda3-analysis-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
import numpy as np
import xarray as xr
def pop_add_cyclic(ds):
nj = ds.TLAT.shape[0]
ni = ds.TLONG.shape[1]
xL = int(ni/2 - 1)
xR = int(xL + ni)
tlon = ds.TLONG.data
tlat = ds.TLAT.data
tlon = np.where(np.greater_equal(tlon, min(tlon[:,0])), tlon-360., tlon)
lon = np.concatenate((tlon, tlon + 360.), 1)
lon = lon[:, xL:xR]
if ni == 320:
lon[367:-3, 0] = lon[367:-3, 0] + 360.
lon = lon - 360.
lon = np.hstack((lon, lon[:, 0:1] + 360.))
if ni == 320:
lon[367:, -1] = lon[367:, -1] - 360.
#-- trick cartopy into doing the right thing:
# it gets confused when the cyclic coords are identical
lon[:, 0] = lon[:, 0] - 1e-8
#-- periodicity
lat = np.concatenate((tlat, tlat), 1)
lat = lat[:, xL:xR]
lat = np.hstack((lat, lat[:,0:1]))
TLAT = xr.DataArray(lat, dims=('nlat', 'nlon'))
TLONG = xr.DataArray(lon, dims=('nlat', 'nlon'))
dso = xr.Dataset({'TLAT': TLAT, 'TLONG': TLONG})
# copy vars
varlist = [v for v in ds.data_vars if v not in ['TLAT', 'TLONG']]
for v in varlist:
v_dims = ds[v].dims
if not ('nlat' in v_dims and 'nlon' in v_dims):
dso[v] = ds[v]
else:
# determine and sort other dimensions
other_dims = set(v_dims) - {'nlat', 'nlon'}
other_dims = tuple([d for d in v_dims if d in other_dims])
lon_dim = ds[v].dims.index('nlon')
field = ds[v].data
field = np.concatenate((field, field), lon_dim)
field = field[..., :, xL:xR]
field = np.concatenate((field, field[..., :, 0:1]), lon_dim)
dso[v] = xr.DataArray(field, dims=other_dims+('nlat', 'nlon'),
attrs=ds[v].attrs)
# copy coords
for v, da in ds.coords.items():
if not ('nlat' in da.dims and 'nlon' in da.dims):
dso = dso.assign_coords(**{v: da})
return dso
@DapengLi17
Copy link

Hi @matt-long:
I am using your util.py to plot sea surface temperature on POP gx1v6 grids. There are strange horizontal lines at the north pole. Please see the attached figure. My jupyter-notebook is available at Could you please let me know how to get rid of the horizontal lines near the north pole. Thanks so much for your help.
Regards!
Dapeng
HMXL_POP_gx1v6_proj_2020Jun30

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment